求解线性最小二乘法的最快方法

Question

在https://math.stackexchange.com/a/2233298/340174中提到，如果通过 LU 分解完成，求解线性方程“M\xc2\xb7x = b”（矩阵 M 是方）会很慢（但甚至更慢）使用 QR 分解）。现在我注意到这numpy.linalg.solve实际上是使用 LU 分解。事实上，我想求解“V\xc2\xb7x = b”以获得最小二乘的非平方范德蒙设计矩阵 V。我不想要正则化。我看到多种方法：

\n\n

\n
1. 用 求解“V\xc2\xb7x = b” numpy.linalg.lstsq，它使用基于 SVD 的 Fortran“xGELSD”。SVD 应该比 LU 分解更慢，但我不需要计算“(V^T\xc2\xb7V)”。
\n
2. 使用 求解“(V^T\xc2\xb7V)\xc2\xb7x = (V^T\xc2\xb7b)” numpy.linalg.solve，它使用 LU 分解。
\n
3. 用 求解“A\xc2\xb7x = b” numpy.linalg.solve，它使用LU分解，但直接根据https://math.stackexchange.com/a/3155891/340174计算“A=xV^T\xc2\xb7V”
\n

\n\n

或者，我可以使用 scipy 的最新版本solve（https://docs.scipy.org/doc/scipy-1.2.1/reference/ generated/scipy.linalg.solve.html），它可以对对称矩阵“A”使用对角旋转“（我猜这比使用 LU 分解更快），但我的 scipy 停留在 1.1.0 上，所以我无法访问它。

\n\n

从/sf/answers/3187486641/看来，它solve比 更快lstsq，包括计算“V^T\xc2\xb7V”，但当我尝试时，lstsq速度更快。也许我做错了什么？

\n\n

解决线性问题的最快方法是什么？

\n\n

---

\n\n

没有实际选择

\n\n

\n
• statsmodels.regression.linear_model.OLS.fit是使用 Moore-Penrose 伪逆或 QR 分解 + np.linalg.inv+ np.linalg.svd+ numpy.linalg.solve，这对我来说似乎不太有效。
\n
• sklearn.linear_model.LinearRegression使用 scipy.linalg.lstsq。
\n
• scipy.linalg.lstsq还使用 xGELSD。
\n
• 我预计计算“(V^T\xc2\xb7V)”的倒数非常昂贵，因此我放弃了“x = (V^T\xc2\xb7V)^-1\xc2\xb7(V ^T\xc2\xb7b)"
\n

\n

Answer 1

尝试scipy.linalg.lstsq()使用lapack_driver=\'gelsy\'！

\n\n

让我们回顾一下求解线性最小二乘的不同例程和方法：

\n\n

\n

• numpy.linalg.lstsq()包装 LAPACK's xGELSD()，如umath_linalg.c.src第 2841+ 行所示。该例程使用分治策略将矩阵 V 简化为双对角形式，并计算该双对角矩阵的 SVD。

\n

• scipyscipy.linalg.lstsq()包装了 LAPACKxGELSD()和：可以修改参数以从一个切换到另一个xGELSY()。根据 LAPACK 的基准测试，比 慢，比 快 5 左右。利用 V 的 QR 分解和列旋转。好消息是这个开关已经在 scipy 1.1.0 中可用！xGELSS()lapack_driverxGELSS()xGELSD()xGELSY()xGELSD()xGELSY()

\n
• LAPACKxGELS()使用矩阵 V 的 QR 分解，但它假设该矩阵具有满秩。根据 LAPACK 的基准， on 预计dgels()会比 快大约 5 倍dgelsd()，但它也更容易受到矩阵条件数的影响，并且可能变得不准确。请参阅C++（LAPACK、sgels）和 Python（Numpy、lstsq）结果之间的差异中的详细信息和进一步参考。xGELS() 在scipy 的 cython-lapack 接口中可用。
\n

\n\n

虽然非常诱人，但计算和使用V^T\xc2\xb7V来求解正规方程可能不是正确的方法。事实上，精度受到该矩阵的条件数的威胁，大约是矩阵 V 的条件数的平方。由于范德蒙德矩阵往往是严重病态的，除了离散傅里叶变换的矩阵之外，它可能会变得危险......最后，您甚至可以继续使用xGELSD()以避免与条件相关的问题。如果您切换到xGELSY()，请考虑估计错误。

\n

Answer 2

我将忽略问题的范德蒙德部分（bubble的评论指出它有一个解析逆）并回答有关其他矩阵的更一般的问题。

我认为这里可能会混淆一些事情，所以我将区分以下几点：

1. V x = b使用LU的精确解
2. V x = b使用QR的精确解
3. V x = b使用QR的最小二乘解
4. V x = b使用SVD的最小二乘解
5. V^T V x = V^T b使用LU的精确解
6. V^T V x = V^T b使用 Cholesky的精确解

您链接到的第一个 maths.stackexchange 答案是关于情况 1 和 2 的。当它说 LU 很慢时，它意味着相对于特定类型矩阵的方法，例如正定矩阵、三角形矩阵、带状矩阵……

但我认为你实际上问的是 3-6。最后一个 stackoverflow 链接指出 6 比 4 快。正如您所说，4 应该比 3 慢，但 4 是唯一适用于rank-deficient 的情况V。一般来说，6 应该比 5 快。

我们应该确保您执行的是 6 而不是 5。要使用 6，您需要使用scipy.linalg.solvewith assume_a="pos"。否则，你最终会做 5 件事。

我还没有找到一个可以执行 3 in numpyor 的函数scipy。Lapack例程是xGELS，它似乎没有在scipy. 您应该能够scupy.linalg.qr_multiply通过后跟来做到这一点scipy.linalg.solve_triangular。