分块矩阵的快速Fourier变换（FFT）在Toeplitz矩阵向量乘法中的应用
题目描述
考虑一个 \(n \times n\) 的实Toeplitz矩阵 \(T\)（即矩阵的每条对角线上的元素相同，\(T_{i,j} = t_{i-j}\)），需要计算矩阵向量乘法 \(y = T x\)，其中 \(x\) 是给定的向量。直接计算需要 \(O(n^2)\) 次浮点运算。但利用Toeplitz矩阵可嵌入为循环矩阵的性质，并结合快速Fourier变换（FFT），可将计算复杂度降至 \(O(n \log n)\)。本题目将详细讲解如何通过分块矩阵的视角，构造扩展的循环矩阵，并利用FFT高效计算 \(T x\)。
解题过程
步骤1：理解Toeplitz矩阵的结构
一个 \(n \times n\) 的Toeplitz矩阵形式如下：
\[T = \begin{bmatrix} t_0 & t_{-1} & t_{-2} & \cdots & t_{-(n-1)} \\ t_1 & t_0 & t_{-1} & \cdots & t_{-(n-2)} \\ t_2 & t_1 & t_0 & \cdots & t_{-(n-3)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ t_{n-1} & t_{n-2} & t_{n-3} & \cdots & t_0 \end{bmatrix} \]
其中 \(t_k\)（\(k = -(n-1), \dots, n-1\)）是给定的常数。
注意：第一列元素是 \(t_0, t_1, \dots, t_{n-1}\)，第一行元素是 \(t_0, t_{-1}, \dots, t_{-(n-1)}\)。
步骤2：将Toeplitz矩阵嵌入循环矩阵
循环矩阵 \(C\) 的每一行是上一行的循环右移。关键观察：任何 \(n \times n\) 的Toeplitz矩阵可以嵌入一个更大的 \(m \times m\) 循环矩阵中，其中 \(m = 2n\)（或至少 \(m \ge 2n-1\)）。
构造一个 \(m \times m\) 循环矩阵 \(C\)（\(m = 2n\)），其第一列 \(c\) 定义为：
\[c = [t_0, t_1, \dots, t_{n-1}, 0, t_{-(n-1)}, t_{-(n-2)}, \dots, t_{-1}]^T \]
即：
- 前 \(n\) 个元素：\(t_0, t_1, \dots, t_{n-1}\)（对应 \(T\) 的第一列）。
- 第 \(n+1\) 个元素：0（填充项，使总长度为 \(2n\)）。
- 后 \(n-1\) 个元素：\(t_{-(n-1)}, t_{-(n-2)}, \dots, t_{-1}\)（对应 \(T\) 第一行的后 \(n-1\) 个元素，顺序倒置并放在末尾）。
验证：循环矩阵 \(C\) 的左上角 \(n \times n\) 子矩阵恰好等于 \(T\)。
步骤3：利用循环矩阵的傅里叶对角化
循环矩阵 \(C\) 可被离散Fourier变换（DFT）矩阵对角化：
\[C = F_m^{-1} \cdot \text{diag}(\hat{c}) \cdot F_m \]
其中：
- \(F_m\) 是 \(m \times m\) 的DFT矩阵，满足 \((F_m)_{j,k} = \omega_m^{jk} / \sqrt{m}\)（通常归一化因子可调整），且 \(F_m^{-1} = F_m^H\)（酉矩阵）。
- \(\hat{c} = F_m c\) 是向量 \(c\) 的DFT。
- 对角矩阵 \(\text{diag}(\hat{c})\) 的对角线元素就是 \(\hat{c}\) 的分量。
因此，对任意向量 \(z \in \mathbb{C}^m\)，计算 \(C z\) 可转换为三步：
- 计算 \(\hat{z} = F_m z\)（DFT，用FFT实现，复杂度 \(O(m \log m)\)）。
- 计算逐点乘法 \(\hat{y} = \hat{c} \odot \hat{z}\)（\(\odot\) 表示逐元素乘，复杂度 \(O(m)\)）。
- 计算逆DFT：\(y = F_m^{-1} \hat{y}\)（用FFT实现，复杂度 \(O(m \log m)\)）。
总复杂度为 \(O(m \log m) = O(n \log n)\)（因为 \(m = 2n\)）。
步骤4：计算 \(T x\) 的嵌入技巧
目标：计算 \(y = T x\)，其中 \(x \in \mathbb{R}^n\)。
构造扩展向量 \(\tilde{x} \in \mathbb{R}^m\)（\(m = 2n\)）：
\[\tilde{x} = [x_0, x_1, \dots, x_{n-1}, 0, 0, \dots, 0]^T \]
即前 \(n\) 个分量等于 \(x\)，后 \(n\) 个分量补零。
计算 \(\tilde{y} = C \tilde{x}\) 利用步骤3的FFT方法，得到长度为 \(m\) 的向量 \(\tilde{y}\)。
取 \(\tilde{y}\) 的前 \(n\) 个分量，即为 \(y = T x\)。
理由：因为 \(C\) 的左上 \(n \times n\) 子矩阵等于 \(T\)，且 \(\tilde{x}\) 的后 \(n\) 个分量为零，所以 \(C\tilde{x}\) 的前 \(n\) 个分量正好等于 \(T x\)。
步骤5：算法伪代码
- 输入：Toeplitz矩阵的第一列 \(a = [t_0, t_1, \dots, t_{n-1}]^T\) 和第一行（去掉第一个元素）\(b = [t_{-1}, t_{-2}, \dots, t_{-(n-1)}]^T\)，以及向量 \(x \in \mathbb{R}^n\)。
- 令 \(m = 2n\)。
- 构造向量 \(c \in \mathbb{R}^m\)：
\[ c = [a^T, 0, \text{reverse}(b)^T]^T \]
其中 \(\text{reverse}(b)\) 表示将 \(b\) 的元素顺序颠倒（即 \(t_{-(n-1)}, t_{-(n-2)}, \dots, t_{-1}\)）。
4. 构造向量 \(\tilde{x} \in \mathbb{R}^m\)：前 \(n\) 个元素为 \(x\)，后 \(n\) 个元素为0。
5. 计算 \(\hat{c} = \text{FFT}(c)\)，\(\hat{\tilde{x}} = \text{FFT}(\tilde{x})\)。
6. 计算逐点乘积 \(\hat{\tilde{y}} = \hat{c} \odot \hat{\tilde{x}}\)。
7. 计算 \(\tilde{y} = \text{IFFT}(\hat{\tilde{y}})\)。
8. 输出 \(y = \tilde{y}[0 : n-1]\)（取前 \(n\) 个分量）。
步骤6：复杂度与注意事项
- 主要计算量：两次长度为 \(m = 2n\) 的FFT和一次IFFT，每次FFT/IFFT需 \(O(m \log m) = O(n \log n)\) 次运算。逐点乘法为 \(O(n)\)。总复杂度 \(O(n \log n)\)，远优于直接矩阵乘法的 \(O(n^2)\)。
- 若 \(T\) 是对称Toeplitz（即 \(t_{-k} = t_k\)），则构造更简单：\(c = [t_0, t_1, \dots, t_{n-1}, 0, t_{n-1}, \dots, t_1]^T\)。
- 实际实现时，常使用复数FFT，即使输入为实数，可利用实数FFT的优化（如使用
numpy.fft.rfft等）。 - 该方法也可推广到分块Toeplitz矩阵（每个块是Toeplitz），此时可利用二维FFT进一步加速，但原理类似：将分块Toeplitz矩阵嵌入分块循环矩阵，再用二维FFT对角化。
总结
本算法通过嵌入技巧将Toeplitz矩阵向量乘法转化为循环矩阵的乘法，并利用FFT实现快速卷积，将计算复杂度从 \(O(n^2)\) 降至 \(O(n \log n)\)。此方法在信号处理、图像卷积、求解Toeplitz线性系统（如Levinson-Durbin算法预处理）中广泛应用。