低秩近似 (low-rank approximation) 是指用一個較低秩的矩陣去近似給定矩陣的過程。更精確地說，它是一個最佳化問題，其中損失函數衡量給定矩陣（資料）與近似矩陣（最佳化變數）之間的擬合程度，並且附帶近似矩陣秩的約束條件。此問題常用於數學模型建構與資料壓縮。秩的約束與對符合資料模型複雜度的限制相關。在應用中，近似矩陣通常還會有其他約束，例如非負性及漢克爾結構。
低秩近似與多種其他技術密切相關，包括主成分分析、因素分析、潛在語義分析、最小全平方法、正交迴歸及動態模態分解。
定義
給定
- 結構映射 \mathcal{S} : \mathbb{R}^{n_p} \to \mathbb{R}^{m\times n},
- 結構參數向量 p\in\mathbb{R}^{n_p},
- 範數 \| \cdot \|，以及
- 希望的秩 r，
:
\min_{\widehat p} \|p - \widehat p\|
\quad \text{subject to} \quad \operatorname{rank}\big(\mathcal{S}(\widehat p)\big) \le r
基本低秩近似
無結構且擬合度以弗羅貝尼烏斯範數衡量的問題，即
:
\min_{\widehat D} \|D - \widehat D\|_{\mathrm{F}}
\quad \text{subject to} \quad \operatorname{rank}(\widehat D) \le r
有解析解，該解可由資料矩陣的奇異值分解得到。此結果稱為矩陣近似引理或Eckart–Young–Mirsky 定理。此問題最初由埃哈德·施密特在無限維積分算子情境中提出（其方法可擴展至希爾伯特空間中任意緊算子），後由C. Eckart與G. Young重新發現。 L. Mirsky將結果推廣至任意單位不變範數。 設
:
D = U\Sigma V^{\top} \in \mathbb{R}^{m\times n}, \quad m \geq n
為 D 的奇異值分解，其中 \Sigma=:\operatorname{diag}(\sigma_1,\ldots,\sigma_r)，r\leq \min\{m,n\}=n，為一個 m\times n 的矩形對角矩陣，具有 r 個非零奇異值，且
\sigma_1 \geq \cdots \geq \sigma_r > \sigma_{r+1} = \cdots = \sigma_n = 0。對於給定的 k \in \{1,\dots,r\}，將 U、\Sigma 和 V 分割為：
:
U =: \begin{bmatrix} U_1 & U_2\end{bmatrix}, \quad
\Sigma =: \begin{bmatrix} \Sigma_1 & 0 \\ 0 & \Sigma_2 \end{bmatrix}, \quad\text{且}\quad
V =: \begin{bmatrix} V_1 & V_2 \end{bmatrix},
其中 U_1 是 m \times k，\Sigma_1 是 k \times k，V_1 是 n \times k。則透過截斷奇異值分解得到的秩為 k 的矩陣為：
:
\widehat D^* = U_1 \Sigma_1 V_1^{\top},
且滿足
:
\|D-\widehat D^*\|_{\text{F}} = \min_{\operatorname{rank}(\widehat D) \leq k} \|D-\widehat D\|_{\text{F}} = \sqrt{\sigma_{k+1}^2 + \cdots + \sigma_r^2}.
當且僅當 \sigma_k > \sigma_{k+1} 時，最小化解 \widehat D^* 唯一。
Eckart–Young–Mirsky 定理的證明（針對譜範數）
設 A \in \mathbb{R}^{m \times n} 為一個實數（可能非方陣）矩陣，且 m \leq n。假設
:
A = U \Sigma V^\top
為 A 的奇異值分解，其中 U 和 V 為正交矩陣，\Sigma 為 m \times n 的對角矩陣，對角線元素為 (\sigma_1, \sigma_2, \ldots, \sigma_m)，且 \sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_m \geq 0。
我們宣稱，在譜範數 \|\cdot\|_2 下，對 A 的最佳秩為 k 近似為：
:
A_k := \sum_{i=1}^k \sigma_i u_i v_i^\top
其中 u_i 和 v_i 分別為 U 與 V 的第 i 欄。
首先，注意
:
\|A - A_k\|_2 = \left\|\sum_{i=1}^n \sigma_i u_i v_i^\top - \sum_{i=1}^k \sigma_i u_i v_i^\top \right\|_2 = \left\|\sum_{i=k+1}^n \sigma_i u_i v_i^\top \right\|_2 = \sigma_{k+1}.
因此，我們要證明若 B_k = XY^\top，其中 X 與 Y 均有 k 欄，則
:
\|A - A_k\|_2 = \sigma_{k+1} \leq \|A - B_k\|_2.
由於 Y 有 k 欄，必存在第一個 k+1 欄的非平凡線性組合
:
w = \gamma_1 v_1 + \cdots + \gamma_{k+1} v_{k+1},
使得 Y^\top w = 0。不失一般性，將 w 正規化，使得 \|w\|_2 = 1 或等價地
\gamma_1^2 + \cdots + \gamma_{k+1}^2 = 1。因此，
:
\|A - B_k\|_2^2 \geq \|(A - B_k) w\|_2^2 = \|A w\|_2^2 = \gamma_1^2 \sigma_1^2 + \cdots + \gamma_{k+1}^2 \sigma_{k+1}^2 \geq \sigma_{k+1}^2.
由上述不等式兩邊取平方根即得證。
Eckart–Young–Mirsky 定理的證明（針對 弗羅貝尼烏斯範數）
設 A \in \mathbb{R}^{m \times n} 為一個實數矩陣（可能為長方形矩陣），且 m \leq n。假設
:
A = U \Sigma V^\top
為 A 的奇異值分解。
我們主張，對於 Frobenius 範數（記為 \|\cdot\|_F），對 A 的最佳秩為 k 的近似矩陣為
:
A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top
其中 u_i 和 v_i 分別為矩陣 U 和 V 的第 i 欄。
首先注意，我們有
:
\|A-A_k\|_F^2 = \left\|\sum_{i=k+1}^n \sigma_i u_iv_i^\top \right\|_F^2 = \sum_{i=k+1}^n \sigma_i^2
因此，我們需要證明：若 B_k = X Y^\top，其中 X 和 Y 均有 k 欄，則
: \|A-A_k\|_F^2=\sum_{i=k+1}^n\sigma_{i}^2\leq \|A-B_k\|_F^2.
利用譜範數的三角不等式，若 A = A' + A*，則
:
\sigma_1(A) \leq \sigma_1(A') + \sigma_1(A*).
設 A'_k 與 A_k 分別為上述奇異值分解所得 A' 與 A 的秩為 k 的近似矩陣。則對任意 i, j \geq 1 有
:
\begin{align}
\sigma_i(A')+\sigma_j(A) &= \sigma_1(A'-A'_{i-1}) + \sigma_1(A-A*_{j-1})\\
&\geq \sigma_1(A-A'_{i-1}-A*_{j-1})\\
&\geq \sigma_1(A-A_{i+j-2})\qquad (\because {\rm rank}(A'_{i-1}+A*_{j-1})\leq i+j-2))\\
&=\sigma_{i+j-1}(A).
\end{align}
因為 \sigma_{k+1}(B_k) = 0，令 A' = A - B_k，A* = B_k，我們得到對所有 i \geq 1 有
:
\sigma_i(A - B_k) \geq \sigma_{k+i}(A).
因此，
:
\|A-B_k\|_F^2 = \sum_{i=1}^n \sigma_i(A-B_k)^2 \geq \sum_{i=k+1}^n\sigma_i(A)^2 = \|A-A_k\|_F^2,
證畢。
加權低秩近似
Frobenius 範數對近似誤差 D - \widehat{D} 的所有元素均等加權。若想考慮誤差分佈的先驗知識，可利用加權低秩近似問題：
:
\min_{\widehat D} \quad
\operatorname{vec}(D - \widehat{D})^\top W \operatorname{vec}(D - \widehat{D})
\quad\text{subject to}\quad \operatorname{rank}(\widehat{D}) \leq r,
其中 \operatorname{vec}(A) 表示將矩陣 A 按欄向量化，而 W 是給定的正（半）定權重矩陣。
一般加權低秩近似問題無法透過奇異值分解求得解析解，需使用局部優化方法，且無法保證能找到全域最優解。
當權重為無相關性時，加權低秩近似問題亦可如下表示： 對非負矩陣 W 及矩陣 A，目標為最小化
:
\sum_{i,j} \big( W_{i,j} (A_{i,j} - B_{i,j}) \big)^2
在秩不超過 r 的矩陣 B 上。
秩約束的像空間與核空間表示
利用等價關係
:
\operatorname{rank}(\widehat D) \leq r
\quad\iff\quad
\exists P \in \mathbb{R}^{m \times r},\ L \in \mathbb{R}^{r \times n} \text{ s.t. } \widehat D = PL
及
:
\operatorname{rank}(\widehat D) \leq r
\quad\iff\quad
\exists R \in \mathbb{R}^{(m - r) \times m},\ \operatorname{rank}(\widehat R) = m - r \text{ s.t. } \ R \widehat D = 0
加權低秩近似問題可轉換為參數優化問題
:
\min_{\widehat D, P, L} \quad
\operatorname{vec}^{\top}(D - \widehat D) W \operatorname{vec}(D - \widehat D)
\quad\text{subject to}\quad \widehat D = PL
及
:
\min_{\widehat D, R} \quad
\operatorname{vec}^{\top}(D - \widehat D) W \operatorname{vec}(D - \widehat D)
\quad\text{subject to}\quad R \widehat D = 0, \quad RR^{\top} = I_r,
其中 I_r 為大小為 r 的單位矩陣。
交替投影演算法
秩約束的像空間表示啟發了一種參數優化方法，即交替固定其中一組變數（P 或 L）最小化目標函數。雖然同時對 P 與 L 最小化是困難的雙凸最佳化問題，但單獨固定一組變數的最小化是線性最小平方法問題，可有效且全域求解。
該演算法稱為交替投影（alternating projections），在加權低秩近似問題中可保證全域收斂且收斂速率為線性，收斂至局部最優解。使用者需給定初始 P（或 L）值，並在滿足預設收斂條件時停止迭代。
Matlab 實作如下：
function [dh, f] = wlra_ap(d, w, p, tol, maxiter)
[m, n] = size(d); r = size(p, 2); f = inf;
for i = 2:maxiter
% minimization over L
bp = kron(eye(n), p);
vl = (bp' w bp) \ bp' w d(:);
l = reshape(vl, r, n);
% minimization over P
bl = kron(l', eye(m));
vp = (bl' w bl) \ bl' w d(:);
p = reshape(vp, m, r);
% check exit condition
dh = p * l; dd = d - dh;
f(i) = dd(:)' w dd(:);
if abs(f(i - 1) - f(i))
變數投影演算法
交替投影演算法利用低秩近似問題（以像空間參數化）對 P 與 L 為雙線性的特性。變數投影法（variable projections）則進一步利用此雙線性結構。
考慮加權低秩近似問題的像空間參數化，對 L 變數（線性最小平方法）求解得關於 P 的近似誤差封閉式
:
f(P) = \sqrt{\operatorname{vec}^{\top}(D)\Big(
W - W (I_n \otimes P) \big( (I_n \otimes P)^{\top} W (I_n \otimes P) \big)^{-1} (I_n \otimes P)^{\top} W
\Big) \operatorname{vec}(D)}.
因此原問題等價於非線性最小平方問題，透過標準優化方法（如萊文伯格-馬夸特方法）可求解 f(P) 的最小值。
Matlab 實作如下：
function [dh, f] = wlra_varpro(d, w, p, tol, maxiter)
prob = optimset(); prob.solver = 'lsqnonlin';
prob.options = optimset('MaxIter', maxiter, 'TolFun', tol);
prob.x0 = p; prob.objective = @(p) cost_fun(p, d, w);
[p, f ] = lsqnonlin(prob);
[f, vl] = cost_fun(p, d, w);
dh = p * reshape(vl, size(p, 2), size(d, 2));
function [f, vl] = cost_fun(p, d, w)
bp = kron(eye(size(d, 2)), p);
vl = (bp' w bp) \ bp' w d(:);
f = d(:)' w (d(:) - bp * vl);
變數投影法亦可應用於核空間參數化的低秩近似問題。當消除變數數量遠大於剩餘優化變數時，該方法尤為有效。此類問題常見於以核形式參數化的系統識別，消除的變數為近似軌跡，剩餘變數為模型參數。在線性時不變系統（LTI）理論中，消除步驟相當於卡尔曼滤波（Kalman smoothing）。
元素逐點的 Lp 低秩近似問題
設 \|A\|_p = \left(\sum_{i,j}|A_{i,j}^p|\right)^{1/p} 。當 p=2 時，最快演算法時間為 nnz(A) + n \cdot \mathrm{poly}(k/\epsilon)。 其中一個重要概念為「Oblivious Subspace Embedding (OSE)」，最早由 Sarlos 提出。
當 p=1 時，元素逐點的 L1 範數在面對離群值時比 Frobenius 範數更具魯棒性，適用於雜訊不服從高斯分布的模型。很自然的會想要最小化 \|B - A\|_1。 對於 p=0 和 p \geq 1，已有部分帶有理論保證的演算法。
距離低秩近似問題
設 P = { p_1 , \ldots, p_m } 與 Q = { q_1, \ldots, q_n } 為任意度量空間中的兩組點集。令 A 表示 m \times n 的矩陣，其中元素定義為 A_{i,j} = \operatorname{dist}(p_i,q_j)。此類距離矩陣常見於軟體套件中，並應用於學習影像流形（image manifolds）、手寫辨識（handwriting recognition）及多維展開（multi-dimensional unfolding）。為了嘗試縮減其描述大小， 可研究此類矩陣的低秩近似。
凸限制低秩近似
通常，我們不僅希望解為低秩，還要符合應用需求的其他凸限制。問題可表述為：
:
\text{minimize} \quad \text{over } \widehat p \quad \|p - \widehat p\| \quad\text{subject to}\quad \operatorname{rank}\big(\mathcal{S}(\widehat p)\big) \leq r \text{ and }
g(\widehat p) \leq 0
此問題有多項實際應用，包括從不精確的（半正定規劃）鬆弛中恢復良好解。若額外限制 g(\widehat p) \leq 0 為線性（如要求元素非負），此問題稱為結構化低秩近似（structured low-rank approximation）。 更一般的形式則稱為凸限制低秩近似（convex-restricted low rank approximation）。
此問題雖具有實用價值，但由於同時涉及凸限制與非凸的低秩限制，因此具有相當的挑戰性。根據不同形式的限制條件 g(\widehat p) \leq 0，已有多種技術被提出。值得注意的是，交替方向乘子法（ADMM）可應用於目標函數為凸、同時具有秩限制與其他凸限制的非凸問題中，，因此特別適合處理上述類型的問題。此外，與一般非凸問題不同，只要對偶變數在迭代過程中收斂，ADMM 即可保證收斂至一個可行解。
應用
- 線性系統識別，此時近似矩陣具漢克爾結構。
- 機器學習，此時近似矩陣具有非線性結構。
- 推薦系統，資料矩陣中存在缺失數據且近似為分类变量。
- 距離矩陣還原，附帶正定性約束。
- 自然語言處理，近似矩陣要求非負矩陣。
- 計算機代數，近似矩陣具西爾維斯特結構。
引用
外部連結
*[https://github.com/slra/slra C++ package for structured-low rank approximation]
评论 (0)