sklearn中LinearRegression关键源码解读

问题的引入

我们知道，线性回归方程的参数，可以用梯度下降法求解，或者用正规方程求解。
那sklearn.linear_model.LinearRegression中，是不是可以指定求解方式呢？能不能从中获取梯度相关信息呢？
下面是线性回归最简单的用法。

from sklearn import  linear_model

# Create linear regression object
regr = linear_model.LinearRegression()

# Train the model using the training sets
regr.fit(x_train, y_train)

# Predict
regr.predict(x_test)

参考LinearRegression文档，发现并不能指定用梯度下降法求解，或者用正规方程求解该线性回归，更不能获取该线性回归梯度相关的信息。

这是为什么呢？sklearn.linear_model.LinearRegression是用什么方法求解线性回归方程（一般用Y=WX+b来表示）参数的呢？

为了回答这些问题，我们需要看一下sklearn.linear_model.LinearRegression的代码，来了解它的实现方式。

sklearn.linear_model.LinearRegression源码解读

在github可以找到sklearn.linear_model.LinearRegression的实现代码。通过寻找LinearRegression类的fit()方法，我们定位到线性回归方程训练阶段的代码。

if sp.issparse(X):#如果是稀疏矩阵
    if y.ndim < 2:
        out = sparse_lsqr(X, y)
        self.coef_ = out[0]
        self._residues = out[3]
    else:
        # sparse_lstsq cannot handle y with shape (M, K)
        outs = Parallel(n_jobs=n_jobs_)(
            delayed(sparse_lsqr)(X, y[:, j].ravel())
            for j in range(y.shape[1]))
        self.coef_ = np.vstack(out[0] for out in outs)
        self._residues = np.vstack(out[3] for out in outs)
else:
    self.coef_, self._residues, self.rank_, self.singular_ = \
        linalg.lstsq(X, y)
    self.coef_ = self.coef_.T

其中，线性回归方程的参数获取，主要靠这两个方法：sparse_lsqr()或linalg.lstsq()。若训练集X是稀疏矩阵，则用sparse_lsqr()，否则用linalg.lstsq()。

这两个方法是怎么求解线性回归方程参数的呢？

linalg.lstsq()方法源码解读

从LinearRegression类的源码注释中，可以发现如下说明。

"""
Notes
-----
From the implementation point of view, this is just plain Ordinary
Least Squares (scipy.linalg.lstsq) wrapped as a predictor object.
"""

也就是说，这里是用scipy.linalg.lstsq()方法来计算线性回归方程参数的，sklearn.linear_model.LinearRegression是用普通的最小二乘（OLS）法求解线性回归方程参数的。
我们跟进去看看scipy.linalg.lstsq()是怎么实现的。

scipy.linalg.lstsq()方法源码解读

从github中找到scipy.linalg.lstsq()方法的代码实现。从它的源码注释中，我们发现，可以通过参数lapack_driver来选择求解线性回归方程参数的方法。

"""
lapack_driver: str, optional
    Which LAPACK driver is used to solve the least-squares problem.
    Options are ``'gelsd'``, ``'gelsy'``, ``'gelss'``. Default
    (``'gelsd'``) is a good choice.  However, ``'gelsy'`` can be slightly
    faster on many problems.  ``'gelss'`` was used historically.  It is
    generally slow but uses less memory.
"""

可见，这里一共提供了三种方法来求解least-squares problem（最小均方问题，这也是线性回归模型的优化目标）：gelsd, gelsy, gelss，它们是什么意思呢，看了这页代码，也看了文档，都没找到解释！不过能看到，这个参数是传给函数get_lapack_funcs()的（如下），那就跟进去看看函数get_lapack_funcs()的代码吧。

driver = lapack_driver
if driver is None:
    driver = lstsq.default_lapack_driver
if driver not in ('gelsd', 'gelsy', 'gelss'):
    raise ValueError('LAPACK driver "%s" is not found' % driver)

lapack_func, lapack_lwork = get_lapack_funcs((driver,
                                             '%s_lwork' % driver),
                                             (a1, b1))

get_lapack_funcs()方法源码解读

找到get_lapack_funcs()方法源码，发现它是调用_get_funcs()实现的（如下）。

"""
In LAPACK, the naming convention is that all functions start with a
type prefix, which depends on the type of the principal
matrix. These can be one of {'s', 'd', 'c', 'z'} for the numpy
types {float32, float64, complex64, complex128} respectevely, and
are stored in attribute `typecode` of the returned functions.
"""
return _get_funcs(names, arrays, dtype,
                  "LAPACK", _flapack, _clapack,
                  "flapack", "clapack", _lapack_alias)

原来这些gelsd, gelsy, gelss，是C函数名，函数是从LAPACK—Linear Algebra PACKage中拿到的。

gelsd：它是用singular value decomposition of A and a divide and conquer method方法来求解线性回归方程参数的。
gelsy：computes the minimum-norm solution to a real/complex linear least squares problem
gelss：Computes the minimum-norm solution to a linear least squares problem using the singular value decomposition of A.

scipy.linalg.lstsq()方法默认选择gelsd。由于sklearn.linear_model.LinearRegression并没有提供参数选择所用解法，所以我们也只能用默认的gelsd。

sparse_lsqr()方法源码解读

当训练集X是稀疏矩阵时，用sparse_lsqr()方法求解线性回归方程的参数。
sparse_lsqr()方法是在这里被导入的。

if sp_version < (0, 15):
    # Backport fix for scikit-learn/scikit-learn#2986 / scipy/scipy#4142
    from ._scipy_sparse_lsqr_backport import lsqr as sparse_lsqr
else:
    from scipy.sparse.linalg import lsqr as sparse_lsqr

原来sparse_lsqr()就是不同版本的scipy.sparse.linalg.lsqr()。它在源码注释中指出，是参考论文C. C. Paige and M. A. Saunders (1982a). "LSQR: An algorithm for sparse linear equations and sparse least squares", ACM TOMS 8(1), 43-71.中的Golub & Kahan双对角线化过程实现的。

结论

sklearn.linear_model.LinearRegression求解线性回归方程参数时，首先判断训练集X是不是稀疏矩阵，如是，就用Golub & Kahan双对角线化过程方法来求解；否则就调用C库LAPACK中的用基于分治法的奇异值分解来求解。

这些解法都跟梯度下降没有半毛钱的关系，可见常见的求解线性回归方程参数理论与实际用代码实现的数值解法还是有很大区别。

读源码时，要特别注意看注释，否则有些数学解法没有解读很难看懂。