A QR Algorithm for Symmetric Tensors

In this paper we focus on symmetric tensors. In particular, a $d$th-order symmetric tensor has all its dimensions being equal, namely, $n_{1}=\cdots=n_{d}=n$, and $\mathcal{A}_{i_{1}\ldots i_{d}}=\mathcal{A}_{\pi(i_{1}\ldots i_{d})}$, where $\pi(i_{1}\ldots i_{d})$ is any permutation of the indices $i_{1}\ldots i_{d}$. To facilitate notation we further introduce some shorthand specific to symmetric tensors. First, we use $\mathbb{S}^{[d,n]}$ to denote the set of all real symmetric $d$th-order $n$-dimensional tensors. When a row vector $x^{T}\in\mathbb{R}^{1\times n}$ is to be multiplied onto $k$ modes of $\mathcal{A}\in\mathbb{S}^{[d,n]}$ (due to symmetry it can be assumed without loss of generality that all but the first $d-k$ modes are multiplied with $x^{T}$), we use the shorthand

(5a)		$\displaystyle\mathcal{A}x$	$\displaystyle\triangleq\mathcal{A}{{}_{\times_{d}}}{x^{T}}\in\mathbb{S}^{[d-1,n]},$
(5b)		$\displaystyle\mathcal{A}x^{2}$	$\displaystyle\triangleq\mathcal{A}{{}_{\times_{d-1}}}{x^{T}}{{}_{\times_{d}}}{x^{T}}\in\mathbb{S}^{[d-2,n]},\cdots$
(5c)		$\displaystyle\mathcal{A}x^{d-2}$	$\displaystyle\triangleq\mathcal{A}{{}_{\times_{3}}}{x^{T}}{{}_{\times_{4}}}{x^{T}}\cdots{{}_{\times_{d}}}x^{T}\in\mathbb{S}^{[2,n]}\in\mathbb{R}^{n\times n},$
(5d)		$\displaystyle\mathcal{A}x^{d-1}$	$\displaystyle\triangleq\mathcal{A}{{}_{\times_{2}}}{x^{T}}{{}_{\times_{3}}}{x^{T}}\cdots{{}_{\times_{d}}}x^{T}\in\mathbb{S}^{[1,n]}=\mathbb{R}^{n},$
(5e)		$\displaystyle\mathcal{A}x^{d}$	$\displaystyle\triangleq\mathcal{A}{{}_{\times_{1}}}{x^{T}}{{}_{\times_{2}}}{x^{T}}\cdots{{}_{\times_{d}}}x^{T}\in\mathbb{S}^{[0,n]}=\mathbb{R}.$

This is a contraction process because the resulting tensors are progressively shrunk wherein we assume any higher order singleton dimensions ($n_{i}=1$) are “squeezed”, e.g., $\mathcal{A}x^{d-1}$ is treated as an $n\times n$ matrix rather than its equivalent $n\times n\times 1\cdots\times 1$ tensor. Besides contraction by a row vector, the notation generalizes to any matrix of appropriate dimensions. Suppose we have a matrix $V\in\mathbb{R}^{n\times p}(p>1)$, then

(6)

$\displaystyle\mathcal{A}V^{d-1}\triangleq\mathcal{A}{{}_{\times 2}}{V^{T}}{{}_{\times 3}}{V^{T}}\cdots{{}_{\times d}}{V^{T}}\in\mathbb{R}^{n\times p\cdots\times p}.$

Now, the real eigenpair of a symmetric tensor $\mathcal{A}\in\mathbb{S}^{[d,n]}$, namely, $(\lambda,x)\in\mathbb{R}\times\mathbb{R}^{n}$, is defined as

(7)

$\displaystyle\mathcal{A}x^{d-1}=\lambda x,\mbox{~~}x^{T}x=1.$

This definition is called the (real) Z-eigenpair [16] or $l^{2}$-eigenpair [14] which is of interest in this paper, apart from other differently defined eigenpairs, e.g., [11]. It follows that

(8)

$\displaystyle x^{T}\mathcal{A}x^{d-1}=\mathcal{A}x^{d}=\lambda.$

For $\mathcal{A}\in\mathbb{S}^{[d,n]}$, depending on whether $d$ is odd or even, we have the following result regarding eigenpairs that occur in pairs [10]:

Subsequently, we do not call two eigenpairs distinct if they follow the relationship in Lemma 1. For even-order symmetric tensors, there exists a symmetric identity tensor counterpart $\mathcal{E}\in\mathbb{S}^{[d,n]}$ such that [16, 10]

(9)

$\displaystyle\mathcal{E}x^{d-1}=x~\mbox{~for all~}~x^{T}x=1.$

The constraint that $||x||=1$ is necessary so it holds for all even $d\geq 2$. Obviously, $(1,x)$ is an eigenpair for $\mathcal{E}$. It is also obvious that such an identity tensor does not exist for an odd $d$ since contracting with $-x$ does not lead to itself. Consequently, for an even $d$, we have similar shifting of eigenvalues when we add a multiple of $\mathcal{E}$ to $\mathcal{A}$, i.e., if $(\lambda,x)$ is an eigenpair of $\mathcal{A}$, then $(\lambda+\alpha,x)$ is an eigenpair of $\mathcal{A}+\alpha\mathcal{E}$.

In this paper, we say that two symmetric tensors $\mathcal{A},\mathcal{B}\in\mathbb{S}^{[d,n]}$ are similar, denoted by $\mathcal{A}\sim\mathcal{B}$, if for some non-singular $P\in\mathbb{R}^{n\times n}$

(10)

$\displaystyle\mathcal{B}=\mathcal{A}P^{d},$

so that the inverse transform of $\mathcal{B}$ by $P^{-1}$ will return $\mathcal{A}$, namely,

(11)

$\displaystyle\mathcal{B}P^{-d}=\left(\mathcal{A}P^{d}\right)P^{-d}=\mathcal{A}{{}_{\times 1}}(P^{-T}P^{T}){{}_{\times 2}}(P^{-T}P^{T})\cdots{{}_{\times d}}(P^{-T}P^{T})=\mathcal{A}.$

Suppose $\mathcal{\tilde{A}}\in\mathbb{S}^{[d,n]}$ is similar to $\mathcal{A}$ through a similarity transform by an orthogonal matrix $\underline{Q}$, i.e., $\underline{Q}^{T}\underline{Q}=I_{n}$ and $\mathcal{\tilde{A}}=\mathcal{A}\underline{Q}^{d}$. It follows that if $(\lambda,e_{1})$ is an eigenpair of $\mathcal{\tilde{A}}$, then $(\lambda,\underline{q}_{1})$ is an eigenpair of $\mathcal{A}$ where $\underline{q}_{1}$ denotes the first column of $\underline{Q}$. This is seen by

	$\displaystyle\mathcal{\tilde{A}}e_{1}^{d-1}=\lambda e_{1}\Leftrightarrow$	$\displaystyle\mathcal{A}{{}_{\times 1}}\underline{Q}^{T}{{}_{\times 2}}(e_{1}^{T}\underline{Q}^{T})\cdots{{}_{\times d}}(e_{1}^{T}\underline{Q}^{T})=\lambda e_{1}$
	$\displaystyle\Leftrightarrow$	$\displaystyle\mathcal{A}{{}_{\times 1}}\underline{Q}^{T}{{}_{\times 2}}(\underline{q}_{1})^{T}\cdots{{}_{\times d}}(\underline{q}_{1})^{T}=\lambda e_{1}$
	$\displaystyle\Leftrightarrow$	$\displaystyle\underline{Q}^{T}\mathcal{A}\underline{q}_{1}^{d-1}=\lambda e_{1}$
(12)		$\displaystyle\Leftrightarrow$	$\displaystyle\mathcal{A}\underline{q}_{1}^{d-1}=\lambda\underline{q}_{1}.$

In fact, we do not need to limit ourselves to $e_{1}$ and the above result easily generalizes to any $e_{i}$, which we present as lemma 2.

The symmetric tensor $\mathcal{A}\in\mathbb{S}^{[3,3]}$ in Figure 2(a) will be used as a running example. We call it a labeling matrix as it labels each distinct entry in a symmetric tensor sequentially with a natural number. Moreover, Figure 2(b) is the matricization of this labeling matrix along the first mode [9]. As a matter of fact, matricization about any mode of $\mathcal{A}$ gives rise to the same matrix due to symmetry.

Now we formally propose the QRST algorithm, summarized in Algorithm 3.1. The core of QRST is a two-step iteration, remarked in Algorithm 3.1 as Steps 1 and 2, for the $i$th slice where $i=1,\cdots,n$. In the Algorithm, it should be noted that Matlab convention is used whereby it can be recognized that $\mathcal{A}_{k}(:,i,i,\cdots,i)=\mathcal{A}_{k}e_{i}^{d-1}$ is a column and $\mathcal{A}_{k}(:,:,i,\cdots,i)=\mathcal{A}_{k}e_{i}^{d-2}$ is a (symmetric) square matrix.

Algorithm 3.1.

(Shifted) QRST
Input: $\mathcal{A}_{0}\in\mathbb{S}^{[d,n]}$, tolerance $\tau$, maximum iteration $k_{max}$, shifts $s_{k}$ from shift scheme
Output: eigenpairs $(\lambda_{i},\underline{q}_{i})$  for $i=1:n$ do # \eqparboxCOMMENTiterate through the $i$th square slice $\mathcal{A}_{0}(:,:,i,\cdots,i)$   $\epsilon\leftarrow||\mathcal{A}_{0}e_{i}^{d-1}-e_{i}||_{2}/||\mathcal{A}_{0}(:,:,i,\cdots,i)||_{2}$   $\underline{Q}\leftarrow I_{n}$   iteration count $k\leftarrow 1$   while $\epsilon>\tau$ or $k\leq k_{max}$ do    Step 1: $Q_{k}R_{k}=\mathcal{A}_{k-1}(:,:,i,\cdots,i)+s_{k-1}I_{n}$   # perform QR factorization    Step 2: $\mathcal{A}_{k}=\mathcal{A}_{k-1}Q_{k}^{d}$   # similarity transform    $\underline{Q}\leftarrow\underline{Q}Q_{k}$    $\epsilon\leftarrow||\mathcal{A}_{k}e_{i}^{d-1}-e_{i}||_{2}/||\mathcal{A}_{k}(:,:,i,\cdots,i)||_{2}$   end while  end for  if converged then   $\lambda_{i}\leftarrow\mathcal{A}_{k}(i,i,i,\cdots,i)$   $\underline{q}_{i}\leftarrow\underline{Q}(:,i)$  end if

Referring to Figure 2(b), if we take $\mathcal{A}_{0}=\mathcal{A}$, $i=1$ and start with $\mathcal{A}_{0}e_{1}^{d-2}=\mathcal{A}_{0}e_{1}$, the first (leftmost) $3\times 3$ “slice” of the $3\times 9$ matrix is extracted for QR factorization to get $Q_{1}R_{1}$, and then all three modes of $\mathcal{A}_{0}$ are multiplied with $Q_{1}^{T}$ so the tensor symmetry is preserved. Using Algorithm 3.1 with the heuristic shift described in Section 3.2, this tensor QR algorithm converges in $152$ iterations when $\mathcal{A}_{k}e_{1}^{2}$ becomes a multiple of $e_{1}$ (i.e., all but the first entry in the $\mathcal{A}_{k}e_{1}^{2}$ column are numerically zero), resulting in

$\displaystyle\left[{\begin{array}[]{*{20}c}-0.14&0.00&0.00&\vline&0.00&-0.19&0.56&\vline&0.00&0.56&-0.40\\ 0.00&-0.19&0.56&\vline&-0.19&-20.69&11.56&\vline&0.56&11.56&-5.56\\ 0.00&0.56&-0.40&\vline&0.56&11.56&-5.56&\vline&-0.40&-5.56&2.46\\ \end{array}}\right],$

for which $e_{1}=[1~0~0]^{T}$ is obviously an eigenvector as its contraction along the $i_{2}$ and $i_{3}$ axes will extract the first column $[-0.14~0~0]^{T}$. Using (12), this corresponds to the eigenpair $(-0.14,[-0.7854~0.6029~-0.1401]^{T})$ of the labeling tensor. In fact, for different $i$’s, we shoot for different converged patterns in the corresponding square slices as depicted in Figure 3, where it is obvious that $\mathcal{A}_{k}e_{i}^{2}=\lambda_{i}e_{i}$ for each $i$.

Since the core computation of QRST lies in the QR factorization ($\mathcal{O}(n^{3})$) and the corresponding similarity transform ($\mathcal{O}(dn^{2}n^{d-1})=\mathcal{O}(dn^{d+1})$), the overall complexity is therefore dominated by the similarity transform process. Note that this high complexity comes about whereby the tensor symmetry is not exploited. When a low-rank decomposition for $\mathcal{A}_{0}$ exists that represents $\mathcal{A}_{0}$ as a finite sum of $r$ rank-$1$ outer factors, the complexity of QRST can be largely reduced to $\mathcal{O}(drn^{2})$. However, we refrain from elaborating as it is beyond our scope of verifying the idea of QRST. Such low-rank implementation will be presented separately in a sequel of this work.

We assume an initial symmetric $\mathcal{A}_{0}=\mathcal{A}\in\mathbb{S}^{[d,n]}$. The vector $e_{i}$ denotes the $i$th column of the identity matrix $I_{n}$ and we define $\underline{Q}_{k}\triangleq Q_{1}Q_{2}\cdots Q_{k}$ with $\underline{Q}_{0}=I_{n}$ and $\underline{Q}_{1}=Q_{1}$ etc. A tensor with higher order singleton dimensions is assumed squeezed, e.g., $\mathcal{A}_{0}e_{i}^{d-2}$ is regarded as an $n\times n$ matrix rather than its equivalent $n\times n\times 1\cdots\times 1$ tensor. Also, we use $\underline{q}_{i,k}$ to denote the $i$th column in $\underline{Q}_{k}$ and $r_{ij,k}$ to stand for the $ij$-entry of $R_{k}$. It will be shown that the unshifted QRST operating on the matrix slice resulting from contracting all but the first two modes $\mathcal{A}_{k}$ by $e_{i}$, i.e., $\mathcal{A}_{k}e_{i}^{d-2}$, is in fact an orthogonal iteration. This can be better visualized by enumerating a few iterations of QRST:

(13a)		$\displaystyle Q_{1}R_{1}$	$\displaystyle=\mathcal{A}_{0}e_{i}^{d-2},$
(13b)		$\displaystyle\mathcal{A}_{1}$	$\displaystyle=\mathcal{A}_{0}Q_{1}^{d}=\mathcal{A}_{0}\underline{Q}_{1}^{d},$
(13c)		$\displaystyle Q_{2}R_{2}$	$\displaystyle=\mathcal{A}_{1}e_{i}^{d-2}=\mathcal{A}_{0}{{}_{\times_{1}}}{\underline{Q}_{1}^{T}}{{}_{\times_{2}}}{\underline{Q}_{1}^{T}}{{}_{\times_{3}}}{\underline{q}_{i,1}^{T}}\cdots{{}_{\times_{d}}}{\underline{q}_{i,1}^{T}}=\underline{Q}_{1}^{T}(\mathcal{A}_{0}\underline{q}_{i,1}^{d-2})\underline{Q}_{1},$
(13d)		$\displaystyle\mathcal{A}_{2}$	$\displaystyle=\mathcal{A}_{1}Q_{2}^{d}=(\mathcal{A}_{0}\underline{Q}_{1}^{d})Q_{2}^{d}=\mathcal{A}_{0}(\underline{Q}_{1}Q_{2})^{d}=\mathcal{A}_{0}\underline{Q}_{2}^{d},$
(13e)		$\displaystyle Q_{3}R_{3}$	$\displaystyle=\mathcal{A}_{2}e_{i}^{d-2}=\mathcal{A}_{0}{{}_{\times_{1}}}{\underline{Q}_{2}^{T}}{{}_{\times_{2}}}{\underline{Q}_{2}^{T}}{{}_{\times_{3}}}{\underline{q}_{i,2}^{T}}\cdots{{}_{\times_{d}}}{\underline{q}_{i,2}^{T}}=\underline{Q}_{2}^{T}(\mathcal{A}_{0}\underline{q}_{i,2}^{d-2})\underline{Q}_{2},$

etc. From (13), we see that in general for $k=1,2,\cdots$,

(14a)		$\displaystyle\mathcal{A}_{k-1}$	$\displaystyle=\mathcal{A}_{0}\underline{Q}_{k-1}^{d},$
(14b)		$\displaystyle Q_{k}R_{k}$	$\displaystyle={\mathcal{A}_{k-1}}{e_{i}^{d-2}}=\underline{Q}_{k-1}^{T}(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})\underline{Q}_{k-1}.$

Multiplying $\underline{Q}_{k-1}$ to both sides of (14b) we get

	$\displaystyle\underline{Q}_{k}R_{k}=(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})\underline{Q}_{k-1}\Leftrightarrow$
(22)		$\displaystyle\left[{\begin{array}[]{*{20}c}|&|&{}\hfil&|\\ {\underline{q}_{1,k}}&{\underline{q}_{2,k}}&\cdots&{\underline{q}_{n,k}}\\ |&|&{}\hfil&|\\ \end{array}}\right]\left[{\begin{array}[]{*{20}c}{r_{11,k}}&{r_{12,k}}&\cdots&{r_{1n,k}}\\ {}\hfil&{r_{22,k}}&{}\hfil&{r_{2n,k}}\\ {}\hfil&{}\hfil&\ddots&\vdots\\ {}\hfil&{}\hfil&{}\hfil&{r_{nn,k}}\\ \end{array}}\right]$
(29)		$\displaystyle~~~~~~~~~~~~~~=\left[{\begin{array}[]{*{20}c}{}\hfil&{}\hfil&{}\hfil\\ {}\hfil&{\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2}}&{}\hfil\\ {}\hfil&{}\hfil&{}\hfil\\ \end{array}}\right]\left[{\begin{array}[]{*{20}c}|&|&{}\hfil&|\\ {\underline{q}_{1,k-1}}&{\underline{q}_{2,k-1}}&\cdots&{\underline{q}_{n,k-1}}\\ |&|&{}\hfil&|\\ \end{array}}\right].$

It should now become obvious that (29) is the tensor generalization of the matrix orthogonal iteration [6] (also called simultaneous iteration or subspace iteration) for finding dominant invariant eigenspace of matrices¹¹1Unlike the matrix case with a constant (static) matrix $A$, the tensor counterpart comes with a dynamic $\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2}$ that is updated in every pass with the $i$th column of $\underline{Q}_{k-1}$.. This result comes at no surprise as the matrix QR algorithm is known to be equivalent to orthogonal iteration applied to a block column vector. However, our work confirms, for the first time, this is also true also for symmetric tensors.

Now, if we set $i=1$ and extract the first columns on both sides of (29), we get

(30)

$\displaystyle\underline{q}_{1,k}r_{11,k}=(\mathcal{A}_{0}\underline{q}_{1,k-1}^{d-2})\underline{q}_{1,k-1}=\mathcal{A}_{0}\underline{q}_{1,k}^{d-1}\mbox{~~or~~}\underline{q}_{1,k}=\frac{\mathcal{A}_{0}\underline{q}_{1,k}^{d-1}}{||\mathcal{A}_{0}\underline{q}_{1,k}^{d-1}||},$

which is just the unshifted S-HOPM method [3, 8] cited as Algorithm 1 in [10]. If we start QRST with $i=2$, then the first two columns of (29) proceed as

(31a)		$\displaystyle\underline{q}_{1,k}r_{11,k}$	$\displaystyle=(\mathcal{A}_{0}\underline{q}_{2,k-1}^{d-2})\underline{q}_{1,k-1},$
(31b)		$\displaystyle\underline{q}_{1,k}r_{12,k}+\underline{q}_{2,k}r_{22,k}$	$\displaystyle=\mathcal{A}_{0}\underline{q}_{2,k-1}^{d-1}.$

Assuming convergence at $k=\infty$,

(32a)		$\displaystyle\underline{q}_{1,\infty}r_{11,\infty}$	$\displaystyle=(\mathcal{A}_{0}\underline{q}_{2,\infty}^{d-2})\underline{q}_{1,\infty},$
(32b)		$\displaystyle\underline{q}_{1,\infty}r_{12,\infty}+\underline{q}_{2,\infty}r_{22,\infty}$	$\displaystyle=(\mathcal{A}_{0}\underline{q}_{2,\infty}^{d-2})\underline{q}_{2,\infty}=\mathcal{A}_{0}\underline{q}_{2,\infty}^{d-1},$

from which we can easily check that $r_{12,\infty}=\underline{q}_{1,\infty}^{T}(\mathcal{A}_{0}\underline{q}_{2,\infty}^{d-2})\underline{q}_{2,\infty}=r_{11,\infty}\underline{q}_{1,\infty}^{T}\underline{q}_{2,\infty}=0$ due to (32a). In fact, for any $i$, if $\underline{Q}_{k}$ and $R_{k}$ converge to $\underline{Q}_{\infty}$ and $R_{\infty}$, respectively, then from (29), we must have

(37)

$\displaystyle R_{\infty}=\left[{\begin{array}[]{*{20}c}{r_{11,\infty}}&0&\cdots&0\\ 0&{r_{22,\infty}}&{}\hfil&\vdots\\ {\vdots}&{}\hfil&\ddots&0\\ 0&{\cdots}&{0}&{r_{nn,\infty}}\\ \end{array}}\right]=\underline{Q}_{\infty}^{T}(\mathcal{A}_{0}\underline{q}_{i,\infty}^{d-2})\underline{Q}_{\infty},$

where $R_{\infty}$ is diagonal due to symmetry on the right of equality. Now if we pre-multiply $\underline{Q}_{\infty}$ onto both sides of (37) and post-multiply with $e_{j}$, then

(38)

$\displaystyle\underline{q}_{j,\infty}r_{jj,\infty}=(\mathcal{A}_{0}\underline{q}_{i,\infty}^{d-2})\underline{q}_{j,\infty}$

which boils down to the eigenpair $(r_{ii,\infty},\underline{q}_{i,\infty})$ when $i=j$. So in effect, (shifted) QRST finds more than just the conventional eigenpairs, but actually all “local” eigenpairs of the square matrix $(\mathcal{A}_{0}\underline{q}_{i,\infty}^{d-2})$, one of them happens to be the tensor eigenpair. This explains why in practice, (shifted) QRST can locate even unstable tensor eigenpairs as it is not restricted to the “global” convexity/concavity constraint imposed on the function $f(x)=\mathcal{A}_{0}x^{d}$ that always leads to positive/negative stable eigenpairs [10].

Till now, it can be seen that except for $i=1$, (shifted) QRST differs from the (S)S-HOPM algorithm presented in [10, 11], and convergence is generally NOT ensured for $i\neq 1$. However, inspired by (32b), the key for convergence starts with the convergence of the first column that reads

(39)

$\displaystyle\underline{q}_{1,k}$

$\displaystyle=\frac{(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})\underline{q}_{1,k-1}}{||(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})\underline{q}_{1,k-1}||}.$

Consequently, applying results in [10, 11], convergence of (39) can be enforced locally around $\underline{q}_{1,k-1}$ provided the scalar function $f(x)=x^{T}(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})x$ is convex (concave), which is “localized” in the sense that the contracted $(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})$ is used instead of the whole $\mathcal{A}_{0}$. Then, for example, convexity can be achieved by adding a shift $s_{k-1}$ that makes $(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})+s_{k-1}I_{n}$ positive semidefinite, namely,

	$\displaystyle s_{k-1}$	$\displaystyle\geq-\lambda_{\min}\left(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2}\right)$
		$\displaystyle=-\lambda_{\min}\left(\underline{Q}_{k-1}^{T}(\mathcal{A}_{0}\underline{q}_{i,k-1}^{d-2})\underline{Q}_{k-1}\right)$
(40)			$\displaystyle=-\lambda_{\min}\left(\mathcal{A}_{k-1}e_{i}^{d-2}\right).$

So the shift step in Algorithm 3.1 is simply $s_{k-1}=-\lambda_{\min}(\mathcal{A}_{k-1}(:,:,i,\cdots,i))+\delta$ where $\delta$ is a small positive quantity which is set to be unity in our tests unless otherwise specified. This is analogous to the core step in SS-HOPM (Algorithm 2 of [10]). In short, the shift strategy in [10, 11] can be locally reused for QRST with the advantage that QRST may further compute the unstable eigenpairs alongside the positive/negative-stable ones, since a group of eigenpairs are generated at a time due to analogy of QRST to orthogonal iteration.

For a symmetric $\mathcal{A}_{0}\in\mathbb{S}^{[d,n]}$, QRST produces at most $n$ eigenpairs in one pass, and fewer than $n$ eigenpairs can result when QRST does not converge or when the same eigenpair is computed repeatedly. To increase the probability of generating distinct eigenpairs, and noting that eigenvectors of similar tensors are just a change of basis, we introduce the heuristic of permutation to scramble the entries in $\mathcal{A}_{0}$ as much as possible. Figure 4 shows how the $\mathbb{S}^{[3,3]}$ labeling tensor is permuted by the $n!$ $(n=3)$ (orthogonal) permutation matrices.

As such, we can run QRST on a tensor subject to different permutations. We call this permuted QRST or PQRST, summarized in Algorithm 3.2.

Algorithm 3.2.

PQRST
Input: $\mathcal{A}_{0}\in\mathbb{S}^{[d,n]}$, preconditioning orthogonal matrices in a cell $P\{p\}$, $p=1,\cdots,\bar{p}$
Output: eigenpairs $(\lambda_{i},\underline{q}_{i})$  for $p=1:\bar{p}$ do   $P\leftarrow P\{p\}$   $\mathcal{A}_{p}\leftarrow\mathcal{A}_{0}P^{d}$   call QRST Algorithm 3.1 with $\mathcal{A}_{p}$ and collect distinct converged eigenpairs  end for

Of course, PQRST requires running QRST $n!$ times which quickly becomes impractical for large $n$. In that case, we should limit the number of permutations by sampling across all possible permutations or use other heuristics for setting the initial pre-transformed $\mathcal{A}_{0}$, which is our on-going research.

As a first example we use the labeling tensor of Figure 2(a). Due to the odd order of the tensor we have that both $(\lambda,x)$ and $(-\lambda,-x)$ are eigenpairs. From solving the polynomial system we learn that there are 4 distinct eigenpairs, listed in Table 1. The stability of each of these eigenpairs was determined by checking the signs of the eigenvalues of the projected Hessian of the Lagrangian as described in [10].

100 runs of the SS-HOPM algorithm were performed, each with a different random initial vector. Choosing the random initial vectors from either a uniform or normal distribution had no effect whatsoever on the final result. The shift was fixed to the “conservative” choice of $(m-1)\sum_{i_{1},\ldots,i_{m}}|a_{i_{1}\cdots i_{m}}|=288$. All 100 runs converged to the eigenpair with the eigenvalue $30.4557$. The median of number of iterations is 169. One run of 169 iterations took about 0.02 seconds. Next, 100 runs of the adaptive SS-HOPM algorithm were performed with a tolerance of $10^{-15}$. When the initial guess for $x$ was a random vector sampled from a uniform distribution, then all 100 runs consistently converged to the $30.4557$ eigenpair with a median number of 14 iterations. This clearly demonstrates the power of the adaptive shift in order to reduce the number of iterations. Choosing the initial vector $x$ from a normal distribution resulted in finding the 3 stable eigenpairs and 2 iterations not converging (Table 2). The error was computed as the average value of $||\mathcal{A}x^{m-1}-\lambda\,x||_{2}$. Similar to the SS-HOPM algorithm, the PQRST algorithm without shifting can only find the “largest” eigenpair with a median of 13 iterations. This took about 0.082 seconds. PQRST with shift $\delta=1$ results in 18 runs that find all eigenpairs (Table 3). From the 18 runs, 8 did not converge.

For the second example we consider the 4th-order tensor from Example 3.5 in [10]. Solving the polynomial system shows that there are 11 distinct eigenpairs (Table 4).

Both a positive and negative shift are used for the (adaptive) SS-HOPM runs. Tables 5 and 6 list the results of applying 200 (100 positive shifts and 100 negative shifts) SS-HOPM and adaptive SS-HOPM runs to this tensor. Using the conservative shift of $37$ and the same convergence criterion as in Example 4.1 resulted in 27 runs not converging for the SS-HOPM algorithm. Again, the main strength of the adaptive SS-HOPM is the reduction of number of iterations required to converge. A faster convergence also implies a smaller error as observed in the numerical results. The same initial random vectors were used for both the SS-HOPM and adaptive SS-HOPM runs. Observe that in the 100 adaptive SS-HOPM runs with a negative shift the eigenpair with eigenvalue $-0.0451$ was never found. The total run time for 100 iterations of one adaptive SS-HOPM algorithm took about 0.025 seconds, which means that 100 runs take about 2.5 seconds to complete.

Similarly, the PQRST algorithm without shifting does not converge in any of the permutations. PQRST however retrieves all eigenpairs except one (Table 7). The total time to run PQRST for all permutations for 100 iterations took about 0.35 seconds.