4.4 Matrix factorization
Often several systems with the same matrix \(A\) and different right-hand sides must be solved. Gaussian elimination repeats the expensive \(O(n^3)\) work on \(A\) for every \(\mathbf b\). Recording the elimination once, as a factorization \(A = LU\), reduces each further solution to \(O(n^2)\) operations.
Theorem 4.4.1 (LU factorization). If Gaussian elimination can be performed on \(A\mathbf x = \mathbf b\) without row interchanges, then \(A = LU\), where \(U\) is the upper triangular matrix produced by the elimination and \(L\) is the unit lower triangular matrix (ones on the diagonal) whose entries below the diagonal are the multipliers: \[L = \begin {pmatrix}1&0&\cdots& 0\\ m_{21}&1&\cdots& 0\\ \vdots&& \ddots& \vdots \\ m_{n1}&m_{n2}&\cdots& 1\end {pmatrix}.\]
Proof. Let \(\mathbf m^{(k)} = (0,\dots ,0,m_{k+1,k},\dots ,m_{nk})^T\), with zeros in its first \(k\) entries, and \(\mathbf e_k\) the \(k\)th unit vector. Stage \(k\) of elimination subtracts \(m_{jk}\) times row \(k\) from row \(j\) for \(j>k\), which is left multiplication by \[M^{(k)} = I-\mathbf m^{(k)}\mathbf e_k^T .\] So \(U = M^{(n-1)}\cdots M^{(2)}M^{(1)}A\). Since \(\mathbf e_k^T\mathbf m^{(k)} = 0\), \[(I-\mathbf m^{(k)}\mathbf e_k^T)(I+\mathbf m^{(k)}\mathbf e_k^T) = I-\mathbf m^{(k)}(\mathbf e_k^T\mathbf m^{(k)})\mathbf e_k^T = I,\] so \(\left (M^{(k)}\right )^{-1} = I+\mathbf m^{(k)}\mathbf e_k^T\), and \[A = \left (M^{(1)}\right )^{-1}\left (M^{(2)}\right )^{-1}\cdots \left (M^{(n-1)}\right )^{-1}U .\] Multiplying out the product, every cross term contains a factor \(\mathbf e_j^T\mathbf m^{(k)}\) with \(j<k\), which is \(0\) because the first \(k\) entries of \(\mathbf m^{(k)}\) are zero. Hence the product equals \(I+\sum _k\mathbf m^{(k)}\mathbf e_k^T\), which is exactly \(L\). □
Once \(A = LU\) is known, \(A\mathbf x = \mathbf b\) is solved in two triangular steps: first \(L\mathbf y = \mathbf b\) by forward substitution, then \(U\mathbf x = \mathbf y\) by backward substitution. Each costs about \(n^2\) operations.
Example 4.4.2. Factor the matrix of Example 4.2.1 and use the factors to solve \(A\mathbf x = (8,7,14,-7)^T\).
Solution. The elimination in Example 4.2.1 used the multipliers \(m_{21} = 2\), \(m_{31} = 3\), \(m_{41} = -1\), \(m_{32} = 4\), \(m_{42} = -3\), \(m_{43} = 0\), so \begin {align*} A &= \begin {pmatrix}1&1&0&3\\2&1&-1&1\\3&-1&-1&2\\-1&2&3&-1\end {pmatrix}\\ &= \begin {pmatrix}1&0&0&0\\2&1&0&0\\3&4&1&0\\-1&-3&0&1\end {pmatrix}\begin {pmatrix}1&1&0&3\\0&-1&-1&-5\\0&0&3&13\\0&0&0&-13\end {pmatrix} = LU . \end {align*}
Forward substitution in \(L\mathbf y = \mathbf b\): \[y_1 = 8,\quad y_2 = 7-2y_1 = -9,\quad y_3 = 14-3y_1-4y_2 = 26,\quad y_4 = -7+y_1+3y_2 = -26 .\] Backward substitution in \(U\mathbf x = \mathbf y\): \(x_4 = 2\), \(x_3 = \frac {26-13x_4}{3} = 0\), \(x_2 = \frac {-9+x_3+5x_4}{-1} = -1\), \(x_1 = 8-x_2-3x_4 = 3\). So \(\mathbf x = (3,-1,0,2)^T\).
Example 4.4.3. Compare the cost of solving \(k\) systems \(A\mathbf x = \mathbf b_j\) (\(j = 1,\dots ,k\)) by repeated Gaussian elimination and by one LU factorization, for \(n = 100\) and \(k = 100\).
Solution. Repeated elimination costs about \(k\cdot \frac {n^3}{3}\) multiplications, that is \(3.3\times 10^7\). One factorization costs about \(\frac {n^3}{3}\) and each pair of triangular solves about \(n^2\), a total of about \(\frac {n^3}{3}+kn^2 = 3.3\times 10^5+10^6 = 1.3\times 10^6\): some \(25\) times less. This is why the inverse matrix is almost never computed: \(A^{-1}\mathbf b\) costs more than the triangular solves and is less accurate.
Row interchanges
If row interchanges are needed, they can be performed in advance. A permutation matrix \(P\) is the identity matrix with its rows reordered; \(PA\) reorders the rows of \(A\) in the same way.
Theorem 4.4.4. If \(A\) is nonsingular, there is a permutation matrix \(P\) such that \(PA = LU\), with \(L\) unit lower triangular and \(U\) upper triangular. To solve \(A\mathbf x = \mathbf b\), solve \(L\mathbf y = P\mathbf b\) and then \(U\mathbf x = \mathbf y\).
The permutation simply records the interchanges made during Gaussian elimination (with or without pivoting); applying them to \(A\) first, elimination needs none.
Example 4.4.5. Find a factorization \(PA = LU\) for the matrix of Example 4.2.2.
Solution. Elimination needed one interchange, of rows \(2\) and \(3\). So take \[P = \begin {pmatrix}1&0&0&0\\0&0&1&0\\0&1&0&0\\0&0&0&1\end {pmatrix},\qquad PA = \begin {pmatrix}1&-1&2&-1\\1&1&1&0\\2&-2&3&-3\\1&-1&4&3\end {pmatrix}.\] Elimination on \(PA\) needs no interchange and gives \[L = \begin {pmatrix}1&0&0&0\\1&1&0&0\\2&0&1&0\\1&0&-2&1\end {pmatrix},\qquad U = \begin {pmatrix}1&-1&2&-1\\0&2&-1&1\\0&0&-1&-1\\0&0&0&2\end {pmatrix}.\]
Direct computation of the factors
The factors can be computed directly from the equations \(a_{ij} = \sum _kl_{ik}u_{kj}\), row by row of \(U\) and column by column of \(L\). With \(l_{ii} = 1\) (Doolittle’s method): \[u_{ij} = a_{ij}-\sum _{k=1}^{i-1}l_{ik}u_{kj}\quad (j\geq i),\qquad l_{ji} = \frac 1{u_{ii}}\left (a_{ji}-\sum _{k=1}^{i-1}l_{jk}u_{ki}\right )\quad (j>i).\] Crout’s method instead takes \(u_{ii} = 1\) and puts the pivots on the diagonal of \(L\). If \(A = LU\) (Doolittle) and \(D\) is the diagonal of \(U\), then \(A = (LD)(D^{-1}U)\) is the Crout factorization.
Example 4.4.6. Find the Doolittle and Crout factorizations of \(A = \begin {pmatrix}2&-1&1\\3&3&9\\3&3&5\end {pmatrix}\).
Solution. Doolittle. Row 1 of \(U\) is row 1 of \(A\): \(u_{11} = 2\), \(u_{12} = -1\), \(u_{13} = 1\). Column 1 of \(L\): \(l_{21} = \frac 32\), \(l_{31} = \frac 32\). Row 2 of \(U\): \(u_{22} = 3-\frac 32(-1) = \frac 92\), \(u_{23} = 9-\frac 32(1) = \frac {15}{2}\). Column 2 of \(L\): \(l_{32} = \frac {3-\frac 32(-1)}{9/2} = 1\). Finally \(u_{33} = 5-\frac 32(1)-1\cdot \frac {15}{2} = -4\). So \[L = \begin {pmatrix}1&0&0\\\frac 32&1&0\\\frac 32&1&1\end {pmatrix},\qquad U = \begin {pmatrix}2&-1&1\\0&\frac 92&\frac {15}2\\0&0&-4\end {pmatrix}.\] Crout. With \(D = \operatorname {diag}(2,\frac 92,-4)\), \[LD = \begin {pmatrix}2&0&0\\3&\frac 92&0\\3&\frac 92&-4\end {pmatrix},\qquad D^{-1}U = \begin {pmatrix}1&-\frac 12&\frac 12\\0&1&\frac 53\\0&0&1\end {pmatrix}.\] In both cases \(\det A = 2\cdot \frac 92\cdot (-4) = -36\), the product of the pivots.
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.