简介围绕 Toeplitz 线性系统快速求解的一份理论讲义 PDF面向计算机科学、数学研究人员及相关领域学者与技术开发者帮助读者在谱估计、线性预测、自回归滤波器设计与纠错码等场景中高效解算大规模方程组。内容聚焦矩阵结构特性的挖掘系统梳理了对称 Toeplitz 情形的 Levinson 与 Durbin 迭代算法、面向一般非对称矩阵的 Trench 算法并针对不可逆 Toeplitz 矩阵讨论基于 Euclid 算法的改进方案对 Berlekamp-Massey 算法及其加速版本也有专门解析说明如何把计算复杂度由传统的 n² 量级压缩至约 n log n log(log n)显著提升运算效率。压缩包内为 1 个 pdf 文件约 351KB篇幅紧凑、推导与结论并重算法脉络清晰便于按章查阅。目前已有 187 人学习适合具备线性代数与信号处理基础、需要深入理解快速求解思路的读者研读。1. 为什么 Toeplitz 系统值得单写一套快速解算法一个 N×N 矩阵如果每条对角线上的元素都取同一个值它只需要 2N-1 个数就能完整描述这就是 Toeplitz 矩阵。信号处理里的自相关矩阵、线性预测编码的 Normal 方程、卷积反演、谱估计最后都会落到解一个 Toeplitz 线性系统 Txb。麻烦在于 T 是稠密的用通用 LU 硬解要 O(N³)N 到几千时在嵌入式 ARM 上基本做不了实时。结构红利恰恰在这里Levinson-Durbin 递推把它压到 O(N²)分治加 FFT 压到 O(N log²N)循环预条件共轭梯度每次迭代只要 O(N log N)。同一个问题四条路常数因子、数值鲁棒性、能否向量化完全不一样。做信号处理、数值计算和嵌入式 DSP 部署的工程师都值得把这套算法从推导到落地完整走一遍。2. Levinson-Durbin 递推把 O(N³) 消元压到 O(N²)2.1 结构红利与四种快速解法的选型边界Toeplitz 矩阵 T 满足 T[i][j] t[i-j]只用到对称情形时T[i][j] t[|i-j|]。整块矩阵完全由 t[-(N-1)] 到 t[N-1] 这 2N-1 个数决定存储从 N² 塌缩到 2N-1。更重要的是所有基于分块消去加边界增量更新的算法都能在 O(N²) 内走完不需要真的去做列主元消元。同样是一个 Toeplitz 系统选哪条路线取决于 T 是否对称正定、N 的量级、右端项是否要反复求解比如多个 b 复用同一个 T。方法时间复杂度适用条件数值鲁棒性可向量化Levinson-DurbinO(N²)对称正定或强对角占优中等N 大时误差累积中Schur 算法O(N²)一般 Toeplitz含非对称优于 Levinson高分治 FFTO(N log²N)N 大、内存充足好高PCG circulant 预条件O(iter · N log N)对称正定、可快速乘向量好迭代次数依赖预条件子高选型的经验是N 小于 500 且 T 对称正定直接上 Levinson实现最短N 在 10³ 到 10⁴ 之间且只解一两次分治 FFT 的裸递归代码量偏大可以考虑 Schur 加上内层 BLASN 再大、或者 T 是病态自相关矩阵PCG 加循环预条件是最可靠的选择代价是要写矩阵向量乘和预条件子迭代次数还不确定。嵌入式部署里如果 N 只有几十LPC 的阶数通常 10 到 20Levinson 几乎是唯一现实的选择循环体短能塞进 I-cache。2.2 从 Yule-Walker 方程到 Durbin 反射系数递推设 r [r_0, r_1, ..., r_{N-1}] 是 T 的第一列T_n 表示它的 n 阶主子矩阵。AR(p) 模型的 Yule-Walker 方程是sum_{i0}^{p} a_i * r[k-i] 0, k 1, ..., p sum_{i0}^{p} a_i * r[i] sigma^2把系数排成 A^(p) [1, a_1^(p), ..., a_p^(p)]。注意上标表示它属于 p 阶模型a_1^(p) 会随阶数变化这一点容易看错。Durbin 递推每一步用前一阶的系数和反射系数 K_p 生成下一阶K_p -(1 / E_{p-1}) * sum_{j0}^{p-1} a_j^(p-1) * r[p-j] a_j^(p) a_j^(p-1) K_p * a_{p-j}^(p-1), 1 j p-1 a_p^(p) K_p E_p E_{p-1} * (1 - K_p^2)E_p 是 p 阶预测误差也是 T 主子矩阵的行列式比。条件 |K_p| 1 等价于 T_{p1} 正定这个判据比去算特征值便宜得多。到这里拿到的只是预测系数要真正解 Txb 还需要把解向量 x 逐维扩展。把 T_{k1} 按 Schur 补写成分块形式 T_{k1} [[T_k, g_k], [g_k^T, r_0]]其中 g_k [r_k, r_{k-1}, ..., r_1]^T。记 y^(k) T_k^{-1} g_k从 Durbin 的系数可以直接读出 y^(k) -[a_k^(k), a_{k-1}^(k), ..., a_1^(k)]这正是预测系数能顺带解线性系统的原因。Schur 补 s_k r_0 - g_k^T y^(k) 恰好等于 E_k于是从 x^(k) 扩到 x^(k1) 的更新是alpha (b[k] - g_k^T x^(k)) / E_k x^(k1)[j] x^(k)[j] alpha * a_{k-j}^(k), j 0, ..., k-1 x^(k1)[k] alpha一个循环里同时更新 a、E 和 x整个算法只有 O(N²) 次乘加。2.3 可复现的最小 Python 实现import numpy as np def levinson_sym(r, b): 求解对称正定 Toeplitz 系统 T x b其中 T[i, j] r[|i - j|]。 r : 长度 n 的数组r[0] 是对角元 b : 长度 n 的右端向量 返回 x : 长度 n 的解向量 n len(b) x np.zeros(n) a np.zeros(n 1) # 当前阶预测系数 a^(k)a[0] 恒为 1 a[0] 1.0 x[0] b[0] / r[0] E float(r[0]) # 0 阶预测误差 E_0 for k in range(1, n): # 1) Durbin 反射系数K_k -(sum_{jk} a_j * r[k-j]) / E_{k-1} num 0.0 for j in range(k): num a[j] * r[k - j] K -num / E # 2) 交叉更新预测系数a_new[j] a[j] K * a[k-j] a_new np.zeros(k 1) a_new[0] 1.0 a_new[k] K for j in range(1, k): a_new[j] a[j] K * a[k - j] a[:k 1] a_new # 3) 预测误差迭代E_k E_{k-1} * (1 - K^2) E * (1.0 - K * K) # 4) 用 Schur 补把 x 从 k 维扩到 k1 维 s 0.0 for j in range(k): s r[k - j] * x[j] alpha (b[k] - s) / E for j in range(k): x[j] alpha * a[k - j] x[k] alpha return x if __name__ __main__: n 128 rng np.random.default_rng(0) r np.exp(-0.1 * np.arange(n)) # 指数衰减自相关保证正定 b rng.standard_normal(n) T np.array([[r[abs(i - j)] for j in range(n)] for i in range(n)]) x_ref np.linalg.solve(T, b) x_lin levinson_sym(r, b) rel np.linalg.norm(x_lin - x_ref) / np.linalg.norm(x_ref) print(relative error:, rel)代码分四块看。第一块的 num 是 Yule-Walker 方程残差 sum_{jk} a_j r[k-j]除以 E_{k-1} 得到反射系数 K_k。第二块做系数交叉更新a_new[j] 依赖上一轮的 a[j] 和 a[k-j]所以必须写进临时数组再拷回原地更新会把后半段污染。第三块 E 的更新是 Durbin 的核心恒等式E 单调递减永远不会变号。第四块的 s 是 g_k 与当前解的內积alpha 是 Schur 消元得到的新分量a[k-j] 来自第二块算出的系数把 x 前 k 维做一次线性修正后补上尾元素。整段代码只有两个长度为 n 的数组空间 O(N)。np.linalg.solve 只用来做正确性对照生产代码里不该出现。2.4 r[0] 归一化、反射系数门限与对角加载Levinson 的两类失败几乎都是数值问题。一是 r 没归一化当 r[0] 很小时例如从归一化信号估出来的自相关K 的分母 E 会先小后大溢出和舍入误差同时出现。工程上先把 r 整体除以 r[0]让 r[0] 1最后再按比例还原 x。二是自相关来自有限长观测、被噪声污染导致 |K_p| 越过 1此时 T 不再正定Levinson 会在某一步发散。参数建议取值说明r 归一化r / r[0]保持 E 在 1 附近抑制溢出反射系数门限abs(K_k) 1 - 1e-6越界即判非正定回退到加对角加载对角加载量 lambda1e-8 ~ 1e-4 倍 r[0]r[0] lambda把 T 拉回正定预测误差停机E_k 1e-12 倍 E_0再迭代下去只剩舍入噪声阶数上限 p信号带宽 / 采样率的经验值LPC 里常取 10 ~ 20对角加载是对付病态自相关最省事的一招它相当于给 T 的每条对角元加一个小常数把最小特征值抬到 lambda 以上代价是解会略微平滑。加了之后反射系数的绝对值会压到 1 以下Levinson 又能顺利跑完。如果在嵌入式上做定点就把 r 先缩放成 Q15 或 Q31 再跑同一套递推只是每个累加器要多留几位保护位否则第 2 块的交叉更新很容易在中途溢出。3. 分治 FFT 与循环预条件把快速解推到准线性3.1 循环嵌入用 FFT 做 O(N log N) 的 Toeplitz 乘向量Toeplitz 矩阵乘向量可以用卷积实现。把 T 嵌入一个更大的循环矩阵循环卷积的前 n 项恰好等于 Tx。做法是取 m 为不小于 2n-1 的 2 的幂构造长度 m 的向量 c前 n 项放 r然后从尾部往前填 r[1:] 的反序。用 FFT 对角化这个循环矩阵一次乘向量就变成两次正变换加一次逆变换。import numpy as np def toeplitz_matvec_fast(r, x, mNone): 计算 y T xT 是由 r 定义的对称 Toeplitz 矩阵复杂度 O(n log n)。 r : 长度 nT[i, j] r[|i - j|] x : 长度 n m : 循环嵌入长度默认取 2n-1 的最小 2 的幂 n len(r) if m is None: m 1 while m 2 * n - 1: m 1 c np.zeros(m) c[:n] r c[m - n 1:] r[1:][::-1] # 嵌入 Toeplitz 的上三角部分 xp np.zeros(m) xp[:n] x y np.fft.ifft(np.fft.fft(c) * np.fft.fft(xp)).real return y[:n]索引 (i-j) mod m 落在 [0, n-1] 时取 r[i-j]落在 [m-n1, m-1] 时取 r[j-i]正好覆盖 T[i][j] 的两种情形。m 取 2 的幂是为了 FFT 长度对齐不取也能跑但速度会掉。numpy 的 ifft 会引入约 1e-15 的虚部残差取 .real 足够。嵌入长度要比 2n-1 大否则循环回绕会污染前 n 项。3.2 CG 加 FFT 矩阵向量乘的完整求解流程有了快速乘向量就可以把 Txb 交给共轭梯度迭代每次迭代只调一次 toeplitz_matvec_fast。def cg_toeplitz(r, b, tol1e-10, maxit1000): 用 CG 求解对称正定 Toeplitz 系统矩阵向量乘走 FFT。 n len(b) x np.zeros(n) r_vec b - toeplitz_matvec_fast(r, x) # 初始残差 p r_vec.copy() rs_old r_vec r_vec for it in range(maxit): Ap toeplitz_matvec_fast(r, p) alpha rs_old / (p Ap) # 精确线搜索步长 x alpha * p r_vec - alpha * Ap rs_new r_vec r_vec if np.sqrt(rs_new) tol * np.linalg.norm(b): return x, it 1 p r_vec (rs_new / rs_old) * p # Fletcher-Reeves 更新 rs_old rs_new return x, maxitr_vec 是残差p 是搜索方向rs_old 和 rs_new 是前后两次残差平方。alpha 取 rs_old / (p·Ap) 是共轭梯度的精确线搜索p 的更新系数取 rs_new / rs_old。没有预条件子时迭代次数大致服从 sqrt(kappa(T))其中 kappa 是条件数。Toeplitz 系统如果来自真实信号的自相关kappa 上千很常见迭代次数会破百到这一步就必须上预条件子。3.3 circulant 预条件子怎么补循环矩阵可以被 FFT 直接对角化所以循环矩阵求逆是 O(N log N) 的。把 T 用一个循环矩阵 C 去逼近M C^{-1} 当预条件子每次迭代多两次 FFT但迭代次数能压到个位数到十几次。预条件子第一列构造适用场景迭代次数典型Strangc [r_0, ..., r_{n/2}, 0, ..., 0, r_{n/2}, ..., r_1]自相关衰减快5 ~ 20Chanc[j] r_j - r_{n-j}首项取 r_0一般对称 Toeplitz5 ~ 15T. Chan 最优最小化循环矩阵与 T 的 Frobenius 距离需要预先估计3 ~ 10不加预条件无条件数小于 100可能上百Strang 预条件子的思路是把 r 的前半段保留、后半段清零后嵌入循环结构对短记忆信号很合适。Chan 预条件子把首行和末行的差当作循环生成元对称性更好工程上更常用。两个都不需要构造 n×n 矩阵只需要长度为 n 的向量加 FFT所以内存开销和 Levinson 一个量级。判断该不该上的标准很简单CG 残差在 50 次迭代内降不到 1e-8就直接换成 Chan。3.4 复杂度、精度与内存的对照路线时间空间相对残差N4096 量级Levinson-DurbinO(N²)O(N)1e-10 以下良态分治 FFTO(N log²N)O(N log N)1e-12 量级CG无预条件O(iter · N log N)O(N)依赖 kappa(T)CG Chan 预条件O(iter · N log N)O(N)1e-10 以下N 在 1000 以内Levinson 的绝对时间往往还更短因为常数小、没有 FFT 的固定开销。N 过 10⁴ 之后O(N²) 的乘加数就到 10⁸ 以上FFT 路线开始反超。内存上 Levinson 只需要两个 n 长数组是嵌入式上几乎无内存压力的选择CG 路线要多留几个 n 长向量给 p、Ap 和残差。4. 把 Toeplitz 快速解部署到嵌入式 Linux4.1 设备树里描述一块 Toeplitz 加速器常见的做法是把 Levinson 或 Schur 的核心循环塞进一块 FPGA 或 DSP 加速器主 CPU 通过 AXI 或 SPI 下发 r 和 b加速器算完回读 x。设备树里要描述寄存器地址、中断号、时钟和 DMA 通道下面是一个最小示例。/ { toeplitz_accel: toeplitz43c00000 { compatible acme,toeplitz-accel-1.0; reg 0x43c00000 0x100 p a hrefhttps://download.csdn.net/download/huanghm88/90473436 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p Toeplitz 快速解：Levinson、分治 FFT 与预条件 CG
简介围绕 Toeplitz 线性系统快速求解的一份理论讲义 PDF面向计算机科学、数学研究人员及相关领域学者与技术开发者帮助读者在谱估计、线性预测、自回归滤波器设计与纠错码等场景中高效解算大规模方程组。内容聚焦矩阵结构特性的挖掘系统梳理了对称 Toeplitz 情形的 Levinson 与 Durbin 迭代算法、面向一般非对称矩阵的 Trench 算法并针对不可逆 Toeplitz 矩阵讨论基于 Euclid 算法的改进方案对 Berlekamp-Massey 算法及其加速版本也有专门解析说明如何把计算复杂度由传统的 n² 量级压缩至约 n log n log(log n)显著提升运算效率。压缩包内为 1 个 pdf 文件约 351KB篇幅紧凑、推导与结论并重算法脉络清晰便于按章查阅。目前已有 187 人学习适合具备线性代数与信号处理基础、需要深入理解快速求解思路的读者研读。1. 为什么 Toeplitz 系统值得单写一套快速解算法一个 N×N 矩阵如果每条对角线上的元素都取同一个值它只需要 2N-1 个数就能完整描述这就是 Toeplitz 矩阵。信号处理里的自相关矩阵、线性预测编码的 Normal 方程、卷积反演、谱估计最后都会落到解一个 Toeplitz 线性系统 Txb。麻烦在于 T 是稠密的用通用 LU 硬解要 O(N³)N 到几千时在嵌入式 ARM 上基本做不了实时。结构红利恰恰在这里Levinson-Durbin 递推把它压到 O(N²)分治加 FFT 压到 O(N log²N)循环预条件共轭梯度每次迭代只要 O(N log N)。同一个问题四条路常数因子、数值鲁棒性、能否向量化完全不一样。做信号处理、数值计算和嵌入式 DSP 部署的工程师都值得把这套算法从推导到落地完整走一遍。2. Levinson-Durbin 递推把 O(N³) 消元压到 O(N²)2.1 结构红利与四种快速解法的选型边界Toeplitz 矩阵 T 满足 T[i][j] t[i-j]只用到对称情形时T[i][j] t[|i-j|]。整块矩阵完全由 t[-(N-1)] 到 t[N-1] 这 2N-1 个数决定存储从 N² 塌缩到 2N-1。更重要的是所有基于分块消去加边界增量更新的算法都能在 O(N²) 内走完不需要真的去做列主元消元。同样是一个 Toeplitz 系统选哪条路线取决于 T 是否对称正定、N 的量级、右端项是否要反复求解比如多个 b 复用同一个 T。方法时间复杂度适用条件数值鲁棒性可向量化Levinson-DurbinO(N²)对称正定或强对角占优中等N 大时误差累积中Schur 算法O(N²)一般 Toeplitz含非对称优于 Levinson高分治 FFTO(N log²N)N 大、内存充足好高PCG circulant 预条件O(iter · N log N)对称正定、可快速乘向量好迭代次数依赖预条件子高选型的经验是N 小于 500 且 T 对称正定直接上 Levinson实现最短N 在 10³ 到 10⁴ 之间且只解一两次分治 FFT 的裸递归代码量偏大可以考虑 Schur 加上内层 BLASN 再大、或者 T 是病态自相关矩阵PCG 加循环预条件是最可靠的选择代价是要写矩阵向量乘和预条件子迭代次数还不确定。嵌入式部署里如果 N 只有几十LPC 的阶数通常 10 到 20Levinson 几乎是唯一现实的选择循环体短能塞进 I-cache。2.2 从 Yule-Walker 方程到 Durbin 反射系数递推设 r [r_0, r_1, ..., r_{N-1}] 是 T 的第一列T_n 表示它的 n 阶主子矩阵。AR(p) 模型的 Yule-Walker 方程是sum_{i0}^{p} a_i * r[k-i] 0, k 1, ..., p sum_{i0}^{p} a_i * r[i] sigma^2把系数排成 A^(p) [1, a_1^(p), ..., a_p^(p)]。注意上标表示它属于 p 阶模型a_1^(p) 会随阶数变化这一点容易看错。Durbin 递推每一步用前一阶的系数和反射系数 K_p 生成下一阶K_p -(1 / E_{p-1}) * sum_{j0}^{p-1} a_j^(p-1) * r[p-j] a_j^(p) a_j^(p-1) K_p * a_{p-j}^(p-1), 1 j p-1 a_p^(p) K_p E_p E_{p-1} * (1 - K_p^2)E_p 是 p 阶预测误差也是 T 主子矩阵的行列式比。条件 |K_p| 1 等价于 T_{p1} 正定这个判据比去算特征值便宜得多。到这里拿到的只是预测系数要真正解 Txb 还需要把解向量 x 逐维扩展。把 T_{k1} 按 Schur 补写成分块形式 T_{k1} [[T_k, g_k], [g_k^T, r_0]]其中 g_k [r_k, r_{k-1}, ..., r_1]^T。记 y^(k) T_k^{-1} g_k从 Durbin 的系数可以直接读出 y^(k) -[a_k^(k), a_{k-1}^(k), ..., a_1^(k)]这正是预测系数能顺带解线性系统的原因。Schur 补 s_k r_0 - g_k^T y^(k) 恰好等于 E_k于是从 x^(k) 扩到 x^(k1) 的更新是alpha (b[k] - g_k^T x^(k)) / E_k x^(k1)[j] x^(k)[j] alpha * a_{k-j}^(k), j 0, ..., k-1 x^(k1)[k] alpha一个循环里同时更新 a、E 和 x整个算法只有 O(N²) 次乘加。2.3 可复现的最小 Python 实现import numpy as np def levinson_sym(r, b): 求解对称正定 Toeplitz 系统 T x b其中 T[i, j] r[|i - j|]。 r : 长度 n 的数组r[0] 是对角元 b : 长度 n 的右端向量 返回 x : 长度 n 的解向量 n len(b) x np.zeros(n) a np.zeros(n 1) # 当前阶预测系数 a^(k)a[0] 恒为 1 a[0] 1.0 x[0] b[0] / r[0] E float(r[0]) # 0 阶预测误差 E_0 for k in range(1, n): # 1) Durbin 反射系数K_k -(sum_{jk} a_j * r[k-j]) / E_{k-1} num 0.0 for j in range(k): num a[j] * r[k - j] K -num / E # 2) 交叉更新预测系数a_new[j] a[j] K * a[k-j] a_new np.zeros(k 1) a_new[0] 1.0 a_new[k] K for j in range(1, k): a_new[j] a[j] K * a[k - j] a[:k 1] a_new # 3) 预测误差迭代E_k E_{k-1} * (1 - K^2) E * (1.0 - K * K) # 4) 用 Schur 补把 x 从 k 维扩到 k1 维 s 0.0 for j in range(k): s r[k - j] * x[j] alpha (b[k] - s) / E for j in range(k): x[j] alpha * a[k - j] x[k] alpha return x if __name__ __main__: n 128 rng np.random.default_rng(0) r np.exp(-0.1 * np.arange(n)) # 指数衰减自相关保证正定 b rng.standard_normal(n) T np.array([[r[abs(i - j)] for j in range(n)] for i in range(n)]) x_ref np.linalg.solve(T, b) x_lin levinson_sym(r, b) rel np.linalg.norm(x_lin - x_ref) / np.linalg.norm(x_ref) print(relative error:, rel)代码分四块看。第一块的 num 是 Yule-Walker 方程残差 sum_{jk} a_j r[k-j]除以 E_{k-1} 得到反射系数 K_k。第二块做系数交叉更新a_new[j] 依赖上一轮的 a[j] 和 a[k-j]所以必须写进临时数组再拷回原地更新会把后半段污染。第三块 E 的更新是 Durbin 的核心恒等式E 单调递减永远不会变号。第四块的 s 是 g_k 与当前解的內积alpha 是 Schur 消元得到的新分量a[k-j] 来自第二块算出的系数把 x 前 k 维做一次线性修正后补上尾元素。整段代码只有两个长度为 n 的数组空间 O(N)。np.linalg.solve 只用来做正确性对照生产代码里不该出现。2.4 r[0] 归一化、反射系数门限与对角加载Levinson 的两类失败几乎都是数值问题。一是 r 没归一化当 r[0] 很小时例如从归一化信号估出来的自相关K 的分母 E 会先小后大溢出和舍入误差同时出现。工程上先把 r 整体除以 r[0]让 r[0] 1最后再按比例还原 x。二是自相关来自有限长观测、被噪声污染导致 |K_p| 越过 1此时 T 不再正定Levinson 会在某一步发散。参数建议取值说明r 归一化r / r[0]保持 E 在 1 附近抑制溢出反射系数门限abs(K_k) 1 - 1e-6越界即判非正定回退到加对角加载对角加载量 lambda1e-8 ~ 1e-4 倍 r[0]r[0] lambda把 T 拉回正定预测误差停机E_k 1e-12 倍 E_0再迭代下去只剩舍入噪声阶数上限 p信号带宽 / 采样率的经验值LPC 里常取 10 ~ 20对角加载是对付病态自相关最省事的一招它相当于给 T 的每条对角元加一个小常数把最小特征值抬到 lambda 以上代价是解会略微平滑。加了之后反射系数的绝对值会压到 1 以下Levinson 又能顺利跑完。如果在嵌入式上做定点就把 r 先缩放成 Q15 或 Q31 再跑同一套递推只是每个累加器要多留几位保护位否则第 2 块的交叉更新很容易在中途溢出。3. 分治 FFT 与循环预条件把快速解推到准线性3.1 循环嵌入用 FFT 做 O(N log N) 的 Toeplitz 乘向量Toeplitz 矩阵乘向量可以用卷积实现。把 T 嵌入一个更大的循环矩阵循环卷积的前 n 项恰好等于 Tx。做法是取 m 为不小于 2n-1 的 2 的幂构造长度 m 的向量 c前 n 项放 r然后从尾部往前填 r[1:] 的反序。用 FFT 对角化这个循环矩阵一次乘向量就变成两次正变换加一次逆变换。import numpy as np def toeplitz_matvec_fast(r, x, mNone): 计算 y T xT 是由 r 定义的对称 Toeplitz 矩阵复杂度 O(n log n)。 r : 长度 nT[i, j] r[|i - j|] x : 长度 n m : 循环嵌入长度默认取 2n-1 的最小 2 的幂 n len(r) if m is None: m 1 while m 2 * n - 1: m 1 c np.zeros(m) c[:n] r c[m - n 1:] r[1:][::-1] # 嵌入 Toeplitz 的上三角部分 xp np.zeros(m) xp[:n] x y np.fft.ifft(np.fft.fft(c) * np.fft.fft(xp)).real return y[:n]索引 (i-j) mod m 落在 [0, n-1] 时取 r[i-j]落在 [m-n1, m-1] 时取 r[j-i]正好覆盖 T[i][j] 的两种情形。m 取 2 的幂是为了 FFT 长度对齐不取也能跑但速度会掉。numpy 的 ifft 会引入约 1e-15 的虚部残差取 .real 足够。嵌入长度要比 2n-1 大否则循环回绕会污染前 n 项。3.2 CG 加 FFT 矩阵向量乘的完整求解流程有了快速乘向量就可以把 Txb 交给共轭梯度迭代每次迭代只调一次 toeplitz_matvec_fast。def cg_toeplitz(r, b, tol1e-10, maxit1000): 用 CG 求解对称正定 Toeplitz 系统矩阵向量乘走 FFT。 n len(b) x np.zeros(n) r_vec b - toeplitz_matvec_fast(r, x) # 初始残差 p r_vec.copy() rs_old r_vec r_vec for it in range(maxit): Ap toeplitz_matvec_fast(r, p) alpha rs_old / (p Ap) # 精确线搜索步长 x alpha * p r_vec - alpha * Ap rs_new r_vec r_vec if np.sqrt(rs_new) tol * np.linalg.norm(b): return x, it 1 p r_vec (rs_new / rs_old) * p # Fletcher-Reeves 更新 rs_old rs_new return x, maxitr_vec 是残差p 是搜索方向rs_old 和 rs_new 是前后两次残差平方。alpha 取 rs_old / (p·Ap) 是共轭梯度的精确线搜索p 的更新系数取 rs_new / rs_old。没有预条件子时迭代次数大致服从 sqrt(kappa(T))其中 kappa 是条件数。Toeplitz 系统如果来自真实信号的自相关kappa 上千很常见迭代次数会破百到这一步就必须上预条件子。3.3 circulant 预条件子怎么补循环矩阵可以被 FFT 直接对角化所以循环矩阵求逆是 O(N log N) 的。把 T 用一个循环矩阵 C 去逼近M C^{-1} 当预条件子每次迭代多两次 FFT但迭代次数能压到个位数到十几次。预条件子第一列构造适用场景迭代次数典型Strangc [r_0, ..., r_{n/2}, 0, ..., 0, r_{n/2}, ..., r_1]自相关衰减快5 ~ 20Chanc[j] r_j - r_{n-j}首项取 r_0一般对称 Toeplitz5 ~ 15T. Chan 最优最小化循环矩阵与 T 的 Frobenius 距离需要预先估计3 ~ 10不加预条件无条件数小于 100可能上百Strang 预条件子的思路是把 r 的前半段保留、后半段清零后嵌入循环结构对短记忆信号很合适。Chan 预条件子把首行和末行的差当作循环生成元对称性更好工程上更常用。两个都不需要构造 n×n 矩阵只需要长度为 n 的向量加 FFT所以内存开销和 Levinson 一个量级。判断该不该上的标准很简单CG 残差在 50 次迭代内降不到 1e-8就直接换成 Chan。3.4 复杂度、精度与内存的对照路线时间空间相对残差N4096 量级Levinson-DurbinO(N²)O(N)1e-10 以下良态分治 FFTO(N log²N)O(N log N)1e-12 量级CG无预条件O(iter · N log N)O(N)依赖 kappa(T)CG Chan 预条件O(iter · N log N)O(N)1e-10 以下N 在 1000 以内Levinson 的绝对时间往往还更短因为常数小、没有 FFT 的固定开销。N 过 10⁴ 之后O(N²) 的乘加数就到 10⁸ 以上FFT 路线开始反超。内存上 Levinson 只需要两个 n 长数组是嵌入式上几乎无内存压力的选择CG 路线要多留几个 n 长向量给 p、Ap 和残差。4. 把 Toeplitz 快速解部署到嵌入式 Linux4.1 设备树里描述一块 Toeplitz 加速器常见的做法是把 Levinson 或 Schur 的核心循环塞进一块 FPGA 或 DSP 加速器主 CPU 通过 AXI 或 SPI 下发 r 和 b加速器算完回读 x。设备树里要描述寄存器地址、中断号、时钟和 DMA 通道下面是一个最小示例。/ { toeplitz_accel: toeplitz43c00000 { compatible acme,toeplitz-accel-1.0; reg 0x43c00000 0x100 p a hrefhttps://download.csdn.net/download/huanghm88/90473436 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p