Question

from sklearn import datasets, linear_model, metrics
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures
import math, scipy, numpy as np
from scipy import linalg

np.set_printoptions(precision=6)

data = datasets.load_diabetes()

feature_names=['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']

trn,test,y_trn,y_test = train_test_split(data.data, data.target, test_size=0.2)

trn.shape, test.shape

# ((353, 10), (89, 10))

def regr_metrics(act, pred):
    return (math.sqrt(metrics.mean_squared_error(act, pred)), 
    metrics.mean_absolute_error(act, pred))

sklearn 如何实现它

sklearn 是如何做到这一点的？ 通过检查 源代码 ，你可以看到在密集的情况下，它调用 scipy.linalg.lstqr ，它调用 LAPACK 方法：

选项是 'gelsd' ， 'gelsy' ， 'gelss' 。默认值 'gelsd' 是个好的选择，但是， 'gelsy' 在许多问题上更快一些。 'gelss' 由于历史原因而使用。它通常更慢但是使用更少内存。

• gelsd ：使用 SVD 和分治方法
• gelsy ：使用 QR 分解
• gelss ：使用 SVD

Scipy 稀疏最小二乘

我们不会详细介绍稀疏版本的最小二乘法。如果你有兴趣，请参考以下信息：

Scipy 稀疏最小二乘使用称为 Golub 和 Kahan 双对角化 的迭代方法。

Scipy 稀疏最小二乘源代码：预处理是减少迭代次数的另一种方法。如果有可能有效地求解相关系统 M*x = b ，其中 M 以某种有用的方式近似 A （例如， M-A 具有低秩或其元素相对于 A 的元素较小），则 LSQR 可以在系统 A*M(inverse)*z = b 更快地收敛。之后可以通过求解 M * x = z 来恢复 x 。

如果 A 是对称的，则不应使用 LSQR！替代方案是对称共轭梯度法（cg）和/或 SYMMLQ。 SYMMLQ 是对称 cg 的一种实现，适用于任何对称的 A ，并且比 LSQR 更快收敛。如果 A 是正定的，则存在对称 cg 的其他实现，每次迭代需要的工作量比 SYMMLQ 少一些（但需要相同的迭代次数）。

linalg.lstqr

sklearn 实现为我们添加一个常数项（因为对于我们正在学习的直线， y 截距可能不是 0）。 我们现在需要手工完成：

trn_int = np.c_[trn, np.ones(trn.shape[0])]
test_int = np.c_[test, np.ones(test.shape[0])]

由于 linalg.lstsq 允许我们指定，我们想要使用哪个 LAPACK 例程，让我们尝试它们并进行一些时间比较：

%timeit coef, _,_,_ = linalg.lstsq(trn_int, y_trn, lapack_driver="gelsd")

# 290 µs ± 9.24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit coef, _,_,_ = linalg.lstsq(trn_int, y_trn, lapack_driver="gelsy")

# 140 µs ± 91.7 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit coef, _,_,_ = linalg.lstsq(trn_int, y_trn, lapack_driver="gelss")

# 199 µs ± 228 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)