Applications of the LU Decomposition

Solving Systems of Linear Equations

Consider a system of \(\,m\,\) linear equations in \(\,n\,\) variables over a field \(\,K,\) \(\\\) with a coefficient matrix \(\,\boldsymbol{A}\in M_{m\times n}(K)\ \) and \(\,\) a column of constants \(\,\boldsymbol{b}\in K^m.\)

Assuming the decomposition \(\,\boldsymbol{P}\,\boldsymbol{A}=\boldsymbol{L}\,\boldsymbol{U},\ \) the system may be transformed as follows:

\[\begin{split}\begin{array}{c} \boldsymbol{A}\boldsymbol{x}\ =\ \boldsymbol{b}\,, \\ \\ (\boldsymbol{P}\,\boldsymbol{A})\,\boldsymbol{x}\ =\ \boldsymbol{P}\,\boldsymbol{b}\,, \\ \\ (\boldsymbol{L}\,\boldsymbol{U})\,\boldsymbol{x}\ =\ \boldsymbol{P}\,\boldsymbol{b}\,, \\ \\ \boldsymbol{L}(\boldsymbol{U}\boldsymbol{x})\ =\ \boldsymbol{P}\,\boldsymbol{b}\,. \end{array}\end{split}\]

Denoting \(\,\boldsymbol{y}\ =\ \boldsymbol{U}\boldsymbol{x}\,\) we get a system of \(\,2m\,\) equations in \(\,m+n\,\) unknowns, \(\\\) separated into two subsystems with triangular coefficient matrices:

(1)\[\begin{split}\left\{\ \ \begin{array}{ll} \boldsymbol{L}\,\boldsymbol{y}\ =\ \boldsymbol{P}\,\boldsymbol{b}\,, \qquad\qquad & \boldsymbol{L}\in M_m(K),\ \ \boldsymbol{y}\in K^m \\ \boldsymbol{U}\,\boldsymbol{x}\ =\ \boldsymbol{y}\,, \qquad\qquad & \boldsymbol{U}\in M_{m\times n}(K),\ \ \boldsymbol{x}\in K^n\,. \end{array}\right.\end{split}\]

The resulting subsystems are easily solved by forward or backward substitution. This approach is especially efficient in the case of several systems with the same matrix \(\,\boldsymbol{A},\ \) but with different columns \(\,\boldsymbol{b},\ \) since the decomposition \(\,\boldsymbol{P}\boldsymbol{A}=\boldsymbol{L}\,\boldsymbol{U}\ \) is then performed only once.

As an example, we shall solve the following system over the rational field \(\,Q:\)

\begin{alignat*}{4} 2\,x_1 & {\,} + {\,} & x_2 & {\,} + {\,} & x_3 & {\;} = {} & 6 \\ 4\,x_1 & {\,} - {\,} & 6\,x_2 & & & {\;} = {} & 14 \\ -2\,x_1 & {\,} + {\,} & 7\,x_2 & {\,} + {\,} & 2\,x_3 & {\;} = {} & -5 \end{alignat*}

The decomposition \(\,\boldsymbol{A}=\boldsymbol{L}\,\boldsymbol{U}\ \) for the coefficient matrix has been done in the previous section:

\[\begin{split}\left[\begin{array}{rrr} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2\end{array}\right]\ =\ \left[\begin{array}{rrr} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1\end{array}\right] \left[\begin{array}{rrr} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & 1\end{array}\right]\,.\end{split}\]

Eqs. (1) give two interconnected linear systems with triangular coefficient matrices:

\begin{alignat*}{4} y_1 & {\,} {\,} & & {\,} {\,} & & {\;} = {} & 6 \\ 2\,y_1 & {\,} + {\,} & y_2 & & & {\;} = {} & 14 \\ -y_1 & {\,} - {\,} & y_2 & {\,} + {\,} & y_3 & {\;} = {} & -5 \end{alignat*}

\begin{alignat*}{4} 2\,x_1 & {\,} + {\,} & x_2 & {\,} + {\,} & x_3 & {\;} = {\ } & y_1 \\ & {\,} - {\,} & 8\,x_2 & {\,} - {\,} & 2\,x_3 & {\;} = {\ } & y_2 \\ & {\,} {\,} & & {\,} {\,} & x_3 & {\;} = {\ } & y_3 \end{alignat*}

The forward and backward substitution in the first and second system, respectively, yields

\[\begin{split}\textstyle \begin{array}{l} y_1\ =\ 6 \\ y_2\ =\ 14\ -\ 2\,y_1\ =\ 2 \\ y_3\ =\ -5\ +\ y_1\ +\ y_2\ =\ 3\,, \end{array} \qquad \begin{array}{l} x_3\ =\ y_3\ =\ 3 \\ x_2\ =\ -{1\over 8}\ (y_2+2\,x_3)\ =\ -1 \\ x_1\ =\ {1\over 2}\ (y_1-x_2-x_3)\ =\ 2\,. \end{array}\end{split}\]

Finally, the solution reads: \(\ x_1=2,\ x_2=-1,\ x_3=3\,.\)

Inversion of a Matrix

Let \(\,\boldsymbol{A}\in M_n(K)\,\) be an invertible matrix. \(\\\) The problem of finding its inverse consists in solving the matrix equation

\[\boldsymbol{A}\,\boldsymbol{X}\ =\ \boldsymbol{I}_n\,,\]

where \(\,\boldsymbol{I}_n\ \) is the identity matrix of size \(\,n.\ \)

Using the column form of matrices

\[\boldsymbol{X}\ =\ \left[\, \boldsymbol{x}_1\,|\, \boldsymbol{x}_2\,|\, \ldots\,|\, \boldsymbol{x}_n\, \right]\,, \qquad \boldsymbol{I}_n\ =\ \left[\, \boldsymbol{e}_1\,|\, \boldsymbol{e}_2\,|\, \ldots\,|\, \boldsymbol{e}_n\, \right]\,,\]

and recalling the Column Rule of Matrix Multiplication, \(\\\) we obtain the set of \(\,n\,\) systems of linear equations, \(\,\) each in \(\,n\,\) unknowns:

\[ \begin{align}\begin{aligned}\boldsymbol{A}\ \left[\,\boldsymbol{x}_1\,|\,\boldsymbol{x}_2\,|\,\ldots\,|\, \boldsymbol{x}_n\,\right]\ =\ \left[\,\boldsymbol{e}_1\,|\,\boldsymbol{e}_2\,|\,\ldots\,|\, \boldsymbol{e}_n\,\right]\,,\\\left[\, \boldsymbol{A}\boldsymbol{x}_1\,|\, \boldsymbol{A}\boldsymbol{x}_2\,|\,\ldots\,|\, \boldsymbol{A}\boldsymbol{x}_n\,\right]\ =\ \left[\,\boldsymbol{e}_1\,|\,\boldsymbol{e}_2\,|\,\ldots\,|\, \boldsymbol{e}_n\,\right]\,,\\\boldsymbol{A}\,\boldsymbol{x}_j\ =\ \boldsymbol{e}_j\,, \qquad j=1,2,\ldots,n.\end{aligned}\end{align} \]

All systems have the same coefficient matrix \(\,\boldsymbol{A},\ \) but differ in columns of constants. As remarked in the previous subsection, in such a situation the (presented therein) method of solution based on the LU decomposition may be a convenient one.

Calculation of Determinants

Suppose we have to calculate the determinant of a matrix \(\,\boldsymbol{A}\in M_n(K)\,\) given in the form:

(2)\[\boldsymbol{A}\ =\ \boldsymbol{P}\,\boldsymbol{L}\,\boldsymbol{U}\,,\]

where \(\ \ \boldsymbol{P}=\boldsymbol{P}_\sigma\,,\ \ \sigma\in S_n\,,\quad \boldsymbol{L}=[l_{ij}]_{n\times n}\,,\quad \boldsymbol{U}=[u_{ij}]_{n\times n}\,.\)

According to the Cauchy’s Theorem on the determinant of a product of matrices, we get

\[\det\boldsymbol{A}\ \,=\ \, \det\boldsymbol{P}_\sigma\,\cdot\,\det\boldsymbol{L}\,\cdot\,\det\boldsymbol{U}\,.\]

Representing the permutation \(\,\sigma\,\) as a product of transpositions, it’s easy to show that

\[\det\boldsymbol{P}_\sigma\,=\,\text{sgn}\,\sigma\,.\]

The determinant of a triangular matrix is equal to the product of its diagonal elements:

\[\det\boldsymbol{L}\ =\ l_{11}\,l_{22}\,\dots\,l_{nn}\,,\qquad \det\boldsymbol{U}\ =\ u_{11}\,u_{22}\,\dots\,u_{nn}\,.\]

Hence the formula for the determinant of a LU-decomposed matrix \(\,\boldsymbol{A}\ \) given by (2) is

(3)\[\det\boldsymbol{A}\ =\ \text{sgn}\,\sigma\,\cdot\, l_{11}\,l_{22}\,\dots\,l_{nn}\,\cdot\, u_{11}\,u_{22}\,\dots\,u_{nn}\,.\]

As an example, let’s take the matrix \(\ \ \boldsymbol{A}\ =\ \left[\begin{array}{rrr} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{array}\right],\ \) \(\\\)\(\\\) for which two decompositions have been already derived:

\[\begin{split}\begin{array}{rl} \left[\begin{array}{rrr} 2 & 1 & 1 \\ 4 & -6 & 0 \\ -2 & 7 & 2 \end{array}\right] & =\ \ \ \left[\begin{array}{rrr} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right]\ \cdot\ \left[\begin{array}{rrr} 1 & 0 & 0 \\ \textstyle{1\over 2} & 1 & 0 \\ -\textstyle{1\over 2} & 1 & 1 \end{array}\right]\ \cdot\ \left[\begin{array}{rrr} 4 & -6 & 0 \\ 0 & 4 & 1 \\ 0 & 0 & 1 \end{array}\right]\,, \\ \\ & =\ \ \ \left[\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right]\ \cdot\ \left[\begin{array}{rrr} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & -1 & 1 \end{array}\right]\ \cdot\ \left[\begin{array}{rrr} 2 & 1 & 1 \\ 0 & -8 & -2 \\ 0 & 0 & 1 \end{array}\right]\,. \end{array}\end{split}\]

In both cases the formula (3) yields the same result: \(\ \det{\boldsymbol{A}}=-16.\)