QR, Normal Equations and Cholesky

Definition Let $A \in \mathbb{F}^{m \times n}$, where $\mathbb{F}$ is $\mathbb{R}$ or $\mathbb{C}$, have full column rank, $\operatorname{rank} A = n$, which forces $m \ge n$. Then all $p = \min(m, n) = n$ singular values are positive, $\sigma_1 \ge \cdots \ge \sigma_n > 0$, the square matrix $A^{*}A$ is invertible, and the pseudoinverse is $A^{+} = (A^{*}A)^{-1}A^{*}$. It satisfies $A^{+}A = I_n$, and its singular values are $1/\sigma_n \ge \cdots \ge 1/\sigma_1$, so $\lVert A^{+} \rVert_2 = 1/\sigma_n$. The two-norm condition number of $A$ is

For a nonsingular square matrix the pseudoinverse is $A^{-1}$, and this is the condition number of the earlier lessons. Because $m \ge n$, the extreme singular values are the largest and the smallest stretch of $A$, that is $\sigma_1 = \max \lVert Ax \rVert_2$ and $\sigma_n = \min \lVert Ax \rVert_2$ over unit vectors $x$, so $\kappa_2(A) \ge 1$. For a right-hand side $b = Ax \ne 0$ in the column space, a change $\delta b$ moves the least-squares solution by $\delta x = A^{+}\delta b$. Since $\lVert \delta x \rVert_2 \le \lVert \delta b \rVert_2 / \sigma_n$ and $\lVert b \rVert_2 \le \sigma_1 \lVert x \rVert_2$,

For a general $b$ the sensitivity of the least-squares solution also depends on the angle between $b$ and the column space, and when $A$ itself is perturbed and the residual is not small it can involve $\kappa_2(A)^2$. So $\kappa_2(A)$ measures the directional scaling of $A$, not every error bound for every data set. The full-column-rank hypothesis is part of the definition. Dividing by the smallest positive singular value of a rank-deficient matrix describes a different problem, and an arbitrarily small change of $A$ can raise its rank and make $\lVert A^{+} \rVert_2$ jump. Here and below a tuple such as $(1, 1, 1)$ stands for a column vector.

Definition Let $A \in \mathbb{F}^{m \times n}$ have full column rank, where $\mathbb{F}$ is $\mathbb{R}$ or $\mathbb{C}$, and let its columns be $a_1, \dots, a_n$. The Gram matrix $A^{*}A$ has entries $(A^{*}A)_{ij} = \langle a_i, a_j \rangle = a_i^{*}a_j$. It is Hermitian positive definite, since $x^{*}A^{*}Ax = \lVert Ax \rVert_2^2 > 0$ for $x \ne 0$. Its eigenvalues, which are also its singular values, are $\sigma_1^2 \ge \cdots \ge \sigma_n^2 > 0$, and

For the proof take a thin SVD $A = U\Sigma V^{*}$, where $U \in \mathbb{F}^{m \times n}$ has orthonormal columns, $\Sigma = \operatorname{diag}(\sigma_1, \dots, \sigma_n)$ and $V \in \mathbb{F}^{n \times n}$ is unitary. Since $U^{*}U = I_n$,

which is a unitary diagonalization with positive diagonal entries and therefore also an SVD of $A^{*}A$. The ratio of its extreme singular values is $\sigma_1^2/\sigma_n^2$. The identity is exact and concerns the matrix of the normal equations. It does not claim that every normal-equations computation loses exactly twice as many digits as a QR computation. It says that the system being solved has condition number $\kappa_2(A)^2$, which exceeds $\kappa_2(A)$ whenever $\kappa_2(A) > 1$, and in floating point the entries of $A^{*}A$ are also rounded before any solve begins.

Definition Let $A \in \mathbb{F}^{m \times n}$ and $b \in \mathbb{F}^m$, where $\mathbb{F}$ is $\mathbb{R}$ or $\mathbb{C}$. A vector $\hat{x}$ minimizes $\lVert Ax - b \rVert_2$ if and only if the residual $\hat{r} = b - A\hat{x}$ is orthogonal to every column of $A$, that is $A^{*}\hat{r} = 0$, or

These are the normal equations. They always have a solution, because the orthogonal projection of $b$ onto the column space has the form $A\hat{x}$. If $A^{*}\hat{r} = 0$, then for every $x$ the vectors $\hat{r}$ and $A(\hat{x} - x)$ are orthogonal, and Pythagoras gives

So $\hat{x}$ is a minimizer, and $x$ is another minimizer exactly when $Ax = A\hat{x}$, in which case $x$ also solves the normal equations. If $A$ has full column rank, then $A^{*}A$ is Hermitian positive definite and the least-squares solution $\hat{x} = (A^{*}A)^{-1}A^{*}b = A^{+}b$ is unique. The normal equations are mathematically exact. Their numerical price is that they replace $A$ by a matrix with condition number $\kappa_2(A)^2$ and round the entries of $A^{*}A$ and $A^{*}b$ along the way.

Definition Let $M \in \mathbb{F}^{n \times n}$ be Hermitian positive definite, where $\mathbb{F}$ is $\mathbb{R}$ or $\mathbb{C}$. For $\mathbb{F} = \mathbb{R}$ this means that $M$ is symmetric and $x^{\mathsf{T}}Mx > 0$ for every $x \ne 0$. Then there is exactly one lower triangular $L \in \mathbb{F}^{n \times n}$ with positive real diagonal entries such that

This is the Cholesky factorization $M = R^{*}R$ of the earlier lesson, written with $L = R^{*}$. A Cholesky solve of $Mx = b$ computes $L$, solves $Ly = b$ by forward substitution and then solves $L^{*}x = y$ by back substitution. The pivots are positive, so no row exchanges are needed, and the factorization takes about half the arithmetic of a general LU factorization. The hypotheses matter. If $M$ is Hermitian but not positive definite, no such factorization exists, and the computation reaches a step where the number under the square root is zero or negative. Cholesky fits a system whose matrix is given as positive definite. Applied to the normal equations of a least-squares problem, it solves a system with condition number $\kappa_2(A)^2$ and inherits the loss described above, however well the triangular steps themselves behave.

Definition Let $A \in \mathbb{F}^{n \times n}$ be nonsingular, where $\mathbb{F}$ is $\mathbb{R}$ or $\mathbb{C}$. Gaussian elimination with partial pivoting works through the columns $k = 1, \dots, n - 1$. At step $k$ it chooses, among the entries of column $k$ on or below the diagonal of the current matrix, one of largest absolute value, swaps its row into row $k$, and subtracts multiples $\ell_{ik} = a_{ik}/a_{kk}$ of row $k$ from the rows below. Because $A$ is nonsingular, a nonzero candidate always exists, and every multiplier satisfies $\lvert \ell_{ik} \rvert \le 1$. The result is

with $P$ a permutation matrix, $L$ unit lower triangular with entries of absolute value at most $1$, and $U$ upper triangular and nonsingular. To solve $Ax = b$, apply the same row swaps to $b$, then solve $Ly = Pb$ and $Ux = y$. Partial pivoting controls the multipliers at each step, but it does not keep the later entries small in general. Entries can grow by a factor up to $2^{n-1}$ relative to the largest entry of $A$, and some matrices reach that factor. It also cannot repair an ill-conditioned matrix, because $\kappa(A)$ belongs to the problem, not to the method.

Definition Rank-aware method choice selects a solver from the stated task, its structure and what is known about its rank, rather than from the shortest formula. For real or complex data the standard choices are these.

• For a general nonsingular square system $Ax = b$, use Gaussian elimination with partial pivoting.
• For a square system whose matrix is given as Hermitian positive definite, use a Cholesky solve.
• For least squares with $A \in \mathbb{F}^{m \times n}$ of full column rank, use Householder QR. A reduced factorization $A = Q_1R$, with $Q_1 \in \mathbb{F}^{m \times n}$ having orthonormal columns and $R \in \mathbb{F}^{n \times n}$ upper triangular and nonsingular, gives the solution of $R\hat{x} = Q_1^{*}b$ by back substitution, and in exact arithmetic the minimal residual satisfies $\lVert b - A\hat{x} \rVert_2^2 = \lVert b \rVert_2^2 - \lVert Q_1^{*}b \rVert_2^2$. In floating point this difference loses digits when the residual is small, and computing the vector $b - Q_1(Q_1^{*}b)$ directly is more accurate. The method works with $A$ itself, whose condition number is $\kappa_2(A)$, and never forms $A^{*}A$.
• When the rank is in doubt, use column-pivoted QR, which can help expose independent columns, or the SVD $A = \sum_{i=1}^{p} \sigma_i u_i v_i^{*}$ with $p = \min(m, n)$, together with a tolerance $\tau$ set by the size of the data errors and the units. The numerical rank is the number of singular values $\sigma_i > \tau$, and directions with $\sigma_i \le \tau$ are dropped or regularized.

One way to drop them is the truncated SVD solution

Even a backward stable method only returns the exact solution of a nearby problem, so its forward error still depends on the conditioning of the problem itself. No factorization recovers information that the data never resolved.