Krylov子空间方法在求解Sylvester方程中的应用
字数 3271 2025-12-07 16:19:19

Krylov子空间方法在求解Sylvester方程中的应用

1. 问题描述
Sylvester方程是矩阵方程的一种常见形式，表示为：

\[A X - X B = C \]

其中 \(A \in \mathbb{R}^{m \times m}\)，\(B \in \mathbb{R}^{n \times n}\)，\(C \in \mathbb{R}^{m \times n}\) 是已知矩阵，而 \(X \in \mathbb{R}^{m \times n}\) 是未知矩阵。当 \(m\) 和 \(n\) 很大时，直接求解（如转化为线性方程组）计算量过大。Krylov子空间方法通过投影到低维空间迭代逼近解，可高效处理大规模稀疏矩阵的情形。

2. 问题转化
将矩阵 \(X\) 和 \(C\) 按列向量化（vec操作），Sylvester方程可等价为线性系统：

\[(I_n \otimes A - B^T \otimes I_m) \text{vec}(X) = \text{vec}(C) \]

其中 \(\otimes\) 表示Kronecker积。但该系统的系数矩阵大小为 \(mn \times mn\)，直接求解不现实。Krylov子空间方法避免显式构造此矩阵，而是利用矩阵与向量的乘积在原始矩阵形式下进行迭代。

3. Krylov子空间构造
对向量形式的线性系统 \(\mathcal{A} x = c\)，其中 \(\mathcal{A} = I_n \otimes A - B^T \otimes I_m\)，\(x = \text{vec}(X)\)，\(c = \text{vec}(C)\)，标准Krylov子空间为：

\[\mathcal{K}_k(\mathcal{A}, c) = \text{span}\{c, \mathcal{A}c, \mathcal{A}^2 c, \dots, \mathcal{A}^{k-1} c\} \]

但 \(\mathcal{A}\) 与向量 \(v\) 的乘积等价于矩阵运算：

\[\mathcal{A} v \quad \Leftrightarrow \quad A V - V B^T \]

其中 \(V\) 是 \(v\) 重塑成的 \(m \times n\) 矩阵。因此，我们可在矩阵层面构造Krylov子空间，避免处理 \(mn\) 维向量。

4. 矩阵Krylov子空间
从初始残差 \(R_0 = C - (A X_0 - X_0 B)\) 开始（初始猜测 \(X_0\) 常取零矩阵），定义块Krylov子空间：

\[\mathcal{K}_k(A, B, R_0) = \text{span} \{ R_0, A R_0 - R_0 B, A^2 R_0 - 2 A R_0 B + R_0 B^2, \dots \} \]

更系统的构造方式是采用Arnoldi或Lanczos过程（当 \(A\) 和 \(B\) 对称时）生成一组正交基。这里以一般矩阵为例，采用全局Arnoldi过程（Global GMRES变体）：

  • 初始残差矩阵 \(R_0 = C - (A X_0 - X_0 B)\)，令 \(V_1 = R_0 / \|R_0\|_F\)。
  • 对 \(j = 1, 2, \dots, k\) 执行：
    1. 计算 \(W = A V_j - V_j B\)。
    2. 对 \(i = 1, \dots, j\)，计算 \(h_{ij} = \langle V_i, W \rangle_F\)（Frobenius内积），并令 \(W := W - h_{ij} V_i\)。
    3. 令 \(h_{j+1,j} = \|W\|_F\)，若 \(h_{j+1,j} \neq 0\)，则 \(V_{j+1} = W / h_{j+1,j}\)。
      此时得到正交基矩阵 \(V_1, \dots, V_{k+1}\) 和上Hessenberg矩阵 \(\bar{H}_k \in \mathbb{R}^{(k+1) \times k}\)。

5. 投影与近似解
设近似解 \(X_k = X_0 + \sum_{i=1}^k y_i V_i\)，写成矩阵形式 \(X_k = X_0 + [V_1, \dots, V_k] y\)，其中 \(y \in \mathbb{R}^k\)。将残差投影到子空间：

\[R_k = C - (A X_k - X_k B) = R_0 - [A V_1 - V_1 B, \dots, A V_k - V_k B] y \]

由Arnoldi过程可知，\([A V_1 - V_1 B, \dots, A V_k - V_k B] = [V_1, \dots, V_{k+1}] \bar{H}_k\)。且 \(R_0 = \|R_0\|_F V_1 = \beta V_1\)（记 \(\beta = \|R_0\|_F\)）。于是：

\[R_k = [V_1, \dots, V_{k+1}] (\beta e_1 - \bar{H}_k y) \]

其中 \(e_1\) 是第一个单位向量。选取 \(y\) 最小化残差范数 \(\|R_k\|_F\)，等价于解最小二乘问题：

\[y_k = \arg\min_y \|\beta e_1 - \bar{H}_k y\|_2 \]

用Givens旋转递推求解，得到近似解 \(X_k = X_0 + \sum_{i=1}^k y_i V_i\)。

6. 算法步骤总结

  1. 选择初始猜测 \(X_0\)，计算 \(R_0 = C - (A X_0 - X_0 B)\)，\(\beta = \|R_0\|_F\)。
  2. 若 \(\beta < \epsilon\) 则停止，否则令 \(V_1 = R_0 / \beta\)。
  3. 对 \(j = 1, 2, \dots, k\) 执行Arnoldi过程：
    • 计算 \(W = A V_j - V_j B\)。
    • 对 \(i = 1, \dots, j\) 正交化：\(h_{ij} = \text{tr}(V_i^T W)\)，\(W := W - h_{ij} V_i\)。
    • \(h_{j+1,j} = \|W\|_F\)，若 \(h_{j+1,j} = 0\) 则终止，否则 \(V_{j+1} = W / h_{j+1,j}\)。
  4. 解最小二乘问题 \(\min_y \|\beta e_1 - \bar{H}_k y\|_2\) 得 \(y_k\)。
  5. 更新近似解 \(X_k = X_0 + [V_1, \dots, V_k] y_k\)。
  6. 若残差满足精度要求则停止，否则可重启或继续扩展子空间。

7. 关键点与适用场景

  • 该方法只需矩阵与矩阵的乘积 \(A V\) 和 \(V B\)，适合 \(A, B\) 稀疏或具有快速矩阵乘法的情况。
  • 当 \(A\) 和 \(B\) 对称时，可利用短递归（如Lanczos）减少存储。
  • 对于大规模问题，常采用重启策略防止子空间过大。
  • 该方法可视为将Krylov子空间方法（如GMRES）应用于Sylvester方程的自然推广，避免了向量化带来的维度爆炸。
Krylov子空间方法在求解Sylvester方程中的应用 1. 问题描述 Sylvester方程是矩阵方程的一种常见形式，表示为： \[ A X - X B = C \] 其中 \( A \in \mathbb{R}^{m \times m} \)，\( B \in \mathbb{R}^{n \times n} \)，\( C \in \mathbb{R}^{m \times n} \) 是已知矩阵，而 \( X \in \mathbb{R}^{m \times n} \) 是未知矩阵。当 \( m \) 和 \( n \) 很大时，直接求解（如转化为线性方程组）计算量过大。Krylov子空间方法通过投影到低维空间迭代逼近解，可高效处理大规模稀疏矩阵的情形。 2. 问题转化 将矩阵 \( X \) 和 \( C \) 按列向量化（vec操作），Sylvester方程可等价为线性系统： \[ (I_ n \otimes A - B^T \otimes I_ m) \text{vec}(X) = \text{vec}(C) \] 其中 \( \otimes \) 表示Kronecker积。但该系统的系数矩阵大小为 \( mn \times mn \)，直接求解不现实。Krylov子空间方法避免显式构造此矩阵，而是利用矩阵与向量的乘积在原始矩阵形式下进行迭代。 3. Krylov子空间构造 对向量形式的线性系统 \( \mathcal{A} x = c \)，其中 \( \mathcal{A} = I_ n \otimes A - B^T \otimes I_ m \)，\( x = \text{vec}(X) \)，\( c = \text{vec}(C) \)，标准Krylov子空间为： \[ \mathcal{K}_ k(\mathcal{A}, c) = \text{span}\{c, \mathcal{A}c, \mathcal{A}^2 c, \dots, \mathcal{A}^{k-1} c\} \] 但 \( \mathcal{A} \) 与向量 \( v \) 的乘积等价于矩阵运算： \[ \mathcal{A} v \quad \Leftrightarrow \quad A V - V B^T \] 其中 \( V \) 是 \( v \) 重塑成的 \( m \times n \) 矩阵。因此，我们可在矩阵层面构造Krylov子空间，避免处理 \( mn \) 维向量。 4. 矩阵Krylov子空间 从初始残差 \( R_ 0 = C - (A X_ 0 - X_ 0 B) \) 开始（初始猜测 \( X_ 0 \) 常取零矩阵），定义块Krylov子空间： \[ \mathcal{K}_ k(A, B, R_ 0) = \text{span} \{ R_ 0, A R_ 0 - R_ 0 B, A^2 R_ 0 - 2 A R_ 0 B + R_ 0 B^2, \dots \} \] 更系统的构造方式是采用Arnoldi或Lanczos过程（当 \( A \) 和 \( B \) 对称时）生成一组正交基。这里以一般矩阵为例，采用 全局Arnoldi过程 （Global GMRES变体）： 初始残差矩阵 \( R_ 0 = C - (A X_ 0 - X_ 0 B) \)，令 \( V_ 1 = R_ 0 / \|R_ 0\|_ F \)。 对 \( j = 1, 2, \dots, k \) 执行： 计算 \( W = A V_ j - V_ j B \)。 对 \( i = 1, \dots, j \)，计算 \( h_ {ij} = \langle V_ i, W \rangle_ F \)（Frobenius内积），并令 \( W := W - h_ {ij} V_ i \)。 令 \( h_ {j+1,j} = \|W\| F \)，若 \( h {j+1,j} \neq 0 \)，则 \( V_ {j+1} = W / h_ {j+1,j} \)。 此时得到正交基矩阵 \( V_ 1, \dots, V_ {k+1} \) 和上Hessenberg矩阵 \( \bar{H}_ k \in \mathbb{R}^{(k+1) \times k} \)。 5. 投影与近似解 设近似解 \( X_ k = X_ 0 + \sum_ {i=1}^k y_ i V_ i \)，写成矩阵形式 \( X_ k = X_ 0 + [ V_ 1, \dots, V_ k ] y \)，其中 \( y \in \mathbb{R}^k \)。将残差投影到子空间： \[ R_ k = C - (A X_ k - X_ k B) = R_ 0 - [ A V_ 1 - V_ 1 B, \dots, A V_ k - V_ k B ] y \] 由Arnoldi过程可知，\( [ A V_ 1 - V_ 1 B, \dots, A V_ k - V_ k B] = [ V_ 1, \dots, V_ {k+1}] \bar{H}_ k \)。且 \( R_ 0 = \|R_ 0\|_ F V_ 1 = \beta V_ 1 \)（记 \( \beta = \|R_ 0\| F \)）。于是： \[ R_ k = [ V_ 1, \dots, V {k+1}] (\beta e_ 1 - \bar{H}_ k y) \] 其中 \( e_ 1 \) 是第一个单位向量。选取 \( y \) 最小化残差范数 \( \|R_ k\|_ F \)，等价于解最小二乘问题： \[ y_ k = \arg\min_ y \|\beta e_ 1 - \bar{H}_ k y\| 2 \] 用Givens旋转递推求解，得到近似解 \( X_ k = X_ 0 + \sum {i=1}^k y_ i V_ i \)。 6. 算法步骤总结 选择初始猜测 \( X_ 0 \)，计算 \( R_ 0 = C - (A X_ 0 - X_ 0 B) \)，\( \beta = \|R_ 0\|_ F \)。 若 \( \beta < \epsilon \) 则停止，否则令 \( V_ 1 = R_ 0 / \beta \)。 对 \( j = 1, 2, \dots, k \) 执行Arnoldi过程： 计算 \( W = A V_ j - V_ j B \)。 对 \( i = 1, \dots, j \) 正交化：\( h_ {ij} = \text{tr}(V_ i^T W) \)，\( W := W - h_ {ij} V_ i \)。 \( h_ {j+1,j} = \|W\| F \)，若 \( h {j+1,j} = 0 \) 则终止，否则 \( V_ {j+1} = W / h_ {j+1,j} \)。 解最小二乘问题 \( \min_ y \|\beta e_ 1 - \bar{H}_ k y\|_ 2 \) 得 \( y_ k \)。 更新近似解 \( X_ k = X_ 0 + [ V_ 1, \dots, V_ k] y_ k \)。 若残差满足精度要求则停止，否则可重启或继续扩展子空间。 7. 关键点与适用场景 该方法只需矩阵与矩阵的乘积 \( A V \) 和 \( V B \)，适合 \( A, B \) 稀疏或具有快速矩阵乘法的情况。 当 \( A \) 和 \( B \) 对称时，可利用短递归（如Lanczos）减少存储。 对于大规模问题，常采用 重启策略 防止子空间过大。 该方法可视为将Krylov子空间方法（如GMRES）应用于Sylvester方程的自然推广，避免了向量化带来的维度爆炸。