An accelerated conjugate gradient method for the Z-eigenvalues of symmetric tensors

1.   Introduction

Since Lim ^[1] and Qi ^[2] presented some theories of the eigenvalues and eigenvectors for higher order tensors, the related research has received much more attention (see ^{[3,4,5,6,7,8,9,10,11,12,13]}, etc). The eigenvalues of symmetric tensors have been applied in blind source separation ^[2] and hypergraph theory ^[9], statistical data analysis ^[14] and high order Markov chains ^[15], etc. Moreover, various definitions of eigenvalues and eigenvectors for tensors have been introduced ^[2,10,16].

There are many works for computing the eigenvalues of tensors, especially for Z-eigenvalues. Qi et al. ^[17] proposed an elimination method for finding all Z-eigenvalues, which is specific to the third-order tensors. Kolda and Mayo ^[6] presented a shifted power method (SPM) for calculating Z-eigenvalues, in which the shifted parameter is crucial. Han ^[18] provided an unconstrained optimization approach for even order symmetric tensors. Hao, Cui and Dai ^[4] found the extreme Z-eigenvalues and corresponding Z-eigenvectors by the sequential subspace projection method. Under certain assumptions, the global convergence and linear convergence were established for symmetric tensors. Hao, Cui and Dai ^[19] proposed a feasible trust region method for finding the extreme Z-eigenvalues of symmetric tensors. The global convergence and local quadratic convergence were established.

Inspired by the idea of improved conjugate parameters proposed in the above works and the application of optimization methods in tensor eigenvalue calculation, our main work is to study the transformation of Z-eigenvalues of symmetric tensors into unconstrained optimization and to propose a new algorithm. The contributions of this article are listed as follows:

For the case of different critical point, we transform the Z-eigenvalues of symmetric tensors into different unconstrained optimization problems which include shifted problem.

We propose a new conjugate gradient method with a new conjugate gradient parameter and an accelerated parameter, which converges to a critical point. The found nonzero critical point is a Z-eigenvector associated with a Z-eigenvalue of symmetric tensors. When the zero critical point is obtained, a shifted problem is solved for finding a Z-eigenvalue.

The global convergence of new method is established. We compare our method with conjugate gradient methods proposed in ^[20,21], for computing the Z-eigenvalues of symmetric tensors. The numerical results show that the proposed method is competitive.

The rest of this paper is organized as follows. In Section 2, we transform the Z-eigenvalues problem into unconstrained optimizations, and propose an accelerated conjugate gradient method for solving it. Global convergence result is established in Section 3. Numerical experiments are shown in Section 4.

2.   New method for the Z-eigenvalues of symmetric tensors

Let $\mathbb{R}$ be the real field, $m, n$ be positive integers and

be an $m$th-order $n$-dimensional real tensor. The set of $m$th-order $n$-dimensional real tensor is denoted by $\mathbb{R}^{[m, n]}$. Tensor $\mathcal{A}$ is symmetric if its entries are invariant under any permutation of their indices. The set of $m$th-order $n$-dimensional real symmetric tensor is denoted by $\mathbb{S}^{[m, n]}$.

If $\mathcal{A} \in \mathbb{R}^{[m, n]}$, an $m$th-degree homogeneous polynomial function with real coefficients is uniquely determined by

For $x = (x_1, \cdots, x_n)^{\rm T} \in \mathbb{R}^n$, $\mathcal{A}x^{m-1}$ denotes a $n$-dimensional column vector, i.e.,

If $\mathcal{A}$ is symmetric, then the gradient of $\mathcal{A}x^{m}$ satisfies

for all $x \in \mathbb{R}^n$. In our work, we consider the following Z-eigenvalues of symmetric tensors.

Definition 1. ^[2] Let $\mathcal{A} \in \mathbb{R}^{[m, n]}$, if there exist $\lambda \in \mathbb{R}$ and a vector $x \in \mathbb{R}^n\setminus \{0\}$ satisfying

then $\lambda$ is called a Z-eigenvalue of $\mathcal{A}$ and $x$ is called the corresponding Z-eigenvector.

Motivated by the work of Auchmuty ^[22], we generalize the unconstrained variational principles to Z-eigenvalues of $\mathcal{A}$. Consider the following unconstrained optimization

The gradient and Hessian of $f(x)$ are listed as follows

where

Obviously, $G(x)$ is a symmetric matrix. In order to research the properties of $f(x)$ in (2.2), we cite the following definition and a nice feature of it.

Definition 2. ^[23] A continuous function $h: \mathbb{R}^n \rightarrow \mathbb{R}$ is called coercive if it satisfies

If $x$ satisfies the equation $\nabla h(x) = 0$, then it is termed as a critical point of $h(x)$.

Theorem 1. ^[23] Let $h: \mathbb{R}^n \rightarrow \mathbb{R}$ be continuous. If $h$ is coercive, then $h$ has at least one global minimizer. In addition, if the first partial derivatives exist on $\mathbb{R}^n$, then $h$ attains its global minimizers at its critical points.

Based on a similar argument, for the Z-eigenvalues of tensors, we have the following result.

Theorem 2. Let $\mathcal{A} \in \mathbb{R}^{[m, n]}$ be symmetric tensors. Assume that $\lambda_{\max}$ is the largest Z-eigenvalue of $\mathcal{A}$. Denote the Z-spectrum of $\mathcal{A}$ by $\sigma_{Z}(\mathcal{A}): = \{\lambda: \lambda \; {\text{is a Z-eigenvalue of}}\; \mathcal{A}\}.$ We have

(ⅰ) $f(x)$ is coercive on $\mathbb{R}^n$.

(ⅱ) The critical points of $f(x)$ are at $x = 0$ and any Z-eigenvector $x\neq 0$ associated with a Z-eigenvalue $\lambda > 0$ of $\mathcal{A}$ satisfying $\lambda = (x^{\rm T}x)^{m-1}$.

(ⅲ) If $\lambda_{\max} > 0$, then $f(x)$ attains its global minimal value

at any Z-eigenvector associated with the Z-eigenvalue $\lambda_{\max}$ such that $\lambda_{\max} = (x^{\rm T}x)^{m-1}$.

(ⅳ) If $\lambda_{\max} \leq 0$, then $x = 0$ is the unique critical point of $f(x)$. Moreover, it is the unique global minimizer of $f(x)$ on $\mathbb{R}^n$.

Proof. (ⅰ) Since

and $\mathcal{A}x^{m}$ is an $mth$-degree homogeneous polynomial function with real coefficients, then

That is, $f(x)$ is coercive on $\mathbb{R}^n$.

(ⅱ) From the definition of the critical point of $f(x)$, we have

It is obvious that $x = 0$ is a critical point of $f(x)$ as $g(0) = 0$. The point $x\in \mathbb{R}^n\backslash\{0\}$ satisfying (2.5) is a Z-eigenvector corresponding to the Z-eigenvalue $\lambda = (x^{\rm T}x)^{m-1} > 0$, which is also a critical point of $f(x)$.

(ⅲ) At the critical point $x\in \mathbb{R}^n\backslash\{0\}$, a Z-eigenvector $x$ associated with a Z-eigenvalue $\lambda$ satisfies $\lambda = (x^{\rm T}x)^{m-1}$ and $\mathcal{A}x^{m} = \lambda x^{\rm T}x.$ Moreover,

By Theorem 2.1 and conclusion (ⅱ), we get the global minimum value $-\frac{1}{2m}(\lambda_{\max})^{\frac{m}{m-1}}$ at any Z-eigenvector $x$ corresponding to the Z-eigenvalue $\lambda_{\max}$ such that $\lambda_{\max} = (x^{\rm T}x)^{m-1}$.

(ⅳ) Since $\lambda_{\max} \leq 0$ implies that $\lambda\leq 0$ for any $\lambda \in \sigma_{Z}(\mathcal{A})$, $\lambda = (x^{\rm T}x)^{m-1}$ does not hold for any Z-eigenvector $x$ associated with a Z-eigenvalue $\lambda$, as $(x^{\rm T}x)^{m-1} > 0$ for any $x\in \mathbb{R}^n\backslash\{0\}$. Therefore, from Theorem 2.1, $x = 0$ is the unique critical point and the unique global minimize of $f(x)$.

Note that, if $\lambda_{\max} \leq 0$, then $x = 0$ is the unique critical point, but it does not result in a Z-eigenvalue. In this case, we solve a following shifted problem

where $t > 0$ is a shifted parameter. It is obvious that, when $t$ is sufficient large, for any $x\neq 0$, we have $f_t(x) < 0$. From $f_t(0) = 0$, we know that $x = 0$ is the unique maximizer of $f_t(x)$. Denote the gradient of $f_t(x)$ as

Obviously, $x = 0$ is also a critical point of $f_t(x)$. The nonzero critical point of problem (2.6) is a Z-eigenvector corresponding to Z-eigenvalue $\lambda_t = (x_t^{\rm T}x_t)^{m-1}-t$. In this case, a suitable descent algorithm for solving shifted problem (2.6) should converge to a nonzero critical point. Therefore, we can get Z-eigenvalues and its associated Z-eigenvectors by solving problem (2.2) or (2.6). The algorithm is described as follows. □

Algorithm 1: Step 0: Given $\mathcal{A} \in \mathbb{S}^{[m, n]}, t\geq 1, \bar{\rho} > 1$.

Step 1: Solving problem (2.2) by using an algorithm to obtain $x_k$ and compute $\lambda_k = (x_k^{\rm T}x_k)^{m-1}$.

Step 2: If $\|x_k\| \leq \varepsilon$, stop, output $x_k$ and compute $\lambda_k = (x_k^{\rm T}x_k)^{m-1}-t$; otherwise, let $t: = \bar{\rho} t$, go to Step 3.

Step 3: Solve problem (2.6) by using an algorithm to obtain $x_k$, go to Step 2.

Remark 1. There is an inner loop between Step 2 and Step 3. Since descent algorithm for solving problem (2.6) will result in a nonzero critical point, then this inner loop can be terminated by finite iteration for sufficient large $t$. Therefore, Algorithm 1 is well defined.

Remark 2. When executing algorithm, it should use the same unconstrained optimization method to solve problems (2.2) and (2.6). We will propose a new accelerated conjugate gradient method, especially for solving problem (2.2) or (2.6).

3.   An accelerated conjugate gradient algorithm and its Convergence

In this section, we devote to giving an accelerated conjugate gradient method for solving unconstrained optimization problem, such as (2.2) or (2.6). Firstly, we consider the iterative formula of nonlinear conjugate gradient algorithm

The stepsize $\alpha_k > 0$ is determined by a line search and the direction $d_k$ are computed by

where $g_{k+1} = g(x_{k+1}), s_k = \alpha_kd_k$.

We first introduce a new conjugate parameter of our conjugate gradient method based on CD nonlinear conjugate gradient method (^[20]). The new conjugate parameter is

is introduced. When using the exact line search, $\beta_{k+1}$ in (3.3) reduces to CD conjugate gradient parameter. An accelerated parameter $\theta_{k+1}$ is obtained by the quasi-Newton direction (^[24,25]). Let $d_{k+1} = -G_{k+1}^{-1}g_{k+1}$, namely

where $G_{k+1}$ satisfies the following secant equation

From (3.4), then

Pre-multiplying at the both sides by $s_k^{\rm T}$, we have

Then, combined with (3.3)–(3.6), we have

The stepsize is generated by the strong Wolfe line search conditions

where $0 < \rho < \sigma < 1.$ We prove the descent property of the direction (3.2) under (3.8) and (3.9) in the following. Multiply both sides of (3.9) by $\alpha_k$, from $s_k = \alpha_kd_k$, we can easily obtain $(g_{k+1}^{\rm T}s_k)^2 \leq \sigma^2(-g_k^{\rm T}s_k)^2$.

Theorem 3. If $\theta_{k+1} \geq \frac{1}{2}+2\sigma^2$, then the direction determined by (3.2) satisfies the sufficient descent condition

Proof. Multiplying (3.2) by $g_{k+1}^{\rm T}$, we obtain

Using the inequality $a^{\rm T}b\leq \frac{1}{2}(\|a\|^2+\|b\|^2)$, where $a, b \in \mathbb{R}^n$, we have

Substituted (3.12) into (3.11), we have

Then from $(g_{k+1}^{\rm T}s_k)^2 \leq \sigma^2(-g_k^{\rm T}s_k)^2$ and $\theta_{k+1}\geq \frac{1}{2}+2\sigma^2$, we have

The proof is completed. □

Now, we describe an accelerated conjugate gradient algorithm (ACG).

ACG algorithm

Step 0: Given $x_0 \in \mathbb{R}^n, \varepsilon \geq 0 \; \mbox{and}\; 0 < \rho < \sigma < 1$. Compute $g_0$, let $d_0 = -g_0$. Set $k: = 0$.

Step 1: If $\|g_k\| \leq \varepsilon$, stop, output $x_k$. Otherwise, calculate $\alpha_k$ from (3.8) and (3.9). Let $x_{k+1} = x_k+\alpha_kd_k$ and $s_k = \alpha_kd_k$.

Step 2: Compute $g_{k+1}, y_k,$ $\beta_{k+1}$ by (3.3) and $\theta_{k+1}$ by (3.7). Let $\theta_{k+1}: = \max\{\theta_{k+1}, 2\sigma^2+\frac{1}{2}\}$.

Step 3: Using (3.2) to obtain $d_{k+1}$. Set $k: = k+1$ and go to step 1.

Now, we establish the convergence result of ACG algorithm. Let the level set $\Omega = \{x\in \mathbb{R}^n| f(x)\leq f(x_0)\}$ be a bounded closed set, i.e., there exists a constant $\gamma > 0$ such that $\|x\|\leq \gamma$ for all $x \in \Omega$. To facilitate analyzing, denote $\mathcal{A}_{i_1} = (a_{{i_1 i_2 \cdots i_m}})_{1\leq i_2, i_3, \cdots, i_m\leq n}$ and $\mathcal{A}_{i_1i_2} = (a_{{i_1 i_2 \cdots i_m}})_{1\leq i_3, i_4, \cdots, i_m\leq n}$.

Lemma 1. Consider tensors $\mathcal{A}\in \mathcal{S}^{[m, n]}$, then $\mathcal{A}x^m$ is Lipschitz continuous on $\Omega$.

Proof. Let $p(x) = \mathcal{A}x^m$, mathematical induction is adopted. If $m = 1$, $p(x) = \sum\limits_{i = 1}^{n}a_ix_i$, utilizing equivalence of vector norm, for all $x$, $y \in \Omega$, we have

Assume that the statements satisfy for all order $m\leq k-1$. When $m = k$, we have

Since $\Omega$ is a bounded close set, then $\|\mathcal{A}_{i_1}x^{k-1}\|$ is bounded on $\Omega$ and $\|z\|\leq P_1$. From $\sum\limits_{i_1 = 1}^{n}|x_{i_1}-y_{i_1}|\leq \|x-y\|$ and $\sum\limits_{i_1 = 1}^{n}|y_{i_1}| \leq \|z\|,$ there exists a positive constant $P_3$ such that

Namely, $\mathcal{A}x^m$ is Lipschitz continuous on $\Omega$. The proof is completed. □

Lemma 2. If $\mathcal{A}\in \mathcal{S}^{[m, n]}$, then $\mathcal{A}x^{m-1}$ and $\mathcal{A}x^{m-2}$ are Lipschitz continuous on $\Omega$.

Proof. For all $x$, $y \in \Omega$, using Lemma 3.1 and equivalence of norm, we have

where $P_4$ depends on tensor $\mathcal{A}$ and set $\Omega$.

Similarly, for all $x$, $y \in \Omega$, using Lemma 3.1 and equivalence of norm, we have

where $P_5$ depends on tensor $\mathcal{A}$ and set $\Omega$.

□

Lemma 3. If $\mathcal{A}\in \mathcal{S}^{[m, n]},$ then $g(x)$ is Lipschitz continuous in a neighbourhood $\mathbb{N}$ of $\Omega$, namely

holds for any $x, y \in \mathbb{N}$, where $L$ is a positive number.

Proof. There are two cases of gradient $g(x)$ to consider. One case is computed by (2.3) and the other case is computed by (2.7).

For (2.3): since $\mathbb{N}$ is a bound closed set, from Lemma 3.1, for all $x, y \in \mathbb{N}$, we have

For (2.7): since $\mathbb{N}$ is a bound closed set, from Lemma 3.1, for all $x, y \in \mathbb{N}$, we have

The proof is completed. □

We can easily get that there exists a constant $M > 0$ such that $\|g(x)\|\leq M.$ The following useful lemma was essentially which proved by Zoutendijk ^[26]. From Theorem 3.1, we can obtain that the sequence $\{d_k\}$ generated by ACG algorithm satisfies the following Lemmas.

Lemma 4. Let the sequences $\{x_k\}$ and $\{d_k\}$ be generated by ACG algorithm, we have

Lemma 5. Let the sequence $\{x_k\}$ be generated by ACG algorithm. If

then

From Theorem 3.1, (3.17) and Lemma 3.4, the result can be proved, which is omitted here.

Theorem 4. Let the sequence $\{x_k\}$ be generated by ACG algorithm. Then, we have

Proof. From Theorem 3.1, there exists a constant $c < 0$ satisfying $g_k^{\rm T}s_k\leq c\|g_k\|\|s_k\| < 0$, i.e., $-g_k^{\rm T}s_k\geq -c\|g_k\|\|s_k\|$. Then, we have

where $\xi = \frac{M(1+\sigma)}{-c}$. According to $|g_{k+1}^{\rm T}d_k| \leq -\sigma g_k^{\rm T}d_k$ and $\|s_k\|\leq 2\gamma$, we have

Because of $\theta_{k+1}\geq 2\sigma^2+\frac{1}{2}$, so $y_k^{\rm T}g_{k+1} > 0$. Without losing generality, let $y_k^{\rm T}g_{k+1} > \kappa > 0$, then we have

Therefore,

We have

From Lemma 3.5, it follows that (3.20) is derived. The proof is completed. □

Theorem 5. Let problems (2.2) and (2.6) are solved by ACG algorithm. ACG algorithm is well defined.

Proof. In ACG algorithm, we obtain the Z-eigenvalues of symmetric tensors by solving problem (2.2) or (2.6). When problem (2.6) is solved, it means that solving problem (2.2) results in a zero critical point. Then ACG algorithm turn to solve problem (2.6) which converges to a nonzero critical point. Moreover, since the convergence of our algorithm has been guaranteed, so the termination criteria condition always holds. That is, ACG algorithm is well defined. The proof is completed. □

4.   Numerical experiments

In this section, we report some numerical performance of ACG algorithm for solving problems (2.2) and (2.6). For convenience, we provide a table of abbreviations for the methods in Table 1.

We compare ACG with HS and PRP, which have been reported to be very efficient for unconstrained optimization. All experiments are done on a PC with CPU 2.40GHz and 2.00GB RAM using MATLAB R2013a. In the implementation of ACG algorithm, we set parameters $\varepsilon = 10^{-5}, \rho = 0.1, \sigma = 0.5, t = 1, \bar{\rho} = 2$. In Table 2, $Ex$ is the number of example, $n$ is the dimension, $k$ is the number of iterations, CPU stands for the time costed by algorithms (in seconds), $\lambda^*$ stands for Z-eigenvalue outputted by algorithms. All algorithms share the same start points and stopping criteria. In the following examples, the tensors $\mathcal{A}$ are originally from ^[27].

Example 1. Let $\mathcal{A} \in S^{[3, n]}$ defined by

Example 2. Let $\mathcal{A} \in S^{[3, n]}$ defined by

Example 3. Let $\mathcal{A} \in S^{[4, n]}$ defined by

Example 4. Let $\mathcal{A} \in S^{[4,3]}$ defined by

Example 5. Let $\mathcal{A} \in S^{[4, n]}$ defined by

Example 6. Let $\mathcal{A} \in S^{[4, n]}$ defined by

In Tables 2 and 3, we compare the numerical results of ACG algorithm and SPM algorithm (in ^[6]), PRP, HS and CD methods. Although the computed Z-eigenvectors $x^*$ associated to $\lambda^*$ are not shown, the values of $\|g_k\|$ are listed which are satisfy the given precision. It implies that $(\lambda^*, x^*)$ are considered as true solutions of problem (2.1). Two methods reach the same Z-eigenvalues for problems with same dimensions. ACG algorithm requires less iterations and CPU time than that of SPM algorithm. That is, we can see that ACG algorithm is competitive for computing Z-eigenvalues of symmetric tensors.

To show the numerical performance of a given optimal method, that the number of iterations ($k$) and CPU time (CPU) are important factors. So, we employ the profiles introduced by Dolan and Mor${\rm \acute{e}}$ ^[28] to analyze the efficiency of ACG, PRP, HS ^[20], QN ^[29] ATTCG and ADL^[30,31] methods, with the following conjugate gradient parameters, respectively,

Let $Y$ and $W$ be the set of methods and test problems, $n_y$, $n_w$ be the number of methods and test problems, respectively. The performance profile $\psi:\mathbb{R} \rightarrow [0, 1]$ is for each $y\in Y$ and $w\in W$ defined that $a_{w, y} > 0$ is $k$ or CPU required to solve problems $w$ by method $y$. Furthermore, the performance profile is obtained by

where $\tau > 0$, $size\{\cdot \}$ is the number of the elements in a set, and $r_{w, y}$ is the performance ratio defined as

In a performance profile plot, the top curve is a method that solved most problems in a time that is within a factor of best time. The horizontal axis gives the percentage ($\tau$) of the test problems for which a method is the fastest (efficiency), while the vertical side gives the percentage ($\psi$) of the test problems that are successfully solved by each of the methods. If program runs failure, or the number of iterations can reach more than 500, it is regarded as failed. And we denote the number of iterations by 500 and CPU time by 200 seconds. In this way, only ACG algorithm can solve all test problems.

As can be seen from Figure 1 shows the CPU time performance of the ACG algorithm and the other algorithms. It can be seen from the figure that when $\tau > 3$, the curves of the ACG algorithm and the ATTCG algorithm are similar, but when $\tau > 3.5$, both the ACG algorithm and ADL algorithm tend to be stable and coincide. Figure 2, the ACG algorithm is better than other algorithms in terms of the number of iterations, especially when $\tau > 2$, the curve of ACG algorithm becomes stable, which indicates that ACG algorithm can solve the problem only with fewer iterations. Therefore, Figures 1 and 2 show that the ACG algorithm proposed in this paper converge to the solution quickly.

5.   Conclusions

We constructed the unconstrained optimization problems with a shifted parameter. Based on the shifted unconstrained optimization problems, we presented an accelerated conjugate gradient method by using the quasi-Newton direction for solving them. Furthermore, we showed the global convergence analysis of the proposed algorithm. Numerical experiments demonstrated that our method has good numerical performance. We further highlight that the proposed algorithm can be used in other fields, such as the symmetric system of nonlinear equations. It is vital to note that some new methods with random technology will be taken into account in our future work.

Acknowledgments

National Natural Science Foundation of China (12061087), the Scientific and Technological Developing Scheme of Jilin Province (YDZJ202101ZYTS167, YDZJ202201ZYTS303), the Yunnan Provincial Ten Thousands Plan Young Top Talents, the Project of education department of Jilin Province (JJKH20210030KJ).

Conflict of interest

The authors declare no conflicts of interest.