分块矩阵的快速Fourier变换加速预处理共轭梯度法求解Toeplitz线性方程组
题目描述
考虑求解大型线性方程组 \(A\mathbf{x} = \mathbf{b}\)，其中 \(A\) 是一个 \(n \times n\) 的对称正定 Toeplitz 矩阵（即矩阵的每条对角线上的元素相同，形式为 \(a_{i,j} = t_{i-j}\)）。当 \(n\) 很大时，直接解法（如 Cholesky 分解）计算复杂度为 \(O(n^3)\)，内存需求为 \(O(n^2)\)，效率低下。本题要求利用 Toeplitz 矩阵的特殊结构，结合快速 Fourier 变换（FFT） 设计一种预处理技术，并嵌入共轭梯度法（CG） 中，实现高效迭代求解。核心思想是：通过循环嵌入将 Toeplitz 矩阵向量乘法复杂度降为 \(O(n \log n)\)，并构造基于循环矩阵的预处理子加速 CG 收敛。
解题过程
步骤1：理解问题结构
Toeplitz 矩阵 \(T_n\) 定义为：
\[T_n = \begin{pmatrix} 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{pmatrix} \]
由于对称正定性，有 \(t_{-k} = t_k\)。直接存储只需 \(2n-1\) 个元素 \(\{t_{-(n-1)}, \dots, t_0, \dots, t_{n-1}\}\)。
步骤2：利用 FFT 加速矩阵向量乘法
关键技巧：将 \(T_n\) 嵌入一个 \(2n \times 2n\) 的循环矩阵 \(C_{2n}\)。循环矩阵完全由其第一列 \(\mathbf{c} = (c_0, c_1, \dots, c_{2n-1})^T\) 决定，且可被 Fourier 矩阵对角化。构造方式：
- 令 \(C_{2n}\) 的第一列为：
\[ \mathbf{c} = (t_0, t_1, \dots, t_{n-1}, 0, t_{-(n-1)}, t_{-(n-2)}, \dots, t_{-1})^T \]
其中 \(t_{-k} = t_k\) 对称时，后半部分即 \((0, t_{n-1}, \dots, t_1)^T\)。
2. 对于任意向量 \(\mathbf{x} \in \mathbb{R}^n\)，补零扩展为 \(\tilde{\mathbf{x}} = (\mathbf{x}^T, \mathbf{0}^T)^T \in \mathbb{R}^{2n}\)。
3. 计算循环矩阵向量积 \(C_{2n} \tilde{\mathbf{x}}\) 等价于 \(\mathbf{c}\) 与 \(\tilde{\mathbf{x}}\) 的循环卷积，可通过 FFT 实现：
\[ C_{2n} \tilde{\mathbf{x}} = \text{IFFT}(\text{FFT}(\mathbf{c}) \odot \text{FFT}(\tilde{\mathbf{x}})) \]
其中 \(\odot\) 表示逐元素相乘。计算复杂度为 \(O(2n \log(2n)) = O(n \log n)\)。
4. 取结果的前 \(n\) 个分量即得 \(T_n \mathbf{x}\)。
步骤3：设计预处理子
共轭梯度法的收敛速度依赖于系数矩阵的条件数。对于 Toeplitz 矩阵，常用循环预处理子 \(M = C_n\)，其中 \(C_n\) 是一个 \(n \times n\) 循环矩阵，其第一列由 \(T_n\) 的部分元素构造，例如取 Strang 预处理子：
\[\mathbf{c}_S = (t_0, t_1, \dots, t_{\lfloor n/2 \rfloor}, t_{-\lfloor n/2 \rfloor}, \dots, t_{-1})^T \]
（对称时即前一半对称元素）。循环矩阵 \(C_n\) 可被 Fourier 矩阵对角化：
\[C_n = F_n^* \Lambda F_n \]
其中 \(F_n\) 是 Fourier 矩阵，\(\Lambda = \text{diag}(\hat{\mathbf{c}}_S)\)，\(\hat{\mathbf{c}}_S = \text{FFT}(\mathbf{c}_S)\)。因此求解预处理子方程 \(M \mathbf{z} = \mathbf{r}\) 只需三次 FFT：
\[\mathbf{z} = F_n^* (\Lambda^\dagger \odot (F_n \mathbf{r})) \]
其中 \(\Lambda^\dagger\) 是 \(\Lambda\) 的伪逆（避免零特征值）。复杂度为 \(O(n \log n)\)。
步骤4：算法实现
结合上述步骤，预处理共轭梯度法（PCG）流程如下：
- 初始化：给定初始解 \(\mathbf{x}_0\)（常取零向量），计算残差 \(\mathbf{r}_0 = \mathbf{b} - T_n \mathbf{x}_0\)（用 FFT 加速矩阵向量乘）。
- 预处理：求解 \(M \mathbf{z}_0 = \mathbf{r}_0\) 得 \(\mathbf{z}_0\)（用 FFT 加速）。
- 设 \(\mathbf{p}_0 = \mathbf{z}_0\)。
- 迭代 对于 \(k=0,1,\dots\) 直到收敛：
a. 计算 \(\alpha_k = \frac{\mathbf{r}_k^T \mathbf{z}_k}{\mathbf{p}_k^T (T_n \mathbf{p}_k)}\)（分母用 FFT 加速）。
b. 更新 \(\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k\)。
c. 更新 \(\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k (T_n \mathbf{p}_k)\)。
d. 若 \(\|\mathbf{r}_{k+1}\|\) 小于容忍误差，则停止。
e. 求解 \(M \mathbf{z}_{k+1} = \mathbf{r}_{k+1}\) 得 \(\mathbf{z}_{k+1}\)（FFT 加速）。
f. 计算 \(\beta_k = \frac{\mathbf{r}_{k+1}^T \mathbf{z}_{k+1}}{\mathbf{r}_k^T \mathbf{z}_k}\)。
g. 更新 \(\mathbf{p}_{k+1} = \mathbf{z}_{k+1} + \beta_k \mathbf{p}_k\)。
步骤5：复杂度与优势
- 每轮迭代：两次 FFT 加速的矩阵向量乘（步骤4a、4c）和一次预处理子求解（步骤4e），总复杂度 \(O(n \log n)\)。
- 存储：仅需存储向量和 Toeplitz 序列，内存 \(O(n)\)。
- 预处理子 \(M\) 近似 \(T_n\) 且易求逆，可显著降低条件数，使 CG 迭代次数从 \(O(n)\) 减至 \(O(\log n)\) 或常数级。
- 整体复杂度从直接法的 \(O(n^3)\) 降为 \(O(n \log n \times \text{迭代次数})\)，适合大规模问题。
总结
本算法通过 FFT 加速 Toeplitz 矩阵的矩阵向量乘法和预处理子求解，将预处理共轭梯度法的每步成本降至 \(O(n \log n)\)，并利用循环预处理子改善收敛性，是求解对称正定 Toeplitz 线性方程组的高效迭代算法。核心在于利用 Toeplitz 矩阵的循环嵌入和 Fourier 变换的快速卷积性质。