张量奇异值分解：从矩阵SVD到t-SVD的原理与NumPy实现-平芜编程栈

做数据处理的人，对“矩阵奇异值分解”应该都不陌生。但当你碰到的对象是一段视频、一组多光谱图像，或者带上下文维度的用户行为数据，强行把这些数据展开成一个大矩阵再去做SVD，总会有一种“信息被压扁”的别扭感：空间关系还在，可时间方向的连续性、不同波段之间的相关性，都被塞进了一个二维平面里。我第一次接触张量奇异值分解框架（tensor Singular Value Decomposition，即t-SVD）是在视频降噪项目里，当时对比了矩阵展开与张量低秩逼近两种路线，效果差异非常明显。后来把整套方法吃透，才发现它并不是矩阵SVD的简单三维版，而是以循环卷积为基本乘法的一套完整分解体系。

这篇文章我从原理讲到代码，再讲到实际项目里的踩坑记录。不涉及复杂推导，也不绕弯子，目标就是让你看完能自己复现一个可用的t-SVD框架，并且知道它适合解决什么问题。

1. 它和矩阵 SVD 的差异：为什么“展开成矩阵”不够

很多人第一次听到t-SVD，第一反应是“把矩阵SVD推广到三维”。这种理解方向是对的，但不能只停留在“多了一个维度”这个层面上。矩阵SVD处理的是一个二维表，t-SVD处理的是一个三阶张量，关键是t-SVD重新定义了张量之间的乘法。这个乘法不是逐元素相乘，也不是矩阵乘法直接外扩，而是在第三维上引入了循环卷积的机制。只有把这个乘法机制搞明白，才能理解为什么t-SVD能把矩阵SVD丢掉的结构信息捡回来。

1.1 数据展开成矩阵到底丢了什么

拿一段灰度视频来举例。假设视频有100帧，每帧分辨率是200×200。把它展开成矩阵有几种常见做法：把每一帧拉成一行，得到100×40000的矩阵；或者把每一帧拉成一列，得到40000×100的矩阵。不管用哪种方式，展开之后矩阵的每一行或每一列都对应某一个完整的帧。

问题出在帧与帧之间的关系上。视频里的物体运动是连续的，同一辆车在第5帧、第6帧、第7帧里只平移了几个像素，这种位移在矩阵展开后会被保留在行与行之间，但矩阵SVD处理的是行与列的相关结构。展开成矩阵后，行索引是帧号，列索引是像素位置，矩阵SVD里出现的奇异向量对应的是“像素方向的模式”和“帧方向的模式”，它没有办法显式地表达“第6帧是第5帧平移两格得到”这种跨帧的位移结构。

如果换一种视角，把视频视为一个整体三阶张量，第一维是高度，第二维是宽度，第三维是帧索引，那么t-SVD的乘法天然包含了循环卷积。循环卷积在信号处理里本身就与平移、位移密切相关。t-SVD处理视频类数据时，能更好地保留相邻帧之间的连续性，这是展开成矩阵后很难做到的。

1.2 t-product 到底是怎样一种乘法

t-SVD里的核心运算是t-product。一个三阶张量可以看成由若干正面切片组成，切片就是沿着第三维方向切出来的多个矩阵。例如形状为n1×n2×n3的张量，它的第k个正面切片就是一个n1×n2的矩阵。

两个张量A和B做t-product，第一步是把A的正面切片按照第三维方向排列成一个块循环矩阵，这个块循环矩阵比A大得多，维度是n1×n3乘n2×n3。第二步再把B的前两个维度方向做unfold操作，也就是把B的正面切片竖着堆成一个更高的大矩阵。第三步执行普通矩阵乘法。最后把乘法结果折叠回三阶张量的形状。整个过程抽象、繁琐，但它背后有一个非常明朗的直觉：t-product在第三维方向上做的是循环卷积，而前两维方向上做的是标准矩阵乘法。

解释卷积这里用一个生活化类比。假设有一个滑动窗口从数据的一头滑到另一头，每次移动一步，每一步都把窗口内对应位置的值做加权求和。如果窗口循环地回到起点继续滑，那么这个过程就是循环卷积。t-product把第三维上的这种循环滑动和矩阵乘法结合在一起，所以它对周期性强或者存在连续位移的数据特别敏感。

1.3 矩阵SVD与t-SVD的关键区别

下面这张表可以帮助快速定位两者的不同。

对比项	矩阵 SVD	t-SVD
输入对象	二维矩阵	三阶张量
基本乘法	普通矩阵乘法	t-product（第三维循环卷积）
奇异值类型	非负标量	沿第三维分布的奇异管子
秩的定义	非零奇异值个数	tubal rank（非零奇异管个数）
计算路径	一次矩阵分解	FFT后逐切片做矩阵分解

矩阵SVD的奇异值是一维排序的一组非负数，t-SVD的“奇异值”在张量里变成了更复杂的结构。它不再是一个数字，而是一段沿着第三维方向存在的序列，常被称为奇异管子。t-SVD的秩也不是数一数有几个非零的奇异值，而是数一数有多少根非零的奇异管子。这些概念上的差别，直接导致了t-SVD在保持张量结构上比矩阵展开更稳。

2. 核心数学：t-SVD 的定义与计算路径

t-SVD的数学表达式与矩阵SVD形式上很像，可以写成A = U * S * V'，其中*表示t-product，V'表示张量转置。不过这个形式背后涉及的块循环矩阵、unfold、fold操作，让第一次接触的人容易头大。我在学习时找到了一条更直接的理解路径：先沿第三维做傅里叶变换。

2.1 从傅里叶域打开突破口

t-product之所以复杂，是因为它牵扯循环卷积在第三维上的纠缠。而卷积定理告诉我们，时域里的卷积对应到频域就是乘法。如果先把张量沿第三维做FFT，那么原本纠缠在一起的第三维关系就会被解耦，每个频率切片上的操作会变成独立的矩阵乘法。

这个视角是整个t-SVD框架最关键的突破口。按照这个思路，t-SVD的计算可以不再直接构造巨大的块循环矩阵，而是转化成三步：沿第三维做FFT；对每个正面切片独立做一次矩阵SVD；把所有切片的结果沿第三维做逆FFT。这三个步骤在数学上是严格等价的，但计算复杂度完全不是一个数量级。

2.2 傅里叶域的三步操作

第一步是对三阶张量A沿第三维做FFT，得到A_hat。这一步只影响第三维，前两维的数值不发生任何变化。在NumPy里，一行代码就能完成：np.fft.fft(A, axis=2)。

第二步是遍历第三维的所有切片，对A_hat[:, :, k]执行一次标准矩阵SVD。这里k对应某频率下的切片索引。每个切片都做矩阵SVD，得到对应频率下的U_hat、S_hat和V_hat。因为各切片之间互相独立，这一步非常容易并行化，也是t-SVD对计算比较友好的一个重要原因。

第三步是把U_hat、S_hat、V_hat沿第三维做逆FFT，得到U、S、V三个因子张量。对实数张量来说，逆变换之后会存在很小的复数残留，直接取实部即可，这不会影响理论正确性。

2.3 截断操作与tubal rank

实际使用中几乎不会去做完整t-SVD，而是只保留前r个主要模式。这个操作叫截断t-SVD。截断后的重构失效很大程度取决于r选在哪里。

截断时要特别注意一点：在傅里叶域里，每个切片保留的奇异值数目必须保持统一。你不能在第一个切片保留前10个奇异值、在第二个切片保留前12个，这样逆变换回去之后，得到的张量在第三维方向上并不具备应有的共轭对称结构，无法对应一个合法的截断低秩张量。实践上所有切片都统一取前r个奇异值，这是最稳妥的做法。

tubal rank的定义依赖于S张量中非零奇异管子的数量。直观理解就是张量在t-product乘法下真正独立的有多少“通道”。某个数据的矩阵展开秩可能很高，但它沿第三维方向很可能存在足够强的一致性，t-SVD能在tubal rank较低的情况下逼近原数据。这正是t-SVD能处理视频、高光谱图像这类数据的底气所在。

3. 手写实现：几十行 NumPy 跑通 t-SVD 框架

如果只停留在看公式和推导，很难建立体感。我自己入门t-SVD时，把三件事亲手写出来之后才算真正掌握：张量转置、t-product、t-SVD分解。下面给一套可以直接复制运行的NumPy实现，没有其他第三方依赖。

3.1 最小可运行的核心三个函数

import numpy as np def t_transpose(A): # 张量转置：先交换前两维，再把第三维从第二块开始整体倒序 n1, n2, n3 = A.shape At = np.zeros((n2, n1, n3), dtype=A.dtype) At[:, :, 0] = A[:, :, 0].T At[:, :, 1:] = np.flip(A.transpose(1, 0, 2)[:, :, 1:], axis=2) return At def t_product(A, B): # t-product：沿第三维做FFT后，逐切片矩阵乘法，再逆FFT n1 = A.shape[0] n2 = B.shape[1] n3 = A.shape[-1] A_hat = np.fft.fft(A, axis=2) B_hat = np.fft.fft(B, axis=2) C_hat = np.zeros((n1, n2, n3), dtype=np.complex128) for k in range(n3): C_hat[:, :, k] = A_hat[:, :, k] @ B_hat[:, :, k] return np.fft.ifft(C_hat, axis=2).real def t_svd(A): n1, n2, n3 = A.shape A_hat = np.fft.fft(A, axis=2) U_hat = np.zeros((n1, n1, n3), dtype=np.complex128) S_hat = np.zeros((n1, n2, n3), dtype=np.complex128) V_hat = np.zeros((n2, n2, n3), dtype=np.complex128) for k in range(n3): u, s, vh = np.linalg.svd(A_hat[:, :, k], full_matrices=True) U_hat[:, :, k] = u S_hat[:, :, k] = np.diag(s) V_hat[:, :, k] = vh.conj().T U = np.fft.ifft(U_hat, axis=2).real S = np.fft.ifft(S_hat, axis=2).real V = np.fft.ifft(V_hat, axis=2).real return U, S, V

t_transpose里的“第三维倒序”可能初看很奇怪。原因在于实张量求转置时，不仅要交换前两维，还要在第三维上把所有非首块切片逆序排列，这样才能保证t-product下的对偶关系成立。我一开始在这里漏掉了翻转，导致重构误差始终无法降到预期量级，后来翻了定义才发现。

3.2 用重构误差验证正确性

写完代码第一件事不是直接拿去用，而是验证是否满足A = U * S * V'这个重构关系。验证方法非常直接。

# 随机生成一个形状为(5, 6, 8)的三阶张量 A = np.random.rand(5, 6, 8) U, S, V = t_svd(A) recon = t_product(U, t_product(S, t_transpose(V))) # 重构误差应接近浮点精度 print(np.linalg.norm(recon - A))

我本地跑出来的误差通常在1e-8量级。如果你跑出误差很大，优先检查两件事：t_transpose里是不是漏掉了第三维翻转；t_product里是不是没有按频率切片逐块做矩阵乘法，而是直接用了整体数组乘法。

还有一个常见的坑：逆FFT之后存在复数残留。上面代码里用.real取实部是合理的，因为输入是实数张量，理论上分解因子也可以是实数张量。但要注意不要直接拿逆变换后的复数结果去算重构误差，否则会得到莫名其妙的多余虚部。

3.3 截断版本的实现思路

如果要做低秩逼近，差分很小。只需要在t_svd的循环内部，对每个频率切片保留前r个奇异值和对应的前r列奇异向量。

def t_svd_truncated(A, r): n1, n2, n3 = A.shape A_hat = np.fft.fft(A, axis=2) U_hat = np.zeros((n1, r, n3), dtype=np.complex128) S_hat = np.zeros((r, r, n3), dtype=np.complex128) V_hat = np.zeros((n2, r, n3), dtype=np.complex128) for k in range(n3): u, s, vh = np.linalg.svd(A_hat[:, :, k], full_matrices=False) U_hat[:, :, k] = u[:, :r] S_hat[:, :, k] = np.diag(s[:r]) V_hat[:, :, k] = vh[:r, :].conj().T U = np.fft.ifft(U_hat, axis=2).real S = np.fft.ifft(S_hat, axis=2).real V = np.fft.ifft(V_hat, axis=2).real return U, S, V

截断后的重构误差自然不可能是浮点精度，具体误差大小取决于数据本身的低秩度。对典型视频数据，r取总维度的10%到30%往往就能获得不错的重构质量，但最优值必须结合奇异值分布来判断。

4. 实战经验：t-SVD 在项目里怎么落地

直接手写t-SVD并不是为了证明数学，而是为了把它用进具体任务。就我的经验来说，t-SVD在视频降噪和缺失数据补全这两个场景下手效果最直接。

4.1 视频低秩去噪：一个可直接套用的流程

一个简单的视频降噪流程可以这样组织：把灰度视频组织成形状为(height, width, frames)的三阶张量，直接做截断t-SVD，再用低秩重构得到降噪结果。

我在一个模拟项目中处理过一段帧数较多的灰度视频，每帧尺寸为240×320，总帧数为60。程序里直接将视频张量送入t_svd_truncated，设定r为12，重构之后峰值信噪比比矩阵展开方案高了2dB以上。这个差距不是偶然现象，因为矩阵展开把时间维度强行抹平，而t-SVD在第三维方向保留了帧序列的连续关系，对随机加性噪声的抑制效果更好。

实际操作中需要把视频分块。整段视频动辄几百上千帧，直接对整个张量做t-SVD，内存压力和计算时间都会非常可观。我习惯用滑动窗口按帧索引切块，比如每32帧一个窗口、步长16帧，窗口与窗口之间重叠一部分，重构后再对重叠区域做等权平均。这样既能控制内存峰值，又能避免窗口边界出现明显的跳变。

4.2 张量补全：面向缺失值的低秩填充

t-SVD在图像修复、推荐系统缺失值填充这些场景里也经常出现。以缺失像素补全为例，思路可以简化为迭代：先把缺失位置用某个初值填充，接着做一次截断t-SVD得到低秩重构，然后用重构值替换缺失位，保留已知像素的原始值，一直迭代到误差收敛。

这里的关键在于更新策略。每次迭代时，已知位置的值必须被固定住，不能被低秩重构覆盖，否则模型的优化目标会受到干扰。实际使用中，每轮迭代后观察重构的变化量，变化量足够小就可以停止。

4.3 什么样的数据适合 t-SVD 框架

t-SVD不是万能的，它适合第三维方向存在连续模式或周期结构的三阶数据。视频帧序列、高光谱图像的波段序列、多时刻传感器阵列信号都在这个范围内。反例也有：如果第三维只是几个互不相关的分组标签，比如乱七八糟的类别编号，t-SVD的循环卷积假设就很难发挥作用，这里用矩阵方法或者简单的矩阵化张量分解可能更合适。

提示：判断数据适不适合t-SVD，可以先看第三维是否存在顺序意义。如果索引顺序换一个排列，数据的语义就完全改变，那么t-SVD大概率能捕捉到有价值的结构；如果换个顺序数据仍然成立，那更像普通的三维表格，t-SVD收益有限。

5. 踩坑速查：常遇到的问题与处理思路

手写实现和实际应用过程里，有一些问题是几乎肯定会碰到的。我把它们整理成速查清单，帮你省下排查时间。

5.1 复数残留与逆变换的实部问题

用FFT实现的t-SVD，逆变换后必然会有浮点精度的复数残留。这是正常现象，不应该让你怀疑代码写错了。问题出在处理方式：如果你直接拿带虚部的复数结果参与后续运算，重构误差会一直比理论值大一点。解决办法是统一在逆FFT之后取实部，后续所有运算都基于实数张量进行。

5.2 内存爆炸与尺寸选择

t-SVD的内存峰值主要来自傅里叶域的奇异向量张量。假设一个张量是300×300×100，如果用full_matrices=True，U_hat就要存300×300×100的复数数组，内存消耗非常可观。我实际跑的时候强烈建议用full_matrices=False，只保留正方形奇异向量部分；或者直接使用截断版本，把r控制在合理范围。第三个维度也不宜选得过大，一般控制在30到200之间。第三维太大时，逐切片SVD的循环次数会明显拖慢速度，可以考虑用批量并行方式处理切片，或者用快速矩阵库来加速。

5.3 截断秩怎么选才有依据

很多初学者把矩阵SVD的奇异值阈值法直接搬来用，但t-SVD的奇异值已经变成奇异管子，不再是简单的一维数。一个比较有效的启发式做法：把所有频率切片上的奇异值平方汇总到对应r上，形成一条沿r方向的能量累积曲线，再选一个能量保持到接近100%的拐点作为r。

另一个更直观的验证方式是把r从1开始逐渐增大，每增加一次就计算一次重构误差。如果误差下降曲线在某个r处出现明显转折，往往对应的就是数据的主要模式数量。需要注意的是，t-SVD的奇异值呈管子状，不同频率上的奇异值变化模式不同，因此观察能量累积曲线时不要只看一两个切片，而是要看整体汇总，不然很容易被单个频段的假信号误导。