3.5 线性最小二乘问题的求解方法

3.5.1 正规方程

这里我们考虑超定线性最小二乘问题的求解

定理3.23 设 $A \in \mathbb{R}^{m \times n}$ ( $m \geq n$ ). 则 $x_{*} \in \mathbb{R}^{n}$ 是线性最小二乘问题 (3.1) 的解当且仅当残量 $r = b - Ax_{*}$ 与 $\operatorname{Ran}(A)$ (值域) 正交, 即 $x_{*}$ 是下面的正规方程的解

A ^ {\mathsf {T}} (b - A x) = 0 \quad \text {或} \quad A ^ {\mathsf {T}} A x = A ^ {\mathsf {T}} b. \qquad \qquad \qquad (3. 2 2)

(板书)

证明. 充分性: 设 $x_*$ 是正规方程 (3.21) 的解. 对任意向量 $y \in \mathbb{R}^n$ , 由 $(b - Ax_*) \bot \operatorname{Ran}(A)$ 可知

\begin{array}{l} \| A y - b \| _ {2} ^ {2} = \| (A x _ {*} - b) + A (y - x _ {*}) \| _ {2} ^ {2} \\ = \| A x _ {*} - b \| _ {2} ^ {2} + \| A (y - x _ {*}) \| _ {2} ^ {2} \\ \geq \left\| A x _ {*} - b \right\| _ {2} ^ {2}. \\ \end{array}

因此, $x_*$ 是线性最小二乘问题 (3.1) 的解

必要性：设 $x_{*}$ 是线性最小二乘问题(3.1)的解．用反证法，假定 $z\triangleq A^{\mathsf{T}}(b - Ax_{*})\neq 0$ 取 $y = x_{*} + \alpha z$ ，其中 $\alpha >0$ ，则有

\| A y - b \| _ {2} ^ {2} = \| A x _ {*} - b + \alpha A z \| _ {2} ^ {2} = \| A x _ {*} - b \| _ {2} ^ {2} - 2 \alpha \| z \| _ {2} ^ {2} + \alpha^ {2} \| A z \| _ {2} ^ {2}.

由于 $\| z\| _2 > 0$ ，当 $\alpha$ 充分小时，有 $2\| z\| _2^2 >\alpha \| Az\| _2^2$ ，即上式右端小于 $\| Ax_{*} - b\|_{2}^{2}$ .这与 $x_{*}$ 是最小二乘解相矛盾.所以 $z = 0$ ,即 $A^{\mathsf{T}}(b - Ax_{*}) = 0.$

由定理3.23可知, 求线性最小二乘问题(3.1)的解等价于求正规方程(3.21)的解. 由于

A ^ {\mathsf {T}} b \in \operatorname {R a n} (A ^ {\mathsf {T}}) = \operatorname {R a n} (A ^ {\mathsf {T}} A),

因此正规方程 $A^{\mathsf{T}}Ax = A^{\mathsf{T}}b$ 是相容 (consistent) 的, 即最小二乘解总是存在的. 当 $A$ 非奇异时, 这个解也是唯一的.

定理3.24 设 $A \in \mathbb{R}^{m \times n}$ ( $m > n$ ). 则 $A^{\mathrm{T}}A$ 对称正定当且仅当 $A$ 是列满秩的, 即 $\operatorname{rank}(A) = n$ . 此时, 线性最小二乘问题 (3.1) 的解是唯一的, 其表达式为

x = \left(A ^ {\mathsf {T}} A\right) ^ {- 1} A ^ {\mathsf {T}} b.

最小二乘解的几何含义

根据定理3.23，我们可以把 $b$ 写成

b = A x _ {*} + r, \quad \text {其 中} \quad r \bot \operatorname {R a n} (A). \tag {3.23}

所以 $Ax_{*}$ 就是 $b$ 在 $\operatorname {Ran}(A)$ 上的正交投影，见图3.2

图3.2.最小二乘解的几何描述

最小二乘解可能并不唯一，但分解(3.22)总是唯一的

3.5.2 Cholesky分解法

当 $A$ 列满秩时, 我们就可以使用 Cholesky 分解来求解正规方程, 总的运算量大约为 $mn^2 + \frac{1}{3} n^3 + mn + \mathcal{O}(n^2)$ , 其中大部分的运算量 $(mn^2)$ 是用来计算 $A^\top A$ (由于 $A^\top A$ 对称, 因此只需计算其下三角部分).

通过直接求解正规方程来求解最小二乘问题, 运算量小, 而且简单直观.
但由于 $A^{\mathsf{T}}A$ 的条件数是 $A$ 的条件数的平方, 因此对于病态情形 (即 $A$ 的条件数比较大), 不建议使用该方法.

例3.6 下面的例子说明, 计算 $A^{\mathsf{T}} A$ 可能会损失计算精度: 设

A = \left[ \begin{array}{c c c} 1 & 1 & 1 \\ \varepsilon & & \\ & \varepsilon & \\ & & \varepsilon \end{array} \right],

则

A ^ {\mathsf {T}} A = \left[ \begin{array}{c c c} 1 + \varepsilon^ {2} & 1 & 1 \\ 1 & 1 + \varepsilon^ {2} & 1 \\ 1 & 1 & 1 + \varepsilon^ {2} \end{array} \right].

记 $\varepsilon_{u}$ 为机器精度, 则当 $\varepsilon_{u} < \varepsilon < \sqrt{\varepsilon_{u}}$ 时有 $\varepsilon^{2} < \varepsilon_{u}$ , 由于舍入误差的原因, 通过浮点运算计算得到的 $A^{\mathsf{T}}A$ 是奇异的. 但我们注意到 $A$ 是满秩的.

3.5.3 QR分解法

这里假定 $A \in \mathbb{R}^{m \times n}$ ( $m \geq n$ ) 是满秩的. 设 $A$ 的 QR 分解为 $A = QR$ , 我们用三种不同的方法来推导线性最小二乘问题的解.

(1) 将 $Q$ 的扩充成一个正交矩阵, 记为 $[Q, \hat{Q}] \in \mathbb{R}^{m \times m}$ . 于是有

\left\| A x - b \right\| _ {2} ^ {2} = \left\| [ Q, \hat {Q} ] ^ {\mathsf {T}} (A x - b) \right\| _ {2} ^ {2} = \left\| [ Q, \hat {Q} ] ^ {\mathsf {T}} (Q R x - b) \right\| _ {2} ^ {2}

= \left\| \left[ \begin{array}{c} R x - Q ^ {\mathsf {T}} b \\ - \hat {Q} ^ {\mathsf {T}} b \end{array} \right] \right\| _ {2} ^ {2} = \left\| R x - Q ^ {\mathsf {T}} b \right\| _ {2} ^ {2} + \left\| \hat {Q} ^ {\mathsf {T}} b \right\| _ {2} ^ {2} \geq \left\| \hat {Q} ^ {\mathsf {T}} b \right\| _ {2} ^ {2},

等号成立当且仅当 $R x = Q^{\mathsf{T}} b$ . 所以最小二乘解为

x _ {*} = R ^ {- 1} Q ^ {\intercal} b.

(2) 由于 $Q$ 的列向量构成 $\operatorname{Ran}(A)$ 的一组单位正交基, 所以 $QQ^{\mathsf{T}}$ 和 $I - QQ^{\mathsf{T}}$ 分别是 $\operatorname{Ran}(A)$ 和 $\operatorname{Ran}(A)^{\perp}$ 上的正交投影矩阵. 又 $Ax \in \operatorname{Ran}(A)$ , 所以

\begin{array}{l} \left\| A x - b \right\| _ {2} ^ {2} = \left\| Q Q ^ {\mathsf {T}} (A x - b) \right\| _ {2} ^ {2} + \left\| (I - Q Q ^ {\mathsf {T}}) (A x - b) \right\| _ {2} ^ {2} \\ = \| Q R x - Q Q ^ {\top} b \| _ {2} ^ {2} + \| (I - Q Q ^ {\top}) b \| _ {2} ^ {2} \\ = \| R x - Q ^ {\mathsf {T}} b \| _ {2} ^ {2} + \| (I - Q Q ^ {\mathsf {T}}) b \| _ {2} ^ {2} \\ \geq \| (I - Q Q ^ {\mathsf {T}}) b \| _ {2} ^ {2}, \\ \end{array}

等号成立当且仅当 $R x = Q^{\mathsf{T}} b$ . 所以最小二乘解为

x _ {*} = R ^ {- 1} Q ^ {\mathsf {T}} b.

事实上，由图3.2可知， $x_{*}$ 满足

A x _ {*} = Q Q ^ {\intercal} b,

即 $QRx_{*} = QQ^{\mathsf{T}}b,$ 由此可得 $x_{*} = R^{-1}Q^{\mathsf{T}}b.$

(3) 解正规方程. 由定理 3.23 可知, 最小二乘解为

x _ {*} = \left(A ^ {\mathsf {T}} A\right) ^ {- 1} A ^ {\mathsf {T}} b = \left(R ^ {\mathsf {T}} Q ^ {\mathsf {T}} Q R\right) ^ {- 1} R ^ {\mathsf {T}} Q ^ {\mathsf {T}} b = \left(R ^ {\mathsf {T}} R\right) ^ {- 1} R ^ {\mathsf {T}} Q ^ {\mathsf {T}} b = R ^ {- 1} Q ^ {\mathsf {T}} b.

用基于MGS的QR分解求解最小二乘问题时的运算量大约为 $2mn^{2}$ ，而基于Householder变换的QR分解则只需 $2mn^{2} - 2n^{3} / 3$ (不需要显式计算 $Q$ ，只需将Householder变换直接作用到 $b$ 上即可计算出 $Q^{\mathsf{T}}b$ )。因此当 $m\gg n$ 时，QR分解法的计算量大约是Cholesky分解的两倍。
因此当 $m \gg n$ 时, QR 分解法的计算量大约是 Cholesky 分解法的两倍. 当 $m = n$ 时, 几乎相同.
通常 QR 算法比较稳定, 是求解最小二乘问题的首选方法, 特别是当 $A$ 条件数较大 (病态)时.

3.5.4 奇异值分解法

设 $A \in \mathbb{R}^{m \times n}$ 列满秩, $A = U_n \Sigma V^\top$ 是 $A$ 的降阶奇异值分解, 则

x _ {*} = \left(A ^ {\mathsf {T}} A\right) ^ {- 1} A ^ {\mathsf {T}} b = \left(V \Sigma U _ {n} ^ {\mathsf {T}} U _ {n} \Sigma V ^ {\mathsf {T}}\right) ^ {- 1} V \Sigma U _ {n} ^ {\mathsf {T}} b = \left(V \Sigma^ {- 2} V ^ {\mathsf {T}}\right) V \Sigma U _ {n} ^ {\mathsf {T}} b = V \Sigma^ {- 1} U _ {n} ^ {\mathsf {T}} b.

相比于Cholesky分解法和QR分解，用SVD求解最小二乘问题具有更高的健壮性，但由于需要计算系数矩阵的SVD,运算量远超Cholesky分解法和QR分解.所以只有当系数矩阵秩亏或者接近秩亏时才使用(此时QR分解法可能会失效).

例3.7 分别用三种方法求解最小二乘问题, 比较运算时间.

(LS_3methods.m)

三种方法的运算时间如下 $(A\in \mathbb{R}^{2n\times n},b\in \mathbb{R}^n$ ，以秒为单位)：