求解双鞍点线性系统的一种改进维数分裂预处理子

期刊菜单

求解双鞍点线性系统的一种改进维数分裂预处理子
An Improved Dimensional Splitting Preconditioner for Solving Double Saddle Point Linear Systems

DOI: 10.12677/aam.2024.137338, PDF, HTML, XML, 下载: 7 浏览: 14
作者: 谷玲倩：温州大学数理学院，浙江温州
关键词: 双鞍点问题；矩阵分裂；谱性质；最优参数；Double Saddle Point Problem； Matrix Splitting； Spectral Property； Optimal Parameters

摘要: 为了提高维数分裂(DS)预处理子和松弛维数分解(RDF)预处理子的性能，针对双鞍点问题，本文提出了一种改进维数分裂(IDS)预处理子，详细分析了预处理矩阵的谱性质并讨论了最优参数。数值实验结果验证了IDS预处理子的有效性。

Abstract: In order to improve the performance of dimensional splitting (DS) preconditioner and relaxed dimensional factorization (RDF) preconditioner, an improved dimensional splitting (IDS) preconditioner is proposed for the double saddle point problem. The spectral properties of the preconditioned matrix are analyzed in detail and the optimal parameters are discussed. The effectiveness of IDS preconditioner is verified by numerical experiments.

文章引用：谷玲倩. 求解双鞍点线性系统的一种改进维数分裂预处理子[J]. 应用数学进展, 2024, 13(7): 3541-3553. https://doi.org/10.12677/aam.2024.137338

1. 引言

考虑下面大型稀疏的3 × 3块线性方程组

$A x = [\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & A_{2} & B_{2}^{T} \\ - B_{1} & - B_{2} & 0 \end{matrix}] [\begin{matrix} u \\ v \\ p \end{matrix}] = [\begin{matrix} f_{1} \\ f_{2} \\ - g \end{matrix}] \equiv b$ (1)

其中， $A_{1} \in ℝ^{n_{1} \times n_{1}}$ 和 $A_{2} \in ℝ^{n_{2} \times n_{2}}$ 是正定的， $B_{1} \in ℝ^{m \times n_{1}}$ 和 $B_{2} \in ℝ^{m \times n_{2}}$ 是两个维数满足 $m \leq n_{1} + n_{2}$ 的矩形矩阵。形式(1)的线性系统常被称为双鞍点线性系统，出现在许多实际问题中，如计算流体力学[1]、约束优化[1]-[3]等。

令 $A = (\begin{matrix} A_{1} & 0 \\ 0 & A_{2} \end{matrix})$ ， $B = (\begin{matrix} B_{1} & B_{2} \end{matrix})$ ， $w = (\begin{matrix} u \\ v \end{matrix})$ ， $f = (\begin{matrix} f_{1} \\ f_{2} \end{matrix})$ ，则上述线性方程组(1)的3 × 3块系统可以等价地改写为

$[\begin{matrix} A & B^{T} \\ - B & 0 \end{matrix}] [\begin{matrix} w \\ p \end{matrix}] = [\begin{matrix} f \\ - g \end{matrix}]$ (2)

近几十年来，许多人对于线性方程组(1)或其等价的2 × 2块形式(2)，提出了许多有效的预处理子，如块对角预处理子[4]，块三角形预处理子[5]-[8]以及约束预处理子[9] [10]；更多细节参见[11]。

2. 一种改进维数分裂预处理子

我们首先引入求解双鞍点问题(1)的维数分裂(Dimensional Splitting，简记DS)预处理子。

记 $A = A_{1} + A_{2}$ ，其中

$A_{1} = [\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & 0 & 0 \\ - B_{1} & 0 & 0 \end{matrix}]$ 和 $A_{2} = [\begin{matrix} 0 & 0 & 0 \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & 0 \end{matrix}]$

基于上述分裂，Benzi和Guo [12]在2011年提出了如下维数分裂(DS)预处理子：

$P_{D S} = \frac{1}{α} (α I + A_{1}) (α I + A_{2}) = \frac{1}{α} [\begin{matrix} α I + A_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}] [\begin{matrix} α I & 0 & 0 \\ 0 & α I + A_{2} & B_{2}^{T} \\ 0 & - B_{2} & α I \end{matrix}]$ (3)

其中， $α > 0$ ， $I$ 是适当阶数的单位矩阵。

另一方面，基于维数分裂(DS)预处理子，Benzi等人[13]进一步提出了下述松弛维数分解(Relaxed Dimensional Factorization，简记RDF)预处理子。

$P_{R D F} = \frac{1}{α} [\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}] [\begin{matrix} α I & 0 & 0 \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & α I \end{matrix}]$ (4)

与维数分裂(DS)预处理子相比，松弛维数分解(RDF)预处理子更接近于系数矩阵 $A$ ，即差值 $R_{R D F} = P_{R D F} - A$ 比差值 $R_{D S} = P_{D S} - A$ 更接近于零矩阵。事实上，上面的分析由如下简单的计算给出：

$R_{D S} = [\begin{matrix} α I & - \frac{1}{α} B_{1}^{T} B_{2} & 0 \\ 0 & α I & 0 \\ 0 & 0 & α I \end{matrix}]$ 和 $R_{R D F} = [\begin{matrix} 0 & - \frac{1}{α} B_{1}^{T} B_{2} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & α I \end{matrix}]$

矩阵 $R_{D S}$ 比 $R_{R D F}$ 多两个非零块，这可能是松弛维数分解(RDF)预处理子优于维数分裂(DS)预处理子的原因。

针对双鞍点问题，2022年，Ren和Chen [14]通过引入一个新的参数 $β$ 以及省略一些项来避免因子 $α$

和 $\frac{1}{α}$ 出现在同一系数矩阵中，基于这一思想进而提出了IAPSS预处理子。本文受此思想启发，通过修改

这个矩阵分解中的某些项，我们可以得到改进维数分裂(IDS)预处理子，事实上，我们知道：

$[\begin{matrix} A_{1} + \frac{1}{α} B_{1}^{T} B_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ 0 & 0 & α I \end{matrix}] [\begin{matrix} I & 0 & 0 \\ 0 & I & 0 \\ - \frac{1}{α} B_{1} & 0 & I \end{matrix}] = [\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}]$ (5)

和

$[\begin{matrix} α I & 0 & 0 \\ 0 & A_{2} + \frac{1}{β} B_{2}^{T} B_{2} & B_{2}^{T} \\ 0 & 0 & β I \end{matrix}] [\begin{matrix} I & 0 & 0 \\ 0 & I & 0 \\ 0 & - \frac{1}{β} B_{2} & I \end{matrix}] = [\begin{matrix} α I & 0 & 0 \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & β I \end{matrix}]$ (6)

令

$P_{I D S} = \frac{1}{α} [\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}] [\begin{matrix} α I & 0 & 0 \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & β I \end{matrix}] = [\begin{matrix} A_{1} & - \frac{1}{α} B_{1}^{T} B_{2} & \frac{β}{α} B_{1}^{T} \\ 0 & A_{2} & B_{2}^{T} \\ - B_{1} & - B_{2} & β I \end{matrix}]$ (7)

其中， $α$ 和 $β$ 是给定的正常数。

$R_{I D S} = P_{I D S} - A = [\begin{matrix} 0 & - \frac{1}{α} B_{1}^{T} B_{2} & (\frac{β}{α} - 1) B_{1}^{T} \\ 0 & 0 & 0 \\ 0 & 0 & β I \end{matrix}]$ (8)

备注1。因为增加了一个参数，改进维数分裂(IDS)预处理子要比松弛维数分解(RDF)预处理子灵活，选取 $β = α$ 时，改进维数分裂(IDS)预处理子就退化为松弛维数分解(RDF)预处理子。

当我们使用改进维数分裂(IDS)预处理子(7)来加速Kryov子空间方法的收敛时，在预处理Kryov子空间方法的每一步都需要求解一个残差线性方程组 $P_{I D S} z = r$ ，其中z和r分别是当前和广义残差向量。令 $z = {[z_{1}^{T}, z_{2}^{T}, z_{3}^{T}]}^{T}$ ， $r = {[r_{1}^{T}, r_{2}^{T}, r_{3}^{T}]}^{T}$ ，引入一个介质向量 $t = {[t_{1}^{T}, t_{2}^{T}, t_{3}^{T}]}^{T}$ ，其中 $z_{1}, r_{1}, t_{1} \in ℝ^{n_{1}}$ ， $z_{2}, r_{2}, t_{2} \in ℝ^{n_{2}}$ 和 $z_{3}, r_{3}, t_{3} \in ℝ^{m}$ ，则线性方程组 $P_{I D S} z = r$ 可以等价地写为以下两个线性方程组：

$[\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}] [\begin{matrix} t_{1} \\ t_{2} \\ t_{3} \end{matrix}] = [\begin{matrix} α r_{1} \\ α r_{2} \\ α r_{3} \end{matrix}]$ (9)

和

$[\begin{matrix} α I & 0 & 0 \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & β I \end{matrix}] [\begin{matrix} z_{1} \\ z_{2} \\ z_{3} \end{matrix}] = [\begin{matrix} t_{1} \\ t_{2} \\ t_{3} \end{matrix}]$ (10)

从式(9)，我们有 $t_{2} = r_{2}$ ，且 $t_{1}$ 、 $t_{3}$ 由下述解得

$[\begin{matrix} A_{1} & B_{1}^{T} \\ - B_{1} & α I \end{matrix}] [\begin{matrix} t_{1} \\ t_{3} \end{matrix}] = [\begin{matrix} α r_{1} \\ α r_{3} \end{matrix}]$ (11)

或等价地

$[\begin{matrix} A_{1} + \frac{1}{α} B_{1}^{T} B_{1} & 0 \\ - B_{1} & α I \end{matrix}] [\begin{matrix} t_{1} \\ t_{3} \end{matrix}] = [\begin{matrix} α r_{1} - B_{1}^{T} r_{3} \\ α r_{3} \end{matrix}]$

使用相似的技术，线性系统(10)可以改写为 $z_{1} = t_{1} / α$ ，且

$[\begin{matrix} A_{2} + \frac{1}{β} B_{2}^{T} B_{2} & 0 \\ - B_{2} & β I \end{matrix}] [\begin{matrix} z_{2} \\ z_{3} \end{matrix}] = [\begin{matrix} t_{2} - \frac{1}{β} B_{2}^{T} t_{3} \\ t_{3} \end{matrix}]$

因此，残差线性方程组 $P_{I D S} z = r$ 可按如下算法1求解。

算法1. 求解 $P_{I D S} z = r$

1) 从以下线性子系统求解 $t_{1} \in ℝ^{n_{1}}$

$(A_{1} + \frac{1}{α} B_{1}^{T} B_{1}) t_{1} = α r_{1} - B_{1}^{T} r_{3}$ ；

2) 通过 $t_{3} = (α r_{3} + B_{1} t_{1}) / α$ 计算 $t_{3} \in ℝ^{m}$ ；

3) 通过 $t_{2} = r_{2}$ 计算 $t_{2} \in ℝ^{n_{2}}$ ；

4) 通过 $z_{1} = t_{1} / α$ 计算 $t_{1} \in ℝ^{n_{1}}$ ；

5) 从以下线性子系统求解 $z_{2} \in ℝ^{n_{2}}$

$(A_{2} + \frac{1}{β} B_{2}^{T} B_{2}) z_{2} = t_{2} - \frac{1}{β} B_{2}^{T} t_{3}$ ；

6) 通过 $z_{3} = (t_{3} + B_{2} z_{2}) / β$ 计算 $z_{3} \in ℝ^{m}$ ；

3. 预处理矩阵 $P_{I D S}^{- 1} A$ 的谱性质

我们首先引入一个引理，它将用于下面分析预处理矩阵 $P_{I D S}^{- 1} A$ 的谱性质。

引理1. 令 $T_{22} = [\begin{matrix} \frac{1}{α β} M_{2}^{- 1} B_{2}^{T} {\hat{S}}_{1} B_{2} & - M_{2}^{- 1} B_{2}^{T} (I - \frac{1}{β} {\hat{S}}_{1}) \\ - \frac{1}{α β} (I - \frac{1}{β} {\hat{S}}_{2}) {\hat{S}}_{1} B_{2} & (I - \frac{1}{β} {\hat{S}}_{2}) (I - \frac{1}{β} {\hat{S}}_{1}) \end{matrix}]$ ，则 $T_{22}$ 的0特征值的重数至少为 $n_{2}$ ，其余的特征值为 $1 - μ_{i}$ ，其中 $μ_{i}$ 是 $m \times m$ 矩阵 $Z = \frac{1}{β} ({\hat{S}}_{1} + {\hat{S}}_{2}) - \frac{α + β}{α β^{2}} {\hat{S}}_{1} {\hat{S}}_{2}$ 的特征值。

证明：首先，我们注意到

$T_{22} = \frac{1}{α β^{2}} [\begin{matrix} - β M_{2}^{- 1} B_{2}^{T} \\ β I - {\hat{S}}_{2} \end{matrix}] [\begin{matrix} - {\hat{S}}_{1} B_{2} & α (β I - {\hat{S}}_{1}) \end{matrix}] = X Y^{T}$ (12)

其中， $X = \frac{1}{α β^{2}} [\begin{matrix} - β M_{2}^{- 1} B_{2}^{T} \\ β I - {\hat{S}}_{2} \end{matrix}] \in ℝ^{(n_{2} + m) \times m}$ 和 $Y^{T} = [\begin{matrix} - {\hat{S}}_{1} B_{2} & α (β I - {\hat{S}}_{1}) \end{matrix}] \in ℝ^{m \times (n_{2} + m)}$ 。

因此，从式(12)我们可以知道， $T_{22}$ 是一个秩不超过m的矩阵，因此， $T_{22}$ 的特征值0的重数至少为 $n_{2}$ ，利用一个众所周知的结果([15]，定理1.3.20)，其剩余特征值就是 $m \times m$ 矩阵的特征值：

$\begin{matrix} Y^{T} X = \frac{1}{α β^{2}} [β {\hat{S}}_{1} B_{2} M_{2}^{- 1} B_{2}^{T} + α (β I - {\hat{S}}_{1}) (β I - {\hat{S}}_{2})] \\ = \frac{1}{α β^{2}} [β {\hat{S}}_{1} {\hat{S}}_{2} + α (β I - {\hat{S}}_{1}) (β I - {\hat{S}}_{2})] \\ = I - \frac{1}{β} ({\hat{S}}_{1} + {\hat{S}}_{2}) + \frac{α + β}{α β^{2}} {\hat{S}}_{1} {\hat{S}}_{2} \\ = I - Z \end{matrix}$

证毕。

基于引理1，我们有以下结果。

定理1. 假设 $A_{1} \in ℝ^{n_{1} \times n_{1}}$ 和 $A_{2} \in ℝ^{n_{2} \times n_{2}}$ 是正定矩阵， $B_{1} \in ℝ^{m \times n_{1}}$ 和 $B_{2} \in ℝ^{m \times n_{2}}$ 是长方阵且维数满足 $m \leq n_{1} + n_{2}$ ，则预处理矩阵 $P_{I D S}^{- 1} A$ 具有特征值1，其重数至少为 $n_{1} + n_{2}$ ，剩余特征值 $μ_{i}, i = 1, 2, \dots, m$ 是广义

特征值问题 $(S_{1} + S_{2}) ξ_{i} = μ_{i} \cdot \frac{1}{α} (α I + S_{1}) (β I + S_{2}) ξ_{i}$ (18)的特征值。

证明：令

$P_{1} = [\begin{matrix} A_{1} & 0 & B_{1}^{T} \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}]$ ， $P_{2} = [\begin{matrix} α I & 0 & 0 \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & β I \end{matrix}]$

易知

$P_{1}^{- 1} = [\begin{matrix} M_{1}^{- 1} & 0 & - \frac{1}{α} M_{1}^{- 1} B_{1}^{T} \\ 0 & \frac{1}{α} I & 0 \\ \frac{1}{α} B_{1} M_{1}^{- 1} & 0 & \frac{1}{α} I - \frac{1}{α^{2}} {\hat{S}}_{1} \end{matrix}]$

和

$P_{2}^{- 1} = [\begin{matrix} \frac{1}{α} I & 0 & 0 \\ 0 & M_{2}^{- 1} & - \frac{1}{β} M_{2}^{- 1} B_{2}^{T} \\ 0 & \frac{1}{β} B_{2} M_{2}^{- 1} & \frac{1}{β} I - \frac{1}{β^{2}} {\hat{S}}_{2} \end{matrix}]$

其中， $M_{1} = A_{1} + \frac{1}{α} B_{1}^{T} B_{1}$ ， $M_{2} = A_{2} + \frac{1}{β} B_{2}^{T} B_{2}$ ， ${\hat{S}}_{1} = B_{1} M_{1}^{- 1} B_{1}^{T}$ ， ${\hat{S}}_{2} = B_{2} M_{2}^{- 1} B_{2}^{T}$ 。

$\begin{matrix} P_{I D S}^{- 1} A = I - P_{I D S}^{- 1} R_{I D S} = I - α P_{2}^{- 1} P_{1}^{- 1} R_{I D S} \\ = I - [\begin{matrix} \frac{1}{α} I & 0 & 0 \\ 0 & M_{2}^{- 1} & - \frac{1}{β} M_{2}^{- 1} B_{2}^{T} \\ 0 & \frac{1}{β} B_{2} M_{2}^{- 1} & \frac{1}{β} I - \frac{1}{β^{2}} {\hat{S}}_{2} \end{matrix}] [\begin{matrix} α M_{1}^{- 1} & 0 & - M_{1}^{- 1} B_{1}^{T} \\ 0 & I & 0 \\ B_{1} M_{1}^{- 1} & 0 & I - \frac{1}{α} {\hat{S}}_{1} \end{matrix}] [\begin{matrix} 0 & - \frac{1}{α} B_{1}^{T} B_{2} & (\frac{β}{α} - 1) B_{1}^{T} \\ 0 & 0 & 0 \\ 0 & 0 & β I \end{matrix}] \\ = I - [\begin{matrix} 0 & - \frac{1}{α} M_{1}^{- 1} B_{1}^{T} B_{2} & - M_{1}^{- 1} B_{1}^{T} \\ 0 & \frac{1}{α β} M_{2}^{- 1} B_{2}^{T} {\hat{S}}_{1} B_{2} & - M_{2}^{- 1} B_{2}^{T} (I - \frac{1}{β} {\hat{S}}_{1}) \\ 0 & - \frac{1}{α β} (I - \frac{1}{β} {\hat{S}}_{2}) {\hat{S}}_{1} B_{2} & (I - \frac{1}{β} {\hat{S}}_{2}) (I - \frac{1}{β} {\hat{S}}_{1}) \end{matrix}] = I - [\begin{matrix} 0 & T_{12} \\ 0 & T_{22} \end{matrix}] \end{matrix}$ (13)

其中，

$T_{12} = [\begin{matrix} - \frac{1}{α} M_{1}^{- 1} B_{1}^{T} B_{2} & - M_{1}^{- 1} B_{1}^{T} \end{matrix}] \in ℝ^{n_{1} \times (n_{2} + m)}$

和

$T_{22} = [\begin{matrix} \frac{1}{α β} M_{2}^{- 1} B_{2}^{T} {\hat{S}}_{1} B_{2} & - M_{2}^{- 1} B_{2}^{T} (I - \frac{1}{β} {\hat{S}}_{1}) \\ - \frac{1}{α β} (I - \frac{1}{β} {\hat{S}}_{2}) {\hat{S}}_{1} B_{2} & (I - \frac{1}{β} {\hat{S}}_{2}) (I - \frac{1}{β} {\hat{S}}_{1}) \end{matrix}] \in ℝ^{(n_{2} + m) \times (n_{2} + m)}$

根据引理1， $T_{22}$ 的特征值由0和 $1 - μ_{i}$ 给出，因此，由式(13)可以看出， $P_{I D S}^{- 1} A$ 的特征值为1，其重数至少为 $n_{1} + n_{2}$ ，剩下的特征值 $μ_{i}$ 是矩阵 $Z$ 的特征值：

$Z = \frac{1}{β} ({\hat{S}}_{1} + {\hat{S}}_{2}) - \frac{α + β}{α β^{2}} {\hat{S}}_{1} {\hat{S}}_{2}$ (14)

现在，我们只需要考虑矩阵 $Z$ 的谱性质，为了方便，令 $S_{k} = B_{k} A_{k}^{- 1} B_{k}^{T}$ ，其中， $k = 1, 2$ 。则

$\begin{matrix} S_{1} = B_{1} A_{1}^{- 1} B_{1}^{T} \\ = B_{1} A_{1}^{- 1} {(I + \frac{1}{α} B_{1}^{T} B_{1} A_{1}^{- 1})}^{- 1} (I + \frac{1}{α} B_{1}^{T} B_{1} A_{1}^{- 1}) B_{1}^{T} \\ = B_{1} {(A_{1} + \frac{1}{α} B_{1}^{T} B_{1})}^{- 1} B_{1}^{T} (I + \frac{1}{α} B_{1} A_{1}^{- 1} B_{1}^{T}) \\ = {\hat{S}}_{1} (I + \frac{1}{α} S_{1}) \end{matrix}$ (15)

通过计算，由(15)式易知

${\hat{S}}_{1} = α S_{1} {(α I + S_{1})}^{- 1}$ (16)

经过类似的推导，我们可以得到

${\hat{S}}_{2} = β S_{2} {(β I + S_{2})}^{- 1}$ (17)

将上述(16)、(17)代入 $Z$ 中得到

$\begin{matrix} Z = \frac{1}{β} (α S_{1} {(α I + S_{1})}^{- 1} + β S_{2} {(β I + S_{2})}^{- 1}) - \frac{α + β}{β} S_{1} {(α I + S_{1})}^{- 1} S_{2} {(β I + S_{2})}^{- 1} \\ = \frac{1}{β} S_{1} {(α I + S_{1})}^{- 1} (α I - (α + β) S_{2} {(β I + S_{2})}^{- 1}) + S_{2} {(β I + S_{2})}^{- 1} \\ = S_{1} {(α I + S_{1})}^{- 1} (α I - S_{2}) {(β I + S_{2})}^{- 1} + S_{2} {(β I + S_{2})}^{- 1} \end{matrix}$

我们注意到 $S_{1} {(α I + S_{1})}^{- 1} = {(α I + S_{1})}^{- 1} S_{1}$ ，矩阵 $Z$ 可以进一步简化为

$\begin{matrix} Z = {(α I + S_{1})}^{- 1} (S_{1} (α I - S_{2}) + (α I + S_{1}) S_{2}) {(β I + S_{2})}^{- 1} \\ = α {(α I + S_{1})}^{- 1} (S_{1} + S_{2}) {(β I + S_{2})}^{- 1} \end{matrix}$

它相似于 $\hat{Z} = α {(β I + S_{2})}^{- 1} {(α I + S_{1})}^{- 1} (S_{1} + S_{2})$ 。因此， $Z$ 的特征值 $μ_{i}, i = 1, 2, \dots, m$ 满足以下广义特征值问题：

$(S_{1} + S_{2}) ξ_{i} = μ_{i} \cdot \frac{1}{α} (α I + S_{1}) (β I + S_{2}) ξ_{i}$ (18)

其中， $ξ_{i}$ 是 $\hat{Z}$ 的特征值 $μ_{i}$ 对应的特征向量。

证毕。

4. 改进维数分裂(IDS)预处理子的最优参数估计

改进维数分裂(IDS)预处理子的计算效率很大程度上取决于 $α$ 和 $β$ 两个参数的选取。为了使当前预处理子的预处理效果更好，我们需要选择合适的参数使得预处理子能够尽可能地接近系数矩阵，从而使预处理矩阵具有更强的谱聚集。在本文中，我们对IDS预处理子中的参数 $α$ 和 $β$ 做了全局最优。

首先，定义

$f (α, β) = {‖ R_{I D S} ‖}_{F}^{2} = t r (R_{I D S}^{T} R_{I D S})$

则我们有

$\begin{matrix} f (α, β) = {‖ - \frac{1}{α} B_{1}^{T} B_{2} ‖}_{F}^{2} + {‖ (\frac{β}{α} - 1) B_{1}^{T} ‖}_{F}^{2} + β^{2} {‖ I ‖}_{F}^{2} \\ = \frac{1}{α^{2}} t r (B_{1}^{T} B_{2} B_{2}^{T} B_{1}) + {(\frac{β}{α} - 1)}^{2} t r (B_{1}^{T} B_{1}) + β^{2} m \end{matrix}$

先对上式求一阶偏导，得到

${\begin{cases} f_{α} = - \frac{2}{α^{3}} t r (B_{1}^{T} B_{2} B_{2}^{T} B_{1}) + 2 (\frac{β}{α} - 1) (- \frac{β}{α^{2}}) t r (B_{1}^{T} B_{1}) \\ f_{β} = 2 (\frac{β}{α} - 1) \frac{1}{α} t r (B_{1}^{T} B_{1}) + 2 β m \end{cases}$ (19)

等价于

${\begin{cases} f_{α} = - \frac{2 t r (B_{1}^{T} B_{2} B_{2}^{T} B_{1})}{α^{3}} - \frac{2 β^{2} t r (B_{1}^{T} B_{1})}{α^{3}} + \frac{2 β t r (B_{1}^{T} B_{1})}{α^{2}} \\ f_{β} = \frac{2 β t r (B_{1}^{T} B_{1})}{α^{2}} - \frac{2 t r (B_{1}^{T} B_{1})}{α} + 2 β m \end{cases}$ (20)

经过简单计算，我们可以得到驻点 $(\frac{1}{\sqrt{\frac{m}{\sqrt{a m}} - \frac{m}{b}}}, \sqrt{\frac{a}{\sqrt{a m}} - \frac{a}{b}})$ ，其中， $a = t r (B_{1}^{T} B_{2} B_{2}^{T} B_{1})$ ， $b = t r (B_{1}^{T} B_{1})$ 。

再对上式(20)求二阶偏导，得到

${\begin{cases} f_{α α} = \frac{6 t r (B_{1}^{T} B_{2} B_{2}^{T} B_{1})}{α^{4}} + \frac{6 β^{2} t r (B_{1}^{T} B_{1})}{α^{4}} - \frac{4 β t r (B_{1}^{T} B_{1})}{α^{3}} \\ f_{α β} = - \frac{4 β t r (B_{1}^{T} B_{1})}{α^{3}} + \frac{2 t r (B_{1}^{T} B_{1})}{α^{2}} \\ f_{β β} = \frac{2 t r (B_{1}^{T} B_{1})}{α^{2}} + 2 m \end{cases}$ (21)

令 $f_{α α} = A$ ， $f_{α β} = B$ ， $f_{β β} = C$ ，其中， $a = t r (B_{1}^{T} B_{2} B_{2}^{T} B_{1})$ ， $b = t r (B_{1}^{T} B_{1})$ ， $α = \frac{1}{\sqrt{\frac{m}{\sqrt{a m}} - \frac{m}{b}}}$ ， $β = \sqrt{\frac{a}{\sqrt{a m}} - \frac{a}{b}}$ 。

则 $A C - B^{2} = (\frac{6 a}{α^{4}} + \frac{6 β^{2} b}{α^{4}} - \frac{4 β b}{α^{3}}) (\frac{2 b}{α^{2}} + 2 m) - {(- \frac{4 β b}{α^{3}} + \frac{2 b}{α^{2}})}^{2}$ ，经过一系列计算，

得到 $A C - B^{2} = \frac{16 a b}{α^{6}}$ 。因为 $a, b > 0$ ， $α^{6} > 0$ ，所以易知 $A C - B^{2} > 0$ 。由此可知， $A C > B^{2} > 0$ ，又因为 $C = \frac{2 b}{α^{2}} + 2 m$ ，把 $α = \frac{1}{\sqrt{\frac{m}{\sqrt{a m}} - \frac{m}{b}}}$ 代入，经过计算可知 $C = \frac{2 m b}{\sqrt{a m}}$ ，因为 $m, b > 0$ ， $\sqrt{a m} > 0$ ，所以 $C > 0$ 。根据 $A C > 0$ ，可知 $A > 0$ 。

根据二元函数极值存在定理， $A C - B^{2} > 0$ ，且 $A > 0$ ，因此函数 $f (α, β)$ 存在极小值，可以获得IDS预处理子中的准最优参数 $α$ 和 $β$ ，其表达式如下：

$α = \frac{1}{\sqrt{\frac{m}{\sqrt{a m}} - \frac{m}{b}}}$ ， $β = \sqrt{\frac{a}{\sqrt{a m}} - \frac{a}{b}}$ (22)

5. 数值实验

在本节中，我们利用计算流体力学中的一个例子来测试IDS预处理子的性能。

在实际计算中，我们分别使用DS、RDF、RSS和IDS预处理来加速广义极小残差(GMRES)方法的收敛。对于DS预处理子(3)，我们利用[16]来计算其最优参数：

$\begin{matrix} α_{D S} = \arg \min α {‖ I ‖}_{F} \cdot α {‖ I ‖}_{F} - ({‖ A_{1} ‖}_{F} + {‖ A_{2} ‖}_{F}) α + {‖ A_{1} ‖}_{F} {‖ A_{2} ‖}_{F} \\ = \frac{\sqrt{t r (A_{1}^{T} A_{1}) + 2 t r (B_{1} B_{1}^{T})} + \sqrt{t r (A_{2}^{T} A_{2}) + 2 t r (B_{2} B_{2}^{T})}}{2 (n_{1} + n_{2} + m)} \end{matrix}$ (23)

利用文献[17]来选择RDF预处理子的最优参数：

$α_{R D F} = \arg \min | \frac{1}{α} t r (S_{1} + S_{2}) - \frac{2}{α^{2}} t r (S_{1} S_{2}) |$ (24)

其中 $S_{i} = B_{i} d i a g {(A_{i})}^{- 1} B_{i}^{T}, i = 1, 2$

针对文献[18]提出的如下形式的RSS预处理子：

$P_{R S S} = \frac{1}{α} [\begin{matrix} A_{1} & 0 & 0 \\ 0 & α I & 0 \\ - B_{1} & 0 & α I \end{matrix}] [\begin{matrix} α I & 0 & B_{1}^{T} \\ 0 & A_{2} & B_{2}^{T} \\ 0 & - B_{2} & α I \end{matrix}] = [\begin{matrix} A_{1} & 0 & \frac{1}{α} A_{1} B_{1}^{T} \\ 0 & A_{2} & B_{2}^{T} \\ - B_{1} & - B_{2} & α I - \frac{1}{α} B_{1} B_{1}^{T} \end{matrix}]$

我们利用Huang提出的代数估计方法[19]，得到了RSS预处理子中的最优参数 $α_{R S S}$ 。

这里，在计算DS、RDF、RSS和IDS预处理子的最优参数时，为了减少工作量，我们使用以下公式来计算矩阵乘积的迹。

$t r (A B) = \sum_{i, j = 1}^{n} {(A \circ B^{T})}_{i j}$ ， $A \in ℝ^{n \times m}$ ， $B \in ℝ^{m \times n}$

其中 $\circ$ 表示哈达玛积。

在下面实验中，我们选取向量 $x^{(0)} = {(0, 0, \dots, 0)}^{T} \in ℝ^{n_{1} + n_{2} + m}$ 为初始向量，利用迭代步数(缩写为IT，单位：次)和算法达到收敛时所需要的总运行时间(缩写为CPU，单位：秒)来比较GMRES和四种预处理GMRES方法的数值结果，一旦当前迭代 $x^{(k)}$ 满足

$\frac{{‖ b - A x^{(k)} ‖}_{2}}{{‖ b - A x^{(0)} ‖}_{2}} \leq 10^{- 6}$

或者迭代次数超过2500，迭代就终止。下面表中“+”表示不收敛或者迭代步数超过2500的情况。

本文的实验都在电脑内存为16.0GB，CPU型号为Intel(R) Core(TM) i5-12500CPU @ 3.10 GHz，操作系统为Windows 10的MATLAB(R2022b)上完成。

例1 [1]考虑不可压缩Navier-Stokes方程的解。

$\begin{array}{l} \frac{\partial u}{\partial t} - v \cdot u + (u \cdot \nabla) u + \nabla p = f, on Ω \times (0, T] . \\ d i v u = 0, on Ω \times [0, T], \\ u = g, on \partial Ω \times [0, T], \\ u (x, 0) = u_{0} (x), on Ω \end{array}$

其中 $Ω \in ℝ^{2}$ 是一个开有界域，时间间隔为 $[0, T]$ ，f为外力场，g和u₀是初始边界数据， $u = u (x, t)$ 为速度场， $p = p (x, t)$ 为压力场，v为运动粘度，与雷诺数成反比。∆、 $\nabla$ 和div符号分别表示拉普拉斯算子、梯度和散度。我们使用IFISS软件包[20]将Ω划分，对二维漏盖驱动空腔问题进行离散化，实验中涉及的网格类型为stretched，选取的有限元分别为：Q2-Q1有限元和Q2-P1有限元，网格大小和粘度常数分别设置为16 × 16，32 × 32，64 × 64，128 × 128和 $v = 10^{- 1}, 10^{- 2}, 10^{- 3}, 10^{- 4}$ 。更多细节可参见文[1]。

针对离散化的3 × 3块线性方程组，根据上面的分析，我们分别计算出DS、RDF、RSS和IDS预处理子的准最优参数值。表1和表2记录了Q2-Q1、Q2-P1有限元离散化的DS、RDF、RSS和IDS预处理子的准最优参数值。

Table 1. The optimal parameter-values of the DS, RDF, RSS and IDS preconditioners for the Q2-Q1 FEM discretization

表1. Q2-Q1有限元离散化的DS、RDF、RSS和IDS预处理子的准最优参数

v	Grid	$α_{D S}$	$α_{R D F}$	$α_{R S S}$	$α_{I D S}$	$β$
10⁻¹	16 × 16	0.0194	0.0523	0.7065	0.2482	0.0589
	32 × 32	0.0121	0.0173	0.9186	0.1295	0.0352
	64 × 64	0.0084	0.0056	1.2298	0.0678	0.0198
	128 × 128	0.0060	0.0018	1.6364	0.0356	0.0108
10⁻²	16 × 16	0.0125	0.5234	0.5280	0.2482	0.0589
	32 × 32	0.0048	0.1734	0.3838	0.1295	0.0352
	64 × 64	0.0019	0.0558	0.3016	0.0678	0.0198
	128 × 128	0.0009	0.0176	0.3602	0.0356	0.0108
10⁻³	16 × 16	0.0124	5.2480	0.5741	0.2482	0.0589
	32 × 32	0.0047	1.7340	0.4641	0.1295	0.0352
	64 × 64	0.0017	0.5580	0.3488	0.0678	0.0198
	128 × 128	0.0006	0.1764	0.2196	0.0356	0.0108
10⁻⁴	16 × 16	0.0124	15.3091	0.5794	0.2482	0.0589
	32 × 32	0.0047	17.3865	0.4763	0.1295	0.0352
	64 × 64	0.0017	5.5803	0.3786	0.0678	0.0198
	128 × 128	0.0006	1.7642	0.2908	0.0356	0.0108

Table 2. The optimal parameter-values of the DS, RDF, RSS and IDS preconditioners for the Q2-P1 FEM discretization

表2. Q2-P1有限元离散化的DS、RDF、RSS和IDS预处理子的准最优参数

v	Grid	$α_{D S}$	$α_{R D F}$	$α_{R S S}$	$α_{I D S}$	$β$
10⁻¹	16 × 16	0.0176	0.1914	1.1738	0.4345	0.1037
	32 × 32	0.0102	0.0605	1.5627	0.2220	0.0588
	64 × 64	0.0069	0.0191	2.2353	0.1148	0.0322
	128 × 128	0.0050	0.0060	3.1209	0.0600	0.0173
10⁻²	16 × 16	0.0122	1.9136	1.1448	0.4345	0.1037
	32 × 32	0.0043	0.6048	0.7645	0.2220	0.0588
	64 × 64	0.0016	0.1912	0.4950	0.1148	0.0322
	128 × 128	0.0007	0.0601	0.5666	0.0600	0.0173
10⁻³	16 × 16	0.0121	19.1486	1.2566	0.4345	0.1037
10⁻³	32 × 32	0.0042	6.0477	1.0090	0.2220	0.0588
10⁻³	64 × 64	0.0015	1.9119	0.7422	0.1148	0.0322
10⁻³	128 × 128	0.0005	0.6015	0.4219	0.0600	0.0173
10⁻⁴	16 × 16	0.0121	200.2231	1.2684	0.4345	0.1037
	32 × 32	0.0042	60.5177	1.0391	0.2220	0.0588
	64 × 64	0.0015	19.1188	0.8219	0.1148	0.0322
	128 × 128	0.0005	6.0148	0.6267	0.0600	0.0173

基于上述这些参数值，我们利用四种预处理GMRES方法对线性方程组(1)的解进行近似。表3和表4分别记录了Q2-Q1、Q2-P1有限元离散化的预处理GMRES方法数值结果。从表3和表4中，我们可以看出，当v = 10⁻¹时，RDF预处理子的效果要比DS、RSS和IDS预处理子好。不过，随着v的减小，我们可以观察到，IDS预处理子的效果是最佳的，因为IDS预处理GMRES方法的迭代步数和CPU时间都比DS、RDF和RSS预处理GMRES方法少得多。

Table 3. Numerical results of the PGMRES methods for the Q2-Q1 FEM discretization

表3. Q2-Q1有限元离散化的预处理GMRES方法数值结果

		16 × 16		32 × 32		64 × 64		128 × 128
v	Method	IT	CPU	IT	CPU	IT	CPU	IT	CPU
10⁻¹	DS	22	0.0052	34	0.0408	88	1.0870	221	21.4641
	RDF	19	0.0054	29	0.0383	48	0.6628	80	8.6092
	RSS	46	0.0127	89	0.0765	181	1.9199	344	22.4499
	IDS	28	0.0105	51	0.0552	100	1.7292	194	19.2940
10⁻²	DS	60	0.0395	94	0.3356	124	2.7392	136	17.3263
	RDF	24	0.0055	37	0.0385	54	0.7256	90	9.4049
	RSS	41	0.0117	83	0.0812	172	1.8069	397	25.5556
	IDS	23	0.0123	33	0.0411	55	0.4100	109	11.2747
10⁻³	DS	190	0.1782	325	1.4786	509	23.8302	690	190.3037
	RDF	77	0.0602	132	0.4574	210	4.3294	273	26.9543
	RSS	57	0.0326	139	0.3519	318	4.3385	618	44.4922
	IDS	42	0.0316	65	0.3076	74	2.0447	98	11.8734
10⁻⁴	DS	428	0.7478	1147	19.1632	1810	188.9457	+	+
	RDF	87	0.0818	318	2.4150	1065	71.5925	1887	654.4964
	RSS	68	0.0503	211	0.9038	677	25.0102	1549	309.1655
	IDS	65	0.0614	140	1.1032	210	11.9102	254	89.7660

Table 4. Numerical results of the PGMRES methods for the Q2-P1 FEM discretization

表4. Q2-P1有限元离散化的预处理GMRES方法数值结果

		16 × 16		32 × 32		64 × 64		128 × 128
v	Method	IT	CPU	IT	CPU	IT	CPU	IT	CPU
10⁻¹	DS	31	0.0082	32	0.0281	50	0.1768	118	2.7674
	RDF	18	0.0037	29	0.0144	44	0.1097	63	1.4456
	RSS	38	0.0098	82	0.0587	150	0.4883	253	6.4549
	IDS	23	0.0071	43	0.0282	79	0.2310	141	3.4027
10⁻²	DS	99	0.0518	135	0.1764	150	0.8612	154	3.8271
	RDF	28	0.0047	43	0.0266	62	0.1614	90	2.0641
	RSS	47	0.0122	89	0.1764	158	0.5270	294	7.8066
	IDS	29	0.0117	29	0.0295	42	0.1514	80	1.8566
10⁻³	DS	283	0.2511	520	1.9303	716	11.9321	858	42.8383
	RDF	147	0.0922	236	0.4707	346	2.9161	490	16.5912
	RSS	92	0.0547	205	0.3420	431	4.0651	853	40.0421
	IDS	69	0.0400	79	0.1226	84	0.4923	82	2.0019
10⁻⁴	DS	413	0.6073	1232	20.7233	+	+	+	+
	RDF	220	0.1827	846	7.3474	2310	238.7753	+	+
	RSS	123	0.0784	430	1.4620	1219	41.9477	+	+
	IDS	131	0.0885	207	0.5668	261	4.1543	244	14.8608

6. 结论

基于3 × 3块系数矩阵的维数分裂(DS)，本文提出了一种新的改进维数分裂(IDS)预处理子，并讨论了IDS预处理矩阵的谱性质。在计算准最优参数之后，通过二维不可压缩Navier-Stokes方程这一个重要算例的数值实验，验证了IDS预处理子与其他预处理子相比具有很大的优越性。

参考文献

[1]	Elman, H.C., Silvester, D.J. and Wathen, A.J. (2014) Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics. Oxford University Press.
[2]	Bai, Z. (2010) Block Preconditioners for Elliptic PDE-Constrained Optimization Problems. Computing, 91, 379-395. https://doi.org/10.1007/s00607-010-0125-9
[3]	Bai, Z. and Lu, K. (2021) Optimal Rotated Block-Diagonal Preconditioning for Discretized Optimal Control Problems Constrained with Fractional Time-Dependent Diffusive Equations. Applied Numerical Mathematics, 163, 126-146. https://doi.org/10.1016/j.apnum.2021.01.011
[4]	de Sturler, E. and Liesen, J. (2005) Block-diagonal and Constraint Preconditioners for Nonsymmetric Indefinite Linear Systems. Part I: Theory. SIAM Journal on Scientific Computing, 26, 1598-1619. https://doi.org/10.1137/s1064827502411006
[5]	Cao, Z. (2007) Positive Stable Block Triangular Preconditioners for Symmetric Saddle Point Problems. Applied Numerical Mathematics, 57, 899-910. https://doi.org/10.1016/j.apnum.2006.08.001
[6]	Zhou, S., Yang, A. and Wu, Y. (2016) A Relaxed Block-Triangular Splitting Preconditioner for Generalized Saddle-Point Problems. International Journal of Computer Mathematics, 94, 1609-1623. https://doi.org/10.1080/00207160.2016.1226500
[7]	Chaparpordi, S.H.A., Beik, F.P.A. and Salkuyeh, D.K. (2018) Block Triangular Preconditioners for Stabilized Saddle Point Problems with Nonsymmetric (1, 1)-Block. Computers & Mathematics with Applications, 76, 1544-1553. https://doi.org/10.1016/j.camwa.2018.07.006
[8]	Aslani, H. and Salkuyeh, D.K. (2023) A Block Triangular Preconditioner for a Class of Three-By-Three Block Saddle Point Problems. Japan Journal of Industrial and Applied Mathematics, 40, 1015-1030. https://doi.org/10.1007/s13160-022-00561-8
[9]	Cao, Z.H. (2006) A Class of Constraint Preconditioners for Nonsymmetric Saddle Point Matrices. Numerische Mathematik, 103, 47-61. https://doi.org/10.1007/s00211-006-0675-0
[10]	Bai, Z., Ng, M.K. and Wang, Z. (2009) Constraint Preconditioners for Symmetric Indefinite Matrices. SIAM Journal on Matrix Analysis and Applications, 31, 410-433. https://doi.org/10.1137/080720243
[11]	Benzi, M., Golub, G.H. and Liesen, J. (2005) Numerical Solution of Saddle Point Problems. Acta Numerica, 14, 1-137. https://doi.org/10.1017/s0962492904000212
[12]	Benzi, M. and Guo, X. (2011) A Dimensional Split Preconditioner for Stokes and Linearized Navier-Stokes Equations. Applied Numerical Mathematics, 61, 66-76. https://doi.org/10.1016/j.apnum.2010.08.005
[13]	Benzi, M., Ng, M., Niu, Q. and Wang, Z. (2011) A Relaxed Dimensional Factorization Preconditioner for the Incompressible Navier-Stokes Equations. Journal of Computational Physics, 230, 6185-6202. https://doi.org/10.1016/j.jcp.2011.04.001
[14]	Ren, B., Chen, F. and Wang, X. (2022) Improved Splitting Preconditioner for Double Saddle Point Problems Arising from Liquid Crystal Director Modeling. Numerical Algorithms, 91, 1363-1379. https://doi.org/10.1007/s11075-022-01305-y
[15]	Horn, R.A. and Johnson, C.R. (2012). Matrix Analysis. 2nd Edition, Cambridge University Press. https://doi.org/10.1017/cbo9781139020411
[16]	Ren, Z. and Cao, Y. (2015) An Alternating Positive-Semidefinite Splitting Preconditioner for Saddle Point Problems from Time-Harmonic Eddy Current Models. IMA Journal of Numerical Analysis, 36, 922-946. https://doi.org/10.1093/imanum/drv014
[17]	Benzi, M., Deparis, S., Grandperrin, G. and Quarteroni, A. (2016) Parameter Estimates for the Relaxed Dimensional Factorization Preconditioner and Application to Hemodynamics. Computer Methods in Applied Mechanics and Engineering, 300, 129-145. https://doi.org/10.1016/j.cma.2015.11.016
[18]	Tan, N., Huang, T. and Hu, Z. (2012) A Relaxed Splitting Preconditioner for the Incompressible Navier‐Stokes Equations. Journal of Applied Mathematics, 2012, Article ID: 402490. https://doi.org/10.1155/2012/402490
[19]	Huang, Y. (2014) A Practical Formula for Computing Optimal Parameters in the HSS Iteration Methods. Journal of Computational and Applied Mathematics, 255, 142-149. https://doi.org/10.1016/j.cam.2013.01.023
[20]	Elman, H.C., Ramage, A. and Silvester, D.J. (2007) Algorithm 866: IFISS, a Matlab Toolbox for Modelling Incompressible Flow. ACM Transactions on Mathematical Software, 33, Article 14. https://doi.org/10.1145/1236463.1236469

为你推荐

友情链接