Novel accelerated methods of tensor splitting iteration for solving multi-systems

1.   Introduction

In this paper, we will discuss the following multi-linear system

where $\mathcal{A} = (a_{i_1\cdots i_m})\in \mathbb{R}^{[m, n]}$ is an order $m$ dimension $n$ tensor, $b \in \mathbb{R}^{n}$ is a dimension $n$ vector.

We know an essential problem in pure and applied mathematics is solving various classes of equations. The rapid and efficient calculation approaches of multi-linear systems ^[1,2,3,4,5] are becoming more and more significant in the field of science and engineering due to their wide applications, especially for the data analysis need of big data era (see ^[6,7,8,9,10]). Normally, it is hard to get the exact solution by means of direct methods even for smaller-scale general linear systems, which greatly promotes the substantial developments of presenting various kinds of iterative strategies. Many research works have been investigated in some literatures on fast solvers for the multi-linear systems (1.1). Ding and Wei ^[11,12] proposed some classical iterative, such as Jocobi, Gauss-Seidel methods, and Newton methods through translating (1.1) into the optimization problem. In general, the computational cost for the Newton method is expensive. Then Han ^[13] investigated an homotopy method by the Euler-Newton prediction-correction technique to solve multi-linear systems with nonsymmetric M-tensors, which is shown a better result than Newton method in the sense of convergence performance. Tensor splitting method and its convergence results have been studied by Liu and Li et al. ^[14]. Further some comparison results for splitting iteration for solving multi-linear systems were investigated widely in ^[15,16], however, we find that some acceleration versions can be introduced further and may probably improve their work. Motivated by ^[15], we propose an tensor alternating splitting iteration scheme for solving multi-linear systems. In practical application, the accelerated overrelaxation method (AOR), symmetric accelerated overrelaxation method (SAOR) and their preconditioned versions are generalized for solving multi-linear systems (1.1).

The remainder of this paper is organized as follows. In Section 2, some basic and useful notations are described simply. In Section 3, we will propose a tensor alternating splitting iteration scheme for solving multi-linear systems. Meanwhile, a preconditioner is introduced to accelerate the novel method. In Section 4, the classical approaches, AOR and SAOR, will be generalized to solve multilinear systems, then some numerical tests are provided to illustrate the superiority of the presented iteration methods. Finally, a concluding remark is given in Section 5.

2.   Preliminaries

Let $A\in \mathbb{R}^{[2, n]}$ and $B\in \mathbb{R}^{[k, n]}$. The matrix-tensor product $\mathcal{C} = A\mathcal{B}\in \mathbb{R}^{[k, n]}$ is defined by

The above formula can be regarded as

where $\mathcal{C}_{(1)}$ and $\mathcal{B}_{(1)}$ are the matrices generated from $\mathcal{C}$ and $\mathcal{B}$ flattened along first index. For more details, see ^[17,18].

Definition 2.1. (^[19]). Let $\mathcal{A} = (a_{i_1i_2\cdots i_m})\in \mathbb{R}^{[m, n]}$. Then the majorization matrix $M(\mathcal{A})$ of $\mathcal{A}$ is the $n\times n$ matrix with the entries

Definition 2.2. (^[15]). Let $\mathcal{A} = (a_{i_1i_2\cdots i_m}) \in \mathbb{R}^{[m, n]}$. If $M(\mathcal{A})$ is a nonsingular matrix and $\mathcal{A} = M(\mathcal{A})\mathcal{I}_m$, we call $M(\mathcal{A})^{-1}$ the order-$2$ left-inverse of tensor $\mathcal{A}$, and $\mathcal{A}$ is called left-nonsingular, where $\mathcal{I}_m$ is identity tensor with all diagonal elements be 1.

Definition 2.3. (^[15]). Let $\mathcal{A}, \mathcal{E}, \mathcal{F } \in \mathbb{R}^{[m, n]}$. $\mathcal{A} = \mathcal{E}-\mathcal{F}$ is named as a splitting of tensor $\mathcal{A}$ if $\mathcal{E}$ is left-nonsingular. A regular splitting of $\mathcal{A}$ if $\mathcal{E}$ is left-nonsingular with $M(\mathcal{E})^{-1}\geq 0$ and $\mathcal{F}\geq 0$ $($here $\leq$ or $\geq$ denotes elementwise$)$. A weak regular splitting of $\mathcal{A}$ if $\mathcal{E}$ is left-nonsingular with $M(\mathcal{E})^{-1}\mathcal{F}\geq 0$. A convergence splitting if spectral radius of $M(\mathcal{E})^{-1}\mathcal{F}$ is less than 1, i.e., $\rho(M(\mathcal{E})^{-1}\mathcal{F}) < 1$.

Definition 2.4. (^[20]). Let $\mathcal{A} = (a_{i_1i_2\cdots i_m})\in \mathbb{R}^{[m, n]}$. A pair $(\lambda, \; x)\in \mathbb{C} \times (\mathbb{C}^n \backslash\{0\})$ is called an eigenvalue-eigenvector of tensor $\mathcal{A}$ if they satisfy the systems

where $x^{[m-1]} = (x_{1}^{m-1}, x_{2}^{m-1}, \cdots, x_{n}^{m-1})^T$. The $(\lambda, x)$ is named as $H$-eigenpair if both $\lambda$ and vector $x$ are real.

Definition 2.5. Let $\rho(A) = \max\{|\lambda| |\lambda \in \sigma(A)\}$ be the spectral radius of $\mathcal{A}$, where $\sigma(\mathcal{A})$ is the set of all eigenvalues of $\mathcal{A}$.

Definition 2.6. (^[21]). Let $\mathcal{A} = (a_{i_1i_2\cdots i_m})\in \mathbb{R}^{[m, n]}$. $\mathcal{A}$ is called a $Z$-tensor if its off-diagonal entries are non-positive. $\mathcal{A}$ is called an $M$-tensor if there exists a nonnegative tensor $\mathcal{B}$ and a positive real number $\eta\geq \rho(B)$ such that

If $\eta > \rho(\mathcal{B})$, then $\mathcal{A}$ is called a strong $M$-tensor.

3.   Tensor alternating splitting iteration

Consider two tensor splittings $\mathcal{A} = \mathcal{E}_1-\mathcal{F}_1 = \mathcal{E}_2-\mathcal{F}_2$. First of all, we describe briefly alternating direction iterative method for solving multi-linear systems $\mathcal{A}x^{m-1} = b$.

By $\mathcal{A} = \mathcal{E}_1-\mathcal{F}_1$, clearly, the above multi-linear systems can be written as

i.e.,

here use the property of order $2$ left-nonsingular of tensor $\mathcal{E}_{1}$. $\mathcal{I}_m$ is an identify tensor with appropriate order. The result leads to Algorithm 3.1.

Algorithm 3.1. (Preconditioned tensor splitting iterative method (PTSI) ^[15])

$\bf{Step\; 1}$ Input a vector $b$, and a preconditioner $P$, a strong $M$-tensor $\mathcal{A}$ with (weak) regular splitting $\mathcal{A}_P: = P\mathcal{A} = \mathcal{E}_{P_1}-\mathcal{F}_{P_1}$. Given a precision $\varepsilon > 0$ and initial vector $x_0$. Set $k: = 1$.

$\bf{Step\; 2}$ If $\|\mathcal{A}_P x_k^{m-1}-b\|_2 < \varepsilon$ stop; otherwise, go to Step $3$.

$\bf{Step\; 3}$

$\bf{Step\; 4}$ Set $k: = k+1$, return to Step $2$.

Based on this, we introduce two-step tensor alternating splitting iteration method, then it gets the iterative scheme

We further consider preconditioned multi-linear systems $P\mathcal{A}x^{m-1} = Pb$. It follows from the splitting $\mathcal{A}_P: = P\mathcal{A} = \mathcal{E}_{P_1}-\mathcal{F}_{P_1} = \mathcal{E}_{P_2}-\mathcal{F}_{P_2}$ that

Set $\mathcal{G}: = M(\mathcal{E}_{P_2})^{-1}\mathcal{F}_{P_2}$. By $\mathcal{I}_m x^{m-1} = x^{[m-1]}$ where $\mathcal{I}_m$ be an identify tensor with appropriate order, we have

Hence,

The above analysis can be described concretely in Algorithm 3.2. Now, let

Next we will show the spectral radius of preconditioned iterative tensor $\rho\big(\mathcal{T}(\mathcal{E}_{P_1}, \mathcal{E}_{P_2})\big) < 1,$ namely, the proof of convergence of Algorithm 3.2.

Algorithm 3.2. (Preconditioned tensor alternating splitting iterative method (PTASI))

$\bf{Step\; 1}$ Input a vector $b$, a preconditioner $P$, a strong $M$-tensor $\mathcal{A}$ with (weak) regular splitting $\mathcal{A}_P: = P\mathcal{A} = \mathcal{E}_{P_1}-\mathcal{F}_{P_1} = \mathcal{E}_{P_2}-\mathcal{F}_{P_2}$. Given a precision $\varepsilon > 0$ and initial vector $x_0$. Set $k: = 1$.

$\bf{Step\; 2}$ If $\|\mathcal{A}_P x_k^{m-1}-b\|_2 < \varepsilon$ stop; otherwise, go to Step $3$.

$\bf{Step\; 3}$

$\bf{Step\; 4}$ Set $k: = k+1$, return to Step $2$.

Lemma 3.1. ^[11]. Let $\mathcal{A}_P = (a_{i_1i_2\cdots i_m})\in \mathbb{R}^{[m, n]}$, and $\mathcal{A}_P = \mathcal{E}_{P_1}-\mathcal{F}_{P_1} = \mathcal{E}_{P_2}-\mathcal{F}_{P_2}$ be a weak regular splitting and a regular splitting, respectively. If $\mathcal{F}_{P_2} < \mathcal{F}_{P_1}$, $\mathcal{F}_{P_2}\neq 0$, and $\rho((\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}) < 1$, then there exists a positive Perron vector $x\in \mathbb{R}^{n}$ such that

where $\rho_k = \rho(M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}),$ $k$ is a positive integer and $\mathcal{S}\in \mathbb{R}^{[m, n]}$ is a positive tensor.

Proof. Let $\mathcal{S}$ be in $\mathbb{R}^{[m, n]}$ whose entries are all equal to $1$. There exists a positive integer $N$ such that $\rho(M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1})\leq \rho(M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}) < 1$ as $k > N$. It is clear to check that $M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}$ is positive and hence is irreducible. Using the strong Perron-Frobenius theorem (see ^[22,23]), $M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}$ has a positive Perron vector $x$ such that

for $k > N$, where $\rho_k = \rho(M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}).$

Hence, it give rises to

By $\mathcal{A}_P = \mathcal{E}_{P_1}-\mathcal{F}_{P_1} = \mathcal{E}_{P_2}-\mathcal{F}_{P_2}$, one gets $M(\mathcal{A}_P) = M(\mathcal{E}_{P_1})-M(\mathcal{F}_{P_1}) = M(\mathcal{E}_{P_2})-M(\mathcal{F}_{P_2}).$

So it generates

Further, it follows from (3.11) that

here, one should notice that the condition $\mathcal{F}_{P_1}\geq\mathcal{F}_{P_2}$ and Definition 2.1, so $M(\mathcal{F}_{P_1})\geq M(\mathcal{F}_{P_2})$, $\mathcal{F}_{P_1}-M(\mathcal{F}_{P_1})\mathcal{I}_m\geq \mathcal{F}_{P_2}-M(\mathcal{F}_{P_2})\mathcal{I}_m$ and $\rho_k < 1$.

By some simple computations, we have

Observe that $M(\mathcal{E}_{P_2})^{-1}\geq0$ and $\mathcal{F}_{P_2}\geq0$ due to the regular splitting of $\mathcal{A}_P = \mathcal{E}_{P_2}-\mathcal{F}_{P_2}$. From (3.14), it gets

Theorem 3.1. Let $\mathcal{A}_P = (a_{i_1i_2\cdots i_m})\in \mathbb{R}^{[m, n]}$, and $\mathcal{A}_P = \mathcal{E}_{P_1}-\mathcal{F}_{P_1} = \mathcal{E}_{P_2}-\mathcal{F}_{P_2}$ be a weak regular splitting and a regular splitting, respectively. If $\mathcal{F}_{P_2} < \mathcal{F}_{P_1}$, $\mathcal{F}_{P_2}\neq 0$ and $\rho\big((\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}\big) < 1$. $\mathcal{G}: = M(\mathcal{E}_{P_2})^{-1}\mathcal{F}_{P_2}$ is order $2$ left-nonsingular, i.e., $\mathcal{G} = M(\mathcal{G})\mathcal{I}_m$, then $\rho\big(\mathcal{T}(\mathcal{E}_{P_1}, \; \mathcal{E}_{P_2})\big) < 1,$ where $\mathcal{T}(\mathcal{E}_{P_1}, \mathcal{E}_{P_2})$ is defined by $(3.7)$.

Proof. First of all, similar to the previous discussion, using the strong Perron-Frobenius theorem, $\exists$ $N > 0$, $M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}$ has a positive Perron vector $x$ such that

for $k > N$, where $\rho_k = \rho(M(\mathcal{E}_{P_1})^{-1}\mathcal{F}_{P_1}+\frac{1}{k}\mathcal{S}).$ That it to say

where the inequality comes from the Lemma 3.1. When $k\rightarrow \infty$, it gets the result. This completes the proof.

4.   Numerical experiments

In this section, some numerical examples are discussed to validate the performance of effectiveness of the proposed preconditioned tensor AOR ('PTAOR') and preconditioned tensor alternating splitting AOR ('PTAAOR') based two-step splitting method for solving the multi-linear systems (see Algorithms 3.1, 3.2). We compare the convergence of preconditioned tensor Jacobi method (denoted as 'PTJb'), preconditioned tensor Gauss-Seidel method (denoted as 'PTGS') and preconditioned tensor SOR method (denoted as 'PTSOR') and unpreconditioned versions by the iteration step (denoted as 'IT'), elapsed CPU time in seconds (denoted as 'CPU'), and residual error (denoted as 'RES').

Now, consider the tensor preconditioned splitting of $(1.1)$.

The layout of splitting description is shown in Table 1, where $\mathcal{D} = D\mathcal{I}_m$, $\mathcal{L} = L\mathcal{I}_m$, $\mathcal{U} = U\mathcal{I}_m$, and $D$, $-L$, $-U$ are the diagonal part, strictly lower and strictly upper triangle part of $M(P\mathcal{A})$. Clearly, the splittings described in Table 1 satisfy the condition of weak regular splitting, i.e., $M(\mathcal{E})^{-1}\mathcal{F}\geq 0$.

A type of preconditioner (as a variant form in ^[24,25]) $P_{\alpha} = I+S_{\alpha}$ is considered for solving $P\mathcal{A}x^{m-1} = Pb$, where

$\alpha_i = 0.01, i = 1, 2, \cdots n-1.$

All the numerical experiments have been carried out by MATLAB R2011b 7.1.3 on a PC equipped with an Intel(R) Core(TM) i7-2670QM, CPU running at 2.20 GHZ with 8 GB of RAM in Windows 7 operating system.

Example 4.1. First consider the multi-linear systems (1.1) with a strong $\mathcal{M}$-tensor $\mathcal{A}$ in difference cases ^[14,15].

Case $1$. $\mathcal{A} = 864.4895 \mathcal{I}_m-\mathcal{B}$, where $\mathcal{B}\in \mathbb{R}^{[3,5]}$ is a nonnegative tensor with $b_{i_1, i_2, i_3} = |tan(i_1+i_2+i_3)|.$

Case $2$. $\mathcal{A} = n^2 \mathcal{I}_3-\mathcal{B}$, where $\mathcal{B}\in \mathbb{R}^{[3,3]}$ is a nonnegative tensor with $b_{i_1, i_2, i_3} = |sin(i_1+i_2+i_3)|.$

Case $3$. $\mathcal{A} = s \mathcal{I}_3-\mathcal{B}$, where $\mathcal{B}$ is generated randomly by MATLAB, and $s = (1+\delta)\max_{1, 2, \cdots, n}(\mathcal{B}{\bf{e}}^2)_i,$ ${\bf{e}} = (1, 1, \cdots, 1)^T$. Let $\delta = 8.$ $b = x_0 = {\bf{e}}$.

We give three different cases for different tensors $\mathcal{A}$, $\mathcal{B}$ with various sizes. Parameter $r$ and $\omega$ are experiential selected according to particular example. In this example, the running is terminated when the current iteration satisfies $\mathrm{RES} = \|\mathcal{A}x^{m-1}-b\|_2 < 10^{-11}$ or if the number of iteration exceeds the prescribed iteration steps $k_{max} = 500$.

The numerical results have been shown in Tables 2, 3 and Figures 1, 2. From the numerical results, we can see that PTAAOR and PTAOR are efficient methods, and both of them overmatch the PTSOR, PTGS, and PTJb in all sides. PTAAOR seems to be a fascinating approach, however the PTAOR is more efficient method than the PTSOR, PTGS, PTJb methods from all aspects due to the flexible selection of parameters. It is clear PTSOR and PTGS are nearly similar efficiency when $\omega$ approaching to $1$. Meanwhile, we find that PTGS is superior to PTSOR from the view of iteration number. It further bears out the conclusions in ^[15], one of which clarifies a fact the spectral radius of PTGS is small than PTSOR's (Corollary 5.6 ^[15]). Further from the residual trend chart with the changing numbers of iteration in Figure 1, one can demonstrably find the desired performance of the proposed methods.

Example 4.2. Consider the following higher-order Markov chain model:

where $\mathcal{B} = (b_{i_1, i_2, \cdots, i_m})$ is an order $m$ tensor representing an $(m-1)th$ order Markov chain, which is called an order $m$ dimension $n$ transiting probability tensor, i.e., $b_{i_1, i_2, \cdots, i_m}\geq 0, \; \; \Sigma_{i_1 = 1}^{n}b_{i_1, i_2, \cdots, i_m} = 1,$ and $x$ is named as a stochastic vector with $x_i\geq0$ and $\Sigma_{i = 1}^{n}x_i = 1$ ^[16].

Observe that $(4.2)$ can be transformed into the following systems

Next we set tensor $\mathcal{A}: = I_m-\beta \mathcal{B}$, $b: = x^{[m-1]}-\beta x$, then use splitting iteration approaches to solve the systems $(4.3)$, where

$\mathcal{I}_m$ is an identity tensor of order $m$ dimension $n$. In this example, the running is terminated when the current iteration satisfies $\mathrm{RES} = \|\mathcal{A}x^{m-1}-b\|_2 < 10^{-6}$ or if the number of iteration exceeds the prescribed iteration steps $k_{max} = 100$.

In this example, all the numerical results are depicted in Table 4 and Figure 3. From the elapsed CPU and numbers of iteration, PTAOR performs always well than all of other approaches, although the optimum parameters are hard to determine, we can just choose them tentatively in some reality application according to experiment initial effects. Preconditioned scheme can modify all these methods to some extent, but the efficiency need to be improved further in some actual applications, which depends closely on the construction of preconditioner. Hence, the research of optimum parameters and preconditioners for PTAOR will be further proceeded in future. PTAAOR can be considered as a novel and efficient approach, however, the selection of splitting tensor $\mathcal{E}_{P_1}$ and $\mathcal{E}_{P_2}$ should be adequately studied discussed in a later work.

5.   Conclusions

In this paper, an tensor alternating splitting iteration scheme is proposed for solving multi-linear systems, and the tensor accelerated overrelaxation and tensor symmetric accelerated overrelaxation splitting iteration strategies are introduced to solve this kind of systems, which can be regarded as the generalizations of AOR and SAOR for linear systems. Meanwhile, some efficient preconditioned techniques are provided to improve the efficiency of solving multi-linear systems. The proposed approaches have been demonstrated to be superior to classical SOR, GS, Jacobi methods under normal conditions, which can be fully validated in our numerical experiments section. Further, we point out our future research work.

Acknowledgments

The project is supported by Fujian Natural Science Foundation (Grant No. 2019J01879) and The New Century Training Plan of Fujian Province University (Grant No. Minjiao2017[54]), JT180586.

Conflict of interest

All authors contributed equally and significantly in writing this article. All authors read and approved the final manuscript.