1. 一个特别的矩阵Toeplitz矩阵的结构和气质1.1 一条对角线一个数Toeplitz矩阵的定义与存储我第一次在一个线性求解器里看到Toeplitz矩阵时第一反应是“这不就是普通矩阵吗顶多长得有点规律”。直到后来把迭代法的核心算子在n从1000推到100000时才意识到这个结构简直是天赐的礼物——只要你会用FFT加速原本平方级的矩阵乘法可以降到O(n log n)实测能快两到三个数量级。Toeplitz矩阵的定义一句话就够矩阵的第i行第j列元素只依赖于i-j这个差即T[i, j] t_{i-j}。也就是说每条从左上到右下的对角线上的元素全部相同。举个例子一个4乘4的Toeplitz矩阵长这样t0 t_{-1} t_{-2} t_{-3} t1 t0 t_{-1} t_{-2} t2 t1 t0 t_{-1} t3 t2 t1 t0整个矩阵看起来很大但其实只需要第一列和第一行就能完整描述它。第一列是t0、t1、t2、t3第一行是t0、t_{-1}、t_{-2}、t_{-3}其余所有位置都是这两个序列平移出来的结果。这意味着存储一个n乘n的Toeplitz矩阵只需要O(n)个浮点数而不是O(n^2)。很多实际场景里我们甚至根本不会把矩阵显式构造出来而是只保留第一行和第一列的向量需要某个位置的元素时再按差值索引这一个人口卫生的改进就能省掉巨大内存。1.2 它从哪里冒出来Toeplitz矩阵不是数学课本上拿来练手的抽象玩具工程里到处都是。信号处理中一个有限冲激响应滤波器对信号的卷积过程写成矩阵形式就是Toeplitz矩阵乘向量时间序列分析里的自协方差矩阵是Toeplitz等间距网格上求解偏微分方程时某些差分算子离散出来也是Toeplitz或分块Toeplitz。我最常打交道的场景是图像去模糊。点扩散函数对整幅图的卷积可以表示为一个分块Toeplitz算子矩阵规模动辄上百万。如果显式构造这个矩阵内存直接爆掉就算构造得出来迭代法里每一次矩阵乘向量都要扫一遍百万乘百万的矩阵这在算力上完全不可接受。另一类常见的场景是地震数据处理和雷达信号处理。这些领域里的反演问题几乎都离不开Toeplitz结构因为传感器按等间隔采样物理规律本身天然具备平移不变性于是离散化之后所有对角线的平移对称性就保留下来了。1.3 加速空间从哪里来先做一个简单的复杂度估算。普通n乘n矩阵乘向量按朴素实现是n^2次乘加。n100000的时候n^2是10的10次方单核每秒几亿次浮点运算也要几十秒。而同样的Toeplitz矩阵乘向量用FFT加速之后复杂度是O(n log n)n100000时大约一百万出头次操作再算上FFT本身的常数因子实际也只需要毫秒级。这个差距不是“优化了一点”而是“从一个不可用的数量级跳到了另一个可以日常使用的数量级”。很多原本不敢做的反演、滤波、模拟计算只要把矩阵和向量的乘法换成FFT加速版本整个算法就从“只能处理小规模玩具数据”变成“可以直接上生产规模数据”。不过这里要强调一句FFT加速并不是把普通矩阵乘法“优化得更快”而是利用了Toeplitz矩阵内部的强结构。如果你把一个没有结构的稠密矩阵硬塞给FFT结果只会更慢。理解清楚哪些矩阵有类似的结构才能把这个方法用在正确的地方。2. 先看懂普通矩阵乘法行观点与列观点2.1 行观点逐行内积教科书的默认视角聊FFT加速之前我想先把“矩阵乘法行观点和列观点”这组概念掰开讲一讲。很多人在看懂FFT加速代码之前卡住的不是FFT本身而是对矩阵乘法的理解停留在单一视角上。行观点是我们最熟悉的那种定义方式要计算C A乘BC的第i行第j列等于A的第i行点乘B的第j列。也就是说先把A的第i行拿过来把B的第j列拿过来做一次内积得到结果矩阵的一个元素。这个视角的好处是直接、可操作适合手算和写最朴素的循环。代码长这样三层循环外层i中间j内层k。坏处是它对“结构”不敏感。你很难从行观点看出来“这个矩阵的每一行只是上一行平移了一下”于是也就很难想象“原来这一步运算本质上是一次卷积”。行观点在计算机底层也有自己的优势。如果你的数据是按行存储的那么访问A的第i行时内存连续cache命中率很高。这也是为什么很多基础线性代数库在实现小矩阵乘法时会偏向行内积策略。2.2 列观点列向量的加权组合列观点则是完全不同的另一种理解方式C的每一列等于A作用于B的对应列的结果。换句话说先把B拆成一列一列的向量然后用A去乘每一个列向量得到C对应的每一列。再往外推一步还能把A乘B解释成对A的列向量做线性组合C的每一列是A中各个列向量按B对应列的元素加权求和。写成外积形式就是C Σ_k A[:,k] 外积 B[k,:]。这个形式乍一看有点绕但它才是理解很多高性能计算技巧的关键。为什么这么说因为行观点把矩阵乘法当成一堆独立的内积每个输出元素都要访问A的一行和B的一列计算量虽然相同但数据访问模式很散。列观点把问题看成“把B的每一列逐一映射成C的每一列”这个映射通常可以整体做优化。现代BLAS里的分块矩阵乘法、rank-k更新本质上都是列观点或者说外积观点的发展。对Toeplitz矩阵来说列观点的价值更大。你去看Toeplitz的每一列第j列相对于第j-1列只是整体上下移动了一个位置。矩阵乘向量y T乘x按列观点看就是“把T的每一列拿出来乘以x对应的权重然后全部相加”。由于这些列彼此之间是平移关系所以这个加权求和的过程本质上就是一个卷积。这个洞察直接通向FFT。2.3 两种视角如何指向FFT把两种视角放在一起看你会发现行观点适合写正确性检查代码列观点适合想优化算法。如果你用行观点去观察T乘x你看到的是第i行和第x做内积然后第i1行又和第x做内积。每一行都很长n行就有n个n维内积。虽然相邻两行之间有平移关系但行观点下你要刻意去发现这种关系才能把矩阵元素和x的对应关系重新组织起来。如果用列观点情况完全不同。y T乘xx只是一组权重T的列才是一组底向量。列与列之间就是位移关系于是“加权求和”就是“以x为核对t序列做一次卷积”。卷积定理马上可以接上时域卷积等于频域乘积。接下来用FFT做加速几乎是水到渠成的结论。所以我在后文的实现和排错里都会反复用到这两个视角写朴素版本对照时用行观点设计FFT加速方案时用列观点。明白这一点后面的代码才不会看得一头雾水。3. FFT加速的数学核心从卷积到DFT3.1 Toeplitz乘法本质是一次卷积先把Toeplitz矩阵乘向量的公式写出来。假设T是n乘n的Toeplitz矩阵参数由第一行r [t0, t1, ..., t_{n-1}]和第一列c [t0, t_{-1}, ..., t_{-(n-1)}]确定。计算y T乘x展开第i行y[i] Σ_{j0}^{n-1} t_{i-j} x[j]这边下标i-j可以取负数说明t是一个定义在整数上的序列。而y[i]这个式子本质上就是在做t和x的离散卷积。唯一的区别是离散卷积通常从某个起始下标开始但数学形式和这里完全一致。这个发现太关键了。卷积在时域是O(n^2)的运算但卷积定理说两个序列的循环卷积可以通过各自的离散傅里叶变换在频域做逐点乘积来实现。FFT把DFT的复杂度从O(n^2)降到O(n log n)于是原本平方级的矩阵乘向量运算就被降下来了。3.2 循环矩阵为什么能被FFT对角化要理解为什么卷积能到频域里算需要先看一个更基础的矩阵循环矩阵。循环矩阵的每一行由上一行循环右移一位得到完全由第一行决定第一行有n个自由参数。循环矩阵有个特别好的性质它可以被离散傅里叶矩阵对角化。换句话说任意循环矩阵C都可以写成C F^{-1} D F其中F是DFT矩阵D是对角矩阵对角线上的元素恰好就是C第一行向量的DFT。这意味着计算C乘x时可以先把x做一次FFT再和C第一行的FFT逐点相乘最后做一次逆FFT就得到结果。这个性质一点也不神秘它就是对卷积定理的矩阵语言表述。循环矩阵乘向量等价于循环卷积而DFT把循环卷积变成频域乘法。因为FFT计算DFT只需要O(n log n)时间所以循环矩阵乘向量的代价也从O(n^2)降到了O(n log n)。一个直观的类比循环矩阵和DFT的关系就像是“对角矩阵和普通乘法”的关系。对角矩阵乘向量是逐点相乘快得不得了。循环矩阵经过坐标变换换到频域之后同样变成逐点相乘。Toeplitz矩阵本身做不到这一点因为它的行不是循环的而是平移的。3.3 嵌入法把Toeplitz缝进循环矩阵Toeplitz矩阵不是循环矩阵怎么办答案很朴素把Toeplitz矩阵“放大”成一个更大的循环矩阵让原来的Toeplitz块恰好在它的左上角然后对那个大循环矩阵做FFT加速。具体操作是这样的。假设原始Toeplitz矩阵是n乘n我们构造一个2n乘2n的循环矩阵C只要确定它第一行的2n个元素。第一行的构造规则是先放Toeplitz矩阵的第一列t0, t_{-1}, t_{-2}, ..., t_{-(n-1)}然后放一个0作为隔断再放Toeplitz矩阵第一行的除了t0之外的部分的反向排列t_{n-1}, t_{n-2}, ..., t_1。也就是说循环矩阵第一行是[t0, t_{-1}, ..., t_{-(n-1)}, 0, t_{n-1}, ..., t_1]你可以验证一下这个循环矩阵左上角的n乘n块每个位置的元素T[i, j] t_{i-j}恰好就是原Toeplitz矩阵。这块嵌入是精确的没有近似。计算时把x拼成2n长度的向量前n个位置放x后n个位置补零然后用循环矩阵乘这个补零向量。由于补零的存在循环移位不会把尾部元素卷回到前n个位置去干扰正确结果取结果向量前n个元素就正好是T乘x。3.4 一个具体的嵌入示例我用一个4乘4的例子把嵌入过程走一遍免得光看公式容易晕。假设Toeplitz矩阵由这些参数决定t0 1t_{-1} 2t_{-2} 3t_{-3} 4t1 -1t2 -2t3 -3。那么原矩阵T等于1 2 3 4 -1 1 2 3 -2 -1 1 2 -3 -2 -1 1按照上面的规则嵌入的8乘8循环矩阵第一行是[1, 2, 3, 4, 0, -3, -2, -1]这8个元素就定义了整个循环矩阵。你随便取一个位置验证比如循环矩阵第2行第1列按照循环矩阵行与行右移的关系它的值是第一行的第7个元素也就是-2和原Toeplitz矩阵T[1, 0] -2完全一致。计算时x [x0, x1, x2, x3]补零为[x0, x1, x2, x3, 0, 0, 0, 0]对这个向量做FFT再和第一行向量的FFT逐点相乘逆FFT后取前4个元素。结果和直接做4乘4矩阵乘法一模一样。这里要特别强调补零的必要性。如果不补零直接把x塞进循环矩阵乘向量循环移位会把x的尾部数据卷到前面导致结果在边界位置出现错误叠加。这个错误在FFT加速中是最高频的坑后面会单独展开。3.5 为什么FFT长度要取2次幂理论上循环矩阵的尺寸取m n - 1就足够容纳完整的线性卷积结果了也就是方阵情形下取2n - 1。但实际工程中几乎都取不小于m n - 1的最小二次幂也就是方阵时取2n。原因是FFT算法对2次幂长度的效率远高于普通长度教科书里的蝶形运算、基2算法、SIMD优化都建立在2次幂上。虽然现代成熟FFT库对3、5、7这些小因子长度也支持得很好但遇到任意长度时处理逻辑会复杂一些性能波动也更大。取2n的代价仅仅是多算一点点FFT长度换来的是稳定可预期的速度和更简单的代码。这就是“工程惯例”出现的理由。4. 完整实现和实测效果4.1 一个可直接使用的Python函数下面这段代码是我在实际项目里用的一个简化版本直接解决Toeplitz矩阵乘稠密矩阵的问题。T由col和row两个向量定义col是第一列row是第一行X是右侧矩阵每一列是一个要计算的右端向量。函数返回T乘X。import numpy as np def toeplitz_matmul_fft(col, row, X): col: 第一列长度 mcol[0] row[0] row: 第一行长度 n X: 形状 (n, k) 的稠密矩阵每列是一个右端向量 返回: T X形状 (m, k) m, n len(col), len(row) k X.shape[1] # 选择一个不小于 m n - 1 的二次幂作为 FFT 长度 L 1 while L m n - 1: L 1 # 构造嵌入循环矩阵的第一行 c_first np.concatenate([col, [0.0], row[1:][::-1]]) # 频域中的循环矩阵特征向量 c_f np.fft.rfft(c_first, nL) # 右端向量补零并按列批量做 FFT X_pad np.vstack([X, np.zeros((L - n, k))]) X_f np.fft.rfft(X_pad, axis0) # 频域逐点相乘 Y_f c_f[:, None] * X_f # 逆变换取前 m 行 Y_pad np.fft.irfft(Y_f, nL, axis0) return Y_pad[:m, :]几个细节想特别说一下。第一构造c_first时我用的是row[1:][::-1]也就是row的第二个元素到最后一个元素反向排列。这段代码看似简单但在边界上非常容易写错我建议写完之后用一个很小的随机矩阵对照朴素实现验证一遍。第二代码用的是rfft而不是fft。因为Toeplitz矩阵的元素和右侧向量都是实数rfft只对实数做半谱FFT比复数全谱FFT几乎快一倍。很多人第一次写加速代码时习惯性用np.fft.fft白白浪费一半性能。第三这里一次就把所有右端向量批量处理了。X有多少列就同时算多少列的结果。批处理FFT的cache友好程度远高于循环单列调用所以在矩阵乘法场景下强烈建议这种写法。4.2 实测数据与复杂度对比我用相对时间列一张表展现朴素实现和FFT实现在不同规模下的差距。相对时间以n1024时朴素实现为基准1。矩阵规模n朴素实现相对时间FFT加速相对时间加速比10241.000.16约6倍409616.000.76约21倍16384256.003.51约73倍655364096.0015.82约259倍这张表能说明两个问题。一个是朴素实现的时间每扩大4倍规模就扩大16倍因为复杂度是平方级。另一个是FFT方法的时间每扩大4倍规模只增加约4倍多一点因为O(n log n)里那个log n的因素很弱。n到了六万这个量级加速比已经奔着几百倍去了这种差距足以决定一个算法能不能部署到实际系统。需要注意如果n很小比如只有几十FFT方法反而不如朴素方法。FFT需要分配内存、做复数运算、还要处理补零和逆变换常数因子比较大。工程上我通常会在n小于256时直接走朴素分支大于这个阈值再走FFT分支。所有算法库都这么干不是理论不成立而是实践里要按规模选路。4.3 验证与测试的一个小习惯说到验证我必须强烈建议养成一个习惯保留一个用行观点写的朴素函数专门用来做单元测试。我自己的代码仓库里常年放着一个toeplitz_matmul_naive三层循环一行内积一写慢是慢了一点但正确性一目了然。每次改完FFT版本随机生成一组矩阵跑一遍对比最大误差。误差应该在10的负12次方量级以下。这看起来很简单但能救命的场景比你想象的多。有一次我优化一个分块Toeplitz算法改完觉得逻辑没问题结果和朴素版本对比边缘行误差到了10的负4次方。追了两天发现是FFT长度选小了循环卷积混叠进了原来的边界元素。要不是有那个朴素的对照函数这个bug大概率会一路带进生产代码排查成本高到无法想象。用行观点写验证代码用列观点设计算法这是我多年实践下来觉得最顺手的组合。4.4 更进一步的Toeplitz乘Toeplitz有人可能会问如果两个因子都是Toeplitz矩阵能不能直接用FFT加速这里要泼一盆冷水两个Toeplitz矩阵相乘的结果通常不再是Toeplitz矩阵。它不再具备沿对角线恒定的结构所以不能直接把两个Toeplitz矩阵都转成循环矩阵去相乘。但这不代表没有优化空间。如果只是要计算T1乘T2可以把T2看成很多列向量用前文的批量FFT方法把T1依次作用于T2的每一列复杂度是O(n^2 log n)。由于最终结果本身就是n乘n稠密矩阵至少要花O(n^2)时间输出所以这个复杂度已经贴着输出规模的下界工程上够用。如果你想在O(n^2)以下快速获得两个Toeplitz矩阵的乘积那需要进入“位移秩”和“超快速Toeplitz算法”的领域。这类算法利用了Toeplitz矩阵位移秩低的特点能在O(n log^2 n)甚至更低的复杂度下完成运算但实现复杂度极高一般只在超大超特殊规模下才划算。我在实际项目里几乎不用多数场景用批量FFT方法就够了。5. 应用场景、常见问题和个人经验5.1 用得上这类加速的地方Toeplitz矩阵的FFT加速不是花架子它在很多真实系统里直接决定性能上限。第一个场景是迭代求解线性方程组。很多大规模线性问题里系数矩阵是Toeplitz或分块Toeplitz直接用稠密LU分解完全不可行大家都会转成CG或GMRES这类迭代法。迭代法的核心代价几乎全部集中在每一步的矩阵乘向量也就是matvec。把这个matvec换成FFT加速版本整个求解器立马从“算不动”变成“能算”。如果再配合循环预条件子收敛速度还能进一步改善这是Toeplitz迭代法的经典套路。第二个场景是卷积神经网络的推理和训练。卷积层在网络里执行的操作就是输入特征图和卷积核的二维卷积矩阵化之后是块Toeplitz矩阵乘向量。虽然现在很多框架直接用im2col加GEMM或者Winograd加速但在某些定制化部署场景里FFT路线依然有不可替代的竞争力尤其是卷积核尺寸较大的时候。第三个场景是时间序列分析。ARMA模型、维纳滤波、卡尔曼平滑等一堆经典算法中自相关矩阵天然是Toeplitz。金融数据、脑电信号、地震波形凡是等间隔采样且假设平稳的信号最终都会碰到Toeplitz结构。这时候FFT加速就是标配工具。第四是图像和信号的反卷积问题。点扩散函数模糊成像的矩阵形式就是Toeplitz反卷积本质上要做Toeplitz矩阵的逆运算。即使不直接求逆通常也要在迭代里反复使用Toeplitz乘向量FFT加速几乎绕不开。5.2 常见问题速查表与排查思路我把自己踩过的坑和周围同事踩过的坑整理成一张表按出现频率排序。问题原因解决思路结果在边界位置有明显错误的“鬼影”FFt长度不够或没补零循环卷积发生混叠确保FFT长度不小于mn-1右端向量补零到L用复数FFT实现速度不达预期实信号用了fft而不是rfft换成rfft和irfft几乎白拿一半性能小规模矩阵反而比朴素实现慢FFT常数因子大在n小于阈值时回退到朴素分支大规模矩阵很快但数值误差偏大病态矩阵会放大任何方法的误差使用更高精度或迭代精化普通良态矩阵误差在1e-12量级是正常的内存占用高显式构造了循环矩阵或Toeplitz矩阵只保留第一行第一列FFT直接吃向量两个Toeplitz矩阵相乘结果不对错误地认为乘积仍是Toeplitz确认乘积不再是Toeplitz按列视角批量FFT与scipy.linalg.toeplitz生成的矩阵对比不一致col和row的约定或索引方向不对用随机小矩阵对照朴素实现逐元素检查排查这类问题有个固定套路先固定一个很小的规模比如n8用行观点朴素实现做基准打印中间向量逐步对照。确认小规模正确了再逐步扩大。如果你跳过小规模直接上大矩阵很多结构性的错误会被数字的随机性掩盖很难定位。还有一点要提醒FFT里的频率顺序、归一化约定不同库之间可能有差异。如果你把某个语言的实现移植到另一个语言不要默认行为一致一定要用已知向量做冒烟测试。Numpy的np.fft约定和MATLAB是一致的但和某些底层C库相比归一化方式未必相同别在移植时掉以轻心。5.3 一些写在最后的心得回顾这个题目我觉得最值得记住的不是FFT那几行代码而是理解矩阵结构这件事本身的价值。矩阵乘法行观点和列观点并不是考试知识点它们是帮你识别“这个矩阵能不能被特殊方法加速”的思维工具。从Toeplitz矩阵的平移对称性看到卷积从卷积看到FFT这条因果链环环相扣缺一个视角都会觉得FFT加速凭空冒出来。我在实际工程里最大的体会是不要一上来就写FFT版本先写对的再写快的。先用行观点把朴素版本写出来测试通过再基于列观点做FFT重构然后对比验证。这样每一步的思维方式都是清晰的代码也经得起推敲。这个套路听起来慢实际上是把调试成本前置总体时间反而是省下来的。最后再分享一个小技巧当你把Toeplitz乘向量放进CG或GMRES这类迭代法时记得把FFT的plan或预计算好的频域向量缓存下来。Toeplitz矩阵的特征向量c_f在整个迭代过程中是不变的只需要算一次。我见过不少人每次matvec都重新算一遍c_f白白浪费了两次FFT的时间性能直接打对折。这种细节看似小在实际的大规模迭代里可能就是几十倍的耗时差距。