Numerical Analysis Lecture (IV): Solving Linear Systems and Matrix Computations Part II

... views ... visitors ... total

中文版

It is best to read Numerical Analysis Lecture (IV): Solving Linear Systems and Matrix Computations Part I first. This part continues with the second half of solving linear systems: the Cholesky factorization of symmetric positive definite matrices, and the effects of condition numbers, perturbations, and rounding errors on computed solutions.

---

4.3 The Cholesky Method

For a general invertible matrix, Gaussian elimination without pivoting may fail; and, for reasons of numerical stability, searching for pivots is usually a sensible choice. For the important class of positive definite matrices, however, Gaussian elimination without pivoting can always be carried out in a stable way.

Definition 4.3.1: A real matrix $A\in\mathbb{R}^{n\times n}$ is called positive definite if

\[A=A^T,\qquad x^TAx>0\quad \forall x\in\mathbb{R}^n\setminus\{0\}.\]

It is called positive semidefinite if

\[A=A^T,\qquad x^TAx\ge 0\quad \forall x\in\mathbb{R}^n.\]

More generally, a complex matrix $A\in\mathbb{C}^{n\times n}$ is called positive definite if

\[A=A^H,\qquad x^HAx>0\quad \forall x\in\mathbb{C}^n\setminus\{0\}.\]

It is called positive semidefinite if

\[A=A^H,\qquad x^HAx\ge 0\quad \forall x\in\mathbb{C}^n.\]

Here $A^H=(\overline a_{ji})_{1\le i\le n,1\le j\le n}$, where the bar denotes complex conjugation.

Positive definite matrices occur frequently in applications, for example in the numerical solution of elliptic partial differential equations such as Laplace’s equation and parabolic partial differential equations such as the heat equation.

A positive definite matrix is necessarily invertible.

The Cholesky method provides an efficient variant of Gaussian elimination for linear systems with positive definite matrices. It is based on the following observation.

Theorem 4.3.2: Let $A\in\mathbb{R}^{n\times n}$ be positive definite. Then there exists a unique lower triangular matrix $L$ with positive diagonal entries, that is, $l_{ii}>0$, such that

\[LL^T=A\]

holds. This is called the Cholesky factorization.

In addition, $A$ has a unique triangular factorization

\[\widetilde L\widetilde R=A,\]

where

\[\widetilde L=LD^{-1},\qquad \widetilde R=DL^T,\qquad D=\operatorname{diag}(l_{11},\ldots,l_{nn}).\]

This factorization is produced by Gaussian elimination without pivoting.

The theorem can be proved by complete induction on $n$; the proof is omitted here.

The Cholesky factorization $LL^T=A$ can be obtained by solving the following $\frac{n(n+1)}2$ equations. By symmetry, it is enough to consider the lower triangular part, including the diagonal:

\[a_{ij}=\sum_{k=1}^j l_{ik}l_{jk}, \qquad j\le i,\quad i=1,\ldots,n. \tag{4.11}\]

The entries of $L$ can therefore be computed column by column, in the order

\[l_{11},\ldots,l_{n1},\ l_{22},\ldots,l_{n2},\ \ldots,\ l_{nn}.\]

For the first column of $L$, setting $j=1$ gives

\[a_{11}=l_{11}^2,\qquad \text{so } l_{11}=\sqrt{a_{11}},\] \[a_{i1}=l_{i1}l_{11},\qquad \text{so } l_{i1}=a_{i1}/l_{11}.\]

Solving successively for $l_{ij}$, $i=j,\ldots,n$, gives the following algorithm.

Algorithm 4.3.3: Cholesky method for computing the factorization $LL^T=A$

For $j=1,\ldots,n$:

\[l_{jj} = \sqrt{ a_{jj}-\sum_{k=1}^{j-1}l_{jk}^2 }.\]

For $i=j+1,\ldots,n$:

\[l_{ij} = \frac{ a_{ij}-\sum_{k=1}^{j-1}l_{ik}l_{jk} }{l_{jj}}.\]

Figure 4-3: Cholesky factorization exploits the symmetric positive definite structure and constructs only a lower-triangular factor.

Positive definiteness keeps the quantity under the square root positive; if a nonpositive quantity appears, the algorithm also provides a valid test for positive definiteness.

Remark 4.3.4 Cholesky’s method has several useful properties:

• Because it exploits symmetry, it requires only $O(n^3/6)$ operations apart from $n$ square roots. This is approximately half the operation count of unstructured Gaussian elimination.

• Equation (4.11) gives

\[|l_{ij}|\le \sqrt{a_{ii}}, \qquad j\le i,\quad i=1,\ldots,n.\]

Thus, the entries of $L$ do not become excessively large. This is an important reason for the numerical stability of the Cholesky method.

• Cholesky’s method is one of the most effective general tests for positive definiteness. Extend Algorithm 4.3.3 as follows:

\[a=a_{jj}-\sum_{k=1}^{j-1}l_{jk}^2.\]

If $a\le 0$, stop: $A$ is not positive definite. Otherwise, set

\[l_{jj}=\sqrt a.\]

4.4 Error Estimates and the Effect of Rounding Errors

When describing direct methods for solving linear systems, we have so far assumed that all input data are exact and that no rounding errors occur during the computation. This is unrealistic, because rounding errors can have a significant effect, especially for large systems.

4.4.1 Error Estimates for Perturbed Systems

We first ask how much the solution of a linear system can change when the matrix and right-hand side are perturbed. Consider

\[Ax=b\]

and the perturbed system

\[(A+\Delta A)\widetilde x=b+\Delta b,\]

where $\Delta A$ and $\Delta b$ are “small.”

How small is $x-\widetilde x$?

This question is important in practice:

• It estimates how sensitive the solution is to perturbations in the matrix and right-hand side.
• A computed approximate solution, for example one obtained from an implementation of Gaussian elimination, is the exact solution of a system such as

\[A\widetilde x=b+\Delta b, \qquad \Delta b=A\widetilde x-b.\]

Therefore, the easily computed residual $\Delta b=A\widetilde x-b$ can be used to derive a bound for the unknown error $|x-\widetilde x|$.

The matrix condition number describes this effect of perturbations.

To measure $x-\widetilde x$, $\Delta b$, and $\Delta A$, we need a notion of length for vectors and matrices.

Definition 4.4.1: A vector norm on $\mathbb{R}^n$ is a map

\[x\in\mathbb{R}^n\mapsto \|x\|\in[0,\infty[\]

that satisfies:

a) $|x|=0$ if and only if $x=0$;

b) for every $\alpha\in\mathbb{R}$ and every $x\in\mathbb{R}^n$,

\[\|\alpha x\|=\lvert\alpha\rvert\,\|x\|;\]

c) for all $x,y\in\mathbb{R}^n$, the triangle inequality

\[\|x+y\|\le \|x\|+\|y\|\]

holds.

We now introduce matrix norms. Let $|\cdot|$ be any norm on $\mathbb{R}^n$. The corresponding matrix norm on $\mathbb{R}^{n\times n}$ is defined by

\[\|A\|:=\sup_{\|x\|=1}\|Ax\| =\sup_{x\ne 0}\frac{\|Ax\|}{\|x\|}, \qquad A\in\mathbb{R}^{n\times n}. \tag{4.12}\]

It is called the matrix norm induced by the vector norm $|\cdot|$.

It satisfies the following properties as well:

a) $|A|=0$ if and only if $A=0$;

b) for every $\alpha\in\mathbb{R}$ and every $A\in\mathbb{R}^{n\times n}$,

\[\|\alpha A\|=\lvert\alpha\rvert\,\|A\|;\]

c) for all $A,B\in\mathbb{R}^{n\times n}$, the triangle inequality

\[\|A+B\|\le \|A\|+\|B\|.\]

In addition, (4.12) gives the useful inequality

d) for all $x\in\mathbb{R}^n$ and all $A\in\mathbb{R}^{n\times n}$,

\[\|Ax\|\le \|A\|\,\|x\|\]

which is called the compatibility condition.

e) for all $A,B\in\mathbb{R}^{n\times n}$,

\[\|AB\|\le \|A\|\,\|B\|\]

which is called submultiplicativity.

Example 4.4.2

\[\|x\|_2=\sqrt{x^Tx} \quad\text{induces}\quad \|A\|_2=\sqrt{\lambda_{\max}(A^TA)}.\] \[\|x\|_1=\sum_{i=1}^n |x_i| \quad\text{induces}\quad \|A\|_1=\max_{j=1,\ldots,n}\sum_{i=1}^n |a_{ij}|\]

This is called the column-sum norm.

\[\|x\|_\infty=\max_{i=1,\ldots,n}|x_i| \quad\text{induces}\quad \|A\|_\infty=\max_{i=1,\ldots,n}\sum_{j=1}^n |a_{ij}|\]

This is called the row-sum norm.

We can now introduce the matrix condition number mentioned above.

Definition 4.4.3: Let $A\in\mathbb{R}^{n\times n}$ be invertible, and let $|\cdot|$ be an induced matrix norm. Then

\[\operatorname{cond}(A)=\|A\|\,\|A^{-1}\|\]

is called the condition number of $A$ with respect to this matrix norm.

The following result can be proved.

Figure 4-4: The condition number is not the algorithmic error itself; it measures the problem's sensitivity to perturbations.

The Hilbert matrix example shows that even a residual caused by tiny rounding errors can be amplified by the condition number into a noticeable solution error.

Theorem 4.4.4 (Effect of perturbations in the matrix and right-hand side)
Let $A\in\mathbb{R}^{n\times n}$ be invertible, let $b,\Delta b\in\mathbb{R}^n$ with $b\ne 0$, and let $\Delta A\in\mathbb{R}^{n\times n}$ satisfy

\[\|\Delta A\|<\frac1{\|A^{-1}\|},\]

where $|\cdot|$ is a matrix norm induced by an arbitrary norm on $\mathbb{R}^n$. If $x$ solves

\[Ax=b\]

and $\widetilde x$ solves

\[(A+\Delta A)\widetilde x=b+\Delta b,\]

then

\[\frac{\|\widetilde x-x\|}{\|x\|} \le \frac{\operatorname{cond}(A)} {1-\operatorname{cond}(A)\|\Delta A\|/\|A\|} \left( \frac{\|\Delta A\|}{\|A\|} +\frac{\|\Delta b\|}{\|b\|} \right).\]

Proof
For simplicity, consider only the case $\Delta A=0$. Subtracting the unperturbed system from the perturbed system gives

\[A(\widetilde x-x)=\Delta b,\]

and therefore

\[\|\widetilde x-x\| =\|A^{-1}\Delta b\| \le \|A^{-1}\|\,\|\Delta b\|.\]

Since

\[\|b\|=\|Ax\|\le \|A\|\,\|x\|,\]

we have

\[\frac1{\|x\|}\le \frac{\|A\|}{\|b\|}.\]

Consequently,

\[\frac{\|\widetilde x-x\|}{\|x\|} \le \|A\|\,\|A^{-1}\|\frac{\|\Delta b\|}{\|b\|}.\]

Thus, the condition number determines how sensitive the solution is to perturbations in the matrix and right-hand side.

4.4.2 Rounding Error Analysis for Gaussian Elimination

A basic, though lengthy, estimate of the amplification of rounding errors in Gaussian elimination gives the following result.

Theorem 4.4.5: Let $A\in\mathbb{R}^{n\times n}$ be invertible. On a computer with machine precision $\mathrm{eps}$, apply Gaussian elimination to $A$ using a pivoting strategy that guarantees $\lvert l_{ij}\rvert\le 1$, such as partial or complete pivoting. Then the computed $\overline L,\overline R$ satisfy

\[\overline L\,\overline R=PAQ+F, \qquad |f_{ij}|\le \frac{2j\,\overline a\,\mathrm{eps}}{1-\mathrm{eps}}.\]

Here $P,Q$ are the permutations generated by the pivoting strategy, and

\[\overline a=\max_k \overline a_k, \qquad \overline a_k=\max_{i,j}|a^{(k)}_{ij}|. \tag{4.13}\]

If an approximate solution $\overline x$ of $Ax=b$ is computed from $\overline L,\overline R$ by forward and backward substitution, then there exists a matrix $E$ such that

\[(A+E)\overline x=b, \qquad |e_{ij}| \le \frac{2(n+1)\mathrm{eps}}{1-n\mathrm{eps}} (|\overline L|\,|\overline R|)_{ij} \le \frac{2(n+1)\mathrm{eps}}{1-n\cdot\mathrm{eps}}\,n\overline a.\]

Here

\[|\overline L|=(|\overline l_{ij}|), \qquad |\overline R|=(|\overline r_{ij}|).\]

Proof
See Stoer [5].

Remark 4.4.6
Theorem 4.4.4 can now also be used to estimate the relative error of the approximate solution $\overline x$.

Effect of Pivoting Strategies

The size of $\overline a$ in (4.13) depends on the pivoting strategy. One can prove:

• Partial pivoting:

\[\overline a_k\le 2^k\max_{i,j}|a_{ij}|.\]

This bound can be attained, but it is usually too pessimistic. In practice, one almost always observes

\[\overline a_k\le 10\max_{i,j}|a_{ij}|.\]

• Partial pivoting for tridiagonal matrices:

\[\overline a_k\le 2\max_{i,j}|a_{ij}|.\]

• Complete pivoting:

\[\overline a_k\le f(k)\max_{i,j}|a_{ij}|, \qquad f(k)=k^{1/2}\left(2\,3^{1/2}\cdots k^{1/(k-1)}\right)^{1/2}.\]

The function $f(n)$ grows very slowly. So far, no example has been found satisfying

\[\overline a_k\ge (k+1)\max_{i,j}|a_{ij}|.\]

Example 4.4.7
Consider the Hilbert matrix $H^n=(h^n_{ij})\in\mathbb{R}^{n\times n}$, where

\[h^n_{ij}=\frac1{i+j-1}, \qquad i,j\in\{1,\ldots,n\}.\]

This matrix is well known to be badly conditioned. For example,

\[\operatorname{cond}(H^5)\approx 9.4\cdot 10^5\]

with respect to $|\cdot|_\infty$, and

\[\|H^5\|_\infty\approx 2.3, \qquad \|(H^5)^{-1}\|_\infty\approx 4.1\cdot 10^5.\]

When Gaussian elimination with partial pivoting is applied to it, $\overline a=1$.

For $n=5$ and $\mathrm{eps}=10^{-16}$, Theorem 4.4.5 gives

\[|e_{ij}| \le \frac{2(n+1)\mathrm{eps}}{1-n\cdot\mathrm{eps}}\,n\overline a = \frac{6\cdot 10^{-15}}{1-5\cdot 10^{-16}} \approx 6\cdot 10^{-15}.\]

Therefore,

\[\|E\|_\infty\le 3\cdot 10^{-14}.\]

Theorem 4.4.4 gives

\[\frac{\|\widetilde x-x\|_\infty}{\|x\|_\infty} \le \frac{\operatorname{cond}(A)} {1-\operatorname{cond}(A)\|E\|_\infty/\|A\|_\infty} \frac{\|E\|_\infty}{\|A\|_\infty} = \frac{\|A^{-1}\|_\infty\|E\|_\infty} {1-\|A^{-1}\|_\infty\|E\|_\infty}\] \[\approx \frac{4.1\cdot 10^5\cdot 3\cdot 10^{-14}} {1-4.1\cdot 10^5\cdot 3\cdot 10^{-14}} \approx 1.23\cdot 10^{-8}.\]

Thus, rounding errors alone “lose” about half of the significant digits. For larger $n$, the rounding error grows rapidly and soon makes the result unusable.

For larger $n$, Theorem 4.4.4 is no longer applicable because

\[\|\Delta A\|>\frac1{\|A^{-1}\|}\]

holds.

---

Return to Numerical Analysis Lecture (IV): Solving Linear Systems and Matrix Computations Part I.

Abbreviations and Notation

• SPD: symmetric positive definite, meaning $A=A^T$ and $x^TAx>0$.
• Cholesky factorization: the factorization $A=LL^T$ of a symmetric positive definite matrix.
• $\operatorname{cond}(A)$: the condition number, which describes how sensitive a linear system is to perturbations in its input.
• pivot: a pivot element; partial pivoting searches rows, while complete pivoting searches both rows and columns.

Source, Copyright, and Usage Notes

This article is organized from Chapter 4 of the locally saved TU Darmstadt 2016 Mathematik 4 ET/3Inf lecture file Skript-Mathe4ET-3Inf-2016-Kap4-5.pdf, with reference to the Chinese translation draft Skript-Mathe4ET-3Inf-2016-Kap4.zh.md in the same local directory. It is published for personal study, translation, and knowledge organization. The English wording, supplementary explanations, and remade figures in this article do not represent the original authors or an official position.

My organization, English wording, supplementary explanations, and remade figures in this article may be used for non-commercial study, discussion, and citation, provided that the author and original material source are credited. Because parts of this article are based on translation and organization of course lecture notes, the original lecture notes and any materials they may contain should still be governed by their original authors, course pages, and relevant authorization statements. For commercial use, systematic republication, publication, or large-scale adaptation, please confirm the authorization status of the original materials first.

If there are omissions or errors in translation, formulas, terminology, or interpretation, or if a relevant rights holder considers the use of any content inappropriate, please contact me and I will handle or remove it promptly.