Fast Randomized Algorithms for t-Product Based Tensor Operations and Decompositions with Applications to Imaging Data

TL;DR

This paper introduces fast randomized algorithms for tensor operations based on the t-product, providing efficient low-rank approximations and decompositions with error guarantees, demonstrated on imaging and signal processing applications.

Contribution

It develops novel randomized algorithms for tensor low-rank approximation and decomposition using the t-product, with theoretical error guarantees and practical efficiency.

Findings

Algorithms achieve significant speed-up over traditional methods.

High accuracy in tensor completion tasks.

Effective in imaging applications like mass spectrometry and video recovery.

Abstract

Tensors of order three or higher have found applications in diverse fields, including image and signal processing, data mining, biomedical engineering and link analysis, to name a few. In many applications that involve for example time series or other ordered data, the corresponding tensor has a distinguishing orientation that exhibits a low tubal structure. This has motivated the introduction of the tubal rank and the corresponding tubal singular value decomposition in the literature. In this work, we develop randomized algorithms for many common tensor operations, including tensor low-rank approximation and decomposition, together with tensor multiplication. The proposed tubal focused algorithms employ a small number of lateral and/or horizontal slices of the underlying 3-rd order tensor, that come with {\em relative error guarantees} for the quality of the obtained solutions. The performance of the proposed algorithms is illustrated on diverse imaging applications, including mass spectrometry data and image and video recovery from incomplete and noisy data. The results show both good computational speed-up vis-a-vis conventional completion algorithms and good accuracy.

Tables

Table 1. Table 1 : Relative errors and running times of tensor multipication algorithms.

	t-product	rt-product
Rank, selected slices		non-uniform sampling			uniform sampling
	time (s)	time(s)	RFE	RSE	time(s)	RFE	RSE
r=100, l=460	9.09	8.13	70e-3	43e-4	1.03	39e-2	11e-2
r=200, l=1000	36.36	14.27	64e-3	25e-4	1.08	12e-2	82e-2
r=300, l=1700	162.99	34.89	56e-3	29e-4	1.39	41e-2	34e-2
r=400, l=2500	321.43	41.04	31e-3	13e-3	1.64	90e-2	38e-3
r=500, l=3100	684.01	67.07	71e-3	12e-4	3.09	23e-2	43e-3

Table 2. Table 2 : RSE of matrix and tensor low rank decomposition relative to the best rank − 5 5 -5 space for identifying important ions in the brain data set.

Number of	CX				t-CX
selected slices	Randomized		Deterministic		Randomized		Deterministic
	RSE	time	RSE	time	RSE	time	RSE	time
25	19.13e-2	4.56	18.84e-2	8.46	13.65e-2	4.13	16.05e-2	4.47
35	17.24e-2	5.12	17.59e-2	8.95	17.43e-2	4.87	13.13e-2	5.99
45	16.35e-2	7.01	16.93e-2	14.99	11.32e-2	6.14	11.01e-2	7.01
55	15.14e-2	8.16	16.52e-2	15.86	16.26e-2	6.63	10.16e-2	8.99

Table 3. Table 3 : RSE of matrix and tensor low rank decomposition relative to the best rank − 5 5 -5 space for finding important pixels in the brain data set.

Number of	CX				t-CX
selected slices	Randomized		Deterministic		Randomized		Deterministic
	RSE	time	RSE	time	RSE	time	RSE	time
25	45.24e-2	3.15	46.03e-2	7.41	37.65e-2	3.01	16.05e-2	4.06
35	35.87e-2	3.81	35.68e-2	7.63	28.00e-2	3.19	13.13e-2	4.55
45	24.13e-2	5.29	24.19e-2	11.89	23.98e-2	4.15	21.01e-2	5.59
55	23.16e-2	5.83	24.39e-2	12.93	15.23e-2	4.23	15.17e-2	5.71

Table 4. Table 4 : RSE of matrix and tensor low rank decomposition relative to the best rank − 15 15 -15 space for finding important ions and pixels in the lung data set.

Number of	Decomposition using ions				Decomposition using ions and pixels
selected slices	CX		t-CX		CUR		t-CUR
	RSE	time	RSE	time	RSE	time	RSE	time
25	40.82e-2	3.31	46.89e-2	3.05	59.76e-2	2.89	24.67e-2	2.14
35	41.98e-2	4.13	29.63e-2	3.24	42.57e-2	3.16	22.71e-2	2.67
45	30.16e-2	5.78	22.35e-2	3.36	12.18e-2	4.25	18.15e-2	3.13
55	30.74e-2	6.67	20.81e-2	4.05	19.00e-2	5.47	18.15e-2	3.26

Table 5. Table 5 : RSE of tensor robust completion methods on Extended Yale B dataset.

Robust Completion Approach	RSE	Time(s)
RPCA NN [6]	0.0056	417
R-TUCKER NN [20]	0.0034	513
RPCA t-NN [35]	0.0021	495
CUR t-NN	0.0026	136

Table 6. Table 6 : RSE of tensor completion results for the Basketball video.

Completion Approach	RSE	Time(s)
EXACT NN [7]	0.1001	687
E-TUCKER NN [33]	0.0900	718
EXACT t-NN [53]	0.0715	695
CUR t-NN	0.0850	205

Equations

$$\mathcal{Z}=\sum _{r_1=1}^{R_1}\sum _{r_2=1}^{R_2}\sum _{r_3=1}^{R_3}{g}_{r_1{r}_2{r}_3}\left(u_{r_1}^{(1)}\circ u_{r_2}^{(2)}\circ u_{r_3}^{(3)}\right)=\mathcal{G}{\times }_1{U}^{(1)}{\times }_2{U}^{(2)}{\times }_3{U}^{(3)},$$

$$\hat{X}={\mathrm{bdiag}}_{k\in [n_3]}({\hat{\mathcal{X}}}_{::k}):=\left(\begin{array}{cccc}{\hat{\mathcal{X}}}_{::1}& & & \\ & {\hat{\mathcal{X}}}_{::2}& & \\ & & \ddots & \\ & & & {\hat{\mathcal{X}}}_{::n_3}\\ & & & \end{array}\right).$$

$$\mathcal{C}=\mathcal{Z}\ast \mathcal{X}=\mathrm{fold}\left(\mathrm{circ}\left(\mathcal{Z}\right)\right).\mathrm{unfold}\left(\mathcal{X}\right),$$

$$\mathrm{circ}\left(\mathcal{Z}\right):=\left(\begin{array}{cccc}{\mathcal{Z}}_{::1}& {\mathcal{Z}}_{::n_3}& \dots & {\mathcal{Z}}_{::2}\\ {\mathcal{Z}}_{::2}& {\mathcal{Z}}_{::1}& \dots & {\mathcal{Z}}_{::3}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathcal{Z}}_{::n_3}& {\mathcal{Z}}_{::n_3-1}& \dots & {\mathcal{Z}}_{::1}\end{array}\right),$$

$$\mathrm{unfold}\left(\mathcal{X}\right):=\left(\begin{array}{c}{X}_{::1}\\ X_{::2}\\ ⋮\\ X_{::n_3}\end{array}\right),\mathrm{fold}\left(\mathrm{unfold}\left(\mathcal{X}\right)\right)=\mathcal{X}.$$

$$\mathcal{X}=\mathcal{U}\ast \Sigma \ast {\mathcal{V}}^T,$$

$$(F_{n_3}\otimes I_{n_1})\cdot \mathrm{circ}\left(\mathcal{X}\right)\cdot (F_{n_3}^{-1}\otimes I_{n_2})=\hat{X},$$

$$r=\text{card}\{i:{\Sigma }_{ii:}\ne 0\}=\underset{i}{\max }{r}_i,$$

$${\mathrm{rank}}_t(\mathcal{X}\ast \mathcal{Z})\le \min ({\mathrm{rank}}_t(\mathcal{X}),{\mathrm{rank}}_t(\mathcal{Z})).$$

$$\Vert \mathcal{X}{\Vert }_{⊛}:=\frac{1}{n_3}\sum _{k=1}^{n_3}\Vert {\hat{\mathcal{X}}}_{::k}{\Vert }_{\ast }.$$

$$\Vert \mathcal{X}{\Vert }_{⊛}$$

$$\Vert \mathcal{X}\Vert :=\Vert \hat{X}{\Vert }_2=\Vert (F_{n_3}\otimes I_{n_1})\cdot \mathrm{circ}\left(\mathcal{X}\right)\cdot (F_{n_3}^{-1}\otimes I_{n_2}){\Vert }_2,$$

$${\mathcal{X}}^{-1}\ast \mathcal{X}=\mathcal{X}\ast {\mathcal{X}}^{-1}=\mathcal{I},$$

$${\rm Y}:=\{{\mathcal{X}}_{:j:}\mid {\mathcal{X}}_{:j:}\in {\mathbb{C}}^{n_1\times 1\times n_3},j\in [n_2]\},$$

$$\sum _{j=1}^m{\mathcal{X}}_{:j:}\ast {\mathcal{C}}_{jj:}=O,$$

$$\sum _{j=1}^n{\mathcal{X}}_{:j:}\ast {\mathcal{C}}_{jj:}=O,$$

$$\underset{\mathcal{L}}{\min }\Vert \mathcal{L}{\Vert }_{⊛}s.t.\Vert P_{\Omega }(\mathcal{X}-\mathcal{L})\Vert \le \Delta ,$$

$$(P_{\Omega }(\mathcal{X}{)}_{ij}$$

$${\mu }_0(\mathcal{V})$$

$${\mu }_1(\mathcal{L}):=\frac{n_1{n}_2{n}_3^2}{r}\Vert \mathcal{U}\ast {\mathcal{V}}^T{\Vert }_{\infty }^2,$$

$${\text{rank}}_t(\mathcal{L})$$

$$\max \{{\mu }_0(\mathcal{U}),{\mu }_0(\mathcal{V})\}$$

$${\mu }_1(\mathcal{L})$$

$$\mathcal{A}\ast \mathcal{B}=\sum _{t=1}^{n_2}{\mathcal{A}}_{:t:}\ast {\mathcal{B}}_{t::},$$

$$\mathcal{A}\ast \mathcal{B}\sim \mathcal{C}\ast \mathcal{R}.$$

$$p_i\ge \beta \frac{\Vert {\mathcal{A}}_{:i:}{\Vert }_F\Vert B_{i::}{\Vert }_F}{{\sum }_{i=1}^{n_2}\Vert {\mathcal{A}}_{:i:}{\Vert }_F\Vert {\mathcal{B}}_{i::}{\Vert }_F},\beta \in (0,1].$$

$$p_i\ge \beta \frac{\Vert {\mathcal{A}}_{:i:}{\Vert }_F^2}{\Vert \mathcal{A}{\Vert }_F^2},\beta \in (0,1].$$

$$\mathbf{E}[\Vert \mathcal{A}\ast \mathcal{B}-\mathcal{C}\ast \mathcal{R}{\Vert }_F]\le \frac{1}{\sqrt{\beta c}}\Vert \mathcal{A}{\Vert }_F\Vert \mathcal{B}{\Vert }_F.$$

$$\Vert \mathcal{A}\ast {\mathcal{A}}^{\top }-\mathcal{C}\ast {\mathcal{C}}^{\top }\Vert$$

$${\Pi }_{\mathcal{C}}(\mathcal{X}):=\mathcal{C}\ast {\mathcal{C}}^{†}\ast \mathcal{X},$$

$${\Pi }_{\mathcal{C}}(\mathcal{X})=\mathcal{C}\ast {\mathcal{C}}^{†}\ast \mathcal{X},$$

Full text

{AMS}

15A69, 15A23, 68W20, 65F30

1 Introduction

Tensors are multi-dimensional arrays that have been used in diverse fields of applications, including chemometrics [45], psychometrics [30], image/video and signal processing [24, 25, 35, 44, 54] and link analysis [29]. They have also been the object of intense mathematical study (see for example the review paper by Kolda and Bader [28] and references therein).

Analogously to the matrix case, a number of tensor decompositions have been proposed in the literature, briefly described next for the case of 3-way tensors. Let $\boldsymbol{\mathscr{{Z}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$. Then, there exists a factorization, called a “Tucker decomposition”of $\boldsymbol{\mathscr{{Z}}}$, of the form

[TABLE]

where $\boldsymbol{\mathscr{{G}}}\in\mathbb{R}^{R_{1}\times R_{2}\times R_{3}}$ is the core tensor associated with this decomposition, $U^{(i)}\in\mathbb{R}^{n_{i}\times R_{i}}$ is the $i$-th factor matrix and $\times_{i}$ is the mode$-i$ product for $i=1,2,3$. The $3$-tuple $(R_{1},R_{2},R_{3})$ is called the Tucker rank or multilinear rank of tensor $\boldsymbol{\mathscr{{Z}}}$ [8, 22]. The conventional Tucker decomposition corresponds to the orthonormal Tucker decomposition, which is also known as the higher-order SVD (HOSVD). De Lathauwer et al. [10] proposed an algorithm to compute a HOSVD decomposition. Soon afterwards they proposed the higher-order orthogonal iteration (HOOI) [11] to provide an inexact Tucker decomposition. The CP decomposition of a tensor is another important notion of tensor-decomposition, which leads to the definition of CP rank. The CP model can be considered as a special case of the Tucker model with a superdiagonal core tensor. Further, the CP rank of a tensor equals that of its Tucker core [21].

1.1 Third order tensor as operator on matrices

While the Tucker-based factorization may be sufficient for many applications, in this paper we consider an entirely different tensor decomposition based on circulant algebra [26]. In this factorization, a tensor in $\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is viewed as a $n_{1}\times n_{2}$ matrix of “tubes”, also known as elements of the ring $\mathbb{R}^{n_{3}}$ where addition is defined as vector addition and multiplication as circular convolution. This “matrix-of-tubes” viewpoint leads to definitions of a new multiplication for tensors (“tubal multiplication”), a new rank for tensors (“tubal rank”), and a new notion of a singular value decomposition (SVD) for tensors (“tubal SVD”). The tubal SVD (t-SVD) provides the “best” tubal rank-$r$ approximation to $\boldsymbol{\mathscr{{Z}}}$, as measured with respect to any unitary invariant tensor norm.

A limitation of the t-SVD decomposition is that it depends directly on the orientation of the tensor, whereas the CP and Tucker decompositions are not. This suggests that the latter decompositions are well suited for data applications where the tensor’s orientation is not critical - e.g. chemometrics [45] and/or psychometrics [30]. However, in applications involving time series or other ordered data, the orientation of the tensor is fixed. Examples include, but not limited to, computed tomography (CT) [44], facial recognition [18] and video compression [54], where the tensor decomposition is dependent on the third dimension. Analogous to compressing a two-dimensional image using the matrix SVD (a classic linear algebra example, with detailed writeup in [23]), the t-SVD decomposition can be used to compress several images taken over time (e.g. successive frames from a video). Since such images do not change substantially from frame to frame, we expect tubal compression strategies to provide better results than performing a matrix SVD on each separate image frame. The former consider the tensor as a whole object, rather than as a collection of distinct images [27, 26, 54, 44, 35]. Further, t-SVD is essentially based on a group theoretic approach, where the multidimensional structure is unraveled by constructing group-rings along the tensor fibers 111We consider the group rings constructed out of cyclic groups, resulting in an algebra of circulants. However, the results presented in this paper hold true for the general group-ring construction.. The advantage of such an approach over existing ones is that the resulting algebra and corresponding analysis enjoys many similar properties to matrix algebra and analysis. For example, it is shown in [20] that recovering a 3-way tensor of length $n$ and Tucker rank $(r,r,r)$ from random measurements requires $O(rn^{2}\log^{2}(n))$ observations under a matrix coherence condition on all mode-$n$ unfoldings. However, the number of samples needed for exact recovery of a 3-way tensor of length $n$ and tensor tubal-rank $r$ is $O(rn^{2}\log(n^{2}))$ under a weaker tensor coherence assumption [53]. Further, consider the decomposition $\boldsymbol{\mathscr{{X}}}=\boldsymbol{\mathscr{{L}}}+\boldsymbol{\mathscr{{E}}}$, where $\boldsymbol{\mathscr{{L}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is low rank and $\boldsymbol{\mathscr{{E}}}$ is sparse. Let $n_{(1)}=\max(n_{1},n_{2})$ and $n_{(2)}=\min(n_{1},n_{2})$. The work in [35] shows that for tensor $\boldsymbol{\mathscr{{L}}}$ with coherence parameter $\mu$, the recovery is guaranteed with high probability for the tubal rank of order $n_{(2)}n_{3}/(\mu(\log n_{(1)}n_{3})^{2})$ and a number of nonzero entries in $\boldsymbol{\mathscr{{E}}}$ of order $O(n_{1}n_{2}n_{3})$. Hence, under the same coherence condition (see, Definitions 2.15 and 2.16), the tubal robust tensor factorization problem perfectly recovers the low-rank and sparse components of the underlying tensor.

A shortcoming of these three classical decompositions is their brittleness with respect to severely corrupted or outlier data entries. To that end, a number of approaches have been developed in the literature to recover a low-rank tensor representation from data subject to noise and corrupted entries. We focus on two instances of the problem based on the t-SVD algorithm: (i) noisy tensor completion, i.e., recovering a low-rank tensor from a small subset of noisy entries, and (ii) noisy robust tensor factorization, i.e., recovering a low-rank tensor from corrupted data by noise and/or outliers of arbitrary magnitude [35, 54]. These two classes of tensor factorization problems have attracted significant interest in the research community [2, 6, 44, 54]. In particular, convex formulations of noisy tensor factorization have been shown to exhibit strong theoretical recovery guarantees and a number of algorithms has been developed for solving them [54, 44, 35].

It is frequently mentioned that (noisy) tensor factorization, despite its numerous advantages, also exhibits a number of drawbacks listed below:

• •

The available methods [4, 27, 26, 54, 44, 35] are inherently sequential and all rely on the repeated and costly computation of t-SVD factors, Discrete Fourier Transform (DFT) and its inverse (IDFT), that limit their scalability.

• •

The basis tensor vectors resulting from t-SVD have little concrete meaning, which makes it challenging for practitioners to interpret the obtained results. For instance, the vector [$(1/2)$ age - $(1/\sqrt{2})$ height + $(1/2)$income], being one of the significant uncorrelated factors from a data set of people’s features is not easily interpretable (see discussion in [15]). Kuruvilla et al. in [31] have also claimed: “it would be interesting to try to find basis vectors for all experiment vectors, using actual experiment vectors and not artificial bases that offer little insight”.

• •

The t-SVD decomposition for sparse tensors does not preserve sparsity in general, which for large size tensors leads to excessive computations and storage requirements. Hence, it is important to compute low-rank tensor factorizations that preserve such structural properties of the original data tensor.

1.2 Main Contributions

In this paper, we study scalable randomized tensor multiplication (rt-product algorithm) and tensor factorization (rt-project) operations and extend the matrix CX and CUR type decompositions [12, 13, 15, 36] to third order tensors using a circulant algebra embedding. To that end, we develop a basic algorithm (t-CX), together with a more general one (t-CUR) based on tensor slice selection that come with relative error guarantees. Finally, we propose a new tensor nuclear norm minimization method, called CUR t-NN, which solves the noisy tensor factorization problem using a small number of lateral and/or horizontal slices of the underlying tensor. Specifically, CUR t-NN uses an adaptive technique to sample slices of the tensor based on a leverage score for them and subsequently solves a convex optimization problem followed by a projection step to recover a low rank approximation to the underlying tensor. Advantages of CUR t-NN include:

• •

Contrary to nuclear norm minimization based approaches, which minimize the sum of the singular values of the underling tensor, our approach only minimizes the sum of singular values of the set of sampled slices corresponding to the largest leverage scores. Thus, we obtain a more accurate and robust approximation to the rank function.

• •

Using subspace sampling, we obtain provable relative-error recovery guarantees for tubal product based tensor factorization. This technique is likely to be useful for other tensor approximation and data analysis problems.

• •

The proposed algorithm for noisy tensor factorization tasks has polynomial time complexity.

1.3 Related work on randomized tensor factorization

Note that there has been prior work on the topic of randomized methods for tensor low rank approximations/decompositions. For example, recent work includes [37, 5]. Our methods are different from these studies, since we rely on the t-product construct in [4, 26, 27] to provide tubal tensor multiplication and low rank approximation. In [37], the proposed tensor-CUR algorithm employs a CUR matrix approximation to one of the unfolding matrix modes (the distinguished mode) providing an approximation based on few tube fibers and few slices. Hence, that algorithm achieves an additive error guarantee in terms of the flattened (unfolded) tensor, rather than the original one. However, the algorithms presented in this work offer relative error guarantees for the obtained approximations, in terms of the original low tubal rank tensor. Further, we propose a new tensor nuclear norm minimization method, which solves the noisy tensor factorization problem in the fully and partially observed setting using a small number of lateral and/or horizontal slices of the underlying tensor. It is worth mentioning that randomized algorithms were used to efficiently solve the robust matrix PCA problem using small sketches constructed from random linear measurements of the low rank matrix [43, 55, 34, 36, 39]. Our proposed nuclear norm minimization approach is different from these studies, since we relay on slice selection and projection (mainly t-CUR factorization). Further, our proposed algorithm uses an adaptive sampling strategy and provides high probability recovery guarantees of the exact low rank tensor approximation under a weaker coherence assumption.

The remainder of the paper is organized as follows. In Section 2, we review some relevant mathematical concepts including the tensor circulant algebra, basic definitions and theorems of tensors, and the t-SVD decomposition based on the t-product concept. In Section 3, we provide the randomized tensor decompositions and provide their relative error guarantees and introduce the concept of slice selection and projection. In Section 4, we introduce, evaluate and analyze our proposed algorithm (CUR t-NN) for large scale noisy tensor decomposition. Experimental results are presented in Section 5 and some concluding remarks are drawn in Section 6.

The detailed proofs of the main results established are delegated to the Appendix.

2 Mathematical preliminaries

Next, we introduce key definitions and concepts used in subsequent developments.

Tensor indexing.

We denote tensors by boldface Euler script letters, e.g., $\mathscr{X}$, matrices by boldface capital letters, e.g., $X$, vectors by lowercase letters, e.g., $x$. The order of a tensor is the number of dimensions (also refereed to as ways or modes). In this work, we focus on 3-way tensors.

Fibers and slices [51].

A fiber of tensor $\mathscr{X}$ is a one-dimensional array defined by fixing two of the indices. Specifically, $\mathscr{X}_{:jk}$ is the (j,k)-th column fiber, $\mathscr{X}_{i:k}$ is the (i,k)-th row fiber, and $\mathscr{X}_{ij:}$ is the (i,j)-th tube fiber. A slice of tensor $\mathscr{X}$ is a two-dimensional array defined by fixing one index only. Specifically, $\mathscr{X}_{i::}$ is the i-th horizontal slice, $\mathscr{X}_{:j:}$ is the j-th lateral slice, and $\mathscr{X}_{::k}$ is the k-th frontal slice. For convenience, $\mathscr{X}_{::k}$ is written as $\mathscr{X}_{k}$. The vectorization of $\mathscr{X}$ is denoted by $vec(\mathscr{X})$. For a 3-way tensor $\mathscr{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, we denote its $(i,j,k)$-th entry as $\mathscr{X}_{ijk}$.

Norms.

We denote the $\ell_{1}$ norm as $\|\boldsymbol{\mathscr{{X}}}\|_{1}:=\sum_{ijk}|x_{ijk}|$, the infinity norm as $\|\boldsymbol{\mathscr{{X}}}\|_{\infty}:=\max_{ijk}|x_{ijk}|$, and the Frobenius norm as $\|\boldsymbol{\mathscr{{X}}}\|_{F}:=\sqrt{\sum_{ijk}|x_{ijk}|^{2}}$. The above norms reduce to the vector or matrix norms if $\boldsymbol{\mathscr{{X}}}$ is a vector or a matrix.

Operators.

For $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, and using the Matlab commands fft,ifft, we denote by $\hat{\boldsymbol{\mathscr{{X}}}}$ the result of applying the Discrete Fourier Transform (DFT) on $\boldsymbol{\mathscr{{X}}}$ along the 3-rd dimension, i.e., $\hat{\boldsymbol{\mathscr{{X}}}}=\texttt{fft}(\boldsymbol{\mathscr{{X}}},[],3)$. Analogously, one can also compute $\boldsymbol{\mathscr{{X}}}$ from $\hat{\boldsymbol{\mathscr{{X}}}}$ via the Inverse Discrete Fourier Transform (IDFT), using $\texttt{ifft}(\hat{\boldsymbol{\mathscr{{X}}}},[],3)$. In particular, we denote by $\hat{X}$ the block diagonal matrix with each blockcorresponding to the frontal slice $\hat{X}_{::k}$ of $\hat{\boldsymbol{\mathscr{{X}}}}$. i.e.,

[TABLE]

2.1 Tensor basics

Next, we review relevant mathematical concepts including the tensor SVD (t-SVD), basic definitions and operations and other technical results of tensors [27, 44, 53], that are used throughout the paper.

Definition 2.1**.**

(Tensor product). Given two tensors $\boldsymbol{\mathscr{{Z}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{2}\times n_{4}\times n_{5}}$, the t-product $\boldsymbol{\mathscr{{Z}}}*\boldsymbol{\mathscr{{X}}}$ is the $n_{1}\times n_{4}\times n_{3}$ tensor,

[TABLE]

where

[TABLE]

and

[TABLE]

Because the circular convolution of two tube fibers can be computed by the DFT, the t-product can be alternatively computed in the Fourier domain, as shown in Algorithm 1.

Definition 2.2**.**

**(Conjugate transpose). The conjugate transpose of a tensor $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is tensor $\boldsymbol{\mathscr{{X}}}^{*}\in\mathbb{R}^{n_{2}\times n_{1}\times n_{3}}$ obtained by conjugate transposing each of the frontal slices and then reversing the order of transposed frontal slices 2 through $n_{3}$. **

Definition 2.3**.**

**(Identity tensor). The identity tensor $\boldsymbol{\mathscr{{I}}}\in\mathbb{R}^{n\times n\times n_{3}}$ is the tensor whose first frontal slice is the $n\times n$ identity matrix, and whose other frontal slices are all zeros. **

Definition 2.4**.**

**(Orthogonal tensor). A tensor $\boldsymbol{\mathscr{{Q}}}\in\mathbb{R}^{n\times n\times n_{3}}$ is orthogonal if it satisfies $\boldsymbol{\mathscr{{Q}}}^{*}*\boldsymbol{\mathscr{{Q}}}=\boldsymbol{\mathscr{{Q}}}*\boldsymbol{\mathscr{{Q}}}^{*}=\boldsymbol{\mathscr{{I}}}.$ **

Definition 2.5**.**

**(F-diagonal Tensor). A tensor is called $f-$diagonal if each of its frontal slices is a diagonal matrix. **

Theorem 2.6**.**

(t-SVD). Let $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$. Then, it can be factored as

[TABLE]

*where $\boldsymbol{\mathscr{{U}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, $\boldsymbol{\mathscr{{V}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ are orthogonal and $\Sigma\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is a f-diagonal tensor. *

Note that t-SVD can be efficiently computed based on the matrix SVD in the Fourier domain. This is based on the key property that the block circulant matrix can be mapped to a block diagonal matrix in the Fourier domain, i.e.

[TABLE]

where $F_{n_{3}}$ denotes the $n_{3}\times n_{3}$ DFT matrix and $\otimes$ denotes the Kronecker product.

Definition 2.7**.**

(Tensor multi and tubal rank). The tensor multi-rank of $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is a vector $\upsilon\in\mathbb{R}^{n_{3}}$ with its $i-th$ entry being the rank of the i-th frontal slice of $\hat{\boldsymbol{\mathscr{{X}}}}$, $r_{i}=\operatorname{rank}(\hat{\boldsymbol{\mathscr{{X}}}}_{i::})$. The tensor tubal rank, denoted by $r=\operatorname{rank}_{t}(\boldsymbol{\mathscr{{X}}})$, is defined as the number of nonzero singular tubes of $\Sigma$, where $\Sigma$ is obtained from the t-SVD of $\boldsymbol{\mathscr{{X}}}=\boldsymbol{\mathscr{{U}}}\ast\Sigma\ast\boldsymbol{\mathscr{{V}}}^{T}$. That is,

[TABLE]

*where card denotes the cardinality of a set. *

The tensor tubal rank shares some properties of the matrix rank, e.g.

[TABLE]

Many tensor completion and decomposition techniques for video and seismic noise reduction rely on a low-rank factorization of a time-frequency transform. Further, certain energy methods broadly used in image processing, e.g., PDEs [1] and belief propagation techniques mainly focus on local relationships. The basic assumption is that the missing entries depend primarily on their neighbors. Hence, the further apart two pixels are, the smaller their dependance is. However, for video and time series of images the value of the missing entry also depends on entries which are relatively far away in the time/sequence dimension. Thus, it is necessary to develop a tool to directly capture such global information in the data. Using (5), $r_{0}$ is the rank of the “mean image” across the video sequence. Meanwhile, $r_{1}$ is the rank of the next frequency’s content across frames, etc. Under a smoothness assumption that captures global information at given pixels across time, the frontal slices (after FFT) for bigger $i$ have smaller singular values.

Definition 2.8**.**

(Tensor nuclear norm). The tensor nuclear norm of a tensor $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, denoted by $\|\boldsymbol{\mathscr{{X}}}\|_{\circledast}$, is defined as the average of the nuclear norm of all frontal slices of $\hat{\boldsymbol{\mathscr{{X}}}}$, i.e.

[TABLE]

The above tensor nuclear norm is defined in the Fourier domain. It is closely related to the nuclear norm of the block circulant matrix in the original domain. Indeed,

[TABLE]

The above relationship gives an equivalent definition of the tensor nuclear norm in the original domain. Thus, the tensor nuclear norm is the nuclear norm (with a factor $1/n_{3}$) of a new matricization (block circulant matrix) of a tensor.

Definition 2.9**.**

(Tensor spectral norm). The tensor spectral norm of $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n\times n\times n_{3}}$, denoted as $\|\boldsymbol{\mathscr{{X}}}\|$, is defined as

[TABLE]

*where $\|\cdot\|_{2}$ denotes the spectral norm of a matrix. *

Definition 2.10**.**

(Inverse of tensor). The inverse of $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n\times n\times n_{3}}$, denoted by $\boldsymbol{\mathscr{{X}}}^{-1}$ satisfies

[TABLE]

*where $\boldsymbol{\mathscr{{I}}}$ is the identity tensor of size $n\times n\times n_{3}$. *

Definition 2.11**.**

**(Standard tensor basis). The lateral basis $\boldsymbol{\mathring{\mathfrak{{e}}}}_{i}$, is of size $n_{1}\times 1\times n_{3}$ with only one entry equal to 1 and the remaining equal to zero, in which the nonzero entry 1 will only appear at the first frontal slice of $\boldsymbol{\mathring{\mathfrak{{e}}}}_{i}$. Normally its transpose $\boldsymbol{\mathring{\mathfrak{{e}}}}_{i}^{\top}$ is called the horizontal basis. The other standard tensor basis is called tube basis $\boldsymbol{\dot{\mathfrak{{e}}}}_{i}$, and corresponds to a tensor of size $1\times 1\times n_{3}$ with one entry equal to 1 and the rest equal to 0. Figure 2 illustrates these bases. **

2.2 Linear algebra with tensors: free submodules

The set of complex numbers $\mathbb{C}$ with standard scalar addition and multiplication is a field and $\mathbb{C}^{n_{3}}$ forms a vector space over this field. However, as pointed out in [25, 27], the set of tubes $\mathbb{C}^{1\times 1\times n_{3}}$ equipped with the tensor product form a ring with unity. A module over a ring can be thought of as a generalization of a vector space over a field, where the corresponding scalars are the elements of the ring. In linear algebra over a ring, the analog of a subspace is a free submodule. Our algorithm relies on submodules, and the following theorem.

Theorem 2.12**.**

The set of slices

[TABLE]

*forms a free module over the set of tubes $\mathbb{C}^{1\times 1\times n_{3}}$. *

Proof 2.13**.**

*A detailed proof is provided in [4]. *

Using Theorem 2.12, $\Upsilon$ has a basis so that any element of it can be written as a “t-linear combination” of elements of a set of basis slices. A “t-linear combination”, is defined as a sum of slices multiplied, based on the t-product, by coefficients from $\mathbb{C}^{1\times 1\times n_{3}}$.

Definition 2.14**.**

(Slice-wise linear independence) The slices in a subset $\Lambda=\{\boldsymbol{\mathscr{{X}}}_{:1:},\dots,\boldsymbol{\mathscr{{X}}}_{:n:}\}$ of $\Upsilon$ are said to be linearly dependent, if there exists a finite number of distinct slices $\boldsymbol{\mathscr{{X}}}_{:1:},\dots,\boldsymbol{\mathscr{{X}}}_{:m:}$ in $\Lambda$, and tubes $\boldsymbol{\mathscr{{C}}}_{11:},\dots,\boldsymbol{\mathscr{{C}}}_{mm:}$, not all zero, such that

[TABLE]

where $O$ denotes the lateral slice comprising of all zeros.

The slices in a subset $\Lambda=\{\boldsymbol{\mathscr{{X}}}_{:1:},\dots,\boldsymbol{\mathscr{{X}}}_{:n:}\}$ of $\Upsilon$ are said to be linearly independent if the equation

[TABLE]

*can only be satisfied by all zero tubes $\boldsymbol{\mathscr{{C}}}_{jj:},~{}j=1,\dots,n$. *

Based on the computable t-SVD, the tensor nuclear norm [44] is used to replace the tubal rank for low-rank tensor recovery (from incomplete/corrupted tensors) by solving the following convex program,

[TABLE]

where $\|\boldsymbol{\mathscr{{L}}}\|_{\circledast}$ denotes the tensor nuclear norm and

[TABLE]

Similarly to the matrix completion problem, recovery of tensor $\boldsymbol{\mathscr{{X}}}$ from its observed entries is essentially infeasible if the large majority of the entries are equal to zero [7]. For the tensor completion case, it is the case that if tensor $\boldsymbol{\mathscr{{X}}}$ only has a few entries which are not equal to zero, in its t-SVD $\boldsymbol{\mathscr{{U}}}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{V}}}^{\top}=\boldsymbol{\mathscr{{X}}}$, the singular tensors $\boldsymbol{\mathscr{{U}}}$ and $\boldsymbol{\mathscr{{V}}}$ will be highly concentrated. Indeed, not all tensors can be recovered from data sets with missing entries and/or large outliers. Hence, in analogy to the main idea in matrix completion [7], tensor slices $\boldsymbol{\mathscr{{U}}}(:,i,:)$ and $\boldsymbol{\mathscr{{V}}}(:,i,:),i=1,2,...,r$, need to be sufficiently spread out, which in turn implies that they should be uncorrelated with the standard tensor basis. Our analysis in Section 4 will focus on noisy tensor completion based on a robust factorization algorithm whose estimation/recovery guarantees are expressed in terms of the coherence of the target low-rank tensor $\boldsymbol{\mathscr{{L}}}$. It establishes that lower values of tensor coherence provide better recovery results. In addition, we propose a randomized approximation algorithm, whose guarantees are also related to the notion of tensor coherence. The following three notions of tensor coherence are defined next.

Definition 2.15**.**

(Tensor $\mu_{0}$-coherence). Let $\boldsymbol{\mathscr{{V}}}\in\mathbb{R}^{n\times r\times n_{3}}$ contain orthonormal lateral basis with $r\leq n$. Then, the $\mu_{0}$-coherence of $\boldsymbol{\mathscr{{V}}}$ is given by:

[TABLE]

*where $\boldsymbol{\mathring{\mathfrak{{e}}}}_{i}$ is standard lateral basis. *

Definition 2.16**.**

(Tensor $\mu_{1}$-coherence). For a tensor $\boldsymbol{\mathscr{{L}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ assume that $\text{rank}_{t}(\boldsymbol{\mathscr{{L}}})=r$. Then, the $\mu_{1}$-coherence of $\boldsymbol{\mathscr{{L}}}$ is defined as

[TABLE]

*where $\boldsymbol{\mathscr{{L}}}$ has the skinny t-SVD $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{V}}}^{\top}$. *

Definition 2.17**.**

(Tensor $(\mu,r)$-coherence). For any $\boldsymbol{\mathscr{{\mu}}}>0$, we call a tensor $\boldsymbol{\mathscr{{L}}}$ $(\mu,r)$-coherent if

[TABLE]

Note that the standard tensor coherence condition is much weaker than the matrix weak coherence one for each frontal slice of $\hat{\boldsymbol{\mathscr{{X}}}}$ [53]. Hence, in the analysis of the proposed randomized algorithms, we will use the standard tensor coherence condition.

3 Approximate Tensor Low Rank Decompositions

We develop extensions of the matrix CX and CUR decompositions [13, 15] to third-order tensors. The proposed algorithms result in computing low-rank tensor approximations that are explicitly expressed in terms of a small number of slices of the input tensor. We start by introducing a basic computational tool -a randomized tensor multiplication procedure- that is used in subsequent developments.

3.1 Approximate tensor multiplication

Next, we present a randomized tensor-product and provide key results on the quality of the resulting approximation. Given two tensors $\boldsymbol{\mathscr{{A}}}$ and $\boldsymbol{\mathscr{{B}}}$, using Definition 2.1, the t-product may be written as follows:

[TABLE]

where $*$ denotes the tensor product.

It can be easily shown that the left hand side of (11) is equivalent to the block multiplication and then summation in the Fourier domain. When tensor multiplication is formulated as (11), we can develop a randomized algorithm to approximate the product $\boldsymbol{\mathscr{{A}}}*\boldsymbol{\mathscr{{B}}}$.

Algorithm 2 takes as input two tensors $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and $\boldsymbol{\mathscr{{B}}}\in\mathbb{R}^{n_{2}\times n_{4}\times n_{3}}$, a positive integer $c\leq n_{2}$, and a probability distribution $\{p_{i}\}_{i=1}^{n_{2}}$ over $[n_{2}]$. It returns as output two tensors $\boldsymbol{\mathscr{{C}}}$ and $\boldsymbol{\mathscr{{R}}}$, where the lateral slices of $\boldsymbol{\mathscr{{C}}}$ correspond to a small number of sampled and rescaled slices of $\boldsymbol{\mathscr{{A}}}$, and similarly the horizontal slices of $\boldsymbol{\mathscr{{R}}}$ constitute a small number of sampled and rescaled slices of $\boldsymbol{\mathscr{{B}}}$. Specifically, consider at most $c$ lateral slices (in expectation) of $\boldsymbol{\mathscr{{A}}}$ selected, with the $i$-th lateral slice of $\boldsymbol{\mathscr{{A}}}$ in $\boldsymbol{\mathscr{{C}}}$ be chosen with probability $p_{i}=\min\{1,cp_{i}\}$. Then, define the sampling tensor $\boldsymbol{\mathscr{{S}}}\in\mathbb{R}^{n_{2}\times n_{2}\times n_{3}}$ to be binary where $S_{ii1}=1$ if the $i$-th slice is selected and $\boldsymbol{\mathscr{{S}}}_{ij2:n_{3}}=0$ otherwise. Define the rescaling tensor $\boldsymbol{\mathscr{{D}}}\in\mathbb{R}^{n_{2}\times c\times n_{3}}$ to be the tensor with $\boldsymbol{\mathscr{{D}}}_{ij1}=1/\sqrt{cp_{j}}$, if $i-1$ of the previous slices have been selected and $\boldsymbol{\mathscr{{D}}}_{ij2:n_{3}}=0$ otherwise. Analogously, the rt-product samples and rescales the corresponding horizontal slices of tensor $\boldsymbol{\mathscr{{B}}}$. In both cases, $\boldsymbol{\mathscr{{C}}}=\boldsymbol{\mathscr{{A}}}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}$ is an $n_{1}\times c\times n_{3}$ tensor consisting of sampled and rescaled copies of the lateral slices of $\boldsymbol{\mathscr{{A}}}$, and $\boldsymbol{\mathscr{{R}}}=(\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}})^{\top}*\boldsymbol{\mathscr{{B}}}=\boldsymbol{\mathscr{{D}}}*\boldsymbol{\mathscr{{S}}}^{\top}*\boldsymbol{\mathscr{{B}}}$ is a $c\times n_{4}\times n_{3}$ tensor comprising of sampled and rescaled copies of the horizontal slices of $\boldsymbol{\mathscr{{B}}}$. For the case of $n_{3}=1$, the algorithm selects column-row pairs in Algorithm 2 and is identical to the algorithm in [13].

3.2 Running time of rt-product

The rt-product is computationally efficient, has small memory requirements, and is well suited for large scale problems. Indeed, the rt-product algorithm can be implemented without storing tensors $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ and $\boldsymbol{\mathscr{{B}}}\in\mathbb{R}^{n_{2}\times n_{4}\times n_{3}}$ and uses $O(\max(n_{1},n_{4})cn_{3}\log(n_{3}))$ flops to transform the input to the Fourier domain and $O(c(n_{1}+n_{2}+n_{4})n_{3})$ flops to construct $\boldsymbol{\mathscr{{C}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$ and $\boldsymbol{\mathscr{{R}}}\in\mathbb{R}^{c\times n_{4}\times n_{3}}$, where

[TABLE]

It is worth mentioning that we do not require storing the sampling and rescaling tensors $\boldsymbol{\mathscr{{S}}}$ and $\boldsymbol{\mathscr{{D}}}$ in our implementation.

Next, we provide the main result on the quality of the approximation obtained by Algorithm 2 that specifies the assumptions under which (12) holds. The most interesting of these assumptions is that the sampling probabilities used to randomly sample the lateral slices of $\boldsymbol{\mathscr{{A}}}$ and the corresponding horizontal slices of $\boldsymbol{\mathscr{{B}}}$ are non-uniform and depend on the product of the norms of the lateral slices of $\boldsymbol{\mathscr{{A}}}$ and/or the corresponding horizontal slices of $\boldsymbol{\mathscr{{B}}}$. Following [13], we consider two examples of nonuniform sampling probabilities:

i.

If we would like to use information from both tensors $\boldsymbol{\mathscr{{A}}}$ and $\boldsymbol{\mathscr{{B}}}$, we consider sampling probabilities $\left\{p_{i}\right\}_{i=1}^{n_{2}}$ such that

[TABLE]

ii.

If only information on $\boldsymbol{\mathscr{{A}}}$ is easily available, then we use sampling probabilities $\left\{p_{i}\right\}_{i=1}^{n_{2}}$ such that

[TABLE]

The following theorem gives the main result for Algorithm 2, and generalizes Theorem 6 in [15] to third order tensors .

Theorem 3.1**.**

Suppose $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, $\boldsymbol{\mathscr{{B}}}\in\mathbb{R}^{n_{2}\times n_{4}\times n_{5}}$, and $c\leq n_{2}$. In Algorithm 2, if the sampling probabilities $\left\{p_{i}\right\}_{i=1}^{n_{2}}$ satisfy (13) or (14), then the following holds

[TABLE]

The following lemma establishes a bound for tensor multiplication with respect to the spectral norm with improved sampling complexity.

Lemma 3.2**.**

Given a tensor $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, choose $c\geq 48r\log(4r/(\beta\delta))/(\beta\epsilon^{2})$. Let $\boldsymbol{\mathscr{{C}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$ and $\boldsymbol{\mathscr{{R}}}=\boldsymbol{\mathscr{{C}}}^{\top}$ be the corresponding subtensors of $\boldsymbol{\mathscr{{A}}}$. If sampling probabilities $\left\{p_{i}\right\}_{i=1}^{n_{2}}$ satisfy (14), then with probability at least $1-\delta$, we have

[TABLE]

*where $k_{\pi}=k\in[n_{3}]$ is the index of the tensor’s frontal slice with the maximum spectral norm $\|\hat{\boldsymbol{\mathscr{{A}}}}_{::k}\hat{\boldsymbol{\mathscr{{A}}}}_{::k}^{\top}-\hat{\boldsymbol{\mathscr{{C}}}}_{::k}\hat{\boldsymbol{\mathscr{{C}}}}_{::k}^{\top}\|_{2}$. *

3.3 Slice-based tensor CX decomposition

Next, using the concept of free submodules and Definition 2.14, we introduce the novel notion of slice selection and projection before formulating a tensor CX decomposition.

Definition 3.3**.**

Let $\boldsymbol{\mathscr{{C}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$ and $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$. Then, the tensor project (t-project) operator for $\boldsymbol{\mathscr{{X}}}$, $\Pi_{\boldsymbol{\mathscr{{C}}}}(\boldsymbol{\mathscr{{X}}})$, is an $n_{1}\times n_{2}\times n_{3}$ tensor obtained as

[TABLE]

*where $\ast$ denotes the t-product, $\Pi_{\boldsymbol{\mathscr{{C}}}}(\boldsymbol{\mathscr{{X}}})$ is the projection of $\boldsymbol{\mathscr{{X}}}$ onto the subspace spanned by the lateral slices of $\boldsymbol{\mathscr{{C}}}$, and $\boldsymbol{\mathscr{{C}}}^{\dagger}$ is the Moore-Penrose generalized inverse of tensor $\boldsymbol{\mathscr{{C}}}$ [26]. *

Because the circulant convolution of two tube fibers can be computed by the DFT, the t-project can be alternatively computed in the Fourier domain, as shown in Algorithm 3.

Definition 3.4**.**

(t-CX). Given an $n_{1}\times n_{2}\times n_{3}$ tensor $\boldsymbol{\mathscr{{X}}}$, let $\boldsymbol{\mathscr{{C}}}$ be an $n_{1}\times c\times n_{3}$ tensor whose lateral slices correspond to $c$ lateral slices from tensor $\boldsymbol{\mathscr{{X}}}$. Then, the $n_{1}\times n_{2}\times n_{3}$ tensor

[TABLE]

*is called a t-CX decomposition of $\boldsymbol{\mathscr{{X}}}$. *

The following remarks are of interest for the previous definition.

• •

The choice for the number of lateral slices $c$ in the t-CX approximation depends on the application under consideration; neverthelees, we are primarily interested in the case $c\ll n_{2}$. For example, $c$ could be constant, independent of the dimension of tensor, or logarithmic in the size of $n_{2}$, or a large constant factor less than $n_{2}$.

• •

The t-CX decomposition expresses each $\boldsymbol{\mathscr{{X}}}$ slice in terms of a linear combination of basis slices (see, Definition  2.14), each of which corresponds to an actual lateral slice of $\boldsymbol{\mathscr{{X}}}$. Hence, t-CX provides a low-rank approximation to the original tensor $\boldsymbol{\mathscr{{X}}}$, even though its structural properties are different than those of the t-SVD.

• •

Given a set of lateral slices $\boldsymbol{\mathscr{{C}}}$, the approximation $\Pi_{\boldsymbol{\mathscr{{C}}}}(\boldsymbol{\mathscr{{X}}})=\boldsymbol{\mathscr{{C}}}*\boldsymbol{\mathscr{{C}}}^{\dagger}*\boldsymbol{\mathscr{{X}}}$ is the ”best” approximation to $\boldsymbol{\mathscr{{X}}}$ in the following sense

[TABLE]

Next, we provide the t-CX decomposition algorithm and provide its relative error. Algorithm 4 takes as input an $n_{1}\times n_{2}\times n_{3}$ tensor $\boldsymbol{\mathscr{{A}}}$, a tubal rank parameter $r$, and an error parameter $\epsilon$. It returns as output an $n_{1}\times c\times n_{3}$ tensor $\boldsymbol{\mathscr{{C}}}$ comprising of a small number of slices of $\boldsymbol{\mathscr{{A}}}$. Let $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$, be a tubal rank-$r$ approximation of tensor $\boldsymbol{\mathscr{{A}}}$. Central to our proposed randomized algorithm is the concept of sampling slices of the tensor based on a leverage score [15], defined by

[TABLE]

where $\hat{\boldsymbol{\mathscr{{V}}}}=\texttt{fft}(\boldsymbol{\mathscr{{V}}},[],3)$. The goal is to select a small number of lateral slices of $\boldsymbol{\mathscr{{A}}}$.

Note that we define rescaling and sampled tensors $\hat{\boldsymbol{\mathscr{{D}}}}$ and $\hat{\boldsymbol{\mathscr{{S}}}}$ in the Fourier domain. This leads to a significant reduction in time complexity of the fast Fourier transform and its inverse. For the case $n_{3}=1$, the algorithm selects column-row pairs in Algorithm 4 and is identical to the algorithm of [15].

Next, we present a lemma of general interest. The derived properties aid in obtaining tensor based randomized $\ell_{2}$ regression and low rank estimation guarantees.

Lemma 3.5**.**

Given a target tensor $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, let $\epsilon\in(0,1]$ and $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$ be a $\operatorname{rank}_{t}$-r approximation of $\boldsymbol{\mathscr{{A}}}$. Let $\Gamma=(\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}})^{\dagger}-(\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}})^{\top}$. If the sampling probabilities $\left\{p_{i}\right\}_{i=1}^{n_{2}}$ satisfy (19), and if $c\geq 48r\log(4r/(\beta\delta))/(\beta\epsilon^{2})$, then with probability at least $1-\delta$, the following hold

[TABLE]

*where $\Sigma_{\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}}$ is an F-diagonal tensor and contains the $r$ non-zero singular tubes of $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}$. *

Note that the equation (20d) shows that in terms of its singular tubes, the tensor $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}$ -i.e., the slice sampled and rescaled version of $\boldsymbol{\mathscr{{V}}}$- is almost an orthogonal tensor. A useful property of an orthogonal tensor $\boldsymbol{\mathscr{{V}}}$ is that $\boldsymbol{\mathscr{{V}}}^{\dagger}=\boldsymbol{\mathscr{{V}}}^{\top}$ (see, Definition 2.4). Equation (20b) shows that although this property does not hold for $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}$, the difference between $(\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}})^{\dagger}$ and $(\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}})^{\top}$ can be bounded.

Using Theorem 3.1 and Lemma 3.5, we provide in proposition 3.6 a sampling complexity for a tensor based randomized “$\ell_{2}$-regression” as in [15]. Indeed, Lemmas 2 and 3 of [15] follow by using the Frobenius norm bound of Theorem 3.1 and applying Markov’s inequality. The claim of Lemma 1 of [15] follows by applying Lemma 3.5. If we consider the failing probability $\delta=0.05$, the claims of all three lemmas hold simultaneously with probability at least $0.85$ and in the following proposition we condition on this event. The proof follows along similar lines to the proof of Theorem 5 in [15].

Proposition 3.6**.**

Given a target tensor $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, let $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$ be a $\operatorname{rank}_{t}$-r approximation of $\boldsymbol{\mathscr{{A}}}$. Choose $c=O(r\log(r/\beta)/(\beta\epsilon^{2}))$ slices of $\boldsymbol{\mathscr{{A}}}$ with their sampling probabilities $\left\{p_{i}\right\}_{i=1}^{n_{2}}$ satisfying (19). Then, with probability at least $0.85$, the following holds

[TABLE]

Let $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$ be a t-SVD of tensor $\boldsymbol{\mathscr{{L}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$. In the remainder of the paper, we use the following two parameters related to the coherence of tensor $\boldsymbol{\mathscr{{L}}}$,

[TABLE]

where $\boldsymbol{\mathscr{{U}}}_{c}$ is the left singular tensor of $\boldsymbol{\mathscr{{C}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$ lateral slices of $\boldsymbol{\mathscr{{L}}}$.

The following theorem shows that projection based on slice sampling leads to near optimal estimation in tensor regression, when the covariate tensor has small coherence, $\mu_{0}(\boldsymbol{\mathscr{{V}}})$.

Theorem 3.7**.**

Given $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, let $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$ be a $\operatorname{rank}_{t}$-r approximation of $\boldsymbol{\mathscr{{A}}}$. Choose $c=O(\varrho\log(\varrho)/\epsilon^{2})$, and let $\boldsymbol{\mathscr{{A}}}_{c}\in\mathbb{R}^{n_{1}\times c\times n_{2}}$ be a tensor of $c$ lateral slices of $\boldsymbol{\mathscr{{A}}}$ sampled uniformly without replacement. Further, let $\boldsymbol{\mathscr{{L}}}_{c}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$ consist of the corresponding slices of $\boldsymbol{\mathscr{{L}}}$. Then, with probability at least $0.85$, the following holds

[TABLE]

Next, by using Lemma 3.5, we provide a low rank estimation bound for the t-CX algorithm, under the coherence assumption. Indeed, we will take advantage of a constant $\beta$ in equation (19) to provide relative error estimation guarantees for the t-CX algorithm under uniform slice sampling, when tensor $\boldsymbol{\mathscr{{A}}}$ under consideration exhibits sufficient incoherence.

Corollary 3.8**.**

Given a target tensor $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, let $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$ be a $\operatorname{rank}_{t}-r$ approximation of $\boldsymbol{\mathscr{{A}}}$. Choose $c=O(\varrho\log(\varrho)\log(1/\delta)/\epsilon^{2})$, and let $\boldsymbol{\mathscr{{C}}}$ be a tensor of $c$ lateral slices of $\boldsymbol{\mathscr{{A}}}$ sampled uniformly without replacement. Then, the following holds

[TABLE]

*with probability at least $1-\delta$. *

3.4 Slice-based tensor CUR decomposition

Similar to the matrix case, the t-CX decomposition suffices when $n_{1}\ll n_{2}$ because $\boldsymbol{\mathscr{{C}}}^{\dagger}*\boldsymbol{\mathscr{{A}}}$ is small in size. However, when $n_{1}$ and $n_{2}$ are almost equal, computing and storing the dense matrix $\boldsymbol{\mathscr{{C}}}^{\dagger}*\boldsymbol{\mathscr{{A}}}$ in memory becomes prohibitive. The CUR decomposition provides a very useful alternative. Next, we introduce a tensor CUR (t-CUR) decomposition based on the rt-product.

Definition 3.9**.**

(t-CUR). Let $\boldsymbol{\mathscr{{X}}}$ be an $n_{1}\times n_{2}\times n_{3}$ tensor. For any given $\boldsymbol{\mathscr{{C}}}$, an $n_{1}\times c\times n_{3}$ tensor whose slices comprise of $c$ lateral slices of tensor $\boldsymbol{\mathscr{{X}}}$, and $\boldsymbol{\mathscr{{R}}}$, an $l\times n_{2}\times n_{3}$ tensor whose slices comprise of $l$ horizontal slices of tensor $\boldsymbol{\mathscr{{X}}}$, the $n_{1}\times n_{2}\times n_{3}$ tensor

[TABLE]

*is a lateral-horizontal-based tensor approximation to $\boldsymbol{\mathscr{{X}}}$ for any $c\times l\times n_{3}$ tensor $\boldsymbol{\mathscr{{U}}}$. *

The following remarks are of interest regarding the previous definition.

• •

The t-CUR decomposition is most appropriate as a data analysis tool, when the data consist of one- and/or two modes that are qualitatively different than the remaining ones. In this case, the t-CUR decomposition approximately expresses the original data tensor in terms of a basis consisting of underlying subtensors that are actual data slices and not artificial bases.

• •

The t-CUR approximation is a t-CX approximation, but one with a very special structure; i.e., every lateral slice of $\boldsymbol{\mathscr{{X}}}$ can be expressed in terms of the basis provided by $\boldsymbol{\mathscr{{C}}}$ (see, Definition 2.14) using only the information contained in a small number of horizontal slices of $\boldsymbol{\mathscr{{X}}}$ and a low-dimensional encoding tensor.

• •

In terms of its t-SVD structure, $\boldsymbol{\mathscr{{U}}}$ contains the ”inverse-of-$\boldsymbol{\mathscr{{X}}}$” information. For the proposed t-CUR decomposition, $\boldsymbol{\mathscr{{U}}}$ will be a generalized inverse of the intersection between selected tensors $\boldsymbol{\mathscr{{C}}}$ and $\boldsymbol{\mathscr{{R}}}$.

Note that the structural simplicity of the t-CUR tensor decomposition becomes apparent in the Fourier domain, as detailed in Algorithm 5. The latter takes as input an $n_{1}\times n_{2}\times n_{3}$ tensor $\boldsymbol{\mathscr{{A}}}$, an $n_{1}\times c\times n_{3}$ tensor $\boldsymbol{\mathscr{{C}}}$ consisting of a small number of lateral slices of $\boldsymbol{\mathscr{{A}}}$, and an error parameter $\epsilon$. Letting $\hat{\boldsymbol{\mathscr{{C}}}}=\hat{\boldsymbol{\mathscr{{U}}}}*\hat{\Sigma}*\hat{\boldsymbol{\mathscr{{V}}}}^{\top},$ Algorithm 5 uses sampling probabilities

[TABLE]

to select a small number of horizontal slices of $\boldsymbol{\mathscr{{A}}}$. It returns an $l\times n_{2}\times n_{3}$ tensor $\boldsymbol{\mathscr{{R}}}$ comprising of a small number of horizontal slices of $\boldsymbol{\mathscr{{A}}}$ and an $c\times l\times n_{3}$ tensor $\boldsymbol{\mathscr{{U}}}$ consisting of the corresponding slices of $\boldsymbol{\mathscr{{C}}}$.

Corollary 3.10**.**

Given $\boldsymbol{\mathscr{{A}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, let $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}$ be a $\operatorname{rank}_{t}-r$ approximation of $\boldsymbol{\mathscr{{A}}}$. In Algorithm 5, choose $c=O(\varrho\log(\varrho)\log(1/\delta)/\epsilon^{2})$, and let $\boldsymbol{\mathscr{{C}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$ be a tensor of $c$ slices of $\boldsymbol{\mathscr{{A}}}$ sampled uniformly without replacement. Further, choose $l=O(\varrho_{c}\log(\varrho_{c})\log(1/\delta)/\epsilon^{2})$, and let $\boldsymbol{\mathscr{{R}}}\in\mathbb{R}^{l\times n_{2}\times n_{3}}$ be a tensor of $l$ horizontal slices of $\boldsymbol{\mathscr{{A}}}$ sampled uniformly without replacement. Then, the following holds

[TABLE]

*with probability at least $1-\delta$. *

Remark 3.11**.**

*In many applications such as tensor completion problems discussed in Section 5.3, it may not feasible to compute the t-SVD of the entire tensor $\boldsymbol{\mathscr{{A}}}$ due to either computational cost or large set of missing values. In these cases, algorithms requiring knowledge of the leverage scores can not be applied. Hence, we may use an estimate where a subset of the lateral slices are chosen uniformly and without replacement, and the horizontal leverage scores (24) are calculated using the top-$r$ right singular slices of this lateral tensor instead of the entire $\boldsymbol{\mathscr{{A}}}$ tensor. *

Remark 3.12**.**

*Note that for any subspace, the smallest $\mu_{0}$ can be is 1. In such case, we are using the optimal subspace sampling with $\beta=1$. Note also that we have provided Corollaries (3.8) and (3.10) based on uniform sampling. However, it can be easily seen that the relative error guarantees hold with nonuniform sampling probabilities (19) and (24), if we set $c=O(r\log(r/\beta)\log(1/\delta)/(\beta\epsilon^{2}))$, and $l=O(c\log(c/\beta)\log(1/\delta)/(\beta\epsilon^{2}))$. *

4 Tensor completion and robust factorization

Next, we provide a CUR tensor nuclear norm minimization (CUR t-NN) procedure that solves a noisy tensor factorization problem, using a small number of lateral and/or horizontal slices of the underlying tensor and exhibits favorable computational complexity and in addition comes with performance guarantees.

4.1 Related work on matrix problems

Low rank plus sparse matrix decomposition. In many image processing and computer vision applications, the given data matrix $X$ can be decomposed as a sum of a low-rank and a sparse component. To that end, Candés et al. proposed Robust PCA [6] to model data matrices generated according to the following mechanism:

[TABLE]

Note that the solution of (25) is an NP-hard problem. It is established in [7] that if $L$ exhibits a certain degree of low-rankness, while $E$ is sparse enough, then the formulation of (25) can be relaxed into a convex problem of the form:

[TABLE]

This model implicitly assumes that the underlying data structure lies in a single low-rank subspace. However, in many applications (e.g. image classification) it is more likely that the data are obtained from a union of multiple subspaces, and hence recovery of the structure based on the above decomposition would be inaccurate. In order to segment the data into their respective subspaces, one needs to compute an affinity matrix that encodes the pairwise affinities between data vectors. Liu [32] proposed a more general rank minimization problem, where the data matrix itself is used as the dictionary, resulting in the following convex optimization problem:

[TABLE]

When the subspaces are globally independent, the data are noiseless and sampling is sufficient, Liu et al. [32] show that the optimal solution, denoted by $L^{*}$, to the problem given by 27 corresponds to the widely used Shape Iteration Matrix (SIM) method [9]. The latter is a ”block-diagonal” affinity matrix that indicates the true segmentation of the data. To handle data corrupted by noise, the popular Low Rank Representation (LRR) introduced in [32] adopts a regularized formulation that introduces an extra penalty term to fit the noise component. Further, after obtaining the self-representation matrix $L$, the affinity matrix $C$ is usually constructed as $C=\frac{1}{2}(|L|+|L^{\top}|),$ where $|\cdot|$ represents the absolute operator. Then, the obtained affinity matrix $C$ will be processed through a spectral clustering algorithm [40] to produce the final clustering result and obtain the corresponding data generating subspaces.

It is established in [52, 47] that combining sparse and low-rank regularization can improve the performance of image classification. The basic objective function of this combination [52] is as follows:

[TABLE]

where $\lambda_{1},\lambda_{2},\lambda_{3}$ are tuning parameters and $\|E\|_{\ell}$ indicates different norms suitable for different types for corrupting the data by noise; for example, the squared Frobenius norm for Gaussian noise and the $\ell_{1}$ norm for random spiked noise. Equation (28) is similar to the objective functions in [47, 16], where a detailed explanation of the formulation given in (28) is also provided.

4.1.1 Low rank tensor decomposition

Lu et al. [35] extended Robust PCA [6] to the third order tensor based on t-SVD and proposed the following convex optimization problem:

[TABLE]

where $\lambda$ is a tuning parameter.

This model implicitly assumes that the underlying data come from a single low-rank subspace. When the data is drawn from a union of multiple subspaces, which is common in image classification, the recovery based on the above formulation may lack in accuracy. To that end, Xie et al. [49] extended LRR based subspace clustering to a multi-view one, by employing the rank sum of different mode unfoldings to constrain the subspace coefficient tensor, resulting in the following convex optimization problem:

[TABLE]

Equation (30) is similar to the objective functions in [25, 42] and a detail explanation of (30) is provided in that paper.

4.2 Proposed algorithm

Next, we propose an algorithm for large scale tensor decomposition of noisy data. Our proposal, called CUR t-NN, extends Algorithm 5 to the noisy tensor factorization. The main steps are outlined in Algorithm 6.

Note that CUR t-NN can be used in combination with an arbitrary optimization algorithm. In this paper, we have chosen to solve the noisy tensor factorization formulations (10), (29) and (30) using an Alternating Direction Method of Multiplier (ADMM) algorithm [3]. ADMM is the most widely used approach for robust tensor PCA in both the fully and partially observed settings. Indeed, ADMM achieves much higher accuracy than (accelerated) proximal gradient algorithm using fewer iterations. It works well across a wide range of problem settings and does not require careful tuning of the regularization parameters. Further, the following empirical finding has been frequently observed: namely, the rank of the iterates often remains bounded by the rank of the initializer, thus enabling efficient computations [6]. This feature is not shared by the block coordinate decent algorithm. We provide a variant of ADMM for solving problem (29) in the Appendix (see, Algorithm 7). With a small modification, Algorithm 7 can be also used to solve the tensor completion problem 10, and tensor subspace clustering problem (30). In the following algorithm, $\textbf{ADMM}(\boldsymbol{\mathscr{{X}}},\lambda)$ denotes the ADMM algorithm for solving regularized tensor nuclear norm minimization problems (10), (29) and (30), where $\boldsymbol{\mathscr{{X}}}$ is the (sampled) data tensor and $\lambda$ is a regularization parameter.

4.3 Running Time of CUR t-NN

Algorithm 6 significantly reduces the per-iteration complexity of nuclear norm minimization problems. Indeed, in each iteration, a base tensor nuclear norm minimization algorithm requires $O(n_{1}n_{2}n_{3}\log(n_{3}))$ flops to transform the tensor to the Fourier domain, $O(n_{1}n_{2}n_{3}\min(n_{1},n_{2}))$ flops for the t-SVD computation and factorization, and $O(n_{1}n_{2}n_{3}\log(n_{3}))$ flops to transform it back to the original domain. On the other hand, Algorithm 6 only requires $O(\max(cn_{1},ln_{2})n_{3}\log(n_{3}))$, $O(\max(n_{1}c,ln_{2})n_{3}r)$ and $O(\max(cn_{1},ln_{2})n_{3}\log(n_{3}))$ flops, for the respective steps. Further, Algorithm 6 can be implemented without storing the data tensor $\boldsymbol{\mathscr{{X}}}$ and can be advantageous when $r\ll\min(n_{1},n_{2})$, which occurs frequently in real data sets.

4.4 Theoretical guarantees

Next, using Lemma 3.5, we establish that Algorithm 6 exhibits high probability recovery guarantees comparable to those of the exact algorithms that use the full data tensor. Our first result bounds the $\mu_{0}$ and $\mu_{1}$-coherence (see, Definitions (2.15) and (2.16)) of a randomized tensor in terms of the coherence of the full tensor.

Lemma 4.1**.**

Let $\boldsymbol{\mathscr{{L}}}_{c}$ be a tensor formed by selecting $c$ slices of a $\operatorname{rank}_{t}-r$ tensor $\boldsymbol{\mathscr{{L}}}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}$ that satisfy the probabilistic conditions in (19). If $c\geq 48r\log(4r/(\beta\delta))/(\beta\epsilon^{2})$, then with probability at least $1-\delta$:

• •

$\mu_{0}(\boldsymbol{\mathscr{{U}}}_{\boldsymbol{\mathscr{{L}}}_{c}})=\mu_{0}(\boldsymbol{\mathscr{{U}}})$,

• •

$\displaystyle\mu_{0}(\boldsymbol{\mathscr{{V}}}_{\boldsymbol{\mathscr{{L}}}_{c}})\leq\frac{1}{1-\epsilon/2}\mu_{0}(\boldsymbol{\mathscr{{V}}})$,

• •

$\displaystyle\mu_{1}(\boldsymbol{\mathscr{{L}}}_{c})\leq\frac{r}{1-\epsilon/2}\mu_{0}(\boldsymbol{\mathscr{{U}}})\mu_{0}(\boldsymbol{\mathscr{{V}}})$,

*where $\epsilon\in(0,1]$. *

Our next theorem provides a bound for the estimation error of the CUR t-NN Algorithm.

Theorem 4.2**.**

Under the notion of Algorithm 6, choose $c=O(\varrho\log(\varrho)\log(1/\delta)/\epsilon^{2})$, and $l=O(\varrho_{c}\log(\varrho_{c})\log(1/\delta)/\epsilon^{2})$. Let $\boldsymbol{\mathscr{{C}}}^{*}$ and $\boldsymbol{\mathscr{{R}}}^{*}$ be the corresponding lateral and horizontal sub-tensors of the exact solution $\boldsymbol{\mathscr{{L}}}^{*}$. If $\boldsymbol{\mathscr{{L}}}^{*}$ is $(\mu,r)$-coherent, then with probability at least $1-\delta$, $\boldsymbol{\mathscr{{C}}}^{*}$ and $\boldsymbol{\mathscr{{R}}}^{*}$ are $(\frac{r\mu^{2}}{1-\epsilon/2},r)$-coherent, and

[TABLE]

*where $\tilde{\boldsymbol{\mathscr{{L}}}}$ is a solution obtained by Algorithm 6. *

5 Experimental Results

Next, we investigate the efficiency of the proposed randomized algorithms on both synthetic and real data sets. The results are organized in the following three sub-sections: in Section 5.1, we compare the performance of the rt-product to that of the t-product using synthetic data. In Section 5.2, we use the proposed t-CX and t-CUR algorithms for finding important ions and positions in two Mass Spectrometry Imaging (MSI) data sets. Finally, in Section 5.3, we apply the proposed tensor factorization CUR t-NN algorithm on data sets related to image and video processing.

All algorithms have been implemented in the MATLAB R2018a environment and run on a Mac machine equipped with a 1.8 GHz Intel Core i5 processor and 8 GB 1600 MHz DDR3 of memory.

5.1 Experimental results for fast tensor multipication

A random data tensor $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ of dimension $n_{1}=10^{5}$, $n_{2}=10^{3}$ and $n_{3}=5$ is generated as follows: the first frontal slice of $\boldsymbol{\mathscr{{X}}}$ is set to $\boldsymbol{\mathscr{{X}}}(:,:,1)=\texttt{sprandn}(n_{1},n_{2},1)$; i.e. comprising of $n_{1}\times n_{2}$ sparse normally distributed random matrix. The remaining frontal slices are generated by $\texttt{sprandn}(\boldsymbol{\mathscr{{X}}}(:,:,1))$, which results in the same sparsity structure as $\boldsymbol{\mathscr{{X}}}(:,:,1)$, with normally distributed random entries with mean 0 and variance 1.

We perform a t-SVD decomposition on $\boldsymbol{\mathscr{{X}}}$ and select the left singular value tensor $\boldsymbol{\mathscr{{U}}}$ for our experiments. We consider both uniform and nonuniform sampling schemes. For nonuniform sampling, we use a randomized approach (similar to Algorithm 4 in [14]) to obtain normalized leverage scores of $\boldsymbol{\mathscr{{U}}}$. Further, the Frobenius and spectral bounds given in (15) and (16), respectively, are used as performance metrics for Algorithm 2. Table (2) shows the Relative Frobenius Error (RFE) -$\|\boldsymbol{\mathscr{{U}}}_{r}^{\top}\boldsymbol{\mathscr{{U}}}_{r}-\tilde{\boldsymbol{\mathscr{{U}}}}_{r}^{\top}\tilde{\boldsymbol{\mathscr{{U}}}}_{r}\|_{F}/\|\boldsymbol{\mathscr{{U}}}_{r}\|_{F}^{2}$- and Relative Spectral Error (RSE) -$\|\boldsymbol{\mathscr{{U}}}_{r}^{\top}\boldsymbol{\mathscr{{U}}}_{r}-\tilde{\boldsymbol{\mathscr{{U}}}}_{r}^{\top}\tilde{\boldsymbol{\mathscr{{U}}}}_{r}\|_{2}/\|\boldsymbol{\mathscr{{U}}}_{r}\|_{2}^{2}$- for the rt-product algorithm to recover $\tilde{\boldsymbol{\mathscr{{U}}}}_{r}^{\top}\tilde{\boldsymbol{\mathscr{{U}}}}_{r}=\boldsymbol{\mathscr{{U}}}_{r}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}*\boldsymbol{\mathscr{{D}}}*\boldsymbol{\mathscr{{S}}}^{\top}*\boldsymbol{\mathscr{{U}}}_{r}$ based on $l\approx r\log(r)$ slices, and its deterministic counterpart. The depicted results are averaged over 10 independent replications.

Very large improvements in computational speed can be seen when using the rt-product, especially when coupled with a non-uniform sampling scheme. Further, the gains become more pronounced for larger number of slices sampled and larger rank.

5.2 Finding Important Ions and Positions in MSI

MSI is used to visualize the spatial distribution of chemical compounds, such as metabolites, peptides or proteins by their molecular masses. The ability of MSI to localize panels of biomolecules in tissue samples has led to a rapid and substantial impact in both clinical and pharmacological research, aided in uncovering biomolecular changes associated with disease and finally provided low cost imaging of drugs. Typical techniques used require finding important ions and positions from a 3D image: ions to be used in fragmentation studies to identify key compounds, and positions for follow up validation measurements using microdissection or other orthogonal techniques. Unfortunately, with modern imaging machines, these must be selected from an overwhelming amount of raw data. Existing popular techniques used to reduce the volume of data, include principal components analysis and non-negative matrix factorization, but return difficult-to-interpret linear combinations of actual data elements. A recent paper [50] shows that CX and CUR matrix decompositions can be used directly to address this selection need. One major shortcoming of CX and/or CUR matrix decompositions is that they can only handle 2-way (matrix) data. However, MSI data form a multi-dimensional array. Hence, in order to use CX/CUR matrix decompositions, one has first to reformulate the multi-way array as a matrix. Such preprocessing usually leads to information loss, which in turn could cause significant performance degradation.

By using instead the t-CX and t-CUR decompositions (Algorithms 4 and 5, respectively) one can obtain good low-rank approximations of the available data, expressed as combinations of actual ions and positions, as opposed to difficult-to-interpret eigen-ions and eigen-positions produced by matrix factorization techniques. We show that this leads to effective prioritization of information for both actual ions and actual positions. In particular, important ions can be discovered by using leverage scores as the importance sampling distribution. Further, selection of important positions from the original tensor can be accomplished based on the random sampling algorithm in [50], since the distribution of the leverage scores of positions is uniform. To this end, we consider the following two ways of computing leverage scores of a given data set:

• •

Deterministic: Compute the normalized tensor leverage scores exactly using probabilities (19) and (24).

• •

Randomized: Compute an approximation to the normalized leverage scores of tensor (mapped to a block diagonal matrix in the Fourier domain) by using Algorithm 4 of [14].

5.2.1 Description and Analysis of MSI data sets

Next, we use the following two data sets for illustration purposes, that are publicly available at the OpenMSI Web gateway222https://openmsi.nersc.gov/openmsi/client/. They represent two diverse acquisition modalities, including one mass spectrometry image of the left coronal hemisphere of a mouse brain (see, Figure 3) acquired using a time-of-flight (TOF) mass analyzer and one MSI data set of a lung acquired using an Orbitrap mass analyzer. These data sets form a $122\times 120\times 80339$ and a $122\times 120\times 500000$ tensor, respectively.

As described in [50], the brain data set is processed using peak-finding to identify the most intense ions. Using this technique, the original brain data is reduced from $122\times 120\times 80339$ values to a data set of size $122\times 120\times 2000$. To compute the CX decomposition, we reshape the MSI data cube into a two-dimensional $(14640\times 2000)$ matrix, where each row corresponds to the spectrum of a pixel in the image, and each column to the intensity of an ion across all pixels, thus describing the distribution of the ion in physical space. No peak-finding was applied to the lung data set.

For the brain data set, we evaluate the quality of the approximation of the leverage scores based on a rank $r=5$ approximation. The distribution of the deterministic leverage scores for the ions and pixels is shown in Figure 4. Table (2) shows the relative square error (RSE) $\|\boldsymbol{\mathscr{{X}}}-\tilde{\boldsymbol{\mathscr{{X}}}}\|_{F}/\|\boldsymbol{\mathscr{{X}}}\|_{F}$ using t-CX decomposition to recover rank $-5$ tensor $\tilde{\boldsymbol{\mathscr{{X}}}}$ for selection of $c=25,35,45$ and, $55$ ions, using both randomized and deterministic leverage scores. Table (3) provides the reconstruction errors to the best rank$-5$ approximation, based on the t-CX decomposition coupled with horizontal slice selection. The results obtained running both randomized and deterministic CX and t-CX algorithms are based on 10 independent replicates and then averaging them. Note that for pixel selection, the deterministic CX decomposition results in larger reconstruction errors than its randomized CX counterpart. The reason for this behavior lies in the distribution of the leverage scores for the pixels, which are fairly uniform (see, Figure 4).

For the lung data set, reconstruction errors to the best rank $r=15$ approximation based on randomized t-CX and t-CUR decompositions are given in Table (4) based on averages over 10 independent replicates. It can be seen that the t-CX and t-CUR match or outperform their matrix variants in terms of accuracy, while improving on computing time.

These results introduce the concept of t-CX/ t-CUR tensor factorizations to MSI, describing their utility and illustrating principled algorithmic approaches to deal with the overwhelming amount of data generated by this technology and their ability to select important and intepretable ions/pixels.

5.3 Robust PCA in the fully and partially observed settings

Many images exhibit an inherent low rank structure and are suitably denoised by low-rank modeling methods, such as robust PCA [6]. In this section, we assess the performance of the CUR t-NN on two popular data sets, and for the typical use cases they represent. We compare the performance of the proposed CUR t-NN algorithm to the following techniques:

• •

EXACT NN, the exact matrix completion [7];

• •

RPCA NN, the robust matrix completion [6];

• •

E-TUCKER NN, the TUCKER based tensor completion [33];

• •

R-TUCKER NN, the robust TUCKER based tensor completion [20];

• •

EXACT t-NN, the t-SVD based tensor completion [53];

• •

RPCA t-NN, the robust t-SVD based tensor completion [35];

All algorithms are terminated either when the relative square error (RSE),

[TABLE]

or the number of iterations and CPU times exceed 1,000 and 20 minutes, respectively.

5.3.1 CUR t-NN on the Extended Yale B data set

We apply the CUR t-NN on the Extended Yale B data set to evaluate the accuracy of the proposed low-rank representations, as well as its computation time. The database consists of 2432 images of 38 individuals, each under 64 different lighting conditions [17]. We used 30 images from each subject and kept them at full resolution. Each image comprises of $192\times 168$ pixels on a grayscale range. The data are organized into a $192\times 1140\times 168$ tensor that exhibits low tubal rank, which is an expected feature due to the spatial correlation within lateral slices [18]. Laplacian (salt and pepper 333This noise can be caused by sharp and sudden disturbances in the image signal. It demonstrates itself as sparsely occurring white and black pixels.) noise was introduced separately in all frontal slices of the observation tensor for 20% of the entries.

We provide visualizations of the reconstructed first image from the first subject at the 20% noise level in Figure 5. We compare the performance of the CUR t-NN, RPCA t-NN, R-TUCKER NN and RPCA NN for solving the robust tensor low rank approximation problem (29). Table 5 shows the accuracy of the CUR t-NN together with that of competing algorithms for the Extended Yale B data set. As shown in Figure 5 and Table 5, CUR t-NN estimates nearly the same face model as the RPCA t-NN requiring only a small fraction of time. On the other hand, all algorithms exhibit small RSE, but the CUR t-NN proves essentially as competitive as the best method RPCA t-NN, but achieves almost the same performance at approximately $1/4$ of the time.

Beyond just speed-ups and/or accuracy improvements of the CUR t- NN algorithm, its output can be directly used in place of the singular slices and tubes that standard methods provide. The latter represent linear combinations of the slices of the tensor, which for an image data set capture an “average eigenface”. On the other hand, CUR t-NN reconstructs the tensor through selection of real faces in the data set, thus giving the opportunity to researchers to inspect them and examine their representativeness. Hence, similar to the original CUR decomposition of matrices, CUR t-NN enhances the interpretability of the tensor decomposition [18, 48, 41].

5.3.2 CUR t-NN on a video data set

Next, we compare the CUR t-NN to the aforementioned listed competing methods for video data representation and compression from randomly missing entries. The video data, henceforth referred to as the Basketball video (source: YouTube) is mapped to a $144\times 256\times 80$ tensor, obtained from with a nonstationary panning camera moving from left to right horizontally following the running players. We randomly sampled 50% entries from the Basketball video. We compare the performance of the CUR t-NN, EXACT t-NN, E-TUCKER NN and EXACT NN for solving the tensor completion problem (10).

The result is shown in Table 6. It can be seen that the CUR t-NN outperforms almost all its competitors in terms of CPU running time and accuracy and essentially matches that of EXACT t-NN.

6 Conclusion

This paper introduced two randomized algorithms for basic tensor operations -rt-product and rt-project- and then used in tensor CX and CUR decomposition algorithms, whose aim is to select informative slices. The randomized tensor operations together with the tensor decompositions algorithms comes with relative error recovery guarantees. These algorithms can be effectively used in the analysis of large scale tensors with small tubal rank.

In addition, we proposed the CUR t-NN algorithm that exploits the advantages of randomization for dimensionality reduction and can be used effectively for large scale noisy tensor decompositions. Indeed, CUR t-NN uses an adaptive technique to sample slices of the tensor based on a leverage score for them and subsequently solves a convex optimization problem followed by a projection step to recover a low rank approximation to the underlying tensor. The proposed algorithm has linear running time, and provably maintains the recovery guarantees of the exact algorithm with full data tensor.

Acknowledgments

The authors would like to thank the Associate Editor and three anonymous referees for many constructive comments and suggestions that improved significantly the structure and readability of the paper.

6.1 Proof of Theorem 3.1

Let $F_{n_{3}}$ denote the Discrete Fourier Matrix. Then, Definition 2.1 implies that

[TABLE]

Using (31) and the unitary invariance of the Frobenius norm, we have

[TABLE]

In Algorithm 2, we define $I_{j}$, $j\in[n_{2}]$ as an indicator variable, which is set to $1$ if both $thej$-th lateral slice of $\boldsymbol{\mathscr{{A}}}$ and the $j$-th horizontal slice of $\boldsymbol{\mathscr{{B}}}$ are selected. In this case, the selected lateral and horizontal slices are scaled by score $1/\sqrt{\min\{1,cp_{j}\}}$. Note that if $\min\{1,cp_{j}\}=1$, then $I_{j}=1$ with probability 1, and $1-I_{j}/\min\{1,cp_{j}\}=0$. Then, taking expectation on both sides of (32) and considering the set $\Upsilon=\{j\in[n_{2}]:cp_{j}<1\}\subseteq[n_{2}]$, we get

[TABLE]

where $\hat{p}_{i_{1}i_{2}j_{1}j_{2}}=\left(1-\frac{I_{j_{1}}}{cp_{j_{1}}}\right)\left(1-\frac{I_{j_{2}}}{cp_{j_{2}}}\right)\hat{\boldsymbol{\mathscr{{A}}}}_{i_{1}j_{1}k}\hat{\boldsymbol{\mathscr{{B}}}}_{j_{1}i_{2}k}\hat{\boldsymbol{\mathscr{{A}}}}_{i_{1}j_{2}k}\hat{\boldsymbol{\mathscr{{B}}}}_{j_{2}i_{2}k}$.

Since for $j\in[\Upsilon]$, $\mathbf{E}[1-I_{j}/cp_{j}]=0$ and $\mathbf{E}[\left(1-I_{j}/cp_{j}\right)^{2}]=\left(1/cp_{j}\right)-1\leq 1/cp_{j}$, we get

[TABLE]

Equation (15) follows from (34) by using Jensen’s inequality and the fact that the sampling probabilities (13) and (14) are defined in the original domain.

6.2 Proof of Lemma 3.2

Using Definitions 2.1 and 2.9, we have that

[TABLE]

where the first equality follows from the unitary invariance of the spectral norm, the second equality from (31) and the third equality follows since the spectral norm of a block matrix is equal to the maximum of block norms.

Let $k_{\pi}=k\in[n_{3}]$ be the index of the tensor’s frontal slice with maximum spectral norm $\|\hat{\boldsymbol{\mathscr{{A}}}}_{::k}\hat{\boldsymbol{\mathscr{{A}}}}_{::k}^{\top}-\hat{\boldsymbol{\mathscr{{C}}}}_{::k}\hat{\boldsymbol{\mathscr{{C}}}}_{::k}^{\top}\|_{2}$. Using (35), and Example 4.3 in [19], we obtain

[TABLE]

with probability at least $1-\delta$.

6.3 Proof of Lemma 3.5

Let $\sigma_{i,k}$ be the $i$-th largest singular value of slice $k$. For all $1\leq i\leq r$ and $1\leq k\leq n_{3}$, we have

[TABLE]

Using Lemma 3.2 and (37), we get

[TABLE]

for all $1\leq i\leq r$ and an index $k_{\pi}\in[n_{3}]$.

Since $\epsilon\in(0,1]$, each tubal singular value of $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}$ is positive, which implies that $\operatorname{rank}_{t}({\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}})=\operatorname{rank}_{t}({\boldsymbol{\mathscr{{V}}}})=\operatorname{rank}_{t}({\boldsymbol{\mathscr{{L}}}})$.

To prove (20b), we use the t-SVD of $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}$ and note that

[TABLE]

The claim follows since $\boldsymbol{\mathscr{{V}}}_{\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}}$ and $\boldsymbol{\mathscr{{U}}}_{\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{D}}}}$ are orthogonal tensors.

To prove (20c), note that

[TABLE]

To remove the pseudoinverse in the above derivations, we use the first part of the Lemma. In this case,

[TABLE]

By combining (40) and (41), we get the result.

To prove (20d), we have that for all $1\leq i\leq r$ and $1\leq k\leq n_{3}$,

[TABLE]

6.4 Proof of Theorem 3.7

From Definition 2.15, and using (19) with $\beta=n_{3}/\mu_{0}(\boldsymbol{\mathscr{{V}}})\in(0,1]$, we have that

[TABLE]

for all $i\in[n_{2}]$.

Using Proposition 3.6, our choice of $c$ implies that

[TABLE]

holds with probability at least $0.85$.

6.5 Proof of Corollary 3.8

Since $\boldsymbol{\mathscr{{C}}}^{+}*\boldsymbol{\mathscr{{A}}}$ minimizes $\|\boldsymbol{\mathscr{{A}}}-\boldsymbol{\mathscr{{C}}}*\boldsymbol{\mathscr{{X}}}\|_{F}$ over all tensors $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$, it follows that

[TABLE]

Now, using (42), we get

[TABLE]

with probability at least $0.85$.

We can trivially boost the success probability to $1-\delta$ by repeating Algorithm 4 $O(\log(1/\delta))$ rounds. Specifically, let $\boldsymbol{\mathscr{{C}}}_{i}$ denote the output of Algorithm 4 at round $i$; using (43) for each $\boldsymbol{\mathscr{{C}}}_{i}$ we have

[TABLE]

with probability at least $0.85$.

Now, let $\boldsymbol{\mathscr{{C}}}$ denote the set of all columns used in the approximation. Since each $\boldsymbol{\mathscr{{C}}}_{i}=\boldsymbol{\mathscr{{C}}}*\boldsymbol{\mathscr{{S}}}_{i}$ for some tensor $\boldsymbol{\mathscr{{S}}}_{i}$ and $\boldsymbol{\mathscr{{C}}}^{+}*\boldsymbol{\mathscr{{A}}}$ minimizes $\|\boldsymbol{\mathscr{{A}}}-\boldsymbol{\mathscr{{C}}}*\boldsymbol{\mathscr{{X}}}\|_{F}$ over all tensors $\boldsymbol{\mathscr{{X}}}\in\mathbb{R}^{n_{1}\times c\times n_{3}}$, it follows that

[TABLE]

for each $i$. Hence, if

[TABLE]

fails to hold, then for each $i$, (44) also fails to hold. Since $0.15<1/e$, the desired conclusion must hold with probability at least $1-(1/e)^{\log(1/\delta)}=1-\delta$.

6.6 Proof of Corollary 3.10

Using (42), it follows that

[TABLE]

with probability at least $0.85$.

Using (43) and (45), our choice of $c$ and $l$ implies the following holds with probabilities at least $0.7$,

[TABLE]

where $\epsilon^{\prime}=3\epsilon$.

The inequality (46) holds with probability at least $1-\delta$ by following the boosting strategy employed in the proof of Corollary 3.8. Indeed, since in each trial inequality (46) fails with probability less than $0.3<1/e$, the claim of Corollary 3.10 will hold with probability greater than $1-(1/e)^{\log(1/\delta)}=1-\delta$.

6.7 Proof of Lemma 4.1

Since from Lemma 3.5 we have $\operatorname{rank}_{t}(\boldsymbol{\mathscr{{L}}}_{c})=\operatorname{rank}_{t}(\boldsymbol{\mathscr{{L}}})$, the first claim follows similarly by using Lemma 1 of [38].

To prove the second claim, using Lemma 3.5, assume that $\boldsymbol{\mathscr{{S}}}^{\top}*\boldsymbol{\mathscr{{V}}}$ consists of the first $c$ horizontal slices of $\boldsymbol{\mathscr{{V}}}$. Then, if $\boldsymbol{\mathscr{{L}}}_{c}=\boldsymbol{\mathscr{{U}}}*\Sigma*\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}$ has $\operatorname{rank}_{t}(\boldsymbol{\mathscr{{L}}}_{c})=\operatorname{rank}_{t}(\boldsymbol{\mathscr{{L}}})=r$, the tensor $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}$ must have full tubal rank. Thus, we can write

[TABLE]

where the second and third equalities follow from $\boldsymbol{\mathscr{{U}}}^{\top}*\boldsymbol{\mathscr{{U}}}=\boldsymbol{\mathscr{{I}}}_{c}$. The fourth and fifth equalities result from $\Sigma$ having full tubal rank and $\boldsymbol{\mathscr{{V}}}$ having full lateral slice rank, and the sixth follows from $\boldsymbol{\mathscr{{S}}}^{\top}*\boldsymbol{\mathscr{{V}}}$ having full horizontal slice rank. Next, denote the right singular vectors of $\boldsymbol{\mathscr{{L}}}_{c}$ by $\boldsymbol{\mathscr{{V}}}_{c}\in\mathbb{R}^{c\times r\times n_{3}}$. Define $\boldsymbol{\mathring{\mathfrak{{e}}}}_{i,c}$ as the $i$-th lateral slice of $\boldsymbol{\mathscr{{I}}}_{c}$ and $\boldsymbol{\mathring{\mathfrak{{e}}}}_{i,n_{2}}$ as the $i$-th lateral slice of $\boldsymbol{\mathscr{{I}}}$. Then we have,

[TABLE]

where $\boldsymbol{\mathscr{{W}}}=\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathscr{{S}}}^{\top}*\boldsymbol{\mathscr{{V}}}$ and the final equality follows from $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathscr{{S}}}*\boldsymbol{\mathring{\mathfrak{{e}}}}_{i,c}=\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathring{\mathfrak{{e}}}}_{i,n_{2}}$ for all $1\leq i\leq c$.

Next, we have

[TABLE]

where the last inequality follows form Hölder’s inequality.

Since $\boldsymbol{\mathscr{{V}}}^{\top}*\boldsymbol{\mathring{\mathfrak{{e}}}}_{i,n_{2}}*\boldsymbol{\mathring{\mathfrak{{e}}}}_{i,n_{2}}^{\top}*\boldsymbol{\mathscr{{V}}}$ has tubal rank one, using the definition of $\mu_{0}$-coherence, we have

[TABLE]

Now, using (38), we have that $\|\boldsymbol{\mathscr{{W}}}^{-1}\|^{2}\leq\frac{n_{2}}{c(1-\epsilon/2)}$. Thus, it follows that

[TABLE]

To prove the last claim under Lemma 3.5, we note that

[TABLE]

6.8 Proof of Theorem 4.2

Define $A(\boldsymbol{\mathscr{{X}}})$ as the event that a tensor $\boldsymbol{\mathscr{{X}}}$ is $(\frac{r\mu^{2}}{1-\epsilon/2},r)$-coherent. To prove Theorem 4.2, let $\tilde{\boldsymbol{\mathscr{{L}}}}$ denote the solution obtained by CUR t-NN and let $\boldsymbol{\mathscr{{L}}}^{*}$ be the exact solution of problems (10), (29) and (30). In Algorithm 6, define $\bar{\boldsymbol{\mathscr{{L}}}}$ as

[TABLE]

where $\tilde{\boldsymbol{\mathscr{{C}}}}=\begin{bmatrix}\tilde{\boldsymbol{\mathscr{{C}}}}_{1}&\tilde{\boldsymbol{\mathscr{{C}}}}_{2}\end{bmatrix}^{\top}$, and $\tilde{\boldsymbol{\mathscr{{R}}}}=\begin{bmatrix}\tilde{\boldsymbol{\mathscr{{R}}}}_{1}&\tilde{\boldsymbol{\mathscr{{R}}}}_{2}\end{bmatrix}$, and $\boldsymbol{\mathscr{{L}}}_{22}^{*}\in\mathbb{R}^{(n_{1}-l)\times(n_{2}-c)\times n_{3}}$ is the bottom right subtensor of $\boldsymbol{\mathscr{{L}}}^{*}$. It can easily be seen that

[TABLE]

Now, define $W(\tilde{\boldsymbol{\mathscr{{L}}}},\bar{\boldsymbol{\mathscr{{L}}}})$ as the event

[TABLE]

If $W(\tilde{\boldsymbol{\mathscr{{L}}}},\bar{\boldsymbol{\mathscr{{L}}}})$ holds, we have

[TABLE]

Next, we consider all three events $W(\tilde{\boldsymbol{\mathscr{{L}}}},\bar{\boldsymbol{\mathscr{{L}}}})$, $A(\boldsymbol{\mathscr{{R}}}^{*})$, and $A(\boldsymbol{\mathscr{{C}}}^{*})$ with $\log(3/\delta)$. Using Lemma 4.1, our choice of $c$ and $l$ implies that $A(\boldsymbol{\mathscr{{C}}}^{*})$ and $A(\boldsymbol{\mathscr{{R}}}^{*})$ holds with probability at least $1-\delta/3$. Since $\tilde{\boldsymbol{\mathscr{{L}}}}$ is a t-CUR approximation to $\bar{\boldsymbol{\mathscr{{L}}}}$, from Corollary 3.10, we get that $W(\tilde{\boldsymbol{\mathscr{{L}}}},\bar{\boldsymbol{\mathscr{{L}}}})$ holds with probability at least $1-\delta/3$.

Using the union bound, it follows that

[TABLE]

6.9 Optimization by ADMM

We provide the optimization and parameter setting details of ADMM used in Algorithm 7 for solving a robust tensor factorization.

As discussed in [54], the updates of $\boldsymbol{\mathscr{{L}}}_{k+1}$ and $\boldsymbol{\mathscr{{E}}}_{k+1}$ have closed form solutions. It is easy to see that the main per-iteration cost of Algorithm 7 is in the update of $\boldsymbol{\mathscr{{L}}}_{k+1}$, which requires computing the fft of $\boldsymbol{\mathscr{{L}}}$ and the SVD of block matrix $\hat{L}=\texttt{fft}(\boldsymbol{\mathscr{{L}}})$ in the Fourier domain.

The next result establishes the global convergence of the ADMM for solving problem (29) (for details on the convergence analysis, see [46]). Note that similar results hold for solving problems (10), and (30).

Theorem 6.1**.**

*The sequence $(\boldsymbol{\mathscr{{L}}}_{k},\boldsymbol{\mathscr{{E}}}_{k},\boldsymbol{\mathscr{{Y}}}_{k})$ generated by Algorithm 7 from any starting point converges to a stationary point of problem (29). *