sklearn的LinearRegression源码理解

Posted 2022-02-13 ybdesire

0. 引入

下面是一个简单的数据集：

x_data = np.array(  [[1,1,1 ],[2,4,8],[3,9,27],[4,16,64]]  )
y_data = np.array( [3,2,0,5] )

用sklearn自带的的LinearRegression，得到的结果（这是正确结果）。

model.coef_ = array([15.16666667, -8.5       ,  1.33333333])
model.intercept_ = −5.00000000000028

也就是说，最终拟合的直线为：y=-5+15.17x_1-8.5x_2+1.33*x_3

但是，用numpy从头实现的LinearRegression（参考2中有详细代码），得到的结果，就与这个大有不同了。

为什么会有不同呢？

1. numpy从头实现的LinearRegression

参考2中有详细代码，下面给出原理注释：

import numpy as np
import pandas as pd
    
class XLinearRegression() :
    def __init__( self, learning_rate, iterations ):
        self.learning_rate = learning_rate# 学习速率
        self.iterations = iterations# 迭代次数
        
    # 训练模型f(x)=wx+b
    def fit( self, X, Y ):
        # self.m是样本数量，self.n是特征数量
        self.m, self.n = X.shape
        # 初始化权重
        self.W = np.zeros( self.n )
        self.b = 0
        self.X = X
        self.Y = Y
          
        # 梯度下降更新权重的过程
        for i in range( self.iterations ):
            self.update_weights()
            
        return self
      
    # 梯度下降法更新权重
    def update_weights( self ):
        Y_pred = self.predict( self.X )#根据上一次求解得到的权重，获取预测值
        # 计算梯度
        #   MSE=Sigma(true-pred)^2， MSE对W和b分别求偏导数即可得到如下梯度
        dW = - ( 2 * ( self.X.T ).dot( self.Y - Y_pred )  ) / self.m
        db = - 2 * np.sum( self.Y - Y_pred ) / self.m 
        # 更新权重（负梯度方向）
        self.W = self.W - self.learning_rate * dW
        self.b = self.b - self.learning_rate * db
        return self
      
    # 根据求得的权重，求待拟合函数f(x)的值
    def predict( self, X ):
        return X.dot( self.W ) + self.b

下面代码使用该线性回归

x_data = np.array(  [[1,1,1 ],[2,4,8],[3,9,27],[4,16,64]]  )
y_data = np.array( [3,2,0,5] )
model = XLinearRegression(learning_rate=0.001, iterations=30)
model.fit(x_data,y_data)
print(model.W,model.b)

得到的权重W为：[ -11305.592508 -41387.85209775 -155378.28247138]
偏移（截距）b为：-3211.3991546315137

结果错的很离谱。

2. sklearn的LinearRegression源码原理

那sklearn的LinearRegression是如何求解权重的呢？我们参考1对sklearn的LinearRegression做个动态调试：
详细步骤见参考1。

1. 找到第三方库所在的位置
   先利用如下Python代码找到sklearn源码位置。我的位置在 /xx/anaconda3/envs/env_test_py38/lib/python3.8/site-packages/sklearn

import sklearn, os
path = os.path.dirname(sklearn.__file__)

2. 在第三方库源码中加断点

首先找到sklearn的LinearRegression的源码，见参考3。
然后在这个位置打开文件，用pdb加上断点

• /xx/anaconda3/envs/env_test_py38/lib/python3.8/site-packages/sklearn/linear_model/_base.py

3. 运行测试程序

测试程序如下，直接调用 sklearn的LinearRegression 来训练模型。

from sklearn.linear_model import LinearRegression
import numpy as np

x_data = np.array(  [[1,1,1],[2,4,8],[3,9,27],[4,16,64]]  )
y_data = np.array( [3,2,0,5] )

model = LinearRegression()
model.fit(x_data,y_data)

4. 源码调试过程

加上断点后，可以单步调试。

从调试中，可以看到，对于如本文这样一般性的数据，使用默认的fit()参数，源码主要做了两件事情

（1）数据预处理，通过如下函数，对数据做了normalization

X, y, X_offset, y_offset, X_scale = self._preprocess_data(
    X,
    y,
    fit_intercept=self.fit_intercept,
    normalize=_normalize,
    copy=self.copy_X,
    sample_weight=sample_weight,
    return_mean=True,
)

在预处理之前，数据值为

(Pdb) X,y
(array([[ 1,  1,  1],
       [ 2,  4,  8],
       [ 3,  9, 27],
       [ 4, 16, 64]]), array([3, 2, 0, 5]))

在预处理之后，数据值变为

(Pdb) X,y
(array([[ -1.5,  -6.5, -24. ],
       [ -0.5,  -3.5, -17. ],
       [  0.5,   1.5,   2. ],
       [  1.5,   8.5,  39. ]]), array([ 0.5, -0.5, -2.5,  2.5]))

注意这里连y值也变了

（2）调用numpy中的linalg.lstsq最小二乘法来求解线性方程

elf.coef_, self._residues, self.rank_, self.singular_ = linalg.lstsq(X, y)

这里是对预处理之后的数据进行求解。

3. sklearn的LinearRegression的数据预处理原理

从源码中，可以找到预处理部分的源码（详见参考4），加上断点调试后，可以发现，最关键的预处理过程源码为：

X_offset = np.average(X, axis=0)
X -= X_offset
y_offset = np.average(y, axis=0)
y = y - y_offset

即对矩阵X按列求平均值，并用每一列的数值减去平均值；对y也是一样，求平均值后相减。

把这个过程用如下代码测试：

import numpy as np
x_data = np.array(  [[1,1,1],[2,4,8],[3,9,27],[4,16,64]]  )
y_data = np.array( [3,2,0,5] )
x_data = x_data-np.average(x_data,axis=0)
y_data = y_data-np.average(y_data)
print(x_data)
print(y_data)

得到的结果为

[[ -1.5  -6.5 -24. ]
 [ -0.5  -3.5 -17. ]
 [  0.5   1.5   2. ]
 [  1.5   8.5  39. ]]
[ 0.5 -0.5 -2.5  2.5]

这个结果与2.4中“在预处理之后，数据值变为”的结果完全一致。

4. 用 numpy.linalg.lstsq 复现sklearn的LinearRegression结果

numpy.linalg.lstsq只能计算y=Wx+b中的W，无法得到b的值。

对sklearn的LinearRegression中计算b（intercept_）的过程进一步动态调试（源码详见参考5），可以发现，b的计算过程为：

-> self.coef_ = self.coef_ / X_scale
(Pdb) X_scale
array([1., 1., 1.])
(Pdb) self.coef_
array([15.16666667, -8.5       ,  1.33333333])
(Pdb) n
-> self.intercept_ = y_offset - np.dot(X_offset, self.coef_.T)
(Pdb) X_offset
array([ 2.5,  7.5, 25. ])
(Pdb) y_offset
2.5
(Pdb) n
--Return--
-> self.intercept_ = y_offset - np.dot(X_offset, self.coef_.T)
(Pdb) self.intercept_
-5.000000000000277

这样得到的截距值-5.00就是正确的了。

下面给出完整版本的使用numpy.linalg.lstsq 复现sklearn的LinearRegression结果的代码

import numpy as np

# setp-01: original dataset
x_data = np.array(  [[1,1,1],[2,4,8],[3,9,27],[4,16,64]]  )
y_data = np.array( [3,2,0,5] )
# step-02: preprocess
x_offset = np.average(x_data,axis=0)
x_data = x_data-x_offset
y_offset = np.average(y_data)
y_data = y_data-y_offset
# step-03: lstsq
coef,_,_,_ = np.linalg.lstsq(x_data,y_data,rcond=None)
x_scale = np.ones(len(coef))
intercept = y_offset - np.dot(x_offset, coef.T)
# result
print(coef)# W = [15.16666667 -8.5         1.33333333]
print(intercept)# b = -5.000000000000277

5. 总结

1. sklearn的LinearRegression源码，会首先对数据做预处理（X,Y都减去列均值），然后再用numpy中lstsq求解线性回归
2. 如果不做这个预处理，得到的值是与直接调用sklearn的LinearRegression得到的不相同的
3. 动态调试结束，一定要删除断点
4. 目前还没有找到靠谱的如原文1中的XLinearRegression一样从头写的多元线性回归源代码

6. 参考

1. 动态调试python第三方库。https://blog.csdn.net/ybdesire/article/details/54649211
2. https://www.geeksforgeeks.org/linear-regression-implementation-from-scratch-using-python/
3. sklearn的LinearRegression源码. https://github.com/scikit-learn/scikit-learn/blob/0d378913b/sklearn/linear_model/_base.py#L507
4. https://github.com/scikit-learn/scikit-learn/blob/0d378913be6d7e485b792ea36e9268be31ed52d0/sklearn/linear_model/_base.py#L213
5. https://github.com/scikit-learn/scikit-learn/blob/0d378913be6d7e485b792ea36e9268be31ed52d0/sklearn/linear_model/_base.py#L366