An iterative algorithm based on strong H -tensors for identifying positive definiteness of irreducible homogeneous polynomial forms

Journal Pre-proof An iterative algorithm based on strong H-tensors for identifying positive definiteness of irreducible homogeneous polynomial forms Qilong Liu, Jianxing Zhao, Chaoqian Li, Yaotang Li

PII: DOI: Reference:

S0377-0427(19)30586-2 https://doi.org/10.1016/j.cam.2019.112581 CAM 112581

To appear in:

Journal of Computational and Applied Mathematics

Received date : 12 January 2018 Revised date : 10 October 2019 Please cite this article as: Q. Liu, J. Zhao, C. Li et al., An iterative algorithm based on strong H-tensors for identifying positive definiteness of irreducible homogeneous polynomial forms, Journal of Computational and Applied Mathematics (2019), doi: https://doi.org/10.1016/j.cam.2019.112581. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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.

© 2019 Elsevier B.V. All rights reserved.

Journal Pre-proof Manuscript Click here to view linked References

of

An iterative algorithm based on strong H-tensors for identifying positive definiteness of irreducible homogeneous polynomial forms

p ro

Qilong Liua , Jianxing Zhaob , Chaoqian Lic , Yaotang Lic,∗ a

Pr e-

School of Mathematical Sciences, Guizhou Normal University, Guiyang, 550025, People’s Republic of China b College of Data Science and Information Engineering, Guizhou Minzu University, Guiyang 550025, People’s Republic of China c School of Mathematics and Statistics, Yunnan University, Kunming 650091, People’s Republic of China

Abstract

urn

al

Identifying the positive definiteness of an even-degree homogeneous polynomial form arises in various applications, such as evolutionary game dynamics, automatical control, magnetic resonance imaging, and spectral hypergraph theory. However, it is difficult to determine whether a given homogeneous polynomial is positive definite or not because the problem is NP-hard. In this paper, an iterative algorithm of identifying the positive definiteness of irreducible homogeneous polynomial forms is proposed by identifying strong H-tensors with weakly irreducibility. The validity of the iterative algorithm is proved theoretically. Experimental results on multilinear systems, high-order Markov chains and symmetric multi-player games are presented to illustrate the applications of the proposed methods.

Jo

Keywords: Strong H-tensors, Weakly irreducible tensors, Positive definiteness, Iterative algorithm MSC: 15A18, 15A69, 65F15, 65F10

∗

Corresponding author Email address: [email protected](Qilong Liu); [email protected](Jianxing Zhao); [email protected](Chaoqian Li); [email protected] (Yaotang Li)

Preprint submitted to Elsevier

October 11, 2019

Journal Pre-proof

1. Introduction

ai1 i2 ···im ∈ R,

ij ∈ [n],

j ∈ [m],

p ro

A = (ai1 i2 ···im ),

of

An mth order n-dimensional real tensor is a multidimensional array of n elements of the form m

Pr e-

where [n] = {1, 2, · · · , n}[1]. A tensor A is called symmetric[2, 3], if its entries are invariant under any permutation of their indices, and A is called nonnegative if each entry is nonnegative. We shall denote the set of all real mth order n-dimensional tensors by R[m,n] , the set of all mth order [m,n] n-dimensional nonnegative tensors by R+ and the set of all mth order [m,n] n-dimensional real symmetric tensors by S . For an mth order n-dimensional tensor A and n-dimensional vector x = (x1 , x2 , · · · , xn )⊤ , we define n-dimensional vectors Axm−1 and x[m−1] , whose i-th component is ∑ (Axm−1 )i := aii2 ···im xi2 · · · xim , i2 ,...,im ∈[n]

and

(x[m−1] )i := xm−1 , i

urn

al

respectively. We let Axm denote the scalar defined by Axm := x⊤ Axm−1 [2, 4, 5]. Consider the following real mth degree homogeneous polynomial of n variables form ∑ f (x) = ai1 i2 ···im xi1 xi2 · · · xim , (1) i1 ,i2 ,...,im ∈[n]

Jo

where x = (x1 , x2 , · · · , xn )⊤ ∈ Rn . The form f (x) is called reducible if there exists a nonempty proper index subset Λ1 ⊆ [n] such that f (x) = f1 (xΛ1 ) + f2 (xΛ2 ),

where Λ2 = [n] \ Λ1 , f1 (xΛ1 ) and f2 (xΛ2 ) are homogeneous polynomials with the following forms ∑ fj (xΛj ) = ai1 i2 ···im xi1 xi2 · · · xim , j = 1, 2, i1 ,i2 ,...,im ∈Λj

2

Journal Pre-proof

p ro

j∈[k]

of

respectively. If f (x) is not reducible, then we call f (x) irreducible[6]. For any mth degree homogeneous polynomial f (x) of n variables, there exists a partition Λ1 , Λ2 , · · · , Λk of [n] such that f (x) can be written as ∑ f (x) = fj (xΛj ), (2) where for each j = 1, 2, · · · , k, fj (xΛj ) is irreducible or a zero polynomial[7]. The form f (x) in (1) is called positive definite if f (x) > 0,

for any x ∈ Rn , x ̸= 0.

urn

al

Pr e-

According to (2), if each irreducible homogeneous polynomial fj (xΛj ), j ∈ [k], is positive definite, then f (x) is also positive definite. Therefore, the problem of identifying positive definiteness of a homogeneous polynomial can be transformed into that of some irreducible homogeneous polynomials with reduced variables. Identifying the positive definiteness of homogeneous polynomials has applications in evolutionary game dynamics [8, 9, 10], which describes how the frequencies of strategies within a population changes in time, according to the strategies’ success. We consider here only symmetric multi-player games [11, 12, 13]. In such games, it is assumed that all players are the same role in the game and moreover the payoff of any player depends only on his strategy and on the number of players playing different types of strategies. The tensor form of symmetric multi-player games could be described by the following three parts[13, 14], (i) The set of players [m]; (ii) Each player has n available pure strategies [n]; (k)

Jo

(iii) For any k ∈ [m], let A(k) = (ai1 i2 ···im ) be the payoff tensor of player k, that is, for any ij ∈ [n] with any j ∈ [m], if player 1 plays i1 th pure strategy, player 2 plays i2 th pure strategy, · · · , player m plays im th pure strategy, then the payoffs of player 1, player 2, · · · , player m are (1) (2) (m) ai1 i2 ···im , ai1 i2 ···im , · · · , ai1 i2 ···im respectively; and for every permutation of the set of players, denoted by σ, satisfies (σ −1 (k))

(k)

ai1 i2 ···im = aiσ(1) iσ(2) ···iσ(m) ,

where σ −1 is the inverse of the permutation σ. 3

Journal Pre-proof

p ro

of

Note that the payoffs in symmetric multi-player games are uniquely determined by the first player’s payoff tensor A(1) . Consider a population consisting of n types (or strategies), and let xi be the frequency of type i. Then the state of the population is given by ∑ x ∈ ∆n := {(x1 , x2 , · · · , xn )⊤ > 0 : xi = 1}. i∈[n]

Pr e-

Assume that the xi is differentiable functions of time t and postulates a law of motion for x(t). If individuals meet randomly and engage in a symmetric multi-player game with payoff tensor A(1) , then (A(1) xm−1 )i is the expected payoff for an individual of type i and A(1) xm is the average payoff in the population state x. According to replicator equation [9, 15], the per capita xi ) i 1 rate of growth, i.e. the logarithmic derivative d(log = dx , is given by dt dt xi the difference between the payoff for type i and the average payoff in the population. Therefore, dxi = xi ((A(1) xm−1 )i − A(1) xm ), dt

i = 1, 2, · · · , n.

(3)

Jo

urn

al

Let int(∆n ) := {x ∈ ∆n : x > 0}. The system (3) is said to be permanent if there exists a compact set Λ ⊆ int(∆n ) with the property that for all x ∈ int(∆n ) there exists a T such that for all t > T , x(t) ∈ Λ. This means roughly that if initially all types (or strategies) are present in the population, then they will be, in the long run, proof against extinction through small, rare random shocks [10]. In order to describe the long-term behaviour of the dynamics, we shall discuss the stability of nonlinear system (3). As Jacobian matrix of the system (3) at an equilibrium might be singular, the stability of the equilibrium point of nonlinear system (3) can not be determined by the symbol of eigenvalues of Jacobian matrix. Therefore, its stability analysis can be reduced, using Lyapunov’s direct method [16], to extend a positive definite function, such that its time derivative along the trajectories of the system is negative. Concretely, if a real-valued function (or multivariate polynomial) V (x) can be found such that V (x) is positive definite and −V˙ (x) = −

∑ ∂V x˙i > 0, ∂xi

i∈[n]

∀x ∈ Rn , x ̸= 0,

then the system (3) is asymptotically stable. Hence, identifying whether a multivariate polynomial is positive definite or not for all real x plays an 4

Journal Pre-proof

p ro

of

important role in the stability of system (3) via Lyapunovs direct method in evolutionary game dynamics. Identifying the positive definiteness of homogeneous polynomials also has important applications in other fields, such as automatical control[17, 18, 19, 20], magnetic resonance imaging[21, 22, 23, 24], spectral hypergraph theory[25, 26, 27, 28, 29], and polynomial problems[30, 31]. For n 6 3, the positive definiteness of the homogeneous polynomial form can be checked by a method based on the Sturm theorem[17]. For n > 3 and m > 4, it is difficult to determine whether a given homogeneous polynomial f (x) is positive definite or not since the problem is NP-hard[32]. For solving this problem, the f (x) in (1) can be represented as

Pr e-

f (x) = Axm ,

(4)

Clearly, if A in (4) is not symmetric, we may replace A by a symmetric tensor A such that f (x) = Axm = Axm .

Jo

urn

al

Without loss of generality, we assume that the tensor A in (4) is symmetric, and call tensor A an underlying tensor associated with homogeneous polynomial form f (x)[33]. In [3, 34], Qi pointed out that the homogeneous polynomial form f (x) is positive definite if and only if the underlying tensor A associated with f (x) is positive definite. In [7], Zhou et al. proved that f (x) is irreducible if and only if the underlying tensor A associated with f (x) is weakly irreducible. Therefore, the positive definiteness of irreducible homogeneous polynomial form is equivalent to that of the underlying tensor with weakly irreducibility. Recently, by introducing the definition of strong H-tensors, Li et al. [35] proposed that an even-order real symmetric strong H-tensor with positive diagonal entries is positive definite. Due to this, we may identify the positive definiteness of irreducible homogeneous polynomial forms via identifying symmetric strong H-tensors with weakly irreducibility. Several criteria for strong H-tensors were given in [35, 36, 37, 38, 39, 40, 41, 42]. However, those criteria are successful only for some special cases and the conditions in them are hard to be verified. In [35], an iterative algorithm for the identification of strong H-tensors was presented and recently improved in [33, 43]. However, the algorithm is parameter-dependent, and an inappropriate choice of the parameter may result in a large amounts of computation. The iterative algorithms in [33, 35] stop in finite iterative steps 5

Journal Pre-proof

2. Preliminaries

Pr e-

p ro

of

only if the original tensor is a strong H-tensor. Some properties of strong H-tensors with weakly irreducibility was introduced in [33, 37]. Unfortunately, to the best of the authors’ knowledge, the study in identifying strong H-tensors with weakly irreducibility is very few. Therefore, finding effective criteria to identify strong H-tensors with weakly irreducibility is interesting. The remainder of this paper is distributed as follows. Some preliminaries about strong H-tensors are presented in Section 2. In Section 3, we first present a parameter-free iterative algorithm, see Algorithm 1, of identifying strong H-tensors for weakly irreducible tensors. The validity of the algorithm is proved in this section. Based on Algorithm 1 and a property of strong Htensors (see Lemma 2.2), we also give an iterative algorithm for identifying the positive definiteness of a homogeneous polynomial form in this section. In Section 4, the numerical experiments and applications are given to show the efficiency of the proposed method. Finally, we draw some conclusions in the last section.

al

We start this section with some fundamental notions and properties on tensors. The mth order n-dimensional identity tensor [44], denoted by I ∈ R[m,n] , is the tensor with entries { 1, if i1 = i2 = · · · = im , δi1 i2 ···im = 0, otherwise.

Jo

urn

Let A ∈ R[m,n] , for i ∈ [n], we call aii2 ···im for ij ∈ [n], j = 2, · · · , m, the entries of A in the ith row, where ai···i is the ith diagonal entry of A, while the other entries are the off-diagonal entries of A in the ith row [45]. We use DA to denote the 2th order n-dimensional diagonal matrix in which the diagonal entries are the same as A. We write |A| to denote the mth order n-dimensional tensor with (i1 , i2 , · · · , im )-th entry being |ai1 i2 ···im |. For a set Λ, |Λ| denotes the cardinality of Λ. Let A ∈ R[m,n] and λ ∈ C. If there is a nonzero vector x ∈ Cn such that Axm−1 = λx[m−1] ,

then λ is said to be an eigenvalue of A and x is said to be an eigenvector of A corresponding to the eigenvalue λ. If the eigenvector x is real, then the eigenvalue is also real. In this case, λ and x are called an H-eigenvalue 6

Journal Pre-proof

p ro

of

and an H-eigenvector of A, respectively. This definition was first introduced and study by Qi [3] and Lim [46], respectively. The largest modulus of the eigenvalues of A, denoted by ρ(A), is called the spectral radius of A, and the set of all eigenvalues of A is said to be the spectrum of A, denoted as σ(A). The following result gives bounds for the spectral radius of a nonnegative tensor. [m,n]

Lemma 2.1. [44] Let A ∈ R+ . Then ∑ min aii2 ···im 6 ρ(A) 6 max i∈[n]

i∈[n]

i2 ,...,im ∈[n]

∑

i2 ,...,im ∈[n]

aii2 ···im .

Pr e-

We now recall the definitions of strong M-tensors, comparison tensors and strong H-tensors. Definition 2.1. [47] A tensor A = (ai1 i2 ···im ) ∈ R[m,n] is called an M-tensor [m,n] if there exists a nonnegative tensor B ∈ R+ and a positive real number η > ρ(B) such that A=ηI − B. If η > ρ(B), then A is called a strong M-tensor.

al

Definition 2.2. [48] Let A = (ai1 i2 ···im ) ∈ R[m,n] . We say that M(A) = (mi1 i2 ···im ) ∈ R[m,n] is the comparison tensor of A if { |ai1 i2 ···im |, if i1 = i2 = · · · = im , mi1 i2 ···im = −|ai1 i2 ···im |, otherwise.

urn

Definition 2.3. [48] Let A = (ai1 i2 ···im ) ∈ R[m,n] . If M(A) is an M-tensor, then A is called an H-tensor. If M(A) is a strong M-tensor, then A is called a strong H-tensor. Moreover, Li et al.[35] also provided the following definition of strong H-tensors, which is equivalent to the Definition 2.3.

Jo

Definition 2.4. [35] A tensor A = (ai1 i2 ···im ) ∈ R[m,n] is called a strong H-tensor if there is an entrywise positive vector x = (x1 , x2 , · · · , xn )⊤ ∈ Rn such that for all i ∈ [n], ∑ |aii2 ···im |xi2 · · · xim . |ai···i |xm−1 > i i2 ,i3 ,...,im ∈[n], δii ···im =0 2

7

Journal Pre-proof

of

The strong H-tensor is closely linked with the positive definiteness of symmetric tensors. Lemma 2.2. [35] Any even order symmetric strong H-tensor with positive diagonal entries is positive definite.

p ro

Remark here that a strong H-tensor with positive diagonal entries is called an H+ -tensor [41]. Therefore, an even order symmetric H+ -tensor is positive definite. The following multiplication of tensors was given in [49].

Pr e-

Definition 2.5. [49] Let A = (ai1 i2 ···im ) ∈ R[m,n] , B = (bi1 i2 ···ik ) ∈ R[k,n] with m > 2, k > 1. Defined the product AB to be the following tensor C of order (m − 1)(k − 1) + 1 and dimension n ∑ ciα1 α2 ···αm−1 = aii2 ···im bi2 α1 · · · bim αm−1 , i2 ,...,im ∈[n]

where i ∈ [n], α1 , · · · , αm−1 ∈ {i2 · · · ik : ij ∈ [n], j = 2, · · · , k}. It was proved in [49] that the product of tensors in Definition 2.5 satisfies the following distributive law and associative law.

al

Lemma 2.3. [49] Let A ∈ R[m+1,n] , B ∈ R[k+1,n] and C ∈ R[l+1,n] . Then the tensor product defined in Definition 2.5 has the following properties:

urn

(i) (The left distributivity): (A + B)C=AC + BC (when A and B have the same order, i.e., m = k); (ii) (The right distributivity in the case of m = 1): A(B + C)=AB + AC (when A is a matrix, B and C have the same order, i.e., k = l. Note that in general when A is not a matrix, then the right distributivity does not hold.);

Jo

(iii) (The associativity): A(BC) = (AB)C. Remark 2.1. Given X, Y ∈ R[2,n] and A, B ∈ R[m,n] . By Lemma 2.3, we have XAY + XBY = X(A + B)Y. Based on the product of tensors, Kannan, et al. established a necessary and sufficient condition for a tensor to be a strong H-tensor in [37]. 8

Journal Pre-proof

−1 ρ(I − DM(A) M(A)) < 1.

of

Lemma 2.4. [37] Let A = (ai1 ···im ) ∈ R[m,n] . Then, A is a strong H-tensor if and only if ai···i ̸= 0 for all i ∈ [n] and

p ro

Lemma 2.5. [37] Let A = (ai1 ···im ) ∈ R[m,n] with ai···i ̸= 0 for all i ∈ [n]. Then, A is an H-tensor if and only if −1 ρ(I − DM(A) M(A)) 6 1.

The following definitions of diagonal similarity and permutation similarity for tensors were proposed by Shao[49].

Pr e-

Definition 2.6. [49] Let A, B ∈ R[m,n] . We say that A and B are diagonal similar, if there exists some invertible diagonal matrix X = diag(x1 , x2 , · · · , xn ) such that B = X 1−m AX, where X 1−m = diag(x1−m , x1−m , · · · , x1−m ). 1 2 n Definition 2.7. [49] Let A, B ∈ R[m,n] . We say that A and B are permutational similar, if there exists some permutation matrix P of dimension n such that B = P AP ⊤ .

al

Lemma 2.6. [49] Let A, B ∈ R[m,n] . If A and B are diagonal similar (or permutational similar, respectively), then σ(A) = σ(B). [m,n]

Jo

urn

Let A = (ai1 ···im ) ∈ R+ , it is associated to a directed graph G(A) = (V, E(A)), where V = [n] and a directed arc (i, j) ∈ E(A) if there exist indices {i2 , · · · , im } such that j ∈ {i2 , · · · , im } and aii2 ···im ̸= 0. A walk in G(A) = (V, E(A)) from vertex i to vertex j is a sequence γ of vertices i = i0 , i1 , · · · , ik = j such that (il , il+1 ) ∈ E(A) is a directed arc for each l = 0, 1, · · · , k−1. If the vertices i0 , i1 , · · · , ik are distinct, then γ is a directed path joining i and j. A directed graph is strongly connected provided that for each pair i and j of distinct vertices, there exists a directed path from i to j. In [50], Friedland et al. introduced the notion of weakly irreducible tensors as follows. [m,n]

Definition 2.8. [50] A nonnegative tensor A ∈ R+ is called weakly irreducible if the associated directed graph G(A) is strongly connected. A tensor A ∈ R[m,n] is said to be weakly irreducible if |A| is weakly irreducible. 9

Journal Pre-proof

of

In [7], Zhou et al. established the relation between the irreducibility of a homogeneous polynomial form f (x) and the weakly irreducibility of underlying tensor associated with f (x).

p ro

Lemma 2.7. [7] Let A ∈ R[m,n] and f (x) = Axm . Then, A is weakly irreducible if and only if f (x) is irreducible. For a tensor A = (ai1 ···im ) ∈ R[m,n] , let ∑ ∑ ri (A) := |aii2 ···im | =

i2 ,...,im ∈[n]

i2 ,...,im ∈[n], δii ···im =0 2

and

R(A) := max ri (A).

Pr e-

r(A) := min ri (A),

|aii2 ···im | − |ai···i |,

i∈[n]

i∈[n]

3. An algorithm of identifying strong H-tensors and an algorithm of identifying positive definiteness for homogeneous polynomial forms

urn

al

In this section, we first present an iterative algorithm of identifying strong H-tensors with weakly irreducibility and prove that this algorithm is convergent for any weakly irreducible tensors. Based on the equivalence of the positive definiteness of the irreducible homogeneous polynomial forms and the positive definiteness of the underlying irreducible tensors, and a property of strong H-tensors (see Lemma 2.2), we propose an iterative algorithm of identifying the positive definiteness for irreducible homogeneous polynomial forms.

Jo

3.1. An algorithm of identifying strong H-tensors for weakly irreducible tensors In this subsection, we first provide a parameter-free iterative algorithm, see Algorithm 1, of identifying whether a weakly irreducible tensor is a strong H-tensor or not, and then we prove that the algorithm always converges. Remark 3.1. If m = 2, then Algorithm 1 reduces to Algorithm AH of [51]. Denote Pn := {x ∈ Rn : xi > 0, i ∈ [n]} and int(Pn ) := {x ∈ Rn : xi > 0, i ∈ [n]}. For any nonnegative vector x ∈ Pn , we define ϕ : Pn → P1 by ∑ ϕ(x) := xi . i∈[n]

10

Journal Pre-proof

Algorithm 1 Identifying strong H-tensors for weakly irreducible tensors

−1 A, X (0) := I(I is an n-by-n identity matrix), q = 1, Set X := I, A(0) := DA (q) Compute X = XX (q−1) , A(q) = (X (q−1) )1−m A(q−1) X (q−1) = (ai1 ···im ), Compute ri (A(q) ), r(A(q) ) and R(A(q) ), If r(A(q) ) > 1, then output “A is not a SHT”, stop; otherwise, If R(A(q) ) < 1, then output “A is a SHT”, stop; otherwise, If r(A(q) ) = R(A(q) ), then output “ρ(|A(0) |) is approximate to 2”, stop; otherwise, (q) (q) (q) Set X (q) := diag(x1 , x2 , · · · , xn ), where

Pr e-

1: 2: 3: 4: 5: 6: 7:

p ro

X = X (0) X (1) · · · X (q−1) .

of

Input: A weakly irreducible tensor A = (ai1 ···im ) ∈ R[m,n] with ai···i ̸= 0 for all i ∈ [n]. Output: “A is not a SHT” if A is not a strong H-tensor or “A is a SHT” if A is a strong H-tensor, and a positive diagonal matrix

(q) xi

8: Set q := q + 1; go to step 2.

[

1 + ri (A(q) ) = 1 + R(A(q) )

1 ] m−1

,

al

In order to prove the convergence of Algorithm 1, let us recall Power Type Algorithm, see Algorithm 2, proposed by Zhou, Qi and Wu[52] for calculating the largest eigenvalue of weakly irreducible nonnegative tensors. It is shown in [52] that the obtained sequences {λq } and {λq } in Algorithm 2 converge to ρ(B). Furthermore, ρ(A) = ρ(B) − 1. (q)

urn

Lemma 3.1. Let A(q) = (ai1 i2 ···im ) be the sequence produced by Algorithm 1. Then (i) A(1) = A(0) ; (q)

(ii) aii···i = 1, for all i ∈ [n], q = 0, 1, · · · ;

Jo

(iii)

∏

[1 + ri (A(k) )] =

k∈[q]

∑

i2 ,...,im ∈[n]

(1)

|aii2 ···im |

for all i ∈ [n], q = 2, 3, · · · .

11

∏

k∈[q−1]

1

1

[1 + ri2 (A(k) )] m−1 · · · [1 + rim (A(k) )] m−1 . (5)

Journal Pre-proof

Algorithm 2 Power Type Algorithm [52] [m,n]

and x(0) ∈

and

(q)

λq = min (q)

xi >0

p ro

of

Input: A weakly irreducible nonnegative tensor A = (ai1 ···im ) ∈ R+ int(Pn ). Output: A scalar λ. 1: Set B := A + I and q := 0, 2: Compute y(q) = B(x(q) )m−1 , (q)

yi (q)

(xi )m−1

,

λq = max (q)

xi >0

yi

(q)

(xi )m−1

,

3: If λq = λq , then let λ := λq and stop, otherwise, compute 1

Pr e-

(y(q) )[ m−1 ]

x(q+1) =

4: Set q := q + 1, go to step 2.

1

ϕ((y(q) )[ m−1 ] )

,

Proof. It follows from Step 2 of Algorithm 1 that

A(1) = (X (0) )1−m A(0) X (0) = I 1−m A(0) I = A(0) ,

al

−1 which implies (i) holds. Since A(0) = DA A, all diagonal entries of A(0) and (1) A are 1. For q = 2, 3, · · · ,

urn

A(q) = (X (q−1) )1−m · · · (X (1) )1−m A(1) X (1) · · · X (q−1) , which implies (q) ai1 i2 ···im

=

Jo

=

(1) ai1 i2 ···im

(1) ai1 i2 ···im

(k) ∏ x(k) i2 · · · xim

k∈[q−1]

(k)

(xi1 )m−1

1 1 ∏ [1 + ri (A(k) )] m−1 · · · [1 + rim (A(k) )] m−1 2 . 1 + ri1 (A(k) )

k∈[q−1]

Therefore, all diagonal entries of A(q) are 1. Hence, (ii) holds. The ith row sum of A(q) are (q)

1 + ri (A ) =

∑

i2 ,...,im ∈[n]

(1) |aii2 ···im |

1 1 ∏ [1 + ri (A(k) )] m−1 · · · [1 + rim (A(k) )] m−1 2 . 1 + ri (A(k) )

k∈[q−1]

12

Journal Pre-proof

∏

Multiplying

(1 + ri (A(k) )) on the both sides of the above equality yields

∏

[1 + ri (A(k) )] =

∑

(1)

i2 ,...,im ∈[n]

∏

1

1

[1 + ri2 (A(k) )] m−1 · · · [1 + rim (A(k) )] m−1 ,

k∈[q−1]

p ro

k∈[q]

|aii2 ···im |

of

k∈[q−1]

which implies (iii) holds. The proof is completed.

By Lemma 3.1 and Power Type Algorithm (Algorithm 2)[52], we prove that Algorithm 1 always converges except in a special case, and its output is correct.

Pr e-

Theorem 3.1. Let A = (ai1 ···im ) ∈ R[m,n] be weakly irreducible with ai···i ̸= 0 for all i ∈ [n]. Then Algorithm 1 always converges except ρ(|A(0) |) = 2. −1 (q) Proof. Let B (q) := I − DM(A ). According to Lemma 3.1, it follows (q) ) M(A that DM(A(q) ) = I, thus

B (q) = I − M(A(q) ) = |A(q) | − I, which implies

(6)

al

|A(q) | = B (q) + I.

For each q = 1, 2, · · · , we have

|A(q+1) | − I |(X (q) )1−m A(q) X (q) | − (X (q) )1−m IX (q) (X (q) )1−m (|A(q) | − I)X (q) (X (q) )1−m B(q) X (q) .

urn

B(q+1) = = = =

Jo

By Definition 2.6, we have B(q+1) and B (q) are diagonal similar, and inductively so are B(q) and B(1) . From Lemma 2.6, we obtain ρ(B (1) ) = ρ(B (2) ) = · · · = ρ(B (q) ) = · · · .

Since A is weakly irreducible so are A(1) and B (1) . Because the diagonal similarity transformation (Step 2) preserves weakly irreducibility, by induction, we have B(q) is nonnegative weakly irreducible. 13

Journal Pre-proof

Let d(1) := (1, 1, · · · , 1)⊤ ∈ int(Pn ). For q = 1, 2, · · · . Compute 1

(q) m−1

= |A |(d )

(q+1)

,

d

and (q)

(q)

(di )m−1

(q)

di >0

1

ϕ((y(q) )[ m−1 ] )

,

(q)

yi

λq = min

=

(y(q) )[ m−1 ]

of

(1)

p ro

y

(q)

,

λq = max (q)

di >0

(1)

(1)

yi

(q)

(di )m−1

.

Pr e-

It is easy to obtain di = 1 and yi = 1 + ri (A(1) ) for i ∈ [n]. Next we will prove that for i ∈ [n] and q = 2, 3, · · · , (q) di

1 ∏ [1 + ri (A(k) )] m−1

=

1

ϕ((y(k) )[ m−1 ] )

k∈[q−1]

and (q)

yi

∏

,

(7)

.

(8)

[1 + ri (A(k) )]

k∈[q]

=

∏

1

[ϕ((y(k) )[ m−1 ] )]m−1

al

k∈[q−1]

We proceed by induction on q. For i ∈ [n] and q = 2, we have (1)

=

1

(yi ) m−1

urn

(2) di

1

ϕ((y(1) )[ m−1 ] )

1

=

[1 + ri (A(1) )] m−1 1

ϕ((y(1) )[ m−1 ] )

.

(9)

By (5), we have (2)

= (|A(1) |(d(2) )m−1 )i ∑ (1) (2) (2) = |aii2 ···im |di2 · · · dim i2 ,...,im ∈[n]

Jo

yi

=

∑

i2 ,...,im ∈[n]

=

∏

[1 (1) |aii2 ···im |

1

[ϕ((y(1) )[ m−1 ] )]m−1

1

1

[ϕ((y(1) )[ m−1 ] )]m−1

[1 + ri (A(k) )]

k∈[2]

1

+ ri2 (A(1) )] m−1 · · · [1 + rim (A(1) )] m−1

, 14

Journal Pre-proof

Pr e-

p ro

of

which together with (9) shows that (7) and (8) hold for i ∈ [n] and q = 2. Assume now that (7) and (8) hold for i ∈ [n] and q = 2, 3, · · · , l. For i ∈ [n], q = l + 1, by the induction hypothesis, we have ( ) 1 (y(l) )[ m−1 ] (l+1) di = 1 ϕ((y(l) )[ m−1 ] ) i 1 1 = [(y(l) )[ m−1 ] ]i 1 [ ] ϕ((y(l) ) m−1 ) ∏ 1 [1 + ri (A(k) )] m−1 1 k∈[l] = · ∏ 1 1 [ ] ϕ((y(l) ) m−1 ) [ϕ((y(k) )[ m−1 ] )] k∈[l−1]

=

1 ∏ [1 + ri (A(k) )] m−1 1

k∈[l]

[ϕ((y(k) )[ m−1 ] )]

,

which, together with (5), yields (l+1)

yi

= (|A(1) |(d(l+1) )m−1 )i ∑ (1) (l+1) (l+1) |aii2 ···im |di2 · · · dim = i2 ,...,im ∈[n]

(1)

∏

al

∑

i2 ,...,im ∈[n]

urn

=

|aii2 ···im |

∏

k∈[l]

∏

∏

1

[ϕ((y(k) )[ m−1 ] )]m−1

k∈[l]

[1 + ri (A(k) )]

k∈[l+1]

=

1

1

[ϕ((y(k) )[ m−1 ] )]m−1

,

k∈[l]

Jo

which means (7) and (8) hold for i ∈ [n] and q = l + 1. Dividing the equalities (7) with (8), then (q)

yi

(q) (di )m−1

= 1 + ri (A(q) ), ∀i ∈ [n], q = 1, 2, · · · ,

which implies

λq = 1 + r(A(q) ),

1

[1 + ri2 (A(k) )] m−1 · · · [1 + rim (A(k) )] m−1

λq = 1 + R(A(q) ), q = 1, 2, · · · . 15

Journal Pre-proof

lim λq = lim λq = ρ(|A(1) |).

q→∞

q→∞

Therefore,

of

According to the Power Type Algorithm (Algorithm 2), it follows that

p ro

lim r(A(q) ) = lim R(A(q) ) = ρ(|A(1) |) − 1 = ρ(B(1) ).

q→∞

q→∞

Consequently, for any ε > 0, there exists a positive integer number K such that q > K implies r(A(q) ), R(A(q) ) ∈ (ρ(B (1) ) − ε, ρ(B (1) ) + ε).

(10)

Pr e-

To prove our claim we have to distinguish the following three cases. Case i If ρ(B(1) ) < 1, then let ε := 1 − ρ(B (1) ). By (10), there exists a positive integer number K1 such that q > K1 implies r(A(q) ), R(A(q) ) ∈ (2ρ(B (1) ) − 1, 1). Hence, for all q > K1 , it will be R(A(q) ) < 1. Case ii If ρ(B(1) ) > 1, then let ε := ρ(B (1) ) − 1. By (10), there exists a positive integer number K2 such that q > K2 implies

al

r(A(q) ), R(A(q) ) ∈ (1, 2ρ(B(1) ) − 1).

urn

Hence, for all q > K2 , it will be r(A(q) ) > 1. Case iii If ρ(B(1) ) = 1, then by (10), there exists a positive integer number K3 such that, for all q > K3 , r(A(q) ) and R(A(q) ) may be very close to 1. The algorithm, theoretically, may converge or not. Computationally, it is inconclusive. According to (6) and the statement (i) of Lemma 3.1, it follows that ρ(B (1) ) = 1 if and only if ρ(|A(0) |) = 2. Therefore, Algorithm 1 always converges except ρ(|A(0) |) = 2.

Jo

Theorem 3.2. Let A ∈ R[m,n] be weakly irreducible with ai···i ̸= 0 for all i ∈ [n]. If Algorithm 1 converges, then its output is correct.

Proof. Assume that Algorithm 1 converges. Then, we have to distinguish cases based on the algorithm’s output as follows. Case i Suppose that the output is “A is not a SHT”. In this case, we have r(A(q) ) > 1. By the proof of Theorem 3.1, we have ρ(B (1) ) = ρ(B (q) ). 16

Journal Pre-proof

i∈[n]

i2 ,...,im ∈[n]

∑

max i∈[n]

i2 ,...,im ∈[n]

p ro

By Lemma 2.1, it follows that

(q)

|bii2 ···im | = R(A(q) ).

of

Since B (q) = |A(q) | − I, we obtain ∑ (q) min |bii2 ···im | = r(A(q) ),

r(A(q) ) 6 ρ(B(q) ) 6 R(A(q) ),

which, together with ρ(B(1) ) = ρ(B(q) ) and r(A(q) ) > 1, yields ρ(B(1) ) > 1.

urn

al

Pr e-

Thus, by Lemma 2.4, A(1) is not a strong H-tensor and so A is not a strong H-tensor. Case ii Suppose that the output is “A is a SHT”. In this case, we have R(A(q) ) < 1. Then by an argument similar to the one used in Case i previously shows that ρ(B (1) ) < 1. Hence, by Lemma 2.4, A is a strong H-tensor. Case iii Suppose that the output is “ρ(|A(0) |) is approximate to 2”, then since Algorithm 1 has just passed from Steps 4 and 5 it is implied that r(A(q) ) 6 1 and R(A(q) ) > 1, respectively. By Step 6 of Algorithm 1, r(A(q) ) = R(A(q) ) = 1, hence by Theorem 3.1, ρ(B (0) ) = ρ(|A(0) | − I) = 1. By virtue of Lemmas 2.4 and 2.5, it follows that A is a H-tensor but not a strong H-tensor. If the equality in Step 6 has been produced computationally, then ρ(|A(0) |) is approximate to 2 and therefore no conclusion whether A is a strong A or not can be draw. Remark 3.2. Each iteration of Algorithm 1, except for initializations, substitutions and comparisons, consists of the following three major steps (q)

(q−1)

(q−1)

(q−1)

Jo

(i) ai1 ···im = ai1 ···im (xi1 )1−m xi2 ∑ (q) (ii) ri (A(q) ) = |aii2 ···im |;

(q−1)

· · · xim

;

i2 ,...,im ∈[n], δii ···im =0 2

(q)

(iii) xi =

[

1+ri (A(q) ) 1+R(A(q) )

1 ] m−1

.

Therefore, Algorithm 1 requires (m+1)nm +3n additions and multiplications per iterations step. 17

Journal Pre-proof

p ro

of

To conclude this section, we present a theorem which guarantees that after at most l(l 6 ⌊ n2 ⌋) iterations it will be [r(A(q+l) ), R(A(q+l) )] ⊂ [r(A(q) ), R(A(q) )], while after at most s(l 6 s 6 n − 1) iterations there will hold r(A(q) ) < r(A(q+s) ) 6 ρ(B(1) ) 6 R(A(q+s) ) < R(A(q) ). Here, ⌊x⌋ is the floor function, that is, it equals the largest integer not greater than x. This theorem is based on the following lemma. Lemma 3.2. Let A(q) be the sequence produced by Algorithm 1. Define S (q) := {i : ri (A(q) ) = R(A(q) )}, and

T (q) := {i : ri (A(q) ) = r(A(q) )}

Pr e-

V (q) := {i : r(A(q) ) < ri (A(q) ) < R(A(q) )}.

Then

(i) for all i ∈ [n], q = 1, 2, · · · , there holds

r(A(q) ) 6 ri (A(q+1) ) 6 R(A(q) );

(11)

(ii) in each iteration there exist at least two indices i and j such that the above inequalities turn strict. Hence, |T (q+1) | 6 |T (q) | − 1

al

|S (q+1) | 6 |S (q) | − 1,

and |V (q+1) | > |V (q) | + 2.

Proof. (i) According to Step 2 of Algorithm 1, it follows that

urn

A(q+1) = (X (q) )1−m A(q) X (q) ,

q = 1, 2, · · · .

Therefore, for i ∈ [n], we have ∑ (q+1) ri (A(q+1) ) = |aii2 ···im | i2 ,...,im ∈[n], δii ···im =0 2

Jo

=

=

∑

|aii2 ···im |(xi )1−m xi2 · · · xim

∑

[1 (q) |aii2 ···im |

i2 ,...,im ∈[n], δii ···im =0 2

i2 ,...,im ∈[n], δii ···im =0 2

(q)

(q)

(q)

(q)

1

1

+ ri2 (A(q) )] m−1 · · · [1 + rim (A(q) )] m−1 . 1 + ri (A(q) ) (12) 18

Journal Pre-proof

Hence, ∑

(q)

i2 ,...,im ∈[n], δii ···im =0 2

|aii2 ···im |

1 + R(A(q) ) 1 + ri (A(q) )

of

ri (A(q+1) ) 6

p ro

ri (A(q) ) + ri (A(q) )R(A(q) ) = 1 + ri (A(q) ) R(A(q) ) + ri (A(q) )R(A(q) ) 6 1 + ri (A(q) ) = R(A(q) ),

Pr e-

and

ri (A(q+1) ) >

∑

(q)

i2 ,...,im ∈[n], δii ···im =0 2

|aii2 ···im |

1 + r(A(q) ) 1 + ri (A(q) )

ri (A(q) ) + ri (A(q) )r(A(q) ) 1 + ri (A(q) ) r(A(q) ) + ri (A(q) )r(A(q) ) > 1 + ri (A(q) ) = r(A(q) ),

al

=

urn

which implies (11) holds. (ii) Assume that

(q)

|S (q) | = n1 , (q)

(q)

|V (q) | = n2 , (q)

(q)

and |T (q) | = n3 .

(q)

Jo

It is easy to obtain n1 + n2 + n3 = n. Suppose also that by a similarity (q) permutation with permutation matrix P (q) , we bring the aforementioned n1 (q) (q) rows first, the n2 rows next, and the n3 rows last, by keeping the same symbol A(q) for P (q) A(q) (P (q) )⊤ . So, the following three series of relationships hold (q)

ri (A(q) ) = R(A(q) ),

i = 1, · · · , n1 , (q)

r(A(q) ) < ri (A(q) ) < R(A(q) ),

(q)

(q)

i = n1 + 1, · · · , n1 + n2 , (q)

ri (A(q) ) = r(A(q) ),

(q)

i = n1 + n2 + 1, · · · , n. 19

Journal Pre-proof

(q)

of

Next, we will prove that there exists at least one index i ∈ {1, · · · , n1 } such that ri (A(q+1) ) < R(A(q) ).

∑

ri (A(q+1) ) <

(q)

|aii2 ···im |

i2 ,...,im ∈[n], δii ···im =0 2

p ro

For q ∈ [n], A(q) is weakly irreducible since A is weakly irreducible. Thus, (q) (q) there exists k ∈ {n1 + 1, · · · , n} such that ai···k···im ̸= 0 for some i ∈ (q) {1, · · · , n1 }, which together with (12), yields 1 + R(A(q) ) = R(A(q) ). 1 + ri (A(q) ) (q)

Pr e-

A similar argument shows that there exists at least one index j ∈ {n1 + (q) n2 + 1, · · · , n} such that rj (A(q+1) ) > r(A(q) ).

(q)

(q)

(q)

Finally, for i ∈ {n1 + 1, · · · , n1 + n2 }, we have ∑

ri (A(q+1) ) 6

and

)>

∑

urn

ri (A

(q+1)

al

i2 ,...,im ∈[n], δii ···im =0 2

(q)

|aii2 ···im |

i2 ,...,im ∈[n], δii ···im =0 2

1 + R(A(q) ) < R(A(q) ), 1 + ri (A(q) )

(q) |aii2 ···im |

1 + r(A(q) ) > r(A(q) ). (q) 1 + ri (A )

The proof is completed.

Jo

Theorem 3.3. Let A(q) be the sequence produced by Algorithm 1. For any given iteration step q there exist two numbers l(l 6 ⌊ n2 ⌋) and s (l 6 s 6 n−1) such that [r(A(q+k) ), R(A(q+k) )] ⊂ [r(A(q) ), R(A(q) )],

k = l, s;

(13)

r(A(q) ) < r(A(q+s) ) 6 ρ(B (1) ) 6 R(A(q+s) ) < R(A(q) ).

(14)

especially for k = s,

20

Journal Pre-proof

p ro

of

Proof. We need to prove that either r(A(q) ) < r(A(q+k) ) or R(A(q+k) ) < R(A(q) ) for k = l, s. Using the notation of Lemma 3.2, it is obvious that the (q) (q) (q) “worst” case is n2 = 0, and (i) for n even n1 = n3 = n2 , while (ii) for n (q) (q) (q) (q) odd either n1 = n+1 and n3 = n−1 or n1 = n−1 and n3 = n+1 . It follows 2 2 2 2 from Lemma 3.2 that in either case, n is even or odd, the maximum number of iterations needed to have one of the two inequalities strict is l = ⌊ n2 ⌋. Finally, for the set of strict inequalities in (14) to be satisfied the “worst” (q) (q) (q) case is again to have n2 = 0 and either n1 = n − 1 or n3 = n − 1. So, combining it with the result of Lemma 3.2 we have Inequalities (14) are both satisfied after at most n − 1 iterations.

Pr e-

3.2. An iterative algorithm of identifying positive definiteness for irreducible homogeneous polynomial forms In this subsection, based on Lemmas 2.2, 2.7 and Algorithm 1, we present an iterative algorithm of identifying the positive definiteness for irreducible homogeneous polynomial forms. Algorithm 3 Identifying positive definiteness for irreducible homogeneous polynomial forms

urn

al

Input: An irreducible homogeneous polynomial f (x) = Axm with ai···i > 0 for all i ∈ [n]. Output: “f (x) is positive definite” or “the positive definiteness of f (x) cannot be determined”. −1 1: Set X := I, A(0) := DA A, X (0) := I, q := 1, (q) 2: Compute X = XX (q−1) , A(q) = (X (q−1) )1−m A(q−1) X (q−1) = (ai1 ···im ), 3: Compute ri (A(q) ), r(A(q) ) and R(A(q) ), 4: If R(A(q) ) < 1, then output “f (x) is positive definite”, stop; otherwise, 5: If r(A(q) ) > 1 or r(A(q) ) = R(A(q) ), then output “the positive definiteness of f (x) cannot be determined”, stop; otherwise, (q) (q) (q) 6: Set X (q) := diag(x1 , x2 , · · · , xn ), where (q) xi

[

1 + ri (A(q) ) = 1 + R(A(q) )

1 ] m−1

,

Jo

7: Set q := q + 1; go to step 2.

4. Numerical experiments and applications In this section, we present numerical examples to demonstrate the theoretical results in Section 3. The numerical experiments will be performed via 21

Journal Pre-proof

p ro

Example 4.1. Let A = (ai1 i2 i3 i4 ) ∈ S[4,4] , where

of

MATLAB R2015a (version: 8.5.0.197613) on a 2.7 GHz Inter computer with 8GB of RAM. In the following, we first report two preliminary numerical results for Algorithms 1 and 3.

a1111 = 4, a2222 = a, a3333 = b, a4444 = 1, a1114 = a1141 = a1411 = a4111 = −1, a1333 = a3133 = a3313 = a3331 = −1, a2223 = a2232 = a2322 = a3222 = −1, a2333 = a3233 = a3323 = a3332 = −2,

Pr e-

and the remaining ai1 i2 i3 i4 = 0, where a, b ∈ R take different values. We use Algorithm 1 to identify whether A is a strong H-tensor or not. The numerical results are showed in Table 1, where q denotes the number of iteration steps, CPU(s) denotes the CPU time in seconds, X denotes the output positively diagonal matrix when the algorithm terminates. Table 1: Numerical results of Algorithm 1 for Example 4.1

b 12 9 10 12 12

q 9 3 1 5 3

CPU(s) 0.3725 0.0177 0.0121 0.0320 0.0177

Output A is not a SHT A is not a SHT (0) ρ(|A |) is approximate to 2 A is a SHT A is a SHT

al

a 3 4 5 4 5

X diag(0.7478,1.0000,0.7720,0.7452) diag(0.9110,1.0000,0.9608,0.9059) diag(1.0000,1.0000,1.0000,1.0000) diag(0.9207,1.0000,0.8707,0.9295) diag(0.9930,0.9806,0.9310,1.0000)

Jo

urn

Example 4.2. Consider the following 6th degree irreducible homogeneous polynomial f (x) = Ax6 where A = (ai1 ···i6 ) ∈ S[6,10] is a real weakly irreducible symmetric tensor, and it is randomly constructed as follows. Ran(i) (i) domly generate three 10-dimensional vectors, denoted by v(i) = (v1 , · · · , v10 )⊤ , i = 1, 2, 3, such that the elements of each vector are generated by uniform distribution [−1, 1]. Let ∑ B= (v(i) )6 , i∈{1,2,3}

where (v(i) )6 , i = 1, 2, 3, are 6th order 10-dimensional rank-one tensors whose entries are (i) (i)

(i)

(v(i) )6i1 i2 ···i6 = vi1 vi2 · · · vi6 ,

i1 , i2 , · · · , i6 ∈ {1, 2, · · · , 10}. 22

Journal Pre-proof

Pr e-

p ro

of

Furthermore, we construct A by taking absolute value for each diagonal entries of B and amplify each diagonal entries 3.5 × 104 times. Therefore A is a 6th order 10-dimensional symmetric tensor with weakly irreducibility, and f (x) is a 6th degree irreducible homogeneous polynomial. Using the above method construct 100 irreducible homogeneous polynomial forms with 6th degree. We identify the positive definiteness of homogeneous polynomial f (x) by using Algorithm 3. The numerical results are showed in Figure 1. The x-axis of Fig.1 refers to the k-th random generated homogeneous polynomial form, and the y-axis refers to the number of iteration steps of Algorithm 3. The plus symbol in blue color denotes f (x) is positive definite. The star symbol in green color denotes the positive definiteness of f (x) cannot be determined. In the 100 generated homogeneous polynomial forms, 75 cases are positive definiteness, the other 25 cases about the positive definiteness of f (x) cannot be determined. 13

f(x) is positive definite

12

the positive definiteness of f(x) cannot be determined

10 9

al

8 7 6 5 4 3 2 1 0 0

urn

The number of iteration steps

11

10

20

30

40

50

60

70

80

90

100

The k-th random generated homogeneous polynomial f(x)

Jo

Figure 1: The randomly generated results of Algorithm 3 for Example 4.2

To illustrate the applications of Algorithms 1 and 3, we give the following three numerical examples, which arise from the study of the multilinear systems[53], the higher-order Markov chains[54] and the symmetric multiplayer games[12]. 23

Journal Pre-proof

Axm−1 = b.

of

4.1. Applications to multilinear systems Let A ∈ R[m,n] and b ∈ Rn . A multilinear system[53] can be expressed as (15)

Pr e-

p ro

Solving the multilinear system (15) is always an important issue in engineering and scientific computing [55, 56]. An interesting result proved by Wang, Che and Wei[57] is that a system of H+ -tensor equations with the positive right-hand side always has a unique positive solution. Therefore, we can use Algorithm 1 to identify the existence and uniqueness of positive solution for the multilinear system (15). In the following, we consider a real application tensor, which arises from the study of the numerical differential equations [55]. Example 4.3. Consider an ordinary differential equation of particle’s movement under the gravitation d2 x(t) GM =− 2 , 2 dt x (t)

in (0, 1)

with Dirichlet’s boundary conditions

x(1) = c1 ,

al

x(0) = c0 ,

urn

where G ≈ 6.67 × 10−11 N m2 /kg 2 is the gravitational constant and M ≈ 5.98 × 1024 kg is the mass of earth. After the discretization, we obtain a system of polynomial equations  x31 = c30 ,  GM 3 2 2 2x − xi xi−1 − xi xi+1 = (n−1) i = 2, 3, · · · , n − 1, 2,  i 3 3 x n = c1 ,

Jo

which can be rewritten into a multilinear equation Ax3 = b,

where A = (aijkl ) is a 4-order n dimensional tensor with elements  a1111 = annnn = 1,    aiiii = 2, i = 2, 3, · · · , n − 1, ai(i−1)ii = aii(i−1)i = aiii(i−1) = − 13 , i = 2, 3, · · · , n − 1,    ai(i+1)ii = aii(i+1)i = aiii(i+1) = − 13 , i = 2, 3, · · · , n − 1, 24

(16)

Journal Pre-proof

p ro

of

and the right-hand side is a vector with elements  b1 = c20 ,  GM bi = (n−1) i = 2, 3, · · · , n − 1, 2,  2 b n = c1 .

Pr e-

We use Algorithm 1 to identify whether A is a strong H-tensor or not. The numerical results are showed in Table 2, where n denotes the dimensional of the tensor A, and q, CPU(s) and X are the same as those used in Table 1. In Table 2, we can see that A is a strong H-tensor for a fix dimensional n. It is easy to see that the diagonal entries of A are positive and the right-hand side vector b is positive. Hence the mulitilinear system (16) has a unique positive solution. Table 2: Numerical results of Algorithm 1 for Example 4.3

q 2 2 3

CPU(s) Output 0.0108 A is a SHT 0.0128 A is a SHT 0.0236 A is a SHT

6

3

0.0253

A is a SHT

7

4

0.0368

A is a SHT

8

4

0.0442

9

5

A is a SHT

urn

10 5

X diag(0.6532,1.0000,0.6532) diag(0.6412,1.0000,1.0000,0.6412) diag(0.5015,0.9565,1.0000,0.9565,0.5015) diag(0.5007,0.9551,1.0000,1.0000,0.9551, 0.5007) diag(0.3969,0.9186,0.9952,1.0000,0.9952, 0.9186,0.3969) diag(0.3969,0.9185,0.9950,1.0000,1.0000, 0.9950,0.9185,0.3969) diag(0.3150,0.8797,0.9889,0.9995,1.0000, 0.9995,0.9889,0.8797,0.3150) diag(0.3150,0.8797,0.9889,0.9995,1.0000, 1.0000,0.9995,0.9889,0.8797,0.3150)

al

n 3 4 5

0.0796

A is a SHT

0.1017

A is a SHT

Jo

Besides, we consider other applications of Algorithm 1 in the higher-order Markov chains[54]. 4.2. Applications to high-order Markov chains Li and Ng [58] gave an approximated high-order Markov chain model as follows: Pxm−1 = x, ∥x∥1 = 1, (17) 25

Journal Pre-proof

p ro

i1 ∈[n]

of

where P = (pi1 i2 ···im ) is an mth order tensor representing an (m − 1)th order Markov chain, which is called an mth order n-dimensional transition probability tensor, i.e., ∑ pi1 i2 ···im > 0, pi1 i2 ···im = 1, and x is called a stochastic vector with xi > 0 and

∑

xi = 1. It was point-

i∈[n]

Pr e-

ed out in [59] that the nonlinear system (17) is equivalent to the following equation: { (I − βP)xm−1 = x[m−1] − βx, (18) ∥x∥1 = 1. If we choose the parameter β such that I − βP is an H+ -tensor, then we may use the tensor splitting method proposed by Liu, Li and Vong [59] to solve the nonlinear equation (18). We consider an example come from DNA sequence data in Tables 6 and 10 of [60] as follows: Example 4.4. Consider by  0.6000  (i) P(:, :, 1) = 0.2000 0.2000

two transition probability tensors [58, 60] defined

and

urn

al

   0.4083 0.4935 0.5217 0.3300 0.4152 0.2568 0.2426, P(:, :, 2) = 0.2232 0.2800 0.2658, 0.3349 0.2639 0.2551 0.3900 0.3190   0.5565 0.3648 0.4500 P(:, :, 3) = 0.2174 0.2742 0.2600, 0.2261 0.3610 0.2900

Jo

    0.5200 0.2986 0.4462 0.6514 0.4300 0.5776 (ii) P(:, :, 1) = 0.2700 0.3930 0.3192, P(:, :, 2) = 0.1970 0.3200 0.2462, 0.2100 0.3084 0.2346 0.1516 0.2500 0.1762   0.5638 0.3424 0.4900 P(:, :, 3) = 0.2408 0.3638 0.2900, 0.1954 0.2938 0.2200 respectively. Their orders m are both 3 and numbers of states n are both 3 by considering three categories ({A/G, C, T }). It was pointed out in [58] 26

Journal Pre-proof

p ro

of

that the model (17) has a unique positive solution for the two transition probability tensors. We will choose different parameters β from 0.200 to 0.400 with the step size 0.02, and determine whether I − βP is a strong H-tensor or not by Algorithm 1. The numerical results are showed in Table 3, where q, CPU(s) and X are the same as those used in Table 1, and TPT denotes the transition probability tensors. We see from Table 3 that the effectiveness of Algorithm 1 for identifying a given tensor whether is a strong H-tensor or not. Meanwhile, we find that I − βP is a strong H-tensor when β 6 0.340. In this case, we can use the tensor splitting method proposed by Liu, Li and Vong [59] to solve the nonlinear equation (18).

Jo

(ii)

CPU(s) 0.0073 0.0068 0.0104 0.0126 0.0104 0.0145 0.0238 0.0409 0.0188 0.0151 0.0120 0.0130 0.0100 0.0135 0.0116 0.0128 0.0181 0.0259 0.0279 0.0306 0.0177 0.0124

Output I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is not a SHT I − βP is not a SHT I − βP is not a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is a SHT I − βP is not a SHT I − βP is not a SHT I − βP is not a SHT

al

q 1 1 1 2 2 3 4 8 3 3 2 1 1 2 2 2 3 4 6 4 3 2

urn

(i)

β 0.200 0.220 0.240 0.260 0.280 0.300 0.320 0.340 0.360 0.380 0.400 0.200 0.220 0.240 0.260 0.280 0.300 0.320 0.340 0.360 0.380 0.400

Pr e-

Table 3: Numerical results of Algorithm 1 for Example 4.4

TPT

27

X diag(1.0000,0.8843,0.9114) diag(1.0000,0.8757,0.9045) diag(1.0000,0.8674,0.8978) diag(1.0000,0.7937,0.8401) diag(1.0000,0.7855,0.8331) diag(1.0000,0.7472,0.8016) diag(1.0000,0.7275,0.7845) diag(1.0000,0.7130,0.7709) diag(1.0000,0.7299,0.7855) diag(1.0000,0.7247,0.7805) diag(1.0000,0.7443,0.7964) diag(1.0000,0.8996,0.8643) diag(1.0000,0.8927,0.8549) diag(1.0000,0.8357,0.7737) diag(1.0000,0.8294,0.7648) diag(1.0000,0.8235,0.7565) diag(1.0000,0.7971,0.7165) diag(1.0000,0.7848,0.6970) diag(1.0000,0.7775,0.6851) diag(1.0000,0.7792,0.6890) diag(1.0000,0.7826,0.6957) diag(1.0000,0.7949,0.7162)

Journal Pre-proof

p ro

of

4.3. Applications to symmetric multi-player games Example 4.5. Consider evolutionary dynamics of symmetric four-players games with three strategies [12]. The first player’s payoff tensor A(1) ∈ R[4,3] given as follows     −1 3.5 1.5 3.5 2.5 1 1 , A(1) (1, 2, :, :) = 2.5 4 −2.5, A(1) (1, 1, :, :) = 3.5 2.5 1.5 1 −3.5 1 −2.5 1.5    1.5 1 −3.5 2.5 3.5 −2 −2.5 1.5 , A(1) (2, 1, :, :) = 3.5 3 0.5, A(1) (1, 3, :, :) =  1 −3.5 1.5 2 −2 0.5 0.5

Pr e-



    −2 0.5 0.5 3.5 3 0.5 5 −5.5, 5 , A(1) (2, 3, :, :) = 0.5 A(1) (2, 2, :, :) =  3 −4 0.5 −5.5 5 0.5 5 −5.5    3 2 2 2 3 −1 A(1) (3, 1, :, :) =  3 2 2 , A(1) (3, 2, :, :) = 2 3 −1, 2 −1 −1 −1 2 −2

al



urn

  −1 2 −2 A(1) (3, 3, :, :) =  2 −1 −1. −2 −1 3

The replicator equation is given by

Jo

dxi = xi ((A(1) xm−1 )i − A(1) xm ), i = 1, 2, 3, dt ∑ where x ∈ ∆3 := {(x1 , x2 , x3 )⊤ > 0 : xi = 1}. (1) m−1

Let xi ((A x

i∈{1,2,3}

(1) m

)i − A x ) = 0,

i = 1, 2, 3, we obatin

∑ 1 1 1 xi = 1} ( , , )⊤ ∈ int(∆3 ) := {(x1 , x2 , x3 )⊤ > 0 : 3 3 3 i∈{1,2,3}

28

(19)

Journal Pre-proof

p ro

of

is an interior equilibrium of the system (19), and Jacobian matrix of the system (19) at the equilibrium point ( 13 , 13 , 13 )⊤ is singular. Therefore, the stability of equilibrium point ( 13 , 13 , 13 )⊤ of the nonlinear system (19) can not be determined by the symbol of eigenvalues of Jacobian matrix. Let ∆ := {(y1 , y2 )⊤ : y1 , y2 > −1}. We could reduce the number of parameters by using the following transformations ψ : int(∆3 ) 7→ ∆ by ψ((x1 , x2 , x3 )⊤ ) = (y1 , y2 )⊤ , where y1 = xx13 − 1, y2 = is rewritten as

x2 x3

− 1, and x33 dt = dτ . The replicator equation (19) i = 1, 2,

Pr e-

dyi = (yi + 1)(By3 )i , dτ

(20)

where y = (y1 , y2 )⊤ , and B is a 4th order 2-dimensional tensor with entries ] ] [ [ 0.5 0.5 −3 0.5 , , B(1, 2, :, :) = B(1, 1, :, :) = 0.5 1 0.5 0.5 ] ] [ 0.5 1 0.5 0.5 . , B(2, 2, :, :) = B(2, 1, :, :) = 1 −7 0.5 1

al

[

urn

Since the transformation ψ is an one to one transformation from the int(∆3 ) of the system (19) to ∆ of the system (20), the stability of the system (19) at ( 31 , 13 , 31 )⊤ remain the same as the system (20) at (0, 0)⊤ . Choosing V (y1 , y2 ) := y1 − log(y1 + 1) + y2 − log(y2 + 1).

Jo

as the Lyapunov function. It can be verified that the function V (y1 , y2 ) is zero at the equilibrium (0, 0)⊤ and is positive for all other values of (y1 , y2 )⊤ ∈ ∆. By calculation, we have −V˙ (y1 , y2 ) = −

∑ ∂V y˙i = (−B)y4 . ∂yi

i∈{1,2}

Algorithm 3 is applied to polynomial −V˙ (y1 , y2 ), it terminates at the fourth step with the output that −V˙ (y1 , y2 ) is positive definite. By Lyapunov’s 29

Journal Pre-proof

of

direct method, the system (20) at (0, 0)⊤ is asymptotically stable. Thus the system (19) at ( 13 , 13 , 31 )⊤ is asymptotically stable. Consequently the system (19) is permanent. By choosing 6 initial values

p ro

(x1 (0), x2 (0), x3 (0)) = {(0.1, 0.5, 0.4), (0.2, 0.4, 0.6), (0.3, 0.2, 0.5), (0.4, 0.2, 0.4), (0.5, 0.3, 0.2), (0.6, 0.1, 0.3)},

5. Conclusions

Pr e-

we investigate the long-term behavior of the solution of the system (19) in the positive quadrant. The trajectories of the system (19) with different initial values are showed in Figure 2, which demonstrates that the equilibrium ( 13 , 13 , 13 )⊤ is asymptotically stable.

In this paper, we propose a parameter-free iterative algorithm of identifying strong H-tensors with weakly irreducibility, and prove that the algorithm always converges. Based on Algorithm 1, we present an algorithm for testing the positive definiteness of irreducible homogeneous polynomial forms. Numerical examples are given to show the effectiveness of the proposed method. Acknowledgements

Jo

urn

al

The authors would like to thank the anonymous referees for their valuable suggestions and comments. Qilong Liu’s work is supported by Doctoral Scientific Research Foundation of Guizhou Normal University in 2017 under grant GZNUD [2017]26. Jianxing Zhao’s work is supported partly by the National Natural Science Foundation of China under grant 11501141; the Foundation of Guizhou Science and Technology Department under grant [2015]2073; the Natural Science Programs of Education Department of Guizhou Province under grant [2016]066 and CAS “Light of West China” Program. Chaoqian Li’s work is supported partly by the Applied Basic Research Programs of Science and Technology Department of Yunnan Province under grant 2018FB001; Program for Excellent Young Talents in Yunnan University; Yunnan Provincial Ten Thousands Plan Young Top Talents. Yaotang Li’s work is supported partly by the National Natural Science Foundations of China under grant 11861077; the Joint Key Project of Yunnan Provincial Science and Technology Department and Yunnan University under grant 2018FY001(-014) and IRTSTYN. 30

Journal Pre-proof

(x1 (0),x2 (0),x3 (0))=(0.1,0.5,0.4) x

0.6

1

(b) x

2

x

(x1 (0),x2 (0),x3 (0))=(0.2,0.4,0.6)

of

(a)

x

0.6

3

0.4

0.2

0.2

p ro

0.4

0

1

x

2

x

3

0

0

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

0

(x1 (0),x2 (0),x3 (0))=(0.3,0.2,0.5)

(c)

(d)

x

0.6

1

x

2

x

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 (x1 (0),x2 (0),x3 (0))=(0.4,0.2,0.4)

0.6

3

1

x

2

x

3

0.4

0.2

Pr e-

0.4

x

0.2

0

0

0

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 (x1 (0),x2 (0),x3 (0))=(0.5,0.3,0.2)

(e)

(f)

x

0.6 0.4

0

1

x

2

x

0.6

3

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 (x1 (0),x2 (0),x3 (0))=(0.6,0.1,0.3) x

1

x

2

x

3

0.4

0.2

0.2

0

0 0

al

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

urn

0

Figure 2: The trajectories of the system (19)with different initial values

References

Jo

[1] L. Qi, Z. Luo, Tensor Analysis, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017. doi:10.1137/1.9781611974751. [2] E. Kofidis, P. A. Regalia, On the best rank-1 approximation of higherorder supersymmetric tensors, SIAM J. Matrix Anal. Appl. 23 (3) (2002) 863–884. [3] L. Qi, Eigenvalues of a real supersymmetric tensor, J. Symb. Comput. 40 (6) (2005) 1302–1324. 31

Journal Pre-proof

of

[4] L. De Lathauwer, B. De Moor, J. Vandewalle, On the best rank-1 and rank-(r1 , r2 , · · · , rn ) approximation of higher-order tensors, SIAM J. Matrix Anal. Appl. 21 (4) (2000) 1324–1342.

p ro

[5] T. Zhang, G. H. Golub, Rank-one approximation to high order tensors, SIAM J. Matrix Anal. Appl. 23 (2) (2001) 534–550. [6] L. Baratchart, M. Berthod, L. Pottier, Optimization of positive generalized polynomials under lp constraints, J. Convex Anal. 5 (2) (1995) 1303–1321.

Pr e-

[7] G. Zhou, L. Qi, S. Y. Wu, On the largest eigenvalue of a symmetric nonnegative tensor, Numer. Linear Algebra 20 (6) (2013) 913–928. [8] J. M. Smith, G. R. Price, The logic of animal conflict, Nature (1973) 15–18. [9] P. D. Taylor, L. B. Jonker, Evolutionarily stable strategies and game dynamics, Math. Biosci. 40 (1-2) (1978) 145–156. [10] K. Sigmund, J. Hofbauer, Evolutionary game dynamics, B. Am. Math. Soc. 40 (4) (2003) 479–519.

al

[11] M. Broom, C. Cannings, G. T. Vickers, Multi-player matrix games, B. Math. Biol. 59 (5) (1997) 931.

urn

[12] M. Broom, J. Rycht´ ar˘, Game-theoretical models in biology, CRC Press, 2013. [13] G. Palm, Evolutionary stable strategies and game dynamics for n-person games, J. Math. Biol. 19 (3) (1984) 329–334.

Jo

[14] Z. H. Huang, L. Qi, Formulating an N-person noncooperative game as a tensor complementarity problem, Comput. Optim. Appl. 66 (3) (2017) 557–576. [15] P. Schuster, K. Sigmund, Replicator dynamics, J. Theor. Biol. 100 (3) (1983) 533–538. [16] T. Myint-U, Ordinary Differential Equations, Elsevier Science Ltd, 1977.

32

Journal Pre-proof

of

[17] N. Bose, A. Modarressi, General procedure for multivariable polynomial positivity test with control applications, IEEE Trans. Automat. Control AC21 21 (5) (1976) 696–701.

p ro

[18] M. A. Hasan, A. A. Hasan, A procedure for the positive definiteness of forms of even order, IEEE Trans. Automat. Control 41 (4) (1996) 615–617. [19] E. Jury, M. Mansour, Positivity and nonnegativity conditions of a quartic equation and related problems, IEEE Trans. Automat. Control 26 (2) (1981) 444–451.

Pr e-

[20] F. Wang, L. Qi, Comments on “Explicit criterion for the positive definiteness of a general quartic form”, IEEE Trans. Automat. Control 50 (3) (2005) 416–418. [21] Y. Chen, Y. Dai, D. Han, W. Sun, Positive semidefinite generalized diffusion tensor imaging via quadratic semidefinite programming, SIAM J. Imaging Sci. 6 (3) (2013) 1531–1552. [22] S. Hu, Z. Huang, H. Ni, L. Qi, Positive definiteness of diffusion kurtosis imaging, Inverse Probl. Imaging 6 (2012) 57–75.

al

[23] L. Qi, G. Yu, E. X. Wu, Higher order positive semidefinite diffusion tensor imaging, SIAM J. Imaging Sci. 3 (3) (2010) 416–433.

urn

[24] L. Qi, G. Yu, Y. Xu, Nonnegative diffusion orientation distribution function, J. Math. Imaging Vis. 45 (2) (2013) 103–113. [25] S. Hu, L. Qi, Algebraic connectivity of an even uniform hypergraph, J. Comb. Optim. 24 (4) (2012) 564–579.

Jo

[26] S. Hu, L. Qi, The eigenvectors associated with the zero eigenvalues of the laplacian and signless laplacian tensors of a uniform hypergraph, Discrete Appl. Math. 169 (2014) 140–151. [27] S. Hu, L. Qi, J. Shao, Cored hypergraphs, power hypergraphs and their Laplacian H-eigenvalues, Linear Algebra Appl. 439 (10) (2013) 2980– 2998.

33

Journal Pre-proof

of

[28] G. Li, L. Qi, G. Yu, The Z-eigenvalues of a symmetric tensor and its application to spectral hypergraph theory, Numerical Linear Algebra Appl. 20 (6) (2013) 1001–1029.

p ro

[29] L. Qi, H+ -Eigenvalues of Laplacian and signless Laplacian tensors, Commun. Math. Sci. 12 (6) (2014) 1045–1064. [30] B. Reznick, Some concrete aspects of Hilbert’s 17th problem, Contemporary Mathematics 253 (2000) 251–272. [31] N. Z. Shor, Nondifferentiable optimization and polynomial problems, Vol. 24, Springer Science and Business Media, 2013.

Pr e-

[32] C. J. Hillar, L. Lim, Most tensor problems are NP-hard, J. ACM 60 (6) (2013) 45. [33] K. Zhang, Y. Wang, An H-tensor based iterative scheme for identifying the positive definiteness of multivariate homogeneous forms, J. Comput. Appl. Math. 305 (2016) 1–10. [34] L. Qi, Eigenvalues and invariants of tensors, J. Math. Anal. Appl. 325 (2) (2007) 1363–1377.

al

[35] C. Li, F. Wang, J. Zhao, Y. Zhu, Y. Li, Criterions for the positive definiteness of real supersymmetric tensors, J. Comput. Appl. Math. 255 (2014) 1–14.

urn

[36] J. Cui, G. Peng, Q. Lu, Z. Huang, New iterative criteria for strong H-tensors and an application, J. Inequal. Appl. (1) (2017) 49. doi: 10.1186/s13660-017-1323-1.

Jo

[37] M. R. Kannan, N. Shaked-Monderer, A. Berman, Some properties of strong H-tensors and general H-tensors, Linear Algebra Appl. 476 (2015) 42–55. [38] Y. Li, Q. Liu, L. Qi, Programmable criteria for strong H-tensors, Numer. Algorithms 74 (1) (2017) 199–221. [39] F. Wang, D. Sun, New criteria for H-tensors and an application, Open Math. 13 (1).

34

Journal Pre-proof

of

[40] F. Wang, D. Sun, New criteria for H-tensors and an application, J. Inequal. Appl. 2016 (1) (2016) 1–12. [41] X. Wang, Y. Wei, H-tensors and nonsingular H-tensors, Front. Math. China 11 (3) (2016) 557–575.

p ro

[42] Y. Wang, G. Zhou, L. Caccetta, Nonsingular H-tensor and its criteria, J. Ind. Manag. Optim 12. [43] Q. Liu, C. Li, Y. Li, On the iterative criterion for strong H-tensors, Comput. Appl. Math. (2016) 1–13.

Pr e-

[44] Y. Yang, Q. Yang, Further results for Perron-Frobenius theorem for nonnegative tensors, SIAM J. Matrix Anal. Appl. 31 (5) (2010) 2517– 2530. [45] L. Qi, Y. Song, An even order symmetric B tensor is positive definite, Linear Algebra Appl. 457 (457) (2014) 303–312. [46] L. Lim, Singular values and eigenvalues of tensors: a variational approach, in: Computational Advances in Multi-Sensor Adaptive Processing, 2005 1st IEEE International Workshop on, IEEE, 2005, pp. 129–132.

al

[47] L. Zhang, L. Qi, G. Zhou, M-tensors and some applications, SIAM J. Matrix Anal. Appl. 35 (2) (2014) 437–452.

urn

[48] W. Ding, L. Qi, Y. Wei, M-tensors and nonsingular M-tensors, Linear Algebra Appl. 439 (10) (2013) 3264–3278. [49] J. Shao, A general product of tensors with applications, Linear Algebra Appl. 439 (8) (2013) 2350–2366.

Jo

[50] S. Friedland, S. Gaubert, L. Han, Perron-Frobenius theorem for nonnegative multilinear forms and extensions, Linear Algebra Appl. 438 (2) (2013) 738–749. [51] M. Alanelli, A. Hadjidimos, A new iterative criterion for H-matrices, SIAM J. Matrix Anal. Appl. 29 (1) (2006) 160–176. [52] G. Zhou, L. Qi, S. Wu, Efficient algorithms for computing the largest eigenvalue of a nonnegative tensor, Front. Math. China 8 (1) (2013) 155–168. 35

Journal Pre-proof

of

[53] Y. Wei, W. Ding, Theory and computation of tensors: multi-dimensional arrays, Academic Press, 2016. [54] W. K. Ching, X. Huang, M. K. Ng, T. K. Siu, Markov Chains: Models, Algorithms and Applications, Springer, 2006.

p ro

[55] W. Ding, Y. Wei, Solving multi-linear systems with H-tensors, J. Sci. Comput. 68 (2) (2016) 689–715. [56] X. Li, M. K. Ng, Solving sparse non-negative tensor equations: algorithms and applications, Front. Math. China 10 (3) (2015) 649–680.

Pr e-

[57] X. Wang, M. Che, Y. Wei, Existence and uniqueness of positive solution for H+ -tensor equations, Appl. Math. Lett. 98 (2019) 191 – 198. doi: https://doi.org/10.1016/j.aml.2019.05.046. [58] W. Li, M. K. Ng, On the limiting probability distribution of a transition probability tensor, Linear Multilinear A. 62 (3) (2014) 362–385. [59] D. Liu, W. Li, S.-W. Vong, The tensor splitting with application to solve multi-linear systems, J. Comput. Appl. Math. 330 (2018) 75–94.

Jo

urn

al

[60] A. Raftery, S. Tavar´e, Estimation and modelling repeated patterns in high order markov chains with the mixture transition distribution model, J. R. Stat. Soc.: Series C (Applied Statistics) 43 (1) (1994) 179–199.

36