Riemannian preconditioned algorithms for tensor completion via tensor ring decomposition

TL;DR

This paper introduces Riemannian preconditioned algorithms for tensor completion using tensor ring decomposition, improving efficiency and accuracy through a novel metric and optimization methods with proven convergence.

Contribution

The paper develops a new Riemannian metric on tensor ring decomposition space and proposes gradient-based algorithms with convergence guarantees for tensor completion.

Findings

Algorithms outperform existing methods in efficiency.

Proposed methods achieve better reconstruction accuracy.

Numerical experiments validate effectiveness on real datasets.

Abstract

We propose Riemannian preconditioned algorithms for the tensor completion problem via tensor ring decomposition. A new Riemannian metric is developed on the product space of the mode-2 unfolding matrices of the core tensors in tensor ring decomposition. The construction of this metric aims to approximate the Hessian of the cost function by its diagonal blocks, paving the way for various Riemannian optimization methods. Specifically, we propose the Riemannian gradient descent and Riemannian conjugate gradient algorithms. We prove that both algorithms globally converge to a stationary point. In the implementation, we exploit the tensor structure and adopt an economical procedure to avoid large matrix formulation and computation in gradients, which significantly reduces the computational cost. Numerical experiments on various synthetic and real-world datasets -- movie ratings, hyperspectral images, and high-dimensional functions -- suggest that the proposed algorithms are more efficient and have better reconstruction ability than other candidates.

Tables

Table 1. Table 1: Training and test errors of noisy case

$\sigma$	Error	TR-RGD (RBB)	TR-RGD (exact)	TR-RCG	TR-GD	TR-ALS
1e-3	Training	8.9142e-04	8.9142e-04	8.9142e-04	1.1003e-03	8.9142e-04
	Test	1.1471e-03	1.1471e-03	1.1472e-03	1.4066e-03	1.1470e-03
1e-4	Training	8.9154e-05	8.9154e-05	8.9154e-05	1.1129e-03	8.9154e-05
	Test	1.1458e-04	1.1459e-04	1.1461e-04	1.6735e-03	1.1458e-04
1e-5	Training	8.9155e-06	8.9155e-06	8.9155e-06	5.7432e-04	8.9155e-06
	Test	1.1457e-05	1.1457e-05	1.1461e-05	9.1646e-04	1.1457e-05
1e-6	Training	8.9155e-07	8.9155e-07	8.9155e-07	2.7753e-03	8.9155e-07
	Test	1.1457e-06	1.1457e-06	1.1458e-06	4.0237e-03	1.1458e-06
1e-7	Training	8.9155e-08	8.9155e-08	8.9155e-08	3.8170e-03	8.9156e-08
	Test	1.1456e-07	1.1456e-07	1.1457e-07	5.4727e-03	1.1462e-07
1e-8	Training	8.9156e-09	8.9156e-09	8.9156e-09	9.5924e-04	8.9289e-09
	Test	1.1450e-08	1.1455e-08	1.1451e-08	1.4597e-03	1.1512e-08

Table 2. Table 2: PSNR and relative errors for completion of two hyperspectral images

$p$	Results	TR-RCG	TR-RGD	TR-ALS	HaLRTC	TT-RCG	CP-WOPT	geomCG
		Ribeira House Shrubs
$0.1$	PSNR	38.4310	37.5599	38.3964	21.3404	19.3895	22.6050	19.9492
	relerr	0.1058	0.1170	0.1062	0.7569	0.9476	0.6544	0.8884
$0.3$	PSNR	39.9493	40.0437	38.8803	23.4366	29.6771	24.3917	35.3738
	relerr	0.0888	0.0879	0.1005	0.5946	0.2899	0.5327	0.1505
$0.5$	PSNR	40.4171	40.4871	38.3704	26.1920	29.8119	25.7736	36.2973
	relerr	0.0842	0.0835	0.1066	0.4330	0.2854	0.4544	0.1353
		Bom Jesus Bush
$0.1$	PSNR	40.5529	40.0065	40.1682	24.7499	21.4320	23.5259	21.9541
	relerr	0.1185	0.1262	0.1239	0.7310	1.0711	0.8417	1.0086
$0.3$	PSNR	40.4514	40.3119	39.9008	27.6014	25.3039	23.1848	35.0618
	relerr	0.1199	0.1219	0.1278	0.5265	0.6859	0.8754	0.2230
$0.5$	PSNR	40.4516	39.9737	39.4380	30.1986	24.9847	23.2231	34.9976
	relerr	0.1199	0.1267	0.1348	0.3904	0.7115	0.8715	0.2247

Table 3. Table 3: Test errors for high-dimensional functions

$p$	$\exp \left(-\Vert 𝐱\Vert \right)$				$1/\Vert 𝐱\Vert$
	TR-RGD	TR-RCG	TR-ALS	TT-RCG	TR-RGD	TR-RCG	TR-ALS	TT-RCG
0.001	8.0884e-2	7.4157e-2	7.4161e-2	1.3445e-1	1.7531e-1	1.8106e-1	1.8081e-1	2.6876e-1
0.005	7.3505e-3	8.7366e-3	9.2121e-3	1.5904e-2	3.4428e-2	2.9218e-2	3.2090e-2	1.2899e-1
0.01	6.2650e-3	9.7247e-4	1.8737e-3	4.1233e-3	2.5230e-2	1.7676e-2	1.8697e-2	3.4675e-2
0.05	3.8862e-4	1.5019e-4	1.8218e-4	2.2991e-4	3.8510e-3	3.6002e-3	5.2173e-3	3.5697e-3
0.1	1.2251e-4	5.9871e-5	6.8898e-5	8.2512e-5	7.7886e-4	2.9423e-4	6.0423e-4	7.4727e-4

Table 4. Table 4: Speedup of efficient gradient computation on MovieLens 1M dataset. “Avg. Speedup”: average speedup

#iter	Time (seconds)		Avg. Speedup	Relerr on $\Omega$ (${\epsilon }_{\Omega }$)		Relerr on $\Gamma$ (${\epsilon }_{\Gamma }$)
	Naive	Proposed		Naive	Proposed	Naive	Proposed
1	266.6065	5.9717	-	4.0151	4.0151	4.0233	4.0233
21	4880.3810	37.7889	145.5539$\times$	0.7636	0.7636	0.8674	0.8674
41	9400.1950	69.2155	144.6929$\times$	0.3979	0.3979	0.5268	0.5268
61	13957.4616	100.5884	144.8814$\times$	0.2796	0.2796	0.4084	0.4084
81	18518.4900	132.6191	144.2527$\times$	0.2612	0.2612	0.3880	0.3880

Table 5. Table 5: Rank selections for different tensor formats in numerical experiments. “#params”: number of parameters

Format	MovieLens 1M				Hyperspectral images		High-dimensional functions
	Rank	#params	Rank	#params	Rank	#params	Rank	#params
CP	36	365112	49	496958	110	67430	-	-
Tucker	(36,36,36)	411768	(60,30,18)	516060	(65,65,7)	67506	-	-
TT	(1,9,9,1)	375822	(1,10,10,1)	457100	(1,15,15,1)	78495	(1,5,5,5,1)	1200
TR	(6,6,6)	365112	(6,10,3)	483660	(7,16,7)	66577	(4,4,4,4)	1280

Equations

$$\begin{array}[]{cc}\min&\ \frac{1}{2}\left\|\operatorname{P}_{\Omega}(\llbracket\mathcal{U}_{1},\mathcal{U}_{2},\dots,\mathcal{U}_{d}\rrbracket)-\operatorname{P}_{\Omega}(\mathcal{A})\right\|_{\mathrm{F}}^{2}\\ \mathrm{s.\,t.}&(\mathcal{U}_{1},\mathcal{U}_{2},\dots,\mathcal{U}_{d})\in\mathcal{M}_{\mathcal{U}},\end{array}$$

$${\mathcal{M}}_{\mathcal{U}}\text{mathchar}58={\mathbb{R}}^{r_1\times n_1\times r_2}\times {\mathbb{R}}^{r_2\times n_2\times r_3}\times \cdots \times {\mathbb{R}}^{r_d\times n_d\times r_1}.$$

$$\mathcal{M}\text{mathchar}58={\mathbb{R}}^{n_1\times r_1{r}_2}\times {\mathbb{R}}^{n_2\times r_2{r}_3}\times \cdots \times {\mathbb{R}}^{n_{d-1}\times r_{d-1}{r}_d}\times {\mathbb{R}}^{n_d\times r_d{r}_1}$$

$${\pi }_k\text{mathchar}58(i_1,\dots ,i_{k-1},i_{k+1},\dots ,i_d)\mapsto 1+\sum _{\ell \ne k,\ell =1}^d(i_{\ell }-1)J_{\ell }$$

$$\mathcal{X}=[[{\mathcal{U}}_1,\dots ,{\mathcal{U}}_d]],$$

$$\mathcal{X}(i_1,\dots ,i_d)=\mathrm{tr}\left({\mathbf{U}}_1(i_1)\cdots {\mathbf{U}}_d(i_d)\right),$$

$$\mathrm{tr}({\mathbf{U}}_1(i_1)\cdots {\mathbf{U}}_d(i_d))=\mathrm{tr}\left({\mathbf{U}}_k(i_k)\cdots {\mathbf{U}}_d(i_d){\mathbf{U}}_1(i_1)\cdots {\mathbf{U}}_{k-1}(i_{k-1})\right)$$

$$=⟨\mathrm{vec}({\mathbf{U}}_k(i_k)),\mathrm{vec}(({\mathbf{U}}_{k+1}(i_{k+1})\cdots {\mathbf{U}}_d(i_d){\mathbf{U}}_1(i_1)\cdots {\mathbf{U}}_{k-1}(i_{k-1}){)}^{\mathsf{T}})⟩,$$

$${\mathrm{ten}}_{(2)}({\mathbf{W}}_k)(i_1,i_2,i_3)\text{mathchar}58={\mathbf{W}}_k(i_2,i_1+(i_3-1)r_k)$$

$${\mathbf{U}}_{\ne k}\left({\pi }_k(i_1,\dots ,i_{k-1},i_{k+1},\dots ,i_d)\right)\text{mathchar}58={\left(\prod _{j=k+1}^d{\mathbf{U}}_j(i_j)\prod _{j=1}^{k-1}{\mathbf{U}}_j(i_j)\right)}^{\mathsf{T}}$$

$${\mathbf{X}}_{(k)}={\mathbf{W}}_k{\mathbf{W}}_{\ne k}^{\mathsf{T}},$$

$$\mathcal{M}={\mathbb{R}}^{n_1\times r_1{r}_2}\times {\mathbb{R}}^{n_2\times r_2{r}_3}\times \cdots \times {\mathbb{R}}^{n_{d-1}\times r_{d-1}{r}_d}\times {\mathbb{R}}^{n_d\times r_d{r}_1}.$$

$$\underset{\vec{\mathbf{W}}=\left({\mathbf{W}}_1,\dots ,{\mathbf{W}}_d\right)\in \mathcal{M}}{\min }$$

$$f_{\Omega }(\vec{\mathbf{W}})\text{mathchar}58=\frac{1}{2p}{\Vert P_{\Omega }(\mathcal{X})-P_{\Omega }(\mathcal{A})\Vert }_{\mathrm{F}}^2=\frac{1}{2p}{\Vert P_{{\Omega }_{(k)}}\left({\mathbf{W}}_k{\mathbf{W}}_{\ne k}^{\mathsf{T}}-{\mathbf{A}}_{(k)}\right)\Vert }_{\mathrm{F}}^2,$$

$$r(\vec{\mathbf{W}})\text{mathchar}58=\frac{\lambda }{2}\Vert \vec{\mathbf{W}}{\Vert }_{\mathrm{F}}^2$$

$$\begin{array}[]{c}(\mathcal{U}_{1},\dots,\mathcal{U}_{d})\in\mathcal{M}_{\mathcal{U}}\\[3.0pt] \mathbb{R}^{r_{1}\times n_{1}\times r_{2}}\times\cdots\times\mathbb{R}^{r_{d}\times n_{d}\times r_{1}}\end{array}\ \begin{array}[]{c}\text{matricization}\\[2.0pt] \underscriptrightarrow{\ \ \ \ \ \ \mathbf{W}_{k}=(\mathcal{U}_{k})_{(2)}}\\[5.0pt] \overscriptleftarrow{\ \ \ \ \mathcal{U}_{k}=\operatorname{ten}_{(2)}(\mathbf{W}_{k})}\\[2.0pt] \text{tensorization}\end{array}\ \begin{array}[]{c}\mathcal{M}\ni(\mathbf{W}_{1},\dots,\mathbf{W}_{d})\\[3.0pt] \mathbb{R}^{n_{1}\times r_{1}r_{2}}\times\cdots\times\mathbb{R}^{n_{d}\times r_{d}r_{1}}\end{array}.$$

$$\tau \text{mathchar}58\mathcal{M}\to {\mathbb{R}}^{n_1\times \cdots \times n_d},\tau (\vec{\mathbf{W}})\text{mathchar}58=[[{\mathrm{ten}}_{(2)}({\mathbf{W}}_1),\dots ,{\mathrm{ten}}_{(2)}({\mathbf{W}}_d)]],$$

$${\partial }_{{\mathbf{W}}_k}f(\vec{\mathbf{W}})={\partial }_{{\mathbf{W}}_k}{f}_{\Omega }(\vec{\mathbf{W}})+{\partial }_{{\mathbf{W}}_k}r(\vec{\mathbf{W}})=\frac{1}{p}{\mathbf{S}}_{(k)}{\mathbf{W}}_{\ne k}+\lambda {\mathbf{W}}_k\text{for }k=1,\dots ,d,$$

$$\nabla f_{\Omega }(\vec{\mathbf{W}})=\left(\frac{1}{p}{\mathbf{S}}_{(1)}{\mathbf{W}}_{\ne 1},\frac{1}{p}{\mathbf{S}}_{(2)}{\mathbf{W}}_{\ne 2},\dots ,\frac{1}{p}{\mathbf{S}}_{(d)}{\mathbf{W}}_{\ne d}\right).$$

$$⟨\mathcal{H}(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}],\vec{\mathbf{\eta }}⟩\approx {\nabla }^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }},\vec{\mathbf{\eta }}]$$

$$T_{\vec{\mathbf{W}}}\mathcal{M}$$

$$={\mathbb{R}}^{n_1\times r_1{r}_2}\times \cdots \times {\mathbb{R}}^{n_{d-1}\times r_{d-1}{r}_d}\times {\mathbb{R}}^{n_d\times r_d{r}_1}.$$

$${\nabla }^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }},\vec{\mathbf{\eta }}]=\sum _{k=1}^d⟨{\partial }_{{\mathbf{W}}_k,{\mathbf{W}}_k}^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}],\vec{\mathbf{\eta }}⟩+\sum _{\ell ,m=1,\ell \ne m}^d⟨{\partial }_{{\mathbf{W}}_{\ell },{\mathbf{W}}_m}^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}],\vec{\mathbf{\eta }}⟩$$

$${\partial }_{{\mathbf{W}}_k,{\mathbf{W}}_k}^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}]=\frac{1}{p}{P}_{{\Omega }_{(k)}}\left({\mathbf{\xi }}_k{\mathbf{W}}_{\ne k}^{\mathsf{T}}\right){\mathbf{W}}_{\ne k}\text{for }k=1,\dots ,d.$$

$${\mathcal{H}}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}]\text{mathchar}58=\left({\partial }_{{\mathbf{W}}_1,{\mathbf{W}}_1}^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}],\dots ,{\partial }_{{\mathbf{W}}_d,{\mathbf{W}}_d}^2{f}_{\Omega }(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}]\right),$$

$$\mathcal{H}(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}]$$

$$=\left({\mathbf{\xi }}_1{\mathbf{W}}_{\ne 1}^{\mathsf{T}}{\mathbf{W}}_{\ne 1},\dots ,{\mathbf{\xi }}_d{\mathbf{W}}_{\ne d}^{\mathsf{T}}{\mathbf{W}}_{\ne d}\right).$$

$$g_{\vec{\mathbf{W}}}(\vec{\mathbf{\xi }},\vec{\mathbf{\eta }})\text{mathchar}58=\sum _{k=1}^d\mathrm{tr}\left({\mathbf{\xi }}_k{\mathbf{H}}_k(\vec{\mathbf{W}}){\left({\mathbf{\eta }}_k\right)}^{\mathsf{T}}\right)$$

$${\mathbf{H}}_k(\vec{\mathbf{W}})\text{mathchar}58={\mathbf{W}}_{\ne k}^{\mathsf{T}}{\mathbf{W}}_{\ne k}+\delta {\mathbf{I}}_{r_k{r}_{k+1}}$$

$$g_{\vec{\mathbf{W}}}\left(\mathrm{grad}f(\vec{\mathbf{W}}),\vec{\mathbf{\xi }}\right)=\mathrm{D}f(\vec{\mathbf{W}})[\vec{\mathbf{\xi }}]\text{mathchar}58=⟨\nabla f(\vec{\mathbf{W}}),\vec{\mathbf{\xi }}⟩$$

$$g_{\vec{\mathbf{W}}}(\vec{\mathbf{\eta }},\vec{\mathbf{\xi }})=\sum _{k=1}^d\mathrm{tr}\left({\mathbf{\eta }}_k{\mathbf{H}}_k(\vec{\mathbf{W}}){\left({\mathbf{\xi }}_k\right)}^{\mathsf{T}}\right)=⟨(\mathcal{H}+\delta \mathcal{I})(\vec{\mathbf{W}})[\vec{\mathbf{\eta }}],\vec{\mathbf{\xi }}⟩,$$

$$\mathrm{grad}f(\vec{\mathbf{W}})=(\mathcal{H}+\delta \mathcal{I}{)}^{-1}(\vec{\mathbf{W}})[\nabla f(\vec{\mathbf{W}})].$$

Full text

1 Introduction

The tensor completion problem, as a natural generalization of the matrix completion problem, is a task of recovering a tensor based on its partially observed entries. In practice, datasets collected from real applications are often assumed to have underlying low-rank structure. Therefore, low-rank matrix decompositions are widely used in matrix completion, which can save the computational cost and storage. In the same spirit, low-rank tensor decompositions play a significant role in tensor completion; applications can be found across various fields, e.g., recommendation systems kasai2016low ; dong2022new , image processing liu2012tensor , and interpolation of high-dimensional functions steinlechner2016riemannian ; khoo2021efficient .

In this paper, we consider tensor ring decomposition and focus on the following tensor completion problem with bounded tensor ring rank $\mathbf{r}\mathrel{\mathop{\mathchar 58\relax}}=(r_{1},\dots,r_{d})$. Given a partially observed $d$-th order tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ on an index set $\Omega\subseteq[n_{1}]\times\cdots\times[n_{d}]$, where $[n_{k}]\mathrel{\mathop{\mathchar 58\relax}}=\{1,2,\dots,n_{k}\}$ for $k=1,\dots,d$ and $d\geq 3$ is a positive interger. The tensor completion problem is formulated as follows,

[TABLE]

where $\llbracket\mathcal{U}_{1},\mathcal{U}_{2},\dots,\mathcal{U}_{d}\rrbracket\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ denotes the tensor ring decomposition with core tensors $\mathcal{U}_{k}\in\mathbb{R}^{r_{k}\times n_{k}\times r_{k+1}}$ for $k\in[d]$; see definitions in section 2. $\operatorname{P}_{\Omega}$ denotes the projection operator onto $\Omega$, namely, $\operatorname{P}_{\Omega}(\mathcal{X})(i_{1},\dots,i_{d})=\mathcal{X}(i_{1},\dots,i_{d})$ if $(i_{1},\dots,i_{d})\in\Omega$, otherwise $\operatorname{P}_{\Omega}(\mathcal{X})(i_{1},\dots,i_{d})=0$. The search space of (1.1) is defined by a product space of core tensors, i.e.,

[TABLE]

Related works and motivation

Tensor completion has several different formulations, one type is based on the nuclear norm minimization. In matrix completion, the nuclear norm of a matrix—a convex relaxation of matrix rank—is minimized. Liu et al. liu2012tensor extended the “nuclear norm” to tensor by calculating the sum of nuclear norms of unfolding matrices of a tensor, and applied the alternating direction method of multipliers algorithm to solve the tensor completion problem. Since the tensor data observed in real world may be perturbed by noise, Zhao et al. zhao2022robust focused on the robust tensor completion problem based on the tubal nuclear norm, and proposed a proximal majorization-minimization algorithm. These algorithms require storing tensors in full size, and the number of parameters is $n_{1}\cdots n_{d}$, which scales exponentially to the dimension $d$.

Instead of working with full-size tensors, tensor decompositions kolda2009tensor allow us to take advantage of the low-rank structure in tensor completion problems, and thus can reduce the number of parameters in search space and save storage. In view of the block structure in tensor decompositions, one can develop alternating minimization methods by updating one block while fixing others. Jain and Oh jain2014provable proposed an alternating minimization method for symmetric third order tensors in canonical polyadic (CP) decomposition. For Tucker decomposition, the alternating minimization, also called the alternating least squares (ALS) algorithm, was considered by Andersson and Bro andersson199893 in which the subproblem is a least squares problem. Tensor train (TT) decomposition oseledets2010tt ; oseledets2011tensor , also known as matrix product states (MPS) verstraete2008matrix ; schollwock2011density in computational physics, decomposes a tensor into $d$ core tensors. Grasedyck et al. grasedyck2015variants presented an alternating direction fitting algorithm for the tensor train completion problem. Recently, tensor ring (TR) decomposition was proposed by Zhao et al. zhao2016tensor as a generalization of tensor train decomposition, and is also known as MPS with periodic boundary conditions in computational physics. The ALS algorithm was developed to solve the TR decomposition problem in zhao2016tensor , and it was further applied to tensor completion wang2017efficient . In general, ALS methods in which the number of parameters scales linearly to the dimension $d$ have been proved to be effective for the tensor completion problem. However, they may suffer from overfitting and require careful initialization in practice chen2020tensor .

More recently, Riemannian optimization working with tensor-based manifolds appears to be prosperous for solving completion problems. Since the search space is a manifold, one can benefit from different geometric tools and develop efficient optimization methods on the manifold where the convergence can be guaranteed; see absil2009optimization ; boumal2023intromanifolds for an overview. Acar et al. acar2011scalable proposed an Euclidean nonlinear conjugate gradient method for tensor completion via CP decomposition. Dong et al. dong2022new proposed Riemannian gradient and Riemannian conjugate gradient methods based on a Riemannian metric on the product space of matrices in polyadic decomposition. For tensor completion in Tucker decomposition, Kressner et al. kressner2014low proposed a Riemannian conjugate gradient method on the manifold of rank-constrained tensors. Kasai and Mishra kasai2016low introduced a preconditioned metric and considered the quotient geometry of the manifold via Tucker decomposition. The corresponding Riemannian conjugate gradient algorithm was proposed. Based on geometric properties of the manifold of tensors with fixed TT rank, Steinlechner steinlechner2016riemannian proposed a Riemannian conjugate gradient algorithm. Furthermore, Cai et al. cai2022tensor investigated the quotient geometry on this manifold, and proposed Riemannian gradient, Riemannian conjugate gradient and Riemannian Gauss–Newton algorithms.

Tensor ring decomposition is a generalization of TT decomposition. On the one hand, it provides a flexible choice of tensor rank. Specifically, the unfolding matrices of core tensors in TR are not necessarily of full rank while the mode-1 and mode-3 unfolding matrices of core tensors in TT have to be of full rank. Moreover, TR decomposition permits a more comprehensive exploration of information along mode-$1$ and mode-$d$ by relaxing the rank constraint on the first and last cores in TT from one to an arbitrary interger, i.e., from TT rank $(1,r_{2},\dots,r_{d},1)$ to TR rank $(r_{1},r_{2},\dots,r_{d})$, enabling higher compressibility and flexibility zhao2019learning . On the other hand, although tensors with bounded TR rank do not form a Riemannian submanifold of $\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$, TR decomposition preserves a similar block structure as TT, which allows us to endow the search space with a manifold structure and to develop a Riemannian metric that has a preconditioning effect. More precisely, we develop Riemannian optimization methods on a product space of unfolding matrices

[TABLE]

for the tensor completion problem via TR decomposition by an equivalent reformulation of (1.1). We refer to the search space $\mathcal{M}$ as a Riemannian manifold by endowing $\mathcal{M}$ with a non-Euclidean metric.

Contributions

We formulate the tensor completion problem via tensor ring decomposition, in which the search space $\mathcal{M}$ is a product space of the mode-2 unfolding matrices of the core tensors in tensor ring decomposition. We design a preconditioned metric on the search space, and propose the Riemannian gradient descent and Riemannian conjugate gradient algorithms to solve the tensor completion problem. We prove that every accumulation point of the sequence generated by the proposed algorithms is a stationary point. To the best of our knowledge, this is the first Riemannian formulation of the tensor completion problem in TR decomposition.

The computation of (Riemannian) gradients involves large matrix formulation and multiplication. Its computational cost is exponential to the dimension $d$, which is unaffordable in practice. In order to improve the efficiency of the proposed algorithms, we adopt an economical procedure, which has polynomial complexity to the dimension $d$, to compute the (Riemannian) gradients by exploiting the tensor structure in TR.

We compare the proposed algorithms with existing methods in various tensor completion tasks on both synthetic and real-world datasets, including movie ratings, hyperspectral images, and high-dimensional functions. The numerical results illustrate that the TR-based algorithms have better or comparable recovery performance than the others. In addition, the proposed algorithms are favorably comparable to the alternating least squares algorithm among TR-based algorithms.

Organization

First, the preliminaries of tensor ring decomposition and the tensor completion problem are introduced in section 2. Next, we develop Riemannian algorithms and present computational details in section 3. The convergence results of the proposed algorithms are shown in section 4. Numerical results are reported in section 5. Finally, the conclusion is given in section 6.

2 Preliminaries

In this section, we first introduce the notation in tensor operations and the definition of tensor ring decomposition. Next, a reformulation of the TR-based tensor completion problem (1.1) is described.

The following notions are involved in tensor computations kolda2009tensor . First, we define an index mapping $\pi_{k}$ by

[TABLE]

with $J_{\ell}=\prod_{m=1,m\neq k}^{\ell-1}n_{m}$ for $k=1,\dots,d$. The mode-$k$ unfolding of a tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ is denoted by a matrix $\mathbf{X}_{(k)}\in\mathbb{R}^{n_{k}\times n_{-k}}$, where $n_{-k}\mathrel{\mathop{\mathchar 58\relax}}=\prod_{i\neq k}n_{i}$. The $(i_{1},i_{2},\dots,i_{d})$-th entry of $\mathcal{X}$ corresponds to the $(i_{k},j)$-th entry of $\mathbf{X}_{(k)}$, where $j=\pi_{k}(i_{1},\dots,i_{k-1},i_{k+1},\dots,i_{d})$. Similarly, the mode-$k$ unfolding of indices in an index set $\Omega$ is defined by $\Omega_{(k)}{\mathrel{\mathop{\mathchar 58\relax}}=}\{(i_{k},\pi_{k}(i_{1},\dots,i_{k-1},i_{k+1},\dots,i_{d}))\mathrel{\mathop{\mathchar 58\relax}}\ (i_{1},\dots,i_{d})\in\Omega\}$. The inner product between two tensors $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ is defined by $\langle\mathcal{X},\mathcal{Y}\rangle\mathrel{\mathop{\mathchar 58\relax}}=\sum_{i_{1}=1}^{n_{1}}\cdots\sum_{i_{d}=1}^{n_{d}}\mathcal{X}({i_{1},\dots,i_{d}})\mathcal{Y}({i_{1},\dots,i_{d}})$. The Frobenius norm of a tensor $\mathcal{X}$ is defined by $\|\mathcal{X}\|_{\mathrm{F}}\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{\langle\mathcal{X},\mathcal{X}\rangle}$.

Definition 1 (tensor ring decomposition zhao2016tensor )

Given a $d$-th order tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$, tensor ring decomposition, denoted by

[TABLE]

decomposes $\mathcal{X}$ into $d$ core tensors $\mathcal{U}_{k}\in\mathbb{R}^{r_{k}\times n_{k}\times r_{k+1}}$ for $k=1,\dots,d$ and $r_{d+1}=r_{1}$. Specifically, the $(i_{1},\dots,i_{d})$-th entry of $\mathcal{X}$ is represented by the trace of products of $d$ matrices, i.e.,

[TABLE]

where $\mathbf{U}_{k}(i_{k})=\mathcal{U}_{k}(\mathrel{\mathop{\mathchar 58\relax}},i_{k},\mathrel{\mathop{\mathchar 58\relax}})\in\mathbb{R}^{r_{k}\times r_{k+1}}$ is the lateral slice matrix of $\mathcal{U}_{k}$ for $i_{k}\in[n_{k}]$.

The tuple $\mathbf{r}=(r_{1},\dots,r_{d})$ is referred to as the tensor ring rank (TR rank). Figure 1 illustrates the TR decomposition of a tensor.

In view of the cyclic symmetry of trace operator, we rewrite the $(i_{1},\dots,i_{d})$-th entry of tensor $\mathcal{X}$ in tensor ring decomposition as follows,

[TABLE]

where $\operatorname{vec}(\cdot)$ denotes the column vectorization of a matrix. In fact, $\operatorname{vec}(\mathbf{U}_{k}(i_{k}))^{\mathsf{T}}$ is the $i_{k}$-th row of the mode-$2$ unfolding matrix $\left(\mathcal{U}_{k}\right)_{(2)}\in\mathbb{R}^{n_{k}\times r_{k}r_{k+1}}$ of the core tensor $\mathcal{U}_{k}$. Additionally, given a matrix $\mathbf{W}_{k}\in\mathbb{R}^{n_{k}\times(r_{k}r_{k+1})}$ for fixed $n_{1},\dots,n_{d}$ and $r_{1},\dots,r_{d}$, the second tensorization operator maps $\mathbf{W}_{k}$ to a tensor $\operatorname{ten}_{(2)}(\mathbf{W}_{k})\in\mathbb{R}^{r_{k}\times n_{k}\times r_{k+1}}$ defined by

[TABLE]

for $(i_{1},i_{2},i_{3})\in[r_{k}]\times[n_{k}]\times[r_{k+1}]$. We observe that $(\operatorname{ten}_{(2)}(\mathbf{W}_{k}))_{(2)}=\mathbf{W}_{k}$ holds. Therefore, the second tensorization operator is invertible. The definition of subchain tensors is given as follows.

Definition 2 (subchain tensor)

The subchain tensor $\mathcal{U}_{\neq k}\in\mathbb{R}^{r_{k}\times n_{-k}\times r_{k+1}}$ is defined by its lateral slice matrices, i.e.,

[TABLE]

for $(i_{1},\dots,i_{k-1},i_{k+1},\dots,i_{d})\in[n_{1}]\times\cdots\times[n_{k-1}]\times[n_{k+1}]\times\cdots\times[n_{d}]$ and $k\in[d]$.

Given the mode-2 unfolding matrix $\left(\mathcal{U}_{k}\right)_{(2)}$ and subchain tensor in Definition 2, the mode-$k$ unfolding of a tensor $\mathcal{X}$ in tensor ring decomposition equals to the product of two smaller matrices (see (zhao2016tensor, , Theorem 3.5)), i.e.,

[TABLE]

where $\mathbf{W}_{k}\mathrel{\mathop{\mathchar 58\relax}}=\left(\mathcal{U}_{k}\right)_{(2)}\in\mathbb{R}^{n_{k}\times r_{k}r_{k+1}}$ and $\mathbf{W}_{\neq k}\mathrel{\mathop{\mathchar 58\relax}}=\left(\mathcal{U}_{\neq k}\right)_{(2)}\in\mathbb{R}^{n_{-k}\times r_{k}r_{k+1}}$ for $k=1,\dots,d$. In the light of this fact, we are able to formulate the TR-based tensor problems on mode-2 unfolding matrices $\mathbf{W}_{1},\dots,\mathbf{W}_{d}$. Specifically, the update of a TR tensor $\mathcal{X}=\llbracket\mathcal{U}_{1},\dots,\mathcal{U}_{d}\rrbracket$ involves the updates of each core tensors $\mathcal{U}_{1},\dots,\mathcal{U}_{d}$, which is equivalent to updating the mode-2 unfolding matrices $\mathbf{W}_{1},\dots,\mathbf{W}_{d}$ by matricizations $\mathbf{W}_{k}=(\mathcal{U}_{k})_{(2)}$. In other words, mode-$2$ unfolding matrices are adequate for all tensor-related computations. It is worth noting that given $\mathbf{W}_{1},\dots,\mathbf{W}_{d}$, the matrices $\mathbf{W}_{k}\mathbf{W}_{\neq k}^{\mathsf{T}}$ and $\mathbf{W}_{j}\mathbf{W}_{\neq j}^{\mathsf{T}}$ represent a same tensor after their respective tensorizations for $j\neq k$. However, we never compute $\mathbf{W}_{k}\mathbf{W}_{\neq k}^{\mathsf{T}}$ explicitly in practice.

Subsequently, the tensor completion problem (1.1) can be reformulated on the product space of the mode-2 unfolding matrices of the core tensors in TR decomposition

[TABLE]

The Frobenius norm on $\mathcal{M}$ is defined by $\|{\vec{\mathbf{W}}}\|_{\mathrm{F}}\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{\sum_{k=1}^{d}\|\mathbf{W}_{k}\|_{\mathrm{F}}^{2}}$ for ${\vec{\mathbf{W}}}=\left(\mathbf{W}_{1},\dots,\mathbf{W}_{d}\right)\in\mathcal{M}$.

In this paper, we focus on the following tensor completion problem combined with a sampling set $\Omega$ and a regularization term $r$,

[TABLE]

The objective function $f$ consists of two parts: one is the cost function

[TABLE]

where $p\mathrel{\mathop{\mathchar 58\relax}}=\frac{\lvert\Omega\rvert}{n_{1}\cdots n_{d}}$ is the sampling rate, $\mathbf{A}_{(k)}$ denotes the mode-$k$ unfolding matrix of $\mathcal{A}$, and $\Omega_{(k)}$ is the mode-$k$ unfolding of the sampling set $\Omega$; the other is a regularization term $r({\vec{\mathbf{W}}})$ and we choose

[TABLE]

in a similar fashion as Maximum-margin Matrix Factorization srebro2004maximum with $\lambda>0$. The regularization term keeps the variable ${\vec{\mathbf{W}}}$ in a compact subset of $\mathcal{M}$ and guarantees the convergence; see analysis in section 4.

Since the second tensorization operator $\operatorname{ten}_{(2)}(\cdot)$ is invertible when $n_{1},\dots,n_{d}$ and $r_{1},\dots,r_{d}$ are fixed, it turns out that the search space in (2.4) is equivalent to the one in the original problem (1.1) and the two search spaces are connected by matricizations $\mathbf{W}_{k}=(\mathcal{U}_{k})_{(2)}$ and tensorizations $\mathcal{U}_{k}=\operatorname{ten}_{(2)}(\mathbf{W}_{k})$ for $(\mathbf{W}_{1},\dots,\mathbf{W}_{d})\in\mathcal{M}$ and $(\mathcal{U}_{1},\dots,\mathcal{U}_{d})\in\mathcal{M}_{\mathcal{U}}$, i.e.,

[TABLE]

Additionally, we observe that an element $\vec{\mathbf{W}}=\left(\mathbf{W}_{1},\dots,\mathbf{W}_{d}\right)$ in $\mathcal{M}$ can also represent an element $\llbracket\operatorname{ten}_{(2)}(\mathbf{W}_{1}),\dots,\operatorname{ten}_{(2)}(\mathbf{W}_{d})\rrbracket\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ by TR decomposition. We define the mapping

[TABLE]

which is adopted to generate systhetic data in section 5.

Given the objective function $f$, the first-order derivative with respect to $\mathbf{W}_{k}$ has the following form,

[TABLE]

where $\mathcal{S}\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{P}_{\Omega}(\tau({\vec{\mathbf{W}}}))-\operatorname{P}_{\Omega}(\mathcal{A})$ is called the residual tensor and $\mathbf{S}_{(k)}$ is the mode-$k$ unfolding matrix of tensor $\mathcal{S}$. Therefore, the Euclidean gradient of $f_{\Omega}$ at ${\vec{\mathbf{W}}}\in\mathcal{M}$ can be computed as follows,

[TABLE]

Moreover, one can develop Euclidean gradient descent algorithms (e.g., yuan2018higher ) to solve the completion problem (2.4). Recently, Riemannian preconditioned algorithms are considered in which the search space is endowed with a non-Euclidean metric. The construction of this metric aims to approximate the Hessian of the cost function by its “diagonal blocks”. These algorithms improve the performance of Euclidean-based algorithms, and are successfully applied to matrix and tensor completion (e.g., mishra2012riemannian ; kasai2016low ; dong2022new ; cai2022tensor ). However, the extension to TR is not straightforward since it involves large matrix formulation and computation. In next section, we consider how to develop efficient preconditioned algorithms for the tensor completion problem (2.4).

3 Tensor completion algorithms

We first develop a preconditioned metric on the manifold $\mathcal{M}$. The corresponding Riemannian gradient is derived under this metric. Next, we propose the Riemannian gradient descent and Riemannian conjugate gradient algorithms. An efficient procedure for computing the Riemannian gradient is proposed in the end.

3.1 A preconditioned metric

The idea of developing a preconditioned metric on $\mathcal{M}$ is to take advantage of the second-order information of the cost function $f_{\Omega}$ and to formulate a search direction that approximates the Newton direction. Specifically, we intend to construct an operator $\mathcal{H}({\vec{\mathbf{W}}})\mathrel{\mathop{\mathchar 58\relax}}\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}\to\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$ such that

[TABLE]

for all $\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}\simeq\mathcal{M}$, where $\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$ denotes the tangent space to $\mathcal{M}$ at ${\vec{\mathbf{W}}}\in\mathcal{M}$ and $\nabla^{2}f_{\Omega}$ denotes the (Euclidean) Hessian of $f_{\Omega}$. Note that

[TABLE]

Therefore, a tangent vector $\vec{\boldsymbol{\xi}}\in\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$ can be expressed by $\vec{\boldsymbol{\xi}}=(\boldsymbol{\xi}_{1},\boldsymbol{\xi}_{2},\dots,\boldsymbol{\xi}_{d})$ with $\boldsymbol{\xi}_{k}\in\mathbb{R}^{n_{k}\times r_{k}r_{k+1}}$.

To this end, we start from computing the explicit form of the Hessian operator

[TABLE]

for all $\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$. By direct calculations, the “diagonal blocks” of $\nabla^{2}f_{\Omega}({\vec{\mathbf{W}}})$ have the following forms,

[TABLE]

Since computing the “off-diagonal blocks” $\partial_{\mathbf{W}_{\ell},\mathbf{W}_{m}}^{2}f_{\Omega}({\vec{\mathbf{W}}})$ is complicated, we consider only the “diagonal blocks” as a trade-off between accuracy and computational cost to form an operator $\mathcal{H}({\vec{\mathbf{W}}})$. An instinctive approach to construct the operator is to directly apply the second-order derivatives, i.e.,

[TABLE]

which relies on a specific sampling set $\Omega$. We intend to design an operator that is applicable to a class of tensor completion problem with a subsampling pattern by taking expectation on (3.8). More accurately, we suppose that the indices in $\Omega$ are i.i.d. samples from the Bernoulli distribution with probability $p$. We eliminate the projection operator $\operatorname{P}_{\Omega}$ by taking expectation on “diagonal blocks” over $\Omega$. Hence, we can define the operator $\mathcal{H}({\vec{\mathbf{W}}})$ by

[TABLE]

Note that both $\mathcal{H}_{\Omega}$ and $\mathcal{H}$ approximate the second-order information of the cost function. The operator $\mathcal{H}$ is an expectation form of $\mathcal{H}_{\Omega}$. We adopt the operator (3.9) for computational convenience but one can consider other subsampling patterns in practice. Consequently, we define a new metric as follows.

Definition 3 (preconditioned metric)

$g$ is an inner product on $\mathcal{M}$ defined by

[TABLE]

for all ${\vec{\mathbf{W}}}\in\mathcal{M}$ and $\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\eta}}\in\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$, where $\mathbf{H}_{k}({\vec{\mathbf{W}}})\in\mathbb{R}^{r_{k}r_{k+1}\times r_{k}r_{k+1}}$ is a matrix defined by

[TABLE]

with a constant parameter $\delta>0$ and the identity matrix $\mathbf{I}_{r_{k}r_{k+1}}\in\mathbb{R}^{r_{k}r_{k+1}\times r_{k}r_{k+1}}$.

Since the matrix $\mathbf{W}_{\neq k}^{\mathsf{T}}\mathbf{W}_{\neq k}$ is not necessarily positive definite, a shifting term $\delta\mathbf{I}_{r_{k}r_{k+1}}$ is added to avoid singularity. Moreover, $g$ is smooth on $\mathcal{M}$, and thus is a well-defined Riemannian metric on $\mathcal{M}$. It turns out that $\mathcal{M}$ is a Riemannian manifold endowed with $g$ and the norm of a tangent vector $\vec{\boldsymbol{\xi}}\in\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$ can be defined by $\|\vec{\boldsymbol{\xi}}\|_{\vec{\mathbf{W}}}\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{g_{\vec{\mathbf{W}}}(\vec{\boldsymbol{\xi}},\vec{\boldsymbol{\xi}})}$. The Riemannian gradient (absil2009optimization, , Sect. 3.6), $\operatorname{grad}\!f({\vec{\mathbf{W}}})$, is the unique element in $\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$ that satisfies

[TABLE]

for all $\vec{\boldsymbol{\xi}}\in\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$. Note that

[TABLE]

where $\mathcal{I}({\vec{\mathbf{W}}})\mathrel{\mathop{\mathchar 58\relax}}\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}\to\operatorname{T}_{\vec{\mathbf{W}}}\mathcal{M}$ is the identity operator. It turns out that

[TABLE]

In view of (2.4) and (3.7), $\operatorname{grad}\!f$ is an approximation of the Newton direction of $f$. As a result, the metric (3.10) has a preconditioning effect on the Euclidean gradient. Therefore, we refer to the new metric (3.10) as a preconditioned metric on $\mathcal{M}$. In summary, the Riemannian gradient can be computed from (2.6) as follows.

Proposition 1 (Riemannian gradient)

The Riemannian gradient of $f$ at ${\vec{\mathbf{W}}}\in\mathcal{M}$ with respect to the metric $g$ is

[TABLE]

3.2 Riemannian preconditioned algorithms

Using the preconditioned metric (3.10) and the Riemannian gradient (3.11), we propose the Riemannian gradient descent and Riemannian conjugate gradient algorithms to solve the tensor completion problem (2.4).

Riemannian gradient descent

The Riemannian gradient descent algorithm is listed in Algorithm 1. Note that the retraction on the manifold $\mathcal{M}$ is the identity map. For the selection of stepsize, we consider the following two strategies: 1) we adopt exact line search by solving the following optimization problem

[TABLE]

Since $h$ is a polynomial of $s$ with degree $2d$, the solution $s^{(t)}_{\mathrm{exact}}$ is the root of the polynomial $h^{\prime}(s)$ with $2d-1$ degree; 2) alternatively, we consider Armijo backtracking line search. Given the initial stepsize $s_{0}^{(t)}>0$, find the smallest integer $\ell$, such that for $s^{(t)}=\rho^{\ell}s_{0}^{(t)}>{s}_{\min}$, the inequality

[TABLE]

holds, where $\rho,a\in(0,1),{{s}_{\min}}>0$ are backtracking parameters. The Riemannian Barzilai–Borwein (RBB) stepsize iannazzo2018riemannian , defined by

[TABLE]

where ${\vec{\mathbf{Z}}}^{(t-1)}\mathrel{\mathop{\mathchar 58\relax}}={\vec{\mathbf{W}}}^{(t)}-{\vec{\mathbf{W}}}^{(t-1)}$ and ${\vec{\mathbf{Y}}}^{(t-1)}\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{grad}\!f({\vec{\mathbf{W}}}^{(t)})-\operatorname{grad}\!f({\vec{\mathbf{W}}}^{(t-1)})$, appears to be favorable in many applications. Therefore, we set RBB to be the initial stepsize $s_{0}^{(t)}$. It is worth noting that the initial stepsize $s_{0}^{(0)}$ can also be generated by exact line search (3.12) in practice.

We illustrate the connection and difference between the alternating least squares algorithm for tensor ring completion wang2017efficient (TR-ALS) and the proposed TR-RGD algorithm in the following remark.

Remark 1

Both TR-ALS and TR-RGD can be viewed as line search methods. In contrast with the TR-ALS, where $\mathbf{W}_{1},\dots,\mathbf{W}_{d}$ are updated sequentially with exact line search, $(\mathbf{W}_{1},\dots,\mathbf{W}_{d})$ are treated as one point ${\vec{\mathbf{W}}}\in\mathcal{M}$ in the proposed TR-RGD and updated by the Riemannian gradient descent algorithm. The proposed TR-RGD benefits from a preconditioned metric (3.10) which is exquisitely tailored for the tensor completion problem (2.4). Moreover, TR-RGD allows a more flexible choice of stepsize rules, e.g., exact line search and Riemannian BB stepsize. Therefore, the TR-RGD method can be potentially competitive in practice.

Riemannian conjugate gradient

The Riemannian conjugate gradient algorithm is given in Algorithm 2. For CG parameter $\beta^{(t)}$, we consider the Riemannian version boumal2014manopt of the modified Hestenes–Stiefel rule (HS+) hestenes1952methods

[TABLE]

The stepsize in Algorithm 2 is determined by Armijo backtracking line search (3.13).

In both algorithms, the computation of Riemannian gradient, involving large matrix formulation and multiplication, dominates the total cost. Since this cost is exponential to $d$, a straightforward implementation is not affordable in practice. To this end, we have to come up with a procedure that can compute the Riemannian gradient efficiently.

3.3 Efficient computation for gradients

We investigate the computational details of the Riemannian gradient (3.11) in this subsection. Generally, the computation of $\operatorname{grad}\!f({\vec{\mathbf{W}}})=(\boldsymbol{\eta}_{1},\dots,\boldsymbol{\eta}_{d})$ involves two steps: the first is computing the Euclidean gradient $\nabla f({\vec{\mathbf{W}}})=(\mathbf{G}_{1},\dots,\mathbf{G}_{d})$ in (2.6); the second is assembling the Riemannian gradient $\operatorname{grad}\!f({\vec{\mathbf{W}}})$, where

[TABLE]

for $k=1,\dots,d$.

There are five operations in straightforward calculating the gradients:

1. form the matrix $\mathbf{W}_{\neq k}$ for $k=1,\dots,d$, which requires $2\sum_{k=1}^{d}n_{-k}R_{k}$ flops, where $R_{k}\mathrel{\mathop{\mathchar 58\relax}}=r_{k}\left(\sum_{j=1,j\neq k}^{d}r_{j}r_{j+1}\right)$; 2) compute the sparse tensor $\mathcal{S}$, requiring $2\lvert\Omega\rvert r_{1}r_{2}$ flops; 3) compute the Euclidean gradient by sparse-dense matrix product which involves $2\lvert\Omega\rvert\bar{r}$ flops, where $\bar{r}\mathrel{\mathop{\mathchar 58\relax}}=\sum_{k=1}^{d}r_{k}r_{k+1}$; 4) compute $\mathbf{H}_{k}({\vec{\mathbf{W}}})$ by a dense-dense matrix multiplication with $2\sum_{k=1}^{d}n_{-k}(r_{k}r_{k+1})^{2}$ flops; 5) compute the Riemannian gradient through Cholesky decomposition for linear systems, which requires $\sum_{k=1}^{d}\left(2n_{k}(r_{k}r_{k+1})^{2}+C_{\mathrm{chol}}(r_{k}r_{k+1})^{3}\right)$ flops. Following above operations, the straightforward computation of the Riemannian gradient totally requires

[TABLE]

flops. If $n_{1}=\cdots=n_{d}=n$ and $r_{1}=\cdots=r_{d}=r$, it boils down to

[TABLE]

flops in total. In practice, the terms with order of $\mathcal{O}(n^{d-1})$ dominate the computational cost.

By using the Kronecker product structure of $\mathbf{W}_{\neq k}^{\mathsf{T}}\mathbf{W}_{\neq k}$, the total cost can be significantly reduced. Algorithm 3 illustrates how to compute $\mathbf{W}_{\neq k}^{\mathsf{T}}\mathbf{W}_{\neq k}$ without forming $\mathbf{W}_{\neq k}$ explicitly. The matrix multiplication is

[TABLE]

where $\tilde{\mathbf{w}}_{k}(\mathbf{i}_{-k})\mathrel{\mathop{\mathchar 58\relax}}=\operatorname{vec}((\prod_{j=k+1}^{d}\mathbf{U}_{j}(i_{j})\prod_{j=1}^{k-1}\mathbf{U}_{j}(i_{j}))^{\mathsf{T}})$ and $\mathbf{i}_{-k}\mathrel{\mathop{\mathchar 58\relax}}=(i_{k+1},\dots,i_{d},$ $i_{1},\dots,i_{k-1})\in[n_{k+1}]\times\cdots\times[n_{d}]\times[n_{1}]\times\cdots\times[n_{k-1}]$. By using $\operatorname{vec}(\mathbf{C}\mathbf{X}\mathbf{B}^{\mathsf{T}})=(\mathbf{B}\otimes\mathbf{C})\operatorname{vec}(\mathbf{X})$ for matrices $\mathbf{B},\mathbf{C},\mathbf{X}$ in appropriate size, we have

[TABLE]

where $\otimes$ denotes the Kronecker product. Taking (3.17) into (3.16), we can compute $\mathbf{W}_{\neq k}^{\mathsf{T}}\mathbf{W}_{\neq k}$ in a recursive way,

[TABLE]

It follows that $\mathbf{W}_{\neq k}^{\mathsf{T}}\mathbf{W}_{\neq k}=\tilde{\mathbf{H}}_{d-1}$. The computation of $\mathbf{W}_{\neq k}^{\mathsf{T}}\mathbf{W}_{\neq k}$ for $k\in[d]$ requires

[TABLE]

flops in Algorithm 3. If $n_{1}=\cdots=n_{d}$ and $r_{1}=\cdots=r_{d}$, computing the Riemannian gradient requires

[TABLE]

flops in total, which has the terms with order of $\mathcal{O}(n)$. In order to verify the improvement of Algorithm 3, we report a numerical comparison on a real dataset in Appendix A.

4 Convergence Analysis

The global convergence of TR-RGD and TR-RCG is analyzed in this section. Let $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ be an infinite sequence generated by Algorithm 1 or Algorithm 2. We prove that every accumulation point of $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ is a stationary point; see Theorem 4.1 and Theorem 4.2.

Lemma 1 ((boumal2019global, , Lemma 2.7))

Let $\mathcal{M}^{\prime}\subseteq\mathcal{M}$ be a compact Riemannian submanifold. Let $\mathcal{R}_{x}\mathrel{\mathop{\mathchar 58\relax}}\operatorname{T}_{x}\mathcal{M}^{\prime}\to\mathcal{M}^{\prime}$ be retraction. $f\mathrel{\mathop{\mathchar 58\relax}}\mathcal{M}^{\prime}\to\mathbb{R}$ has Lipschitz continuous gradient on the convex hull of $\mathcal{M}^{\prime}$. Then, there exists a constant $L>0$, such that

[TABLE]

Since the search space $\mathcal{M}$ is flat, the retraction map in Lemma 1 is chosen as the identity map. We observe the coercivity of $f$ in (2.4) from the regularization term $\frac{\lambda}{2}\|{\vec{\mathbf{W}}}\|_{\mathrm{F}}^{2}$. Hence, the level set $\mathcal{L}\mathrel{\mathop{\mathchar 58\relax}}=\{{\vec{\mathbf{W}}}\mathrel{\mathop{\mathchar 58\relax}}\ f({\vec{\mathbf{W}}})\leq f({\vec{\mathbf{W}}}^{(0)})\}$ is compact. The sequence of function values $\{f({\vec{\mathbf{W}}}^{(t)})\}_{t\geq 0}$, obtained from (3.12) and (3.13), is monotonically decreasing. It holds that $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ is a bounded sequence in $\mathcal{L}$ with

[TABLE]

Therefore, there exist accumulation points for the sequence $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$. Moreover, the objective function $f$ has Lipschitz continuous gradient. By using Lemma 1 and (dong2022new, , Proposition 4.3), we can prove the global convergence of the proposed Riemannian gradient descent algorithm in a same fashion.

Theorem 4.1

Let $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ be an infinite sequence generated by Algorithm 1. Then, there exists $C>0$ such that $f({\vec{\mathbf{W}}}^{(t)})-f({\vec{\mathbf{W}}}^{(t+1)})>C\|\operatorname{grad}\!f({\vec{\mathbf{W}}}^{(t)})\|_{{\vec{\mathbf{W}}}^{(t)}}$. Furthermore, we have: 1) every accumulation point of $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ is a stationary point of $f$; 2) the algorithm returns ${\vec{\mathbf{W}}}\in\mathcal{M}$ satisfying $\|\operatorname{grad}\!f({\vec{\mathbf{W}}})\|_{{\vec{\mathbf{W}}}}<\epsilon$ after $\left\lceil{f({\vec{\mathbf{W}}}^{(0)})}/{(C\epsilon^{2})}\right\rceil$ iterations at most.

Now, we discuss the global convergence of the proposed Riemannian conjugate gradient algorithm in Algorithm 2; interested readers are referred to absil2009optimization ; sato2022riemannian for the convergence of RCG on general manifolds. Here, we follow the convergence analysis in absil2009optimization . To fulfill the basic assumptions in  (absil2009optimization, , Theorem 4.3.1), i.e., the sequence of search directions $\{\vec{\boldsymbol{\eta}}^{(t)}\}_{t\geq 0}$ is gradient-related to $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$, we enforce $\vec{\boldsymbol{\eta}}^{(t)}$ to be

[TABLE]

Therefore, $\{\vec{\boldsymbol{\eta}}^{(t)}\}_{t\geq 0}$ are descent directions with

[TABLE]

In addition, it follows from Algorithm 2 that

[TABLE]

is uniformly bounded. In view of (4.19) and (4.20), $\{\vec{\boldsymbol{\eta}}^{(t)}\}_{t\geq 0}$ is gradient-related to $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$. By using (absil2009optimization, , Theorem 4.3.1), we have the following result.

Theorem 4.2

Let $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ be an infinite sequence generated by Algorithm 2 with modified search direction (4.18). Then, every accumulation point of $\{{\vec{\mathbf{W}}}^{(t)}\}_{t\geq 0}$ is a stationary point of $f$.

5 Numerical Experiments

In this section, we numerically compare TR-RGD (Algorithm 1) and TR-RCG (Algorithm 2) with existing algorithms based on different tensor decompositions on synthetic and real-world datasets, including movie ratings, hyperspectral images, and high-dimensional functions. First, we introduce all the compared algorithms and default settings.

We consider the Riemannian conjugate gradient algorithm111TTeMPS toolbox: https://www.epfl.ch/labs/anchp/index-html/software/ttemps/. in steinlechner2016riemannian based on tensor train decomposition, which is denoted by “TT-RCG”. For CP-based algorithm, we choose “CP-WOPT” acar2011scalable in Tensor-Toolbox222Tensor-Toolbox v3.4: http://www.tensortoolbox.org/. by Bader and Kolda bader2008efficient with limited-memory BFGS algorithm333Available from https://github.com/stephenbeckr/L-BFGS-B-C. recommended by the authors. “GeomCG” kressner2014low is a Riemannian conjugate gradient algorithm444GeomCG toolbox: https://www.epfl.ch/labs/anchp/index-html/software/geomcg/. for Tucker-based tensor completion problem. In addition, we consider a nuclear-norm-based algorithm555Available from https://github.com/andrewssobral/mctc4bmi/tree/master/algs_tc/LRTC. “HaLRTC” liu2012tensor ; sobral_sbmi_prl_2016 . If not specified, we adopt default settings for all the compared algorithm.

Note that there are different ways to choose stepsize. In the preliminary numerical experiments, TR-RGD with stepsize RBB2 in (3.14) performs better than the algorithm with RBB1. Therefore, we only consider stepsize RBB2 in the following comparisons. “TR-RGD (exact)” and “TR-RGD (RBB)” denote Algorithm 1 with exact line search (3.12) and Algorithm 1 with Armijo backtracking (3.13) and RBB2, respectively. Moreover, the vanilla Euclidean gradient descent algorithm, denoted by “TR-GD”, is implemented for comparison. The stepsize for TR-GD is based on Armijo backtracking and the standard BB bb1988 . The default settings of line search parameters are $\rho=0.4$, $a=10^{-5}$, and $s_{\min}=10^{-10}$. In addition, we implement the alternating least squares algorithm wang2017efficient ; zhao2019learning called “TR-ALS” for the tensor completion problem (2.4).

All algorithms are initialized from $\mathcal{X}^{(0)}=\tau({\vec{\mathbf{W}}})$ in which $\tau$ is defined in (2.5) and ${\vec{\mathbf{W}}}$ is randomly generated. We define the relative error on a sampling set $\Omega$

[TABLE]

We refer to the training error as $\varepsilon_{\Omega}({\vec{\mathbf{W}}})$. Furthermore, we evaluate the test error $\varepsilon_{\Gamma}({\vec{\mathbf{W}}})$ on a test set $\Gamma$ different from $\Omega$. The default setting of $\lvert\Gamma\rvert$ is $100$. For stopping criteria, we terminate the algorithms once one of the following criteria is reached: 1) the relative error $\varepsilon_{\Omega}({\vec{\mathbf{W}}}^{(t)})<10^{-12}$; 2) the relative change $\lvert{(\varepsilon_{\Omega}({\vec{\mathbf{W}}}^{(t)})-\varepsilon_{\Omega}({\vec{\mathbf{W}}}^{(t-1)}))}/{\varepsilon_{\Omega}({\vec{\mathbf{W}}}^{(t-1)})}\rvert<\varepsilon$; 3) the gradient norm $\|\operatorname{grad}\!f({\vec{\mathbf{W}}})\|_{\mathrm{F}}<\varepsilon$. This criterion is only activated for Riemannian methods; 4) maximum iteration number; 5) time budget. The tolerance $\varepsilon$ is chosen as $10^{-8}$. All experiments are performed on a MacBook Pro 2019 with MacOS Ventura 13.1, 2.4 GHz 8 core Intel Core i9 processor, 32GB memory, and Matlab R2020b. The codes of TR-RGD and TR-RCG can be downloaded from https://github.com/JimmyPeng1998.

5.1 Synthetic data

In this subsection, we investigate the numerical performance and the reconstruction ability of tensor completion algorithms in TR decomposition on noiseless and noisy observations.

Noiseless observations

We consider a synthetic low-rank tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}}$ in (2.4) generated by

[TABLE]

where $\tau$ is defined in (2.5), each entry of ${\vec{\mathbf{W}}}\in\mathcal{M}$ is uniformly sampled from $[0,1]$, $d=3$, $n_{1}=n_{2}=n_{3}=100$, and TR rank $\mathbf{r}^{*}=(6,6,6)$. Given the sampling rate $p$, we formulate the sampling set $\Omega$ by randomly selecting $pn_{1}\cdots n_{d}$ samples from $[n_{1}]\times\cdots\times[n_{d}]$. In order to get an unbiased recovery result, we choose the regularization parameter $\lambda=0$, the sampling rate $p=0.3$, and $\mathbf{r}=(6,6,6)$. The time budget is $1800$s and the maximum iteration number is $10000$.

Figure 2 presents the numerical results of noiseless case. First, we observe that all algorithms successfully recover the true tensor within time budget. TR-RGD (RBB) and TR-RCG perform better than the alternating least squares algorithm (TR-ALS) in terms of training error and test error, and TR-RGD (RBB) is comparable to TR-RCG. Second, TR-RGD with RBB stepsize is more efficient than the algorithm with exact line search, since computing exact line seach needs to find the roots for a polynomial of degree $2d-1$. However, computing the coefficients of such polynomial is expensive in practice. Third, it is worth noting that the proposed Riemannian gradient algorithms outperform the Euclidean gradient algorithm with BB stepsizes, implying that the new metric have a preconditioning effect indeed.

The question “how many samples are required to recover a low-rank data” is interesting but challenging. In the matrix case $(d=2)$, around $\mathcal{O}(nr\log n)$ samples are necessary to recover a low-rank matrix keshavan2009matrix ; candes2010power . We aim to numerically investigate this question by exploring the relationship between sampling rate and tensor size $n$ for third-order tensors with TR rank $\mathbf{r}^{*}=(3,3,3)$. To this end, we randomly generate third order synthetic tensors $\mathcal{A}$ with TR rank $\mathbf{r}=(3,3,3)$ in the same fashion as above, and select tensor size $n\mathrel{\mathop{\mathchar 58\relax}}=n_{1}=n_{2}=n_{3}$ from $\{60,70,\dots,180\}$ and the sample size $\lvert\Omega\rvert$ from $\{1000,1500,\dots,20000\}$. We run TR-RGD, TR-RCG, and TR-ALS algorithms five times for each combination of $n$ and $\lvert\Omega\rvert$. One algorithm is regarded to successfully recover a tensor if the test error $\varepsilon_{\Gamma}<10^{-4}$ within maximum iteration number 250. Figure 3 illustrates the phase plots of recovery results, in which the gray level of each block represents the number of successful recovery. White blocks imply successful recovery in all five runs. The red line represents $\mathcal{O}(n\log n)$ when $n$ scales. The phase plots suggest a similar behavior to the matrix case, which is consistent to other numerical experiments for $d=3$, e.g., kressner2014low .

Noisy observations

Furthermore, we investigate the reconstruction ability of TR-based algorithms under different noise levels. We consider a synthetic noisy tensor

[TABLE]

where $\hat{\mathcal{A}}$ is generated in the same rule as the noiseless case with TR rank $\mathbf{r}^{*}=(3,3,3)$ and $n_{1}=n_{2}=n_{3}=100$. $\mathcal{E}$ is a tensor with its entries being sampled from normal distribution $\mathcal{N}(0,1)$. The parameter $\sigma$ measures the noise level of a tensor. Ideally, the relative errors $\varepsilon_{\Omega}$ and $\varepsilon_{\Gamma}$ are supposed to be the noise level $\sigma$ when the algorithm terminates. We set $\sigma=10^{-3},10^{-4},10^{-5},10^{-6},10^{-7},10^{-8}$, and the regularization parameter $\lambda=10^{-12}$ to get unbiased recovery results. We observe from Fig. 3 that it requires at least $8000$ samples to recover the noiseless tensor, i.e., $p$ should be larger than $8000/100^{3}=0.008$. Therefore, we choose the sampling rate $p=0.05$. The time budget is $120$s and the maximum iteration number is $1000$.

Figure 4 shows the recovery performance of TR-based algorithms under different noise levels. All algorithms successfully recover the underlying low-rank tensor within time budget except TR-GD. Since the TR rank parameter $\mathbf{r}$ is exactly the true TR rank $\mathbf{r}^{*}$ of data tensor $\mathcal{A}$, we observe from Table 1 that the test errors are comparable to and slightly larger than the training errors.

5.2 Experiments on MovieLens 1M dataset

We consider a real-world tensor completion problem on the MovieLens 1M666Available from https://grouplens.org/datasets/movielens/1m/. dataset including one million movie ratings from 6040 users on 3952 movies from September 19th, 1997 to April 22nd, 1998. We formulate the movie ratings as a third order tensor $\mathcal{A}$ of size $6040\times 3952\times 150$ by choosing one week as a period. We randomly select $80\%$ of the known ratings as training set $\Omega$ and the rest $20\%$ ratings are test set $\Gamma$. We compared the proposed algorithms TR-RGD and TR-RCG with other tensor completion algorithms including TR-GD, TR-ALS, TT-RCG, CP-WOPT, and geomCG. We choose the parameters: TR rank $\mathbf{r}=(6,10,3)$ and $\mathbf{r}=(6,6,6)$, regularization parameter $\lambda=1$, and the time budget 500s. To ensure a close number of parameters in different search spaces, we choose the TT rank as $(1,10,10,1)$ and $(1,9,9,1)$, the Tucker rank as $(60,30,18)$ and $(36,36,36)$, and the CP rank as $49$ and $36$; see Appendix B for details.

The test error under two different rank selections are shown in Fig. 5. We observe that the proposed TR-RGD (RBB) and TR-RCG are comparable with CP-WOPT since CP-WOPT also benefits from the second-order information. They have faster convergence with lower test errors than the others. Moreover, TR-RGD (exact) spends more time than TR-RGD (RBB). For the sake of brevity, we only consider TR-RGD (RBB) in the following experiments and simplify “TR-RGD (RBB)” to “TR-RGD”.

5.3 Experiments on hyperspectral images

In this experiment, we consider the completion task on hyperspectral images. A hyperspectral image is formulated as a third order tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$. Mode three of $\mathcal{A}$ represents $n_{3}$ wavelength values of light. Mode one and two represents reflectance level of light under different wavelengths. We select two images777Image source: hsi_34.mat and hsi_46.mat from https://figshare.manchester.ac.uk/articles/dataset/Fifty_hyperspectral_reflectance_images_of_outdoor_scenes/14877285. from “50 reduced hyperspectral reflectance images” by Foster foster2022colour : “Ribeira Houses Shrubs”, abbreviated as “Ribeira”, with size $249\times 330\times 33$; and “Bom Jesus Bush”, abbreviated as “Bush”, with size $250\times 330\times 33$. We transform hyperspectral images to RGB images by following the tutorial888https://personalpages.manchester.ac.uk/staff/david.foster/Tutorial_HSI2RGB/Tutorial_HSI2RGB.html.. Figure 6 shows these two hyperspectral images represented as RGB images.

To evaluate the recovery performance of image completion, peak signal-to-noise ratio (PSNR) is used to measure the similarity of two images, defined by

[TABLE]

where $\max(\mathcal{A})$ denotes the highest pixel value of $\mathcal{A}$, and $\mathrm{MSE}$ is the mean square error defined by $\mathrm{MSE}\mathrel{\mathop{\mathchar 58\relax}}={\|\mathcal{X}-\mathcal{A}\|_{\mathrm{F}}^{2}}/(n_{1}n_{2}n_{3})$. The relative error

[TABLE]

is also reported. In our experiments, the PSNR between a recovered image and original image will be reported after being transformed into RGB images. Given sampling rate $p$, we formulate the sampling set $\Omega$ in the same fashion as subsection 3.3. The initial guess $\mathcal{X}^{(0)}$ is randomly generated with normal distribution $\mathcal{N}(0,1)$. We set the TR rank as $\mathbf{r}=(7,16,7)$, the Tucker rank as (65,65,7) following kressner2014low , the TT rank as $(1,15,15,1)$, and the CP rank as $110$. The search spaces are of a similar size for different tensor formats; see Appendix B. Moreover, the nuclear-norm-based algorithm, HaLRTC liu2012tensor , is compared. The maximum iteration number is 200.

Numerical results for the completion of two hyperspectral images are shown in Fig. 7 and Table 2. One difficulty of recovering the “Ribeira” image is shrubs around the buildings, which triggers different recovery quality of different algorithms. The recovery quality of “Bush” image depends on the details of leaves. Figure 7 displays the recovery results under different sampling rates $p=0.1,0.3,0.5$. The images recovered by TR-based algorithms (TR-RGD, TR-RCG, TR-ALS) depict more details than other candidates, e.g., the shrubs and the bushes. In fact, this observation can be quantified by the PSNR and relative errors; see Table 2. The proposed methods are favorably comparable to TR-ALS. Specifically, we observe that the proposed TR-RCG algorithm reaches the highest PSNR and the lowest relative error among all compared methods in most cases.

5.4 Experiments on high-dimensional functions

The discretization of high-dimensional functions ${h\mathrel{\mathop{\mathchar 58\relax}}[0,1]^{d}\to\mathbb{R}}$ requires storing a large tensor $\mathcal{A}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$, which is unfavorable in practice. To this end, one can apply tensor decomposition and tensor completion to recover the tensor $\mathcal{A}$ by partially observed entries; see applications in glau2020low . In this subsection, we compare TR-based algorithms (TR-RGD, TR-RCG, TR-ALS) with TT-RCG on completion of the data tensor $\mathcal{A}$ generated from $h$ by evenly dividing $[0,1]$ in dimension $k$ of $[0,1]^{d}$ into $n_{k}-1$ intervals for $k\in[d]$. Specifically,

[TABLE]

We consider the following two functions (steinlechner2016riemannian, , Sect. 5.4),

[TABLE]

We set $d=4,\ n_{1}=n_{2}=n_{3}=n_{4}=20$, and the sampling rate $p=0.001,0.005,$ $0.01,0.05,0.1$. We follow the same way in subsection 5.1 to create the sampling set $\Omega$. The rank-increasing strategy is implemented to both TT and TR algorithms by following steinlechner2016riemannian . In order to choose a search space with similar size, we set the maximum rank of TT-RCG as $(1,5,5,5,1)$, and the maximum rank for TR-based algorithms as $(4,4,4,4)$; see Appendix B for details. We adopt the stopping criteria in section 5. Additionally, due to the rank-increasing strategy, an algorithm is also terminated if: 1) maximum iteration number 50 is reached in each fixed-rank searching; 2) the maximum rank is achieved; 3) there is no acceptance of rank increase along any mode for a given point.

The numerical results for recovering high-dimensional functions $h_{1}$ and $h_{2}$ are reported in Table 3. It illustrates that all algorithms have comparable performance. TR-based algorithms perform favorably comparable to TT-RCG. Among all TR-based algorithms, TR-RCG has a better performance in most experiments.

6 Conclusion and perspectives

We have developed Riemannian preconditioned algorithms for the tensor completion problem based on tensor ring decomposition. The preconditioning effect stems from a metric defined on the product space of matrices generated from the mode-2 unfolding of core tensors. However, the straightforward calculation of Riemannian gradient requires large matrix multiplications that are unaffordable in practice. To this end, we have adopted a procedure that efficiently computes the Riemannian gradient without forming large matrices explicitly. The proposed algorithms enjoy global convergence results. Numerical comparisons on both synthetic and real-world datasets present promising results.

Since the rank parameter has to be fixed a priori, we are interested in rank-adaptive strategies that may result in more accurate completion.

Appendix A Speedup of efficient gradient computation

In contrast with naive computation, Algorithm 3 provides an efficient way to compute the gradients. Table 4 reports the speedup results on MovieLens 1M dataset. “Naive” denotes the non-optimized implementation, i.e., explicitly forming $\mathbf{W}_{\neq k}$ and large matrix multiplication. The results show that the proposed algorithm is of great potential to be applied in large-scale problems.

Appendix B The rank selection for numerical experiments

Table 5 introduces the rank selections in numerical experiments. For instance, we initially choose the Tucker rank $(65,65,7)$ following kressner2014low in experiments on hyperspectral images with size $250\times 330\times 33$. In order to ensure a fair comparison, the search spaces for different tensor formats are of a similar size. To this end, we compute the number of parameters by

[TABLE]

For TT decomposition, we select the rank as $(1,15,15,1)$ with

[TABLE]

parameters. For TR decomposition, the rank is selected by $(7,16,7)$. Then, the number of parameters is

[TABLE]

For CP decomposition, the rank is selected by

[TABLE]

Therefore, the search spaces appear to have comparable sizes across various tensor formats. The ranks are chosen using the same procedure in other experiments.

Acknowledgements

We would like to thank the two anonymous reviewers for helpful comments. We acknowledge Shuyu Dong for helpful discussions on the preconditioned metric.

Declaration

The authors declare that the data supporting the findings of this study are available within the paper. The authors have no competing interests to declare that are relevant to the content of this article.