10.解析解方法推导线性回归——不容小觑的线性回归算法

引言

线性回归是许多复杂机器学习模型的基础。作为一种基本的机器学习方法，线性回归提供了清晰的思路和工具，通过理解其推导过程，可以更好地掌握机器学习的基本原理和模型设计。

通过阅读本篇博客，你可以：

1.学会如何用解析解的方法推导线性回归的最优解

2.了解如何判定损失函数是凸函数或非凸

一、解析解的推导

通过上一篇9.深入线性回归推导出MSE——不容小觑的线性回归算法-CSDN博客的讲解，我们已经得到了线性回归的损失函数形式，也明确了目标就是最小化损失函数，那么问题就变成了 $\theta$ 什么时候可以使得损失函数最小。

1.最小二乘形式变化

我们已知损失函数公式为 :

$J(\theta ) = \frac{1}{2}\sum_{m}^{i=1}(h_{\theta }x_{i}-y_{i})^{2}$

我们将损失函数变化一个形式，变为以下：

$\frac{1}{2}(X\theta - y)^{T}(X\theta -y)$

其中 $X$ 为变量集，即所有 $x_{i}$ 样本行程的 $m$ 行 $n$ 列的样本矩阵。 $\theta$ 即我们要求的最优解，它是一个 $m$ 行 1 列的矩阵。这个公式是如何变化而来的呢？

首先原损失函数中的 $h_{\theta }$ 是线性回归模型中的预测函数， $h_{\theta }x_{i}$ 用来表示预测值 $\hat{y_{i}}$ 。所以我们可以得出以下结论：

$\hat{y_{i}} = h_{\theta }x_{i} = x_{i}\theta$

$\hat{y} = h_{\theta }X = X\theta$

得到这个结论之后，我们回归到公式本身:

$J(\theta ) = \frac{1}{2}\sum_{m}^{i=1}(h_{\theta }x_{i}-y_{i})^{2}$

$\Rightarrow J(\theta ) = \frac{1}{2} \sum_{m}^{i=1}(h_{\theta }x_{i}-y_{i})(h_{\theta }x_{i}-y_{i})$

将上述公式代入，又由于矩阵的性质，我们需要将其中一项转置，这里就相当于一个长度为 $m$ 的向量乘以它自己，说白了就是对应位置相乘相加。

所以我们的公式变为:

$J(\theta ) = \frac{1}{2}(X\theta - y)^{T}(X\theta - y)$

由矩阵运算的基本性质

可继续推出公式：

$J(\theta) = \frac{1}{2}((X\theta)^{T} - y^{T} )(X\theta - y)$

$\Rightarrow J(\theta) = \frac{1}{2} (\theta^{T}X^{T} - y^{T})(X\theta - y)$

最终，我们得到：

$J(\theta) = \frac{1}{2}(\theta^{T}X^{T}X\theta - \theta^{T}X^{T}y-y^{T}X\theta+y^{T}y)$

2.推导出模型的解析解形式

假使我们开着小车，从下图中寻找最优解。为了便于理解，我们假设存在横轴表示 $\theta$ ，存在纵轴表示 loss损失，曲线是 loss function。

我们把最小二乘看成是一个函数曲线，最优解一定是驻点中某个极小值(驻点顾名思义就是小车可以停驻的点)。从图中我们可以看出，驻点的特定是梯度全为0(梯度：函数在某点上的切线的斜率)。

所以要求出 $\theta$ 的解析解形式，我们就可以通过把函数的一阶导函数推导出来，再使其的值为0以求出 $\theta$ 。依据以下求导公式：

我们能将公式进行推导：

${J}'(\theta) = \frac{1}{2}\left [ {(\theta^{T}X^{T}X\theta)}' - {(\theta^{T}X^{T}y)}'-{(y^{T}X\theta)}' + {(y^{T}y)}'\right ]$

由于 $X$ 和 $y$ 是已知的， $\theta$ 是我们要求的答案。所以和 $\theta$ 没关系的部分在求导时可以忽略不计，继续推导为以下公式：

$\Rightarrow {J}'(\theta) = \frac{1}{2}[2X^{T}X\theta-X^{T}y-(y^{T}X)^T]$

$\Rightarrow {J}'(\theta) = \frac{1}{2}[2X^{T}X\theta-2X^{T}y]$

$\Rightarrow {J}'(\theta) = X^{T}X\theta - X^{T}y$

然后我们设置导函数为0，去进一步解出来驻点对应的 $\theta$ 值为多少：

$0 = X^{T}X\theta - X^{T}y$

$\Rightarrow X^{T}X\theta = X^{T}y$

由于矩阵与逆矩阵相乘可以得到单位矩阵，所以我们最终可以求出 $\theta$ 的解析解形式(解析解为方程的解析式，是方程的精确解，能在任意精度下满足方程)：

$\theta = (X^{T}X)^{-1}X^{T}y$

这样，我们有数据集 $X$ ， $y$ 时，就可以将数据代入上面解析解公式，去直接求出对应的 $\theta$ 值了。比如我们可以设想 $X$ 为 $m$ 行 $n$ 列的矩阵， $y$ 为 $m$ 行 1 列的列向量。 $X^{T}$ 是 $n$ 行 $m$ 列的，所以 $X^{T}X$ 就是 $n$ 行 $n$ 列的矩阵。又因为矩阵求逆形状不变，再次乘以 $X^{T}$ 后变为 $n$ 行 $m$ 列的矩阵。最后乘以 $y$ ，结果 $\theta$ 就是 $n$ 行 1 列的列向量！

二、判断损失函数是否为凸函数

对于求解最优解而言，判断一个损失函数是否为凸函数是极其重要的。如果一个损失函数是凸函数，那么局部最优解即为全局最优解，这是因为在凸函数上没有局部极小值的存在，所有的局部极小值都位于全局最小值处。

如上图所示，左上和右下为非凸函数，左下和右上为凸函数。在非凸函数中，有很多条极值点，我们无法直接得到最优解。对于二次可微的函数，我们可以通过判断黑塞矩阵(hessian matrix)是否为半正定的来进行判断，所以我们要对目标函数在点 $x$ 处的二阶偏导数进行求解。

对于我们的式子来说，就是在导函数的基础上再次对 $\theta$ 进行求偏导，由于 $X^{T}y$ 对 $\theta$ 的导数为0，所以再次求偏导后的答案为 $X^{T}X$ 。所谓的正定就是答案的特征值全为正数，而半正定无非就是特征值大于等于0。这里我们对损失函数求二阶导的黑塞矩阵是 $X^{T}X$ ，自己和自己做点乘，所以答案一定是半正定的。

在此处我们不用深入去讨论数学推导的证明。在机器学习中，损失函数往往是凸函数，在深度学习中的损失函数往往是非凸函数。并且在实际应用当中，我们并不要求找到全局最优解，只要模型适用。机器学习的特点就是不强调模型 100% 正确，而是有价值的，堪用的。

三、代码实战求解线性回归算法模型

经过大片的理论讲解，相信大家已经对线性回归模型的实现有了深刻的认识，接下来我们就要通过代码的形式来实战求解线性回归模型。

1.导入需要使用的库

import numpy as np
import matplotlib.pyplot as plt

我们需要使用numpy模块进行矩阵之间的运算，最后用matplotlib模块中的绘图功能绘制 $X$ 与 $y$ 的关系图。

2.定义样本集

# 回归，有监督监督机器学习，X，y
X = 2 * np.random.rand(100,1)
y = 5 * 4 * X + np.random.randn(100,1)

这里的 $X$ 是所有 $x_{i}$ 组成的100行1列的矩阵。 $y$ 是真实值，5是偏置(截距)，4是 $X$ 的权重，后面的 np.random.randn(100,1) 是100行1列以正态分布形成的误差矩阵。

3.实现解析解公式求解模型

# 为了求解W0截距项,我们给X矩阵加上一列全为1的X0
X_b = np.c_[np.ones((100,1)),X]

上述代码是通过 np.c_ 的方式将截距项恒为1的权重拼接到 $X$ 矩阵中。

# 实现解析解的公式来求解θ
θ = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print(θ)
"""[[5.21509616]
    [3.77011339]]
"""

我们将得到的 $\theta$ 的解析解公式通过代码去实现，np.linalg.inv() 是numpy模块中用来计算逆矩阵的函数，X_b.T 是变量 X_b 的转置，.dot() 是进行点乘运算(对numpy模块的讲解在专栏前面的文章当中7.科学计算模块Numpy（4）ndarray数组的常用操作（二）_ndarray逐元素相加-CSDN博客)。我们通过以上代码就可以表示公式：

$\theta = (X^{T}X)^{-1}X^{T}y$

输出 $\theta$ 之后我们可以看到，截距项与权重都是相当接近真实情况，但由于误差的存在，我们不可能得到真实值，只能拟合数据得到最优解。

4.使用模型去预测

X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
print(X_new_b)
"""[[1. 0.]
    [1. 2.]]
"""

我们定义一个2行1列的矩阵作为要预测样本的自变量。再加上截距恒为1的项，就会形成一个2行2列的矩阵。

y_predict = X_new_b.dot(θ)
print(y_predict)
"""[[ 5.21509616]
    [12.75532293]]
"""

随后，我们计算预测值，也就是 $\hat{y}$ 。

y_predict = X_new_b.dot(θ)
print(y_predict)
"""[[ 5.21509616]
    [12.75532293]]
"""

只要使用刚刚拼接完的矩阵点乘 $\theta$ 就可以得到预测值了。

5.绘图

plt.plot(X_new, y_predict, 'r-')
plt.plot(X, y, 'b.')
plt.axis([0, 2, 0, 15])
plt.show()

这边我们使用到了matplotlib模块中的绘图功能，成功绘制了下方的坐标图。其中，横轴表示了输入 $x$ 的值，纵轴表示了 $y$ 的值，红线代表了整体的函数，蓝色的点则是真实值的分布情况。我们可以发现，红色的直线尽可能地穿过了蓝色的点，这就是我们一直说的线性回归模型。

总结

这篇博客讲述了模型解析解的推导原理以及代码实现。希望可以对大家起到作用，谢谢。

关注我，内容持续更新（后续内容在作者专栏《从零基础到AI算法工程师》）！！！