Comparison results for splitting iterations for solving multi-linear systems

Accepted Manuscript Comparison results for splitting iterations for solving multi-linear systems

Wen Li, Dongdong Liu, Seak-Weng Vong

PII: DOI: Reference:

S0168-9274(18)30157-0 https://doi.org/10.1016/j.apnum.2018.07.009 APNUM 3400

To appear in:

Applied Numerical Mathematics

Received date: Revised date: Accepted date:

1 February 2018 10 May 2018 20 July 2018

Please cite this article in press as: W. Li et al., Comparison results for splitting iterations for solving multi-linear systems, Appl. Numer. Math. (2018), https://doi.org/10.1016/j.apnum.2018.07.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Comparison Results for Splitting Iterations for Solving Multi-linear Systems Wen Li∗, Dongdong Liu†, Seak-Weng Vong‡ July 23, 2018

Abstract It is known that the spectral radius of the iterative tensor can be seen as an approximate convergence rate for solving multi-linear systems by tensor splitting iterative methods. So in this paper, first we give some spectral radius comparisons between two different iterative tensors. Then, we propose the preconditioned tensor splitting method for solving multi-linear systems, which provides an alternative algorithm with the choice of a preconditioner. In particular, also we give some spectral radius comparisons between the preconditioned iterative tensor and the original one. Numerical examples are given to demonstrate the efficiency of the proposed preconditioned methods.

Keywords. Comparison theorem, tensor splitting, spectral radius, multi-linear systems, preconditioned methods. MSC 15A48 15A69 65F10 65H10

1

Introduction

Recently, because of some applications, the following multi-linear system: Axm−1 = b,

(1.1)

has attracted more and more attention (e.g., see [4],[12],[14],[23]), where A = (aii2 ···im ) is an order m dimension n tensor, b is an n dimensional vector, and the tensor-vector product Axm−1 is defined by (Ax

m−1

)i =

n 

aii2 ···im xi2 · · · xim , i = 1, 2, ..., n,

i2 ,··· ,im =1 ∗

School of Mathematical Sciences, South China Normal University, China Department of Mathematics, University of Macau, Macau, China ‡ Department of Mathematics, University of Macau, Macau, China. †

1

(1.2)

where xi denotes the i-th component of x. Some theoretical analysis and algorithms for solving the equation (1.1) are also studied (see [4],[12],[14],[23]). In [12], the authors proposed tensor splitting A = E − F, and then established a general tensor splitting iterative method for solving (1.1) as follows: for a given initial vector x0 , 1

−1 [ m−1 ] xk = [M (E)−1 Fxm−1 , k = 1, 2, · · · . k−1 + M (E) b]

(1.3)

The tensor M (E)−1 F is called the iterative tensor of the splitting method. In particular, they discussed the convergence rate for the splitting iterative method, and showed that the spectral radius ρ(M (E)−1 F) can be seen as an approximate convergence rate of the iteration (1.3). Numerical examples given in [12] showed that ρ(M (E)−1 F) performs well as the proposed approximate convergence rate of the iteration (1.3). For matrix splitting methods, it is well-known that for the nonsingular M-matrix case, the Guass-Seidel method demonstrates faster convergence than the Jacobi method, that is to say, the spectral radius of the iterative matrix of the Guass-Seidel method is not larger than the one of the Jacobi method (see [17]). Naturally, we may ask: ”is this comparison property true for the Jacobi type and Guass-Seidel type splitting methods of a strong M-tensor?” As one of motivations in this paper, we study the comparison of the proposed approximate convergence rate for these different tensor splitting iterative methods. As we know, the preconditioned technique for solving linear systems is very important, which can be used to improve the rate of convergence of the iterations when a suitable preconditioner is chosen. However, this technique has not been employed to solve the tensor system (1.1). So in this paper, we explore the preconditioned tensor splitting methods, and develop the preconditioned Gauss-Seidel type and the SOR type iterative methods for solving (1.1), which can be used as alternative methods for solving (1.1). The main contributions of this paper are • to present some comparison theorems for spectral radii of the iterative tensors for the tensor splittings. • to propose some preconditioned iterative methods for solving the equation (1.1). The rest of this paper is organized as follows. In Section 2 we introduce some definitions and some related lemmas which will be used in the sequel. In Section 3, we present some comparison theorems for spectral radii of the iterative tensors for different tensor splittings. In Section 4, we explore a preconditioned splitting iterative method, and a special preconditioner is proposed. The corresponding theoretical analysis for the preconditioned iteration is given in Section 5. In Section 6, we proposed the preconditioned SOR type method. In Section 7, some numerical examples are given to show the efficiency of the proposed preconditioned iterations. The final section is the concluding remark. 2

2

Preliminaries

Let n be a positive integer, by n we denote the set {1, · · · , n}. A tensor A consists of n1 × · · · × nm entries in the real field R: A = (ai1 ···im ), ai1 ···im ∈ R, ij ∈ nj , j = 1, · · · , m. If n1 = · · · = nm = n, A is called an order m dimension n tensor. We denote the set of all order m dimension n tensors by R[m,n] . When m = 2, R[2,n] denotes the set of all n × n real matrices. When m = 1, R[1,n] is simplified as Rn , which is the set of all n-dimension real vectors. Similarly, the above notions can be generalized to the complex number field C. Let R+ be the nonnegative real field. If each entry of A is nonnegative, we call A a nonnegative [m,n] tensor, and the set of all the order m dimension n nonnegative tensors is denoted by R+ . Let 0, O and O denote a zero vector, a zero matrix and a zero tensor, respectively. Let A and B be two tensors (vectors or matrices) with the same sizes. The order A ≥ B(> B) means that each entry of A is no less than (larger than) corresponding one of B. Let A ∈ R[2,n] and B ∈ R[k,n] . The matrix-tensor product C = AB ∈ R[k,n] is defined by cji2 ···ik =

n 

ajj2 bj2 i2 ···ik .

(2.1)

j2 =1

The formula (2.1) can be written as follows (see [9]): C(1) = (AB)(1) = AB(1) ,

(2.2)

where C(1) and B(1) are the matrices obtained from C and B flattened along the first index (see [9, 10]). For example, if B = (bijk ) ∈ C[3,n] , then ⎡ B(1)

⎢ ⎢ =⎢ ⎢ ⎣

b111 b211 .. .

··· ··· .. .

bn11 · · ·

b1n1 b2n1 .. .

b112 b212 .. .

··· ··· .. .

bnn1 bn12 · · ·

b1n2 b2n2 .. .

··· ··· .. .

bnn2 · · ·

b11n b21n .. .

··· ··· .. .

b1nn b2nn .. .

bn1n · · ·

bnnn

⎤ ⎥ ⎥ ⎥. ⎥ ⎦

Next we introduce the following definition of tensor eigenvalues and eigenvectors, which were given by Qi [19] and Lim [15] independently. Definition 2.1. [15, 19] Let A ∈ R[m,n] . A pair (λ, x) ∈ C×(Cn \{0}) is called an eigenvalueeigenvector(or simply eigenpair) of A if they satisfy the equation Axm−1 = λx[m−1] ,

(2.3)

where x[m−1] = (xm−1 , ..., xm−1 )T . We call (λ, x) an H-eigenpair if both λ and x are real. n 1 3

Let ρ(A) = max{|λ| |λ ∈ σ(A)} be the spectral radius of A, where σ(A) is the set of all eigenvalues of A. An unit tensor Ik = (δi1 ···ik ) ∈ C[k,n] for a positive integer k ≥ 2 is given by: 1, i1 = · · · = i k , δi1 ···ik = 0, else. We call aii···i a diagonal entry of A = (ai1 ···im ) ∈ R[m,n] for ∀i ∈ n . In particular, a tensor A is called diagonal if its entries are null except the diagonal entries. In [15], Lim generalized the definition of an irreducible matrix into an irreducible tensor. Definition 2.2. [15] An A ∈ C[m,n] is called reducible if there exists a nonempty proper index subset I ⊆ {1, 2, ..., n} such that ai1 i2 ···im = 0, ∀i1 ∈ I, ∀i2 , ..., im ∈ / I. If A is not reducible, then we call A irreducible. The definitions of Z-tensors, M-tensors and strong M-tensors are introduced as follows. Definition 2.3. [3, 26] Let A ∈ R[m,n] . A is called a Z-tensor if its off-diagonal entries are non-positive. A is called an M-tensor if there exists a nonnegative tensor B and a positive real number η ≥ ρ(B) such that A = ηIm − B. If η > ρ(B), then A is called a strong M-tensor. We recall the following definitions and lemmas for the completeness of our presentation. Definition 2.4. [18] Let A ∈ C[m,n] . The majorization matrix M (A) of A is defined as an n × n matrix with its entries M (A)ij = aij...j

i, j = 1, .., n.

In [13], Liu and Li gave the left-inverse of a tensor. Definition 2.5. [13] Let A ∈ R[m,n] . If M (A) is a nonsingular matrix and A = M (A)Im , we call M (A)−1 the order 2 left-inverse of A. Based on Definition 2.5, Liu, Li and Vong [12] gave the definitions of a left-nonsingular tensor and the tensor splitting. Definition 2.6. [12] Let A ∈ R[m,n] . If A has an order 2 left-inverse, A is called a leftinvertible tensor or a left-nonsingular tensor. Definition 2.7. [12] Let A, E, F ∈ R[m,n] . A = E −F is said to be a splitting of A if E is leftnonsingular; a regular splitting of A if E is left-nonsingular with M (E)−1 ≥ O, and F ≥ O; a weak regular splitting of A if E is left-nonsingular with M (E)−1 ≥ O and M (E)−1 F ≥ O; a convergent splitting if ρ(M (E)−1 F) < 1. 4

Furthermore, they gave some equivalent conditions based on the definition of the tensor splitting. Some of these conditions are listed as follows. Lemma 2.8. If A is a Z-tensor, then the following conditions are equivalent: (1) A is a strong M-tensor. (2) There exist an inverse-positive Z-matrix B and a semi-positive Z-tensor C with A = BC. (3) A has a convergent (weak) regular splitting. (4) All (weak) regular splittings of A are convergent. For the multi-linear system (1.1), Ding and Wei [4] gave Lemma 2.9. [4] If A is a strong M-tensor, then for every positive vector b the multi-linear system Axm−1 = b has a unique positive solution.

3

The comparison theorem

In order to give comparison theorems of spectral radii of iterative tensors for the tensor splittings, we first provide some lemmas as follows. [m,n]

Lemma 3.1. [24][25] Suppose that A ∈ R+

. Then

(Axm−1 )i . x≥0,x=0 xi >0 xm−1 i

ρ(A) = max min

(3.1)

Based on Lemma 3.1, we get [m,n]

Lemma 3.2. For A ∈ R+

,

μx[m−1] ≤ (<)Axm−1 , x ≥ 0, x = 0 impies μ ≤ (<)ρ(A),

(3.2)

Axm−1 ≤ νx[m−1] , x > 0, impies ρ(A) ≤ ν.

(3.3)

and Proof. First we prove (3.2). From μx[m−1] ≤ (<)Axm−1 , x ≥ 0 we have μ ≤ (<) minxi >0 Then by (3.1), we get μ ≤ (<)ρ(A), which proves (3.2). The inequality (3.3) follows from Lemma 3.22 of [20]. [m,n]

Lemma 3.3. [20][24] Let A ∈ R+ eigenvalue λ, then λ = ρ(A).

(Axm−1 )i . xm−1 i

. If A has a positive eigenvector x corresponding to an

Next, we give the comparison theorem of spectral radii of iterative tensors for tensor splittings. 5

Theorem 3.4. Let A ∈ R[m,n] and A = E1 − F1 = E2 − F2 be a weak regular splitting and a regular splitting, respectively. If F2 ≤ F1 , F2 = O, then one of the following statements holds: (1) ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ) < 1. (2) ρ(M (E2 )−1 F2 ) ≥ ρ(M (E1 )−1 F1 ) ≥ 1. Furthermore, if F2 < F1 , F2 = O and ρ(M (E1 )−1 F1 ) > 1, the first inequality is strict. Proof. (1): Let G be in ∈ R[m,n] whose entries are all equal to 1. Let Ck = M (E1 )−1 F1 + k1 G. Because A = E1 − F1 is a weak regular splitting, M (E1 )−1 F1 ≥ O. By the proof of Theorem 2.3 in [24], we have lim ρ(M (E1 )−1 F1 + k1 G) = ρ(M (E1 )−1 F1 ) and there exists a positive k→∞

integer N such that ρ(M (E1 )−1 F1 ) ≤ ρ(M (E1 )−1 F1 + k1 G) < 1 as k > N . It is easy to check that M (E1 )−1 F1 + k1 G is positive and hence is irreducible. By the strong Perron-Frobenius theorem (see [2],[24] or [20]), for any given k > N , M (E1 )−1 F1 + k1 G has a positive Perron vector x such that 1 (M (E1 )−1 F1 + G)xm−1 = ρk x[m−1] , k where ρk = ρ(M (E1 )−1 F1 + k1 G), and then 1 M (E1 )(ρk Im − G)xm−1 = F1 xm−1 . k

(3.4)

As the assumption A = E1 − F1 , we have M (E1 ) = M (A) + M (F1 ) which together with (3.4) gives 1 1 M (A)(ρk Im − G)xm−1 = F1 xm−1 − M (F1 )(ρk Im − G)xm−1 k k 1 m−1 = (1 − ρk )M (F1 )Im x + M (F1 )Gxm−1 + (F1 − M (F1 )Im )xm−1 . k (3.5) Because F1 ≥ F2 , we have M (F1 ) ≥ M (F2 ) and F1 − M (F1 )Im ≥ F2 − M (F2 )Im . Then by (3.5) and M (A) = M (E2 ) − M (F2 ), we get 1 (M (E2 ) − M (F2 ))(ρk Im − G)xm−1 k 1 m−1 ≥(1 − ρk )M (F2 )Im x + M (F2 )Gxm−1 + (F2 − M (F2 )Im )xm−1 , k which implies that 1 M (E2 )(ρk Im − G)xm−1 ≥ F2 xm−1 . k 6

(3.6)

Because A = E2 − F2 is a regular splitting, M (E2 )−1 ≥ O. By (3.6), we have 1 0 ≤ M (E2 )−1 F2 xm−1 ≤ (ρk Im − G)xm−1 ≤ ρk x[m−1] . k

(3.7)

Then by (3.3) of Lemma 3.2, we have ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 + k1 G). For k → ∞, we get ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ) < 1. (2) Since M (E1 )−1 F1 ≥ O, there is a Perron vector x such that M (E1 )−1 F1 xm−1 = ρ(M (E1 )−1 F1 )x[m−1] , which implies: F1 xm−1 = ρ(M (E1 )−1 F1 )M (E1 )x[m−1] = ρ(M (E1 )−1 F1 )E1 xm−1 .

(3.8)

By E1 = A + F1 and (3.8), we get ρ(M (E1 )−1 F1 )Axm−1 = (1 − ρ(M (E1 )−1 F1 ))F1 xm−1 . If ρ(M (E1 )−1 F1 ) ≥ 1, by F2 ≤ F1 , it is easy to check that ρ(M (E1 )−1 F1 )Axm−1 ≤ (1 − ρ(M (E1 )−1 F1 ))F2 xm−1 .

(3.9)

Noting that A = E2 − F2 . Then by (3.9) we have ρ(M (E1 )−1 F1 )E2 xm−1 ≤ F2 xm−1 .

(3.10)

ρ(M (E1 )−1 F1 )x[m−1] ≤ M (E2 )−1 F2 xm−1 .

(3.11)

By M (E2 )−1 ≥ O, we have

By (3.2) of Lemma 3.2, we have ρ(M (E2 )−1 F2 ) ≥ ρ(M (E1 )−1 F1 ) ≥ 1. Moreover, if F2 < F1 , F2 = O and ρ(M (E1 )−1 F1 ) > 1, it is easy to check that the inequalities in (3.9)-(3.11) are strict. By (3.2) of Lemma 3.2, the first inequality is strict. The proof is completed. Remark 3.5. The first inequality in (1) of Theorem 3.4 can be strict, e.g., taking an order 3 dimension 2 tensor A = (aijk ) with a111 = 2, a222 = 2 and all other entries aijk = − 14 . Clearly, A is a strong M-tensor. Let E1 = D and E2 = E1 − L, where D = DI3 , L = LI3 and D, −L are the diagonal matrix and the strictly lower triangle matrix of M (A), respectively. It is easy to check that F1 = E1 − A ≥ F2 = E2 − A. Then ρ(M (E2 )−1 F2 ) = 0.3273 and ρ(M (E1 )−1 F1 ) = 0.3750. A Z-tensor A is said to be a non-strong M-tensor if A is an M-tensor but not a strong M-tensor. For a Z-tensor A, we have the following corollary. 7

Corollary 3.6. Let A ∈ R[m,n] be an irreducible Z-tensor, and let A = E1 − F1 = E2 − F2 be a weak regular splitting and a regular splitting, respectively, and F2 ≤ F1 , F2 = O. Then (1) ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ) < 1 if and only if A is a strong M-tensor. (2) ρ(M (E2 )−1 F2 ) = ρ(M (E1 )−1 F1 ) = 1 if and only if A is a non-strong M-tensor. (3) ρ(M (E2 )−1 F2 ) ≥ ρ(M (E1 )−1 F1 ) > 1 if and only if A is not an M-tensor. In particular, if F2 < F1 , F2 = O, the first inequality is strict. Proof. (1).” ⇐ ” If A is a strong M-tensor, by Lemma 2.8, we have that ρ(M (E2 )−1 F2 ) < 1 and ρ(M (E1 )−1 F1 ) < 1. By Theorem 3.4, we get ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ). ” ⇒ ” Since ρ(M (E1 )−1 F1 ) < 1, by Lemma 2.8, A is a strong M-tensor. [m,n] (2).” ⇐ ” Since A is a non-strong M-tensor, there exist s ∈ R+ and B ∈ R+ such that A = sIm − B with ρ(B) = s. Note that A is irreducible, so is B. By the strong PerronFrobenius theorem, B has a positive Perron eigenvector with Bxm−1 = ρ(B)x[m−1] = sx[m−1] . Then we have Axm−1 = 0. By A = E1 − F1 = E2 − F2 , we have M (E1 )−1 F1 xm−1 = x[m−1] and M (E2 )−1 F2 xm−1 = x[m−1] . By Lemma 3.3, we get ρ(M (E2 )−1 F2 ) = ρ(M (E1 )−1 F1 ) = 1. ” ⇒ ” Because ρ(M (E1 )−1 F1 ) = 1 and M (E1 )−1 F1 ≥ O , by the weak Perron-Frobenius theorem (see [2],[24] or [20]), there exists a nonnegative eigenvector x such that M (E1 )−1 F1 xm−1 = [m,n] x[m−1] . Then we get Axm−1 = 0. For any splitting A = sIm − B with s ∈ R+ and B ∈ R+ , it is easy to get Bxm−1 = sx[m−1] . Since A is irreducible, so is B. By the strong PerronFrobenius theorem, ρ(B) = s, and then A is a non-strong M-tensor. (3). ” ⇐ ” Since A is not an M-tensor, by the proof of (1) and (2), we know that ρ(M (E1 )−1 F1 ) > 1. By Theorem 3.4, we get ρ(M (E2 )−1 F2 ) ≥ ρ(M (E1 )−1 F1 ). ” ⇒ ” By (1) and (2), it is easy to check that A is not an M-tensor. Remark 3.7. It is easy to see that the conclusion (1) of Corollary 3.6 holds without the assumption of irreducibility of A. Immediately, by Corollary 3.6, we obtain the following Corollaries 3.8 and 3.9. Corollary 3.8. Let A ∈ R[m,n] be an irreducible Z-tensor, and let A = E1 − F1 = E2 − F2 be a weak regular splitting and a regular splitting, respectively, and F2 < F1 , F2 = O. Then one of the following statements holds: (1) ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ) < 1. (2) ρ(M (E2 )−1 F2 ) = ρ(M (E1 )−1 F1 ) = 1. (3) ρ(M (E2 )−1 F2 ) > ρ(M (E1 )−1 F1 ) > 1. Corollary 3.9. Let A ∈ R[m,n] be an irreducible Z-tensor. If A = Im − L − F, one of the following statements holds: 8

(1) ρ(M (Im − L)−1 F) ≤ ρ(L + F) < 1. (2) ρ(M (Im − L)−1 F) > ρ(L + F) > 1. (3) ρ(M (Im − L)−1 F) = ρ(L + F) = 1. Here L = LIm , and −L is the strictly lower triangle part of M (A). Corollary 3.10. Let A ∈ R[m,n] be an irreducible Z-tensor. If A = Im − L − F and 0 < ω1 < ω2 ≤ 1, then one of the following statements holds: (1) ρ(M (Im − ω2 L)−1 (ω2 F + (1 − ω2 )Im )) ≤ ρ(M (Im − ω1 L)−1 (ω1 F + (1 − ω1 )Im )) < 1. (2) ρ(M (Im − ω2 L)−1 (ω2 F + (1 − ω2 )Im )) > ρ(M (Im − ω1 L)−1 (ω1 F + (1 − ω1 )Im )) > 1. (3) ρ(M (Im − ω2 L)−1 (ω2 F + (1 − ω2 )Im )) = ρ(M (Im − ω1 L)−1 (ω1 F + (1 − ω1 )Im )) = 1. Proof. It is easy to check that A = ω1i (Im − ωi L) − ω1i (ωi F + (1 − ωi )Im ) (i = 1, 2), and 1 1 ω2 (ω2 F + (1 − ω2 )Im ) < ω1 (ω1 F + (1 − ω1 )Im ). By Corollary 3.8, the proof is completed.

4

Preconditioning techniques for solving multi-linear systems

As we introduced in Section 1, for solving the multi-linear system (1.1) with a strong Mtensor and a positive vector b in [12] the authors proposed tensor splitting iterative (TSI) method as follows: TSI Method 1) Given a positive vector b, a strong M-tensor A with its (weak) regular splitting A = E − F, a maximal number of iterations kmax , a machine precision ε and a positive initial vector x0 . Initialize k = 1. 2) while k < kmax 1

−1 [ m−1 ] , k = 1, 2, ... 3) xk = (M (E)−1 Fxm−1 k−1 + M (E) b)

4) If ||Axm−1 − b||2 < ε, break and output xk . k 5) k = k + 1, back to 2).

Also they gave the following convergence analysis: Theorem 4.1. [12] Let A be a strong M-tensor and b be a positive vector. For every (weak) regular splitting A = E − F, the sequence {xk } generated by the tensor splitting iterative (TSI) method converges to a unique positive solution of the multi-linear system (1.1) for any ≤ b. positive vector x0 > 0 with 0 < Axm−1 0 9

Now taking A=D−L−F

(4.1)

and E = D and E = D − L, we get the Jacob type method and the Gauss-Seidel type method in [12], where D = DIm , L = LIm , and D and −L are the diagonal part and strictly lower triangle part of M (A). In [12], Liu, Li and Vong demonstrated that the spectral radius ρ((M (E)−1 F) can serve as an approximate convergent rate for the tensor splitting iteration (1.3). Therefore, in this section, based on the idea of preconditioning linear systems we consider solving multi-linear system (1.1) by the preconditioned technique for improving the spectral radius of the iterative tensor. Since A is a strong M-tensor, aii···i > 0 for any i ∈ n(see [3] or [11]). Taking D = diag(a11 ···1 , · · · , ann···n ), by Lemma 2.8, D−1 A is also a strong M-tensor whose the diagonal entries are all 1, and D−1 b is also positive vector. Therefore, the equation D−1 Axm−1 = D−1 b is equivalent to the equation (1.1). Without loss of generality, we always assume that each diagonal entry of the tensor A in (1.1) is 1. Let P be a nonsingular matrix. Then we consider the following preconditioned multilinear system (4.2) which is equivalent to the equation (1.1) for a strong M-tensor A, P Axm−1 = P b.

(4.2)

The matrix P is called a preconditioner of solving (1.1). Let P A = Ep − Fp be a splitting of P A. Then the corresponding preconditioned tensor splitting iterative (PTSI) method is given as follows: PTSI Method 1) Given a positive vector b, a preconditioner P , a strong M-tensor A with its (weak) regular splitting P A = Ep − Fp , a maximal number of iterations kmax , a machine precision ε and a positive initial vector x0 . Initialize k = 1. 2) while k < kmax 3)

1

[ m−1 ] −1 . xk = (M (Ep )−1 Fp xm−1 k−1 + M (Ep ) P b)

(4.3)

4) If ||Axm−1 − b||2 < ε, break and output xk . k 5) k = k + 1, back to 2).

It is noted that when m = 2 the preconditioned multi-linear system reduces to the preconditioned linear system P Ax = P b, which has been studied extensively, and has powerful 10

methods to solve linear system by taking a suitable preconditioner P . For examples, when the coefficient matrix A is a nonsingular M-matrix, we may take the preconditioner P as P = I + Gα (see [5]), where α = (αi ) and αi is a parameter, i = 1, ..., n − 1, and ⎡ ⎤ 0 ··· 0 0 −α1 a12 ⎢ 0 ⎥ 0 −α2 a23 · · · 0 ⎢ ⎥ ⎢ . ⎥ . . . . ⎥. .. .. .. .. Gα = ⎢ ⎢ .. ⎥ ⎢ ⎥ ⎣ 0 0 0 · · · −αn−1 an−1,n ⎦ 0 0 0 ··· 0 If αi = 1, i = 1, ..., n − 1, the preconditioner P = I + G1 was given in [8](or see [11]). By an analogous technique, we consider the preconditioner Pα = I + Sα for solving the multi-linear system (4.2), where ⎡ ⎤ 0 ··· 0 0 −α1 a12···2 ⎢ 0 ⎥ 0 −α2 a23···3 · · · 0 ⎢ ⎥ ⎢ . ⎥ . . . . ⎥. .. .. .. .. Sα = ⎢ ⎢ .. ⎥ ⎢ ⎥ ⎣ 0 0 0 · · · −αn−1 an−1,n···n ⎦ 0 0 0 ··· 0 Remark 4.2. Let A be a strong M-tensor and b be a positive vector. Taking preconditioner Pα = I + Sα , then by Theorem 4.1 it is easy to know that for every (weak) regular splitting P A = Ep −Fp , the sequence {xk } generated by the PTSI method converges to a unique positive ≤ b. solution of the multi-linear system (1.1) for any positive vector x0 > 0 with 0 < Axm−1 0 That is to say, the PTSI method also converges when the TSI method converges. Let A = (I + Sα )A and

− L − F,

A = D

(4.4)

= DI

m , L = LI

m , and D,

−L

are the diagonal part, the strictly lower triangle part where D

respectively. Since A = Im − L − F, and L = LIm , −L is the strictly lower triangle of M (A), part of M (A), we have A = Im − L − Sα L − F − Sα F + Sα Im .

(4.5)

Then, we propose the following preconditioned Gauss-Seidel type method: 1

[ m−1 ] xk = (Tα xm−1 , k−1 + b)

where Tα is the iterative tensor of the method (4.6), Tα = M (Eα )−1 Fα , b = M (Eα )−1 (I + Sα )b, 11

(4.6)

and Eα = Im − L − Sα L, Fα = F + Sα F − Sα Im .

(4.7)

Remark 4.3. If α1 = · · · = αn−1 = 0, Tα = T0 , and the preconditioned Gauss-Seidel type iterative methods (4.6) reduces to the Gauss-Seidel type iterative method in [12].

we can get the corresponding preconditioned Jacob type method. Remark 4.4. If Ep = D, By Corollary 3.9, we know that the spectral radius of the iterative tensor for the Gauss-Seidel type method is equal or lesser than the corresponding one of Jacob type method. Therefore, we only consider the theoretical analysis of the preconditioned Gauss-Seidel type method for solving the equation (1.1).

5

The monotonicity of approximate convergence rate of the preconditioned Gauss-Seidel type method (4.6)

In this section, we discuss the relationship between the parameter α and ρ(Tα ). Firstly, we give the following lemma. Lemma 5.1. Let A ∈ R[m,n] and A = Im −L−F, F ≥ O, where L = LIm , −L is the strictly lower part of M (A). If A is a strong M-tensor, then for all αi ∈ [0, 1], i = 1, 2, ..., n − 1, (I + Sα )A is a strong M-tensor. Proof. Since

aji2 ···im =

aji2 ···im − αj ajj+1···j+1 aj+1i2 ···im , 1 ≤ j < n − 1, ani2 ···im , j = n,

(5.1)

aji2 ···im ≤ 0 for αi ∈ [0, 1], i.e., A is a Z-tensor. for (j, i2 , · · · , im ) = (j, j, · · · , j), we have

Taking (5.2) Eα = (I + Sα )(Im − L), Fα = (I + Sα )F, it is easy to check that A = Eα − Fα and Fα ≥ O. Since M (Eα )−1 Fα = (I − L)−1 (I + Sα )−1 (I + Sα )F = (I − L)−1 F ≥ O,

(5.3)

A = Eα − Fα is a weak regular splitting. Because A is a strong M-tensor and by (5.3)  = ρ((I − L)−1 F) < 1. By Lemma 2.8, A is a strong M-tensor.  −1 F) ρ(M (E) Since (I + Sα )b ≥ b > 0 with any αi ∈ [0, 1], i = 1, 2, ..., n − 1, by Lemmas 2.9 and 5.1, we obtain Theorem 5.2. For P = I + Sα with any αi ∈ [0, 1], i = 1, 2, ..., n − 1, the preconditioned multi-linear system (4.2) has the same unique positive solution as in (1.1). 12

Lemma 5.3. Let A ∈ R[m,n] be a strong M-tensor, and let A = E1 − F1 = E2 − F2 be two weak regular splittings with M (E2 )−1 ≥ M (E1 )−1 . If the Perron vector x of M (E2 )−1 F2 satisfies Axm−1 ≥ 0, then ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ). Proof. By M (E2 )−1 F2 xm−1 = ρ(M (E2 )−1 F2 )x[m−1] , we get F2 xm−1 = ρ(M (E2 )−1 F2 ))M (E2 )x[m−1] , which gives Axm−1 = (1 − ρ(M (E2 )−1 F2 ))M (E2 )x[m−1] .

(5.4)

Because A is a strong M-tensor and Axm−1 ≥ 0, by (4) of Lemma 2.8, we obtain M (E2 )x[m−1] ≥ 0.

(5.5)

By (5.4), we get E1 xm−1 = F1 xm−1 + (1 − ρ(M (E2 )−1 F2 ))M (E2 )xm−1 . Thus, x[m−1] = M (E1 )−1 F1 xm−1 + (1 − ρ(M (E2 )−1 F2 ))M (E1 )−1 M (E2 )x[m−1] , which together with M (E2 )−1 ≥ M (E1 )−1 ≥ O and (5.5) gives x[m−1] ≤ M (E1 )−1 F1 xm−1 + (1 − ρ(M (E2 )−1 F2 ))M (E2 )−1 M (E2 )x[m−1] .

(5.6)

By (5.6), we obtain M (E1 )−1 F1 xm−1 ≥ ρ(M (E2 )−1 F2 )x[m−1] .

(5.7)

By (3.2) of Lemma 3.2 and M (E1 )−1 F1 ≥ O, we obtain ρ(M (E2 )−1 F2 ) ≤ ρ(M (E1 )−1 F1 ). Note that A = Pα A = Eα − Fα , we consider the following splitting of A A = EPα − FPα ,

(5.8)

EPα = Pα−1 Eα = Pα−1 Im − L, FPα = Pα−1 Fα = Pα−1 Im − Im + F.

(5.9)

M (EPα )−1 FPα = (Pα−1 M (Eα ))−1 Pα−1 Fα = Tα ≥ O,

(5.10)

where

Thus, together with M (EPα )−1 = (Pα−1 M (Eα ))−1 = M (Eα )−1 Pα ≥ O gives that A = EPα − FPα is a weak regular splitting. Lemma 5.4. Let A ∈ R[m,n] be a strong M-tensor. If αi ∈ [0, 1], i = 1, 2, ..., n − 1, and ≥ 0. (ρα , xα ) is Perron eigenpair of Tα , Axm−1 α 13

Proof. If ρα = 0, there is nothing to do. Now we assume that ρα > 0. Since (ρα , xα ) is Perron eigenpair of Tα , by (5.10) [m−1]

ρα M (EPα )xα

= FPα xm−1 . α

(5.11)

2 )−1 (I − S ) = (I + S 2 (I − S 2 )−1 )(I − S ), we get Notice that Pα−1 = (I − Sα α α α α 2 2 −1 (I − Sα ) )(I − Sα )Im − L, EPα = (I + Sα

FPα = F − Sα Im +

2 Sα (I

−

2 −1 Sα ) (I

(5.12)

− Sα )Im .

(5.13)

Because F ≥ Sα Im , xα ≥ 0, and L ≥ O, by (5.11)-(5.13), we obtain [m−1]

ρα (I + Qα )(I − Sα )xα

[m−1]

≥ Qα (I − Sα )xα

,

(5.14)

2 (I − S 2 )−1 . Furthermore, where Qα = Sα α

ρα (I −

1 − ρα [m−1] Qα )(I − Sα )xα ≥ 0. ρα

(5.15)

2 )−1 = I + Q , it is easy to check that Q is a strict lower triangle matrix, which By (I − Sα α α 1−ρα together with 0 < ρα < 1 gives that I − ρα Qα is a nonsingular M-matrix. Thus, we have [m−1]

(I − Sα )xα Axm−1 = α

≥ 0. Thus, by (5.11),

1 − ρα 2 2 −1 FPα xm−1 = (F − Sα Im )xm−1 + Sα (I − Sα ) (I − Sα )Im xm−1 ≥ 0. (5.16) α α α ρα

The proof is completed. Let α1 = (α1,i ) and α2 = (α2,i ). Then we give the following theorem. Theorem 5.5. Let A ∈ R[m,n] be a strong M-tensor. If α1,i , α2,i ∈ [0, 1], i = 1, ..., n − 1 and α2 ≤ α1 , ρ(Tα1 ) ≤ ρ(Tα2 ). Proof. If ρ(Tα1 ) = 0, ρ(Tα2 ) ≥ 0 holds. We assume that ρ(Tα1 ) > 0. Since α1,i , α2,i ∈ [0, 1], i = 1, ..., n − 1 and α2 ≤ α1 , we have Pα2 ≤ Pα1 . Furthermore, Eα1 ≤ Eα2 . Because both M (Eα1 ) and M (Eα2 ) are nonsingular M-matrices, we have M (Eα1 )−1 ≥ M (Eα2 )−1 , and then M (EPα1 )−1 ≥ M (EP α2 )−1 . If (ρ(Tα1 ), xα1 ) is an eigenpair of Tα1 , by Lemma 5.4, we have Axm−1 α1 ≥ 0. Furthermore, by Lemma 5.3, we get ρ(Tα1 ) ≤ ρ(Tα2 ). By Theorem 5.5, immediately, we give the comparison result between the preconditioned Gauss-Seidel type iterative method (4.6) and the Gauss-Seidel type iterative method in [12]. Corollary 5.6. Let A ∈ R[m,n] and A = Im − L − F, F ≥ O, where L = LIm with −L being the strictly lower part of M (A). If A is a strong M-tensor, then for all αi ∈ [0, 1], i = 1, 2, ..., n − 1, ρ(Tα ) ≤ ρ(T0 ) < 1. Remark 5.7. From Theorems 5.5, we know that ρ(Tα ) have the smallest value as α = 1. 14

6

The preconditioned SOR type method

In [12], the SOR type method for solving (1.1) is proposed by taking E = (1.3):

1 υ (Im

1

−1 [ m−1 ] xk = (M (Im − υL)−1 ((1 − υ)Im + υF)xm−1 . k−1 + υM (Im − υL) b)

− υL) in

(6.1)

In this section, we consider the preconditioned SOR type method: 1

[ m−1 ] xk = (Rα (ω)xm−1 , k−1 + b)

(6.2)

where Rα (ω) is the iterative tensor of the method (6.2), i.e., Rα (ω) = (RE α (ω))−1 RF α (ω), b = M (RE α (ω))−1 (I + Sα )b, and RE α (ω) =

1 1 (Im − ω(L + Sα L)), RF α (ω) = ((1 − ω)Im + ω(F + Sα F − Sα Im )). (6.3) ω ω

By Corollary 3.10, we obtain Theorem 6.1. Let A ∈ R[m,n] be a strong M-tensor. If A = Im −L−F and 0 < ω1 < ω2 ≤ 1, ρ(Rα (ω2 )) ≤ ρ(Rα (ω1 )) < 1. If ω = 1, the preconditioned SOR type method (6.2) reduces to the Gauss-Seidel type method (4.6). Thus, by Theorem 6.1, we have Corollary 6.2. If A ∈ R[m,n] is a strong M-tensor, for any ω ∈ (0, 1], ρ(Tα ) ≤ ρ(Rα (ω)) < 1. Remark 6.3. For the preconditioned method, although theoretically we only consider the parameter αi ∈ [0, 1] in the preconditioner Pα , in numerical experiments we can search the optimal parameter when αi > 1(i = 1, 2, ..., n); see the next section. Remark 6.4. As shown in [12], the convergence rate of the tensor splitting iterative method is only a theoretical result, which can’t be computed. In this case, the spectral radius of the iterative tensor can be taken as an approximate convergent rate. Hence we apply the preconditioned technique to solve multi-linear systems in order to improve the spectral radii of the iterative tensors by tensor splitting methods, which can be an alternative algorithm for solving multi-linear systems. 15

7

Numerical examples

In this section, we will compare preconditioned tensor splitting methods with those without preconditioning. All tests will be done in MATLAB R2014a with the configuration: Intel(R) Core(TM)i7-2600 CPU 3.40GHz and 16.00G RAM. In the following examples, two aspects are given to check the efficiency of the proposed algorithms: the number of iteration steps (denoted by IT), the CPU time in seconds (denoted by CPU(s)). We set the maximum iterative number as 1000 and the stopping criteria is ε = 10−11 . In numerical experiments, we take Sα = αS1 . In [4], the authors proposed the SOR-like type method (also see [12] or [6]) to accelerate their Jacob type and Gauss-Seidel type method as follows: 1

xk+1 = M (E − Im )−1 ((F − Im )xm−1 + b)[ m−1 ] . k

(7.1)

We also use this method to accelerate the preconditioned Jacobi and Gauss-Seidel type methods for solving (1.1). All the optimized parameters are chosen by numerical experiments. The product Axm−1 defined in (1.2) is computed by transforming into the following matrix-vector product: Axm−1 = A(1) (x · · ⊗ x), (7.2) ⊗ ·  m−1

where ⊗ is the Kronecker product. The product BA with a matrix B and a tensor A is computed by (2.2). All the tested methods are stopped if ||Axm−1 − b||2 < ε.

(7.3)

Let P A = Ep − Fp . We use the same notations as in (4.4)-(4.7) and (6.1)-(6.3). Furthermore, we use simple notations to denote the proposed methods as given in Table 1. Table 1: Abbreviations of the preconditioned iterative methods and the corresponding splitting Ep . Abbreviations of the preconditioned methods Pα Jacob Pα GS Pα SOR Pα Jacobsorlike Pα GSsorlike

Ep Ep Ep Ep Ep Ep

=D

− L

=D

− ω L)

= ω1 (D

= D − Im

− L − Im =D

Remark 7.1. If α = 0, Pα = P0 = I and P0 Jacob, P0 GS, P0 SOR, P0 Jacobsorlike , P0 GSsorlike reduce to the corresponding methods without preconditioning in [12]. The following examples are taken from [12] and [23], respectively. In numerical experiments, we take the value of α from 0.01 to 2 in the interval of 0.01. For the SOR type or SOR-like methods, we search the optimal parameter from 0.01 to 2 in the interval of 0.01 to research the optimal parameter. We take initial value x0 = 0 and the right-hand side b = 1. 16

Example 7.2. Let B ∈ R[3,5] be a nonnegative tensor with bi1 i2 i3 = |tan(i1 + i2 + i3 )|. We can compute ρ(B) = 364.4895. Then A = 864.4895Im − B is a strong M-tensor. Example 7.3. Let B ∈ R[3,3] be a nonnegative tensor with bi1 i2 i3 = |sin(i1 + i2 + i3 )|. By [23], we know A = n2 I3 − B is a strong M-tensor. We report the numerical results in Tables 2 and 3, respectively. Besides, we plot the Figures 1-6 to show the relationship between the CPU time and the parameter α. The CPU times are taken as the base-10 logarithm in the corresponding figures. Remark 7.4. From Figures 1 -4 and Tables 2-3, we know that the Pα Jacob, Pα GS and Pα SOR methods converge faster than the corresponding ones without preconditioning for some suitable parameters. For Example 7.2, the optimal parameter α is taken in (1.4, 1.6) from Figures 1,3 and 5. For Example 7.3, the optimal parameter α is 2.

Table 2: Numerical results of the preconditioned tensor splitting methods for Example 7.2. α 0.00 0.01 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48 1.50 1.55 1.60 1.70 1.80 1.89 1.90 2.00

CPU(s) 0.026975 0.002533 0.001986 0.002049 0.001964 0.001951 0.001904 0.001801 0.001763 0.001720 0.001619 0.001633 0.001586 0.001457 0.001398 0.001312 0.001276 0.001189 0.001194 0.001171 0.001135 0.001188 0.001193 0.001190 0.001204 0.001199 0.001181 0.001307 0.001291 0.001362 0.001409 0.001447

Pα Jacob IT 24 24 23 23 22 22 22 21 20 20 19 19 18 17 16 15 15 14 14 14 13 14 14 14 14 14 14 15 15 16 16 17

CPU(s) 0.002077 0.001743 0.001967 0.001984 0.001951 0.001926 0.001800 0.001791 0.001706 0.001714 0.001613 0.001627 0.001576 0.001448 0.001405 0.001298 0.001181 0.001176 0.001181 0.001078 0.001125 0.001177 0.001182 0.001219 0.001186 0.001228 0.001258 0.001298 0.001320 0.001354 0.001391 0.001437

Pα GS IT 24 24 23 23 22 22 21 21 20 20 19 19 18 17 16 15 14 14 14 13 13 14 14 14 14 14 14 15 15 16 16 17

CPU(s) 0.001594 0.001363 0.001648 0.001518 0.001575 0.001444 0.001436 0.001377 0.001273 0.001292 0.001188 0.001110 0.001105 0.001016 0.000987 0.000933 0.000833 0.000836 0.000860 0.000911 0.000932 0.000927 0.000819 0.000824 0.000921 0.000937 0.000929 0.001021 0.001079 0.001209 0.001172 0.001209

17

Pα SOR IT 19 19 19 18 18 17 17 16 15 15 14 13 13 12 12 11 10 10 10 11 11 11 10 10 11 11 11 12 13 14 14 14

CPU(s) 0.003133 0.001694 0.001597 0.001524 0.001521 0.001428 0.001489 0.001351 0.001258 0.001279 0.001189 0.001221 0.001106 0.001015 0.001017 0.000910 0.000974 0.000917 0.000822 0.000906 0.000915 0.000910 0.000921 0.000921 0.000947 0.000992 0.001011 0.001092 0.001074 0.001180 0.001169 0.001263

Pα Jacobsorlike IT 20 20 19 18 18 17 17 16 15 15 14 14 13 12 12 11 11 11 10 11 11 11 11 11 11 12 12 13 13 14 14 15

CPU(s) 0.001665 0.001694 0.001606 0.001507 0.001519 0.001438 0.001485 0.001362 0.001261 0.001273 0.001184 0.001271 0.001108 0.001009 0.001039 0.000916 0.000905 0.000821 0.000907 0.000907 0.000916 0.000904 0.000915 0.000932 0.000989 0.000996 0.000995 0.001092 0.001091 0.001231 0.001232 0.001260

Pα GSsorlike IT 20 20 19 18 18 17 17 16 15 15 14 15 13 12 12 11 11 10 11 11 11 11 11 11 11 12 12 13 13 14 14 15

Table 3: Numerical results of the preconditioned tensor splitting methods for Example 7.3. α 0.00 0.01 0.07 0.10 0.19 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.15 1.19 1.20 1.30 1.40 1.49 1.60 1.70 1.79 1.80 1.86 1.87 1.88 1.89 1.90 1.91 1.92 1.93 1.94 1.95 1.98 1.99 2.00

CPU(s) 0.005181 0.004756 0.004755 0.004737 0.004753 0.004674 0.004712 0.004573 0.004514 0.004468 0.004449 0.004366 0.004398 0.004297 0.004313 0.004338 0.004270 0.004311 0.004244 0.004145 0.004102 0.004123 0.004152 0.004055 0.003999 0.004031 0.004043 0.004047 0.004172 0.004044 0.004003 0.004298 0.004138 0.003964 0.003951 0.003999 0.004073

Pα Jacob IT 54 54 53 53 53 52 52 51 51 50 50 49 49 48 48 48 48 47 47 46 46 45 45 45 45 45 45 45 45 45 45 45 45 44 44 44 44

CPU(s) 0.004145 0.004191 0.004198 0.004261 0.004265 0.004275 0.004097 0.004054 0.004007 0.004013 0.003939 0.003974 0.003839 0.003836 0.003787 0.003739 0.003724 0.003741 0.003652 0.003635 0.003660 0.003604 0.003574 0.003583 0.003538 0.003577 0.003485 0.003455 0.003584 0.003472 0.003683 0.003476 0.003480 0.003476 0.003507 0.003440 0.003493

Pα GS IT 47 47 47 47 47 46 46 45 45 44 44 43 43 43 42 42 42 42 41 41 41 40 40 40 40 40 39 39 39 39 39 39 39 39 39 39 39

CPU(s) 0.002284 0.002306 0.002214 0.002277 0.002261 0.002118 0.002129 0.002099 0.002017 0.002030 0.002019 0.002016 0.001949 0.001931 0.001942 0.001945 0.001950 0.001951 0.001869 0.001924 0.001850 0.001866 0.001831 0.001842 0.001761 0.001793 0.001764 0.001811 0.001757 0.001773 0.001763 0.001764 0.001769 0.001782 0.001768 0.001753 0.001758

18

Pα SOR IT 27 27 26 26 26 25 25 25 24 24 24 24 23 23 23 23 23 23 22 22 22 22 21 21 21 21 21 21 21 21 21 21 21 21 21 21 21

CPU(s) 0.003014 0.003102 0.003115 0.003186 0.003051 0.003139 0.002980 0.002930 0.003045 0.002945 0.003151 0.002913 0.002758 0.002756 0.002747 0.002720 0.002669 0.002720 0.002652 0.002581 0.002558 0.002584 0.002469 0.002491 0.002482 0.002539 0.002497 0.002473 0.002499 0.002477 0.002458 0.002531 0.002481 0.002477 0.002527 0.002469 0.002487

Pα Jacobsorlike IT 35 35 35 35 34 34 33 33 33 32 32 32 31 31 31 30 30 30 30 29 29 29 28 28 28 28 28 28 28 28 28 28 28 28 28 28 28

CPU(s) 0.002380 0.002471 0.002462 0.002505 0.002546 0.002383 0.002385 0.002396 0.002295 0.002360 0.002283 0.002326 0.002207 0.002190 0.002207 0.002257 0.002462 0.002148 0.002299 0.002130 0.002019 0.002027 0.002019 0.002122 0.002090 0.002019 0.002032 0.001996 0.001935 0.001940 0.001944 0.001937 0.001942 0.001936 0.001922 0.001927 0.001984

Pα GSsorlike IT 28 28 28 28 28 27 27 27 26 26 26 25 25 25 25 25 25 24 24 24 23 23 23 23 23 23 23 22 22 22 22 22 22 22 22 22 22

Table 4: The corresponding values of the optimal parameters (a) The corresponding values of the optimal parameters for the proposed methods in Table 2. α 0.00 0.01 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48 1.50 1.55 1.60 1.70 1.80 1.89 1.90 2.00

Pα SOR ω 1.14 1.16 1.18 1.16 1.15 1.18 1.15 1.20 1.19 1.21 1.18 1.20 1.19 1.17 1.17 1.13 1.13 1.12 1.09 1.09 1.10 1.14 1.13 1.12 1.16 1.15 1.13 1.17 1.18 1.14 1.23 1.22

Pα Jacobsorlike 0.04 0.04 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.04 0.04 0.04 0.04 0.03 0.04 0.04 0.04 0.04 0.04 0.04 0.05 0.05 0.05 0.05 0.05 0.06

(b) The corresponding values of the optimal parameters for the proposed methods in Table 3.

Pα GSsorlike 0.04 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.06 0.05 0.05 0.05 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.05 0.05 0.03 0.05 0.05 0.05 0.03

α 0.00 0.01 0.07 0.10 0.19 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.15 1.19 1.20 1.30 1.40 1.49 1.60 1.70 1.79 1.80 1.86 1.87 1.88 1.89 1.90 1.91 1.92 1.93 1.94 1.95 1.98 1.99 2.00

Pα SOR ω 1.45 1.45 1.45 1.45 1.46 1.48 1.47 1.48 1.45 1.45 1.46 1.46 1.46 1.44 1.43 1.43 1.44 1.44 1.44 1.46 1.42 1.43 1.44 1.44 1.44 1.44 1.44 1.44 1.44 1.44 1.45 1.43 1.45 1.44 1.45 1.42 1.43

Pα Jacobsorlike 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07

Pα GSsorlike 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07

−2.22

−1.5 Jacob PαJacob

Jacob PαJacob

−2.24

−2.26

−2.28 −2

CPU(s)

CPU(s)

−2.3

−2.32

−2.34 −2.5 −2.36

−2.38

−2.4

−3

0

0.2

0.4

0.6

0.8

1 α

1.2

1.4

1.6

1.8

−2.42

2

Figure 1: The relationship between CPU(s) and α for the Jacob type method applied to Example 7.2.

0

0.2

0.4

0.6

0.8

1 α

1.2

1.4

1.6

1.8

2

Figure 2: The relationship between CPU(s) and α for the Jacob type method applied to Example 7.3.

19

−2.5

−2.3 GS P GS

GS P GS

α

−2.55

α

−2.32

−2.6

−2.34

−2.65 −2.36

CPU(s)

CPU(s)

−2.7

−2.75

−2.38

−2.4

−2.8 −2.42 −2.85 −2.44

−2.9

−2.46

−2.95

−3

0

0.2

0.4

0.6

0.8

1 α

1.2

1.4

1.6

1.8

−2.48

2

Figure 3: The relationship between CPU(s) and α for the Gauss-Seidel type method applied to Example 7.2.

0

0.2

0.4

0.6

0.8

1 α

1.2

1.4

1.6

−2.6 SOR PαSOR

SOR P SOR α

−2.75

−2.62

−2.8

−2.64

−2.85

−2.66 CPU(s)

CPU(s)

2

Figure 4: The relationship between CPU(s) and α for the Gauss-Seidel type method applied to Example 7.3.

−2.7

−2.9

−2.68

−2.95

−2.7

−3

−2.72

−3.05

−2.74

−3.1

1.8

0

0.2

0.4

0.6

0.8

1 α

1.2

1.4

1.6

1.8

−2.76

2

Figure 5: The relationship between CPU(s) and α for the SOR type method applied to Example 7.2.

0

0.2

0.4

0.6

0.8

1 α

1.2

1.4

1.6

1.8

2

Figure 6: The relationship between CPU(s) and α for the SOR type method applied to Example 7.3.

20

It is expensive to compute matrix-tensor products. Therefore, in the following Example 7.5 arising from an ordinary differential equation with Dirichlet’s boundary conditions, we consider sparse tensor structure. We take the value of α from 0.1 to 3 in the interval of 0.1. For the SOR type or SOR-like methods, we search the optimal parameter from 0.1 to 4 with the interval of 0.1 to research the optimal parameter. We choose initial value x0 = 1. Example 7.5. Let A ∈ R[3,n] and b ∈ Rn×n with ⎧ a111 = annn = 1, ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ aiii = 2, i = 2, 3, · · · , n − 1, aii−1i = −1/2, i = 2, 3, · · · , n − 1, ⎪ ⎪ ⎪ a ii−1i−1 = −1/2, i = 2, 3, · · · , n − 1, ⎪ ⎪ ⎩ aii+1i+1 = −1/2, i = 2, 3, · · · , n − 1. and

⎧ 2 ⎪ ⎨ b1 = c0 , a bi = (n−1) 2 , i = 2, 3, · · · , n − 1, ⎪ ⎩ 2 bn = c 1 .

Taking c0 = 1/2, c1 = 1/3 and a = 2, we solve the Ax2 = b by the proposed splitting methods, and report the numerical results in Table 5(a) and the corresponding parameters in 5(b). Remark 7.6. In Example 7.5 it is known that the preconditioned methods perform better than those without preconditioning.

21

Table 5: The preconditioned methods for Example 7.5. (a) Numerical results for the preconditioned tensor splitting methods applied to Example 7.5. n 10(α = 0) 10(preconditioning) 20(α = 0) 20(preconditioning) 50(α = 0) 50(preconditioning) 90(α = 0) 90(preconditioning) 100(α = 0) 100(preconditioning) 150(α = 0) 150(preconditioning) 200(α = 0) 200(preconditioning) 300(α = 0) 300(preconditioning)

CPU(s) 0.005591 0.002603 0.006466 0.003628 0.008497 0.004744 0.017382 0.006175 0.018719 0.006728 0.017970 0.009941 0.024957 0.014201 0.056784 0.032174

Pα Jacob IT 74 39 85 46 91 50 92 51 92 51 93 52 94 53 95 54

CPU(s) 0.003908 0.001598 0.004816 0.001946 0.005897 0.002582 0.007994 0.003851 0.008996 0.004124 0.012885 0.005558 0.018254 0.008588 0.040287 0.018222

Pα GS IT 55 24 62 28 65 29 66 30 66 30 67 30 67 31 68 31

CPU(s) 0.002535 0.001395 0.003008 0.001561 0.003730 0.002295 0.005293 0.003019 0.005674 0.003309 0.007892 0.005297 0.011434 0.006753 0.026092 0.014680

Pα SOR IT 35 20 39 22 41 23 41 23 42 24 42 24 42 24 43 25

CPU(s) 0.004457 0.002539 0.005696 0.003072 0.010665 0.004254 0.015412 0.005900 0.016490 0.006458 0.015735 0.009183 0.022289 0.014069 0.051869 0.030766

Pα Jacobsorlike IT 65 37 75 44 80 48 82 48 82 49 83 49 83 49 84 50

CPU(s) 0.003384 0.001442 0.004011 0.001713 0.004867 0.002334 0.007233 0.003362 0.007549 0.003457 0.010523 0.005038 0.014946 0.007456 0.033526 0.015874

Pα GSsorlike IT 46 21 52 24 55 25 55 26 56 26 56 26 56 27 57 27

(b) The corresponding values of the optimal parameters for the proposed methods in Table 5(a). n 10 20 50 90 100 150 200 300

P0 SOR ν 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3

P0 Jacobsorlike 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

P0 GSsorlike 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

Pα Jacob α 1.6 1.7 1.8 1.8 1.8 1.8 1.8 1.7

22

Pα GS α 2.0 2.0 2.1 2.0 2.0 2.0 2.0 2.0

Pα SOR (ω, α) (1.2, 1.5) (1.2, 1.6) (1.3, 1.2) (1.3, 1.2) (1.3, 1.1) (1.3, 1.2) (1.3, 1.2) (1.3, 1.1)

Pα Jacobsorlike ( , α) (0.1, 1.4) (0.1, 1.4) (0.1, 1.5) (0.1, 1.5) (0.1, 1.5) (0.1, 1.5) (0.1, 1.5) (0.1, 1.5)

Pα GSsorlike ( , α) (0.1, 1.7) (0.1, 1.7) (0.1, 1.8) (0.1, 1.7) (0.1, 1.7) (0.1, 1.7) (0.1, 1.7) (0.1, 1.7)

8

Concluding Remarks

In this paper, we compare spectral radii of iterative tensors for two tensor splitting methods, and propose a preconditioned tensor splitting iterative method for solving multi-linear system (1.1). The theoretical analysis for the preconditioned iterative method is also given. A type of preconditioners is proposed when the system tensor is a strong M-tensor. We theoretically show that this preconditioned technique can improve approximate rate of convergence. Numerical examples are given to show the efficiency of the proposed methods. Furthermore, we will study the computable rate of convergence for the tensor splitting method, and then develop more efficient preconditioned iterative methods for solving the multi-linear system (1.1).

9

Acknowledgement

The authors would like to thank the editor and reviewers for their kind comments. The first author was supported by National Natural Science Foundation of China (Grant Nos. 11671185, 11771159), and Major Project (Grant No.2016KZDXM025) and Innovation Team Project (Grant No. 2015KCXTD007) of Guangdong Provincial Universities, the third author was supported by University of Macau (Grant No. MYRG2017-00098-FST) and the Macao Science and Technology Development Fund (050/2017/A).

References [1] A. Berman, J.R. Plemmons, Nonnegative matrices in the Mathematical Science, SIAM, Philadelphia. 1994. [2] K. C. Chang, K. Pearson, T. Zhang, Perron-Frobenius theorem for nonnegative tensors, Commun. Math. Sci., 6(2008): 507-520. [3] W. Ding, L. Qi, Y. Wei, M-tensors and nonsingular M-tensors, Linear Algebra and Its Applications, 439(2013): 3264-3278. [4] W. Ding, Y. Wei, Solving multi-linear system with M-tensors, Journal of Scientific Computing, 68(2)(2016): 689-715. [5] A. D. Gunawardena, S. K. Jain, L. Snyder, Modified iterative methods for consistent linear systems, Linear Algebra and Its Applications, 154(1991): 123-143. [6] G. H. Golub, C.F. Van Loan, Matrix computations, 4th edition, Johns Hopkins University Press, Baltimore, 2013. 23

[7] A. Hadjidimos, D. Noutsos, M. Tzoumas, More on modifications and improvements of classical iterative schemes for M-matrices. Linear Algebra and Its Applications, 364(2003): 253-279. [8] T. Kohno , H. Kotakemori, H. Niki, et al, Improving the modified Gauss-Seidel method for Z-matrices, Linear Algebra and its Applications, 267(1997): 113-123. [9] T.G. Kolda, B.W. Bader, Tensor decompositions and applications, SIAM Review, 51(2009): 455-500. [10] T.G. Kolda, Multilinear operators for higher-order decompositions, Technical report SAND2006-2081, Sandia National Laboratories, Albuquerque, NM and Livermore, CA. [11] W. Li, W. Sun, Modified Gauss-Seidel type methods and Jacobi type methods for Zmatrices. Linear Algebra and Its Applications, 317(1-3)(2000): 227-240. [12] D. D. Liu, W. Li, S. W. Vong, The tensor splitting with application to solve multi-linear systems, Journal of Computational and Applied Mathematics, 330 (2018), 75-94. [13] W. Liu, W. Li, On the inverse of a tensor, Linear Algebra and its Applications, 495(2016): 199-205. [14] D. H. Li, S. Xie , H. R. Xu, Splitting methods for tensor equations, Numerical Linear Algebra with Applications, DOI: 10.1002/nla.2102, 2017, online. [15] L.H. Lim, Singular values and eigenvalues of tensors: a variational approach, in: Proceedings of the IEEE International Workshop on ComputationalAdvances in Multi-Sensor Adaptive Processing, CAMSAP 05, vol.1, IEEE Computer Society Press, Piscataway, NJ. 2005, pp.129-132. [16] J. P. Milaszewicz, Improving jacobi and gauss-seidel iterations, Linear Algebra and Its Applications, 93(1987): 161-170. [17] R.S. Varga, Matrix iterative analysis, Springer Science and Business Media, Berlin Heidelberg, 2009. [18] K. Pearson, Essentially positive tensors, Int. J. Algebra, 4(2010): 421-427. [19] L. Qi, Eigenvalues of a real supersymmetric tensor, Journal of Symbolic Computation, 40(2005): 1302-1324. [20] L. Qi, Z. Luo, Tensor analysis: Spectral theory and special tensors. Society for Industrial and Applied Mathematics, 2017. [21] A.E. Raftery, A model for high-order Markov chains, Journal of the Royal Statistical Society. Series B (Methodological), 1985: 528-539. 24

[22] J. Shao, L. You, On some properties of three different types of triangular blocked tensors, Linear Algebra and its Applications, 511(2016): 110-140. [23] Z. J. Xie, X. Q. Jin, Y. M. Wei, Tensor methods for solving symmetric M-tensor systems, Journal of Scientific Computing, 2017: 1-14. [24] Y. Yang, Q. Yang, Further results for Perron-Frobenius theorem for nonnegative tensors, SIAM Journal on Matrix Analysis and Applications, 31(2010): 2517-2530. [25] Y. Yang, Q. Yang, A study on eigenvalues of higher-order tensors and related polynomial optimization problems, Science Press, Beijing, 2015. [26] L. Zhang, L. Qi, G. Zhou, M-tensors and some applications, SIAM Journal on Matrix Analysis and Applications, 35(2014): 437-452.

25