基于预处理分裂迭代算法求解Mz张量方程组

Solving the Multi-Linear System with Mz-Tensor via Preconditioned Tensor Splitting Iteration

摘要: 偶数阶张量的正定性和半定性由于在谱超图理论、自动控制、多项式理论、随机过程、磁共振成像等众多领域均有着广泛的应用，近年来得到了深入而广泛的研究。研究表明，M张量、B张量、H张量、希尔伯特张量以及随机张量在特定条件下可以是正定的。然而，仍有许多正定张量无法通过上述准则来确定。在本文中，我们提出了一类全新的正定张量，其非对角项相对于强M张量可以为正，我们将其称为强Mz张量。这类强Mz张量可从离散微分方程中产生，原因在于它是基于Z特征值而非传统的H特征值。由于Z特征值定义中涉及的非齐次性通常比H特征值更难处理相应的问题，因此Z特征值的计算结果远远小于h本征值的计算结果。定义在Z特征值基础上的Mz张量在结构上比已被广泛研究的M张量更复杂。对Mz张量的研究几乎是一片空白。到目前为止，还没有人提出求解强Mz方程的算法。作为本文的动机之一，我们创新性地利用张量迭代给出了求解系数张量为Mz张量的多线性系统的张量分裂迭代算法。其次，给出了相应的预处理分裂迭代算法以及理论分析。最后，通过数值算例说明了所提算法的有效性。这在图像处理以及物理领域等多维信息分析方面，乃至高阶非线性系统的稳定性研究等领域均具有广泛的应用前景。

Abstract: The positive definiteness and semi-definiteness of even-order tensors have received extensive and in-depth research in recent years due to their wide applications in various fields such as spectral hypergraph theory, automatic control, polynomial theory, stochastic processes, and magnetic resonance imaging. Studies have shown that M-tensors, B-tensors, H-tensors, Hilbert tensors, and random tensors can be positive definite under specific conditions. However, there are still many positive definite tensors that cannot be determined by the above criteria. In this paper, we propose a new class of positive definite tensors whose off-diagonal terms can be positive relative to strong M-tensors. We call them strong Mz-tensor. This kind of strong Mz-tensor can be generated from discrete differential equations because it is based on Z-eigenvalues rather than the traditional H-eigenvalues. Due to the non-homogeneity involved in the definition of Z-eigenvalues, which is usually more difficult to handle corresponding problems than H-eigenvalues, the computational results of Z-eigenvalues are much smaller than those of h-eigenvalues. The Mz-tensor defined on the basis of Z-eigenvalues is structurally more complex than the widely studied M-tensor. The research on Mz-tensor is almost a blank. So far, no one has proposed an algorithm for solving strong Mz equations. As one of the motivations of this paper, we innovatively use tensor iteration to give a tensor splitting iterative algorithm for solving multilinear systems with coefficient tensors being Mz-tensor. Secondly, the corresponding preconditioned splitting iterative algorithm and theoretical analysis are given. Finally, numerical examples illustrate the effectiveness of the proposed algorithm. This has broad application prospects in multi-dimensional information analysis such as image processing and the physical field, and even in the stability research of high-order nonlinear systems.

1. 引言

“张量”一词最初是由威廉·罗恩·汉密尔顿于1846年提出的，用来表示现代意义上的“模块”概念。从数学上讲，张量的形式本质上是一个多维数组，可以看作是向量和矩阵的高阶推广。张量的研究方向主要包括：张量的秩、张量的行列式、张量分解、张量本征值问题、张量的低秩近似、结构化张量、张量互补问题、张量特征值互补问题、强M-张量的确定和张量方程[1]-[5]等。

求解各种方程是纯数学和应用数学中的一个基本问题。例如[4]中的线性方程和Riccati方程。

由于张量是矩阵的高阶推广，并且由于张量的应用价值和良好性质，求解张量方程的问题引起了许多学者的关注和研究。讨论了一些特殊的张量，如非负张量、M张量、Mz-张量。

最近，Mo等人[6]提出了一个多线性系统

$A{x}^{m-1}=b,$

其中$A=\left(a_{i{i}_2\cdots i_m}\right)$ 为m阶n维Mz张量，b为n维向量，左手边定义为

$A{x}^{m-1}={\displaystyle \sum }_{i_2,\cdots ,i_m=1}^n{a}_{i{i}_2\cdots i_m}{x}_{i_2}\cdots x_{i_m},i=1,2,\cdots ,n,$

其中$x_i$ 表示x的第i个分量。

Mo等人在参考文献[7]中最初提出了一类新的张量，称为(强) Mz张量，它基于Z特征值。他们证明偶数阶(强) M张量总是(强) Mz张量，而反过来通常不是真的。随后，Mo [1]证明了具有强Mz张量的多线性系统在特定条件下总是具有非负解。

由于Z特征值定义中涉及的非齐次性通常比H特征值更难处理相应的问题，因此Z特征值的计算结果远远小于h本征值的计算结果。定义在Z特征值基础上的Mz张量在结构上比已被广泛研究的M张量更复杂。对Mz张量的研究几乎是一片空白。到目前为止，还没有人提出求解强Mz方程的算法。作为本文的动机之一，我们给出了用Mz张量求解多线性系统的张量分裂算法。

众所周知，求解多线性系统的预处理技术是非常重要的。然而，这种技术还没有被用来解决多线性系统系数张量A是Mz张量的问题。因此，在本文中，我们探讨了相应的预条件张量分裂方法。

本文的主要贡献是：

• 提出了求解强Mz方程的张量分裂算法。

• 相应的预处理张量分裂迭代。

• 对提出的算法进行了理论分析，并给出了具体的最有效预处理子的收敛性分析。

2. 预备知识

定义1 [6] 我们称$\lambda$ 为张量$A\in R^{\left[m,n\right]}$ 的Z特征值，如果存在一个非零实向量x，使得

$A{x}^{m-1}=\lambda x,x^{\text{T}}x=1$ 。

这样的x称为与$\lambda$ 相关的Z特征向量，称为Z特征对。

张量的Z谱半径定义为${\rho }_z\left(B\right)=\sup \left\{\left|\lambda \right|:\lambda \text{is}\text{z-eigenvalue}\text{of}B\right\}$ 。

利用Z谱半径，回顾了Mz张量的定义。

定义2 [7] 假设m是偶数。张量$A\in R^{\left[m,n\right]}$ 称为Mz张量，如果存在一个非负张量B和一个正实数$s\ge {\rho }_z\left(B\right)$ 使得

$A=s{I}_z-B$ ，

其中$I_z=e_{i_1\cdots i_m}\in R^{\left[m,n\right]}$ 定义为非负张量的Z单位张量，这样

$I_z{x}^{m-1}=x$

对于所有$x\in R^{\left[m,n\right]}$ ，$x^{\text{T}}x=1$ 。另外，如果$s>{\rho }_z\left(B\right)$ ，则A被称为强Mz张量。

定义3 [8] 让$A\in R^{\left[m,n\right]}$ 。则A的主化矩阵$J\left(A\right)$ 是带有条目的n级矩阵

$J{\left(A\right)}_{ij}=a_{ij\cdots j},i,j=1,\cdots ,n$ 。

定义4 [9] 让A是Mz张量，$M,N\in R^{\left[m,n\right]}$ 。如果M是左非奇异的，则称$A=M-N$ 为A的分裂；如果M是左非奇异的，且$J{\left(M\right)}^{-1}\ge 0$ 以及$J{\left(M\right)}^{-1}N\ge 0$ ，则称之为弱正则分裂；如果M是左非奇异的，且$J{\left(M\right)}^{-1}\ge 0$ 以及$N\ge 0$ ，则称之为正则分裂；如果${\rho }_z\left(J{\left(M\right)}^{-1}N\right)<1$ ，则称之为收敛分裂。

3. 张量的分裂

对于一般张量，我们将系数张量A分解为$A=M-N$ ，使得具有系数张量M的方程易于求解且N是非负的。特别地，对于一个强Mz张量和一个正向量b，我们给出了一个分裂迭代方法来解决多线性系统如下。

算法1 张量分裂法

(1) 给定一个正向量b，一个具有(弱)正则分裂的强Mz张量，最大迭代步数$k_{\max }$ ，一个机器精度$\epsilon$ 和一个正初始向量x₀。初始化k = 1。

(2) 当$k<k_{\max }$ 。

(3) $x_k={\left(J{\left(M\right)}^{-1}N{x}_{k-1}^{m-1}+J{\left(M\right)}^{-1}b\right)}^{\left[\frac{1}{m-1}\right]}$ 。

(4) 如果${\Vert A{x}^{m-1}-b\Vert }_2<ϵ$ ，断开并输出$x_k$ 。

(5) $k=k+1$ ，回到(2)。

定理1 设A是强Mz张量，则对于其每一个正则分裂，算法1线性收敛。

证明 我们首先证明算法生成的序列是单调的，再证明其有界，于是其收敛。再通过证明迭代雅可比矩阵的谱半径小于1证明其是线性收敛的。

4. 预处理张量分裂法

在本节中，我们考虑用预处理方法求解具有Mz张量的多线性系统。设P是一个非奇异矩阵。然后我们考虑下面的预处理多线性系统，它等价于$PA{x}^{m-1}=Pb$ 。矩阵P称为求解的预条件。设$PA=M_P-N_P$ 为$PA$ 的一个分裂。

算法2 预处理张量分裂法

(1) 给定一个正向量b，一个具有(弱)正则分裂的强Mz张量，最大迭代步数$k_{\max }$ ，一个机器精度$\epsilon$ 和一个正初始向量x₀。初始化k = 1。

(2) 当$k<k_{\max }$ 。

(3) $x_k={\left(J{\left(M_P\right)}^{-1}{N}_P{x}_{k-1}^{m-1}+J{\left(M_P\right)}^{-1}Pb\right)}^{\left[\frac{1}{m-1}\right]}$ 。

(4) 如果${\Vert PA{x}^{m-1}-Pb\Vert }_2<ϵ$ ，断开并输出$x_k$ 。

(5) $k=k+1$ ，回到(2)。

定理2 设A是强Mz张量，$PA=M_P-N_P$ 是一个收敛的正则分裂。

证明 我们首先证明对于具体的预处理矩阵来说$PA$ 是强Mz张量，再证明$PA=M_P-N_P$ 是一个正则分裂，于是由定理1，证明完成。

5. 数值例子

在这一部分，我们将给出一些数值例子来检验提出的算法，并比较预处理张量分裂方法和那些没有预处理的方法。所有测试将在MATLAB R2016a中完成，配置为Intel(R) Core(TM)i3-8100 CPU 3.60 GHz 和8.00 GB内存。在下面的例子中，给出了两个方面来检验提出的算法的效率：迭代步骤的数量(由IT表示)，CPU时间(以秒表示)。在数值实验中，我们取$s_{\alpha }=\alpha s_1$ 。$A{x}^{m-1}$ 是通过转换成以下矩阵向量积来计算的：

$A{x}^{m-1}=A_{\left(1\right)}\left(x\otimes \cdots \otimes x\right)$

其中$\otimes$ 是克罗内克内积，x出现$m-1$ 次。如果${\Vert A{x}^{m-1}-b\Vert }_2<ϵ$ ，所有被测试的方法都停止。

注1 如果预处理子$\alpha =0$ ，那么$P_{\alpha }=P_0=I$ 和P₀ Jacobi，P₀ GS，P₀ SOR在没有预处理的情况下还原为相应的方法。

由于强Mz张量必须是偶数阶的(见[7])，我们选择m = 4阶的强Mz张量进行数值实验。我们在0.1的区间内将α或β的值从0.1取到2。对于SOR类型方法，我们在0.1的区间内搜索0.1至2的最佳参数以研究最佳参数。例1被认为是测试提出的预处理方法。我们在表1和表2中用粗体显示了CPU的最佳结果。下面的例子可以在[1] [2]中类似地获得。

例1 由[6]中的(2.2)给出的Z单位张量$I_Z$ 是弱对称的。假设$B1\in R^{[4,2]}$ 是一个所有条目都为1的张量，设$B=I_Z+B1$ ，则b是弱对称且不可约的。可以计算${\rho }_Z\left(B\right)=5$ 。那么$A=5.2I_Z-B$ 是一个强Mz张量。

注2 在表1、表2和图1中，停机准则为$\epsilon ={10}^{-13}$ ，右侧向量$b=\left(0,1\right)$ 。数值结果见表1和表2，结果表明，预处理$P_{\alpha }$ 或$P_{\beta }$ 的预处理方法在CPU时间内比没有预处理的方法表现得更好。此外，预处理剂$P_{\alpha }$ 的效果通常优于$P_{\beta }$ 。

Table 1. Numerical results of the preconditioned tensor splitting methods for Example 1

表1. 例1的预处理张量分裂方法的数值结果

	PaJacobi			PaGS			PaSOR
Alpha	CPU(s)	IT	err	CPU(s)	IT	err	CPU(s)	IT	err
0	0.1135	235	8.95E−14	0.0094132	199	9.65E−14	0.0035793	75	7.02E−14
0.01	0.0104059	233	9.18E−14	0.0090086	198	9.11E−14	0.0032484	74	6.75E−14
0.02	0.010516	231	9.39E−14	0.0089181	197	8.62E−14	0.0031261	73	8.20E−14
0.1	0.0097803	217	8.73E−14	0.0080842	186	9.08E−14	0.0029742	69	9.65E−14
0.2	0.0087005	198	9.92E−14	0.0077507	173	8.39E−14	0.0027678	65	8.39E−14
0.3	0.0079942	180	8.92E−14	0.0073992	159	8.98E−14	0.0025636	61	7.29E−14
0.4	0.0072657	161	8.63E−14	0.0066482	145	8.87E−14	0.0028706	57	5.94E−14
0.5	0.0067549	141	8.14E−14	0.0057596	130	9.40E−14	0.0022554	52	7.94E−14
0.6	0.0054432	119	7.72E−14	0.0052237	115	7.75E−14	0.0020183	48	5.29E−14
0.7	0.0041602	91	8.39E−14	0.0044313	97	7.20E−14	0.0017928	43	5.39E−14
0.8	0.0035717	75	8.97E−14	0.0030115	68	9.66E−14	0.0016882	37	8.62E−14
0.9	0.0037026	86	6.87E−14	0.0031694	72	6.80E−14	0.0014438	29	5.48E−14
1	0.0044173	100	9.28E−14	0.0034502	78	7.89E−14	0.0009204	22	2.14E−14
1.1	0.0058616	128	7.57E−14	0.0040607	89	9.30E−14	0.0012547	25	9.43E−14
1.2	0.0108985	166	5.90E−14	0.0047465	99	9.60E−14	0.0012483	30	6.79E−14
1.3	0.0104176	238	7.93E−14	0.0050708	115	5.80E−14	0.0015859	34	6.42E−14
1.4	0.0189009	423	8.62E−14	0.0063244	136	7.58E−14	0.0016338	39	6.44E−14

Table 2. Numerical results of the preconditioned tensor splitting methods for Example 1

表2. 例1的预处理张量分裂方法的数值结果

	PbJacobi		PbGS		PbSOR
Beta	CPU(s)	IT	CPU(s)	IT	CPU(s)	IT
0	0.1135	235	0.0094132	199	0.0035793	75
0.01	0.0106548	233	0.0094925	199	0.0032394	75
0.02	0.0111039	232	0.0093884	199	0.0036244	76
0.1	0.010139	220	0.0088232	196	0.0036678	80
0.2	0.009392	206	0.0089294	193	0.0046896	85
0.3	0.0087654	191	0.0086217	190	0.0055055	89
0.4	0.0080678	176	0.0092442	187	0.0052891	93
0.5	0.0079344	161	0.0084815	184	0.0061409	102
0.6	0.0070149	145	0.0090109	181	0.0089577	111
0.7	0.0065293	129	0.0084287	178	0.0064621	117
0.8	0.0057004	110	0.008424	175	0.0070495	124
0.9	0.0045892	87	0.0081532	172	0.0083511	129
1	0.0032896	72	0.007638	169	0.0074936	133
1.1	0.0035227	81	0.0077751	166	0.0069604	138
1.2	0.0044582	92	0.0074952	163	0.0071947	141
1.3	0.0054573	106	0.0070962	160	0.006949	144
1.4	0.0058115	128	0.0069804	157	0.0080971	147

Figure 1. The relationship between errors and iteration steps in Example 1

图1. 例1中误差和迭代步骤之间的关系

下面给出了来自微分方程离散化的一个例子。

例2 设$A\in R^{\left[4,2\right]}$ 且所有条目均为0，除了

$\begin{array}{ll}A\left(i,i,i,i\right)=1,\hfill & \text{for}i=1,\cdots ,n,\hfill \\ A\left(i,i,i+1,i+1\right)=1\hfill & \text{for}i=1,\cdots ,n-1,\hfill \\ A\left(i,i,i,i+1\right)=-2\hfill & \text{for}i=1,\cdots ,n-1.\hfill \end{array}$

让$s>1$ 且$B=s{I}_Z-A$ ，则B是非负的。这个四阶张量A是一个强Mz张量；有关更多细节，请参见[1]。

注3 我们取$b=\left(0,1\right)$ ，并在表2中为示例1重新报告了运行时间(CPU(s))和迭代步骤(IT)的数值结果。在图1中，我们取$M=J\left(A\right)I$ ，并显示示例1中误差的对数与迭代步骤之间的关系。从表3和图2，图3中，我们取$b=\left(1,1\right)$ 和初始向量$x_0=\left(0,0\right)$ 。显然，即使数值规模非常大，算法1仍然有一个良好的性能。我们注意到在图2到图3中，当n变大时，拐点向右移动。这是因为我们提出的算法是局部线性收敛的。也就是说，迭代向量越接近真值，收敛速度就越快。

Table 3. The running time (CPU(s)) and iterative steps (IT) of different scales for Example 2

表3. 例2中不同规模的运行时间(CPU)和迭代步骤(IT)的关系

SIZE	CPU(s)	IT	err
2*2*2*2	0.0040000	75	6.84E−14
10*10*10*10	0.0416766	702	9.58E−12
30*30*30*30	3.6406317	3693	9.55E−12
50*50*50*50	58.4339360	7981	9.66E−10
60*60*60*60	1.44E+02	9927	9.94E−07
70*70*70*70	3.64E+02	12754	9.73E−06

Figure 2. The relationship between errors and iteration steps when n = 30 in Example 2

图2. 例2中，当n = 30时误差与迭代步长之间的关系

Figure 3. The relationship between errors and iteration steps when n = 70 in Example 2

图3. 例2中，当n = 70时误差与迭代步长之间的关系

例3 让$A=s{I}_Z-B\in R^{\left[4,2\right]}$ 是一个对称张量，其中$I_Z\in R^{\left[4,2\right]}$ 是一个z单位张量且

$I_z\left(:,:,1,1\right)=\left(\begin{array}{cc}1& 0\\ 0& 1/3\end{array}\right),I_z\left(:,:,2,1\right)=\left(\begin{array}{cc}0& 1/3\\ 1/3& 0\end{array}\right)$

$I_z\left(:,:,1,2\right)=\left(\begin{array}{cc}0& 1/3\\ 1/3& 0\end{array}\right),I_z\left(:,:,2,2\right)=\left(\begin{array}{cc}1/3& 0\\ 0& 1\end{array}\right),$

和$B\in R^{\left[4,2\right]}$ 是具有如下条目的非负张量

$B\left(:,:,1,1\right)=\left(\begin{array}{cc}6& 3.75\\ 3.75& 5\end{array}\right),B\left(:,:,2,1\right)=\left(\begin{array}{cc}3.75& 5\\ 5& 7.5\end{array}\right),$

$B\left(:,:,1,2\right)=\left(\begin{array}{cc}3.75& 5\\ 5& 7.5\end{array}\right),B\left(:,:,2,2\right)=\left(\begin{array}{cc}5& 7.5\\ 7.5& 5\end{array}\right).$

我们用Matlab计算${\rho }_Z\left(B\right)=21.75$ 。所以A是一个强Mz张量。

注4 我们在表4中显示了示例3的数值结果。这表明了我们提出的求解Mz张量的算法的有效性。值得一提的是，当s < 21.75时，算法不收敛。也就是说，在例3中若A不是Mz张量，算法1不收敛。

Table 4. The running time (CPU(s)) and IT of Example 3 with different x₀ and b

表4. 例3中不同x₀和b的运行时间(CPU(s))和迭代步骤(IT)的关系

x0	b	CPU(s)	IT	err
(0,0)	(1,0)	0.0264562	563	9.75E−12
(1,1)	(1,0)	0.0239597	501	9.92E−12
(1,1)	(1,1)	0.0259110	547	9.61E−12
(0,0)	(1,1)	0.0269045	561	9.93E−12

6. 结束语

本文给出了求解强Mz方程的张量分裂算法，包括相应的预条件分裂迭代算法。然后给出了算法的收敛性和收敛速度的证明。最后，通过数值算例说明了所提算法的有效性。在未来，我们将探索强Mz张量的一些特征，以便更好地理解(强) Mz张量。

致 谢

作者在此感谢编辑和审稿人提出的有益意见和建议，使本文的呈现有了很大的改进。