Decomposition into Gaussian Triangular Factors

We know that Gaussian Elimination achieves the factorization of the matrix of coefficients, but we want to look at other ways of doing this. First we give a uniqueness result.

Theorem 2.2.1   Let L₁U₁ and L₂U₂ both be LU decompositions of nonsingular $A\in \hbox{{\sf I}\kern-.4em\hbox{\sf C}}^{n\times n}$, where L₁,L₂ each have unit diagonal entries. Then L₁=L₂ and U₁=U₂.

Proof.   A=L₁U₁=L₂U₂, so L₂^-1L₁=U₂U₁^-1. The left hand member is lower unit triangular and the right hand member is upper triangular, so both must equal I, which gives the result.

Gaussian elimination produces triangular factors via a number of intermediate stages. Now that we know the factors are unique, why not obtain them more directly, that is, compute each element once and for all (rather than as in GE, where elements are successively reduced).

Let L=[m_ij], U=[u_ij] and A=[a_ij]. Then LU=A becomes

$$\sum _{t=1}^{\min i,j}m_{it}u_{tj}=a_{ij}, \quad i,j=1,2,\ldots ,n.$$ (2.1)

We shall choose m_kk=1 for the moment to conform with Gaussian Elimination, but we shall use other choices later.

It is convenient to change i to k to find the elements of U:

For $i=k=1, ~ j=1,2,\ldots ,n,$ we have u_1j=a_1j.

For $i=k>1, j=i,i+1,\ldots ,n,$ (5.1) becomes $\sum _{t=1}^{k-1}m{_kt}u_{tj}+m_{kk}u_{kj}=a_{kj}$, and solving this gives

$$u_{kj}=a_{kj}-\sum _{t=1}^{k-1}m_{kt}u_{tj}.$$ (2.2)

For the elements of L, we change j to k:

For $i=2,3,\ldots ,n,~ j=k=1$ we have $m_{i1}=a_{i1}/u_{11}\mbox{ if }u_{11}\not= 0$.

For $j=k>1,~ i=j+1,j+2,\ldots ,n,$ (5.1) becomes $\sum _{t=1}^{k-1}m_{it}u_{tk}+m_{ik}u_{kk}=a_{ik}$, which yields

$$m_{ik}=\frac{1}{u_{kk}}(a_{ik}-\sum _{t=1}^{k-1}m_{it}u_{tk}) \mbox{ if }u_{kk}\not= 0.$$ (2.3)

These equations give the elements of U and L in terms of previously computed elements, provided we use them in an appropriate order. This is not unique, but one possibility is to use (5.2) to compute the kth row of U for $k=1,2,\ldots ,n,$ followed by (5.3) to compute the kth column of L. This procedure is illustrated below. The right hand side or sides can be reduced at the same time using f₁⁽ⁿ⁾=f₁ and

$$f_k^{(n)}=f_k-\sum _{k=1}^{k-1}m_{kt}f_t^{(n)},\quad k=2,3,\ldots ,n,$$

and the usual back-substitution is used for $\mbox{\boldmath$\space x $ }$.

Example.  Solve the equations represented by the augmented matrix

$$\left[ \begin{array}{rrrr\vert r}2 & 1 & 1 & 1 & 5\\ 1 & -2 &... ...\\ 1 & 1 & 1 & 1 & 4\\ 1 & -1 & 2 & -1 & 1\end{array} \right].$$

The following matrix shows the order in which the solution is computed:

$$\left[ \begin{array}{cccc\vert c}(1) & (1) & (1) & (1) & (1)\... ... & (5) & (5)\\ (2) & (4) & (6) & (7) & (7)\end{array} \right].$$

The lower and upper factors are put together into one matrix, in which the diagonal elements of U are shown, since those of L are all 1:

$$\left[ \begin{array}{rrrr\vert r}2.0 & 1.0 & 1.0 & 1.0 & 5.0\... ...& -0.4 & 0.6\\ 0.5 & 0.6 & 0.0 & 1.2 & 1.2\end{array} \right].$$

Interchanges can, and indeed should, be incorporated, but it is somewhat complicated, so we omit the details. Note that each entry is computed in a single operation, which can conveniently be done in double precision arithmetic.