A block IDR(s) method for nonsymmetric linear systems with multiple right-hand sides

基于诱导降维（IDR）定理的IDR（s）是一类新的求解大型非对称线性系统的有效算法。本文将IDR（s）推广到求解具有多个右边的大型非对称线性系统。

本文首先给出了IDR定理的一个变体，然后在变体IDR（s）定理的基础上提出了块IDR（s）。

1. Introduction

考虑求解具有多个右侧边的大的稀疏线性系统
$$AX=B\text{\hspace{1em}(1)}$$

系数矩阵A是n阶非奇异实数矩阵，, $X=[x_1,x_2,...,x_m]$ and $B=[b_1,b_2,...,b_m]\in R^{n\times m}$(m ≪ n)

本文的其余部分组织如下。在下一节中，我们将回顾块Krylov子空间和IDR(s)。然后，我们将给出IDR定理的一个变体，它是IDR定理的推广和我们的块IDR(s)算法的理论基础，并在第3节中分析块方法。在第4节中，一些数值结果证明了该方法的有效性。最后，我们在第5节中作了一些总结性。

符号。$M_{n\times n}$，$N_{n\times m}$表示矩阵。如果一个矩阵的维数从上下文中是明显的，并且没有混淆，我们将去掉索引，用$M$表示$M_{n\times n}$。$M^H$是M的厄米特转置。

$\Vert M{\Vert }_F^2=Tr(M^{H}M)$，其中$\Vert \cdot {\Vert }_F$和$Tr(\cdot )$分别表示方阵的F范数和迹。除非另有说明，‖·‖是贯穿整个论文的欧几里得范数。
$\mathcal{N}(M)$表示矩阵M的空空间。

2. The block Krylov subspace and IDR(s)

2.1. Block Krylov subspace

设$X_0\in R^{n\times m}$为初始猜测，$R_0=B-A{X}_0$为对应的块残差矩阵。

Definition 2.1 由$A$产生的子空间${\mathcal{K}}_k(A，R_0)$和A应用于$R_0$的升幂

被称为块Krylov子空间。

矩阵Krylov子空间，其定义如下。
Definition 2.2 由$A$产生的子空间${\mathbb{K}}_k(A，R_0)$和A应用于$R_0$的升幂

Definition 2.3 正整数$v:=v(R_0,A)$，定义为

对于A，被称为$R_0$的块级。

Corollary 2.4 设$X_{\ast }$是$AX=B$的精确块解

用于寻找近似解的方法 $X_k\in X_0+{\mathcal{K}}_k\left(A,R_0\right)$ 被称为块方法， 选择 $X_k\in X_0+{\mathbb{K}}_k\left(A,R_0\right)$ 导致所谓的全局方法. 如果我们选择 ${\gamma }_i={\alpha }_i{I}_{m\times m}(i=0,1,\dots ,k-1),{\mathbb{K}}_k$ and ${\mathcal{K}}_k$ 可以是同一个子空间。 由此看来，矩阵 Krylov 子空间${\mathbb{K}}_k\left(A,R_0\right)$ 可以认为是块 Krylov 子空间的子空间${\mathcal{K}}_k\left(A,R_0\right)$, i.e. ${\mathbb{K}}_k\left(A,R_0\right)\subset {\mathcal{K}}_k\left(A,R_0\right)$.

对于块求解器，让 $Z=\left[z_1,z_2,\dots ,z_m\right]\in {\mathcal{K}}_k$, where $z_i\in {\mathbb{R}}^n(i=1,\dots ,m)$. F从(2)的定义来看，有${\gamma }_j{}^{\prime }s\in {\mathbb{R}}^{m\times m}(j=0,\dots ,k-1)$ 满足
$$Z=\sum _{j=0}^{k-1}{A}^j{R}_0{\gamma }_j$$
这意味着
$$z_i=\sum _{l=1}^m\sum _{j=0}^{k-1}{\gamma }_j(l,i)A^j{R}_0(:,l)\in {\mathcal{B}}_k\left(A,R_0\right)$$
where ${\mathcal{B}}_k\left(A,R_0\right):={\mathcal{K}}_k\left(A,R_0(:,1)\right)+\cdots +{\mathcal{K}}_k\left(A,R_0(:,m)\right)$.
所以 $X_k=X_0+Z\in {\mathcal{K}}_k$ 的列对应于 $m$ 单右手线性系统的近似解。但是，与标准 Krylov 求解器不同，块 Krylov 求解器在每个右侧的搜索空间要大得多，即,近似解 $X_k(:,l)=X_0(:,l)+{\mathcal{B}}_k\left(A,R_0\right)$ are updated instead of $X_k(:,l)=X_0(:,l)+{\mathcal{K}}_k\left(A,R_0(:,l)\right)$. 这是使用块方法的主要原因。

2.2. IDR(s)

Theorem 2.5 (IDR). 设 $A\in {\mathbb{C}}^{n\times n},v_0\in {\mathbb{C}}^n$, and $g_0={\mathcal{K}}_v\left(A,v_0\right)$. 设 $S\subset {\mathbb{C}}^n$ 并定义递归子空间 $g_j$ 为
$$g_j=\left(I-{\omega }_{j}A\right)\left(g_{j-1}\cap S\right),\left({\omega }_j\ne 0\right)\in \mathbb{C},j=1,2,\dots$$
如果$g_0\cap S$ 不包含 $A$,的任何特征向量，则保持如下：
(a) $g_j\subset g_{j-1},\forall j>0$.
(b) g j = { 0 } \mathcal{g}_{j}=\{0\} gj​={0} for some $j\le n$.

从该定理中，我们知道级数嵌套子空间的维数随着序列子空间$g_j$的缩小而减小。如果所有的剩余$r_{i}s$都可以在嵌套的子空间$g_j$中构造，我们可以在有限步长中得到近似解。在IDR(s)方法[23]的一般情况下，最多需要$n+n/s$个矩阵向量积。

s+2项IDR(s)算法可以作为IDR定理的推导如下。

要初始化$r_1,\dots ,r_s$，可以使用任何相当迭代的方法，例如，BiCGStab

然后，从方程(3)(4)，近似解$x_{i+1}$可以更新为

为了确定$s$个变量${\gamma }_j$，可以选择空间$S$为某些$n\times s$矩阵的左空空间，$P=[p_1,...,p_s]$, i.e., $S=N(P^H)$,这可以随机生成，因为空间$S\cap g_0$包含A的一些特征向量(s)的概率为零[23]。然后可以从方程中求解${\gamma }_j$
$$P^H{v}_i=0.$$

在形成整个算法之前，还需要计算一个参数${\omega }_{j+1}$，[23]建议通过每$s+1$步最小化残差$r_{i+1}$的范数来选择$\omega$。

3. The block IDR(s)

在本节中，我们考虑具有多个右侧的非对称线性系统。为了提出IDR(s)的块版本，我们首先给出了IDR定理的一个变体，它是在第3.1节中对IDR定理的一个扩展。为了从变体IDR定理中给出一个可执行的过程，我们分别在第3.2节和第3.3节中讨论了实现细节，并阐述了所提出的算法。然后，在第3.4节中，我们比较了IDR(s)和块IDR(s)之间的计算成本和内存需求，并分析了通过矩阵向量运算计算的收敛性。最后，在第3.5节中给出了一个块IDR(s)算法的预处理版本。

3.1. A variant of the IDR theorem

从块初始残差r0和系数矩阵A中，可以生成块Krylov子空间。基于块Krylov子空间，可以自然地得到IDR定理的一个变体。在本节中，我们首先给出了变体IDR定理，这是我们所提出的方法的理论基础。

3.2. Implementation details

因为块IDR(s)是IDR(s)的自然扩展，所以s+2项块IDR(s)方法可以推导出类似于等式。（3）-（5）作为变体IDR定理的翻译。

3.5. Preconditioning