最小二乘计算时用QR分解的目的是什么?

在进行最小二乘算法的计算时，有时候会用到QR分解，你知道为什么吗？最小二乘不是有闭式解吗，QR分解的意义何在？本文告诉你！

1. 最小二乘算法
最小二乘算法是一种以最小化误差平方和为目标来拟合一组系数的方法。设有$N$个目标值$y_1,...,y_N$，每个目标值由$M$个变量$x_1,...,x_M$决定，即$y_i={\sum }_m{\omega }_m{x}_{i,m}$。在已知$y_1,...,y_N$和$x_1,...,x_M$的情况下，${\omega }_m$该如何取值才能使$y_i-{\sum }_m{\omega }_m{x}_{i,m}$误差最小？
我们将上面的问题用矩阵表示：
$$\{\begin{array}{c}{y}_1={\omega }_1{x}_{1,1}+{\omega }_2{x}_{1,2}+...+{\omega }_M{x}_{1,M}\\ y_2={\omega }_1{x}_{2,1}+{\omega }_2{x}_{2,2}+...+{\omega }_M{x}_{2,M}\\ ⋮\\ y_N={\omega }_1{x}_{N,1}+{\omega }_2{x}_{N,2}+...+{\omega }_M{x}_{N,M}\end{array}\to \mathbf{y}=\mathbf{X}\omega$$
其中，$\mathbf{y}=[y_1,...,y_N{]}^T$表示目标向量，$\omega =[{\omega }_1,...,{\omega }_N{]}^T$表示系数向量，
$$\mathbf{X}=\left[\begin{array}{ccc}{x}_{1,1}& ...& x_{1,M}\\ x_{2,1}& ...& x_{2,M}\\ ⋮& ...& ⋮\\ x_{N,1}& ...& x_{N,M}\end{array}\right]$$
当$N\ge M$时，方程数目大于未知参数数目，上述方程为超定方程，可以用最小二乘法求解。将上述问题转化为
$${\omega }^{\ast }=\arg \underset{\omega }{\min }||\mathbf{y}-\mathbf{X}\omega |{|}^2$$
上式的求解过程如下:
$$\begin{gathered} ||\mathbf{y}-\mathbf{X}\omega |{|}^2=(\mathbf{y}-\mathbf{X}\omega {)}^T(\mathbf{y}-\mathbf{X}\omega )={\mathbf{y}}^{\mathbf{T}}\mathbf{y}-2{\mathbf{y}}^{\mathbf{T}}\mathbf{X}\omega +{\omega }^{\mathbf{T}}{\mathbf{X}}^{\mathbf{T}}\mathbf{X}\omega =J(\omega ) \\ \downarrow \\ \frac{dJ(\omega )}{d\omega }=-2{\mathbf{X}}^{\mathbf{T}}\mathbf{y}+2{\mathbf{X}}^{\mathbf{T}}\mathbf{X}\omega =0 \\ \downarrow \\ {\omega }^{\ast }=({\mathbf{X}}^{\mathbf{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathbf{T}}\mathbf{y} \end{gathered}$$
可见，最小二乘法有闭式解，可以根据上式直接求得最优系数。

2. QR分解
QR分解是将一个矩阵$\mathbf{A}$分解为一个正交矩阵$\mathbf{Q}$和一个上三角矩阵$\mathbf{R}$的乘积，即$\mathbf{A}=\mathbf{QR}$。其中，正交矩阵满足$\mathbf{Q}{\mathbf{Q}}^{\mathbf{T}}=\mathbf{I}$，上三角矩阵$\mathbf{R}$的形式如下，它在对角线下边的元素全为0。
$$\mathbf{R}=\left[\begin{array}{cccc}{r}_{1,1}& r_{1,2}& ...& r_{1,M}\\ 0& r_{2,2}& ...& r_{2,M}\\ ⋮& ...& ⋮& ⋮\\ 0& 0& ...& r_{M,M}\end{array}\right]$$
QR分解可以很方便地进行线性方程组的求解，说明如下：
$$\begin{gathered} \mathbf{Ax}=\mathbf{b}\to \mathbf{QRx}=\mathbf{b}\to \mathbf{Rx}={\mathbf{Q}}^{-1}\mathbf{b}={\mathbf{Q}}^{\mathbf{T}}\mathbf{b}=\mathbf{c} \\ \downarrow \\ \left[\begin{array}{cccc}{r}_{1,1}& r_{1,2}& ...& r_{1,M}\\ 0& r_{2,2}& ...& r_{2,M}\\ ⋮& ...& ⋮& ⋮\\ 0& 0& ...& r_{M,M}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_M\end{array}\right]=\left[\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_M\end{array}\right] \\ \downarrow \\ \{\begin{array}{c}{r}_{1,1}{x}_1+r_{1,2}{x}_2+...+r_{1,M}{x}_M=c_1\\ r_{2,2}{x}_2+...+r_{2,M}{x}_M=c_2\\ ⋮\\ r_{M,M}{x}_M=c_M\end{array} \end{gathered}$$

3. 条件数
条件数常用于衡量线性方程组的稳定性，条件数越大，矩阵稳定性越差，越接近于病态(奇异矩阵)。线性方程组$\mathbf{Ax}=\mathbf{b}$中，我们称$\mathbf{b}$为观测向量，$\mathbf{A}$为系数矩阵，$\mathbf{x}$为待求向量。我们期望观测向量$\mathbf{b}$和系数矩阵$\mathbf{A}$的微小变化不会引起解向量$\mathbf{x}$的巨大变化，否则若观测向量中有些许噪声就会导致解向量的极大误差。换言之，我们希望线性方程是具备一定的抗噪能力的。首先考察待求向量$\mathbf{x}$受观测向量$\mathbf{b}$的波动的影响。
$$\begin{gathered} \mathbf{Ax}=\mathbf{b}\to \mathbf{A}(\mathbf{x}+\mathbf{\Delta x})=\mathbf{b}+\mathbf{\Delta b}\to \mathbf{A\Delta x}=\mathbf{\Delta b}\to \mathbf{\Delta x}={\mathbf{A}}^{-1}\mathbf{b}\to ||\mathbf{\Delta x}||\le ||{\mathbf{A}}^{-1}|||\mathbf{b}||| \\ \mathbf{Ax}=\mathbf{b}\to ||\mathbf{b}||\le ||\mathbf{A}||||\mathbf{x}|| \\ \downarrow \\ \frac{||\mathbf{\Delta x}||}{||\mathbf{x}||}\le (||\mathbf{A}||||{\mathbf{A}}^{-1}||)\frac{||\mathbf{\Delta b}||}{||\mathbf{b}||} \end{gathered}$$
下面考虑待求向量$\mathbf{x}$受系数矩阵$\mathbf{A}$的波动的影响：
$$\begin{gathered} (\mathbf{A}+\mathbf{\Delta A})(\mathbf{x}+\mathbf{\Delta x})=b\to \mathbf{\Delta x}=-{\mathbf{A}}^{-1}\mathbf{\Delta A}(\mathbf{x}+\mathbf{\Delta x})\to ||\Delta x||\le ||{\mathbf{A}}^{-1}||||\mathbf{\Delta A}||||\mathbf{x}+\mathbf{\Delta x}|| \\ \downarrow \\ \frac{||\mathbf{\Delta x}||}{||\mathbf{x}+\mathbf{\Delta x}||}\le (||\mathbf{A}||||{\mathbf{A}}^{-1}||)\frac{||\mathbf{\Delta A}||}{||\mathbf{A}||} \end{gathered}$$
在上面的表达式中，$\frac{||\mathbf{\Delta b}||}{||\mathbf{b}||}$表示观测向量的变化率，$\frac{||\mathbf{\Delta A}||}{||\mathbf{A}||}$表示系数向量的变化率, $\frac{||\mathbf{\Delta x}||}{||\mathbf{x}||}$表示待求向量的变化率。条件数定义为$cond(\mathbf{A})=||\mathbf{A}||||{\mathbf{A}}^{-1}||$，显然条件数越大，待求向量的变化受观测向量和系数矩阵变化的影响越大，方程越不稳定。
当上述推导过程中的范数||||取$l_2$范数时，条件数$cond(\mathbf{A})=\frac{{\sigma }_{max}}{{\sigma }_{min}}$，其中${\sigma }_{max}$和${\sigma }_{min}$分别表示矩阵$\mathbf{A}$的最大和最小奇异值。显然，$cond(\mathbf{A})$是一个大于或等于1的正数，当矩阵奇异时，${\sigma }_{min}=0$，此时条件数无穷大。当条件数不是无穷大，但很大时，矩阵接近奇异，此时矩阵的行向量或列向量的线性相关性很强。此外，由条件数的定义可知，正交矩阵或酉矩阵的条件数为1。
对矩阵$\mathbf{A}$进行奇异值分解有$\mathbf{A}=\mathbf{V\Sigma }{\mathbf{V}}^{\mathbf{H}}$，对${\mathbf{A}}^{\mathbf{H}}\mathbf{A}$进行奇异值分解有${\mathbf{A}}^{\mathbf{H}}\mathbf{A}=\mathbf{V}{\mathbf{\Sigma }}^2{\mathbf{V}}^{\mathbf{H}}$，所以${\mathbf{A}}^{\mathbf{H}}\mathbf{A}$的最大最小奇异值分别是矩阵$\mathbf{A}$的最大和最小奇异值的平方，所以$cond(\mathbf{A}=\mathbf{V\Sigma }{\mathbf{V}}^{\mathbf{H}})=\frac{{\sigma }_{max}^2}{{\sigma }_{min}^2}=cond(\mathbf{A}{)}^2$，即矩阵${\mathbf{A}}^{\mathbf{H}}\mathbf{A}$的条件数是矩阵$\mathbf{A}$的条件数的平方。

4. QR分解在最小二乘计算中的应用
从最小二乘的表达式$\mathbf{x}=({\mathbf{A}}^{\mathbf{T}}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\mathbf{T}}\mathbf{b}$中，我们可以认为待求向量是方程${\mathbf{A}}^{\mathbf{T}}\mathbf{Ax}={\mathbf{A}}^{\mathbf{T}}\mathbf{b}$的解。从上面的分析中，我们知道矩阵${\mathbf{A}}^{\mathbf{H}}\mathbf{A}$的条件数是矩阵$\mathbf{A}$的条件数的平方，所以以矩阵${\mathbf{A}}^{\mathbf{H}}\mathbf{A}$为系数矩阵的方程相比以$\mathbf{A}$为系数矩阵的方程稳定性更差。如果将矩阵$\mathbf{A}$进行QR分解，并带入上述线性方程中有
$$(\mathbf{QR}{)}^{\mathbf{T}}(\mathbf{QR})\mathbf{x}=(\mathbf{QR}{)}^{\mathbf{T}}\mathbf{b}\to {\mathbf{R}}^{\mathbf{T}}{\mathbf{Q}}^{\mathbf{T}}\mathbf{QRx}={\mathbf{R}}^{\mathbf{T}}{\mathbf{Q}}^{\mathbf{T}}\mathbf{b}\to {\mathbf{R}}^{\mathbf{T}}\mathbf{Rx}={\mathbf{R}}^{\mathbf{T}}{\mathbf{Q}}^{\mathbf{T}}\mathbf{b}\to \mathbf{Rx}={\mathbf{Q}}^{\mathbf{T}}\mathbf{b}\to \mathbf{x}={\mathbf{R}}^{-1}{\mathbf{Q}}^{\mathbf{T}}\mathbf{b}$$
即此时将原方程中的系数矩阵${\mathbf{A}}^{\mathbf{H}}\mathbf{A}$转化成了系数矩阵$\mathbf{R}$。由于$cond(\mathbf{Q})=1$,$cond(\mathbf{A})=cond({\mathbf{Q}}^{\mathbf{T}}\mathbf{A})=cond({\mathbf{Q}}^{\mathbf{T}}\mathbf{QR})=cond(\mathbf{R})$，所以经过QR分解后，将线性方程组的条件数由$cond({\mathbf{A}}^{\mathbf{T}}\mathbf{A})$转化为了$cond(\mathbf{A})$，可提高方程的稳定性。

5. 计算量对比

6. 仿真验证