低秩近似

低秩近似 (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 實作如下：

```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)) < tol, break, end
endfor
```

### 變數投影演演算法

交替投影演算法利用低秩近似問題（以像空間參數化）對 $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 實作如下：

```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 即可保證收斂至一個可行解。

## 應用

- 線性系統識別，此時近似矩陣具漢克爾結構。
- 機器學習，此時近似矩陣具有非線性結構。
- 推薦系統，資料矩陣中存在缺失數據且近似為分類變量。
- 距離矩陣還原，附帶正定性約束。
- 自然語言處理，近似矩陣要求非負矩陣。
- 電腦代數，近似矩陣具西爾維斯特結構。

## 外部連結

- C++ package for structured-low rank approximation