NumPy 中最小二乘算法的高效计算答案

【问题标题】：Efficient computation of the least-squares algorithm in NumPyNumPy 中最小二乘算法的高效计算
【发布时间】：2017-01-15 03:31:18
【问题描述】：

我需要解决大量线性系统，在最小二乘意义上。我无法理解numpy.linalg.lstsq(a, b)、np.dot(np.linalg.pinv(a), b) 和数学实现的计算效率差异。

我使用以下矩阵：

h=np.random.random((50000,100))
a=h[:,:-1].copy()
b=-h[:,-1].copy()

算法的结果是：

# mathematical implementation
%%timeit
np.dot(np.dot(np.linalg.inv(np.dot(a.T,a)),a.T),b)

10 次循环，3 次取胜：36.3 毫秒 每个循环

# numpy.linalg.lstsq implementation
%%timeit
np.linalg.lstsq(a, b)[0]

10 个循环，3 个循环中的最佳：103 毫秒每个循环

%%timeit
np.dot(np.linalg.pinv(a), b)

1 次循环，3 次取胜：216 毫秒 每个循环

为什么会有差异？

【问题讨论】：

为什么它们应该是一样的？
在数学上它们返回相同的结果
@blaz 您是否希望进一步优化它们，或者您对这些方法是否满意并只是想了解时间差异？
我只是想先了解一下:)。

标签： python numpy least-squares

【解决方案1】：

例程lstsq 处理任何系统：过度确定、未确定或完全确定。它的输出是你从 pinv(a)*b 得到的，但它比计算伪逆要快。原因如下：

一般建议：除非确实需要，否则不要计算逆矩阵。求解一个特定右手边的系统比反转它的矩阵要快。

然而，您通过求解 a^Ta = a^Tb 的方法更快，即使您正在反转矩阵。是什么赋予了？问题是，反转 a^Ta 仅在 a 具有完整列排名时才有效。因此，您已将问题限制在这种特殊情况下，并获得了速度，以换取较少的通用性，如下所示，为了降低安全性。

但是对矩阵求逆仍然是低效的。如果您知道 a 具有完整的列排名，则以下操作比您的三个尝试中的任何一个都快：

np.linalg.solve(np.dot(a.T, a), np.dot(a.T, b))

也就是说，在处理条件不佳的矩阵时，lstsq 仍然比上述方法更可取。形成产品 a^Ta 基本上是条件数的平方，因此您更有可能得到无意义的结果。这是一个警示示例，使用 SciPy 的 linalg 模块（本质上等同于 NumPy，但有更多方法）：

import numpy as np
import scipy.linalg as sl
a = sl.hilbert(10)    # a poorly conditioned square matrix of size 10
b = np.arange(10)     # right hand side
sol1 = sl.solve(a, b)
sol2 = sl.lstsq(a, b)[0]
sol3 = sl.solve(np.dot(a.T, a), np.dot(a.T, b))

这里lstsq 提供与solve 几乎相同的输出（该系统的独特解决方案）。然而，sol3 是完全错误的，因为数字问题（你甚至不会被警告）。

sol1：

  [ -9.89821788e+02,   9.70047434e+04,  -2.30439738e+06,
     2.30601241e+07,  -1.19805858e+08,   3.55637424e+08,
    -6.25523002e+08,   6.44058066e+08,  -3.58346765e+08,
     8.31333426e+07]

sol2：

  [ -9.89864366e+02,   9.70082635e+04,  -2.30446978e+06,
     2.30607638e+07,  -1.19808838e+08,   3.55645452e+08,
    -6.25535946e+08,   6.44070387e+08,  -3.58353147e+08,
     8.31347297e+07]

sol3：

  [  1.06913852e+03,  -4.61691763e+04,   4.83968833e+05,
    -2.08929571e+06,   4.55280530e+06,  -5.88433943e+06,
     5.92025910e+06,  -5.56507455e+06,   3.62262620e+06,
    -9.94523917e+05]

【讨论】：