线性因子模型：从生成模型原理到EM算法与代码实现

1. 为什么单独拎出线性因子模型来讲如果你翻过《深度学习》这本书的第十三章大概率会有一种“这章怎么突然变简单了”的错觉。前面第十二章还在讲各种复杂的正则化技巧第十四章马上要进入自编码器的世界夹在中间的线性因子模型看起来像是暴风雨前的宁静。但恰恰是这一章把生成模型最核心的几个思想用最朴素的方式讲透了。线性因子模型要解决的问题很具体你手头有一堆高维数据比如一张张图片、一段段音频每个样本的维度可能上千甚至上万但真正决定这些数据变化的“幕后推手”其实没几个。就像你看一个人每天穿什么衣服表面上有无数种搭配组合但真正起作用的变量可能只有季节、场合、心情这三四个。线性因子模型就是要把这几个隐藏的“推手”找出来并且假设它们对观测数据的影响是线性的。这个假设听起来很强但它的价值在于可解性。非线性模型虽然灵活但训练困难、解释性差而线性因子模型在数学上有一整套漂亮的闭式解或者简单的迭代算法。更重要的是它是理解PCA、独立成分分析、稀疏编码、甚至后面自编码器的基础。你如果跳过这一章直接去看自编码器会发现很多设计选择莫名其妙但回头补上这一课一切就顺了。这篇文章适合两类人一类是正在啃花书、被概率图模型绕晕的读者另一类是已经用过PCA或因子分析但没搞明白背后统一框架的从业者。我会从模型定义讲到参数估计再落到代码实现和踩坑经验尽量把书上的公式翻译成能跑的东西。2. 线性因子模型到底在建模什么2.1 从生成过程的角度理解模型线性因子模型的核心是一个生成故事。假设你有一个隐藏变量z它的维度是k通常远小于观测变量x的维度d。生成过程分两步先从某个先验分布里采样出z然后通过一个线性变换加上噪声得到x。用公式写就是x Wz b ε其中W是一个d×k的矩阵叫做因子载荷矩阵b是偏置项ε是噪声项通常假设为高斯噪声。这个式子看起来简单但它包含了线性因子模型的全部精髓观测数据的高维变化被分解成了少数几个因子的线性组合再加上每个维度独立的噪声。为什么这个分解有意义因为如果k选得合适W的每一列就代表一个“因子”对各个观测维度的影响模式。比如在人脸数据里第一个因子可能对应“光照方向”第二个因子对应“是否戴眼镜”第三个对应“性别”。这些因子是潜在的、不可直接观测的但它们通过W决定了你看到的那张脸长什么样。2.2 概率视角下的模型形式从概率角度写线性因子模型定义了一个联合分布p(z) N(z|0,I)p(x|z) N(x|Wzb,Ψ)这里Ψ通常是对角矩阵表示每个观测维度的噪声方差可以不同但维度之间噪声不相关。这个假设很关键它意味着给定z之后x的各个维度是条件独立的。如果你把Ψ限制成标量乘以单位矩阵那就得到了概率PCA如果进一步让噪声趋近于零就退化成普通PCA。边缘化z之后x的边际分布也是高斯p(x) N(x|b,WWᵀΨ)这个协方差矩阵WWᵀΨ是理解线性因子模型的关键。WWᵀ是低秩的秩最多为k而Ψ是对角的。所以模型在说观测数据之间的相关性全部由少数几个因子通过W来解释剩下的部分被认为是各维度独立的噪声。2.3 和PCA、因子分析的关系很多人学到这里会问这不就是PCA吗区别在于噪声假设。PCA完全不做概率假设它只是找一个低维线性子空间使得重构误差最小。而线性因子模型明确假设了噪声结构。当Ψ σ²I时最大似然解对应的子空间和PCA的主子空间一致但尺度不同。当Ψ允许每个维度不同时模型会更关注噪声小的维度这在某些数据上比PCA更合理。因子分析是线性因子模型的一个特例它通常假设z的先验是标准正态Ψ是对角矩阵。而独立成分分析则换了一个思路它不假设z是高斯而是假设z的各个分量相互独立且非高斯通过最大化非高斯性来估计W。稀疏编码又把先验换成了稀疏先验让z的大部分分量为零。所以你看线性因子模型是一个框架换不同的先验和噪声假设就得到不同的具体算法。花书第十三章的价值就在于把这个统一视角讲清楚了。3. 参数估计从最大似然到EM算法3.1 最大似然估计的困难给定数据矩阵X大小为n×d我们想估计W、b和Ψ。最直接的想法是最大化对数似然log p(X) Σᵢ log N(xᵢ |b,WWᵀΨ)但这个优化问题不好解因为W出现在协方差矩阵里而且W和Ψ耦合在一起。对W求导会得到一个复杂的矩阵方程没有闭式解。更麻烦的是模型存在旋转不确定性如果W是最优解那么WR对于任意正交矩阵R也是最优解因为WR(WR)ᵀWWᵀ。这意味着因子本身不可唯一识别只能识别它们张成的子空间。3.2 EM算法的推导思路EM算法是解决这个问题的标准工具。核心思想是如果z已知估计W就很简单如果W已知推断z也很简单。交替进行这两步就能收敛到局部最优。E步需要计算z的后验分布。给定x和当前参数后验是p(z|x) N(z|μ_z,Σ_z)其中Σ_z (IWᵀΨ⁻¹W)⁻¹μ_zΣ_z WᵀΨ⁻¹(x-b)这两个式子看着吓人但实际计算时因为Ψ是对角矩阵求逆就是逐元素取倒数复杂度可控。Σ_z对所有样本都一样只需要算一次。M步更新参数W_new [Σᵢ (xᵢ-b) E[zᵢ]ᵀ] [Σᵢ E[zᵢzᵢᵀ]]⁻¹b_new (1/n) Σᵢ (xᵢ-WE[zᵢ])Ψ_new diag((1/n) Σᵢ [ (xᵢ-b)(xᵢ-b)ᵀ -WE[zᵢ](xᵢ-b)ᵀ ])这里的 E[zᵢ] 就是μ_zE[zᵢzᵢᵀ] 是Σ_zμ_z μ_zᵀ。实际写代码时把这些期望用E步的结果代入就行。3.3 实际操作中的参数选择k怎么选这是线性因子模型最实际的问题。理论上可以用交叉验证但在无监督场景下没有标签交叉验证不好做。常见的做法是看重构误差随k的变化曲线找拐点。或者看Ψ的估计值如果某些维度的噪声方差已经很小再增加k意义不大。另一个经验是如果数据维度d很大但样本量n不大k不要超过n/2否则容易过拟合。我试过在n500、d2000 的数据上把k设到100结果EM算法收敛后W的很多列几乎为零说明模型容量过剩了。还有一个坑是初始化。EM算法对初始值敏感W初始化太大会导致Σ_z接近奇异数值不稳定。通常用小的随机高斯初始化WΨ初始化为数据每个维度的方差。这样起步比较稳。4. 代码实现从零写一个线性因子模型4.1 数据准备和预处理我用一个模拟数据集来演示这样你能清楚看到因子恢复的效果。生成方式设定真实的W_true为d×kz_true为n×k的标准正态然后Xz_true W_trueᵀ 噪声。import numpy as np np.random.seed(42) n, d, k 500, 50, 5 W_true np.random.randn(d, k) * 0.5 Z_true np.random.randn(n, k) noise np.random.randn(n, d) * 0.3 X Z_true W_true.T noise预处理很重要。线性因子模型假设数据是零均值的所以要先减去均值。如果各维度量纲差异大还要标准化。但注意标准化会改变噪声结构如果原始数据的噪声方差本来就不一样标准化可能破坏模型假设。我的做法是只做中心化不做标准化除非确实有必要。X_mean X.mean(axis0) X_centered X - X_mean4.2 EM算法主循环下面是一个完整的EM实现。关键点在于E步的数值稳定性计算Σ_z时用np.linalg.solve而不是直接求逆避免数值误差。def linear_factor_model_em(X, k, max_iter200, tol1e-6): n, d X.shape # 初始化 W np.random.randn(d, k) * 0.1 b X.mean(axis0) Psi np.var(X, axis0) log_likelihood_old -np.inf for iteration in range(max_iter): # E步 Psi_inv 1.0 / Psi # 计算 Sigma_z (I W^T Psi^{-1} W)^{-1} M np.eye(k) (W.T * Psi_inv) W Sigma_z np.linalg.inv(M) # 计算 mu_z for all samples X_centered X - b # mu_z Sigma_z W^T Psi^{-1} (x - b) mu_z (X_centered * Psi_inv) W Sigma_z # 计算 E[z z^T] Sigma_z mu_z mu_z^T E_zz Sigma_z mu_z.T mu_z / n # M步 # W_new (X_centered^T mu_z) (E_zz)^{-1} W_new (X_centered.T mu_z) np.linalg.inv(E_zz * n) # b_new b_new X.mean(axis0) - W_new mu_z.mean(axis0) # Psi_new X_recon mu_z W_new.T b_new residual X - X_recon Psi_new np.mean(residual**2, axis0) Psi_new np.maximum(Psi_new, 1e-6) # 防止为零 # 计算对数似然 cov W_new W_new.T np.diag(Psi_new) sign, logdet np.linalg.slogdet(cov) log_likelihood -0.5 * (n * (d * np.log(2*np.pi) logdet) np.sum((X_centered) np.linalg.inv(cov) * X_centered)) # 检查收敛 if abs(log_likelihood - log_likelihood_old) tol: break W, b, Psi W_new, b_new, Psi_new log_likelihood_old log_likelihood return W, b, Psi, mu_z这段代码里有个细节计算对数似然时用了np.linalg.inv(cov)当d很大时这会很慢。实际应用中可以用Cholesky分解来加速但作为教学示例这样写更直观。4.3 结果验证和因子旋转跑完EM之后W的列和真实的W_true的列不会一一对应因为存在旋转不确定性。要比较它们需要做Procrustes对齐找一个正交矩阵R使得W R最接近W_true。from scipy.linalg import orthogonal_procrustes R, _ orthogonal_procrustes(W, W_true) W_aligned W R对齐之后可以计算因子之间的相关性看看恢复得怎么样。在我的模拟里前三个因子通常恢复得很好相关系数在0.9以上后两个因子因为对应的奇异值小容易被噪声淹没相关系数可能只有0.6左右。这提示我们如果真实因子的影响强度差异大弱因子很难被准确估计。另一个验证方式是看重构误差。用估计的W和mu_z重构X计算均方误差。如果k选得合适重构误差应该接近真实噪声方差。如果远大于噪声方差说明k太小如果远小于说明过拟合了。5. 常见问题与排查技巧实录5.1 EM算法不收敛或收敛到奇怪的结果这是最常见的问题。表现是迭代很多次对数似然还在震荡或者Ψ的某些元素趋近于零导致数值爆炸。原因通常有三个初始化不好、k太大、数据没有中心化。排查步骤先检查数据均值是否接近零如果不是手动减去均值。然后看Ψ的最小值如果小于1e-8说明模型在试图用因子解释噪声这时候要么减小k要么给Ψ加一个下限。我通常设下限为数据方差中位数的1e-3倍这样既不会限制模型又能保证数值稳定。还有一个隐蔽的坑如果数据中存在常数维度方差为零Ψ的初始值就是零E步计算Psi_inv会直接报错。处理方法是先删掉这些维度或者在初始化时给Ψ加一个很小的正数。5.2 因子个数k怎么定除了前面说的看重构误差拐点还有一个更统计的做法计算每个因子对应的特征值看累计方差贡献率。但注意线性因子模型的特征值和PCA的特征值不完全一样因为噪声结构不同。我的经验是先用PCA快速看一下特征值谱如果前几个特征值明显大于后面的k就取那个拐点位置。然后再用线性因子模型验证如果估计出的Ψ在各个维度上比较均匀说明k合适如果某些维度的Ψ特别小说明这些维度几乎被因子完全解释了可能需要增加k。下面这个表总结了几种k选择方法的适用场景方法适用场景优点缺点重构误差拐点数据量大、噪声均匀直观、计算快拐点有时不明显特征值累计贡献快速初筛有现成工具忽略噪声结构交叉验证似然有验证集统计严谨无监督场景难划分噪声方差均匀性精细调参直接反映模型拟合需要多次训练5.3 因子解释性差怎么办线性因子模型估计出的W列向量很多时候看起来像噪声没有明显的物理意义。这通常是因为旋转不确定性任何正交旋转都给出相同的似然但旋转后的因子可能更难解释。解决办法是加一个旋转准则比如Varimax旋转让每个因子只在少数几个维度上有大载荷这样更容易解释。from scipy.stats import ortho_group def varimax(W, gamma1.0, max_iter100, tol1e-6): d, k W.shape R np.eye(k) dW W.copy() for _ in range(max_iter): R_old R.copy() # 计算旋转后的载荷 L dW R # Varimax准则的梯度 u, s, vh np.linalg.svd( dW.T (L**3 - (gamma/d) * L np.diag(np.sum(L**2, axis0))) ) R u vh if np.max(np.abs(R - R_old)) tol: break return dW R, RVarimax旋转之后W的每一列会趋向于稀疏只有少数几个维度上的值比较大。这时候再去看哪些维度载荷大往往能对应到数据中可解释的模式。我在一个用户行为数据上试过旋转前因子看不出什么意思旋转后第一个因子明显对应“高频使用时段”第二个对应“功能偏好”解释性提升很大。5.4 和PCA结果对不上有人会发现自己写的线性因子模型跑出来的子空间和PCA不一样。这很正常因为当Ψ不是标量矩阵时线性因子模型的极大似然解本来就不等于PCA。PCA最小化的是重构误差而线性因子模型最大化的是似然两者在噪声结构不均匀时会有差异。如果你希望结果接近PCA可以把Ψ限制为所有维度相同即Ψ σ²I然后优化 σ²。这样得到的解和PCA的主子空间一致只是尺度不同。6. 从线性因子模型到更复杂的模型线性因子模型虽然简单但它是很多现代生成模型的起点。理解了这个框架再看变分自编码器就会觉得顺理成章VAE无非是把线性映射Wz换成了神经网络把先验和近似后验换成了更灵活的形式但EM算法的影子还在只是E步变成了神经网络的编码器M步变成了解码器的训练。稀疏编码是另一个方向把z的先验从高斯换成拉普拉斯就得到了稀疏的因子。独立成分分析则是假设z非高斯通过最大化非高斯性来估计W。这些变体在信号处理、神经科学、金融数据分析里都有大量应用。如果你想把线性因子模型用到实际项目里我的建议是先跑PCA看数据有没有低维结构如果有再用线性因子模型精细建模噪声。如果数据量很大EM算法可能慢可以考虑随机优化或者用PCA初始化再微调。另外如果数据明显是非线性的别硬套线性因子模型直接上自编码器或者核方法。最后分享一个我在实际项目中总结的小技巧当d远大于n时直接对d×d的协方差矩阵做特征分解很慢可以先对数据矩阵做SVD用右奇异向量来初始化W这样能大大加速EM的收敛。具体做法是取X的前k个右奇异向量乘以对应的奇异值作为W的初始值。这个初始化比随机初始化稳定得多通常迭代十几次就能收敛。