Solving systems of linear equations

Solving a set of simultaneous linear equations is a fundamental problem that occurs in diverse applications. A linear system can be expressed as a matrix equation in which each matrix or vector element belongs to a field, typically the real numbers R. This section discusses how to solve a system of linear equations using a method called LUP decomposition.

We start with a set of linear equations in n unknowns x₁, x₂, . . . , x_n:

a x + a₁₂x₂ + + a_1nx_n = b₁,
a x + a₂₂x₂ + + a_2nx_n = b₂,

a_n x + a_n2x₂ + + a_nnx_n = b_n.

A set of values for x₁, x₂, . . . , x_n that satisfy all of the equations (17) simultaneously is said to be a solution to these equations. In this section, we only treat the case in which there are exactly n equations in n unknowns.

We can conveniently rewrite equations (17) as the matrix-vector equation

or, equivalently, letting A = (a_ij), x = (x_j), and b = (b_i), as

Ax = b .

If A is nonsingular, it possesses an inverse A^-1, and

x = A^-1 b

is the solution vector. We can prove that x is the unique solution to equation (18) as follows. If there are two solutions, x and x', then Ax = Ax' = b and

x (A^-1 A)x
A^-1(Ax)
A^-1(Ax')
(A^-1 A)x'
x'.

In this section, we shall be concerned predominantly with the case in which A is nonsingular or, equivalently (by Theorem 1), the rank of A is equal to the number n of unknowns. There are other possibilities, however, which merit a brief discussion. If the number of equations is less than the number n of unknowns--or, more generally, if the rank of A is less than n--then the system is underdetermined. An underdetermined system typically has infinitely many solutions (see Exercise 4-9), although it may have no solutions at all if the equations are inconsistent. If the number of equations exceeds the number n of unknowns, the system is overdetermined, and there may not exist any solutions. Finding good approximate solutions to overdetermined systems of linear equations is an important problem that is addressed in Section 6.

Let us return to our problem of solving the system Ax = b of n equations in n unknowns. One approach is to compute A^-1 and then multiply both sides by A^-1, yielding A^-1Ax = A^-1b, or x = A^-1b. This approach suffers in practice from numerical instability: round-off errors tend to accumulate unduly when floating-point number representations are used instead of ideal real numbers. There is, fortunately, another approach--LUP decomposition--that is numerically stable and has the further advantage of being about a factor of 3 faster.

Overview of LUP decomposition

The idea behind LUP decomposition is to find three n n matrices L, U, and P such that

PA = LU ,

where

L is a unit lower-triangular matrix,

U is an upper-triangular matrix, and

P is a permutation matrix.

We call matrices L, U, and P satisfying equation (20) an LUP decomposition of the matrix A. We shall show that every nonsingular matrix A possesses such a decomposition.

The advantage of computing an LUP decomposition for the matrix A is that linear systems can be solved more readily when they are triangular, as is the case for both matrices L and U. Having found an LUP decomposition for A, we can solve the equation (18) Ax = b by solving only triangular linear systems, as follows. Multiplying both sides of Ax = b by P yields the equivalent equation P Ax = Pb, which by Exercise 1-2 amounts to permuting the equations (17). Using our decomposition (20), we obtain

LUx = Pb .

We can now solve this equation by solving two triangular linear systems. Let us define y = Ux, where x is the desired solution vector. First, we solve the lower-triangular system

Ly = Pb

for the unknown vector y by a method called 'forward substitution.' Having solved for y, we then solve the upper-triangular system

Ux = y

for the unknown x by a method called 'back substitution.' The vector x is our solution to Ax = b, since the permutation matrix P is invertible (Exercise 1 -2):

Ax = P^- LUx
= P^-1 Ly
= P^-1 Pb
= b .

Our next step is to show how forward and back substitution work and then attack the problem of computing the LUP decomposition itself.

Forward and back substitution

Forward substitution can solve the lower-triangular system (21) in (n²) time, given L, P, and b. For convenience, we represent the permutation P compactly by an array [1 . . n]. For i = 1, 2, . . . , n, the entry [i] indicates that P_i,[i] = 1 and P_ij = 0 for j [i]. Thus, PA has a[i],j in row i and column j, and Pb has b[i] as its ith element. Since L is unit lower-triangular, equation (21) can be rewritten as

y = b_[1],
l y + y₂ = b_[2],
l y + l₃₂y₂ + y₃ = b_[3],

l_n y + l_n2y₂ + l_n3y₃ + + y_n = b_[n].

Quite apparently, we can solve for y₁ directly, since the first equation tells us that y₁ = b_[1]. Having solved for y₁, we can substitute it into the second equation, yielding

y = b_[2] - l₂₁y₁.

Now, we can substitute both y₁ and y₂ into the third equation, obtaining

y = b_[3] - (l₃₁y₁ + l₃₂y₂).

In general, we substitute y₁, y₂, . . . ,y_i-1 'forward' into the ith equation to solve for y_i:

Back substitution is similar to forward substitution. Given U and y, we solve the nth equation first and work backward to to the first equation. Like forward substitution, this process runs in (n²) time. Since U is upper-triangular, we can rewrite the system (22) as

u x + u₁₂x₂ + + u_1,n-2x_n-2 + u_1,n-1x_n-1 + u_1nx_n = y₁,
u x + + u_2,n-2x_n-2 + u_2,n-1x_n-1 + u_2nx_n = y₂,

un-2,n-2xn-2 + un-2,n-1xn-1 + un-2,nxn = yn-2,
un-1,n-1xn-1 + un-1,nxn = yn-1,
u_n,nx_n = y_n .

Thus, we can solve for x_n, x_n-1, . . . , x₁ successively as follows:

x_n = y_n/u_nn ,
x_n- (y_n-1 - u_n-1,nx_n)/u_n-1,n-1 ,
x_n (y_n-2 - (u_n-2,n-1x_n-1 + u_n-2,nx_n))/u_{n-2,n-2 ,}

or, in general,

Given P, L, U, and b, the procedure LUP-SOLVE solves for x by combining forward and back substitution. The pseudocode assumes that the dimension n appears in the attribute rows[L] and that the permutation matrix P is represented by the array .

Procedure LUP-SOLVE solves for y using forward substitution in lines 2-3, and then it solves for x using backward substitution in lines 4-5. Since there is an implicit loop in the summations within each of the for loops, the running time is (n²).

As an example of these methods, consider the system of linear equations defined by

where

and we wish to solve for the unknown x. The LUP decomposition is

(The reader can verify that PA = LU.) Using forward substitution, we solve Ly = Pb for y:

obtaining

by computing first y₁, then y₂, and finally y₃. Using back substitution, we solve Ux = y for x:

thereby obtaining the desired answer

by computing first x₃, then x₂, and finally x₁.

Computing an LU decomposition

We have now shown that if an LUP decomposition can be computed for a nonsingular matrix A, forward and back substitution can be used to solve the system Ax = b of linear equations. It remains to show how an LUP decomposition for A can be found efficiently. We start with the case in which A is an n n nonsingular matrix and P is absent (or, equivalently, P = I_n). In this case, we must find a factorization A = LU. We call the two matrices L and U an LU decomposition of A.

The process by which we perform LU decomposition is called Gaussian elimination. We start by subtracting multiples of the first equation from the other equations so that the first variable is removed from those equations. Then, we subtract multiples of the second equation from the third and subsequent equations so that now the first and second variables are removed from them. We continue this process until the system that is left has an upper-triangular form--in fact, it is the matrix U. The matrix L is made up of the row multipliers that cause variables to be eliminated.

Our algorithm to implement this strategy is recursive. We wish to construct an LU decomposition for an n n nonsingular matrix A. If n = 1, then we're done, since we can choose L = I₁ and U = A. For n > 1, we break A into four parts:

where v is a size-(n - 1) column vector, w^T is a size-(n - 1) row vector, and A' is an (n - 1) (n - 1) matrix. Then, using matrix algebra (verify the equations by simply multiplying through), we can factor A as

The 0's in the first and second matrices of the factorization are row and column vectors, respectively, of size n - 1. The term vw^T/a₁₁, formed by taking the outer product of v and w and dividing each element of the result by a₁₁, is an (n - 1) (n - 1) matrix, which conforms in size to the matrix A' from which it is subtracted. The resulting (n - 1) (n - 1) matrix

A' - vw^T/a₁₁

is called the Schur complement of A with respect to a₁₁.

We now recursively find an LU decomposition of the Schur complement. Let us say that

A' - vw^T/a₁₁ = L'U' ,

where L' is unit lower-triangular and U' is upper-triangular. Then, using matrix algebra, we have

thereby providing our LU decomposition. (Note that because L' is unit lower-triangular, so is L, and because U'is upper-triangular, so is U.)

Of course, if a₁₁ = 0, this method doesn't work, because it divides by 0. It also doesn't work if the upper leftmost entry of the Schur complement A' - vw^T/a₁₁ is 0, since we divide by it in the next step of the recursion. The elements by which we divide during LU decomposition are called pivots, and they occupy the diagonal elements of the matrix U. The reason we include a permutation matrix P during LUP decomposition is that it allows us to avoid dividing by zero elements. Using permutations to avoid division by 0 (or by small numbers) is called pivoting.

An important class of matrices for which LU decomposition always works correctly is the class of symmetric positive-definite matrices. Such matrices require no pivoting, and thus the recursive strategy outlined above can be employed without fear of dividing by 0. We shall prove this result, as well as several others, in Section 6.

Our code for LU decomposition of a matrix A follows the recursive strategy, except that an iteration loop replaces the recursion. (This transformation is a standard optimization for a 'tail-recursive' procedure--one whose last operation is a recursive call to itself.) It assumes that the dimension of A is kept in the attribute rows[A]. Since we know that the output matrix U has 0's below the diagonal, and since LU-SOLVE does not look at these entries, the code does not bother to fill them in. Likewise, because the output matrix L has 1's on its diagonal and 0's above the diagonal, these entries are not filled in either. Thus, the code computes only the 'significant' entries of L and U.

The outer for loop beginning in line 2 iterates once for each recursive step. Within this loop, the pivot is determined to be u_kk = a_kk in line 3. Within the for loop in lines 4-6 (which does not execute when k = n), the v and w^T vectors are used to update L and U. The elements of the v vector are determined in line 5, where v_i is stored in l_ik, and the elements of the w^T vector are determined in line 6, where w_i^T is stored in u_ki. Finally, the elements of the Schur complement are computed in lines 7-9 and stored back in the matrix A. Because line 9 is triply nested, LU-DECOMPOSITION runs in time (n³).

Figure 1 The operation of LU-DECOMPOSITION. (a) The matrix A. (b) The element a₁₁ = 2 in black is the pivot, the shaded column is v/a₁₁, and the shaded row is w^T. The elements of U computed thus far are above the horizontal line, and the elements of L are to the left of the vertical line. The Schur complement matrix A' - vw^T/a₁₁ occupies the lower right. (c) We now operate on the Schur complement matrix produced from part (b). The element a₂₂ = 4 in black is the pivot, and the shaded column and row are v/a₂₂ and w^T (in the partitioning of the Schur complement), respectively. Lines divide the matrix into the elements of U computed so far (above), the elements of L computed so far (left), and the new Schur complement (lower right). (d) The next step completes the factorization. (The element 3 in the new Schur complement becomes part of U when the recursion terminates.) (e) The factorization A = LU.

Figure 1 illustrates the operation of LU-DECOMPOSITION. It shows a standard optimization of the procedure in which the significant elements of L and U are stored 'in place' in the matrix A. That is, we can set up a correspondence between each element a_ij and either l_ij (if i > j) or u_ij (if i j) and update the matrix A so that it holds both L and U when the procedure terminates. The pseudocode for this optimization is obtained from the above pseudocode merely by replacing each reference to l or u by a; it is not difficult to verify that this transformation preserves correctness.

Computing an LUP decomposition

Generally, in solving a system of linear equations Ax = b, we must pivot on off-diagonal elements of A to avoid dividing by 0. Not only is division by 0 undesirable, so is division by any small value, even if A is nonsingular, because numerical instabilities can result in the computation. We therefore try to pivot on a large value.

The mathematics behind LUP decomposition is similar to that of LU decomposition. Recall that we are given an n n nonsingular matrix A and wish to find a permutation matrix P, a unit lower-triangular matrix L, and an upper-triangular matrix U such that PA = LU. Before we partition the matrix A, as we did for LU decomposition, we move a nonzero element, say a_k1, from the first column to the (1,1) position of the matrix. (If the first column contains only 0's, then A is singular, because its determinant is 0, by Theorems 4 and 5.) In order to preserve the set of equations, we exchange row 1 with row k, which is equivalent to multiplying A by a permutation matrix Q on the left (Exercise 1-2). Thus, we can write QA as

where v = (a₂₁, a₃₁, . . . , a_n1)^T, except that a₁₁ replaces a_k1; w^T = (a_k2, a_k3, . . . , a_kn); and A' is an (n - 1) (n - 1) matrix. Since a_k1 0, we can now perform much the same linear algebra as for LU decomposition, but now guaranteeing that we do not divide by 0:

The Schur complement A' - vw^T/a_k1 is nonsingular, because otherwise the second matrix in the last equation has determinant 0, and thus the determinant of matrix A is 0; but this means that A is singular, which contradicts our assumption that A is nonsingular. Consequently, we can inductively find an LUP decomposition for the Schur complement, with unit lower-triangular matrix L', upper-triangular matrix U', and permutation matrix P', such that

P'(A' - vw^T/a_k1) = L'U' .

Define

which is a permutation matrix, since it is the product of two permutation matrices (Exercise 1-2). We now have

yielding the LUP decomposition. Because L' is unit lower-triangular, so is L, and because U' is upper-triangular, so is U.

Notice that in this derivation, unlike the one for LU decomposition, both the column vector v/a_k1 and the Schur complement A' - vw^T/a_k1 must be multiplied by the permutation matrix P'.

Like LU-DECOMPOSITION, our pseudocode for LUP decomposition replaces the recursion with an iteration loop. As an improvement over a direct implementation of the recursion, we dynamically maintain the permutation matrix P as an array , where [i] = j means that the ith row of P contains a 1 in column j. We also implement the code to compute L and U 'in place' in the matrix A. Thus, when the procedure terminates,

LUP-DECOMPOSITION(A)
1 n rows[A]
2 for i 1 to n
do [i] i
4 for k 1 to n - 1
do p 0
6 for i k to n
do if |a_ik| > p
then p |a_ik|
k' i
if p = 0
then error 'singular matrix'
exchange [k] [k']
for i 1 to n
do exchange a_ki a_k'_i
for i k + 1 to n
do a_ik a_ik/a_kk
for j k + 1 to n
do a_ij a_ij - a_ika_kj

Figure 2 illustrates how LUP-DECOMPOSITION factors a matrix. The array is initialized by lines 2-3 to represent the identity permutation. The outer for loop beginning in line 4 implements the recursion. Each time through the outer loop, lines 5-9 determine the element a_k'k with largest absolute value of those in the current first column (column k) of the (n - k + l ) (n k + 1) matrix whose LU decomposition must be found.

Figure 2 The operation of LUP-DECOMPOSITION. (a) The input matrix A with the identity permutation of the rows on the left. The first step of the algorithm determines that the element 5 in black in the third row is the pivot for the first column. (b) Rows 1 and 3 are swapped and the permutation is updated. The shaded column and row represent v and w^T. (c) The vector v is replaced by v/5, and the the lower right of the matrix is updated with the Schur complement. Lines divide the matrix into three regions: elements of U (above), elements of L (left), and elements of the Schur complement (lower right).(d)-(f) The second step.(g)-(i) The third step finishes the algorithm. (j) The LUP decomposition PA = LU.

If all elements in the current first column are zero, lines 10-11 report that the matrix is singular. To pivot, we exchange [k'] with [k] in line 12 and exchange the kth and k'th rows of A in lines 13-14, thereby making the pivot element a_kk. (The entire rows are swapped because in the derivation of the method above, not only is A' - vw^T/a_k1 multiplied by P', but so is v/a_k1.) Finally, the Schur complement is computed by lines 15-18 in much the same way as it is computed by lines 4-9 of LU-DECOMPOSITION, except that here the operation is written to work 'in place.'

Because of its triply nested loop structure, the running time of LUP-DECOMPOSITION is (n³), the same as that of LU-DECOMPOSITION. Thus, pivoting costs us at most a constant factor in time.

Exercises

Solve the equation

by using forward substitution.

Find an LU decomposition of the matrix

Why does the for loop in line 4 of LUP-DECOMPOSITION run only up to n - 1, whereas the corresponding for loop in line 2 of LU-DECOMPOSITION runs all the way to n?

Solve the equation

by using an LUP decomposition.

Describe the LUP decomposition of a diagonal matrix.

Describe the LUP decomposition of a permutation matrix A, and prove that it is unique.

Show that for all n 1, there exist singular n n matrices that have LU decompositions.

Show how we can efficiently solve a set of equations of the form Ax = b over the boolean quasiring (,

Suppose that A is an m n real matrix of rank m, where m < n. Show how to find a size-n vector x₀ and an m (n - m) matrix B of rank n m such that every vector of the form x₀ + By, for y R^n-m, is a solution to the underdetermined equation Ax = b.