【机器学习】用QR分解求最小二乘法的最优闭式解

写在前面

今天刷知乎，看到张皓在面试官如何判断面试者的机器学习水平？的回答里面讲到了关于用矩阵的QR分解求解最小二乘法闭式解的问题，碰巧前几天《矩阵分析》课堂上刚讲到QR分解，觉得挺有意思，值得深究，遂产生写本篇博文的动机。另外本人知识水平理解能力有限，如有错漏请各位大侠批评指正。

QR分解

先来看看维基百科上怎么说:

QR分解法是三种将矩阵分解的方式之一。这种方式，把矩阵分解成一个半正交矩阵与一个上三角矩阵的积。QR分解经常用来解线性最小二乘法问题。QR分解也是特定特征值算法即QR算法的基础。

定义

实数矩阵的QR分解是把矩阵A分解为：

这里的Q是正交矩阵（意味着Q^TQ = E）而R是上三角矩阵(即矩阵R的对角线下面的元素全为0)。为了便于理解，这里引用一下daniel-D的图：

这个将A分解成这样的矩阵Q和R的过程就是QR分解，QR分解可用于求解矩阵A的特征值、A的逆等问题。具体参考血影雪梦大神的博文。

QR的求解

QR 分解的实际计算有很多方法，例如 Givens 旋转、Householder 变换，以及 Gram-Schmidt 正交化等等。每一种方法都有其优点和不足。对于QR分解的求解这里不再展开论述，在python中可以用numpy库实现：

import numpy as np
Q, R = np.linalg.qr(A)

线性回归模型

给定数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x m , y m ) } D=\{(\boldsymbol x_1,y_1),(\boldsymbol x_2,y_2),...,(\boldsymbol x_m,y_m)\} D={(x1​,y1​),(x2​,y2​),...,(xm​,ym​)}其中${\mathbf{x}}_i=(x_{i1};x_{i2};...;x_{id}),y_i\in \mathbb{R}$.即样本由$d$个属性描述。“线性回归”(linear regression)试图学得一个线性模型以尽可能准确地预测实值输出标记。
线性回归试图学得的模型为：$$f({\mathbf{x}}_i)={\mathbf{w}}^T{\mathbf{x}}_i+b,使得f({\mathbf{x}}_i)\simeq y_i$$
这称为“多元线性回归”(multivariate linear regression)
可以利用最小二乘法来对参数$\mathbf{w}和b$进行估计。我们知道均方误差是回归任务中最常用的性能度量，即通常所说的损失函数，这里的损失函数可以表示为：
$$J(\mathbf{w},b)=\frac{1}{2}\sum _{i=1}^m(f({\mathbf{x}}_i)-y_i{)}^2=\frac{1}{2}\sum _{i=1}^m(y_i-\mathbf{w}{x}_i-b{)}^2$$
因此我们可试图让均方误差最小化，求出最优闭式解，即
$$({\mathbf{w}}^{\ast },b^{\ast })=\arg mi{n}_{(\mathbf{w},b)}J(\mathbf{w},b)$$
在处理大型数据集时，我们一般都把参数和训练数据进行向量化(vectorization)，即吸入矩阵中，再对矩阵进行计算，数据的向量化能给计算速度带来非常大的提升。具体做法是：对于训练参数，把 $\mathbf{w}和b吸收入向量形式\hat{\mathbf{w}}=(b;\mathbf{w})$，即在列向量$\mathbf{w}$的第一个元素前插入偏置项$b$，此时的参数$\hat{\mathbf{w}}$变成$(d+1)\times 1$的列向量；即
$$\hat{\mathbf{w}}=\left[\begin{array}{c}b\\ w_1\\ ⋮\\ w_d\end{array}\right]$$
对于训练数据集$D$，我们把$D$吸入一个$m\times (d+1)$大小的矩阵$\mathbf{X}$中，矩阵$\mathbf{X}$的第一列用全为1的元素填充，即
$$\mathbf{X}=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{1d}\\ 1& x_{21}& \cdots & x_{2d}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{m1}& \cdots & x_{md}\end{array}\right]=\left[\begin{array}{cc}1& {\mathbf{x}}_1^T\\ 1& {\mathbf{x}}_2^T\\ ⋮& ⋮\\ 1& {\mathbf{x}}_m^T\end{array}\right]$$
相应地，我们也把标记(label)写成向量形式：
$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$$
至此，所有数据已经向量化完。经过向量化后，线性回归模型的损失函数可表示为：
$$J(\mathbf{w},b)=J(\hat{\mathbf{w}})=\frac{1}{2}(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y}{)}^T(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y})$$
将$J(\hat{\mathbf{w}})$对参数$\hat{\mathbf{w}}$进行求导：
$$\frac{\partial J(\hat{\mathbf{w}})}{\partial \hat{\mathbf{w}}}={\mathbf{X}}^T(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y})$$
令导数为零即可得到最优闭式解
$${\hat{\mathbf{w}}}^{\ast }=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$$
重点来了，那么我们怎样才能高效求解这个 ${\hat{\mathbf{w}}}^{\ast }$ 呢？

下面介绍的内容才是本篇文章的核心

用QR分解求解最优闭式解

本文的前面已经介绍过QR分解的定义与求法了，这里我们直接应用。
根据QR分解的定义，我们可以将训练集矩阵$\mathbf{X}$分解为
$$\mathbf{X}=QR$$
代入上面的最优闭式解的式子得到
$$\begin{array}{l}{\hat{\mathbf{w}}}^{\ast }=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}\\ \\ \Rightarrow {\mathbf{X}}^T\mathbf{X}{\hat{\mathbf{w}}}^{\ast }={\mathbf{X}}^T\mathbf{y}\\ \\ \Rightarrow R^T{Q}^{T}QR{\hat{\mathbf{w}}}^{\ast }=R^T{Q}^T\mathbf{y}\\ \\ \Rightarrow R^{T}R{\hat{\mathbf{w}}}^{\ast }=R^T{Q}^T\mathbf{y}\\ \\ \Rightarrow R{\hat{\mathbf{w}}}^{\ast }=Q^T\mathbf{y}\\ \\ \Rightarrow {\hat{\mathbf{w}}}^{\ast }=R^{-1}{Q}^T\mathbf{y}\end{array}$$

至此，我们已经用QR分解求得了最优闭式解的表达式。可以看到，仅需知道Q和R就能求出最小二乘法的最优闭式解，so easy!

看到这里可能有的人就要问了，为什么不直接用${\hat{\mathbf{w}}}^{\ast }=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$进行最优闭式解的求解呢？用上面这种解法有什么好处？
解答这个问题，要从矩阵的条件数说起。

矩阵的条件数

先说一下矩阵的奇异值分解
我们知道对于任意 $m\times n$ 的矩阵$\mathbf{X}$，均可以进行奇异值分解： $\mathbf{X}=U\Sigma V^H$，其中$U$是$m\times m$阶的酉矩阵，$\Sigma$是 $m\times n$阶非负实数对角矩阵，而$V^H$，即$V$的共轭转置，是$n\times n$阶的酉矩阵。$\Sigma$对角线上的元素${\sigma }_i$就是矩阵$\mathbf{X}$的奇异值。

对于条件数，维基百科这样说:

“数值分析中，一个问题的条件数是该数量在数值计算中的容易程度的衡量，也就是该问题的适定性。一个低条件数的问题称为良置的，而高条件数的问题称为病态（或者说非良置）的”

上面这句话可能不太好理解，我们只要记住这一点就行：通常情况下，矩阵的条件数越低越好。由于计算机的计算不是精确的，条件数越高，计算精度的误差对解的影响也越大。
矩阵$\mathbf{X}$的条件数定义为:
$$\mathbf{\kappa }(\mathbf{X})=\Vert {\mathbf{X}}^{-1}\Vert \cdot \Vert \mathbf{X}\Vert$$
若$\Vert \cdot \Vert$是矩阵2范数，则$$\mathbf{\kappa }(\mathbf{X})=\frac{{\mathbf{\sigma }}_{\max }(\mathbf{X})}{{\mathbf{\sigma }}_{\min }(\mathbf{X})}$$
其中${\mathbf{\sigma }}_{\max }(\mathbf{X})$和${\mathbf{\sigma }}_{\min }(\mathbf{X})$分别是矩阵$\mathbf{X}$的极大极小奇异值。
回到上面最优闭式解的求解的问题，如果我们直接拿${\hat{\mathbf{w}}}^{\ast }=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$进行求解，其条件数：
$$\mathbf{\kappa }({\mathbf{X}}^T\mathbf{X})=\mathbf{\kappa }(V{\Sigma }^T{U}^{T}U\Sigma V^H)=\mathbf{\kappa }(V{\Sigma }^2{V}^H)=\mathbf{\kappa }(\mathbf{X}{)}^2$$
也就是说，矩阵${\mathbf{X}}^T\mathbf{X}$的条件数是矩阵$\mathbf{X}$条件数的平方！这和我们想要条件数尽量低的愿望相违背。而使用QR分解将闭式解的求解等价为解$R{\hat{\mathbf{w}}}^{\ast }=Q^T\mathbf{y}$，这就巧妙的绕过了条件数被平方的操作。

实验

我们将这种解法应用到吴恩达的《机器学习》课程上的房价预测的数据集上。这里直接放出代码

import numpy as np
import matplotlib.pyplot as plt



if __name__ == '__main__':

    data = np.loadtxt('ex1data1.txt', delimiter=',')  # load data
    X = data[:, 0]
    y = data[:, 1]  # 注意这里得到的y的shape为(m, ),要把它reshape为(m,1)
    y = np.reshape(y, (y.shape[0], 1))
    
    plt.scatter(X, y, marker='x', c='r')
    # 在X的第一列之前插入元素全为1的列
    X = np.column_stack(
        (np.ones((X.shape[0], 1)), np.reshape(X, (X.shape[0], 1))))
    #QR分解
    Q, R = np.linalg.qr(X)
    w = np.linalg.linalg.inv(R).dot(Q.T).dot(y)
    print(w)
    plt.plot(X[:, 1], np.dot(X, w), '-')
    plt.show()

运行结果

w = [-3.895, 1.193]

完。

Refrences:
《机器学习》周志华
https://www.cnblogs.com/daniel-D/p/3208534.html
https://www.zhihu.com/people/hao-zhang-0214/activities
https://blog.csdn.net/u014046022/article/details/78877369?utm_source=blogxgwz3
https://zh.wikipedia.org/wiki/%E6%9D%A1%E4%BB%B6%E6%95%B0
https://zh.wikipedia.org/wiki/%E9%85%89%E7%9F%A9%E9%98%B5
https://zh.wikipedia.org/wiki/QR%E5%88%86%E8%A7%A3