Direct Matrix Factorisations: LU, Cholesky and QR

Forward substitution. For $i=0$, [L1] gives $x_0=b_0\mathrm{/}{\mathrm{\ell }}_{00}$, so the first row equation is satisfied. Assume the first $i$ rows are satisfied. Then the $(i+1)$-st row of a lower-triangular system is ${\sum }_{j\le i+1}{\mathrm{\ell }}_{i+1,j}{x}_j=b_{i+1}$, and solving for $x_{i+1}$ gives exactly the recursion in [L1]. Because ${\mathrm{\ell }}_{i+1,i+1}\mathrm{\ne }0$, this value is forced and is unique. Induction proves that the forward recursion solves $Lx=b$ and that no other vector can differ in any coordinate.

Backward substitution is the same argument from the last row upward: the last row determines $x_{n-1}=b_{n-1}\mathrm{/}{u}_{n-1,n-1}$, and once $x_{i+1},\dots ,x_{n-1}$ are fixed, the $i$-th row forces $x_i=u_{ii}^{-1}(b_i-{\sum }_{j>i}{u}_{ij}{x}_j)$. Hence the backward recursion solves $Ux=b$ and the solution is unique.

In either triangular case, the $i$-th step forms one sum with $i$ products and $i$ additions in the lower-triangular solve, or with $n-1-i$ products and $n-1-i$ additions in the upper-triangular solve, followed by one division. Summing over all rows gives $n$ divisions and ${\sum }_{i=0}^{n-1}i={\sum }_{i=0}^{n-1}(n-1-i)=n(n-1)\mathrm{/}2$ multiplications and the same number of additions.

Step 1.1 proves claim 1, step 1.2 proves claim 2, and step 2.1 gives the operation counts.

For claim 1, the first row of $A$ equals the first row of both $U$ and $U^{\mathrm{\prime }}$, because the first rows of $L$ and $L^{\mathrm{\prime }}$ are $(1,0,\dots ,0)$. Hence the first rows of $U$ and $U^{\mathrm{\prime }}$ agree, in particular $u_{00}=u_{00}^{\mathrm{\prime }}\mathrm{\ne }0$. Then the first-column equations are $a_{i0}={\mathrm{\ell }}_{i0}{u}_{00}={\mathrm{\ell }}_{i0}^{\mathrm{\prime }}{u}_{00}^{\mathrm{\prime }}$ for $i>0$, so ${\mathrm{\ell }}_{i0}={\mathrm{\ell }}_{i0}^{\mathrm{\prime }}$. Subtracting the common first-column outer product leaves two factorizations of the same trailing principal block by unit lower and upper triangular matrices of size $n-1$ with nonzero diagonal. Induction gives equality of the trailing blocks, hence $L=L^{\mathrm{\prime }}$ and $U=U^{\mathrm{\prime }}$.

For claim 2, put $V:=DU$ and $V^{\mathrm{\prime }}:=D^{\mathrm{\prime }}{U}^{\mathrm{\prime }}$. Because $D,D^{\mathrm{\prime }}$ are diagonal and $U,U^{\mathrm{\prime }}$ are unit upper triangular, $V,V^{\mathrm{\prime }}$ are upper triangular and their diagonals are exactly the diagonals of $D,D^{\mathrm{\prime }}$, so they are nonzero. The factorisations $A=LV=L^{\mathrm{\prime }}{V}^{\mathrm{\prime }}$ are therefore of the kind handled in step 1.1, which gives $L=L^{\mathrm{\prime }}$ and $V=V^{\mathrm{\prime }}$. Taking diagonals of $V=V^{\mathrm{\prime }}$ yields $D=D^{\mathrm{\prime }}$, and then $DU=D^{\mathrm{\prime }}{U}^{\mathrm{\prime }}$ forces $U=U^{\mathrm{\prime }}$.

Steps 1.1 and 2.1 prove claims 1 and 2. The nonzero-pivot hypothesis is essential: without it, a zero row of the upper factor leaves lower-triangular multipliers undetermined.

By [L3], $A_{11}^{-1}$ exists, and by [L4] the block product $\left(\begin{array}{cc}I& 0\\ A_{21}{A}_{11}^{-1}& I\end{array}\right)\left(\begin{array}{cc}{A}_{11}& A_{12}\\ 0& S\end{array}\right)$ has upper-left block $A_{11}$, upper-right block $A_{12}$, lower-left block $A_{21}{A}_{11}^{-1}{A}_{11}=A_{21}$, and lower-right block $A_{21}{A}_{11}^{-1}{A}_{12}+S=A_{22}$. Hence it equals $A$.

If $S=\tilde{L}\tilde{U}$, substitute this into step 1.1 and absorb the factorisation of $S$ into the lower-right block. The resulting left factor is block lower triangular and the right factor is block upper triangular, so this is a block LU factorisation in the sense of [L1].

The off-diagonal block $A_{21}{A}_{11}^{-1}$ is exactly the block of multipliers obtained by solving $X{A}_{11}=A_{21}$. If $A_{11}$ itself has a triangular factorisation, [L2] computes these multipliers by triangular solves.

Steps 1.1 and 2.1 prove claims 1 and 2, and step 2.2 identifies the solve interpretation of the multiplier block.

Suppose $A=LU$ as in [L1] and that the diagonal of $U$ is nonzero. For each $k$, the leading block also factors as $A_k=L_k{U}_k$, where $L_k$ is unit lower triangular and $U_k$ is upper triangular with diagonal $u_{00},\dots ,u_{k-1,k-1}$. By [L3] and [L4], ${\mathrm{\Delta }}_k=\det (A_k)=\det (L_k)\det (U_k)=1\cdot {\prod }_{i<k}{u}_{ii}\mathrm{\ne }0.$

For the converse, use induction on $n$. When $n=1$, the condition is ${\mathrm{\Delta }}_1=a_{00}\mathrm{\ne }0$, and $A=(1)(a_{00})$ is already a normalised LU factorisation with nonzero diagonal.

Assume $n>1$ and the converse for size $n-1$. Write $A=\left(\begin{array}{cc}a& r^{\mathsf{T}}\\ c& C\end{array}\right)\mathrm{.}$ Then $a={\mathrm{\Delta }}_1\mathrm{\ne }0$. Set $S:=C-c{a}^{-1}{r}^{\mathsf{T}}$. For each $1\le j<n$, apply [L2] to the leading $(j+1)\times (j+1)$ block of $A$: its determinant is $a\det (S_j)$, where $S_j$ is the leading $j\times j$ block of $S$. Hence $\det (S_j)={\mathrm{\Delta }}_{j+1}\mathrm{/}a\mathrm{\ne }0.$ By the induction hypothesis, $S$ has a normalised LU factorisation $S=\tilde{L}\tilde{U}$ with nonzero diagonal.

Step 1.3 and [L2] give $A=\left(\begin{array}{cc}1& 0\\ c{a}^{-1}& \tilde{L}\end{array}\right)\left(\begin{array}{cc}a& r^{\mathsf{T}}\\ 0& \tilde{U}\end{array}\right),$ which is a normalised LU factorisation of $A$. Its upper diagonal is $a$ followed by the diagonal of $\tilde{U}$, so it is nonzero.

Step 1.1 proves the forward implication, and steps 1.2-2.1 prove the reverse implication for every $n\ge 1$. The final sentence is the case $k=n$.

At stage $k$, the matrix $E_k$ is unit lower triangular and differs from $I$ only in column $k$ below the diagonal. Its inverse is $E_k^{-1}=I+{\sum }_{i>k}{m}_{ik}{e}_i{e}_k^{\mathsf{T}}$, because all cross terms vanish: $e_i{e}_k^{\mathsf{T}}{e}_j{e}_k^{\mathsf{T}}=0$ for $i,j>k$. Thus $E_k^{-1}$ is unit lower triangular with strict lower entries exactly the multipliers from column $k$.

Since $U=E_{n-2}\cdots E_{0}A$, one has $A=E_0^{-1}\cdots E_{n-2}^{-1}U$. The product of unit lower triangular matrices is unit lower triangular, and the first time a strict lower entry in column $k$ appears is in $E_k^{-1}$, so later factors do not change that column above the already-created zeros. Hence the strict lower entries of $L:=E_0^{-1}\cdots E_{n-2}^{-1}$ are exactly the stored multipliers $m_{ik}$.

The matrix $U$ is upper triangular by construction, so step 2.1 is a normalised LU factorisation in the sense of [L1]. This is the scalar-entry form of the same elimination mechanism that [L2] packages blockwise. Equivalently, the row operations of elimination store their multipliers in $L$ and leave the final reduced matrix as $U$.

Steps 2.1 and 3.1 prove the claim.

For $n=1$, either $A=[a]$ with $a\mathrm{\ne }0$ or $a=0$. In both cases $I_{1}A=(1)[a]$ is a PLU factorisation, and the zero-pivot case is exactly the singular case.

Assume $n>1$. After an initial row swap, write $PA=\left(\begin{array}{cc}a& r^{\mathsf{T}}\\ c& C\end{array}\right),$ where $\mathrm{\mid }a\mathrm{\mid }$ is maximal in the first column below the swap position. If $a\mathrm{\ne }0$, perform the first elimination step. This produces $\left(\begin{array}{cc}a& r^{\mathsf{T}}\\ 0& S\end{array}\right)$ with multipliers $c\mathrm{/}a$. By the induction hypothesis, the trailing block $S$ has a PLU factorisation. Reassembling the first elimination step and the recursive factorisation gives $P^{\mathrm{\prime }}A=LU$ for suitable $P^{\mathrm{\prime }},L,U$. If $a=0$, then the whole first column is zero, so $PA$ already has the block form $\left(\begin{array}{cc}0& r^{\mathsf{T}}\\ 0& C\end{array}\right)$; by the induction hypothesis $C$ has a PLU factorisation, and adjoining the zero first column yields a PLU factorisation of $PA$, hence of $A$.

If the first failed pivot occurs at stage $k$, then after the previous row swaps and eliminations the active matrix has the form $\left(\begin{array}{cc}\ast & \ast \\ 0& T\end{array}\right),$ where the first column of $T$ is zero. Therefore $T$ is singular because $T{e}_0=0$ for a nonzero standard basis vector $e_0$. The whole active matrix is then singular, hence so are all earlier row permutations and elimination preimages, including $A$.

Claim 1 is steps 1.1-1.2. Step 2.1 proves claim 2, and its contrapositive is the final sentence: an invertible matrix cannot reach a zero pivot under partial pivoting.

In exact arithmetic, each computed upper-triangular entry satisfies a formula of the form ${\widehat{u}}_{ij}=\mathrm{fl}((\widehat{P}A{)}_{ij}-{\sum }_{k<i}{\widehat{\mathrm{\ell }}}_{ik}{\widehat{u}}_{kj}),$ and each computed strict-lower entry satisfies ${\widehat{\mathrm{\ell }}}_{ij}=\mathrm{fl}\left(\frac{(\widehat{P}A{)}_{ij}-{\sum }_{k<j}{\widehat{\mathrm{\ell }}}_{ik}{\widehat{u}}_{kj}}{{\widehat{u}}_{jj}}\right)\mathrm{.}$ By [L2]-[L4], each such floating-point combination is the exact relation for a perturbed entry of $\widehat{P}A$, with perturbation bounded by ${\gamma }_n$ times the corresponding dot-product sum. Collecting those entrywise perturbations into a matrix $E$ gives $\widehat{P}A+E=\widehat{L}\widehat{U}$ and $\mathrm{\mid }{E}_{ij}\mathrm{\mid }\le {\gamma }_n{\sum }_{k\le \min (i,j)}\mathrm{\mid }{\widehat{\mathrm{\ell }}}_{ik}\mathrm{\mid }\mathrm{\mid }{\widehat{u}}_{kj}\mathrm{\mid }\mathrm{.}$

Set $\mathrm{\Delta }A:={\widehat{P}}^{\mathsf{T}}E$. Because $\widehat{P}$ is a permutation matrix, $\mathrm{\mid }\mathrm{\Delta }A\mathrm{\mid }={\widehat{P}}^{\mathsf{T}}\mathrm{\mid }E\mathrm{\mid }\le {\gamma }_n{\widehat{P}}^{\mathsf{T}}(\mathrm{\mid }\widehat{L}\mathrm{\mid }\mathrm{\mid }\widehat{U}\mathrm{\mid })$ by step 1.1. Multiplying the identity of step 1.1 by ${\widehat{P}}^{\mathsf{T}}$ on the left gives $\widehat{P}(A+\mathrm{\Delta }A)=\widehat{L}\widehat{U},$ which is the nearby-factorisation statement.

By [L1], every row of $\mathrm{\mid }\widehat{L}\mathrm{\mid }$ has one diagonal $1$ and at most $n-1$ strict-lower entries of modulus at most $1$, so its row sum is at most $n$. Every entry of $\widehat{U}$ is bounded by $\rho (A){\max }_{p,q}\mathrm{\mid }{a}_{pq}\mathrm{\mid }$ from the definition of pivot growth. Therefore each row sum of $\mathrm{\mid }\widehat{L}\mathrm{\mid }\mathrm{\mid }\widehat{U}\mathrm{\mid }$ is at most $n^2\rho (A){\max }_{p,q}\mathrm{\mid }{a}_{pq}\mathrm{\mid }$, and step 1.1 yields $\mathrm{\parallel }E{\mathrm{\parallel }}_{\mathrm{\infty }}\le n^2{\gamma }_n\rho (A){\max }_{p,q}\mathrm{\mid }{a}_{pq}\mathrm{\mid }\mathrm{.}$ Because row permutations preserve the infinity norm, the same bound holds for $\mathrm{\Delta }A$.

Step 2.1 proves the backward-error identity, and step 2.2 gives the pivot-growth dependent bound.

Suppose $A=L{L}^{\ast }$ as in [L1]. Then $A^{\ast }=L{L}^{\ast }=A$, so $A$ is Hermitian. For any nonzero $x$, one has $x^{\ast }Ax=x^{\ast }L{L}^{\ast }x=(L^{\ast }x{)}^{\ast }(L^{\ast }x)=\mathrm{\parallel }{L}^{\ast }x{\mathrm{\parallel }}_2^2\mathrm{.}$ The diagonal of $L$ is positive, hence nonzero; [L2] and [L3] give $\det L\mathrm{\ne }0$, so $L$ is invertible and $L^{\ast }x\mathrm{\ne }0$ when $x\mathrm{\ne }0$. Therefore $x^{\ast }Ax>0$.

For the converse, use induction on $n$. The case $n=1$ is immediate: $A=[a]$ is Hermitian positive definite exactly when $a>0$, in which case $A=[\sqrt{a}][\sqrt{a}{]}^{\ast }$.

Assume $n>1$ and the converse for size $n-1$. Write $A=\left(\begin{array}{cc}a& r^{\ast }\\ r& C\end{array}\right)\mathrm{.}$ Positive definiteness of $A$ at $x=e_0$ gives $a>0$. Set $\mathrm{\ell }:=r\mathrm{/}\sqrt{a}$ and $S:=C-\mathrm{\ell }{\mathrm{\ell }}^{\ast }$. For nonzero $y\in {\mathbb{C}}^{n-1}$, apply positive definiteness to $x=\left(\genfrac{}{}{0px}{}{-a^{-1}{r}^{\ast }y}{y}\right)$: $0<x^{\ast }Ax=y^{\ast }Sy\mathrm{.}$ Thus $S$ is Hermitian positive definite. By the induction hypothesis, $S=\tilde{L}{\tilde{L}}^{\ast }$ with $\tilde{L}$ lower triangular and positive diagonal.

Define $L:=\left(\begin{array}{cc}\sqrt{a}& 0\\ \mathrm{\ell }& \tilde{L}\end{array}\right)\mathrm{.}$ Direct multiplication gives $L{L}^{\ast }=\left(\begin{array}{cc}a& r^{\ast }\\ r& \mathrm{\ell }{\mathrm{\ell }}^{\ast }+\tilde{L}{\tilde{L}}^{\ast }\end{array}\right)=\left(\begin{array}{cc}a& r^{\ast }\\ r& C\end{array}\right)=A\mathrm{.}$ The matrix $L$ is lower triangular and its diagonal is positive real, so this is a Cholesky factorisation.

For uniqueness, suppose $A=L_1{L}_1^{\ast }=L_2{L}_2^{\ast }$ are two such factorizations. Then $M:=L_2^{-1}{L}_1$ is lower triangular, and $M{M}^{\ast }=I$ because $L_{2}M{M}^{\ast }{L}_2^{\ast }=L_1{L}_1^{\ast }=L_2{L}_2^{\ast }$. A lower-triangular unitary matrix has zero strict-lower entries and diagonal entries of modulus $1$; since the diagonal of $M$ is the ratio of two positive real diagonals, each diagonal entry is also positive real, hence equals $1$. Therefore $M=I$ and $L_1=L_2$.

Step 1.1 proves that a Cholesky factorisation implies Hermitian positive definiteness. Steps 1.2-2.1 prove the converse by induction, and step 3.1 gives uniqueness.

By [L1] and [L2], condition 1 is equivalent to positive definiteness, and positive definiteness is equivalent to existence of a Cholesky factorisation with positive diagonal. Hence conditions 1 and 2 are equivalent.

Suppose $A=L{L}^{\mathsf{T}}$ with $L$ lower triangular and positive diagonal. Its leading $k\times k$ principal block factors as $A_k=L_k{L}_k^{\mathsf{T}}$, where $L_k$ is the leading $k\times k$ block of $L$. By [L3], ${\mathrm{\Delta }}_k=\det (A_k)=\det (L_k{)}^2={\left({\prod }_{i<k}{\mathrm{\ell }}_{ii}\right)}^2={\prod }_{i<k}{\mathrm{\ell }}_{ii}^2\mathrm{.}$ Since each ${\mathrm{\ell }}_{ii}>0$, every ${\mathrm{\Delta }}_k$ is positive.

Step 1.1 proves the equivalence, and step 2.1 identifies the leading principal minors with the squared pivot products.

By [L2] and [L3], the equation $Ax=b$ is $L{L}^{\ast }x=b$. Solve $Ly=b$ by forward substitution and then $L^{\ast }x=y$ by backward substitution. Stepwise uniqueness in [L1] makes both solves unique, and substituting $y=L^{\ast }x$ back into the first equation gives $Ax=b$.

The solve cost is the sum of two triangular solves, so [L1] gives $2(n+\frac{n(n-1)}{2}+\frac{n(n-1)}{2})=2n^2+O(n)$ scalar operations, or $n^2+O(n)$ multiply-add pairs.

In the dense factorisation, column $k$ computes one square root, $n-1-k$ divisions in the column below the pivot, and updates the trailing symmetric submatrix of size $n-1-k$ by a rank-one outer product. That update touches $\frac{(n-1-k)(n-k)}{2}$ stored entries, and each one costs one multiplication and one subtraction, so the update work is $(n-1-k)(n-k)\mathrm{.}$ Summing with $m=n-1-k$ gives ${\sum }_{m=0}^{n-1}m(m+1)=\frac{1}{3}n(n-1)(n+1)=\frac{1}{3}{n}^3+O(n^2),$ which dominates the cost. LU performs the same style of update on the full trailing matrix, giving $2{\sum }_{m=0}^{n-1}{m}^2=\frac{1}{3}n(n-1)(2n-1)=\frac{2}{3}{n}^3+O(n^2)$ scalar operations. Thus Cholesky uses about half the factorisation work.

Step 1.1 proves the solve statement, while steps 1.2 and 2.1 give the operation counts.

For $H_v=I-2v{v}^{\ast }$ with $v^{\ast }v=1$, one has $H_v^{\ast }=I-2v{v}^{\ast }=H_v$ and $H_v^2=I-4v{v}^{\ast }+4v(v^{\ast }v)v^{\ast }=I$. Hence $H_v^{\ast }{H}_v=I$, so $H_v$ is orthogonal or unitary.

In the real case, for $(a,b)\mathrm{\ne }(0,0)$ set $r=(a^2+b^2{)}^{1\mathrm{/}2}$, $c=a\mathrm{/}r$, and $s=b\mathrm{/}r$. Then $c^2+s^2=1$ and $\left(\begin{array}{cc}c& s\\ -s& c\end{array}\right)\left(\genfrac{}{}{0px}{}{a}{b}\right)=\left(\genfrac{}{}{0px}{}{r}{0}\right)\mathrm{.}$ In the complex case, if $a=0$ take $c=0$ and $s=b\mathrm{/}\mathrm{\mid }b\mathrm{\mid }$; otherwise write $a=\varphi \mathrm{\mid }a\mathrm{\mid }$ with $\mathrm{\mid }\varphi \mathrm{\mid }=1$, set $r=(\mathrm{\mid }a{\mathrm{\mid }}^2+\mathrm{\mid }b{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}$, $c=\mathrm{\mid }a\mathrm{\mid }\mathrm{/}r$, and $s=\varphi \bar{b}\mathrm{/}r$. Then $c^2+\mathrm{\mid }s{\mathrm{\mid }}^2=1$ and $\left(\begin{array}{cc}c& s\\ -\bar{s}& c\end{array}\right)\left(\genfrac{}{}{0px}{}{a}{b}\right)=\left(\genfrac{}{}{0px}{}{\varphi r}{0}\right)\mathrm{.}$ Thus a Givens transformation annihilates the second entry while preserving the Euclidean norm.