國立高雄大學應用數學學系
碩士論文
多維陣列上的
Perron-Frobenius 定理及計算
Perron-Frobenius theorem and computation on
multidimensional arrays
研究生：林冠宇 撰
指導教授：郭岳承
致謝
首先，我要感謝指導教授郭岳承老師，在我的碩士期間教導
我許多矩陣理論及計算的相關知識，也在我遇到困難時給予適當
的建議及鼓勵我參與許多校外活動。感謝吳宗芳老師大學及碩士
期間的照顧，能在大學時期上你的課以及擔任你兩年的大學入門
助教可以說是十分的幸運。感謝系辦雅鳳姐及佩君姊這些年的行
政上的各種協助。感謝我的家人一路上精神的支持讓我得以完成
碩士的學位。感謝學弟妹平常與我聊天及分享，陪我度過在高大
應數最後的時光。最後感謝系上的一切，讓我在這裡留下許多美
好的回憶。
林冠宇 謹於
國立高雄大學應用數學系
Perron-Frobenius theorem and
computation on multidimensional
arrays
by
Guan-Yu Lin
Advisor
Yueh-Cheng Kuo
Department of Applied Mathematics,
National University of Kaohsiung
Kaohsiung, Taiwan 811, R.O.C.
June 2018
Contents
1 Introduction 1
2 Preliminaries for matries 2
2.1 Irreducible matries . . . 2
2.2 Gershgorin’s disk theorem . . . 2
2.3 Implicit function theorem . . . 3
3 Perron-Frobenius theorem 5 3.1 Analysis of nonnegative irreducible matrices . . . 5
3.2 Proof of Theorem 3.1 . . . 7
4 Numerical experiments 9 4.1 Homotopy Continuation method . . . 9
4.2 Other numerical method . . . 12
4.3 Compare the numerical results . . . 13
5 Preliminaries for tensors 16 5.1 Tensor notations and operations . . . 16
5.2 Eigenvalue of tensors . . . 18
5.3 Symmetric and semi-symmetric tensors . . . 19
5.4 Irreducible tensors . . . 20
5.5 Rank-1 tensors . . . 20
6 Perron-Frobenius theorem on tensors 21 6.1 Analysis of nonnegative irreducible tensors . . . 21
6.2 Proof of Theorem 6.1 . . . 22
7 Numerical experiments on tensors 25 7.1 Jacobian matrix . . . 26
7.2 Other algorithms for tensors . . . 26
7.3 Compare the numerical results on tensors . . . 27
多維陣列上的
Perron-Frobenius 定理及計算
指導教授：郭岳承 教授 國立高雄大學應用數學系 學生：林冠宇 國立高雄大學應用數學系 摘要 Perron-Frobenius 定理在非負不可約矩陣上展示了一些重要的結果。此定理目前 已有許多的應用及推廣。在本論文中，我們著重在非負不可約矩陣上。在建構出一線 性同倫後，我們分析解曲線並證明任意一個非負不可約矩陣有一正特徵對。這個是 Frobenius 定理主要的結果。此證明技巧可被推廣用來證明張量上的 Perron-Frobenius 定理。除此之外，我們也呈現了延拓法在非負不可約矩陣及張量上的數值結 果。 關鍵字： Perron-Frobenius 定理、不可約矩陣及張量、延拓法 、隱函數定理Perron-Frobenius theorem and computation on
multidimensional arrays
Advisor: Dr. Yueh-Cheng Kuo Department of Applied Mathematics
National University of Kaohsiung
Student: Guan-Yu Lin
Department of Applied Mathematics National University of Kaohsiung
ABSTRACT
The Perron-Frobenius theorem shows some important results on nonnegative irreducible matrices. This theorem has various applications and extensions. In this paper, we focus on the nonnegative irreducible matrices. After constructing a linear homotopy, we analyze the solution curve of the linear homotopy and prove that any nonnegative irreducible matrix has the positive eigenpair. This is the main result of the Perron-Frobenius theorem. The skill of the proof can be extended to prove the Perron-Frobenius theorem on tensors. Furthermore, the numerical results computed by the homotopy continuation method on nonnegative irreducible matrices and tensors are presented.
Keywords: Perron-Frobenius theorem, Nonnegative irreducible matrices and tensors, Continuation method, Implicit function theorem
1
Introduction
Perron-Frobenius theorem is a fundamental and useful theorem for nonnegative irreducible matrices. It guarantees that the spectral radius is a simple eigenvalue corresponding to a positive eigenvector for any nonnegative irreducible matrix. The theorem has various applications, e.g., Markov chains, graph theories, and population models. The result of the theorem is elegent but the proof of the theorem is not familiar to the beginner of learning linear algebra. The origin proof has lots of matrix limits and it can be found in [10]. So, how to give a simpler proof is the first problem that what we concern and we will re-prove it with a new way.
The multidimensional arrays are often called the tensors. In fact, there are many properties of matrices that have been extended on tensors, including the eigenvalue problem. There are many researches about the spectral theories and computation of the eigenvalues on tensors, including the extension of the Perron-Frobenius theorem. W. Ding and Y. Wei [15] organized the knowledges of tensors, including the notations and applications. Chang, Pearson, and Zhang used the Brower’s fixed point theorem to extend the Perron-Frobenius theorem on tensors [5]. They even surveyed the spectral theories and applications of tensors in [4]. Kuo, Lin, and Liu [14] used the continuation method to compute the positive eigenpair for nonnegative weakly irreducible tensors. Kuo’s research almost proved the Perron-Frobenius theorem on tensors, except the spectral radius. So, we will continue the work of [14] and use the skill of the proof in [14] to prove the Perron-Frobenius theorem on tensors.
The main contributions are that re-proving the Perron-Frobenius theorem by a new way and extending the skill of the proof on tensors. This paper is organized as follow. In section 2 is about the preliminaries, including the irreducible matrices, the Gershgorin’s disk theorem and the implicit function theorem. In section 3, we analyze the nonnegative irreducible matrices first and re-prove the Perron-Frobenius theorem by obtaining the positive eigenpair. In section 4, the homotopy continuation method is presented to compute the positive eigenpair for nonnegative irreducible matrices. Given an initial condition at t = 0 and construct a linear
homotopy F1(x, λ, t) = 0. We compute the predict vector of possible solution and
use the Newton’s method to correct the solution. Furthermore, we compare the numerical results with other algorithms. Section 5 is about the preliminaries for tensors. In section 6, we analyze the nonnegative irreducible tensors first and prove the Perron-Frobenius theorem on tensors by obtaining the positive eigenpair of any nonnegative irreducible tensor. Section 7 is the numerical results on tensors. We
construct a linear homotopy F2(x, λ, t) = 0 for tensor and use the same homotopy
continuation method in section 4 to compute the positive eigenpair on tensors. Final section is about the conclusion.
2
Preliminaries for matries
In this section, some useful theorems will be introduced. For the consistency of symbols, we use the capital letters to denote matrices and lowercase letters to denote the scalars, respectively. Bold letters are denoted as the vectors. Let Rn×n>0 (R
n×n
≥0 ) is the set of all positive(nonnegative) matrices and ρ(.) is the spectral
radius of a matrix.
2.1
Irreducible matries
Definition 2.1. A ∈ Cn×n is reducible if there is a permutation matrix P ∈ Cn×n
such that
PTAP = B C
O D
, (2.1)
where B ∈ Cr×r, C ∈ Cr×(n−r), and D ∈ C(n−r)×(n−r). If A is not reducible, then
A is irreducible.
On matrices, the reducibility is an important property. We can find that if a matrix A is reducible then the eigenvalues of A are the collection of the eigenvalues of B and D and det(A) = det(B)det(D). Furthermore, the less computations will be required if a matrix of linear system or eigenvalue problem is reducible.
2.2
Gershgorin’s disk theorem
The Gershgorin’s disk theorem is a fundamental theorem for the upper-bound of eigenvalues of matrices. The statement of the theorem is as follows.
Definition 2.2. Let A ∈ Cn×n and ρi(A) =Pnj=1|Ai,j| for i = 1, 2, ...n. Then the
ith Gerschgorin’s disk Ci is defined
Ci = {z ∈ C : |z − Aii| ≤ ρi(A) − |Aii|}. (2.2)
Theorem 2.1. (Gershgorin’s Disk Theorem) Let A ∈ Cn×n. Then every
eigen-value of A is contained in a Gershgorin disk. We give an example to explain this theorem.
Example 1. Let A = 2 + 2i 2
1 −3i
, then the eigenvalues λ1, λ2 are λ1 =
2.1648 + 1.6463i, λ2 = −0.1648 − 2.6463i. We can see that there are two
-1 0 1 2 3 4 real axis -4 -3 -2 -1 0 1 2 3 4 imaginary axis 2.1648 + 1.6463i -0.1648 - 2.6463i
Figure 1: Eigenvalues of A in each Gershgorin Disk
2.3
Implicit function theorem
In this subsection, we will introduce the implicit function theorem. The proofs can be found in [3]. Following definitions are about the critical value and regular value.
Definition 2.3. Let F : Rm → Rn be continuously differentiable and let
C = {x ∈ Rm | rank(DxF (x)) < min(m, n)},
then
(a). Let x ∈ C and x is called the singular point of G.
(b). Let x ∈ Rm\C and x is called the regular point of G.
(c). F(C) is the set of critical values of F.
(d). Rn\F (C) is the set of regular values of F.
We give a simple example to explain.
Example 2. Let F (x, t) = (x − 2)(2x2 − t), where x, t ∈ R. Then
DF (x, t) = [6x2− 8x − t, 2 − x],
1. (2,8) is the critical point.
2. R2\{(2, 8)} is the set of regular point.
3. 0 is the critical value.
4. R\{0} is the set of regular value.
Next, we introduce the implicit function theorem
Theorem 2.2. (Implicit function theorem) Let F : Rn× R → Rn and it satisfies
(a). F (u0, t0) = 0, where u0 ∈ Rn, t ∈ R.
(b). DuF (u0, t0) is invertible.
(c). F and DuF (u, t) are smooth near u0.
Then there is a smooth function u(t),u : Bρ1(t0) → Bρ2(u0) such that
1. F (u(t), t) = 0, for all t ∈ Bρ1(t0) and u(t0) = u0,
2. The solution of F (u, t) = 0 in Bρ2(u0) is only u(t).
We see the example 2 and obtain that there is no smooth function near (2, 8)(See the figure 2). Theorem 2.2 says that if a function F satisfies the conditions in The-orem 2.2 and there will be a unique solution curve that can be parametrized by t.
-5 -4 -3 -2 -1 0 1 2 3 4 5 x -20 0 20 40 60 80 100 t (2,8)
3
Perron-Frobenius theorem
3.1
Analysis of nonnegative irreducible matrices
First, we introduce the Perron-Frobenius theorem.
Theorem 3.1. (Perron-Frobenius Theorem) A matrix A ∈ Rn×n≥0 is irreducible,
then the following statements are hold (a). ρ(A) > 0 is a simple eigenvalue.
(b). There is a positive eigenvector x ∈ Rn
>0 such that Ax = ρ(A)x .
(c). If λ is an eigenvalue with a nonnegative eigenvector, then λ = ρ(A).
Our goal is to prove the Theorem 3.1. The most noticeable feature of any nonnegative irreducible matrices is the positive eigenpair. In fact, it is sufficient to prove the Perron-Frobenius theorem by obtaining a positive eigenpair for any nonnegative irreducible matrix. First, there is no nonnegative eigenpair (λ, x) ∈
∂(Rn+1≥0 ) for any nonnegative irreducible matrix, where ∂(Rn+1≥0 ) is the boundary of
Rn+1≥0 .
Lemma 3.1. Let A ∈ Rn×n≥0 be irreducible. Suppose x ∈ Rn≥0 is an eigenvector of
A coeersponding to λ, then x ∈ Rn
>0 and λ > 0.
Proof. Suppose x ≥ 0, x is nonzero and there is a nonempty set S ⊂ {1, 2, ..., n}
such that xi = 0, where i ∈ S. Then there is a permutation matrix P ∈ Rn×n such
that PTx = x∗
0
∈ Rn
≥0, where x∗ ∈ Rk>0, 0 < k < n. And we obtain
λPTx = PTλx = PTAx = PTAP PTx. (λ, PTx) is an eigenpair of (PTAP ). So, PTAP = B C D E x∗ 0 = λ x∗ 0 , where B ∈ Rk×k≥0 , C ∈ R k×(n−k) ≥0 , D ∈ R (n−k)×(n−k) ≥0 , and E ∈ R (n−k)×(n−k) ≥0 . Hence,
Dx∗ = 0. Since x∗ > 0, we get D is a zero matrix and A is a reducible matrix,
a contradiction. Since A is nonnegative irreducible matrix and x ∈ Rn
>0, λx =
Ax > 0. So, λ > 0.
By Lemma 3.1, if there is an eigenvector x ∈ Rn≥0 of an irreducible matrix
A ∈ Rn×n≥0 , then x must be positive. And the following lemma can be obtained if
Lemma 3.2. Let A ∈ Rn×n≥0 be irreducible. Suppose x ∈ Rn>0 is an eigenvector of
A corresponding to λ ∈ R>0, then
(a). nullity(A − λI) = 1.
(b). The left eigenvector is also positive. (c). λ is a simple eigenvalue.
(d). λ = ρ(A).
(e). The matrix M = A − λI −x
2xT 0
∈ R(n+1)×(n+1) is invertible.
Proof. (a). Suppose (λ, y) ∈ Rn+1 is also an eigenpair of A. Then there exists
β ∈ R such that
v = x − βy ≥ 0.
If y ∈ Rn
>0, β = min1≤i≤n(xi/yi). Otherwise, β = max1≤i≤n(xi/yi) < 0. Vector
v is a nonnegative vector with at least an element is 0. If v 6= 0, then (λ, v) is an eigenpair of A. But it is a contradiction to Lemma 3.1. So, v = 0 and nullity(A − λI) = 1.
(b). Let a nonzero y ∈ Rn be a left eigenvector of A such that yTA = λyT. So,
λ | yT |=| λyT |=| yTA |≤| yT || A |=| yT | A,
and let zT =| yT | (A − λI) ≥ 0. Since (λ, x) is an eigenpair of A. If z 6= 0,
then zTx =| yT | (A − λI)x = 0. But it is a contradiction, because z ≥ 0, x > 0
, zTx 6= 0. So, z = 0, and (λ, | y |) is the left eigenpair of A. We know the left
eigenvector of A is the eigenvector of AT. By Lemma 3.1, the elements of | y | are
all positive .
(c). Suppose λ is not simple, then there is a x1 ∈ Rn such that
(A − λI)x1 = x ∈ Rn>0.
There is also a left positive eigenpair (λ, y) ∈ Rn+1>0 of A by Lemma 3.2(b). Then
yT(A − λI)x
1 = 0 = yTx, a contradiction. So, λ is simple.
(d). λ ≤ ρ(A) is clearly. Suppose u is an eigenvalue of A and |u| = ρ(A). Then
there is a Gershgorin’s disk Ci such that ρ(A) ∈ Ci with some i ∈ {1, 2, ...n}. And
|u| ≤ |u − Aii| + |Aii| ≤Pnj=1Aij. Then
ρ(A) = |u| ≤ max
1≤i≤n n
X
j=1
Since A is similar to Q−1AQ, where Q is the diagonal matrix with Qii = vi and
vi > 0 for i = 1, 2, ...n. It implies
ρ(A) = |u| ≤ max
1≤i≤n 1 vi n X j=1 Aijvj (3.2)
for any positive vector v ∈ Rn. Take (λ, x) in the inequality (3.2), and we get
ρ(A) ≤ λ. So, λ = ρ(A).
(e). Suppose M is singular, then there exist a nonzero vector zT = zT1, z2
T ∈ Rn+1. such that M z = A − λI −x 2xT 0 z1 z2 = 0.
By (b), there is a left eigenvector y ∈ Rn>0 such that
yT(A − λI)z1 = z2yTx = 0.
It implies
z2 = 0, z1 = αx,
where α ∈ R and α 6= 0. And then z1x = 2αxTx 6= 0, a contradiction.
Let A ∈ Rn×n≥0 be irreducible. By Lemma 3.2, if we can obtain a positive
eigenpair (λ, x) ∈ R>0× Rn>0, then λ = ρ(A) is a simple eigenvalue. And x is the
unique positive eigenvector. By Lemma 3.2(a)-(d), the Perron-Frobenius theorem can be proved by obtaining a positive eigenvector for any nonnegative irreducible matrix. Lemma 3.2(e) is an useful result for later analysis.
3.2
Proof of Theorem 3.1
Let A0 = x1xT2 be a rank-1 matrix, where x1, x2 ∈ Rn>0. And we let
λ0 = xT1x2, x0 =
x1
||x1||
. (3.3)
By Lemma 3.2, we know that (λ0, x0) is the unique positive eigenpair of A0. Now,
we give an irreducible matrix A ∈ Rn×n≥0 and let
A(t) = (1 − t)A0+ tA, t ∈ [0, 1]. (3.4)
Let F1 : Rn+1 × R → Rn+1 be continuous and differentiable. We consider the
following linear homotopy
F1(x, λ, t) =
A(t)x − λx
xTx − 1
From equation (3.5), (x0, λ0, 0) is the unique positive solution of F1(x, λ, t) = 0
and (λ0, x0) is the positive eigenpair of A0. Let (x∗, λ∗, 1) be the solution of
F1(x, λ, 1) = 0 as t = 1 and it is clearly to see that (λ∗, x∗) is the eigenpair of A.
By Lemma 3.2, if (λ∗, x∗) is the positive eigenpair of A, then the proof of Theorem
3.1 will be finished. We differentiate F1(x, λ, t) and get the Jacobian matrix of
F1(x, λ, t), DF1(x, λ, t) = [Dx,λF1(x, λ, t)|DtF1(x, λ, t)] ∈ R(n+1)×(n+2), where Dx,λF1(x, λ, t) = A(t) − λI −x 2xT 0 ∈ R(n+1)×(n+1), (3.6a) DtF1(x, λ, t) = (A − A0)x 0 ∈ Rn+1. (3.6b)
So, the following lemma can be obtained.
Lemma 3.3. If (x∗, λ∗, t∗) ∈ Rn+1≥0 × [0, 1) is an solution of F1(x, λ, t) = 0, then
the Jacobian matrix Dx,λF1(x∗, λ∗, t∗) is invertible.
Proof. For t∗ ∈ [0, 1), A(t∗) ∈ Rn×n>0 is irreducible. (λ∗, x∗) ∈ Rn+1>0 is an eigenpair
of A(t∗) by Lemma 3.1. By Lemma 3.2(e), Dx,λF1(x∗, λ∗, t∗) is invertible.
Since A0 is a positive rank-1 matrix and (λ0, x0) that defined in (3.3) is the
eigenpair of A. By Lemma 3.2(d) and the implicit function theorem, there is a
solution curve s1(t) = (x(t), λ(t), t) ∈ Rn+1≥0 × [0, 1) of F1(x, λ, t) = 0 that is
parametrized by t, where s1(0) = (x0, λ0, 0). For t ∈ [0, 1), we can get s(t) ∈
Rn+1>0 × [0, 1) and (λ(t), x(t)) /∈ ∂(R n+1
≥0 ) by Lemma 3.1. By Lemma 3.2(d), λ(t) =
ρ(A(t)) and we get ||(λ(t), x(t))|| is bounded. We can obtain the nonnegative
solution set of F1(x, λ, t) = 0,
C1 = {(x(t), λ(t), t)|F1(x(t), λ(t), t) = 0} ⊆ Rn+1≥0 × [0, 1), (3.7)
where (x0, λ0, 0) ∈ C1. So, the following theorem can be obtained.
Theorem 3.2. Let A ∈ Rn×n≥0 be irreducible. Then the nonnegative solution of
F1(x, λ, 1) = 0 is isolated.
Proof. Let (x∗, λ∗, 1) be the accumulation point of C1. (λ∗, x∗) is the eigenpair of
A. Since A is irreducible. So, (λ∗, x) ∈ Rn+1>0 . λ∗ is the root of the characteristic
polynomial of A and λ∗ is isolated. By Lemma 3.2(c), we get the nonnegative
solution of F1(x, λ, 1) = 0 is isolated.
Theorem 3.2 shows that any irreducible matrix A ∈ Rn×n≥0 has the positive
eigenpair (λ∗, x∗) of A by tracing the solution curve s1(t) of F1(x, λ, t) = 0 in (3.5)
with initial (λ0, x0) in (3.3). By Lemma 3.2(a)-(d), the Theorem 3.1 has been
4
Numerical experiments
In this section, the homotopy continuation method will be introduced for
comput-ing the positive solution of F1(x, λ, 1) = 0 which is defined in (3.5) first. The other
numerical methods will be introduced and the experiments results will be shown.
4.1
Homotopy Continuation method
In section 3, the Perron-Frobenius theorem has been proved and we know that the positive eigenpair of nonnegative irreducible matrices can be obtained by tracing the solution curve of the linear homotopy in (3.5). So, we use the homotopy continuation method to experiment. The homotopy continuation method includes
the predictor process and the correction process for each iteration. Let x0, λ0 in
(3.3) be the initial condition of F1(x, λ, t) = 0 in (3.5) as t = 0. Let (xTi , λi, ti)T
be the i-th step which satisfies F1(x, λ, t) = 0 for i = 1, 2, ... . Following is about
the prediction process and the correct process.
Predict process
We use the tangent vector as the predict vector. We differentiate the equation
F1(x, λ, t) = 0 at (xi, λi, ti) and get Dx,λF1(xi, λi, ti) ˙ xi ˙ λi + DtF1(xi, λi, ti) ˙ti = 0.
Let ˙ti = ∆ti, where ∆ti is the step length and ti+1= ti+ ∆ti ≤ 1. And we get
˙ xi ˙λi = −Dx,λF1(xi, λi, ti)−1DtF1(xi, λi, ti). So, h x(1)Ti+1, λ(1)i+1 iT
is the predict vector, where " x(1)i+1 λ(1)i+1 # = xi λi + ∆ti ˙ xi ˙ λi .
Correct process
For computing the approximate solution of F1(xi+1, λi+1, ti+1) = 0, we use the
Newton method with the initial x(1)i+1, λ(1)i+1 and ti+1 is fixed. Each iteration of
Newton method satisfies " x(l+1)i+1 λ(l+1)i+1 # = " x(l)i+1 λ(l)i+1 # + δl, where δl satisfies Dx,λF1(x (l) i+1, λ (l) i+1, ti+1)δl = −F1(x (l) i+1, λ (l)
i+1, ti+1) for l = 2, 3, . . ..
We use the Newton method until that the (x(l+1)i+1 , λ(l+1)i+1 , ti+1) converges to (xi+1,
λi+1, ti+1) and satisfies F1(xi+1, λi+1, ti+1) = 0 approximately. The convergent
behavior of ||F1(x
(l+1) i+1 , λ
(l+1)
i+1 , ti+1)|| will be quadratic. After we get (xi+1, λi+1,
ti+1), we continue the predict-correct process until t = 1.
Adjustment of the step length
An iteration is composed by a predict and a correct process. ”∆t” is the step length for the predict process and ”times” is the numbers of the successful iterations. After the predict process, we execute the correct process. If the correct process
is not convergent to the possible solution of F1(x, λ, t) = 0, then ∆t will be cut
half. If ∆t < 10−16, the homotopy continuation method will be stop. If we get
a successful iteration, times will be added 1. While the iteration successes twice, the step length will be enlarged to accelerate the iteration. Then we set times = 0 and continue the iteration. Let the initial ”acc”= 1 and following is the algorithm of the homotopy continuation method for computing the nonnegative eigenpair of a nonnegative irreducible matrix.
Algorithm 1 Homotopy continuation method
Step 1. Set initial condition A0, x0, λ0, t0 = 0, times = 0
Set count = the number of possible solution
Step 2. Compute the tangent vector ( ˙xi, ˙λi)
Compute the predict vector (x(1)i+1, λ(1)i+1) = (xi, λi) + ∆t( ˙xi, ˙λi)
ti+1= ti+ ∆t Step 3. Compute δ = −Dx,λF1(x (1) i+1, λ (1) i+1, ti+1)−1F1(x (1) i+1, λ (1) i+1, ti+1)
(x(2)i+1, λ(2)i+1) = (x(1)i+1, λ(1)i+1) + δ
If (||F1(x (1) i+1, λ (1) i+1, ti+1)||/||F1(x (2) i+1, λ (2)
i+1, ti+1)|| <10 and
||F1(x
(2) i+1, λ
(2)
i+1, ti+1)|| > 10−8acc and ||δ|| > 10−7)
∆t = ∆t/2, back to (xi, λi, ti), times=0
if (∆t < 10−16)
stop the algorithm end go to step 2 else go to step 4 end Step 4. While(||F1(x (l) i+1, λ (l)
i+1, ti+1)|| > 10−8acc)
Compute δ = −Dx,λF1(x (l) i+1, λ (l) i+1) −1F 1(x (l) i+1, λ (l) i+1)
(x(l+1)i+1 , λ(l+1)i+1 , ti+1) = (x(l)i+1, λ (l) i+1, ti+1) + δ If(||F1(x (l) i+1, λ (l) i+1, ti+1)||/||F1(x (l+1) i+1 , λ (l+1) i+1 , ti+1)|| < 5) If ||δ|| < 10−7
break the while loop and go to step 5 else
∆t = ∆t/2 ,back to (xi, λi, ti),times=0
if (∆t < 10−16)
stop the algorithm end go to step 2 end end end Step 5. If (times≥ 2) times=0, ∆t = 2∆t else times = times+1 end
(xi+1, λi+1, ti+1) = (x
(l) i+1, λ
(l) i+1, ti+1)
acc = p||Dx,λF1(xi+1, λi+1, ti+1)||||(xi+1, λi+1)||
If (acc < 1) acc=1 end If t = 1
stop the algorithm
4.2
Other numerical method
In this subsection, other algorithms [8] that compute the positive eigenpair of nonnegative irreducible matrices will be introduced.
Power method
Power method is a simple iteration method to compute the eigenpair of matrices.
It is developed by the matrix limits. Let A ∈ Cn×n and (λ
i, xi) ∈ C × Cn is the
eigenpair of A, where i = 1, 2, . . . , n and |λ1| > |λ2| ≥ |λ3| ≥ . . . ≥ |λn|. We can
use the following algorithm to compute (λ1, x1). And the convergent rate is base
on |λ2/λ1|.
Algorithm 2 Power method
1. Set q(0) ∈ Cn×n 2. For i=1,2,. . . z(k)= Aq(k−1) q(k)= z(k)/||z(k)|| λ(k) = [q(k)]HAq(k) end
Noda’s algorithm
The main concept of Noda’s algorithm [8] is the inverse iteration. For a
nonnega-tive irreducible matrix A ∈ Rn×n, we let λ(0) be an approximation eigenvalue of A
with λ(0) > ρ(A) and x(0) ∈ Rn
>0be a unit random vector. Let (λ(i), x(i)) be the
ap-proximation eigenpair of A. For each iteration, Noda’s algorithm keeps λ(i)> ρ(A)
and λ(i) is decreasing. Otherwise, (λ(i)I − A)−1x keeps positive. Finally, (λ(i), x(i))
Algorithm 3 Noda’s algorithm
1. Set initial eigenpair (λ(0), x0) ∈ R × Rn, where λ(0) > ρ(A). And let ¯λ(0) = λ(0).
2. For k = 0,1,2,. . . x∗(n+1) = (¯λ(n)I − A)−1x(n) τ(n+1) =< x∗(n+1), v > / < x(n), v > x(n+1) = x∗(n+1)/τ(n+1) ¯ τ = maxk((x∗(n+1))k/x (n) k ), τ = mink((x∗(n+1))k/x (n) k ) ¯ λ(n+1) = ¯λ(n)− 1/¯τ(n+1), λ(n+1) = ¯λ(n)− 1/τ(n+1) λ(n+1) = ¯λ(n)− 1/τ(n+1)
if ¯λ(n+1)− λ(n+1) < stop the algorithm end
Prakash’s algorithm
Prakash Chanchana presented the algorithm in 2007 [8]. The main concept of the algorithm are Collatz’s formula and the inverse iteration. The most interesting is that Prakash’s algorithm used the LU-factorization to solve the linear system twice.
Algorithm 4 Prakash’s algorithm
1. Let x(0)∈ Rn
>0, λ(0) = maxi(
Pn
j=1aij) and B(i) = (λ(0)I − A)−1
2. For i = 1,2,. . .
(λ(i)I − A) = L(i)U(i)
L(i)U(i)xˆ(i) = x(i)
L(i)U(i)x(n+1) = ˆx(i) ¯
λ(i) = λ(i)− minj(
ˆ x(i)j B(i)ˆx(i)
j
), λ(i) = λ(i)− maxj(
ˆ x(i)j B(i)ˆx(i) j ) λ(n+1) = ¯λ(n)
If ¯λ(i)− λ(i) < stop the algorithm
end
4.3
Compare the numerical results
In this subsection, we will present the numerical results of homotopy continuation method, including the behavior of the residual and the numbers of the positive
elements of x of F1(x, λ, t) = 0 for t ∈ [0, 1]. Each iteration includes a
predic-tion process and a correct process. The step length is ∆t = 0.1. Each iterapredic-tion
is required to stop by the stop criteria with tolerance = 10−10. Otherwise, we
also compare the numerical results with other algorithms in section 4.2. For con-venience, ”Continuation” is the homotopy continuation method, ”Power” is the
power method, ”Noda” is the Noda’s algorithm, and ”PC” is the Prakash’s al-gorithm. ”Count” is the number of solving linear system. Following numerical results are experimented by Matlab 2015b on an Intel(R) Core(TM) i5-8250 CPU @1.80 GHz with 4GB RAM and NVIDIA GeForce MX150.
We will construct the table to compare the different algorithms. The first row is the number of iterations. If the number of iterations is bigger than 2000, than we will stop the algorithm. The second row is the number of computing linear system.
The third row is the residual of F1(x, λ, 1). The fourth row is the subtraction with
the solution of each algorithm and the solution that be computed by the command ”eig” in Matlab.
Example 3. Let A, A0 ∈ Rn×n≥0 be random matrix, where A0 = x1xT2 with x1, x2 ∈
Rn>0. A is the case dwt 992 in [2] and A is a sparse symmetric matrix with
n = 992 and 16744 nonnzero elements. In Figure 3(a), it shows that the behavior
of the residual of F1(x, λ, t) = 0 for t ∈ [0, 1] and the homotopy continuation
method can compute the positive eigenvector correctly. In Figure 3(b), it shows that the numbers of the positive elements of the eigenvector for each iteration.
Figure 3: Compute for case in dwt 992
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 3 3.5 Residual of F1 ×10-14 (a) Residual of F1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 991 993
numbers of positive elements
(b) Numbers of positive elements of eigenvec-tor.
We can see that the homotopy continuation method can get the positive
eigen-pair successfully with the tolerance 10−10 but it need more computations in solving
linear systems. The power method presents the poor behavior of convergence
be-cause the rate of |λ2/λ1| = 0.9904.
Table 1: Compare for case in dwt 992
Continuation Power Noda PC
Iteration 5 816 8 3
Count 22 0 8 3
Residual 3.2041e−15 9.8707e−11 4.3471e−11 6.6069e−13
Difference 1.7091e−14 2.1754e−10 4.3478e−11 1.1316e−12
a p-cycle block matrix if it satisfies
A = 0 A1 0 . . . 0 0 0 A2 . . . 0 .. . ... ... . .. . . . 0 0 0 . . . Ap−1 Ap 0 0 . . . 0 , (4.1)
where Ai ∈ Cn×n for i = 1 . . . p. Let λ be the eigenvalue of A and
βk = ei2kπ/p, k = 0, 1, . . . p. (4.2)
In [6], it says that βkλ is also the eigenvalue of A. So the eigenvalue of A is
peri-odic with period 2π/p. It is clearly to see that the power method can not compute
the eigenpair of A. Now we let Ai ∈ Rn×n>0 for i = 1, 2 . . . , p and A be constructed in
(4.1). A is irreducible, there is a positive eigenpair of A by the Perron-Frobenius theorem. We let n = 50, p = 20 and get the numerical result like in example 3. Excluding the power method, each iteration can compute the positive eigenpair suc-cessfully and steadily.
Table 2: Compare for block-cyclic matrix
Continuation Power Noda PC
Iteration 5 > 2000 7 4
Count 17 0 7 4
Residual 2.0944e−15 0.3298 3.8447e−15 3.2600e−15
Difference 2.1936e−14 0.0647 2.2235e−14 2.2136e−14
Remark 4.1. The eigenvalue of A in (4.1) is the eigenvalue of Qp
i=1Ai with a
Figure 4: Compute for block-cyclic matrix 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 0 0.5 1 1.5 2 2.5 3 3.5 4 Residual of F1 ×10-13 (a) Residual of F1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t 999 1001
numbers of positive elements
(b) Numbers of positive elements of eigen-vector.
So, the homotopy continuation method can compute the nonnegative eigen-pair of nonnegative irreducible matrices and the solution of each iteration keeps positive.
5
Preliminaries for tensors
Here we use calligraphic letters, like A, to denote tensors. A mth-order tensor
A ∈ Cn1×n2×...×nm is a multidimensional array. In fact, a vector and a matrix are
the first-order tensor and second-order tensor, respectively. If m ≥ 3, A is called
the higher-order tensor. If n1 = n2 = ... = nm = n, we called A is a mth-order,
n-dimensional tensor and denote the set of all mth-order, n-dimensional tensor
over real(complex) field R(C) by R[m,n](C[m,n]). We let A ∈ Cn1×n2×...×nm be a
tensor and the element (i1, i2, ..., im) of A is denoted by ai1i2...im.
5.1
Tensor notations and operations
Both matrices and tensors can operate with vectors. However, the computation of tensors is more complex than matrix. So, we usually transform a tensor into a matrix and the following definition is very useful for us [13].
Definition 5.1. Let A ∈ Cn1×n2×...×nm be a tensor. The mode-k fiber is that fixing
every index of ai1i2...im except the k-th element. The mode-k unfolding is denoted
by A(k) ∈ Cnk×(
Qm
i6=kni) and arranges the k-fibers to be the columns of the matrix
Example 5. We use the example in [13]. Let a tensor A ∈ R3×4×2 exist, where A::1 = 1 4 7 10 2 5 8 11 3 6 9 12 , A::2 = 13 16 19 22 14 17 20 23 15 18 21 24 .
Following are the mode-k unfoldings for k = 1, 2, 3.
A(1)= 1 4 7 10 13 16 19 22 2 5 8 11 14 17 20 23 3 6 9 12 15 18 21 24 A(2) = 1 2 3 13 14 15 4 5 6 16 17 18 7 8 9 19 20 21 10 11 12 22 23 24 A(3) = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
Next, we introduce the Kronecker product. Kronecker is an useful tool to compute on tensors [13].
Definition 5.2. Let A ∈ Cm×n, B ∈ Cs×t, then the Kronecker product of A and
B is A ⊗ B = a11B a12B . . . a1nB a21B a22B . . . a2nB .. . ... . .. ... am1B am2B . . . amnB ∈ Cms×nt (5.1)
Next, it is about the k-mode product with vector of tensors.
Definition 5.3. For a tensor A ∈ Cn1×n2×...×nm, the k-mode product of A with a
vector x ∈ Cnk is denoted by A × kx ∈ Cn1×n2×...×nk−1×nk+1...×nm, and we have (A ×kx)i1...ik−1ik+1...im = nk X ik=1 ai1i2...imxik. (5.2)
Example 6. We use the tensor A in example 5. Let x = [1, 2, 3, 4]T. Then
A ×2x = 70 190 80 200 90 210 = A(1)(I2⊗ x),
5.2
Eigenvalue of tensors
Let A ∈ C[m,n] and x ∈ Cn. We denote Axm−1 = A ×
2x ×3x . . . ×mx. Following
definition is about the eigenvector and eigenvalue on tensors [12][15].
Definition 5.4. Let A ∈ C[m,n]. If the pair (λ, x) ∈ C × Cn\ {0} satisfies
Axm−1 = λx[m−1], (5.3)
then λ is called the eigenvalue of A corresponding to eigenvector x. Where x =
[x1, x2, ...xn]T and x[m−1] = [x1m−1, xm−12 , ..., xm−1n ]T. If (λ, x) ∈ R × Rn, then we
call (λ, x) is H-eigenpair.
Next definition is about the spectral radius of tensors.
Definition 5.5. Let A ∈ C[m,n]. We let
%(A) = max{ |λ| |λ is an eigenvalue of A.} (5.4)
and %(A) is called the spectral radius of A.
Solving an eigenvalue problem on tensors is a homogeneous problem. In [12], it defined the characteristic polynomial of a tensor A as det(A − λI) by using the resultant, where each element of I is
ii1i2...im =
1, if i1 = i2 = ...im
0, others. In [12], it also showed the following theorem.
Theorem 5.1. Let A ∈ C[m,n] be a tensor. And the det(A − λI) can be written
as
det(A − λI) = Y
λi∈σ(A)
(λ − λi)mi, (5.5)
where σ(A) = {λ|λ is the eigenvalue of A.} and mi is the algebraic multiplicities
of λi.
In addition to algebraic multiplicity, the geometric multiplicity on tensors is also defined. Next definition is about the geometric multiplicity of tensors [5].
Definition 5.6. Let A ∈ C[m,n] and λ be the eigenvalue of A. Then the
geomet-ric multiplicity of λ is the maximum number of linear independent eigenvector corresponding to λ.
Example 7. Let A ∈ R[3,2]≥0 , where a111 = a222 = 1, a122 = a211 = for 0 < < 1
and other aijk = 0. The eigenvalue of A will satisfy
x21+ x22 = λx21 x21+ x22 = λx22
We obtain det(A − λI) = (λ − 1 + )2(λ − 1 − )2. And we get the following results:
(a). λ = 1 + with eigenvector x1 = [1, 1]T, x2 = [1, −1]T. The geometric
multi-plicity of λ is 2. And x1, x2 are the H-eigenvectors corresponding to λ.
(b). λ = 1 − with eigenvector x1 = [1, i]T, x2 = [1, −i]T. The geometric
multi-plicity of λ is 2. But x1, x2 are complex eigenvectors.
5.3
Symmetric and semi-symmetric tensors
Symmetric tensors are very useful for computing the eigenvalue problems on ten-sors [15].
Definition 5.7. Let A ∈ C[m,n] be a tensor and A is symmetric if it satisfies
ai1i2...im = aj1j2...jm, (5.6)
where (j1, j2, . . . , jm) ∈ Θ(i1, i2, . . . , im) and Θ(i1, i2, . . . , im) is the set of all the
permutation of (i1, i2, . . . , im).
Following definition is the semi-symmetric tensor. The difference between sym-metric tensor and semi-symsym-metric tensor is about the first index.
Definition 5.8. Let A ∈ C[m,n] be a tensor and A is semi-symmetric if it satisfies
ai1i2...im = ai1j2...jm, (5.7)
where (j2, j3, . . . , jm) ∈ Θ(i2, i3, . . . , im) and Θ(i2, . . . , im) is the set of all the
per-mutation of (i2, . . . , im)
Remark 5.1. Let A ∈ C[m,n] be tensor. If A is symmetric then A is
semi-symmetric.
In [9] [1], let A ∈ C[m,n] and there is a semi-symmetric tensor A
s∈ C[m,n] such that Axm−1 = A sxm−1. It is efficient to compute Dx(Axm−1) = m X k=2 A ×2x... ×k−1x ×k+1... ×mx ∈ Rn×n. (5.8)
We can get Dx(Axm−1) = Dx(Asxm−1) and Dx(Asxm−1) = m X k=2 As×2x... ×k−1x ×k+1... ×mx (5.9) = (m − 1)As×3x ×4 x . . . ×mx ∈ Rn×n. (5.10)
The equation (5.8) is equivalent to equation (5.10).
5.4
Irreducible tensors
Following definition is about the irreducible tensors[15] [5].
Definition 5.9. A tensor A ∈ R[m,n] is called reducible if there is a nonempty
proper subindex I ⊂ {1, 2, ...n} such that
a1i2...im = 0, ∀i1 ∈ I, ∀i2, i3...im ∈ I./
If A is not reducible, then A is irreducible.
The tensor A in Example 7 is an irreducible tensor. If m = 2, the definition of irreducible tensors is equivalent to the definition of irreducible matrices.
5.5
Rank-1 tensors
Next definition is about the rank-1 tensors. Rank-1 tensors play an important role in tensor decomposition and rank computing [13] .
Definition 5.10. A tensor A ∈ Rn1×n2×...×nm is rank-1 if it can be written as the
outer product of m nonzeros vectors xk ∈ Rnk for k = 1, 2, ...m, and it is denoted
by A = x1◦ x2... ◦ xm. Each element of A is
ai1,i2,...im = x1,i1x2,i2...xm,im, (5.11)
where xk,ik is the ik-th element of vector xk.
In [14], it recorded the theorem about the H-eigenvalues of rank-1 tensors.
Lemma 5.1. Let A0 = x1 ◦ x2◦ ... ◦ xm ∈ R [m,n] ≥0 , where x1, x2, ...xm ∈ Rn≥0. Let λ0 = Qm k=2x T kx [1/(m−1)] 1 and x0 = x[1/(m−1)]1
kx[1/(m−1)]1 k. Then (λ0, x0) and (0, w) with w ∈
Sm
k=2span{xk}
⊥ are H-eigenpairs. In addition, if x
1 > 0, then the eigenvalue
λ0 > 0 with positive eigenvector x0. Furthermore, if x1, x2, ...xm ∈ Rn>0, then
cx ∈ Rn
>0 with c > 0 is the unique nonnegative eigenvector of A0.
So, we can find an unique positive eigenpair of a positive rank-1 tensor. To construct the linear homotopy like section 3, it is useful for us to set the initial condition.
6
Perron-Frobenius theorem on tensors
In this section, we will prove the Perron-Frobenius theorem on tensors. Let A ∈ R[m,n]≥0 be an irreducible tensor. And we focus on finding the positive H-eigenpair
of A. Here we let R[m,n]>0 (R
[m,n]
≥0 ) be the set of all real positive(nonnegative)
mth-order, n-dimension tensors.
6.1
Analysis of nonnegative irreducible tensors
Following theorem is the Perron-Frobenius theorem on tensors [5].
Theorem 6.1. (Perron-Frobenius Theorem on tensors) Let A ∈ R[m,n] be
irre-ducible, then
(a). There is a H-eigenvalue λ0 = %(A) > 0.
(b). There is a H-eigenvector x ∈ Rn
>0 corresponding to λ0.
(c). If λ is an eigenvalue with nonnegative eigenvector y, then λ = λ0. Moreover,
y = kx, where k is a constant.
Remark 6.1. Some different of Perron-Frobenius theorem on tensors is that λ0 in
Theorem 6.1 may not be simple. In Example 7, A is irreducible and we can find the λ = 1 + is the %(A) with the algebraic multiplicity 2 and two H-eigenvectors. No matter the Perron-Frobenius theorem on matrices or tensors, the positive eigen-vector is unique up to a constant.
In fact, we can prove Theorem 6.1 by obtaining the positive eigenpair of any
nonnegative irreducible tensors. First, there is no H-eigenpair in ∂(Rn+1≥0 ), where
∂(Rn+1≥0 ) is the boundary of Rn+1≥0 .
Lemma 6.1. Let A ∈ R[m,n]≥0 be irreducible. If x ∈ Rn≥0 is a H-eigenvector of A,
then x ∈ Rn>0. Furthermore, there is a λ > 0 and (λ, x) is the H-eigenpair of A.
Proof. Let x ∈ Rn
≥0 be a H-eigenvector with some subset I ⊂ {1, 2, ...n} such that
xi = 0 ∀i ∈ I and xj 6= 0 ∀j /∈ I. We have
n
X
i2,i3,...im=1
ai1i2...imxi2xi3...xim = 0 ∀i1 ∈ I.
Since xi 6= 0 ∀ i /∈ I and it implies ai1,i2,...im = 0 with ∀ i1 ∈ I. It implies
ai1,i2,...im = 0 with ∀i1 ∈ I, ∀ i2. . . im ∈ I. We obtain A is reducible, a contradic-/
tion. So, x ∈ Rn
>0and we can get that there is a H-eigenvalue λ > 0 corresponding
Next, we discuss the H-eigenvalue of A. Here we use a theorem in [5] and the proof can also be found in [5].
Theorem 6.2. Let A ∈ R[m,n]≥0 be irreducible. If (λ, x), (µ, y) ∈ Rn+1≥0 satisfies
Axm−1 = λx[m−1] and Aym−1 ≥ µy[m−1](Aym−1 ≤ µy[m−1]), then λ ≥ µ(λ ≤ µ).
So, the following theorem can be obtained.
Theorem 6.3. Let A ∈ R[m,n]≥0 be irreducible. If x ∈ Rn>0 is a H-eigenvector of A
corresponding to λ, then λ = %(A) > 0 is a H-eigenvalue of A.
Proof. Let (u, y) ∈ Rn+1 be an H-eigenpair of A, where |u| = %(A) and y is a
nonzero vector. Then
%(A)|y|[m−1]= |uy[m−1]| = |Aym−1| ≤ A|y|m−1.
So, by Theorem 6.2, %(A) ≤ λ. Since λ ≤ %(A) is clearly to see, and we have done λ = %(A).
Now, we know that there is no H-eigenpair on ∂(Rn+1≥0 ) for nonnegative
ir-reducible tensors. If there is a H-eigenvector x ∈ Rn
>0 of an irreducible tensor
A ∈ R[m,n]≥0 , then λ = %(A) is the H-eigenvalue corresponding to x. The
unique-ness of positive H-eigenvector will be discussed in later subsection by the implicit function theorem.
6.2
Proof of Theorem 6.1
We use the skill in [14]. Let A0 = x1◦ x2◦ ... ◦ xm ∈ R
[m,n] >0 be a rank-1 tensor, where x1, x2, ..., xm ∈ Rn>0. Let x0 = x[1/(m−1)] ||x[1/(m−1)]||, λ0 = m Y k=2 (xTx1/(m−1)1 ). (6.1)
By Lemma 5.1, (λ0, x0) is the unique positive H-eigenpair of A0. Now, we let
A(t) = (1 − t)A0+ tA, ∀t ∈ [0, 1]. (6.2)
We consider the following linear homotopy
F2(x, λ, t) =
A(t)xm−1 − λx[m−1]
xTx − 1
= 0, ∀t ∈ [0, 1]. (6.3)
Differentiating the F2(x, λ, t), and the Jacobian matrix of F2(x, λ, t) can be
ob-tained. The Jacobian matrix is denoted as DF2(x, λ, t) = [Dx,λF2(x,λ,t)|DtF2(x, λ, t)]
Dx,λF2(x, λ, t) = At− (m − 1)λ[|x[m−2]|] −x[m−1] 2xT 0 ∈ R(n+1)×(n+1) (6.4a) DtF2(x, λ, t) = (A − A0)xm−1 0 ∈ Rn+1. (6.4b)
Here [|x|] is the square diagonal matrix with the elements of x on the main diagonal and At ≡ Dx(A(t)xm−1) = m X k=2 A(t) ×2x... ×k−1x ×k+1... ×mx ∈ Rn×n. (6.5)
Suppose (x∗, λ∗, t∗) be the solution of F2(x, λ, t) = 0. The fact is
At∗x∗ = (m − 1)A(t∗)x m−1 ∗ = (m − 1)λ∗x[m−1]∗ , and (At∗− (m − 1)λ∗[|x [m−2] ∗ |])x∗ = 0. (6.6)
So, the following lemma can be obtained.
Lemma 6.2. Let (x∗, λ∗, t∗) ∈ Rn+1≥0 × [0, 1) be a solution F2(x, λ, t) = 0. Then
the Jacobian matrix Dx,λF2(x∗, λ∗, t∗) is invertible.
Proof. For t∗ ∈ [0, 1), A(t∗) > 0 is irreducible. So, (λ∗, x∗) must be positive. And
the Jacobian matrix Dx,λF2(x, λ, t) can be get as follow
Dx,λF2(x∗, λ∗, t∗) = At∗− (m − 1)λ[|x [m−2] ∗ |] −x[m−1]∗ 2xT∗ 0 ∈ R(n+1)×(n+1).
Since A∗ > 0 and x∗ > 0, we can get At∗ > 0 directly. Let Q =
1 m−1[|x [m−2] ∗ |]−1 And Q 0 0T 1 At∗− (m − 1)λ[|x [m−2] ∗ |] −x[m−1]∗ 2xT ∗ 0 = QAt∗ − λ∗I − 1 m−1x∗ 2xT ∗ 0 = T ∈ R(n+1)×(n+1).
Q is the diagonal matrix and At∗ is positive matrix. So, QAt∗ > 0 and (λ∗, x∗) is
the positive eigenpair of QAt∗. By Lemma (3.2)(e), T is invertible.
Q 0
0T 1
is
For t = 0, A(0) = A0 and (x0, λ0, 0) is the positive solution of F2(x, λ, t) = 0.
By Lemma 6.2 and the implicit function theorem, there is a solution curve s2(t) =
(x(t), λ(t), t) ∈ Rn+1≥0 × [0, 1) that is parametrized by t for t ∈ [0, 1). We let the
following set
C2 = {(x(t), λ(t), t)|t ∈ [0, 1)}, (x(0), λ(0), 0) = (x0, λ0, 0) (6.7)
be the set of s2(t). The following theorem can be obtained.
Theorem 6.4. C2 is the solution set of H2(x, λ, t) = 0 in Rn+1≥0 × [0, 1).
Proof. If (x∗, λ∗, t∗) ∈ Rn+1≥0 ×[0, 1) is a solution of F1(x, λ, t) = 0 and (x∗, λ∗, t∗) 6=
(x(t∗), λ(t∗), t∗) ∈ C2. By Lemma 6.2 and the implicit function theorem, there is
a solution curve s∗2 = (x∗(t), λ∗(t), t) can be parameterized by t for t ∈ [0, t∗] with
(x∗(t∗), λ∗(t∗), t∗) = (x∗, λ∗, t∗). Let C2∗ be the set of s1∗(t), where
C2∗ = {(x∗(t), λ∗(t), t)|t ∈ [0, t∗]} with (x∗(t∗), λ∗(t∗), t∗) = (x∗, λ∗, t∗).
We can obtain C2∩ C2∗ = ∅ by Lemma 6.2 and the implicit function theorem. By
Lemma 6.1, C2 ⊆ Rn+1≥0 ×[0, 1). Since (x∗(0), λ∗(0), 0) is the solution F2(x, λ, t) = 0
and the rank-1 matrix A0 ∈ Rn×n>0 is irreducible, we obtained (λ∗(0), x∗(0)) is the
eigenpair of A0 and (λ∗(0), x∗(0)) = (λ0, x0) by Lemma 5.1. It is a contradiction
because C2 ∩ C2∗ = ∅.
For t ∈ [0, 1), ||x(t)|| = 1 and A(t) is bounded. Since A(t) ∈ R[m,n]≥0 is
irre-ducible and (λ(t), x(t)) ∈ Rn+1≥0 is the eigenpair of A(t) for t ∈ [0, 1). By Lemma
6.3, λ(t) = %(A(t)) and we get ||(λ(t), x(t))|| is bounded for t ∈ [0, 1). Suppose
(x∗, λ∗, 1) is an accumulation point for C2 and we know (λ∗, x∗) is the H-eigenpait
of A. The H-eigenvalues of tensors had been shown are the roots of nonzero
polynomial by Theorem 5.1. So, λ∗ is isolated and limt→1−λ(t) = λ∗. Let
Γ2 = {x ∈ Rn≥0|(x∗, λ∗, 1) is an accumulation point of C2.}. (6.8)
We get Γ2 is connected and the following theorem.
Theorem 6.5. Let A ∈ R[m,n]≥0 be irreducible. Then the nonnegative solution of
F2(x, λ, 1) = 0 is isolated.
Proof. Let x∗ ∈ Γ2 and we know that (x∗, λ∗, 1) is the nonnegative solution of
F2(x, λ, t) = 0. Since A is irreducible, we obtain (λ∗, x∗) ∈ Rn+1>0 . The Jacobian
matrix Dx,λF2(x∗, λ∗, 1) is Dx,λF2(x∗, λ∗, 1) = A1− (m − 1)λ[|x [m−2] ∗ |] −x[m−1]∗ 2xT∗ 0 ∈ R(n+1)×(n+1).
Now, we prove A1is irreducible. Suppose A1is reducible , then there is a nonempty
subset S ⊂ {1, 2, . . . , n} such that (A1)i,j = 0, ∀i ∈ S, ∀j /∈ S. We obtain
(A1)i,j = ( m X k=2 A ×2x∗. . . ×k−1x∗×k+1x∗. . . ×mx∗)ij = m X k=2 n X i2,i3,...,im=1
aii2...iik−1jiik+1...imxi2. . . xik−1xik+1. . . xim = 0,
∀ i ∈ S, ∀j /∈ S. Since x∗ ∈ Rn>0, we get
aii2...im = 0, ∀ i ∈ S.
And it implies
aii2...im = 0, ∀ i ∈ S, ∀ i2, . . . im ∈ S./
So, A is reducible, a contradiction. Next, we let Q = m−11 [|x[m−2]∗ |]−1 and we can
know Q 0 0T 1 A1− (m − 1)λ[|x [m−2] ∗ |] −x[m−1]∗ 2xT ∗ 0 = QA1− λ∗I − 1 m−1x∗ 2xT∗ 0 = T ∈ R(n+1)×(n+1).
Since Q is diagonal matrix and A1 is irreducible matrix. We obtain QA1 is
irre-ducible and (λ∗, x∗) is the positive eigenpair of QA1. By Lemma (3.2)(e), T is
invertible. Since
Q 0
0T 1
is invertible, and it implies Dx,λF2(x∗, λ∗, 1) is
invert-ible. We know that the nonnegative solution of F2(x, λ, 1) = 0 is isolated. And
limt→1−(λ(t), x(t)) = (λ∗, x∗).
Theorem 6.5 shows that given a irreducible tensor A ∈ R[m,n]≥0 and we can obtain
the unique positive H-eigenpair (λ∗, x∗) of A by tracing the solution curve with
the initial rank-1 tensor A0 and the (λ0, x0) that is defined in (6.1). So, we have
proven the Theorem 6.1.
7
Numerical experiments on tensors
In this section, the homotopy continuation method will be used to compute the positive eigenpair on the nonnegative irreducible tensors. We use the algorithm
4.1 and replace F1(x, λ, t) by F2(x, λ, t) that is defined in (6.3). The equipment of
7.1
Jacobian matrix
The difficulty is to present the Jacobian matrices of F2(x, λ, t) that are defined
in (6.4a) and (6.4b). In fact, we always translate a tensor into matrix. Not only is is more familiar to understand, but also it is more efficient to execute tensor operations. So, we can use the Kronecker product to get
Axm−1 = A
(1)(x ⊗x . . . ⊗
| {z }
m−2 times
x) ∈ Rn. (7.1)
At= Dx(A(t)xm−1) in (6.5) and we know there is a semi-symmetric As(t) ∈ R
[m,n] ≥0
such that A(t)xm−1 = As(t)xm−1. From (5.10), we can get
At= Dx(A(t)xm−1) = (m − 1)As(t) ×3x ×4 x . . . ×mx (7.2)
= (m − 1)As(1)(t)(I ⊗x . . . ⊗
| {z }
m−2 times
x). (7.3)
There is a simple way to get the semi-symmetric tensor As(t) [1]. Each element
of As(t) can be obtained by as(t)i1,i2,...,im = 1 (m − 1)! X j2,j3,...jm a(t)i1j2,...,jm, (7.4)
where (j2, j3, . . . , jm) is any permutation of (i2, i3, . . . , im).
7.2
Other algorithms for tensors
NQZ algorithm
The NQZ was developed by M. Ng, L. Qi and G. Zhou in 2009 [7]. It is like a power method with Collatz formula on tensors and it guarantees to compute the positive eigenpair for nonnegative irreducible tensors. The convergence of NQZ algorithm is linear.
Newton-Noda iteration (NNI)
We let r(x, λ) = λx[m−1]− Axm−1, f (x, λ) = −r(x, λ) 1 2(1 − x Tx) . (7.5)
The idea of NNI [1] is that using the Newton method to solve f (x, λ) = 0 and keeping the positive of each solution like Noda’s algorithm. The Jacobian matrix of f (x, λ) is Df (x, λ) = − Dxr(x, λ) x [m−1] xT 0 , (7.6)
Algorithm 5 NQZ algorithm
Step 1. For a irreducible tensor A ∈ R[m,n]≥0 , choosing x0 ∈ Rn>0. Let y0 = Axm−1
Step 2. For k = 0,1,... Compute xk+1 = ||yyk k|| yk+1= Axm−1 λ = minxi>0( yi xi) ¯ λ = maxxi>0( yi xi) Step 3. If ¯λ − λ <
stop the algorithm else
go to step 2. end
where
Dxr(x, λ) = (m − 1)λ[|x[m−2]|] − Dx(Axm−1) (7.7)
and Dx(Axm−1) is defined in (6.5). The key of the algorithm is the chosen of θ ∈
(0, 1] and we let θ = 1 in the following tests. It guarantees the global convergence and the convergence of NNI is quadratic near the solution.
Algorithm 6 NNI
Step 1. Given a unit random vector x0 ∈ Rn>0, ¯λ0 = max(
Axm−10 x[m−1]0 ), and = 10 −10 Step 2. For k = 0, 1, 2... Solve Dxr(xk, ¯λk)wk= x [m−1] k yk= wk/||wk|| Choose θk > 0 ˆ xk+1= (m − 2)xk+ θkyk xk+1= ˆxk+1/||ˆxk+1|| ¯ λk+1 = max( Axm−1 k+1 x[m−1]k+1 ), λk+1 = min( Axm−1 k+1 x[m−1]k+1 ) Step 3. If ||¯λk+1− λk+1||/¯λk+1 <
stop the algorithm else
go to Step 2. end
7.3
Compare the numerical results on tensors
Following table is the numerical results for each size of tensor. (m,n) is the order and the dimension of tensor A, respectively. ”Continuation” is the homotopy
Table 3: Random tensor
Size Continuation NQZ NNI
m n Iter Num Res Num Res Num Res
3 50 5 20 6.4874e−11 6 1.2648e−12 4 6.0991e−12
60 5 20 1.8723e−11 6 1.1783e−12 4 1.7212e−11
70 5 20 2.3022e−11 6 1.1715e−12 4 1.5664e−12
80 5 20 8.5636e−11 6 8.5676e−13 4 3.0801e−11
90 5 20 1.7804e−11 6 5.8902e−13 4 3.1113e−12
100 5 20 6.6102e−11 6 6.6187e−13 4 4.7704e−12
4 50 5 19 2.2976e−11 6 1.5369e−12 4 3.2836e−11
60 5 19 5.2278e−11 6 1.832e−12 4 4.7681e−11
70 5 19 9.1175e−11 6 1.9992e−12 4 8.5072e−11
80 5 19 1.4789e−11 6 2.5718e−12 4 9.9697e−11
90 5 19 2.2007e−11 6 3.1720e−12 4 2.9086e−11
100 5 19 2.9606e−11 6 3.1736e−12 4 2.9932e−11
continuation method. ”NQZ” is the NQZ algorithm. ”NNI” is the Newton-Noda iteration. ”Iteration” is the numbers of iteration. ”Res” is the residual. ”Num” is
the number of execution of Axm−1. If Num≥2000, we will stop the algorithm.
Example 8. Let A ∈ R[m,n] be a random tensor. We examine the numerical
results with m = 3, 4 and n = 50, 60, 70, 80, 90, 100. We choose A0 ∈ R
[m,n] >0
with all the elements of A0 are 1, x0 = √1n[1, 1, . . . 1]T ∈ Rn and λ0 = nm−1 as
the initial condition for each algorithm. Each algorithm can compute the positive eigenpair successfully. In Table 3, the homotopy continuation method can compute
the solution accurately but it costs more evaluations of Axm−1 than NQZ and NNI.
Because there are many predict and correct processes in homotopy continuation method.
Example 9. Let A ∈ R[m,n] = D + L be a signless Laplacian tensor with the
edge set E = {{i − m + 2, i − m + 3, . . . , i, i + 1}, i = m − 1, . . . , n} [11] [1] and
n+1 is defined with 1. The elements of D are zeros except the di1i2...im 6= 0 for
i1 = i2 = . . . = im. So, A is a sparse tensor. We choose A0 = nm−11 E ∈ R
[m,n] >0 ,
where all the elements of E ∈ R[m,n]>0 are 1, x0 = √1n[1, 1, . . . 1]T ∈ Rn and λ0 = 1
as the initial condition for each algorithm. We examine the size of tensor by m = 3, 4 and n = 50, 60, 70, 80, 90, 100. In Table 4, it is clear to see that the number of iterations of NQZ is more than the homotopy continuation method and NNI. The convergence of NQZ is poor. The NNI and the homotopy continuation method can compute to the positive solution successfully.
Table 4: A = D + L
Size Continuation NQZ NNI
m n Iter Num Res Num Res Num Res
3 50 11 64 1.7402e−16 1402 9.9359e−11 8 2.6009e−16
60 13 79 2.0582e−16 1971 9.9329e−11 8 4.3628e−16
70 13 79 2.6795e−16 2000 1.3644e−8 8 5.4009e−16
80 15 93 1.6959e−16 2000 3.8976e−7 8 1.6211e−12
90 15 93 1.8075e−16 2000 3.7315e−6 8 1.9874e−142
100 15 94 1.1344e−15 2000 1.8135e−6 9 7.4063e−16
4 50 9 54 2.8566e−11 757 9.881e−11 9 1.8967e−15
60 12 71 1.2655e−16 1069 9.9096e−11 9 1.2254e−14
70 12 72 6.2383e−16 1431 9.9514e−11 9 1.1820e−13
80 14 86 1.2461e−16 1842 9.9996e−11 9 4.1863e−14
90 13 82 5.1522e−12 2000 1.1576e−8 9 7.1314e−15
100 13 82 3.4720e−11 2000 2.0444e−7 10 1.5562e−15
8
Conclusion
We have proven the Perron-Frobenius theorem by a new way. By constructing
a linear homotopy with a nonnegative irreducible matrix A ∈ Rn×n, we solve it
to obtain the positive eigenpair of A and the positive eigenpair of A is sufficient to prove the Perron-Frobenius theorem. Furthermore, the skill of the proof has been extended to prove the Perron-Frobenius theorem on tensors. Both numerical results on matrices and on tensors are presented. Whether on matrices or on tensors, we can compute the positive eigenpair by the homotopy continuation method.
References
[1] C.S. Liu, C.H. Guo, and W. W. Lin, NewtonNoda iteration for finding the Perron pair of a weakly irreducible nonnegative tensor, preprint, (2017). [2] DIMACS10: DIMACS10 test set and the University of Florida Sparse Matrix
Collection. https : //sparse.tamu.edu/DIM ACS10
[3] H.B. Keller, Lectures on Numerical Methods in Bifurcation problems, 1987. [4] K.C. Chang, K. Pearson, and T. Zhang, A survey on the spectral theory of
[5] K.C. Chang, K. Pearson, and T. Zhang, Perron-Frobenius theorem for non-negative tensors, Commum. Math. Sci., 6(2008) no. 1, 507-520.
[6] K Kontovasilis, RJ Plemmons and WJ Stewart, Block Cyclic SOR for Markov Chains With p-cyclic Infinitesimal Generator, 1991.
[7] M. Ng, L. Qi, and G. Zhou, Finding the largest eigenvalue of a nonnegative tensor, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 10901099.
[8] P. Chanchana, An algorithm for computing the Perron root of a nonnegative irreducible matrices, 2007.
[9] Q. Ni and L. Qi, A quadratically convergent algorithm for finding the largest eigenvalue of a nonnegative homogeneous polynomial map, J. Global Optim., 61 (2015) pp. 627641
[10] R. A. Horn, and C. R. Johnson, (2012)Matrix analysis, Cambridge University Press.
[11] S. Hu, L. Qi and J. Xie, The largest Laplacian and signless Laplacian H-eigenvalues of a uniform hypergraph, Linear Algebra Appl., 469 (2015), pp. 127.
[12] S. Hu, Z.-H. Huang, C. Ling, and L. Qi, On determinants and eigenvalue theory of tensors, J. Symb. Comput., 50 (2013), pp. 508531.
[13] T.G. Kolda and B.W. Bader, Tensor Decompositions and Applications, SIAM Review, vol. 51, no. 3, pp. 455-500, September 2009.
[14] Y.C. Kuo, W.W. Lin and C.S. Liu, Continuation Method For Computing Z-/H-Eigenpairs of Nonnegative Tensors, arXiv preprint arXiv:1702.05841.
[15] Y. Wei, and W. Ding, Theory and Computation of Tensors :
