# Low Rank Tensor Decompositions and Approximations

Jiawang Nie, Li Wang, and Zequn Zheng  Jiawang Nie and Zequn Zheng Department
of Mathematics, University of California San Diego, 9500 Gilman Drive, La
Jolla, CA, USA, 92093.  Li Wang Department of Mathematics, University of Texas
at Arlington, 411 South Nedderman Drive, Arlington, TX, 76019.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

There exist linear relations among tensor entries of low rank tensors. These
linear relations can be expressed by multi-linear polynomials, which are
called generating polynomials. We use generating polynomials to compute tensor
rank decompositions and low rank tensor approximations. We prove that this
gives a quasi-optimal low rank tensor approximation if the given tensor is
sufficiently close to a low rank one.

###### Key words and phrases:

tensor, decomposition, rank, approximation, generating polynomial

###### 2010 Mathematics Subject Classification:

15A69,65F99

##  1\. Introduction

Let $m$ and $n_{1},\ldots,n_{m}$ be positive integers. A tensor $\mathcal{F}$
of order $m$ and dimension $n_{1}\times\cdots\times n_{m}$ can be labelled
such that

$\mathcal{F}=(\mathcal{F}_{i_{1},\ldots,i_{m}})_{1\leq i_{1}\leq
n_{1},\ldots,1\leq i_{m}\leq n_{m}}.$

Let $\mathbb{F}$ be a field (either the real field $\mathbb{R}$ or the complex
field $\mathbb{C}$). The space of all tensors of order $m$ and dimension
$n_{1}\times\cdots\times n_{m}$ with entries in $\mathbb{F}$, is denoted as
$\mathbb{F}^{n_{1}\times\cdots\times n_{m}}$. For vectors
$v_{i}\in\mathbb{F}^{n_{i}}$, $i=1,\ldots,m$, their outer product
$v_{1}\otimes\ldots\otimes v_{m}$ is the tensor in
$\mathbb{F}^{n_{1}\times\cdots\times n_{m}}$ such that

$(v_{1}\otimes\cdots\otimes
v_{m})_{i_{1},\ldots,i_{m}}\,=\,(v_{1})_{i_{1}}\cdots(v_{m})_{i_{m}}$

for all labels in the corresponding range. A tensor like
$v_{1}\otimes\ldots\otimes v_{m}$ is called a rank-$1$ tensor. For every
tensor $\mathcal{F}\in\mathbb{F}^{n_{1}\times\cdots\times n_{m}}$, there exist
vector tuples $(v^{s,1},\ldots,v^{s,m})$, $s=1,\ldots,r$,
$v^{s,j}\in\mathbb{F}^{n_{j}}$, such that

(1.1) $\displaystyle\mathcal{F}=\sum_{s=1}^{r}v^{s,1}\otimes\cdots\otimes
v^{s,m}.$

The smallest such $r$ is called the $\mathbb{F}$-rank of $\mathcal{F}$, for
which we denote $\text{rank}_{\mathbb{F}}(\mathcal{F})$. When $r$ is minimum,
the equation (1.1) is called a rank-$r$ decomposition over the field
$\mathbb{F}$. In the literature, this rank is sometimes referenced as the
candecomp-parafac (CP) rank. We refer to [22, 24] for various notions of
tensor ranks. Recent work for tensor decompositions can be found in [1, 9, 10,
13, 26, 36, 38]. Tensors are closely related to polynomial optimization [7,
14, 29, 31, 32]. Various applications of tensors can be found in [18, 20, 30].
Throughout the paper, we use the Hilbert-Schmidt norm for tensors:

$\|\mathcal{F}\|=\sqrt{\sum_{\begin{subarray}{c}1\leq i_{j}\leq n_{j},1\leq
j\leq m\end{subarray}}|\mathcal{F}_{i_{1},\ldots,i_{m}}|^{2}}.$

The low rank tensor approximation (LRTA) problem is to find a low rank tensor
that is close to a given one. This is equivalent to a nonlinear least square
optimization problem. For a given tensor
$\mathcal{F}\in\mathbb{F}^{n_{1}\times\cdots\times n_{m}}$ and a given rank
$r$, we look for $r$ vector tuples $v^{(s)}:=(v^{s,1},\ldots,v^{s,m})$,
$s=1,\ldots,r$, such that

$\mathcal{F}\,\approx\,\sum\limits_{s=1}^{r}v^{s,1}\otimes\cdots\otimes
v^{s,m},\quad v^{s,j}\in\mathbb{F}^{n_{j}}.$

This requires to solve the following nonlinear least squares optimization

(1.2)
$\displaystyle\min\limits_{v^{(s,j)}\in\mathbb{F}^{n_{j}},\,j=1,\ldots,m}\big{\|}\mathcal{F}-\sum\limits_{s=1}^{r}v^{s,1}\otimes\cdots\otimes
v^{s,m}\big{\|}^{2}.$

When $r=1$, the best rank-$1$ approximating tensor always exists and it is
equivalent to computing the spectral norm [16, 19]. When $r>1$, the best
rank-$r$ tensor approximation may not exist [11]. Classically used methods for
solving low rank tensor approximations are the alternating least squares (ALS)
method [5, 25, 40], higher order power iterations [8], semidefinite
relaxations [7, 29], SVD based methods [17], optimization based methods [15,
36]. We refer to [1, 27] for recent work on low rank tensor approximations.

### Contributions

In this paper, we extend the generating polynomial method in [26, 27] to
compute tensor rank decompositions and low rank tensor approximations for
nonsymmetric tensors. First, we estimate generating polynomials by solving
linear least squares. Second, we find their approximately common zeros, which
can be done by computing eigenvalue decompositions. Third, we get a tensor
decomposition from their common zeros, by solving linear least squares. To
find a low rank tensor approximation, we first apply the decomposition method
to obtain a low rank approximating tensor and then use nonlinear optimization
methods to improve the approximation. Our major conclusion is that if the
tensor to be approximated is sufficiently close to a low rank one, then the
obtained low rank tensor is a quasi-optimal low rank approximation. The proof
is based on perturbation analysis of linear least squares and eigenvalue
decompositions.

The paper is organized as follows. In Section 2, we review some basic results
about tensors. In Section 3, we introduce the concept of generating
polynomials and study their relations to tensor decompositions. In Section 4,
we give an algorithm for computing tensor rank decompositions for low rank
tensors. In Section 5, we give an algorithm for computing low rank
approximations. The approximation error analysis is also given. In Section 6,
we present numerical experiments. Some conclusions are made in Section 7.

##  2\. preliminary

### Notation

The symbol $\mathbb{N}$ (resp., $\mathbb{R}$, $\mathbb{C}$) denotes the set of
nonnegative integers (resp., real, complex numbers). For an integer $r>0$,
denote the set $[r]\coloneqq\\{1,\ldots,r\\}$. Uppercase letters (e.g., $A$)
denote matrices, $A_{ij}$ denotes the $(i,j)$th entry of the matrix $A$, and
Curl letters (e.g., $\mathcal{F}$) denote tensors,
$\mathcal{F}_{i_{1},...,i_{m}}$ denotes the $(i_{1},...,i_{m})$th entry of the
tensor $\mathcal{F}$. For a complex matrix $A$, $A^{T}$ denotes its transpose
and $A^{*}$ denotes its conjugate transpose. The Kruskal rank of $A$, for
which we denote $\kappa_{A}$, is the largest number $k$ such that every set of
$k$ columns of $A$ is linearly independent. For a vector $v$, the $(v)_{i}$
denotes its $i$th entry and $\mbox{diag}(v)$ denotes the square diagonal
matrix whose diagonal entries are given by the entries of $v$. The subscript
$v_{s:t}$ denotes the subvector of $v$ whose label is from $s$ to $t$. For a
matrix $A$, the subscript notation $A_{:,j}$ and $A_{i,:}$ respectively denote
its $j$th column and $i$th row. Similar subscript notation is used for
tensors. For two matrices $A,B$, their classical Kronecker product is denoted
as $A\boxtimes B$. For a set $S$, its cardinality is denoted as $|S|$.

For a tensor decomposition for $\mathcal{F}$ such that

(2.1) $\mathcal{F}\,=\,\sum_{s=1}^{r}u^{s,1}\otimes u^{s,2}\otimes...\otimes
u^{s,m}.$

we denote the matrices

$U^{(j)}\,=\,[u^{1,j},...,u^{r,j}],\quad\,j=1,\ldots,m.$

The $U^{(j)}$ is called the $j$th decomposing matrix for $\mathcal{F}$. For
convenience of notation, we denote that

$\displaystyle U^{(1)}\circ\cdots\circ
U^{(m)}=\sum_{i=1}^{r}(U^{(1)})_{:,i}\otimes\ldots\otimes(U^{(m)})_{:,i}.$

Then the above tensor decomposition is equivalent to
$\mathcal{F}=U^{(1)}\circ\cdots\circ U^{(m)}$.

For a matrix $V\in\mathbb{C}^{p\times n_{t}}$, define the matrix-tensor
product

$\mathcal{A}\coloneqq V\times_{t}\mathcal{F}$

is a tensor in $\mathbb{C}^{n_{1}\times...\times n_{t-1}\times p\times
n_{t+1}\times...\times n_{m}}$ such that the $i$th slice of $\mathcal{A}$ is

$\mathcal{A}_{i_{1},...,i_{t-1},:,i_{t+1},...,i_{m}}=V\mathcal{F}_{i_{1},...,i_{t-1},:,i_{t+1},...,i_{m}}.$

###  2.1. Flattening matrices

We partition the dimensions $n_{1},n_{2},...,n_{m}$ into two disjoint groups
$I_{1}$ and $I_{2}$ such that the difference

$\Big{|}\prod_{i\in I_{1}}n_{i}-\prod_{j\in I_{2}}n_{j}\Big{|}$

is minimum. Up to a permutation, we write that
$I_{1}=\\{n_{1},\ldots,n_{k}\\}$, $I_{2}=\\{n_{k+1},\ldots,n_{m}\\}.$ For
convenience, denote that

$\begin{array}[]{l}I=\left\\{\left(\imath_{1},\ldots,\imath_{k}\right):1\leq\imath_{j}\leq
n_{j},j=1,\ldots,k\right\\},\\\
J=\left\\{\left(\imath_{k+1},\ldots,\imath_{m}\right):1\leq\imath_{j}\leq
n_{j},j=k+1,\ldots,m\right\\}.\end{array}$

For a tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times\ldots\times n_{m}}$, the
above partition gives the flattening matrix

(2.2)
$\displaystyle\mbox{Flat}(\mathcal{F})\coloneqq\left(\mathcal{F}_{\imath,\jmath}\right)_{\imath\in
I,\jmath\in J}.$

This gives the most square flattening matrix for $\mathcal{F}$. Let
$\sigma_{r}$ denote the closure of all rank-$r$ tensors in
$\mathbb{C}^{n_{1}\times\ldots\times n_{m}}$, under the Zariski topology (see
[6]). The set $\sigma_{r}$ is an irreducible variety of
$\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$. For a given tensor
$\mathcal{F}\in\sigma_{r}$, it is possible that
$\operatorname{rank}(\mathcal{F})>r.$ This fact motivates the notion of border
rank:

(2.3)
$\operatorname{rank}_{B}(\mathcal{F})=\min\left\\{r:\mathcal{F}\in\sigma_{r}\right\\},$

For every tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times\ldots\times n_{m}}$,
one can show that

(2.4)
$\operatorname{rank}\mbox{Flat}(\mathcal{F})\leq\operatorname{rank}_{B}(\mathcal{F})\leq\operatorname{rank}(\mathcal{F}).$

A property $\mathtt{P}$ is said to hold $generically$ on $\sigma_{r}$ if
$\mathtt{P}$ holds on a Zariski open subset $T$ of $\sigma_{r}$. For such a
property $\mathtt{P}$, each $u\in T$ is called a generic point. Interestingly,
the above three ranks are equal for generic points of $\sigma_{r}$ for a range
of values of $r$.

######  Lemma 2.1.

Let $s$ be the smaller dimension of the matrix $\mbox{Flat}(\mathcal{F})$. For
every $r\leq s$, the equalities

(2.5)
$\operatorname{rank}\mbox{Flat}(\mathcal{F})\,=\,\operatorname{rank}_{B}(\mathcal{F})\,=\,\operatorname{rank}(\mathcal{F})$

hold for tensors $\mathcal{F}$ in a Zariski open subset of $\sigma_{r}$.

###### Proof.

Let $\phi_{1},\ldots,\phi_{\ell}$ be the $r\times r$ minors of the matrix

(2.6) $\displaystyle\mbox{Flat}(\sum_{i=1}^{r}x^{i,1}\otimes\cdots\otimes
x^{i,m}).$

They are homogeneous polynomials in $x^{i,j}(i=1,\ldots,r,j=1,\ldots,m).$ Let
$x$ denote the tuple $\left(x^{1,1},x^{1,2},\ldots,x^{r,m}\right).$ Define the
projective variety in $\mathbb{P}^{r\left(n_{1}+\cdots+n_{m}\right)-1}$

(2.7) $\displaystyle
Z=\left\\{x:\phi_{1}(x)=\cdots=\phi_{\ell}(x)=0\right\\}.$

Then $Y:=\mathbb{P}^{r\left(n_{1}+\cdots+n_{m}\right)-1}\backslash Z$ is a
Zariski open subset of full dimension. Consider the polynomial mapping
$\pi:Y\rightarrow\sigma_{r}$,

(2.8)
$\displaystyle\left(x^{1,1},x^{1,2},\ldots,x^{r,m}\right)\mapsto\sum_{i=1}^{r}\left(x^{i,1}\right)\otimes\cdots\otimes\left(x^{i,m}\right).$

The image $\pi(Y)$ is dense in the irreducible variety $\sigma_{r}.$ So,
$\pi(Y)$ contains a Zariski open subset $\mathscr{Y}$ of $\sigma_{r}$ (see
[35]). For each $\mathcal{F}\in\mathscr{Y}$, there exists $u\in Y$ such that
$\mathcal{F}=\pi(u)$. Because $u\notin Z$, at least one of
$\phi_{1}(u),\ldots,\phi_{\ell}(u)$ is nonzero, and hence
$\operatorname{rank}\mbox{Flat}(\mathcal{F})\geq r.$ By (2.4), we know (2.5)
holds for all $\mathcal{F}\in\mathscr{Y}$ since
$\operatorname{rank}(\mathcal{F})\leq r.$ Since $\mathscr{Y}$ is a Zariski
open subset of $\sigma_{r}$, the lemma holds. ∎

By Lemma 2.1, if $r\leq s$ and $\mathcal{F}$ is a generic tensor in
$\sigma_{r}$, we can use $\operatorname{rank}\mbox{Flat}(\mathcal{F})$ to
estimate $\operatorname{rank}(\mathcal{F})$. However, for a generic
$\mathcal{F}\in\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$ such that
$\operatorname{rank}\mbox{Flat}(\mathcal{F})=r$, we cannot conclude
$\mathcal{F}\in\sigma_{r}$.

###  2.2. Reshaping of tensor decompositions

A tensor $\mathcal{F}$ of order greater than $3$ can be reshaped to another
tensor $\widehat{\mathcal{F}}$ of order $3$. A tensor decomposition of
$\widehat{\mathcal{F}}$ can be converted to a decomposition for $\mathcal{F}$
under certain conditions. In the following, we assume a given tensor
$\mathcal{F}$ has the decomposition (2.1). Suppose the set $\\{1,\ldots,m\\}$
is partitioned into $3$ disjoint subsets

$\\{1,\ldots,m\\}\,=\,I_{1}\cup I_{2}\cup I_{3}.$

Let $p_{i}=|I_{i}|\mathrm{~{}for~{}}i=1,2,3$. For the reshaped vectors

(2.9)
$\left\\{\begin{array}[]{rcll}w^{s,1}&=&u^{s,i_{1}}\boxtimes\cdots\boxtimes
u^{s,i_{p_{1}}}&\mathrm{~{}for~{}}I_{1}\,=\,\\{i_{1},\ldots,i_{p_{1}}\\},\\\
w^{s,2}&=&u^{s,j_{1}}\boxtimes\cdots\boxtimes
u^{s,j_{p_{2}}}&\mathrm{~{}for~{}}I_{2}\,=\,\\{j_{1},\ldots,j_{p_{2}}\\},\\\
w^{s,3}&=&u^{s,k_{1}}\boxtimes\cdots\boxtimes
u^{s,k_{p_{3}}}&\mathrm{~{}for~{}}I_{3}\,=\,\\{k_{1},\ldots,k_{p_{3}}\\},\end{array}\right.$

we get the following tensor decomposition

(2.10) $\widehat{\mathcal{F}}\,=\,\sum_{s=1}^{r}w^{s,1}\otimes w^{s,2}\otimes
w^{s,3}.$

Conversely, for a decomposition like (2.10) for $\widehat{\mathcal{F}}$, if
all $w^{s,1},w^{s,2},w^{s,3}$ can be expressed as rank-$1$ products as in
(2.9), then the equation (2.10) can be reshaped to a tensor decomposition for
$\mathcal{F}$ as in (2.1). When the flattened tensor $\widehat{\mathcal{F}}$
satisfies some conditions, the tensor decomposition of $\widehat{\mathcal{F}}$
is unique. For such a case, we can obtain a tensor decomposition for
$\mathcal{F}$ through the decomposition (2.10). A classical result about the
uniqueness is the Kruskal’s criterion [21].

######  Theorem 2.2.

(Kruskal’s Criterion, [21]) Let $\mathcal{F}=U^{(1)}\circ U^{(2)}\circ
U^{(3)}$ be a tensor with each $U^{(i)}\in\mathbb{C}^{n_{i}\times r}$. Let
$\kappa_{i}$ be the Kruskal rank of $U^{(i)}$, for $i=1,2,3$. If

$2r+2\leq\kappa_{1}+\kappa_{2}+\kappa_{3},$

then $\mathcal{F}$ has a unique rank-$r$ tensor decomposition.

The Kruskal’s Criterion can be generalized for more range of $r$ as in [4].
Assume the dimension $n_{1}\geq n_{2}\geq n_{3}\geq 2$ and the rank $r$ is
such that

$2r+2\leq\mathrm{min}(n_{1},r)+\mathrm{min}(n_{2},r)+\mathrm{min}(n_{3},r),$

or equivalently, for $\delta=n_{2}+n_{3}-n_{1}-2$, $r$ is such that

$r\leq n_{1}+\mathrm{min}\\{\frac{1}{2}\delta,\delta\\}.$

If $\mathcal{F}$ is a generic tensor of rank $r$ as above in the space
$\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$, then $\mathcal{F}$ has a unique
rank-$r$ decomposition. There following is a uniqueness result for reshaped
tensor decompositions.

######  Theorem 2.3.

(Reshaped Kruskal Criterion, [4, Theorem 4.6]) For the tensor space
$\mathbb{C}^{n_{1}\times n_{2}\times\cdots\times n_{m}}$ with $m\geq 3$, let
$I_{1}\cup I_{2}\cup I_{3}=\\{1,2,...,m\\}$ be a union of disjoint sets and
let

$p_{1}=\prod_{i\in I_{1}}n_{i},\quad p_{2}=\prod_{j\in I_{2}}n_{j},\quad
p_{3}=\prod_{k\in I_{3}}n_{k}.$

Suppose $p_{1}\geq p_{2}\geq p_{3}$ and let $\delta=p_{2}+p_{3}-p_{1}-2$.
Assume

(2.11) $r\leq p_{1}+\mathrm{min}\\{\frac{1}{2}\delta,\delta\\}.$

If $\mathcal{F}$ is a generic tensor of rank $r$ in $\mathbb{C}^{n_{1}\times
n_{2}\times\cdots\times n_{m}}$, then the reshaped tensor
$\widehat{\mathcal{F}}\in\mathbb{C}^{p_{1}\times p_{2}\times p_{3}}$ as in
(2.10) has a unique rank-$r$ decomposition.

##  3\. Generating Polynomials

Generating polynomials can be used to compute tensor rank decompositions. We
consider tensors whose ranks $r$ are not bigger than the highest dimension,
say, $r\leq n_{1}$ where $n_{1}$ is the biggest of $n_{1},\ldots,n_{m}$.
Denote the indeterminate vector variables

$\mathbf{x_{1}}=(x_{1,1},...x_{1,n_{1}}),\,\ldots,\,\mathbf{x_{m}}=(x_{m,1},...x_{m,n_{m}}).$

A tensor in $\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$ can be labelled by
monomials like $x_{1,i_{1}}x_{2,i_{2}}...x_{m,i_{m}}$. Let

(3.1)
$\begin{array}[]{rcl}\mathbb{M}&\coloneqq&\big{\\{}x_{1,i_{1}}...x_{m,i_{m}}\;|\;1\leq
i_{j}\leq n_{j}\big{\\}}.\end{array}$

For a subset $J\subseteq\\{1,2,...,m\\}$, denote that

(3.2) $\boxed{\begin{array}[]{rcl}J^{c}&\coloneqq&\\{1,2,...,m\\}\backslash
J,\\\
\mathbb{M}_{J}&\coloneqq&\big{\\{}x_{1,i_{1}}...x_{m,i_{m}}\;|\;x_{j,i_{j}}=1\,~{}\forall\,j\in
J^{c}\big{\\}},\\\ \mathcal{M}_{J}&\coloneqq&{\rm
span}\\{\mathbb{M}_{J}\\}.\end{array}}$

A label tuple $(i_{1},\ldots,i_{m})$ is uniquely determined by the monomial
$x_{1,i_{1}}\cdots x_{m,i_{m}}$. So a tensor
$\mathcal{F}\in\mathbb{C}^{n_{1}\times\ldots\times n_{m}}$ can be equivalently
labelled by monomials such that

(3.3) $\displaystyle\mathcal{F}_{x_{1,i_{1}}\ldots
x_{m,i_{m}}}\coloneqq\mathcal{F}_{i_{1},\ldots,i_{m}}.$

With the new labelling by monomials, define the bi-linear product

(3.4)
$\displaystyle\langle\sum_{\mu\in\mathbb{M}}c_{\mu}\mu,\mathcal{F}\rangle\,\coloneqq\,\sum_{\mu\in\mathbb{M}}c_{\mu}\mathcal{F}_{\mu}.$

In the above, each $c_{\mu}$ is a scalar and $\mathcal{F}$ is labelled by
monomials as in (3.3).

######  Definition 3.1.

([28, 33]) For a subset $J\subseteq\\{1,2,...,m\\}$ and a tensor
$\mathcal{F}\in\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$, a polynomial
$p\in\mathcal{M}_{J}$ is called a generating polynomial for $\mathcal{F}$ if

(3.5) $\langle pq,\mathcal{F}\rangle=0\quad\mbox{for
all}\,\,q\in\mathbb{M}_{J^{c}}.$

The following is an example of generating polynomials.

######  Example 3.2.

Consider the cubic order tensor $\mathcal{F}\in\mathbb{C}^{3\times 3\times 3}$
given as

$\displaystyle\begin{bmatrix}\mathcal{F}_{:,:,1}&\vline&\mathcal{F}_{:,:,2}&\vline&\mathcal{F}_{:,:,3}\\\
\end{bmatrix}=\begin{bmatrix}\begin{matrix}[r]11&20&10\\\ 7&10&5\\\ 12&18&9\\\
\end{matrix}&\vline&\begin{matrix}[r]15&24&12\\\ 15&18&9\\\ 24&30&15\\\
\end{matrix}&\vline&\begin{matrix}[r]7&10&5\\\ 9&10&5\\\ 14&16&8\\\
\end{matrix}\\\ \end{bmatrix}.$

The following is a generating polynomial for $\mathcal{F}$:

$p\coloneqq(2x_{1,1}-x_{1,2})(2x_{2,1}-x_{2,2}).$

Note that $p\in\mathcal{M}_{\\{1,2\\}}$ and for each $i_{3}=1,2,3$

$\displaystyle p\cdot
x_{3,i_{3}}=(2x_{1,1}-x_{1,2})(2x_{2,1}-x_{2,2})x_{3,i_{3}}.$

One can check that for each $i_{3}=1,2,3$,

$4\mathcal{F}_{1,1,i_{3}}-2\mathcal{F}_{1,2,i_{3}}-2\mathcal{F}_{2,1,i_{3}}+\mathcal{F}_{2,2,i_{3}}=0.$

This is because $\begin{bmatrix}4&-2&-2&1\end{bmatrix}$ is orthogonal to

$\begin{bmatrix}\mathcal{F}_{1,1,i_{3}}&\mathcal{F}_{1,2,i_{3}}&\mathcal{F}_{2,1,i_{3}}&\mathcal{F}_{2,2,i_{3}}\end{bmatrix}$

for each $i_{3}=1,2,3$.

Suppose the rank $r\leq n_{1}$ is given. For convenience of notation, denote
the label set

(3.6) $\displaystyle J\coloneqq\\{(i,j,k):1\leq i\leq r,~{}2\leq j\leq
m,~{}2\leq k\leq n_{j}\\}.$

For a matrix $G\in\mathbb{C}^{[r]\times J}$ and a triple $\tau=(i,j,k)\in J$,
define the bi-linear polynomial

(3.7)
$\displaystyle\phi[G,\tau](x)\coloneqq\sum_{\ell=1}^{r}G(\ell,\tau)x_{1,\ell}x_{j,1}-x_{1,i}x_{j,k}\,\,\in\,\mathcal{M}_{\\{1,j\\}}.$

The rows of $G$ are labelled by $\ell=1,2,...,r$ and the columns of $G$ are
labelled by $\tau\in J$. We are interested in $G$ such that $\phi[G,\tau]$ is
a generating polynomial for a tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times
n_{2}\times\ldots\times n_{m}}$. This requires that

$\displaystyle\langle\phi[G,\tau]\cdot\mu,\mathcal{F}\rangle=0\,\quad\mbox{for
all}\,\,\mu\in\mathbb{M}_{\\{1,j\\}^{c}}.$

The above is equivalent to the equation ($\mathcal{F}$ is labelled as in
(3.3))

(3.8)
$\sum_{\ell=1}^{r}G(\ell,\tau)\mathcal{F}_{x_{1,\ell}\cdot\mu}=\mathcal{F}_{x_{1,i}x_{j,k}\cdot\mu}.$

######  Definition 3.3.

([28, 33]) When (3.8) holds for all $\tau\in J$, the $G$ is called a
generating matrix for $\mathcal{F}$.

For given $G$, $j\in\\{2,\ldots,m\\}$ and $k\in\\{2,\ldots,n_{j}\\}$, we
denote the matrix

(3.9)
$M^{j,k}[G]\coloneqq\begin{bmatrix}G(1,(1,j,k))&G(2,(1,j,k))&\dots&G(r,(1,j,k))\\\
G(1,(2,j,k))&G(2,(2,j,k))&\dots&G(r,(2,j,k))\\\ \vdots&\vdots&\ddots&\vdots\\\
G(1,(r,j,k))&G(2,(r,j,k))&\dots&G(r,(r,j,k))\\\ \end{bmatrix}.$

For each $(j,k)$, define the matrix/vector

(3.10)
$\left\\{\begin{array}[]{rcl}A[\mathcal{F},j]&\coloneqq&\Big{(}\mathcal{F}_{x_{1,\ell}\cdot\mu}\Big{)}_{\mu\in\mathbb{M}_{\\{1,j\\}^{c}},1\leq\ell\leq
r},\\\ b[\mathcal{F},j,k]&\coloneqq&\Big{(}\mathcal{F}_{x_{1,\ell}\cdot
x_{j,k}\cdot\mu}\Big{)}_{\mu\in\mathbb{M}_{\\{1,j\\}^{c}},1\leq\ell\leq
r}.\end{array}\right.$

The equation (3.8) is then equivalent to

(3.11) $A[\mathcal{F},j](M^{j,k}[G])^{T}\,=\,b[\mathcal{F},j,k].$

The following is a useful property for the matrices $M^{j,k}[G]$.

######  Theorem 3.4.

([28, 33]) Suppose $\mathcal{F}=\sum_{s=1}^{r}u^{s,1}\otimes...\otimes
u^{s,m}$ for vectors $u^{s,j}\in\mathbb{C}^{n_{j}}$. If $r\leq n_{1}$,
$(u^{s,2})_{1}...(u^{s,m})_{1}\neq 0$, and the first $r$ rows of the first
decomposing matrix

$U^{(1)}\coloneqq[u^{1,1}\,\,\cdots\,\,u^{r,1}]$

are linearly independent, then there exists a $G$ satisfying (3.11) and
satisfying (for all $j\in\\{2,\ldots,m\\}$, $k\in\\{2,\ldots,n_{j}\\}$ and
$s=1,\ldots,r$)

(3.12) $\displaystyle
M^{j,k}[G]\cdot(u^{s,1})_{1:r}=(u^{s,j})_{k}\cdot(u^{s,1})_{1:r}.$

##  4\. Low rank tensor decompositions

Without loss of generality, assume the dimensions are decreasing as $n_{1}\geq
n_{2}\geq\cdots\geq n_{m}.$ We discuss how to compute tensor decomposition for
a tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$ when the
rank $r$ is not bigger than the highest dimension, i.e., $r\leq n_{1}$. As in
Theorem 3.4, the decomposing vectors $(u^{s,1})_{1:r}$ are common eigenvectors
of the matrices $M^{j,k}[G]$, with $(u^{s,j})_{k}$ being the eigenvalues
respectively. This implies that the matrices $M^{j,k}[G]$ are simultaneously
diagonalizable. This property can be used to compute tensor decompositions.

Suppose $G$ is a matrix such that (3.11) holds and $M^{j,k}[G]$ are
simultaneously diagonalizable. That is, there is an invertible matrix
$P\in\mathbb{C}^{r\times r}$ such that all the products $P^{-1}M^{j,k}[G]P$
are diagonal for all $j=2,\ldots,m$ and for all $k=2,\ldots,n_{j}$. Suppose
$M^{j,k}[G]$ are diagonalized such that

(4.1)
$P^{-1}M^{j,k}[G]P\,=\,\mbox{diag}[\lambda_{j,k,1},\lambda_{j,k,2},\ldots,\lambda_{j,k,r}]$

with the eigenvalues $\lambda_{j,k,s}$. For each $s=1,\ldots,r$ and
$j=2,\ldots,m$, denote the vectors

(4.2) $u^{s,j}\,\coloneqq\,(1,\lambda_{j,2,s},\ldots,\lambda_{j,n_{j},s}).$

When $\mathcal{F}$ is rank-$r$, there exist vectors
$u^{1,1},\ldots,u^{r,1}\in\mathbb{C}^{n_{1}}$ such that

(4.3) $\mathcal{F}=\sum_{s=1}^{r}u^{s,1}\otimes u^{s,2}\otimes\cdots\otimes
u^{s,m}.$

The vectors $u^{s,1}$ can be found by solving linear equations after $u^{s,j}$
are obtained for $j=2,\ldots,m$ and $s=1,\ldots,r$. The existence of vectors
$u^{s,1}$ satisfying the tensor decomposition (4.3) is shown in the following
theorem.

######  Theorem 4.1.

Let $\mathcal{F}=V^{(1)}\circ...\circ V^{(m)}$ be a rank-r tensor, for
matrices $V^{(i)}\in\mathbb{C}^{n_{i}\times r}$, such that the first $r$ rows
of $V^{(1)}$ are linearly independent. Suppose $G$ is a matrix satisfying
(3.11) and $P\in\mathbb{C}^{r\times r}$ is an invertible matrix such that all
matrix products $P^{-1}\cdot M^{j,k}[G]\cdot P$ are simultaneously
diagonalized as in (4.1). For $j=2,\ldots,m$ and $s=1,\ldots,r$, let $u^{s,j}$
be vectors given as in (4.2). Then, there must exist vectors
$u^{1,1},\ldots,u^{r,1}\in\mathbb{C}^{n_{1}}$ such that the tensor
decomposition (4.3) holds.

###### Proof.

Since the matrix $P=\begin{bmatrix}p_{1}&\cdots&p_{r}\end{bmatrix}$ is
invertible, there exist scalars $c_{1},\ldots,c_{r}\in\mathbb{C}$ such that

(4.4) $\mathcal{F}_{1:r,1,...,1}=c_{1}p_{1}+c_{2}p_{2}+\cdots+c_{r}p_{r}.$

Consider the new tensor

$\mathcal{H}\,\coloneqq\,\sum_{s=1}^{r}c_{s}p_{s}\otimes
u^{s,2}\otimes\cdots\otimes u^{s,m}.$

In the following, we show that $\mathcal{F}_{1:r,:,...,:}=\mathcal{H}$ and
there exist vectors $u^{1,1},\ldots,u^{r,1}\in\mathbb{C}^{n_{1}}$ satisfying
the equation (4.3).

By Theorem 3.4, one can see that the generating matrix $G$ for $\mathcal{F}$
is also a generating matrix for $\mathcal{H}$, so it holds that

(4.5)
$\langle\phi[G,\tau]p,\mathcal{F}\rangle=\langle\phi[G,\tau]p,\mathcal{H}\rangle=0,\quad\mbox{for
all}\,\,p\in\mathbb{M}_{\\{1,j\\}^{c}}.$

Therefore, we have

(4.6) $\langle\phi[G,\tau]p,\mathcal{H}-\mathcal{F}\rangle=0,\quad\mbox{for
all}\,\,p\in\mathbb{M}_{\\{1,j\\}^{c}}.$

By (4.4), one can see that

(4.7) $\displaystyle(\mathcal{H}-\mathcal{F})_{1:r,1,\ldots,1}=0.$

In (4.6), for each $\tau=(i,2,k)\in J$ and $p=1$, we can get

$(\mathcal{H}-\mathcal{F})_{1:r,:,1,\ldots,1}=0.$

Similarly, for $\tau=(i,2,k)\in J$ and $p=x_{3,j_{3}}$, we can get

$(\mathcal{H}-\mathcal{F})_{1:r,:,:,1,\ldots,1}=0.$

Continuing this, we can eventually get
$\mathcal{H}=\mathcal{F}_{1:r,:,:,\ldots,:}$. Since the matrix
$(V^{(1)})_{1:r,:}$ is invertible, there exists a matrix
$W\in\mathbb{C}^{n_{1}\times r}$ such that $V^{(1)}=W(V^{(1)})_{1:r,:}.$
Observe that

$\mathcal{F}=W\times_{1}\mathcal{F}_{1:r,:,:,\ldots,:}=W\times_{1}\mathcal{H}.$

Let $u^{s,1}=W\cdot(c_{s}p_{s})$ for $s=1,\ldots,r$. Then the tensor
decomposition (4.3) holds. ∎

###  4.1. An algorithm for computing tensor decompositions

Consider a tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times
n_{2}\times\cdots\times n_{m}}$ with a given rank $r$. Recall that the
dimensions are ordered such that $n_{1}\geq n_{2}\geq\cdots\geq n_{m}$. We
discuss how to compute the rank-$r$ tensor decomposition for $\mathcal{F}$.
Recall $A[\mathcal{F},j]$, $b[\mathcal{F},j,k]$ as in (3.10), for $j>1$. Note
that $A[\mathcal{F},j]$ has the dimension $N_{j}\times r$, where

(4.8) $N_{j}\,\coloneqq\,\frac{n_{2}\cdots n_{m}}{n_{j}}.$

If $r\leq N_{j}$, then the matrices $A[\mathcal{F},j]$ have full column rank
for generic cases. For instance, $r\leq N_{3}$ if $m=3$ and $r\leq n_{2}$.
Since $N_{2}$ is the smallest, we often use the matrices $A[\mathcal{F},j]$
for $j\geq 3$. For convenience, denote the label set

(4.9) $\Upsilon\,\coloneqq\,\\{(j,k):3\leq j\leq m,2\leq k\leq n_{j}\\}.$

In the following, we consider the case that $r\leq N_{3}$. For each pair
$(j,k)\in\Upsilon$, the linear system (3.11) has a unique solution, for which
we denote

$Y^{j,k}\,=\,M^{j,k}[G].$

For $j=2$, the equation (3.11) may not have a unique solution if $r>N_{2}$. In
the following, we show how to get the tensor decomposition without using the
matrices $M^{2,k}[G]$. By Theorem 3.4, the matrices $Y^{j,k}$ are
simultaneously diagonalizable, that is, there is an invertible matrix
$P\in\mathbb{C}^{r\times r}$ such that all products $P^{-1}Y^{j,k}P$ are
diagonal for every $(j,k)\in\Upsilon$. Suppose they are diagonalized as

(4.10)
$P^{-1}Y^{j,k}P\,=\,\mbox{diag}[\lambda_{j,k,1},\lambda_{j,k,2},\ldots,\lambda_{j,k,r}]$

with the eigenvalues $\lambda_{j,k,s}$. Write $P$ in the column form

$P\,=\,\begin{bmatrix}p_{1}&\cdots&p_{r}\end{bmatrix}.$

For each $s=1,\ldots,r$ and $j=3,\ldots,m$, let

(4.11) $v^{s,j}\,\coloneqq\,(1,\lambda_{j,2,s},\ldots,\lambda_{j,n_{j},s}).$

Suppose $\mathcal{F}$ has the rank-$r$ decomposition

$\mathcal{F}\,=\,\sum_{s=1}^{r}u^{s,1}\otimes...\otimes u^{s,m}.$

Under the assumptions of Theorem 3.4, the linear system (3.11) has a unique
solution for each pair $(j,k)\in\Upsilon$. For every $j\in\\{3,...,m\\}$,
there exist scalars $c_{s,j},c_{s,1}$ such that

$u^{s,j}\,=\,c_{s,j}v^{s,j},\quad u^{s,1}\,=\,c_{s,1}p_{s}.$

Then, we consider the sub-tensor equation in the vector variables
$y_{1},\ldots,y_{r}\in\mathbb{C}^{n_{2}}$

(4.12) $\mathcal{F}_{1:r,:,\ldots,:}\,=\,\sum_{s=1}^{r}p_{s}\otimes
y_{s}\otimes v^{s,3}\otimes\cdots\otimes v^{s,m}.$

There are $rn_{2}\cdots n_{m}$ equations and $rn_{2}$ unknowns. This
overdetermined linear system has solutions such that

$y_{s}=c_{s,2}u^{s,2},~{}\mathrm{for~{}some}~{}c_{s,2}\in\mathbb{C}.$

After all $y_{s}$ are obtained, we solve the linear equation in
$z_{1},\ldots,z_{r}\in\mathbb{C}^{n_{1}-r}$

(4.13) $\mathcal{F}_{r+1:n_{1},:,\ldots,:}=\sum_{s=1}^{r}z_{s}\otimes
y_{s}\otimes v^{s,3}\otimes\cdots\otimes v^{s,m}.$

After all $y_{s},z_{s}$ are obtained, we choose the vectors ($s=1,\ldots,r$)

$v^{s,1}=\begin{bmatrix}p_{s}\\\ z_{s}\end{bmatrix},\quad v^{s,2}=y_{s}.$

Then we get the tensor decomposition

(4.14) $\mathcal{F}\,=\,\sum_{s=1}^{r}v^{s,1}\otimes
v^{s,2}\otimes\cdots\otimes v^{s,m}.$

Summarizing the above, we get the following algorithm for computing tensor
decompositions when $r\leq n_{1}$ and $r\leq N_{3}$. Suppose the dimensions
are ordered such that $n_{1}\geq n_{2}\geq\cdots\geq n_{m}$.

######  Algorithm 4.2.

(Rank-$r$ tensor decomposition.)

  * Input:

A tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times...\times n_{m}}$ with rank
$r\leq\min(n_{1},N_{3})$.

  * Step 1

For each pair $(j,k)\in\Upsilon$, solve the matrix equation for the solution
$Y^{j,k}$:

(4.15) $A[\mathcal{F},j]Y^{j,k}\,=\,b[\mathcal{F},j,k].$

  * Step 2

Choose generic scalars $\xi_{j,k}$. Then compute the eigenvalue decomposition
$P^{-1}YP=D$ for the matrix

$Y\,\coloneqq\,\frac{1}{\sum\limits_{(j,k)\in\Upsilon}\xi_{j,k}}\sum\limits_{(j,k)\in\Upsilon}\xi_{j,k}Y^{j,k}.$

  * Step 3

For $s=1,\ldots,r$ and $j\geq 3$, let $v^{s,j}$ be the vectors as in (4.11).

  * Step 4

Solve the linear system (4.12) for vectors $y_{1},\ldots,y_{r}$.

  * Step 5

Solve the linear system (4.13) for vectors $z_{1},\ldots,z_{r}$.

  * Step 6

For each $s=1,\ldots,r$, let $v^{s,1}=\begin{bmatrix}p_{s}\\\
z_{s}\end{bmatrix}$ and $v^{s,2}=y_{s}$.

  * Output:

The tensor rank-$r$ decomposition as in (4.14).

The correctness of Algorithm 4.2 is justified as follows.

######  Theorem 4.3.

Suppose $n_{1}\geq n_{2}\geq\cdots\geq n_{m}$ and $r\leq\min(n_{1},N_{3})$ as
in (4.8). For a generic tensor $\mathcal{F}$ of rank-$r$, Algorithm 4.2
produces a rank-$r$ tensor decomposition for $\mathcal{F}$.

###### Proof.

This can be implied by Theorem 4.1. ∎

###  4.2. Tensor decompositions via reshaping

A tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$ can be
reshaped as a cubic order tensor $\widehat{\mathcal{F}}$ as in (2.10). One can
apply Algorithm 4.2 to compute the tensor decomposition (2.10) for
$\widehat{\mathcal{F}}$. If the decomposing vectors $w^{s,1},w^{s,2},w^{s,3}$
can be reshaped to rank-$1$ tensors, then we can convert (2.10) to a tensor
decomposition for $\mathcal{F}$. This is justified by Theorem 2.3, under some
assumptions. A benefit for doing this is that we may be able to compute tensor
decompositions for the case that

$N_{3}<r\leq p_{2},$

with the dimension $p_{2}$ as in Theorem 2.3. This leads to the following
algorithm for computing tensor decompositions.

######  Algorithm 4.4.

(Tensor decompositions via reshaping.) Let $p_{1},p_{2},p_{3}$ be dimensions
as in Theorem 2.3.

  * Input:

A tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times\cdots\times n_{m}}$ with rank
$r\leq p_{2}$.

  * Step 1

Reshape the tensor $\mathcal{F}$ to a cubic tensor
$\widehat{\mathcal{F}}\in\mathbb{C}^{p_{1}\times p_{2}\times p_{3}}$ as in
(2.10).

  * Step 2

Use Algorithm 4.2 to compute the tensor decomposition

(4.16) $\widehat{\mathcal{F}}\,=\,\sum_{s=1}^{r}w^{s,1}\otimes w^{s,2}\otimes
w^{s,3}.$

  * Step 3

If all $w^{s,1},w^{s,2},w^{s,3}$ can be expressed as outer products of
rank-$1$ tensors as in (2.9), then output the tensor decomposition as in
(2.1). If one of $w^{s,1},w^{s,2},w^{s,3}$ cannot be expressed as in (2.9),
then the reshaping does not produce a tensor decomposition for $\mathcal{F}$.

  * Output:

A tensor decomposition for $\mathcal{F}$ as in (2.1).

For Algorithm 4.4, we have a similar conclusion like Theorem 4.3. For
cleanness of the paper, we do not repeat it here.

##  5\. low rank Tensor Approximations

When a tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times...\times n_{m}}$ has the
rank bigger than $r$, the linear systems in Algorithm 4.2 may not be
consistent. However, we can find linear least squares solutions for them. This
gives an algorithm for computing low rank tensor approximations. Recall the
label set $\Upsilon$ as in (4.9). The following is the algorithm.

######  Algorithm 5.1.

(Rank-$r$ tensor approximation.)

  * Input:

A tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times n_{2}\times...\times n_{m}}$
and a rank $r\leq\min(n_{1},N_{3})$.

  * Step 1

For each pair $(j,k)\in\Upsilon$, solve the linear least squares problem

(5.1) $\min\limits_{Y^{j,k}\in\mathbb{C}^{r\times
r}}\quad\bigg{\lVert}A[\mathcal{F},j](Y^{j,k})^{T}-b[\mathcal{F},j,k]\bigg{\rVert}^{2}.$

Let $\hat{Y}^{j,k}$ be an optimizer.

  * Step 2

Choose generic scalars $\xi_{j,k}$ and let

$\hat{Y}[\xi]\,=\,\frac{1}{\sum\limits_{(j,k)\in\Upsilon}\xi_{j,k}}\sum\limits_{(j,k)\in\Upsilon}\xi_{j,k}\hat{Y}^{j,k}.$

Compute the eigenvalue decomposition $\hat{P}^{-1}\hat{Y}[\xi]\hat{P}=\Lambda$
such that $\hat{P}=\begin{bmatrix}\hat{p}_{1}&\cdots&\hat{p}_{r}\end{bmatrix}$
is invertible and $\Lambda$ is diagonal.

  * Step 3

For each pair $(j,k)\in\Upsilon$, select the diagonal entries

$\mbox{diag}[\hat{\lambda}_{j,k,1}\,\,\hat{\lambda}_{j,k,2}\,\,\ldots\,\,\hat{\lambda}_{j,k,r}]\,=\,\mbox{diag}(\hat{P}^{-1}\hat{Y}^{j,k}\hat{P}).$

For each $s=1,\ldots,r$ and $j=3,\ldots,m$, let

$\hat{v}^{s,j}=(1,\hat{\lambda}_{j,2,2},\ldots,\hat{\lambda}_{j,n_{j},s}).$

  * Step 4

Let $(\hat{y}_{1},\ldots,\hat{y}_{r})$ be an optimizer for the following least
squares:

(5.2)
$\min\limits_{(y_{1},\ldots,y_{r})}\quad\bigg{\lVert}\mathcal{F}_{1:r,:,\ldots,:}-\sum_{s=1}^{r}\hat{p}_{s}\otimes
y_{s}\otimes\hat{v}^{s,3}\otimes\cdots\otimes\hat{v}^{s,m}\bigg{\rVert}^{2}.$

  * Step 5

Let $(\hat{z}_{1},\ldots,\hat{z}_{r})$ be an optimizer for the following least
squares:

(5.3)
$\min\limits_{(z_{1},\ldots,z_{r})}\quad\bigg{\lVert}\mathcal{F}_{r+1:n_{1},:,\ldots,:}-\sum_{s=1}^{r}z_{s}\otimes\hat{y}_{s}\otimes\hat{v}^{s,3}\otimes\cdots\otimes\hat{v}^{s,m}\bigg{\rVert}^{2}.$

  * Step 6

Let $\hat{v}^{s,1}=\begin{bmatrix}\hat{p}_{s}\\\ \hat{z}_{s}\end{bmatrix}$ and
$\hat{v}^{s,2}=\hat{y}_{s}$ for each $s=1,\ldots,r$.

  * Output:

The rank-$r$ approximation tensor

(5.4)
$\mathcal{X}^{gp}\,\coloneqq\,\sum_{s=1}^{r}\hat{v}^{s,1}\otimes\hat{v}^{s,2}\otimes...\otimes\hat{v}^{s,m}.$

If $\mathcal{F}$ is sufficiently close to a rank-$r$ tensor, then
$\mathcal{X}^{gp}$ is expected to be a good rank-$r$ approximation.
Mathematically, the tensor $\mathcal{X}^{gp}$ produced by Algorithm 5.1 may
not be the best rank-$r$ approximation. However, in computational practice, we
can use (5.4) as a starting point to solve the nonlinear least squares
optimization

(5.5)
$\min\limits_{(u^{s,1},\ldots,u^{s,m})}\bigg{\lVert}\mathcal{F}-\sum_{s=1}^{r}u^{s,1}\otimes
u^{s,2}\otimes...\otimes u^{s,m}\bigg{\rVert}^{2}.$

to improve the approximation quality. Let $\mathcal{X}^{opt}$ be the rank-$r$
approximation tensor

(5.6) $\mathcal{X}^{opt}\,\coloneqq\,\sum_{s=1}^{r}u^{s,1}\otimes
u^{s,2}\otimes...\otimes u^{s,m}$

which is an optimizer to (5.5) obtained by nonlinear optimization methods with
$\mathcal{X}^{opt}$ as the initial point.

###  5.1. Approximation error analysis

Suppose the tensor $\mathcal{F}$ has the best (or nearly best) rank-$r$
approximation

(5.7)
$\mathcal{X}^{bs}\,\coloneqq\,\sum_{s=1}^{r}(x^{s,1})\otimes(x^{s,2})\otimes...\otimes(x^{s,m}).$

Let $\mathcal{E}$ be the tensor such that

(5.8) $\displaystyle\mathcal{F}=\mathcal{X}^{bs}+\mathcal{E}.$

We analyze the approximation performance of $\mathcal{X}^{gp}$ when the
distance $\epsilon=\|\mathcal{E}\|$ is small. For a generating matrix $G$ and
a generic $\xi=(\xi_{j,k})_{(j,k)\in\Upsilon}$, denote that

(5.9)
$M[\xi,G]\,\coloneqq\,\frac{1}{\sum\limits_{(j,k)\in\Upsilon}\xi_{j,k}}\sum\limits_{(j,k)\in\Upsilon}\xi_{j,k}M^{j,k}[G].$

Recall the $A[\mathcal{F},j]$, $b[\mathcal{F},j,k]$ as in (3.10). Note that

(5.10)
$\displaystyle\begin{split}A[\mathcal{F},j]=&~{}A[\mathcal{X}^{bs},j]+A[\mathcal{\mathcal{E}},j],\\\
b[\mathcal{F},j,k]=&~{}b[\mathcal{X}^{bs},j,k]+b[\mathcal{\mathcal{E}},j,k].\end{split}$

Suppose $\left(x^{s,j}\right)_{1}\neq 0$ for $j=2,\ldots,m$.

######  Theorem 5.2.

Let $\mathcal{X}^{gp}$ be produced by Algorithm 5.1. Let
$\mathcal{F},\mathcal{X}^{bs},\mathcal{X}^{opt},\mathcal{E},x^{s,j},\xi_{j,k}$
be as above. Assume the following conditions hold:

  1. (i)

The subvectors $(x^{1,1})_{1:r},\ldots,(x^{r,1})_{1:r}$ are linearly
independent.

  2. (ii)

All matrices $A[\mathcal{F},j]$ and $A[\mathcal{X}^{bs},j]$ ($3\leq j\leq m$)
have full column rank.

  3. (iii)

The first entry $\left(x^{s,j}\right)_{1}\neq 0$ for all $j=2,\ldots,m$.

  4. (iv)

The following scalars are pairwisely distinct

(5.11)
$\sum_{(j,k)\in\Upsilon}\xi_{j,k}(x^{1,j})_{k},...,\sum_{(j,k)\in\Upsilon}\xi_{j,k}(x^{r,j})_{k}.$

If the distance $\epsilon=\lVert\mathcal{F}-\mathcal{X}^{bs}\rVert$ is
sufficiently small, then

(5.12)
$\lVert\mathcal{X}^{bs}-\mathcal{X}^{gp}\rVert=O(\epsilon)\quad\quad\quad\lVert\mathcal{F}-\mathcal{X}^{gp}\rVert=O(\epsilon).$

where the constants in the above $O(\cdot)$ only depend on $\mathcal{F}$ and
$\xi$.

###### Proof.

By conditions (i) and (iii) and by Theorem 3.4, there exists a generating
matrix $G^{bs}$ for $\mathcal{X}^{bs}$ such that

(5.13)
$A[\mathcal{X}^{bs},j](M^{j,k}[G^{bs}])^{T}\,=\,b[\mathcal{X}^{bs},j,k]$

for all $j\in\\{2,\ldots,m\\}$ and $k\in\\{2,\ldots,n_{j}\\}$. Note that
$Y^{j,k}$ is the least squares solution to (5.1), so for each
$(j,k)\in\Upsilon$,

$Y^{j,k}=A[\mathcal{F},j]^{\dagger}\cdot b[\mathcal{F},j,k],\quad
M^{j,k}[G^{bs}_{0}]=A[\mathcal{X}^{bs},j]^{\dagger}\cdot
b[\mathcal{X}^{bs},j,k].$

(The super script † denotes the Pseudo-inverse of a matrix.) By (5.8), for
$j=2,\ldots,m$, we have

(5.14)
$\begin{array}[]{l}\left\|A[\mathcal{F},j]-A[\mathcal{X}^{bs},j]\right\|_{F}\leq\left\|\mathcal{F}-\mathcal{X}^{bs}\right\|\leq\epsilon,\\\
\left\|b[\mathcal{F},j,k]-b[\mathcal{X}^{bs},j,k]\right\|_{F}\leq\left\|\mathcal{F}-\mathcal{X}^{bs}\right\|\leq\epsilon.\end{array}$

Hence, by the condition (ii), if $\epsilon>0$ is small enough, we have

(5.15) $\left\|Y^{j,k}-M^{j,k}[G^{bs}]\right\|=O(\epsilon).$

for all $(j,k)\in\Upsilon$. This follows from perturbation analysis for linear
least squares (see [12, Theorem 3.4]).

By (5.7) and Theorem 3.4, for $s=1,\ldots,r$ and $(j,k)\in\Upsilon$, it holds
that

$M^{j,k}[G^{bs}]\left(x^{s,1}\right)_{1:r}\quad=\quad\left(x^{s,j}\right)_{k}\left(x^{s,1}\right)_{1:r}.$

This means that each $\left(x^{s,1}\right)_{1:r}$ is an eigenvector of
$M^{j,k}[G^{bs}]$, associated to the eigenvalue $\left(x^{s,j}\right)_{k}$,
for each $s=$ $1,\ldots,r$. The matrices $M^{j,k}[G^{bs}]$ are simultaneously
diagonalizable, by the condition (i). So $M[\xi,G^{bs}]$ is also
diagonalizable. Note the eigenvalues of $M[\xi,G^{bs}]$ are the sums in
(5.11). They are distinct from each other, by the condition (iv). When
$\epsilon>0$ is small enough, $M[\xi,G^{bs}]$ also has distinct eigenvalues.
Write that

$Q=\begin{bmatrix}(x^{1,1})_{1:r}&\cdots&(x^{r,1})_{1:r}\end{bmatrix}.$

Note that $Q^{-1}M[\xi,G^{bs}]Q=D$ is an eigenvalue decomposition. Up to a
scaling on $\hat{P}$ in algorithm 5.1, it holds that

(5.16)
$\|\hat{p}_{s}-x^{s,1}\|_{2}=O(\epsilon),\quad\|D-\Lambda\|_{F}=O(\epsilon).$

We refer to [2] for the perturbation bounds in (5.16). The constants in the
above $O(\cdot)$ eventually only depend on $\mathcal{F},\xi$.

Note that $(\hat{y}_{s},\ldots,\hat{y}_{r})$ is the least squares solution to
(5.2) and

(5.17) $\mathcal{X}^{bs}_{1:r,:,\ldots,:}=\sum_{s=1}^{r}x^{s,1}\otimes
x^{s,2}\otimes x^{s,3}\otimes\cdots\otimes x^{s,m}.$

Due to perturbation analysis of linear least squares, we also have

(5.18) $\|\hat{y}_{s}-x^{s,2}\|_{2}=O(\epsilon).$

Note that the subvectors $(x^{s,1})_{r+1:n_{1}}$ satisfy the equation

(5.19)
$\mathcal{X}^{bs}_{r+1:n_{1},:,\ldots,:}=\sum_{s=1}^{r}(x^{s,1})_{r+1:n_{1}}\otimes
x^{s,2}\otimes\cdots\otimes x^{s,m}.$

Recall that $(\hat{z}_{1},\ldots,\hat{z}_{r})$ is the least squares solution
to (5.3). Due to perturbation analysis of linear least squares, we further
have the error bound

(5.20) $\|(x^{s,1})_{r+1:n_{1}}-\hat{z}_{s}\|_{2}=O(\epsilon).$

Summarizing the above, we eventually get
$\|\mathcal{X}^{gp}-\mathcal{X}^{bs}\|=O(\epsilon)$, so

$\left\|\mathcal{F}-\mathcal{X}^{gp}\right\|\leq\left\|\mathcal{F}-\mathcal{X}^{bs}\right\|+\left\|\mathcal{X}^{bs}-\mathcal{X}^{gp}\right\|=O(\epsilon).$

The constant for the above $O(\cdot)$ eventually only depends on
$\mathcal{F}$, $\xi$. ∎

###  5.2. Reshaping for low rank approximations

Similar to tensor decompositions, the reshaping trick as in Section 4.2 can
also be used for computing low rank tensor approximations. For $m>3$, a tensor
$\mathcal{F}\in\mathbb{C}^{n_{1}\times n_{2}\times\cdots\times n_{m}}$ can be
reshaped as a cubic tensor $\widehat{\mathcal{F}}\in\mathbb{C}^{p_{1}\times
p_{2}\times p_{3}}$ as in (2.10). Similarly, Algorithm 5.1 can be used to
compute low rank tensor approximations. Suppose the computed rank-$r$
approximating tensor for $\widehat{\mathcal{F}}$ is

(5.21)
$\widehat{\mathcal{X}}^{gp}:=\sum_{s=1}^{r}\hat{w}^{s,1}\otimes\hat{w}^{s,2}\otimes\hat{w}^{s,3}.$

Typically, the decomposing vectors $\hat{w}^{s,1},\hat{w}^{s,2},\hat{w}^{s,3}$
may not be reshaped to rank-$1$ tensors. Suppose the reshaping is such that
$I_{1}\cup I_{2}\cup I_{3}=\\{1,2,...,m\\}$ is a union of disjoint label sets
and the reshaped dimensions are

$p_{1}=\prod_{i\in I_{1}}n_{i},\quad p_{2}=\prod_{i\in I_{2}}n_{i},\quad
p_{3}=\prod_{i\in I_{3}}n_{i}.$

Let $m_{i}=|I_{i}|$ for $i=1,2,3$. By the reshaping, the vectors
$\hat{w}^{s,i}$ can be reshaped back to a tensor $\hat{W}^{s,i}$ of order
$m_{i}$, for each $i=1,2,3$. If $m_{i}=1$, $\hat{W}^{s,i}$ is a vector. If
$m_{i}=2$, we can find a best rank-$1$ matrix approximation for
$\hat{W}^{s,i}$. If $m_{i}\geq 3$, we can apply Algorithm 5.1 with $r=1$ to
get a rank-$1$ approximation for $\hat{W}^{s,i}$. In application, we are
mostly interested in reshaping such that all $m_{i}\leq 2$. Finally, this
produces a rank-$r$ approximation for $\mathcal{F}$.

The following is a low rank tensor approximation algorithm via reshaping
tensors.

######  Algorithm 5.3.

(low rank tensor approximations via reshaping.)

  * Input

A tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times n_{2}\times...\times n_{m}}$
and a rank $r$.

  * Step 1

Reshape $\mathcal{F}$ to a cubic order tensor
$\widehat{\mathcal{F}}\in\mathbb{C}^{p_{1}\times p_{2}\times p_{3}}$.

  * Step 2

Use Algorithm 5.1 to compute a rank-$r$ approximating tensor
$\widehat{\mathcal{X}}^{gp}$ as in (5.21) for $\widehat{\mathcal{F}}$.

  * Step 3

For each $i=1,2,3$, reshape each vector $\hat{w}^{s,i}$ back to a tensor
$\widehat{W}^{s,i}$ of order $m_{i}$ as above.

  * Step 4

For each $i=1,2,3$, compute a rank-$1$ approximating tensor
$\widehat{X}^{s,i}$ for $\widehat{W}^{s,i}$ of order $m_{i}$ as above.

  * Output

Reshape the sum
$\sum\limits_{s=1}^{r}\widehat{X}^{s,1}\otimes\widehat{X}^{s,2}\otimes\widehat{X}^{s,3}$
to a tensor in $\mathbb{C}^{n_{1}\times n_{2}\times...\times n_{m}}$, which is
a rank-$r$ approximation for $\mathcal{F}$.

We can do a similar approximation analysis for Algorithm 5.3 as for Theorem
5.2. For cleanness of the paper, we do not repeat that.

##  6\. Numerical Experiments

In this section, we apply Algorithms 4.2 and 5.1 to compute tensor
decompositions and low rank tensor approximations. We implement these
algorithms in MATLAB 2020b on a workstation with Ubuntu 20.04.2 LTS, Intel®
Xeon(R) Gold 6248R CPU @ 3.00GHz and memory 1TB. For computing low rank tensor
approximations, we use the function $cpd\\_nls$ provided in Tensorlab 3.0 [39]
to solve the nonlinear least squares optimization (5.5). The
$\mathcal{X}^{gp}$ denotes the approximating tensor returned by Algorithm 5.1
and $\mathcal{X}^{opt}$ denotes the approximating tensor obtained by solving
(5.5), with $\mathcal{X}^{gp}$ as the initial point. In our numerical
experiments, if the rank $r$ is unknown, we use the most square flattening
matrix to estimate $r$ as in (2.2) and Lemma 2.1.

######  Example 6.1.

Consider the tensor $\mathcal{F}\in\mathbb{C}^{4\times 4\times 3}$ whose
slices $\mathcal{F}_{:,:,1},\mathcal{F}_{:,:,2},\mathcal{F}_{:,:,3}$ are
respectively

$\displaystyle\begin{bmatrix}\begin{matrix}[r]27&25&35&42\\\ 48&68&80&80\\\
26&24&34&40\\\ 33&41&49&66\\\
\end{matrix}&\vline&\begin{matrix}[r]44&32&52&56\\\ 68&76&100&96\\\
42&30&50&52\\\ 46&46&62&76\\\
\end{matrix}&\vline&\begin{matrix}[r]42&26&48&45\\\ 64&60&88&76\\\
47&27&53&47\\\ 45&37&57&60\\\ \end{matrix}\\\ \end{bmatrix}.$

By Lemma 2.1, the estimated rank is $r=4$. Applying Algorithm 4.2 with $r=4$,
we get the rank-$4$ decomposition $\mathcal{F}=U^{(1)}\circ U^{(2)}\circ
U^{(3)}$, with

$U^{(1)}=\begin{bmatrix}[r]8&6&4&9\\\ 8&12&16&12\\\ 4&6&4&12\\\ 4&12&8&9\\\
\end{bmatrix},~{}U^{(2)}=\begin{bmatrix}[r]1&1&1&1\\\
\frac{1}{2}&1&3&\frac{1}{3}\\\ 1&1&3&1\\\ 1&4&1&\frac{2}{3}\\\
\end{bmatrix},~{}U^{(3)}=\begin{bmatrix}[r]1&1&1&1\\\ 2&1&1&2\\\
1&\frac{2}{3}&\frac{3}{4}&3\\\ \end{bmatrix}.$

######  Example 6.2.

Consider the tensor in $\mathbb{C}^{5\times 4\times 3\times 3}$

$\displaystyle\mathcal{F}=V^{(1)}\circ V^{(2)}\circ V^{(3)}\circ V^{(4)},$

where the matrices $V^{(i)}$ are

$\displaystyle V^{(1)}=\begin{bmatrix}[r]10&5&-9&-5&7\\\ 8&6&-3&-9&7\\\
-9&-1&7&-3&-1\\\ 9&-7&-8&8&-5\\\ -1&10&7&-3&10\\\ \end{bmatrix},\quad
V^{(2)}=\begin{bmatrix}[r]-1&9&-8&8&2\\\ 0&-1&-4&6&8\\\ 7&-7&-2&2&10\\\
2&10&-3&-1&-3\\\ \end{bmatrix},$ $\displaystyle
V^{(3)}=\begin{bmatrix}[r]5&2&-2&-7&3\\\ 9&-3&-7&7&-2\\\ 0&-10&10&6&10\\\
\end{bmatrix},\quad V^{(4)}=\begin{bmatrix}[r]8&2&-7&10&-5\\\ 4&-8&4&-6&-10\\\
5&0&7&-1&-2\\\ \end{bmatrix}.$

By Lemma 2.1, the estimated rank $r=5$. Applying Algorithm 4.2 with $r=5$, we
get the rank-$5$ tensor decomposition $\mathcal{F}=U^{(1)}\circ U^{(2)}\circ
U^{(3)}\circ U^{(4)}$, where the computed matrices $U^{(i)}$ are

$U^{(1)}=\begin{bmatrix}[r]-400&180&1008&2800&-210\\\
-320&216&336&5040&-210\\\ 360&-36&-784&1680&30\\\ -360&-252&896&-4480&150\\\
40&360&-784&1680&-300\\\ \end{bmatrix},\quad
U^{(2)}=\begin{bmatrix}[r]1&1&1&1&1\\\
0&-\frac{1}{9}&\frac{1}{2}&\frac{3}{4}&4\\\
-7&-\frac{7}{9}&\frac{1}{4}&\frac{1}{4}&5\\\
-2&\frac{10}{9}&\frac{3}{8}&-\frac{1}{8}&-\frac{3}{2}\\\ \end{bmatrix},$
$U^{(3)}=\begin{bmatrix}[r]1&1&1&1&1\\\
\frac{9}{5}&-\frac{3}{2}&\frac{7}{2}&-1&-\frac{2}{3}\\\
0&-5&-5&-\frac{6}{7}&\frac{10}{3}\\\ \end{bmatrix},\quad
U^{(4)}=\begin{bmatrix}[r]1&1&1&1&1\\\
\frac{1}{2}&-4&-\frac{4}{7}&-\frac{3}{5}&2\\\
\frac{5}{8}&0&-1&-\frac{1}{10}&\frac{2}{5}\\\ \end{bmatrix}.$

######  Example 6.3.

Consider the tensor $\mathcal{F}\in\mathbb{C}^{5\times 5\times 4}$ such that

$\displaystyle\mathcal{F}_{i_{1},i_{2},i_{3}}=i_{1}+\frac{i_{2}}{2}+\frac{i_{3}}{3}+\sqrt{i_{1}^{2}+i_{2}^{2}+i_{3}^{2}}$

for all $i_{1},i_{2},i_{3}$ in the corresponding range. The $5$ biggest
singular values of the flattening matrix $\mbox{Flat}(\mathcal{F})$ are

$\displaystyle 109.7393,~{}~{}5.2500,~{}~{}0.1068,~{}~{}8.325\times
10^{-3},~{}~{}3.401\times 10^{-4}.$

Applying Algorithm 5.1 with rank $r=2,3,4,5$, we get the approximation errors

$r$ | 2 | 3 | 4 | 5
---|---|---|---|---
$\|\mathcal{F}-\mathcal{X}^{gp}\|$ | $5.1237\times 10^{-1}$ | $6.8647\times 10^{-2}$ | $1.0558\times 10^{-2}$ | $9.9449\times 10^{-3}$
$\|\mathcal{F}-\mathcal{X}^{opt}\|$ | $1.5410\times 10^{-1}$ | $1.3754\times 10^{-2}$ | $2.6625\times 10^{-3}$ | $4.9002\times 10^{-4}$

For the case $r=3$, the computed approximating tensor by Algorithm 5.1 and by
solving (5.5) is $U^{(1)}\circ U^{(2)}\circ U^{(3)}$, with

$\displaystyle U^{(1)}=$
$\displaystyle\begin{bmatrix}[r]-0.4973&-7.6813&11.7465\\\
-0.2525&-6.9651&12.4970\\\ -0.0872&-6.0497&13.2858\\\
-0.0132&-5.0521&14.1423\\\ -0.0010&-4.0469&15.0771\\\ \end{bmatrix},\quad
U^{(2)}=\begin{bmatrix}[r]1.0000&1.0000&1.0000\\\ 0.5058&0.9211&1.0306\\\
0.1713&0.8167&1.0649\\\ 0.0262&0.7003&1.1042\\\ 0.0136&0.5807&1.1490\\\
\end{bmatrix},$ $\displaystyle U^{(3)}=$
$\displaystyle\begin{bmatrix}[r]1.0000&1.0000&1.0000\\\
0.5075&0.9289&1.0216\\\ 0.1756&0.8323&1.0469\\\ 0.0399&0.7231&1.0771\\\
\end{bmatrix}.$

######  Example 6.4.

Consider the tensor $\mathcal{F}\in\mathbb{C}^{6\times 6\times 6\times 5\times
4}$ such that

$\mathcal{F}_{i_{1},i_{2},i_{3},i_{4},i_{5}}=\
arctan(i_{1}+2i_{2}+3i_{3}+4i_{4}+5i_{5}),$

for all $i_{1},i_{2},i_{3},i_{4},i_{5}$ in the corresponding range. The $5$
biggest singular values of the flattening matrix $\mbox{Flat}(\mathcal{F})$
are

$\displaystyle 101.71,~{}~{}~{}7.7529\times 10^{-2},~{}~{}~{}2.2870\times
10^{-3},~{}~{}~{}7.2294\times 10^{-5},~{}~{}~{}2.0633\times 10^{-6}.$

Applying Algorithm 5.1 with rank $r=2,3,4,5$, we get the approximation errors
as follows:

$r$ | 2 | 3 | 4 | 5
---|---|---|---|---
$\|\mathcal{F}-\mathcal{X}^{gp}\|$ | $9.8148\times 10^{-3}$ | $3.1987\times 10^{-3}$ | $5.7945\times 10^{-3}$ | $1.0121\times 10^{-5}$
$\|\mathcal{F}-\mathcal{X}^{opt}\|$ | $5.3111\times 10^{-3}$ | $2.2623\times 10^{-4}$ | $3.0889\times 10^{-5}$ | $1.7523\times 10^{-6}$

For the case $r=3$, the computed approximating tensor by Algorithm 5.1 and by
solving (5.5) is $U^{(1)}\circ U^{(2)}\circ U^{(3)}\circ U^{(4)}\circ
U^{(5)}$, with

$\displaystyle U^{(1)}=$
$\displaystyle\begin{bmatrix}[r]-0.0134&-0.0347&1.5524\\\
-0.0112&-0.0329&1.5525\\\ -0.0094&-0.0312&1.5526\\\ -0.0079&-0.0295&1.5527\\\
-0.0066&-0.0280&1.5528\\\ -0.0056&-0.0265&1.5529\\\ \end{bmatrix},\quad
U^{(2)}=\begin{bmatrix}[r]1.0000&1.0000&1.0000\\\ 0.7011&0.8992&1.0001\\\
0.4939&0.8080&1.0003\\\ 0.3485&0.7260&1.0004\\\ 0.2459&0.6523&1.0006\\\
0.1734&0.5861&1.0007\\\ \end{bmatrix},$ $\displaystyle U^{(3)}=$
$\displaystyle\begin{bmatrix}[r]1.0000&1.0000&1.0000\\\
0.5886&0.8523&1.0002\\\ 0.3490&0.7258&1.0004\\\ 0.2064&0.6183&1.0006\\\
0.1214&0.5269&1.0008\\\ 0.0715&0.4489&1.0011\\\ \end{bmatrix},\quad
U^{(4)}=\begin{bmatrix}[r]1.0000&1.0000&1.0000\\\ 0.4949&0.8078&1.0003\\\
0.2463&0.6521&1.0006\\\ 0.1211&0.5269&1.0008\\\ 0.0596&0.4256&1.0011\\\
\end{bmatrix},$ $\displaystyle U^{(5)}=$
$\displaystyle\begin{bmatrix}[r]1.0000&1.0000&1.0000\\\
0.4161&0.7656&1.0003\\\ 0.1730&0.5862&1.0007\\\ 0.0711&0.4489&1.0011\\\
\end{bmatrix}.$

######  Example 6.5.

As in Theorem 5.2, we have shown that if the tensor to be approximated is
sufficiently close to a rank-$r$ tensor, then the computed rank-$r$
approximation $\mathcal{X}^{gp}$ is quasi-optimal. It can be further improved
to a better approximation $\mathcal{X}^{opt}$ by solving the nonlinear
optimization (5.5). In this example, we explore the numerical performance of
Algorithms 5.1 and 5.3 for computing low rank tensor approximations. For the
given dimensions $n_{1},\ldots,n_{m}$, we generate the tensor

$\mathcal{R}\,=\,\sum_{s=1}^{r}u^{s,1}\otimes u^{s,2}\otimes\cdots\otimes
u^{s,m},$

where each $u^{s,j}\in\mathbb{C}^{n_{j}}$ is a complex vector whose real and
imaginary parts are generated randomly, obeying the Gaussian distribution. We
perturb $\mathcal{R}$ by another tensor $\mathcal{E}$, whose entries are also
generated with the Gaussian distribution. We scale the perturbing tensor
$\mathcal{E}$ to have a desired norm $\epsilon$. The tensor $\mathcal{F}$ is
then generated as

$\mathcal{F}\,=\,\mathcal{R}+\mathcal{E}.$

We choose $\epsilon$ to be one of $10^{-2},10^{-4},10^{-6}$, and use the
relative errors

$\text{$\rho$\\_gp}=\frac{\left\|\mathcal{F}-\mathcal{X}^{gp}\right\|}{\|\mathcal{E}\|},\quad\text{$\rho$\\_opt
}=\frac{\left\|\mathcal{F}-\mathcal{X}^{opt}\right\|}{\|\mathcal{E}\|}$

to measure the approximation quality of $\mathcal{X}^{\text{gp }}$,
$\mathcal{X}^{\text{opt }}$ respectively. For each case of
$(n_{1},\ldots,n_{m}),r$ and $\epsilon$, we generate $10$ random instances of
$\mathcal{R},\mathcal{F},\mathcal{E}$. For the case
$(n_{1},\ldots,n_{m})=(20,20,20,20,10)$, Algorithm 5.3 is used to compute
$\mathcal{X}^{gp}$. All other cases are solved by Algorithm 5.1. The
computational results are reported in Tables 1. For each case of
$(n_{1},\ldots,n_{m})$ and $r$, we also list the median of above relative
errors and the average CPU time (in seconds). The t_gp and t_opt denote the
average CPU time (in seconds) for Algorithms 5.1/5.3 and for solving (5.5)
respectively.

Table 1. Computational performance of Algorithms 5.1 and 5.3 and of nonlinear
optimization (5.5).

$r$ | $\epsilon$ |  $\rho$_gp | t_gp |  $\rho$_opt | t_opt | $r$ | $\epsilon$ |  $\rho$_gp | t_gp |  $\rho$_opt | t_opt
---|---|---|---|---|---|---|---|---|---|---|---
($n_{1}$,$n_{2}$,$n_{3}$)=(50,50,50) |  ($n_{1}$,$n_{2}$,$n_{3}$)=(60, 50, 40)
10 | $10^{-2}$ | 1.63 | 0.08 | 0.99 | 1.57 | 15 | $10^{-2}$ | 17.49 | 0.19 | 0.99 | 2.17
$10^{-4}$ | 6.32 | 0.10 | 0.99 | 1.16 | $10^{-4}$ | 10.80 | 0.15 | 0.99 | 1.36
$10^{-6}$ | 3.84 | 0.09 | 0.99 | 0.83 | $10^{-6}$ | 5.16 | 0.20 | 0.99 | 1.10
20 | $10^{-2}$ | 25.83 | 0.29 | 0.99 | 2.99 | 30 | $10^{-2}$ | 28.70 | 0.40 | 0.98 | 6.95
$10^{-4}$ | 5.41 | 0.28 | 0.99 | 1.99 | $10^{-4}$ | 15.77 | 0.37 | 0.98 | 3.61
$10^{-6}$ | 30.41 | 0.29 | 0.99 | 1.49 | $10^{-6}$ | 50.96 | 0.37 | 0.98 | 2.27
30 | $10^{-2}$ | 27.91 | 0.50 | 0.98 | 7.08 | 45 | $10^{-2}$ | 35.48 | 0.61 | 0.97 | 25.73
$10^{-4}$ | 213.82 | 0.43 | 0.98 | 3.73 | $10^{-4}$ | 35.03 | 0.63 | 0.97 | 8.08
$10^{-6}$ | 17.97 | 0.47 | 0.98 | 2.20 | $10^{-6}$ | 34.67 | 0.61 | 0.97 | 5.69
($n_{1}$,$n_{2}$,$n_{3}$)=(100,100,100) |  ($n_{1}$,$n_{2}$,$n_{3}$)=(150,150,150)
20 | $10^{-2}$ | 11.21 | 0.86 | 1.00 | 6.36 | 30 | $10^{-2}$ | 8.59 | 2.92 | 1.00 | 17.17
$10^{-4}$ | 3.48 | 0.85 | 1.00 | 4.24 | $10^{-4}$ | 3.18 | 3.05 | 1.00 | 11.20
$10^{-6}$ | 3.88 | 0.83 | 1.00 | 3.20 | $10^{-6}$ | 4.24 | 3.42 | 1.00 | 11.75
40 | $10^{-2}$ | 24.17 | 1.76 | 0.99 | 17.80 | 60 | $10^{-2}$ | 49.80 | 6.04 | 1.00 | 87.31
$10^{-4}$ | 11.60 | 1.65 | 0.99 | 11.02 | $10^{-4}$ | 13.77 | 5.89 | 1.00 | 24.96
$10^{-6}$ | 11.09 | 1.61 | 0.99 | 7.97 | $10^{-6}$ | 17.49 | 6.07 | 1.00 | 18.81
60 | $10^{-2}$ | 18.71 | 3.40 | 0.99 | 28.16 | 90 | $10^{-2}$ | 29.44 | 10.64 | 0.99 | 98.78
$10^{-4}$ | 26.28 | 3.41 | 0.99 | 17.25 | $10^{-4}$ | 152.49 | 10.53 | 0.99 | 43.58
$10^{-6}$ | 19.12 | 3.49 | 0.99 | 13.14 | $10^{-6}$ | 17.01 | 10.06 | 0.99 | 26.98
($n_{1}$,$n_{2}$,$n_{3}$,$n_{4}$)=(20,20,20,20,10) |  ($n_{1}$,$n_{2}$,$n_{3}$,$n_{4}$)=(60,50,40,30)
24 | $10^{-2}$ | 37.93 | 0.88 | 1.00 | 45.56 | 20 | $10^{-2}$ | 31.42 | 2.78 | 1.00 | 31.16
$10^{-4}$ | 9.10 | 0.92 | 1.00 | 15.86 | $10^{-4}$ | 1.17 | 2.76 | 1.00 | 9.39
$10^{-6}$ | 715.16 | 0.91 | 1.00 | 15.63 | $10^{-6}$ | 4.14 | 2.79 | 1.00 | 9.48
48 | $10^{-2}$ | 166.00 | 1.95 | 1.00 | 270.56 | 40 | $10^{-2}$ | 6.99 | 7.52 | 1.00 | 31.81
$10^{-4}$ | 161.62 | 1.93 | 1.00 | 40.63 | $10^{-4}$ | 2.58 | 7.32 | 1.00 | 20.07
$10^{-6}$ | 52.01 | 1.93 | 1.00 | 21.71 | $10^{-6}$ | 2.49 | 7.22 | 1.00 | 20.22
72 | $10^{-2}$ | 73.70 | 3.10 | 1.00 | 102.90 | 60 | $10^{-2}$ | 11.48 | 9.83 | 1.00 | 48.08
$10^{-4}$ | 113.13 | 3.06 | 1.00 | 70.13 | $10^{-4}$ | 6.38 | 9.80 | 1.00 | 38.97
$10^{-6}$ | 34.28 | 3.03 | 1.00 | 36.72 | $10^{-6}$ | 16.35 | 9.76 | 1.00 | 30.38

In the following, we give a comparison with the generalized eigenvalue
decomposition (GEVD) method, which is a classical one for computing tensor
decompositions when the rank $r\leq n_{2}$. We refer to [23, 34] for the work
about the GEVD method. Consider a cubic order tensor
$\mathcal{F}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$ with $n_{1}\geq
n_{2}\geq n_{3}$. Suppose $\mathcal{F}=U^{(1)}\circ U^{(2)}\circ U^{(3)}$ is a
rank-$r$ decomposition and $r\leq n_{2}$. Assume its first and second
decomposing matrices $U^{(1)},U^{(2)}$ have full column ranks and the third
decomposing matrix $U^{(3)}$ does not have colinear columns. Denote the slice
matrices

(6.1) $F_{1}\,\coloneqq\,\mathcal{F}_{1:r,1:r,1},\quad
F_{2}\,\coloneqq\,\mathcal{F}_{1:r,1:r,2}.$

One can show that

(6.2) $\displaystyle
F_{1}=U^{(1)}_{1:r,:}\cdot\mbox{diag}(U^{(3)}_{1,:})\cdot(U^{(2)}_{1:r,:})^{T},\quad
F_{2}=U^{(1)}_{1:r,:}\cdot\mbox{diag}(U^{(3)}_{2,:})\cdot(U^{(2)}_{1:r,:})^{T}.$

This implies that the columns of $(U^{(1)}_{1:r,r})^{-T}$ are generalized
eigenvectors of the matrix pair $(F_{1}^{T},F_{2}^{T})$. Consider the
transformed tensor

(6.3)
$\hat{\mathcal{F}}\,=\,(U^{(1)}_{1:r,r})^{-1}\times_{1}\mathcal{F}_{1:r,:,:},$

For each $s=1,\ldots,r$, the slice
$\hat{\mathcal{F}}_{s,:,:}=U^{(2)}_{:,s}\cdot(U^{(3)}_{:,s})^{T}$ is a
rank-$1$ matrix. The matrices $U^{(2)}$, $U^{(3)}$ can be obtained by
computing rank-$1$ decompositions for the slices $\hat{\mathcal{F}}_{s,:,:}$.
After this is done, we can solve the linear system

(6.4) $U^{(1)}\circ U^{(2)}\circ U^{(3)}=\mathcal{F}$

to get the matrix $U^{(1)}$. The following is the GEVD method for computing
cubic order tensor decompositions when the rank $r\leq n_{2}$.

######  Algorithm 6.6.

(The GEVD method.)

  * Input:

A tensor $\mathcal{F}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$ with the
rank $r\leq n_{2}$.

  * 1.

Formulate the tensor $\hat{\mathcal{F}}$ as in (6.3).

  * 2.

For $s=1,\ldots,r,$ compute $U^{(2)}_{:,s}$, $U^{(3)}_{:,s}$ from the rank-$1$
decomposition of the matrix $\hat{\mathcal{F}}_{s,:,:}$.

  * 3.

Solve the linear system (6.4) to get $U^{(1)}$.

  * Output:

The decomposing matrices $U^{(1)},U^{(2)},U^{(3)}.$

We compare the performance of Algorithm 4.2 and Algorithm 6.6 for randomly
generated tensors with the rank $r\leq n_{2}$. We generate
$\mathcal{F}=U^{(1)}\circ U^{(2)}\circ U^{(3)}$ such that each
$U^{(i)}\in\mathbb{C}^{n_{i}\times r}$. The entries of $U^{(i)}$ are randomly
generated complex numbers. Their real and imaginary parts are randomly
generated, obeying the Gaussian distribution. For each case of
$(n_{1},...,n_{m})$ and $r$, we generate $20$ random instances of
$\mathcal{F}$. Algorithm 6.6 is implemented by the function $cpd\\_gevd$ in
the software Tensorlab. All the tensor decompositions are computed correctly
by both methods. The average CPU time (in seconds) for Algorithm 4.2 is
denoted as time-gp, while the average CPU time for the GEVD method is denoted
as time-gevd. The computational results are reported in Table 2. The numerical
experiments show that Algorithm 4.2 is more computationally efficient than
Algorithm 6.6.

Table 2. A comparison for the performance of Algorithms 4.2 and 6.6. ($n_{1}$,$n_{2}$,$n_{3}$) | $r$ | time-gevd | time-gp
---|---|---|---
(40,30,30) | 30 | 0.91 | 0.29
(50,50,50) | 50 | 4.77 | 0.85
(100,100,100) | 80 | 12.17 | 5.54
(150,150,150) | 100 | 79.85 | 13.30
(200,200,200) | 120 | 161.83 | 25.71
(250,250,250) | 140 | 285.03 | 55.71
(300,300,300) | 100 | 306.64 | 61.38
(400,400,400) | 180 | 934.15 | 271.21
(500,500,500) | 200 | 1688.98 | 539.75

##  7\. Conclusions

.

This paper gives computational methods for computing low rank tensor
decompositions and approximations. The proposed methods are based on
generating polynomials. For a generic tensor of rank
$r\leq\mathrm{min}(n_{1},N_{3})$, its tensor decomposition can be obtained by
Algorithm 4.2 . Under some general assumptions, we show that if a tensor is
sufficiently close to a low rank one, then the low rank approximating tensor
produced by Algorithm 5.1 is quasi-optimal. Numerical experiments are
presented to show the efficiency of the proposed methods.

Acknowledgement Jiawang Nie is partially supported by the NSF grant
DMS-2110780. Li Wang is partially supported by the NSF grant DMS-2009689.

## References

  * [1] P. Breiding and N. Vannieuwenhoven, A Riemannian trust region method for the canonical tensor rank approximation problem, SIAM J. Optim. 28(3), 2435–2465, 2018.
  * [2] F. Chatelin, Eigenvalues of Matrices: revised edition, SIAM, 2012.
  * [3] L. Chiantini, G. Ottaviani, and N. Vannieuwenhoven, An algorithm for generic and low rank specific identifiability of complex tensors, SIAM J. Matrix Anal. Appl. 35(4), 1265–1287, 2014.
  * [4] L. Chiantini, G. Ottaviani, and N. Vannieuwenhoven, Effective criteria for specific identifiability of tensors and forms, SIAM J. Matrix Anal. Appl. 38(2), 656–681, 2017.
  * [5] P. Comon, X. Luciani, and A. L. F. de Almeida, Tensor decompositions, alternating least squares and other tales, J. of Chemometrics 23, 393–405, 2009.
  * [6] D. Cox, J. Little, and D. O’Shea, Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra, Springer Science & Business Media, 2013.
  * [7] C-F. Cui, Y-H. Dai and J. Nie, All real eigenvalues of symmetric tensors, SIAM J. Matrix Anal. Appl. 35(4), 1582–1601, 2014.
  * [8] L. de Lathauwer, B. de Moor and J. Vandewalle, On the best rank-$1$ and rank-$(R_{1},R_{2},\ldots,R_{n})$ approximation of higher-order tensors, SIAM J. Matrix Anal. Appl. 21(4), 1324–1342, 2000.
  * [9] L. de Lathauwer, B. de Moor and J. Vandewalle, Computation of the canonical decomposition by means of a simultaneous generalized Schur decomposition, SIAM J. Matrix Anal. Appl. 26(2), 295–327, 2004.
  * [10] L. de Lathauwer, A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization, SIAM J. Matrix Anal. Appl. 28(3), 642–666, 2006.
  * [11] V. de Silva and L.-H. Lim, Tensor rank and the ill-posedness of the best low rank approximation problem, SIAM J. Matrix Anal. Appl. 30(3), 1084–1127, 2008.
  * [12] J. Demmel, Applied Numerical Linear Algebra, SIAM, 1997.
  * [13] I. Domanov and L. de Lathauwer, Canonical polyadic decomposition of third-order tensors: Reduction to generalized eigenvalue decomposition, SIAM J. Matrix Anal. Appl. 35(2), 636–660, 2014.
  * [14] J. Fan, J. Nie and A. Zhou, Tensor eigenvalue complementarity problems, _Mathematical Programming_ , 170(2), 507–539, 2018.
  * [15] S. Friedland and V. Tammali, Low-rank approximation of tensors, Numerical algebra, matrix theory, differential-algebraic equations and control theory, 377-411, Springer, 2015.
  * [16] S. Friedland and L. Wang, Spectral norm of a symmetric tensor and its computation, Mathematics of Computation 89, 2175–2215, 2020.
  * [17] Y. Guan, M. T. Chu, and D. Chu, Convergence analysis of an svd-based algorithm for the best rank-$1$ tensor approximation, Linear Algebra and its Applications 555, 53–69, 2018.
  * [18] B. Guo, J. Nie, and Z. Yang, Learning diagonal gaussian mixture models and incomplete tensor decompositions, Vietnam J. Math. (2021). doi.org/10.1007/s10013-021-00534-3
  * [19] L. Qi and S. Hu, Spectral norm and nuclear norm of a third order tensor, Preprint, 2019. arXiv:1909.01529
  * [20] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Review 51(3):455–500, 2009.
  * [21] J. B. Kruskal, Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics, Linear Algebra Appl 18(2), 95–138, 1977.
  * [22] J. Landsberg, Tensors: Geometry and Applications, Grad. Stud. Math., 128, AMS, Providence, RI, 2012.
  * [23] S.E. Leurgans, R.T. Ross and R.B. Abel, A decomposition for three-way arrays, SIAM J. Matrix Anal. Appl. 14, 1064–1083, 1993.
  * [24] L.-H. Lim, Tensors and hypermatrices, in: L. Hogben (Ed.), Handbook of linear algebra (2nd Ed.), CRC Press, Boca Raton, 2013.
  * [25] X. Mao, G. Yuan and Y. Yang, A self-adaptive regularized alternating least squares method for tensor decomposition problems, _Analysis and Applications_ 18 (01), 129–147, 2020.
  * [26] J. Nie, Generating polynomials and symmetric tensor decompositions, Found. Comput. Math. 17(2), 423–465, 2017.
  * [27] J. Nie, Low rank symmetric tensor approximations, SIAM J. Matrix Anal. Appl. 38(4), 1517–1540, 2017.
  * [28] J. Nie, Nearly Low Rank Tensors and Their Approximations, Preprint, 2014. arXiv:1412.7270.
  * [29] J. Nie and L. Wang, Semidefinite relaxations for best rank-$1$ tensor approximations, SIAM J. Matrix Anal. Appl. 35(3), 1155–1179, 2014.
  * [30] J. Nie and Z. Yang, Hermitian tensor decompositions, SIAM J. Matrix Anal. Appl. 41(3), 1115–1144, 2020.
  * [31] J. Nie and X. Zhang, Real eigenvalues of nonsymmetric tensors, Comp. Opt. and Appl., 70(1), 1–32, 2018.
  * [32] J. Nie, Z. Yang and X. Zhang, A complete semidefinite algorithm for detecting copositive matrices and tensors, _SIAM J. Optim._ 28(4), 2902–2921, 2018.
  * [33] J. Nie, L. Wang, and Z. Zheng. Higher Order Correlation Analysis for Multi-View Learning, Pacific Journal of Optimization, to appear, 2022.
  * [34] E. Sanchez and B. R. Kowalski, Tensorial resolution: A direct trilinear decomposition, J. Chemometrics 4, 29–45, 1990.
  * [35] I. Schafarevich, _Basic Algebraic Geometry I: Varieties in Projective Space_ , Springer-Verlag, 1988.
  * [36] L. Sorber, M. van Barel, and L. de Lathauwer,  Optimization-based algorithms for tensor decompositions: Canonical polyadic decomposition, decomposition in rank-$(L_{r},L_{r},1)$ terms, and a new generalization, SIAM J. Optim. 23(2), 695–720, 2013.
  * [37] N. Vervliet, O. Debals, L. Sorber, M. van Barel, and L. de Lathauwer, Tensorlab 3.0, March 2016.
  * [38] S. Telen and N. Vannieuwenhoven, Normal forms for tensor rank decomposition, Preprint, 2021. arXiv:2103.07411
  * [39] N. Vervliet, O. Debals and L. de Lathauwer, Tensorlab 3.0 — numerical optimization strategies for large-scale constrained and coupled matrix/tensor factorization, In 2016 50th Asilomar Conference on Signals, Systems and Computers, 1733–1738, 2016.
  * [40] Y. Yang, The epsilon-alternating least squares for orthogonal low-rank tensor approximation and its global convergence, SIAM J. Matrix Anal. Appl. 41(4), 1797–1825, 2020.