A class of low-rank short recurrences for nonsymmetric linear matrix equations

技术摘要：一类用于非对称线性矩阵方程的低秩短递推方法

问题陈述
本文探讨了如下形式的大规模多术语线性矩阵方程的数值解：
$$A_1XB_1 + \dots + A_pXB_p = CD^T$$
其中 $A_i \in \mathbb{R}^{n_A \times n_A}$ 和 $B_i \in \mathbb{R}^{n_B \times n_B}$ 是大型、稀疏且非对称的矩阵，而 $C, D$ 是秩为 $q \ll n_A, n_B$ 的瘦高型满秩矩阵。解 $X$ 通常为矩形矩阵。虽然该方程可通过克罗内克积向量化为标准的线性系统 $Ax=b$，但由此产生的维度$n_A n_B$使得受限于内存约束的全向量迭代方法变得不可行。作者聚焦于$p > 2$的情况（将$p=2$ 时的 Sylvester 方程推广），并假设解 $X$ 具有低秩近似性质，这一性质常见于偏微分方程（PDE）的离散化、随机 Galerkin 方法以及控制系统中。

方法论
作者提出了一类新的面向矩阵的迭代方法族，将基于向量的短递推方法（特别是极小残差法和广义共轭残差法）推广至非对称矩阵情形。核心策略是在整个迭代过程中将解和残差保持为因式分解的低秩形式，以最小化内存使用。

1. 短递推公式：
   该方法不将问题转化为向量系统，而是利用低秩修正直接更新近似解 $X_k$：
   $$X_{k+1} = X_k + P^{(l)}_k \alpha_k (P^{(r)}_k)^T$$
   其中 $P^{(l)}_k$ 和 $P^{(r)}_k$ 是低秩因子矩阵，$\alpha_k$ 是通过求解局部极小化问题确定的小型稠密矩阵。由此导出了两种具体算法：

   • ss–mr（子空间极小残差法）： 向量 MR 迭代的矩阵对应形式，在每一步最小化残差的 Frobenius 范数。
   • ss–gcr(1)（子空间广义共轭残差法）： 向量 GCR 方法的矩阵对应形式，采用短递推（$\ell=1$），通过对搜索方向施加正交条件以加速收敛。

2. 局部优化：
   系数矩阵 $\alpha_k$（以及在 ss–gcr(1) 中的更新系数 $\beta_k$）通过求解一个缩减的多术语矩阵方程获得。这涉及将算子 $L$ 投影到由当前因子张成的子空间上。根据秩的大小，这些缩减系统要么通过直接 Cholesky 分解求解（针对小秩），要么通过预处理共轭梯度法（PCG）求解（针对较大秩）。

3. 秩截断与随机化：
   为防止秩增长并控制内存，作者采用了低秩截断策略。关键在于，他们引入了随机化草图技术（使用随机子采样三角变换），用于计算残差的 Frobenius 范数并执行秩截断，而无需显式构建大型稠密残差矩阵。这使得该方法能够处理 $n_A$ 或 $n_B$ 过大以至于无法存储完整残差的情况。

4. 预处理：
   该框架支持通过将算子分解为 $L = P - N$ 来进行预处理。作者讨论了一阶预处理（求逆单个 $A_i, B_i$ 对）和二阶预处理（通过 ADI 迭代近似求解 Sylvester 方程）。预处理后的残差用于生成下一个搜索方向。

主要贡献

• 推广至非对称算子： 本文将此前仅适用于对称/正定算子的短递推矩阵方法推广到了通用的非对称情形。
• 内存效率： 通过将短递推与低秩因式分解表示及随机化草图技术相结合，该方法避免了存储长度为 $n_A n_B$ 的完整向量，使其在向量化方法失效的问题中变得可行。
• 算法框架： 推导了具体的算法（ss–mr 和 ss–gcr(1)），在多术语方程的计算成本、内存和收敛性质之间取得了平衡。
• 随机 Galerkin 应用： 将这些方法应用于一个具有挑战性的参数化 PDE 问题（具有随机输入的达西流），该问题通过随机 Galerkin 混合有限元方法（SG-MFEM）离散化，形成了一个具有多项的非对称矩阵方程。

结果
数值实验在两个基准问题上进行：

1. 对流 - 扩散问题： 使用有限差分法离散化的二维稳态对流 - 扩散方程。
   • 新方法（ss–mr 和 ss–gcr(1)）在 CPU 时间上比最先进的低秩 GMRES（lr-pgmres）收敛快得多，通常快一个数量级。
   • 它们所需的内存远少于 lr-pgmres。
   • 与基于向量的 BiCGStab(2) 在克罗内克形式上的表现相比，矩阵方法在时间和内存方面均具有巨大优势，向量方法在更大的网格尺寸下要么无法收敛，要么耗尽内存。
2. 随机达西流： 一个具有不确定渗透率的参数化鞍点问题。
   • 该方法利用适度的内存资源成功求解了自由度超过 $10^9$ 的问题。
   • 迭代次数与空间网格尺寸和随机参数的数量无关。
   • 在“慢衰减”情况（需要更高秩）下，由于无需计算 $\beta_k$ 更新带来的更低成本，ss–mr 在 CPU 时间上优于 ss–gcr(1)，而两者均优于基于向量的 MINRES 求解器，后者在最大规模实例上因内存/时间限制而失败。

意义与主张
本文声称，所提出的方法论为求解大规模非对称多术语矩阵方程提供了一种稳健且高效的替代方案，特别是在解本身具有低秩特性的场景中。作者强调，他们的方法有效地弥合了矩阵形式（保持结构）的理论优势与低秩近似的实际需求之间的差距。

其意义在于能够解决由于克罗内克积引入的“维数灾难”而导致标准基于向量的迭代求解器目前无法处理的问题。随机化策略的整合被强调为实现无需显式矩阵构建即可计算残差范数和截断秩的关键赋能因素。作者指出，虽然这些方法行之有效，但其性能在很大程度上取决于良好预处理器的可用性以及解的低秩可近似性。这项工作被视为此类问题算法发展的前进一步，此类问题在历史上一直滞后于其在偏微分方程和数据科学中的应用。