Drug delivery from engineered organisms and nanocarriers as monitored by multimodal imaging technologies

1.   Introduction

In this work, we consider the following multi-linear system:

where $\mathcal{A}\in\mathbb{R}^{[m, n]}$ and ${\mathbf y}, {\mathbf b}\in\mathbb{R}^{n}$. Here, $\mathbb{R}^{[m, n]}$ and $\mathbb{R}^{n}$ denote the sets of the $m$-order $n$-dimensional tensors and $n$-dimensional real vectors, respectively. Besides, the tensor-vector product $\mathcal{A}{\mathbf y}^{m-1}$ in (1.1) is defined by ^[1]:

The multi-linear system (1.1) arises in many scientific and engineering fields, such as the tensor complementarity problems, image processing, numerical partial differential equations and so on ^{[2,3,4,5,6,7,8]}. Moreover, it is demanded to compute the solutions of multi-linear system (1.1) in many practical problems, thus solving the multi-linear system (1.1) as one of the topics of very active research in computational mathematics. So in this paper, we focus on the numerical solution of the multi-linear system (1.1).

It is noteworthy that the coefficient tensors of the multi-linear systems have special structures and they are $\mathcal{M}$-tensors in lots of practical applications. Thus, it is inevitable to solve the multi-linear systems with $\mathcal{M}$-tensors in many practical problems. In view of this fact, designing the efficient methods for solving the multi-linear systems with $\mathcal{M}$-tensors has been a hot research topic. Many researchers have established many results and numerous methods for the multi-linear systems with $\mathcal{M}$-tensors recently. For instance, Ding and Wei ^[8] proved that the multi-linear system (1.1) has a unique positive solution when $\mathcal{A}$ is a strong $\mathcal{M}$-tensor and ${\mathbf b} > {\mathbf 0}$, and they also designed the Newton, Jacobi and Gauss-Seidel methods for (1.1). Han ^[9] proposed the homotopy method for the multi-linear systems with $\mathcal{M}$-tensors. Xie et al. ^[10] derived the tensor methods for solving symmetric $\mathcal{M}$-tensor systems. Wang et al. ^[11] proposed the neural networks based method for (1.1). He et al. ^[12] developed a Newton-type method to solve multi-linear systems with $\mathcal{M}$-tensors. Besides, tensor splitting method ^[13,14] has been presented for solving (1.1) recently. Among these methods, the tensor splitting method is one of the most popular methods, because it is efficient and easy to be implemented in the practical calculation. To solve multi-linear systems effectively, Liu et al. ^[13] proposed the tensor splitting $\mathcal{A} = \mathcal{E}-\mathcal{F}$, and then established a general tensor splitting iterative method in the following form:

where $M(\mathcal{E})^{-1}\mathcal{F}$ is called the iteration tensor of the tensor splitting method (1.2), and $M(\mathcal{E})$ stands for the majorization matrix of the tensor $\mathcal{E}$ (see Definition 3). As stated in ^[13], the spectral radius of the iteration tensor can be served as an approximate convergence rate for the tensor splitting iterative method (1.2).

It is well-known that preconditioning technique can decrease the convergence factors of the iteration methods so as to accelerate their convergence speeds. Thus, it has been widely applied to the tensor splitting methods and many efficient preconditioners have been designed for (1.1) in recent years. In ^[14], Li et al. constructed the preconditioner $P_{\alpha} = I+S_{\alpha}$ with

and $\alpha_{q}\in{\mathbb{R}}$ $(q = 1, \cdots, n-1)$ are parameters. They proposed the related preconditioned Gauss-Seidel type and successive over-relaxation (SOR) type methods for (1.1). Then, Cui et al. ^[15] put forward a preconditioner $P_{\max} = I+S_{\max}$ with

where $k_{i} = \min\{j|\max\limits_{j}|a_{ij\cdots j}|, i < j, i < n\}$. It is shown in ^[15] that the preconditioner $P_{\max}$ outperforms $P_{\alpha}$ from the perspective of both the theoretical and numerical results. Subsequently, Zhang et al. ^[16] derived a new preconditioner $\tilde{P}_{\alpha} = I+G_{\alpha}$, where

and $\alpha_{w}\in{\mathbb{R}}$ $(w = 2, \cdots, n)$ are parameters, and proposed three preconditioned Jacobi type methods. Then, by making use of the preconditioner $\tilde{P}_{\alpha}$, Liu et al. ^[17] presented a new preconditioned SOR method for the multi-linear system (1.1), and a new preconditioned Gauss-Seidel type method can be obtained as the relaxation factor $\tau = 1$. Very recently, Chen and Li ^[18] constructed a new preconditioner $P_{\gamma\beta} = I+S_{\gamma\beta}$ with

where $\gamma_{j} > 0$ and $\beta_{j}$ $(j = 2, \cdots, n)$ are constants. Besides, Cui et al. ^[19] introduced a matrix $R_{\beta}$ into the preconditioner $P_{\alpha} = I+S_{\alpha}$ ^[14] and proposed a more effective preconditioner $\hat{P}_{\beta} = I+S_{\alpha}+R_{\beta}$ with

and $\beta_{q}\in\mathbb{R}$ $(j = 2, \cdots, n)$. Xie and Miao ^[20] stated that the preconditioning effect of $P_{\alpha}$ is not observed on the last horizontal slice of the system tensor $\mathcal{A}$ and derived a new preconditioner $\hat{P}_{\gamma} = I+S_{\alpha}+R_{\gamma}$ and the corresponding preconditioned Gauss-Seidel method, where

and $\gamma_{t}\in \mathbb{R}$ ($t = 1, \cdots, n-1$) are parameters. Comparison theorems and test results in ^[20] indicate that the preconditioned Gauss-Seidel method in ^[20] is more efficient than the one in ^[14]. In the sequel, Nobakht-Kooshkghazi and Najafi-Kalyani ^[21] presented a preconditioner $\tilde{P} = I+S_{\alpha}+R_{\gamma}+C_{\alpha}$, where $(C_{\alpha})_{kn} = -\alpha_{k}a_{kn\cdots n}$ $(k = 1, \cdots, n-1)$ and the remaining elements in $C_{\alpha}$ are zero. Based on it, the authors in ^[21] established a new preconditioned Gauss-Seidel method for solving (1.1) and compared it with the one in ^[20]. Very recently, An et al. ^[22] replaced the matrix $S_{\alpha}$ in $\tilde{P} = I+S_{\alpha}+R_{\gamma}+C_{\alpha}$ ^[21] by $S_{\max}$ and proposed a new preconditioner $\hat{P} = I+S_{\max}+R_{\gamma}+C_{\alpha}$. It is proved in ^[22] that the corresponding preconditioned Gauss-Seidel method has faster convergence rate than the one in ^[21].

In addition to the preconditioners for multi-linear systems with $\mathcal{M}$-tensors, some researchers have extended the related results to the multi-linear systems with other or more general tensors. For example, Wang et al. ^[23] investigated and developed a preconditioned tensor splitting accelerated over-relaxation (AOR) for solving the multi-linear systems with nonsingular $\mathcal{H}$-tensors. Moreover, Cui et al. ^[24] provided a general preconditioner that contains several existing preconditioners. And they constructed two different preconditioned SOR-type splitting methods for the multi-linear systems with more general $\mathcal{Z}$-tensors. For more preconditioners of the multi-linear system (1.1), we can refer to ^[25,26] and the references therein.

It can be seen that the preconditioners developed in ^{[14,15,16,17,18,19,24]} have preconditioning effect on some but not all horizontal slices of the coefficient tensor $\mathcal{A}$ in (1.1) ^[27]. That is, the preconditioners in ^{[14,15,19,24]} have no preconditioning effect on the last horizontal slice of the system tensor $\mathcal{A}$, and the preconditioners in ^[16,17,18] have no preconditioning effect on the first horizontal slice of $\mathcal{A}$. Thus, in this work, we will construct a new preconditioner for the multi-linear system (1.1), which has preconditioning effect on all horizontal slices of $\mathcal{A}$ and ameliorates the numerical performances of the preconditioners $P_{\alpha}$ ^[14] and $\tilde{P}_{\alpha}$ ^[16,17]. Inspired by the constructions of the preconditioners in ^{[19,20,21,22]}, and based on the fact that utilizing more entries of the majorization matrix of coefficient tensor in the preconditioner can make it more efficient, we design a new preconditioner $P_{\alpha, \beta} = I+G_{\alpha}+R_{\beta}$. It is constructed by introducing some elements of the first row of the majorization matrix of the system tensor $\mathcal{A}$ into the preconditioner $\tilde{P}_{\alpha} = I+G_{\alpha}$ ^[16], where

and ${\beta} = (\beta_{t}), {{\alpha}} = (\alpha_{t})\in{\mathbb{R}^{n-1}}$, $\alpha_{t}, \beta_{t}\in{\mathbb{R}}$ $(t = 2, \cdots, n)$. The new preconditioner $P_{\alpha, \beta} = I+G_{\alpha}+R_{\beta}$ has the form

For the preconditioned multi-linear system with the new preconditioner $P_{\alpha, \beta}$, we propose three Gauss-Seidel type splittings of preconditioned system tensor, then establish three new preconditioned Gauss-Seidel type methods. We first prove that the preconditioned coefficient tensor is a strong $\mathcal{M}$-tensor, then investigate the convergence properties of the three proposed preconditioned Gauss-Seidel type methods to show that they are convergent. In the meanwhile, we establish several comparison theorems, which indicate that the convergence factors of the proposed preconditioned Gauss-Seidel type methods are monotonically decreasing and the last one is the most efficient. Also, the third proposed preconditioned Gauss-Seidel type method has a faster convergence rate than the original Gauss-Seidel method and the one in ^[17]. In addition, we also propose the monotonic theorem to illustrate that for a given parameter $\alpha$, the convergence factor of the third proposed preconditioned Gauss-Seidel type method decreases with the parameter $\beta$ in $[0, 1]$, i.e., its convergence rate becomes faster with the parameter $\beta$ in $[0, 1]$. Finally, some numerical examples are given to validate the correctness of the theoretical results, and illustrate that the proposed preconditioned Gauss-Seidel type methods are efficient, and outperform the Gauss-Seidel method and preconditioned Gauss-Seidel ones in ^[14,17]. We also show the feasibilities of our preconditioner and methods by applying them to solve practical problems on higher-order Markov chain and directed graph.

The remainder of this paper is organized as follows. In Section 2, we review some basic notations, definitions, and conclusions, which will be used in the subsequent sections. In Section 3, we propose a new preconditioner and present the corresponding preconditioned Gauss-Seidel type methods for solving multi-linear system (1.1). We prove that the proposed methods are convergent in this section as well. In Section 4, we establish some comparison theorems and the monotonic theorem for the proposed preconditioned Gauss-Seidel type methods. In Section 5, numerical experiments are performed to show the efficiencies and superiorities of the proposed preconditioner and iteration methods. As applications, we also apply the preconditioned Gauss-Seidel type methods to solve the practical problems on higher-order Markov chain and directed graph in Section 6. Finally, in Section 7, some conclusions are given to end this paper.

2.   Preliminaries

In this section, we list some notations, definitions, and lemmas, which will be used throughout this paper.

Let ${\mathbf 0}$, $O$, and $\mathcal{O}$ denote a zero vector, a zero matrix, and a zero tensor, respectively. Let $\mathcal{A}, \mathcal{B}\in{\mathbb{R}^{[m, n]}}$, $\mathcal{A}\geq\mathcal{B}$ mean that each entry of $\mathcal{A}$ is no less than the corresponding one of $\mathcal{A}$. If all elements of $\mathcal{B}$ are nonnegative, then $\mathcal{B}$ is called a nonnegative tensor. The set of all $m$-order $n$-dimensional nonnegative tensors is denoted by $\mathbb{R}_{+}^{[m, n]}$. Let ${\rho(\mathcal{B})}$ = $\max\{|\lambda||\lambda\in{\sigma(\mathcal{B})}\}$ be the spectral radius of $\mathcal{B}$, where $\sigma(\mathcal{B})$ is the set of all eigenvalues of $\mathcal{B}$. Let $m\geq 2$ be an integer, and a unit tensor $\mathcal{I}_{m}$ = ($\mu_{j_{1}\cdots j_{m}}$)$\in{\mathbb{C}^{[m, n]}}$ is given by

Below some needed definitions are listed.

Definition 1. ^[28] Let $\mathcal{T}\in{\mathbb{R}^{[m, n]}}$. A pair ($\zeta$, ${\mathbf z}$)$\in{\mathbb{C}\times(\mathbb{C}^{n}\setminus{{\mathbf 0}}})$ is called an eigenpair of $\mathcal{T}$ if they satisfy the equation

where ${\mathbf z}^{[m-1]} = (z_{1}^{m-1}, \cdots, z_{n}^{m-1})^{T}$. Further, ($\zeta$, ${\mathbf z}$) is called an $H$-eigenpair of $\mathcal{T}$ if both $\zeta$ and ${\mathbf z}$ are real.

Definition 2. ^[7] Let $\mathcal{U}\in{\mathbb{R}^{[m, n]}}$. $\mathcal{U}$ is called a $\mathcal{Z}$-tensor if its non-diagonal entries are nonpositive, $\mathcal{U}$ is called an $\mathcal{M}$-tensor if there exist a tensor $\mathcal{Q}\geq \mathcal{O}$ and a positive constant $\delta\geq\rho(\mathcal{Q})$ such that

In addition, if $\delta > \rho(\mathcal{Q})$, then $\mathcal{U}$ is called a strong $\mathcal{M}$-tensor.

Definition 3. ^[13] Let $\mathcal{U}\in{\mathbb{R}^{[m, n]}}$. The majorization matrix $M(\mathcal{U})$ of $\mathcal{U}$ is defined as an $n\times n$ matrix with its entries

Definition 4. ^[34] Let $F\in{\mathbb{R}^{[2, n]}}$ and $\mathcal{U} = (u_{i_{1}\cdots i_{m}})\in{\mathbb{R}^{[m, n]}}$. The matrix-tensor product $\mathcal{T} = F\mathcal{U}\in{\mathbb{R}^{[m, n]}}$ is defined by

The above equation can be written as $\mathcal{T}_{(1)} = (F\mathcal{U})_{(1)} = F\mathcal{U}_{(1)}$, where $\mathcal{T}_{(1)}$ and $\mathcal{U}_{(1)}$ are the matrices obtained from $\mathcal{T}$ and $\mathcal{U}$ flattened along the first index ^[27].

Definition 5. ^[29] Let $\mathcal{W}\in{\mathbb{R}^{[m, n]}}$. If $M(\mathcal{W})$ is nonsingular and $\mathcal{W} = M (\mathcal{W})\mathcal{I}_{m}$, we call $M(\mathcal{W})^{-1}$ the order 2 left-inverse of $\mathcal{W}$.

Definition 6. ^[13] Let $\mathcal{W}\in{\mathbb{R}^{[m, n]}}$. If $\mathcal{W}$ has an order 2 left-inverse, then $\mathcal{W}$ is called a left-invertible tensor.

Definition 7. ^[13] Let $\mathcal{W}$, $\mathcal{E}$, $\mathcal{F}\in{\mathbb{R}^{[m, n]}}$. $\mathcal{W} = \mathcal{G}-\mathcal{H}$ is said to be a splitting of $\mathcal{W}$ if $\mathcal{G}$ is left-invertible; a convergent splitting if $\rho( M (\mathcal{G})^{-1}\mathcal{H})\leq1$; a regular splitting of $\mathcal{W}$ if $M(\mathcal{G})^{-1}\geq O$ and $\mathcal{H}\geq\mathcal{O}$; and a weak regular splitting if $M(\mathcal{G})^{-1}\geq O$ and $M(\mathcal{G})^{-1}\mathcal{H}\geq\mathcal{O}$.

Definition 8. ^[7] A tensor $\mathcal{A}\in{{\mathbb{R}^{[m, n]}}}$ is said to be a semi-positive tensor, if there exists a vector ${\mathbf x} > {\mathbf 0}$ such that $\mathcal{A}{\mathbf x}^{m-1} > {\mathbf 0}$.

Additionally, some significant lemmas are reviewed in the following.

Lemma 1. ^[13] If $\mathcal{W}\in{\mathbb{R}^{[m, n]}}$ is a $\mathcal{Z}$-tensor, then the following conditions are equivalent:

(1) $\mathcal{W}$ is a strong $\mathcal{M}$-tensor;

(2) There exists ${\mathbf y}\geq {\mathbf 0}$ such that $\mathcal{W}{\mathbf y}^{m-1} > {\mathbf 0}$;

(3) $\mathcal{W}$ has a convergent (weak) regular splitting;

(4) All (weak) regular splittings of $\mathcal{W}$ are convergent.

Lemma 2. ^[14] Let $\mathcal{T}\in{\mathbb{R}_{+}^{[m, n]}}$, and it holds that

and

Lemma 3. ^[8] If $\mathcal{T}\in{\mathbb{R}^{[m, n]}}$ is a strong $\mathcal{M}$-tensor, then for every ${\mathbf b} > {\mathbf 0}$, the multi-linear system $\mathcal{T}{\mathbf z}^{m-1} = {\mathbf b}$ has a unique positive solution.

Lemma 4. ^[19] Let $\mathcal{T}\in{\mathbb{R}^{[m, n]}}$ be a strong $\mathcal{M}$-tensor. If $\mathcal{T} = \mathcal{G}-\mathcal{H}$ is a splitting with $\mathcal{G}$ being a $\mathcal{Z}$-tensor and $\mathcal{H}$ being nonnegative, then the splitting $\mathcal{T} = \mathcal{G}-\mathcal{H}$ is a convergence regular splitting.

Lemma 5. ^[14] Let $\mathcal{T}\in{\mathbb{R}^{[m, n]}}$ be a strong $\mathcal{M}$-tensor and $\mathcal{T} = \mathcal{G}_{1}-\mathcal{H}_{1} = \mathcal{G}_{2}-\mathcal{H}_{2}$ be two weak regular splittings with $M(\mathcal{G}_{2})^{-1}\geq M(\mathcal{G}_{1})^{-1}$. If the Perron vector ${\mathbf z}$ of $M(\mathcal{G}_{2})^{-1}\mathcal{H}_{2}$ satisfies ${\mathcal{T}{\mathbf z}^{m-1}}\geq{\mathbf 0}$, then

Lemma 6. ^[14] Let $\mathcal{T}\in{\mathbb{R}}^{[m, n]}$ and $\mathcal{T} = \mathcal{U}_{1}-\mathcal{V}_{1} = \mathcal{U}_{2}-\mathcal{V}_{2}$ be a weak regular splitting and a regular splitting, respectively. If $\mathcal{F}_{2}\leq \mathcal{F}_{1}$, $\mathcal{F}_{2}\neq \mathcal{O}$, then one of the following statements holds:

(1) $\rho(M(\mathcal{U}_{2})^{-1}\mathcal{V}_{2})\leq \rho(M(\mathcal{U}_{1})^{-1}\mathcal{V}_{1}) < 1$.

(2) $\rho(M(\mathcal{U}_{2})^{-1}\mathcal{V}_{2})\geq \rho(M(\mathcal{U}_{1})^{-1}\mathcal{V}_{1})\geq1$. Furthermore, if $\mathcal{U}_{2} < \mathcal{V}_{1}$, $\mathcal{U}_{2}\neq \mathcal{O}$ and $\rho(M(\mathcal{U}_{1})^{-1}\mathcal{V}_{1}) > 1$, the first inequality is strict.

Lemma 7. ^[7] A semi-positive $\mathcal{Z}$-tensor has all positive diagonal entries.

Lemma 8. ^[13] If $\mathcal{A}$ is a strong $\mathcal{M}$-tensor, then $M(\mathcal{A})$ is a nonsingular $M$-matrix.

Lemma 9. ^[37] Let $A$ and $B$ be two nonsingular $M$-matrices. If $B\geq A$, then $A^{-1}\geq B^{-1}\geq O$.

3.   Three new preconditioned Gauss-Seidel type methods and their convergence analyses

In this section, we present a new preconditioner $P_{\alpha\beta}$, then establish three new preconditioned Gauss-Seidel type methods and their convergence properties.

We consider a new preconditioner $P_{\alpha\beta} = I+G_{\alpha}+R_{\beta}$, where

and $\alpha = (\alpha_{j})$, $\beta = (\beta_{j})$ $(j = 2, \cdots, n)$ are parameters. By premultiplying multi-linear system (1.1) by the new preconditioner $P_{\alpha\beta} = I+G_{\alpha}+R_{\beta}$, we transform the original system (1.1) into the preconditioned form

with $\mathcal{A}_{\alpha\beta} = P_{\alpha\beta}\mathcal{A} = (I+G_{\alpha}+R_{\beta})\mathcal{A}$, ${\mathbf b}_{\alpha\beta} = P_{\alpha\beta}{\mathbf b}$. Without loss of generality, we assume that $a_{j\cdots j} = 1$ $(j = 1, \cdots, n)$ as in ^{[13,14,15,16,18,19,20,21,22,23,24]}.

Let $\mathcal{A} = \mathcal{I}_{m}-\mathcal{L}-\mathcal{F}$, where $\mathcal{L} = L\mathcal{I}_{m}$, $-L$ is the strictly lower triangular part of $M(\mathcal{A})$ and $\mathcal{I}_{m}$ is a unit tensor. According to $G_{\alpha}\mathcal{L} = \mathcal{O}$, by some computations, we can obtain

Let $G_{\alpha}\mathcal{F} = \mathcal{D}_{1}-\mathcal{L}_{1}-\mathcal{F}_{1}$ and $R_{\beta}\mathcal{L} = \mathcal{D}_{2}-\mathcal{F}_{2}$, where $\mathcal{D}_{1} = D_{1}\mathcal{I}_{m}$, $\mathcal{D}_{2} = D_{2}\mathcal{I}_{m}$ and $\mathcal{L}_{1} = L_{1}\mathcal{I}_{m}$. Here, $D_{1}$ and $D_{2}$ are the diagonal parts of $M(G_{\alpha}\mathcal{F})$ and $M(R_{\beta}\mathcal{L})$, respectively, and $-L_{1}$ is the strictly lower triangular part of $M(G_{\alpha}\mathcal{F})$.

In what follows, we will study the convergence of the new preconditioned Gauss-Seidel type methods. To this end, we need the following lemma, which shows that $\mathcal{A}_{\alpha\beta}$ is a strong $\mathcal{M}$-tensor.

Lemma 10. Let $\mathcal{A}\in{\mathbb{R}}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $\alpha_{j}, \beta_{j}\in[0, 1]$ $(j = 2, \cdots, n)$, then $\mathcal{A}_{\alpha\beta}$ is a strong $\mathcal{M}$-tensor.

Proof. We first prove that $\mathcal{A}_{\alpha\beta}$ is a $\mathcal{Z}$-tensor. We denote $\mathcal{A}_{\alpha\beta} = (\widetilde{a}_{ji_{2}, \cdots, i_{m}})$. By $\mathcal{A}_{\alpha\beta} = (I+G_{\alpha}+R_{\beta})\mathcal{A}$,

Now we prove that the non-diagonal entries of $\mathcal{A}_{\alpha\beta}$ are nonpositive in the following cases:

(1) When $j = 1$, $(i_{2}, \cdots, i_{m}) = (k, \cdots.k), k\neq j, k\in\{2, \cdots, n\}$, it holds that

For $(i_{2}, \cdots, i_{m})\neq(k, \cdots.k)$, $k\in\{1, 2, \cdots, n\}$, we have

(2) When $j = 2, \cdots, n$ and $(i_{2}, \cdots, i_{m}) = (1, \cdots, 1)$, it follows that

For $(i_{2}, \cdots, i_{m}) = (k, \cdots, k), k\in\{2, \cdots, n\}$, it has

For $(i_{2}, \cdots, i_{m})\neq(k, \cdots.k), k\in\{1, 2, \cdots, n\}$, we have

Thus, for $(j, i_{2}, \cdots, i_{m})\neq(j, j, \cdots, j)$, when $0\leq\alpha_{j}\leq1$, and $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$, we have $\widetilde{a}_{ji_{2}, \cdots, i_{m}}\leq0$. Therefore, $\mathcal{A}_{\alpha\beta}$ is a $\mathcal{Z}$-tensor. Now we prove that the $\mathcal{Z}$-tensor $\mathcal{A}_{\alpha\beta}$ is a strong $\mathcal{M}$-tensor. Having in mind that $\mathcal{A}$ is a strong $\mathcal{M}$-tensor, according to Lemma 1, there exists ${\mathbf x}\geq{\mathbf 0}$ such that $\mathcal{A}{\mathbf x}^{m-1} > 0$. Since $G_{\alpha}\geq O$ and $R_{\beta}\geq O$, we can conclude that $\mathcal{A}_{\alpha\beta}{\mathbf x}^{m-1} = (I+G_{\alpha}+R_{\beta})\mathcal{A}{\mathbf x}^{m-1} = \mathcal{A}{\mathbf x}^{m-1}+G_{\alpha}\mathcal{A}{\mathbf x}^{m-1}+R_{\beta}\mathcal{A}{\mathbf x}^{m-1} > {\mathbf 0}$. Therefore, again using Lemma 1, $\mathcal{A}_{\alpha\beta}$ is a strong $\mathcal{M}$-tensor. □

Theorem 1. For $P_{\alpha\beta} = I+G_{\alpha}+R_{\beta}$ with $\alpha_{j}, \beta_{j}\in[0, 1]$ $(j = 2, \cdots, n)$, the preconditioned multi-linear system $P_{\alpha\beta}\mathcal{A}{\mathbf x}^{m-1} = P_{\alpha\beta}{\mathbf b}$ has the same unique positive solution as in $\mathcal{A}{\mathbf x}^{m-1} = {\mathbf b}.$

Proof. The conclusion of this theorem is obtained by Lemma 3 and Lemma 10. □

Below, we construct three new preconditioned Gauss-Seidel type methods

in terms of the following three Gauss-Seidel type splittings of $\mathcal{A}_{\alpha\beta}$

respectively. According to (3.2), the three preconditioned Gauss-Seidel type iteration tensors are:

In ^[17], when the relaxation factor is $\tau = 1$, the Gauss-Seidel type splitting of $\widetilde{\mathcal{A}} = (I+G_{\alpha})\mathcal{A}$ is defined as:

If $M(\check{\mathcal{E}}_{4})^{-1}$ exists, then we use $\check{\mathcal{T}} = M(\check{\mathcal{E}}_{4})^{-1}\check{\mathcal{F}}_{4}$ to denote the corresponding iteration tensor.

Next, we prove that the three Gauss-Seidel type splittings in (3.2) are convergent.

Theorem 2. Let $\mathcal{A}\in\mathbb{R}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$, then $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{1}-\check{\mathcal{F}}_{1} = \check{\mathcal{E}}_{2}-\check{\mathcal{F}}_{2} = \check{\mathcal{E}}_{3}-\check{\mathcal{F}}_{3}$ are convergence regular splittings.

Proof. We divide our proof into three cases.

(1) It is not difficult to verify that $\mathcal{L}-G_{\alpha}\mathcal{I}_{m}\geq \mathcal{O}$ and $\mathcal{F}-R_{\beta}\mathcal{I}_{m}\geq \mathcal{O}$ under the conditions $\alpha_{j}, \beta_{j}\in[0, 1]$ $(j = 2, \cdots, n)$, which means that the off-diagonal entries of $\check{\mathcal{E}}_{1} = \mathcal{I}_{m}-(\mathcal{L}-G_{\alpha}\mathcal{I}_{m})$ are non-positive, and, therefore, $\check{\mathcal{E}}_{1}$ is a $\mathcal{Z}$-tensor. Besides, $\check{\mathcal{F}}_{1} = \mathcal{F}+G_{\alpha}\mathcal{F}-R_{\beta}\mathcal{I}_{m}+R_{\beta}\mathcal{L}+R_{\beta}\mathcal{F}$ is a nonnegative tensor in view of $\mathcal{F}-R_{\beta}\mathcal{I}_{m}\geq \mathcal{O}$. Then, according to Lemma 4, $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{1}-\check{\mathcal{F}}_{1}$ is a convergence regular splitting.

(2) Inasmuch as $-L_{1}$ is the strictly lower triangular part of $M({G}_{\alpha}\mathcal{F})$ and ${G}_{\alpha}\mathcal{F}\geq\mathcal{O}$, it has $-\mathcal{L}_{1}\geq\mathcal{O}$. It follows from $\mathcal{L}-G_{\alpha}\mathcal{I}_{m}\geq \mathcal{O}$ that the off-diagonal entries of $\check{\mathcal{E}}_{2} = (\mathcal{I}_{m}-\mathcal{D}_{1})-(\mathcal{L}-G_{\alpha}\mathcal{I}_{m}-\mathcal{L}_{1})$ are non-positive, so $\check{\mathcal{E}}_{2}$ is a $\mathcal{Z}$-tensor. Owing to $\mathcal{F}-R_{\beta}\mathcal{I}_{m}\geq \mathcal{O}$, it holds that $\check{\mathcal{F}}_{2} = \mathcal{F}-\mathcal{F}_{1}-R_{\beta}\mathcal{I}_{m}+R_{\beta}\mathcal{L}+R_{\beta}\mathcal{F}\geq \mathcal{O}$. By making use of Lemma 4, we have that $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{2}-\check{\mathcal{F}}_{2}$ is a convergence regular splitting.

(3) In the same way, we can prove that $\check{\mathcal{E}}_{3} = (\mathcal{I}_{m}-\mathcal{D}_{1}-\mathcal{D}_{2})-(\mathcal{L}-G_{\alpha}\mathcal{I}_{m}-\mathcal{L}_{1})$ is a $\mathcal{Z}$-tensor, and $\check{\mathcal{F}}_{3} = \mathcal{F}-\mathcal{F}_{1}-R_{\beta}\mathcal{I}_{m}-\mathcal{F}_{2}+R_{\beta}\mathcal{F}\geq \mathcal{O}$, then $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{3}-\check{\mathcal{F}}_{3}$ is a convergence regular splitting by Lemma 4. □

4.   Comparison theorem

In this section, we first prove that the third proposed Gauss-Seidel splitting of $\mathcal{A}_{\alpha\beta}$ (the corresponding method is denoted by "$P_{\check{\mathcal{T}}_{3}}$G-S" throughout this paper) is optimal. Then, we present several comparison theorems, which show that the $P_{\check{\mathcal{T}}_{3}}$G-S method has faster convergence rate than the original Gauss-Seidel method and the one in ^[17]. Lastly, we prove that for a given $\alpha$, the convergence speed of the $P_{\alpha\beta}$G-S method increases gradually with the parameters $\beta_{j}$ ($j = 2, \cdots, n$) in the interval [0, 1].

Theorem 3. Let $\mathcal{A}\in\mathbb{R}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $0\leq\alpha_{j}\leq1$ $(j = 2, \cdots, n)$, then $\mathcal{A}_{\alpha} = \mathcal{E}_{4}-\mathcal{F}_{4}$ is a convergence regular splitting.

Proof. According to Theorem 2, we see that $\check{\mathcal{E}}_{4} = (\mathcal{I}_{m}-\mathcal{D}_{1})-(\mathcal{L}-G_{\alpha}\mathcal{I}_{m}-\mathcal{L}_{1})$ is a $\mathcal{Z}$-tensor. It is easy to verify that $\mathcal{F}_{4} = \mathcal{F}-\mathcal{F}_{1}$ is a nonnegative tensor. Then, by Lemma 4, we have that $\mathcal{A}_{\alpha} = \mathcal{E}_{4}-\mathcal{F}_{4}$ is a convergence regular splitting. □

Next, we will prove that the $P_{\mathcal{\check{T}}_{3}}$G-S method is optimal among the three proposed methods, and compare it with the preconditioned Gauss-Seidel type method according to the splitting $\widetilde{\mathcal{A}} = \check{\mathcal{E}}_{4}-\check{\mathcal{F}}_{4}$ in ^[17].

Theorem 4. Let $\mathcal{A}\in\mathbb{R}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$, then

Proof. It can be seen from Theorem 2 that $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{1}-\check{\mathcal{F}}_{1} = \check{\mathcal{E}}_{2}-\check{\mathcal{F}}_{2} = \check{\mathcal{E}}_{3}-\check{\mathcal{F}}_{3}$ are three convergence regular splittings. By straightforward computations, it holds that

i.e., $\check{\mathcal{F}}_{2}\leq\check{\mathcal{F}}_{1}$, $\check{\mathcal{F}}_{3}\leq\check{\mathcal{F}}_{1}$, and $\check{\mathcal{F}}_{3}\leq\check{\mathcal{F}}_{2}$. Below, we distinguish two cases to discuss.

1) If $\check{\mathcal{F}}_{2}\neq\mathcal{O}$, then it follows from Lemma 6 that $\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2})\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1}) < 1$. If $\check{\mathcal{F}}_{3} = \mathcal{O}$, then $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3}) = 0$, which together with $\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2})\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1})$ yields that $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})\leq\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2})\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1})$. If $\mathcal{F}_{3}\neq\mathcal{O}$, then it holds that $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})\leq\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2}) < 1$ in terms of Lemma 6, which combines with $\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2})\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1})$ and leads to $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})\leq\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2})\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1})$.

2) If $\check{\mathcal{F}}_{2} = \mathcal{O}$, then it follows from $\mathcal{O}\leq\check{\mathcal{F}}_{3}\leq\check{\mathcal{F}}_{2}$ that $\check{\mathcal{F}}_{3} = \mathcal{O}$, and, hence, $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3}) = \rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2}) = 0\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1})$.

By summarizing the above discussions, we obtain the conclusion of this theorem. □

Remark 4.1 Theorem 4 shows that $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})\leq\rho(M(\check{\mathcal{E}}_{2})^{-1}\check{\mathcal{F}}_{2})\leq\rho(M(\check{\mathcal{E}}_{1})^{-1}\check{\mathcal{F}}_{1})$, which implies that the convergence factors of the three proposed preconditioned Gauss-Seidel type methods are monotonically decreasing, and the $P_{\mathcal{\check{T}}_{3}}$G-S method has the fastest convergence speed among the proposed methods.

Note that the three proposed different preconditioned Gauss-Seidel type methods are proved to be similar. From the results of numerical experiments in Sections 5 and 6, we see that there is not much difference between the three methods. Thus, in the following analyses and theorems, we mainly discuss the properties and results of the third proposed preconditioned Gauss-Seidel type method, which is denoted by "$P_{\check{\mathcal{T}}_{3}}$G-S method".

Now, we give a comparison theorem to compare the convergence factors between the original Gauss-Seidel method and the $P_{\check{\mathcal{T}}_{3}}$G-S method.

Theorem 5. Let $\mathcal{A}\in\mathbb{R}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$. Assume that $M(\mathcal{I}_{m}-\mathcal{L})^{-1}\mathcal{F}$ is irreducible, then we have

Proof. By assumptions, we can deduce that $\mathcal{A} = (\mathcal{I}_{m}-\mathcal{L})-\mathcal{F}$ is a convergent regular splitting. Since $M(\mathcal{I}_{m}-\mathcal{L})^{-1}\mathcal{F}$ is irreducible, by the strong Perron-Frobenius theorem ^[35], there exists a positive Perron vector ${\mathbf x}$ of $M(\mathcal{I}_{m}-\mathcal{L})^{-1}\mathcal{F}$ such that

where $0 < \lambda = \rho(M(\mathcal{I}_{m}-\mathcal{L})^{-1}\mathcal{F}) < 1$. Then it holds that

and

Since $\mathcal{A} = (\mathcal{I}_{m}-\mathcal{L})-\mathcal{F}$ and $M(\mathcal{I}_{m}-\mathcal{L})^{-1}\mathcal{F}{\mathbf x}^{m-1} = \lambda{\mathbf x}^{[m-1]}$, it has $\frac{1}{\lambda}\mathcal{F}{\mathbf x}^{m-1} = (\mathcal{I}_{m}-\mathcal{L}){\mathbf x}^{m-1}$ and, therefore,

By pre-multiplying equation (4.3) by $R_{\beta}$, we have

Hence the combination of (4.1)–(4.4) results in

By Theorem 2, we know that $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{3}-\check{\mathcal{F}}_{3}$ is a convergence regular splitting, so $M(\check{\mathcal{E}}_{3})^{-1}\geq O$. Since $\mathcal{D}_{1}, \mathcal{D}_{2}, -\mathcal{L}_{1}, R_{\beta}$, and $\mathcal{F}$ are nonnegative, it follows from (4.5) that $\check{\mathcal{T}}_{3}{\mathbf x}^{m-1}\leq\lambda{\mathbf x}^{m-1}$. By Lemma 2, we have $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})\leq\rho(M(\mathcal{I}_{m}-\mathcal{L})^{-1}\mathcal{F})$, which completes the proof of this theorem. □

Remark 4.2. Theorem 5 implies that the convergence rate of the $P_{\check{\mathcal{T}}_{3}}$G-S method is faster than the original Gauss-Seidel method, which indicates that the new preconditioner $P_{\alpha, \beta}$ can improve the computational efficiency of the original Gauss-Seidel method.

In the sequel, we prove that the convergence rate of the $P_{\check{\mathcal{T}}_{3}}$G-S method is faster than that of the preconditioned Gauss-Seidel type method in ^[17]. To this end, we start with a useful lemma below.

Lemma 11. Let $\mathcal{A}\in{\mathbb{R}}^{[m, n]}$ be a strong $\mathcal{M}$-tensor. If $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$, then

Proof. According to Lemma 10, we know that $\mathcal{A}_{\alpha\beta}$ is a strong $\mathcal{M}$-tensor. Then, by making use of Lemma 1, there exists ${\mathbf x}\geq{\mathbf 0}$ such that $\mathcal{A}_{\alpha\beta}{\mathbf x}^{m-1} > {\mathbf 0}$. Since $\check{\mathcal{F}}_{3}\geq\mathcal{O}$, we can get $\check{\mathcal{E}}_{3}\geq \mathcal{A}_{\alpha\beta}$, and, hence,

It is known from Theorem 4 that $\check{\mathcal{E}}_{3}$ is a $\mathcal{Z}$-tensor, thus, we can conclude from Lemma 1 that $\check{\mathcal{E}}_{3}$ is a strong $\mathcal{M}$-tensor. In the same way, we can prove that $\check{\mathcal{E}}_{4}$ is a strong $\mathcal{M}$-tensor. According to Lemma 8, both $M(\check{\mathcal{E}}_{3})$ and $M(\check{\mathcal{E}}_{4})$ are nonsingular $M$-matrices.

Besides, it follows from $\check{\mathcal{E}}_{3} = \mathcal{I}_{m}-\mathcal{L}+G_{\alpha}\mathcal{I}_{m}-\mathcal{D}_{1}+\mathcal{L}_{1}-\mathcal{D}_{2}$, $\check{\mathcal{E}}_{4} = \mathcal{I}_{m}-\mathcal{L}+G_{\alpha}\mathcal{I}_{m}-\mathcal{D}_{1}+\mathcal{L}_{1}$ that $\check{\mathcal{E}}_{4}\geq\check{\mathcal{E}}_{3}$, which leads to $M(\check{\mathcal{E}}_{4})\geq M(\check{\mathcal{E}}_{3})$. Having in mind that $M(\check{\mathcal{E}}_{3})$ and $M(\check{\mathcal{E}}_{4})$ are nonsingular $M$-matrices, then

in view of Lemma 9. The proof is completed. □

Theorem 6. Let $\mathcal{A}\in\mathbb{R}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$. Suppose that $\overline{{\mathbf x}}$ is the Perron vector of $M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3}$, and if $\overline{{\mathbf x}}$ satisfies $\mathcal{A}\overline{{\mathbf x}}^{m-1}\geq{\mathbf 0}$, then

Proof. For $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{3}-\check{\mathcal{F}}_{3}$ and $\widetilde{\mathcal{A}} = \mathcal{E}_{4}-\mathcal{F}_{4}$, we consider the following two splittings of the tensor $\mathcal{A}$:

where $\mathcal{E}_{3}^{'} = P_{\alpha, \beta}^{-1}\check{\mathcal{E}}_{3}$, ${\mathcal{F}}_{3}^{'} = P_{\alpha, \beta}^{-1}\check{\mathcal{F}}_{3}$, $\mathcal{E}_{4}^{'} = P_{\alpha}^{-1}\check{\mathcal{E}}_{4}$, and ${\mathcal{F}}_{4}^{'} = P_{\alpha}^{-1}\check{\mathcal{F}}_{4}$. Then, we can deduce that

and

Theorems 2 and 3 prove that both $\mathcal{A}_{\alpha\beta} = \check{\mathcal{E}}_{3}-\check{\mathcal{F}}_{3}$ and $\widetilde{\mathcal{A}} = \check{\mathcal{E}}_{4}-\check{\mathcal{F}}_{4}$ are convergence regular splittings, then it follows that

which shows that $\mathcal{A} = \mathcal{E}_{3}^{'}-\mathcal{F}_{3}^{'} = \mathcal{E}_{4}^{'}-\mathcal{F}_{4}^{'}$ are two weak regular splittings of the strong $\mathcal{M}$-tensor $\mathcal{A}$.

From Lemma 11, we get $M(\check{\mathcal{E}}_{3})^{-1}\geq M(\check{\mathcal{E}}_{4})^{-1}$, which together with $P_{\alpha, \beta}\geq P_{\alpha}$ leads to

Inasmuch as $\mathcal{A} = \mathcal{E}_{3}^{'}-\mathcal{F}_{3}^{'} = \mathcal{E}_{4}^{'}-\mathcal{F}_{4}^{'}$ are two convergence weak regular splittings, $M(\mathcal{E}_{3}^{'})^{-1}\geq M(\mathcal{E}_{4}^{'})^{-1}$, and the Perron vector $\overline{{\mathbf x}}$ of $M(\mathcal{E}_{3}^{'})^{-1}\mathcal{F}_{3}^{'}$ satisfies $\mathcal{A}\overline{{\mathbf x}}^{m-1}\geq0$, we have

in view of Lemma 5, which is equivalent to

□

Finally, we will give a conclusion that for a given $\alpha$, the spectral radius of the iteration tensor of the $P_{\check{\mathcal{T}}_{3}}$G-S method decreases with $\beta$.

Theorem 7. Let $\mathcal{A}\in\mathbb{R}^{[m, n]}$ be a strong $\mathcal{M}$-tensor and $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$. Suppose that $\overline{{\mathbf x}}$ is the Perron vector of $M(\mathcal{E}_{3_{\alpha\beta^{2}}})^{-1}\mathcal{F}_{3_{\alpha\beta^{2}}}$, if $\overline{{\mathbf x}}$ satisfies $\mathcal{A}\overline{{\mathbf x}}^{m-1}\geq{\mathbf 0}$, $\alpha_{j}, \beta_{j}^{1}, \beta_{j}^{2}\in[0, 1]$ $(j = 2, \cdots, n)$ and $\beta^{1} = (\beta_{2}^{1}, \cdots, \beta_{n}^{1})\leq\beta^{2} = (\beta_{2}^{2}, \cdots, \beta_{n}^{2})$, then the following inequality holds

Proof. Since $\check{\mathcal{E}}_{3} = \mathcal{I}_{m}-\mathcal{L}+G_{\alpha}\mathcal{I}_{m}-\mathcal{D}_{1}+\mathcal{L}_{1}-\mathcal{D}_{2}$ and $\beta^{1}\leq\beta^{2 }$, we can deduce that $\check{\mathcal{E}}_{3_{\alpha\beta^{2}}}\leq\check{\mathcal{E}}_{3_{\alpha\beta^{1}}}$, which yields that $M(\check{\mathcal{E}}_{3_{\alpha\beta^{2}}})\leq M(\check{\mathcal{E}}_{3_{\alpha\beta^{1}}})$. According to Lemma 11, we see that $\check{\mathcal{E}}_{3}$ is a strong $\mathcal{M}$-tensor when $0\leq\alpha_{j}\leq1$, $0\leq\beta_{j}\leq1$ $(j = 2, \cdots, n)$. Then, from Lemma 8, we obtain that $M(\check{\mathcal{E}}_{3_{\alpha\beta^{2}}})$ and $M(\check{\mathcal{E}}_{3_{\alpha\beta^{1}}})$ are nonsingular $M$-matrices. Therefore, $M(\check{\mathcal{E}}_{3_{\alpha\beta^{2}}})^{-1}\geq M(\check{\mathcal{E}}_{3_{\alpha\beta^{1}}})^{-1}\geq O$ in terms of Lemma 9.

Besides, we derive

which means that $\mathcal{A}_{\alpha\beta^{2}}\overline{{\mathbf x}}^{m-1}\geq\mathcal{A}_{\alpha\beta^{1}}\overline{{\mathbf x}}^{m-1}\geq{\mathbf 0}$. This combines with $M(\check{\mathcal{E}}_{3_{\alpha\beta^{2}}})^{-1}\geq M(\check{\mathcal{E}}_{3_{\alpha\beta^{1}}})^{-1}$ and results in

i.e.,

So we can get $M(\check{\mathcal{E}}_{3_{\alpha\beta^{1}}})^{-1}\check{\mathcal{F}}_{3_{\alpha\beta^{1}}}\overline{{\mathbf x}}^{m-1}\geq M(\check{\mathcal{E}}_{3_{\alpha\beta^{2}}})^{-1}\check{\mathcal{F}}_{3_{\alpha\beta^{2}}}\overline{{\mathbf x}}^{m-1}$. Inasmuch as $\overline{{\mathbf x}}\neq {\mathbf 0}$ is the nonnegative Perron vector of the nonnegative tensor $M(\check{\mathcal{E}}_{3_{\alpha\beta^{2}}})^{-1}\check{\mathcal{F}}_{3_{\alpha\beta^{2}}}$, it has

According to Lemma 5, we have

□

Remark 4.3 From Theorem 3.5 of ^[17], we observe that $\rho(M(\check{\mathcal{E}}_{4})^{-1}\check{\mathcal{F}}_{4})$ has the smallest value as $\alpha_{j} = 1$ $(j = 2, \cdots, n)$. Thus we adopt $\alpha_{j} = 1$ $(j = 2, \cdots, n)$ in $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})$. Besides, it can be concluded from Theorem 7 that for given $\alpha_{i}$ $(i = 2, \cdots, n)$, $\rho(M(\check{\mathcal{E}}_{3})^{-1}\check{\mathcal{F}}_{3})$ is monotonically decreasing with $\beta_{j}$ ($j = 2, \cdots, n$) in the interval [0, 1]. Thus, we can choose $\alpha_{j} = \beta_{j} = 1$ ($j = 2, \cdots, n$) for the $P_{\check{\mathcal{T}}_{3}}$G-S method in practical computations.

5.   Numerical experiments

This section provides three numerical examples to illustrate the correctness of the theoretical results and the efficiencies of the proposed methods. To show the advantages of the proposed preconditioned Gauss-Seidel type methods, we also compare the convergence behaviors of the proposed methods with those of the original Gauss-Seidel, and the existing preconditioned Gauss-Seidel ones in ^[14,17], with respect to the number of iteration steps (IT) and the computational time in seconds (CPU).

All numerical results are computed in MATLAB (version R2016b) on a personal computer with Intel (R) Pentium (R) CPU G3240T 2.870 GHz, 16.0 GB memory and Windows 10 system. We use the power method in ^[33] to compute the spectral radii of all tested preconditioned Gauss-Seidel type methods. Besides, all iterations are terminated once the residual(RES) satisfies $\mathrm{RES}: = \|\mathcal{A}{\mathbf y}_{k}^{m-1}-{\mathbf b}\|_{2}\leq10^{-10}$, or IT exceeds the maximal number of iteration steps $10,000$. According to Remark 3.7 in ^[17] and Remark 5.7 in ^[14], we see that $\rho(\mathcal{T}_{\alpha})$ and $\rho(\check{\mathcal{T}})$ attain minimum values as $\alpha_{j} = 1$ ($j = 2, \ldots, n$). Moreover, based on Theorem 7 and Remark 4.3, we can choose $\alpha_{l} = 1$ and $\beta_{l} = 1$ ($l = 2, \ldots, n$) in the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method. Then, we compare the numerical results of the preconditioned Gauss-Seidel type methods with all parameters being 1. Also, we list the numerical experimental results of the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method when the parameters $\alpha$ and $\beta$ are not equal, and give a reasonable suggestion of the values of the parameters in the $P_{\mathcal{\check{T}}_{3}}$G-S method. The product $\mathcal{A}{\mathbf x}^{m-1}$ is computed by transforming into the matrix-vector product $\mathcal{A}{\mathbf x}^{m-1} = \mathcal{A}_{(1)}(\underbrace{{\mathbf x}\otimes\cdots\otimes{\mathbf x}}_{m})$, where $\otimes$ is the Kronecker product. The product of a matrix $U$ and a tensor $\mathcal{V}$ is computed by $(U\mathcal{V})_{(1)} = U\mathcal{V}_{(1)}$. In our computations, the preconditioned Gauss-Seidel type methods in ^[14] and ^[17] are denoted by $P_{\mathcal{T}_{\alpha}}$G-S and $P_{\check{\mathcal{T}}}$G-S, respectively. The three proposed preconditioned Gauss-Seidel type methods are denoted by $P_{\mathcal{\check{T}}_{1}}$G-S, $P_{\mathcal{\check{T}}_{2}}$G-S, and $P_{\mathcal{\check{T}}_{3}}$G-S, respectively.

Example 1. ^[32] Let $\mathcal{A}\in \mathbb{R}^{[3, n]}$ and ${\mathbf b}\in \mathbb{R}^{n}$, where

and

This example comes from the high order bellman equations in ^[32] with few modifications. It can be verified that $\mathcal{A}$ is a strong $\mathcal{M}$-tensor. Both the righthand side vector ${\mathbf b}$ and the initial vector ${\mathbf y}_{0}$ are chosen to be $ones(n, 1)$.

Example 2. ^[8,22] Consider the multi-linear system (1.1) with coefficient tensor $\mathcal{A} = s\mathcal{I}_{3}-\mathcal{B}$ and righthand side vector ${\mathbf b} = (1, 1, \cdots, 1)^{T}$, where $\mathcal{B}$ is generated randomly by Matlab, and $s = 1.01\max\limits_{1\leq j\leq n}(B{\mathbf e}^{2})_{j}$ with ${\mathbf e} = (1, 1, \cdots, 1)^{T}$.

Example 3. ^[23] Let $\mathcal{A}\in \mathbb{R}^{[3, n]}$ and ${\mathbf b}\in \mathbb{R}^{n}$ with

and $b_{i} = 1$ $(i = 1, \cdots, n)$.

This example comes from ^[23] with few modifications. We can verify that $\mathcal{A}$ is a strictly diagonally dominant $\mathcal{Z}$-tensor with positive diagonal elements, then $\mathcal{A}$ is a strong $\mathcal{M}$-tensor. The initial vector is taken as ${\mathbf y}_{0} = ones(n, 1)$.

We compare the convergence factors of the tested Gauss-Seidel and preconditioned Gauss-Seidel type methods for Examples 1, 2, and 3 in Tables 1, 5, and 9, respectively. As observed in Tables 1, 5, and 9, the convergence factors of the preconditioned Gauss-Seidel type methods are smaller than that of the original Gauss-Seidel one, which reveals that applying the preconditioning technique can accelerate the tensor splitting methods and is in accordance with the result in Theorem 5. Also, the convergence factors of the three proposed preconditioned Gauss-Seidel type methods satisfy $\rho(\mathcal{\check{T}}_{3})\leq\rho(\mathcal{\check{T}}_{2})\leq\rho(\mathcal{\check{T}}_{1})$ for all cases. This coincides with the conclusion in Theorem 4. Moreover, $\rho(\mathcal{\check{T}}_{2})$ and $\rho(\mathcal{\check{T}}_{3})$ are smaller than $\rho(\mathcal{T}_{\alpha})$ ^[14] and $\rho(\check{\mathcal{T}})$ ^[17], which demonstrates the correctness of Theorem 6. We conclude that the proposed $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods have faster convergence rates than the existing ones in ^[14,17].

To reflect the changing of $\rho(\mathcal{\check{T}}_{3})$ with respect to parameters, in Tables 2, 6, and 10, we report the convergence factor $\rho(\mathcal{\check{T}}_{3})$ and the corresponding IT of the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method with $\alpha = \beta$ from $0.1$ to $1$ with step size $0.1$ for Examples 1–3, respectively. The corresponding IT of the $P_{\mathcal{\check{T}}_{3}}$G-S method are listed in brackets. Numerical results in Tables 2, 6, and 10 show that the convergence factor and IT of the $P_{\mathcal{\check{T}}_{3}}$G-S method decrease with increasing of $\alpha = \beta$ in [0, 1] for Example 2. However, $\rho(\mathcal{\check{T}}_{3})$ increases first with respect to $\alpha = \beta$ in [0, 0.9], and then decreases with respect to $\alpha = \beta$ in [0.9, 1]. The reason is that the assumption $\mathcal{A}\overline{{\mathbf x}}^{m-1}\geq{\mathbf 0}$ in Theorem 7 is invalid. Besides, the IT and convergence factor of the $P_{\mathcal{\check{T}}_{3}}$G-S method are monotonically increasing with the problem size $n$ for Example 1, while it is opposite for Example 3.

Also, we list the convergence factor $\rho(\mathcal{\check{T}}_{3})$ and IT of the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method with different values of $(\alpha, \beta)$ and $\alpha, \beta$ from $0.1$ to $1$ with step size $0.1$ for Examples 1–3 in Tables 3, 7, and 11, respectively. The best results are marked in bold. It can be observed in Tables 3, 7, and 11 that for a given $\alpha$, $\rho(\mathcal{\check{T}}_{3})$ and IT of the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method become smaller with $\beta$, which coincides with the result in Theorem 7. While for a given $\beta$, $\rho(\mathcal{\check{T}}_{3})$ and IT of the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method are not monotonically decreasing with $\alpha$ for all cases. In addition, the $P_{\mathcal{\check{T}}_{3}}$G-S method has the highest computational efficiency as $\alpha = 0.1$ and $\beta = 1$. This reveals that we can adopt small values of $\alpha$ and large values of $\beta$ for the $P_{\mathcal{\check{T}}_{3}}$G-S method in numerical experiments and practical computations.

Additionally, numerical results of the tested Gauss-Seidel and preconditioned Gauss-Seidel methods for Examples 1, 2, and 3 with different problem sizes are reported in Tables 4, 8, and 12, respectively. From Tables 4, 8, and 12, we observe that all tested methods can successfully compute approximate solutions satisfying the prescribed stopping criterion, and their IT and CPU time increase with increasing of the problem size $n$. For Examples 1 and 3, the proposed $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods have advantages over the other ones from the point of view of IT and CPU time, and the $P_{\mathcal{\check{T}}_{3}}$G-S method performs the best among the tested methods. For Example 2, the three proposed methods outperform the other ones except for the case of $n = 2$, and the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method is the most effective method in all tested methods, which is in accordance with the results in Theorems 5 and 6. Another observation is that for Examples 1 and 3, the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method is more stable than the other ones since the variation amplitude of its IT is the smallest among the tested methods.

In order to further show the efficiencies of the proposed preconditioned Gauss-Seidel type algorithms, we also apply the tested algorithms to solve the multi-linear system in Example 3 with larger scales $n = 100,150,200,250,300$, and the obtained numerical results are listed in Table 13. From Table 13, it is observed that when $n$ increases, the IT and CPU time of the tested methods increase. The proposed $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods return better numerical performances than the other ones in terms of IT and CPU time, and the $P_{\mathcal{\check{T}}_{3}}$G-S method performs the best among the tested methods. The advantages of the $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods become more pronounced as the problem size $n$ becomes larger. This shows that the proposed methods are feasible and efficient for solving multi-linear systems with large scales.

To investigate the dependence of the tested methods on the problem size and parameter, the comparisons for the convergence factors of the tested methods and the changes of $\rho(\mathcal{\check{T}}_{3})$ with $\alpha = \beta$ for Examples 1–3 are depicted in Figures 2, 4, and 6, respectively. We plot the RES(log10) curves of the tested algorithms for Examples 1–3 in Figures 1, 3, and 5, respectively, to show the advantages of the proposed methods. Here, "PGS($\mathrm{S}\alpha$)", "PGS($\mathrm{G}\alpha$)", and "PGS-$j$" $(j = 1, 2, 3)$ denote the $P_{\mathcal{T}_{\alpha}}$G-S, $P_{\check{\mathcal{T}}}$G-S, and $P_{\mathcal{\check{T}}_{j}}$G-S $(j = 1, 2, 3)$ methods, respectively. Their convergence factors are represented by $\rho(\mathrm{S}\alpha)$, $\rho(\mathrm{G}\alpha)$, and $\rho j$ $(j = 1, 2, 3)$, respectively. It can be seen from Figures 1, 3 and 5 that the $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods are superior to the other ones in view of IT, and the $P_{\mathcal{\check{T}}_{1}}$G-S and $P_{\mathcal{T}_{\alpha}}$G-S ^[14] methods are comparable for Examples 1 and 3. In addition, for Example 2, the three proposed preconditioned Gauss-Seidel type methods outperform the other ones. Also, the $P_{\mathcal{T}_{\alpha}}$G-S method performs better than the $P_{\check{\mathcal{T}}}$G-S one, while it is opposite for Examples 1 and 3. This is in consistent with the results in Tables 4, 8 and 12. On the other hand, Figures 2, 4 and 6 clearly show that $\rho(\mathcal{\check{T}}_{2})$ and $\rho(\mathcal{\check{T}}_{3})$ are less than the other ones, and $\rho(\mathcal{\check{T}}_{3})$ is the minimum one among all convergence factors except for $n = 2$ in Example 2. Meanwhile, $\rho(\mathcal{T}_{\alpha})$ and $\rho(\mathcal{\check{T}}_{1})$, and $\rho(\check{\mathcal{T}})$ and $\rho(\mathcal{\check{T}}_{2})$ are close for Example 1. Furthermore, all convergence factors become closer with $n$ for Example 2. Lastly, $\rho(\mathcal{\check{T}}_{3})$ decreases first and then increases as $\alpha = \beta$ increase in [0, 1] for Examples 1 and 3, while it decreases with the increasing of $\alpha = \beta$ in [0, 1] for Example 2. The reason may be that the assumption $\mathcal{A}\overline{{\mathbf x}}^{m-1}\geq{\mathbf 0}$ in Theorem 7 is invalid.

6.   Applications in the higher-order Markov chain and directed graph

As applications of the proposed preconditioned Gauss-Seidel type methods for solving the multi-linear system (1), we consider the higher-order Markov chain ^[4] and directed graph ^[5,6]. Recently, Li and Ng ^[4] presented an approximate tensor model for higher-order Markov chain. The approximate tensor model for higher-order Markov chain is given as follows:

where $\mathcal{P} = (p_{j_{1}j_{2}\cdots j_{m}})$ is an $m$-th order $n$-dimensional stochastic tensor (or transition probability tensor), that is:

and ${\mathbf z} = (z_{j})$ is called a stochastic vector with $z_{j}\geq0$, $\sum\limits_{j = 1}^{n}z_{j} = 1$. The following example relating to the higher-order Markov chain comes from ^[36].

Example 4. ^[36] Let

This example comes from occupational mobility of physicists date given in ^[36]. Let $\mathcal{A} = \mathcal{I}_{m}-\theta \mathcal{P}$, with a positive real number $\theta$. In ^[13], Liu et al. gave an equivalent equation with (6.1):

We test and compare the preconditioned Gauss-Seidel type methods in this paper and in ^[14,17] for this example. According to the numerical experiments in ^[19,24], the initial vector is adopted to be ${\mathbf z}_{0} = \frac{1}{n}ones(n, 1)$, and all runs stopped once the RES satisfies $\mathrm{RES}: = \|\mathcal{P}{\mathbf z}^{m-1}-{\mathbf z}\|_{2}\leq10^{-10}$ or the IT exceeds the maximal number of iteration steps $10,000$.

Example 5. ^[5,6] Let $\Gamma$ be a directed graph with $\mathbb{V} = \{1, 2, 3, 4, 5\}$, where $\mathbb{P} = \{1, 2\}$ and $\mathbb{D} = \{3, 4, 5\}$. Then, by definition, we can get the stochastic tensor $\hat{\mathcal{Q}} = (\hat{q}_{i_{1}i_{2}i_{3}})\in\mathbb{R}^{[3,5]}$, where $\hat{q}_{221} = \hat{q}_{421} = \hat{q}_{521} = \hat{q}_{312} = \hat{q}_{412} = \hat{q}_{512} = 0$, $\hat{q}_{122} = \hat{q}_{133} = \hat{q}_{144} = \hat{q}_{155} = \hat{q}_{121} = \hat{q}_{212}$ = $\hat{q}_{211} = \hat{q}_{311} = \hat{q}_{411} = \hat{q}_{321} = \hat{q}_{112} = \hat{q}_{511} = 1$, and other elements are 1/6. We apply the tested preconditioned Gauss-Seidel type methods to solve the following multi-linear PageRank:

where $\mathcal{V} = (v_{i_{1}i_{2}\cdots i_{m}})$ with $v_{i_{1}i_{2}\cdots i_{m}} = v_{i_{1}}$, $\forall i_{2}, \cdots, i_{m}$, and ${\mathbf v} = (v_{i})$ is a stochastic vector, $\xi\in(0, 1)$ is a parameter. This example comes from ^[5,6,24] with few modifications. As in Example 4, we take the initial vector as ${\mathbf z}_{0} = \frac{1}{n}ones(n, 1)$, and all runs stopped once the RES satisfies $\mathrm{RES}: = \|\mathcal{P}{\mathbf z}^{m-1}-{\mathbf z}\|_{2}\leq10^{-10}$ or the IT exceeds the maximal number of iteration steps $10,000$.

In Tables 14–16, we list the numerical results of the tested methods for Example 4 with five different values of $\theta$. For Example 5, in Tables 17–19, we list the same items as those in Tables 14–16. In Figures 7 and 9, we plot the RES(log10) curves of the tested methods for Example 4 and Example 5, respectively. The comparisons for convergence factors of the tested methods and the variation curves of $\rho(\mathcal{\check{T}}_{3})$ with $\alpha = \beta$ for Examples 4 and 5 are displayed in Figures 8 and 10, respectively. As in Section 5, we use "PGS($\mathrm{S}\alpha$)", "PGS($\mathrm{G}\alpha$)", and "PGS-$j$" $(j = 1, 2, 3)$ to denote the $P_{\mathcal{T}_{\alpha}}$G-S, $P_{\check{\mathcal{T}}}$G-S and $P_{\mathcal{\check{T}}_{j}}$G-S $(j = 1, 2, 3)$ methods, respectively. Their convergence factors are represented by $\rho(\mathrm{S}\alpha)$, $\rho(\mathrm{G}\alpha)$ and $\rho j$ $(j = 1, 2, 3)$, respectively.

By the numerical results presented in Tables 14–19, we can conclude some observations: First, all tested algorithms are convergent for all cases. Second, the IT and CPU time of all tested algorithms increase with $\theta(\xi)$. Third, the $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods have higher computational efficiencies than the other ones, except the $P_{\mathcal{T}_{\alpha}}$G-S method is slightly better than the $P_{\mathcal{\check{T}}_{2}}$G-S one in Example 4, and the advantages become more pronounced with $\theta(\xi)$. Fourth, numerical performances of the $P_{\mathcal{T}_{\alpha}}$G-S and the $P_{\mathcal{\check{T}}_{1}}$G-S methods are comparable, and they are more efficient than the Gauss-Seidel and $P_{\check{\mathcal{T}}}$G-S ones with respect to IT and CPU time. Fifth, all convergence factors increase with $\theta(\xi)$, and $\rho(\mathcal{\check{T}}_{3})$ is less than the other ones. Also, $\rho(\mathcal{\check{T}}_{3}) < \rho(\mathcal{\check{T}}_{2}) < \rho(\mathcal{\check{T}}_{1})$ holds for all cases, which coincides with the conclusion in Theorem 4. Finally, $\rho(\mathcal{\check{T}}_{3})$ decreases gradually with increasing of $\alpha = \beta$.

From Figures 7 and 9, we observe that $P_{\mathcal{\check{T}}_{2}}$G-S and $P_{\mathcal{\check{T}}_{3}}$G-S methods have faster convergence speeds than the other ones. The convergence behaviors of the $P_{\mathcal{T}_{\alpha}}$G-S and $P_{\mathcal{\check{T}}_{1}}$G-S methods are comparable, and they perform better than the remaining two in view of IT. Additionally, it follows from Figures 8 and 10 that all convergence factors become larger and closer with $\theta(\xi)$, and $\rho(\mathcal{\check{T}}_{3})$ is the smallest in all convergence factors. Finally, $\rho(\mathcal{\check{T}}_{3})$ increases with $\theta(\xi)$ and decreases with $\alpha = \beta$. The above results are in consistent with Theorems 5–7 and Remark 4.3.

7.   Conclusions

In this paper, by introducing some elements in the first row of the majorization matrix associated with the coefficient tensor into the preconditioner $\tilde{P}_{\alpha}$ ^[16,17], we propose a new preconditioner for the multi-linear system (1.1). Compared with preconditioners in ^[14,17], the new preconditioner has preconditioning effect on all horizontal slices of the coefficient tensor and improves their computational efficiencies. Then we construct three new preconditioned Gauss-Seidel type methods for (1.1) according to three different Gauss-Seidel type tensor splittings of the preconditioned system tensor. Convergence and comparison theorems have been given to show that the proposed methods are convergent. Meanwhile, the last proposed method has a faster convergence rate than the original Gauss-Seidel method and the preconditioned Gauss-Seidel method in ^[17]. Also, it is proved that for a given $\alpha$, the convergence factor of the proposed $P_{\mathcal{\check{T}}_{3}}$G-S method decreases monotonically with $\beta$ in [0, 1]. Finally, numerical experiments are performed to illustrate that the theoretical results are correct, and the proposed methods are efficient and have better numerical performances than the Gauss-Seidel method and preconditioned Gauss-Seidel methods in ^[14,17]. As applications, we apply the proposed methods to solve the multi-linear system (1.1) arising from the higher-order Markov chain and directed graph to further show their feasibilities and superiorities.

It can be observed that comparison and monotonic theorems (Theorems 6 and 7) are valid under the condition $\mathcal{A}\overline{{\mathbf x}}^{m-1}\geq{\mathbf 0}$, which is not easy to be verified in the practical calculation; thus, we will work to weaken this assumption in the future. Last but not least, it also makes sense to study more efficient preconditioners and the related preconditioned tensor splitting methods in our future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 12361078), and the Guangxi Natural Science Foundation (No. 2024GXNSFAA010498). We would like to express our sincere thanks to the editor and the anonymous reviewers for their valuable suggestions and constructive comments which greatly improved the presentation of this paper.

Conflict of interest

The authors declare there is no conflicts of interest.