The process known as Gaussian elimination is equivalent to computing the LU factorization of the matrix \(A \in \Cmxm \)
\begin{equation*}
:
A = L U,
\end{equation*}
where \(L \) is a unit lower trianguular matrix and \(U\) is an upper triangular matrix.
Definition5.6.2.1.
Given a matrix \(A \in \C^{m \times n} \) with \(m \geq n \text{,}\) its LU factorization is given by \(A = L U \) where \(L \in \C^{m\times n}\) is unit lower trapezoidal and \(U \in \C^{n \times n} \) is upper triangular with nonzeroes on its diagonal.
Definition5.6.2.2. Principal leading submatrix.
For \(k \leq n \text{,}\) the \(k \times k \) principal leading submatrix of a matrix \(A \) is defined to be the square matrix \(A_{TL} \in \C^{k \times k} \) such that \(A =
\FlaTwoByTwo{A_{TL}}{A_{TR}}{A_{BL}}{A_{BR}} \text{.}\)
Lemma5.6.2.3.
Let \(L \in \Cnxn \) be a unit lower triangular matrix and \(U \in \Cnxn\) be an upper triangular matrix. Then \(A = L U \) is nonsingular if and only if \(U \) has no zeroes on its diagonal.
Theorem5.6.2.4. Existence of the LU factorization.
Let \(A \in \C^{m \times n} \) and \(m \geq n \) have linearly independent columns. Then \(A \) has a (unique) LU factorization if and only if all its principal leading submatrices are nonsingular.
Figure5.6.2.6. Left-looking LU factorization algorithm. \(L_{00} \) is the unit lower triangular matrix stored in the strictly lower triangular part of \(A_{00} \) (with the diagonal implicitly stored).
Solving \(A x = b \) via LU factorization:
Compute the LU factorization \(A = L U \text{.}\)
Solve \(L z = b \text{.}\)
Solve \(U x = z \text{.}\)
Cost of LU factorization: Starting with an \(m \times n \) matrix \(A \text{,}\) LU factorization requires approximately \(m n^2 - \frac{1}{3} n^3 \) flops. If \(m = n \) this becomes
where \(\{ \pi_0, \pi_1, \ldots ,
\pi_{n-1} \} \) is a permutation (rearrangement) of the integers \(\{ 0 , 1, \ldots , n-1 \} \text{,}\) we define the permutation matrix \(P( p )\) by
\begin{equation*}
\newcommand{\routinename}{\mbox{Solve } U x = z,
\mbox{ overwriting } z \mbox{ with } x \mbox{ (Variant 1)}}
\newcommand{\guard}{
n( U_{BR} ) \lt n( U )
}
\newcommand{\partitionings}{
U \rightarrow
\FlaTwoByTwo{U_{TL}}{U_{TR}}
{U_{BL}}{U_{BR}},
z \rightarrow
\FlaTwoByOne{z_T}{z_B}
}
\newcommand{\partitionsizes}{
U_{BR} {\rm ~is~} 0 \times 0 \mbox{ and }
z_B \mbox{ has } 0 \mbox{ elements}
}
\newcommand{\repartitionings}{
\FlaTwoByTwo{U_{TL}}{U_{TR}}
{U_{BL}}{U_{BR}}
\rightarrow
\FlaThreeByThreeTL{U_{00}}{u_{01}}{U_{02}}
{u_{10}^T}{\upsilon_{11}}{u_{12}^T}
{U_{20}}{u_{21}}{U_{22}},
\FlaTwoByOne{z_T}{z_B}
\rightarrow
\FlaThreeByOneT{z_0}{\zeta_1}{z_2}
}
\newcommand{\moveboundaries}{
\FlaTwoByTwo{U_{TL}}{U_{TR}}
{U_{BL}}{U_{BR}}
\leftarrow
\FlaThreeByThreeBR{U_{00}}{u_{01}}{U_{02}}
{u_{10}^T}{\upsilon_{11}}{u_{12}^T}
{U_{20}}{u_{21}}{U_{22}},
\FlaTwoByOne{z_T}{z_B}
\leftarrow
\FlaThreeByOneB{z_0}{\zeta_1}{z_2}
}
\newcommand{\update}{
\begin{array}{l}
\zeta_1 := \zeta_1 / \upsilon_{11} \\
z_0 := z_0 - \zeta_1 u_{01}
\end{array}
}
\FlaAlgorithm
\end{equation*}
Figure5.6.2.13. Upper triangular solve Variant 1
\begin{equation*}
\newcommand{\routinename}{\mbox{Solve } U x = z,
\mbox{ overwriting } z \mbox{ with } x \mbox{ (Variant 2)}}
\newcommand{\guard}{
n( U_{BR} ) \lt n( U )
}
\newcommand{\partitionings}{
U \rightarrow
\FlaTwoByTwo{U_{TL}}{U_{TR}}
{U_{BL}}{U_{BR}},
z \rightarrow
\FlaTwoByOne{z_T}{z_B}
}
\newcommand{\partitionsizes}{
U_{BR} {\rm ~is~} 0 \times 0 \mbox{ and }
z_B \mbox{ has } 0 \mbox{ elements}
}
\newcommand{\repartitionings}{
\FlaTwoByTwo{U_{TL}}{U_{TR}}
{U_{BL}}{U_{BR}}
\rightarrow
\FlaThreeByThreeTL{U_{00}}{u_{01}}{U_{02}}
{u_{10}^T}{\upsilon_{11}}{u_{12}^T}
{U_{20}}{u_{21}}{U_{22}},
\FlaTwoByOne{z_T}{z_B}
\rightarrow
\FlaThreeByOneT{z_0}{\zeta_1}{z_2}
}
\newcommand{\moveboundaries}{
\FlaTwoByTwo{U_{TL}}{U_{TR}}
{U_{BL}}{U_{BR}}
\leftarrow
\FlaThreeByThreeBR{U_{00}}{u_{01}}{U_{02}}
{u_{10}^T}{\upsilon_{11}}{u_{12}^T}
{U_{20}}{u_{21}}{U_{22}},
\FlaTwoByOne{z_T}{z_B}
\leftarrow
\FlaThreeByOneB{z_0}{\zeta_1}{z_2}
}
\newcommand{\update}{
\begin{array}{l}
\zeta_1 := \zeta_1 - u_{12}^T z_2 \\
\zeta_1 := \zeta_1 / \upsilon_{11}
\end{array}
}
\FlaAlgorithm
\end{equation*}
Figure5.6.2.14. Upper triangular solve Variant 2
Cost of triangular solve Starting with an \(n \times n \) (upper or lower) triangular matrix \(T \text{,}\) solving \(T x = b \) requires approximately \(n^2 \) flops.
Provided the solution of \(A x = b \) yields some accuracy in the solution, that accuracy can be improved through a process known as iterative refinement.
Let \(\hat x \) is an approximate solution to \(A x = b
\text{.}\)
Let \(\widehat \delta\! x \) is an approximate solution to \(A \delta\! x = b - A \widehat x \text{,}\)
Then \(\widehat x + \widehat \delta\! x \text{,}\) is an improved approximation.
This process can be repeated until the accuracy in the computed solution is as good as warranted by the conditioning of \(A\) and the accuracy in \(b \text{.}\)
A matrix \(A \in \C^{n \times n} \) is Hermitian positive definite (HPD) if and only if it is Hermitian (\(A^H = A\)) and for all nonzero vectors \(x \in \C^n \) it is the case that \(x ^H A x \gt 0 \text{.}\) If in addition \(A \in \R^{n \times n} \) then \(A \) is said to be symmetric positive definite (SPD).
Some insights regarding HPD matrices:
\(B \) has linearly independent columns if and only if \(A = B^H B \) is HPD.
A diagonal matrix has only positive values on its diagonal if and only if it is HPD.
If \(A \) is HPD, then its diagonal elements are all real-valued and positive.
If \(A =
\left( \begin{array}{c | c }
A_{TL} \amp A_{TR} \\ \hline
A_{BL} \amp A_{BR}
\end{array} \right), \) where \(A_{TL} \) is square, is HPD, then \(A_{TL} \) and \(A_{BR} \) are HPD.
Theorem5.6.2.16. Cholesky Factorization Theorem.
Given an HPD matrix \(A \) there exists a lower triangular matrix \(L \) such that \(A = L L^H \text{.}\) If the diagonal elements of \(L \) are restricted to be positive, \(L \) is unique.
Figure5.6.2.17. Cholesky factorization algorithm (right-looking variant). The operation "syr" refers to "symmetric rank-1 update", which performs a rank-1 update, updating only the lower triangular part of the matrix in this algorithm.
Lemma5.6.2.18.
Let \(A
=
\left( \begin{array}{c | c }
\alpha_{11} \amp a_{21}^H \\ \hline
a_{21} \amp A_{22}
\end{array} \right),
\in \C^{n \times n} \) be HPD and \(l_{21} = a_{21} / \sqrt{\alpha_{11}} \text{.}\) Then \(A_{22} - l_{21} l_{21}^H \) is HPD.
Let \(\hat x \in \C^n \) equal the solution to the linear least-squares (LLS) problem
\begin{equation}
\| b - A \widehat x \|_2 = \min_{x \in \C^n} \| b - A x \|_2 ,\label{chapter05-cholesky-lls-eqn-1-summary}\tag{5.6.1}
\end{equation}
where \(A \) has linearly independent columns, equals the solution to the normal equations
\begin{equation*}
\begin{array}[t]{c}
\underbrace{ A^H A }\\
B
\end{array}
\widehat x =
\begin{array}[t]{c}
\underbrace{A^H b} \\
y
\end{array} .
\end{equation*}
This solution can be computed via the steps
Form \(B = A^H A \text{.}\) Cost: approximately \(m n^2 \) flops.
Factor \(B = L L^H \) (Cholesky factorization). Cost: approximately \(n^3/3 \) flops.
Compute \(y = A^H b \text{.}\) Cost: approximately \(2 m n\) flops.
Solve \(L z = y \text{.}\) Cost: approximately \(n^2 \) flops.
Solve \(L^H \hat x = z \text{.}\) Cost: approximately \(n^2 \) flops.
for a total of, approximately, \(m n^2 + n^3/3 \) flops.
