Computing tensor Z-eigenpairs via an alternating direction method

Problem formulation

An mth order n-dimensional real tensor consisting of n^m entries in ℝ:

$\mathcal{A}$ is called symmetric if the value of a_i1i2⋯im is invariant under any permutation of its indices i₁, i₂, …, i_m. For convenience, we use ${\mathbb{S}}^{\left[m,n\right]}$ to denote the set of all m th order n-dimensional real symmetric tensors.

Throughout this article, we use $\mathcal{A}{x}^{m-k}\left(0\le k\le m\right)$ to simply denote a k th order n-dimensional tensor defined by (1)${\displaystyle {\left(\mathcal{A}{x}^{m-k}\right)}_{i_1{i}_2\cdots i_k}=\sum _{i_{k+1},\dots ,i_m=1}^n{a}_{i_1{i}_2\cdots i_k{i}_{k+1}\cdots i_m}{x}_{i_{k+1}}\cdots x_{i_m},}$

for all 1 ≤ i₁, i₂, …, i_k ≤ n and 1 ≤ k ≤ m.

Obviously, $\mathcal{A}{x}^m$ is a scalar, $\mathcal{A}{x}^{m-1}$ is a vector, and $\mathcal{A}{x}^{m-2}$ is a matrix. For brevity, let $\mathcal{A}{x}^p{y}^q$ denote the result of $\left(\mathcal{A}{x}^p\right)y^q=\left(\mathcal{A}{x}^{m-\left(m-p\right)}\right)y^{m-p-\left(m-p-q\right)}$, where 0 ≤ p, q, p + q ≤ m, and $\mathcal{A}{x}_1^{p_1}{x}_2^{p_2}\cdots x_2^{p_k}$ can be computed in a similar way, where 0 ≤ p₁, p₂, …, p_k ≤ m, and k is an arbitrary integer such that 0 ≤ p₁ + p₂ + ⋯ + 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 $\mathcal{A}$, i.e., $f\left(x\right)=\mathcal{A}{x}^m$. We call $\mathcal{A}$ positive definite tensor if $\mathcal{A}{x}^m>0$ for all x ∈ ℝⁿ∖{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 $\mathcal{A}$ be an m th order n-dimensional symmetric real tensor. If there exists a nonzero vector x ∈ ℝⁿ and a scalar λ ∈ ℝ satisfying (2)${\displaystyle \left\{\begin{array}{c}\mathcal{A}{x}^{m-1}=\lambda x,\hfill \\ x^{T}x=1,\hfill \end{array}\right.}$

then we call the scalar λ a Z-eigenvalue of $\mathcal{A}$, and the vector x a Z-eigenvector associated with the Z-eigenvalue λ. We also say the pair (λ, x) is a Z-eigenpair of $\mathcal{A}$.

Iterative algorithms to find the largest or smallest eigenvalues and corresponding eigenvectors are usually designed to solve a nonlinearly constrained optimization problem

${\displaystyle maxf\left(x\right)=\mathcal{A}{x}^m}$ (3)${\displaystyle s.t.x\in {\mathbb{S}}^{n-1},}$

where 𝕊ⁿ⁻¹ denotes the unit sphere in the Euclidean norm, i.e., ${\mathbb{S}}^{n-1}=\left\{x\in {\mathbb{R}}^n|{∥x∥}^2=1\right\}$. We can determine the gradient and Hessian of the objective function of Eq. (3) through some simple calculations, as follows: (4)${\displaystyle g\left(x\right)\equiv \nabla f\left(x\right)=m\mathcal{A}{x}^{m-1}}$

and (5)${\displaystyle H\left(x\right)\equiv {\nabla }^{2}f\left(x\right)=m\left(m-1\right)\mathcal{A}{x}^{m-2}}$

Some existing methods for Z-eigenproblems

In this subsection, we introduce some typical methods for computing Z-eigenpairs by solving the problem Eq. (3) or its variants. From Theorem 3.2 of Kolda & Mayo (2011), we know that $\left(\lambda ,x\right)$ is a Z-eigenpair of $\mathcal{A}$ if and only if x is a constrained stationary point of Eq. (3) and $\lambda =\mathcal{A}{x}^m/∥x∥$. Based on the theorem, De Lathauwer, De Moor & Vandewalle (2000) proposed the S-HOPM method for solving the problem Eq. (3) to find the best symmetric rank-1 approximation of a symmetric tensor $\mathcal{A}\in {\mathbb{S}}^{\left[m,n\right]}$, which is equivalent to finding the largest Z-eigenvalue of $\mathcal{A}$ (Qi, 2005). The main step of the S-HOPM algorithm is (6)${\displaystyle x_{k+1}=\frac{\mathcal{A}{x}_k^{m-1}}{∥\mathcal{A}{x}_k^{m-1}∥},{\lambda }_{k+1}=\mathcal{A}{x}_{k+1}^m.}$

Under the assumption of convexity on $\mathcal{A}{x}^m$, S-HOPM could be convergent for even-order 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 to (7)${\displaystyle \hat{f}\left(x\right)=\mathcal{A}{x}^m+\alpha {∥x∥}^m,}$

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 $\mathcal{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.

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 $\mathcal{A}{x}^m$ with even number m is rewritten as $\mathcal{A}{x}_1^2{x}_2^2\cdots x_p^2$, 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:

${\displaystyle max\tilde{f}\left(x\right)=\mathcal{A}{x}_1^2{x}_2^2\cdots x_p^2}$ (8)${\displaystyle s.t.x_i=x_j,i,j=1,2,\dots ,p}$ ${\displaystyle x_i\in {\mathbb{S}}^{n-1},i=1,2,\dots ,p.}$

When $\mathcal{A}$ is symmetric and conditions x_i = x_j = x(i, j = 1, 2, …, p) hold, we can obtain $\mathcal{A}{x}_1^2{x}_2^2\cdots x_p^2=\mathcal{A}{x}^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 $\mathcal{A}{x}_1^2\cdots x_{i-1}^2{x}_{i+1}^2\cdots x_p^2$.

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 $\mathcal{A}{x}_1^2\cdots x_{i-1}^2{x}_{i+1}^2\cdots x_p^2$. 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.

The main computational cost lies in tensor-vector multiplications $\mathcal{A}{x}_1^2\cdots x_{i-1}^2{x}_{i+1}^2\cdots x_p^2$ and matrix eigenvalue computations. For an mth order n-dimensional symmetric tensor $\mathcal{A}$, it costs O(mn^m) operations to compute the matrix $\mathbf{A}=\mathcal{A}{x}_1^2\cdots x_{i-1}^2{x}_{i+1}^2\cdots x_p^2$. The cost of computing the largest or smallest eigenvalue of the matrix A is (4/3)n³, 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₁ and x₂, and the inner iteration can be omitted because p = 1. According to the symmetry property of $\mathcal{A}$, we also have $\mathcal{A}{x}_1^2{x}_2^2=\mathcal{A}{x}_2^2{x}_1^2$. Therefore, it is not necessary to explicitly write out the variable x₂, 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.

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 $\mathcal{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 $\mathcal{A}$, x_k+1 is the largest eigenvalue of the matrix $\mathcal{A}{x}_k^2$. Therefore, the quadratic function $q\left(y\right)=y^T\left(\mathcal{A}{x}_k^2\right)y=\mathcal{A}{x}_k^2{y}^2$ reaches a maximum value λ_k+1 at the point y = x_k+1 over the unit sphere 𝕊ⁿ, i.e., ${\lambda }_{k+1}=\mathcal{A}{x}_k^2{x}_{k+1}^2\ge \mathcal{A}{x}_k^2{y}^2$ for all y ∈ 𝕊ⁿ. At the same time, x_k is the largest eigenvalue of the matrix $\mathcal{A}{x}_{k-1}^2$. These results give (9)${\displaystyle {\lambda }_{k+1}=\mathcal{A}{x}_k^2{x}_{k+1}^2\ge \mathcal{A}{x}_k^2{x}_{k-1}^2=\mathcal{A}{x}_{k-1}^2{x}_k^2={\lambda }_k.}$

Here, the second equality holds because of the symmetric property of $\mathcal{A}$. From Eq. (9), we know that the sequence λ_k generated by Algorithm 3 is nondecreasing. On the other side, λ_k is computed by ${\lambda }_k=\mathcal{A}{x}_{k-1}^2{x}_k^2$, where x_k ∈ 𝕊ⁿ. Due to the compactness of the unit sphere 𝕊ⁿ, 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 $\mathcal{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 $\mathcal{A}{x}_k^2{x}_{k+1}={\lambda }_k{x}_{k+1}$, $\mathcal{A}{x}_{k+1}^2{x}_k={\lambda }_{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.