张量分解系列第二篇咱们把CP分解彻底聊透。上一篇文章讲了张量分解的整体框架理解了张量就是高维数组分解就是把多维数据拆成更紧凑的成分结构。这一篇的主角是CP分解全名CANDECOMP/PARAFAC一个双名字的分解方法CANDECOMP来自Carroll和Chang在1970年提出的“规范分解”思路PARAFAC则来自Harshman同年提出的“平行因子分析”。名字虽然拗口核心思想一句话就能讲完把一张高阶张量写成若干个“秩一张量”的和。这个直觉非常朴素就像把一个三维数据体切成一叠一层一层叠加的薄片每层都能单独解释。CP分解在脑电分析、化学光谱、推荐系统、交通数据补全里都极其常见也是所有张量分解算法里最值得先掌握的一个。这篇内容会比较硬核但是我会把每一步都尽量拆开讲CP分解的数学定义、核心运算、ALS求解原理再到用NumPy手写一个简化版、用Tensorly跑正式实验以及选秩、调参、排查问题的实战经验。适合刚接触张量分解的读者也适合已经会用SVD、想把矩阵方法推广到高维的算法工程师以及想在自己的数据处理流程里引入分解降维的研究人员。1. 从矩阵到张量CP分解解决的是哪一类问题1.1 为什么我们要拆张量先把场景想清楚。矩阵分解处理的是二维表格比如用户-物品评分矩阵做法是拆成两个低秩因子的乘积也就是常说的隐语义模型。但现实数据往往不只二维一个用户行为序列有用户、时间、地点三个维度一段脑电信号有电极、时间窗、频率三个维度一个路段交通数据有日期、时段、路口三个维度。如果你把这些高维数据压平成矩阵再去做矩阵分解就等于强行丢弃了维度之间的联合结构。举例来说交通数据压平成“时间×路口”之后日期信息就隐藏在哪一列或哪一行对应的编码里跨日期的周期性规律很难被矩阵分解单独提取。但如果我们保留三维张量可以直接分解出“路口模式”“时段模式”“星期模式”每个模式都有实际业务含义。这就是张量分解相比矩阵分解最本质的优势它不破坏原始维度结构拆出来的因子矩阵天然对应着不同模态的解释。1.2 CP分解的一页纸定义CP分解的数学形式非常简单。一个三阶张量X维度为I×J×KCP分解把它近似为R个秩一张量之和X ≈ Σ_{r1}^{R} a_r ∘ b_r ∘ c_r其中a_r是长度为I的向量b_r长度为Jc_r长度为K符号∘表示外积。所谓“秩一张量”就是三维数组可以写成三个向量外积的形式比如矩阵里的秩一矩阵可以写成两个向量的外积升级到三维就是三个向量的外积。把R个这样的成分加起来就得到完整的近似张量。这里的R就是“CP秩”。一个张量的CP秩定义为能精确表示它所需的最少秩一成分个数。注意这和矩阵秩有本质区别矩阵秩在实数域上有成熟算法CP秩的计算却是NP难问题实际使用中几乎没人去精确求CP秩而是把它当作一个超参数通过实验确定。这个后面会专门讲。1.3 适合看这篇文章的人这篇博文的定位是“理论够用、实操可复现”。我假设你已经有最基本的张量概念知道张量是一个多维数组会基本的矩阵运算、SVD概念不需要懂太多数学证明。如果你完全没接触过张量分解建议先把张量的mode-n展开和n-mode乘积了解清楚再看这篇会轻松很多。在阅读过程中只会用到几样工具NumPy、Tensorly、Python脚本。所有代码我都尽量给出完整可运行的版本。你最好配一台有8GB以上内存的电脑因为我们做合成数据实验时张量维度虽然不大但中间变量的展开矩阵可能会让你体会到什么叫内存膨胀。2. 数学基础三个运算一个都不能少2.1 张量展开mode-n matricizationCP分解的推导绕不开张量展开也就是把高维张量按某个维度重排成二维矩阵。三阶张量X ∈ R^{I×J×K}按第一个维度展开得到的矩阵X_{(1)}尺寸是I×(JK)意思是把每个固定第一维下标i的“切片”沿着另外两维拉平然后横向拼接。比如一个2×3×4的张量按第一维展开后是2×12的矩阵。展开的顺序有约定不同的库实现略有差异但原理一致固定目标维度下标剩下维度按固定次序枚举。Tensorly和NumPy的实现都遵循column-major的观念但我们在实际使用中只需要保证展开和运算自洽即可最终结果不受影响。理解张量展开是理解ALS求解的关键因为CP分解的每次更新都会把张量按某个模态展开进而转化成一个标准的最小二乘问题。2.2 Khatri-Rao积和Hadamard积CP分解的推导还会用到两个特殊矩阵运算如果你看过相关论文一定见过这两个符号。第一个叫Khatri-Rao积用符号⊙表示。它本质上是“列对应的Kronecker积”如果有矩阵A ∈ R^{I×R}和B ∈ R^{J×R}它们的Khatri-Rao积是(IJ)×R的矩阵其中每一列都是A的第r列和B的第r列做Kronecker积后拼接成长向量A⊙B [vec(b1 ⊗ a1), vec(b2 ⊗ a2), ..., vec(bR ⊗ aR)]注意这里要用克罗内克积的顺序。之所以出现这个运算是因为按张量展开后线性模型的形式会自然落到它身上。第二个叫Hadamard积用符号*表示就是两个同尺寸矩阵逐元素相乘。在CP分解的ALS更新里需要计算C^T C和B^T B然后做Hadamard积得到一个小尺寸R×R矩阵从而避免直接构造巨大的Khatri-Rao积矩阵。这是工程实现里一个非常重要的性能优化点。2.3 CP分解的标准形式把因子矩阵写出来理解更清晰。设A [a_1, ..., a_R] ∈ R^{I×R}B [b_1, ..., b_R] ∈ R^{J×R}C [c_1, ..., c_R] ∈ R^{K×R}。CP分解近似式等价于X ≈ Σ_{r1}^{R} a_r ∘ b_r ∘ c_r写成张量展开形式则是X_{(1)} ≈ A (C⊙B)^T。这个公式是ALS推导的起点。在实际库的实现里CP分解通常还会返回一个权重向量λ表示每个成分的整体幅度因子矩阵的列会被归一化成单位向量。李式表达为X ≈ Σ_{r1}^{R} λ_r (a_r/||a_r||) ∘ (b_r/||b_r||) ∘ (c_r/||c_r||)权重λ在推荐系统里可以理解为每个成分的重要性。这个归一化不是数学必需但对数值稳定性和结果可解释性很有帮助否则分解结果会被尺度差异淹没。3. CP分解的“个性”优势、局限与选型判断3.1 CP分解 vs Tucker分解很多初学者在第N次搞懂Tucker分解之后又开始纠结CP分解和Tucker到底有什么区别。这里我说一个最直接的理解角度Tucker分解是“张量的SVD推广”核心张量是一张完整的高阶小张量每个因子矩阵对应一个模态的主成分CP分解则可以看成Tucker分解的特例核心张量被强制约束为超对角张量也就是只有R个对角元素非零其余位置全为零。换句话说CP分解不允许不同成分在各模态之间“交叉混合”。一个成分在第一模态编号为1在第二模态编号为1这一组合有能量但编号组合(1,2)在CP里一定为零。Tucker分解则允许这种组合存在。所以CP分解得到的成分结构更稀疏、更符号化适合做“独立成分”式解释Tucker分解更灵活适合做压缩表示比如用较小核心张量保留主要信息然后去重构。如果你的目标是从数据里找“到底有哪些潜在模式”优先考虑CP如果目标是压缩存储或预训练特征提取Tucker通常效果更好。3.2 典型应用从脑电数据到化学光谱CP分解最有名的应用场景之一是脑电EEG数据分析。一次实验记录到的信号是电极×时间×频域的三维张量CP分解出的每个成分对应一个潜在神经过程某个电极分布模式、某个时间波形、某个频率谱。研究者不需要手动指定这些模式分解结果就能直接给出可解释成分。化学计量学领域CP分解还有个更经典的名字PARAFAC常用于荧光光谱数据分析。样品-发射波长-激发波长三阶张量分解出的每个成分对应一种荧光物质配合唯一性理论可以直接做定量分析。这套方法在环境、食品、生命科学仪器分析里用了四十多年稳定性非常高。推荐系统领域用户×物品×上下文的三阶张量CP分解可以让每个成分对应一个用户群体、一类物品、一种场景偏好。实际做推荐类项目时CP尤其适合数据稀疏的场景因为它的参数数量是O((IJK)R)而矩阵分解是O((IJ)R)张量分解不会因增加一个维度而导致参数爆炸。3.3 唯一性与退化文献里经常一笔带过CP分解有个让数学家非常兴奋但又容易让使用者踩坑的性质在特定条件下分解结果是“本质唯一”的。Kruskal提出过严格的条件跟因子矩阵的Kruskal秩有关。通俗理解只要每个因子矩阵的列足够独立且成分数不太大CP分解的解在最坏情况下也只在列交换和列缩放两种尺度上有歧义没有旋转问题。这和矩阵分解截然不同。矩阵SVD里虽然奇异值唯一但左奇异向量和右奇异向量可以乘以任意旋转矩阵而不改变结果CP分解在一定条件下没有这个自由度因此每个成分可以直接对应真实物理或业务意义。这是CP分解能用于“盲源分离”的理论基础。但坏消息是CP分解存在“退化”问题某些张量在试图逼近某个低秩近似时因子矩阵中某些列的范数会趋向无穷大但它们的线性组合却刚好抵消导致数值过程极不稳定。这在实际数据里宁可增加少量正则项比如在目标函数加一个λ(||A||^2||B||^2||C||^2)也不要放任ALS去裸奔。4. 求解CP分解ALS推导到代码实现4.1 目标函数与最优解推导求解CP分解最常见的方法是交替最小二乘ALS核心思路就一句话初始化A、B、C后每次固定其中两个矩阵更新第三个矩阵然后循环直到收敛。为什么能这么操作因为固定B和C后更新A对应的问题是一个凸的最小二乘问题有解析解。以三阶张量X为例展开表达式X_{(1)} ≈ A (C⊙B)^T。固定B和C时最小化||X_{(1)} - A (C⊙B)^T||_F^2对A求导令其为零得到A X_{(1)} (C⊙B) [ (C^T C) * (B^T B) ]^{-1}这个式子里核心矩阵(C⊙B)^T(C⊙B)恰好等于(C^T C) * (B^T B)大小是R×RR也就是成分数。这让计算变得非常轻量不需要真的构造(C⊙B)这种尺寸可能巨大的矩阵只要算两个R×R的小矩阵做Hadamard积。实际代码里会用线性方程组求解比如np.linalg.solve而不是显式求逆速度和稳定性都会更好。同时对R×R矩阵加一个很小的单位阵例如1e-10倍防止数值误差导致不可逆。4.2 用NumPy从零实现一个mini CP分解理论再漂亮不如自己敲一遍。这里我写了一个尽量简短但逻辑完整的NumPy实现先定义Khatri-Rao积、mode-1展开然后实现ALS更新。import numpy as np def khatri_rao(A, B): # A: I x R, B: J x R - (IJ) x R I, R A.shape J B.shape[0] C np.zeros((I * J, R)) for r in range(R): C[:, r] np.kron(A[:, r], B[:, r]) return C def mode1_unfold(X): # X: I x J x K - I x (J*K) I, J, K X.shape U np.zeros((I, J * K)) idx 0 for j in range(J): for k in range(K): U[:, idx] X[:, j, k] idx 1 return U def cp_als_mini(X, rank, max_iter200, tol1e-6): I, J, K X.shape rng np.random.default_rng(42) A rng.normal(size(I, rank)) B rng.normal(size(J, rank)) C rng.normal(size(K, rank)) # 单位列归一化 def normalize(M): scale np.linalg.norm(M, axis0) scale[scale 0] 1.0 return M / scale, scale A, _ normalize(A) B, _ normalize(B) C, _ normalize(C) X1 mode1_unfold(X) prev_err np.inf for it in range(max_iter): # 更新 A V (C.T C) * (B.T B) # (C⊙B)^T X1^T 形状 R x (I) M khatri_rao(C, B).T X1.T A np.linalg.solve(V 1e-10 * np.eye(rank), M).T # 这里未展开写 B、C 的更新实际需构造 X2、X3 # 完整实现请使用 Tensorly 或展开 mode2/mode3 # 计算重构误差 if it % 10 0: recon cp_to_tensor(A, B, C) err np.linalg.norm(X - recon) / np.linalg.norm(X) if abs(prev_err - err) tol: break prev_err err return A, B, C def cp_to_tensor(A, B, C): I A.shape[0] J B.shape[0] K C.shape[0] R A.shape[1] T np.zeros((I, J, K)) for r in range(R): T np.einsum(i,j,k-ijk, A[:, r], B[:, r], C[:, r]) return T这个mini实现里我特意只写了更新A的片段B和C的更新是一模一样的逻辑把张量展开到对应模式再套同一个最小二乘公式。真实的工程实现还要处理展开顺序、不同模态之间的变量循环、收敛判断等。写成教学代码没问题但真用于实验我更推荐直接用Tensorly。这里解释一下为什么更新B时要展开到mode-2。把X_{(2)}写成B (C⊙A)^T的近似形式所以更新公式为B X_{(2)} (C⊙A) [ (C^T C) * (A^T A) ]^{-1}更新C时类似。三个更新循环做下去就是完整的ALS。4.3 生产环境直接上Tensorly谈到正式实验Tensorly是Python生态里最顺手的张量库封装了CP、Tucker、张量回归、稀疏操作等一堆功能而且接口风格像sklearn学习成本低。安装很简单pip install tensorly使用CP分解的典型流程如下import numpy as np import tensorly as tl from tensorly.decomposition import parafac from tensorly import cp_to_tensor tl.set_backend(numpy) # 合成一个低秩张量方便验证 I, J, K, R 30, 40, 50, 4 A np.abs(np.random.randn(I, R)) B np.abs(np.random.randn(J, R)) C np.abs(np.random.randn(K, R)) X_true cp_to_tensor((np.ones(R), [A, B, C])) # 低秩张量 # 加一点噪声 X X_true 0.05 * np.random.randn(I, J, K) # CP分解返回 (weights, factors) weights, factors parafac(X, rankR, initsvd, tol1e-6, n_iter_max500) # 重构误差 X_hat cp_to_tensor((weights, factors)) err np.linalg.norm(X - X_hat) / np.linalg.norm(X) print(重构相对误差:, err) # 恢复的因子矩阵 A_est, B_est, C_est factors print(因子A形状:, A_est.shape)这里initsvd表示用张量按每个模态展开后的SVD来初始化因子矩阵比随机初始化稳定很多。Tucker初始化内部做了多模态SVD作为CP分解的起点效果也很好。每次调用parafac前我都建议固定随机种子这样实验结果可复现。4.4 初始化、收敛和超参数的经验ALS有一个老生常谈但实测影响极大的坑初始化。随机初始化可能导致不同轨道上的局部最优解质量差异巨大。我的习惯是先用SVD初始化然后结合多次随机初始化取重构误差最低的解。如果数据量能接受跑20个随机初始化一点都不浪费因为ALS的一次迭代成本很低。收敛判断也有讲究。常用的标准是相对重构误差变化量小于tol比如1e-6同时配合最大迭代次数兜底。但有些情况下误差下降很慢需要耐心跑几千次迭代尤其是成分数R偏大时。还有一点不要把收敛容忍度设得太紧否则可能陷在数值震荡里白白浪费算力。1e-5到1e-6已经够工程使用了。成分数R的选择更关键我见过不少人不管三七二十一直接取R10结果分解出一堆“混合成分”解释起来非常痛苦。后面第5.2节专门讲选秩方法。5. 实操案例验证算法、选秩、提取成分5.1 合成数据验证无论自己写ALS还是用库第一次跑都应该用合成数据验证“算法有没有写对”。方法很简单先手工生成一组已知的A、B、C用它们构造一个低秩张量再加少量噪声然后做CP分解看恢复出来的因子和原始因子之间能否对应上。由于CP分解存在列顺序和缩放歧义直接比较因子矩阵的数值意义不大。正确的比较方式有两种。一种是算因子矩阵列空间的重叠度比如对每列估计相关系数另一种是直接比较重构张量X_hat与X_true的误差。我一般两个都做。如果重构误差能压到噪声水平以下说明分解算法本身没问题再看恢复因子与真实因子的相关矩阵如果是接近置换矩阵的结构说明唯一性也保持了。我在合成实验里还经常做一个扰动测试把数据按行随机打乱再做分解。如果重构误差基本不变说明算法稳定性不错。相反如果误差波动巨大那可能是成分数选得不对或者初始化质量太差。5.2 如何选择CP秩CP秩选择没有像矩阵SVD那样一锤定音的“肘部法则”但最常用的刻度有两个。第一个是解释方差explained varianceEV 1 - ||X - X_hat||_F^2 / ||X||_F^2你从R1开始逐一测试每次增加R画出EV-R曲线。曲线变平时继续增加R往往只是开始拟合噪声这个拐点就是合理秩的候选。这个方法简单直接但病是“拐点”未必明显尤其在噪声大的数据上。第二个更接近张量领域原生的办法是CORCONDIA核心一致性诊断。原理上它检查计算出的低秩结构是否和“完美超对角核心”足够一致。Tensorly里可以这样粗略实现思路先做CP分解得到因子矩阵再估计核心张量然后计算一致度。一致性大于90%说明成分数合适低于50%则通常说明R取大了。对于工程实践我的建议是同时看EV曲线和CORCONDIA并配合业务可解释性综合判断不要纯靠数字说话。另外还有个简单粗暴但有效的校验法把数据随机分成两部分分别做同样R的CP分解看两组因子是否高度一致。R取太大会导致成分分裂、互相抵消结果不稳定R取太小时误差又太大稳定但拟合不足。5.3 真实场景案例时间-空间-频率张量的结构分解用一个具体案例把这些内容串起来。假设我们有若干个传感器记录了一段时间的振动信号每段信号取20个时间窗每个时间窗算32个频带的功率谱。于是数据是传感器×时间窗×频带的三阶张量尺寸是8×20×32。直接做CP分解设置R3得到三个成分。第一个成分在高频带上有明显峰值传感器权重集中在一号和三号传感器时间上在某个时段激活这很可能对应某个高频振动源。第二个成分低频宽带激活传感器分布均匀可能对应环境噪声。第三个成分中频窄带时间上与某个操作周期同步可能对应某台设备的周期性运转。这个例子想说明的是CP分解的意义不是给你一个漂亮的近似张量而是把一个混杂的多模态信号拆成几个独立可解释的物理过程。做完分解后别急着收工一定要结合业务去看每个成分如果某个成分从物理上说不通那可能是秩选大了或者数据预处理不够干净需要回头调整。6. 排查与避坑ALS不收敛、负成分、内存爆炸6.1 迭代曲线震荡或发散ALS最常见的临床表现是误差曲线震荡或者长期不下降。遇到这种情况排查顺序如下。先检查数据是否规范化。如果三个维度量纲差距过大比如第一维数值在1万量级第二维在0.01量级ALS求解最小二乘时会出现严重的条件数问题矩阵V接近奇异。我的习惯是分解前先做标准化比如每个模态按均值方差标准化或者至少把张量整体缩放到量级接近1。其次检查初始化。随机初始化的方差如果太大比如直接用np.random.randn因子矩阵列范数可能一开始就很大导致最小二乘解来回震荡。建议在初始化后马上对每个因子矩阵做列归一化并考虑用小方差随机初始化。再检查更新顺序和展开是否自洽这个主要影响自己写的实现。尤其是mode-1展开的列顺序必须和你构造Khatri-Rao积时的列顺序保持一致否则更新公式虽然形式上正确数值结果却完全是乱的。用Tensorly的话不会遇到这种问题但自己手写就容易踩。6.2 分解出负成分CP分解本身不保证因子非负。很多实际场景要求物理可解释比如光谱、医学影像、计数数据负的成分不仅没有意义还会污染后续分析。解决思路是加非负约束也就是非负CP分解NCP。Tensorly提供了相关功能底层用非负最小二乘替代普通最小二乘。如果不想引入太多依赖也可以在自己实现的ALS里每次更新因子后把负值截断为0然后用优化库做投影。不过直接截断不是一个严格收敛的方法更稳妥的做法是使用多层投影或者NNLS求解器。经验上如果数据本身物理上应为非负比如功率谱、计数、浓度直接跑无约束CP得到负成分的概率不低。所以这类数据一开始就应该走非负分解路线而不是等发现了负成分再补救。6.3 数据规模大怎么办CP分解在大规模数据上最大的瓶颈是内存mode-n展开矩阵的尺寸是I×(J×K)一旦每个维度都是几万展开矩阵直接就超出内存。工程上通常用稀疏张量格式配合特殊计算来解决。Tensorly支持sparse后端CP分解的核心计算可以只在非零元上进行但前提是数据确实稀疏比如推荐系统的“用户×物品×上下文”交互张量。如果数据是稠密的那我的建议是降维或者采样后再分解。先对每个模态做SVD降维把每个维度的有效尺寸压到几百甚至几十再做CP分解最后再映射回原始空间。这个方法虽然有点“绕”但在很多工程场景里是性价比最高的方案。别忘了CP分解的参数数量是O((IJK)R)主要开销其实在线性方程组的构建也就是计算各个模态的Gram矩阵这个开销只和维度相关因此降维后速度提升非常明显。6.4 结果校验三板斧调参调到怀疑人生时我会用这三个方法验证分解结果到底可不可信。第一重构误差曲线。给不同R分别画误差曲线如果某个R比前一个R误差下降非常有限那大概率是秩选大了。第二因子一致性。用不同随机种子跑多次看因子之间的相关系数矩阵是否稳定。稳定说明落入了一个有意义的局部最优不稳定说明分解被局部解带跑了。第三残差分析。把X - X_hat的残差画出来如果残差存在明显结构比如某块区域特别大说明当前秩或模型结构不足以刻画这个区域需要针对性修整。这三板斧合起来基本能覆盖“结果靠不靠谱”的判断。我相信多做几次这样的检查你也会形成属于自己的分解直觉。最后再分享一个心得遇到CP分解结果不理想先别急着怪算法先检查数据预处理。张量分解对尺度异常值非常敏感一个量纲失控的维度会让整个分解崩塌。先标准化、再分解、再做归一化因子反映尺度这个套路比任何高级调参都管用。搞明白这一点你手里的CP分解就超过大多数照搬文档的选手了。