首頁猿問解決線性最小二乘法的最快方法

解決線性最小二乘法的最快方法

Python

慕妹3146593 2021-12-17 10:40:22

在https://math.stackexchange.com/a/2233298/340174中提到，如果通過 LU 分解完成線性方程“M·x = b”（矩陣 M 是平方）的求解速度很慢（但使用 QR 甚至更慢）分解）?，F(xiàn)在我注意到numpy.linalg.solve實(shí)際上是在使用 LU 分解。事實(shí)上，我想為最小二乘的非平方 Vandermonde 設(shè)計(jì)矩陣 V 求解“V·x = b”。我不想正則化。我看到多種方法：numpy.linalg.lstsq使用基于 SVD 的 Fortran "xGELSD"求解 "V·x = b" 。SVD 應(yīng)該比 LU 分解還要慢，但我不需要計(jì)算“(V^T·V)”。numpy.linalg.solve使用 LU 分解求解 "(V^T·V)·x = (V^T·b)" 。用求解“A·x = b” numpy.linalg.solve，使用LU分解，但直接根據(jù)https://math.stackexchange.com/a/3155891/340174計(jì)算“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，包括計(jì)算“V^T·V”，但是當(dāng)我嘗試它時(shí)，lstsq速度更快。也許我做錯(cuò)了什么？解決線性問題的最快方法是什么？沒有真正的選擇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。我預(yù)計(jì)計(jì)算“(V^T·V)”的倒數(shù)會(huì)非常昂貴，所以我放棄了“x = (V^T·V)^-1·(V^T·b)”的直接計(jì)算

查看完整描述

2 回答

猛跑小豬

TA貢獻(xiàn)1858條經(jīng)驗(yàn) 獲得超8個(gè)贊

我將忽略問題的 Vandermonde 部分（bubble的評論指出它有一個(gè)解析逆），而是回答關(guān)于其他矩陣的更一般的問題。

我認(rèn)為這里可能會(huì)混淆一些事情，因此我將區(qū)分以下幾點(diǎn)：

V x = b使用LU的精確解
V x = b使用QR的確切解決方案
V x = b使用 QR 的最小二乘解
V x = b使用SVD的最小二乘解
V^T V x = V^T b使用LU的精確解
V^T V x = V^T b使用 Cholesky 的精確解

您鏈接到的第一個(gè) maths.stackexchange 答案是關(guān)于案例 1 和 2。當(dāng)它說 LU 很慢時(shí)，這意味著相對于特定類型矩陣的方法，例如正定矩陣、三角矩陣、帶狀矩陣……

但我認(rèn)為你實(shí)際上是在問 3-6。最后一個(gè) stackoverflow 鏈接指出 6 比 4 快。正如您所說，4 應(yīng)該比 3 慢，但 4 是唯一適用于 rank-deficient 的V。一般來說，6 應(yīng)該比 5 快。

我們應(yīng)該確保你做了 6 而不是 5。要使用 6，你需要使用scipy.linalg.solvewith assume_a="pos"。否則，你最終會(huì)做 5。

我還沒有找到在numpy或中執(zhí)行 3 的單個(gè)函數(shù)scipy。Lapack 例程是 xGELS，它似乎沒有在scipy. 你應(yīng)該可以通過scupy.linalg.qr_multiply后跟來做到這一點(diǎn)scipy.linalg.solve_triangular。

反對回復(fù) 2021-12-17

qq_花開花謝_0

TA貢獻(xiàn)1835條經(jīng)驗(yàn) 獲得超7個(gè)贊

嘗試scipy.linalg.lstsq()使用lapack_driver='gelsy'！

讓我們回顧一下求解線性最小二乘法的不同例程和方法：

numpy.linalg.lstsq()包裝 LAPACK's xGELSD()，如第2841+行的 umath_linalg.c.src 所示。該例程使用分而治之的策略將矩陣 V 簡化為雙對角形式，并計(jì)算該雙對角矩陣的 SVD。
scipy'sscipy.linalg.lstsq()包裝 LAPACK's xGELSD(), xGELSY()and xGELSS():lapack_driver可以修改參數(shù)以從一個(gè)切換到另一個(gè)。據(jù)LAPACK的標(biāo)桿，xGELSS()是慢xGELSD()和xGELSY()大約5比快xGELSD()。xGELSY()利用列旋轉(zhuǎn)使用 V 的 QR 分解。好消息是這個(gè)開關(guān)已經(jīng)在 scipy 1.1.0 中可用！
LAPACKxGELS()使用矩陣 V 的 QR 分解，但它假設(shè)該矩陣具有滿秩。根據(jù) LAPACK 的基準(zhǔn)，on 可以預(yù)期dgels()比快約 5 倍dgelsd()，但它也更容易受到矩陣條件數(shù)的影響并且可能變得不準(zhǔn)確。請參閱C++ (LAPACK, sgels) 和 Python (Numpy, lstsq) 結(jié)果之間的區(qū)別中的詳細(xì)信息和進(jìn)一步參考。xGELS() 在scipy 的 cython-lapack 接口中可用。

雖然非常誘人，但計(jì)算和使用V^T·V來求解正規(guī)方程可能不是要走的路。實(shí)際上，精度受到該矩陣的條件數(shù)的威脅，大約是矩陣 V 的條件數(shù)的平方。由于Vandermonde 矩陣往往是嚴(yán)重病態(tài)的，除了離散傅立葉變換的矩陣，它可能變得危險(xiǎn)......最后，您甚至可以繼續(xù)使用xGELSD()以避免與條件相關(guān)的問題。如果切換到xGELSY()，請考慮估計(jì)錯(cuò)誤。