低秩近似

低秩近似 (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