Stable ALS Approximation in the TT-Format for Rank-Adaptive Tensor Completion

TL;DR

This paper introduces SALSA, a stable ALS-based method with regularization for tensor completion, improving reconstruction quality and stability over existing ALS and Riemannian approaches.

Contribution

It develops a singular value regularization for ALS, ensuring stability and enabling effective rank adaptation in tensor completion.

Findings

SALSA significantly improves reconstruction accuracy.

The method is stable under small perturbations.

Numerical experiments demonstrate orders of magnitude better results.

Abstract

Low rank tensor completion is a highly ill-posed inverse problem, particularly when the data model is not accurate, and some sort of regularization is required in order to solve it. In this article we focus on the calibration of the data model. For alternating optimization, we observe that existing rank adaption methods do not enable a continuous transition between manifolds of different ranks. We denote this characteristic as $\textit{instability (under truncation)}$. As a consequence of this property, arbitrarily small changes in the iterate can have arbitrarily large influence on the further reconstruction. We therefore introduce a singular value based regularization to the standard alternating least squares (ALS), which is motivated by averaging in microsteps. We prove its $\textit{stability}$ and derive a natural semi-implicit rank adaption strategy. We further prove that the standard ALS microsteps for completion problems are only stable on manifolds of fixed ranks, and only around points that have what we define as $\textit{internal tensor restricted isometry property, iTRIP}$. In conclusion, numerical experiments are provided that show improvements of the reconstruction quality up to orders of magnitude in the new Stable ALS Approximation (SALSA) compared to standard ALS and the well known Riemannian optimization RTTC.

Tables

Table 1. Table 1: Results for Subsection 9.5 (with arithmetic and geometric variances in brackets) using a (pure) Matlab implementation. For ALS, exact least squares solution are computed, whereas for SALSA, coarse CG is used. Note that most iterations are performed while the rank is not at its maximum yet.

$n=12$		ALS				SALSA
$d$	$c_{\mathrm{sf}}$	${⟨R_C/\Vert M_C\Vert ⟩}_{\mathrm{geo}}$	${⟨R_P/\Vert M_P\Vert ⟩}_{\mathrm{geo}}$	${⟨\text{time}⟩}_{\mathrm{ar}}$	${⟨\text{iter}⟩}_{\mathrm{ar}}$	${⟨R_C/\Vert M_C\Vert ⟩}_{\mathrm{geo}}$	${⟨R_P/\Vert M_P\Vert ⟩}_{\mathrm{geo}}$	${⟨\text{time}⟩}_{\mathrm{ar}}$	${⟨\text{iter}⟩}_{\mathrm{ar}}$
6	2	6.9e-04(2.1)	1.1e-04(2.9)	85(57)	467(151)	3.1e-04(1.1)	6.8e-06(3.7)	81(10)	412(33)
	4	2.0e-04(1.9)	3.9e-05(3.8)	131(86)	527(154)	5.8e-05(1.5)	1.4e-06(2.0)	170(29)	546(58)
	6	4.2e-05(1.8)	4.7e-06(2.8)	215(76)	624(108)	2.3e-05(1.4)	1.5e-06(1.4)	209(23)	581(41)
9	2	7.6e-04(2.1)	2.3e-04(4.4)	145(88)	606(198)	1.8e-04(1.1)	6.9e-06(2.8)	129(9)	455(23)
	4	8.0e-05(1.1)	2.0e-05(1.1)	276(57)	804(81)	2.7e-05(1.2)	7.4e-07(2.0)	298(34)	644(40)
	6	6.6e-05(1.2)	2.5e-05(1.2)	355(94)	834(105)	9.6e-06(1.2)	5.5e-07(1.5)	457(50)	732(48)
15	2	6.6e-04(1.5)	3.3e-04(2.1)	384(113)	951(159)	1.4e-04(1.4)	1.1e-05(1.8)	266(20)	498(31)
	4	3.2e-05(1.1)	6.8e-06(1.1)	961(62)	1396(44)	3.0e-05(1.1)	4.3e-06(1.6)	540(35)	642(32)
	6	2.7e-05(1.1)	8.3e-06(1.1)	1278(85)	1416(47)	5.5e-06(1.8)	2.9e-07(1.6)	1133(92)	827(48)

Table 2. Table 2: Results for Subsection 9.6 (with arithmetic and geometric variances in brackets) using a (pure) Matlab implementation. For ALS, exact least squares solution are computed, whereas for SALSA, coarse CG is used. Note that most iterations are performed while the rank is not at its maximum yet.

$n=12$		ALS				SALSA
$d$	$c_{\mathrm{sf}}$	${⟨R_C/\Vert M_C\Vert ⟩}_{\mathrm{geo}}$	${⟨R_P/\Vert M_P\Vert ⟩}_{\mathrm{geo}}$	${⟨\text{time}⟩}_{\mathrm{ar}}$	${⟨\text{iter}⟩}_{\mathrm{ar}}$	${⟨R_C/\Vert M_C\Vert ⟩}_{\mathrm{geo}}$	${⟨R_P/\Vert M_P\Vert ⟩}_{\mathrm{geo}}$	${⟨\text{time}⟩}_{\mathrm{ar}}$	${⟨\text{iter}⟩}_{\mathrm{ar}}$
	2	6.5e-03(3.2)	5.3e-03(3.4)	72(46)	423(139)	6.5e-05(1.8)	1.1e-06(5.8)	69(16)	359(66)
	4	4.4e-03(1.2)	4.1e-03(1.3)	34(34)	258(114)	3.7e-06(1.8)	1.9e-08(5.5)	133(32)	497(87)
	6	3.8e-03(1.2)	3.6e-03(1.2)	47(40)	288(122)	4.7e-07(2.3)	4.5e-09(4.0)	157(34)	520(88)

Equations

$$M{|}_P=\{M_p{\}}_{p\in P}\in {\mathbb{R}}^P\cong {\mathbb{R}}^{|P|}\text{mbox}forP\subset \mathcal{I}:=\{1,\dots ,n\}\times \{1,\dots ,m\},$$

$$\mathrm{minimize}$$

$$\mathrm{subject}\mathrm{to}$$

$${\tau }_r:(X,Y)\mapsto A=XY\in {\mathbb{R}}^{n\times m}\text{mbox}for(X,Y)\in {\mathcal{D}}_r:={\mathbb{R}}^{n\times r}\times {\mathbb{R}}^{r\times m}.$$

$${\mathcal{T}}_r:=\{A\mid \mathrm{rank}(A)=r\},r=0,\dots ,\min (n,m).$$

$${\mathcal{M}}_r^{(1)}(X,Y)$$

$${\mathcal{M}}_r^{(2)}(X,Y)$$

$$f_{\mathcal{M}}:{\mathbb{R}}^{\mathcal{I}}\to {\mathbb{R}}^{\mathcal{I}},f_{\mathcal{M}}(A):={\tau }_{r(A)}\circ {\mathcal{M}}_{r(A)}\circ {\tau }_{r(A)}^{-1}(A),$$

$$M:=\left(\begin{array}{ccc}\text{framebox}?& 1.1& 0.9\\ 1& 1& 1.1\\ 1.1& 1& 1\end{array}\right),A(\epsilon ):=\left(\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right)+\epsilon \left(\begin{array}{ccc}0.5+a& 0.5+a& -a\\ 1+a& 1+a& -1-a\\ 1-a& 1-a& -1+a\end{array}\right),$$

$$A(0)=\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right)\left(\begin{array}{ccc}1& 1& 1\end{array}\right),A(\epsilon )=\left(\begin{array}{cc}1& 0.5+a\\ 1& 1+a\\ 1& 1-a\end{array}\right)\left(\begin{array}{ccc}1& 1& 1\\ \epsilon & \epsilon & -\epsilon \end{array}\right)\text{mbox}if\epsilon >0.$$

$$B(0)=\left(\begin{array}{ccc}1.05& \ast & \ast \\ 1.05& \ast & \ast \\ 1.05& \ast & \ast \end{array}\right),B(\epsilon )=\left(\begin{array}{ccc}1+\frac{1}{40a}& \ast & \ast \\ 1.0& \ast & \ast \\ 1.1& \ast & \ast \end{array}\right)\text{mbox}if\epsilon >0.(\ast \text{mbox}issomevalue)$$

$$f_{{\mathcal{M}}^{(2)}}(A)={\tau }_r({\mathcal{M}}_r^{(2)}(U,\Sigma V^T))=U{\mathrm{argmin}}_Y\Vert UY-M{\Vert }_P=\underset{\epsilon \searrow 0}{\lim }{f}_{{\mathcal{M}}^{(2)}}(A{|}_{{\sigma }_r=\epsilon })$$

$$f_{{\mathcal{M}}^{(2)}}(A{|}_{{\sigma }_r=0})={\tau }_{r-1}({\mathcal{M}}_{r-1}^{(2)}(U_c,{\Sigma }_c{V}_c^T))=U_c{\mathrm{argmin}}_Y\Vert U_{c}Y-M{\Vert }_P$$

$$f_{\mathcal{M}}^{\ast }(A)$$

$${\mathbb{V}}_{A,\omega }$$

$${\hat{f}}_{\mathcal{M}}(A):=\frac{1}{|V_{\omega }|}{\int }_{V_{\omega }}{f}_{\mathcal{M}}\left((u_1\mid u_2+\Delta u_2)\Sigma V^T\right)\mathrm{d}\Delta u_2$$

$$=\underset{\text{mbox}optimization}{\underbrace{u_1{u}_1^{T}M}}+\underset{\text{mbox}regularization}{\underbrace{(1-{\alpha }_{\omega }{)}^2{u}_2{u}_2^{T}M}}+\underset{\text{mbox}replenishment}{\underbrace{\frac{2{\alpha }_{\omega }-{\alpha }_{\omega }^2}{m-2}(I_m-u_1{u}_1^T-u_2{u}_2^T)M}}$$

$$V_{\omega }:=\left\{\Delta u_2\mid \Vert (u_1\mid u_2+\Delta u_2)\Sigma V^T-A{\Vert }_F=\omega ,(u_1\mid u_2+\Delta u_2)\text{mbox}hasorthonormalcolumns\right\}$$

$$\frac{1}{|V_{\omega }|}{\int }_{V_{\omega }}{f}_{\mathcal{M}}\left((u_1\mid \Delta u_2)\Sigma V^T\right)\mathrm{d}\Delta u_2=u_1{u}_1^{T}M+\frac{1}{m-1}(I_m-u_1{u}_1^T)M,$$

$$\text{mbox}wherehere{V}_{\omega }:=\{\Delta u_2\mid (u_1\mid \Delta u_2)\text{mbox}hasorthonormalcolumns\}.$$

$$f_{\mathcal{M}}((u_1\mid u_2+\Delta u_2)\Sigma V^T)=(u_1\mid u_2+\Delta u_2)(u_1\mid u_2+\Delta u_2{)}^{T}M.$$

$${\int }_{V_{\omega }}{f}_{\mathcal{M}}={\int }_{V_{\omega }}{u}_1{u}_1^{T}M+{\int }_{V_{\omega }}(1-{\alpha }_{\omega }{)}^2{u}_2{u}_2^{T}M+{\int }_{V_{\omega }}\Delta u_2^{\perp }{\Delta u_2^{\perp }}^{T}M$$

$${\int }_{V_{\omega }}\Delta u_2^{\perp }{\Delta u_2^{\perp }}^{T}M={\int }_{{\beta }_{\omega }{S}^{m-2}}(Hx)(Hx{)}^{T}M\mathrm{d}x=H{H}^T\frac{2{\alpha }_{\omega }-{\alpha }_{\omega }^2}{m-2}|V_{\omega }|M$$

$${\mathrm{argmin}}_{\tilde{V}}\frac{1}{|V_{\omega }|}{\int }_{V_{\omega }}\Vert (u_1\mid u_2+\Delta u_2)\tilde{V}-M{\Vert }_F^2\mathrm{d}\Delta u_2=(u_1^{T}M\mid (1-{\alpha }_{\omega })u_2^{T}M)$$

$$|V_{\omega }|\tilde{V}$$

$$=\left(u_1^{T}M|V_{\omega }|\mid u_2^{T}M|V_{\omega }|+{\int }_{V_{\omega }}-{\alpha }_{\omega }{u}_2^{T}M+\Delta u_2^{\perp }M\text{mbox}d\Delta u_2\right)$$

$$=|V_{\omega }|\left(u_1^{T}M\mid (1-{\alpha }_{\omega })u_2^{T}M\right).$$

$$V_{\omega }(A):=\{\Delta A\mid A+\Delta A\in {\mathcal{T}}_r,\Vert \Delta A{\Vert }_F\le \omega \},r=\text{mbox}rank(A)$$

$$\Vert \Delta A{\Vert }_F^2$$

$$=\Vert \omega (\Delta XY+X\Delta Y){\Vert }_F^2+{\left(\mathcal{O}({\omega }^2)\right)}^2$$

$$C(\tilde{X},\tilde{Y})$$

$${\mathbb{V}}_{\omega }(X,Y)$$

$$({\mathcal{M}}^{(1)}{)}^{\ast }(X,Y)$$

$$({\mathcal{M}}^{(2)}{)}^{\ast }(X,Y)$$

$$V_{\omega }^{(n,m)}:=\{X\in {\mathbb{R}}^{n\times m}\mid \Vert X{\Vert }_F=\omega \}.$$

$${\int }_{V_{\omega }^{(n,m)}}{X}^{T}HX\mathrm{d}X=\frac{{\omega }^2|V_{\omega }^{(n,m)}|}{nm}\text{mbox}trace(H)I_m,|V_{\omega }^{(n,m)}|:={\int }_{V_{\omega }^{(n,m)}}1.$$

Full text

Keywords. tensor completion, MPS, tensor train, TT, hierarchical Tucker, HT, alternating optimization, ALS, high-dimensional, low rank, SVD, ill-posedness, stability

AMS subject classifications. 15A18, 15A69, 65F22, 90C06, 90C31

Introduction

Low rank tensor completion is a highly ill-posed inverse problem. In order for any recovery to succeed, regularity assumptions are required. This is either achieved by adding certain penalty terms, or by using an explicit reduction of degrees of freedom, which in the context of high-dimensional tensors can be obtained by using low rank representations, cf. [15, 17, 19]. The main goal of this work is to derive a rank adaptive method for tensor completion and moreover to discuss the benefits of such as well as the reasons why heuristics tend to be insufficient. Since the concept of stability in the sense of Definition 1.1 has not yet been considered in literature, the initial part is dedicated to the simpler matrix case (Sections 2 and 3) in order to provide an easier access. The subsequent analysis will focus on the importance of these concepts to least squares tensor completion where the calibration of model complexity is more challenging.

Introduction to Stability for Ill-posed Inverse Problems through the Example of Matrix Completion

In the setting of low rank matrix completion, the target of recovery is a matrix $M\in\mathbb{R}^{n\times m}$ which is only observable at points

[TABLE]

where $P$ is a given, fixed sampling set, which we hence can not enlarge. One very strict regularity assumption is given by $\mathrm{rank}(M)=r$ for some sufficiently small $r\in\mathbb{N}$, which leads to the minimization problem

[TABLE]

$\mbox{where }\|B\|^{2}_{\rm{P}}:=\sum_{i\in P}B^{2}_{i}$ for matrices $B$. A favorable data model for this task is the low rank representation, i.e. a function

[TABLE]

Since every matrix has a unique rank, we can partition the target space $\mathbb{R}^{n\times m}$ into the disjoint subsets

[TABLE]

With $\mathrm{image}(\tau_{r})=\bigcup_{\widetilde{r}\leq r}\mathcal{T}_{\widetilde{r}}$ in mind, the optimization is performed on the representation or data space $\mathcal{D}_{r}$. In an alternating least squares (ALS) method for example, one then applies two optimization methods $\mathcal{M}^{(1)}$, $\mathcal{M}^{(2)}$,

[TABLE]

which in this context are called microsteps. Note that the $\mathrm{argmin}$ is not necessarily unique, and we choose the element minimizing the Frobenius norm $\|XY\|_{F}$. Formally, for each value of the matrix rank $r$, every single $\mathcal{M}_{r}^{(1)}$, $\mathcal{M}_{r}^{(2)}$ is a different function.

For most realistic applications, it is more reasonable to relax the regularity assumption to $M$ being nearly rank $r$, which means that after $r$ entries, the singular values of $M$ become sufficiently smaller. Ultimately, if no assumptions are made, the appropriate model complexity is a matter of the quality and magnitude of $P$ with respect to $M$. Yet in the general case, the missing structure of given data hardly allows to obtain knowledge about this relation. Therefor, since overestimating the model complexity (i.e. the rank) ultimately leads to flawed results, a cautious learning process is required to adapt such, as rank increasing strategies already suggest. Due to the difficult nature of the problem, we do not expect to be able to find the global minimizer, but instead focus on single aspects that are likely to improve the approximation quality.

Each adaption of the rank during the optimization will cause the algorithm to change between data spaces $\mathcal{D}_{r}$. Intuitively, considering that the generated spaces $\mathcal{T}_{r}$ have pairwise distance [math] within $\mathbb{R}^{n\times m}$, one would want that a change of rank does not have large impact, given the problematic nature of overfitting. This, however, is not true for both $\mathcal{M}^{(1)},\mathcal{M}^{(2)}$, while arbitrarily small perturbations of the iterate $A=\tau_{r}(X,Y)$ may change its rank. For ill-posed inverse problems, we hence propose the following concept:

Definition 1.1** (Stability).**

Let $\mathcal{M}$ be a method that maps any rank $r$ to a function $\mathcal{M}_{r}:\mathcal{D}_{r}\rightarrow\mathcal{D}_{r}$ (the optimization method for fixed rank). We define the following properties:

• •

$\mathcal{M}$* is called representation independent, if $\tau_{r}(\mathcal{M}_{r}(G))=\tau_{r}(\mathcal{M}_{r}(\widetilde{G}))$ for all $r$ and $G,\widetilde{G}\in\mathcal{D}_{r}$ with $\tau_{r}(G)=\tau_{r}(\widetilde{G})$. We then define $\tau_{r}^{-1}$ to map to one possible representation (we want to circumvent the use of equivalence classes).*

• •

$\mathcal{M}$* is called fixed-rank stable, if it is representation independent and for any fixed rank $r$, the map $\tau_{r}\circ\mathcal{M}_{r}\circ\tau^{-1}_{r}:\mathcal{T}_{r}\rightarrow\mathbb{R}^{\mathcal{I}}$ is continuous.*

• •

$\mathcal{M}$* is called stable, if it is representation independent and the function*

[TABLE]

where $r(A)$ is the rank of $A$, is continuous.

This definition of stability is not restricted to matrix completion. Data spaces, the rank (possibly generalized to any model complexity) and the method can be replaced by any appropriate type, in particular tensor completion in hierarchical tensor formats (cf. Section 4). The only assumption should usually be a nestedness of spaces, i.e. $\mathcal{T}_{r}\subset\mathcal{T}_{\widetilde{r}}$ whenever $r\preceq\widetilde{r}$ (entry-wise inequality).

Properly calibrating the rank $r$ for unstable methods in the context of ill-posed inverse problems may lead to complications. Most of the operators applied to representations are stable, e.g. truncations based on matrix singular values. The situation however changes if we apply the partial optimization (or micro-) step $\mathcal{M}^{(1)}$ or $\mathcal{M}^{(2)}$ on a low rank representation:

Example 1.2** (Instability of alternating least squares matrix completion steps).**

Let $a\in\mathbb{R}\setminus\{0,1\}$ be a possibly very small parameter. We consider the target matrix $M$ and an $\varepsilon$-dependent initial approximation $A=A(\varepsilon)$

[TABLE]

where the entry $M_{1,1}$ (the question mark above) is not known or given. The matrix $M$ is of rank $3$ and $A(\varepsilon)$ is of rank $r=1$ for $\varepsilon=0$ and of rank $r=2$ otherwise. We seek a best approximation of (at most) rank $2$ in the least squares sense for the known entries of $M$. In a single ALS step, as defined by \tagform@1.2, we replace $Y(\varepsilon)$ of the low rank representation $A(\varepsilon)=X(\varepsilon)Y(\varepsilon)$ by the local minimizer, where in this case

[TABLE]

This optimization yields a new matrix, $B(\varepsilon)=f_{\mathcal{M}^{(2)}}(A(\varepsilon))=\tau_{r}\circ\mathcal{M}^{(2)}_{r}\circ\tau^{-1}_{r}(A(\varepsilon))$ (independently of the chosen representation), given by

[TABLE]

Now let $a$ be fixed and let $\varepsilon$ tend to zero so that the initial guess $A(\varepsilon)\to A(0)$. However, $B(\varepsilon)\nrightarrow B(0)$, thus violating the stability. Furthermore, the rank two approximation $B(\varepsilon)$, given an arbitrary, fixed $\varepsilon>0$, diverges as $a\to 0$, in particular it is not convergent although the initial guess $A(\varepsilon)$ converges to a rank two matrix as $a\to 0$. Thus, the microstep is not even stable for fixed rank. We want to stress that the initial guess is bounded for all $\varepsilon,a\in(0,1)$, but the difference between $B(0)$ and $B(\varepsilon)$ is unbounded for $a\to 0$ (cf. Definition 5.10). The unboundedness can be remedied by adding a regularization term in the least squares functional, e.g. $+\|XY\|$, but the ALS step remains unstable.

This example likewise demonstrates that ALS for tensor completion is not stable and thus, as discussed before, problematic when adapting the rank (cf. Section 4). We will further show that this is not a marginal phenomenon, but occurs systematically during any rank change (cf. Example 2.1).

Relation to Other Matrix and Tensor Methods

Whenever a tensor is point-wise available, algorithms such as the TT-SVD [33] can just establish the exact rank based on its very definition or a reliable rank estimate as well as representation can be obtained through cross-approximation methods, a setting in which the subset of used entries can be chosen freely [32, 4].

If only indirectly given, adapting the rank of the sought low rank tensor can still be straight-forward, e.g. when the rank has to be limited only due to computational complexity, while in principle the exact solution is desired [5, 2]. Here, an optimal regulation of thresholding parameters becomes most important. This mainly includes classical problems that have been transferred to large scales. These may for example be solved with iterative methods [3, 1, 29], which naturally increase the rank and rely on subsequent reductions, or also by rank preservative optimization, such as alternating optimization [21, 9, 37, 11], possibly combined with a separate rank adaption.

Provided that the tensor restricted isometry property holds, the task may be interpreted as distance minimization with respect to a norm that is sufficiently similar to the Frobenius norm and analyzed based on compressed sensing [35]. Black box tensor completion for a fixed sampling set, however, requires a certain solution to a positive-semi definite linear system. Hence neither an exact solution is reasonable nor does any norm equivalence hold. Thus, the available data is easily misinterpreted, the more so if the rank is overestimated, and truncation based algorithms, including DMRG [21, 24], are misled.

Nuclear norm minimization, being closely related to compressed sensing as well, has a very strong theoretical background [7, 8, 16, 36] for the matrix case. These approaches rely on a direct adaption of the target function, that is convex relaxation. Yet it appears that they are outperformed in practice by alternating least squares approaches [23] and the simplifications required for an adaption to tensors [28, 38, 12] do not seem to allow for an appropriate generalization [31]. Also the approaches which themselves retreat to alternating least squares [20] treat the iterations as necessity for the minimization of an objective function with regularization term. The penalty term is, as usual, based on the singular values of the output of a microstep (a posteriori), as it is also the case for the work [40] on tensor completion through Riemannian optimization. Although their term may appear similar to the term we derive (cf. Theorem 3.5), the stability property, on the contrary, requires that the penalty term depends on the current singular values before the microstep (a priori), which is an essential difference: we explicitly allow and exploit small singular values instead of penalizing such. In that sense, we treat each update and adaption as part of a learning progress, where the magnitudes of singular values indicate in some respects an uncertainty of approximation. To the best of our knowledge, this point of view has not yet been considered in literature and hence a concept of stability as we define it has not been investigated.

For fixed or uniform rank, there have been least squares based proposals in hierarchical tensor formats [26, 39] as well. The essential adaption of the rank however, including similar matrix approaches, is rarely considered, all the less in numerical tests, and remains an open problem in this setting. A mentionable approach so far is the rank increasing strategy [44, 14] and its regularization properties are a first starting point for this article.

The rest of the article is organized as follows: In Section 2, we further investigate instability and exemplarily analyze approaches towards it in the matrix case. Based on this insight, we motivate a variational residual function, derive its minimizer and present a stable algorithm for matrix completion in Section 3. In Section 4, we begin to generalize former results to high dimensional tensors, yet essentially work in three dimensions. In the main Section 5, we then derive its minimizer and prove stability (Theorem 5.14) for the thereby obtained regularized microsteps, further analyzing these in Section 6. Subsequently, in Section 7, these results are transferred back to arbitrarily dimensional tensors. Section 8 finishes with the necessary details for the algorithm, including its rank adaption as it is naturally given through stable alternating least squares. Comprehensive numerical tests (exclusively for unknown ranks) are provided in Section 9. Detailed tables with the numerical values shown in figures can be found in the appendix.

Instability and Approaches to Resolve the Problem

The following example shows that instability can be observed systematically during rank changes in ALS, or more general, in any such range based optimization. In that sense, the implied complications for ill-posed, inverse problems may frequently occur:

Example 2.1** (ALS for ill-posed, inverse problems is unstable).**

Consider the microstep $\mathcal{M}^{{(2)}}$ as in \tagform@1.2. Let $U\in\mathbb{R}^{n\times r}$, $V\in\mathbb{R}^{m\times r}$ be orthogonal, such that $U\Sigma V^{T}$ is a truncated SVD of a rank $r$ matrix $A=\tau_{r}(U,\Sigma V^{T})\in\mathbb{R}^{n\times m}$. We now let $\sigma_{r}\rightarrow 0$, $\sigma_{r}>0$, such that in the limit $A^{\ast}:=A|_{\sigma_{r}=0}$ has rank $r-1$. The update is independent of this last singular value though:

[TABLE]

However, if $\sigma_{r}=0$, then $A|_{\sigma_{r}=0}$ has rank $r-1$ and a truncated SVD $U_{c}\Sigma_{c}V_{c}^{T}$. Hence, the update

[TABLE]

is in general different from the limit \tagform@2.1, given that the range of $U$ is different from $U_{c}$. The same holds for an analogous update $\mathcal{M}^{{(1)}}$ of $X=U\Sigma$. Note that these updates are indeed representation independent.

The microsteps of ALS in the tensor case behave in the same way. The only difference is that there are two tuples of singular values $\sigma^{(\mu-1)}$ and $\sigma^{(\mu)}$ adjacent to a core $G_{\mu}$ (cf. Lemma 5.2). There may be many ways to stabilize the microsteps. However, we aim for an as little distorting as practical way to do this. A quite successful approach for completion has been the rank increasing strategy, e.g. [44]. Therein, the model complexity is slowly increased, step by step attempting a more distinct approximation. Thereby, local minima that correspond to overfitting are avoided.

A similar kind of effect can be achieved by assuming an uncertainty on the current iterate, or, equivalently, averaging the tensor update function. That way, the level of regularization can be adapted continuously and is less dependent on the discrete rank but the more meaningful, real-valued singular values. We will first view this in a minimal fashion for the matrix case and the method $\mathcal{M}^{{(2)}}$ defined by \tagform@1.2 (for the remainder of this section called $\mathcal{M}$ for simplicity). With this approach, we can motivate an algorithm that is stable under truncation and allows to straightforwardly adapt ranks. While other methods in fact keep distance from the border of one manifold of a fixed rank, we aim to optimize, in a certain sense, continuously between manifolds of different ranks.

Assuming local integrability of $f_{\mathcal{M}}$ (as defined in \tagform@1.3), we obtain that the variational function

[TABLE]

is continuous within $image(\tau_{\widetilde{r}})$, where $\|\cdot\|_{F}$ is the Frobenius norm and $\widetilde{r}$ may be considered an upper bound to the rank (cf. Figure 2).

However, this function does not preserve low rank structure, apart from appearing to be to complicated to evaluate, and therefore we cannot find a method $\mathcal{M}^{\ast}$ for which $f_{\mathcal{M}}^{\ast}=f_{\mathcal{M}^{\ast}}$. Still, all simplifications which we will make remain subject to this motivation.

Investigations into the Connection between Averaging and Stability

Although the following examples are heavily narrowed down, we observe which kind of simplifications we may apply in order to obtain a feasible, but still stable method (cf. Section 3). First, we consider a scenario in which we limit the perturbation that the left singular vectors $U$ receive due to the variation of $A$ to only one component as it is approximately the case if $\sigma_{1}\gg\sigma_{2}\approx\omega$.

Lemma 2.2** (Variational low rank matrix approximation).**

Let $\mathcal{M}$ be defined by \tagform@1.2 and $P=\mathcal{I}$ (full sampling). Let further $A=U\Sigma V^{T}\in\mathbb{R}^{n\times m}$ be of rank two, given by its SVD components $U=(u_{1}\mid u_{2})$, $\Sigma=diag(\sigma_{1},\sigma_{2})$ and $V$ as well as $M\in\mathbb{R}^{m\times m}$ arbitrary and $0<\omega<\sqrt{2}\sigma_{2}$. Then

[TABLE]

for $\alpha_{\omega}=\frac{\omega^{2}}{2\sigma_{2}^{2}}$, such that $\alpha_{\omega}\rightarrow 1$ if $\omega\rightarrow\sqrt{2}\sigma_{2}$. Alternatively, considering complete uncertainty concerning the second singular vector, we obtain

[TABLE]

Proof.

We parameterize $V_{\omega}$. First, $\omega=\|(u_{1}\mid u_{2}+\Delta u_{2})\Sigma V^{T}-A\|_{F}=\|\Delta u_{2}\|_{F}\ \sigma_{2}$ and hence $\|\Delta u_{2}\|_{F}=\frac{\omega}{\sigma_{2}}$. By orthogonality conditions, we obtain $\Delta u_{2}=-\alpha_{\omega}u_{2}+\Delta u^{\perp}_{2}$ with $\Delta u^{\perp}_{2}\perp\mbox{range}(U)$ for a fixed $\alpha_{\omega}=\frac{\omega^{2}}{2\sigma_{2}^{2}}$. Hence, $V_{\omega}$ is an $(m-3)$-sphere of radius $\beta_{\omega}=\sqrt{\frac{\omega^{2}}{\sigma_{2}^{2}}-\alpha_{\omega}^{2}}$, that is $\beta_{\omega}S^{m-2}$. The update for each instance of $\Delta u^{\perp}_{2}$ is given by

[TABLE]

We integrate this over $V_{\omega}$ and obtain

[TABLE]

since all integrals of summands which contain $\Delta u^{\perp}_{2}$ exactly once vanish due to symmetry. We can simplify the last summand with Lemma 3.4 to

[TABLE]

for a linear, orthonormal map $H$ that maps $x\in\beta_{\omega}S^{m-2}$ to $\Delta u^{\perp}_{2}$, that is, embeds it into $\mathbb{R}^{m}$. One can then conclude that $HH^{T}=I_{m}-u_{1}u_{1}^{T}-u_{2}u_{2}^{T}$, since the rank of $H$ is $m-2$ and $range(H)\perp range(U)$. The division by $|V_{\omega}|$ then finishes the first part. The second part is analogous. ∎

We can observe that, in this case, choosing $\omega$ close to $\sigma_{2}$, or in that sense a small $\sigma_{2}$, will filter out influence of $u_{2}$. This is indeed in agreement to the update which the rank $1$ best-approximation to $A$ would yield111Note that we fixed $\|\Delta u_{2}\|_{F}=\omega$ for simplicity as well as that for $\omega>\sqrt{2}\sigma$, Example 2.2 does not make sense. Allowing perturbations up to a magnitude $\omega$ will prohibit that the influence of $u_{2}$ vanishes completely, hence $u_{2}$ is never actually truncated.. More importantly, the result $\widehat{f}_{\mathcal{M}}(A)$ in \tagform@2.2 is not low rank, yet close to the rank $2$ approximation $U(u_{1}^{T}M\mid(1-\alpha_{\omega})^{2}u_{2}^{T}M)$, in which the first component $U$ has remained the same. While the variational model as in \tagform@2.2 remains the basic idea, it is too complicated to be used for the derivation of a stable method $\mathcal{M}^{\ast}$. We instead consider a slightly modified approach (which will be used in the following Sections 3 and 5).

Lemma 2.3** (Low rank matrix approximation using a variational residual function).**

In the situation of Lemma 2.2, we have

[TABLE]

Proof.

Let $\widetilde{V}$ be the corresponding minimizer. With the same derivation as in Lemma 2.2, we obtain

[TABLE]

∎

Comparing this to the rank $2$ approximation of the previous result \tagform@2.2, we observe that solely $(1-\alpha_{\omega})^{2}$ has been replaced by $1-\alpha_{\omega}$. For our purpose, these terms are sufficiently similar for small $\alpha_{\omega}\geq 0$.

Remark 2.4** (Replenishment and lower limit).**

The so called replenishment term in \tagform@2.2 points at an important aspect which we analyze in Section 6. We later bypass related problems through an additional manipulation of singular values, the intensity of which is proportional to the current residual.

We here refer to a Matlab implementation of a (superficially random) Monte Carlo approach to the unsimplified variational microstep $f^{\ast}_{\mathcal{M}}$ as in \tagform@2.2 for matrix completion, which is linked on the personal webpage of the author Sebastian Krämer222by the time the paper is written, the address is www.igpm.rwth-aachen.de/team/kraemer. Likewise, an implementation of the final algorithm SALSA (Algorithms 2 and 4), which is developed from the idea in Lemma 2.3, can be found for the matrix case as well as for the tensor case.

Stable Alternating Least Squares Microsteps for Matrix Completion

In this section, we adapt the target function of each microstep $\mathcal{M}$ in order to obtain a stable method $\mathcal{M}^{\ast}$. When performing ALS, instead of just one specific iteration point $A\in\mathcal{T}_{r}$, the results in Section 2.1 suggest to instead consider a set $V_{\omega}(A)$ of variations, or uncertainty, $\Delta A$ along the manifold $\mathcal{T}_{r}$ of rank $r$ matrices:

[TABLE]

The initial idea is slightly similar to gradient sampling (e.g. [6]), but is then pursued differently. Let $A=\tau_{r}(X,Y)$ and $(\Delta X,\Delta Y)$ such that $A+\Delta A=\tau_{r}(X+\omega\Delta X,Y+\omega\Delta Y)$. Then

[TABLE]

The term $\|\Delta XY+X\Delta Y\|^{2}_{F}$ can be approximated, assuming the angles between the three summands are small, by $\|\Delta XY\|^{2}_{F}+\|X\Delta Y\|^{2}_{F}$. This and Lemma 2.3 then motivate the following definition.

Definition 3.1** (Variational residual function).**

Let $\omega\geq 0$, $M$ the target matrix, $P$ the sampling set and $A=XY$ the current iterate. We define the variational residual function $C:=C_{M,P,X,Y}:\mathcal{D}_{r}\rightarrow\mathbb{R}$ by

[TABLE]

Minimizer of the Variational Residual Function for Matrices

We define our modified methods as

[TABLE]

with $C=C_{M,P,X,Y}$ as in \tagform@3.2. It should further be noted that $\mathbb{V}_{\omega}$ does not depend on the unknown $\widetilde{X}$, $\widetilde{Y}$, respectively, but on the current iterate. We will later see that the minimizers are unique, but for formality we again use the minimization of Frobenius norm of the iterate as secondary criterion.

Lemma 3.2**.**

The two methods $(\mathcal{M}^{(1)})^{\ast}$ \tagform@3.3 and $(\mathcal{M}^{(1)})^{\ast}$ \tagform@3.4 are representation independent.

Proof.

We later prove this for the generalized tensor case (cf. Lemma 5.4). ∎

Here, and throughout the remainder of the paper, we use the following, convenient notations, since we often have to reshape, restrict or project objects.

Definition 3.3** (Restrictions).**

For any object $A\in\mathbb{R}^{I}$ and index set $S\subset I$, we use $A|_{S}\in\mathbb{R}^{S}$ as restriction. For a matrix $M$, let $M_{:,i}$ be its $i$-th column and $M_{i,:}$ be its $i$-th row. Furthermore, whenever we apply a restriction to an object or reshape it, we also use the same notation to correspondingly modify index sets (cf. Theorem 3.5).

Lemma 3.4** (Integral over all variations).**

Let $n,m\in\mathbb{N}$, $\omega\geq 0$ and $H\in\mathbb{R}^{n\times n}$ be a matrix as well as

[TABLE]

Then

[TABLE]

Proof.

Let $Y$ be the result of the above integral. Then

[TABLE]

Due to symmetry, for some $\alpha\in\mathbb{R}$, we have

[TABLE]

Since the second term is the trace of the first one, it follows that $\alpha=\omega^{2}|V_{\omega}^{(n,m)}|/(nm)$. We can hence simplify

[TABLE]

which is the to be proven statement. ∎

We now derive the minimizer of the variational residual function for matrices \tagform@3.2. One can use an SVD of the current iterate $A=U\Sigma V^{T}$ for simplification. In this case, $\mathbb{V}_{\omega}$ as in Definition 3.1 takes the easier forms

[TABLE]

Theorem 3.5** (Minimizer of the ALS variational residual function for matrices).**

*Let $A\in\mathbb{R}^{n\times m}$ be the current iterate with $r=\mathrm{rank}(A)$ and let $U\Sigma V^{T}$ be a truncated SVD, $U\in\mathbb{R}^{n\times r}$, $\Sigma\in\mathbb{R}^{r\times r}$, $V^{T}\in\mathbb{R}^{r\times m}$, of $A$. Let further $C$ be the variational residual function as in \tagform@3.2. The minimizer $X^{+}$ of $\widetilde{X}\mapsto C_{M,P,U\Sigma,V^{T}}(\widetilde{X},V^{T})$ is given by *

[TABLE]

*where $P_{i,:}:=\{p_{2}^{(k)}\mid p_{1}^{(k)}=i,\ k=1,\ldots,|P|\}$ is the corresponding part of the index set $P$. The minimizer $Y^{+}$ of $\widetilde{Y}\mapsto C_{M,P,U,\Sigma V^{T}}(U,\widetilde{Y})$ is given by *

[TABLE]

where $P_{:,j}:=\{p_{1}^{(k)}\mid p_{2}^{(k)}=j,\ k=1,\ldots,|P|\}$. The constants $\zeta_{1}$ and $\zeta_{2}$ only depend on the proportions of the representation and sampling set (cf. Remark 3.6).

The factors $\frac{|{P_{i,:}}|}{m}$ and $\frac{|P_{:,j}|}{n}$ normalize the penalty terms to the particular magnitudes of the corresponding shares of the sampling set $P$ and hence the standard ALS part. The factors $\zeta_{1},\zeta_{2}\in(0,1)$ are in turn independent of $i,j$, respectively. The ratio of both terms equals the ratio of the mode sizes $n,m$. The reason for the latter scaling will become apparent in the tensor case (see comments to Theorem 5.6).

Proof.

We search for $Y^{+}:=\mathop{\rm argmin}_{\widetilde{Y}}C_{M,P,U,\Sigma V^{T}}(U,\widetilde{Y})$ (the counterpart for $X^{+}$ is analogous). Substituting

[TABLE]

we can (up to a constant Jacobi determinant) restate $C$ as

[TABLE]

Each of the independent columns in the minimizer is restated as

[TABLE]

Let $j$ be arbitrary but fixed from now on. For any vector $x$, it is $\|x\|_{\mbox{vec}(P_{:,j})}=\|H(j)x\|_{F}=x^{T}H(j)x$ for a diagonal, square matrix $H(j)\in\mathbb{R}^{n\times n}$ with

$H(j)_{(s),(s)}=\delta_{s\in P_{:,j}}$ (hence $H(j)^{2}=H(j)$). Using the normal equation, we obtain $Y^{+}_{:,j}=W^{-1}b$, where

[TABLE]

and

[TABLE]

In both $W$ and $b$, any perturbation that appears only one-sided in the expanded product vanishes due to symmetry of $\mathbb{V}_{\omega}$. Thus $b=|\mathbb{V}_{\omega}|\ U_{P_{:,j},:}^{T}\ M_{P_{:,j},j}$ and

[TABLE]

Now, let $\ell=\#_{X}=nr,\ k=\#_{Y}=rm$. Since $\mathbb{V}$ is a version of the $(\ell+k)$-sphere, we can use the following integration formula: For any $n,m,k\in\mathbb{N}$ let $f:\mathbb{R}^{n+m}\rightarrow\mathbb{R}^{k}$ be a sufficiently smooth function and $S^{v-1}_{\omega}$ be the $v$-sphere of radius $\omega$. Then

[TABLE]

For a function $f$ we then obtain

[TABLE]

When $f$ is independent of $\Delta Y$, this simplifies to

[TABLE]

We further use the identity

[TABLE]

We apply these and Lemma 3.4 for different $f_{0}=U^{T}\ H(j)\ U$ ($\delta=0$) and $f_{1}=(\Delta X\Sigma^{-1})^{T}\ H(j)\ (\Delta X\Sigma^{-1})$ ($\delta=1$). For the summands $W(0)+W(1)=W$ this then yields

[TABLE]

for $c=\frac{4\pi^{(k+\ell)/2}}{\Gamma(\ell/2)\Gamma(k/2)}$. The constant matrices $C_{H}$ are given by

[TABLE]

Furthermore, it is $|\mathbb{V}_{\omega}|=c\frac{\omega^{\ell+k}}{\ell+k}\nu(\ell-1,k-1)$. Factoring out this base volume in $W=|\mathbb{V}_{\omega}|\widetilde{W}$ by using properties of the $\Gamma$ function, one derives:

[TABLE]

Restating the result again as a least squares problem finishes the proof. ∎

Remark 3.6** (Specification of constants).**

Let $\#_{X}:=nr,\ \#_{Y}:=rm$ be the sizes of the components in the matrix decomposition. The constants in Theorem 3.5 are given by

[TABLE]

However, when changing the rank, this would impose a slight offset in continuity of both $f_{{\mathcal{M}^{(1)}}^{\ast}}$ and $f_{{\mathcal{M}^{(2)}}^{\ast}}$. This problem is simply resolved by substituting $\omega$ by $\widetilde{\omega}$ properly for each value $r$, such that we can set and normalize

[TABLE]

We will still just write $\omega$ although we replace $\zeta_{1}$ and $\zeta_{2}$ by the adapted values.

Manipulating $\omega$ does not change the representation independency of the two methods $(\mathcal{M}^{(1)})^{\ast}$ and $(\mathcal{M}^{(2)})^{\ast}$, since for fixed $r$, the value $\omega>0$ just remains an arbitrary constant (cf. Lemma 5.4). We arrive at the main result for the matrix case:

Theorem 3.7** (Stability of the variational matrix methods).**

The methods $(\mathcal{M}^{(1)})^{\ast}$ \tagform@3.3 and $(\mathcal{M}^{(2)})^{\ast}$ \tagform@3.4 are stable.

Proof.

Follows as special case from the proof of the tensor version, Theorem 5.14. ∎

Despite the technicalities involved in the proof, the simplicity of the idea becomes apparent by setting $P=\mathcal{I}$ (being analogous to the argumentation in section 2.1). For a certain constant $c\in\mathbb{R}$, in the setting of Theorem 3.5, we then have

[TABLE]

If now $\sigma_{r}\rightarrow 0$, then also $(1+c\sigma_{r}^{-2})^{-1}\rightarrow 0$ and we obtain the same result as if we would have truncated the representation $(X,Y)=(U,\Sigma V^{T})$ to rank $r-1$ beforehand. We also denote these additional factors as filter, as they filter out influence corresponding to low singular values. Note that obtaining small singular values is not penalized, but using components corresponding to small ones is.

Algorithm 1 performs one stable ALS approximation (SALSA) sweep, that is it applies the two stable methods from Theorem 3.7. Before each update, any current singular value $\sigma_{i}<\sigma_{\mathrm{min}}$ is replaced by $\sigma_{\mathrm{min}}$, which in turn is set as fraction $f_{\min}\ll 1$ of the current residual (Algorithm 2). Although the influence of this manipulation on the subsequent step is thereby marginal, it is necessary since otherwise $\sigma_{i}$ may irreversibly converge to zero (cf. Remark 2.4, for more details, see Section 6).

Rank Adaption

The key aspect of stability is that it rendered explicit rank adaption near unnecessary, since the optimization relies on the magnitude of singular values, as \tagform@3.7 suggests. Starting from an initial representation and a large value $\omega$ proportional to the norm of the iterate $A$, the parameter $\omega$ is decreased after each iteration. The singular values are decided into two types:

Definition 3.8** (Stabilized rank and minor singular values).**

A singular value $\sigma_{i}$ is called stabilized, if it is larger than a certain fixed fraction of $\omega$ (meaning any corresponding terms have an increased influence, cf. \tagform@3.7). Otherwise, it is called minor (as a removal of such does not notably change the next steps). The stabilized rank only counts the number of stabilized singular values.

The actual rank is only modified as to make sure that there is always a fixed, small amount of minor singular values, i.e.

[TABLE]

for constants $f_{\mathrm{minor}}<1$ and $k_{\mathrm{minor}}\in\mathbb{N}$. This states the basic concept of implicit rank adaption and we will provide a more detailed discussion in the later section 8 for the elaborate tensor case. For performance evaluation, a validation set may be used:

Definition 3.9** (Validation set).**

For a given $P$, the sampling or training set, we define $P_{2}\subset P$ as validation set. This set may be chosen randomly or specifically distributed. The actual set used for the optimization is replaced by $P\leftarrow P\setminus P_{2}$ (keeping the same symbol).

The stopping criteria in Algorithm 2 may depend on the behavior of $P_{2}$, or may simply be based on a rank bound, e.g. $r\leq|P|/(n+m)$. The latter criterion, however, only suffices in the matrix case.

Remark 3.10** (On convergence estimates).**

*Due to the specific dependency of the regularization on the current singular values of the iterate, it may be impossible to restate the problem as minimization of a modified, global cost function. Furthermore, the iterate does not remain on a fixed manifold of low rank. Given that also rank deficient points may be included in the optimization due to the stability, theoretical convergence results so far remain subject to future work. *

Generalization to High-Dimensional Tensors

For the rest of this article we consider the problem of (approximately) reconstructing a tensor $M\in\mathbb{R}^{\mathcal{I}}$ from a given data set $M|_{P}=\{M_{p}\}_{p\in P}$, where now

[TABLE]

We further generalize the representation map and data space $\tau_{r}:\mathcal{D}_{r}\rightarrow\mathcal{D}_{r}$ as well as the rank $r$ together with the spaces $\mathcal{T}_{r}$ to the tensor train (TT-)format (Definition 4.2, [34, 33], also called Matrix Product States (MPS) [46, 41] or interpreted as a special case of the Hierarchical Tucker format [18, 13]). Stability for tensor methods is thereby defined through Definition 1.1 as well. For the underlying tensor $M$ it is now assumed that its TT-singular values (cf. Definition 4.1) decline sufficiently fast in order to yield an approximation $M_{\varepsilon}\in\mathcal{T}_{r}$.

Definition 4.1** (TT-singular values and TT-rank).**

Analogously to the TT-ranks, we define the TT-singular values $\sigma=(\sigma^{(1)},\ldots,\sigma^{(d-1)})$ of a tensor $A$ as the unique singular values of the corresponding matricizations $A^{(\{1,\ldots,\mu\})}\in\mathbb{R}^{(\mathcal{I}_{1}\times\cdots\times\mathcal{I}_{\mu})\times(\mathcal{I}_{\mu+1}\times\cdots\times\mathcal{I}_{d})}$ with entries

[TABLE]

of $A$, such that $\sigma^{(\mu)}$ contains the ordered singular values of $A^{(\{1,\ldots,\mu\})}$, $\mu=1,\ldots,d-1$. The TT-rank $r_{\mu}$ is the number of nonzero TT-singular values in $\sigma^{(\mu)}$. We also call $\sigma^{(\mu)}$ the $\mu$-th singular values.

Definition 4.2** (TT-format).**

*A tensor $A\in\mathbb{R}^{\mathcal{I}}$ is in the set $\mathcal{T}_{r}$, often also denoted $\mathrm{TT}(r)$, $r\in\mathbb{N}^{d-1}$, if for $\mu=1,\ldots,d$ and $i_{\mu}\in{\mathcal{I}}_{\mu}$ there exist $G_{\mu}(i_{\mu})\in\mathbb{R}^{r_{\mu-1}\times r_{\mu}}$  ($r_{0}=r_{d}=1$) such that *

[TABLE]

The representation of $A$ in this form is shortly called the TT format. The single $G_{\mu}$, as well as similarly structured objects, are called cores (or sometimes nodes).

The TT-SVD [33] provides that if $r$ is the TT-rank of $A$, then it follows $A\in\mathcal{T}_{r}$. So $G$ is also called low-rank representation of $A$. The function $\tau_{r}$ will appear in different forms, each denoting the mapping from a type of representation of rank $r$ to the represented object.

Definition 4.3** (Unfoldings).**

For a core $H$ (possibly a product of smaller cores in the TT representation) with $H(i)\in\mathbb{R}^{k_{1}\times k_{2}}$, $i=1,\ldots,n$, we denote the left unfolding $\mathfrak{L}(H)\in\mathbb{R}^{k_{1}\cdot n\times k_{2}}$, in which the matrices contained in the core are stacked below each other, and right unfolding $\mathfrak{R}(H)\in\mathbb{R}^{k_{1}\times k_{2}\cdot n}$, in which they are stacked side by side, by

[TABLE]

for $1\leq j\leq n$, $1\leq\ell\leq k_{1}$ and $1\leq q\leq k_{2}$.

Our goal is to determine the ranks adaptively. We will demonstrate why this can be even more troublesome than in the matrix case in the following.

For many tensors that stem from practical applications one observes exponentially decaying singular values, yet the rate of decay may strongly vary for the $d-1$ different matricizations. Theoretically, there is no non-trivial limitation to the shape of the TT-singular values as described in following lemma.

Lemma 4.4** (Feasibility of TT-singular values [25]).**

Let $\sigma=(\sigma^{(1)},\ldots,\sigma^{(d-1)})$ be a $(d-1)$-tuple of weakly decreasing $r_{\mu}$-tuples, $\mu=1,\ldots,d-1$, for which $\|\sigma^{(i)}\|_{2}=\|\sigma^{(j)}\|_{2}$ for all $i,j=1,\ldots,d-1$. Then there exist mode sizes $n_{1},\ldots,n_{d}\in\mathbb{N}$ and a tensor $A\in\mathbb{R}^{\mathcal{I}}$ such that this tensor has TT-singular values $\sigma$.

The following exemplary tensor emphasizes the problems heuristics encounter in the rank adaption of a tensor with inconvenient singular values.

Example 4.5** (Rank adaption test tensor).**

For $k\in\mathbb{N}$, let $Q\in\mathbb{R}^{n_{1}\times\ldots\times n_{4}}$ be an orthogonally decomposable $4$-dimensional Tensor with TT-rank $(k,k,k)$ and uniform singular values $\sigma^{(1)}=\sigma^{(2)}=\sigma^{(3)}=(\alpha,\alpha,\ldots)$ as well as $B\in\mathbb{R}^{n_{5}\times n_{6}}$ be a rank $2k$ matrix with exponentially decaying singular values $\sigma^{(5)}\propto(\beta^{-1},\beta^{-2},\ldots)$ for some $\alpha,\beta>0$. Then the separable tensor $M\in\mathbb{R}^{n_{1}\times\ldots\times n_{6}}$ defined by $M(i)=Q(i_{1},\ldots,i_{4})\cdot B(i_{5},i_{6})$ has singular values $\sigma$ and rank $r^{(M)}=(k,k,k,1,2k)$.

By its definition, $M$ is separable into a $4$- and a $2$-dimensional tensor ($Q$, $B$). Knowing this would of course drastically simplify the problem. We now consider the performance of two very basic rank adaption ideas.

Greedy, single rank increase: We test for maximal improvement by increase of one of the ranks $r_{\mu}$ ($\mu=1,\ldots,d-1$) of the iterate starting from $r\equiv 1$. Solely increasing either of $r_{2}$, $r_{3}$ or $r_{4}$ will give close to no improvement. As further shown in [11], the approximation of orthogonally decomposable tensors with lower rank can be problematic. In numerical tests (cf. Section 9.8), we can observe that $r_{5}$ is often increased to a maximum first. Thereby, extremely small singular values are involved that lie far beneath the current approximation error, although the rank is not actually overestimated. 2. 2.

Uniform rank increase and coarsening: We increase every rank $r_{\mu}$ ($\mu=1,\ldots,d-1$) starting from $r\equiv 1$ and decrease ranks when the corresponding singular values are below a threshold. The problem with this strategy is quite plain, namely that for the target tensor $M$, it holds $r^{(M)}_{5}=1$. If this rank is overestimated, the observed sampling points will be misinterpreted (overfitting) and it does not matter how small corresponding singular values become (see Lemma 2.1).

These indicated difficulties gain more importance with high dimension, but for one microstep at a time, can be resolved by regarding only three components of a tensor. We describe this in the following Section 4.1.

The tensor from Example 4.5 can also be constructed explicitly. We define a representation $G$ for $A=\tau_{r}(G)$ via left and right unfoldings by

[TABLE]

for (column-) orthogonal matrices $Q_{1}\in\mathbb{R}^{n_{1}\times r_{1}},\ Q_{4}\in\mathbb{R}^{n_{4}r_{4}\times r_{3}},\ Q_{5}\in\mathbb{R}^{r_{4}n_{5}\times r_{5}}\ Q_{6}\in\mathbb{R}^{n_{6}\times r_{5}}$ and $(\sigma_{5})_{i}\propto\beta^{-i}$, $\beta>1$. This tensor has exactly the properties postulated in the example.

Notations and Reduction to Three Dimensions

For necessary simplicity, we reduce the $d$ dimensional setting to a three dimensional one:

Definition 4.6** (Interface matrices).**

Let $\mu\in\{1,\ldots,d\}$. For a representation $G$, we define the left interface matrix $G^{<\mu}\in\mathbb{R}^{n_{1}\ldots n_{\mu-1}\times r_{\mu-1}}$ (cf. Definitions 4.2) via

[TABLE]

as well as the right interface matrix $G^{>\mu}\in\mathbb{R}^{r_{\mu}\times n_{\mu+1}\ldots n_{d}}$ via

[TABLE]

We further define the core $A_{(\mu)}\in\left(\mathbb{R}^{n_{1}\ldots n_{\mu-1}\times n_{\mu+1}\ldots n_{d}}\right)^{\mathcal{I}_{\mu}}$ as core unfolding with respect to mode $\mu$ of a tensor $A$ by

[TABLE]

For any representation it hence holds $(\tau_{r}(G))_{(\mu)}(j)=G^{<\mu}\ G_{\mu}(j)\ G^{>\mu}$, $j=1,\ldots,n_{\mu}$. Multiplication of a core $H$ with a matrix $B$ yields again cores, $HB$ and $BH$, given by

[TABLE]

respectively, for all possible $j$. The previous notations then allow us to compactly write

[TABLE]

This relation is displayed in Figure 3.

The interface matrices equal left and right unfoldings, respectively:

[TABLE]

In terms of Definition 3.3, for $P=\{p^{(i)}\mid i=1,\ldots,|P|\}$ and the operation $(\cdot)_{(\mu)}$ \tagform@4.2, which is used to combine components $s=1,\ldots,\mu-1$ as well as $s=\mu+1,\ldots,d$, let

[TABLE]

Thereby, $A|_{P}$ contains the same entries as $(A_{(\mu)})|_{P_{(\mu)}}$. For the selection of one slice, $(\cdot)_{(\mu)}(j)$, we denote

[TABLE]

Likewise, the vectorization of an index set $S\subset\mathbb{R}^{n\times m}$ is defined by $\mbox{vec}(S)=\{s_{1}+n(s_{2}-1)\in\mathbb{R}\mid s\in\mathbb{R}^{2}\}$.

Without loss of generality we can restrict our consideration to three dimensional tensors that correspond to the left and right interface matrices as well as the respective intermediate cores (cf. Definition 4.3):

Remark 4.7** (Reduction to three dimensions).**

When $\noindent\mu\in\{1,\ldots,d\}$ is fixed, we will only use the short notations

• •

$(L,\ \mathbb{N},\ R)=(G^{<\mu},\ G_{\mu},\ G^{>\mu})$**

• •

$(n_{L},\ n_{N},\ n_{R})=(n_{1}\cdot\ldots\cdot n_{\mu-1},\ n_{\mu},\ n_{\mu+1}\cdot\ldots\cdot n_{d})$**

• •

$(\gamma,\ \theta)=(\sigma^{(\mu-1)},\ \sigma^{(\mu)})$* and $(\Gamma,\ \Theta)=(\Sigma^{(\mu-1)},\ \Sigma^{(\mu)})=(\operatorname{diag}(\sigma^{(\mu-1)}),\ \operatorname{diag}(\sigma^{(\mu)}))$*

• •

$(r_{\gamma},\ r_{\theta})=(r_{\mu-1},\ r_{\mu})$**

• •

$B=M_{(\mu)}$* and $S=P_{(\mu)}$*

Hence, the important variables are $L$ (left part), $R$ (right part), $N$ (new part), $B$ (right hand side) and $S$ (sampling). The microsteps $\mathcal{M}^{(1)},\ldots,\mathcal{M}^{(d)}$ of ALS for the tensor train format only change the respective $G_{\mu}$ and are given by

[TABLE]

or equivalently $G^{+}_{\mu}(j)=\mathop{\rm argmin}_{{\widetilde{N}}(j)}\|L\cdot{\widetilde{N}}(j)\cdot R-B(j)\|_{S(j)}$ — an equation in which only matrices are involved. We only need to consider three-dimensional tensors $A\in\mathbb{R}^{n_{L}\times n_{N}\times n_{R}}$, with TT-rank $(r_{\gamma},r_{\theta})$ and TT-singular values $(\gamma,\theta)$. For simplicity, we redefine $\tau_{r}$ for this case via $A=\tau_{r}(L,\mathbb{N},R)$.

Stable Alternating Least Squares Microsteps for Tensor Completion

With the derivation in Section 4.1, we have seen that it is sufficient to only consider ALS for order $3$ tensors.

Definition 5.1** (Variational residual function (cf. Definition 3.1)).**

Let $\omega\geq 0,\ s_{1},s_{2}>0$ and $B,S,L,\mathbb{N},R$ as in Definition 4.7. We define the variational residual function $C:=C_{B,S,L,\mathbb{N},R}$ for $\mathbb{V}_{\omega}:=\mathbb{V}_{\omega}(L,\mathbb{N},R)$ by

[TABLE]

where $s_{1},s_{2}$ are scalings that only depend on the proportions of the representation, to be specified by Lemma 7.1.

As in the matrix case, $\Delta\mathbb{N}$ does not influence the minimizer, so we omit it from now on. It should further be noted that $\mathbb{V}_{\omega}$ does not depend on the unknown $\widetilde{\mathbb{N}}$. The scalings $s_{1}$, $s_{2}$ will become relevant in Section 7 in order to achieve a similar effect as in Remark 3.6.

Standard Representation of a TT-Tensor

A representation $G=(L,\mathbb{N},R)$ can be changed without changing the generated tensor $A=\tau_{r}(G)$ [37, 22], more specifically

[TABLE]

for two regular matrices $T_{1}\in\mathbb{R}^{r_{\gamma}\times r_{\gamma}},T_{2}\in\mathbb{R}^{r_{\theta}\times r_{\theta}}$. One can define an extended standard representation that explicitly contains the unique TT-singular values $(\gamma,\theta)$ which is essentially unique (in terms of uniqueness of the truncated matrix SVD333Both $U\Sigma V^{T}$ and $\widetilde{U}\Sigma\widetilde{V}^{T}$ are truncated SVDs of $A$ if and only if there exists an orthogonal matrix $W$ that commutes with $\Sigma$ and for which $\widetilde{U}=UW$ and $\widetilde{V}=VW$. For any subset of pairwise distinct nonzero singular values, the corresponding submatrix of $W$ needs to be diagonal with entries in $\{-1,1\}$.). For the construction, a slightly modified TT-SVD [33] is used444Although called TT-SVD, unlike the matrix SVD, decompositions constructed by the algorithm do not explicitly contain the singular values.. An analogous decomposition also appeared earlier in [42] and is, within the quantum computing community, sometimes referred to as canonical MPS form (not to be confused with canonical polyadic decomposition). This normalization is needed to obtain the same simplification as in Theorem 3.5.

Lemma 5.2** (Standard representation).**

*Let $A\in\mathbb{R}^{n_{L}\times n_{N}\times n_{R}}$ be a tensor.

There exists an essentially unique (extended) representation*

[TABLE]

for which $A=\tau_{r}(\mathcal{L},\ \Gamma\mathcal{N}\Theta,\ \mathcal{R})$ as well as $\mathcal{L}\ \Gamma\ \mathfrak{R}(\mathcal{N}\Theta\mathcal{R})$ and $\mathfrak{L}(\mathcal{L}\Gamma\mathcal{N})\ \Theta\ \mathcal{R}$ are (truncated) SVDs of $A^{(\{1\})}$ and $A^{(\{1,2\})}$, respectively. This in turn implies that $\mathcal{L}$ and $\mathfrak{L}(\Gamma\mathcal{N})$ are column orthogonal, as well as $\mathcal{R}$ and $\mathfrak{R}(\mathcal{N}\Theta)$ are row orthogonal.

Proof.

$1$. Uniqueness: Let there be two such representations $\widetilde{\mathcal{G}}$ and $\mathcal{G}$. Since the left-singular vectors of $A^{(\{1\})}$ are essentially unique, we conclude $\widetilde{\mathcal{L}}=\mathcal{L}W_{1}$ for an orthogonal matrix $W_{1}$ that commutes with $\Gamma$. Via an SVD of $A^{(\{1,2\})}$ it follows that $\widetilde{\mathcal{R}}=W_{2}^{T}\mathcal{R}$ for an orthogonal matrix $W_{2}$ that commutes with $\Theta$. Furthermore $\mathfrak{L}(\mathcal{L}\Gamma W_{1}\widetilde{\mathcal{N}})=\mathfrak{L}(\widetilde{\mathcal{L}}\Gamma\widetilde{\mathcal{N}})=\mathfrak{L}(\mathcal{L}\Gamma\mathcal{N})W_{2}$. The map $x\mapsto\mathfrak{L}(\mathcal{L}\Gamma x)$ is linear and, in this case, of full rank. This implies $\widetilde{\mathcal{N}}=W_{1}^{T}\mathcal{N}W_{2}$.

$2.$* Existence (constructive):* Let $A=\tau_{r}(\widetilde{L},\widetilde{\mathbb{N}},\widetilde{R})$ where $\mathfrak{R}(\widetilde{\mathbb{N}})$ and $\widetilde{R}$ are column orthogonal (this can always be achieved using \tagform@5.2). An SVD of $\widetilde{L}$ yields $\widetilde{L}=\mathcal{L}\ \Gamma\ V_{1}^{T}$, since $\mathcal{L}\ \Gamma\ \mathfrak{R}(V_{1}^{T}\widetilde{\mathbb{N}}\widetilde{R})$ is a truncated SVD of $A^{(\{1\})}$. A subsequent SVD of $\mathfrak{L}(\Gamma V_{1}^{T}\widetilde{\mathbb{N}})$ yields $\Gamma V_{1}^{T}\widetilde{\mathbb{N}}=\widehat{\mathbb{N}}\ \Theta\ V_{2}^{T}$, since $\mathfrak{L}(\mathcal{L}\widehat{\mathbb{N}})\ \Theta\ (V_{2}^{T}\widetilde{R})$ is a truncated SVD of $A^{(\{1,2\})}$. We can finish the proof defining $\mathcal{N}:=\Gamma^{-1}\widehat{\mathbb{N}}$ and $\mathcal{R}=V_{2}^{T}\widetilde{R}$. Note that, by construction, $\mathfrak{L}(\Gamma\mathcal{N})$ is column-orthogonal.

$3.$* Implied orthogonality:* Using the essential uniqueness, it follows that $\mathfrak{L}(\Gamma\mathcal{N})$ must indeed be column-orthogonal. By analogously constructing the extended representation from right to left we would obtain that $\mathfrak{R}(\mathcal{N}\Theta)$ is row-orthogonal. By uniqueness it follows again that this is always the case. ∎

Remark 5.3** (Conventional form of standard representation).**

Throughout the rest of the article, the standard representation will mostly appear in form of a specific, conventional representation

[TABLE]

hence with interface matrices $\mathcal{L}$ and $\mathcal{R}$ given by corresponding singular vectors.

Minimizer of the Variational Residual Function

We define (from now on) our method as

[TABLE]

with $C=C_{B,S,L,\mathbb{N},R}$ as in \tagform@5.1. Although Theorem 5.6, or more specifically the regularity of $Y(j)$ given by \tagform@5.10, later provide the uniqueness of the minimizer, we up to that point formally use the minimization of $\|\tau_{r}(L,\widetilde{\mathbb{N}},R)\|_{F}$ as secondary and representation independent criterion. The special cases $\mu\in\{1,d\}$ can be derived from the general case (Remark 7.2) and comply with the matrix case.

Lemma 5.4** (Representation independent).**

The method $\mathcal{M}^{\ast}$ is representation independent.

Proof.

Let $\mathbb{N}^{+}:=\mathop{\rm argmin}_{\widetilde{N}}C$, $C=C_{B,S,L,\mathbb{N},R}(\widetilde{N})$ and $\widehat{\mathbb{N}}^{+}:=\mathop{\rm argmin}_{\widetilde{N}}\widehat{C}$,

$\widehat{C}=C_{B,S,\widehat{L},\widehat{\mathbb{N}},\widehat{R}}(\widetilde{N})$ for representations $\tau_{r}(L,\mathbb{N},R)=\tau_{r}(\widehat{L},\widehat{\mathbb{N}},\widehat{R})$ as well as $\widehat{\mathbb{V}}_{\omega}=\mathbb{V}_{\omega}(\widehat{L},\widehat{\mathbb{N}},\widehat{R})$ and $\mathbb{V}_{\omega}=\mathbb{V}_{\omega}(L,\mathbb{N},R)$. According to \tagform@5.2, there exist two matrices $T_{1},T_{2}$ such that

[TABLE]

Hence

[TABLE]

The substitution $(\Delta\widehat{L},\Delta\widehat{\mathbb{N}},\Delta\widehat{R})\overset{\iota}{\rightarrow}(\Delta LT_{1},T_{1}^{-1}\Delta\mathbb{N}T_{2},T_{2}^{-1}\Delta R)$ introduces a constant Jacobi Determinant $|\det(J_{\iota})|=1$. We obtain

[TABLE]

Hence $\widehat{\mathbb{N}}^{+}=T_{1}^{-1}\mathbb{N}^{+}T_{2}$. This is the same relation given for $\mathbb{N}$ and $\widehat{\mathbb{N}}$ and therefore $\tau_{r}(L,\mathbb{N}^{+},R)=\tau_{r}(\widehat{L},\widehat{\mathbb{N}}^{+},\widehat{R})$ (which is a set equality if the minimizer is not assumed to be unique). ∎

Corollary 5.5** (Integral over Kronecker product).**

*Let $\omega_{1}>0$.

Further, let $H\in\mathbb{R}^{(n_{X}n_{Y})\times(n_{X}n_{Y})}$ as well as $Y\in\mathbb{R}^{n_{Y}\times n_{Y}}$ be matrices and*

[TABLE]

Then

[TABLE]

for $H^{\ast}=\mathrm{tr}_{1}(H)\in\mathbb{R}^{n_{Y}\times n_{Y}}$.555$\mathrm{tr}_{1}$ is the partial trace, i.e. $(H^{\ast})_{i,j}=\mbox{tr}(h_{i,j}),\ H=\sum_{i,j}h_{i,j}\otimes e_{i}e_{j}^{T},\ h_{i,j}\in\mathbb{R}^{n_{X}\times n_{X}}$ For an analog $V_{\omega_{2}}^{(n_{Y},m_{Y})}$, $\omega_{2}>0$, we further have

[TABLE]

Proof.

Using the splitting $H=\sum_{i,j}h_{i,j}\otimes e_{i}e_{j}^{T}$, $h_{i,j}\in\mathbb{R}^{n_{X}\times n_{X}}$, Lemma 3.4 can be applied to each summand, separately for $X$ and $Y$. ∎

We now derive the minimizer of the variational residual function \tagform@5.1. Due to Lemma 5.4, we can use the standard representation in form of Remark 5.3 for simplification. In this case, $\mathbb{V}_{\omega}$ takes the convenient form

[TABLE]

Theorem 5.6** (Minimizer of the ALS variational residual function).**

*Let $(\mathcal{L},\ \Gamma,\ \mathcal{N},\ \Theta,\ \mathcal{R})$ be the standard representation \tagform@5.3 for a tensor $A$. The minimizer $\mathbb{N}^{+}$ of the residual function $C_{B,S,\mathcal{L},\Gamma\mathcal{N}\Theta,R}$ as in \tagform@5.1 is given by *

[TABLE]

with $S(j)_{u}=(x^{(1)}_{u},x^{(2)}_{u},\ldots)$, $u=1,2$, for $S(j)=\{x^{(1)},x^{(2)},\ldots,x^{(|S(j)|)}\}\subset\mathbb{N}^{2}$. The constants $\zeta\noindent$ only depend on the proportions of the representation and sampling set (cf. Remark 5.7) as well as the constant scalings $s_{1}$, $s_{2}$.

The regularization term is more straightforward than it might appear. The computational complexity for the calculation of the minimizer is of the same order (with near same constant) as for standard ALS, for which the matrices $\mathcal{L}_{S(j)_{1},:}\in\mathbb{R}^{\noindent|S(j)|\times r_{\gamma}}$ and $\mathcal{R}_{:,S(j)_{2}}\in\mathbb{R}^{r_{\theta}\times\noindent|S(j)|}$ are required anyhow (for further explanation, see \tagform@7.2,\tagform@7.3). The scalings $n_{L}^{-1},n_{R}^{-1}$ go in hand with these to adjust the regularization terms to the magnitudes of the corresponding shares of the sampling set and disappear for the approximation of a fully available tensor (Corollary 5.8). The constants $\zeta$ are weighted according to the proportions of the left and right interface matrices. For example, for $\mu=1$ (the position of the core being updated), we have $\zeta_{1}=0$ and accordingly for $\mu=2$, $\zeta_{1}\ll\zeta_{2}$. In practice, instead of evaluating the minimizer exactly, a few steps (that is for very coarse tolerance) of a suitably computed, preconditioned cg are performed as described in Section 5.3. This can be achieved without changing results.

Proof.

We omit the scalings $s_{1},s_{2}$ for simplicity since they only have to be carried along the lines. We search for $\mathbb{N}^{+}:=\mathop{\rm argmin}_{\widetilde{\mathbb{N}}}C_{B,S,\mathcal{L},\Gamma\mathcal{N}\Theta,\mathcal{N}}(\widetilde{\mathbb{N}})$. Substituting

[TABLE]

we can (up to a constant factor) restate $C$ as

[TABLE]

Each of the independent matrices of the minimizing core is restated as

[TABLE]

where $\otimes_{K}$ is the matrix Kronecker product. Let $j$ be arbitrary but fixed from now on. For any $x$, it is $\|x\|_{\mbox{vec}(S(j))}=\|H(j)x\|_{F}=x^{T}H(j)x$ for a diagonal, square matrix $H(j)\in\mathbb{R}^{|\mathcal{I}|/n_{N}\times|\mathcal{I}|/n_{N}}$ with $H(j)_{(s),(s)}=\delta_{s\in S(j)}$ (hence $H(j)^{2}=H(j)$). Using the normal equation, we obtain $\mathbb{N}^{+}(j)=Y^{-1}b$, where

[TABLE]

and

[TABLE]

In both $Y$ and $b$, any perturbation that appears only one-sided vanishes due to symmetry of $\mathbb{V}_{\omega}$. Hence $b=|\mathbb{V}_{\omega}|\ {({\mathcal{R}}^{T}\otimes_{K}\mathcal{L})_{\mbox{vec}(S(j)),:}}^{T}\ \mbox{vec}(B(j))_{\mbox{vec}(S(j))}$ and for $\mathop{}\!\mathrm{d}\delta:=\mathop{}\!\mathrm{d}\Delta\mathcal{L}\mathop{}\!\mathrm{d}\Delta\mathcal{N}\mathop{}\!\mathrm{d}\Delta\mathcal{R}$

[TABLE]

Now, let $\ell=\#_{\mathcal{R}},\ n=\#_{\mathcal{N}},\ k=\#_{\mathcal{L}}$. Since $\mathbb{V}$ is a version of the $(\ell+n+k)$-sphere, we can use the following integration formula: Let $f:\mathbb{R}^{n+m}\rightarrow\mathbb{R}^{k}$ be a sufficiently smooth function and $S^{v-1}_{\omega}$ be the $v$-sphere of radius $\omega$. Then

[TABLE]

We use it twice and thereby split the integral. For a function $f$ we then obtain

[TABLE]

If $f$ is independent of $\Delta\mathcal{N}$, this then simplifies to

[TABLE]

We further use the identity (where the function $\Gamma(\cdot)$ is not to be confused with the given diagonal matrix $\Gamma$)

[TABLE]

We apply these and Corollary 5.5 for different $f=(X\otimes_{K}Y)^{T}H(j)(X\otimes_{K}Y)$. For $\delta_{1},\delta_{2}\in\{0,1\}$ we set $X$ as ${\mathcal{R}}^{T}$ ($\delta_{1}=0$) or $\Delta{\mathcal{R}}^{T}\Theta^{-1}$ ($\delta_{1}=1$) and analogously $Y$ as $\mathcal{L}$ ($\delta_{2}=0$) or $\Delta\mathcal{L}\Gamma^{-1}$ ($\delta_{2}=1$). For the summands $Y(0,0)+Y(1,0)+Y(0,1)+Y(1,1)=Y$ this then yields

[TABLE]

for $c=\frac{8\pi^{(n+k+\ell)/2}}{\Gamma(n/2)\Gamma(\ell/2)\Gamma(k/2)}$. The constant matrices $C_{H}$ are given by

[TABLE]

Furthermore, it is $|\mathbb{V}_{\omega}|=c\frac{\omega^{n+\ell+k}}{n+\ell+k}\nu(\ell+k-1,n-1)\ \nu(k-1,\ell-1)$. Factoring out this base volume in $Y=|\mathbb{V}_{\omega}|\widetilde{Y}$ by using properties of the $\Gamma$ function, one derives:

[TABLE]

where the constants $s_{1}$ and $s_{2}$ have been added again. Restating the result again as a least squares problem finishes the proof. ∎

Remark 5.7** (Specification of constants).**

*Let $\#_{\mathcal{R}}:=size(\mathcal{R}),\ \#_{\mathcal{N}}:=size(\mathcal{N}),\ \#_{\mathcal{L}}:=size(\mathcal{L})$ be the sizes of the tensor components. The constants in Theorem 5.6 are given by *

[TABLE]

Evaluation with Coarse Conjugate Gradient

For each slice $j$, the solution to the least squares problem in Theorem 5.6 is described through the normal equation

[TABLE]

with

[TABLE]

In practice, since within each microstep we do not benefit from an exact solution of this transitory system, it is much more economic to perform a preconditioned conjugate gradient method and terminate when a coarse, relative tolerance (e.g. $\mathrm{tol}=10^{-2}$) is reached. This tolerance is empirically chosen such that the number of required cg steps is minimized, however under the condition that the approximation quality does not notably change in either direction — such that neither loss of accuracy nor additional regularization can be observed (cf. Section 8). With the following consideration, we can construct a preconditioner (cf. \tagform@3.7).

Corollary 5.8** (Filter properties).**

For full sampling, $P=\mathcal{I}$, the update is given by the so called filter (a diagonal matrix)

[TABLE]

Proof.

From $P=\mathcal{I}$, it follows that $\mathcal{R}_{:,S(j)_{2}}$ is an $n_{L}$-order copy of $\mathcal{R}$ and $\mathcal{L}_{S(j)_{1},:}$ is an $n_{R}$-order copy of $\mathcal{L}$ (cf. \tagform@5.10). Hence the regularization terms are the same for all $j=1,\ldots,n_{N}$. The minimizer $N^{+}(j)$ is given by

[TABLE]

The factors $\sqrt{n_{L}^{-1}}$ and $\sqrt{n_{R}^{-1}}$ vanish in $Z(j)^{T}Z(j)$ due to the multiple rows involving the orthogonal matrices $\mathcal{R}$ and $\mathcal{L}$. Furthermore, $(Z(j)^{T}Z(j))^{-1}$ is diagonal. ∎

Remark 5.9** (Application of cg algorithm and order of computational complexity).**

The matrix $\mathcal{F}$ serves as excellent preconditioner in the sense that $n_{L}^{-1}n_{R}^{-1}|S(j)|\mathcal{F}^{-1}\approx Z(j)^{T}Z(j)$, which holds as equality for $P=\mathcal{I}$ (as in the previous Corollary 5.8). This relation can as well be quantified through upper bounds on the condition number of $Z(j)\mathcal{F}^{1/2}$ as in Lemma 5.15 by which we expect (and observe in practice) very few iterations (and at most $r_{\gamma}r_{\theta}$) to be sufficient to reach a given coarse tolerance (for standard ALS, this is however not necessarily the case, cf. Lemma 5.11). Each single cg step then has complexity $\mathcal{O}(r_{\gamma}r_{\theta}|P|)$, while the full least squares problem has complexity $\mathcal{O}(r_{\gamma}^{2}r_{\theta}^{2}|P|)$ (per slice).

Stability and Restricted Isometry Properties

Before we derive the central theoretical statement of this paper, Theorem 5.14, some preparation is necessary. The tensor restricted isometry property (e.g. [35]) does not hold for any non trivial sampling set $P\varsubsetneq\mathcal{I}$. We however only need to work with a modified version as follows, in which left and right interface matrices are fixed. Apart from that, the shape is exactly the same.

Definition 5.10** (Internal tensor restricted isometry property (iTRIP)).**

We say a rank $r$ tensor $A=\tau_{r}(L,\mathbb{N},R)$ has the internal tensor restricted isometry property for the sampling set $S=P_{(\mu)}$, if there exist $0\leq c<1$ and $\rho>0$ with

[TABLE]

for all $\widetilde{A}\in\mathcal{A}(L,R):=\{\tau_{r}(L,\widetilde{\mathbb{N}},R)\mid\widetilde{\mathbb{N}}\mbox{ arbitrary }\}$.

Given a tensor $A$, if the iTRIP does not hold, then we can not expect the next update to have good completion properties, since changes on the sampling subset are unrelated to changes on the whole space. As indicated in Example 1.2, if $c\rightarrow 1$, then the next iterate can be arbitrarily bad. Note that the constants only depend on the tensor $A$, and not on its representation, and that this property is easy to check. In particular (as proven by Lemma 5.15 for $\omega=0$), the iTRIP with constant $c$ is equivalent to

[TABLE]

where $(\mathcal{L},N,\mathcal{R})$ corresponds to the standard representation as in \tagform@5.4. Hence, the property holds as long as that matrix has full rank (where $\kappa_{2}$ is the condition number regarding the spectral norm $\|\cdot\|_{2}$) or equivalently $\mathbb{N}\mapsto(L\mathbb{N}R)_{S}$ is injective. Note that, by a slice wise consideration, the condition number of a single $(\mathcal{R}^{T}\otimes\mathcal{L})|_{\mathrm{vec}(S(j)),:}$ can be improved since different magnitudes of sampling for each slice can be compensated (after all, the slices are solved independently).

Lemma 5.11** (Likelihood of the iTRIP).**

Let $\mathcal{T}_{r}$ be the subset of $3$ dimensional tensors with rank $r=(r_{\gamma},r_{\theta})$. Let $P$ be a (random) sampling that fulfills $|S(j)|\geq r_{\gamma}r_{\theta}$ for all $j=1,\ldots,n_{N}$. Then, for (only) almost every representation $(L,N,R)\in\mathbb{R}^{n_{L}\times r_{\gamma}}\times\mathbb{R}^{r_{\gamma}\times n_{N}\times r_{\theta}}\times\mathbb{R}^{r_{\theta}\times n_{R}}$ (with respect to the Lebesgue measure), the tensor $A=\tau_{r}(L,N,R)\in\mathcal{T}_{r}$ has the iTRIP. If for one $j$, $|S(j)|<r_{\gamma}r_{\theta}$, then no $A\in\mathcal{T}_{r}$ has the iTRIP.

Proof.

A tensor $A=\tau_{r}(L,\mathbb{N},R)$ has the iTRIP (for some valid constants) if and only if the linear map $\mathbb{N}\mapsto(L\mathbb{N}R)_{S}$ is injective, or equivalently, $(R^{T}\otimes L)_{\mbox{vec}(S(j)),:}$ has full rank for each $j$. Due to the provided slice density of $P$, each matrix $(\mathcal{R}^{T}\otimes\mathcal{L})_{\mbox{vec}(S(j)),:}$ is of size $|S(j)|\times r_{\gamma}r_{\theta}$. Hence generically, it is of full rank. If $|S(j)|<r_{\gamma}r_{\theta}$, then the matrix cannot have full rank. ∎

Tensors themselves that do not have the iTRIP, assuming sufficient sampling, pose just a marginal phenomenon for high dimension $d$.

A quite simple construction however shows that the iTRIP does not behave well under perturbation:

Lemma 5.12** (iTRIP under perturbation).**

Let $B\in\mathbb{R}^{n_{L}\times n_{N}\times n_{R}}$ with singular values $(\gamma^{(B)},\theta^{(B)})$. Assume further that for one $j$ it holds $|S(j)_{1}|<n_{L}$. Then for every $\sigma^{\ast}>0$, there exists a tensor $A$ with rank $(r_{\gamma},r_{\theta})$ and $\|A-B\|^{2}_{F}\leq\sum_{i=r_{\gamma}}^{\infty}(\gamma^{(B)}_{i})^{2}+2\sum_{i=r_{\theta}+1}^{\infty}(\theta^{(B)}_{i})^{2}+(\sigma^{\ast})^{2}$ such that $A$ does not have the iTRIP (and $\gamma^{(A)}_{r_{\gamma}}=\sigma^{\ast}$). If $B$ already has rank $(r_{\gamma},r_{\theta})$, then $\|A-B\|^{2}_{F}\leq(\gamma^{(B)}_{r_{\gamma}})^{2}+(\sigma^{\ast})^{2}$ suffices.

Proof.

Truncation of $B$ yields a tensor $\widetilde{A}$ with rank $(\ast,r_{\theta})$ and $\|\widetilde{A}-B\|^{2}_{F}\leq\delta:=\sum_{i=r_{\theta}+1}^{\infty}(\theta^{(B)}_{i})^{2}$. The tensor $\widetilde{A}$ hence has perturbed singular values such that $\|\gamma^{(\widetilde{A})}-\gamma^{(B)}\|_{2}\leq\delta$ (Mirsky’s Theorem [30]). Let $(\widetilde{\mathcal{L}},\Gamma^{(\widetilde{A})},\widetilde{\mathcal{N}},\Theta^{(\widetilde{A})},\widetilde{\mathcal{R}})$ be the standard representation of $\widetilde{A}$. Without loss of generality, we may assume that $j=1$ and that only points in the first $k:=|S(1)|$ rows of $B(1)$ are contained in the sampling $S(1)$. Let now

[TABLE]

If $X$ is already singular, then we may choose $A=\widetilde{A}$. Otherwise, then we may choose $x=\alpha Xv$, $y=\alpha\widehat{y}$ for $\alpha=\|(Xv;\widehat{y})\|_{2}^{-1}$ and $v=-(X^{T}X)^{-1}Y^{T}\widehat{y}$, for an arbitrary vector $\widehat{y}\neq 0$. In all three cases, $\mathcal{L}$ is orthogonal and for $A:=\tau_{r}(\mathcal{L},\mathrm{diag}(\gamma^{(\widetilde{A})}_{1},\ldots,\gamma^{(\widetilde{A})}_{r_{\gamma}-1},\sigma^{\ast})\widetilde{\mathcal{N}}\Theta^{(\widetilde{A})},\widetilde{\mathcal{R}})$ it holds $\|A-B\|_{F}\leq\|A-\widetilde{A}\|_{F}+\|\widetilde{A}-B\|_{F}\leq\sum_{i=r_{\gamma}}^{\infty}(\sigma^{(\widetilde{A})}_{\gamma})_{i}^{2}+(\sigma^{\ast})^{2}+\delta\leq\sum_{i=r_{\gamma}}^{\infty}(\sigma^{(B)}_{\gamma})_{i}^{2}+\delta+(\sigma^{\ast})^{2}+\delta$. Yet $(\mathcal{R}^{T}\otimes\mathcal{L})|_{\mathrm{vec}(S(1)),:}$ is a singular matrix, since $\mathcal{L}_{\{1,\ldots,k\},:}$ is already singular. ∎

The statement analogously holds true for $r_{\theta}$ and can easily be transferred to matrix completion as well. Tensors that do not have the iTRIP are hence densely scattered depending on $\gamma^{(B)}_{r_{\gamma}}$, as are, more importantly, surroundings in which the constant $c$ is close to $1$ and overfitting becomes more likely (cf. Example 1.2).

In case of the regularized update, the additional term $Y(j)$ \tagform@5.10 does not allow the condition number to change that easily (cf. Lemma 5.16).

Lemma 5.13** (Partial matrix inverse by divergent parts).**

*We partition $\{1,\ldots,n\}=\omega_{j}\cup\omega_{j}^{c}$ $(\omega_{j}^{c}=\{1,\ldots,n\}\setminus\omega_{j})$, $j=1,2$ and define $\Omega:=\omega_{1}\times\omega_{2}$, $\widetilde{\Omega}:=\omega^{c}_{1}\times\omega^{c}_{2}$. Let $\{A^{(k)}\}_{k},\ \{J^{(k)}\}_{k}\subset\mathbb{R}^{n\times n}$ be series of symmetric matrices, $supp(J^{(k)})\subset\Omega$.

If $\mbox{lim}_{k\rightarrow\infty}A^{(k)}|_{\widetilde{\Omega}}=A|_{\widetilde{\Omega}}$, $A|_{\widetilde{\Omega}}$ s.p.d, and $\sigma_{\mathrm{min}}(J^{(k)}|_{\Omega})\rightarrow\infty$, then $V:=\mbox{lim}_{k\rightarrow\infty}(A^{(k)}+J^{(k)})^{-1}$ exists and we have $V|_{\widetilde{\Omega}}=(A|_{\widetilde{\Omega}})^{-1}$ and $V|_{\widetilde{\Omega}^{c}}=0$ $(\widetilde{\Omega}^{c}=\{1,\ldots,n\}^{2}\setminus\widetilde{\Omega})$.*

Proof.

First, w.l.o.g., let $\Omega=\{m+1,\ldots,n\}^{2}$. Otherwise we can apply permutations. Further, let $V^{(k)}:=A^{(k)}+J^{(k)}$. We partition our (symmetric) matrices $M$ for $M_{1,1}\in\mathbb{R}^{m\times m}$ block-wise as

[TABLE]

Note that $J^{(k)}_{1,1},J^{(k)}_{1,2}\equiv 0$. Since $A^{(k)}_{1,1}=V^{(k)}_{1,1}$ and $A_{1,1}=A|_{\widetilde{\Omega}}$ is s.p.d, $A^{(k)}_{1,1}$ is invertible for all $k>K$ for some $K$ and hence $\mbox{lim}_{k\rightarrow\infty}(V^{(k)}_{1,1})^{-1}=A_{1,1}^{-1}$. Further, $\sigma_{\mathrm{min}}(B_{2,2}^{(k)})>\sigma_{\mathrm{min}}(J_{2,2}^{(k)})-\sigma_{\mbox{max}}(A_{2,2}^{(k)})\rightarrow\infty$ and hence $\|(V^{(k)}_{2,2})^{-1}\|\rightarrow 0$. Therefore, for $k>\widetilde{K}$ and $H^{(k)}:=V^{(k)}_{1,1}-V^{(k)}_{1,2}(V^{(k)}_{2,2})^{-1}(V^{(k)}_{1,2})^{T}$, it is $\sigma_{\mathrm{min}}(H^{(k)})>\sigma_{\mathrm{min}}(A_{1,1})/2$. By block-wise inversion of $V^{(k)}$, it then follows $((V^{(k)})^{-1})_{1,1}=(H^{(k)})^{-1}\rightarrow(A^{(k)}_{1,1})^{-1}$. Similarly, $((V^{(k)})^{-1})|_{\Omega}\rightarrow 0$. ∎

One last step remains, since we cannot allow $\zeta$ to depend on the rank $r$. For now, we redefine the method $\mathcal{M}^{\ast}$ to directly yield the result in Theorem 5.6 for arbitrary constants $\zeta$, i.e.

[TABLE]

We explain in Section 7 and Lemma 7.1 how the scalings $s_{1},s_{2}$ as well as $\omega$ are used to obtain one specific $\mathcal{M}_{\zeta}^{\ast}$ from $\mathcal{M}^{\ast}$, for which $\zeta$ is indeed independent of $r$.

Theorem 5.14** (Stability of the method $\mathcal{M}_{\zeta}^{\ast}$).**

Let $B$ be the target tensor, $S\ \subsetneq\mathcal{I}$ the sampling set, arbitrary but fixed, and $\mathcal{M}_{\zeta}^{\ast}$ as in \tagform@5.16.

• •

The regularized method $\mathcal{M}_{\zeta}^{\ast}$ $(\omega>0)$ as defined by \tagform@5.16 (for $\zeta_{1},\zeta_{2}\geq 0$ and $\zeta_{(1,2)}>0$ that do not depend on $r$) is stable at all points $A^{\ast}$ (and hence also fixed-rank stable).

• •

The unregularized method $(\omega=0)$ \tagform@4.6 provides stability only for fixed rank (cf. Example 2.1), and only at those points $A^{\ast}$ that have the iTRIP (cf. Def. 5.10).

Proof.

Let $A^{\ast}$ be a fixed tensor with TT-ranks $r^{\ast}$.

1. Fixed-rank stability: We first show that $\mathcal{M}^{\ast}$ is stable for fixed rank. Let $A_{i}$ be a sequence with $\mbox{rank}(A_{i})=r^{\ast}$ and $A_{i}\rightarrow A^{\ast}$. Let $\mathcal{G}^{\ast}=(\mathcal{L}^{\ast},\Gamma^{\ast},\mathcal{N}^{\ast},\Theta^{\ast},\mathcal{R}^{\ast})$ be the standard representation of $A^{\ast}$ as well as $\mathcal{G}_{i}$ correspond to $A_{i}$. We partition the indices for $\gamma^{\ast}$ and $\theta^{\ast}$ by $k$ and $\ell$ according to equality of entries, such that $\gamma^{\ast}_{1}=\ldots=\gamma^{\ast}_{k_{1}}>\gamma^{\ast}_{k_{1}+1}=\ldots=\gamma^{\ast}_{k_{2}}>\ldots>\gamma^{\ast}_{k_{K-1}+1}=\ldots=\gamma^{\ast}_{k_{K}}>0$ and likewise for $\ell_{1},\ldots,\ell_{L}$. Since $A_{i}\rightarrow A^{\ast}$, their singular values also converge (e.g. [45]). We can hence conclude from [10, 43] that there exist sequences of block diagonal, orthogonal matrices $W_{i}$ and $M_{i}$ with block sizes $k_{1},k_{2}-k_{1},\ldots,k_{K}-k_{K-1}$ and $\ell_{1},\ell_{2}-\ell_{1},\ldots,\ell_{L}-\ell_{L-1}$, respectively, such that

[TABLE]

since the standard representation includes left and right singular vectors. We have to show that the tensors $H_{i}=\tau_{r}(\mathcal{L}_{i},\mathbb{N}_{i},\mathcal{R}_{i})=\tau_{r}(\mathcal{M}^{\ast}(\mathcal{L}_{i},\Gamma_{i}\mathcal{N}_{i}\Theta_{i},\mathcal{R}_{i}))$ converge to the analogously defined $H^{\ast}$. For fixed $j$, we define for each single $\mathcal{G}_{i}$ the matrix $Y_{i}=Y(j)$ (cf. Theorem 5.6, Remark 5.14) and $z_{i}:=(\mathcal{R}_{i}^{T}\otimes\mathcal{L}_{i})$ such that

[TABLE]

We define the shifted matrices

[TABLE]

Due to \tagform@5.17, it holds $(z^{M,W}_{i})_{\mbox{vec}(S(j)),:}\rightarrow z^{\ast}_{\mbox{vec}(S(j)),:}$. Inserting $I=(M_{i}^{T}\otimes W_{i})(M_{i}^{T}\otimes W_{i})^{T}$ into \tagform@5.18, we obtain

[TABLE]

Since $W_{i}^{T}\Gamma^{\ast}W_{i}=\Gamma^{\ast}$ for all $i$, it follows $W_{i}^{T}\Gamma_{i}W_{i}\rightarrow\Gamma^{\ast}$. Likewise $M_{i}^{T}\Theta_{i}M_{i}\rightarrow\Theta^{\ast}$ and thereby also ${Y^{M,W}_{i}}^{T}\ Y^{M,W}_{i}\rightarrow{{Y}^{\ast}}^{T}\ {Y}^{\ast}$. We treat the cases $\omega=0$ and $\omega>0$ separately:

(i) $\omega=0$: In this case, $Y^{M,W}_{i}=0=Y^{\ast}$. If the iTRIP holds for $A^{\ast}$, then

$\sigma_{\mathrm{min}}(z^{\ast}_{\mbox{vec}(S(j)),:})>0$ and therefore

[TABLE]

This directly yields convergence of $(H_{i}(j))\rightarrow(H^{\ast}(j))$ since all involved factors converge.

(ii) $\omega>0$: Here, we use that $\sigma_{\mathrm{min}}(Y^{\ast})>0$ and $\sigma_{\mathrm{min}}(z^{\ast}_{\mbox{vec}(S(j)),:})\geq 0$. We then obtain convergence since

[TABLE]

This proves fixed-rank stability.

2. Stability: Let now $A_{i}$ have arbitrary ranks. Without loss of generality by consideration of a finite amount of infinite subsequences, we can assume that $\mbox{rank}(A_{i})\equiv r$ for all $i$. Then, since $TT(r^{\ast})$ is a manifold, it follows $\gamma\geq\gamma^{\ast}$ and $\theta\geq\theta^{\ast}$. We can therefore have singular values $(\gamma_{i})_{k_{K}+1},\ldots,(\gamma_{i})_{k_{K+1}}\rightarrow 0$ as well as $(\theta_{i})_{\ell_{L}+1},\ldots,(\theta_{i})_{\ell_{L+1}}\rightarrow 0$. We expand the matrices $W_{i}$ and $M_{i}$ by identities of appropriate sizes to account for the vanishing singular values: $W_{i}\leftarrow diag(W_{i},I_{k_{K+1}-k_{K}})$, $M_{i}\leftarrow diag(M_{i},I_{\ell_{L+1}-\ell_{L}})$. In regard of Proposition 5.13, let $\Omega$ be the smallest cross product set, such that $({Y^{M,W}_{i}}^{T}\ Y^{M,W}_{i})|_{\widetilde{\Omega}}$ converges (which is the set that corresponds to vanishing singular values). Then, due to the definition of $Y^{M,W}_{i}$, $\sigma_{\mathrm{min}}(({Y^{M,W}_{i}}^{T}\ Y^{M,W}_{i})|_{\Omega})\rightarrow\infty$. We can conclude that

[TABLE]

and

[TABLE]

Because of this restriction, we in turn again get convergence to the limit $(H_{i}(j))\rightarrow(H^{\ast}(j))$, since all parts that correspond to vanishing singular values, also vanish within the update. This finishes the proof. ∎

Definition 5.15** (Stabilized internal tensor restricted isometry property (siTRIP)).**

We say a rank $r$ tensor $A=\tau_{r}(L,\mathbb{N},R)$ has the stable internal tensor restricted isometry property for the sampling set $S$, if there exist $0\leq c<1$ and $\rho>0$ such that for all $\widetilde{N}$ holds

[TABLE]

where $\widetilde{A}_{\Delta}:=\tau_{r}(L+\Delta L,\widetilde{\mathbb{N}}+\Delta\mathbb{N},R+\Delta R)|^{2}_{F}$ and $\mathop{}\!\mathrm{d}\Delta=\mathop{}\!\mathrm{d}\Delta L\mathop{}\!\mathrm{d}\Delta\mathbb{N}\mathop{}\!\mathrm{d}\Delta R$.

The constants are independent of the specific representation (cf. Lemma 5.4).

Lemma 5.16**.**

Let $Z(j)=\begin{pmatrix}(\mathcal{R}^{T}\otimes\mathcal{L})|_{\mathrm{vec}(S(j)),:}\\ Y(j)\end{pmatrix}$ as in \tagform@5.14. The siTRIP with constant $c$ for $A\in\mathbb{R}^{\mathcal{I}}$ is equivalent to

[TABLE]

*where $A=\tau_{r}(\mathcal{L},\mathbb{N},\mathcal{R})$ is its standard representation and $\mathcal{F}$ is as in Lemma 5.8. *

Proof.

Let $Z=\operatorname{diag}(Z(1),\ldots,Z(n_{N}))$ and $\mathrm{vec}(\widetilde{\mathbb{N}})^{T}=(\mathrm{vec}(N(1))^{T},\ldots,\mathrm{vec}(N(n_{N}))^{T})$. We abbreviate the siTRIP \tagform@5.19 as $(1-c)\beta\leq\rho\xi\leq(1+c)\beta$ ($\beta=\beta(\widetilde{\mathbb{N}}),\ \xi=\xi(\widetilde{\mathbb{N}})$). Let $a:=\max_{x\neq 0}\frac{\|Zx\|^{2}}{\|(I\otimes\mathcal{F}^{-1/2})x\|^{2}}$ and $b:=\min_{x\neq 0}\frac{\|Zx\|^{2}}{\|(I\otimes\mathcal{F}^{-1/2})x\|^{2}}$. Since the perturbation $\Delta\mathbb{N}$ is independent of $\widetilde{N}$, this term can be neglected in consideration of that $\|\widetilde{N}\|$ is not bounded. With $P=\mathcal{I}$ it holds $\beta=\sum_{j=1}^{n_{N}}\mathrm{vec}(\widetilde{\mathbb{N}}(j))^{T}\mathcal{F}^{-1}\mathrm{vec}(\widetilde{\mathbb{N}}(j))=\|(I\otimes\mathcal{F}^{-1/2})\mathrm{vec}(\widetilde{\mathbb{N}})\|^{2}_{F}$ (cf. proof of Theorem 5.6). For the actual sampling $P$, we have $\xi=\|Z\mathrm{vec}(\widetilde{\mathbb{N}})\|^{2}_{2}$. Thereby $a=\max_{\widetilde{\mathbb{N}}\mid\beta=1}\xi$ and $b=\min_{\widetilde{\mathbb{N}}\mid\beta=1}\xi$. Now, given the siTRIP, it follows

[TABLE]

For the opposite implication, define $\rho=\frac{1-c}{b}$. Then

[TABLE]

∎

The siTRIP holds for every tensor (possibly with $c$ close to $1$). For $\omega\rightarrow 0$, the constant $c$ converges to the one of the iTRIP and for $\omega\rightarrow\infty$, $c\rightarrow 0$. As well as for Definition 5.10, a slice wise consideration yields a better condition number for just $Z(j)\mathcal{F}^{1/2}$. Note that this value not only bounds the number of steps required for the cg method, but appears to be important for the reconstruction quality obtained through one microstep. However, a further investigation into the siTRIP, how exactly it behaves under perturbation and if it requires modifications, is a matter of future research.

As for the matrix case, we have to limit the singular values from below by a decreasing value proportional to the current residual, cf. Section 6. This leads to a slight complication, which is resolved through the following, simple corollary to Theorem 5.6.

Corollary 5.17**.**

Let $A$ and $\widetilde{A}$ be tensors with standard representations $(\mathcal{L},\Gamma,\mathcal{N},\Theta,\mathcal{R})$ and $(\mathcal{L},\Gamma,\widetilde{\mathcal{N}},\Theta,\mathcal{R})$, respectively. Then both yield the same update $N^{+}$.

Hence, if we want to modify the singular values $\gamma$ and $\theta$, we may do so without knowledge about an appropriate core $\widetilde{\mathcal{N}}$ (for example in the sense of some unknown best approximation). In the least squares problem in Theorem 5.6, we therefor simply set $\gamma_{i}:=\max(\gamma_{i},\sigma_{\mathrm{min}})$, $i=1,\ldots,r_{\mu-1}$ and $\theta_{i}:=\max(\theta_{i},\sigma_{\mathrm{min}})$, $i=1,\ldots,r_{\mu}$ (and can thereby also ignore that the combination of the new $\gamma$ and $\theta$ might not be feasible, cf. [25]

Behavior of the SALSA Filter

A deeper understanding of the regularization utilized by SALSA and the reason for the lower bound $\sigma_{\min}$ is provided by the filter as indicated by \tagform@3.7 for the matrix case and as defined by Corollary 5.8 for tensors. Throughout this section, we assume that the sampling is such that for the minimizer in Theorem 5.6 it (approximately) holds

[TABLE]

which is true at last for $P=\mathcal{I}$ (cf. Corollary 5.8). Since $\zeta_{1}\zeta_{2}=\zeta_{(1,2)}$ (see the later equation \tagform@7.1), we can rewrite

[TABLE]

We are interested in the fixpoints of this update, i.e. we postulate $\mathbb{N}^{+}=\Gamma\ \mathcal{N}\ \Theta$. Then, since $\mathfrak{R}(\mathcal{N}\ \Theta)$ is row-orthogonal (cf. Lemma 5.2), it holds

[TABLE]

where $Z=:\operatorname{diag}(\sigma^{(Z)})$ is necessarily a diagonal matrix (certainly, an analogous argument holds for $\Theta$ as well). Because \tagform@6.2 can only hold if $d_{\sigma^{(Z)},\zeta_{1}}(\gamma_{i})=\gamma_{i}$ for all $i$, the focus of our analysis lies on the fixpoints of the function $d_{\sigma^{(Z)},c}:\sigma\mapsto(1+c\sigma^{-2})^{-1}\sigma^{(Z)}$. For each pair $(\sigma^{(Z)},c)$, the only attractive fixpoint (if existent) is given by $f_{\mathrm{stab}}=\frac{1}{2}\sigma^{(Z)}+\frac{1}{2}\sqrt{(\sigma^{(Z)})^{2}-4c}$ and the repelling one by $f_{\mathrm{rep}}=\frac{1}{2}\sigma^{(Z)}-\frac{1}{2}\sqrt{(\sigma^{(Z)})^{2}-4c}$. At the point where $f_{\mathrm{stab}}=f_{\mathrm{rep}}$, it holds $\sigma=\sqrt{c}=\frac{1}{2}\sigma^{(Z)}$. The minimal value which the term $(1+c\sigma^{-2})^{-1}$ can hence take in any attractive fixpoint, is $F=1/2$. This behavior is shown in Figure 4.

A stabilized singular value corresponds to some attractive fixpoint of $d_{\sigma^{(Z)},c}$ (cf. Definition 3.8), and it necessarily holds $(D_{\omega^{2}\zeta_{1}}({\Gamma}))_{i,i}>0.5\Leftrightarrow\gamma_{i}>\omega\sqrt{\zeta_{1}}$. This explains why the lower limit $\sigma_{\min}$ is necessary. As displayed in Figure 4 (left), for any fixed $\sigma^{(Z)}$, a singular value $\sigma$ must be above a certain threshold (that corresponds to the repelling fixpoint) to be increased by the microstep. It therefore must not converge to zero.

Results Transferred Back to a d-Dimensional Tensor

In this Section, we transfer the previous results for $S=P_{(\mu)}$ and $B=M_{(\mu)}$ to the $d$-dimensional setting. In Remark 5.7, we have $\#_{\mathcal{R}}=size(\mathcal{R})=r_{\theta}\prod_{i=s+1}^{d}n_{i}$, $\#_{\mathcal{N}}=size(\mathcal{N})=r_{\gamma}n_{s}r_{\theta}$, $\#_{\mathcal{L}}=size(\mathcal{L})=r_{\gamma}\prod_{i=1}^{s-1}n_{i}$. By combining modes (cf. Definition 4.7), the sizes of the left as well as right side have been distorted, considering that the degrees of freedom of $\mathcal{L}=\mathcal{G}^{<s}$ and $\mathcal{R}=\mathcal{G}^{>s}$ are given by a sum, not a product, of the degrees of freedom of the single modes (ignoring minor gauge conditions). We choose one of the few remaining options through which the method becomes stable. We artificially set

[TABLE]

using appropriate scalings $s_{1}=s_{1}^{(\mu)}$, $s_{2}=s_{2}^{(\mu)}$ (differently for each mode $\mu$). Otherwise, we will not obtain a stable microstep. Furthermore, the near common parts of the denominators, $\#_{\mathcal{R}}+\#_{\mathcal{N}}+\#_{\mathcal{R}}+2(+2)$, can be incorporated into $\omega^{2}$, so we omit them in the following sense:

Lemma 7.1** (Rescaled target function).**

*The previously discussed rescaling is achieved by choosing *

[TABLE]

Thereby,

[TABLE]

Proof.

First, $\zeta^{(\mu)}_{1}=s_{1}^{2}\frac{\prod_{s=1}^{\mu-1}n_{s}}{E}=\frac{\sum_{s=1}^{\mu-1}n_{s}}{\sum_{s=1}^{d}n_{s}}$, with an analog result for $\zeta^{(\mu)}_{2}$. For the mixed term, we have $\zeta^{(\mu)}_{(1,2)}=\zeta^{(\mu)}_{1}\zeta^{(\mu)}_{2}\frac{E}{E+2}=\zeta^{(\mu)}_{1}\zeta^{(\mu)}_{2}(1-\frac{2}{E+2})$. ∎

The value $E^{-1}$ is in general far below machine accuracy, such that we (from now on) ignore the factor $(1+\mathcal{O}(E^{-1}))$. There might be a more suitable realization of this result and it should be remarked that the exact scalings are not important for the validity of Theorem 5.14. In this context, for fixed $\mu$, the matrices $\mathcal{L}_{S(j)_{1},\cdot\cdot}\in\mathbb{R}^{\noindent|P_{(\mu)}(j)|\times r_{\gamma}}$ and $\mathcal{R}_{:,S(j)_{2}}\in\mathbb{R}^{r_{\theta}\times\noindent|P_{(\mu)}(j)|}$ (cf. \tagform@4.5), are given by

[TABLE]

for $p^{(i_{\ell})}\in P_{(\mu)}(j)$, $\ell=1,\ldots,\noindent|P_{(\mu)}(j)|$ and a representation $G$ for which $\mathcal{L}=G^{<s}$ and $\mathcal{R}=G^{>s}$ (cf. \tagform@4.4).

Remark 7.2** (Case $\mu=1,d$).**

For $\mu=1,d$ in Theorem 5.6, the same formula can be used by formally setting $G^{<1}=\mathcal{L}=1$, $G^{>d}=\mathcal{R}=1$ and $\zeta^{(1)}_{1}=0$, $\zeta^{(d)}_{2}=0$, $\zeta^{(1)}_{(1,2)},\zeta^{(d)}_{(1,2)}=0$, respectively. These comply with the result in the matrix case (cf. Remark 3.6).

Since all microsteps $\mathcal{M}^{\ast}$ are stable, we call this regularized ALS method stable - hence the name SALSA (Stable ALS Approximation). We summarize in Algorithm 3 one full left sweep $\mu=1\rightarrow d$ of SALSA for fixed rank $r$. Note that in practice, the complexity is reduced to the minimal necessary order in the optimal case (cf. Remark 5.9). The simpler matrix case ($d=2$) is carried out in Algorithm 1.

Semi Implicit and Non Uniform Rank Adaption

The rank adaption for tensor completion is carried out analogously to the matrix case, Section 3.2. The exact choices of the following parameters are not important, such that we only indicate them roughly. The specific values which we used in all numerical tests are provided in Section 9.4. The number of minor singular values (cf. Definition 3.8) for each matricization is kept constant, such that

[TABLE]

for each $\mu=1,\ldots,d-1$ for certain constants $f_{\mathrm{minor}}<1$ and $k_{\mathrm{minor}}\in\mathbb{N}$, subject to the theoretical bound $r_{\mu}\leq\min(n_{\mu}r_{\mu-1},n_{\mu+1}r_{\mu+1})$. Furthermore, a common, upper limit $r_{\mu}\leq r_{\mathrm{lim}}$ is applied, which is chosen large enough, but likewise in order to avoid unnecessary computation time. The factor $f_{\mathrm{minor}}$ is related to the filter in Corollary 5.8 and the analysis in Section 6, but has ultimately been chosen empirically. Whenever necessary, then the rank $r_{\mu}$ is decreased through a simple truncation $\sigma^{(\mu)}_{r_{\mu}}\leftarrow 0$, while it is increased using a minor singular value $0<\sigma^{(\mu)}_{r_{\mu}+1}\ll\sigma_{\mathrm{min}}$. In the latter case, the required, corresponding singular vectors can for example be chosen randomly.

As in the matrix case, the lower limit $\sigma_{\min}$ is a fraction $f_{\sigma_{\min}}\ll 1$ of the residual on the sampling set (cf. Algorithm 4). The parameter $\omega>0$ is reduced by a factor $f_{\omega}$ in each iteration, slowly reducing the magnitude of regularization. The factor $1<f_{\omega}\in(f^{(\min)}_{\omega},f^{(\max)}_{\omega})$, in turn, is increased or decreased after each iteration through a simple heuristic, in such a way that

[TABLE]

or rather, that it stays close to this fixed value. The tensor $A^{({\tt iter})}=\tau_{r^{({\tt iter})}}(G^{({\tt iter})})$ is the iterate at iteration number ${\tt iter}$. This adaption ensures that $f_{\omega}$ is not too large as to impair the approximation quality, but neither so small that the required runtime becomes unreasonable.

The algorithm will terminate if one of the following stopping criteria is fulfilled:

• •

stagnation: $\omega\ll\sigma_{\min}$ and $f_{\omega}=f^{(\max)}_{\omega}$

• •

convergence: $\omega\rightarrow 0$ or $\|A^{({\tt iter})}-M\|_{P}/\|M\|_{P}\rightarrow 0$

• •

early stop: $\|A^{({\tt iter})}-M\|_{P_{2}}\gg\min_{i<{\tt iter}}\|A^{(i)}-M\|_{P_{2}}$

We have neglected minor implementation details and practical tweaks in this subsection to focus on the essence of the above criteria, such that we refer to the Matlab code for remaining parts.

The SALSA Algorithm

SALSA (Stable ALS Approximation) for tensors is summarized in Algorithm 4. For the (recommended) choices of tuning parameters also used in the numerical tests, see Subsection 9.4. The Matlab implementation as well as a video showing the rank adaption by means of plotting the singular values during a runtime can be found on the personal webpage of the author Sebastian Krämer, along with all sources that were used to create the presented results.666by the time the paper is written, the address is www.igpm.rwth-aachen.de/team/kraemer The notion of stability we address in this paper does however not mean a low sensibility to roundoff errors accumulated over multiple iterations. We encountered that even different processor architectures or rearrangement of brackets with regard to associativity in products can change the intermediate approximations. As in almost all cases the algorithm does not find the global minimum, this also holds for the final output, both in case of ALS and SALSA. For once, this is not an actual drawback, since even exact arithmetic would not consistently cause the algorithm to find better local minima, but it should be kept in mind when reconstructing results. The order of computational complexity does not exceed $\mathcal{O}(dr^{2}\noindent|P|)$ using coarse cg, where $r=\max_{\mu}r_{\mu}$. Note that the computational complexity per sweep can actually be lower, since not all ranks are kept equal, but some are lower than others.

Numerical Experiments

We consider the following three algorithms:

• •

ALS (modified Algorithm 3 for $\omega\equiv 0$)

• •

SALSA (Algorithm 4, the algorithm proposed in this work)

• •

RTTC (Riemannian cg for tensor train completion [40])

We explain how ranks are adapted for ALS in Section 9.1, shortly present the idea behind RTTC in Section 9.2, give details for data acquisition and measurements in Section 9.3 as well as tuning parameters in Section 9.4. We analyze the results in the latter Section 9.9. For each test, we give a (too large) upper bound $r_{\mathrm{lim}}$ for the maximal rank of the iterates, in order to rule out excessive computation times (although this bound is seldomly reached). We would like to emphasize that, in contrast to rank adaption itself, the one dimensional problem of choosing such a bound is easily controlled for example based on the validation set. For simplicity, we use a common mode size $n=n_{1}=\ldots=n_{d}$.

Rank Adaption for Standard ALS

Since ALS itself is not rank adaptive, the (so far) most promising approach, that is greedy rank adaption, is chosen. When the progress stagnates, the algorithm searches for the highest (new) singular value $\sigma^{(\mu)}_{+}$ which any of the rank increases may yield. These values are estimated as follows. Let $\mu$ be fixed and $G$ be a representation for which $G^{<\mu-1}$ is column-orthogonal and $G^{>\mu}$ is row-orthogonal. Further, let

[TABLE]

We define the core $H(\cdot,\cdot)$, $H(i_{\mu-1},i_{\mu})=\alpha_{i_{\mu-1},i_{\mu}}T(i_{\mu-1},i_{\mu})\in\mathbb{R}^{r_{\mu-2}\times r_{\mu}}$ and stack its entries to form the matrix $\mathfrak{H}\in\mathbb{R}^{r_{\mu-2}n_{\mu-1}\times r_{\mu}n_{\mu}}$. This yields the candidate $\sigma^{(\mu)}_{+}:=\|\mathfrak{H}\|_{2}$, the largest singular value of $\mathfrak{H}$. This approach is very similar to the two-fold microsteps as defined in [21] and the rank adaption in AMEn [9], which are both based on DMRG. It however prevents overfitting, since it is equivalent to performing only one, preconditioned CG step (similarly to a Landweber iteration). Furthermore, our experiments suggest that it is more reliable. The corresponding rank $\mu=\mathop{\rm argmin}_{\widetilde{\mu}}\sigma^{(\widetilde{\mu})}_{+}$ is increased by $1$, using a rank $1$ approximation of $\mathfrak{H}$. Since ALS works differently than SALSA, only some stopping criteria can be overtaken, while additional ones are introduced in order to prevent premature termination but also to avoid unnecessary runtime. No rank decreases are proceeded since this involves tremendous difficulties, of which the most important one is the sheer incapability to decide when and which rank actually to decrease.

The RTTC Algorithm

The article [40], in which RTTC is derived and explained in detail, focuses exclusively on tensor completion using the tensor train format as well (but it can likewise be assumed that it is generalizable to other problems). Instead of alternating optimization, RTTC provides a nonlinear conjugate gradient scheme based on Riemannian optimization, which has comparable computational complexity per sweep. Naturally, the problem of rank adaption also poses a challenge in that setting. Therefore, a heuristic rank adaption is introduced (Algorithm 3 in [40]) which successively tests if a single rank increase yields a tolerable change of the residual on the validation set (based on a parameter $\rho\geq 0$). If so, it continues normally with the next test; otherwise, the algorithm priorly resets to the iterate with previous rank.

We observed however that for example the choice $\rho\neq 0$ worked better given the assignments in Section 9.6, but choices other than $\rho=0$ caused the algorithm to not recover a single instance in case of the rank adaption test tensor in Section 9.8. We therefore used three different choices $\rho\in\{0,0.2,1\}$ ($1$ is the default in [40]) and granted RTTC, as opposed to SALSA, the advantage to choose the best result (based on the test set) for each of the following problem classes.

Some minor modification to RTTC were necessary in order to provide fair tests. The allowed number of iterations per rank increase test was increased, since $10$ turned out to be too few. In exchange, the relative improvement parameter was raised to $10^{-3}$, as lower values did not yield improvements. Since RTTC does not have an actual stopping criterion, we stopped whenever the normalized validation residual increased and at the same time was $1000$ times higher than the normalized sampling residual, or whenever the current validation residual was much higher than any previously obtained one (as described in Section 8). When terminating, each time the iterate with lowest validation residual was chosen (as described in Algorithm 4). These adaptations were done carefully in order to obtain only improvements in the quality of approximation. We did further not compare the required number of iterations of RTTC, and the algorithm was granted as much time as necessary.

Data Acquisition and Measurements

Sampling: In order to obtain a sufficient sampling for each slice of $M$, we generate the set $P$ in a quasi-random way as follows: For each direction $\mu=1,\ldots,d$ and each index $i_{\mu}\in{\mathcal{I}}_{\mu}$ we pick $c_{\mathrm{sf}}\cdot r_{P}^{2}$ indices $i_{1},\ldots,i_{\mu-1},i_{\mu+1},\ldots,i_{d}$ at random (uniformly). This gives in total $|P|\noindent=c_{\mathrm{sf}}\cdot dnr_{P}^{2}$ samples (excluding duplicate samples). The rank $r_{P}$ is artificial, such that $c_{\mathrm{sf}}$ can be interpreted as sampling factor since the number of degrees of freedom of a TT-tensor of common rank $r$ is slightly less than $dnr^{2}$.

Testing: As a test set $C$, we use a set of the same cardinality as $P$ that is generated in the same way. Of course, neither this set nor the values $M|_{C}$ are known by the algorithm. The residuals are then measured with respect to the chosen iterate which had the lowest validation residual (cf. Algorithm 4).

Order of optimization: Instead of the sweep we gave before ($\mu=1,\ldots,d$) for simplicity, we alternate between two sweeps ($\mu=1,\ldots,h,\quad\mu=d,\ldots,h,\quad h=\lfloor d/2\rfloor$) to enhance symmetry.

Averaging: With $\langle\cdot\rangle_{\operatorname{ar}}$ we denote the arithmetic mean and by $\langle\cdot\rangle_{\operatorname{geo}}$ the geometric mean which we use for logarithmic scales.

Implementation Details and Tuning Parameters

All tests for ALS and SALSA were done using a (pure) Matlab implementation. This includes the toolbox multiprod [27], which allows a reasonably swift evaluation of products between arrays of matrices, $H(i)J(i)$, $i=1,\ldots,k$, and is much faster than a plain loop. In contrast, some subfunctions of RTTC are based on .mex routines.

Instead of solving full problems in each microstep, both ALS and SALSA use coarse cg (cf. Remark 5.9), for which the tolerance was empirically chosen low enough such that it did not influence the quality of approximation . Note that the cg steps of RTTC are not comparable, since they are performed on low rank manifolds and used to update all cores at once.

We only list time performances in the appendix, which should be interpreted carefully, while the iteration numbers may provide a clearer picture due to similar computational complexities. All parameters have been chosen equally for all experiments (except $\rho$ for RTTC) with respect to best results, not speed, and could be relaxed for easier problems (or in practice for first trials) to reduce timing considerably. Straightening the tolerances for ALS or RTTC, hence allowing more iterations, did however not lead to notable improvements.

The parameter choices (cf. Section 8) for SALSA are given by $\varepsilon_{\mathrm{progr}}=5\cdot 10^{-3}$, $f_{\mathrm{minor}}=0.5$, $k_{\mathrm{minor}}=2$, $f^{(\min)}_{\omega}=1+5\cdot 10^{-4}$, $f^{(\max)}_{\omega}=1.1$, $f_{\sigma_{\min}}=0.1$. The size of the validation set is $|P_{2}|/(|P|+|P_{2}|)=1/20$. These have in parts been chosen empirically and are recommendable for other problems. We observed that any reasonably close values work as well, the more so for larger sampling sets. The performance is in that sense not based on how close the parameters are to some unknown optimal choices. We also refer to the implementation for all details.

Approximation of a Tensor with Near Uniform Singular Spectrum

At first, we consider the completion of the following tensor:

[TABLE]

This tensor is not low rank, but has well ordered modes and uniformly exponentially decaying singular values. It can therefore very well be approximated with uniform ranks (for a black box, rank adaptive algorithm however, this is not trivial to recognize) and the low variance of results suggests that mostly a near best approximation is found. Hence, standard ALS can barely be outperformed. The results are plotted in Figure 5 (see Appendix C for Table 1).

Approximation of Three Generic Tensors with non Uniform Singular Spectrum

We want to demonstrate how different results can be through proper rank adaption, considering the following three tensors, generated by generic functions:

[TABLE]

In contrast to the tensor in Section 9.5, the modes are not (and hardly can be) ordered in accordance with the TT format. A different ordering may of course yield other results, but we cannot assume to find a better ordering if the approximation fails in the general case. The results are plotted in Figure 6 (see Appendix C for Table 2).

Recovery of Random Tensors with Exact Low Rank

We next consider the recovery of quasi-random tensors with exact low ranks. Although this in practice will never occur, it is a very neutral test777Note that in some papers, uniform distributions on $[0,1]$ are used such that all entries of the target tensor are positive, causing each first singular value to be huge compared to all following ones. This leads to a tremendous simplification of the completion problem. There seems to be no indication yet that the sampling required for the completion of a random tensor is in general close to what is stated for the matrix case [8].. The ranks are generated randomly, but it is ensured that $\langle r\rangle_{\operatorname{ar}}\geq 2/3k$ and $\max(r)\leq k$ for some bound $k\in\mathbb{N}$.

Each of these is generated via a TT representation $A=\tau_{r}(G)$ where we assign to each entry of each block $G_{1},\ldots,G_{d}$ a uniformly distributed random value in $[-0.5,0.5]$. Subsequently, the singular values $\Sigma^{(1)},\ldots,\Sigma^{(d-1)}$ are forced to take uniformly distributed random values in $[0,1]$ (up to scaling). This is achieved by successive replacements of the current values in $G$.

As results, we plot the number of successful recoveries ($\|A-M\|_{C}/\|M\|_{C}<10^{-5}$) for different mode sizes $n$ (each single tuple uniform), dimensions $d$ and maximal ranks $k$ of the target tensor (Figures 7, 8).

Recovery of the Rank Adaption Test Tensor

Last but not least, we consider the recovery of tensors as in Example 4.5, for which $Q_{1},Q_{4},Q_{5}$ and $Q_{6}$ are generated quasi-randomly for each trial. For an explanation of the results in Figure 9, we refer to Section 9.7.

Analysis of Results

SALSA is superior in nearly all observed cases. For tensors which could as well be approximated with uniform ranks, the differences are marginal as is to be expected. The three generic functions show that the residuals can be multiple orders of magnitude better, and although the functions were chosen quite randomly, we do of course not want to over-interpret these specific results. Finally, for the more neutral test of random tensor recovery, the required sampling seems to be overall $4$ to $8$ times lower. For the rank adaption test tensor, the performance of SALSA becomes even better for larger rank $k$ (this is due to the larger total sampling), while greedy ALS runs into the predicted trouble, as does RTTC. In general, RTTC performs slightly better than greedy ALS in the approximation of tensors with exponentially declining singular values, while the latter is slightly better in the random recovery tests. This difference is unlikely because of their optimization techniques, which are known to achieve similar results, but rather due to the different rank adaption heuristics. This insight is also based on closer inspections of single tests, which suggest that the rank adaption of RTTC is inferior even to the one of ALS in the majority of cases.

Conclusions

In this article, we have demonstrated that the most successful completion algorithms are very sensitive to rank changes and that existing rank adaption methods suffer from this.

In order to correct this, as proven for SALSA, we suggested a regularization motivated by averaged microsteps in order to uncouple the optimization of a discrete, technical rank. While there is likely room for improvements and rigorous convergence bounds remain subject to future work, we take the noteworthy numerical results as indication that stability (under truncation) is a worthwhile property. The computational complexity of SALSA is further reduced to the minimal order through use of a coarse cg method. Although we focused on tensor completion (with possibly small sampling sets), the derivations given in this paper allow for a generalization to other semi-elliptic problems. Furthermore, it may be possible to adapt the presented ideas to manifold based methods such as RTTC.

Appendix (Experimental Data)

Following are the precise values for Figures 5 and 6, for $R_{C}:=\|A-M\|_{C}$ and $R_{P}:=\|A-M\|_{P}$.