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

写在前面
QR分解
- 定义
- QR的求解
线性回归模型
用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={(x1,y1),(x2,y2),...,(xm,ym)}D=\{(\boldsymbol x_1,y_1),(\boldsymbol x_2,y_2),...,(\boldsymbol x_m,y_m)\}D={(x1,y1),(x2,y2),...,(xm,ym)}其中xi=(xi1;xi2;...;xid),yi∈R\boldsymbol x_i=(x_{i1};x_{i2};...;x_{id}),y_i\in\mathbb Rxi=(xi1;xi2;...;xid),yi∈R.即样本由ddd个属性描述。“线性回归”(linear regression)试图学得一个线性模型以尽可能准确地预测实值输出标记。
线性回归试图学得的模型为：f(xi)=wTxi+b,使得f(xi)≃yif(\boldsymbol x_i)=\boldsymbol w^T\boldsymbol x_i+b, 使得f(\boldsymbol x_i)\simeq y_if(xi)=wTxi+b,使得f(xi)≃yi
这称为“多元线性回归”(multivariate linear regression)
可以利用最小二乘法来对参数w和b\boldsymbol w和bw和b进行估计。我们知道均方误差是回归任务中最常用的性能度量，即通常所说的损失函数，这里的损失函数可以表示为：
J(w,b)=12∑i=1m(f(xi)−yi)2=12∑i=1m(yi−wxi−b)2J(\boldsymbol w,b)=\frac{1}{2} \sum_{i=1}^m (f(\boldsymbol x_i)-y_i)^2 = \frac{1}{2} \sum_{i=1}^m(y_i-\boldsymbol wx_i-b)^2 J(w,b)=21i=1∑m(f(xi)−yi)2=21i=1∑m(yi−wxi−b)2
因此我们可试图让均方误差最小化，求出最优闭式解，即
(w∗,b∗)=arg⁡min(w,b)J(w,b)(\boldsymbol w^*,b^*) = \arg min_{(\boldsymbol w,b)} J(\boldsymbol w,b) (w∗,b∗)=argmin(w,b)J(w,b)
在处理大型数据集时，我们一般都把参数和训练数据进行向量化(vectorization)，即吸入矩阵中，再对矩阵进行计算，数据的向量化能给计算速度带来非常大的提升。具体做法是：对于训练参数，把 w和b吸收入向量形式w^=(b;w)\boldsymbol w和b吸收入向量形式\hat \boldsymbol w=(b; \boldsymbol w)w和b吸收入向量形式w^=(b;w)，即在列向量w\boldsymbol ww的第一个元素前插入偏置项bbb，此时的参数w^\hat \boldsymbol ww^变成(d+1)×1(d+1)\times 1(d+1)×1的列向量；即
w^=[bw1⋮wd]\hat \boldsymbol w= \begin{bmatrix} b \\ w_1 \\ \vdots \\ w_d \end{bmatrix} w^=⎣⎢⎢⎢⎡bw1⋮wd⎦⎥⎥⎥⎤
对于训练数据集DDD，我们把DDD吸入一个m×(d+1)m\times(d+1)m×(d+1)大小的矩阵X\boldsymbol XX中，矩阵X\boldsymbol XX的第一列用全为1的元素填充，即
X=[1x11⋯x1d1x21⋯x2d⋮⋮⋱⋮1xm1⋯xmd]=[1x1T1x2T⋮⋮1xmT]\boldsymbol X=\begin{bmatrix} 1 &x_{11} & \cdots & x_{1d} \\ 1 &x_{21} & \cdots & x_{2d} \\ \vdots & \vdots &\ddots & \vdots \\ 1 &x_{m1} & \cdots & x_{md} \end{bmatrix}= \begin{bmatrix} 1 &\boldsymbol x_1^T \\ 1 &\boldsymbol x_2^T \\ \vdots & \vdots \\ 1 &\boldsymbol x_m^T \end{bmatrix} X=⎣⎢⎢⎢⎡11⋮1x11x21⋮xm1⋯⋯⋱⋯x1dx2d⋮xmd⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡11⋮1x1Tx2T⋮xmT⎦⎥⎥⎥⎤
相应地，我们也把标记(label)写成向量形式：
y=[y1y2⋮ym]\boldsymbol y= \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_m \end{bmatrix} y=⎣⎢⎢⎢⎡y1y2⋮ym⎦⎥⎥⎥⎤
至此，所有数据已经向量化完。经过向量化后，线性回归模型的损失函数可表示为：
J(w,b)=J(w^)=12(Xw^−y)T(Xw^−y)J(\boldsymbol w, b) =J(\hat \boldsymbol w) =\frac{1}{2}(\boldsymbol X \hat \boldsymbol w-\boldsymbol y)^T(\boldsymbol X \hat \boldsymbol w - \boldsymbol y) J(w,b)=J(w^)=21(Xw^−y)T(Xw^−y)
将J(w^)J(\hat\boldsymbol w)J(w^)对参数w^\hat\boldsymbol ww^进行求导：
∂J(w^)∂w^=XT(Xw^−y)\frac{\partial J(\hat\boldsymbol w)}{\partial \hat\boldsymbol w}=\boldsymbol X^T(\boldsymbol X \hat\boldsymbol w - \boldsymbol y) ∂w^∂J(w^)=XT(Xw^−y)
令导数为零即可得到最优闭式解
w^∗=(XTX)−1XTy\hat\boldsymbol w^* = (\boldsymbol X^T \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol y w^∗=(XTX)−1XTy
重点来了，那么我们怎样才能高效求解这个 w^∗\hat\boldsymbol w^*w^∗ 呢？

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

用QR分解求解最优闭式解

本文的前面已经介绍过QR分解的定义与求法了，这里我们直接应用。
根据QR分解的定义，我们可以将训练集矩阵X\boldsymbol XX分解为
X=QR\boldsymbol X = QR X=QR
代入上面的最优闭式解的式子得到
w^∗=(XTX)−1XTy⇒XTXw^∗=XTy⇒RTQTQRw^∗=RTQTy⇒RTRw^∗=RTQTy⇒Rw^∗=QTy⇒w^∗=R−1QTy\begin{array}{lcl} \quad \quad \hat\boldsymbol w^* = (\boldsymbol X^T \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol y\\ \\ \Rightarrow \quad \boldsymbol X^T \boldsymbol X \hat\boldsymbol w^* = \boldsymbol X^T \boldsymbol y\\ \\ \Rightarrow \quad R^T Q^T Q R \hat\boldsymbol w^* = R^T Q^T \boldsymbol y\\ \\ \Rightarrow \quad R^T R \hat \boldsymbol w^*=R^T Q^T \boldsymbol y\\ \\ \Rightarrow \quad R \hat \boldsymbol w^* = Q^T \boldsymbol y \\ \\ \Rightarrow \quad \hat\boldsymbol w^* = R^{-1} Q^T \boldsymbol y \end{array}w^∗=(XTX)−1XTy⇒XTXw^∗=XTy⇒RTQTQRw^∗=RTQTy⇒RTRw^∗=RTQTy⇒Rw^∗=QTy⇒w^∗=R−1QTy

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

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

矩阵的条件数

先说一下矩阵的奇异值分解
我们知道对于任意 m×nm\times nm×n 的矩阵X\boldsymbol XX，均可以进行奇异值分解： X=UΣVH\boldsymbol X =U \Sigma V^HX=UΣVH，其中UUU是m×mm \times mm×m阶的酉矩阵，Σ\SigmaΣ是 m×nm \times nm×n阶非负实数对角矩阵，而VHV^HVH，即VVV的共轭转置，是n×nn \times nn×n阶的酉矩阵。Σ\SigmaΣ对角线上的元素σi\sigma_iσi就是矩阵X\boldsymbol XX的奇异值。

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

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

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

实验

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

import numpy as np
import matplotlib.pyplot as pltif __name__ == '__main__':data = np.loadtxt('ex1data1.txt', delimiter=',')  # load dataX = 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