简介矩阵分解中的SVD是处理二维数据的经典工具但当数据扩展到三维乃至更高维时如图像的通道、时间序列的批次传统的矩阵方法便难以刻画多个维度间的耦合结构。张量分解作为矩阵分解的高维推广通过多个因子矩阵与核心张量同时压缩各维度信息能够保留数据内在的多方向结构。HOSVD作为Tucker分解的一种高效实现路径分别对每个模式展开矩阵做SVD再重构核心张量计算简单且逼近误差有理论保证被广泛用于图像压缩、信号处理、推荐系统与科学计算等领域。本文从三阶张量的基本概念出发详解HOSVD的数学原理、完整NumPy代码实现、多线性秩的选择策略以及常见工程陷阱帮助读者从原理到实践系统掌握这一高维数据分析利器。 我们拿到一批三维数组模样的数据时第一直觉往往会去试PCA、SVD这类矩阵工具但一旦数据本身是“立体的”比如一张灰度图像天然就是行乘列的两维再加上通道、时间、批次这些维度SVD就撑不住了。不是不能压平了算而是强行把三维数组拉成矩阵等于把原本沿不同方向变化的规律全部搅在一起丢失了张量作为整体时那种张成的结构信息。这正是张量分解、尤其是HOSVD和Tucker分解派上用场的场景。这篇东西想聊的就是HOSVD高阶奇异值分解在三阶张量上的完整落地过程从数学原理到NumPy手写实现再到多线性秩怎么选、奇异值怎么看、哪些坑我踩过。适合理工科背景、准备把张量分解用到信号处理、图像压缩、推荐系统或科学计算里的朋友。如果你只是想调包跑个结果可以直接跳到第三节如果你想搞清楚Tucker分解和HOSVD到底什么关系、为什么HOSVD有“准最优”的说法建议从头看起。1. 三阶张量分解到底在解决什么问题1.1 从矩阵SVD不够用说起先想一个最简单的场景你有100个用户每人看了50部电影评分构成一个100乘50的矩阵。SVD能做的是把这个矩阵分解成用户因子矩阵、奇异值对角阵和电影因子矩阵从而把用户和电影各自嵌到一个低维空间里。这个逻辑成立的前提是你手里的数据确实可以被“两个方向”描述干净。但现实里多数数据没那么乖。同样是评分数据你可能还想加入时间维度比如每个用户在不同月份对电影的评分那数据就从矩阵变成三阶张量用户乘电影乘时间。此时如果仍然压平成矩阵就有两种做法一是把所有时间片叠在一起变成一个大矩阵但这等于默认时间维度和用户维度在同一个语义层面显然不合理二是对每个时间片单独跑SVD但这样不同时间片的分解结果又对不上号没法对齐比较。无论是哪种都把三个维度之间的耦合关系破坏掉了。张量分解的核心思路是既然数据是一个三阶对象那就用三个因子矩阵或者叫因子向量组从三个方向同时挤压它再配一个“核心张量”去描述维度之间的交互。这跟矩阵SVD把A拆成U乘S乘V^T是同一个逻辑的不同版本只不过原来矩阵U和V各自只有一个方向现在Tucker/HOSVD里有U1、U2、U3三个方向核心张量也不要求是对角阵。1.2 Tucker分解与HOSVD的定位差异很多人会把HOSVD、Tucker分解混为一谈在不少论文里也确实没做严格区分但实际落地时两者差异还挺重要。Tucker分解的核心公式长这样X ≈ G ×₁ A₁ ×₂ A₂ ×₃ A₃其中G是核心张量尺寸是r1×r2×r3A1、A2、A3分别是三个方向的因子矩阵也可以理解成三个方向上的“主成分方向”。这个形式表达的是原始大张量可以被压缩成一个小核心张量外加三组低维投影方向。Tucker分解本身不限定怎么求A1、A2、A3甚至不限定G是什么它只是描述了一种“多方向线性子空间逼近”的框架。HOSVD则是求解Tucker分解的一种具体算法而且名字里带“SVD”也说明它是从矩阵SVD推广过来的——分别在三个mode上做矩阵化再做普通SVD取左奇异向量作为因子矩阵最后用因子矩阵反推核心张量。所以在数学结构上HOSVD输出的结果天然就是一个Tucker分解但反过来Tucker分解的求解路径并不一定是HOSVD还有更高阶的优化算法如HOOI。这里值得一提的是HOSVD往往被叫做“Tucker分解的准最优解”。原因是Tucker的最佳低秩逼近问题在数学上是个非凸优化问题全局最优极难求而HOSVD走的路线是“分别对每个方向做截断SVD再拼接”每一步都有显式解所以计算上非常快在小规模问题上效果也很接近最优。你可以把HOSVD理解成一种贪心策略而HOOI这类算法是在HOSVD给出的基础上再迭代精调。1.3 三阶张量的本质模式、纤维与切片要真正理解HOSVD还要把三阶张量的几个基础概念过一遍不然代码里的mode-1、mode-2、mode-3很容易把人绕晕。三阶张量X的尺寸如果是I×J×K那么沿着第一维row方向固定索引i得到的是J×K的矩阵这叫第i个水平切片沿着第二维固定索引j得到I×K的矩阵沿着第三维固定索引k得到I×J的矩阵。切片对应的是“在某个维度上切一刀”是观察张量的窗口。更底层的是纤维fiber概念。固定其他维度索引只让一个维度变化得到的一维数组就是纤维。在矩阵里行和列就是两种纤维在三阶张量里纤维有三种列纤维沿第一维、行纤维沿第二维、管纤维沿第三维。SVD中对矩阵做分解本质上是在找两种纤维共享的低维结构而HOSVD就是把它推广成三种纤维的共享子空间逼近。如果把三阶张量X沿第一维展开成矩阵X(1)尺寸是I×(JK)每一列就是一条“沿第一维变化的纤维”被摊平后的结果。HOSVD中的第一步——对每个mode做矩阵展开并做SVD——其实就是在问在这个方向上、把张量看作一大撮纤维的集合时它们主要沿着哪些正交方向分布这样一想HOSVD的过程就非常直观先在三个方向上分别做PCA再把三个方向的主成分坐标拼成核心张量。2. HOSVD的核心计算路径与原理拆解2.1 三个关键步骤展开、SVD、重构核心张量HOSVD的完整算法只有三步而且每一步都有清晰的线性代数背景。第一步mode-k展开。对三阶张量X做mode-1展开得到矩阵X(1)尺寸为I×(JK)展开规则其实就是在代码里用reshape加permute实现。展开不是数据搬家而是把张量按特定顺序重排让“沿mode-1变化的元素”排成一列这样做的目的是让SVD能吃到一个标准二维矩阵。mode-2和mode-3同理只是排列顺序不同。第二步对每个展开矩阵做截断SVD。对X(1)做SVD得到左奇异向量矩阵U1尺寸I×I按奇异值从大到小排序后取前r1列作为因子矩阵A1。对X(2)和X(3)做同样的操作分别得到A2和A3。注意这里SVD不用算V因为我们只需要左奇异向量而且截断后V是在后续重构核心张量时才隐式参与的。第三步求解核心张量。有了A1、A2、A3之后核心张量的估计公式是G X ×₁ A₁ᵀ ×₂ A₂ᵀ ×₃ A₃ᵀ这个操作在物理意义上可以理解为把原始张量沿着三个方向分别投影到低维子空间里投影后得到的“坐标张量”就是核心张量G。如果你对矩阵的SVD熟悉这会让你想到X UΣVᵀ里的Σ就是原始矩阵在U和V方向上的投影坐标只是张量情形下坐标不是对角阵而是一个实心的小张量。这也是Tucker分解和矩阵SVD最直观的差异矩阵SVD的“核心”是对角的Tucker分解的“核心”是稠密的。2.2 为什么HOSVD是“准最优”而不是“最优”这里需要多解释几句因为很多人在对比HOSVD和HOOI时会把准最优理解成“HOSVD逼近效果不如HOOI”这话大体对但容易忽略一个前提HOSVD的出发点根本就不是求最优逼近。我们把HOSVD的数学本质拆开看它对X(1)、X(2)、X(3)分别做SVD这其实是在三个不同的矩阵空间中求三个正交子空间。这组子空间本身的性质很好——每个子空间都是对应展开矩阵的前r个主左奇异方向因此各自都是“最佳r维子空间”。但是三个子空间各自最佳并不代表它们张成的张量积空间是最佳的三维子空间组合。打个比方三个人各自选了对自己最优的出行路线但三个人拼成一辆车之后整体的最优路线并不等于这三条路线的笛卡尔积。这是HOSVD准最优的根本原因。不过HOSVD的误差其实有理论保证。如果原始张量X的Tucker秩是(r1,r2,r3)的“有效秩”也就是说忽略掉小奇异值分量后残余很小那么HOSVD截断后的逼近误差和理论最优之间的差距不超过最优误差的√3倍在三阶情形下更精确地说差距有界于√r1√r2√r3的量级。这个结果保证了HOSVD不是瞎猜方向而是有界的近似最优。如果追求真正的最优Tucker分解通常用的是HOOI高阶正交迭代。HOOI的思路是把Tucker分解当成一个交替优化问题固定A2、A3优化A1固定A1、A3优化A2不断迭代直到目标函数收敛。每轮迭代的核心操作其实就是在利用当前的A2、A3把张量压缩成一个小张量再对这个压缩结果做SVD提取新的A1。HOOI每一步都在提升拟合度所以理论上能找到更好的局部最优代价是慢。实践中HOSVD经常被用来给HOOI提供初始值这个组合拳几乎能对标全局最优。2.3 HOSVD、Tucker、高阶PCA之间的关系我在看论文和做实现时经常发现这些术语在不同文章里指的东西略有差别这里整理一下方便你避坑HOSVD是广义叫法有的文献也把高阶奇异值分解称为Tucker分解的HOSVD形式特指“通过逐mode SVD求Tucker分解”的算法过程。Tucker分解强调的是分解形式本身即三个因子矩阵加一个核心张量的结构不限定算法。任何能求出这个结构的算法都可以称为在做Tucker分解。higher-order SVD和Tucker-1有时是同义词在早期文献里Tucker-1指只压缩一个方向的分解Tucker-2是压缩两个方向Tucker-3是三个方向都压缩。HOSVD通常等价于Tucker-3这也是为什么大家默认HOSVD就是“全方向的Tucker”。还有一个容易混淆的概念是CP分解CANDECOMP/PARAFAC它的形式是把张量分解成若干个秩一张量之和X ≈ Σ λᵣ aᵣ ∘ bᵣ ∘ cᵣ。CP和Tucker的最大区别是核心张量是否对角——CP的核心张量强制对角Tucker不强制。因此CP的参数更少、解释更强硬但求解稳定性更差对噪声敏感Tucker则灵活得多。一句话概括HOSVD是求Tucker分解的一条具体路径而Tucker分解是你要的那个结构本身。3. 手写NumPy实现HOSVD从零搭建完整流程3.1 环境与基础工具接下来上代码我这边的实现环境是Python 3.10 NumPy 1.24不需要额外装Tensorly或PyTorchHOSVD最核心的数学操作其实就两个张量展开reshape transpose和SVDnumpy.linalg.svd。这两个操作对中小规模张量完全够用上了大规模再考虑稀疏分解或随机化SVD。先定义一个三阶张量的mode-k展开函数。这里要注意NumPy里没有原生的“张量展开”函数需要自己写而且很容易在维度顺序上出错。一个稳妥的做法是先用np.moveaxis把目标维度挪到最前再做reshape这样逻辑最清晰不容易错import numpy as np def mode_k_unfold(X, k): 将三阶张量X按第k个模式展开为矩阵 k0,1,2 对应第一、第二、第三维度 返回形状为 (X.shape[k], 其它维度乘积) 的二维数组 # 把第k维移到最前面其他地方保持原顺序 moved np.moveaxis(X, k, 0) shape moved.shape # 第一维是第k维后面所有维度合并 return moved.reshape(shape[0], -1)这个函数后面所有步骤都要复用所以花点时间确认逻辑是对的X被moveaxis后第0轴就是我们要的mode方向剩下的轴原封不动按顺序排。reshape成(shape[0], -1)后每一列恰好对应一条沿mode方向变化的纤维。3.2 HOSVD主函数实现接下来实现HOSVD主函数。核心流程是三步走三个方向展开、三个SVD取左奇异向量、投影求核心张量。我写的版本要求用户传入截断秩列表ranks [r1, r2, r3]如果你不确定怎么选秩可以先跑一遍完整SVD看奇异值分布再回来填这个参数。def hosvd(X, ranks): 对三阶张量X执行HOSVD截断版本 X: numpy.ndarray形状 (I, J, K) ranks: list或tuple[r1, r2, r3] 返回: G: 核心张量形状 (r1, r2, r3) factors: 因子矩阵列表 [U1, U2, U3] ndim X.ndim if ndim ! 3: raise ValueError(本实现仅支持三阶张量) factors [] for k in range(ndim): # 第k模式展开 X_k mode_k_unfold(X, k) # 对展开矩阵做SVD只需要左奇异向量 U, _, _ np.linalg.svd(X_k, full_matricesFalse) # 截断到前r个主方向 r ranks[k] U_r U[:, :r] factors.append(U_r) # 核心张量X ×1 U1^T ×2 U2^T ×3 U3^T G X for k in range(ndim): G np.tensordot(factors[k].T, G, axes([1], [k])) # tensordot之后维度顺序可能变化需要调整回标准顺序 # 这里采用连续tensordot后手动确认shape # 一个更稳妥的写法是循环moveaxis但上面已经够用验证下 # 实际运行中每次tensordot会消掉第k个维度并把它放到最前 # 由于我们循环k0,1,2最终G的维度顺序是 (r1, r2, r3) return G, factors这段代码有几个细节值得注意。第一np.linalg.svd默认返回的U、S、Vh三个量对展开矩阵X_k我们只关心US和Vh都可以直接丢弃。第二我用了full_matricesFalse这样当展开矩阵的列数小于行数时U的列数自动与S的长度对齐不会出现多余的全零奇异值列截断时更省心。第三tensordot的axes传参容易写错我这里在注释里已经说明第一次tensordot会消掉X的第0维把因子矩阵的转置乘上去得到的新张量维度顺序为(r1, J, K)第二次tensordot消掉第1维得到(r1, r2, K)第三次消掉第2维得到(r1, r2, r3)。由于我们恰好是按顺序消的维度顺序不用额外调整但如果你在循环里跳着消一定要记得最后修正维度顺序。3.3 从HOSVD到Tucker重构的验证光有分解还不够还得能验证重构误差。Tucker/HOSVD的重构公式是X_approx G ×₁ U1 ×₂ U2 ×₃ U3也就是把核心张量沿着三个方向用因子矩阵映射回原始空间。用一个完整的测试脚本来验证整个流程# 构造一个随机三阶张量并加入一点低秩结构便于观察 np.random.seed(42) I, J, K 30, 20, 25 R1, R2, R3 5, 4, 6 # 先构造一个低秩张量A1(30x5), A2(20x4), A3(25x6)核心张量(5,4,6) A1 np.random.randn(I, R1) A2 np.random.randn(J, R2) A3 np.random.randn(K, R3) G_true np.random.randn(R1, R2, R3) X np.einsum(ir,js,kt,rst-ijk, A1, A2, A3, G_true) X 0.05 * np.random.randn(I, J, K) # 加一点噪声 # 执行HOSVD截断秩设为真实秩 G_hat, factors hosvd(X, [R1, R2, R3]) U1, U2, U3 factors # 重构 X_approx np.einsum(rst,ir,js,kt-ijk, G_hat, U1, U2, U3) # 计算相对误差 rel_err np.linalg.norm(X - X_approx) / np.linalg.norm(X) print(f相对重构误差: {rel_err:.6f})跑完这段代码你会得到一个相对误差约0.00x量级的结果具体数值取决于噪声水平。这个验证过程非常关键如果结果误差很大大概率是维度顺序写错了而不是算法本身的问题。我在调通第一版实现时就因为tensordot的维度顺序搞错误差一直是0.7左右排查很久才发现是核心张量重构时维度错位了。3.4 用Tensorly对比验证实现正确性如果你还是担心自己的实现有bug可以拿Tensorly这个库来交叉验证。Tensorly的高阶SVD实现很成熟在学术界和工业界都被广泛使用用它作为基准对比是最快的办法。import tensorly as tl from tensorly.decomposition import tucker # Tensorly的tucker直接支持指定秩 X_tensor tl.tensor(X) core_tl, factors_tl tucker(X_tensor, rank[R1, R2, R3], initsvd, n_iter_max0) # initsvd n_iter_max0 意味着只做HOSVD初始化不做HOOI迭代注意上面这段里initsvd且n_iter_max0意思就是“只用HOSVD结果不继续迭代优化”这正好对应我们手写的HOSVD。对比一下手写版本的因子矩阵与Tensorly输出的因子矩阵会发现它们的方向完全一致全局符号可能相反不影响语义核心张量也几乎相同。符号相反这件事我们放到第五节细聊它本身是个典型的坑。4. 多线性秩怎么选奇异值谱怎么读4.1 三种方向、三种奇异值谱HOSVD对你透露信息最丰富的时刻其实发生在选择秩之前的那个完整SVD计算。对三个模式展开矩阵分别做SVD后你会得到三组奇异值s1长度min(I, JK)、s2长度min(J, IK)、s3长度min(K, IJ)。这三组奇异值刻画的是张量沿三个方向各自的“能量分布”。矩阵SVD里大家习惯用奇异值累积能量占比来决定主成分个数这个习惯在张量情形下依然有效只是要分别对三个方向做判断。比如s1的前几个值迅速衰减、后面基本平了那就说明第一维方向上有效子空间维度很低s2衰减慢说明第二维方向信息冗余较少需要更多秩去刻画。画图时我建议把三个谱画在同一张图上横轴都是索引纵轴是奇异值或者归一化到最大奇异值。看图时重点观察两个地方第一个是“拐点”在哪里也就是奇异值从大变小趋势明显变缓的位置这个点往往是有效秩的候选第二个是一开始的几个值是不是占了绝对主导如果前两个奇异值占比超过90%那么这个方向取r2基本就够了。4.2 常用选秩方法对比选秩没有统一答案不同任务有不同的考虑方向。我把常见的几种方法连同适用场景一起列出来方法思路适用场景注意事项奇异值累积能量保留总能量90%/95%对应的前r个奇异值压缩、降噪需要对三组谱分别判断取max或逐方向取拐点/差值最大找奇异值二阶差分最大的位置低噪声、结构明显的数据对噪声敏感适合肉眼辅助判断交叉验证把张量元素随机屏蔽一部分用不同秩重构选误差最小的秩组合缺失值填充、推荐系统计算量大要小心随机屏蔽对空间结构的破坏任务驱动不追求最优重构而是按下游任务效果选秩分类、回归、聚类和分解目标分离容易过拟合需要控制验证集贝叶斯/自动秩选择用概率模型估计秩后验分布复杂噪声、分布式数据实现复杂常规任务不需要上这种手段实操中我常用的办法是先用能量占比确定每个方向秩的粗范围再在这个小范围内做交叉验证。因为三阶张量的秩组合空间往往是r1乘r2乘r3的笛卡尔积穷举很不现实但如果你把三个方向的秩分别限制在2到3个候选值里组合数也就十几个跑一轮交叉验证完全可行。4.3 压缩率与重构误差的权衡实操选秩本质上是在回答“我打算牺牲多少信息换取多大压缩”这个问题。为了量化这个权衡我习惯在选定秩之前先把“压缩率-重构误差曲线”画出来。压缩率的定义可以这样算原始张量需要存储I×J×K个floatTucker分解后需要存储r1×r2×r3核心张量加上(I×r1 J×r2 K×r3)因子矩阵所以压缩率是ratio (r1×r2×r3 I×r1 J×r2 K×r3) / (I×J×K)比如I100, J100, K100秩取(10,10,10)那么存储量是1000 1000 1000 1000 4000原始是1000000压缩率高达250倍。如果取(50,50,50)存储量就是125000 5000 5000 5000 140000压缩率约7.1倍。可见秩对压缩率的影响是非线性的尤其三个方向同时增大时核心张量的体积按三次方增长。有了压缩率的概念后我的选秩流程就很机械先看奇异值谱得到候选秩集合然后对每个候选秩组合算一遍重构误差和压缩率画成二维散点或者平行坐标图最后根据业务对“压缩率优先”还是“误差优先”的偏好定夺。在图像压缩任务里往往压到20倍以内肉眼几乎无差异这时候多花一倍的存储换0.1%的精度改善就没必要了。5. 实操中绕不开的坑与排查技巧5.1 奇异值分解的符号不稳定与左奇异向量漂移我在最初跑HOSVD时踩的第一个坑就是“两次运行结果看起来不一致”。后来发现HOSVD对因子矩阵做了SVD分解而SVD本身在符号上存在固有的模糊性如果U的一列乘以-1那么对应的右奇异向量V的同一列也乘以-1等式依然成立。所以两次运行或者不同实现得到的因子矩阵某些列可能是完全反向的。这个符号问题在矩阵SVD里通常不构成困扰因为你可以约定奇异值非负并且常用的算法库会给出稳定的结果。但换到张量场景由于展开矩阵往往是“非方阵且列远多于行”的细长矩阵SVD数值行为更复杂不同库或不同版本之间更容易出现符号差异。另一个类似的场景是PCA大家第一次接触PCA时也总会疑惑为什么不同库跑出来的主成分符号不一样。这个问题的通用解法是事后对齐在输出因子矩阵时统一把每列的第一个非零元素置为正号。def fix_sign(U): 将因子矩阵每列的符号统一到第一个绝对值最大元素为正 for i in range(U.shape[1]): # 找到该列绝对值最大的元素索引 max_idx np.argmax(np.abs(U[:, i])) if U[max_idx, i] 0: U[:, i] * -1 return U注意符号翻转不会改变重构结果但会影响因子矩阵的可解释性。在Tensorly里如果你用tucker函数做分解输出的因子矩阵符号也可能和你手工实现不同跨库比较时一定要先做符号对齐否则你会以为某个因子方向反了。5.2 秩取太大导致过拟合取太小导致欠拟合前文说过HOSVD是一个线性压缩过程所以它的“过拟合”不像深度学习那样严重但仍然存在。当秩取得太大时核心张量的参数数量会膨胀噪声部分也被精细建模这在有噪声的数据上尤其明显。举个例子我测过一个信噪比约10dB的三阶张量真实秩是(4,4,4)。当我取秩(4,4,4)时重构相对误差约0.08而取秩(15,15,15)时重构误差反而升到0.12。原因就在于秩增大后核心张量的许多位置开始记忆噪声细节这些细节在重构时并不能提高泛化能力。判断是否过拟合有个简单的经验法则对比训练集和验证集上的重构误差。在张量分解里“训练”和“验证”可以通过随机屏蔽一部分元素来模拟。如果你在屏蔽比例为10%的验证集上重构误差远大于完整数据上的误差说明秩可能过大了反之如果验证误差和训练误差很接近说明还有增加秩的空间。这个逻辑和机器学习里的泛化误差分析一致。5.3 数据归一化与发散维度的处理HOSVD对数据的尺度非常敏感。如果在构造三阶张量时某个维度的取值范围远大于其他维度比如第一个维度上的数值在0到1第二个维度上的数值在0到10000那么展开矩阵X(1)的SVD奇异值会被第二维主导导致第一个方向的有效子空间被严重低估。这其实是所有基于SVD的方法的通病也是为什么做PCA前要先做标准化。张量分解同理如果你的数据不同维度之间物理意义不同、量纲不同务必先做归一化。常用的归一化有两种一是按每个mode的纤维做标准化即将每个纤维向量减去均值除以标准差二是按整体张量的最大绝对值缩放到[-1, 1]。哪种更好取决于下游任务——如果是聚类或分类按纤维标准化更能保留形态差异如果是压缩整体的缩放就够用了。我自己的习惯是在拿到原始数据后先做一次快速探索性的HOSVD画出三个方向的第一奇异值占比观察是否存在某个方向第一奇异值几乎占了99%以上的情况。如果有多半是这个方向上有极端离群点或者量纲问题先处理数据比调分解参数有效得多。5.4 tensordot维度顺序错乱与排查再回来说技术细节里最容易出错的tensordot维度顺序。HOSVD的核心张量计算需要连续三个tensordot每一步都可能把维度顺序打乱。假如你对G的计算顺序比较随意比如先乘U3的转置、再乘U1的转置、最后乘U2的转置那么最终得到的核心张量维度顺序是乱的和后续重构的einsum对不上结果就全错了。排查手段有两个。第一个是检查核心张量的shape是否完全等于(r1, r2, r3)如果维度数量对但顺序乱shape会暴露出来——比如你明明传的ranks是[5,4,6]结果核心张量shape是(4,5,6)那顺序肯定不对。第二个是重构验证完全舍去秩截断即ranks取完整维度时HOSVD重构出的结果应该在数值误差范围内等于原始张量。这个“无损等价性”是检验维度顺序对不对的黄金标准。# 无损验证秩取满时重构误差应为接近0 full_ranks X.shape G_full, factors_full hosvd(X, full_ranks) X_recon np.einsum(rst,ir,js,kt-ijk, G_full, *factors_full) print(无损重构误差:, np.linalg.norm(X - X_recon) / np.linalg.norm(X))如果这个误差大于1e-10基本可以断定展开或重构过程中有维度顺序问题。我见过不少人把错误归咎于浮点精度其实正常SVD下误差应该到1e-13甚至更低才对。6. 应用场景与扩展方向HOSVD还能用在哪儿6.1 图像与视频压缩三阶张量的天然场景一张RGB图像天然就是三阶张量高乘宽乘通道。传统JPEG压缩是对每个通道分别做DCT再用量化表丢弃高频分量这个流程并没有利用三个通道之间的跨通道相关性。而HOSVD则可以一次性对三个维度做联合压缩把跨通道的冗余也利用起来。具体做法是把图像放进三阶张量X尺寸(H, W, 3)然后跑HOSVD取秩(r1, r2, r3)。如果r3取2甚至1就相当于在通道维度上也做了压缩这在实际场景中往往很有效因为RGB三个通道之间存在强相关性尤其在自然图像中。实测下来对512×512×3的彩色图像取秩(128, 128, 3)时压缩率约4倍PSNR还在32dB以上进一步取秩(100, 100, 2)时压缩率约8倍PSNR大概跌到28dB视觉上依然可接受。这个思路比单独压缩通道然后再合并的传统方案在低比特率下有更平滑的画质退化曲线。6.2 推荐系统与多路数据的联合建模推荐系统里的用户-物品-上下文数据是最典型的三阶张量场景。比如用户u在时间t对物品i的评分或点击行为可以组织成用户乘物品乘时间的张量然后用HOSVD同时学习三组隐向量。这样得到的结果相比矩阵分解的好处是用户和物品的隐向量不仅能刻画彼此关系还能把时间上下文的因素解耦出来。我在一个简化版的推荐场景里试过把评分张量做HOSVD分解取核心张量里每个切片固定时间维作为该时段的“交互模式矩阵”可以看到不同时段的用户偏好结构有明显差异。这种分解方式的最大优点是可解释性——Tucker分解的核心张量是稠密的你可以从核心张量的取值判断哪些“用户因子维度”和“物品因子维度”在哪个时间维度上强耦合这对运营分析很有帮助。代价是Tucker分解不能像矩阵分解那样直接对未知评分做内积预测需要额外设计重构策略一般是先补全缺失值再分解迭代进行。6.3 脑电与多通道信号处理中的HOSVD实践EEG、MEG这类多通道时间序列数据天然是三维的通道乘时间乘试验次数或频带。传统方法通常对通道和时间先做矩阵化再提取特征但这样做在试验维度上的跨样本结构就丢了。HOSVD可以同时把通道、时间、试验三个方向压缩尤其适合做群体级别的信号模式提取。一个经典的流程是把多受试者的EEG数据组织成通道×时间×受试者张量跑HOSVD后取通道因子矩阵和受试者因子矩阵做后续分析。受试者因子矩阵的行向量其实就是每个受试者在共享子空间里的“投影坐标”可以直接拿去做群体差异检验。这个方法在神经科学里叫作tensor PCA或tensor decomposition-based feature extraction它的核心优势是用一个统一模型同时估计了共享的时间模式、空间模式和受试者间变异比逐受试者提取特征再统计要稳健得多。6.4 从HOSVD到HOOI、增量分解与随机化加速最后一个想展开的方向是你在把HOSVD用于大规模数据时会遇到的性能瓶颈。HOSVD对每个展开矩阵做完整SVD当展开矩阵的尺寸很大时比如I10000, J5000, K5000X(1)就是10000×25,000,000的巨型矩阵完整SVD几乎不可能直接算。这时有两类加速手段一是使用随机化SVDrandomized SVD思路是先对展开矩阵做随机线性投影把矩阵降到低维空间再在低维空间做SVD得到近似主方向。这个思路对HOSVD尤其适用因为HOSVD本身只求每个模式的最优子空间并不需要极其精确的奇异向量随机化带来的微小误差完全在可接受范围。二是在线/增量HOSVD。如果数据是流式到达的比如监控视频按帧输入每次做全量HOSVD成本太高可以维护一套已有的因子矩阵和核心张量新数据到来时只在增量子空间上做更新。这类方法在文献里叫incremental tensor decomposition实现比离线版本复杂不少但能把每帧的处理时间降到毫秒级。HOOI作为更高阶的优化方法则是在HOSVD基础上做交替迭代精调。如果你有足够算力推荐流程是先用HOSVD得到初始因子矩阵再跑几步HOOI迭代修正。这个组合在几乎所有公开数据集上都能稳定优于单独HOSVD而且迭代轮数不需要多三到五轮就足够收敛。我个人在实际项目里用过的最顺手的组合拳是HOSVD做初始化交叉验证选秩HOOI精调三到五轮最后用张量化后的下游任务做效果评估。这套流程兼顾了速度、稳定性和可解释性你可以在自己的数据上复现一遍大概率会觉得比直接调包更有掌控感。本文还有配套的精品资源点击获取 HOSVD高阶奇异值分解：三阶张量分解原理与NumPy实现
简介矩阵分解中的SVD是处理二维数据的经典工具但当数据扩展到三维乃至更高维时如图像的通道、时间序列的批次传统的矩阵方法便难以刻画多个维度间的耦合结构。张量分解作为矩阵分解的高维推广通过多个因子矩阵与核心张量同时压缩各维度信息能够保留数据内在的多方向结构。HOSVD作为Tucker分解的一种高效实现路径分别对每个模式展开矩阵做SVD再重构核心张量计算简单且逼近误差有理论保证被广泛用于图像压缩、信号处理、推荐系统与科学计算等领域。本文从三阶张量的基本概念出发详解HOSVD的数学原理、完整NumPy代码实现、多线性秩的选择策略以及常见工程陷阱帮助读者从原理到实践系统掌握这一高维数据分析利器。 我们拿到一批三维数组模样的数据时第一直觉往往会去试PCA、SVD这类矩阵工具但一旦数据本身是“立体的”比如一张灰度图像天然就是行乘列的两维再加上通道、时间、批次这些维度SVD就撑不住了。不是不能压平了算而是强行把三维数组拉成矩阵等于把原本沿不同方向变化的规律全部搅在一起丢失了张量作为整体时那种张成的结构信息。这正是张量分解、尤其是HOSVD和Tucker分解派上用场的场景。这篇东西想聊的就是HOSVD高阶奇异值分解在三阶张量上的完整落地过程从数学原理到NumPy手写实现再到多线性秩怎么选、奇异值怎么看、哪些坑我踩过。适合理工科背景、准备把张量分解用到信号处理、图像压缩、推荐系统或科学计算里的朋友。如果你只是想调包跑个结果可以直接跳到第三节如果你想搞清楚Tucker分解和HOSVD到底什么关系、为什么HOSVD有“准最优”的说法建议从头看起。1. 三阶张量分解到底在解决什么问题1.1 从矩阵SVD不够用说起先想一个最简单的场景你有100个用户每人看了50部电影评分构成一个100乘50的矩阵。SVD能做的是把这个矩阵分解成用户因子矩阵、奇异值对角阵和电影因子矩阵从而把用户和电影各自嵌到一个低维空间里。这个逻辑成立的前提是你手里的数据确实可以被“两个方向”描述干净。但现实里多数数据没那么乖。同样是评分数据你可能还想加入时间维度比如每个用户在不同月份对电影的评分那数据就从矩阵变成三阶张量用户乘电影乘时间。此时如果仍然压平成矩阵就有两种做法一是把所有时间片叠在一起变成一个大矩阵但这等于默认时间维度和用户维度在同一个语义层面显然不合理二是对每个时间片单独跑SVD但这样不同时间片的分解结果又对不上号没法对齐比较。无论是哪种都把三个维度之间的耦合关系破坏掉了。张量分解的核心思路是既然数据是一个三阶对象那就用三个因子矩阵或者叫因子向量组从三个方向同时挤压它再配一个“核心张量”去描述维度之间的交互。这跟矩阵SVD把A拆成U乘S乘V^T是同一个逻辑的不同版本只不过原来矩阵U和V各自只有一个方向现在Tucker/HOSVD里有U1、U2、U3三个方向核心张量也不要求是对角阵。1.2 Tucker分解与HOSVD的定位差异很多人会把HOSVD、Tucker分解混为一谈在不少论文里也确实没做严格区分但实际落地时两者差异还挺重要。Tucker分解的核心公式长这样X ≈ G ×₁ A₁ ×₂ A₂ ×₃ A₃其中G是核心张量尺寸是r1×r2×r3A1、A2、A3分别是三个方向的因子矩阵也可以理解成三个方向上的“主成分方向”。这个形式表达的是原始大张量可以被压缩成一个小核心张量外加三组低维投影方向。Tucker分解本身不限定怎么求A1、A2、A3甚至不限定G是什么它只是描述了一种“多方向线性子空间逼近”的框架。HOSVD则是求解Tucker分解的一种具体算法而且名字里带“SVD”也说明它是从矩阵SVD推广过来的——分别在三个mode上做矩阵化再做普通SVD取左奇异向量作为因子矩阵最后用因子矩阵反推核心张量。所以在数学结构上HOSVD输出的结果天然就是一个Tucker分解但反过来Tucker分解的求解路径并不一定是HOSVD还有更高阶的优化算法如HOOI。这里值得一提的是HOSVD往往被叫做“Tucker分解的准最优解”。原因是Tucker的最佳低秩逼近问题在数学上是个非凸优化问题全局最优极难求而HOSVD走的路线是“分别对每个方向做截断SVD再拼接”每一步都有显式解所以计算上非常快在小规模问题上效果也很接近最优。你可以把HOSVD理解成一种贪心策略而HOOI这类算法是在HOSVD给出的基础上再迭代精调。1.3 三阶张量的本质模式、纤维与切片要真正理解HOSVD还要把三阶张量的几个基础概念过一遍不然代码里的mode-1、mode-2、mode-3很容易把人绕晕。三阶张量X的尺寸如果是I×J×K那么沿着第一维row方向固定索引i得到的是J×K的矩阵这叫第i个水平切片沿着第二维固定索引j得到I×K的矩阵沿着第三维固定索引k得到I×J的矩阵。切片对应的是“在某个维度上切一刀”是观察张量的窗口。更底层的是纤维fiber概念。固定其他维度索引只让一个维度变化得到的一维数组就是纤维。在矩阵里行和列就是两种纤维在三阶张量里纤维有三种列纤维沿第一维、行纤维沿第二维、管纤维沿第三维。SVD中对矩阵做分解本质上是在找两种纤维共享的低维结构而HOSVD就是把它推广成三种纤维的共享子空间逼近。如果把三阶张量X沿第一维展开成矩阵X(1)尺寸是I×(JK)每一列就是一条“沿第一维变化的纤维”被摊平后的结果。HOSVD中的第一步——对每个mode做矩阵展开并做SVD——其实就是在问在这个方向上、把张量看作一大撮纤维的集合时它们主要沿着哪些正交方向分布这样一想HOSVD的过程就非常直观先在三个方向上分别做PCA再把三个方向的主成分坐标拼成核心张量。2. HOSVD的核心计算路径与原理拆解2.1 三个关键步骤展开、SVD、重构核心张量HOSVD的完整算法只有三步而且每一步都有清晰的线性代数背景。第一步mode-k展开。对三阶张量X做mode-1展开得到矩阵X(1)尺寸为I×(JK)展开规则其实就是在代码里用reshape加permute实现。展开不是数据搬家而是把张量按特定顺序重排让“沿mode-1变化的元素”排成一列这样做的目的是让SVD能吃到一个标准二维矩阵。mode-2和mode-3同理只是排列顺序不同。第二步对每个展开矩阵做截断SVD。对X(1)做SVD得到左奇异向量矩阵U1尺寸I×I按奇异值从大到小排序后取前r1列作为因子矩阵A1。对X(2)和X(3)做同样的操作分别得到A2和A3。注意这里SVD不用算V因为我们只需要左奇异向量而且截断后V是在后续重构核心张量时才隐式参与的。第三步求解核心张量。有了A1、A2、A3之后核心张量的估计公式是G X ×₁ A₁ᵀ ×₂ A₂ᵀ ×₃ A₃ᵀ这个操作在物理意义上可以理解为把原始张量沿着三个方向分别投影到低维子空间里投影后得到的“坐标张量”就是核心张量G。如果你对矩阵的SVD熟悉这会让你想到X UΣVᵀ里的Σ就是原始矩阵在U和V方向上的投影坐标只是张量情形下坐标不是对角阵而是一个实心的小张量。这也是Tucker分解和矩阵SVD最直观的差异矩阵SVD的“核心”是对角的Tucker分解的“核心”是稠密的。2.2 为什么HOSVD是“准最优”而不是“最优”这里需要多解释几句因为很多人在对比HOSVD和HOOI时会把准最优理解成“HOSVD逼近效果不如HOOI”这话大体对但容易忽略一个前提HOSVD的出发点根本就不是求最优逼近。我们把HOSVD的数学本质拆开看它对X(1)、X(2)、X(3)分别做SVD这其实是在三个不同的矩阵空间中求三个正交子空间。这组子空间本身的性质很好——每个子空间都是对应展开矩阵的前r个主左奇异方向因此各自都是“最佳r维子空间”。但是三个子空间各自最佳并不代表它们张成的张量积空间是最佳的三维子空间组合。打个比方三个人各自选了对自己最优的出行路线但三个人拼成一辆车之后整体的最优路线并不等于这三条路线的笛卡尔积。这是HOSVD准最优的根本原因。不过HOSVD的误差其实有理论保证。如果原始张量X的Tucker秩是(r1,r2,r3)的“有效秩”也就是说忽略掉小奇异值分量后残余很小那么HOSVD截断后的逼近误差和理论最优之间的差距不超过最优误差的√3倍在三阶情形下更精确地说差距有界于√r1√r2√r3的量级。这个结果保证了HOSVD不是瞎猜方向而是有界的近似最优。如果追求真正的最优Tucker分解通常用的是HOOI高阶正交迭代。HOOI的思路是把Tucker分解当成一个交替优化问题固定A2、A3优化A1固定A1、A3优化A2不断迭代直到目标函数收敛。每轮迭代的核心操作其实就是在利用当前的A2、A3把张量压缩成一个小张量再对这个压缩结果做SVD提取新的A1。HOOI每一步都在提升拟合度所以理论上能找到更好的局部最优代价是慢。实践中HOSVD经常被用来给HOOI提供初始值这个组合拳几乎能对标全局最优。2.3 HOSVD、Tucker、高阶PCA之间的关系我在看论文和做实现时经常发现这些术语在不同文章里指的东西略有差别这里整理一下方便你避坑HOSVD是广义叫法有的文献也把高阶奇异值分解称为Tucker分解的HOSVD形式特指“通过逐mode SVD求Tucker分解”的算法过程。Tucker分解强调的是分解形式本身即三个因子矩阵加一个核心张量的结构不限定算法。任何能求出这个结构的算法都可以称为在做Tucker分解。higher-order SVD和Tucker-1有时是同义词在早期文献里Tucker-1指只压缩一个方向的分解Tucker-2是压缩两个方向Tucker-3是三个方向都压缩。HOSVD通常等价于Tucker-3这也是为什么大家默认HOSVD就是“全方向的Tucker”。还有一个容易混淆的概念是CP分解CANDECOMP/PARAFAC它的形式是把张量分解成若干个秩一张量之和X ≈ Σ λᵣ aᵣ ∘ bᵣ ∘ cᵣ。CP和Tucker的最大区别是核心张量是否对角——CP的核心张量强制对角Tucker不强制。因此CP的参数更少、解释更强硬但求解稳定性更差对噪声敏感Tucker则灵活得多。一句话概括HOSVD是求Tucker分解的一条具体路径而Tucker分解是你要的那个结构本身。3. 手写NumPy实现HOSVD从零搭建完整流程3.1 环境与基础工具接下来上代码我这边的实现环境是Python 3.10 NumPy 1.24不需要额外装Tensorly或PyTorchHOSVD最核心的数学操作其实就两个张量展开reshape transpose和SVDnumpy.linalg.svd。这两个操作对中小规模张量完全够用上了大规模再考虑稀疏分解或随机化SVD。先定义一个三阶张量的mode-k展开函数。这里要注意NumPy里没有原生的“张量展开”函数需要自己写而且很容易在维度顺序上出错。一个稳妥的做法是先用np.moveaxis把目标维度挪到最前再做reshape这样逻辑最清晰不容易错import numpy as np def mode_k_unfold(X, k): 将三阶张量X按第k个模式展开为矩阵 k0,1,2 对应第一、第二、第三维度 返回形状为 (X.shape[k], 其它维度乘积) 的二维数组 # 把第k维移到最前面其他地方保持原顺序 moved np.moveaxis(X, k, 0) shape moved.shape # 第一维是第k维后面所有维度合并 return moved.reshape(shape[0], -1)这个函数后面所有步骤都要复用所以花点时间确认逻辑是对的X被moveaxis后第0轴就是我们要的mode方向剩下的轴原封不动按顺序排。reshape成(shape[0], -1)后每一列恰好对应一条沿mode方向变化的纤维。3.2 HOSVD主函数实现接下来实现HOSVD主函数。核心流程是三步走三个方向展开、三个SVD取左奇异向量、投影求核心张量。我写的版本要求用户传入截断秩列表ranks [r1, r2, r3]如果你不确定怎么选秩可以先跑一遍完整SVD看奇异值分布再回来填这个参数。def hosvd(X, ranks): 对三阶张量X执行HOSVD截断版本 X: numpy.ndarray形状 (I, J, K) ranks: list或tuple[r1, r2, r3] 返回: G: 核心张量形状 (r1, r2, r3) factors: 因子矩阵列表 [U1, U2, U3] ndim X.ndim if ndim ! 3: raise ValueError(本实现仅支持三阶张量) factors [] for k in range(ndim): # 第k模式展开 X_k mode_k_unfold(X, k) # 对展开矩阵做SVD只需要左奇异向量 U, _, _ np.linalg.svd(X_k, full_matricesFalse) # 截断到前r个主方向 r ranks[k] U_r U[:, :r] factors.append(U_r) # 核心张量X ×1 U1^T ×2 U2^T ×3 U3^T G X for k in range(ndim): G np.tensordot(factors[k].T, G, axes([1], [k])) # tensordot之后维度顺序可能变化需要调整回标准顺序 # 这里采用连续tensordot后手动确认shape # 一个更稳妥的写法是循环moveaxis但上面已经够用验证下 # 实际运行中每次tensordot会消掉第k个维度并把它放到最前 # 由于我们循环k0,1,2最终G的维度顺序是 (r1, r2, r3) return G, factors这段代码有几个细节值得注意。第一np.linalg.svd默认返回的U、S、Vh三个量对展开矩阵X_k我们只关心US和Vh都可以直接丢弃。第二我用了full_matricesFalse这样当展开矩阵的列数小于行数时U的列数自动与S的长度对齐不会出现多余的全零奇异值列截断时更省心。第三tensordot的axes传参容易写错我这里在注释里已经说明第一次tensordot会消掉X的第0维把因子矩阵的转置乘上去得到的新张量维度顺序为(r1, J, K)第二次tensordot消掉第1维得到(r1, r2, K)第三次消掉第2维得到(r1, r2, r3)。由于我们恰好是按顺序消的维度顺序不用额外调整但如果你在循环里跳着消一定要记得最后修正维度顺序。3.3 从HOSVD到Tucker重构的验证光有分解还不够还得能验证重构误差。Tucker/HOSVD的重构公式是X_approx G ×₁ U1 ×₂ U2 ×₃ U3也就是把核心张量沿着三个方向用因子矩阵映射回原始空间。用一个完整的测试脚本来验证整个流程# 构造一个随机三阶张量并加入一点低秩结构便于观察 np.random.seed(42) I, J, K 30, 20, 25 R1, R2, R3 5, 4, 6 # 先构造一个低秩张量A1(30x5), A2(20x4), A3(25x6)核心张量(5,4,6) A1 np.random.randn(I, R1) A2 np.random.randn(J, R2) A3 np.random.randn(K, R3) G_true np.random.randn(R1, R2, R3) X np.einsum(ir,js,kt,rst-ijk, A1, A2, A3, G_true) X 0.05 * np.random.randn(I, J, K) # 加一点噪声 # 执行HOSVD截断秩设为真实秩 G_hat, factors hosvd(X, [R1, R2, R3]) U1, U2, U3 factors # 重构 X_approx np.einsum(rst,ir,js,kt-ijk, G_hat, U1, U2, U3) # 计算相对误差 rel_err np.linalg.norm(X - X_approx) / np.linalg.norm(X) print(f相对重构误差: {rel_err:.6f})跑完这段代码你会得到一个相对误差约0.00x量级的结果具体数值取决于噪声水平。这个验证过程非常关键如果结果误差很大大概率是维度顺序写错了而不是算法本身的问题。我在调通第一版实现时就因为tensordot的维度顺序搞错误差一直是0.7左右排查很久才发现是核心张量重构时维度错位了。3.4 用Tensorly对比验证实现正确性如果你还是担心自己的实现有bug可以拿Tensorly这个库来交叉验证。Tensorly的高阶SVD实现很成熟在学术界和工业界都被广泛使用用它作为基准对比是最快的办法。import tensorly as tl from tensorly.decomposition import tucker # Tensorly的tucker直接支持指定秩 X_tensor tl.tensor(X) core_tl, factors_tl tucker(X_tensor, rank[R1, R2, R3], initsvd, n_iter_max0) # initsvd n_iter_max0 意味着只做HOSVD初始化不做HOOI迭代注意上面这段里initsvd且n_iter_max0意思就是“只用HOSVD结果不继续迭代优化”这正好对应我们手写的HOSVD。对比一下手写版本的因子矩阵与Tensorly输出的因子矩阵会发现它们的方向完全一致全局符号可能相反不影响语义核心张量也几乎相同。符号相反这件事我们放到第五节细聊它本身是个典型的坑。4. 多线性秩怎么选奇异值谱怎么读4.1 三种方向、三种奇异值谱HOSVD对你透露信息最丰富的时刻其实发生在选择秩之前的那个完整SVD计算。对三个模式展开矩阵分别做SVD后你会得到三组奇异值s1长度min(I, JK)、s2长度min(J, IK)、s3长度min(K, IJ)。这三组奇异值刻画的是张量沿三个方向各自的“能量分布”。矩阵SVD里大家习惯用奇异值累积能量占比来决定主成分个数这个习惯在张量情形下依然有效只是要分别对三个方向做判断。比如s1的前几个值迅速衰减、后面基本平了那就说明第一维方向上有效子空间维度很低s2衰减慢说明第二维方向信息冗余较少需要更多秩去刻画。画图时我建议把三个谱画在同一张图上横轴都是索引纵轴是奇异值或者归一化到最大奇异值。看图时重点观察两个地方第一个是“拐点”在哪里也就是奇异值从大变小趋势明显变缓的位置这个点往往是有效秩的候选第二个是一开始的几个值是不是占了绝对主导如果前两个奇异值占比超过90%那么这个方向取r2基本就够了。4.2 常用选秩方法对比选秩没有统一答案不同任务有不同的考虑方向。我把常见的几种方法连同适用场景一起列出来方法思路适用场景注意事项奇异值累积能量保留总能量90%/95%对应的前r个奇异值压缩、降噪需要对三组谱分别判断取max或逐方向取拐点/差值最大找奇异值二阶差分最大的位置低噪声、结构明显的数据对噪声敏感适合肉眼辅助判断交叉验证把张量元素随机屏蔽一部分用不同秩重构选误差最小的秩组合缺失值填充、推荐系统计算量大要小心随机屏蔽对空间结构的破坏任务驱动不追求最优重构而是按下游任务效果选秩分类、回归、聚类和分解目标分离容易过拟合需要控制验证集贝叶斯/自动秩选择用概率模型估计秩后验分布复杂噪声、分布式数据实现复杂常规任务不需要上这种手段实操中我常用的办法是先用能量占比确定每个方向秩的粗范围再在这个小范围内做交叉验证。因为三阶张量的秩组合空间往往是r1乘r2乘r3的笛卡尔积穷举很不现实但如果你把三个方向的秩分别限制在2到3个候选值里组合数也就十几个跑一轮交叉验证完全可行。4.3 压缩率与重构误差的权衡实操选秩本质上是在回答“我打算牺牲多少信息换取多大压缩”这个问题。为了量化这个权衡我习惯在选定秩之前先把“压缩率-重构误差曲线”画出来。压缩率的定义可以这样算原始张量需要存储I×J×K个floatTucker分解后需要存储r1×r2×r3核心张量加上(I×r1 J×r2 K×r3)因子矩阵所以压缩率是ratio (r1×r2×r3 I×r1 J×r2 K×r3) / (I×J×K)比如I100, J100, K100秩取(10,10,10)那么存储量是1000 1000 1000 1000 4000原始是1000000压缩率高达250倍。如果取(50,50,50)存储量就是125000 5000 5000 5000 140000压缩率约7.1倍。可见秩对压缩率的影响是非线性的尤其三个方向同时增大时核心张量的体积按三次方增长。有了压缩率的概念后我的选秩流程就很机械先看奇异值谱得到候选秩集合然后对每个候选秩组合算一遍重构误差和压缩率画成二维散点或者平行坐标图最后根据业务对“压缩率优先”还是“误差优先”的偏好定夺。在图像压缩任务里往往压到20倍以内肉眼几乎无差异这时候多花一倍的存储换0.1%的精度改善就没必要了。5. 实操中绕不开的坑与排查技巧5.1 奇异值分解的符号不稳定与左奇异向量漂移我在最初跑HOSVD时踩的第一个坑就是“两次运行结果看起来不一致”。后来发现HOSVD对因子矩阵做了SVD分解而SVD本身在符号上存在固有的模糊性如果U的一列乘以-1那么对应的右奇异向量V的同一列也乘以-1等式依然成立。所以两次运行或者不同实现得到的因子矩阵某些列可能是完全反向的。这个符号问题在矩阵SVD里通常不构成困扰因为你可以约定奇异值非负并且常用的算法库会给出稳定的结果。但换到张量场景由于展开矩阵往往是“非方阵且列远多于行”的细长矩阵SVD数值行为更复杂不同库或不同版本之间更容易出现符号差异。另一个类似的场景是PCA大家第一次接触PCA时也总会疑惑为什么不同库跑出来的主成分符号不一样。这个问题的通用解法是事后对齐在输出因子矩阵时统一把每列的第一个非零元素置为正号。def fix_sign(U): 将因子矩阵每列的符号统一到第一个绝对值最大元素为正 for i in range(U.shape[1]): # 找到该列绝对值最大的元素索引 max_idx np.argmax(np.abs(U[:, i])) if U[max_idx, i] 0: U[:, i] * -1 return U注意符号翻转不会改变重构结果但会影响因子矩阵的可解释性。在Tensorly里如果你用tucker函数做分解输出的因子矩阵符号也可能和你手工实现不同跨库比较时一定要先做符号对齐否则你会以为某个因子方向反了。5.2 秩取太大导致过拟合取太小导致欠拟合前文说过HOSVD是一个线性压缩过程所以它的“过拟合”不像深度学习那样严重但仍然存在。当秩取得太大时核心张量的参数数量会膨胀噪声部分也被精细建模这在有噪声的数据上尤其明显。举个例子我测过一个信噪比约10dB的三阶张量真实秩是(4,4,4)。当我取秩(4,4,4)时重构相对误差约0.08而取秩(15,15,15)时重构误差反而升到0.12。原因就在于秩增大后核心张量的许多位置开始记忆噪声细节这些细节在重构时并不能提高泛化能力。判断是否过拟合有个简单的经验法则对比训练集和验证集上的重构误差。在张量分解里“训练”和“验证”可以通过随机屏蔽一部分元素来模拟。如果你在屏蔽比例为10%的验证集上重构误差远大于完整数据上的误差说明秩可能过大了反之如果验证误差和训练误差很接近说明还有增加秩的空间。这个逻辑和机器学习里的泛化误差分析一致。5.3 数据归一化与发散维度的处理HOSVD对数据的尺度非常敏感。如果在构造三阶张量时某个维度的取值范围远大于其他维度比如第一个维度上的数值在0到1第二个维度上的数值在0到10000那么展开矩阵X(1)的SVD奇异值会被第二维主导导致第一个方向的有效子空间被严重低估。这其实是所有基于SVD的方法的通病也是为什么做PCA前要先做标准化。张量分解同理如果你的数据不同维度之间物理意义不同、量纲不同务必先做归一化。常用的归一化有两种一是按每个mode的纤维做标准化即将每个纤维向量减去均值除以标准差二是按整体张量的最大绝对值缩放到[-1, 1]。哪种更好取决于下游任务——如果是聚类或分类按纤维标准化更能保留形态差异如果是压缩整体的缩放就够用了。我自己的习惯是在拿到原始数据后先做一次快速探索性的HOSVD画出三个方向的第一奇异值占比观察是否存在某个方向第一奇异值几乎占了99%以上的情况。如果有多半是这个方向上有极端离群点或者量纲问题先处理数据比调分解参数有效得多。5.4 tensordot维度顺序错乱与排查再回来说技术细节里最容易出错的tensordot维度顺序。HOSVD的核心张量计算需要连续三个tensordot每一步都可能把维度顺序打乱。假如你对G的计算顺序比较随意比如先乘U3的转置、再乘U1的转置、最后乘U2的转置那么最终得到的核心张量维度顺序是乱的和后续重构的einsum对不上结果就全错了。排查手段有两个。第一个是检查核心张量的shape是否完全等于(r1, r2, r3)如果维度数量对但顺序乱shape会暴露出来——比如你明明传的ranks是[5,4,6]结果核心张量shape是(4,5,6)那顺序肯定不对。第二个是重构验证完全舍去秩截断即ranks取完整维度时HOSVD重构出的结果应该在数值误差范围内等于原始张量。这个“无损等价性”是检验维度顺序对不对的黄金标准。# 无损验证秩取满时重构误差应为接近0 full_ranks X.shape G_full, factors_full hosvd(X, full_ranks) X_recon np.einsum(rst,ir,js,kt-ijk, G_full, *factors_full) print(无损重构误差:, np.linalg.norm(X - X_recon) / np.linalg.norm(X))如果这个误差大于1e-10基本可以断定展开或重构过程中有维度顺序问题。我见过不少人把错误归咎于浮点精度其实正常SVD下误差应该到1e-13甚至更低才对。6. 应用场景与扩展方向HOSVD还能用在哪儿6.1 图像与视频压缩三阶张量的天然场景一张RGB图像天然就是三阶张量高乘宽乘通道。传统JPEG压缩是对每个通道分别做DCT再用量化表丢弃高频分量这个流程并没有利用三个通道之间的跨通道相关性。而HOSVD则可以一次性对三个维度做联合压缩把跨通道的冗余也利用起来。具体做法是把图像放进三阶张量X尺寸(H, W, 3)然后跑HOSVD取秩(r1, r2, r3)。如果r3取2甚至1就相当于在通道维度上也做了压缩这在实际场景中往往很有效因为RGB三个通道之间存在强相关性尤其在自然图像中。实测下来对512×512×3的彩色图像取秩(128, 128, 3)时压缩率约4倍PSNR还在32dB以上进一步取秩(100, 100, 2)时压缩率约8倍PSNR大概跌到28dB视觉上依然可接受。这个思路比单独压缩通道然后再合并的传统方案在低比特率下有更平滑的画质退化曲线。6.2 推荐系统与多路数据的联合建模推荐系统里的用户-物品-上下文数据是最典型的三阶张量场景。比如用户u在时间t对物品i的评分或点击行为可以组织成用户乘物品乘时间的张量然后用HOSVD同时学习三组隐向量。这样得到的结果相比矩阵分解的好处是用户和物品的隐向量不仅能刻画彼此关系还能把时间上下文的因素解耦出来。我在一个简化版的推荐场景里试过把评分张量做HOSVD分解取核心张量里每个切片固定时间维作为该时段的“交互模式矩阵”可以看到不同时段的用户偏好结构有明显差异。这种分解方式的最大优点是可解释性——Tucker分解的核心张量是稠密的你可以从核心张量的取值判断哪些“用户因子维度”和“物品因子维度”在哪个时间维度上强耦合这对运营分析很有帮助。代价是Tucker分解不能像矩阵分解那样直接对未知评分做内积预测需要额外设计重构策略一般是先补全缺失值再分解迭代进行。6.3 脑电与多通道信号处理中的HOSVD实践EEG、MEG这类多通道时间序列数据天然是三维的通道乘时间乘试验次数或频带。传统方法通常对通道和时间先做矩阵化再提取特征但这样做在试验维度上的跨样本结构就丢了。HOSVD可以同时把通道、时间、试验三个方向压缩尤其适合做群体级别的信号模式提取。一个经典的流程是把多受试者的EEG数据组织成通道×时间×受试者张量跑HOSVD后取通道因子矩阵和受试者因子矩阵做后续分析。受试者因子矩阵的行向量其实就是每个受试者在共享子空间里的“投影坐标”可以直接拿去做群体差异检验。这个方法在神经科学里叫作tensor PCA或tensor decomposition-based feature extraction它的核心优势是用一个统一模型同时估计了共享的时间模式、空间模式和受试者间变异比逐受试者提取特征再统计要稳健得多。6.4 从HOSVD到HOOI、增量分解与随机化加速最后一个想展开的方向是你在把HOSVD用于大规模数据时会遇到的性能瓶颈。HOSVD对每个展开矩阵做完整SVD当展开矩阵的尺寸很大时比如I10000, J5000, K5000X(1)就是10000×25,000,000的巨型矩阵完整SVD几乎不可能直接算。这时有两类加速手段一是使用随机化SVDrandomized SVD思路是先对展开矩阵做随机线性投影把矩阵降到低维空间再在低维空间做SVD得到近似主方向。这个思路对HOSVD尤其适用因为HOSVD本身只求每个模式的最优子空间并不需要极其精确的奇异向量随机化带来的微小误差完全在可接受范围。二是在线/增量HOSVD。如果数据是流式到达的比如监控视频按帧输入每次做全量HOSVD成本太高可以维护一套已有的因子矩阵和核心张量新数据到来时只在增量子空间上做更新。这类方法在文献里叫incremental tensor decomposition实现比离线版本复杂不少但能把每帧的处理时间降到毫秒级。HOOI作为更高阶的优化方法则是在HOSVD基础上做交替迭代精调。如果你有足够算力推荐流程是先用HOSVD得到初始因子矩阵再跑几步HOOI迭代修正。这个组合在几乎所有公开数据集上都能稳定优于单独HOSVD而且迭代轮数不需要多三到五轮就足够收敛。我个人在实际项目里用过的最顺手的组合拳是HOSVD做初始化交叉验证选秩HOOI精调三到五轮最后用张量化后的下游任务做效果评估。这套流程兼顾了速度、稳定性和可解释性你可以在自己的数据上复现一遍大概率会觉得比直接调包更有掌控感。本文还有配套的精品资源点击获取