阅读本文需要的背景知识点：矩阵求导、一丢丢编程知识

一、引言

前面介绍了两种二元分类算法——感知器算法、口袋算法，这些算法解决的都是分类的问题，但是现实中更多的是例如预测某一地区的房价、银行该给某个人多少额度的信用卡、今天应该买卖多少股票等等这种最后得到一个具体数值结果的问题，这种类型的问题在机器学习中统一被称为回归问题。

回归分析在统计学中是研究多组变量间关系的方法，在机器学习中也是应用广泛，下面介绍其中一种算法模型 - 线性回归¹（Linear Regression）

二、模型介绍

在介绍模型前，我们先来看一个例子，假设有某地统计了不同工作年限人群的收入情况（模拟数据），如下表所示：

工作年限	平均月工资
1年	1598元
2年	3898元
3年	6220元
4年	7799元
5年	10510元

由上面表格中的数据可以画出每隔一年当地的平均月收入的图像：

由上图可以直观的看到这些点好像都在一条直线上或在其附近，我们根据直觉可以判断出当地的平均月收入与工作年限似乎有一定的线性关系，我们只需要找到一条直线，那么我们就能预测出当地有 6 年工资年限的人的平均月收入了。寻找这样一条直线的过程就称为线性回归分析。

定义：给定一些随机样本点 { x1, y1 }，{ x2, y2 }，…，找到一个超平面（单变量时为一条直线，两变量时为一个平面）去拟合这些样本点。线性方程如下：

（1）超平面函数方程一般形式

（2）同感知器算法一样，将 b 看作第零个 w，将超平面函数方程简写成两个向量 w 和 x 的点积的形式

$\begin{array}{rcc} y & =b+w_{1} x_{1}+w_{2} x_{2}+\cdots+w_{M} x_{M} \\ & =\quad w^{T} x \end{array}$

单变量为直线、两变量为平面、多变量为超平面

三、算法步骤

那么如何从一堆样本点中找到最佳的那个超平面呢？在上面工作年限与平均月工资的例子中，好像可以画很多直线去拟合这些点。

就像上面两条直线一样，哪条直线直观上拟合的更好呢？应该可以通过图中虚线看到，每个样本点与直线 B 的距离相对直线 A 来说要远，直线 A 明显是比直线 B 拟合的更好的，这说明可以通过样本点与直线的距离来判断拟合的好坏情况。

假设有 N 个样本点{ x, y }，每个样本点的自变量有 M 个{ x1, x2, … }，则可以定义所有样本点与这个超平面的距离之和为拟合这些样本点的代价函数，其中 w 为 M 维列向量，x 也为M 维列向量，y 为实数：

$\operatorname{Cost}(w)=\sum_{i=1}^{N}\left|w^{T} x_{i}-y_{i}\right|$

由于函数中存在绝对值，将其改写成平方的形式，这在几何中被称为欧几里得距离²。

$\operatorname{Cost}(w)=\sum_{i=1}^{N}\left(w^{T} x_{i}-y_{i}\right)^{2}$

直观上我们只需要让代价函数的值最小，也就是所有样本点到超平面的欧几里得距离之和最小，其对应的 w 则为这个超平面的权重系数。

$w=\underset{w}{\operatorname{argmin}}\left(\sum_{i=1}^{N}\left(w^{T} x_{i}-y_{i}\right)^{2}\right)$

[argmin](https://en.wikipedia.org/wiki/Arg_max#Arg_min)³ 函数表示当括号内的函数方程值最小时返回此时的变量

基于上面函数的最小化来进行求解的方法被称为最小二乘法，由于代价函数是一个凸函数⁴，根据凸函数的性质可知其局部最小值即全局最小值，可以直接求得 w 的最佳解析解，其中 X 为 N x M 矩阵，y 为 N 维列向量：

$w=\left(X^{T} X\right)^{-1} X^{T} y$

$X=\left[\begin{array}{c} x_{1}^{T} \\ x_{2}^{T} \\ \vdots \\ x_{N}^{T} \end{array}\right]=\left[\begin{array}{cccc} X_{11} & X_{12} & \cdots & X_{1 M} \\ X_{21} & X_{22} & \cdots & X_{2 M} \\ \vdots & \vdots & \ddots & \vdots \\ X_{N 1} & X_{N 2} & \cdots & X_{N M} \end{array}\right] \quad y=\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{N} \end{array}\right)$

四、原理证明

线性回归代价函数为凸函数

凸函数是一个定义在某个向量空间的凸子集C上的实值函数f，而且对于凸子集C中任意两个向量x1，x2满足下式：

$f\left(\frac{x_{1}+x_{2}}{2}\right) \leq \frac{f\left(x_{1}\right)+f\left(x_{2}\right)}{2}$

不等式左边：

$\operatorname{Cost}\left(\frac{w_{1}+w_{2}}{2}\right)=\sum_{i=1}^{N}\left[\left(\frac{w_{1}+w_{2}}{2}\right)^{T} x_{i}-y_{i}\right]^{2}$

不等式右边：

$\frac{\operatorname{Cost}\left(w_{1}\right)+\operatorname{Cost}\left(w_{2}\right)}{2}=\frac{\sum_{i=1}^{N}\left(w_{1}^{T} x_{i}-y_{i}\right)^{2}+\sum_{i=1}^{N}\left(w_{2}^{T} x_{i}-y_{i}\right)^{2}}{2}$

（1）不等式左边乘以 2

（2）将 2 移入连加运算内，括号写出三项

（3）将平方展开成 6 项

$\begin{aligned} 2 \operatorname{Cost}\left(\frac{w_{1}+w_{2}}{2}\right) &=2 \sum_{i=1}^{N}\left[\left(\frac{w_{1}+w_{2}}{2}\right)^{T} x_{i}-y_{i}\right]^{2} & (1) \\ &=\sum_{i=1}^{N} 2\left(\frac{w_{1}^{T} x_{i}}{2}+\frac{w_{2}^{T} x_{i}}{2}-y_{i}\right)^{2} & (2) \\ &=\sum_{i=1}^{N}\left(\frac{w_{1}^{T} x_{i} x_{i}^{T} w_{1}}{2}+\frac{w_{2}^{T} x_{i} x_{i}^{T} w_{2}}{2}+w_{1}^{T} x_{i} w_{2}^{T} x_{i}-2 w_{1}^{T} x_{i} y_{i}-2 w_{2}^{T} x_{i} y_{i}+2 y_{i}^{2}\right) & (3) \end{aligned}$

（1）不等式右边也乘以 2

（2）两项连加运算的和写成两项和的连加运算

（3）将中括号内的平方项展开

（4）移项至与上面保持一致

$\begin{aligned} \operatorname{Cost}\left(w_{1}\right)+\operatorname{Cost}\left(w_{2}\right) &=\sum_{i=1}^{N}\left(w_{1}^{T} x_{i}-y_{i}\right)^{2}+\sum_{i=1}^{N}\left(w_{2}^{T} x_{i}-y_{i}\right)^{2} & (1) \\ &=\sum_{i=1}^{N}\left[\left(w_{1}^{T} x_{i}-y_{i}\right)^{2}+\left(w_{2}^{T} x_{i}-y_{i}\right)^{2}\right] & (2)\\ &=\sum_{i=1}^{N}\left(w_{1}^{T} x_{i} x_{i}^{T} w_{1}-2 w_{1}^{T} x_{i} y_{i}^{2}+y_{i}^{2}+w_{2}^{T} x_{i} x_{i}^{T} w_{2}-2 w_{2}^{T} x_{i} y_{i}+y_{i}^{2}\right) & (3)\\ &=\sum_{i=1}^{N}\left(w_{1}^{T} x_{i} x_{i}^{T} w_{1}+w_{2}^{T} x_{i} x_{i}^{T} w_{2}-2 w_{1}^{T} x_{i} y_{i}-2 w_{2}^{T} x_{i} y_{i}+2 y_{i}^{2}\right) & (4) \end{aligned}$

将不等式右边的项减去不等式左边的项，将得到的差值记为Δ，现在只需要证明这两项的差值大于等于零。

（1）观察两项的结果，最后 3 项是一样的，相减后在化简

（2）将二分之一提出来

（3）观察括号里面会发现为一个平方式

$\begin{aligned} \Delta &=\sum_{i=1}^{N}\left(\frac{w_{1}^{T} x_{i} x_{i}^{T} w_{1}}{2}+\frac{w_{2}^{T} x_{i} x_{i}^{T} w_{2}}{2}-w_{1}^{T} x_{i} w_{2}^{T} x_{i}\right) & (1) \\ &=\frac{1}{2} \sum_{i=1}^{N}\left(w_{1}^{T} x_{i} x_{i}^{T} w_{1}+w_{2}^{T} x_{i} x_{i}^{T} w_{2}-2 w_{1}^{T} x_{i} w_{2}^{T} x_{i}\right) & (2) \\ &=\frac{1}{2} \sum_{i=1}^{N}\left(w_{1}^{T} x_{i}-w_{2}^{T} x_{i}\right)^{2} & (3) \end{aligned}$

不等式右边减去不等式左边的差值为平方式的连加运算，在实数范围内其最后的结果必然是大于等于零，证毕。

线性回归代价函数的解析解

（1）线性回归的代价函数

（2）可以将其改写成两个 N 维向量的点积的形式

（3）使用 A 表示第一个 N 维行向量，则代价函数其实是 A 向量乘以 A 向量的转置

$\begin{aligned} \operatorname{Cost}(w) &=\sum_{i=1}^{N}\left(w^{T} x_{i}-y_{i}\right)^{2} & (1)\\ &=\left(w^{T} x_{1}-y_{1} \ldots w^{T} x_{N}-y_{N}\right)\left(\begin{array}{c} w^{T} x_{1}-y_{1} \\ \cdots \\ w^{T} x_{N}-y_{N} \end{array}\right) & (2)\\ &=\quad A A^{T} & (3) \end{aligned}$

（1）向量 A 的定义

（2）向量 A 可以写成两个 N 维行向量相减

（3）第一个行向量可以将 w 的转置提出来，写成一个 M 维向量和一个 M x N 矩阵相乘（注意 x 原本就是 M 维列向量，N 个列向量 x 就组成了 M x N 矩阵）

（4）定义一个 N x M 矩阵 X ，N 行对应 N 个样本数，M 列对应 M 个变量，定义一个 N 维列向量 y ，N 维对应 N 个样本数。可以看到前面一堆列向量 x 的组合就是矩阵 X 的转置，后面一堆实数 y 的组合就是列向量 y 的转置

$\begin{aligned} A &=\left(w^{T} x_{1}-y_{1} \ldots w^{T} x_{N}-y_{N}\right) & (1)\\ &=\left(w^{T} x_{1} \ldots w^{T} x_{N}\right)-\left(y_{1} \ldots y_{N}\right) & (2) \\ &=w^{T}\left(x_{1} \ldots x_{N}\right)-\left(y_{1} \ldots y_{N}\right) & (3) \\ &=\quad w^{T} X^{T}-y^{T} & (4) \end{aligned}$

（1）代价函数写成两个向量的点积的形式

（2）将上面的简化后的 A 带入到代价函数中

（3）根据转置的性质⁵，可以将后面的转置改写

（4）将乘法展开

（5）观察中间两项发现他们互为转置，又因为这两项最后的结果都是实数，所以这两个项的结果必然是相等的，所以可以合并这两项

$\begin{aligned} \operatorname{Cost}(w) &=A A^{T} & (1)\\ &= \left(w^{T} X^{T}-y^{T}\right)\left(w^{T} X^{T}-y^{T}\right)^{T} & (2) \\ &= \left(w^{T} X^{T}-y^{T}\right)(X w-y) & (3) \\ &= w^{T} X^{T} X w-w^{T} X^{T} y-y^{T} X w+y^{T} y & (4) \\ &= w^{T} X^{T} X w-2 w^{T} X^{T} y+y^{T} y & (5) \end{aligned}$

（1）代价函数对 w 求偏导数，根据向量求导公式，只有第一项和第二项与 w 有关，最后一项为常数，又因为代价函数是个凸函数，当对 W 的偏导数为 0 向量时，代价函数为最小值。

（2）将第二项移项后同时除以2，再两边同时在前面乘以一个逆矩阵，等式左边的矩阵和逆矩阵乘后为单位矩阵，所以只剩下 w 向量。

$\begin{aligned} \frac{\partial \operatorname{Cost}(w)}{\partial w} &=2 X^{T} X w-2 X^{T} y=0 & (1)\\ w &=\left(X^{T} X\right)^{-1} X^{T} y & (2) \end{aligned}$

这样就求出了 w 的解析解，可以看到其中需要求一个矩阵的逆矩阵，当括号内 N x N 矩阵不是满秩矩阵⁶时没有逆矩阵，其本质是自变量 x 之间存在多重共线性⁷（Multicollinearity），下一节将介绍线性回归中的多重共线性与如何解决这种多重共线性的问题。下式中 y 向量前面的部分被称为矩阵 X 的伪逆矩阵，可以通过奇异值分解⁸（SVD）求得矩阵的伪逆矩阵。

$w=\underbrace{\left(X^{T} X\right)^{-1} X^{T}}_{X^{+}} y$

可以看到 X 为 N x M 矩阵，y 为 N 维列向量。在上面工作年限与平均月工资的例子中，X 为一个 5 x 2 的矩阵，y为一个 5 维列向量，最后可以算得 w 为一个 2 维列向量，则这个例子的线性方程为 y = 2172.5 * x - 512.5。

$\begin{bmatrix} 1 & 1\\ 1 & 2\\ 1 & 3\\ 1 & 4\\ 1 & 5 \end{bmatrix} \quad y = \begin{pmatrix} 1598\\ 3898\\ 6220\\ 7799\\ 10510 \end{pmatrix}$

$w=\left(X^{T} X\right)^{-1} X^{T} y=\left(\begin{array}{c} -512.5 \\ 2172.5 \end{array}\right)$

线性回归代价函数解析解的几何解释

线性方程的矩阵形式，其中 y 为 N 维列向量，X 为 N x M 矩阵，w 为 M 维列向量：

$y = X w$

先来看下图，由 M 个黑色 N 维列向量 X 组成了灰色超平面 S，红色 N 维列向量 y 为实际样本点的 y 值，在工作年限与平均月工资的例子中 x1 = (1, 1, 1, 1, 1)、x2 = (1, 2, 3, 4, 5)、y = (1598, 3898, 6220, 7799, 10510)。现在要找由向量 x 经过线性组合后的 y’ ，使得 y - y’ 最小，也就是图中紫色的向量。由图中可以看出，当 y - y’ 是超平面的法向量时最短，也等价于 y - y’ 与每一个向量 x 垂直。

几何解释

（1） y - y’ 与每一个向量 x 垂直，注意结果为 M 维零向量

（2）替换 y’ 的线性方程

（3）展开括号

（4）移项后可解得 w

$\begin{array}{l} X^{T}\left(y-y^{\prime}\right)=0 & (1) \\ X^{T}(y-X w)=0 & (2)\\ X^{T} y-X^{T} X w=0 & (3) \\ w=\left(X^{T} X\right)^{-1} X^{T} y & (4) \end{array}$

可以看到几何解释与通过求导的方式的结果是一致的

五、代码实现

使用 Python 实现线性回归算法：

def linear(X, y):
    """
    线性回归
    args:
        X - 训练数据集
        y - 目标标签值
    return:
        w - 权重系数
    """
   # pinv 函数直接求矩阵的伪逆矩阵
   return np.linalg.pinv(X).dot(y)

六、第三方库实现

scikit-learn⁹ 实现：

from sklearn.linear_model import LinearRegression

# 初始化线性回归器
lin = LinearRegression(fit_intercept=False)
# 拟合线性模型
lin.fit(X, y)
# 权重系数
w = lin.coef_

七、动画演示

当线性方程取不同权重系数时，其对应的代价函数的变化，可以看到代价函数是先减小到一个最小值后然后逐渐增大。

八、思维导图

完整演示请点击这里

注：本文力求准确并通俗易懂，但由于笔者也是初学者，水平有限，如文中存在错误或遗漏之处，恳请读者通过留言的方式批评指正

机器学习算法系列（三）- 标准线性回归算法（Standard Linear Regression Algorithm）