分块矩阵的Kronecker积求解Lyapunov方程算法
1. 问题描述
Lyapunov方程 是控制论、系统理论和数值线性代数中的一个基本矩阵方程，形式通常为：
\[A X + X A^T = -C \]
其中 \(A, C \in \mathbb{R}^{n \times n}\) 是已知矩阵，且 \(C\) 通常对称（常为半正定），\(X \in \mathbb{R}^{n \times n}\) 是未知矩阵，需要求解。这个方程出现在线性时不变系统的稳定性分析、模型降阶、最优控制等场景中。
当矩阵规模 \(n\) 很大时，直接求解 \(n^2\) 维的线性系统计算量极大。本题目介绍一种利用 Kronecker 积和分块矩阵技术，将 Lyapunov 方程转化为结构化线性系统，并设计高效求解算法的过程。
2. 将 Lyapunov 方程转化为线性系统
2.1 Kronecker 积与向量化
- 定义 向量化 操作 \(\text{vec}(X)\)：将矩阵 \(X\) 按列堆叠成一个 \(n^2 \times 1\) 的列向量。
- Kronecker 积 性质：对于适当维数的矩阵 \(A, B, X\)，有
\[ \text{vec}(AXB) = (B^T \otimes A) \text{vec}(X) \]
- 对 Lyapunov 方程 \(A X + X A^T = -C\) 两边同时向量化：
\[ \text{vec}(AXI) + \text{vec}(I X A^T) = -\text{vec}(C) \]
利用上述性质：
\[ (I \otimes A) \text{vec}(X) + (A \otimes I) \text{vec}(X) = -\text{vec}(C) \]
合并得：
\[ \big[ (I \otimes A) + (A \otimes I) \big] \text{vec}(X) = -\text{vec}(C) \]
记：
\[ \mathcal{A} = (I \otimes A) + (A \otimes I) \in \mathbb{R}^{n^2 \times n^2}, \quad b = -\text{vec}(C) \in \mathbb{R}^{n^2} \]
则原方程化为大型线性系统：
\[ \mathcal{A} x = b, \quad x = \text{vec}(X) \]
3. 分块结构分析与求解难点
3.1 矩阵 \(\mathcal{A}\) 的分块特性
- 因为 \(I \otimes A = \text{blockdiag}(A, A, \dots, A)\)（由 \(n\) 个 \(A\) 块构成的分块对角矩阵）。
- 而 \(A \otimes I\) 是一个块矩阵，其第 \((i,j)\) 块是 \(a_{ij} I\)（其中 \(a_{ij}\) 是 \(A\) 的元素）。
- 因此 \(\mathcal{A}\) 是一个稀疏分块矩阵，但维度是 \(n^2 \times n^2\)，直接存储和求解（如高斯消元）的复杂度为 \(O(n^6)\)，不可接受。
3.2 利用 \(A\) 的特殊结构
如果 \(A\) 可对角化或可相似化为较简单形式（如 Schur 分解），则可大幅简化。
4. 基于 Schur 分解的 Bartels–Stewart 类算法思想
实际上，标准的快速 Lyapunov 求解算法（如 Bartels–Stewart 算法）是先对 \(A\) 进行实 Schur 分解，然后通过回代求解。但这里我们重点讲解如何利用分块 Kronecker 结构设计迭代法或直接法。
4.1 转化为块三角系统（当 \(A\) 为上三角时）
- 假设 \(A\) 已经是上三角矩阵（可通过 Schur 分解实现，但这里先假设，以展示分块求解思路）。
- 设 \(A = (a_{ij})\)，\(a_{ij} = 0\) 当 \(i > j\)。
- 将未知矩阵 \(X\) 也分块按列考虑：记 \(X = [x_1, x_2, \dots, x_n]\)，\(C = [c_1, c_2, \dots, c_n]\)。
- Lyapunov 方程第 \(j\) 列为：
\[ A x_j + \sum_{i=1}^n x_i a_{ij} = -c_j \]
这里用到了 \(X A^T\) 的第 \(j\) 列为 \(\sum_{i=1}^n x_i a_{ij}\)。
- 由于 \(A\) 上三角，\(a_{ij} = 0\) 对 \(i > j\)，因此上式可重写为：
\[ A x_j + \sum_{i=1}^j x_i a_{ij} = -c_j, \quad j=1,\dots,n \]
移项：
\[ (A + a_{jj} I) x_j = -c_j - \sum_{i=1}^{j-1} a_{ij} x_i \]
- 这是一个块上三角递归系统：当按 \(j=1,2,\dots,n\) 求解时，右边只依赖已求出的 \(x_1,\dots,x_{j-1}\)。
- 因此，只需顺序解 \(n\) 个 \(n\times n\) 线性系统：
\[ (A + a_{jj} I) x_j = \text{已知向量} \]
复杂度从 \(O(n^6)\) 降至 \(O(n^4)\)（因为每个系统解需 \(O(n^3)\)，共 \(n\) 个系统）。
5. 利用 Kronecker 积的分块迭代求解（当 \(A\) 为一般矩阵时）
如果 \(A\) 不是上三角，我们可以用分裂迭代法，利用 \(\mathcal{A} = I\otimes A + A\otimes I\) 的分裂。
5.1 分裂为
\[\mathcal{A} = M - N \]
其中一种自然选择：
\[M = I \otimes A, \quad N = -A \otimes I \]
但这样不收敛性好。更好的选择是对称分裂，例如当 \(A\) 稳定（特征值实部小于0）时，可用：
\[M = I \otimes A + A \otimes I_{\text{low}} \]
这里更实用的是用 Krylov 子空间方法（如 GMRES）直接求解 \(\mathcal{A}x = b\)，并利用 \(\mathcal{A}\) 的 Kronecker 积结构快速计算矩阵-向量乘积，而不显式构造 \(\mathcal{A}\)。
5.2 快速矩阵-向量乘积
对任意向量 \(y \in \mathbb{R}^{n^2}\)，要计算 \(z = \mathcal{A} y\)：
- 将 \(y\) 重组为 \(n\times n\) 矩阵 \(Y\)，使得 \(y = \text{vec}(Y)\)。
- 由 \(\mathcal{A} = I\otimes A + A\otimes I\) 的性质：
\[ z = \text{vec}(A Y + Y A^T) \]
- 因此计算步骤为：
- 从 \(y\) 得到矩阵 \(Y\)（reshape）。
- 计算 \(W_1 = A Y\)（矩阵乘，\(O(n^3)\) 但可优化如果 \(A\) 稀疏）。
- 计算 \(W_2 = Y A^T\)。
- 计算 \(Z = W_1 + W_2\)。
- 将 \(Z\) 向量化得 \(z\)。
- 这样，每次矩阵-向量乘只需求两个 \(n\times n\) 矩阵乘法，无需构造 \(n^2\times n^2\) 大矩阵。
6. 算法步骤总结
结合上述思想，一个实用的基于分块 Kronecker 结构和 Krylov 子空间的算法步骤如下：
- 输入：矩阵 \(A, C \in \mathbb{R}^{n \times n}\)，\(C\) 对称。
- 预处理：对 \(A\) 进行实 Schur 分解 \(A = Q T Q^T\)，其中 \(Q\) 正交，\(T\) 上三角。将方程转换为：
\[ T \tilde{X} + \tilde{X} T^T = -\tilde{C}, \quad \tilde{C} = Q^T C Q, \quad \tilde{X} = Q^T X Q \]
求解后再回代 \(X = Q \tilde{X} Q^T\)。
3. 解三角 Lyapunov 方程（此时 \(T\) 上三角）：
- 对 \(j=1\) 到 \(n\)：
\[ (T + t_{jj} I) \tilde{x}_j = -\tilde{c}_j - \sum_{i=1}^{j-1} t_{ij} \tilde{x}_i \]
解这个 $n$ 阶线性系统得 $\tilde{x}_j$。
- 得到 \(\tilde{X} = [\tilde{x}_1, \dots, \tilde{x}_n]\)。
- 回代：\(X = Q \tilde{X} Q^T\)。
- 输出：解矩阵 \(X\)。
7. 复杂度与注意事项
- 实 Schur 分解：\(O(n^3)\)。
- 解三角 Lyapunov 方程：每个 \(j\) 需解一个 \(n\times n\) 系统，用直接法（如 LU 分解）需 \(O(n^3)\)，但注意系数矩阵 \(T + t_{jj}I\) 是上三角，可用前代回代法在 \(O(n^2)\) 内解出。因此总复杂度为 \(O(n^3)\)。
- 若 \(A\) 是大型稀疏矩阵，可采用 Krylov 方法 解 \(\mathcal{A}x = b\)，配合前述快速矩阵-向量乘积，并利用预处理技术（例如基于 \(A\) 的近似 Schur 分解的预处理子）加速收敛。
8. 核心思想
本算法的核心是：
- 通过 Kronecker 积的向量化技巧 将 Lyapunov 方程转为大线性系统。
- 利用 矩阵乘法结构 快速实现矩阵-向量乘，避免大矩阵存储。
- 通过 Schur 分解 将一般问题转化为上三角问题，然后利用分块递归法高效求解。
- 整个算法将计算复杂度从 \(O(n^6)\) 降至 \(O(n^3)\)，适合中等规模（\(n\) 几千以内）的稠密 Lyapunov 方程。对于更大规模或稀疏 \(A\)，则采用迭代法（如 Krylov 方法）配合预处理。