Least Squares and QR

Overdetermined Systems

$Ax = b$ with $A \in \mathbb{R}^{m \times n}$, $m > n$, has no exact solution in general. E.g. fitting $m$ points with fewer than $m$ basis functions instead of interpolating (to avoid overfitting),

$$f(t) = c_1 f_1(t) + \cdots + c_n f_n(t), \qquad V = \begin{bmatrix} f_1(t_1) & \cdots & f_n(t_1) \\ \vdots & & \vdots \\ f_1(t_m) & \cdots & f_n(t_m) \end{bmatrix} \in \mathbb{R}^{m \times n}, \qquad Vc \approx y$$

where the fit depends linearly on $c$. Least squares minimizes the sum of squares of the residuals $y_i - f(t_i)$,

$$R(c_1, \ldots, c_n) = \sum_{i=1}^m [y_i - f(t_i)]^2 = r^T r = \|r\|_2^2, \qquad r = y - Vc$$

i.e. in general $\mathrm{argmin}_{x \in \mathbb{R}^n} \|b - Ax\|_2^2$.

Some nonlinear fits become linear after a transformation, minimizing the residuals of $\log y$ rather than of $y$:

• $g(t) = a_1 e^{a_2 t}$: $\log y \approx \log a_1 + a_2 t = c_1 + c_2 t$,
• $g(t) = a_1 t^{a_2}$: $\log y \approx \log a_1 + a_2 \log t$, so fit $z = \log y$ against $s = \log t$.

Normal Equations

$x$ solves the linear least squares problem iff

$$A^T(Ax - b) = 0 \iff A^T A x = A^T b$$

(expand $\|A(x + y) - b\|_2^2 = \|Ax - b\|_2^2 + 2y^T A^T(Ax - b) + \|Ay\|_2^2$). This turns the $m \times n$ problem into an $n \times n$ one. $A^T A$ is symmetric, and if $A^T A z = 0$ then $0 = z^T A^T A z = \|Az\|_2^2$, so $A^T A$ is singular only if the columns of $A$ are linearly dependent, and otherwise it is positive definite. Then

$$x = A^+ b, \qquad A^+ = (A^T A)^{-1} A^T$$

Algorithm $N = A^T A$, $z = A^T b$, solve $Nx = z$ by Cholesky.

With $\kappa(A) = \|A\| \cdot \|A^+\|$ ($\infty$ if $\mathrm{rank}\, A < n$), $\kappa(A^T A) = \kappa(A)^2$ in the 2-norm, so the normal equations square the condition number and are very sensitive for ill-conditioned $A$.

Orthogonality

$u, v$ are orthogonal if $u^T v = 0$, and $q_1, \ldots, q_k$ are orthonormal if also $q_i^T q_i = 1$. For orthogonal $q_1, q_2$,

$$\|q_1 - q_2\|^2 = (q_1 - q_2)^T(q_1 - q_2) = \|q_1\|^2 + \|q_2\|^2$$

so adding and subtracting them is well-conditioned. For a matrix $Q$, $Q^T Q$ is diagonal iff its columns are orthogonal, and $Q^T Q = I$ iff they are orthonormal (an ONC matrix). For $Q \in \mathbb{R}^{n \times k}$ ONC and $x \in \mathbb{R}^k$,

$$Q^T Q = I, \qquad \|Qx\|_2 = \|x\|_2, \qquad \|Q\|_2 = 1$$

A square ONC matrix is orthogonal: $Q^{-1} = Q^T$, $Q^T$ is orthogonal, $\kappa_2(Q) = 1$, $\|AQ\|_2 = \|A\|_2$, and $QU$ is orthogonal for $U$ orthogonal.

QR Factorization

$A = QR$ with $Q \in \mathbb{R}^{m \times m}$ orthogonal and $R \in \mathbb{R}^{m \times n}$ upper triangular. Since the last $m - n$ rows of $R$ are zero, keeping only the first $n$ columns of $Q$ gives the thin QR,

$$A = \hat{Q}\hat{R}, \qquad \hat{Q} \in \mathbb{R}^{m \times n} \text{ ONC}, \quad \hat{R} \in \mathbb{R}^{n \times n} \text{ upper triangular}$$

If $A$ is not rank deficient, the normal equations become

$$A^T A x = A^T b \implies \hat{R}^T \hat{R} x = \hat{R}^T \hat{Q}^T b \implies \hat{R} x = \hat{Q}^T b$$

an $n \times n$ triangular system, stable since it does not square $\kappa$. Gram-Schmidt computes $\hat{Q}$ but is not stable.

Householder Reflections

Like LU, introduce zeros below the diagonal column by column, but with orthogonal transformations. For a vector $z$, find $P$ with $Pz$ nonzero only in its first entry,

$$Pz = \pm \|z\|_2 e_1, \qquad v = \|z\|_2 e_1 - z, \qquad P = I - 2\frac{vv^T}{v^T v}$$

($P = I$ if $v = 0$). $P$ is symmetric and orthogonal.

Algorithm Take $z$ as the first column of $A$, find $P$, and set $A \leftarrow PA$. Then ignore the first row and column, take $z$ as the first column of the remaining submatrix, and repeat. At step $k$,

$$Q_k = \begin{bmatrix} I_{k-1} & 0 \\ 0 & P_k \end{bmatrix}, \qquad P_k \in \mathbb{R}^{(m-k+1) \times (m-k+1)}$$

is orthogonal, and $Q_n Q_{n-1} \cdots Q_1 A = R$, so $Q = (Q_n \cdots Q_1)^T$. $Q^T = Q_n(Q_{n-1}(\cdots(Q_1 I)))$ is built by applying the same steps to $I$, without forming the $P$'s. About $2mn^2 - \frac{2}{3}n^3$ flops.