EE5120 Linear Algebra - Linear system of equations

Take a single equation, say \(5x + 6y + 3z = 7\). You know from school that this is a plane in three dimensions. Let us derive that fact properly, because the derivation generalizes to \(n\) dimensions without any change.

A plane is fixed by two pieces of data: the direction perpendicular to it, given by a unit normal \(\hat{n}\), and how far it sits from the origin, given by a distance \(d\). Drop a perpendicular from the origin \(O\) onto the plane and call the foot \(N\). Then \(\vec{ON} = d\,\hat{n}\). Now take any point \(P\) on the plane with position vector \(\vec{r}\). The segment \(\vec{NP} = \vec{r} - d\,\hat{n}\) lies entirely within the plane, so it must be perpendicular to \(\hat{n}\).

Writing that perpendicularity as a dot product,

For \(d \neq 0\) this gives \(\vec{r}\cdot\hat{n} = d\). (When \(d=0\) the plane passes through the origin and \(\vec{r}\cdot\hat{n}=0\) holds directly, so the relation \(\vec{r}\cdot\hat{n}=d\) covers every case.) With \(\vec{r}=(x,y,z)\) and \(\hat{n}=(n_x,n_y,n_z)\) this is

Two points about this equation. First, an equation and any nonzero multiple of it describe the same plane, so the coefficients by themselves do not fix \(\hat{n}\) and \(d\) — you must normalize before you read the right hand side as a distance. For \(5x+6y+3z=7\) we have \(\hat{n} = (5,6,3)/\sqrt{70}\) and \(d = 7/\sqrt{70}\). Second, if you want \(d\) to be a genuine (non-negative) distance and the arithmetic hands you \(d<0\), simply absorb the sign into \(\hat{n}\); the normal then points the other way.

Nothing in the argument used the number three. In \(\mathbb{R}^n\), with \(\vec{r}=(x_1,\ldots,x_n)\) and \(\hat{n}=(a_1,\ldots,a_n)\), the same steps give

and the solution set is called a hyperplane: a flat object of dimension \(n-1\) sitting inside \(n\) dimensions. One scalar equation removes exactly one degree of freedom.

With this, \(Ax=b\) with \(A\) of size \(n\times n\) is a set of \(n\) hyperplanes in \(\mathbb{R}^n\), row \(i\) giving the hyperplane \(\sum_j a_{ij}x_j = b_i\). A solution is a point lying on all of them at once. So solving the system means intersecting \(n\) hyperplanes.

In two dimensions this is two lines in a plane, and there are only three things that can happen: they cross at one point, they are parallel and distinct, or they are the same line.

In three dimensions the arrangements are richer — three planes can meet at a point, meet along a common line, form a triangular prism with no common point, be parallel, and so on — but the outcomes are still only three in number, and there is a reason for that.

Claim 1.1. For a system \(Ax=b\) over \(\mathbb{R}\), the number of solutions is either zero, exactly one, or infinite. It can never be, say, exactly two.

The cases with no solution or with infinitely many are lumped together and called the singular case. The word will be given a sharper meaning shortly, once we have elimination and pivots in hand.

A remark on why the case of no solution is the common one in practice. In an experiment you typically measure more quantities than you have unknowns, so \(A\) is tall (\(m > n\)), and each measurement is corrupted by noise. The hyperplanes then fail to have a common point by a small amount. Insisting on an exact solution is the wrong thing to ask for; what one wants is the point that misses all of them by as little as possible. That is the least squares problem, and we do it properly in the chapter on orthogonality.

Take the system

Subtract \(2\times\) row 1 from row 2, and \((-1)\times\) row 1 from row 3, so as to create zeros below the leading \(2\). Then subtract \((-1)\times\) the new row 2 from the new row 3. Carrying the right hand side along, the augmented array goes

The entries \(2, -8, 1\) that sit at the head of each row at the moment it is used are the pivots, and this stage is called forward elimination. The matrix is now upper triangular, and the system can be read off from the bottom up: \(w=2\), then \(-8v-2(2)=-12\) gives \(v=1\), then \(2u+1+2=5\) gives \(u=1\). This stage is back substitution, and the solution is \((1,1,2)\).

In general, at step \(k\) we use row \(k\) to clear column \(k\) in all the rows below it. Denoting by \(a^{(k)}_{ij}\) the entries at the start of step \(k\), the algorithm is

for \(k = 1,\ldots,n-1\). The numbers \(\ell_{ik}\) are called the multipliers, and \(a^{(k)}_{kk}\) is the \(k\)-th pivot. Keep the multipliers in mind — they are about to reappear, arranged in a matrix.

The algorithm divides by the pivot, so a zero in the pivot position stops it. Two rather different things can be happening.

The first is a nuisance rather than a real problem. In

the second pivot is zero, but there is a nonzero entry below it in the same column. Exchanging rows 2 and 3 restores the elimination, and the system is perfectly well behaved. The zero was an artefact of the order in which the equations happened to be written down.

The second is genuine trouble. If in the same column every entry from the pivot position downwards is zero, no exchange can help. Then the matrix is singular in the sense of the previous section, and depending on the right hand side we get either no solution or infinitely many. Let us see this happen inside the elimination.

The two look alike in the middle of a calculation, so we state the test that separates them carefully.

For a \(3\times 3\) system nobody cares about cost. For \(n = 10^6\), which is a routine size in engineering computation, cost is the only thing one cares about. Let us count.

At step \(k\) there are \(n-k\) rows below the pivot row. For each of them we compute one multiplier (a division), and then update the entries in columns \(k+1\) through \(n\) — that is \(n-k\) entries, each needing one multiplication and one subtraction. The entry in column \(k\) itself is set to zero and costs nothing. So step \(k\) costs \((n-k)\) divisions and \((n-k)^2\) multiply-subtract pairs, and the total over all steps is

multiply-subtract pairs, plus \(n(n-1)/2 \approx n^2/2\) divisions. The leading behaviour is what matters: elimination on the matrix costs \(O(n^3)\), and more precisely about \(n^3/3\) operations.

The right hand side is far cheaper. Carrying \(b\) through forward elimination costs \(n-k\) operations at step \(k\), i.e. \(n(n-1)/2\) in all, and back substitution costs about the same. So

Is \(n^3\) the end of the story? Ordinary matrix multiplication is also an \(O(n^3)\) operation, and Strassen showed in 1969 that it can be done in \(O(n^{\log_2 7}) = O(n^{2.81})\); later theoretical constructions have pushed the exponent below \(2.38\), though those algorithms have enormous constants and are never used in practice. Solving a linear system can be made to run at the same exponent as matrix multiplication. In practice, for dense matrices, \(n^3/3\) is what you get. The real savings come from a different direction — from the structure of \(A\), which is the subject of the last section of this chapter.

Start from \(IA = A\) and modify \(I\) a little. Define, for \(i \neq k\),

which is the identity with the single extra entry \(-\ell\) in position \((i,k)\). Then

and since \(e_i\) times a row vector places that row in position \(i\) and zeros elsewhere, the effect is precisely \(\text{row}_i \leftarrow \text{row}_i - \ell\,\text{row}_k\). These are called elementary or elimination matrices. Note that \(i>k\) for our purposes (we always subtract an earlier row from a later one), so \(E_{ik}\) is lower triangular with ones on the diagonal — unit lower triangular, in the standard phrase.

Elimination on a \(3\times3\) matrix uses \(E_{21}\), \(E_{31}\) and then \(E_{32}\). It is natural to ask whether the order matters. The outer product form answers this in one line. Using \((e_a e_b^T)(e_c e_d^T) = (e_b^Te_c)\, e_a e_d^T\) and \(e_b^Te_c = \delta_{bc}\),

So the two commute exactly when \(\delta_{kj} = \delta_{mi} = 0\), i.e. when neither operation touches a row that the other one uses. \(E_{21}\) and \(E_{31}\) commute, because both of them modify rows using row 1 and neither of them modifies row 1. But \(E_{32}\) and \(E_{21}\) do not commute, because \(E_{21}\) changes row 2 which \(E_{32}\) then uses. This is exactly the intuition you would have from doing elimination by hand, now made precise.

The inverse of \(E_{ik}(\ell)\) must undo the row operation, i.e. it must add back what was subtracted. So we expect \(E_{ik}(\ell)^{-1} = E_{ik}(-\ell) = I + \ell\,e_ie_k^T\). The algebra confirms it:

since \(e_k^Te_i = 0\) for \(i \neq k\). So to invert an elementary matrix we only flip the sign of its off-diagonal entry.

Two facts we will keep using.

Lemma 1.3. The product of two lower triangular matrices is lower triangular. If both have unit diagonal, so does the product.

Lemma 1.4. A lower triangular matrix with nonzero diagonal entries is invertible, and its inverse is lower triangular, with unit diagonal if the original had one.

Taken together, the unit lower triangular matrices are closed under products and inverses. Everything we are about to do stays inside this family.

Given \(A = LU\), the system \(Ax=b\) becomes \(L(Ux)=b\), and is solved in two triangular sweeps:

Each sweep costs \(O(n^2)\). Now recall the cost gap noted earlier. The factorization is done once, at \(\sim n^3/3\). After that, every new right hand side costs only \(\sim n^2\). In engineering this is the common situation: the same physical system (same \(A\)) is driven by many different sources (many \(b\)). Factor once, solve many times.

The pivots sit on the diagonal of \(U\). Pulling them out into a diagonal matrix \(D = \mathrm{diag}(d_1,\ldots,d_n)\) and dividing each row of \(U\) by its pivot leaves an upper triangular matrix with ones on the diagonal:

This gives \(A = LDU\), now with both triangular factors carrying unit diagonals and all the scaling collected in \(D\). The symmetry of the notation is the point: \(L\) and \(U\) are on the same footing, which is what makes the next result possible.

Theorem 1.5. Let \(A\) be nonsingular and suppose \(A = L_1D_1U_1 = L_2D_2U_2\), where the \(L_i\) are unit lower triangular, the \(U_i\) are unit upper triangular and the \(D_i\) are diagonal with no zero on the diagonal. Then \(L_1=L_2\), \(D_1=D_2\) and \(U_1=U_2\).

The book states this and moves on. Here is the proof; it is short, and we will need it in a moment for symmetric matrices.

We assumed no row exchanges were needed. That assumption is not free, and it can be characterized exactly. Write \(A_k\) for the leading principal submatrix of \(A\), i.e. the top-left \(k\times k\) block.

Theorem 1.6. Gaussian elimination on \(A\) produces \(n\) nonzero pivots without any row exchanges if and only if \(A_k\) is nonsingular for every \(k=1,\ldots,n\). In that case \(A=LU\) exists and, by Theorem 1.5, is unique.

Using determinants, which we will bring in as and when they are needed, this can be sharpened into a formula for the \(k\)-th pivot:

which also shows immediately that \(\det A = d_1 d_2\cdots d_n\).

Row exchanges are also matrix multiplications. The matrix \(P_{ij}\) obtained by swapping rows \(i\) and \(j\) of the identity does the job:

More generally, a permutation matrix \(P\) has exactly one 1 in each row and each column and zeros elsewhere; its rows are \(e_{\sigma(1)}^T,\ldots,e_{\sigma(n)}^T\) for some rearrangement \(\sigma\) of \(1,\ldots,n\). There are \(n!\) of them, which for \(n=3\) gives the six matrices you can write down by hand.

Permutation matrices have a pleasant property that we use again in the chapter on orthogonality.

Claim 1.7. \(P^TP = I\), i.e. \(P^{-1}=P^T\).

The row exchanges are discovered in the course of elimination, not by inspecting \(A\) beforehand. But once elimination has finished we know which exchanges were made, and it is a fact — which we state here without proof — that they can be collected in hindsight into a single permutation applied at the start.

Theorem 1.8. For any nonsingular \(A\) there is a permutation matrix \(P\) such that \(PA = LU\).

To use it: from \(Ax=b\), form \(PAx = Pb\), then factor \(PA = LU\) and do the two triangular solves against the reordered right hand side \(Pb\). Note that \(x\) itself is not permuted — we have only reordered the equations, not the unknowns.

So far pivoting was forced on us by exact zeros. In floating point arithmetic we must pivot for a second, subtler reason: a pivot that is merely small is almost as bad as one that is zero.

Suppose our machine keeps only three significant digits. Then \(0.123400 + 0.005567 = 0.128967\) is stored as \(0.129\); the tail is simply lost. Note that it is significant digits and not decimal places that a floating point machine keeps — that is what lets it hold \(10^{-4}\) and \(10^{4}\) with equal relative accuracy, and it is the reason the failure below is a failure of relative precision. Now take

whose exact solution is \(u = 1/0.9999 = 1.00010\ldots\) and \(v = 0.99989\ldots\); to three significant digits, \((1.00,\,1.00)\).

Eliminate without pivoting. The multiplier is \(1/0.0001 = 10^4\), and row 2 becomes \((0 \quad 1-10^4)\) with right hand side \(2-10^4\). Exactly, that is \(-9999\,v = -9998\); but our machine keeps three digits, so it stores both numbers as \(-1.00\times 10^4\) and computes \(v = 1.00\). Back substitution then gives

The computed answer is \((0,1)\) instead of \((1,1)\), which is not a small error. What went wrong is that the huge multiplier \(10^4\) swamped the original row 2, so everything that equation had to say about \(u\) was rounded away before it could be used; and then \(u\) was recovered by dividing the surviving rounding error by the tiny pivot, which blew it up.

Now exchange the two rows first, so that the pivot is \(1\) and the multiplier is \(0.0001\). Row 2 becomes \((1-0.0001)v = 1-0.0002\), i.e. \(0.9999\,v = 0.9998\), which the machine stores as \(1.00\,v = 1.00\). So \(v = 1.00\) and \(u = 2 - 1.00 = 1.00\). Correct to the working precision.

The escape is a narrow one. Had the machine kept four digits instead of three, \(-9999\) and \(-9998\) would have been stored exactly, \(v\) would have come out as \(0.9999\), and back substitution would have given \(u=1\). The unpivoted elimination fails not because the arithmetic is bad but because the algorithm arranged for the answer to depend on digits the machine had already thrown away.

This is partial pivoting: at step \(k\), search column \(k\) from the diagonal downwards, and exchange rows so that the entry of largest magnitude becomes the pivot. The reason it works is one line, and we state it explicitly.

Two caveats. First, "at most double" applied over \(n-1\) steps only bounds the growth by \(2^{n-1}\), and that bound is attained — Wilkinson constructed matrices for which the entries do grow that fast. So partial pivoting does not come with a good worst-case guarantee. What saves it is that such matrices essentially never arise in practice; observed growth is small, and partial pivoting is the standard choice in every serious library. Second, partial pivoting compares entries in a column without regard to the scale of the rows, and that scale is ours to choose.

For a square matrix \(A\), an inverse is a matrix \(A^{-1}\) with \(AA^{-1} = A^{-1}A = I\). Two immediate points.

Uniqueness. If \(B\) and \(C\) are both inverses of \(A\), then \(B = B(AC) = (BA)C = C\). So "the" inverse is justified: if it exists, there is only one.

A test for non-existence. If \(Ax=0\) for some \(x \neq 0\), then \(A\) cannot have an inverse, because \(A^{-1}Ax = A^{-1}0\) would give \(x=0\), a contradiction.

The converse is also true, and elimination is what proves it. Collecting the threads of this chapter:

Theorem 1.9. For a square matrix \(A\) the following are equivalent.

Note in passing that for a square matrix, a one-sided inverse is automatically two-sided; for rectangular matrices this fails, and one-sided inverses become an interesting topic in their own right (Chapter 2).

For invertible \(A,B\) of the same size:

For \(2\times2\) there is a closed form to remember:

The inverse is computed by the Gauss-Jordan method, which is elimination applied to the \(n\) systems \(Ax = e_j\) simultaneously: augment \(A\) with \(I\), eliminate downwards, then eliminate upwards, then divide each row by its pivot, and \(I \,|\, A^{-1}\) is left behind. The cost is again \(O(n^3)\), but with a bigger constant — roughly three times the work of one LU factorization.

Given that, here is the practical advice, which students often get wrong:

The third point is the most important one in practice and is not intuitive, so here is a concrete instance of it.

The point to take away is that knowing whether \(A\) is invertible is important; actually computing \(A^{-1}\) rarely is.

Symmetry ought to buy us something in the factorization, and it does. Here the uniqueness theorem, Theorem 1.5, does some real work.

Theorem 1.10. If \(A\) is symmetric and nonsingular and admits the factorization \(A=LDU\), then \(U = L^T\), so that

The practical gain is a factor of two: only \(L\) and \(D\) need to be computed and stored. If in addition all the pivots are positive, we can split \(D = D^{1/2}D^{1/2}\) and write \(A = RR^T\) with \(R = LD^{1/2}\). This is the Cholesky factorization, and the matrices for which it works are exactly the positive definite matrices — the subject of a later chapter.

Poisson’s equation for the potential \(V\) due to a charge density \(\rho\) is \(\nabla^2 V = -\rho/\epsilon_0\). Take the one-dimensional version on the interval \(0 \le x \le 1\) with grounded conductors at both ends:

Two boundary conditions are exactly the right number, because if \(V(x)\) is a solution then so is \(V(x)+bx+c\) for any constants; the two conditions pin down \(b\) and \(c\). For a general \(\rho(x)\) there is no analytical solution, so we discretize.

Divide \([0,1]\) into \(n+1\) equal segments of length \(h = 1/(n+1)\), with interior grid points \(x_j = jh\) for \(j=1,\ldots,n\), and write \(V_j = V(x_j)\). The boundary conditions say \(V_0 = V_{n+1} = 0\). Approximating the first derivative by a difference of half-step values and repeating,

Substituting and multiplying through by \(-h^2\), the equation at grid point \(j\) is

with \(V_0\) and \(V_{n+1}\) set to zero. In matrix form, for \(n=5\),

A differential equation has become a matrix equation, and look at what the matrix has inherited. It is tridiagonal: only \(a_{i,i-1}, a_{ii}, a_{i,i+1}\) are nonzero, because the difference formula couples each point only to its immediate neighbours. It is symmetric, so \(A = LDL^T\) applies. And its pivots turn out to be positive, which makes it positive definite — the property that the later chapter on positive definite matrices is built around.

Now run elimination on the tridiagonal matrix and watch the cost. At step \(k\) only one row has a nonzero entry in column \(k\) (namely row \(k+1\)), and row \(k\) has only one nonzero to the right of the pivot. So a step costs a fixed number of operations rather than \((n-k)^2\), and no new nonzeros are created outside the three diagonals. The entire factorization costs \(O(n)\) instead of \(O(n^3)\).

This generalizes. Call \(w\) the half bandwidth of \(A\), meaning \(a_{ij}=0\) whenever \(|i-j| > w\). Tridiagonal is \(w=1\).

Lemma 1.11. If \(A\) has half bandwidth \(w\) and can be factored without row exchanges, then \(L\) and \(U\) have the same half bandwidth, and the factorization costs \(O(w^2 n)\).

The practical gain is large. For the discretized Poisson problem, \(w=1\) and we solve the system in \(O(n)\) work, both for the factorization and for each right hand side. Had we computed \(A^{-1}\) instead, it would have been dense (as shown earlier), needing \(O(n^2)\) storage and \(O(n^2)\) work for every solve.

Two practical remarks. First, partial pivoting spoils this a little: row exchanges can widen the band of \(U\) to \(2w\). The cost stays \(O(w^2n)\) in order of magnitude, so the loss is acceptable in exchange for stability. Second, the happy \(w=1\) situation is special to one dimension. Discretizing Poisson’s equation on a two-dimensional \(N\times N\) grid gives \(n=N^2\) unknowns with half bandwidth \(w=N\), so direct elimination costs \(O(w^2n) = O(N^4) = O(n^2)\). In three dimensions it is worse still. This is precisely the point at which one abandons direct elimination and turns to iterative methods such as the conjugate gradient method; those are covered in the optimization notes.