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

Question

在https://math.stackexchange.com/a/2233298/340174中提到，如果通过 LU 分解完成线性方程“M·x = b”（矩阵 M 是平方）的求解速度很慢（但使用 QR 甚至更慢）分解）。现在我注意到numpy.linalg.solve实际上是在使用 LU 分解。事实上，我想为最小二乘的非平方 Vandermonde 设计矩阵 V 求解“V·x = b”。我不想正则化。我看到多种方法：numpy.linalg.lstsq使用基于 SVD 的 Fortran "xGELSD"求解 "V·x = b" 。SVD 应该比 LU 分解还要慢，但我不需要计算“(V^T·V)”。numpy.linalg.solve使用 LU 分解求解 "(V^T·V)·x = (V^T·b)" 。用 求解“A·x = b” numpy.linalg.solve，使用LU分解，但直接根据https://math.stackexchange.com/a/3155891/340174计算“A=xV^T·V”或者，我可以使用 scipy 的最新solve版本（https://docs.scipy.org/doc/scipy-1.2.1/reference/generated/scipy.linalg.solve.html），它可以对对称矩阵“A "（我猜这比使用 LU 分解要快），但是我的 scipy 停留在 1.1.0 上，所以我无法访问它。从https://stackoverflow.com/a/45535523/4533188看来，它似乎solve比 快lstsq，包括计算“V^T·V”，但是当我尝试它时，lstsq速度更快。也许我做错了什么？解决线性问题的最快方法是什么？没有真正的选择statsmodels.regression.linear_model.OLS.fit是使用 Moore-Penrose 伪逆或 QR-factorization + np.linalg.inv+ np.linalg.svd+ numpy.linalg.solve，这对我来说似乎不太有效。sklearn.linear_model.LinearRegression 使用 scipy.linalg.lstsq。scipy.linalg.lstsq 也使用 xGELSD。我预计计算“(V^T·V)”的倒数会非常昂贵，所以我放弃了“x = (V^T·V)^-1·(V^T·b)”的直接计算

Answer 1

我将忽略问题的 Vandermonde 部分（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。

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

Answer 2

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

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

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

• scipy'sscipy.linalg.lstsq()包装 LAPACK's xGELSD(), xGELSY()and xGELSS():lapack_driver可以修改参数以从一个切换到另一个。据LAPACK的标杆，xGELSS()是慢xGELSD()和xGELSY()大约5比快xGELSD()。xGELSY()利用列旋转使用 V 的 QR 分解。好消息是这个开关已经在 scipy 1.1.0 中可用！

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

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