机器算法验证 - 为什么不使用“正规方程”来找到简单的最小二乘系数？ - 吾爱随笔录

为什么不使用“正规方程”来找到简单的最小二乘系数？

机器算法验证回归最小二乘 scikit-学习

2022-02-03 14:30:50

我在这里看到了这个列表，不敢相信有这么多方法可以解决最小二乘问题。维基百科上的“正规方程”似乎是一种相当直接的方式：

\begin{aligned} \hat{α} & = \bar{y} - \hat{β} \bar{x}, \\ \hat{β} & = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x}) (y_{i} - \bar{y})}{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2}} \end{aligned}

${\begin{aligned}{\hat {\alpha }}&={\bar {y}}-{\hat {\beta }}\,{\bar {x}},\\{\hat {\beta }}&={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}\end{aligned}}$

那么为什么不直接使用它们呢？我认为一定存在计算或精度问题，因为在上面的第一个链接中，Mark L. Stone 提到 SVD 或 QR 是统计软件中的流行方法，并且正规方程“从可靠性和数值精度的角度来看是可怕的”。然而，在下面的代码中，与三个流行的 python 函数相比，正规方程的精度可以精确到小数点后 12 位：numpy 的polyfit；scipy 的线性回归；和 scikit-learn 的LinearRegression。

更有趣的是，当 n = 100000000 时，正规方程方法是最快的。我的计算时间是：2.5s for linregress；polyfit 12.9s；线性回归 4.2s；正常方程为 1.8s。

代码：

import numpy as np
from sklearn.linear_model import LinearRegression
from scipy.stats import linregress
import timeit

b0 = 0
b1 = 1
n = 100000000
x = np.linspace(-5, 5, n)
np.random.seed(42)
e = np.random.randn(n)
y = b0 + b1*x + e

# scipy                                                                                                                                     
start = timeit.default_timer()
print(str.format('{0:.30f}', linregress(x, y)[0]))
stop = timeit.default_timer()
print(stop - start)

# numpy                                                                                                                                      
start = timeit.default_timer()
print(str.format('{0:.30f}', np.polyfit(x, y, 1)[0]))
stop = timeit.default_timer()
print(stop - start)

# sklearn                                                                                                                                    
clf = LinearRegression()
start = timeit.default_timer()
clf.fit(x.reshape(-1, 1), y.reshape(-1, 1))
stop = timeit.default_timer()
print(str.format('{0:.30f}', clf.coef_[0, 0]))
print(stop - start)

# normal equation                                                                                                                            
start = timeit.default_timer()
slope = np.sum((x-x.mean())*(y-y.mean()))/np.sum((x-x.mean())**2)
stop = timeit.default_timer()
print(str.format('{0:.30f}', slope))
print(stop - start)

3个回答

对于问题，形成正规方程通过形成的条件数。粗略地说，是如果一切顺利，您在计算中丢失的位数。的倒数并没有任何关系。无论是如何解决的，你已经失去了位的准确性。即，形成正态方程使准确度损失的位数增加了一倍，马上就可以了。 $Ax \approx b$ $A$ $A^TA$ $log_{10}(cond)$ $A^TA$ $A^TAx = A^Tb$ $log_{10}(cond(A^TA)) = 2*log_{10}(cond(A))$

如果条件数很小（一个是最好的），则无关紧要。如果条件数 =并且您使用 QR 或 SVD 等稳定方法，则双精度可能有 8 位左右的精度。如果您形成正态方程，则您已将条件数平方为，并且您的答案基本上没有准确性。 $10^8$ $10^{16}$

有时您可以摆脱正常方程，有时则不会。

如果您只需要解决这一变量问题，请继续使用公式。它没有任何问题。例如，我可以看到您在 ASM 中为嵌入式设备编写了几行代码。事实上，我在某些情况下使用了这种解决方案。当然，您不需要为了解决这个小问题而拖拽大型统计库。

数值不稳定性和性能是较大问题和一般设置的问题。如果您解决多元最小二乘等问题。对于一般问题，您当然不会使用它。

没有现代统计软件包可以用正态方程求解线性回归。正规方程只存在于统计书籍中。

不应该使用正规方程，因为计算矩阵的逆是非常有问题的。

当有封闭式数学解决方案可用时，为什么要使用梯度下降进行线性回归？

...尽管可以使用直接法线方程。请注意，在正规方程中，必须反转矩阵。现在反转矩阵的计算成本为 O(N3)，其中 N 是 X 矩阵中的行数，即观察值。此外，如果 X 是病态的，那么它将在估计中产生计算错误......

其它你可能感兴趣的问题

上一篇在贝叶斯统计中，数据被认为是非随机的，但可以有概率或有条件。如何？下一篇为什么我们不使用有效数字？