# Recovering orthogonal tensors under arbitrarily strong, but locally
correlated, noise

Oscar Mickelin  Department of Mathematics, Massachusetts Institute of
Technology, Massachusetts, USA<EMAIL_ADDRESS>and  Sertac Karaman
Department of Aeronautics and Astronautics, Massachusetts Institute of
Technology, Massachusetts, USA<EMAIL_ADDRESS>
###### Abstract.

We consider the problem of recovering an orthogonally decomposable tensor with
a subset of elements distorted by noise with arbitrarily large magnitude. We
focus on the particular case where each mode in the decomposition is corrupted
by noise vectors with components that are correlated locally, i.e., with
nearby components. We show that this deterministic tensor completion problem
has the unusual property that it can be solved in polynomial time if the rank
of the tensor is sufficiently large. This is the polar opposite of the low-
rank assumptions of typical low-rank tensor and matrix completion settings. We
show that our problem can be solved through a system of coupled Sylvester-like
equations and show how to accelerate their solution by an alternating solver.
This enables recovery even with a substantial number of missing entries, for
instance for $n$-dimensional tensors of rank $n$ with up to $40\%$ missing
entries.

###### Key words and phrases:

Tensors, tensor completion, compressed tensor formats, canonical
decomposition, orthogonally decomposable tensors.

###### 2010 Mathematics Subject Classification:

Primary 65F99, 15A69

##  1\. Introduction

We consider the problem of reconstructing a tensor $T$ in
$\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ of the form

(1.1) $T=\sum_{k=1}^{r}v_{k1}\otimes v_{k2}\otimes\cdots\otimes v_{kd},$

from observations where its elements have been corrupted by noise. We will be
interested in the case of orthogonally decomposable tensors $T$, i.e., one or
more of the sets $\\{v_{1j},v_{2j},\ldots,v_{rj}\\}$ consists of orthogonal
vectors, for some index $j$ in $1,\ldots,d$. These tensors arise for instance
in independent component analysis [7] and learning latent variable models [1].

The tractability of this problem strongly depends on the specific kind of
noise considered. For general noise terms $E$, the problem becomes that of
finding an optimal orthogonally decomposable approximation to the perturbed
tensor $T+E$. This is in general NP-hard [13, 20], and several different
approximation techniques have been proposed in the recent literature.
Approaches using perturbations of tensor power iteration [31, 1, 21, 22] enjoy
useful stability properties and can therefore cope with noise terms $E$ of
bounded magnitude. Jacobi-type algorithms [19, 16, 28] and algorithms based on
polar decompositions [6, 14], alternating least-squares [29] or singular
value-decompositions [11] have also seen increased interest in recent years. A
different setting includes tensor completion problems, where the location of
the non-zero elements of $E$ are typically chosen probabilistically.
Reconstruction guarantees are then provided under the assumption that the
tensor $T$ is of sufficiently low rank [8, 17, 9, 15, 30, 10, 26, 2].

The goal of this article is to highlight how the structure of an orthogonally
decomposable tensor enables reconstruction under arbitrarily strong, non-
sparse noise, provided that the non-zero entries of the noise term are
sufficiently structured, and that the rank of the tensor is above a certain
threshold. Our results are meant to complement the tensor completion
literature, since high-rank assumptions are required instead of the typical
low-rank assumptions, and we achieve recovery for deterministic sampling
patterns that cannot be efficiently tackled by e.g., nuclear-norm based
techniques. Although simple, the techniques are surprisingly powerful, and are
able to recover $n$-dimensional tensors of rank $n$ with up to $40\%$ missing
entries. The outlook of the article is most closely related to the problem of
decomposing a matrix into the sum of a matrix of low-rank and a diagonal
matrix [23, 24, 25], but relies on different techniques. We would also like to
mention recent articles on deterministic tensor completion [4, 3, 5, 27]. The
tensor structure in (1.1) encodes significant redundancy, which we will show
enables a simple algebraic approach to exact recovery of the tensor $T$ and
its factors $v_{kj}$. Focusing on a particular locally correlated noise model,
we show how to construct coupled Sylvester-type equations that reconstruct $T$
even with a large number of unknown entries. We also show how an alternating
linear solver then recovers $T$ efficiently. Our approach also easily extends
to different sampling patterns.

The remainder of this article is organized as follows. Section 2 introduces
the notation, and section 3 the problem statement. Section 4 details the
proposed algorithm and section 5 shows numerical results. Implementations of
the algorithms in this paper are publicly available
online.111https://github.com/oscarmickelin/locally-correlated-recovery

##  2\. Notation

For two integers $n$ and $m$ with $m\geq n$, we denote the range of integers
contained (inclusively) between them by

(2.1) $\llbracket n,m\rrbracket:=\\{n,n+1,\ldots,m-1,m\\}.$

We will denote the cardinality of a set of integers $\mathcal{I}$ by
$\lvert\mathcal{I}\rvert$.

The symbols $e_{k}$ for $k=1,\ldots,n$ will be reserved for the standard basis
vectors in $\mathbb{R}^{n}$.

For vectors $v_{1},\ldots,v_{d}$ with $v_{k}$ in $\mathbb{R}^{n_{k}}$, their
tensor product $v_{1}\otimes v_{2}\otimes\cdots\otimes v_{d}$ is an element of
$\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$ defined by

(2.2) $\left(v_{1}\otimes v_{2}\otimes\cdots\otimes
v_{d}\right)(i_{1},\ldots,i_{d}):=v_{1}(i_{1})v_{2}(i_{2})\cdots
v_{d}(i_{d}).$

For two tensors $T$ and $S$ in $\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$,
we define their Euclidean inner product by

(2.3) $\langle
T,S\rangle=\sum_{i_{1}=1}^{n_{1}}\cdots\sum_{i_{d}=1}^{n_{d}}T(i_{1},\ldots,i_{d})S(i_{1},\ldots,i_{d}).$

The Euclidean norm of $T$ is then defined by $\|T\|=\sqrt{\langle
T,T\rangle}$.

For a tensor $T$ in $\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$, and an index
$\ell\in\\{1,\ldots,n_{d}\\}$, the $\ell$th slice of $T$ refers to the tensor
$S$ in $\mathbb{R}^{n_{1}\times\cdots\times n_{d-1}}$ with elements

(2.4) $S(i_{1},\ldots,i_{d-1})=T(i_{1},\ldots,i_{d-1},\ell),$

for $i_{k}=1,\ldots,n_{k}$, where $k=1,\ldots,d-1$.

For a matrix $M$ in $\mathbb{R}^{n\times m}$, denote the $\ell$th row of $M$
by $[M]_{\ell}$. We will also denote the set of non-zero elements of $M$ by

(2.5) $\text{supp}(M)=\Big{\\{}(i,j)\in\llbracket
1,n\rrbracket\times\llbracket 1,m\rrbracket:M(i,j)\neq 0\Big{\\}}.$

For a subset $\mathcal{I}\subseteq\llbracket 1,n\rrbracket\times\llbracket
1,m\rrbracket$, we will denote the set of elements of $M$ contained in the
index set $\mathcal{I}$ by

(2.6)
$M\Big{|}_{\mathcal{I}}:=\Big{\\{}M(i,j)\in\mathbb{R}:(i,j)\in\mathcal{I}\Big{\\}}.$

When accessing matrix elements, we will denote by a colon all the elements in
the corresponding mode, e.g., $M(:,m)$ or $M(n,:)$.

##  3\. Problem statement

###  3.1. Motivation

As a motivation for the sampling patterns treated in this article, consider
the problem of recovering a tensor $T=\sum_{k=1}^{r}a_{k}\otimes b_{k}\otimes
c_{k}$ in $\mathbb{R}^{n\times n\times n}$ from a set corrupted samples of the
form $\sum_{k=1}^{r}u_{k}(t)\otimes v_{k}(t)\otimes w_{k}(t)$, indexed by the
variable $t$. Here

(3.1) $\displaystyle u_{k}(t)$ $\displaystyle=a_{k}+n_{1,k}(t),$ (3.2)
$\displaystyle v_{k}(t)$ $\displaystyle=b_{k}+n_{2,k}(t),$ (3.3)
$\displaystyle w_{k}(t)$ $\displaystyle=c_{k}+n_{3,k}(t),$

where $n_{j,k}(t)$ are noise terms with distributions symmetric around the
origin. Our results allow for terms $n_{j,k}(t)$ with arbitrary magnitude,
provided that their covariance terms obey a certain structure. To illustrate
this statement, denote the covariance matrix of $n_{j,k}(t)$ and $n_{i,k}(t)$
by $Q_{ijk}:=\mathbb{E}_{t}\left[n_{i,k}(t)n_{j,k}(t)^{T}\right]$. Since the
distributions of the $n_{j,k}(t)$ are symmetric around the origin, the
corrupted tensor can then be expanded in expectation as

(3.4) $\begin{split}&\mathbb{E}_{t}\\!\left[\sum_{k=1}^{r}\\!u_{k}(t)\otimes
v_{k}(t)\otimes w_{k}(t)\\!\right]\\!\\!=\\!\\!\sum_{k=1}^{r}\\!a_{k}\otimes
b_{k}\otimes c_{k}+\mathbb{E}_{t}\\!\left[\sum_{k=1}^{r}\\!n_{1,k}(t)\otimes
n_{2,k}(t)\otimes c_{k}\\!\right]\\\
&\qquad\qquad\quad\,+\mathbb{E}_{t}\left[\sum_{k=1}^{r}n_{1,k}(t)\otimes
b_{k}\otimes n_{3,k}(t)\right]+\mathbb{E}_{t}\left[\sum_{k=1}^{r}a_{k}\otimes
n_{2,k}(t)\otimes n_{3,k}(t)\right]\\\ &\qquad\qquad\quad\,=T+E.\end{split}$

Here, we treat $E$ as a fully unknown tensor, with sparsity pattern determined
by the sparsity patterns of the matrices $Q_{ijk}$. Each slice of the first
term

(3.5) $\mathbb{E}_{t}\left[\sum_{k=1}^{r}n_{1,k}(t)\otimes n_{2,k}(t)\otimes
c_{k}\right]=\sum_{k=1}^{r}\mathbb{E}_{t}\left[n_{1,k}(t)\otimes
n_{2,k}(t)\right]\otimes c_{k}=\sum_{k=1}^{r}Q_{12k}\otimes c_{k},$

has (in the general case) non-zero elements in the positions given by
$\bigcup_{k=1}^{r}\text{supp}(Q_{12k})$. In other words, this first term
$\mathbb{E}_{t}\left[\sum_{k=1}^{r}n_{1,k}(t)\otimes n_{2,k}(t)\otimes
c_{k}\right]$ has support contained in

(3.6) $\left(\bigcup_{k=1}^{r}\text{supp}(Q_{12k})\times\llbracket
1,n\rrbracket\right).$

Similarly, the second two terms in (LABEL:eq:motivaton_unknown) have support
contained in permutations of the index-sets

(3.7) $\left(\bigcup_{k=1}^{r}\text{supp}(Q_{13k})\times\llbracket
1,n\rrbracket\right),\qquad\left(\bigcup_{k=1}^{r}\text{supp}(Q_{23k})\times\llbracket
1,n\rrbracket\right),$

respectively. Taken together, this specifies the pattern of the allowable non-
zero elements of $E$. The problem of recovering the tensor $T$ can therefore
be cast as a tensor completion problem with a deterministic sparsity pattern
encoded by the support of $E$.

In particular, this article will focus on the case of locally correlated noise
terms $n_{j,k}(t)$. By this, we mean terms for which the covariance terms
$Q_{ijk}$ are band-diagonal, for all $i$, $j$, and $k$. In detail, we will
specify a band-width $b$ and write the support of the covariance matrices as

(3.8) $\text{supp}(Q_{ijk})=\Big{\\{}(i_{1},i_{2})\in\llbracket
1,n\rrbracket^{2}:\lvert i_{1}-i_{2}\rvert\leq b\Big{\\}}.$

One can then verify that the resulting sparsity pattern for $E$ is of the form

(3.9) $\Big{\\{}(i_{1},i_{2},i_{3})\in\llbracket 1,n\rrbracket^{3}:\lvert
i_{k}-i_{\ell}\rvert\leq b,\text{ for some }k,\ell\in\\{1,2,3\\}\text{ with
}k\neq\ell\Big{\\}}.$

An illustration of this sparsity pattern is shown in figure 1.

Figure 1. Illustration of the sparsity pattern in (3.9). Left: first slice of
$E$. Middle: $\ell$th slice of $E$, for $b\leq\ell\leq n-b$. Right: $n$th
slice of $E$.

However, we would like to emphasize that the techniques of this article are
not by any means restricted to this specific sparsity pattern of $E$. Our
approach enables recovery of $T$ for a variety of different sampling patterns
encoded by $E$. In the three-dimensional case, two (non-exhaustive) examples
are illustrated in figure 2.

Figure 2. Two illustrations of sparsity patterns of $E$ for which $T$ in
(LABEL:eq:motivaton_unknown) can be recovered from knowledge only of the
corrupted tensor. Each slice of $E$ has the sparsity pattern indicated either
in $(1)$ or $(2)$.

###  3.2. Problem statement

Consider an unknown tensor $T=\sum_{k=1}^{r}\otimes_{j=1}^{d}v_{kj}$ of rank
$r$, which we aim to recover by knowledge of those elements of $T$ contained
in a specific sampling pattern. We will consider a fixed sampling pattern,
encoded by a tensor $E$. Although unknown, $E$ does, however, have a known
sparsity pattern. We therefore seek to reconstruct all elements of $T$ from
knowing only a corrupted version of $T$, denoted by $S$ and defined by

(3.10) $S=T+E.$

We will in the following focus on the sparsity pattern in (3.9), and in
particular, study three-dimensional $n\times n\times n$ tensors. We therefore
write $T=\sum_{k=1}^{r}a_{k}\otimes b_{k}\otimes c_{k}$.

Note that there clearly are examples of tensors $T$ that cannot be recovered
from knowledge of $S$ alone. One set of examples consists of tensors $T$ with
the same sparsity pattern as $E$. However, we will show that examples such as
these are degenerate cases, in that almost all tensors $T$ can be recovered,
under a few natural assumptions. Our results enable recovery under the
assumption

######  Assumption 1.

The two sets of vectors $\\{a_{k}\\}_{k=1}^{r}$ and $\\{b_{k}\\}_{k=1}^{r}$
both consist of pairwise orthogonal vectors.

This assumption is reasonable to impose, since general tensor decomposition
problems are NP-hard without similar assumptions [13]. Without loss of
generality, we will assume that the vectors $a_{k}$ and $b_{k}$ are unit
vectors, by absorbing their magnitudes into the vectors $c_{k}$.

Additionally, we will require a number of technical non-degeneracy conditions
during the course of our recovery algorithm. They hold true for almost all
collections of
$\\{a_{k}\\}_{k=1}^{r},\\{b_{k}\\}_{k=1}^{r},\\{c_{k}\\}_{k=1}^{r}$, so for
instance with probability one when these vectors are chosen at random. For
ease of exposition, we will therefore state our results in terms of generic
tensors. However, importantly, the degeneracy conditions can be explicitly
verified during a run of the algorithm. Should they not hold in a specific
use-case, this will be discovered and the user can remedy the deficiency by
generating more data, which corresponds to reducing the bandwidth $b$. Our
main result is the following.

######  Theorem 3.1.

Assume that Assumption 1 holds, and let $T$ be a tensor of the form in (3.10)
with unknown elements given by (3.9). For generic matrices $A,B,C$, there are
then constants $\alpha_{1},\alpha_{2},\beta_{1},\beta_{2}>0$ such that $T$ can
be uniquely recovered in polynomial time if

(3.11) $n\geq\alpha_{1}b+\beta_{1},$

and

(3.12) $r\geq\alpha_{2}b+\beta_{2}.$

An iterative recovery algorithm converges linearly to the true tensor, with
cost $O\left(\text{max}(n^{3},n^{2}r^{2},n^{2}b^{2}N)\right)$, where $N$ is
the number of required iterations.

The iterative recovery algorithm is presented in the next section. We would
again like to emphasize that the rank bound in (3.12) enforces a high-rank
assumption, which is principally the opposite of typical low-rank assumptions
for tensor and matrix completion.

The proof of the theorem is constructive, and upper bounds on the constants
$\alpha_{1},\beta_{1},\alpha_{2},\beta_{2}$ can be read out from the proof. We
also note lower bounds of the kind in theorem 3.1 are required for our
approach. In fact, we will show that tensors $T$ cannot be recovered with this
approach without similar bounds, which leads to corresponding lower bounds.
More precisely, the proof results in the bounds

(3.13)
$\begin{split}8(4b+2)&\geq\alpha_{1}b+\beta_{1}\geq\text{max}\\!\left(\\!12b+7,\frac{93b-69+\sqrt{1419b^{2}-2704b+921}}{10}\right),\\\
3(4b+2)&\geq\alpha_{2}b+\beta_{2}\geq 4b+1.\end{split}$

For large $b$, the lower bound on $n$ in (3.13) becomes

(3.14)
$n\geq\left(\frac{93}{10}+\frac{\sqrt{1419}}{10}\right)b-\frac{1352}{10\sqrt{1419}}-\frac{69}{10}+O\left(\frac{1}{b}\right)\approx
13.07b-10.49+O\left(\frac{1}{b}\right).$

However, the upper bounds in (3.13) seem pessimistic empirically, and
numerical experiments show that our techniques can in fact recover a sizable
percentage of unknown elements, as shown in table 1.

Table 1. The minimally admissible value of $n$ from (3.11) for different values of $b$, together with the percentage of unknown elements in the resulting tensor. $b$ | 1 | 3 | 5 | 7 | 10
---|---|---|---|---|---
Minimally admissible $n$ | 19 | 43 | 68 | 94 | 134
Unknown elements ($\%$) | 40.5 | 41.4 | 41.1 | 40.7 | 40.0
Minimally admissible $r$ | 7 | 15 | 23 | 30 | 42

##  4\. Algorithm

###  4.1. Overview

This section presents an algorithm to recover the tensor in theorem 3.1. The
steps in the algorithm are as follows:

  1. (1)

For a positive integer $m$, choose $m$ distinguished slices of $S$ (equation
(4.6)). Using these slices, construct a set of $\binom{m}{2}$ coupled
Sylvester-like equations (equation (4.10)).

  2. (2)

Solve the linear equations in the previous step, using an alternating least-
squares procedure to accelerate convergence. This recovers parts of the $m$
chosen slices.

  3. (3)

Construct a second set of coupled Sylvester-like equations (equation (4.20))
to recover the remaining entries of a subset of the $m$ distinguished slices.

  4. (4)

Use the reconstructed slices to recover the vectors $a_{k}$ and $b_{k}$ for
$k=1,\ldots,r$.

  5. (5)

Recover the vectors $c_{k}$ by solving one linear system per slice.

Sections 4.2.1-4.2.5 below present the details of each step. We first
introduce some notation. Write $A\in\mathbb{R}^{n\times r}$ for the matrix
obtained by stacking the $r$ vectors $a_{k}$, $k=1,\ldots,r$, side-by-side,
and similarly for $B\in\mathbb{R}^{n\times r}$ and $C\in\mathbb{R}^{n\times
r}$, i.e.,

(4.1) $\displaystyle A$ $\displaystyle=\begin{bmatrix}|&|&&|\\\
a_{1}&a_{2}&\ldots&a_{r}\\\
|&|&&|\end{bmatrix}\\!,\,\,\,B=\begin{bmatrix}|&|&&|\\\
b_{1}&b_{2}&\ldots&b_{r}\\\
|&|&&|\end{bmatrix}\\!,\,\,\,C=\begin{bmatrix}|&|&&|\\\
c_{1}&c_{2}&\ldots&c_{r}\\\ |&|&&|\end{bmatrix}.$

Define also a set of diagonal matrices $D^{c}_{\ell}\in\mathbb{R}^{r\times r}$
for $\ell=1,\ldots,n$ by

(4.2) $D^{c}_{\ell}=\begin{bmatrix}c_{1}(\ell)&\\\ &c_{2}(\ell)&\\\
&&\ddots&\\\ &&&c_{r}(\ell)\end{bmatrix}.$

We next define two index sets. Define $\mathcal{I}_{k}(b)$, for
$k=1,\ldots,n$, by

(4.3)
$\mathcal{I}_{k}(b)=\llbracket\text{max}(1,k-b),\text{min}(n,k+b)\rrbracket.$

Define also $\mathcal{D}(b)$ to be the indices contained in a band around the
diagonal of width $b$ in the $x$\- and $y$-directions, i.e.,

(4.4) $\mathcal{D}(b):=\Big{\\{}(i,j)\in\llbracket 1,n\rrbracket^{2}:\lvert
i-j\rvert\leq b\Big{\\}}.$

For each $k=1,\ldots,n$, we let $\mathcal{U}(k)$ denote the set of matrices
with support contained in the set of unknown elements of the $k$th slice of
$T$. More precisely, let

(4.5) $\mathcal{U}(k):=\Big{\\{}M\in\mathbb{R}^{n\times
n}:\text{supp}(M)\subseteq\mathcal{D}(b)\cup\Big{(}\mathcal{I}_{k}(b)\times\llbracket
1,n\rrbracket\Big{)}\cup\Big{(}\llbracket
1,n\rrbracket\times\mathcal{I}_{k}(b)\Big{)}\Big{\\}}.$

###  4.2. Details of algorithm

This section presents the details of the algorithm in section 4.1.

####  4.2.1. Step 1: recover parts of $m$ distinguished slices of $T$

We consider $m$ slices $k_{1},\ldots,k_{m}$ of $S$, chosen so that the
corresponding index sets $\mathcal{I}_{k_{1}}(b)$, $\mathcal{I}_{k_{2}}(b)$,
$\ldots$, $\mathcal{I}_{k_{m}}(b)$ are disjoint. Explicitly, we take

(4.6) $k_{i}=(2b+1)(i-1)+1,\quad\text{ for }i\in\llbracket 1,m\rrbracket.$

Note that this is possible when

(4.7) $(2b+1)(m-1)+1\leq n.$

We will see below that recovery is guaranteed for any $m\geq 5$. Nonetheless,
a higher value of $m$ provides more information to the algorithm and one
therefore expects this to lead to possible recovery for larger values of $b$,
compared to $n$ and $r$. However, the disjointness condition in (4.7)
introduces a bound on the highest possible value of $m$ in terms of $n$ and
$b$, meaning that increasing $m$ only increases recovery performance up to a
point.

From slices $k_{1},\ldots,k_{m}$ of $S$, we have access to the matrices

(4.8) $\displaystyle M_{i}$
$\displaystyle:=AD^{c}_{k_{i}}B^{T}+X_{i},\quad\text{ for }i\in\llbracket
1,m\rrbracket,$

where $X_{i}\in\mathcal{U}(k_{i})$, $i\in\llbracket 1,m\rrbracket$ are unknown
matrices. The first step of the algorithm recovers parts of the $X_{i}$ from a
set of Sylvester-type equations that we now construct.

For each pair $M_{i},M_{j}$ for $i,j\in\llbracket 1,m\rrbracket,i\neq j$, note
that

(4.9)
$\begin{split}(M_{i}-X_{i})\cdot(M_{j}-X_{j})^{T}-&(M_{j}-X_{j})\cdot(M_{i}-X_{i})^{T}\\\
&=AD^{c}_{k_{i}}D^{c}_{k_{j}}A^{T}-AD^{c}_{k_{j}}D^{c}_{k_{i}}A^{T}=0.\end{split}$

Expanding the left hand side therefore results in the set of equations

(4.10)
$\begin{split}X_{i}M_{j}^{T}-M_{j}X_{i}^{T}+M_{i}X_{j}^{T}&-X_{j}M_{i}^{T}+X_{j}X_{i}^{T}-X_{i}X_{j}^{T}\\\
&=M_{i}M_{j}^{T}-M_{j}M_{i}^{T},\quad i,j\in\llbracket
1,m\rrbracket.\end{split}$

We will next restrict these sets of equations to the complement of the support
of the quadratic terms $X_{j}X_{i}^{T}-X_{i}X_{j}^{T}$. The resulting linear
equations will then be used to determine parts of the $X_{i}$ uniquely. We
first write down the support of the quadratic terms, as follows.

######  Lemma 4.1.

Define the index set

(4.11)
$\begin{split}\mathcal{I}_{ij}(b):=\mathcal{D}(2b)&\cup\Big{(}\mathcal{I}_{i}(2b)\times\llbracket
1,n\rrbracket\Big{)}\cup\Big{(}\llbracket
1,n\rrbracket\times\mathcal{I}_{i}(2b)\Big{)}\\\
&\cup\Big{(}\mathcal{I}_{j}(2b)\times\llbracket
1,n\rrbracket\Big{)}\cup\Big{(}\llbracket
1,n\rrbracket\times\mathcal{I}_{j}(2b)\Big{)}.\end{split}$

For any two $N_{i}\in\mathcal{U}(i)$, $N_{j}\in\mathcal{U}(j)$ with
$\mathcal{I}_{i}(b)\cap\mathcal{I}_{j}(b)=\emptyset$, we have

(4.12)
$\text{supp}(N_{j}N_{i}^{T}-N_{i}N_{j}^{T})\subseteq\mathcal{I}_{ij}(b).$

###### Proof.

Write

(4.13) $\displaystyle
N_{i}=D_{i}+\sum_{k\in\mathcal{I}_{i}(b)}\left(x_{k}\otimes e_{k}+e_{k}\otimes
y_{k}\right),\,\,\,N_{j}=D_{j}+\sum_{\ell\in\mathcal{I}_{j}(b)}\left(z_{\ell}\otimes
e_{\ell}+e_{\ell}\otimes w_{\ell}\right),$

for vectors $x_{k},y_{k},z_{\ell},w_{\ell}$ in $\mathbb{R}^{n}$ and where
$D_{i},D_{j}$ have support contained in $\mathcal{D}(b)$. We start with the
term $N_{i}N_{j}^{T}$. In the expansion of the product
$N_{i}N_{j}^{T}-N_{j}N_{i}^{T}$, the cross terms $D_{i}D_{j}^{T}$ and
$D_{j}D_{i}^{T}$ have support contained in $\mathcal{D}(2b)$. The terms
$D_{i}\left(z_{\ell}\otimes e_{\ell}+e_{\ell}\otimes w_{\ell}\right)$ have
support in the rows in $\mathcal{I}_{j}(2b)$. Likewise, the terms
$(x_{k}\otimes e_{k})D_{j}^{T}$ have support in the columns with indices in
$\llbracket\text{min}(\mathcal{I}_{i}(b))-b,\text{max}(\mathcal{I}_{i}(b))+b\rrbracket$,
i.e., in $\mathcal{I}_{i}(2b)$. It follows that the terms $(x_{k}\otimes
e_{k})N_{j}^{T}$ have support in the columns with indices in
$\mathcal{I}_{i}(2b)\cup\mathcal{I}_{j}(b)$, since $\mathcal{I}_{i}(b)$ and
$\mathcal{I}_{j}(b)$ have empty intersection. Similarly, the terms
$(e_{k}\otimes y_{k})N_{j}^{T}$ have support contained in the rows with
indices in $\mathcal{I}_{i}(b)$.

For the term $N_{j}N_{i}^{T}$, transposing the support of the term
$N_{i}N_{j}^{T}$ and exchanging $i$ and $j$ concludes the proof. ∎

Restricting (4.10) to only the entries contained in the complement of
$\mathcal{I}_{k_{i}k_{j}}(b)$ therefore results in a linear equation in the
elements of $X_{i}$ and $X_{j}$. Note that only the rows of $X_{i}$ with
indices not contained in $\mathcal{I}_{k_{i}}(2b)$ appear in any of these
equations. We therefore use the collection of these equations for all
$i,j\in\llbracket 1,m\rrbracket$ to determine the rows of $X_{i}$ with indices
not contained in $\mathcal{I}_{k_{i}}(2b)$, as in the following lemma.

######  Lemma 4.2.

Let $X_{i}$ be matrices in $\mathcal{U}(k_{i})$, for $i=1,\ldots,m$ and define
$M_{i}$ in $\mathbb{R}^{n\times n}$ as in (4.8). Assume $m\geq 5$.

  1. (1)

For generically chosen matrices $A,B,C$ in $\mathbb{R}^{n\times r}$ satisfying
assumption 1, the equations

(4.14)
$\begin{split}X_{i}M_{j}^{T}-M_{j}X_{i}^{T}&+M_{i}X_{j}^{T}-X_{j}M_{i}^{T}\Big{|}_{\mathcal{I}_{k_{i}k_{j}}(b)^{\complement}}\\\
&=M_{i}M_{j}^{T}-M_{j}M_{i}^{T}\Big{|}_{\mathcal{I}_{k_{i}k_{j}}(b)^{\complement}},\quad
i,j\in\llbracket 1,m\rrbracket,\end{split}$

uniquely determine the rows of $X_{i}$ with indices not contained in
$\mathcal{I}_{k_{i}}(2b)$, provided

(4.15) $\begin{split}n&\geq 8(4b+2),\\\ r&\geq 3(4b+2).\end{split}$

  2. (2)

If either

(4.16)
$\begin{split}n&\leq\text{max}\left(12b+7,\frac{93}{10}b+\frac{\sqrt{1419b^{2}-2704b+921}}{10}-\frac{69}{10}\right)\quad\text{
or }\\\ r&\leq 4b+1,\end{split}$

then (4.14) has a non-trivial kernel when restricted to the rows of $X_{i}$
with indices not contained in $\mathcal{I}_{k_{i}}(2b)$.

A proof is given in appendix A. Note that the second part of lemma 4.2 shows
that high-rank assumptions of the type in (4.15) are indeed required in our
approach.

####  4.2.2. Step 2: solve the linear systems using alternating least-squares

The coupled linear systems in (4.14) constitute an overdetermined system of
$O(n^{2})$ equations in $O(nb)$ unknowns. A direct solution would therefore
have complexity $O(n^{4}b^{2})$. However, because of the particular structure
of the system, convergence can be accelerated through an alternating
minimization strategy. The least-squares solution of the linear system amounts
to minimizing the Euclidean norm of the residual. The proof of lemma 4.2 will
show that the associated linear system has full column rank, for generic
choices of $A,B,C$, so this cost function is strongly convex. The minimizer of
this problem can therefore be found by alternatingly minimizing the residual
with respect to a subset of the variables [18, Proposition 3.4]. Iterating
this procedure results in convergence to the solution of the linear system.

In detail, we alternate over each matrix $X_{i}$ for $i\in\llbracket
1,m\rrbracket$ in turn. Denote the $t$th iterate by $X_{i}^{(t)}$. For each
$i$ and $t$, we in turn update each row, $[X_{i}^{(t)}]_{\ell}$, of
$X_{i}^{(t)}$ for $\ell$ not contained in $\mathcal{I}_{k_{i}}(2b)$. Denote
the result of updating rows $\llbracket
1,\ell-1\rrbracket\smallsetminus\mathcal{I}_{k_{i}}(2b)$ of $X_{i}^{(t)}$ by
$X_{i}^{(t,\ell-1)}$. The update equation of the $\ell$th row reads as

(4.17)
$\begin{split}[X_{i}^{(t+1)}]_{\ell}=&\operatorname*{argmin\,}_{[X_{i}^{(t)}]_{\ell}}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{m}\Big{\|}X_{i}^{(t,\ell-1)}M_{j}^{T}-M_{j}X_{i}^{(t,\ell-1)T}\\\
&+M_{i}X_{j}^{(t,\ell-1)T}-X_{j}^{(t,\ell-1)}M_{i}^{T}-\left(M_{i}M_{j}^{T}-M_{j}M_{i}^{T}\right)\Big{|}_{\mathcal{I}_{k_{i}k_{j}}{(b)}^{\complement}}\Big{\|}^{2}.\end{split}$

The variables in $[X_{i}^{(t,\ell-1)}]_{\ell}$ only appear in the $\ell$th row
and column of the matrix

$X_{i}^{(t,\ell-1)}M_{j}^{T}-M_{j}X_{i}^{(t,\ell-1)T}+M_{i}X_{j}^{(t,\ell-1)T}-X_{j}^{(t,\ell-1)}M_{i}^{T}-\left(M_{i}M_{j}^{T}-M_{j}M_{i}^{T}\right)\Big{|}_{\mathcal{I}_{k_{i}k_{j}}{(b)}^{\complement}}.$

The update equation (4.17) is therefore a least squares problem in $O(b)$
unknowns and $O(n)$ equations. It can be solved with cost $O(nb^{2})$. The
total complexity of the alternating procedure therefore becomes
$O(n^{2}b^{2}N)$, where $N$ is the number of iterations used.

####  4.2.3. Step 3: recover remaining entries of a subset of the
distinguished slices of $T$

Next, we recover the remaining entries of $X_{i}$, for $i$ contained in a
subset $I$ of $\llbracket 1,m\rrbracket$. We will require that $I$ satisfies

(4.18) $\lvert I\rvert\geq 3$, (4.19) For all $i,j\in I$ with $i\neq j$,
$\mathcal{I}_{k_{i}}(2b)\cap\mathcal{I}_{k_{j}}(2b)=\emptyset$.

Note the multiplicative factor $2$ in the second requirement. As an example,
these requirements are necessarily satisfied for $I=\\{1,3,5\\}$, if $m\geq
5$. A larger value of $\lvert I\rvert$ will guarantee recovery under more
beneficial bounds of $n$ and $r$ in terms of $b$, but we will phrase our
results in terms of the limiting case $\lvert I\rvert=3$.

At this stage, the $i$th slice $M_{i}$ has all entries known except for
elements in the rows with indices contained in $\mathcal{I}_{k_{i}}(2b)$. This
is illustrated in figure 3.

Figure 3. Illustration of the remaining unknown elements in section 4.2.3.
Compare to figure 1. Left: first slice of $E$. Middle: $\ell$th slice of $E$,
for $b\leq\ell\leq n-b$. Right: $n$th slice of $E$.

We will recover these rows, for slices $k_{i}$ with $i$ contained in $I$. For
this purpose, we construct the system of equations

(4.20)
$(M_{i}-X_{i})^{T}\cdot(M_{j}-X_{j})-(M_{j}-X_{j})^{T}\cdot(M_{i}-X_{i})=0,\quad
i,j\in I,$

where each $X_{i}$ is contained in the subspace

(4.21)
$\Big{\\{}X_{i}\in\mathcal{U}(k_{i}):\text{supp}(X_{i})\subseteq\mathcal{I}_{k_{i}}(2b)\times\llbracket
1,n\rrbracket\Big{\\}}.$

Note that the order of the terms within each product in (4.20) is reversed, as
compared to (4.9). Moreover, the quadratic cross terms in (4.20) satisfy

(4.22) $\displaystyle X_{i}^{T}X_{j}=0=X_{j}^{T}X_{i},\quad i,j\in I,i\neq j,$

since the sets $\mathcal{I}_{k_{i}}(2b)$ and $\mathcal{I}_{k_{j}}(2b)$ are
disjoint by the second assumption on the set $I$. Equation (4.20) is therefore
a set of coupled linear systems for the unknowns $X_{i}$. We will show that
they determine the $X_{i}$ uniquely.

We next present the main result of this section.

######  Lemma 4.3.

For generic $A$ and $B$ and $I$ satisfying (4.18)–(4.19), equation (4.20) has
a unique solution, provided that

(4.23) $\begin{split}n&\geq 16b+4\\\ r&\geq 8b+4\end{split}$

The bound on $r$ in (4.23) is slightly pessimistic and we observe a looser
bound in the numerical experiments below. Equation (4.20) can be solved in an
alternating fashion similar to (4.17) and we omit the details of this step.

####  4.2.4. Step 4: use the recovered slices to recover $a_{k}$, $b_{k}$

From the $\lvert I\rvert$ distinguished slices recovered in the previous step,
we can use Jennrich’s algorithm [12] to recover the vectors $a_{k}$, $b_{k}$.
In detail, write $\widehat{T}\in\mathbb{R}^{n\times n\times\lvert I\rvert}$ as
the tensor obtained by stacking the $\lvert I\rvert$ recovered slices
$k_{i},i\in I$ of $T$, i.e.,

(4.24) $\widehat{T}(i_{1},i_{2},i_{3})=T(i_{1},i_{2},s_{i_{3}}),$

where the elements of $I$ are enumerated as $s_{1},\ldots,s_{\lvert I\rvert}$.
Writing $\widehat{c}_{\ell}$ as the restriction of $c_{\ell}$ to the elements
with indices in $I$, $\widehat{T}$ has the decomposition

(4.25) $\widehat{T}=\sum_{\ell=1}^{r}a_{\ell}\otimes
b_{\ell}\otimes\widehat{c}_{\ell},$

and the vectors $a_{k}$, $b_{k}$ for $k=1,\ldots,r$ can be recovered by an
application of Jennrich’s algorithm.

####  4.2.5. Step 5: recover the $c_{k}$

Lastly, we use the $\ell$th slice of $S$, for $\ell=1,\ldots,n$, to recover
the elements $c_{1}(\ell),\ldots,c_{r}(\ell)$. From each slice, we have access
to

(4.26) $M_{\ell}=AD_{\ell}^{c}B^{T}+N_{\ell},$

for $N_{\ell}$ an unknown matrix in $\mathcal{U}(\ell)$. Since $A$ and $B$ are
known from the preceding step, restricting (4.26) to the complement of the
support of $\mathcal{U}(\ell)$ results in a linear equation in
$c_{1}(\ell),\ldots,c_{r}(\ell)$. We therefore construct the linear system

(4.27)
$M_{\ell}(i,j)=\sum_{k=1}^{r}A(i,k)B(j,k)c_{k}(\ell),\quad(i,j)\in\text{supp}(\mathcal{U}(\ell))^{\complement}.$

######  Lemma 4.4.

For generic $A$ and $B$, the system in (4.27) has a unique solution, provided
that $n\geq 8b+4$.

###### Proof.

Let $i_{0}=\text{max}\left(\mathcal{I}_{\ell}(b)\right)+1$, if $\ell<n-b$, and
$i_{0}=\text{min}\left(\mathcal{I}_{n}(b)\right)-1$ if $\ell\geq n-b$.
Restrict the system (4.27) to the indices $(i_{0},j_{1})$ for
$(i_{0},j_{1})\in\text{supp}(\mathcal{U}(\ell))^{\complement}$ and
$(j_{2},i_{0})$ where
$(j_{2},i_{0})\in\text{supp}(\mathcal{U}(\ell))^{\complement}$. If we
enumerate the possible values of $j_{1}$ as $s_{1},\ldots,s_{N}$ and the
values of $j_{2}$ as $t_{1},\ldots,t_{M}$, the restricted system can be
written in the form

(4.28)

$\begin{bmatrix}A(i_{0},1)B(s_{1},1)&A(i_{0},2)B(s_{1},2)&\ldots&A(i_{0},r)B(s_{1},r)\\\
\vdots&\vdots&\ddots&\vdots\\\
A(i_{0},1)B(s_{N},1)&A(i_{0},2)B(s_{N},2)&\ldots&A(i_{0},r)B(s_{N},r)\\\
A(t_{1},1)B(i_{0},1)&A(t_{1},2)B(i_{0},2)&\ldots&A(t_{1},r)B(i_{0},r)\\\
\vdots&\\\
A(t_{M},1)B(i_{0},1)&A(t_{M},2)B(i_{0},2)&\ldots&A(t_{M},r)B(i_{0},r)\end{bmatrix}\begin{bmatrix}c_{1}(\ell)\\\
c_{2}(\ell)\\\ \vdots\\\
c_{r}(\ell)\end{bmatrix}=\begin{bmatrix}M_{\ell}(i_{0},s_{1})\\\
M_{\ell}(i_{0},s_{2})\\\ \vdots\\\ M_{\ell}(i_{0},s_{N})\\\
M_{\ell}(t_{1},i_{0})\\\ M_{\ell}(t_{2},i_{0})\\\ \vdots\\\
M_{\ell}(t_{M},i_{0})\end{bmatrix}.$

Note that $i_{0}$ is not contained in either of the lists $s_{1},\ldots,s_{N}$
or $t_{1},\ldots,t_{M}$ so the matrix in the left hand side has full column
rank for generic $A$ and $B$ if $M+N\geq r$. Since $M+N\geq 2(n-2(2b+1))$ and
$n\geq r$, the conclusion follows if we impose $2(n-2(2b+1))\geq n$. This is
equivalent to $n\geq 8b+4$, which concludes the proof. ∎

The system in (4.27) therefore determines $c_{1}(\ell),\ldots,c_{r}(\ell)$
uniquely. Repeating this procedure for $\ell=1,\ldots,n$ therefore recovers
the vectors $c_{1},\ldots,c_{r}$.

##  5\. Numerical results

This section details different numerical results of the algorithm of the
article. The tensors used in the simulations are all of the form
$\sum_{k=1}^{r}a_{k}\otimes b_{k}\otimes c_{k}$, with $c_{k}$ having
independently generated normal entries, and the $a_{k},b_{k}$ randomly
generated orthogonal vectors. All computations were carried out on a MacBook
Pro with a 3.1 GHz Intel Core i5 processor and 16 GB of memory.

As stopping criterion for the iterative solutions of equations (4.17) and
(4.20), we terminate the iterations once the relative improvement from one
iteration to the next is below a threshold denoted by
$\varepsilon_{\text{tol}}$.

Figure 4. Numerical results for the experiments in sections 5.1, 5.2 and 5.3.
respectively. (a) Successfully recovered tensors as function of $b$ and the
number of slices used. Green color indicates recovery and red non-recovery.
(b) Number of iterations required for convergence of (4.17) as a function of
$b$ and $r$. (c) Effect of entrywise noise on recovery error, as function of
the relative perturbation magnitude $\frac{\|E_{\text{noise}}\|}{\|T\|}$ and
$r$.

###  5.1. Effect of number of slices on highest possible $b$

For a fixed value $n=100=r$, we show the effect of recovery by the value of
$b$ and the number of slices used. We used $\varepsilon_{\text{tol}}=10^{-6}$.
The results are shown in figure 4(a). Adding additional slices increases the
maximum possible recoverable $b$ only up to a certain point, since our
algorithm requires the slices to have disjoint index sets
$\mathcal{I}_{k_{i}}(b)$ for all slices $i=1,\ldots,m$ which imposes the bound
on $b$ in terms of $m$ and $n$ in (4.7). The use of $m=7$ slices is clearly
optimal in this case, and the maximum value of recoverable $b$ coincides with
that of the lower bound on $n$ in (3.13).

###  5.2. Number of iterations required for convergence as function $b$

For a fixed value $n=100$, we show the number of iterations needed for
convergence of the alternating least squares procedure in (4.17). For each
value of $b$, we run $10$ randomly generated trials, and report the average
number of iterations required for convergence. We used
$\varepsilon_{\text{tol}}=10^{-7}$. The results are shown in figure 4(b). Note
that the number of required iterations remains relatively low for $b$ below
the maximum admissible $b$ for this value of $n$. However, the number of
required iterations increases for $b$ approaching its maximally recoverable
value, and for decreasing $r$.

###  5.3. Effect of entrywise noise on recovery error

For a fixed value $n=65$, we study the effect of an additional entrywise noise
term on the recovery error. In detail, we generate random tensors $T$ of the
form $\sum_{k=1}^{r}a_{k}\otimes b_{k}\otimes c_{k}$ with $a_{k},b_{k}$
orthogonal vectors for $k=1,\ldots,r$, and apply the algorithm of the article
to the tensor $T+E+E_{\text{noise}}$, where $E$ is the structured noise tensor
discussed in the algorithm with $b=5$ and $E_{\text{noise}}$ is a tensor with
independent normal entries. The algorithm of the article is then run on
$T+E+E_{\text{noise}}$ with resulting tensor $T_{0}$. For each $b$, we run
$100$ randomly generated trials with $\varepsilon_{\text{tol}}=10^{-12}$. We
show the resulting average recovery error $\frac{\|T-{T_{0}}\|}{\|T\|}$ as a
function of $\frac{\|E_{\text{noise}}\|}{\|T\|}$ and $r$ in figure 4(c).

##  6\. Conclusion

We presented an algorithm to recover orthogonally decomposable tensors
corrupted by arbitrarily strong, but structured noise. The problem can be seen
as a deterministic tensor completion problem, and we have shown how this can
be solved provided the tensor dimension and rank are larger than an affine
function of the corruption bandwidth. Notably, this enables recovery under a
high-rank assumption on the tensor, as opposed to low-rank assumptions
commonly required for completion problems. The techniques in the article are
not limited to the specific pattern of unknown elements treated. Future work
therefore includes studying the possibility of bridging the gap between high-
and low-rank completion techniques to develop algorithms with weaker
conditions on the tensor rank. Another avenue for future work lies in relaxing
the orthogonality condition in assumption 1 to soft-orthogonality constraints
encoded instead by incoherent tensor components.

##  Appendix A Proof of lemma 4.2

We start with the first part of lemma 4.2. For each slice $X_{i}$, only the
variables not contained in the rows with indices in $\mathcal{I}_{k_{i}}(2b)$
enter into any of the equations in (4.14). We therefore need to show that the
linear operator in (4.14) has trivial kernel under the assumptions in (4.15),
for generic matrices $A,B,C$. This amounts to ensuring that the minors of the
operator in (4.14) are non-zero. Since these minors are polynomial expressions
in the elements of $A,B,C$, this will be true generically, provided that we
can show that these polynomials are not identically vanishing. This will hold
if we can exhibit one example of the operator in (4.14) with full column rank,
under the bounds in (4.15). We first prove an auxiliary result.

######  Lemma A.1.

Let $A_{1},A_{2},A_{3}$ be generic matrices in $\mathbb{R}^{m\times n}$. Let
$I_{1},I_{2},I_{3}$ be disjoint subsets of $\llbracket 1,n\rrbracket$ and let
$Y_{1},Y_{2},Y_{3}$ be matrices in $\mathbb{R}^{n\times p}$. If

  1. (1)

$Y_{i}$ has non-zero elements only in the rows with indices in the set
$I_{i}$, for $i=1,2,3$

  2. (2)

each column of $Y_{i}$ has at most $m$ non-zero elements, for $i=1,2,3$

then the system

(A.1) $\displaystyle A_{2}Y_{1}$ $\displaystyle=-(A_{1}Y_{2})^{T}$ (A.2)
$\displaystyle A_{3}Y_{1}$ $\displaystyle=-(A_{1}Y_{3})^{T}$ (A.3)
$\displaystyle A_{3}Y_{2}$ $\displaystyle=-(A_{2}Y_{3})^{T}$

has the unique solution $Y_{1}=Y_{2}=Y_{3}=0$.

###### Proof.

Just as in the opening paragraph of this section, we need only find one
example of matrices $A_{1},A_{2},A_{3}$ for which the conclusion holds, to
ensure that it holds generically. To do this, let the columns of $A_{1}$ with
indices in $I_{2}$ be zero. It follows that $A_{1}Y_{2}=0$, so, from (A.1),
$A_{2}Y_{1}=0$. In the first column of this equation, discarding the zero
elements of the first column of $Y_{1}$ gives a matrix equation for the non-
zero elements in the first column of $Y_{1}$. From the second assumption in
the lemma, the corresponding matrix has full column rank generically. It
follows that the first column of $Y_{1}$ is zero, and similarly the remaining
columns are as well. Inserting this into (A.2) similarly gives $Y_{3}=0$,
since $I_{2}$ and $I_{3}$ are disjoint. Inserted into (A.3), this gives also
$Y_{2}=0$, which concludes the proof. ∎

We now proceed to construct one example of the operator in (4.14) with full
column rank. For any choice of three slices, we will show that the resulting
variables $X_{i}$ are zero. which will suffice to prove the lemma. We first
consider slices $i=1,3,5$. Let therefore $X_{1},X_{3},X_{5}$ be in the kernel
of (4.14). These matrices then satisfy

(A.4)
$X_{i}M_{j}^{T}-M_{j}X_{i}^{T}+M_{i}X_{j}^{T}-X_{j}M_{i}^{T}\Big{|}_{\mathcal{I}_{k_{i}k_{j}}(b)^{\complement}}=0,\quad
i,j\in\\{1,3,5\\}.$

We will construct $M_{1},M_{3},M_{5}$ with support contained in disjoint rows.
In detail, let $J_{1},J_{3},J_{5}$ be disjoint subsets of $\llbracket
1,n\rrbracket$ such that

(A.5) $\displaystyle J_{i}\cap\mathcal{I}_{\ell}(2b)=\emptyset,\text{ for
}i,\ell\in\\{1,3,5\\},$ (A.6) $\displaystyle\lvert J_{i}\rvert=4b+2,\text{ for
}i\in\\{1,3,5\\},$ (A.7)
$\displaystyle\text{min}(J_{3})\geq\text{max}(J_{1})+4b+2,$ (A.8)
$\displaystyle\text{min}(J_{5})\geq\text{max}(J_{3})+4b+2.$

Let $M_{i}$ have support in the rows contained in $J_{i}$. Note that this is
possible when $r\geq\lvert J_{1}\rvert+\lvert J_{2}\rvert+\lvert
J_{3}\rvert=3(4b+2)$ and $n\geq 8(4b+2)$, by letting the first $\lvert
J_{1}\rvert$ columns of $A$ have support contained in the rows with indices in
$J_{1}$, the subsequent $\lvert J_{3}\rvert$ columns of $A$ have support
contained in the rows with indices in $J_{3}$ and the last $\lvert
J_{5}\rvert$ columns of $A$ have support in rows in $J_{5}$. If
$D_{k_{1}}^{c}$ has support in the first $\lvert J_{1}\rvert$ diagonal
elements, $D_{k_{3}}^{c}$ in the subsequent $\lvert J_{3}\rvert$ diagonal
entries and $D_{k_{5}}^{c}$ in the subsequent $\lvert J_{5}\rvert$ ones, the
$M_{i}$ have the desired support.

We now show that this choice of $M_{i}$ enforces $X_{1}=X_{3}=X_{5}=0$. The
matrix

(A.9) $M_{i}X_{j}^{T}-M_{j}X_{i}^{T},$

has support contained in the rows with indices in $J_{i}\cup J_{j}$. From
(A.4), it follows that the entries of $M_{i}X_{j}^{T}-M_{j}X_{i}^{T}$ not
contained in $J_{i}\times J_{i}$, $J_{i}\times J_{j}$, $J_{j}\times J_{i}$,
$J_{j}\times J_{j}$ or $\mathcal{I}_{k_{i}k_{j}}(b)$ are zero. Therefore,
elements with rows in $J_{j}$ and columns not contained in

(A.10) $\llbracket\text{min}(J_{j})-2b,\text{max}(J_{j})+2b\rrbracket\cup
J_{i}\cup\mathcal{I}_{k_{i}}(2b)\cup\mathcal{I}_{k_{j}}(2b),$

are zero. We now claim that also the columns of $X_{i}^{T}$ not contained in
(A.10) are zero. To see this, observe that the non-zero elements $x$ of any
such column appears in (A.9) as an equation of the form $Mx=0$, where $M$ is a
generic matrix with $\lvert J_{j}\rvert$ rows. Since each column of
$X_{i}^{T}$ has at most $4b+2$ non-zero elements, (A.6) enforces $x=0$.

Next, intersect the sets in (A.10) for all $j\neq i$. Equations (A.7)–(A.8)
show that the columns of $X_{i}^{T}$ not contained in
$J_{i}\cup\mathcal{I}_{k_{i}}(2b)$ are zero, for each $i=1,3,5$. We lastly
show that the remaining columns of $X_{i}^{T}$ with indices in $J_{i}$ are
zero as well. Equation (A.4) now reads as

(A.11) $\displaystyle M_{3}(J_{3},:)X_{1}^{T}(:,J_{1})$
$\displaystyle=-\Big{(}M_{1}(J_{1},:)X_{3}^{T}(:,J_{3})\Big{)}^{T},$ (A.12)
$\displaystyle M_{5}(J_{5},:)X_{1}^{T}(:,J_{1})$
$\displaystyle=-\Big{(}M_{1}(J_{1},:)X_{5}^{T}(:,J_{5})\Big{)}^{T},$ (A.13)
$\displaystyle M_{5}(J_{5},:)X_{3}^{T}(:,J_{3})$
$\displaystyle=-\Big{(}M_{3}(J_{3},:)X_{5}^{T}(:,J_{5})\Big{)}^{T}.$

By (A.5)–(A.8), this system of equations satisfies the conditions in lemma
A.1. It follows that $X_{1}=X_{3}=X_{5}=0$ for generic matrices $A,B,C$
satisfying assumption 1.

Next, for any remaining slice $X_{i}$, for $i\neq 1,3,5$, repeat the above
argument with slices $1,3,i$ to conclude that also $X_{i}=0$. This concludes
the proof of the first part of lemma 4.2.

We next prove the second part of lemma 4.2 and start with the bound on $n$ in
(4.16). This comes from counting the number of equations present in (4.14) and
enforcing that this equals at least the number of unknowns. Write

(A.14) $\displaystyle n_{i}$
$\displaystyle=\text{dim}\left(\mathcal{U}(k_{i})\Big{|}_{\mathcal{I}_{k_{i}}(2b)^{\complement}\times\llbracket
1,n\rrbracket}\right).$

We first count the number of unknowns, i.e., $\sum_{i=1}^{m}n_{i}$. One can
verify that

(A.15) $\displaystyle n_{1}$
$\displaystyle=(n-2(b+1))(3b+2)-\frac{b(b+1)}{2},$ (A.16) $\displaystyle
n_{2}$ $\displaystyle=(2b+1)(2n-8b-3)+b+1-\frac{b(b+1)}{2},$ (A.17)
$\displaystyle n_{\ell}$
$\displaystyle=(2b+1)(2n-8b-2)-b(b+1),\quad\ell\not\in\\{1,2,m\\}.$

For $n_{m}$, we write $\Delta=n-m\cdot(2b+1)$ and distinguish the two cases
$\Delta\geq 0$ and $\Delta<0$. We have

(A.18) $\displaystyle n_{m}$
$\displaystyle=\begin{cases}(2b+1)(2n-8b-2)-b(b+1),\quad\text{ if }\Delta\geq
0,\\\
(n-\Delta_{r})(2b+1)+\Delta_{c}(n-\Delta_{r})-\frac{b(b+1)}{2},\quad\text{ if
}\Delta<0,\end{cases}$

where

(A.19) $\displaystyle\Delta_{r}$ $\displaystyle=n-(m-2)(2b+1)-2,$ (A.20)
$\displaystyle\Delta_{c}$ $\displaystyle=\begin{cases}2b+1,\quad\text{ if
}(m-1)(2b+1)+b+1\leq n,\\\ n-((m-2)(2b+1)+b+1),\text{ if }(m-1)(2b+1)+b+1\geq
n.\end{cases}$

Next, we count the number of equations present in (4.14). Figure 5 illustrates
the sparsity patterns of the equations retained in (4.14).

Figure 5. Illustration of the equations retained in (4.14). Orange color
indicates that the equation in the corresponding element is retained, and
white color that it is not retained. The top row shows the entire matrix and
the bottom row removes the white space induced by the deleted rows and
columns. $(a),(b)$ and $(c)$ illustrate $n_{58}$, $n_{68}$, and $n_{78}$,
respectively.

One can verify that the number of equations present in (4.14) for a specific
choice of $i,j$ is

(A.21) $\frac{(\widehat{n}-2b)(\widehat{n}-2b+1)}{2}+n_{ij},$

where
$\widehat{n}=n-\lvert\mathcal{I}_{k_{i}}(2b)\cup\mathcal{I}_{k_{j}}(2b)\rvert$
and the $n_{ij}$ are real numbers with $n_{ij}=n_{ji}$. With $i<j$, the
$n_{ij}$ are determined by

Case 1: $\Delta>2b$. If $i,j\not\in\\{1,2\\}$, then

(A.22) $n_{ij}=\begin{cases}b(2b+1),\quad\lvert i-j\rvert=1,\\\
b(2b+1)+2b,\quad\lvert i-j\rvert=2,\\\ 2b(2b+1),\quad\lvert i-j\rvert\geq
3.\end{cases}$

If $i=1$, then

(A.23) $n_{1j}=\begin{cases}0,\quad j=2,\\\ 2b,\quad j=3,\\\ b(2b+1),\quad
j\geq 4.\end{cases}$

If $i=2$, then

(A.24) $n_{2j}=\begin{cases}2b,\quad j=3,\\\ 4b,\quad j=4,\\\ 2b+b(2b+1),\quad
j\geq 5.\end{cases}$

Case 2: $\Delta\leq 2b$. Same as case 1, except subtract
$\frac{(2b-\max(0,\Delta))(2b-\max(0,\Delta)+1)}{2}$ from $n_{im}$.

Case 3: $\delta=n-(m-1)\cdot(2b+1)\leq 2b$. Same as case 2, except subtract
$\frac{(2b-\delta)(2b-\delta+1)}{2}$ from $n_{i,m-1}$.

Ensuring that the number of unknowns is at most the number of available
equations with these expressions results in a lower bound for admissible $n$
in terms of $b$. This bound is shown in figure 6 for a few different values of
$m$. The figure also shows the line $L(m)$ determined by $n=(2b+1)(m-1)+1$,
which is the additional lower bound from (4.7), since the $m$ chosen slices
$k_{1},\ldots,k_{m}$ were required to be disjoint.

Figure 6. Lower bounds of admissible $n$ as function of $b$.

Since the lower bound from using seven slices is contained above $L(7)$, but
beneath $L(8)$, it follows that the optimal choice of $m$ is $m=7$. For this
choice, explicitly writing out the bound results in

(A.25)
$n\geq\text{max}\left(12b+7,\frac{93}{10}b+\frac{\sqrt{1419b^{2}-2704b+921}}{10}-\frac{69}{10}\right),$

which is precisely the bound on $n$ in (4.16).

We lastly treat the bound on $r$ in (4.16). If $r<4b+2$, study a column of
$X_{2}$ with $4b+2$ non-zero elements. For each row of $B^{T}$, pick out the
columns with same sparsity pattern as the chosen column of $X_{2}$. The
resulting matrix necessarily has linearly dependent columns. Take therefore
$v$ in $\mathbb{R}^{4b+2}$ as a non-zero vector in the kernel of this matrix
and define $X_{2}$ by inserting $v$ into the non-zero elements of the chosen
column. Let all other entries of $X_{2}$ and $X_{i}$ for $i\neq 2$ be zero.
Clearly, $X_{i}\in\mathcal{U}(k_{i})$, for $i\in\llbracket 1,m\rrbracket$ and
the $X_{i}$ are non-trivial elements of the kernel of (4.14). This concludes
the proof of lemma 4.2.

##  Appendix B Proof of lemma 4.3

It is possible to prove lemma 4.3 in a similar fashion to the proof of lemma
4.2. However, we present a different type of argument that allows for tighter
bounds on $n$ and $r$ in terms of $b$. For ease of reference, we first prove
the following two auxiliary results.

######  Lemma B.1.

Let $M$ in $\mathbb{R}^{p\times q}$ have linearly independent columns and let
$X$ be a matrix in $\mathbb{R}^{q\times p}$. The statements

  1. (1)

$MX$ is a symmetric $p\times p$-matrix

  2. (2)

$X^{T}=MS$, for $S$ a symmetric $q\times q$-matrix

are then equivalent.

###### Proof.

It is clear that the second statement implies the first. For the reverse
direction, complete the column vectors of $M$ to a basis $m_{1},\ldots,m_{p}$,
and let $n_{1},\ldots,n_{p}$ be a dual basis. Since $MX$ is symmetric, it can
be written in the form

(B.1) $MX=\sum_{i,j=1}^{p}\alpha_{ij}m_{i}\otimes m_{j},$

where $\alpha_{ij}=\alpha_{ji}$. For a fixed $\ell=q+1,\ldots,p$, acting on
this equation with $n_{\ell}$ from the left yields

(B.2) $0=n_{\ell}^{T}MX=\sum_{j=1}^{p}\alpha_{\ell j}m_{j}.$

By linear independence of the $m_{j}$, it follows that $\alpha_{\ell
j}=0=\alpha_{j\ell}$, so we can write $MX=MSM^{T}$, where $S$ is the symmetric
$q\times q$-matrix defined by $S_{ij}=\alpha_{ij}$. Acting on this equation by
the pseudoinverse $M^{\dagger}$ from the left results in $X=SM^{T}$, so
$X^{T}=MS$, which concludes the proof. ∎

######  Lemma B.2.

Let $I,J,K\subseteq\llbracket 1,n\rrbracket$ be index sets and
$A,B,C\in\mathbb{R}^{n\times r}$ generic matrices satisfying assumption 1. Let
$i,j,k$ be distinct integers in $\llbracket 1,n\rrbracket$. Write

(B.3) $M_{i}=BD_{i}^{c}A^{T},\quad M_{j}=BD_{j}^{c}A^{T},\quad
M_{k}=BD_{k}^{c}A^{T}.$

Then

  1. (1)

The matrix

(B.4) $\begin{bmatrix}M_{j}(I,J),&-M_{i}(I,K)\end{bmatrix}$

has full column rank if $J$ and $K$ are disjoint, $\lvert I\rvert\geq\lvert
J\rvert+\lvert K\rvert$, and $r\geq\lvert J\rvert+\lvert K\rvert$.

  2. (2)

The matrix

(B.5)
$\begin{bmatrix}M_{i}(:,J),&-M_{j}(:,I),&M_{i}(:,K),&-M_{k}(:,I)\end{bmatrix}$

has full column rank if $I,J$ and $K$ are disjoint, and $r\geq 2\lvert
I\rvert+\lvert J\rvert+\lvert K\rvert$.

###### Proof.

Just as in the first paragraph of appendix A, we need only produce examples of
the matrices in statements $1$ and $2$ with full rank, under the given
assumptions.

For statement $1$, choose $D_{i}^{c}$ to have non-zero diagonal elements for
the indices contained in $J$, and $D_{j}^{c}$ in $K$, which is possible since
$r\geq\lvert J\rvert+\lvert K\rvert$. Since $J$ and $K$ were assumed disjoint
and $\lvert I\rvert\geq\lvert J\rvert+\lvert K\rvert$, the entries of $A$ and
$B$ with row indices in $I$ and column indices in $J$ and $K$ can be chosen
separately, so that the matrix in (B.4) equals e.g., the $\lvert
I\rvert\times(\lvert J\rvert+\lvert K\rvert)$ identity matrix. This has full
column rank, which concludes the proof of the first statement.

For statement $2$, let $D_{i}^{c}$, $D_{j}^{c}$, $D_{k}^{c}$ have non-zero
elements only in the diagonal elements with indices in $I,J,K$, respectively.
Since $I,J,K$ are disjoint, we can choose the columns of $A,B,C$ contained in
$I,J,K$ separately to ensure the matrix in (B.5) equals the $n\times(2\lvert
I\rvert+\lvert J\rvert+\lvert K\rvert)$ identity matrix, with full column
rank. ∎

We can now present the proof of lemma 4.3.

###### Proof of lemma 4.3.

We show that the kernel of (4.20) is trivial, so take $X_{i}$ from this
kernel, for $i\in I$. We first show that the rows of each $X_{i}$ contained in
$\mathcal{I}_{k_{i}}(2b)\smallsetminus\mathcal{I}_{k_{i}}(b)$ are zero. To see
this, we can express the fact that the $X_{i}$ are in the kernel of
$\eqref{eq:solve_N2nd}$ by

(B.6) $M_{i}^{T}X_{j}-M_{j}^{T}X_{i}=S_{ij},\quad i,j\in I,i\neq j,$

where each $S_{ij}$ is a symmetric $n\times n$-matrix. Fix a pair $i,j$ in $I$
with $\mathcal{I}_{k_{j}}(3b)\cap\mathcal{I}_{k_{i}}(3b)=\emptyset$, which is
possible under the assumptions on $I$ in (4.18)–(4.19). We restrict (B.6) to
rows and columns with indices not contained in $J:=\mathcal{I}_{k_{j}}(3b)$.
The columns of $X_{i}$ and $X_{j}$ not contained in $J$ have non-zero elements
in rows with indices $\mathcal{I}_{k_{i}}(2b)$ and $\mathcal{I}_{k_{j}}(b)$,
respectively. Taking this zero structure into account, (B.6) means that

(B.7)
$\left[M_{j}^{T}(J^{\complement},\mathcal{I}_{k_{i}}(2b)),-M_{i}^{T}(J^{\complement},\mathcal{I}_{k_{j}}(b))\right]\begin{bmatrix}X_{i}(\mathcal{I}_{k_{i}}(2b),J^{\complement})\\\
X_{j}(\mathcal{I}_{k_{j}}(b),J^{\complement})\end{bmatrix},$

is a symmetric matrix. By lemma B.2, the left matrix has full column rank
provided $\lvert J^{\complement}\rvert=n-(6b+1)\geq 6b+2$, i.e., $n\geq 12b+3$
and $r\geq 6b+2$. By lemma B.1, this means that

(B.8) $\begin{bmatrix}X_{i}(\mathcal{I}_{k_{i}}(2b),J^{\complement})\\\
X_{j}(\mathcal{I}_{k_{j}}(b),J^{\complement})\end{bmatrix}^{T}=\left[M_{j}^{T}(J^{\complement},\mathcal{I}_{k_{i}}(2b)),-M_{i}^{T}(J^{\complement},\mathcal{I}_{k_{j}}(b))\right]S,$

for some symmetric matrix $S$.

Next, the columns of $X_{i}(\mathcal{I}_{k_{i}}(2b),J^{\complement})^{T}$ with
indices in $\mathcal{I}_{k_{i}}(2b)\smallsetminus\mathcal{I}_{k_{i}}(b)$ have
by construction at least $n-(4b+1)-\lvert J\rvert=n-(10b+2)$ zero elements.
The corresponding rows of

(B.9)
$\left[M_{j}^{T}(J^{\complement},\mathcal{I}_{k_{i}}(2b)),-M_{i}^{T}(J^{\complement},\mathcal{I}_{k_{j}}(b))\right],$

are linearly independent by lemma B.2, provided it has at least $n-(10b+2)$
columns, i.e., provided $n\geq 16b+4$, and $r\geq 6b+2$. It follows that the
columns of $S$ with indices in
$\mathcal{I}_{k_{i}}(2b)\smallsetminus\mathcal{I}_{k_{i}}(b)$ are zero. The
same columns of $X_{i}^{T}$ are therefore zero as well.

We next show that the remaining elements of the $X_{i}$ are zero. Equation
(B.6) now says that

(B.10)
$\left[M_{i}^{T}(:,\mathcal{I}_{k_{j}}(b)),-M_{j}^{T}(\mathcal{I}_{k_{i}}(b))\right]\begin{bmatrix}X_{j}(\mathcal{I}_{k_{j}}(b),:)\\\
X_{i}(\mathcal{I}_{k_{i}}(b),:)\end{bmatrix}$

is a symmetric matrix. Lemma B.1 shows that

(B.11)
$\begin{bmatrix}X_{j}(\mathcal{I}_{k_{j}}(b),:)^{T},X_{i}(\mathcal{I}_{k_{i}}(b),:)^{T}\end{bmatrix}=\left[M_{i}^{T}(:,\mathcal{I}_{k_{j}}(b)),-M_{j}^{T}(:,\mathcal{I}_{k_{i}}(b))\right]T_{ij},$

for $T_{ij}$ symmetric matrices. Every column of
$X_{i}(\mathcal{I}_{k_{i}}(b),:)^{T}$ is therefore contained in the
intersection of the ranges of the matrices
$\left[M_{i}^{T}(:,\mathcal{I}_{k_{j}}(b)),-M_{j}^{T}(:,\mathcal{I}_{k_{i}}(b))\right]$,
for $j$ in $I$, $j\neq i$. We now claim that this intersection consists of
only the zero vector. Enumerating three distinct elements of $I$ as
$i,j,\ell$, a vector in the intersection of these spaces is contained in the
null space of the matrix

(B.12)
$\begin{bmatrix}M_{i}^{T}(:,\mathcal{I}_{k_{j}}(b)),&-M_{j}^{T}(:,\mathcal{I}_{k_{i}}(b)),&M_{i}^{T}(:,\mathcal{I}_{k_{\ell}}(b)),&-M_{\ell}^{T}(:,\mathcal{I}_{k_{i}}(b))\end{bmatrix}.$

By lemma B.2, this matrix has null space consisting of the zero vector if
$r\geq 8b+4$. It follows that $X_{i}=0$, which concludes the proof.

∎

## References

  * [1] A. Anandkumar, R. Ge, D. Hsu, S. M. Kakade, and M. Telgarsky, Tensor decompositions for learning latent variable models, J. Mach. Learn. Res., 15 (2014), pp. 2773–2832.
  * [2] A. Anandkumar, P. Jain, Y. Shi, and U. N. Niranjan, Tensor vs. matrix methods: Robust tensor decomposition under block sparse perturbations, in Artificial Intelligence and Statistics, 2016, pp. 268–276.
  * [3] M. Ashraphijuo, V. Aggarwal, and X. Wang, Deterministic and probabilistic conditions for finite completability of low-tucker-rank tensor, IEEE Trans. Inf. Theory, 65 (2019), pp. 5380–5400.
  * [4] M. Ashraphijuo and X. Wang, Fundamental conditions for low-cp-rank tensor completion, J. Mach. Learn. Res., 18 (2017), pp. 2116–2145.
  * [5] M. Ashraphijuo and X. Wang, Characterization of sampling patterns for low-tt-rank tensor retrieval, Ann. Math. Artif. Intell., (2020), pp. 1–28.
  * [6] J. Chen and Y. Saad, On the tensor SVD and the optimal low rank orthogonal approximation of tensors, SIAM J. Matrix Anal. Appl., 30 (2009), pp. 1709–1734.
  * [7] P. Comon, Independent component analysis, a new concept?, Signal processing, 36 (1994), pp. 287–314.
  * [8] S. Gandy, B. Recht, and I. Yamada, Tensor completion and low-n-rank tensor recovery via convex optimization, Inverse Probl., 27 (2011), p. 025010.
  * [9] D. Goldfarb and Z. Qin, Robust low-rank tensor recovery: Models and algorithms, SIAM J. Matrix Anal. Appl., 35 (2014), pp. 225–253.
  * [10] Q. Gu, H. Gui, and J. Han, Robust tensor decomposition with gross corruption, Advances in Neural Information Processing Systems, 27 (2014), pp. 1422–1430.
  * [11] Y. Guan and D. Chu, Numerical computation for orthogonal low-rank approximation of tensors, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 1047–1065.
  * [12] R. A. Harshman, Foundations of the parafac procedure: Models and conditions for an “explanatory” multimodal factor analysis, UCLA Working Papers in Phonetics, 16 (1970), pp. 1–84.
  * [13] C. J. Hillar and L.-H. Lim, Most tensor problems are NP-hard, J. ACM, 60 (2013), p. 45.
  * [14] S. Hu and K. Ye, Linear convergence of an alternating polar decomposition method for low rank orthogonal tensor approximations, arXiv preprint arXiv:1912.04085, (2019).
  * [15] D. Kressner, M. Steinlechner, and B. Vandereycken, Low-rank tensor completion by riemannian optimization, BIT Numer. Math., 54 (2014), pp. 447–468.
  * [16] J. Li, K. Usevich, and P. Comon, Jacobi-type algorithm for low rank orthogonal approximation of symmetric tensors and its convergence analysis, arXiv preprint arXiv:1911.00659, (2019).
  * [17] J. Liu, P. Musialski, P. Wonka, and J. Ye, Tensor completion for estimating missing values in visual data, IEEE Trans. Pattern Anal. Mach. Intell., 35 (2012), pp. 208–220.
  * [18] Z.-Q. Luo and P. Tseng, Error bounds and convergence analysis of feasible descent methods: a general approach, Ann. Oper. Res., 46 (1993), pp. 157–178.
  * [19] C. D. M. Martin and C. F. Van Loan, A jacobi-type method for computing orthogonal tensor decompositions, SIAM J. Matrix Anal. Appl., 30 (2008), pp. 1219–1232.
  * [20] O. Mickelin and S. Karaman, Optimal orthogonal approximations to symmetric tensors cannot always be chosen symmetric, arXiv preprint arXiv:1906.06407, (2019).
  * [21] C. Mu, D. Hsu, and D. Goldfarb, Successive rank-one approximations for nearly orthogonally decomposable symmetric tensors, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 1638–1659.
  * [22] C. Mu, D. Hsu, and D. Goldfarb, Greedy approaches to symmetric orthogonal tensor decomposition, SIAM J. Matrix Anal. Appl., 38 (2017), pp. 1210–1226.
  * [23] I. Oseledets and E. Tyrtyshnikov, A unifying approach to the construction of circulant preconditioners, Linear Algebra Appl., 418 (2006), pp. 435–449.
  * [24] J. Saunderson, V. Chandrasekaran, P. A. Parrilo, and A. S. Willsky, Diagonal and low-rank matrix decompositions, correlation matrices, and ellipsoid fitting, SIAM J. Matrix Anal. Appl., 33 (2012), pp. 1395–1416.
  * [25] J. Saunderson, P. A. Parrilo, and A. S. Willsky, Diagonal and low-rank decompositions and fitting ellipsoids to random points, in Proc. IEEE Conf. Decis. Control, IEEE, 2013, pp. 6031–6036.
  * [26] P. Shah, N. Rao, and G. Tang, Sparse and low-rank tensor decomposition, Advances in Neural Information Processing Systems, 28 (2015), pp. 2548–2556.
  * [27] M. Sørensen and L. De Lathauwer, Fiber sampling approach to canonical polyadic decomposition and application to tensor completion, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 888–917.
  * [28] K. Usevich, J. Li, and P. Comon, Approximate matrix and tensor diagonalization by unitary transformations: convergence of jacobi-type algorithms, arXiv preprint arXiv:1905.12295, (2019).
  * [29] L. Wang, M. T. Chu, and B. Yu, Orthogonal low rank tensor approximation: Alternating least squares method and its global convergence, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 1–19.
  * [30] M. Yuan and C.-H. Zhang, On tensor completion via nuclear norm minimization, Found. Comput. Math., 16 (2016), pp. 1031–1068.
  * [31] T. Zhang and G. H. Golub, Rank-one approximation to high order tensors, SIAM J. Matrix Anal. Appl., 23 (2001), pp. 534–550.