考虑如下形式的 Stein 张量方程

其中矩阵 $A_i\in \mathbf{R}^{I_i\times I_i}$($i=1,2,\cdots,n$) 和张量 $\mathcal{F}$ 是给定的, 张量 $\mathcal{X}\in \mathbf{R}^{I_1\times I_2\times\cdots\times I_n}$ 是待求张量. 当 $\mathcal{X}\in \mathbf{R}^{I_1 \times I_2}$ 时, 方程(1.1)退化为

上式显然是 Stein 矩阵方程, 因此方程(1.1)称为 Stein 张量方程.

为了建立求解 Stein 张量方程(1.1)的有效算法, 先回顾一些求解线性方程组的子空间方法. 对于求解下列大型稀疏线性方程组$A{x}={b},$式中矩阵 $A\in \mathbf{R}^{n \times n}$ 和向量 ${b} \in \mathbf{R}^n$ 是已知的, ${x} \in \mathbf{R}^n$ 是未知的.许多学者已经提出一些有效的算法用于求解该方程组^{[1⇓⇓⇓⇓⇓⇓⇓-9]}. 众所周知,使用最广泛的方法是 Krylov 子空间方法. Krylov 子空间方法最显著的特点是只涉及矩阵与向量的乘积, 并且在所提出的算法中可以充分利用 $A$ 的特殊结构. 当 $A$ 是对称正定矩阵时, 共轭梯度法 (CG) 是最行之有效的 Krylov 子空间方法之一. 当 $A$ 是对称半正定矩阵时, 极小残差法 (MINRES) 和共轭残差法 (CR) 也有非常好的效果. 对于 $A$ 不是对称矩阵的情形, 可由 CG 法、MINRES 法和 CR 法推广得到 BiCG 法、GMRES 法和 GCR 法. 由 BiCG 算法短循环的特点, 其有效的变形, 即双共轭梯度稳定化方法 (BCGSTAB) 被提出. 此外, BiCG 算法的一些推广也逐渐被提出^[10⇓-12].

为求解Stein 矩阵方程, 许多理论分析及算法已经被提出^{[13⇓⇓⇓⇓-18]}. Zhou 等人在文献 [13] 讨论了Smith 迭代算法的计算复杂性并给出其变形, 并且由他们提出的新的 Smith 加速迭代算法在求解 Stein 矩阵方程时有更好的效果. Fuhrmann^[15]运用张量积理论的函数模型来分析Stein 矩阵方程. Goureia^[16]运用广义逆和非交换环的理论讨论矩阵方程 $AXB-X=C$ 的连续性. Chiang^[18]提出一类特殊的 Stein 矩阵方程, 并通过收缩子空间来求解. 此外, 他们还给出了解的存在性的充分必要条件.

Xu 和 Wang 等^[19⇓-21]推广 RGI、BiCG 和 BiCR 算法用于求解 Stein 张量方程, 并研究了所提出的迭代算法的收敛性. 此外, 他们给出与式(1.1)等价的线性方程组并讨论在一定条件下方程(1.1)的解. 在他们的数值实验中, 将所提出的算法与张量格式的 CGNR 算法和 CGNE 算法进行了对比, 显然张量格式的 BiCG 和 BiCR 算法更具优势. 受其启发, 本文尝试构造 BCGSTAB 算法的张量格式, 并通过数值实验验证了算法的可行性和有效性.

在本节, 将提出求解 Stein 张量方程的张量格式的 BCGSTAB 算法. 为分析 Stein 张量方程解的存在性, 先给出一些基本引理.

引理 2.1 设 $A\in \mathbf{R}^{n\times n}$, $B\in \mathbf{R}^{m\times m}$, “$\otimes$” 为 Kronecker 积, 则

(1) $\|A\otimes B\|_2=\|A\|_2\cdot\|B\|_2$;

(2) $\rho(A\otimes B)=\rho(A)\cdot\rho(B)$.

证 (1) 首先证明第一个等式. 事实上,

(3) 设 $\lambda_1, \lambda_2, \cdots, \lambda_n$ 是矩阵 $A$ 的特征值, $\mu_1, \mu_2, \cdots, \mu_m$ 是 $B$ 的特征值. 不失一般性, 设 $|\lambda_1|\geq|\lambda_2|\geq\cdots\geq|\lambda_n|$, $|\mu_1|\geq|\mu_2|\geq\cdots\geq|\mu_m|$. 则 $A\otimes B$ 的特征值是 $\lambda_l\mu_r$, 其中 $l=1,2,\cdots,n, r=1,2,\cdots,m$.由于

故

证毕.

引理 2.2^[19] 设 $A\in \mathbf{R}^{n\times n}$, 则

(1) 矩阵序列 $\sum\limits_{k=0}^\infty A^k$ 收敛当且仅当 $\rho(A)<1$;

(2) 若矩阵序列 $\sum\limits_{k=0}^\infty A^k$ 收敛, 则 $\sum\limits_{k=0}^\infty A^k=(I-A)^{-1}$.

张量 $\mathcal{X}\in \mathbf{R}^{I_1\times I_2\times\cdots\times I_n}$ 的 $k$-模式矩阵化可记为 $\mathcal{X}(k)$, 并且 $k$-模式纤维用来表示所得到的矩阵的列. 算子 "$vec$" 表示将矩阵或张量的列拉直, 从而将其转换成一个向量. 对于矩阵 $A\in \mathbf{R}^{m\times n}$ 和张量 $\mathcal{X}\in \mathbf{R}^{I_1\times I_2\times\cdots\times I_n}$, 经过列拉直后得到的向量 ${vec}(A)$ 和 ${vec}(\mathcal{X})$ 分别定义为

式中 $a_k$ 是矩阵 $A$ 的第 $k\,(k=1,2,\cdots,n)$ 列且 $\mathcal{X}(1)$ 表示张量 $\mathcal{X}$ 的 1-模式矩阵化. 由于

若将算"${vec}$"作用于等式(1.1)的两边,即

(2.1) 式中 $I$ 是阶数为 $I_1I_2\cdots I_n$ 的单位矩阵. 令

易知 $I-E$ 的谱集为

如果方程组(2.1)有唯一解, 则它的系数矩阵 $I-E$ 满秩, 这表明系数矩阵的特征值非零, 可以立即得到下面的引理.

引理 2.3 Stein张量方程(1.1)有唯一解当且仅当对任意 $\lambda_{i_k}\in\lambda(A_k)$, $i_k\in{1,2,\cdots,n_k}$, $k=1,2,\cdots,n$, 有

基于矩阵谱条件数的定义, 即对于矩阵 $A$, 其谱条件数定义为

根据此定义，可以得到系数矩阵 $I-E$ 的谱条件数的上界和下界.一方面,

由于谱半径小于或等于任何范数, 因此可知当 $\|E\|_2<1$ 时, $\rho(E)<1$. 此外, 由 $\|E\|_2<1$ 可知 $I-E$ 是非奇异的且有 $\|(I-E)^{-1}\|_2\leqslant\frac{1}{1-\|E\|_2}$. 另一方面,

显然, $\mathrm{cond}_2(I-E)$ 的上界为

由 (2.1) 式以及引理2.2, 可知

由此, 得到下面引理.

引理 2.4^[19] Stein 张量方程 (1.1) 有唯一解 $\mathcal{X}$, 且 $\mathcal{X}$ 有序列表示

当且仅当 $\rho(A_1)\rho(A_2)\cdots\rho(A_n)<1$.

根据上述引理, 也可用部分和序列来估计Stein张量方程(1.1)的解, 即

或者表示为下面的迭代格式

若使用 (2.2) 式的部分和的形式来表示(1.1)的解, 则有如下误差估计.

引理 2.5 令 $\tau=\|E\|_2=\prod\limits_{i=1}^n\|A_i\|_2$. 若 $\tau<1$ 且 $\mathcal{X}^*$ 是 Stein 张量方程 (1.1) 的唯一解, 则不等式

成立, 式中 $\mathcal{X}_m=\sum\limits_{k=0}^m\mathcal{F}\times_1A_1^k\times_2A_2^k\times\cdots\times_nA_n^k$.

证 由 $\mathcal{X}^*$ 和 $\mathcal{X}_l$ 的序列表示, 可得

证毕.

引理 2.6 令 $\tau=\prod_{i=1}^n\|A_i\|_2$. 若 $\tau<1$ 且 $\mathcal{X}^*$ 是Stein张量方程的唯一解, 则

(1) $\|\mathcal{X}^*-\mathcal{X}_k\|_{\rm F}\leqslant\tau^k\|\mathcal{X}^*-\mathcal{X}_0\|_{\rm F}$;

(2) $\|\mathcal{X}^*-\mathcal{X}_k\|_{\rm F}\leqslant\frac{\tau}{1-\tau}\|\mathcal{X}^*-\mathcal{X}_{k-1}\|_{\rm F}$;

(3) $\|\mathcal{X}^*-\mathcal{X}_k\|_{\rm F}\leqslant\frac{\tau^k}{1-\tau}\|\mathcal{X}_1-\mathcal{X}_0\|_{\rm F}$.

证 (1) 注意到

(2) 由 ${vec}(\mathcal{X}_k)=E {vec}(\mathcal{X}_{k-1})+ {vec}(\mathcal{F})$ 及 ${vec}(\mathcal{X}^*)=E {vec}(\mathcal{X}^*)+ {vec}(\mathcal{F})$, 可得

(3) 因为

所以,

证毕.

在Krylov子空间方法中, BiCG 算法具有双正交性, 这种性质同样可以推广至张量情形中.因此文献 [19] 中提出了张量格式的 BiCG 算法, 具体如下

算法 1 (张量格式的 BiCG 算法)

步骤 1 输入初始张量 $\mathcal{X}_0$, $\varepsilon>0$, 计算 $\mathcal{R}_0=\mathcal{F}-\mathcal{X}_0+\mathcal{X}_0\times_1A_1\times_2A_2\times\cdots \times_nA_n$;

步骤 2 选取 $\widetilde{\mathcal{R}}_0$ (可选取 $\widetilde{\mathcal{R}}_0=\mathcal{R}_0$). 令 $\mathcal{P}_0=\mathcal{R}_0,\widetilde{\mathcal{P}}_0=\widetilde{\mathcal{R}}_0$.置 $k=:0$;

步骤 3 计算

步骤 4 置$k=:k+1$. 若 Err $\leqslant\varepsilon$ 停算, 输出迭代次数 $k$ 和误差Err. 否则转步骤 3.

定义两个算子

根据张量和矩阵模式积的定义, 有下面的等量关系

下面的定理是文献 [19] 中张量格式的 BiCG 算法的双正交性质, 其具体证明过程如下

定理 2.1 设张量序列 ${\mathcal{R}_i}$, ${\widetilde{\mathcal{R}}_i}$, ${\mathcal{P}_i}$ 和 ${\widetilde{\mathcal{P}}_i}(i=0,1,\cdots,k)$ 由算法 1 生成, 则成立

(2.5) 式中 $\mathcal{L}$ 由 (2.4) 式所定义.

证 使用归纳法证明当 $i\neq j$ 时 $\langle\mathcal{R}_i,\widetilde{\mathcal{R}}_j\rangle=0$. 首先, 证明 $j<i$ 的情形. 当 $i=1$ 时, 若 $j=0$, 那么显然成立

假设当 $j<i\leqslant k$ 时, $\langle\mathcal{R}_i,\widetilde{\mathcal{R}}_j\rangle=0$. 则当 $i=k+1$ 时, 若 $j<k$, 那么

若 $j=k$, 则

因此, 当 $i=k+1$ 时结论成立. 用同样的方法, 也可证明当 $j>i$ 时的情形.

下面, 需证明 $\langle\mathcal{L}(\mathcal{P}_i),\widetilde{\mathcal{P}}_j\rangle=0 (i\neq j)$. 首先, 证明 $j<i$ 的情形.

当 $i=1$, $j=0$ 时,

假设当 $j<i\leqslant k$ 时, $\langle\mathcal{L}(\mathcal{P}_i),\widetilde{\mathcal{P}}_j\rangle=0$. 则当 $i=k+1$ 时, 若 $j<k$, 那么

若 $j=k$, 则

因此, 当 $i=k+1$ 时, 结论成立.用同样的方法, 可证明 $j>i$ 的情形.

推论 2.1 设张量序列 ${\mathcal{R}_i}$, ${\widetilde{\mathcal{R}}_i}$, ${\mathcal{P}_i}$ 和 ${\widetilde{\mathcal{P}}_i}(i=0,1,\cdots,k)$ 由算法1 生成, 则成立

式中 $\mathcal{L}$ 由(2.4) 式所定义.

定理 2.2^[19] 假设Stein张量方程 (1.1) 是相容的, 在不考虑舍入误差的情况下, 对任意的初始张量 $\mathcal{X}_0\in \mathbf{R}^{I_1\times I_2\times\cdots\times I_n}$, 由算法1 生成的序列 $\{\mathcal{X}_k\}$ 至多 N 步迭代收敛到 (1.1)式的精确解.

受张量格式 BiCG 算法的启发, 给出求解 Stein 张量方程的张量格式的 BCGSTAB 迭代算法如下.

算法 2 (张量格式的 BCGSTAB 算法)

步骤 1 输入初始张量 $\mathcal{X}_0$, $0\leqslant\varepsilon\ll1$, 计算 $\mathcal{R}_0=\mathcal{F}-\mathcal{X}_0+\mathcal{X}_0\times_1A_1\times_2A_2\times\cdots \times_n A_n$;

步骤 2 选取 $\widetilde{\mathcal{R}}_0$ (可选取 $\widetilde{\mathcal{R}}_0=\mathcal{R}_0$);令 $\mathcal{P}_0=\mathcal{R}_0$, 计算 $\rho_0=\langle\widetilde{\mathcal{R}}_0,\mathcal{R}_0\rangle$; 置 $k:=0$;

步骤 3 依次计算

步骤 4 置 $k:=k+1$. 若 ${\rm Err}\leqslant\varepsilon$, 则停止迭代, 输出近似解 $\mathcal{X}_k$ 和迭代次数 $k$ 及残差 Err. 否则, 转步骤 3.

同样, 类似于文献 [19] 的证明, 容易得到算法2 的收敛性定理.

定理 2.3 假设 Stein 张量方程 (1.1) 是相容的, 在不考虑舍入误差的情况下, 对任意的初始张量 $\mathcal{X}_0\in \mathbf{R}^{I_1\times I_2\times\cdots\times I_n}$, 由算法2 生成的序列 $\{\mathcal{X}_k\}$ 至多 N 步迭代收敛到 (1.1) 式的精确解.

本节, 将给出几个数值例子来说明我们提出算法的有效性, 并且将提出的算法与张量格式的 CGNR (关于法方程残差的共轭梯度) 算法和 CGNE (关于法方程误差的共轭梯度) 算法进行比较. 所有实验都是在如下环境中实现的: Windows 10 系统, 处理器Intel(R) Core(TM) i5-6500 CPU 3.20GHz, 8GB RAM. 利用 Matlab R2014b 软件环境以及 Bader 和 Kolda 开发的 MATLAB 张量工具箱执行实验. 在实验过程中, 记录了所需的迭代次数("IT")、 所消耗的 CPU 时间(以秒为单位, 记为 "CPU"),在数值实验中, 初始张量取为 $\mathcal{X}_0=\mathcal{O}$, 且终止准则为

下面两个算法分别是张量格式的CGNR算法和CGNE算法.

算法 3 (张量格式的 CGNR 算法)

步骤 1 输入初始张量 $\mathcal{X}_0$, $0\leqslant \varepsilon\ll 1$, 计算 $\mathcal{R}_0=\mathcal{F}-\mathcal{X}_0+\mathcal{X}_0\times_1A_1\times_2A_2\times\cdots \times_nA_n$;

步骤 2 计算 $\mathcal{Z}_0=\mathcal{R}_0-\mathcal{R}_0\times_1A_1^{\rm T}\times_2A_2^{\rm T}\times\cdots \times_nA_n^{\rm T}$; 令 $\mathcal{P}_0=\mathcal{Z}_0$; 置 $k:=0$;

步骤 3 依次计算

步骤 4 置 $k:=k+1$. 若 ${\rm Err}\leqslant\varepsilon$, 则停止迭代, 输出近似解 $\mathcal{X}_k$ 和迭代次数 $k$ 及残差 Err. 否则, 转第 3 步.

算法 4 (张量格式的 CGNE 算法)

步骤 1 输入初始张量 $\mathcal{X}_0$, $0\leqslant \varepsilon\ll 1$, 计算 $\mathcal{R}_0=\mathcal{F}-\mathcal{X}_0+\mathcal{X}_0\times_1A_1\times_2A_2\times\cdots \times_nA_n$;

步骤 2 计算 $\mathcal{Z}_0=\mathcal{R}_0-\mathcal{R}_0\times_1A_1^{\rm T}\times_2A_2^{\rm T}\times\cdots \times_nA_n^{\rm T}$; 令 $\mathcal{P}_0=\mathcal{Z}_0$; 置 $k:=0$;

步骤 3 依次计算

步骤 4 置 $k:=k+1$. 若 ${\rm Err}\leqslant\varepsilon$, 则停止迭代, 输出近似解 $\mathcal{X}_k$ 和迭代次数 $k$ 及残差 Err. 否则, 转步骤 3.

例 3.1 求解下面Stein张量方程的解:$\mathcal{X}-\mathcal{X}\times_1A_1\times_2A_2\times_3A_3=\mathcal{F},$ 其中 $A_i (i=1,2,3)$ 由 MATLAB 随机生成.

下面列出随机生成的 $A_1\in \mathbf{R}^{6\times 6}, A_2\in \mathbf{R}^{5\times 5}, A_3\in \mathbf{R}^{4\times 4}$, 并且构造张量 $\mathcal{F}$, 使得元素都为1 的张量 $\mathcal{X}^*$ 为精确解.实验结果可见图1, 显然从图上可知张量格式的 BCGSTAB 算法有更佳的收敛效果, 它的迭代次数远远少于张量格式的CGNR 和CGNE 算法.

$\begin{equation*} A_2=\left(\begin{array}{lllll} 0.3267 & 0.5928 & 0.2002 & 0.1851 & 0.5579 \\ 0.2798 & 0.0685 & 0.8435 & 0.0817 & 0.6388 \\ 0.9318 & 0.2052 & 0.4237 & 0.4641 & 0.0342 \\ 0.3997 & 0.7236 & 0.5448 & 0.0306 & 0.7099 \\ 0.3794 & 0.5751 & 0.5280 & 0.4350 & 0.1693 \end{array}\right), \end{equation*}$$\begin{equation*} A_3=\left(\begin{array}{llll} 0.5934 & 0.8547 & 0.5554 & 0.6302 \\ 0.6081 & 0.3832 & 0.2954 & 0.6644 \\ 0.7724 & 0.3996 & 0.3661 & 0.9921 \\ 0.0563 & 0.3254 & 0.3490 & 0.9444 \end{array}\right), \end{equation*}$

例 3.2 求解下面Stein张量方程的解

其中系数矩阵 $A_i (i=1,2,3)$ 为稀疏矩阵, 令

初始张量 $\mathcal{X}_0=\mathcal{O}$, 构造的右端项 $\mathcal{F}$ 使得 $\mathcal{X}^*=\mathrm{ones}(n,n,n)$ 是(3.1)式的精确解.

分别令 $n=20,50,100$, 计算结果见图2、图3、图4 及表1. 尽管当维数增加到 100 时, 张量格式的 CGNR 和 CGNE 算法迭代步数大大增加, 但张量格式的 BCGSTAB 算法仍保持着较好的收敛效果.

本文针对下列 Stein 张量方程方程$\mathcal{X}-\mathcal{X}\times_1A_1\times_2A_2\times\cdots\times_nA_n=\mathcal{F},$设计了基于张量形式求解该方程的双共轭梯度稳定化 (BCGSTAB) 方法. 从两个数值算例的实验结果来看, BCGSTAB 方法较 CGNR 和 CGNE 有较快的收敛速度, 不仅表现在迭代次数大大较少, CPU 时间也缩短了许多.