40 | 线性回归(中):如何使用最小二乘法进行直线拟合?

你好,我是黄申。

上一节,我提到了,求解线性回归和普通的线性方程组最大的不同在于误差ε。在求解线性方程组的时候,我们并不考虑误差的存在,因此存在无解的可能。而线性回归允许误差ε的存在,我们要做的就是尽量把ε最小化,并控制在一定范围之内。这样我们就可以求方程的近似解。而这种近似解对于海量的大数据分析来说是非常重要的。

但是现实中的数据一定存在由于各种各样原因所导致的误差,因此即使自变量和因变量之间存在线性关系,也基本上不可能完美符合这种线性关系。总的来说,线性回归分析并不一定需要100%精确,而误差ε的存在可以帮助我们降低对精度的要求。通常,多元线性回归会写作:

$y=b_0+b_1·x_1+b_2·x_2+…+$
$b_{n-1}·x_{n-1}+b_n·x_n+ε$

这里的$x_1,x_2,…,x_n$是自变量,$y$是因变量,$b_0$是截距,$b_1$,$b_2$,…,$b_n$是自变量的系数,$ε$是随机误差。

在线性回归中,为了实现最小化$ε$的目标,我们可以使用最小二乘法进行直线的拟合。最小二乘法通过最小化误差的平方和,来寻找和观测数据匹配的最佳函数。由于这些内容有些抽象,下面我会结合一些例子来解释最小二乘法的核心思想,以及如何使用这种方法进行求解。

使用观测值拟合

在详细阐述最小二乘法之前,我们先来回顾一下第32讲介绍的模型拟合。在监督式学习中,拟合模型其实是指通过模型的假设和训练样本,推导出具体参数的过程。有了这些参数,我们就能对新的数据进行预测。而在线性回归中,我们需要找到观测数据之间的线性关系。

假设我们有两个观测数据,对应于二维空间中的两个点,这两个点可以确定唯一的一条直线,两者呈现线性关系。你可以参考下面这张图。

之后,我们又加入了一个点。这个点不在原来的那条直线上。

这个时候,从线性方程的角度来看,就不存在精确解了。因为没有哪条直线能同时穿过这三个点。这张图片也体现了线性回归分析和求解线性方程组是不一样的,线性回归并不需要求精确解。

如果我们加入更多的观察点,就更是如此了。比如下面这张图。

从上图中你应该可以看出,这根直线不是完全精准地穿过这些点,而只是经过了其中两个,大部分点和这根直线有一定距离。这个时候,线性回归就有用武之地了。

由于我们假设ε的存在,因此在线性回归中,我们允许某条直线只穿过其中少量的点。不过,既然我们允许这种情况发生,那么就存在无穷多这样的直线。比如下面我随便画了几条,都是可以的。

当然,我们从直觉出发,一定不会选取那些远离这些点的直线,而是会选取尽可能靠近这些点的那些线。比如下面这张图里展示的这两条。

好了,既然这样,我们就需要定义哪根线是最优的,以及在给出了最优的定义之后,如何能求解出这条最优的直线呢?最小二乘法可以回答这两个问题,下面我们具体来看。

最小二乘法

最小二乘法的主要思想就是求解未知参数,使得理论值与观测值之差(即误差,或者说残差)的平方和达到最小。我们可以使用下面这个公式来描述。

其中,$y_i$表示来自数据样本的观测值,而$y$^是假设的函数的理论值,$ε$就是我们之前提到的误差,在机器学习中也常被称为损失函数,它是观测值和真实值之差的平方和。最小二乘法里的“二乘”就是指的平方操作。有了这个公式,我们的目标就很清楚了,就是要发现使ε最小化时候的参数。

那么最小二乘法是如何利用最小化$ε$的这个条件来求解的呢?让我们从矩阵的角度出发来理解整个过程。

有了上面的定义之后,我们就可以写出最小二乘问题的矩阵形式。

$min||XB-Y||_{2}^{2}$

其中$B$为系数矩阵,$X$为自变量矩阵,$Y$为因变量矩阵。换句话说,我们要在向量空间中,找到一个$B$,使向量$XB$与$Y$之间欧氏距离的平方数最小的$B$。

结合之前所讲的矩阵点乘知识,我们把上述式子改写为:

$||XB-Y||_{2}^{2}=tr((XB-Y)'(XB-Y))$

其中$(XB-Y)'$表示矩阵$(XB-Y)$的转置。而$tr()$函数表示取对角线上所有元素的和,对于某个矩阵$A$来说,$tr(A)$的值计算如下:

进一步,根据矩阵的运算法则,我们有:

$tr((XB-Y)‘(XB-Y))$
$=tr(B’X’-Y’)(XB-Y)$
$=tr(B’X’XB-B’X’Y-Y’XB+Y’Y)$

因此我们可以得到:

$||XB-Y||_{2}^{2}$
$=tr((XB-Y)‘(XB-Y))$
$=tr(B’X’-Y’)(XB-Y)$
$=tr(B’X’XB-B’X’Y-Y’XB+Y’Y)$

我们知道,求最极值问题直接对应的就是导数为0,因此我对上述的矩阵形式进行求导,得到如下的式子:

$\frac{d||XB-Y||_{2}^{2}}{dB}$
$=\frac{d(tr(B’X’XB-B’X’Y-Y’XB+Y’Y))}{dB}$
$=X’XB+X’XB-X’Y-X’Y$
$=2X’XB-2X’Y$

如果要$||XB-Y||_{2}^{2}$最小,就要满足两个条件。

第一个条件是$\frac{d||XB-Y||_{2}^{2}}{dB}$为0,也就是$2X’XB-2X’Y=0$。

第二个条件是$\frac{d(2X’XB-2X’Y)}{dB}>0$。

由于$\frac{d(2X’XB-2X’Y)}{dB}=2X’X>0$,所以,第二个条件是满足的。只要$2X’XB=2X’Y$。

我们就能获得$ε$的最小值。从这个条件出发,我们就能求出矩阵$B$:

$2X’XB=2X’Y$
$X’XB=X’Y$
$(X’X)^{-1}X’XB=(X’X)^{-1}X’Y$
$IB=(X’X)^{-1}X’Y$
$B=(X’X)^{-1}X’Y$

其中$I$为单位矩阵。而$(X’X)^{-1}$表示$X’X$的逆矩阵。所以,最终系数矩阵为:

$B=(X’X)^{-1}X’Y$

补充证明和解释

为了保持推导的连贯性,在上述的推导过程中,我跳过了几个步骤的证明。下面我会给出详细的解释,供你更深入的学习和研究。

步骤a:

$(XB)‘=B’X’$

证明:

对于$XB$中的每个元素$xb_{i,j}$,有:

而对于$(XB)‘$中的每个元素$xb’_{i,j}$,有:

对于$B’$中的每个元素有:

$b’{i,k}=b{k,i}$

$X’$中的每个元素有:

$x’{k,j}=x{j,k}$

那么,对于$B’X’$中的每个元素$b’x’_{i,j}$,就有:

所以有$(XB)’ = B’X’$。

步骤b:

$(XB-Y)‘=B’X’-Y’$

证明:

和步骤a类似,对于$XB-Y$中的每个元素 $xb-y’_{i,j}$有:

步骤c:

$\frac{d(tr(B’X’Y))}{dB}=X’Y$

证明:

同理,可以证明:

$\frac{d(tr(Y’XB))}{dB}=(Y’X)'=X’Y$

步骤d:

$\frac{d(tr(B’X’XB))}{dB}=2X’XB$

证明:

$\frac{d(tr(B’X’XB))}{dB}$
$=\frac{d(tr(B’(X’XB)))}{dB}+\frac{d(tr((B’X’X)B))}{dB}$
$=(X’XB)+(B’X’X)'$
$=X’XB+X’XB$
$=2X’XB$

步骤e:

常量对于变量求导为0,例如:

$\frac{d(Y’Y)}{dB}=0$

好了,弄明白了这些细节上的证明,你就能更好地理解最小二乘法中的推导步骤。不过,你可能还是会奇怪,为什么最终要对矩阵求导数来求ε的最小值。最后,我们就聊聊如何使用求导获取极小值。

极值是一个函数的极大值或极小值。如果一个函数在一点的某个邻域内每个地方都有确定的值,而该点所对应的值是最大(小)的,那么这函数在该点的值就是一个极大(小)值。而函数的极值可以通过它的一阶和二阶导数来确定。

对于一元可微函数$f(x)$,它在某点$x_0$有极值的充分必要条件是$f(x)$在$x_0$的邻域上一阶可导,在$x_0$处二阶可导,且一阶导数$f’(x_0)=0$,二阶导数$f’‘(x_0)≠0$。其中$f’$和$f’'$分别表示一阶导数和二阶导数。

在一阶导数$f’(x0)=0$的情况下,如果$f’‘(x0)<0$,则$f$在$x_0$取得极大值;如果$f’'(x0)>0$,则$f$在$x_0$取得极小值。这就是为什么在求矩阵$B$的时候,我们要求$2X’XB-2X’Y$为$0$,并且$2X’XB-2X’Y$的导数要大于$0$,这样我们才能确保求得极小值。

总结

今天我们探讨了为什么简单的线性方程组无法满足线性函数拟合的需求,最主要的原因就是现实的观测数据往往不是精确的线性关系,存在一定的误差。我们所要做的就是,在允许一定范围的误差前提下,找到一种线性关系,尽量的满足观察数据,使得我们所定义的误差最小。

最小二乘法通过向量空间的欧氏距离之平方,定义了预测值和真实值之间的误差。在给定自变量和因变量的观测值之后,最小二乘法可以帮助我们推导出所有自变量的系数,并最小化误差。我使用矩阵的形式,为你推导了整个过程。

不过,到目前为止,我们都只是从理论上理解最小二乘法,可能你还没有太深的感触。下一节,我会通过一个具体的例子来逐步进行演算,并使用Python代码对最终的结果进行验证。

思考题

还记得在29讲的线性回归案例吗?我们使用了Boston Housing的数据,拟合出了十多个自变量的系数。请使用这些系数,计算train.csv中所有样本因变量预测值和真实值之间的误差。你可以使用Python代码来实现一下。

欢迎留言和我分享,也欢迎你在留言区写下今天的学习笔记。你可以点击“请朋友读”,把今天的内容分享给你的好友,和他一起精进。

精选留言

  • Peng

    2019-03-18 08:52:16

    老师,线性代数这部分开始很难看懂了,是不是需要先复习一遍线代?请老师指点,明确课程前置条件。
    作者回复

    线性代数基础应该就够了,或者是根据我课程的内容有针对性的补一下,这样效率比较高

    2019-03-19 00:18:08

  • 罗耀龙@坐忘

    2020-05-04 18:33:12

    茶艺师学编程

    刚打开课文,“哇,全是公式,AWSL。”

    但真的学下来,脑子居然蹦出来说:“啊,我能跟上!”(手斜眼看着脑)

    这节课我是这么消化的:

    最小二乘(平方)法,就是在考虑误差的情况下,求解一组方程组的系数矩阵B,并且确保误差最小。按照这样的要求,系数矩阵B=(X’X)^(-1)X’Y。

    这节课就是*花大力气*在推导B是怎么来的。

    在推导过程中,矩阵操作是主线,tr()对角线加和、矩阵点乘的转置/逆。而在确保误差最小时,线性函数的求导操作起大作用,一阶导数为0表示函数为“极值”,二阶导数>0表示该“极值”为“最小值”(反之为“最大值”)。


    感觉回到了大学,但在大学没能现在理解如此清晰···············不愧是黄老师!
    作者回复

    嗯,一步步来,总能理解

    2020-05-05 10:21:41

  • zzz

    2019-04-08 23:10:32

    不太明白为什么要取tr()
    作者回复

    你可以看下L2范式的定义,就能明白为什么我们只关心对角线上的值

    2019-04-09 00:17:15

  • howhigh

    2019-03-21 06:01:04

    黄老师,关于矩阵求导的步骤我依然没有看懂,有没有矩阵求导的资料推荐?
    作者回复

    其实如果你理解了一般函数的求导,矩阵求导并不难理解。你可把矩阵先简化为向量,也就是单个方程式来看,然后整个矩阵就是不同向量的集合。

    2019-03-22 02:28:04

  • 米饭

    2020-08-24 15:02:00

    分享一下国外的网站“数学乐“,里面对于一些数学概念,运算法则讲的很形象具体
    就比如这节课的重点就在如何求ε的最小值,这里用到了导数的概念,可以参考用导数求极大值极小值(https://www.shuxuele.com/calculus/maxima-minima.html)
    对于我这种数学渣也能快速理解
    作者回复

    感谢推荐,我也去看了,不过你确定是国外的网站吗?都是中文的内容,不过节约阅读时间了:)

    2020-08-28 10:36:08

  • Paul Shan

    2019-10-07 09:29:21

    tr(B'X'XB) 对 矩阵 B 的求导,拆成了两部分 tr(B'(X'XB)) 、tr((B'X'X)B) 分别对 B, 这一步没有看懂,主要是因为tr函数,这个函数只有在矩阵中出现而不在其他单因变量的函数中出现,我的疑问是求导的乘法法则可以直接应用在tr函数的里面而不用管tr函数对求导的影响吗?
  • qinggeouye

    2019-03-30 13:10:40

    补充证明的步骤 d ,tr(B'X'XB) 对 矩阵 B 的求导,拆成了两部分 tr(B'(X'XB)) 、tr((B'X'X)B) 分别对 B 求导,求和;意思是 B' 和 B 实际上都与矩阵B相关,第一部分求导相当于固定 B 再对 B 求导,第二部分固定 B‘ 后对 B 求导,而后两部分利用了补充证明步骤 c 的结论得出证明结果?
    作者回复

    是的

    2019-04-01 02:23:58

  • 建强

    2020-10-18 15:59:36

    用python代码计算了一下因变量预测值和真实值之间的误差,误差的平均值和方差都很小,输出结果:
    误差平均值= 8.05648566542149e-17
    误差的方差= 0.2594756827028548

    python程序代码如下:

    # 用线性回归模型预测房价
    import numpy as np
    import pandas as pd
    from sklearn.linear_model import LinearRegression
    from sklearn.preprocessing import StandardScaler
    from sklearn.preprocessing import MinMaxScaler

    df = pd.read_csv("G:/开发实验区/实验数据/房价预测样本数据/boston_housing_data.csv") #读取Boston Housing中的train.csv
    #删除目标值是nan的样本
    nanindex = list(df[df['MEDV'].isna()].index)
    for row in nanindex:
    df.drop(index=row, inplace=True)

    df_features = df.drop(['MEDV'], axis=1) #Dataframe中除了最后一列,其余列都是特征,或者说自变量
    df_targets = df['MEDV'] #Dataframe最后一列是目标变量,或者说因变量

    # 标准化处理后,再进行回归分析
    standardScaler = StandardScaler() #基于Z分数的标准化
    standardScaler.fit(df)
    df_standardized = standardScaler.transform(df) #对原始数据进行标准化,包括特征值和目标变量
    df_features_standardized = df_standardized[:, 0:-1] #获取标准化之后的特征值
    df_targets_standardized = df_standardized[:, -1] #获取标准化之后的目标值

    # 线性回归
    regression_standardized = LinearRegression().fit(df_features_standardized, df_targets_standardized)

    # 计算预测误测

    # 设置系数矩阵B
    B=regression_standardized.coef_.reshape(len(regression_standardized.coef_),1)

    # 设置特征矩阵X
    X = df_features_standardized

    # 设置目标矩阵Y
    Y = df_targets_standardized.reshape(len(df_targets_standardized),1)

    # 计算预测误差绝对值
    E = np.dot(X,B) - Y

    # 输出预测误差的平均值和方差
    print('误差平均值=', np.mean(E))
    print('误差的方差=', np.var(E))
  • 骑行的掌柜J

    2020-07-02 15:47:47

    感觉零基础的我看懂了七七八八。。。但是在求导那里有点懵 还是下来好好补一下线代
    思考题这个“计算 train.csv 中所有样本因变量预测值和真实值之间的误差。”是指的MSE吗?黄老师
    如果是 那么我用traIn_test_split 把训练集和测试集 按默认的测试集占比33%划分数据集后
    最后拟合得到的MSE 是 23.4644515655
    但是按8: 2划分后 MSE 是 23.4148125597
    发现多运行几次 会有不一样的MSE 。。。这是怎么回事呢 黄老师 谢谢
    作者回复

    MSE是指Mean Square Error,常用语线性分析中的误差分析。至于你说的MSE会变化的问题,我认为是和数据划分有关的,具体你可以参考第32讲,有关机器学习中的拟合。即使是同一个标注样本,如果训练和测试的划分不同,也会导致训练集中的数据分布和测试集的不一样,这就会产生新的MSE。除非你能确保每次的训练集和测试集完全一样,MSE才会一样。不过,人们通常都希望尝试不同的训练、测试划分,例如交叉验证,这样能避免误差计算的偶然性

    2020-07-10 00:50:15

  • liying

    2020-04-28 02:22:43

    老师,课程中好像缺步骤b的证明公式
    作者回复

    我的原文是有的,可能是系统显示问题,我让编辑同学看看,谢谢提醒

    2020-05-05 10:38:43

  • 逐风随想

    2019-03-31 00:24:38

    黄老师,越到后面越吃力,怎么办。当年家里条件不好,读到初二就辍学出去工作了。是否需要从新把初中到大学的数学知识补一补呢?
    作者回复

    咱们这个专栏还是讲究实用性,从应用和需求出发。如果觉得有些知识需要补,我的建议是根据不理解知识点出发,逐步倒推,缺啥补啥,一定会有收获

    2019-03-31 02:24:01

  • 那时刻

    2019-03-18 22:51:58

    重温下久违的线性代数
  • 拉欧

    2019-03-18 09:29:17

    今天的数学推导看的有点吃力,不过还好
    作者回复

    加油加油 一步步来

    2019-03-19 00:18:24

  • Geek_62f62c

    2025-02-03 10:03:07

    老师,关于文稿中的证明我有一个疑问,文稿中
    (X'XB)+ (B'X'X)'
    = X'XB+X'XB
    =2X'XB

    我的理解是否应该是
    (X'XB)+ (B'X'X)'
    = X'XB+(X'X)'B
    =X'XB+XX'B

    请问XX‘B和X'XB相等吗
  • 米饭

    2020-08-24 17:29:49

    补充一下步骤c的同理证明
    d(ΣΣ(y'x)_j,i * b_i,j)/db_i,j
    = d(ΣΣ(y'x)'_i,j * b_i,j)/db_i,j
    =(y'x)'_i,j
    =(Y'X)'
    =X'Y
    老师对吗?
    作者回复

    是的

    2020-08-28 10:37:57

  • 猛仔

    2020-04-01 12:08:05

    谢谢老师,之前对最小二乘法的理解一直是记公式,看了老师的讲解终于明白了
    作者回复

    很高兴对你有价值

    2020-04-06 01:52:31

  • Temme

    2019-10-29 13:03:11

    有个疑问,tr()是取对角线,但这里实际tr当中只是一个一行一列的矩阵,感觉这里应该会有拓展,实际应用中有没有可能tr当中是多行多列的矩阵,比如多个自变量对应多个因变量需要去拟合?用最小二乘也应该适用?
  • Paul Shan

    2019-10-03 06:10:55

    最小二乘法是最小化误差作为目标求出系数矩阵。
    虽然以前学过微积分和线性代数,但是矩阵函数的求导还是不太熟悉,老师在这方面有什么文章书籍推荐吗,多谢!
    作者回复

    我觉得矩阵求导,可以细分为单个方程式的求导,如果单个方程式求导理解了,矩阵形式就不难了

    2019-10-05 02:15:25

  • Minis-Shu

    2019-06-19 06:17:04

    老师 我说的非正交矩阵就是行列式不为1的矩阵,我其实是用您讲的最小二乘法求解刚体变换来着,把矩阵看成空间坐标系的一组基向量,矩阵非正交的话,基向量之间的外积就不为0,而利用非正交矩阵进行刚体变换时,向量变换之后,其长度及夹角会发生变化。我就想着,如何将非正交矩阵正交化呢?
  • Minis-Shu

    2019-06-18 16:44:20

    老师,您好!我想问一下,我通过最小二乘求出的系数矩阵是非正交矩阵,有没有什么办法,使得结果是正交阵呢?
    作者回复

    系数矩阵确切的说,应该是个向量,你能否解释一下你这里所说的“非正交矩阵”是指什么?

    2019-06-19 00:43:53