2.6 低秩矩阵补全

动机

部分观测的数据表包含两种信息：测得的数值，以及尚未测量的位置。矩阵补全利用行与列之间的关系，估计这些缺值。低秩模型用少量共同潜在因子描述关系；同一批因子连接已观测与未观测项，使表内不同位置的信息能够相互影响。

缺值分布与测量数量同样重要。三项观测可以确定一个小型秩一矩阵；换一种观测分布，则可能存在很多补全方式。以下先准确指定哪些观测进入损失，再用凸的奇异值惩罚鼓励较简单的拟合矩阵。算法结合观测残差的梯度步，以及作用于整个矩阵的奇异值软阈值。

定义

设 $P\in\mathbb R^{n\times d}$ 存储可用数据，$\Omega\subseteq\{1,\ldots,n\}\times\{1,\ldots,d\}$ 记录已观测位置。$P$ 在 $\Omega$ 以外的数字只是存储选择，不是拟合目标。观测到的零仍是测量，属于 $\Omega$。

定义

观测投影与矩阵范数

观测算子只保留已观测坐标：

$$(\mathcal P_\Omega X)_{ij}= \begin{cases}X_{ij},&(i,j)\in\Omega,\\0,&(i,j)\notin\Omega.\end{cases} \qquad \mathcal P_{\Omega^c}=I-\mathcal P_\Omega.$$

同形矩阵的 Frobenius 内积与范数是 $\langle X,Y\rangle_F=\sum_{ij}X_{ij}Y_{ij}$、$\|X\|_F^2=\sum_{ij}X_{ij}^2$。若 $X$ 的奇异值为 $\sigma_1,\ldots,\sigma_p$，$p=\min(n,d)$，核范数定义为

$$\|X\|_*=\sum_{j=1}^p\sigma_j,\qquad \|X\|_F^2=\sum_{j=1}^p\sigma_j^2.$$

两个范数描述奇异谱的不同方面。Frobenius 范数量度总平方大小，核范数则合计奇异值大小；秩计算的是正奇异值数量。核范数惩罚通过阈值步把部分奇异值变成零，无须直接搜索离散的秩。

观测几何与低秩因子

定理

观测掩码是正交投影

在 Frobenius 内积下，$\mathcal P_\Omega$ 自伴且幂等。特别地，

$$\|X\|_F^2=\|\mathcal P_\Omega X\|_F^2+ \|\mathcal P_{\Omega^c}X\|_F^2,\qquad \|\mathcal P_\Omega X\|_F\le\|X\|_F.$$

证明

投影到已测量坐标

$\langle\mathcal P_\Omega X,Y\rangle_F$ 与 $\langle X,\mathcal P_\Omega Y\rangle_F$ 均等于 $\sum_{(i,j)\in\Omega}X_{ij}Y_{ij}$，所以算子自伴。同一掩码作用两次不变，所以幂等。已观测部分与未观测部分的非零位置不重叠，内积为零；勾股定理给出范数等式及不等式。

秩不超过 $r$ 的模型可写成 $X=LR^T$，其中 $L\in\mathbb R^{n\times r}$、$R\in\mathbb R^{d\times r}$。$X_{ij}$ 是 $L$ 第 $i$ 行与 $R$ 第 $j$ 行的内积。这些共同坐标让缺值依赖解释观测的同一批因子。直接搜索两组因子是非凸问题，因为它们相互相乘。

因子坐标本身并不唯一。任取可逆 $Q\in\mathbb R^{r\times r}$，把因子换成 $LQ$ 和 $RQ^{-T}$，乘积 $LR^T$ 不变。完整矩阵也可能无法由观测唯一确定。讨论恢复时，需要区分因子坐标自由与缺值是否唯一：低秩提供结构假设，观测分布仍须有足够联系。

例题

1. 确定的缺值与含糊的观测分布

设已知项构成

$$X=\begin{bmatrix}2&b\\3&6\end{bmatrix}.$$

秩一要求 $\det X=12-3b=0$，所以唯一得到 $b=4$。完成矩阵可分解为 $(2,3)^T(1,2)$。

若只观测两个对角项，则全部矩阵

$$X_t=\begin{bmatrix}2&t\\12/t&6\end{bmatrix},\qquad t\ne0,$$

都是秩一且符合观测。例如 $t=2$ 与 $t=4$ 给出不同缺值。仅靠低秩无法选出真实数据表；正则化可以表达对拟合的偏好，却不会增加观测信息。

只计算观测项的损失

允许测量含噪声时，考虑

$$\Phi(X)=f(X)+\lambda\|X\|_*,\qquad f(X)=\frac12\|\mathcal P_\Omega(X-P)\|_F^2,\qquad\lambda>0.$$

只有已测量残差进入 $f$。改写 $P$ 未使用的项不改变目标。核范数作用于整个候选矩阵，用奇异方向连接各坐标。两项都凸，但不完整观测下，整体目标未必严格凸。

参数控制拟合与复杂度的权衡。较大参数更重视缩小奇异值，常会移除更多奇异方向；观测位置的拟合值也可偏离测量。精确测量可另建约束模型：在 $\mathcal P_\Omega X=\mathcal P_\Omega P$ 下最小化 $\|X\|_*$。本节推导的是上述惩罚模型，观测残差仍参与最优化。

核范数是凸矩阵范数。它也可写成：在谱范数不超过一的矩阵中，与 $X$ 作内积的最大值。奇异值迹不等式把内积界于奇异值总和，对齐奇异方向即可达到上界；线性函数的最大值是凸函数。这把惩罚的凸性联系到以下同一套奇异向量几何。

奇异值软阈值

定义

对奇异值作软阈值

设 $Y=U\operatorname{diag}(s_j)V^T$、$\tau\ge0$，定义

$$\operatorname{SVT}_\tau(Y) =U\operatorname{diag}\big((s_j-\tau)_+\big)V^T, \qquad (a)_+=\max(a,0).$$

定理

核范数的近端问题

对任意完整指定的矩阵 $Y$ 及 $\tau\ge0$，

$$\frac12\|X-Y\|_F^2+\tau\|X\|_*$$

的唯一最小化矩阵是 $\operatorname{SVT}_\tau(Y)$。

证明

把矩阵问题化为逐个奇异值的阈值

设 $X$、$Y$ 已排序奇异值分别为 $t_j$、$s_j$。von Neumann 迹不等式给出 $\langle X,Y\rangle_F\le\sum_jt_js_j$。展开平方距离，目标因此至少是

$$\sum_j\left[\frac12(t_j-s_j)^2+\tau t_j\right],\qquad t_j\ge0.$$

每项在 $t_j=(s_j-\tau)_+$ 最小；此运算保持排序。把 $X$ 的奇异向量与 $Y$ 对齐，迹不等式取等号，并同时达到全部标量最小值。平方距离严格凸，所以最小化矩阵唯一，即使奇异向量基未必唯一。

软阈值从每个正奇异值减去相同数量，把小于或等于阈值的奇异值设成零；秩截断则把选中的大奇异值原样保留。两者都可能产生较低秩矩阵，但解不同目标。遇到重复奇异值，数值 SVD 可换用该子空间的另一正交归一基；在整个子空间施加同一阈值，重建矩阵仍相同。

例题

2. 完整观测的阈值计算

取 $Y=\operatorname{diag}(3,1)$、$\tau=1/2$。奇异方向为坐标轴，解是 $X=\operatorname{diag}(5/2,1/2)$。平方距离项为 $\tfrac12(1/4+1/4)=1/4$，惩罚为 $\tfrac12(5/2+1/2)=3/2$，目标总值是 $7/4$。保持 $Y$ 不变则得到 $2$，所以精确符合数据并非此惩罚问题的最小值；两个观测奇异值都被缩小。

带掩码的近端梯度步

观测损失的梯度为

$$\nabla f(X)=\mathcal P_\Omega(X-P),\qquad \|\nabla f(X)-\nabla f(Y)\|_F\le\|X-Y\|_F.$$

不等式来自投影定理，所以梯度的 Lipschitz 常数可取一。步长 $0<\eta\le1$ 给出

$$X^+=\operatorname{SVT}_{\eta\lambda} \left(X-\eta\mathcal P_\Omega(X-P)\right).$$

单位步长时，化为

$$X^+=\operatorname{SVT}_\lambda \left(\mathcal P_\Omega P+\mathcal P_{\Omega^c}X\right).$$

梯度步的中间矩阵在观测处采用测量，其余保留当前估计，然后对整个矩阵施加阈值。结果可以改变观测项并填入缺值，符合残差与核范数的权衡。缺值估计会留到下一个中间矩阵，影响下一次 SVD，而不是每步都被视为测得的零。

为证明步长合理，对增量 $D$ 展开二次损失：

$$f(X+D)=f(X)+\langle\nabla f(X),D\rangle_F +\tfrac12\|\mathcal P_\Omega D\|_F^2 \le f(X)+\langle\nabla f(X),D\rangle_F+\tfrac12\|D\|_F^2.$$

加上核范数后，得到在 $X$ 处等于 $\Phi$ 的二次上界模型。由近端定理，其最小值正是单位步更新，所以更新不增加 $\Phi$。较小步长把二次项系数换成 $1/(2\eta)$，仍是有效上界。步长限制直接来自掩码的范数。

惩罚大小与零解的精确判据

令 $Y=\mathcal P_\Omega P$，用 $\|Y\|_2$ 表示其最大奇异值。迹不等式给出

$$\Phi(X)-\Phi(0) =\tfrac12\|\mathcal P_\Omega X\|_F^2-\langle Y,X\rangle_F +\lambda\|X\|_* \ge\tfrac12\|\mathcal P_\Omega X\|_F^2 +(\lambda-\|Y\|_2)\|X\|_*.$$

因此，参数至少等于此最大奇异值时，零矩阵是全局最小值。这精确解释了为什么迭代会移除全部奇异方向：模型对矩阵大小所收取的代价，已超过任何观测拟合改善。这个尺度来自已观测表及其掩码，不依赖想象中的完整矩阵。比较一系列参数时，它提供明确端点；逐步减小参数，才允许更显著的拟合结构。

不动点证明什么

非光滑惩罚目标的最优条件不只看梯度。核范数在 $X$ 的次梯度 $G$，意指对所有 $Z$ 都有 $\|Z\|_*\ge\|X\|_*+\langle G,Z-X\rangle_F$。单位步更新满足

$$\frac{X-X^+-\nabla f(X)}{\lambda}\in\partial\|X^+\|_*.$$

若更新返回相同矩阵，这条支撑不等式就把观测梯度与惩罚平衡。加上凸函数 $f$ 的一阶不等式，即证明任意候选的目标均不小于当前值；不动点是指定凸目标的全局最小值。此时观测残差仍可非零，因为其梯度由核范数次梯度平衡。

所以，除了目标值，也可检查近端更新的大小。变化很小度量的是整个惩罚问题的近似驻定性，不要求观测梯度本身为零，也不等于已经核对缺值的预测准确度。

例题

3. 一次迭代填入缺值

沿用首例的三项观测，存储 $P=\begin{bmatrix}2&0\\3&6\end{bmatrix}$，取 $X_0=0$、$\lambda=1$。右上缺值不在 $\Omega$ 中。中间矩阵 $Y=P$ 的奇异值为 $(\sqrt{73}+5)/2$、$(\sqrt{73}-5)/2$，两者均大于一。其正交极分解因子为

$$UV^T=\frac1{\sqrt{73}}\begin{bmatrix}8&-3\\3&8\end{bmatrix},\qquad X_1=\begin{bmatrix} 2-8/\sqrt{73}&3/\sqrt{73}\\ 3-3/\sqrt{73}&6-8/\sqrt{73} \end{bmatrix}.$$

缺值由存储的零变成约 $0.351123$，观测项也移动。目标从 $49/2=24.5$ 降至 $\sqrt{73}-2+137/146\approx7.482360$。精确秩一补全的缺值为 $4$；首次惩罚迭代给出的是拟合与惩罚准则下的中间估计。

可执行计算

布尔掩码独立记录观测状态。以下实现使用单位步长，以及刚推导的目标。

import numpy as np

observed = np.array([[True, False], [True, True]])
P = np.array([[2., 0.], [3., 6.]])
lam = 1.

def svt(Y, threshold):
    U, s, Vt = np.linalg.svd(Y, full_matrices=False)
    return (U * np.maximum(s - threshold, 0.)) @ Vt

def objective(X):
    residual = np.where(observed, X - P, 0.)
    return 0.5 * np.sum(residual**2) + lam * np.linalg.svd(
        X, compute_uv=False
    ).sum()

X = np.zeros_like(P)
first = svt(P, lam)
history = [objective(X)]
for _ in range(100):
    gradient = np.where(observed, X - P, 0.)
    X = svt(X - gradient, lam)
    history.append(objective(X))

assert np.all(np.diff(history) <= 1e-10)
print(np.round(first, 6))
print(round(history[0], 6), round(history[-1], 6))

首个输出矩阵约为 $\begin{bmatrix}1.063671&0.351123\\2.648877&5.063671\end{bmatrix}$；目标值输出为 $24.5$ 和约 $7.073338$。固定迭代次数使例子可重现。较大型应用可同时记录目标与迭代变化量，决定何时停止；未提供数值的位置，其预测准确度需要额外验证数据。

参数尺度与拟合评估

参数数值受数据及损失的尺度影响。把测量乘以正常数 $c$ 时，也把参数乘以同一常数；候选矩阵同样缩放后，残差项与惩罚都乘以 $c^2$，全部候选目标的排序不变。这说明更换测量单位后等价模型怎样变化，而不是给某个参数数值一个绝对意义。

若把观测平方损失总和改成均值，要同时把惩罚系数除以观测数量，才得到相同最优解。只平均损失却保留原系数，会提高惩罚的相对影响。比较不同实现或观测数量不同的数据表前，先说明损失的归一化方式。

评估预测时，可留出部分真正测得的项作验证，用其余项建立拟合掩码。梯度步不使用验证位置，再用预测值与留出测量比较，就度量未参与拟合的观测之预测表现。这也改变了拟合的观测分布，因此划分时要保留模型所需的行列信息。最优化诊断检查是否已最小化选定目标，预测诊断则评估估计值是否有用。

常见错误

常见错误

用数值是否为零判断观测状态

观测零也产生平方残差。缺值则不产生残差，与存储的数字无关。掩码必须由观测状态建立。更新顺序也按推导：先作梯度步的中间矩阵，再作奇异值软阈值，才能对应指定的惩罚目标与步长论证。

总结

观测算子是正交投影，决定损失和梯度。低秩结构连接缺值，可辨识性则取决于测量位置。核范数给出凸奇异值惩罚，其近端映射是奇异值软阈值。把它与观测梯度结合，就得到可由二次上界直接解释单位步长的算法。

练习

思考检查

1. 秩不超过一的表已观测 $X_{11}=0$、$X_{21}=2$、$X_{22}=4$。求缺值 $X_{12}$，并解释观测零为什么属于掩码。

使用行列式及观测定义。

思考检查

2. 完整观测 $Y=\operatorname{diag}(5,2)$，惩罚参数 $\tau=3$。求核范数近端解、其秩及最小目标值。

分别计算平方残差项和核范数惩罚。

思考检查

3. 在掩码例子中，把惩罚改为 $\lambda=7$ 并从零开始。首次近端梯度更新是什么？此目标是否强制精确符合三项观测？

比较阈值与已列出的较大奇异值。

解答

解答 · 1. 零测量也提供信息

行列式为 $0\cdot4-2X_{12}$，所以秩不超过一要求 $X_{12}=0$。第二行非零，完成矩阵秩为一。若漏掉 $X_{11}=0$，就移除了一条实际方程，允许其他秩一补全。数值为零不等于缺失。

解答 · 2. 缩小奇异值

解为 $\operatorname{diag}(2,0)$，秩一。平方距离项为 $\tfrac12(3^2+2^2)=13/2$，惩罚为 $3\cdot2=6$，最小目标值是 $25/2=12.5$。

解答 · 3. 首步移除全部奇异值

中间矩阵的最大奇异值为 $(\sqrt{73}+5)/2<7$，阈值七给出零矩阵。模型对观测残差施加惩罚，允许拟合与测量不同；观测相等并非其约束。

相关笔记

• 奇异值分解与主成分
• 最小二乘与残差
• 二次函数的梯度下降
• 课程读物