本文主要考虑Sylvester矩阵方程

$ \boldsymbol{A} \boldsymbol{X}+\boldsymbol{X} \boldsymbol{B}=\boldsymbol{C}, $

(1)

其中A∈R^n×n，B∈R^m×m，C∈R^n×m，且系数矩阵满足下列条件之一：

(Ⅰ)A, B是非奇异M-矩阵，C≥0；

(Ⅱ)A, B是不可约M-矩阵且其中至少有一个是非奇异矩阵，C≥0。

方程(1)称为M-矩阵Sylvester方程(MSE)^[1]。当系数矩阵没有结构要求时，称为一般的Sylvester矩阵方程；当A=B^T时，称为Lyapunov矩阵方程。此类广泛应用控制理论领域^[2]，海洋平台减振系统中反馈增益模型需要求解Sylvester方程^[3]。

求解方程(1)的常用数值算法有不动点迭代法^{[1, 4]}、牛顿法^[2]、交替方向加倍算法^[5]、子空间迭代法^[6]和保结构加倍算法^[2]等。当系数矩阵A, B阶数n, m很大时，MSE的求解速度会非常慢，且占用大量内存。P. Benner等针对大型矩阵代数Riccati方程，结合牛顿法和交替方向隐式迭代的思想提出了低秩牛顿交替方向法^[7]；D. Kressner等针对线性矩阵方程在低秩情况下提出了简单而快速的数值方法，基于这种思想，对系数矩阵具有特定的低秩分级结构提出了分而治之方法^[8]。

本文将考虑系数矩阵具有分级非对角低秩(Hierarchically off-diagonal low-rank，HODLR)结构M-矩阵Sylvester方程，提出了分而治之方法，把方程化为低阶MSE和低秩MSE分别求解。具体思路如下：

设矩阵A, B, C可以分裂为

$ \boldsymbol{A}=\boldsymbol{A}_{0}+\boldsymbol{\delta} \boldsymbol{A}, \boldsymbol{B}=\boldsymbol{B}_{0}+\boldsymbol{\delta} \boldsymbol{B}, \boldsymbol{C}=\boldsymbol{C}_{0}+\boldsymbol{\delta} \boldsymbol{C}, $

(2)

其中δA, δB, δC是低秩矩阵。假设X₀是矩阵方程

$ \boldsymbol{A}_{0} \boldsymbol{X}_{0}+\boldsymbol{X}_{0} \boldsymbol{B}_{0}=\boldsymbol{C}_{0} $

(3)

的解，求δX使得X=X₀+δX是方程(4)的解。

$ \begin{gathered} \left(\boldsymbol{A}_{0}+\boldsymbol{\delta} \boldsymbol{A}\right)\left(\boldsymbol{X}_{0}+\boldsymbol{\delta} \boldsymbol{X}\right)+ \\ \left(\boldsymbol{X}_{0}+\boldsymbol{\delta} \boldsymbol{X}\right)\left(\boldsymbol{B}_{0}+\boldsymbol{\delta B}\right)=\left(\boldsymbol{C}_{0}+\boldsymbol{\delta} \boldsymbol{C}\right) \end{gathered} $

(4)

由式(3)和(4)得

$ \left(\boldsymbol{A}_{0}+\boldsymbol{\delta} \boldsymbol{A}\right) \boldsymbol{\delta} \boldsymbol{X}+\boldsymbol{\delta} \boldsymbol{X}\left(\boldsymbol{B}_{0}+\boldsymbol{\delta} \boldsymbol{B}\right)=\overset\frown{\boldsymbol{C}}, $

(5)

其中$\hat{\mathit{\boldsymbol{C}}}=\mathit{\boldsymbol{\delta}} \mathit{\boldsymbol{C}}-\mathit{\boldsymbol{\delta}} \mathit{\boldsymbol{A}} \mathit{\boldsymbol{X}}_{0}-\mathit{\boldsymbol{X}}_{0} \mathit{\boldsymbol{\delta}} \mathit{\boldsymbol{B}}$，易知

$ {\rm{rank}}\left( {\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over C} }}} \right) \le {\rm{rank}}\left( {\mathit{\boldsymbol{\delta C}}} \right) + {\rm{rank}}\left( {\mathit{\boldsymbol{\delta A}}} \right) + {\rm{rank}}\left( {\mathit{\boldsymbol{\delta B}}} \right), $

故式(5)是右端矩阵为低秩矩阵的方程，当式(5)是M-矩阵Sylvester方程时，我们简称为更新MSE问题。当式(3)中已知矩阵都是块对角阵时，式(3)将转化为多个低阶MSE问题的求解。

注1    求解方程(1)转化为求解方程(3)和(5)。

1 M-矩阵及相关性质

本节主要介绍M-矩阵的定义与性质，以及对应的矩阵方程的可解性。

定义1^[1]    设A∈R^n×n，如果所有非对角元素都小于等于0即对∀i≠j, a_ij≤0，称A是Z-矩阵。如果存在数s和非负矩阵N，使得A=sI-N且s≥ρ(N)，就称矩阵A是M-矩阵。当s>ρ(N)时，称A是非奇异M-矩阵；当s=ρ(N)时，称A是奇异M-矩阵。

关于M-矩阵有如下性质：

引理1^[1]    (Ⅰ)设A是Z-矩阵，且存在v≥0，使得Av≥0，则A是M-矩阵；

(Ⅱ)M-矩阵的主子矩阵仍然是M-矩阵；

(Ⅲ)如果A是非奇异M-矩阵，且Z-矩阵B满足B≥A，则B是非奇异M-矩阵；

(Ⅳ)如果A是不可约奇异M-矩阵，且Z-矩阵B满足B≥A，当B≠A时，则B是非奇异M-矩阵；

(Ⅴ)如果A是非奇异的M-矩阵或不可约奇异M-矩阵，分块如下：

$ \boldsymbol{A}=\left[\begin{array}{ll} \boldsymbol{A}_{11} & \boldsymbol{A}_{12} \\ \boldsymbol{A}_{21} & \boldsymbol{A}_{22} \end{array}\right], $

其中A₁₁和A₂₂是方阵，则A₁₁和A₂₂都是非奇异的M-矩阵。

引理2^[1]    设A∈R^n×n是Z-矩阵，则下列几个命题等价：

(Ⅰ)A是非奇异M-矩阵；

(Ⅱ)A^-1≥0；

(Ⅲ)存在v∈Rⁿ满足v>0, 使得Av>0；

(Ⅳ)A的全部特征值都有正的实部。

关于矩阵方程(1)有如下定理。

定理1    (Ⅰ)设A和B都是非奇异M-矩阵，且C≥0，则方程(1)存在唯一非负解。

(Ⅱ)设A和B是不可约M-矩阵且至少一个是非奇异矩阵，C≥0，则方程(1)存在唯一非负解。

证明    由A和B是非奇异M-矩阵或是不可约M-矩阵且至少一个是非奇异矩阵，可得

$ \boldsymbol{P}=I_{m} \otimes \boldsymbol{A}+\boldsymbol{B}^{\mathrm{T}} \otimes I_{n} $

是非奇异M-矩阵，因此方程(1)有唯一解。

由引理2中(Ⅱ)得P^-1≥0，注意到C≥0，所以解是非负的。

定义2^[8]    若A∈R^n×n，分裂p层记为Γ_p, 则：

(Ⅰ)对于k∈N, 称A为Γ_p, k)-HODLR矩阵当且仅当每一个非对角块A(I_i^l, I_j^l)和A(I_j^l, I_i^l)的秩不超过k，其中(i, j)={(2, 1), (4, 3), …, (2^l, 2^l-1)}，l=1, 2, …, p。

(Ⅱ)对于Γ_p，矩阵A的HODLR秩指：最小正整数k使得A是Γ_p, k)-HODLR矩阵。

2 低秩更新

本节主要考虑矩阵C是非负低秩时，MSE的求解方法。对于系数矩阵是小型矩阵或稠密矩阵情形，已有多种求解算法，例如矩阵分解法、Hoskins迭代法、交替方向法等^[2]。文献[4]中提出了Smith方法，但Smith方法中的参数不易选择，文献[9]提出了求解MSE方程的交替方向Smith方法(Alternating directional Smith method，ADSM)，迭代序列是非负单调上升的，且平方收敛。ADSM是Smith方法的改进。具体算法如下：

算法 1[9]：ADSM(A, B, C) % 求解AX+XB=C
1.令$\alpha=\max _{1 \leqslant i \leqslant n} a_{i i}, \beta=\max _{1 \leqslant i \leqslant m} b_{i i}$。 2.计算     E0=(B-βI)(B+αI)-1,     F0=(A+βI)-1(A-αI),     X0=(α+β)(A+βI)-1C(B+αI)-1。 3.对k=1, 2, …, 重复下列步骤直到矩阵序列Xk收敛     Xk+1=Xk+FkXkEk,     Ek+1=Ek2, Fk+1=Fk2。

对于方程(5)，给定如下形式的低秩分解

$ \boldsymbol{\delta} \boldsymbol{A}=\boldsymbol{U}_{\boldsymbol{A}} \boldsymbol{V}_{\boldsymbol{A}}^{\mathrm{T}}, \boldsymbol{\delta} \boldsymbol{B}=\boldsymbol{U}_{\boldsymbol{B}} \boldsymbol{V}_{\boldsymbol{B}}^{\mathrm{T}}, \boldsymbol{\delta} \boldsymbol{C}=\boldsymbol{U}_{\boldsymbol{C}} \boldsymbol{V}_{\boldsymbol{C}}^{\mathrm{T}}。$

则

$ \begin{gathered} \overset\frown{\boldsymbol{C}}=\boldsymbol{\delta} \boldsymbol{C}-\boldsymbol{\delta} \boldsymbol{A} \boldsymbol{X}_{0}-\boldsymbol{X}_{0} \boldsymbol{\delta} \boldsymbol{B}= \\ {\left[\boldsymbol{U}_{\boldsymbol{C}},-\boldsymbol{U}_{\boldsymbol{A}},-\boldsymbol{X}_{0} \boldsymbol{U}_{\boldsymbol{B}}\right]\left[\boldsymbol{V}_{\boldsymbol{C}}, \boldsymbol{X}_{0}^{\mathrm{T}} \boldsymbol{V}_{\boldsymbol{A}}, \boldsymbol{V}_{\boldsymbol{B}}\right]^{\mathrm{T}}=\boldsymbol{U} \boldsymbol{V}^{\mathrm{T}},} \end{gathered} $

式中U, V均为低秩且令

$ s=\operatorname{rank}(\boldsymbol{\delta} \boldsymbol{A})+\operatorname{rank}(\boldsymbol{\delta} \boldsymbol{B})+\operatorname{rank}(\boldsymbol{\delta} \boldsymbol{C})。$

则方程(5)可以表示为:

$ \boldsymbol{A} \boldsymbol{\delta} \boldsymbol{X}+\boldsymbol{\delta} \boldsymbol{X} \boldsymbol{B}=\boldsymbol{U} \boldsymbol{V}^{\mathrm{T}} \text { 。} $

(6)

式中：A=A₀+δA, B=B₀+δB, U∈R^n×s, V∈R^m×s, s$ \ll $min{n, m}。低秩近似解X_k=Z_k(Y_k)^T可以直接利用低秩ADI方法求得^{[5, 10]}，但迭代中的参数需要事先给定，且最优参数不易获得。本文将使用文献[8] 中的扩展Krylov子空间方法，此类算法不需要事先选定参数。基本思路是利用块Arnoldi过程得到扩展Krylov子空间：

$ \begin{aligned} \mathscr{U}_{t}: &=\varepsilon K_{t}(\boldsymbol{A}, \boldsymbol{U})=\\ &\left\{\boldsymbol{U}, \boldsymbol{A}^{-1} \boldsymbol{U}, \boldsymbol{A} \boldsymbol{U}, \boldsymbol{A}^{-2} \boldsymbol{U}, \cdots, \boldsymbol{A}^{t-1} \boldsymbol{U}, \boldsymbol{A}^{-t} \boldsymbol{U}\right\}, \end{aligned} $

$ \begin{aligned} &\ \ \ \ \ \ \ \ \ \ \mathscr{V}_{t}:=\varepsilon K_{t}\left(\boldsymbol{B}^{\mathrm{T}}, \boldsymbol{V}\right)= \\ &\left\{\boldsymbol{V},\left(\boldsymbol{B}^{\mathrm{T}}\right)^{-1} \boldsymbol{V}, \boldsymbol{B}^{\mathrm{T}} \boldsymbol{V},\left(\boldsymbol{B}^{\mathrm{T}}\right)^{-2} \boldsymbol{V}, \cdots,\left(\boldsymbol{B}^{\mathrm{T}}\right)^{t-1} \boldsymbol{V},\left(\boldsymbol{B}^{\mathrm{T}}\right)^{-t} \boldsymbol{V}\right\}, \end{aligned} $

式中2ts < min{m, n}。分别取$\mathscr{U}_t $和$ \mathscr{V}_t$的标准正交基矩阵为U_t和V_t，得到低阶矩阵方程：

$ \widetilde{\boldsymbol{A}} \boldsymbol{X}_{t}+\boldsymbol{X}_{t} \widetilde{\boldsymbol{B}}=\widetilde{\boldsymbol{C}}, $

其中$\tilde{\mathit{\boldsymbol{A}}}=\mathit{\boldsymbol{U}}_{t}^{\mathrm{T}} \mathit{\boldsymbol{A}} \mathit{\boldsymbol{U}}_{t}, \tilde{\mathit{\boldsymbol{B}}}=\mathit{\boldsymbol{V}}_{t}^{\mathrm{T}} \mathit{\boldsymbol{B}} \mathit{\boldsymbol{V}}_{t}, \tilde{\mathit{\boldsymbol{C}}}=\mathit{\boldsymbol{U}}_{t}^{\mathrm{T}} \mathit{\boldsymbol{C}} \mathit{\boldsymbol{V}}_{t}$。

具体算法如下：

算法 2 [8]：LRMSE(A，B，U，V) %求解AδX+δXB=UVT

1.$\left[ {{U_1}, - } \right] = {\mathop{\rm qr}\nolimits} \left( {\left[ {U, {A^{ - 1}}} \right]} \right), \left[ {{V_1}, - } \right] = {\mathop{\rm qr}\nolimits} \left( {\left[ {V, {B^{ - 1}}} \right]} \right)$。 2.对k=1, 2, …, 计算     $\widetilde{A}=U_{k}^{\mathrm{T}} A U_{k}, \widetilde{B}=V_{k}^{\mathrm{T}} B V_{k}, \widetilde{C}=U_{k}^{\mathrm{T}} U V^{\mathrm{T}} V_{k}$     Xk←利用稠密方法计算$\widetilde{A} X_{k}+X_{k} \widetilde{B}=\widetilde{C}$的解，     判断Xk是否收敛，如果收敛，令     $\widetilde{X}=U_{k}^{\mathrm{T}} X_{k} V_{k}$，停止。 否则，令     $U_{k}=\left[U^{0}, U^{+}, U^{-}\right], V_{k}=\left[\mathit{\boldsymbol{V}}^{0}, V^{+}, V^{-}\right], $ 计算     $\tilde{U}=\left[A U^{+}, A^{-1} U^{-}\right], \tilde{V}=\left[B V^{+}, B^{-1} V^{-}\right], $     $\widetilde{U} \leftarrow \widetilde{U}-U_{k} U_{k}^{\mathrm{T}} \widetilde{U}, \tilde{V} \leftarrow \widetilde{U}-V_{k} V_{k}^{\mathrm{T}} \tilde{V}, $     $[\widetilde{U}, -]=\operatorname{qr}(\widetilde{U}), [\tilde{V}, -]=\operatorname{qr}(\widetilde{V}), $     $U_{k+1}=\left[U_{k}, \tilde{U}\right], V_{k+1}=\left[V_{k}, \widetilde{V}\right]_{。}$

注2    算法2中利用稠密方法求解的矩阵方程$\mathit{\boldsymbol{\tilde A}}{\mathit{\boldsymbol{X}}_k} + {\mathit{\boldsymbol{X}}_k}\mathit{\boldsymbol{\tilde B}} = \mathit{\boldsymbol{\tilde C}}$一般不是MSE形式，但仍可采用ADSM方法求解。

3 分而治之法

设方程(1)的系数矩阵具有分级非对角低秩矩阵结构，首先进行如下分裂：

$ \boldsymbol{A}=\boldsymbol{A}_{0}+\boldsymbol{\delta} \boldsymbol{A}, \boldsymbol{B}=\boldsymbol{B}_{0}+\boldsymbol{\delta} \boldsymbol{B}, \boldsymbol{C}=\boldsymbol{C}_{0}+\boldsymbol{\delta} \boldsymbol{C}。$

(7)

式中: A₀, B₀, C₀是块对角矩阵；δA, δB, δC是低秩矩阵；此时方程(1)可化为方程(3)和(5)。当对角块矩阵仍具HODLR结构时，将递归地进行分裂。具体递归分裂的层数将根据矩阵的HODLR结构、计算机的存储能力和计算速度确定，相关内容见文献[11]。图 1展示了HODLR结构, 对角块黑色部分为低阶矩阵, 其余为低秩矩阵。具体如下：

$ \boldsymbol{A}=\left[\begin{array}{ll} \boldsymbol{A}_{11} & \boldsymbol{A}_{12} \\ \boldsymbol{A}_{21} & \boldsymbol{A}_{22} \end{array}\right]。$

式中: A₁₂, A₂₁为低秩矩阵，A₁₁, A₂₂仍具有HODLR结构特征，对主对角块矩阵继续分裂，直至最小块阶数满足要求后停止。

设方程(1)的系数矩阵A, B, C都是HODLR矩阵，则类似文献[8]，得出MSE的分而治之方法。对系数矩阵进行如下分裂：

$ \begin{array}{c} \boldsymbol{A}=\left[\begin{array}{cc} \boldsymbol{A}_{11} & 0 \\ 0 & \boldsymbol{A}_{22} \end{array}\right]+\left[\begin{array}{cc} 0 & \boldsymbol{U}_{1} \boldsymbol{V}_{2}^{\mathrm{T}} \\ \boldsymbol{U}_{2} \boldsymbol{V}_{1}^{\mathrm{T}} & 0 \end{array}\right]=\\ \left[\begin{array}{cc} \boldsymbol{A}_{11} & 0 \\ 0 & \boldsymbol{A}_{22} \end{array}\right]+\boldsymbol{U}_{\boldsymbol{A}} \boldsymbol{V}_{\boldsymbol{A}}^{\mathrm{T}}。\end{array} $

(8)

式中：

$ \boldsymbol{U}_{\boldsymbol{A}}=\left[\begin{array}{ll} \boldsymbol{U}_{1} & \\ & \boldsymbol{U}_{2} \end{array}\right], \boldsymbol{V}_{\boldsymbol{A}}=\left[\begin{array}{ll} & \boldsymbol{V}_{1} \\ \boldsymbol{V}_{2} & \end{array}\right]。$

类似地分裂B和C：

$ \boldsymbol{B} =\left[\begin{array}{cc} \boldsymbol{B}_{11} & 0 \\ 0 & \boldsymbol{B}_{22} \end{array}\right]+\boldsymbol{U}_{\boldsymbol{B}} \boldsymbol{V}_{\boldsymbol{B}}^{\mathrm{T}}, $

(9)

$ \boldsymbol{C} =\left[\begin{array}{cc} \boldsymbol{C}_{11} & 0 \\ 0 & \boldsymbol{C}_{22} \end{array}\right]+\boldsymbol{U}_{\boldsymbol{C}} \boldsymbol{V}_{\boldsymbol{C}}^{\mathrm{T}}。$

(10)

注3    分块后的子矩阵需要保证对应的矩阵方程有意义。若对角块A_ii, B_ii, C_ii, i=1, 2仍是HODLR矩阵，则继续进行分裂。低秩部分对应的更新MSE的近似解将用低秩算法求解，如Krylov子空间方法、低秩Smith法、低秩ADSM法等。若每一递归层均有唯一的非负解，则方程(1)有唯一的非负解。

下列定理可以保证解存在唯一性。

定理2    设A和B是非奇异M-矩阵或是不可约M-矩阵且至少一个是非奇异矩阵，C≥0，矩阵分裂由方程(8)、(9)和(10)得到，则下列结论成立：

(Ⅰ)矩阵方程

$ \boldsymbol{A}_{i i} \boldsymbol{X}_{i i}+\boldsymbol{X}_{i i} \boldsymbol{B}_{i i}=\boldsymbol{C}_{i i}, i=1,2 $

(11)

是M-矩阵Sylvester方程，且有唯一非负解X₁₁, X₂₂，即X₀=diag(X₁₁, X₂₂)≥0是式(3)的非负解；

(Ⅱ)此时得到的矩阵方程(5)是M-矩阵Sylvester方程，且有唯一非负解δX；

(Ⅲ)方程(1)的非负解为X=X₀+δX；

(Ⅳ)分裂k次时解为

$ \boldsymbol{X}^{(k)}=\boldsymbol{X}_{0}^{(k)}+\boldsymbol{\delta} \boldsymbol{X}^{(k)}+\boldsymbol{\delta} \boldsymbol{X}^{(k-1)}+\cdots+\boldsymbol{\delta} \boldsymbol{X}^{(1)}, $

其中，

$ \boldsymbol{X}_{0}^{(k)}=\left[\begin{array}{llll} \boldsymbol{X}_{11}^{(k)} & & & \\ & \boldsymbol{X}_{22}^{(k)} & & \\ & & \ddots & \\ & & & \boldsymbol{X}_{2 ^k, 2^ k}^{(k)} \end{array}\right], $

$ \boldsymbol{\delta} \boldsymbol{X}^{(k)}=\left[\begin{array}{ccc} &\boldsymbol{\delta} \boldsymbol{X}_{12}^{(k)} &&&&&& \\ \boldsymbol{\delta} \boldsymbol{X}_{21}^{(k)} &&&&&&&\\ &&& \boldsymbol{\delta} \boldsymbol{X}_{34}^{(k)} &&& &\\ && \boldsymbol{X} _ { 4 3 } ^ { ( k ) } &&&&& \\ &&&&& \ddots && \\ &&&& \ddots &&& \\ &&&&&&& \boldsymbol{\delta} \boldsymbol{X}_{2 ^{k}-1,2^k}^{(k)}\\ &&&&&&\boldsymbol{\delta} \boldsymbol{X}_{2 ^k, 2 ^{k}-1}^{(k)}& \end{array}\right]。$

证明     (Ⅰ)由引理1中(Ⅴ)可知，A₁₁, B₁₁, A₂₂, B₂₂都是非奇异M-矩阵。由分裂方程(10)可知C₁₁, C₂₂都是非负矩阵，所以方程(11)是MSE，故有唯一非负解。

(Ⅱ)只需要证明$\mathit{\boldsymbol{\tilde C}}$是非负矩阵即可。注意A和B是M-矩阵，C≥0，由矩阵分裂方程(8)、(9)和(10)可知：A₁₂, A₂₁, B₁₂, B₂₁是非正矩阵，C_ij都是非负矩阵，即δA≤0, δB≤0, δC≥0。由Ⅰ知X₀≥0，从而

$ \widetilde{\boldsymbol{C}}=\boldsymbol{\delta} \boldsymbol{C}-\boldsymbol{\delta} \boldsymbol{A} \boldsymbol{X}_{0}-\boldsymbol{X}_{0} \boldsymbol{\delta} \boldsymbol{B} \geqslant 0。$

故公式(5)是M-矩阵Sylvester方程，且有唯一非负解δX。

(Ⅲ)由(Ⅰ)、(Ⅱ)证明可知(Ⅲ)成立。

(Ⅳ)由(Ⅰ)~(Ⅲ)知，当k=1时成立。

k=2时，原方程的解为

$ \begin{array}{c} \boldsymbol{X}^{(2)}=\boldsymbol{X}_{0}^{(2)}+\boldsymbol{\delta} \boldsymbol{X}^{(2)}+\boldsymbol{\delta} \boldsymbol{X}^{(1)}=\\ \left[\begin{array}{llll} \boldsymbol{X}_{11}^{(2)} & & & \\ & \boldsymbol{X}_{22}^{(2)} & & \\ & & \boldsymbol{X}_{33}^{(2)} & \\ & & & \boldsymbol{X}_{44}^{(2)} \end{array}\right]+\\ \left[\begin{array}{llll} & \boldsymbol{\delta} \boldsymbol{X}_{12}^{(2)} & & \\ \boldsymbol{\delta} \boldsymbol{X}_{21}^{(2)} & & & \\ & & & \boldsymbol{\delta} \boldsymbol{X}_{34}^{(2)}\\ &&\boldsymbol{\delta} \boldsymbol{X}_{43}^{(2)}& \end{array}\right]+\\ \left[\begin{array}{ll} & \boldsymbol{\delta} \boldsymbol{X}_{12}^{(1)} \\ \boldsymbol{\delta} \boldsymbol{X}_{21}^{(1)} & \end{array}\right]。\end{array} $

设第k-1次分裂时，原方程的解也满足

$ \boldsymbol{X}^{(k-1)}=\boldsymbol{X}_{0}^{(k-1)}+\boldsymbol{\delta} \boldsymbol{X}^{(k-1)}+\cdots+\boldsymbol{\delta} \boldsymbol{X}^{(2)}+\boldsymbol{\delta} \boldsymbol{X}^{(1)}。$

第k层次分裂时，由于X₀^(k-1)=X₀^(k)+δX^(k), 所以原方程的解为

$ \begin{array}{c} \boldsymbol{X}^{(k)}=\boldsymbol{X}_{0}^{(k)}+\boldsymbol{\delta} \boldsymbol{X}^{(k)}+\cdots+\boldsymbol{\delta} \boldsymbol{X}^{(2)}+\boldsymbol{\delta} \boldsymbol{X}^{(1)}=\\ \left[\begin{array}{cccc} \boldsymbol{X}_{11}^{(k)} & & & \\ & \boldsymbol{X}_{22}^{(k)} & & \\ & & \ddots & \\ & & & \boldsymbol{X}_{2^{k}, 2^{k}}^{(k)} \end{array}\right]+\\ \left[\begin{array}{ccc} &\boldsymbol{\delta} \boldsymbol{X}_{12}^{(k)} &&&&&& \\ \boldsymbol{\delta} \boldsymbol{X}_{21}^{(k)} &&&&&&&\\ &&& \boldsymbol{\delta} \boldsymbol{X}_{34}^{(k)} &&& &\\ && \boldsymbol{X} _ { 4 3 } ^ { ( k ) } &&&&& \\ &&&&& \ddots && \\ &&&& \ddots &&& \\ &&&&&&& \boldsymbol{\delta} \boldsymbol{X}_{2 ^{k}-1,2^k}^{(k)}\\ &&&&&&\boldsymbol{\delta} \boldsymbol{X}_{2 ^k, 2 ^{k}-1}^{(k)}& \end{array}\right]+\cdots+\\ \left[\begin{array}{llll} & \boldsymbol{\delta} \boldsymbol{X}_{12}^{(2)} & & \\ \boldsymbol{\delta} \boldsymbol{X}_{21}^{(2)} & & & \\ & & & \boldsymbol{\delta} \boldsymbol{X}_{34}^{(2)}\\ &&\boldsymbol{\delta} \boldsymbol{X}_{43}^{(2)}& \end{array}\right]+ \left[\begin{array}{ll} & \boldsymbol{\delta} \boldsymbol{X}_{12}^{(1)} \\ \boldsymbol{\delta} \boldsymbol{X}_{21}^{(1)} & \end{array}\right]。\end{array} $

定理3    在定理2条件下，且A, B, C具有HOLDR结构，则递归层对应的矩阵方程都是M-矩阵Sylvester方程。

证明    由引理1中(Ⅱ)可知，A和B的对角块子矩阵都是非奇异M-矩阵。从而A和B的非对角块子矩阵都是非正矩阵，C≥0，由定理2可得结论成立。

下面给出求解MSE的分而治之方法。

算法3    DACMSE(A，B，C) % 求解AX+XB=C

1.给定A, B, C。 2.当A, B, C具有HOLDR结构时，分块。 3.$\begin{aligned} &A=\left[\begin{array}{cc} A_{11} & 0 \\ 0 & A_{22} \end{array}\right]+U_{A} V_{A}^{\mathrm{T}}, \\ &B=\left[\begin{array}{cc} B_{11} & 0 \\ 0 & B_{22} \end{array}\right]+U_{B} V_{B}^{\mathrm{T}}, C=\left[\begin{array}{cc} C_{11} & 0 \\ 0 & C_{22} \end{array}\right]+U_{C} V_{C}^{\mathrm{T}} 。\end{aligned}$ 4.$X_{i i} \leftarrow \mathrm{DAC}_{\mathrm{MSE}}\left(A_{i i}, B_{i i}, C_{i i}\right), i=1, 2 。$ 5.令$X_{0}=\left[\begin{array}{ll} X_{11} & \\ & X_{22} \end{array}\right]$。 6.令$U=\left[U_{C}, -U_{A}, -X_{0} U_{B}\right], V=\left[V_{C}, X_{0}^{\mathrm{T}} V_{A}, V_{B}\right]$。 7.利用低秩算法求解$\delta X \leftarrow \mathrm{LR}_{\mathrm{MSE}}(A, B, U, V)$。 8.返回：X0+δX。

4 数值算例

在本节中，我们通过数值算例来比较两种不同分而治之法求解MSE的效果。基于工具包hm-toolbox^[11]在MATLAB2019b运行所有算法。相对误差如下:

$ \operatorname{Res}=\frac{\|\boldsymbol{A} \boldsymbol{X}+\boldsymbol{X} \boldsymbol{B}-\boldsymbol{C}\|_{2}}{\|\boldsymbol{A} \boldsymbol{X}\|_{2}+\|\boldsymbol{X} \boldsymbol{B}\|_{2}+\|\boldsymbol{C}\|_{2}}。$

例1    系数矩阵分别为：

$ \boldsymbol{A}=\left[\begin{array}{ccccc} 3 & -1 & & & \\ -1 & 3 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 3 & -1 \\ & & & -1 & 3 \end{array}\right]_{n \times n}, $

$ \boldsymbol{B}=2 \boldsymbol{A}, \boldsymbol{C}=\left[\begin{array}{llll} 5 & & & \\ & 1 & & \\ & & & \\ & & \ddots & \\ 2 & & & 1 \end{array}\right]_{n \times n}。$

我们设置最小块阶数不超过100，分别使用拓展Krylov分而治之方法(dac-ek)和拓展Krylov ADSM分而治之方法(dac-ek ADSM)求解，从表 1可知dac-ek ADSM在时间、精度方面均得到改进，随着阶数的增大，计算时间节省。图 2是求解不同阶数MSE的计算时间。

例2    系数矩阵为随机带状M-矩阵，

$ x=\operatorname{rand}(n-1,1), $

$ \boldsymbol{A}=2 \operatorname{eye}(n)-\operatorname{diag}(x,-1)-\operatorname{diag}(x, 1), $

$ \boldsymbol{B}=2 \boldsymbol{A}^{\mathrm{T}}, $

$ \boldsymbol{C}_{0}=\operatorname{diag}(\operatorname{rand}(n, 1))+\operatorname{diag}(x,-1)+\operatorname{diag}(x, 1), $

$ \boldsymbol{U}_{\boldsymbol{C}}=\operatorname{rand}(n, 2), \boldsymbol{V}_{\boldsymbol{C}}=\operatorname{eye}(n, 2), $

$ \boldsymbol{C}=\boldsymbol{C}_{0}+\boldsymbol{U}_{\boldsymbol{C}} \boldsymbol{V}_{\boldsymbol{C}}^{\mathrm{T}} \text {。} $

图 3表明拓展Krylov ADSM分而治之法比拓展Krylov分而治之法计算速度更快。

5 结语

本文给出求解大型M-矩阵Sylvester方程化为低阶方程和低秩方程分别求解的分而治之算法，证明了递归层对应的低阶方程和低秩矩阵方程都是M-矩阵Sylvester方程。数值实验表明本文提出的算法对大型M-矩阵Sylvester方程有效。