低秩矩阵补全（2）：广义似然模型

什么是广义似然下的低秩矩阵补全

第一篇讨论的低秩矩阵补全主要对应平方损失，即观测满足 $Y_{ij}=M^\ast_{ij}+\varepsilon_{ij}$，并用 $\frac12(Y_{ij}-M_{ij})^2$ 衡量拟合误差。这一模型适合连续型、近似高斯噪声的数据，但并不适合二元点击、计数、评分等级、曝光转化等非高斯观测。广义似然模型（generalized likelihood model）的目标，是把低秩结构与更合适的观测分布结合起来。

设 $\Theta^\ast\in\mathbb R^{n_1\times n_2}$ 为潜在自然参数矩阵（natural parameter matrix），观测集合为 $\Omega\subset[n_1]\times[n_2]$。对每个 $(i,j)\in\Omega$，假设 $Y_{ij}$ 条件独立，并且其条件分布属于一参数指数族：

$$p(y\mid \theta)=h(y)\exp\{y\theta-b(\theta)\}.$$

这里 $b$ 是 log-partition function，$\mathbb E(Y_{ij}\mid\Theta^\ast_{ij})=b'(\Theta^\ast_{ij})$，$\operatorname{Var}(Y_{ij}\mid\Theta^\ast_{ij})=b''(\Theta^\ast_{ij})$。因此，矩阵补全的对象可以是自然参数矩阵 $\Theta^\ast$，也可以是均值矩阵 $\mathbf M^\ast=b'(\Theta^\ast)$。本文默认低秩结构施加在 $\Theta^\ast$ 上，即 $\operatorname{rank}(\Theta^\ast)=r\ll\min(n_1,n_2)$。

这一框架包含若干常见模型：高斯模型对应 $b(\theta)=\theta^2/2$，负对数似然等价于平方损失；Bernoulli-logistic 模型对应 $b(\theta)=\log(1+e^\theta)$，适合二元观测；Poisson 模型对应 $b(\theta)=e^\theta$，适合计数观测。更一般地，广义低秩模型（generalized low-rank model, GLRM）允许不同条目使用不同损失，从而处理混合类型数据。

原始优化问题及其困难

在指数族模型下，忽略与 $\Theta$ 无关的常数，观测负对数似然为 $\ell_\Omega(\Theta)=\sum_{(i,j)\in\Omega}\{b(\Theta_{ij})-Y_{ij}\Theta_{ij}\}$。若直接要求低秩，最自然的估计问题是

$$\min_{\Theta}\ \ell_\Omega(\Theta) \qquad \text{s.t.}\qquad \operatorname{rank}(\Theta)\le r,\ \|\Theta\|_\infty\le \alpha.$$

其中 $\|\Theta\|_\infty\le\alpha$ 是常见的有界性约束。它不是形式上的附加条件，而是非高斯模型中非常重要的识别和曲率条件。例如在 logistic 模型中，若 $|\Theta_{ij}|$ 过大，则概率 $b'(\Theta_{ij})$ 接近 $0$ 或 $1$，似然曲率 $b''(\Theta_{ij})$ 接近 $0$，有限样本下很难稳定恢复自然参数。

该问题仍然是非凸的。困难来自两个方面：第一，秩约束本身导致组合型搜索，与第一篇中的秩最小化困难相同；第二，似然损失虽通常关于 $\Theta$ 凸，但与秩约束结合后不再是凸优化问题。因此，广义似然矩阵补全通常采用两条路线：核范数正则化的凸估计，或低秩因子化的非凸最大似然估计。

核范数正则化似然

核范数路线把秩约束替换为核范数惩罚，得到广义似然下的凸估计：

$$\min_{\Theta:\ \|\Theta\|_\infty\le\alpha} \ell_\Omega(\Theta)+\lambda\|\Theta\|_\ast.$$

当 $b$ 为凸函数时，$\ell_\Omega(\Theta)$ 关于 $\Theta$ 是凸函数；核范数 $\|\Theta\|_\ast$ 也是凸函数；无穷范数约束 $\|\Theta\|_\infty\le\alpha$ 也是凸约束。因此，这一问题是凸优化问题。高斯平方损失只是其中的特例：当 $b(\theta)=\theta^2/2$ 时，$\ell_\Omega(\Theta)$ 与 $\frac12\|\mathcal P_\Omega(\mathbf Y-\Theta)\|_F^2$ 只差一个常数。

理论分析通常需要以下条件：观测位置足够分散或采样分布可控；$\Theta^\ast$ 低秩且不应过度尖峰化；参数空间有界，使得 $b''(\theta)$ 在 $[-\alpha,\alpha]$ 上有上下界；正则化参数 $\lambda$ 能控制随机梯度 $\|\nabla\ell_\Omega(\Theta^\ast)\|_2$。在这些条件下，可以得到关于自然参数矩阵的 Frobenius 风险界，或关于分布的 KL 风险界。需要注意，非高斯模型中“恢复矩阵”不一定等价于逐项恢复观测值，而更常见的是恢复自然参数、均值或条件分布。

算法

一、凸松弛：核范数正则化最大似然

凸松弛方法直接求解 $\ell_\Omega(\Theta)+\lambda\|\Theta\|_\ast$ 及其带约束版本。其优势是目标函数凸，便于给出全局最优性和统计误差分析；主要代价是每一步通常涉及 SVD，规模很大时计算较重。

• 近端梯度法（proximal gradient）：梯度为 $\nabla\ell_\Omega(\Theta)=\mathcal P_\Omega(b'(\Theta)-\mathbf Y)$。若暂不考虑无穷范数约束，一步近端更新为 $\Theta^{t+1}=\mathcal D_{\eta\lambda}(\Theta^t-\eta\nabla\ell_\Omega(\Theta^t))$，其中 $\mathcal D_{\tau}$ 是奇异值软阈值算子。带 $\|\Theta\|_\infty\le\alpha$ 约束时，可结合投影、Dykstra 迭代或其他复合近端方法。
• ADMM / ALM：引入变量分裂 $\Theta=\mathbf Z$，把似然项和核范数项分开处理。$\Theta$ 子问题通常按条目分解，可用 Newton 或 IRLS 求解；$\mathbf Z$ 子问题是核范数近端步，即奇异值软阈值化；乘子更新用于逐步满足约束 $\Theta=\mathbf Z$。
• Frank-Wolfe / 条件梯度：当用核范数球约束替代核范数惩罚时，每步只需计算主奇异向量，而不必做完整 SVD。这类方法适合大规模问题，但收敛速度和精度控制需要单独分析。

二、非凸分解法：低秩广义线性矩阵分解

非凸路线直接设 $\Theta=\mathbf U\mathbf V^T$，其中 $\mathbf U\in\mathbb R^{n_1\times r}$ 且 $\mathbf V\in\mathbb R^{n_2\times r}$，并求解

$$\min_{\mathbf U,\mathbf V} \sum_{(i,j)\in\Omega}\{b(\mathbf u_i^T\mathbf v_j)-Y_{ij}\mathbf u_i^T\mathbf v_j\} +\frac{\lambda}{2}(\|\mathbf U\|_F^2+\|\mathbf V\|_F^2).$$

该问题关于 $(\mathbf U,\mathbf V)$ 非联合凸，但变量维度从 $n_1n_2$ 降到约 $r(n_1+n_2)$，适合大规模稀疏观测。

• 交替广义线性模型 / IRLS：固定 $\mathbf V$ 时，每一行 $\mathbf u_i$ 的更新是一个低维 GLM 问题；固定 $\mathbf U$ 时，每一列 $\mathbf v_j$ 同理。高斯情形下该算法退化为第一篇中的 ALS；logistic 或 Poisson 情形下通常用 Newton 或 IRLS 解子问题。
• 梯度下降与随机梯度下降：对每个观测 $(i,j)$，梯度残差为 $b'(\mathbf u_i^T\mathbf v_j)-Y_{ij}$。因此可用全量梯度、小批量梯度或在线更新处理超大规模推荐数据。实际实现中需要正则化、步长控制和参数截断，以避免因子范数不稳定。
• 黎曼梯度法：也可以直接在固定秩矩阵流形上优化 $\ell_\Omega(\Theta)$。其基本步骤是计算欧氏梯度 $\mathcal P_\Omega(b'(\Theta)-\mathbf Y)$，投影到当前秩 $r$ 流形的切空间，再通过截断 SVD 或其他 retraction 回到低秩集合。

概括地说，凸方法更接近统计理论中的标准估计量，非凸分解方法更接近大规模实践。广义似然设定下，算法与第一篇的主要区别不在低秩结构，而在数据拟合项从平方损失变成了由 $b$ 决定的负对数似然。

它和第一篇是什么关系

第一篇的平方损失模型可以看作本文的高斯特例。若 $Y_{ij}\mid\Theta^\ast_{ij}\sim N(\Theta^\ast_{ij},1)$，则 $b(\theta)=\theta^2/2$，负对数似然 $\ell_\Omega(\Theta)$ 等价于 $\frac12\|\mathcal P_\Omega(\mathbf Y-\Theta)\|_F^2$。因此，第一篇中的核范数最小化、SVT、ALM、ALS、黎曼梯度下降都可以看作本文算法在高斯似然下的特例。

广义似然模型的新增内容是：不同数据类型应使用不同的似然；估计对象往往是自然参数或均值，而不是直接把观测值当作连续数值；理论上必须控制似然曲率，尤其要避免 logistic、Poisson 等模型在极端参数区域出现退化。

几个容易混淆的点

• 低秩施加在哪个矩阵上必须说清楚。本文默认 $\Theta^\ast$ 低秩，而均值矩阵 $b'(\Theta^\ast)$ 一般不再严格低秩。
• 二元或计数矩阵不应简单套用平方损失。平方损失会隐含高斯误差，可能给出不合适的概率或方差结构。
• 1-bit matrix completion 估计的是生成二元观测的潜在概率或自然参数，不是从一个 bit 中恢复任意实数条目。
• 非高斯模型更依赖有界性和曲率条件。若 $\Theta^\ast$ 太大，似然可能接近饱和，矩阵即使低秩也难以稳定估计。
• 若缺失机制依赖于未观测值本身，即 MNAR，则仅靠上述似然仍不足够，需要进一步建模采样机制或选择偏差。

一句话总结

广义似然下的低秩矩阵补全，是把第一篇的平方损失替换为适合数据类型的负对数似然；在自然参数矩阵低秩、采样可控、参数有界且似然曲率良好的条件下，可以用核范数正则化最大似然或低秩因子化最大似然估计完整矩阵或其生成分布。

References

1. M. Collins, S. Dasgupta, and R. E. Schapire, A Generalization of Principal Components Analysis to the Exponential Family, NeurIPS, 2001.
2. N. Srebro and T. Jaakkola, Weighted Low-Rank Approximations, ICML, 2003.
3. N. Srebro, J. D. M. Rennie, and T. S. Jaakkola, Maximum-Margin Matrix Factorization, NeurIPS, 2004.
4. M. A. Davenport, Y. Plan, E. van den Berg, and M. Wootters, 1-Bit Matrix Completion, Information and Inference, 2014.
5. J. Lafond, Low Rank Matrix Completion with Exponential Family Noise, COLT / PMLR, 2015.
6. S. Gunasekar, P. Ravikumar, and J. Ghosh, Exponential Family Matrix Completion under Structural Constraints, arXiv, 2015.
7. M. Udell, C. Horn, R. Zadeh, and S. Boyd, Generalized Low Rank Models, Foundations and Trends in Machine Learning, 2016.