Computing tensor Z-eigenpairs via an alternating direction method

Tensor eigenproblems have wide applications in blind source separation, magnetic resonance imaging, and molecular conformation. In this study, we explore an alternating direction method for computing the largest or smallest Z-eigenvalue and corresponding eigenvector of an even-order symmetric tensor. The method decomposes a tensor Z-eigenproblem into a series of matrix eigenproblems that can be readily solved using off-the-shelf matrix eigenvalue algorithms. Our numerical results show that, in most cases, the proposed method converges over two times faster and could determine extreme Z-eigenvalues with 20–50% higher probability than a classical power method-based approach.

Unlike for matrices, there are several definitions of tensor eigenvalues and corresponding eigenvectors. For example, Qi (2005) proposed the definition of an H-eigenvalue and Z-eigenvalue as being equivalent to the l m -eigenvalue and l 2 -eigenvalue in Lim (2005), respectively. In Chang, Pearson & Zhang (2009), these definitions were unified by employing a positive definite tensor B while m was even. In this work, we mainly focus on computing Z-eigenvalues of symmetric tensors.
In general, the calculation of all eigenvalues of a higher-order tensor is very difficult due to the NP-hardness of deciding tensor eigenvalues over R (Hillar & Lim, 2013). Fortunately, one only needs to compute the largest or smallest eigenvalue of a tensor in certain scenarios. For instance, to guarantee the positive definiteness of the diffusivity function in higher-order diffusion tensor imaging, we just need to compute the smallest Z-eigenvalue of the tensor and make sure it is nonnegative. In automatic control (Ni, The main contributions of this work are listed as follows: • We reformulate a typical tensor Z-eigenvalue problem as an equivalent multi-variable linearly constrained problem using the variable splitting method, which has a structure similar to the matrix eigenvalue problem.
• We design an efficient algorithm to solve the reformulated problem, in which the tensor Z-eigenproblem is decomposed into a series of matrix eigenproblems that has been extensively studied and can be readily solved using off-the-shelf matrix eigenvalue algorithms.
The remainder of the article is organized as follows. In the next section, we formulate the tensor eigenproblem and introduce some classical methods for solving it. In Section 3, we propose a simple and efficient algorithm for the problem and analyze its convergence property. Section 4 reports some experimental results to show the efficiency of our proposed method. Finally, we conclude this article in Section 5.

Problem formulation
An mth order n-dimensional real tensor consisting of n m entries in R: A is called symmetric if the value of a i 1 i 2 ···i m is invariant under any permutation of its indices i 1 ,i 2 ,...,i m . For convenience, we use S [m,n] to denote the set of all m th order n-dimensional real symmetric tensors.
Throughout this article, we use Ax m−k (0 ≤ k ≤ m) to simply denote a k th order n-dimensional tensor defined by for all 1 ≤ i 1 ,i 2 ,...,i k ≤ n and 1 ≤ k ≤ m. Obviously, Ax m is a scalar, Ax m−1 is a vector, and Ax m−2 is a matrix. For brevity, let Ax p y q denote the result of can be computed in a similar way, where 0 ≤ p 1 ,p 2 ,...,p k ≤ m, and k is an arbitrary integer such that 0 ≤ p 1 + p 2 + ··· + p k ≤ m.
Using the definition of Eq. (1), an mth degree homogeneous polynomial function f (x) with real coefficients can be represented by a symmetric tensor A, i.e., f (x) = Ax m . We call A positive definite tensor if Ax m > 0 for all x ∈ R n \ {0}. It is easy to understand that m must be even in this case.
As mentioned before, there are several definitions of tensor eigenvalues and corresponding eigenvectors. In this work, we mainly focus on computing Z-eigenvalues of symmetric tensors defined as follows.

Definition 1 (Qi, 2005). Let
A be an m th order n-dimensional symmetric real tensor. If there exists a nonzero vector x ∈ R n and a scalar λ ∈ R satisfying then we call the scalar λ a Z-eigenvalue of A, and the vector x a Z-eigenvector associated with the Z-eigenvalue λ. We also say the pair (λ,x) is a Z-eigenpair of A. Iterative algorithms to find the largest or smallest eigenvalues and corresponding eigenvectors are usually designed to solve a nonlinearly constrained optimization problem where S n−1 denotes the unit sphere in the Euclidean norm, i.e., S n−1 = {x ∈ R n | x 2 = 1}. We can determine the gradient and Hessian of the objective function of Eq. (3) through some simple calculations, as follows: and

Some existing methods for Z-eigenproblems
In this subsection, we introduce some typical methods for computing Z-eigenpairs by solving the problem Eq.  (2000) proposed the S-HOPM method for solving the problem Eq.
(3) to find the best symmetric rank-1 approximation of a symmetric tensor A ∈ S [m,n] , which is equivalent to finding the largest Z-eigenvalue of A (Qi, 2005). The main step of the S-HOPM algorithm is Under the assumption of convexity on Ax m , S-HOPM could be convergent for evenorder tensors. However, it has been pointed out that S-HOPM can not guarantee to converge globally (Kofidis & Regalia, 2002). To address this issue, Kolda & Mayo (2011) modified the objective function tô and proposed the SS-HOPM for solving Eq. (3) with the objective function Eq. (7). SS-HOPM has a similar iterative scheme to S-HOPM, but at the same time has a shortcoming in the choice of the shift α. To overcome the limitation, the same authors proposed an adaptive method, called GEAP, which is monotonically convergent and much faster than the SS-HOPM method due to its adaptive shift choice of the shift. GEAP is originally designed to calculate generalized eigenvalues (Chang, Pearson & Zhang, 2009) with a positive definite tensor B. The authors also presented a specialization of the method to the Z-eigenvalue problem, which is equivalent to SS-HOPM except for the adaptive shift. The details of the GEAP specialization are briefly summarized in Algorithm 1.
Algorithm 1. GEAP method for the problem Eq.
(3) with the objective function Eq. (7) Initialization: Given a tensor A ∈ S [m,n] , an initial vector x 0 ∈ R n , and a tolerance > 0. Let β = 1 if we want to compute the largest Z-eigenvalue, and let β = −1 if we want to compute the smallest Z-eigenvalue. Let τ be the tolerance on being positive/negative definite. 1: GEAP is a simple and effective approach for computing Z-eigenvalues of a symmetric tensor, but it is not guaranteed to determine the largest eigenvalue or the smallest one, which is exactly the goal in some applications. To obtain these extreme eigenvalues with a higher probability, we propose to reformulate the problem Eq. (3) to make its structure similar to a matrix eigenproblem, so that it can be solved by existing methods for matrices that always converge to extreme eigenvalues.

An alternating direction method for Z-eigenproblems
Motivated by that algorithms for solving matrix eigenproblem can always converge to the largest or smallest eigenvalue, we propose to compute extreme eigenvalues by combining the method of solving the matrix eigenproblem and tensor optimization techniques. To this end, we adopt a variable splitting strategy in which we introduce some superfluous variables and equality constraints over these variables. Specifically, the term Ax m with even number m is rewritten as Ax 2 1 x 2 2 ···x 2 p , where p = m/2, with the equality constraints x i = x j (i,j = 1,2,...,p). Therefore, problem Eq. (3) is transformed into the following model: When A is symmetric and conditions x i = x j = x(i,j = 1,2,...,p) hold, we can obtain Ax 2 1 x 2 2 ···x 2 p = Ax m . Using this fact, the equivalence between problems Eqs. (3) and (8) can be easily checked. It is also worthwhile to note that if all variables except x i are available and those equality constraints are not considered, the problem Eq. (8) reduces to the standard matrix eigenproblem for the matrix Ax 2 1 ···x 2 i−1 x 2 i+1 ···x 2 p . Directly solving the problem Eq. (8) may be inefficient because its special structure is not considered, and in doing so, it is easy to converge to a locally optimal point, thus the largest or smallest eigenvalue could not be determined. On the other hand, it is comparatively easy to compute extreme eigenvalues for the matrix cases. In Eq. (8), if all variables except x i are known and those equality constraints are not considered, solving Eq. (8) can exactly get the largest eigenvalue and the corresponding eigenvector of the matrix Ax 2 1 ···x 2 i−1 x 2 i+1 ···x 2 p . Following this observation and the fact that there are many efficient algorithms available for tackling matrix eigenproblems, we propose a simple alternating direction scheme between solving different matrix eigenproblems for the problem Eq. (8). The details of this method are given in Algorithm 2.
Algorithm 2 Alternating direction method (ADM) for Eq. (8) Initialization: Given an even-order tensor A ∈ S [m,n] , initial unit vectors x i ∈ R n ,i =1 ,2,...,p, where p = m/2, and > 0 is the tolerance. Set x = x p , λ = Ax m , and δ as the absolute difference between successive values of λ.

3: Update the variable
End for 4: Set x = x p and λ = Ax m . End while Output: Z-eigenvalue λ and its associated Z-eigenvector x.
The main computational cost lies in tensor-vector multiplications Ax 2 1 ···x 2 i−1 x 2 i+1 ···x 2 p and matrix eigenvalue computations. For an mth order n-dimensional symmetric tensor A, it costs O(mn m ) operations to compute the matrix A = Ax 2 1 ···x 2 i−1 x 2 i+1 ···x 2 p . The cost of computing the largest or smallest eigenvalue of the matrix A is (4/3)n 3 , which is much less than the products for large n. Compared with GEAP, which has similar computations to our method but needs to compute tensor multiplication at least twice in each iteration, our method can save about half of the time because it only needs to calculate tensor multiplication once in each iteration. This can be verified in the numerical experiments in Section 4.

Specialization of ADM to fourth-order tensors
The proposed ADM transforms the tensor eigenvalue problem Eq. (3) into a series of matrix eigenvalue problems that are easy to solve. For fourth-order tensors, there are two related variables of x 1 and x 2 , and the inner iteration can be omitted because p = 1. According to the symmetry property of A, we also have Ax 2 1 x 2 2 = Ax 2 2 x 2 1 . Therefore, it is not necessary to explicitly write out the variable x 2 , and the procedure of Algorithm 2 can be simply described, as shown in Algorithm 3, for fourth-order tensors. To better describe the iterative steps, we use x k to denote a kth iterate in Algorithm 3, rather than the splitting variable as in Algorithm 2.
Algorithm 3 Specialization of the ADM to fourth-order tensors Initialization: Given a tensor A ∈ S [4,n] , initial unit vectors x 0 ∈ R n , and > 0 is the tolerance. Set λ 0 = Ax m 0 , k := 0, and δ as the absolute difference between successive values of λ.
For k = 0,1,2··· do 1: Compute the matrix A k = Ax 2 k . 2: Find the largest or smallest eigenvalueλ and the corresponding unit eigenvector v of A k using any eigenvalue algorithm for matrices.

Convergence analysis
As shown in the main steps of Algorithm 2 and Algorithm 3, the equality constraints in Eq. (8) are not considered in the process of calculation. A natural question arises about whether the algorithms can converge, and furthermore, whether the algorithms can converge to a Z-eigenvalue of A. In this subsection, we handle these issues using properties of extreme eigenvalues and corresponding eigenvectors of matrices. For simplicity, only the convergence property of Algorithm 3 is analyzed. The convergence property of Algorithm 2 can be analyzed in a similar way.
Let x k denote the kth iterate generated by Algorithm 3. According to Steps 2 and 3, in the case of computing the largest Z-eigenvalue of A, x k+1 is the largest eigenvalue of the matrix Ax 2 k . Therefore, the quadratic function q y = y T Ax 2 k y = Ax 2 k y 2 reaches a maximum value λ k+1 at the point y = x k+1 over the unit sphere S n , i.e., λ k+1 = Ax 2 k x 2 k+1 ≥ Ax 2 k y 2 for all y ∈ S n . At the same time, x k is the largest eigenvalue of the matrix Ax 2 k−1 . These results give Here, the second equality holds because of the symmetric property of A. From Eq. (9), we know that the sequence λ k generated by Algorithm 3 is nondecreasing. On the other side, λ k is computed by λ k = Ax 2 k−1 x 2 k , where x k ∈ S n . Due to the compactness of the unit sphere S n , we also know that the sequence λ k is bounded above. Consequently, λ k has a unique limit, and we can readily conclude by posing this as a theorem.
Theorem 1. Let λ k be a sequence generated by Algorithm 3. Then the sequence λ k is nonincreasing and there exists λ * such that λ k → λ * .
While Theorem 1 ensures that Algorithm 3 always terminates in finitely many iterations, theoretically, it cannot ensure that the sequence λ k converges to a Z-eigenvalue of A because the equality constraints in Eq. (8) are omitted in the implementation of Algorithm 3. One possible result is the occurrence of cyclic solutions, that is, two consecutive iterates x k and x k+1 that satisfy Ax 2 k x k+1 = λ k x k+1 , Ax 2 k+1 x k = λ k+1 x k , and λ k = λ k+1 . However, this situation is rarely encountered in the numerical experiments presented in the next section. Additionally, how to theoretically avoid this situation is the subject for future research.

NUMERICAL EXPERIMENTS
In this section, we present some numerical results of the ADM for computing the largest or smallest Z-eigenvalues of tensors. The proposed ADM is compared with the GEAP method, which is an adaptive shifted power method first proposed by Kolda & Mayo (2014), and the AG method (Yu et al., 2016), which is an adaptive gradient method with inexact stepsize. All experiments are performed in MATLAB R2017a and the Tensor Toolbox (Bader & Kolda, 2012) under a Windows 10 operating system on a laptop with an Intel(R) Core (TM) i7-10510U CPU and 12 GB RAM. In all numerical experiments, we terminate the computation when the absolute difference between successive eigenvalues is less than 10 −10 , i.e., |λ k+1 − λ k | ≤ 10 −10 , or the number of iterations exceeded the maximum number 500.
In our experiments, we use some typical examples from references (Cui, Dai & Nie, 2014;Kofidis & Regalia, 2002;Nie & Wang, 2014) to assess the performance of the proposed method in finding the largest or smallest Z-eigenvalue of a symmetric tensor. All of the largest or smallest Z-eigenvalues in these examples are given in the original literature. Therefore, those desired values are known in advance. The largest and smallest Z-eigenvalue of the tensor A are respectively We first test the convergence performance of the proposed ADM in comparison to GEAP. Figure 1 shows the convergence trajectories of the two methods for computing extreme Zeigenvalues of A from Example 1, with the starting point x 0 = (0.0417,−0.5618,0.6848) T . As shown on the left in Fig. 1, both GEAP and ADM can find the largest Z-eigenvalue 0.8893. Until the stopping criterion is met, GEAP runs 63 iterations in 0.2188 s, while the proposed ADM runs only 26 iterations in 0.0313 s. When computing the smallest Z-eigenvalue with the same starting point (right in Fig. 1), although ADM runs longer than GEAP, the ADM find the desired value of −1.0954, while GEAP fail.
To test the performance of the proposed algorithm in finding extreme eigenvalues, 1,000 random starting guesses drawn uniformly from [-1, 1] are employed in the experiments. All three methods are implemented 1,000 times, each with the same initial point. We list the number of occurrences of extreme eigenvalues (#Occu.), the average number of iterations (Iter ave ), and the average running time in seconds (CPU ave ) for the two types of Z-eigenvalues in Tables 1, 2, 3, 4, 5 and 6. As shown in Tables 1-6, for the cases of computing the largest eigenvalues, the proposed ADM runs a similar or slightly larger number of iterations but in a similar or shorter time compared to both GEAP and AG. For the case of computing the smallest Z-eigenvalue, ADM runs slower for Example 1 but runs much faster for the other examples. It can also be seen from the fourth columns of Tables 1-6 that GEAP can only obtain the largest Z-eigenvalues with a probability of about 0.55, and the smallest Z-eigenvalue with a probability of about 0.65. AG performs clearly better than GEAP. The proposed ADM performs best in almost all cases, which can reach In the case of n = 5, the largest and smallest Z-eigenvalues of the tensor A are respectively Example 3 (Cui, Dai & Nie, 2014). Consider the symmetric tensor A ∈ S [4,n] such that a ijkl = tan(i) + tan j + tan(k) + tan(l),1 ≤ i,j,k,l ≤ n.
In the case of n = 5, we can obtain the largest and smallest Z-eigenvalues of the tensor A from the reference as follows. Example 4 (Nie & Wang, 2014). Let A ∈ S [4,n] be a symmetric tensor with a ijkl = arctan (−1) i i n + arctan (−1) j j n + arctan (−1) k k n + arctan (−1) l l n .   In the case of n = 5, the largest and smallest Z-eigenvalues of the tensor A are respectively Example 5 (Nie & Wang, 2014). Let A ∈ S [4,n] be a symmetric tensor with  Besides the point λ 1 = 9.5821 of the largest Z-eigenvalue, the tensor A has another stable eigenvalue λ 2 = 0. As shown in Figure 2 and Table 5, GEAP falls into the latter point with a probability of around 0.4, while ADM can always converge to the previous point.
Example 6 (Sheng & Ni, 2021 and a i 1 ···i 6 =0 if (i 1 ···i 6 ) is not a permutation of an index in the above. We can get the largest Z-eigenvalue λ max = 1 and smallest Z-eigenvalue λ min = 0, and these corresponding eigenvectors are not unique. The comparison results are shown in Table 6, from which we find that GEAP converges very slowly when computing the smallest eigenvalue of A.
In contrast, ADM reaches the smallest eigenvalue with only about 12 iterations for each execution.

DISCUSSION
In general, algorithms for computing the largest or smallest eigenvalue of a higher-order tensor are prone to getting stuck in local extrema and then converging to an arbitrary eigenvalue of the tensor depending on the initial conditions. However, the counterpart for symmetric matrices can always converge to the largest or smallest one. Motivated by this, we propose combining algorithms for matrix eigenproblem and tensor optimization techniques in order to obtain extreme eigenvalues. Specifically, the tensor eigenproblem is split into a series of matrix eigenvalue problems using a variable splitting method, and then an alternating scheme is proposed to solve the problem.
To solve the tensor eigenproblems, many algorithms for matrix eigenproblems are extended to the tensor case. However, these generalizations cannot guarantee the low complexity of these algorithms, and the global convergence to the extreme eigenvalues is also not ensured. In this article, the tensor eigenvalue problem is directly transformed into a series of matrix eigenvalue problems so that its algorithms can be directly used to solve the original tensor eigenvalue problem. This method not only overcomes local minima problems existing in direct generalizations, but also has great potential to speed up the convergence. The experimental results verify the effectiveness and advancement of the proposed algorithm, which converges rapidly in most cases and reaches extreme Z-eigenvalues with a significantly higher probability. In many cases, we determine the extreme eigenvalue with a probability of 1, indicating that we can obtain the extreme eigenvalue under any given initial value in these cases. This demonstrates the significant robustness of the proposed method.
However, the proposed algorithm cannot guarantee global convergence for each type of tensor. For Examples 1 and 6, we could only obtain the largest eigenvalue with a probability of about 0.6. In future research, the question of why this kind of tensor cannot obtain the extreme eigenvalue under any initial point will be discussed in more detail.

CONCLUSION
In this article, we transform a tensor Z-eigenvalue problem into a series of matrix eigenvalue problems using a variable splitting method and propose an alternating scheme for computing the largest or smallest Z-eigenvalue of symmetric tensors. Just like the classical power method, which constantly uses the intermediate iterates to construct a vector, the proposed algorithm uses them to construct a matrix and computes the eigenvalues and corresponding eigenvectors of the matrix. We analyze the convergence properties of the proposed method which is verified in the numerical experiments. The limitations of this work are twofold. First, as the authors note themselves, it can only ensure the convergence of eigenvalues, but not the convergence of eigenvectors due to the possible existence of cyclic solutions. Second, as can be seen from the numerical examples, it cannot obtain the largest eigenvalues with a 100% probability. The numerical results are reported for some testing examples, which showed that the proposed method converged much faster than both GEAP and AG in most cases and could reach the extreme Z-eigenvalues with a significantly higher probability than GEAP.