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

Question

在https://math.stackexchange.com/a/2233298/340174中提到，如果使用LU分解求解线性方程“M·x = b”（矩阵M为方阵），则速度较慢（使用QR分解更慢）。然而，我注意到numpy.linalg.solve实际上是使用LU分解的。事实上，我想对于非方阵Vandermonde设计矩阵V最小二乘地解决“V·x = b”。我不需要正则化。我看到多种方法：

1. 使用numpy.linalg.lstsq求解“V·x = b”，它使用基于SVD的Fortran“xGELSD”。SVD应该比LU分解更慢，但我不需要计算“(V^T·V)”。
2. 使用numpy.linalg.solve求解“(V^T·V)·x = (V^T·b)”，它使用LU分解。
3. 使用numpy.linalg.solve求解“A·x = b”，它使用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://dev59.com/9qTia4cB1Zd3GeqP8hkE#45535523，似乎solve比lstsq更快，包括计算“V^T·V”，但是当我尝试时，lstsq更快。也许我做错了什么？ 如何最快地解决我的线性问题？

• statsmodels.regression.linear_model.OLS.fit 使用Moore-Penrose伪逆或QR分解+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)”。

- Make42

---

范德蒙矩阵具有解析表达式的逆矩阵。我没有深入研究，在这种情况下不知道浮点误差，但这种方法应该是最快的。此外，您可以查看np.linalg.pinv，这几乎是"(V^T·V)^-1·V^T"。 - bubble

@bubble 这是一个不错的资源，但是这个证明只适用于方阵。OP提到了一个非方的Vandermonde设计矩阵。 - Brenlla

Answer 1

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

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

• numpy.linalg.lstsq() 封装了 LAPACK 的 xGELSD()，如 umath_linalg.c.src 第 2841 行所示。该例程使用分治策略将矩阵 V 缩减为双对角形式，并计算该双对角矩阵的奇异值分解。

• scipy 的 scipy.linalg.lstsq() 封装了 LAPACK 的 xGELSD()、xGELSY() 和 xGELSS()：参数 lapack_driver 可以修改以在这些例程之间切换。根据 LAPACK 的基准测试，xGELSS() 比 xGELSD() 更慢，而 xGELSY() 大约比 xGELSD() 快 5 倍。xGELSY() 利用具有列置换的 V 的 QR 分解。好消息是，这个切换在 scipy 1.1.0 中已经 可用了！

• LAPACK 的 xGELS() 利用矩阵 V 的 QR 分解，但假定该矩阵具有满秩。根据 LAPACK 的基准测试，可以预期 dgels() 比 dgelsd() 快约 5 倍，但它也更容易受到矩阵条件数的影响而变得不精确。有关详细信息和进一步参考，请参见 C++（LAPACK，sgels）和 Python（Numpy，lstsq）结果之间的差异。xGELS() 可在 scipy 的 cython-lapack 接口中使用。

虽然计算和使用 V^T·V 来解决正规方程非常诱人，但很可能不是正确的方法。实际上，该矩阵的条件数会危及精度，大约为 矩阵 V 的条件数的平方。由于 范德蒙矩阵通常情况下具有严重的病态特性，除了离散傅里叶变换的矩阵外，这可能会带来风险... 最后，您甚至可以继续使用 xGELSD() 以避免与条件相关的问题。如果您切换到 xGELSY()，请考虑 估计误差。

- francis

---

我进行了一些基准测试。看看新的答案。 - lovetl2002

Answer 2

我将忽略问题中的Vandermonde部分（评论者bubble指出它具有解析逆），并回答关于其他矩阵的更一般的问题。

我认为这里可能混淆了几件事，因此我将区分以下内容：

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

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

但我认为您实际上在询问3-6。最后一个stackoverflow链接说明6比4更快。正如您所说，4应该比3慢，但4是唯一适用于秩缺失的V的方法。6通常比5更快。

我们应该确保您使用的是6而不是5。要使用6，您需要使用scipy.linalg.solve和assume_a="pos"。否则，您将最终执行5。

我没有在numpy或scipy中找到一个单独执行3的函数。 Lapack例程是xGELS，似乎没有在scipy中公开。您应该能够使用scupy.linalg.qr_multiply，然后使用scipy.linalg.solve_triangular来执行它。

- Subhaneil Lahiri

Answer 3

根据 @francis 的回答，以下是一些基准测试结果。为了达到最佳性能，我关闭了 check_finite。

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from numpy.linalg import lstsq as np_lstsq
from scipy.linalg import lstsq as sp_lstsq
from scipy.linalg import solve

#custom OLS by LU factorization
def ord_lu(X, y):
    A = X.T @ X
    b = X.T @ y
    beta = solve(A, b, overwrite_a=True, overwrite_b=True,
                 check_finite=False)
    return beta

#load iris dataset for testing
data = load_iris()
iris_df = pd.DataFrame(data=data.data, columns= data.feature_names)
species = pd.get_dummies(data.target)
y = iris_df.pop('sepal length (cm)').to_numpy()
X = np.hstack([iris_df.to_numpy(), species.to_numpy()])

%timeit -n10000 -r10 result = np_lstsq(X, y)
#33.7 µs ± 1.19 µs per loop (mean ± std. dev. of 10 runs, 10000 loops each)

%timeit -n10000 -r10 result = sp_lstsq(X, y, lapack_driver='gelsd', check_finite=False)
#44.9 µs ± 1.05 µs per loop (mean ± std. dev. of 10 runs, 10000 loops each)

%timeit -n10000 -r10 result = sp_lstsq(X, y, lapack_driver='gelsy', check_finite=False)
#13.5 µs ± 661 ns per loop (mean ± std. dev. of 10 runs, 10000 loops each)

%timeit -n10000 -r10 result = sp_lstsq(X, y, lapack_driver='gelss', check_finite=False)
#43.5 µs ± 2.11 µs per loop (mean ± std. dev. of 10 runs, 10000 loops each)

%timeit -n10000 -r10 result = ord_lu(X, y)
#18 µs ± 360 ns per loop (mean ± std. dev. of 10 runs, 10000 loops each)

根据结果，gelsy 是最快的最小二乘算法。由于未知原因，在 SciPy 中，gelsd 比 gelss 还要慢，这不应该发生。但是 NumPy 的 lstsq（也使用 gelsd）表现正常，并且比 SciPy 的 gelss 快得多。使用 LU 分解的自定义函数非常快。但正如 @francis 所说，它并不安全。

- lovetl2002