Egwald Mathematics: Linear Algebra

Systems of Linear Equations

by

Elmer G. Wiens

Egwald's popular web pages are provided without cost to users.
Please show your support by joining Egwald Web Services as a Facebook Fan:
Follow Elmer Wiens on Twitter:

linear equations | matrix representation | non-homogeneous system | homogeneous system
general solution | inconsistent system | echelon algorithm | singular value algorithm
overdetermined systems | least squares solution | normal equations | svd least squares
square A matrix | A inverse solution | P*A = L*U solution | A = Q*R solution | svd solution | references

Systems of Linear Equations.

A system of linear equations can be expressed as:

Example.

Matrix Representation.

The same problem expressed in matrix and vector form is:

where A = [a_{i, j}] is the matrix of coefficients, x^T = (x₁, x₂, . . . . x_n) is the vector of unknowns, and b^T = (b₁, b₁, . . . . b_m) is the vector of RHS (right hand side) coefficients.

If all the coefficients of the vector b are zero, the system A * x = b is called a homogeneous system.

Otherwise, the system is called nonhomogeneous.

Example: Nonhomogeneous System.

The columns (and rows) of the matrix A are linearly independent, as confirmed by the following diagram, where I use the standard (x, y, z) coordinate system.

Please note that I distinguish between the axis x, and the vector x^T = (x₁, x₂, x₃). Also, note that while the vector x^T is a row vector, and the vector x is a column vector, the diagrams below do not distinguish between column and row vectors.

Graphically, the matrix A, its column vectors, and the column vector b can be represented

Since the columns of A span R³, one can express any three dimensional vector, such as b in the example, as a linear combination of the columns of A. The solution to the linear equation problem, the vector

x^T = (x₁, x₂, x₃) = ( 1, 1, -1)

provides the coefficients of the linear combination of the columns of A that yield b.

Example: Homogeneous System.

Since the columns of A span R³, with the RHS vector b^T = (0, 0, 0) the only solution is x^T = (0, 0, 0), called the zero solution.

Any homogeneous system, (ie a homogeneous system with arbitrary matrix A and b^T = o^T the zero vector), always has at least one solution, the zero or trivial solution where x^T = o^T.

However, if the columns of A are not linearly independent, its homogeneous system has nontrivial solutions.

General Solution.

If x_p is a particular solution of the nonhomogeneous system, and x_h is a solution of the homogeneous system, then:

A * (x_p + x_h) = b + o = b,

proving that (x_p + x_h) is also a solution of the nonhomogeneous system.

The Image of a Matrix.

Let A be an m by n matrix. Write A = {a₁, a₂, ... a_n} as the set of m by 1 column vectors of A, ie A = [a₁, a₂, ... a_n].

The image of A, im(A), is the span of the columns of A. Formally:

im(A) = span(A) = {y in R^m | y = A * x,  for some x in Rⁿ.

Sometimes this is referred to as the column space of the matrix A, denoted by R(A).

The Kernel of a Matrix A.

The kernel of A, ker(A), is the set of all vectors x in Rⁿ for which A * x = o, the zero vector in R^m.  Formally:

ker(A) = {x in Rⁿ | A * x = 0  in R^m}

Sometimes this is referred to as the null space of the matrix A, denoted by N(A).

The matrix A, and its transpose matrix A^T (represented by A' in the diagrams below) generate four important subspaces: im(A), ker(A), im(A^T), and ker(A^T).

Example: General Solution with Linearly Dependent Columns (and Rows).

The matrix of the next linear equation problem has linearly dependent columns (and rows).

The following diagram displays this linear equation problem. The image of the matrix A, Im(A), is the two dimensional yellow subspace of R³ generated by the columns of A.

Because the vector b^T = (4, -6, -3) lies within the Im(A), it can be expressed as a linear combination of the columns of A. The particular solution to the linear equation problems is x_p^T = (1, -2, 0), since:

The ker(A) is the line (also a subspace) through the vector x_h^T = (0, -6, 4), since x_h is a solution of the homogeneous system for the matrix A:

Now let us look at A^T, the transpose matrix of the matrix A.

The image of the matrix A^T, Im(A^T), is the two dimensional yellow subspace of R³ generated by the rows of A (the columns of A^T).

As suggested by the diagram, ker(A) — the line through the vector x_h — is perpendicular to the Im(A^T) subspace.

Similarly, the kernel of the matrix A^T, ker(A^T) is the line through the vector y_h^T = (0, -3, 6), since y_h is a solution of the homogeneous system for the matrix A^T:

The subspace ker(A^T) is perpendicular to the subspace Im(A), as one can see two diagrams above.

Example: General Solution.

The general solution to the linear equation problem A * x = b is x = x_p + s * x_h for any scalar s. If x is defined by:

then

A * x = A * x_p + s * A * x_h = b + s * o = b.

Example: Inconsistent System.

Consider the system with the same matrix A as in the example above, but with a RHS vector b that lies outside of the Im(A).

The following diagram displays this linear equation problem. The image of the matrix A, Im(A), is the same two dimensional yellow subspace of R³ generated by the columns of A.

However, the vector b^T = (4, -6, 3) lies outside the Im(A) subspace. Therefore, it cannot be expressed as a linear combination of the columns of A. Thus, this linear equation problem has no particular solution, although its homogeneous system has solutions consisting of each vector on the line through the vector x_h^T = (0, -6, 4).

Solutions: Inconsistent System.

In the diagram above, I projected the vector b^T = (4, -6, 3) onto the Im(A) subspace, yielding the vector p^T = (4, -3.6, -1.8). Recall that the point of the vector p is the closest point on Im(A) to the point of the vector b. Consequently, the vector (b - p) is perpendicular to the vector p. To verify this, note that the dot product p^T * (b - p) = (4, -3.6, -1.8)^T * (0, -2.4, 4.8) = (0, 0, 0)^T.

A case can be made that the vector p^T = (4, -3.6, -1.8) provides a "solution" to the inconsistent system A * x = b, if one solves the system A * x = p, instead.

In least-squares problems, with more equations than unknowns, the solution provided by the algorithms is the vector x such that A * x = p. That is A * x equals the projection of the vector b onto Im(A).

Writing the matrix A = [a₁, a₂, a₃], one can see that a₃ = 1.5 * a₂, ie a₂ and a₃ are linearly dependent. Dropping the vector a₂ yields the least-squares problem B * z = b given by:

The least-squares solution is z^T = (0.8286, -0.8762) since

A least-squares problem always has a solution for B * z = p, since p is in the Im(B). Furthermore, if the columns of B are linearly independent, the solution vector z is unique.

To get back to the inconsistent system, write w^T = (0.82857, 0, -0.87619). Then, w is a solution of A * x = p since

However, the vector x = w + s * x_h, where s is a scalar and x_h is a solution of A * x = o, given by:

is also a solution to A * x = p,  because  A * (w + s * x_h) = A * w + s * A * x_h = p + s * o = p.

If small changes in the coefficients of the matrix A of a least-squares problem produce a matrix  with linearly dependent columns, the columns of A are "almost" linearly dependent. The estimates of the coefficients of the solution vector are tentative, since small changes in the data matrix A result in a system with a matrix  having many solutions.

Solution: Inconsistent System - Minimum Length.

One might choose as "optimal" the solution vector x^T = (x₁, x₂, x₃)^T = w + s * x_h with minimum length. The length (norm) of x = norm(w + s * x_h) is:

norm(x) = sqrt(x²₁ + x²₂ + x²₃) = sqrt[841/1225+36*s^2+(-92/105+4*s)^2],

having a minimum for s⁺ = 0.0674. The solution of minimum length to the inconsistent system is:

x⁺ = w + s⁺ * x_h = (0.8286, -0.4044, -0.6066).

The singular value decomposition algorithm provides x⁺ = A⁺ * b = (0.8286, -0.4044, -0.6066), where A⁺ is the pseudoinverse of A.

Echelon Algorithm.

The Echelon Algorithm is similar to the Gauss-Jordan method and the Gaussian method. The Echelon Algorithm permits one to calculate the general solution of a system of linear equations having more columns than rows. Begin by writing the matrix A and the vector b side by side as shown in Tableau⁰ below. This matrix arrangement is called the augmented matrix.

Using the following row operations on the tableaux of the augmented matrix, convert the left hand matrix to echelon form.

Begin in the upper left hand corner of the left hand matrix. Proceed along the diagonal converting each column into the appropriate column of the identity matrix as shown in the sequence of tableaux. Exchange rows with the lower row that has the largest entry (in absolute value) for the column. If all entries on and below the diagonal of the column are zeros, proceed to the next column.

The last tableau contains the information one needs to obtain the general solution to the original linear equation problem A * x = b.

The solution to the Echelon Tableau yields the particular solution, x_p, and the entire set of solutions of the associated homogeneous system, ie all vectors in the Ker(A).

Vectors in the Ker(A) are scalar multiples of the basis vectors, x_h1, x_h2, of the Ker(A) subspace.

Use the form below to solve some other linear equation problem.

Specify the dimensions of A (m <= n, 2 <= m,n) and fill in / change the numbers in the upper left hand corner, and then click "Solve".

Linear Equations Echelon Algorithm.

If you add a constant to an entry of the b vector, the system might be inconsistent. If you add a constant to an entry of a basis vector of Ker(A), you can create a matrix with more linearly independent columns.

Singular Value Decomposition Algorithm.

Consider the same linear equation problem A*x = b solved with the Echelon Algorithm.

Using singular value decomposition yields the orthogonal matrices U and V, and the diagonal matrix S of singular values of A:

A * x = U * S * V^T * x = b,

Since the matrix S has 2 positive singular values, the rank of the matrix is 2. Furthermore, the columns of U and V provide orthonormal bases for the following subspaces:

The Pseudoinverse is A⁺ = V * S⁺ * U^T

Both the columns and rows of A are linearly dependent. Setting x⁺ = A⁺ * b yields p = A * x⁺, the projection of b onto the Im(A) subspace.

x⁺ = A⁺ * b

p = A * x⁺

The RHS vector, b, lies inside the Im(A). Consequently, the system A * x = b is consistent. Therefore, p = b.

Overdetermined Systems
With # Rows m >= # Columns n

Suppose one wants to derive a linear equation model of the form A * x = b from a data set with a design (data) matrix A having more rows than columns. The columns of A and the column vector b represent measurable factors of the phenomena of interest. Each row of A represents one observation of the measured factors of the phenomena.

Since the matrix A has more rows than columns, one expects that A * x = b will be inconsistent, ie that b will lie outside the Im(A). Consequently, b will not be represented exactly as a linear combination of the columns of A.

I discuss least squares methods, also called multiple regression methods, on the Statistical web pages. Here I derive methods that statisticians and econometricians use to evaluate the least square (multiple regression) estimates for overdetermined systems of linear equations.

Least Squares (Multiple Regression) Solution.

With an inconsistent system, the least squares solution obtains a vector x⁺ from the data such that the vector p given by:

p = A * x⁺

is the projection of the vector b onto Im(A).

To be the projection of the vector b onto Im(A), p must be the vector in Im(A) that minimizes the distance between p and b, ie:

minimize sqrt{ (b₁ - p₁)² + (b₂ - p₂)² + .   .   . + (b_n - p_n)² }

giving the name "least squares" to this solution.

As the number of observations increases, one hopes that ones understanding of the phenomena improves. This could be the case if the least squares estimated values of each parameter of the vector x⁺ converge as the number of observations increases.

Furthermore, one might deduce that a model of the form:

z = a₁* x⁺ ₁ + a₂* x⁺ ₂ +   .   .   + a_n*x⁺ _n,

for a given instance of factors represented by the vector

a^T = (a₁, a₂, .   .   .   a_n)

could be used to predict the measured value of the factor represented in the b column, if the a vector was observed.

Normal Equations.

The normal equations provide a direct method for obtaining the least squares parameters:

x^+T = (x⁺₁, x⁺₂, .   .   .   x⁺_n)

to the linear equations system A * x = b.

To reduce the dimensions of the system, multiply the equation A * x = b by the matrix A^T, the transpose of A. The resulting equations are called the normal equations:

A^T * A * x = A^T * b.

Example.

As an experimenter, I want to summarize the data set with an equation in R³ space with standard coordinates (x, y, z):

z = x⁺₁ * x + x⁺₂ * y

described by the linear equations system A * x = b. The matrix A has two columns representing the independent variables (factors) of my model, x and y. The RHS vector b represents the dependent variable, z. I want to estimate the coefficients x⁺₁ and x⁺₂ using least squares. Perhaps, I want to use the equation to predict values of the dependent variable z for as yet unknown values of the independent variables x and y.

The normal equations system A^T * A * x = A^T * b is:

Multiplying by the matrix A^T yields the reduced system with 2 equations and 2 unknowns:

The matrix (A^T*A) is symmetric. If the data matrix A has linearly independent columns, (A^T*A) has an inverse. The normal equations provide a system with a square coefficient matrix, solvable by standard methods.

Proceeding with the inverse of (A^T*A):

The least squares solution using the normal equations is x⁺ = ( 0.6439, 0.5625 ).

The projection of b onto Im(A) is p = A * x⁺. Comparing p and b and computing b - p:

The distance from Im(A) to b — the length of the vector (b-p) — is 4.134, found by computing the square root of the sums of the squares of the (b-p) column. (Im(A) is a subspace of R⁸, while b and p are vectors in R⁸)

For each row i of the data matrix A = [a₁, a₂], the triples (a_i,1, a_i,2, b_i) and (a_i,1, a_i,2, p_i) are points in R³ with standard coordinates (x, y, z). Furthermore, the least squares solution can be represented by the plane:

z = 0.6439 * x + 0.5625 * y

The plane is displayed below. If the point (a_i,1, a_i,2, b_i) lies above the plane, a blue line joins it to the point (a_i,1, a_i,2, p_i) on the plane. If (a_i,1, a_i,2, b_i) lies below the plane, a red line joins it to (a_i,1, a_i,2, p_i).

By clicking the refresh button on your browser, you can generate a new system of linear equations, with the new solution with least squares parameters x⁺ = (x⁺₁, x⁺₂) displayed above.

The functional form,

z = x⁺₁ * x + x⁺₂ * y,

ensures that the solution plane passes through the origin.

In the example of the next section, the estimated solution plane is not restricted to pass through the origin.

Singular Value Decomposition and Least Squares.

The "normal equations" least squares method of solving systems of linear equations computes the matrix (A^T * A). For systems with many observations, the data matrix A has many rows. Multiplying two matrices using a computer incurs the risk of roundoff errors distorting the components of the matrix (A^T * A) from their true values. Furthermore, the matrix (A^T * A) may be ill-conditioned if the columns of A are nearly linearly dependent. Experimenters often include as many factors as possible in their design matrix A, hoping to improve the "explanatory power" of their model. The probability that subtle inter dependencies among the columns of A will render A nearly linearly dependent increases with the number of included factors.

The singular value decomposition method controls for such errors. Firstly, it controls roundoff errors by decomposing the data matrix A using Householder and Jacobi transformations, known for their stability. This decomposition is obtained without computing (A^T * A). Secondly, singular value decomposition controls for linear dependence (and roundoff error) by setting very small singular values to zero.

Example.

The matrix A has two columns representing the independent variables (factors) of my model, x and y. The RHS vector b represents the dependent variable, z. As an experimenter, I want to summarize the data set with an equation in R³ space with standard coordinates (x, y, z):

z = x⁺₁ + x⁺₂ * x + x⁺₃ * y

described by the linear equations system A * x = b. I need to estimate three coefficients x⁺₁, x⁺₂, and x⁺₃ since I do not require the solution plane to pass through the origin. To obtain a z-intercept coefficient given by x⁺₁, I append a column of ones to the front of the matrix A

Singular value decomposition yields the orthogonal matrices U and V, and the diagonal matrix S of singular values of A:

A * x = U * S * V^T * x = b,

Since the matrix S has 3 positive singular values, the rank of the matrix A is 3. Computing A⁺ = V * S⁺ * U^T, the least squares solution to A * x = b is:

x⁺ = A⁺ * b.

The least squares solution using singular value decomposition is x⁺ = ( 1.8413, 0.4334, 0.4609 ).

The projection of b onto Im(A) is p = A * x⁺. Comparing p and b and computing b - p:

The distance from Im(A) to b — the length of the vector (b-p) — is 1.677, found by computing the square root of the sums of the squares of the (b-p) column. (Im(A) is a subspace of R⁸, while b and p are vectors in R⁸)

For each row i of the data matrix A = [a₁, a₂, a₃], the triples (a_i,2, a_i,3, b_i) and (a_i,2, a_i,3, p_i) are points in R³ with standard coordinates (x, y, z). Furthermore, the least squares solution can be represented by the plane:

z = 1.8413 + 0.4334 * x + 0.4609 * y

The plane is displayed below. If the point (a_i,2, a_i,3, b_i) lies above the plane, a blue line joins it to the point (a_i,2, a_i,3, p_i) on the plane. If (a_i,2, a_i,3, b_i) lies below the plane, a red line joins it to (a_i,2, a_i,3, p_i).

By clicking the refresh button on your browser, you can generate a new system of linear equations, with the new solution with least squares parameters x⁺ = (x⁺₁, x⁺₂, x⁺₃) displayed above.

A Square Matrix.

If A is a square matrix of dimension n, the system of linear equations problem has of n equations in n unknowns. If the rank of A, rank(A), equals n, its columns and rows are linearly independent and the problem has a unique solution.

The variables x = [x_i] that solve the system of linear equations A * x = b can be obtained by various methods.

You can adjust the matrix A and the RHS vector b in the form below to see how the different methods solve the problem. The available methods are based on the techniques for decomposing a matrix that I discussed on the matrices web page. These methods are the A inverse solution using the Gauss-Jordan method,   P*A = L*U solution using the Gaussian method, A = Q*R solution using Householder transformations, and the svd solution using singular value decomposition.

The Gauss-Jordan Inverse Method

The Gausss-Jordan inverse method uses row operations on the initial tableau [ A  |  I ] to convert the RHS matrix to the identity matrix I and the LHS matrix to the inverse matrix A^-1.

Solve A * x = b with A^-1.

The Gaussian Decomposition Method

The Gaussian decomposition method factors the matrix A into the product of a lower triangular matrix L and an upper triangular matrix U, keeping track of row switches to construct a permutation matrix P so that P * A = L * U.

Solve A * x = b using P * A * x = L * U * x = P * b.

The solution vector x is obtained in two steps. First, solve for y in L * y = P * b by forward substitution:

Next, solve for x in U * x =y by backward substitution:

The Householder QR Method

The Householder decomposition method factors the matrix A into the product of an orthonormal matrix Q and an upper triangular matrix R, such that A = Q * R.

Solve A * x = b using A * x = Q * R * x = b.

Multiplying by Q^T gives R * x = Q^T * b.

These equations are solved for x by back substitution yielding:

The Singular Value Decomposition Method

The singular value decomposition method factors the matrix A into the product A = U * S* V^T, where U consists of the eigenvectors of A * A^T, V consists of the eigenvectors of A^T * A, and S consists of the square roots of the eigenvalues of A * A^T (which equal the eigenvalues of A^T * A).

Solve A * x = b using A * x = U * S * V^T * x = b.

Since the matrix S has 5 positive singular values, the rank of the matrix A is 5. Computing A⁺ = V * S⁺ * U^T, the solution to A * x = b is:

References.

• Ayres, Frank Jr. Matrices. New York: Schaum McGraw-Hill, 1962.
• Ayres, Frank Jr. Modern Algebra. New York: Schaum McGraw-Hill 1965.
• Bretscher, Otto. Linear Algebra with Applications. Upper Saddle River: Prentice Hall, 1997.
• Burden, Richard L. and J. Douglas Faires. Numerical Analysis. 6th ed. Pacific Grove: Brooks/Cole, 1997.
• Cohn, P. M. Linear Equations. London: Routledge, 1964.
• Demmel, James W. Applied Numerical Linear Algebra. Philadelphia: Siam, 1997.
• Dowling, Edward T. Mathematics for Economists. New York: Schaum McGraw-Hill, 1980.
• Lipschutz, Seymour. Linear Algebra. New York: Schaum McGraw-Hill, 1968.
• Mathews, John H. and Kurtis D. Fink. Numerical Methods Using MATLAB. 3rd ed. Upper Saddle River: Prentice Hall, 1999.
• Press, William H., Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling. Numerical Recipes: The Art of Scientific Computing. Cambridge: Cambridge UP, 1989.
• Strang, Gilbert. Linear Algebra and Its Applications. 3d ed. San Diego: Harcourt, 1976.
• Varah, James. Numerical Linear Algebra: Computer Science 402 Lecture Notes. Vancouver: University of B.C., 2000.
• Watkins, David S. Fundamentals of Matrix Computations. New York: John Wiley, 1991.