使用Numpy的最小二乘法不会产生残差

Question

我正在尝试使用Numpy来计算最小二乘问题（即简单回归的普通最小二乘法），以找到相应的R²值。但是，在某些情况下，Numpy会返回空列表作为残差。以下是一个过度确定的示例（即方程数大于未知数）, 说明了这个问题:

（注：没有常数因子（即拦截器）（即所有1的初始列向量），因此将使用未居中的总平方和（TSS）。）

import numpy as np

A = np.array([[6, 6, 3], [40, 40, 20]]).T
y = np.array([0.5, 0.2, 0.6])

model_parameters, residuals, rank, singular_values = np.linalg.lstsq(A, y, rcond=None)

# No Intercept, therefore use Uncentered Total Sum of Squares (TSS)
uncentered_tss = np.sum((y)**2)  
numpy_r2 = 1.0 - residuals / uncentered_tss

print("Numpy Model Parameter(s): " + str(model_parameters))
print("Numpy Sum of Squared Residuals (SSR): " + str(residuals))
print("Numpy R²: " + str(numpy_r2))

以下代码会产生以下输出：

Numpy Model Parameter(s): [0.00162999 0.01086661]
Numpy Sum of Squared Residuals (SSR): []
Numpy R²: []

根据numpy文档：

当方程是欠定的或者是完全确定的时，残差为空；当方程是超定的时，残差会有返回值。 然而，这个问题显然是超定的（3个方程和2个未知数）。我甚至可以通过计算statsmodels的OLS函数给出的回归结果来证明残差（从而证明平方残差和（SSR））是存在的。

import statsmodels.api as sm

A = np.array([[6, 6, 3], [40, 40, 20]]).T
y = np.array([0.5, 0.2, 0.6])

statsmodel_model = sm.OLS(y, A)
regression_results = statsmodels_model.fit()

calculated_r_squared = 1.0 - regression_results.ssr / np.sum((y)**2)

print("Parameters: " + str(regression_results.params))
print("Residuals: " + str(regression_results.resid))
print("Statsmodels R²: " + str(regression_results.rsquared))
print("Manually Calculated R²: " + str(calculated_r_squared))

以下代码会产生以下输出：

Parameters: [0.00162999 0.01086661]
Residuals: [ 0.05555556 -0.24444444  0.37777778]
Statsmodels R²: 0.6837606837606838
Manually Calculated R²: 0.6837606837606838

(正如您所看到的，Statsmodels和Numpy模型具有相同的参数。) 为什么使用以下示例时Numpy会返回一个空的SSR数组？这是numpy.linalg.lstsq的错误吗？ 如果这不是一个错误，那么为什么Statsmodels能够计算平方残差和（SSR）而numpy不能呢？我们也可以根据最佳拟合平面手动计算残差：

- Code Doggo

Answer 1

从numpy.linalg.lstsq()文档中：

残差：{（），（1），（K，）} ndarray

……如果a的秩< N或M <= N，则这是一个空数组。…

您的矩阵的秩为1。

---

注意：您认为“缺失”的残差也可以使用 numpy 找到（无需其他软件包）：

residuals = y - np.dot(A, model_parameters)

- AGN Gazer

---

关于您更新的答案。我仍然不明白为什么Numpy选择不在最小二乘调用中给出残差。除非您还有其他需要纠正的地方（如果可以，请告诉我），否则我将提交一个新问题以了解背后的原因。 - Code Doggo

@DanHoynoski 请参考http://web.gps.caltech.edu/classes/ge193.old/lectures/Lecture2.pdf，第7页。这就是你所拥有的。想一想，你有同一个方程（“数量”）等于许多值：`3*b0+20*b1=0.25; 3*b0+20*b1=0.1; 3*b0+20*b1=0.6。“许多值”（0.25、0.1、0.6）并不意味着它是一个超定系统。使你的系统“欠定”的原因是你无法解决“两个”未知数（b0和b1`），因为你所有的方程都是线性相关的。因此，你的解决方案是不确定的。Numpy只是给出了一个可能的解决方案，而这个解决方案有无限多种可能。 - AGN Gazer

“Numpy返回了所有两个参数。因此，你关于“只有一个未知数可以被确定”的说法是不正确的。” Numpy返回具有最小范数的“a”解。这并不意味着你的解决方案已经确定。 - AGN Gazer

"...而不意识到它是垃圾。" - 为什么许多解决方案中的一个会是垃圾？因为Numpy选择了最小化范数的解决方案，所以这是在当前数据下最好的解决方案。 - Code Doggo

@DanHoynoski 一个最小范数的解并不意味着它是最小化残差的解。不同的情况需要不同的解决方案，你认为对于你特定问题的好解决方案可能不是其他问题的最佳解决方案（而且很可能不是，因为有无限多个解决方案）。对我来说，拥有一个最小范数的解没有任何意义（一般而言）：我想要“最佳”解决方案——最小化残差的解决方案，但在欠定系统中我无法得到它。 - AGN Gazer