Inexact Newton Method for M-Tensor Equations
aa r X i v : . [ m a t h . O C ] J u l Inexact Newton Method for M-Tensor Equations ∗ Dong-Hui Li † Hong-Bo Guan ‡ Jie-Feng Xu § August 7, 2020
Abstract
We first investigate properties of M-tensor equations. In particular, we show that ifthe constant term of the equation is nonnegative, then finding a nonnegative solution ofthe equation can be done by finding a positive solution of a lower dimensional M-tensorequation. We then propose an inexact Newton method to find a positive solution tothe lower dimensional equation and establish its global convergence. We also show thatthe convergence rate of the method is quadratic. At last, we do numerical experimentsto test the proposed Newton method. The results show that the proposed Newtonmethod has a very good numerical performance.
Keywords
M-tensor equation, inexact Newton method, global convergence, quadraticconvergence
AMS
Tensor equation is a special system of nonlinear equations. It is also called multilinearequation. Tensor equation can be expressed as F ( x ) = A x m − − b = 0 , (1.1)where x, b ∈ R n and A is an m th-order n -dimensional tensor that takes the form A = ( a i i ...i m ) , a i i ...i m ∈ R , ≤ i , i , · · · , i m ≤ n, and A x m − ∈ R n with elements( A x m − ) i = X i ,...,i m a ii ...i m x i · · · x i m , i = 1 , , . . . , n. ∗ Supported by the NSF of China grant 11771157 and the NSF of Guangdong Province grantNo.2020B1515310013. † School of Mathematical Sciences, South China Normal University, Guangzhou, 510631, China, [email protected] ‡ School of Mathematical Sciences, South China Normal University, Guangzhou, 510631, China, [email protected]. School of Mathematics, Physics and Energy Engineering, Hunan Institute of Tech-nology, Hengyang, 421002, China. § School of Mathematical Sciences, South China Normal University, Guangzhou, 510631, China,[email protected]. A x m will denote the homogenous polynomial of degree m , i.e., A x m = x T A x m − = X i ,...,i m a i ...i m x i · · · x i m . For convenience of presentation, we introduce some concepts and notations, which willbe used throughout the paper. We denote the set of all m th-order n -dimensional tensors by T ( m, n ). We first introduce the concepts of Z-matrix and M-matrix. Definition 1.1. [3]
A matrix A is called a Z-matrix if all its off-diagonal entries are non-positive. It is apparent that a Z-matrix A can be written as A = sI − B, where B is a nonnegative matrix ( B ≥
0) and s >
0; When s ≥ ρ ( B ), we call A is anM-matrix; And further when s > ρ ( B ), we call A as a nonsingular M-matrix.The concept of M-tensor is an extension of the definition of M-matrix. Now we introducethe definition of M-tensor and other structure tensors that will be involved in this paper. Definition 1.2. [6, 7, 8, 19, 25, 26, 27, 36] Let
A ∈ T ( m, n ). • A is called a non-negative tensor, denoted by A ≥
0, if all its elements are non-negative,i.e., a i i ...i m ≥ ∀ i , . . . , i m ∈ [ n ], where [ n ] = { , · · · , n } . • A is called a symmetric tensor, if its elements a i i ...i m are invariant under any per-mutation of their indices. In particular, for every index i ∈ [ n ], if an (m-1)th ordern-dimensional square tensor A i := ( a ii ...i m ) ≤ i ,...,i m ≤ n is symmetric, then A is calledsemi-symmetric tensor with respect to the indices { i , . . . , i m } . The set of all m th-order n -dimensional symmetric tensors is denoted by ST ( m, n ). • A is called the identity tensor, denoted by I , if its diagonal elements are all ones andother elements are zeros, i.e., all a i i ...i m = 0 except a ii...i = 1, ∀ i ∈ [ n ]. • If a real number λ and a nonzero real vector x ∈ R n satisfy A x m − = λx [ m − , then λ is called an H-eigenvalue of A and x is called an H-eigenvector of A associatedwith λ . • A is called an M-tensor, if it can be written as A = s I − B , B ≥ , s ≥ ρ ( B ) , (1.2)where ρ ( B ) is the spectral radius of tensor B , that is ρ ( B ) = max {| λ | : λ is an eigenvalue of B} . If s > ρ ( B ), then A is called a strong or nonsingular M-tensor.2 A is called a lower triangular tensor, if its possibly nonzero elements are a i i ...i m with i = 1 , , . . . , n and i , . . . , i m ≤ i and all other elements of A are zeros. A is calleda strictly lower triangular tensor, if its possibly nonzero elements are a i i ...i m with i = 1 , , . . . , n and i , . . . , i m < i and all other elements of A are zeros. • A is called reducible it there is an index set I ⊂ [ n ] such that the elements of A satisfy a i i ...i m = 0 , ∀ i ∈ I, ∀ i , . . . , i m / ∈ I. If A is not reducible, then we call A irreducible.In the case A ∈ ST ( m, n ), the derivative of the homogeneous polynomial A x m can beexpressed as ∇ ( A x m ) = m A x m − .In the definition of reducible tensor, the index set I ⊂ [ n ] can be arbitrary. In ourpaper, we will need some special reducible tensor where the index set I is contained in somespecified set. For the sake of convenience, we make a slight extension to the definition ofreducible tensors. Definition 1.3.
Tensor
A ∈ T ( m, n ) is called reducible respect to I ⊂ [ n ] if its elementssatisfies a i i ...i m = 0 , ∀ i ∈ I, ∀ i , . . . , i m / ∈ I. It is easy to see that tensor
A ∈ T ( m, n ) is reducible if and only if there is an index I ⊂ [ n ] such that it is reducible respect to I .We call index pair ( I, I c ) a partition to [ n ] if I, I c ⊂ [ n ] and I ∪ I n = [ n ].For x, y ∈ R n and α ∈ R , the notations x ◦ y and x [ α ] are vectors in R n defined by x ◦ y = ( x y , · · · , x n y n ) T and x [ α ] = ( x α , . . . , x αn ) T respectively.We use R n + and R n ++ to denote the sets of all nonnegative vectors and positive vectorsin R n . That is, R n + = { x ∈ R n | x ≥ } and R n ++ = { x ∈ R n | x > } . If A is an M-tensor, we call the tensor equation an M-tensor equation and abbreviate itas M-Teq.The following theorem comes from [3, 8, 11]. Theorem 1.4.
Let
A ∈ ST ( m, n ) . • ([8]) If A is a strong M-tensor and b ∈ R n ++ , then the M-Teq (1.1) has a unique positivesolution. • ([11]) If A is a strong M-tensor and b ∈ R n + , then the M-Teq (1.1) has a nonnegativesolution. • ([3]) For a Z-matrix A ∈ R n × n , the following statements are equivalent. A is a nonsingular M-matrix. (ii) Av ∈ R n ++ for some vector v ∈ R n ++ . (iii) All the principal minors of A are positive.
Tensor equation is also called multilinear equation. It appears in many practical fieldsincluding data mining and numerical partial differential equations [5, 8, 9, 10, 14, 15, 16, 32].The study in numerical methods for solving tensor equations has begun only a few yearsago. Most existing methods focus on solving the M-Teq under the restriction b ∈ R n ++ or b ∈ R n + . Such as the iterative methods in [8], the homotopy method in [12], the tensorsplitting method in [20], the Newton-type method in [13], the continuous time neural networkmethod in [28], the preconditioned tensor splitting method in [22], the preconditioned SORmethod in [21], the preconditioned Jacobi type method in [35], the nonnegativity preservingalgorithm in [2]. There are also a few methods that can solve M-Teq (1.1) without restriction b ∈ R n ++ or that A is an M tensor. Those methods include the splitting method by Li, Guanand Wang [18], and Li, Xie and Xu [15], the alternating projection method by Li, Dai andGao [17], the alternating iterative methods by Liang, Zheng and Zhao [23] etc.. Relatedworks can also be found in [4, 5, 16, 24, 29, 30, 31, 32, 33, 34].Newton’s method is a well-known efficient method for solving nonlinear equations. Anattractive property of the method is its quadratic convergence rate. However, in many cases,the standard Newton method may fail to work or loss its quadratic convergence propertywhen applied to solve tensor equation (1.1). We refer to [18] for details.Recently, He, Ling, Qi and Zhou [13] proposed a Newton type method for solving theM-Teq (1.1) with b ∈ R n ++ . Unlike the standard Newton method for solving nonlinearequations, by utilizing the special structure of the equation (1.1), the authors transformedthe equation into an equivalent form through a variable transformation y = x [ m ] . Startingfrom some positive initial point, the method generates a sequence of positive iterates. Anattractive property of the method is that the Jacobian matrices of the equation at the iteratesare nonsingular. As a result, the method is well defined and retains the global and quadraticconvergence. The reported numerical results in [13] confirmed the quadratic convergenceproperty of that method.It should be pointed out that the positivity of b plays an important role in the Newtonmethod by He, Ling, Qi and Zhou [13]. It is not known if the method in [13] is still welldefined and reserves quadratic convergence property if there is some i satisfying b i = 0. Thepurpose of this paper is to develop a Newton method to find the a nonnegative solutionof the equation (1.1) with b ∈ R n + . Our idea is similar to but different from that of themethod in [13]. Specifically, we will reformulate the equation via the variable transformation y = x [ m − . Such an idea comes from the following observation. Consider a vary specialtensor equation Ax [ m − − b = 0 , corresponding to the tensor equation (1.1) where the only nonzero elements of A are a ij...j = a ij , i, j = 1 , , . . . , n . For that special equation, the tensor equation is equivalent to thesystem of linear equation Ay − b = 0 with y = x [ m − . As a result, the corresponding Newtonmethod terminates at a solution of the equation within one iteration. Another differencebetween our method and the method in [13] is that we will consider the equation (1.1) with b ∈ R n + . The case where b has zero elements cause the problem be much more difficult.Existing techniques that deals with equation (1.1) with b ∈ R n ++ are no longer available.To overcome that difficult, we will propose a criterion that can identify the zero elementsin a nonnegative solution of the M-tensor equation. From computational view point, the4riterion is easy to implement. By the use of that criterion, we can get a nonnegative solutionof the M-tensor equation (1.1) by finding a positive solution to a lower dimensional M-tensorequation with nonnegative constant term.Based on that criterion, we propose a Newton method for finding a positive solution ofthe M-Teq with b ∈ R n + and establish its global and quadratic convergence.The remainder of the paper is organized as follows. In the next section, we investigatesome nice properties of the M-tensor equation (1.1). In particular, we propose a criterion todistinguish zero and nonzero elements of a nonnegative solution of the equation. In Section3, we propose a Newton method to get a positive solution to the M-Teq (1.1) with b ∈ R n ++ and establish its global and quadratic convergence. In Section 4, we extend the methodproposed in Section 3 to the M-Teq (1.1) with b ∈ R n + and show its global and quadraticconvergence. At last, we do numerical experiments to test the proposed method in SectionSection 5. Throughout this section, we suppose that tensor
A ∈ T ( m, n ) is a strong M-tensor.The following lemma was proved by Li, Guan and Wang [18]. Lemma 2.1. If A is a strong M-tensor, and the feasible set S defined by S △ = { x ∈ R n + | F ( x ) = A x m − − b ≤ } is not empty, then S has a largest element that is the largest nonnegative solution to theM-tensor equation F ( x ) = A x m − − b = 0 . As an application of the last lemma, we have the following proposition.
Proposition 2.2.
Let A be a strong M-tensor and b (1) , b (2) ∈ R n satisfy b (2) ≥ b (1) . Supposethat the M-tensor equation A x m − − b (1) = 0 (2.1) has a nonnegative solution x (1) . Then the M-tensor equation A x m − − b (2) = 0 (2.2) has a nonnegative solution x (2) satisfying x (2) ≥ x (1) . In particular, if b (1) > , then theunique positive solution ¯ x (1) of (2.1) and the unique positive solution ¯ x (2) of (2.2) satisfies ¯ x (2) ≥ ¯ x (1) .Proof. Define S △ = { x ∈ R n + | A x m − − b (1) ≤ } and S △ = { x ∈ R n + | A x m − − b (2) ≤ } . Since b (1) ≤ b (2) , we obviously have S ⊆ S .
5y the assumption that (2.1) has a nonnegative solution, we claim from Lemma 2.1 that theset S is nonempty and has a largest element ¯ x (1) that is a solution to the equation (2.1).Consequently, the set S is nonempty and has a largest element x (2) that is a solution tothe equation (2.2). It is clear that x (2) ≥ ¯ x (1) ≥ x (1) . If b (1) >
0, then the unique positive solution ¯ x (1) is the largest element of S and ¯ x (2) is thelargest element of S . As a result, we have ¯ x (2) ≥ ¯ x (1) . The proof is complete. Theorem 2.3.
Suppose that A is a strong M-tensor. Then the following statements aretrue. (i) The tensor equation A x m − = 0 (2.3) has a unique solution x = 0 . (ii) If − b ∈ R n + \{ } , then the M-Teq (1.1) has no nonnegative solutions. (iii) The following relationship holds x ◦ A x m − = 0 ⇐⇒ x = 0 . (iv) It holds that lim k x k→∞ kA x m − k = + ∞ . (v) For any b ∈ R n , the solution set of the M-tensor equation (1.1) , if not empty, isbounded.Proof. Conclusion (i) is trivial because zero is not an eigenvalue of any strong M-tensor.(ii) Suppose for some b ≤ b = 0, the M-Teq (1.1) has a nonnegative solution ¯ x = 0.Clearly, ¯ x ≥
0. Denote I = { i : ¯ x > } . Let D be a diagonal tensor whose diagonals are d i ··· i = − b i ¯ x − ( m − i , ∀ i ∈ I , and d i ··· i = 0, ∀ i / ∈ I . Let ¯ A = A + D . It is obvious that ¯ A is astrong M-tensor. Clearly, ¯ A I is a strong M-tensor too. However, it holds that ¯ A I ¯ x m − I = 0,which yields a contradiction.Conclusion (iii) follows from (i) directly because any principal subtensor of a strongM-tensor is a strong M-tensor.(iv) Suppose on the contrary that there is some sequence { x k } satisfyinglim k →∞ k x k k = + ∞ such that the sequence {kA x m − k k} is bounded. Then we havelim k →∞ kA x m − k kk x k k m − = 0 . Let y k = x k / k x k k and ¯ y be an accumulation point of the sequence { y k } . It is easy to seethat ¯ y = 0 but A ¯ y m − = 0, which contradicts with (i).The conclusion (v) is a direct corollary of the conclusion (iv).6he latter part of this section focuses on the M-Teq (1.1) with b ∈ R n + . We denote I + ( b ) = { i : b i > } , and I ( b ) = { i : b i = 0 } . We first show the following theorem.
Theorem 2.4.
Suppose that A is irreducible and is a strong M-tensor. Then every non-negative solution of the M-Teq (1.1) with b ∈ R n + must be positive.Proof. Suppose on the contrary that the M-Teq (1.1) with b ∈ R n + has a nonnegative solution¯ x satisfying I = { i | ¯ x i = 0 } 6 = ∅ . We have for any i ∈ I ,0 = X i ,...,i m a ii ...i m ¯ x i · · · ¯ x i m − b i = X { i ,...,i m }⊆ I c a ii ...i m ¯ x i · · · ¯ x i m − b i ≤ . Since ¯ x j > ∀ j ∈ I c , the last inequality yields b i = 0 and a ii ...i m = 0 , ∀ i , . . . , i m I. It shows that tensor A is reducible with respect to I , which yields a contradiction.By the proof of the last theorem, we have the following corollary. Corollary 2.5.
Suppose that A is a strong M-tensor. If the M-Teq (1.1) with b ∈ R n + hasa nonnegative ¯ x with zero elements, then A is reducible with respective some I ⊆ I ( b ) . The following theorem characterizes a nonnegative solution of the M-Teq (1.1) with b ∈ R n + . Theorem 2.6.
Suppose that A is a strong M-tensor and b ∈ R n + . Then the M-Teq (1.1) hasa nonnegative solution with zero elements if and only if A is reducible with respect to some I ⊆ I ( b ) . Moreover, for a nonnegative solution ¯ x of the M-Teq (1.1) , ¯ x i = 0 iff i ∈ I .Proof. The “only if” part follows from Corollary 2.5 directly.Suppose that tensor A is reducible with respect to some I ⊆ I ( b ). It is easy to see thatthe M-tensor equation (1.1) has a solution ¯ x with ¯ x I = 0. Denote I c = [ n ] \ I . Considerthe lower dimension M-tensor equation A I c x m − I c − b I c = 0 . (2.4)Since b I c ≥
0, the last equation has a nonnegative solution ¯ x I c , which together with ¯ x I = 0forms a nonnegative solution to the M-tensor equation (1.1).If A I c is irreducible, then I = I is the desired index set. Otherwise, A I c is reducible withrespect to some I ⊂ I c satisfying I ⊆ b I c . We consider the lower dimensional M-tensorequation (2.4). Following a similar discussion to the above process, we can get another lowerdimensional tensor equation whose nonnegative solution together with some zeros forms anonnegative solution to (1.1). Continuing this process finitely many times, we can get adesirable index set I ⊂ I ( b ). 7 emark 2.7. The above theorem provides a way to reduce the size of an M-tensor equationwith b ∈ R n + . Specifically, in the case tensor A is reducible with respect to some I ⊆ I ( b ),we can get a solution to (1.1) by finding a positive solution to the lower dimensional tensorequation A I c x m − I c − b I c = 0 , where I c = [ n ] \ I .As a direct corollary of the last theorem, we have the following results, which gives anecessary and sufficient condition for the M-tensor equation (1.1) with b ∈ R n + to have apositive solution. Corollary 2.8.
Suppose that A is a strong M-tensor and b ∈ R n + . Then there is an indexset I ⊆ I ( b ) such that every nonnegative solution to the following lower dimensional tensorequation with I c = [ n ] \ I A I c x m − I c − b I c = 0 is positive. Moreover, the positive solution x I c of the last equation together with x I = 0 forms a nonnegative solution to the M-Teq (1.1) . The following lemma gives another interesting property for an M-Teq.
Lemma 2.9.
Suppose that A is a strong M-tensor and b ∈ R n + . Suppose further that everynonnegative solution of the M-Teq (1.1) is positive. Then there is an index set J ⊆ I ( b ) such that for each i ∈ J , there are at least one i j J , ≤ j ≤ m such that a ii ...i m = 0 .Proof. Let J = I ( b ). It is easy to see that there must be at least one i ∈ J and at leastone i j J , 2 ≤ j ≤ m such that a ii ...i m = 0. Otherwise, the M-Teq has a nonnegativesolution ˜ x with ˜ x i = 0, ∀ i ∈ J , with yields a contradiction.If J does not meet the requirement, we get an index set J ⊂ J consisting of all indices i ∈ J that does not meet the requirement. If J still is not the desired index set, we canfurther get a small index set J ⊂ J . We proceed the process. At last, we get the desiredindex J .Based on the last lemma, we can show the nonsingularity property of the Jacobian F ′ at the positive solutions. Theorem 2.10.
Suppose that the M-Teq (1.1) with a strong M-tensor A and b ∈ R n + hasa positive solution ¯ x . Then the Jacobian F ′ (¯ x ) is a nonsingular M-matrix. Let f ( y ) = F ( y [ m − ] ) and ¯ y = ¯ x [ m − . Then f ′ (¯ y ) is also a nonsingular M-matrix.Proof. It is easy to derive for any y > f ′ ( y ) = A ( y [ m − ] ) m − diag( y [ m − − ) . It shows that the nonsingularity of f ′ (¯ y ) is the same as the nonsingularity of F ′ (¯ x ).8et J ⊆ I ( b ) be the index set specified by Lemma 2.9 and I = [ n ] \ J . Write the Jacobianmatrix F ′ (¯ x ) as the block form F ′ (¯ x ) = ( m − A ¯ x m − = A II A IJ A JI A JJ ! . Since ¯ x is a positive solution of (1.1), it follows from Lemma 2.9 that A JI has no zero rows.That ¯ x is a solution to (1.1) yields F ′ (¯ x )¯ x = ( m − b . Writing it as block form, we get ( A II ¯ x I + A IJ ¯ x J = ( m − b I ,A JI ¯ x I + A JJ ¯ x J = ( m − b J . (2.5)It follows from the last equality of the above system that A JJ ¯ x J = ( m − b J − A JI ¯ x I = − A JI ¯ x I > . Since A JJ is a Z-matrix, the last inequality implies that A JJ is a nonsingular M-matrix.It then suffices to show that the Schur complement A II − A IJ A − JJ A JI is a nonsingularM-matrix.If J = I ( b ), then I = I + ( b ). We get from the first equality of (2.5), (cid:0) A II − A IJ A − JJ A JI (cid:1) ¯ x I = ( m − b I > . Clearly, matrix A II − A IJ A − JJ A JI is a Z-matrix. Consequently, the last inequality showsthat the Schur complement of A JJ is a nonsingular M-matrix too. Therefore, A = F ′ (¯ x ) isa nonsingular M-matrix.In the case J ⊂ I ( b ),we denote J = J , I = I and A = A I I − A I J A − J J A J I . Thento show F ′ (¯ x ) is a nonsingular M-matrix is equivalent to show that the lower dimensionalZ-matrix A is a nonsingular M-matrix. It is clear that ¯ x I satisfies the lower dimensionalsystem of linear equations A ¯ x I = ( m − b I . Similar to above arguments, we can get a partition ( I , J ) to the index set I thatpossesses the same properties as ( I , J ). Repeat the process finitely many times, we canget J t = I ( b I t − ). As a result, we can verify that F ′ (¯ x ) is a nonsingular M-matrix. b ∈ R n ++In this section, we propose a Newton method to find the unique positive solution to (1.1)with b ∈ R n ++ . Throughout this section, without specification, we always suppose that thefollowing assumption holds. Assumption 3.1.
Tensor A is a semi-symmetric and strong M-tensor, and b ∈ R n ++ .9ecently, He, Ling, Qi and Zhou [13] developed a Newton method for solving the M-Teq(1.1) with b ∈ R n ++ . By making a variable transformation x = y [ m ] , they formulated theequation to the following equivalent nonlinear equation: W ( y ) = D ( y ) · F ( y [ m ] ) = D ( y ) · A (cid:16) y [ m ] (cid:17) m − − D ( y ) · b = 0 , where D ( y ) = diag (cid:16) y m − i (cid:17) . The above equation has some nice properties such as thenonsingularity of the Jacobian W ′ ( y ) for any y >
0. In the case where A is symmetric, thetensor equation (1.1) is the stationary equation of the minimization problemmin ¯ f ( y ) = 1 m A (cid:16) y [ m ] (cid:17) m − b T (cid:16) y [ m ] (cid:17) because the gradient of ¯ f ( y ) is ∇ ¯ f ( y ) = 1 m W ( y ) = 1 m D ( y ) · ∇ f ( y [ m ] ) . In what follows, we propose a Newton method for finding the unique positive solution ofthe M-Teq (1.1). Our idea to develop the Newton method is similar to but different fromthat in [13]. Details are given below.Since our purpose is to get a positive solution of the M-Teq (1.1), we restrict x ∈ R n ++ .Making a variable transformation y = x [ m − , we formulate the M-Teq (1.1) as f ( y ) = F ( y [ m − ] (cid:17) = A (cid:16) y [ m − ] (cid:17) m − − b = 0 . (3.1)A direct computation gives f ′ ( y ) = A (cid:16) y [ m − ] (cid:17) m − diag (cid:16) y [ m − − (cid:17) . It follows that f ′ ( y ) y = A (cid:16) y [ m − ] (cid:17) m − diag (cid:16) y [ m − − i (cid:17) y = A (cid:16) y [ m − ] (cid:17) m − = f ( y ) + b. Clearly, the positive solutions of the M-Teq (1.1) are positive solutions of the followingnonlinear equation: E ( y ) △ = diag (cid:16) y [ − (cid:17) f ( y ) = (cid:16) y [ − ◦ f ( y ) (cid:17) = 0 . (3.2)The Jacobian of E ( y ) is E ′ ( y ) = diag (cid:16) y [ − (cid:17) f ′ ( y ) − diag ( f ( y ))diag (cid:16) y [ − (cid:17) = diag (cid:16) y [ − (cid:17)h f ′ ( y ) − diag ( f ( y ))diag (cid:16) y [ − (cid:17)i . It is a non-symmetric Z-matrix. For any y >
0, it holds that E ′ ( y ) y = diag (cid:16) y [ − (cid:17)h f ′ ( y ) y − diag ( f ( y )) (cid:16) y [ − (cid:17) y i = diag (cid:16) y [ − (cid:17)h A (cid:16) y [ m − ] (cid:17) m − − f ( y ) i = diag (cid:16) y [ − (cid:17) b > . Consequently, we have got the following proposition.10 roposition 3.2.
Let E : R n ++ → R be defined by (3.2) . For any y > , the Jacobian E ′ ( y ) is an M-matrix. Moreover, the equation (3.2) has a unique positive solution that isthe unique positive solution to the M-Teq (1.1) . We are going to develop a Newton method for solving the nonlinear equation (3.2) inwhich the Newton direction d k is the solution to the system of linear equations E ′ ( y k ) d + E ( y k ) = 0 , i.e., diag (cid:16) y [ − k (cid:17)h f ′ ( y k ) − diag (cid:16) f ( y k ) y k (cid:17)i d + diag (cid:16) y [ − k (cid:17) f ( y k ) = 0 , or equivalently h f ′ ( y k ) − diag (cid:16) f ( y k ) y k (cid:17)i d + f ( y k ) = 0 (3.3)Here diag (cid:16) f ( y k ) y k (cid:17) is a diagonal matrix whose diagonals are f i ( y k )( y k ) i , i = 1 , , . . . , n . Wecan regard d k as an inexact Newton method for solving the equation f ( y ) = 0 because theNewton equation (3.3) can be written as f ′ ( y k ) d k + f ( y k ) = r k , r k = diag (cid:16) f ( y k ) y k (cid:17) d k = O ( k f ( y k ) k k d k k ) , if y k > y k ( α ) = y k + αd k . Then y k ( α ) satisfies h f ′ ( y k ) − diag (cid:16) f ( y k ) y k (cid:17)i y k ( α ) = h f ′ ( y k ) − diag (cid:16) f ( y k ) y k (cid:17)i y k − αf ( y k ) = b − αf ( y k ) . Since the Jacobian E ′ ( y k ) = diag (cid:16) y [ − k (cid:17)h f ′ ( y k ) − diag (cid:16) f ( y k ) y k (cid:17)i is an M-matrix and y k >
0, it is clear that the matrix f ′ ( y k ) − diag (cid:16) f ( y k ) y k (cid:17) is an M-matrix too. Therefore, the inequality y k ( α ) > b − αf ( y k ) > . (3.4)Let ¯ α max k = min n b i f i ( y k ) : f i ( y k ) > o . (3.5)It is clear that y k + αd k > , ∀ α ∈ (0 , ¯ α max k ) . The iterative process of the Newton method is stated as follows.
Algorithm 3.3. ( Newton’s Method ) 11 nitial.
Given constant σ, ρ ∈ (0 ,
1) and ǫ >
0. Select an initial point x >
0. suchthat y = x [ m − satisfies f ( y ) < b . Let k = 0. Step 1.
Stop if k E ( y k ) k < ǫ . Step 2.
Solve the system of linear equations (3.3) to get d k . Step 3.
For given constant σ ∈ (0 , α k = max { ρ i : i = 0 , , . . . } such that y k + α k d k > k E ( y k + α k d k ) k ≤ (1 − σα k ) k E ( y k ) k , σ ∈ (0 , . (3.6)is satisfied. Step 3.
Let y k +1 = y k + α k d k . Go to Step 1. Remark 3.4.
It is easy to see that the inequality (3.4) is guaranteed if f ( y k ) < b . So, at thebeginning, we select y > f ( y ) < b and at each iteration, we let y k +1 = y k + α k d k such that f ( y k +1 ) < b . In this way, the inequalities f ( y k ) < b for all k . Lemma 3.5.
Let { y k } be generated by Algorithm . Then there is a positive constant c such that y k ≥ c e , ∀ k ≥ , (3.7) where e = (1 , , . . . , T .Proof. It is clear that the sequence of the function evaluations {k E ( y k ) k} is decreasing andhence bounded by some constant M >
0, i.e., k E ( y k ) k ≤ M .
Since A is an M-tensor, there is a constant s > B ≥ A = s I − B , where I stands for the identity tensor whose diagonals are all ones and allother elements are zeros.By the definition of E ( y ), we have E ( y ) = s e − diag ( y − ) B (cid:16) y [ m − ] (cid:17) m − − b ◦ ( y [ − ) . Since
B ≥
0, the last inequality implies for any y > i ∈ [ n ] | E i ( y ) | ≥ b i y i + y − i (cid:16) B (cid:16) y [ m − ] (cid:17) m − (cid:17) i − s ≥ b i y i − s. Suppose there is an index i and an infinite set K such that lim k ∈ K, k →∞ ( y k ) i = 0. We have M ≥ lim k ∈ K, k →∞ | E i ( y k ) | ≥ lim k ∈ K, k →∞ b i ( y k ) i − s = + ∞ , which yields a contradiction. The contradiction shows that the inequality in (3.5) is satisfiedwith some positive constant c . 12he following theorem establishes the global convergence of the proposed method. Theorem 3.6.
Suppose that the sequence { y k } generated by Algorithm is bounded. Then { y k } converges to the unique positive solution to the M-Teq (1.1) .Proof. We first show that the maximum step length ¯ α max k satisfying (3.5) can be boundedaway from zero. That is, there is a constant ¯ α such that¯ α max k ≥ ¯ α, ∀ k ≥ . (3.8)Indeed, it follows from the last lemma that M ≥ | E i ( y k ) | = | f i ( y k ) | ( y k ) i . Since { y k } is bounded, the last inequality implies that for each i , {| f i ( y k ) |} is bounded too.By the definition of ¯ α max k , it is bounded away from some constant ¯ α . Consequently, theinequality (3.8) is satisfied for all k ≥ y of { y k } that is a positive solutionto (1.1).Suppose { y k } K → ¯ y . By Lemma 3.5, it is clear that ¯ y >
0. Consequently, E ′ (¯ y ) is anM-matrix. Moreover, lim k ∈ K, k →∞ d k = − E ′ (¯ y ) − E (¯ y ) △ = ¯ d. Without loss of generality, we let lim k ∈ K, k →∞ α k = ˜ α .If ˜ α >
0, then the inequality (3.6) shows that {k E ( y k +1 ) k} k → α = 0, then when k ∈ K is sufficiently large, the inequality (3.6) is not satisfied with α ′ k = ρ − α k , i.e., k E ( y k + α ′ k d k ) k − k E ( y k ) k > − σα ′ k k E ( y k ) k , σ ∈ (0 , . Dividing both sizes of the inequality by α ′ k and then taking limits as k → ∞ with k ∈ K ,we get −k E (¯ y ) k = E (¯ y ) T E ′ (¯ y ) ¯ d ≥ − σ k E (¯ y ) k , which implies E (¯ y ) = 0.Since {k E ( y k ) k} converges, it follows from Lemma 3.5 that every accumulation pointof { y k } is a positive solution to (1.1). However, the positive solution of (1.1) is unique.Consequently, the whole sequence { y k } converges to the unique positive solution to (1.1).By a standard argument, it is not difficult to show that the convergence rate of { y k } isquadratic. Theorem 3.7.
Let the conditions in Theorem hold. Then the convergence rate of { y k } is quadratic. An Extension
In this section, we extend the Newton method proposed in the last section to the M-Teq (1.1)with b ∈ R n + . In the case b has zero elements, the M-Teq may have multiple nonnegativeor positive solutions. Our purpose is to find one nonnegative or positive solution of theequation.We see from the definition of E ( y ) that the function E ( y ) and its Jacobian are not welldefined at a point with zero elements. Therefore, the Newton method proposed in the lastsection can not be applied to find a nonnegative solution with zero elements. Fortunately,from Corollary 2.8, we can get a nonnegative solution of (1.1) by finding a positive solutionto a lower dimensional M-Teq.Without loss of generality, we make the following assumption. Assumption 4.1.
Suppose that tensor A is a semi-symmetric and strong M-tensor, and b ∈ R n + . Moreover, every nonnegative solution of the M-Teq (1.1) is positive.Similar to the Newton method by He, Ling, Qi and Zhou [13], we propose another Newtonmethod, which we call a regularized Newton method, such that the method is still globallyand quadratically convergent without assuming the boundedness of the generated sequenceof iterates.It is easy to see that the M-Teq (1.1) is equivalent to the following nonlinear equation E ( t, y ) △ = (cid:18) ty [ − ◦ f ( y ) + ty (cid:19) = (cid:18) tE ( t, y ) (cid:19) = 0 , (4.1)where E ( t, y ) = E ( y ) + ty = y [ − ◦ f ( y ) + ty. The Jacobian of E ( t, y ) is E ′ ( t, y ) = (cid:18) y E ′ y ( t, y ) (cid:19) , where E ′ y ( t, y ) = E ′ ( y ) + tI satisfying E ′ y ( t, y ) y = E ′ ( y ) y + ty = y [ − ◦ b + ty > , ∀ y ∈ R n ++ , ∀ t > . Since E ′ y ( t, y ) is a Z-matrix, the last inequality shows that it is a nonsingular M-matrix. Asa result, for any t > y ∈ R n ++ , the Jacobian E ′ ( t, y ) is nonsingular.Now, we propose a Newton method for solving the equivalent nonlinear equation (4.1)to the M-Teq (1.1). The idea is similar to the Newton method by He, Ling, Qi and Zhou[13]. Details are given below.Given constant γ ∈ (0 , θ ( t, y ) = 12 k E ( t, y ) k , β ( t, y ) = γ min { , k E ( t k , y k ) k } . The subproblem of the method is the following system of linear equations: E ′ ( t k , y k ) d k + E ( t k , y k ) = β ( t k , y k ) e , (4.2)14here e = (1 , , . . . , T ∈ R n +1 . Let d k = ( d tk , d yk ).Suppose t k ≤ ¯ t with ¯ t satisfying ¯ tγ <
1. Then the Newton direction d k satisfies ∇ θ ( t k , y k ) T d k = E ( t k , y k ) T E ′ ( t k , y k ) d k = −k E ( t k , y k ) k + β ( t k , y k ) ≤ − (1 − γ ¯ t ) k E ( t k , y k ) k . (4.3)As a result, for given constant σ ∈ (0 , θ ( t k + α k d tk , y k + α k d yk ) ≤ [1 − σ (1 − γ ¯ t ) α k ] θ ( t k , y k ) (4.4)is satisfied for all α k > Algorithm 4.2. Regularized Newton MethodInitial.
Given constants γ, σ, ρ ∈ (0 , ǫ > t > tγ <
1. Giveninitial point x > t = ¯ t . Let y = x [ m − and k = 0. Step 1.
Stop if k E ( t k , y k ) k ≤ ǫ . Step 2.
Solve the system of linear equations (4.2) to get d k . Step 3.
Find α k = max { ρ i : i = 0 , , . . . } such that y k + ρ i d yk > α k = ρ i . Step 4.
Let y k +1 = y k + α k d yk and t k +1 = t k + α k d tk . Step 5.
Let k := k + 1. Go to Step 1.Following a similar argument as the proof of Lemma 3.2 of [13], it is not difficult to getthe following proposition. It particularly shows that the above algorithm is well-defined. Proposition 4.3.
Suppose that A is a strong M-tensor and b ∈ R n + . Then the sequence ofiterates { ( t k , y k ) } generated by Algorithm satisfies < t k +1 ≤ t k ≤ ¯ t and t k > ¯ tβ ( t k , y k ) . In addition, the sequence of function evaluations { θ ( t k , y k ) } is decreasing. Since A is an M-tenor, there are a constant s > B = ( b i ...i m )such that A = s I − B , where I is the identity tensor whose diagonal entities are all onesand all other elements are zeros. By the definition of E ( y ), it is easy to get E ( y ) = s e − y [ − ◦ B (cid:16) y [ m − ] (cid:17) m − − y [ − ◦ b. Lemma 4.4.
Suppose that A is a strong M-tensor and b ∈ R n + . Then the sequence ofiterates { y k } generated by Algorithm is bounded away from zero. In other words, thereis a constant η > such that ( y k ) i ≥ η, ∀ k ≥ , ∀ i = 1 , , . . . , n. roof. Suppose that there is an index i and a subsequence { y k } K such that lim k →∞ , k ∈ K ( y k ) i =0. Without loss of generality, we suppose { y k } K → ¯ y , where some elements of ¯ y may be+ ∞ . Denote I = { i : ¯ y i = 0 } and I c = [ n ] \ I c . Since Let { θ ( t k , y k ) } is decreasing, itis bounded and so is the sequence {k E ( t k , y k ) k} . Let C > {k E ( t k , y k ) k} .For each i ∈ I , it holds that C ≥ | E i ( t k , y k ) | = (cid:12)(cid:12)(cid:12) y k ) i X i ,...,i m a ii ...i m (cid:16) ( y k ) m − i · · · ( y k ) m − i m (cid:17) − b i ( y k ) i + t k ( y k ) i (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) s − y k ) i X i ,...,i m b ii ...i m (cid:16) ( y k ) m − i · · · ( y k ) m − i m (cid:17) − b i ( y k ) i + t k ( y k ) i (cid:12)(cid:12)(cid:12) ≥ X i ,...,i m b ii ...i m (cid:18) ( y k ) i ( y k ) i · · · ( y k ) i m ( y k ) i (cid:19) m − + b i ( y k ) i − t k ( y k ) i − s ≥ X i ,...,i m ∈ I c b ii ...i m (cid:18) ( y k ) i ( y k ) i · · · ( y k ) i m ( y k ) i (cid:19) m − + b i ( y k ) i − t k ( y k ) i − s. Notice that for any i ∈ I c , ¯ y i >
0. Since t k ≤ ¯ t and ( y k ) i →
0, as k → ∞ with k ∈ K , thelast inequality implies b i = 0 and a ii ...i m = b ii ...i m = 0, ∀ i , . . . , i m ∈ I c . It means thattensor A is reducible with respect to index set I . It then follows from Theorem 2.6 that theM-Teq (1.1) has a nonnegative solution that has zero elements. It is a contradiction. Thecontradiction shows that { y k } is bounded away from zero. Lemma 4.5.
Suppose that A is a strong M-tensor and b ∈ R n + . If there is a ˜ t > such that t ≥ ˜ t , then the sequence of iterates { y k } generated by Algorithm is bounded.Proof. Denote by i k the index satisfying ( y k ) i k = k y k k ∞ . Since { θ ( t k , y k ) } has an upperbound, so is {k E ( t k , y k ) k} . Let C be an upper bound of {k E ( t k , y k ) k} . It is clear that (cid:12)(cid:12)(cid:12) X i ,...,i m a i k i ...i m (cid:18) ( y k ) i ( y k ) i k · · · ( y k ) i m ( y k ) i k (cid:19) m − (cid:12)(cid:12)(cid:12) ≤ X i ,...,i m | a i k i ...i m | △ = ˜ a i k is bounded. Therefore, we obtain C ≥ k E ( t k , y k ) k≥ (cid:12)(cid:12)(cid:12) y k ) i k X i ,...,i m a i k i ...i m (cid:16) ( y k ) m − i · · · ( y k ) m − i m (cid:17) − b i k ( y k ) i k + t k ( y k ) i k (cid:12)(cid:12)(cid:12) ≥ t k ( y k ) i k − ˜ a i k − b i k ( y k ) i k . The last inequality together with t k ≥ ˜ t implies that {k y k k} is bounded.The following theorem establishes the global convergence of Algorithm 4.2. Theorem 4.6.
Suppose that A is a strong M-tensor and b ∈ R n + . Then every accumulationpoint of the sequence of iterates { ( t k , y k ) } generated by Algorithm is a positive solutionto the M-Teq (1.1) . roof. It suffices to show that the sequence { θ ( t k , y k ) } converges to zero by contradiction.Suppose on the contrary that there is a constant δ > θ ( t k , y k ) ≥ δ , ∀ k ≥ t △ = lim t →∞ t k ≥ ¯ t lim t →∞ β ( t k , y k ) ≥ ¯ tγ min { , δ } > . By Lemma 4.5, { y k } is bounded. Let the subsequence { y k } K converges to some point ¯ y .Lemma 4.4 ensures ¯ y >
0. It is easy to show that the Jacobian E ′ (˜ t, ¯ y ) is a nonsingularM-matrix. Consequently, { d k } K is bounded. Without loss of generality, we suppose { d k } K converges to some ¯ d . Since ¯ y >
0, there is a constant α min > y k + α k d k > ∀ α k ∈ (0 , α min ). Let ¯ α = lim inf k →∞ , k ∈ K α k . If ¯ α >
0, the line search condition (4.4)implies θ (¯ y, ˜ t ) = 0. If ¯ α = 0, then when k is sufficiently large, the inequality (4.4) is notsatisfied with α ′ k = α k ρ − , i.e., θ ( t k + α ′ k d tk , y k + α ′ k d yk ) − θ ( t k , y k ) ≥ − σ (1 − γ ¯ t ) α ′ k θ ( t k , y k ) . Dividing both sizes of the last inequality by α ′ k and then taking limits as k → ∞ with k ∈ K ,we get ∇ θ (˜ t, ¯ y ) T ¯ d ≥ − σ (1 − γ ¯ t ) θ (˜ t, ¯ y ) . On the other hand, by taking limits in both sizes of (4.3) as k → ∞ with k ∈ K , we obtain ∇ θ (˜ t, ¯ y ) T ¯ d ≤ − − γ ¯ t ) θ (˜ t, ¯ y ) . Since σ ∈ (0 , θ (¯ y, ˜ t ) = 0, which yields a contradiction.As a result, we claim that { θ ( t k , y k ) } converges to zero. The proof is complete.The last theorem has shown that every accumulation is a positive solution to the M-Teq (1.1). However, it does not the existence of the accumulation point. The followingtheorem shows that the sequence { y k } is bounded. As a result, it ensure the existence ofthe accumulation point. Theorem 4.7.
Suppose that A is a strong M-tensor and b ∈ R n + . Then the sequence { y k } generated by Algorithm is bounded.Proof. First, similar to the proof of Lemma 4.5, it is not difficult to show that the sequence { t k y k } is bounded.Case (i), { t k y k } →
0. Since { θ ( y k , t k ) } →
0, we immediately have { E ( y k ) } →
0. Denote µ k = k y k k ∞ , ˜ y k = µ − k y k and ˜ b k = µ − k b . Clearly, the sequence { ˜ y k } is bounded. If { y k } isunbounded, then there is a subsequence { µ k } K → ∞ , and hence { ˜ b k } K →
0. Without lossof generality, we suppose that the subsequence { ˜ y k } K converges to some ˜ y ≥
0. Denote by J the set of indices i satisfying ˜ y i >
0. Obviously, J = ∅ .For some i ∈ J , satisfies y i = k y k k ∞ , we have | E i ( y k , t k ) | = (cid:12)(cid:12)(cid:12) y k ) i X i ,...,i m a ii ...i m (cid:16) ( y k ) m − i · · · ( y k ) m − i m (cid:17) − b i ( y k ) i + t k ( y k ) i (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) X i ,...,i m a ii ...i m (cid:16) (˜ y k ) m − i · · · (˜ y k ) m − i m (cid:17) − (˜ b k ) i + t k ( y k ) i (cid:12)(cid:12)(cid:12) . k → ∞ with k ∈ K yields0 = X i ,...,i m a ii ...i m (cid:16) ˜ y m − i · · · ˜ y m − i m (cid:17) = X i ,...,i m ∈ J a ii ...i m (cid:16) ˜ y m − i · · · ˜ y m − i m (cid:17) , Let A J be the principal subtensor of A with elements a i i ...i m , ∀ i , i , . . . , i m ∈ J . It is astrong M-tensor but A J (cid:16) ˜ y [ m − ] (cid:17) m − J = 0 with ˜ y = 0. It is a contradiction. Consequently, { y k } is bounded.Case (ii), there are at least one i such that lim inf k →∞ t k ( y k ) i >
0. In other words, thereis a subsequence { t k y k } K → ˜ y ≥ y i > i . Again, denote by J the set of indices for satisfying ˜ y i >
0. Since { t k } →
0, it is easy to see thatlim k →∞ , k ∈ K ( y k ) i = + ∞ , ∀ i ∈ J. Denote ˜ y k = t k y k . Similar to Case (i), we can get We derive for any i ∈ J | E i ( y k , t k ) | = (cid:12)(cid:12)(cid:12) y k ) i X i ,...,i m a ii ...i m (cid:16) ( y k ) m − i · · · ( y k ) m − i m (cid:17) − b i ( y k ) i + t k ( y k ) i (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) y k ) i X i ,...,i m a ii ...i m (cid:16) (˜ y k ) m − i · · · (˜ y k ) m − i m (cid:17) − b i ( y k ) i + (˜ y k ) i (cid:12)(cid:12)(cid:12) . Taking limits in both sizes of the equality as k → ∞ with k ∈ K yields X i ,...,i m a ii ...i m (cid:16) ˜ y m − i · · · ˜ y m − i m (cid:17) + ˜ y i = X i ,...,i m ∈ J a ii ...i m (cid:16) ˜ y m − i · · · ˜ y m − i m (cid:17) + ˜ y i , ∀ i ∈ J. It contradicts Theorem 2.3 (ii).The proof is complete.Similar to theorem 3.3 of [13], we have the following theorem.
Theorem 4.8.
Let the conditions in Assumption hold, then the sequence of iterates { t k , y k } generated by Algorithm converges to a positive solution of the equation . Andthe convergence rate is quadratic. In this section, we do numerical experiments to test the effectiveness of the proposed meth-ods. We implemented our methods in Matlab R2015b and ran the codes on a personalcomputer with 2.30 GHz CPU and 8.0 GB RAM. We used a tensor toolbox [1] to proceedtensor computation.While do numerical experiments, similar to [12, 13], we solved the tensor equationˆ F ( x ) = ˆ A x m − − ˆ b = 0instead of the tensor equation (1.1), where ˆ A := A /ω and ˆ b := b/ω with ω is the largestvalue among the absolute values of components of A and b . The stopping criterion is set to k ˆ F ( x k ) k ≤ − .
18r the number of iteration reaches to 300. The latter case means that the method is failurefor the problem.
Problem 1. [8] We solve tensor equation (1.1) where A is a symmetric strong M-tensorof order m ( m = 3 , ,
5) in the form A = s I − B , where tensor B is symmetric whose entriesare uniformly distributed in (0 , s = (1 + 0 . · max i =1 , ,...,n ( B e m − ) i , where e = (1 , , . . . , T . Problem 2. [32] We solve tensor equation (1.1) where A is a symmetric strong M-tensorof order m ( m = 3 , ,
5) in the form A = s I − B , and tensor B is a nonnegative tensor with b i i ...i m = | sin( i + i + . . . + i m ) | , and s = n m − . Problem 3. [8] Consider the ordinary differential equation d x ( t ) dt = − GMx ( t ) , t ∈ (0 , , with Dirichlet’s boundary conditions x (0) = c , x (1) = c , where G ≈ . × − N m /kg and M ≈ . × is the gravitational constant andthe mass of the earth.Discretize the above equation, we have x = c , x i − x i x i − − x i x i +1 = GM ( n − , i = 2 , , · · · , n − ,x n = c . It is a tensor equation, i.e., A x = b, where A is a 4-th order M tensor whose entries are a = a nnnn = 1 ,a iiii = 2 , i = 2 , , · · · , n − ,a i ( i − ii = a ii ( i − i = a iii ( i − = − / , i = 2 , , · · · , n − ,a i ( i +1) ii = a ii ( i +1) i = a iii ( i +1) = − / , i = 2 , , · · · , n − , and b is a positive vector with b = c ,b i = GM ( n − , i = 2 , , · · · , n − ,b n = c . Problem 4. [18] We solve tensor equation (1.1) where A is a non-symmetric strongM-tensor of order m ( m = 3 , ,
5) in the form A = s I − B , and tensor B is nonnegativetensor whose entries are uniformly distributed in (0 , s is set to s = (1 + 0 . · max i =1 , ,...,n ( B e m − ) i . roblem 5. We solve tensor equation (1.1) where A is a lower triangle strong M-tensorof order m ( m = 3 , ,
5) in the form A = s I − B , and tensor B is a strictly lower triangularnonnegative tensor whose entries are uniformly distributed in (0 , s is setto s = (1 − . · max i =1 , ,...,n ( B e m − ) i . For Problem 4 and 5, we need to semi-symmetrize the tensor A , i.e., find a semi-symmetric tensor ˜ A such that A x m − = ˜ A x m − . The time of semi-symmetrize the tensor is not included in CPU time.We first test the performance of the Inexact Newton method. We set the start point x = ε e , where parameter ε is selected to satisfy f ( y ) < b . We set the parameter σ = 0 . ρ = 0 .
5. And b is uniformly distributed in (0 ,
1) except the b in the problem 3.For the stability of numerical results, we test the problems of different sizes. For eachpair ( m, n ), we randomly generate 100 tensors A and b . In order to test the effectiveness ofthe proposed method, we compare Inexact Newton method with the QCA method in [13].We take parameters δ = 0 . , γ = 0 . , σ = 0 . , ¯ t = 2 / (5 γ ) as the same as in [13]. The resultsare listed in Tables 1, whereIR = the number of iteration steps of the Inexact Newton methodthe number of iteration steps of the QCA methodand TR = the CPU time used by the Inexact Newton methodthe CPU time used by the QCA method . Table 1:
Comparison between Inexact Newton method and QCA method with b ∈ R n + . ( m, n ) (3,10) (3,100) (3,300) (3,500) (4,10) (4,50) (4,100) (5,10) (5,30)IR Problem 1 89.2% 91.5% 91.5% 91.0% 93.0% 93.7% 94.3% 95.2% 96.3%Problem 2 91.0% 91.4% 90.2% 90.5% 95.7% 94.8% 93.1% 97.2% 96.2%Problem 3 - - - - 11.1% 9.1% 8.3% - -Problem 4 91.8% 91.2% 91.3% 90.5% 94.4% 94.7% 93.2% 97.1% 93.9%Problem 5 89.8% 90.4% 89.6% 89.3% 95.2% 91.5% 93.0% 95.1% 95.6%TR Problem 1 48.0% 66.6% 87.7% 88.2% 67.3% 92.3% 94.0% 80.0% 97.0%Problem 2 50.0% 73.8% 88.8% 88.8% 67.4% 94.0% 93.7% 79.0% 96.6%Problem 3 - - - - 20.3% 15.4% 14.2% - -Problem 4 54.1% 73.3% 89.6% 89.6% 66.0% 94.7% 94.1% 74.6% 95.4%Problem 5 45.7% 74.1% 87.4% 88.5% 59.1% 92.3% 95.0% 74.6% 99.4% We then test the effectiveness of the Regularized Newton method. We set the initialpoint x = 0 . ∗ e and b ∈ R n + has 0 zero elements except the problem 3. We first generatea vector b ∈ R n whose elements are uniformly distributed in (0 , b i = (cid:26) b i , if b i ≤ . , , if b i > . . to get a vector b ∈ R n + . In order to get the positive solution of the problem 5, the firstcomponent of vector b can’t be equal to 0, so we set the first component b = 0 . .
20e compare the Regularized Newton Method with QCA method. We take the pa-rameters σ = 0 . , ρ = 0 . , γ = 0 . t = 0 .
01 in Regularized Newton Method and theparameters in QCA method is the same as above. The results are listed in Tables 2, whereIR = the number of iteration steps of the Regularized Newton methodthe number of iteration steps of the QCA methodand TR = the CPU time used by the Regularized Newton methodthe CPU time used by the QCA method . Table 2:
Comparison between Regularized Newton method and QCA method with b ∈ R n ++ . ( m, n ) (3,10) (3,100) (3,300) (3,500) (4,10) (4,50) (4,100) (5,10) (5,30)IR Problem 1 92.4% 59.7% 71.6% 67.4% 93.2% 67.2% 59.5% 95.7% 78.4%Problem 2 83.9% 61.5% 50.3% 49.4% 89.1% 56.0% 59.3% 87.0% 59.4%Problem 3 - - - - 83.3% 80.0% 81.0% -Problem 4 94.3% 65.8% 58.1% 60.9% 95.1% 67.7% 53.8% 95.5% 76.5%Problem 5 80.0% 81.0% 81.1% 81.4% 75.4% 77.6% 76.4% 72.4% 74.0%TR Problem 1 80.0% 78.7% 89.2% 82.4% 93.2% 77.2% 71.6% 86.8% 97.8%Problem 2 72.7% 72.9% 61.3% 60.4% 87.5% 61.4% 64.0% 94.1% 65.3%Problem 3 - - - - 76.5% 97.6% 98.1% - -Problem 4 80.6% 81.3% 65.7% 76.7% 92.3% 78.1% 65.0% 97.4% 99.6%Problem 5 72.1% 94.9% 89.6% 91.0% 86.9% 79.3% 77.7% 82.5% 80.8% The datas in Table 1 and 2 show that for all test problems the Inexact Newton methodand the Regularized Newton method are better than QCA method in terms of the numberof iterations and CPU time. It is worth noting that although the QCA method in [13] doesnot established the convergence property in the case of b ∈ R n + , we find that in the case of b ∈ R n + , the QCA method can still find the solution of the problem successfully. For theconvenience of readers, we only list the relative results. More detailed numerical results canbe found in the Appendix. References [1]
B.W. Bader, T.G. Kolda and others , MATLAB Tensor Toolbox Version 2.6,2015.[2]
X. Bai, H. He, C. Ling and G. Zhou , An efficient nonnegativity preserving algorithmfor multilinear systems with nonsingular M-tensors , arXiv:1811.09917v1, 2018.[3]
A. Berman, R. Plemmons , Nonnegative Matrices in the Mathematical Sciences ,SIAM, Philadephia, 1994.[4]
H. Bozorgmanesh, M. Hajarian and A.T. Chronopoulos , Interval Tensors andtheir application in solving multi-linear systems of equations , Comput. Math. Appl., 79(2020), pp. 697–715. 215]
M. Brazell, N, Li, C. Navasca and C. Tamon , Solving multilinear systems viatensor inversion , SIAM J. Matrix Anal. Appl., 34 (2013), pp. 542–570.[6]
K.C. Chang, K. Pearson and T. Zhang , Primitivity, the convergence of the NQZmethod, and the largest eigenvalue for nonnegative tensors . SIAM J. Matrix Anal. Appl.,32 (2011), pp. 806–819.[7]
W. Ding, L. Qi and Y. Wei , M-tensors and nonsingular M-tensors , Linear AlgebraAppl., 10 (2013), pp. 3264–3278.[8]
W. Ding and Y. Wei , Solving multi-linear systems with M-tensors , J. Sci. Comput.,68 (2016), pp. 1–27.[9]
F. Facchinei and C. Kanzow , Beyond mootonicity in regularization methods fornonlinear complementarity problems , SIAM J. Control Optim., 37 (1999), pp. 150–1161.[10]
H.Y. Fan, L. Zhang, E.K. Chu and Y. Wei , Numerical solution to a linear equationwith tensor product structure , Numeri. Linear Algebra Appl., 24 (2017), pp. 1–22.[11]
M.S. Gowda, Z. Luo, L. Qi and N. Xiu , Z-tensor and complementarity problems ,arXiv:1510.07933, 2015.[12]
L. Han , A homotopy method for solving multilinear systems with M-tensors , Appl.Math. Lett., 69 (2017), pp. 49–54.[13]
H. He, C. Ling, L. Qi and G. Zhou , A globally and quadratically convergent al-gorithm for solving multilinear systems with M-tensors , J. Sci. Comput.,, 76 (2018),pp. 1718–1741.[14]
D. Kressner and C. Tobler , Krylov subspace methods for linear systems with tensorproduct structure , SIAM J. Matrix Anal. Appl., 31 (2010), pp. 1688–1714.[15]
D.H. Li, S. Xie and H.R. Xu , Splitting methods for tensor equations , Numeri. LinearAlgebra Appl., https://doi.org/10.1002/nla.2102, 2017.[16]
X. Li and M.K. Ng , Solving sparse non-negative tensor equations: algorithms andapplications , Front. Math. China., 10 (2015), pp. 649–680.[17]
Z. Li, Y.-H. Dai and H. Gao , Alternating projection method for a class of tensorequations , J. Comput. Appl. Math., 346 (2019), pp. 490–504.[18]
D.H. Li, H.B. Guan and X.Z. Wang , Finding a nonnegative solution to an M-tensorequation , arXiv: 1811.11343v1, 2018.[19]
L.H. Lim , Singular values and eigenvalues of tensors, a variational approach , In: Pro-ceedings of the 1st IEEE International Workshop on Computational Advances of Multi-tensor Adaptive Processing, 1 (2005), pp. 129–132.[20]
D. Liu, W. Li and S.W. Vong , The tensor splitting with application to solve multi-linear systems , J. Comput. Appl. Math., 330 (2018), pp. 75–94.[21] D. Liu, W. Li and S.W. Vong, A new preconditioned SOR method for solving multilinearsystems with an M-tensor. Calcolo, http://doi.org/10.1007/s10092-020-00364-8, 2020.[22]
W. Li Wen, D. Liu and S.W. Vong , Comparison results for splitting iterations forsolving multi-linear systems , Appl. Numer. Math., 134 (2018), pp. 105–121.2223]
M. Liang, B. Zheng, R. Zhao , Alternating iterative methods for solving tensor equa-tions with applications , Numer. Algorithms, 80(4) (2019), pp. 1437–1465.[24]
C.Q. Lv and C.F. Ma , A Levenberg-Marquardt method for solving semi-symmetrictensor equations , J. Comput. Appl. Math., 332 (2018), pp. 13–25.[25]
L. Qi , Eigenvalues of a real supersymmeytic tensor . J. Symbolic Comput., 40 (2005),pp. 1302–1324.[26]
L. Qi, H. Chen and Y. Chen , Tensor Eigenvalues and Their Applications , Springer,2018.[27]
L. Qi and Z. Luo , Tensor Analysis, Spectral Theory and Special Tensors , SIAM,Philadelphia, 2017.[28]
X. Wang, M. Che and Y. Wei , Neural networks based approach solving multi-linearsystems with M-tensors , Neurocomputing, 351 (2019), pp. 33–42.[29]
X. Wang, M. Che and Y. Wei , Neural network approach for solv-ing nonsingular multi-linear tensor systems , J. Comput. Appl. Math.,http://doi.org/10.1016/j.cam.2019.112569, 2020.[30]
X. Wang, M. Che and Y. Wei , Existence and uniqueness of positive solution for H + -tensor equations , Appl. Math. Lett., 98 (2019), pp. 191–198.[31] Z.J. Xie, X.Q. Jin and Y.M. Wei , A fast algorithm for solving circulant tensorsystems , Linear Multilinear Algebra., 65 (2017), pp. 1894–1904.[32]
Z.J. Xie, X.Q. Jin and Y.M. Wei , Tensor methods for solving symmetric M-tensorsystems , J. Sci. Comput., 74 (2018), pp. 412–425.[33]
Y. Xu, W. Gu and Z.H. Huang , Properties of the nonnegative solution set of multi-linear equations , Pac. J. Optim., 15 (2019), pp. 441–456.[34]
W. Yan, C. Ling, L. Ling and H. He , Generalized tensor equations with leadingstructured tensors , Appl. Math. Comput., 361 (2019), pp. 311–324.[35]
Y. Zhang, Q. Liu, Z. Chen , Preconditioned Jacobi type methodfor solving multi-linear systems with M-tensors , Appl. Math. Lett.,http://doi.org/10.1016/j.aml.2020.106287, 2020.[36]
L. Zhang, L. Qi and G. Zhou , M-tensor and some applications , SIAM J. MatrixAnal. Appl., 35 (2014), pp. 437–452.
A Detailed Numerical Results
In this section, we list the detailed numerical results of the proposed methods comparedwith QCA method. The results are listed in Tables 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, wherethe columns ‘Iter’, ‘Time’, ’Res’ and ’Ls-iter’ stand for the total number of iterations, thecomputational time (in second) used for the method, the residual k ˆ A x ( m − k − ˆ b k and thetotal number of iterations of linear search. 23able 4: Comparison between Inexact Newton method and QCA method on Problem 2.
Inexact Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 7.1 0.00020 1.0E-11 0 7.8 0.00040 4.7E-12 0(3,100) 9.6 0.00899 5.7E-12 0 10.5 0.01218 9.1E-12 0(3,300) 11.9 0.32718 8.6E-12 0 13.2 0.36855 1.1E-11 0(3,500) 12.4 1.49714 8.0E-12 0 13.7 1.68625 8.3E-12 0(4,10) 6.7 0.00029 8.6E-12 0 7.0 0.00043 9.8E-12 0(4,50) 9.1 0.04654 7.8E-12 0 9.6 0.04950 1.4E-11 0(4,100) 9.5 0.74050 1.5E-11 0 10.2 0.79047 1.7E-11 0(5,10) 6.9 0.00049 6.5E-12 0 7.1 0.00062 1.3E-11 0(5,30) 7.6 0.15073 1.0E-11 0 7.9 0.15600 1.4E-11 0Table 5:
Comparison between Inexact Newton method and QCA method on Problem 3.
Inexact Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(4,10) 1.0 0.00012 9.2E-15 1.0 9.0 0.00059 6.4E-12 1.0(4,50) 1.0 0.00947 2.0E-15 1.0 11.0 0.06154 4.2E-14 1.0(4,100) 1.0 0.14564 2.1E-15 1.0 12.0 1.02232 1.9E-15 1.0Table 3:
Comparison between Inexact Newton method and QCA method on Problem 1 .
Inexact Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 6.6 0.00024 8.5E-12 0 7.4 0.00050 6.1E-12 0(3,100) 9.7 0.00829 9.0E-12 0 10.6 0.01244 1.0E-11 0(3,300) 11.9 0.32238 9.8E-12 0 13.0 0.36766 1.7E-11 0(3,500) 12.1 1.44961 5.1E-12 0 13.3 1.64275 7.4E-12 0(4,10) 6.6 0.00033 5.0E-12 0 7.1 0.00049 9.2E-12 0(4,50) 8.9 0.04552 1.2E-11 0 9.5 0.04931 1.4E-11 0(4,100) 10.0 0.77240 1.3E-11 0 10.6 0.82138 6.7E-12 0(5,10) 6.0 0.00052 1.3E-11 0 6.3 0.00065 1.3E-11 0(5,30) 7.9 0.15599 8.9E-12 0 8.2 0.16078 1.2E-11 0From the data in the Tables 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, we can see that theproposed methods are effective for all test problems. In terms of the number of iterationsand CPU time, Inexact Newton method and Regularized Newton method are better thanQCA method, and the number of linear search of the Regularized Newton method are farless than that of the QCA method. 24able 6:
Comparison between Inexact Newton method and QCA method on Problem 4.
Inexact Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 6.7 0.00020 8.8E-12 0 7.3 0.00037 9.6E-12 0(3,100) 10.3 0.00934 7.9E-12 0 11.3 0.01274 1.1E-11 0(3,300) 11.6 0.31909 1.2E-11 0 12.7 0.35600 1.3E-11 0(3,500) 12.4 1.50356 7.9E-12 0 13.7 1.67812 8.5E-12 0(4,10) 6.8 0.00031 3.7E-12 0 7.2 0.00047 8.9E-12 0(4,50) 8.9 0.04571 1.4E-11 0 9.4 0.04826 1.1E-11 0(4,100) 9.6 0.74759 1.4E-11 0 10.3 0.79482 1.4E-11 0(5,10) 6.6 0.00047 5.4E-12 0 6.8 0.00063 1.0E-11 0(5,30) 7.7 0.15334 1.5E-11 0 8.2 0.16067 1.5E-11 0Table 7:
Comparison between Inexact Newton method and QCA method on Problem 5.
Inexact Newton method QCA( m, nm, n
Inexact Newton method QCA( m, nm, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 7.9 0.00016 7.2E-12 0.4 8.8 0.00035 4.7E-12 0.6(3,100) 10.3 0.00758 1.2E-11 0.1 11.4 0.01023 9.5E-12 0.6(3,300) 12.1 0.27135 1.2E-11 0 13.5 0.31054 8.8E-12 0.5(3,500) 12.5 1.44273 1.4E-11 0 14.0 1.63077 1.7E-11 0.4(4,10) 8.0 0.00026 4.2E-12 0.5 8.4 0.00044 1.1E-11 0.7(4,50) 9.7 0.04921 8.8E-12 0.2 10.6 0.05329 8.2E-12 0.5(4,100) 10.6 0.82689 7.0E-12 0.2 11.4 0.87036 1.5E-11 0.7(5,10) 7.7 0.00047 6.3E-12 0.5 8.1 0.00063 1.0E-11 0.7(5,30) 8.6 0.17388 6.0E-12 0.4 9.0 0.17493 1.4E-11 0.6Table 8:
Comparison between Regularized Newton method and QCA method on Problem 1.
Regularized Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 7.3 0.00036 4.3E-12 0.0 7.9 0.00045 7.1E-12 0.0(3,100) 4.6 0.00734 1.2E-11 0.6 7.7 0.00933 7.5E-12 4.8(3,300) 5.3 0.19312 9.8E-12 0.8 7.4 0.21646 1.1E-11 15.2(3,500) 6.0 0.97474 1.4E-11 0.9 8.9 1.18350 2.1E-11 18.0(4,10) 8.2 0.00055 7.4E-12 0.0 8.8 0.00059 7.9E-12 0.0(4,50) 4.3 0.02707 8.4E-12 0.6 6.4 0.03507 1.1E-11 3.6(4,100) 5.0 0.47484 1.3E-11 0.8 8.4 0.66361 3.1E-11 14.3(5,10) 8.9 0.00079 8.3E-12 0.0 9.3 0.00091 8.8E-12 0.0(5,30) 4.0 0.11252 1.5E-12 1.0 5.1 0.11508 1.8E-11 1.525able 9:
Comparison between Regularized Newton method and QCA method on Problem 2.
Regularized Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 5.2 0.00024 6.7E-12 0.8 6.2 0.00033 6.8E-12 3.6(3,100) 6.7 0.00894 7.9E-12 0.8 10.9 0.01227 1.1E-11 35.4(3,300) 7.2 0.25393 1.0E-11 1.0 14.3 0.41420 2.2E-12 69.1(3,500) 7.9 1.22181 1.0E-11 0.8 16.0 2.02200 1.0E-12 87.0(4,10) 4.9 0.00035 1.2E-11 0.7 5.5 0.00040 1.6E-11 3.8(4,50) 6.5 0.03804 1.2E-11 0.8 11.6 0.06193 8.4E-13 42.3(4,100) 7.0 0.61673 1.2E-11 0.8 11.8 0.96383 1.3E-12 50.2(5,10) 4.7 0.00048 1.3E-11 0.7 5.4 0.00051 7.0E-12 3.5(5,30) 6.0 0.13580 1.3E-11 0.7 10.1 0.20807 1.7E-12 31.6Table 10:
Comparison between Regularized Newton method and QCA method on Problem 3.
Regularized Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(4,10) 15.0 0.00101 6.9E-16 0 18.0 0.00132 1.6E-12 0(4,50) 16.0 0.09657 3.7E-11 0 20.0 0.09892 3.9E-12 0(4,100) 17.0 1.58356 9.0E-12 0 21.0 1.61399 8.8E-14 0Table 11:
Comparison between Regularized Newton method and QCA method on Problem 4.
Regularized Newton method QCA( m, n ) Iter Time Res Ls-iter Iter Time Res Ls-iter(3,10) 6.6 0.00029 7.4E-12 0.0 7.0 0.00036 6.7E-12 0.0(3,100) 4.8 0.00669 1.4E-11 0.6 7.3 0.00823 4.6E-12 5.8(3,300) 5.4 0.18935 1.2E-11 0.5 9.3 0.28820 5.9E-11 17.9(3,500) 5.6 0.91665 9.6E-12 0.8 9.2 1.19580 2.0E-11 20.0(4,10) 7.8 0.00048 8.1E-12 0.0 8.2 0.00052 6.9E-12 0.0(4,50) 4.4 0.02734 9.8E-12 0.6 6.5 0.03502 1.0E-11 4.0(4,100) 5.0 0.46779 1.4E-11 0.7 9.3 0.72021 3.7E-11 17.0(5,10) 8.5 0.00075 8.8E-12 0.0 8.9 0.00077 1.0E-11 0.0(5,30) 3.9 0.10700 1.6E-11 1.2 5.1 0.10740 1.2E-11 1.726able 12:
Comparison between Regularized Newton method and QCA method on Problem 5.
Regularized Newton method QCA( m, nm, n