NumPy 如何解決欠定係統的最小二乘法?
假設我們有形狀為 (2, 5) 的 X 和形狀為 (2,)
的 y
這有效:
np.linalg.lstsq(X, y)
我們希望這僅在 X 的形狀為 (N,5) 其中 N>=5 時才有效。但是為什麼以及如何?
我們確實按預期返回了 5 個權重,但是這個問題是如何解決的呢?
是不是就像我們有 2 個方程和 5 個未知數?
numpy 怎麼能解決這個問題?
它必須做一些像插值這樣的事情來創建更多的人工方程?…
我的理解是numpy.linalg.lstsq依賴於LAPACK例程dgelsd。
問題是要解決:
$$ \text{minimize} (\text{over} ; \mathbf{x}) \quad | A\mathbf{x} - \mathbf{b} |_2 $$
當然,這對於秩小於向量長度的矩陣 A沒有唯一解 $ \mathbf{b} $ . 在未確定係統的情況下,
dgelsd
提供解決方案 $ \mathbf{z} $ 這樣:
- $ A\mathbf{z} = \mathbf{b} $
- $ | \mathbf{z} |_2 \leq |\mathbf{x} |_2 $ 對所有人 $ \mathbf{x} $ 滿足 $ A\mathbf{x} = \mathbf{b} $ . (IE $ \mathbf{z} $ 是待定係統的最小範數解。
例如,如果系統是 $ x + y = 1 $ , numpy.linalg.lstsq 將返回 $ x = .5, y = .5 $ .
dgelsd 是如何工作的?
該例程
dgelsd
計算 A 的奇異值分解(SVD)。我將概述使用 SVD 解決線性系統的想法。奇異值分解是因式分解 $ U \Sigma V' = A $ 在哪裡 $ U $ 和 $ V $ 是正交矩陣和 $ \Sigma $ 是一個對角矩陣,其中對角元素被稱為奇異值。
矩陣的有效秩 $ A $ 將是有效非零的奇異值的數量(即,相對於機器精度等與零有很大差異……)。讓 $ S $ 是非零奇異值的對角矩陣。因此,SVD 為:
$$ A = U \begin{bmatrix} S & 0 \ 0 & 0 \end{bmatrix} V' $$
的偽逆 $ A $ 是(誰)給的:
$$ A^\dagger = V \begin{bmatrix} S^{-1} & 0 \ 0 & 0 \end{bmatrix} U' $$
考慮解決方案 $ \mathbf{x} = A^\dagger \mathbf{b} $ . 然後:
$$ \begin{align*} A\mathbf{x} - \mathbf{b} &= U \begin{bmatrix} S & 0 \ 0 & 0 \end{bmatrix} V' V \begin{bmatrix} S^{-1} & 0 \ 0 & 0 \end{bmatrix} U' \mathbf{b} - \mathbf{b} \ &= U \begin{bmatrix} I & 0 \ 0 & 0 \end{bmatrix} U' \mathbf{b} - \mathbf{b}\ \end{align*} $$
這里基本上有兩種情況:
- 非零奇異值的個數(即矩陣的大小 $ I $ ) 小於的長度 $ \mathbf{b} $ . 這裡的解決方案並不准確;我們將以最小二乘的方式求解線性系統。
- $ A\mathbf{x} - \mathbf{b} = \mathbf{0} $
最後一部分有點棘手……需要跟踪矩陣尺寸並使用它 $ U $ 是一個正交矩陣。
偽逆等價
什麼時候 $ A $ 具有線性獨立的行,(例如,我們有一個胖矩陣),然後: $$ A^\dagger = A'\left(AA' \right)^{-1} $$
對於未定係統,您可以證明偽逆為您提供了最小範數解。
什麼時候 $ A $ 具有線性獨立的列,(例如,我們有一個瘦矩陣),然後: $$ A^\dagger = \left(A’A \right)^{-1}A' $$