LU decomposition







Numerical Methods

LU decomposition

Matrix $A$ can be written as the product $A=LU$ where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix of the form,

\[\begin{equation} \left[ \begin{array}{cccc} A_{11} & A_{12} & A_{13} & A_{14} \\ A_{21} & A_{22} & A_{23} & A_{24} \\ A_{31} & A_{32} & A_{33} & A_{34} \\ A_{41} & A_{42} & A_{43} & A_{44} \end{array} \right]= \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ L_{21} & 1 & 0 & 0 \\ L_{31} & L_{32} & 1 & 0 \\ L_{41} & L_{42} & L_{43} & 1 \end{array} \right]\left[ \begin{array}{cccc} U_{11} & U_{12} & U_{13} & U_{14} \\ 0 & U_{22} & U_{23} & U_{24} \\ 0 & 0 & U_{33} & U_{34} \\ 0 & 0 & 0 & U_{44} \end{array} \right]. \end{equation}\]

Matrix multiplication yields,

\[\begin{equation} \begin{array}{y} A_{11} = U_{11} \\ A_{21} = L_{21}U_{11} \\ A_{31} = L_{31}U_{11}\\ A_{41} = L_{41}U_{11}\\ \end{array} \qquad \begin{array}{y} A_{12} = U_{12} \\ A_{22} = L_{21}U_{12}+U_{22} \\ A_{32} = L_{31}U_{12}+L_{32}U_{22} \\ A_{42} = L_{41}U_{12}+L_{42}U_{22} \\ \end{array} \qquad \begin{array}{y} A_{13} = U_{13} \\ A_{23} = L_{21}U_{13}+U_{23} \\ A_{33} = L_{31}U_{13}+L_{32}U_{23} +U_{33}\\ A_{43} = L_{41}U_{13}+L_{42}U_{23} +L_{43}U_{33}\\ \end{array} \qquad \begin{array}{y} A_{14} = U_{14} \\ A_{24} = L_{21}U_{14}+U_{24} \\ A_{34} = L_{31}U_{14}+L_{32}U_{24} + U_{34}\\ A_{44} = L_{41}U_{14}+L_{42}U_{24} + L_{43}U_{34}+U_{44}\\ \end{array} \end{equation}\]

These equations can be solved starting in the upper-left and solving $U_{1j}$ in the first row then for $L_{i1}$ in the first column, then for $U_{2j}$ in the second row and for $L_{i2}$ in the second column, etc. Code that performs LU decomposition is shown below.

<script> var A = [[2,7,6,2],[9,5,1,3],[4,3,8,4],[5,6,7,8]]; var N = A.length; document.write("A = <br \/>"); for (i=0; i<N; i++) { for (j=0; j<N; j++) { document.write(A[i][j]+" "); } document.write("<br \/>"); } document.write("<br \/>"); var L = []; var U = []; function identity(N) { // returns an NxN identity matrix var I = []; for (j=0; j<N; j++) { I[j] = new Array(N); for (k=0; k<N; k++) { I[j][k] = 0; } I[j][j] = 1; } return I; } L = identity(N); U = identity(N); for (i=0; i<N; i++) { j = i; while (j<N) { //Determine U across row i q = 0; U[i][j] = A[i][j]; while (q<i) { U[i][j] -= L[i][q]U[q][j]; q++; } j++; } j = i+1; while (j<N) { //Determine L down column i q = 0; L[j][i] = A[j][i]; while (q<i) { L[j][i] -= L[j][q]U[q][i]; q++; } L[j][i] = L[j][i]/U[i][i]; j++; } } document.write("L = <br \/>"); for (i=0; i<N; i++) { for (j=0; j<N; j++) { document.write(L[i][j]+" "); } document.write("<br \/>"); } document.write("<br>U = <br>"); for (i=0; i<N; i++) { for (j=0; j<N; j++) { document.write(U[i][j]+" "); } document.write("<br \/>"); } </script>

Since the matrix elements $A_{ij}$ are only used once in the calculation, it possible to store the components of $L_{ij}$ and $U_{ij}$ in the same memory locations where $A_{ij}$ was stored. Code that generates the $LU$ matrix is shown below. It is more compact than the version above and uses less memory but is harder to read.

<script> function LU_matrix(A) { //returns the LU matrix of matrix A N = A.length; for (i=0; i<N; i++) { j = i; while (j<N) { //Determine U across row i q = 0; while (q<i) { A[i][j] -= A[i][q]A[q][j]; q++; } j++; } j = i+1; while (j<N) { //Determine L down column i q = 0; while (q<i) { A[j][i] -= A[j][q]A[q][i]; q++; } A[j][i] = A[j][i]/A[i][i]; j++; } } return A; } function write_matrix(A) { rows = A.length; cols = A[0].length; for (i=0; i<rows; i++) { for (j=0; j<cols; j++) { document.write(A[i][j]+" "); } document.write("<\/br> "); } document.write("<\/br> "); } var A = [[2,7,6,2],[9,5,1,3],[4,3,8,4],[5,6,7,8]]; document.write("A = <br \/>"); write_matrix(A); LU = LU_matrix(A); document.write("LU = <br \/>"); write_matrix(LU); </script>

LU-decomposition

Optimization of the rounding error through partial pivoting

The numerical results obtained are not affected by methodological error (this is a direct method). All deviations from the 'true' results, therefore, must come from rounding errors. We give now a concrete example. Let us consider the following set of three equations:

\begin{eqnarray}\label{lin14} x_{1} + 5923181 x_{2} + 1608 x_{3} = 5924790 \nonumber\\ 5923181 x_{1} + 337116 x_{2} - 7 x_{3} = 6260290 \\ 6114 x_{1} + 2 x_{2} + 9101372 x_{3} = 9107488 \nonumber \end{eqnarray}

The exact solution of this system is: $x_{1}=x_{2}=x_{3}=1$. This can easily be checked by substituting this solution into the equations.

The numerical solution yields the following LU-Matrix and vector:

\begin{equation} LU =\begin{bmatrix}1&5923181&1608\\5923181&-35084072821645&-9524475055\\6114&0.001032215638591928&9101372.101149714\end{bmatrix}, \nonumber \end{equation} \begin{equation} \vec{x} =\begin{bmatrix}0.9999999998590283\\1\\1.0000000000000877\end{bmatrix}. \nonumber \end{equation}

We can see the effect of the rounding error. If we now change the order of the equations to 2/3/1, which in principle should have no influence on the result, the same program returns:

\begin{equation} LU =\begin{bmatrix}6114&2&9101372\\968.7898266274125&335178.4203467452&-8817316608.951586\\0.00016355904481517829&17.671725385975915&155817197873.98508\end{bmatrix}, \nonumber \end{equation} \begin{equation} \vec{x} =\begin{bmatrix}0.9999999999993907\\1.0000000000070886\\1.0000000000000004\end{bmatrix}. \nonumber \end{equation}

The result is now different because of the rounding errors. If we re-order the lines with the sequence 2/1/3 we obtain the correct solution up to machine precision:

\begin{equation} LU =\begin{bmatrix}5923181&337116&-7\\1.6882820227847166e-7&5923180.943085312&1608.0000011817974\\0.0010322156287305756&-0.000058410574861642146&9101372.101149714\end{bmatrix}, \nonumber \end{equation} \begin{equation} \vec{x} =\begin{bmatrix}1\\1\\1\end{bmatrix}. \nonumber \end{equation}

This means that the order in which the equations appear affect the size of the rounding error and on the quality of the numerical solution in a decisive way.

What is the 'right' order? Round off errors occur when large, nearly-equal numbers are subtracted from each other. In the construction of the LU-matrix, we divide by the diagonal elements $U_{ii}$ so large elements will appear if $U_{ii}$ is small. Even worse than the round-off error is the fact that if the upper-left element of the matrix $A$ is zero, our LU-decomposition routine divides by zero and fails. To try to prevent this from happening we rearrange the rows of the equations so that the matrix has as large as possible elements along the diagonal. This is known as pivoting. It is possible to rearrange the rows and the columns. If two rows are swapped, this just changes the order in which the equations appear but has no effect on the solution. If two columns are swapped, it relabels the variables. For instance, if column 1 and column 3 are swapped, the variable that was named $x_1$ is now nammed $x_3$ and the old $x_3$ is now named $x_1$. This also has no effect on the solutions. Often just the rows are exchanged during pivoting. This is known as partial pivoting. The function pivot(A) defined in the code below puts the largest elements of matrix $A$ along the diagonal. It looks down the first column for the element with the largest absolute value and swaps that row to the first row. Then it looks down the second column starting from the second row for the element with the largest absolute value and swaps that with the second row etc.

<script> function pivot(A) { //determines pivot matrix for A N = A.length; P = identity(N); for (j = 0; j<N; j++) { row = j; Amax = 0; for (i = j; i <N; i++) {//look down column j for largest element absA = Math.abs(A[i][j]); if (absA>Amax) { Amax = absA; row = i; } } if (j!=row) {//swap rows for (q=0; q<N; q++) { tmp = P[j][q]; P[j][q]=P[row][q]; P[row][q]=tmp; tmp = A[j][q]; A[j][q]=A[row][q]; A[row][q]=tmp; } } } return {P:P,PA:A};//returns permutation matrix P and PA } function identity(N) { // returns an NxN identity matrix var I = []; for (j=0; j<N; j++) { I[j] = new Array(N); for (k=0; k<N; k++) { I[j][k] = 0; } I[j][j]=1; } return I; } function write_matrix(A) { rows = A.length; cols = A[0].length; for (i=0; i<rows; i++) { for (j=0; j<cols; j++) { document.write(A[i][j]+" "); } document.write("<\/br> "); } document.write("<\/br> "); } var A = [[2,7,6,2],[9,5,1,3],[4,3,8,4],[5,6,7,8]]; write_matrix(A); document.write("The permutation matrix P that rearranges the rows of A is:<br \/>"); Pvt = pivot(A); document.write("P = <br />"); write_matrix(Pvt.P); document.write("The permutation matrix P times A is:<br \/>"); document.write("PA = <br \/>"); write_matrix(Pvt.PA); </script>

The function pivot(A) accepts a square matrix as its argument and returns an object that contains two square matrices called $P$ and $PA$. $P$ is the permutation matrix that rearranges the matrix $A$. A permutation matrix has the property that one element equals $1$ in every row and every column. All other elements are zero. The other matrix that pivot(A) returns is the matrix product of $P$ an $A$. This matrix is passed to the LU-decomposition routine.

In general, expecially for higher-order systems, there is no chance to get a correct result without using some kind of pivoting. There are however also sets of equations for which the LU-decomposition is stable also without pivoting. This happens if the matrix is symmetric, tridiagonal, cyclically tridiagonal, diagonal-dominant or positive definite.