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

Question

在https://math.stackexchange.com/a/2233298/340174 中提到，如果通过 LU 分解求解线性方程“M·x = b”（矩阵 M 为正方形）会很慢（但使用 QR 分解会更慢）。现在我注意到numpy.linalg.solve 实际上是在使用 LU 分解。事实上，我想为最小二乘的非平方范德蒙德设计矩阵 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 (https://docs.scipy.org/doc/scipy-1.2.1/reference/generated/scipy.linalg.solve.html) 的最新 solve，它可以对对称矩阵“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 分解 + 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使用二维码的精确解
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.solve 和 assume_a="pos"。否则，你最终会做 5 个。

我还没有找到在numpy 或scipy 中执行 3 的单个函数。 Lapack 例程是 xGELS，在scipy 中似乎没有暴露。您应该可以使用scupy.linalg.qr_multiply 后跟scipy.linalg.solve_triangular 来完成此操作。

Answer 2

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

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

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

• scipy 的 scipy.linalg.lstsq() 包装了 LAPACK 的 xGELSD()、xGELSY() 和 xGELSS()：参数 lapack_driver 可以修改为从一个切换到另一个。根据 LAPACK 的基准，xGELSS() 比 xGELSD() 慢，xGELSY() 比 xGELSD() 快大约 5。 xGELSY() 使用带有列旋转的 V 的 QR 分解。好消息是这个开关是already available in scipy 1.1.0！

• LAPACK 的xGELS() 利用矩阵 V 的 QR 分解，但它假设该矩阵具有满秩。根据 LAPACK 的 benchmark，on 可以预期 dgels() 比 dgelsd() 快 5 倍左右，但它也更容易受到矩阵条件数的影响，可能会变得不准确。请参阅The difference between C++ (LAPACK, sgels) and Python (Numpy, lstsq) results 中的详细信息和更多参考资料。 xGELS() 在cython-lapack interface of scipy 中可用。

虽然很诱人，但计算和使用V^T·V 求解正规方程可能不是可行的方法。实际上，精度受到该矩阵的条件数的影响，大约the square of the condition number of the matrix V。由于Vandermonde matrices tend to be badly ill-conditioned, except for the matrices of the discrete Fourier transform，它可能会变得危险...... 最后，您甚至可以继续使用xGELSD() 以避免与调节相关的问题。如果您切换到xGELSY()，请考虑estimating the error。

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 所说，这并不安全。