\section{Introduction}\label{sec:introduction}}

\IEEEPARstart{L}{ow}-rank tensor completion (LRTC) \cite{gandy2011tensor,acar2011scalable,liu2012tensor,kressner2014low,yuan2016tensor}, as a high-order generalization of low-rank matrix completion (LRMC)\cite{srebro2005rank,CandesRecht2009,wangxu2012,MC_TNNR_PAMI2013,MC_ShattenP_AAAI125165, MC_MultiLlabel2015PAMI,Fan_2019_CVPR,NIPS2019_8754}, aims to recover the missing entries of a low-rank tensor. LRTC has recently drawn the attention of many researchers owing to the following reasons. First, nowadays, with the development of information science and technologies, especially data acquisition techniques, real data in many areas such as computer vision, biomedical engineering, and e-commerce are often organized in tensors (multi-way arrays) of low rank or approximately low rank. Often caused by failure or high expense in data acquisition, the presence of missing data in tensors leads to difficulties for data analysis and decision making. Second, some important properties hold for matrices but not for tensors \cite{hillar2013most}, making tensor completion more challenging than matrix completion both numerically and theoretically. For instance, it is NP-hard to compute the (CANDECOMP/PARAFAC, or CP \cite{kolda2009tensor}) rank of a tensor, while the rank of a matrix can be computed efficiently by singular value decomposition (SVD) \cite{sidiropoulos2017tensor}. Finally, directly recovering the missing entries in a tensor is more reasonable and reliable than recovering these entries in the unfolding of this tensor by matrix completion algorithms \cite{yuan2016tensor}.

For LRMC, many algorithms have been proposed in the past decade. Some popular ones are nuclear norm (a convex relaxation of rank) minimization \cite{CandesRecht2009}, Schatten-$p$ quasi-norm minimization \cite{MC_ShattenP_AAAI125165,Nie2015}, and low-rank factorization \cite{srebro2005rank,NIPS2019_8754}. It is natural to extend these ideas to tensor completion \cite{liu2012tensor,gandy2011tensor,bazerque2013rank,zhang2016exact,kong2018t,ghadermarzy2019near}. Given an incomplete tensor, as mentioned before, one may recover the missing entries by performing matricization \cite{kolda2009tensor} along every mode and minimizing the nuclear norms of these matrices, which we call tensor unfoldings. Liu et al. \cite{liu2012tensor} proposed to minimize a weighted sum of the nuclear norms of the $d$ unfoldings. Tomioka et al. \cite{tomioka2013convex} proposed to minimize the sum of Schatten-$p$ norms of the unfoldings to recover the missing entries and provided theoretical guarantee for the recovery error. Yuan and Zhang \cite{yuan2016tensor} pointed out that these unfolding-based tensor completion algorithms require much more observed entries for exact recovery than direct tensor completion without matricization. They proposed to minimize the tensor nuclear norm defined as the dual norm of the tensor spectral norm, and showed that exact recovery of a size $n \times n\times n$ tensor can be obtained if the number of uniformly observed entries satisfy $\vert\Omega\vert\geq C(3r^2n+\sqrt{rn^3})\textup{polylog}(n)$, assuming that the ranks of the unfoldings are all $r$. Obviously, this is much less than that required by unfolding-based tensor completion algorithms ($Crn^2\log^2(n)$) \cite{yuan2016tensor}. Since tensor nuclear norm is hard to compute in practice, Barak and Moitra \cite{barak2016noisy} and Foster and Risteski \cite{foster2019sum} used the sum-of-squares hierarchy as a tractable relaxation of the tensor nuclear norm. Particularly, Foster and Risteski \cite{foster2019sum} showed that for noisy completion of a rank-$r$ $n\times n\times n$ tensor, the generalization bound is $\tilde{O}(r^2n^{3/2}/\vert\Omega\vert)$.
Kong et al. \cite{kong2018t} used the t-Schatten-$p$ quasi-norm ($0<p<1$) for low-rank tensor (3-order) recovery and provided a multi-tensor-Schatten-$p$ quasi-norm surrogate to transform the nonconvex objective into a sum of several convex objectives. They provided an error bound for robust tensor recovery. It is not clear whether the method can be extended to higher-order tensors. In addition, the generalization bound for tensor completion with multi-tensor-Schatten-$p$ quasi-norm surrogate remains unknown.

Another category of LRTC algorithms is based on tensor factorization or decomposition  \cite{karatzoglou2010multiverse,acar2011scalable,xu2013parallel,xu2013block,jain2014provable,kressner2014low,zhao2015bayesian,kasai2016low} and they have relatively lower computational costs than the nuclear norm minimization based algorithms and hence are suitable for big tensors. For instance, Karatzoglou et al. \cite{karatzoglou2010multiverse} extended the regularized low-rank matrix factorization to regularized low-rank tensor decomposition for collaborative filtering. Jain and Oh \cite{jain2014provable} proposed to learn a CP (CANDECOMP/PARAFAC) decomposition \cite{kolda2009tensor} from incomplete tensor via alternating minimization and proved that an $n\times n\times n$ tensor of rank $r$ can be recovered from $O(n^{3/2}r^5 \log^4(n))$ randomly sampled entries. Nevertheless, such a number of required entries is significantly larger than the intrinsic dimension of the tensor, say $3rn$. Ghadermarzy et al. \cite{ghadermarzy2019near} proposed to use their defined tensor max-qnorm and M-norm as a rank proxy for tensor completion, and proved that $O(r^{1.5d}dn)$ observed entries are sufficient to recover a $d$th-order tensor of rank $r$ and size $n$ with high probability. More recent works on tensor completion can be found in \cite{huang2015provable,barak2016noisy,zhou2017tensor,8000407,yuan2019tensor,liu2020tensor,jiang2020multi,wang2020robust}.

In decomposition-based LRTC methods (e.g. \cite{acar2011scalable,xu2013parallel,jain2014provable}), one has to determine or estimate the rank (CP or multilinear rank \cite{kolda2009tensor}) in advance. A too large or too small rank may lead to a large recovery error in tensor completion. Xu et al.\cite{xu2013parallel} suggested adjusting the rank adaptively and provided a parallelizable algorithm to further improve the computational efficiency of tensor completion. The algorithm may have high recovery error when the number of missing entries is large, which will be shown in Section \ref{sec_exp_LRTC}. Zhao et al. \cite{zhao2015bayesian} proposed a Bayesian CP decomposition method, in which the rank is estimated automatically. The method has high computational cost when the rank is high.

In LRMC, regularized matrix factorizations have been widely studied \cite{srebro2005rank,MC_ShattenP_AAAI125165, NIPS2019_8754}. Similar ideas can be applied to LRTC. A few researchers have considered adding rank regularization to the CP decomposition model. For instance, for 3rd-order tensor completion, Bazerque et al. \cite{bazerque2013rank} used the sum of squared Frobenius norms as regularizer and showed that it is related to the $\ell_{2/3}$ norm of the weights in CP decomposition. Yang et al. \cite{yang2016tensor} applied group-sparse regularization to LRTC and showed that the regularizer is related to the $\ell_{1/3}$ norm, which is shaper than the $\ell_{2/3}$ norm \cite{bazerque2013rank}. Shi et al. \cite{shi2017tensor} proposed to directly minimize the $\ell_1$ norm of the weights in CP decomposition for LRTC. In fact, the regularizers presented in the above three works are closely related to the tensor Schatten-$p$ (quasi) norm with $p=1$, $2/3$, or $1/3$. Notice that these works are purely empirical-motivated and have no theoretical guarantee on the tensor completion performance.  Moreover, the regularizations are limited in the discrete values $\lbrace 1, 2/3,1/3\rbrace$ and only for 3rd-order tensors. One may expect to exploit Schatten-$p$ quasi-norm with arbitrary $p$ on arbitrary-order tensor and have theoretical guarantees for recovery performance.

Besides LRTC, tensor robust principal component analysis (TRPCA) is another important problem of tensor recovery. TRPCA is a generalization of robust PCA \cite{de2001robust,ding2006r,RPCAcandes11,shang2017bilinear,haeffele2019structured} and aims to decompose a noisy tensor to the sum of a low-rank tensor and a sparse tensor.  Many recent works on TRPCA can be found in \cite{anandkumar2016tensor,lu2016tensor,8000407,zheng2019mixed,lu2019tensor,wang2020robust}. For example, Lu et al. \cite{lu2019tensor} defined a new tensor nuclear norm based on the $t$-product \cite{kilmer2011factorization} of tensors and provided sufficient conditions for exact recovery. Notice that these TRPCA algorithms have high computational costs on large-scale data. If we take advantage of factorization model (e.g. \cite{anandkumar2016tensor}) for fast TRPCA, we need to estimate the rank. In addition, one may use sharper rank proxies like the Schatten-$p$ quasi-norm \cite{MC_ShattenP_AAAI125165} to achieve higher recovery accuracy if it can be handled efficiently.

In this paper, we focus on fast, accurate LRTC and TRPCA. Our contributions are as follows.
\begin{itemize}
\item We propose a new class of regularizers as a tensor rank proxy. These are based on the Euclidean norms of the component vectors in the form of CP decomposition. We show that the regularizers are monotonic transformations of Schatten-$p$ quasi-norms, where $p$ could be any positive value and the tensor could have arbitrary order. Hence, the regularizers proposed in \cite{bazerque2013rank,yang2016tensor,shi2017tensor} are special cases of our framework.
\item We also provide asymmetric variational forms for the tensor Schatten-$p$ quasi-norm with discrete $p$. Each of the asymmetric regularizers could have only one nonconvex and nonsmooth term such that the optimization has fewer nonconvex and nonsmooth subproblems than the case of the symmetric regularizers for small $p$.
\item The regularizers enable us to minimize Schatten-$p$ quasi-norms efficiently without performing SVD and hence solve large tensor recovery problems more efficiently. The regularizers also provide high recovery accuracy when $p$ is much less than $1$. When applying the regularizers to LRTC and TRPCA, the recovery performance is robust to the choice of initial rank.
\item We provide a generalization error bound for LRTC, which is tighter than that of Foster and Risteski \cite{foster2019sum}. The bound applies to any CP-decomposition-based LRTC models and indicates that one can obtain a tighter bound by reducing the number and Euclidean norms of CP component vectors.
\item We prove that when the number of observed entries is large enough, a smaller $p$ (but not too small) in the tensor Schatten-$p$ quasi-norm leads to a tighter generalization error bound. This guides the practical selection of sharper regularizers in LRTC from a theoretical perspective.
\item We prove that a smaller $p$ (but not too small) in the tensor Schatten-$p$ quasi-norm leads to a tighter error bound for TRPCA. This guarantees the superiority of sharper rank regularizers in recovering low-rank tensors corrupted by sparse noises.
\end{itemize}
The numerical results of LRTC and TRPCA on synthetic data, image inpainting, and image denoising corroborate the effectiveness and superiority of our methods over state-of-the-art baseline methods.

This paper is organized as follows. Section \ref{sec:CVEN} introduces the proposed regularizers and their properties. Sections \ref{sec:app_ENR_LRTC} and  \ref{sec:app_ENR_TRPCA} are the applications of the regularizers in LRTC and TRPCA, respectively, as well as the optimization strategies. Section \ref{sec:bound} analyzes the LRTC error bound. Section \ref{sec:bound_TRPCA} provides a recovery guarantee for TRPCA. Sections \ref{sec_exp_LRTC} and \ref{sec_exp_TRPCA} are the respective experiments of LRTC and TRPCA. Section \ref{sec:conclusion} draws a conclusion for this paper. Throughout this paper, we use the notations listed in Table \ref{table_notation}.

\begin{table}[h]
\centering
\caption{Mathematical notations}
\vspace{3pt}
\begin{tabular}{l|l||l|l}
\hline
$\mathbb{R}$ & real number & $\mathbb{N}$ & natural number   \\ \hline
$\bm{X}$ & matrix &$\bm{\mathcal{X}}$ & tensor\\ \hline
$\bm{I}$& identity matrix 	&$\bm{\mathcal{I}}$ & identity tensor\\ \hline
$\bm{x}$ & column vector &  $[d]$ & $\lbrace 1,2,\ldots,d\rbrace$ \\ \hline
$P[\cdot]$ & probability & $\mathcal{P}_{\Omega}(\cdot)$  & linear projection onto $\Omega$\\ \hline
$\Omega$ & index set & $\vert\Omega\vert$  & cardinality of $\Omega$\\ \hline
$\mathcal{L}$ & loss function & $\mathcal{S}$ & matrix or tensor sets \\ \hline
$\bm{x}_i^{(j)}$ & \multicolumn{3}{l}{column vector with index $i$ and $j$}\\ \hline
$\bm{X}_{(j)}$ & \multicolumn{3}{l}{mode-$j$ matricization of $\bm{\mathcal{X}}$}\\ \hline
$\circ$ & \multicolumn{3}{l}{outer product of vectors} \\ \hline
$\ast$ & \multicolumn{3}{l}{Hadamard product of vectors, matrices, or tensors} \\ \hline
$\odot$ & \multicolumn{3}{l}{Khatri-Rao product  of vectors or matrices} \\ \hline
$\langle\cdot,\cdot\rangle$& \multicolumn{3}{l}{inner product of vectors, matrices, or tensors}\\ \hline
$\times_j$ & \multicolumn{3}{l}{mode-$j$ product of a tensor and a matrix} \\ \hline
$\Vert \cdot\Vert$& \multicolumn{3}{l}{Euclidean norm or $\ell_2$ norm of vector}\\ \hline
$\Vert \cdot\Vert_F$& \multicolumn{3}{l}{Frobenius norm of matrix or tensor}\\ \hline
$\Vert \cdot\Vert_1$& \multicolumn{3}{l}{$\ell_1$ norm of vector, matrix, or tensor}\\ \hline
$\Vert \cdot\Vert_2$& \multicolumn{3}{l}{spectral norm of matrix or tensor}\\ \hline
$\Vert \cdot\Vert_\ast$& \multicolumn{3}{l}{nuclear norm of matrix or tensor}\\ \hline
$\Vert \cdot\Vert_{Sp}$& \multicolumn{3}{l}{Schatten-$p$ norm or quasi-norm of matrix or tensor}\\ \hline
$\Vert \cdot\Vert_\infty$& \multicolumn{3}{l}{the maximal absolute entry of matrix or tensor}\\ \hline
$\mathcal{N}(\mu,\sigma^2)$ & \multicolumn{3}{l}{Gaussian distribution with mean $u$ and variance $\sigma^2$} \\ \hline
\end{tabular}\label{table_notation}
\end{table}

\section{Euclidean Regularization}\label{sec:CVEN}
We introduce the following definitions.
\begin{definition}
Tensor CANDECOMP/PARAFAC (CP) decomposition factorizes a tensor $\bm{\mathcal{X}}\in\mathcal{R}^{n_1\times n_2\ldots\times n_d}$ into a sum of component rank-one tensors:
$$\bm{\mathcal{X}}=\sum_{i=1}^r\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)},$$
where $\bm{x}_i^{(j)}\in\mathbb{R}^{n_j\times 1}$, $i\in[r]$, $j\in[d]$.
See Figure \ref{fig_cp}.
\begin{figure}[h]
\centering
\includegraphics[width=8.5cm]{fig_cp}
\caption{CP decomposition of a 3rd-order tensor.}\label{fig_cp}
\end{figure}

\end{definition}
\begin{definition}
The rank  of a tensor $\bm{\mathcal{X}}$ is
defined as the minimum number of rank-one tensors that sum to $\bm{\mathcal{X}}$:
$$
\textup{rank}(\bm{\mathcal{X}})=\min\left\lbrace r\in\mathbb{N}:\ \bm{\mathcal{X}}=\sum_{i=1}^r\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)}\right\rbrace.
$$
\end{definition}
In contrast to matrix, it is NP-hard to determine the rank of a tensor. In practice, the rank of a tensor is numerically determined or estimated by fitting various CP rank-$r$ decomposition\footnote{An exact CP decomposition with $r = \textup{rank}(\bm{\mathcal{X}})$ components is called \textit{rank decomposition} \cite{kolda2009tensor}.} models such as alternating least squares \cite{carroll1970analysis,harshman1970foundations}. Another major difference between tensors and matrices is the rank decompositions of tensors are often unique but the rank decompositions of matrices are not unique \cite{kolda2009tensor}.
\begin{definition}
The Frobenius norm of a tensor $\bm{\mathcal{X}}$ is defined as the root of squared sum of its elements:
\begin{equation*}
\begin{aligned}
\Vert \bm{\mathcal{X}}\Vert_F=\sqrt{\left\langle \bm{\mathcal{X}}, \bm{\mathcal{X}}\right\rangle}=\left(\sum_{i_1=1}^{n_1}\ldots\sum_{i_d=1}^{n_d}x^2_{i_1\ldots i_d}\right)^{1/2}.
\end{aligned}
\end{equation*}
\end{definition}
\begin{definition}\label{def_tnn}
The nuclear norm \cite{friedland2018nuclear} of a tensor $\bm{\mathcal{X}}$ is defined as
\begin{equation*}
\begin{aligned}
\Vert \bm{\mathcal{X}}\Vert_\ast=\inf\Bigg\lbrace&\sum_{i=1}^r\vert s_i\vert:\ \bm{\mathcal{X}}=\sum_{i=1}^r s_i\bm{u}_i^{(1)}\circ\bm{u}_i^{(2)}\ldots\circ\bm{u}_i^{(d)},\\
& \Vert\bm{u}_i^{(j)}\Vert=1, \ r\in\mathbb{N}\Bigg\rbrace.
\end{aligned}
\end{equation*}
\end{definition}
Notice that this definition of the tensor nuclear norm is commonly used in mathematics and machine learning, though there is an alternate definition based on the sum of nuclear norms of matricizations of $\bm{\mathcal{X}}$ \cite{liu2012tensor}. In this paper, Definition \ref{def_tnn} defines the tensor nuclear norm. It is also NP-hard to compute the tensor nuclear norm \cite{friedland2018nuclear}.

Using Definition \ref{def_tnn}, we see the tensor nuclear norm is a convex relaxation of tensor rank. Similarly, we may define a  tensor Schatten-$p$ quasi-nor
that extends matrix Schatten-$p$ quasi-norm. The $p$-th power of the tensor Schatten-$p$ quasi-norm is a nonconvex relaxation of tensor rank.
\begin{definition}\label{def_tspn}
The tensor Schatten-$p$ quasi-norm ($0<p<1$) of a tensor $\bm{\mathcal{X}}$ is defined as
\begin{small}
\begin{equation*}
\begin{aligned}
\Vert \bm{\mathcal{X}}\Vert_{S_p}=\inf\Bigg\lbrace&\left(\sum_{i=1}^r\vert s_i\vert^p\right)^{\tfrac{1}{p}}: \ \bm{\mathcal{X}}=\sum_{i=1}^r s_i\bm{u}_i^{(1)}\circ\bm{u}_i^{(2)}\ldots\circ\bm{u}_i^{(d)},\\
& \Vert\bm{u}_i^{(j)}\Vert=1, \ r\in\mathbb{N}, \ 0<p<1\Bigg\rbrace.
\end{aligned}
\end{equation*}
\end{small}
\end{definition}

According to Definition \ref{def_tnn}, we see that, unlike matrix norms, it is difficult to directly compute the tensor nuclear norm. To tackle this problem, Barak and Moitra \cite{barak2016noisy} and Foster and Risteski \cite{foster2019sum} proposed to use the sum-of-squares hierarchy as a tractable relaxation of tensor nuclear norm for tensor completion. The definition and computation of the sum-of-squares hierarchy are complex and rely on a hyper-parameter. Moreover, tensor Schatten-$p$ quasi-norm (Definitions \ref{def_tspn}) remains intractable.

In this paper, for convenience, we define an operator that forms a tensor $\bm{\mathcal{X}}$ from a collection of vectors $\bm{x}$.
\begin{definition}
Let $k,d\in\mathbb{N}$. Define the inverse CP operator as
$$\mathcal{CP}_k^d(\bm{x})=\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\cdots\circ\bm{x}_i^{(d)}.$$
\end{definition}

We provide the following variational form of the tensor Schatten-$p$ (quasi) norm ($0<p\leq 1$):
\begin{theorem}[Symmetric Regularizer]\label{theorem_sp}
The tensor Schatten-$p$ (quasi) norm can be represented as a function of the Euclidean norms of the CP components:
\begin{equation*}
\begin{aligned}
\Vert \bm{\mathcal{X}}\Vert_{S_p}^p=\inf\Bigg\lbrace  &\dfrac{1}{d}\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^{pd}:   \bm{\mathcal{X}}=\mathcal{CP}_k^d(\bm{x}),\ 0<p\leq 1\Bigg\rbrace.
\end{aligned}
\end{equation*}
\end{theorem}
Compared to the definition (Definition \ref{def_tspn}) of the tensor Schatten-$p$ quasi-norm, $\Vert \bm{\mathcal{X}}\Vert_{S_p}$ given by Theorem \ref{theorem_sp} internalizes the weight factors $\bm{s}$ into the vectors $\bm{u}$ and removes the constraints on $\bm{u}$. This reduces decision variables and the computational cost of calculating the gradient in optimization. In addition, for some choice of $p$, there is no nonsmooth functions on the factors in Theorem \ref{theorem_sp}, while there are always nonsmooth functions in Definition \ref{def_tspn}. Hence, the formulation of Schatten-$p$ quasi-norm in  Theorem \ref{theorem_sp} is more tractable than that in Definition \ref{def_tspn}.
For instance, in Theorem \ref{theorem_sp}, letting $p=1$, $\tfrac{2}{d}$, or $\tfrac{1}{d}$, we have the following rank regularizers:
\begin{corollary}\label{corollary_sp_sym}
Characterizations of the tensor nuclear norm and Schatten-$p$ quasi-norm:
\begin{equation*}
\begin{aligned}[r]
\textit{(1)}
\ \Vert \bm{\mathcal{X}}\Vert_\ast=\inf_{ \bm{\mathcal{X}}=\mathcal{CP}_k^d(\bm{x})}\  &\dfrac{1}{d}\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^d;\qquad \qquad \qquad \qquad
\end{aligned}
\end{equation*}
\begin{equation*}
\begin{aligned}
\textit{(2)}
\ \Vert \bm{\mathcal{X}}\Vert_{S_{2/d}}^{2/d}=\inf_{ \bm{\mathcal{X}}=\mathcal{CP}_k^d(\bm{x})}\ & \dfrac{1}{d}\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^2;\qquad \qquad \qquad \qquad
\end{aligned}
\end{equation*}
\begin{equation*}
\begin{aligned}
\textit{(3)}
\ \Vert \bm{\mathcal{X}}\Vert_{S_{1/d}}^{1/d}=\inf_{ \bm{\mathcal{X}}=\mathcal{CP}_k^d(\bm{x})}\  &\dfrac{1}{d}\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert.\qquad \qquad \qquad \qquad
\end{aligned}
\end{equation*}
\end{corollary}
These three special cases given by Corollary \ref{corollary_sp_sym} are based on convex functions of Euclidean norms of component vectors, and hence are possibly easier to handle in optimization.

According to Theorem~\ref{theorem_sp}, for a relatively low-order tensor, when we want to obtain a sharp enough regularizer, we have to use a sufficiently small $p$ for all $i$ and $j$. Thus, every term in $\lbrace \Vert\bm{x}_i^{(j)}\Vert^{pd} \rbrace_{i\in[k]}^{j\in[d]}$ is nonconvex and nonsmooth, making it difficult to solve the optimization problem of low-rank tensor recovery. The following theorem provides a class of asymmetric regularizers that have fewer nonconvex  and nonsmooth terms than those in Theorem \ref {theorem_sp}.

\begin{theorem}[Asymmetric Regularizer]\label{theorem_sp_asym}
Suppose $q\in\lbrace 1,1/2,1/3,1/4,\ldots\rbrace$. Let $p_1=q/(1+qd-q)$ and $p_2=2q/(2+qd-q)$. We have:
\begin{equation*}
\textit{(a)}\ \Vert \bm{\mathcal{X}}\Vert_{S_{p_1}}^{p_1}
=\inf_{ \bm{\mathcal{X}}=\mathcal{CP}_k^d(\bm{x})}\   p_1\sum_{i=1}^k\Big(\dfrac{1}{q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert\Big);\qquad \qquad \qquad \qquad
\end{equation*}
\begin{equation*}
\begin{aligned}
\textit{(b)}\ \Vert \bm{\mathcal{X}}\Vert_{S_{p_2}}^{p_2}
=\inf_{ \bm{\mathcal{X}}=\mathcal{CP}_k^d(\bm{x})}\  p_2\sum_{i=1}^k\Big(\dfrac{2}{q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert^2\Big).\qquad \qquad \qquad \qquad
\end{aligned}
\end{equation*}
\end{theorem}
In Theorem  \ref{theorem_sp_asym} (b), the terms of $j\geq 2$ are convex and smooth while those in Theorem  \ref{theorem_sp_asym} (a) are nonsmooth. Therefore, the optimization related to Theorem  \ref{theorem_sp_asym} (b) is easier.

Since $3$rd-order tensors (e.g. color images and videos) are more prevalent than tensors with other orders, we list their symmetric and asymmetric regularizers with only convex terms in Table \ref{tab_reg3} for convenience. Note that the last two regularizers in the table are not the consequences of Theorem  \ref{theorem_sp_asym}. The derivations are in Appendix \ref{proof_S25asym}. The first three regularizers in Table \ref{tab_reg3} are equivalent to the regularizers in  \cite{shi2017tensor}, \cite{bazerque2013rank}, and  \cite{yang2016tensor} respectively. But our motivation is different: we generalize the matrix Schatten-$p$ (quasi) norms ($p=1$, $2/3$, or $1/3$) to tensors while they aimed to penalize the decomposition factors and did not consider Schatten-$p$ (quasi) norm.

\begin{table}[h]
\centering
\caption{Symmetric and asymmetric regularizers with only convex terms for 3rd-order tensor ($\bm{\mathcal{X}}=\mathcal{CP}_k^3(\bm{x})$).}
\begin{tabular}{c|c|c}
\hline
& Schatten regularizer  & Characterization based on Euclidean norm\\ \hline
\multirow{3}{*}{\rotatebox{90}{Symmetric}} &$\Vert \bm{\mathcal{X}}\Vert_\ast$&$\tfrac{1}{3}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^3+\Vert\bm{x}_i^{(2)}\Vert^3+\Vert\bm{x}_i^{(3)}\Vert^3\right)$\\
&$\Vert \bm{\mathcal{X}}\Vert_{S_{2/3}}^{2/3}$&$\tfrac{1}{3}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^2+\Vert\bm{x}_i^{(2)}\Vert^2+\Vert\bm{x}_i^{(3)}\Vert^2\right)$\\
&$\Vert \bm{\mathcal{X}}\Vert_{S_{1/3}}^{1/3}$&$\tfrac{1}{3}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert+\Vert\bm{x}_i^{(2)}\Vert+\Vert\bm{x}_i^{(3)}\Vert\right)$\\ \hline
\multirow{2}{*}{\rotatebox{90}{Asymmetric}} & $\Vert \bm{\mathcal{X}}\Vert_{S_{1/2}}^{1/2}$&$\tfrac{\sqrt{2}}{4}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^{2}+\Vert\bm{x}_i^{(2)}\Vert^{2}+\Vert\bm{x}_i^{(3)}\Vert\right)$\\
&$\Vert \bm{\mathcal{X}}\Vert_{S_{2/5}}^{2/5}$&$\tfrac{16^{1/5}}{5}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^{2}+\Vert\bm{x}_i^{(2)}\Vert+\Vert\bm{x}_i^{(3)}\Vert\right)$\\
&$\Vert \bm{\mathcal{X}}\Vert_{S_{3/7}}^{3/7}$&$\tfrac{81^{1/7}}{7}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^{3}+\Vert\bm{x}_i^{(2)}\Vert+\Vert\bm{x}_i^{(3)}\Vert\right)$\\ \hline
\end{tabular}\label{tab_reg3}
\end{table}

In sum, the regularizers we proposed in this section covers all Schatten-$p$ (quasi) norms with any $0<p\leq 1$ for any-order tensors. Some of the regularizers, especially those asymmetric ones, are based on convex or/and smooth functions on the tensor factors, which provide convenience for application and optimization. These regularizers can be applied to low-rank tensor completion and robust tensor PCA that enjoy a variety of real applications in many areas such machine learning, computer vision, and neuroscience.

It is worth noting that the tensor regularizers we presented in this paper are closely related to the variational forms of matrix nuclear norm \cite{srebro2005maximum} and Schatten-$p$ quasi-norms \cite{shang2016tractable,shang2017bilinear,NIPS2019_8754}. For instance, in \cite{srebro2005maximum}, the squared sum of Frobenius norms of two factors of a matrix is lower bounded by the nuclear norm of the matrix, which is a special case of our Theorem \ref{theorem_sp} with $d=2$ and $p=1$. In \cite{shang2017bilinear}, the authors showed that the sum of nuclear norms of two factors of a matrix is a monotonic transformation of matrix Schatten-1/2 quasi-norm. In \cite{NIPS2019_8754}, the authors provided a class of SVD-free variational forms of the matrix Schatten-$p$ quasi-norm with discrete $p$ that can be arbitrarily small. Our tensor regularizers can be regarded as a generalization of the variational form of Schatten-$p$ quasi-norm from matrix to tensor. Nevertheless, theoretical guarantees about these tensor regularizers in low-rank tensor recovery are more difficult to derive than their matrix counterparts. To tackle the problem, in Section \ref{sec:bound}, we provide theoretical guarantees for the regularizers in LRTC.

\section{Low-Rank Tensor Completion}\label{sec:app_ENR_LRTC}
\subsection{LRTC-ENR Model}\label{sec:LRTC}
Let $\bm{\mathcal{T}}\in\mathbb{R}^{n_1\times n_2\ldots\times n_d}$ be a rank-$r$ tensor. Suppose we observed a few noisy entries of $\bm{\mathcal{T}}$ randomly (without replacement):
\begin{equation}\label{eq_LRTC_ob}
	[\bm{\mathcal{D}}]_{j_1 j_2 \ldots j_d}=[\bm{\mathcal{T}}]_{j_1 j_2 \ldots j_d}+[\bm{\mathcal{N}}]_{j_1 j_2 \ldots j_d}, \quad (j_1, j_2, \ldots, j_d) \in \Omega
\end{equation}
where $\Omega$ consists of the locations of the observed entries and each entry of the noise tensor $\bm{\mathcal{N}}$ is drawn from $\mathcal{N}(0,\sigma^2)$. To recover $\bm{\mathcal{T}}$ from $\bm{\mathcal{D}}$, one may solve
\begin{align}\label{problem_LRTC_Sp}
\mathop{\text{minimize}}_{\bm{\mathcal{X}}}\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}\left(\bm{\mathcal{D}}-\bm{\mathcal{X}}\right)\Vert_F^2+\lambda\Vert \bm{\mathcal{X}}\Vert_{S_p}^p,
\end{align}
of which the solution is an estimate of $\bm{\mathcal{T}}$. However, Problem  \eqref{problem_LRTC_Sp} is intractable.  Instead, based on the analysis in Section \ref{sec:CVEN}, we propose to solve
\begin{equation}\label{problem_LRTC_CVEN}
\begin{aligned}
\mathop{\text{minimize}}_{\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}}\ &\dfrac{1}{2}\left\Vert \bm{\mathcal{M}}\ast\left(\bm{\mathcal{D}}-\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)}\right)\right\Vert_F^2\\
&+\lambda \mathcal{R}\left(\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}\right),
\end{aligned}
\end{equation}
where $\bm{\mathcal{M}}$ is a binary tensor with $0$ for missing entries and $1$ for observed entries and the sign ``$\ast$" denotes the element-wise product of tensors.
$\mathcal{R}\left(\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}\right)$ denotes a certain regularizer in Theorem \ref{theorem_sp}, Corollary \ref{corollary_sp_sym}, Theorem \ref{theorem_sp_asym}, or Table \ref{tab_reg3}. For example,
\begin{equation}\label{reg_exp}
\mathcal{R}\left(\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}\right)=\lambda\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^{pd},
\end{equation}
in which we omit the constant for simplicity. Note that even finding a stationary point of \eqref{problem_LRTC_Sp} is
hard due to computation of the Schatten-$p$ quasi-norm while finding
a stationary point of (3) is computationally feasible.
For convenience, we call problem \eqref{problem_LRTC_CVEN} Low-Rank Tensor Completion with Euclidean Norm Regularization (\textbf{LRTC-ENR}).

\subsection{Optimization for LRTC-ENR}
\subsubsection{Block Coordinate Descent with Extrapolation}\label{sec_LRTC_opt1}
There are $d$ blocks of decision variables in Problem \eqref{problem_LRTC_CVEN}, i.e. $\lbrace \bm{X}^{(j)}\rbrace^{j\in[d]}$, where $\bm{X}^{(j)}=[\bm{x}_1^{(j)},\bm{x}_2^{(j)},\ldots,\bm{x}_k^{(j)}]$. We propose to find a critical point of  \eqref{problem_LRTC_CVEN} by Block Coordinate Descent (BCD) with Extrapolation (BCDE for short) \cite{xu2013block}, which is more efficient than BCD. For simplicity, here we only present the optimization of \eqref{problem_LRTC_CVEN} with the regularizer shown in \eqref{reg_exp}. The optimization can be easily extended to \eqref{problem_LRTC_CVEN} with other regularizers we proposed in Section \ref{sec:CVEN}.

Let
\begin{align*}
&\mathcal{L}\left(\lbrace \bm{X}^{(j)}\rbrace^{j\in[d]}\right)\\
\triangleq&\dfrac{1}{2}\left\Vert \bm{\mathcal{M}}\ast\left(\bm{\mathcal{D}}-\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)}\right)\right\Vert_F^2\\
=&\dfrac{1}{2}\left\Vert \bm{M}_{(j)}\ast\left(\bm{D}_{(j)}-\bm{X}^{(j)}[(\bm{X}^{(i)})^{\odot_{i\neq j}}]^\top\right)\right\Vert_F^2
\end{align*}
and
$$\mathcal{R}_{pd}(\bm{X}^{(j)})\triangleq\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^{pd},$$
where $[(\bm{X}^{(i)})^{\odot_{i\neq j}}]=\bm{X}^{(d)}\odot\ldots\odot\bm{X}^{(j+1)}\odot\bm{X}^{(j-1)}\odot\ldots\odot\bm{X}^{(1)}$.
We initialize $\lbrace \bm{X}^{(j)}\rbrace^{j\in[d]}$ randomly and rescale $\bm{x}$ to have unit Euclidean norm:
\begin{equation}\label{eq_x_ini}
\bm{x}_i^{(j)}\sim\mathcal{N}(0,1), \quad \bm{x}_i^{(j)}\leftarrow \bm{x}_i^{(j)}/\Vert \bm{x}_i^{(j)}\Vert, \quad i\in[k],\ j\in[d].
\end{equation}
Then at iteration $t$, for $j=1,2,\ldots,d$, we first perform an extrapolation
\begin{equation}\label{eq_extrap}
 \hat{\bm{X}}^{(j)}_{t-1}=\bm{X}^{(j)}_{t-1}+\omega_{j,t-1}(\bm{X}^{(j)}_{t-1}-\bm{X}^{(j)}_{t-2}),
\end{equation}
where $\omega_{j,t-1}\geq 0$ controls the size of the extrapolation at iteration $t$. Then we perform a proximal step as
\begin{equation}\label{eq_prox}
\begin{aligned}
\bm{X}^{(j)}_t=&\mathop{\textup{argmin}}_{\bm{X}^{(j)}}\ \left\langle \bm{\nabla}_{\hat{\bm{X}}_{t-1}^{(j)}}\mathcal{L},\ \bm{X}^{(j)}-\hat{\bm{X}}^{(j)}_{t-1}\right\rangle\\
&\qquad\quad\ \ +\dfrac{\hat{L}_{j,t}}{2}\Vert\bm{X}^{(j)}-\hat{\bm{X}}^{(j)}_{t-1} \Vert_F^2+\lambda\mathcal{R}_{pd}(\bm{X}^{(j)})\\
:=&\textup{prox}_{pd}^{\lambda}\Big(\hat{\bm{X}}^{(j)}_{t-1}-\hat{L}_{j,t}^{-1}\bm{\nabla}_{\hat{\bm{X}}_{t-1}^{(j)}}\mathcal{L}\Big).
\end{aligned}
\end{equation}
In \eqref{eq_prox}, $\bm{\nabla}_{\hat{\bm{X}}_{t-1}^{(j)}}\mathcal{L}=\left(\bm{M}_{(j)}\ast\big(\bm{D}_{(j)}-\hat{\bm{X}}_{t-1}^{(j)}[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]^\top\big)\right)\big(-[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\big)$, $[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]=\bm{X}_{t-1}^{(d)}\odot\cdots\odot\bm{X}_{t-1}^{(j+1)}\odot\bm{X}_t^{(j-1)}\odot\cdots\odot\bm{X}_t^{(1)}$, and $\hat{L}_{j,t}$ is larger than the Lipschitz constant of $\bm{\nabla}_{\hat{\bm{X}}_{t-1}^{(j)}}\mathcal{L}$. We approximate it by
\begin{equation}\label{eq_estimateL}
\hat{L}_{j,t}=\varrho\sqrt{\dfrac{\vert\Omega\vert}{\prod_{i=1}^dn_i}}\left\Vert [(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\right\Vert_2^2,
\end{equation}
where $\varrho=0.5$, $1$, or $2$. The derivation of \eqref{eq_estimateL} is in Appendix \ref{appendix_Lips}.

In \eqref{eq_prox}, when $p=1/d$, $\textup{prox}_{1}^{\lambda}\big(\bm{Y}\big)=\Phi_\lambda(\bm{Y})$, where $\Phi$ is the column-wise soft-thresholding operator \cite{parikh2014proximal}  defined by
\begin{equation*}
\Phi_\lambda(\bm{y})=\left\{
\begin{array}{ll}
\tfrac{(\Vert\bm{y}\Vert-\lambda)\bm{y}}{\Vert\bm{y}\Vert}, &\textup{if}\ \Vert\bm{y}\Vert>\lambda;\\
\bm{0}, &\textup{otherwise.}
\end{array}
\right.
\end{equation*}
When $p=2/d$, $\textup{prox}_{2}^{\lambda}\big(\bm{Y}\big)=\dfrac{\hat{L}_{j,t}}{\hat{L}_{j,t}+2\lambda}\bm{Y}$. When $p>2/d$, we can estimate $\bm{X}^{(j)}_t$ by gradient descent. When $p<1/d$, \eqref{eq_prox} is nonconvex and nonsmooth. Then we update $\bm{X}^{(j)}$ with iteratively reweighted method \cite{lu2014iterative}, which is given by Algorithm \ref{alg_q_IRLS}.

\begin{algorithm}[h]
\caption{$\mathop{\textup{min}}_{\bm{Y}}\tfrac{1}{2}\Vert \bm{Y}-\bm{G}\Vert_F^2+\tilde{\lambda}\sum_{i=1}^k\Vert\bm{y}_i^{(1)}\Vert^{q}$}
\label{alg_q_IRLS}
\begin{algorithmic}[1]
\Require
$\bm{G}$, $q$, $\tilde{\lambda}$, $t_q$, $\epsilon$.
\State $\bm{Y}\leftarrow\bm{G}$.
	  \For{$t=1,2,\ldots,t_{q}$}
	\State $\bm{W}=\textup{diag}\big((\Vert\bm{y}_1^{(1)}\Vert+\epsilon)^{\tfrac{q-2}{2}},\ldots,(\Vert\bm{y}_k^{(1)}\Vert+\epsilon)^{\tfrac{q-2}{2}}\big)$.
	\State $\bm{Y}=\bm{G}(\bm{I}+2\tilde{\lambda}\bm{W}\bm{W}^T)^{-1}$.
      \EndFor
\Ensure $\bm{Y}$.
\end{algorithmic}
\end{algorithm}

In \eqref{eq_extrap}, the parameter $\omega_{jt}$ is determined as
\begin{equation}\label{eq_omega_jt}
\omega_{j,t-1}=\delta\sqrt{\hat{L}_{j,t-2}/\hat{L}_{j,t-1}},
\end{equation}
where $\delta<1$. We set $\delta=0.95$ for simplicity. The whole procedure is summarized in Algorithm \ref{alg_LRTC}. The convergence analysis can be found in \cite{xu2013block}.

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
\begin{algorithm}[h]
\caption{solve LRTC-ENR by BCDE}
\label{alg_LRTC}
\begin{algorithmic}[1]
\Require
$\bm{\mathcal{D}}$, $\bm{\mathcal{M}}$, $k$, $\lambda$, $t_\mathrm{max}$.
\State Initialize $\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}$ with \eqref{eq_x_ini}, let $t=0$.
\Repeat
\State $t\leftarrow t+1$.
	  \For{$j=1,2,\ldots,d$}
	  \If{$t<=2$}
	\State $\omega_{j,t-1}=0$.
	\Else
	\State Compute $\omega_{j,t-1}$ by \eqref{eq_omega_jt}.
	\EndIf
	\State Compute $\hat{\bm{X}}^{(j)}_{t-1}$ by \eqref{eq_extrap}.
	\State Compute $L_{j,t-1}$ by \eqref{eq_estimateL}.
	\State Compute $\bm{X}^{(j)}_{t}$ by \eqref{eq_prox}.
      \EndFor
\State Remove the zero columns of $\bm{X}^{(j)}_{t}$, $j\in[d]$.
\Until{converged or $t=t_\mathrm{max}$}
\Ensure $\hat{\bm{\mathcal{T}}}=\bm{\mathcal{I}}\times_{1}\bm{X}_t^{(1)}\times_{2}\bm{X}_t^{(2)}\ldots\times_{d}\bm{X}_t^{(d)}$.
\end{algorithmic}
\end{algorithm}

\subsubsection{LBFGS}
Though faster than BCD, the computational cost per iteration of BCDE is high when $d$ is not small because we need to compute $\bm{X}^{(j)}[(\bm{X}^{(i)})^{\odot_{i\neq j}}]^\top$ for $d$ times in every iteration.  One may consider the Jacobi-type iteration (cheaper computation but slower convergence) rather than the Gauss-Seidel iteration in BCD and BCDE. In practice, we can use quasi-Newton methods such as L-BFGS \cite{liu1989limited} to solve problem \eqref{problem_LRTC_CVEN} even though the objective function is nonsmooth. Particularly, we can drop the columns of $\bm{X}^{(j)}_{t}$ ($j\in[d]$) with nearly-zero Euclidean norms for acceleration. The corresponding algorithm\footnote{The implementation in this paper is based on the \textit{minFunc} MATLAB toolbox of M. Schmidt: \url{http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html}.} is shown in Algorithm \ref{alg_LRTC_lbfgs}.

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
\begin{algorithm}[h]
\caption{solve LRTC-ENR by LBFGS}
\label{alg_LRTC_lbfgs}
\begin{algorithmic}[1]
\Require
$\bm{\mathcal{D}}$, $\bm{\mathcal{M}}$, $k$, $\lambda$, $t_\mathrm{max}$.
\State Initialize $\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}$ with \eqref{eq_x_ini}, let $t=0$.
\Repeat
\State $t\leftarrow t+1$.
\State Compute the search directions by LBFGS.
\State Use line search to determine the step size.
\State Update $\bm{X}^{(j)}_{t}$, $j\in[d]$.
\State For $j\in[d]$, remove the columns of $\bm{X}^{(j)}_{t}$ with Euclidean norms less than a small threshold, e.g. $10^{-5}$.
\Until{converged or $t=t_\mathrm{max}$}
\Ensure $\hat{\bm{\mathcal{T}}}=\bm{\mathcal{I}}\times_{1}\bm{X}_t^{(1)}\times_{2}\bm{X}_t^{(2)}\ldots\times_{d}\bm{X}_t^{(d)}$.
\end{algorithmic}
\end{algorithm}

\subsubsection{LRTC model based on the definition of $\Vert \bm{\mathcal{X}}\Vert_{S_p}^p$?}
It is worth mentioning that according to Definition \ref{def_tspn}, one may consider the following LRTC problem
\begin{equation}\label{problem_LRTC_Sp_min}\tag{LRTC-Schatten-$p$}
\begin{aligned}
\mathop{\text{minimize}}_{\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}}\ \ &\dfrac{1}{2}\left\Vert \bm{\mathcal{M}}\ast\left(\bm{\mathcal{D}}-\sum_{i=1}^k s_i\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)}\right)\right\Vert_F^2\\
&+\lambda\sum_{i=1}^k\vert s_i\vert^{p},\\
\textup{subject to}\ \ & \Vert\bm{x}_i^{(j)}\Vert=1, \forall i\in[k], j\in[d].
\end{aligned}
\end{equation}
Nevertheless, it is more difficult to solve  \ref{problem_LRTC_Sp_min} than LRTC-ENR due to the following reasons. First, \ref{problem_LRTC_Sp_min} has one more block of variables $\bm{s}$, of which the gradient computation is costly. Second, the constraint in  \ref{problem_LRTC_Sp_min} hinders the application of quasi-Newton methods such as LBFGS. Heuristically, we suggest using BCDE with iteratively reweighted minimization embedded to solve \ref{problem_LRTC_Sp_min}. The procedures are similar to Algorithm \ref{alg_LRTC} and will not be detailed in this paper. \ref{problem_LRTC_Sp_min} is compared with our LRTC-ENR in Section \ref{sec_compare_p_syn}.

\section{Tensor Robust PCA}\label{sec:app_ENR_TRPCA}
\subsection{TRPCA-ENR Model}
Let $\bm{\mathcal{T}}\in\mathbb{R}^{n_1\times n_2 \times \cdots \times n_d}$ be a rank-$r$ tensor. Suppose $\bm{\mathcal{T}}$ is corrupted as
\begin{equation}\label{eq.TRPCA_model}
\bm{\mathcal{D}}=\bm{\mathcal{T}}+\bm{\mathcal{N}}+\bm{\mathcal{E}},
\end{equation}
where $\bm{\mathcal{N}}$ is a dense noise tensor drawn from $\mathcal{N}(0,\sigma^2)$ and $\bm{\mathcal{E}}$ is a sparse noise tensor with randomly distributed nonzero entries.
To recover $\bm{\mathcal{T}}$ from $\bm{\mathcal{D}}$, we wish to solve
\begin{align}\label{problem_RTPCA_Sp}
\mathop{\text{minimize}}_{\bm{\mathcal{X}},\bm{\mathcal{E}}}\dfrac{1}{2}\Vert \bm{\mathcal{D}}-\bm{\mathcal{X}}-\bm{\mathcal{E}}\Vert_F^2+\lambda_x\Vert \bm{\mathcal{X}}\Vert_{S_p}^p+\lambda_e\Vert\bm{\mathcal{E}}\Vert_1,
\end{align}
but the problem is intractable.  Then, based on the analysis in Section \ref{sec:CVEN}, we propose to solve
\begin{equation}\label{problem_RTPCA_CVEN}
\begin{aligned}
\mathop{\text{minimize}}_{\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]},\ \bm{\mathcal{E}}}\ &\dfrac{1}{2}\left\Vert \bm{\mathcal{D}}-\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)}-\bm{\mathcal{E}}\right\Vert_F^2\\
&+\lambda_x\mathcal{R}\left(\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}\right)+\lambda_e\Vert\bm{\mathcal{E}}\Vert_1,
\end{aligned}
\end{equation}
where $\lambda_e$ is a penalty parameter for the sparse tensor $\bm{\mathcal{E}}$.
For convenience, we call \eqref{problem_RTPCA_CVEN} \textbf{TRPCA-ENR}.

\subsection{Optimization for TRPCA-ENR}
When $p=2/d$, we use alternating minimization to solve the optimization of TRPCA-ENR because every subproblem has a closed-form solution.
When $p=1/d$, we may, like Section \ref{sec_LRTC_opt1}, use block coordinate descent with extrapolation \cite{xu2013block}. However, the additional variable $\bm{\mathcal{E}}$ further slows down the convergence. We hence propose to solve \eqref{problem_RTPCA_CVEN} by the (nonconvex) alternating direction method of multipliers (ADMM) \cite{wang2015global}. Specifically, by adding auxiliary variables $\lbrace\bm{Y}^{(j)}\rbrace_{j=1}^d$, we reformulate \eqref{problem_RTPCA_CVEN} as
\begin{align*}
\mathop{\textup{minimize}}_{\lbrace\bm{X}^{(j)},\bm{Y}^{(j)}\rbrace_{j=1}^d,\bm{\mathcal{E}}}\ \ &\dfrac{1}{2}\left\Vert \bm{D}_{(j)}-\bm{X}^{(j)}[(\bm{X}^{(i)})^{\odot_{i\neq j}}]^\top-\bm{E}_{(j)}\right\Vert_F^2\\
&+\lambda_x\mathcal{R}\left(\lbrace\bm{y}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}\right)+\lambda_e\Vert\bm{\mathcal{E}}\Vert_1\\
\textup{subject to\quad\ }\ \ & \bm{Y}^{(j)}=\bm{X}^{(j)}, j=1,2,\ldots,d.
\end{align*}
Let $\lbrace\bm{Z}^{(j)}\rbrace_{j=1}^d$ be Lagrange multipliers and solve
\begin{equation}\label{eq.admm_f}
\begin{aligned}
&\mathop{\textup{minimize}}_{\lbrace\bm{X}^{(j)},\bm{Y}^{(j)}\rbrace_{j=1}^d,\bm{\mathcal{E}}}\ \ \dfrac{1}{2}\left\Vert \bm{D}_{(j)}-\bm{X}^{(j)}[(\bm{X}^{(i)})^{\odot_{i\neq j}}]^\top-\bm{E}_{(j)}\right\Vert_F^2\\
&\quad+\lambda_x\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{y}_i^{(j)}\Vert^{pd}+\lambda_e\Vert\bm{\mathcal{E}}\Vert_1\\
&\quad+\sum_{j=1}^d\left\langle \bm{Y}^{(j)}-\bm{X}^{(j)},\bm{Z}^{(j)}\right\rangle+\dfrac{\mu}{2}\Vert \bm{Y}^{(j)}-\bm{X}^{(j)}\Vert_F^2,
\end{aligned}
\end{equation}
where $\mu$ is the augmented Lagrange penalty parameter. Then update $\lbrace\bm{X}^{(j)},\bm{Y}^{(j)}\rbrace_{j=1}^d$ and $\bm{\mathcal{E}}$ sequentially to minimize \eqref{eq.admm_f} and update $\lbrace\bm{Z}^{(j)}\rbrace_{j=1}^d$ lastly. The procedures are summarized into Algorithm \ref{alg_RTPCA}, where
\begin{equation}\label{eq_trpca_x}
\begin{aligned}
\bm{X}^{(j)}_t=&\Big((\bm{D}_{(j)}-\bm{E}_{(j)})[(\bm{X}^{(i)})^{\odot_{i\neq j}}]+\mu\bm{Y}^{(j)}_{t-1}-\bm{Z}^{(j)}\Big)\\
&\times \Big([(\bm{X}^{(i)})^{\odot_{i\neq j}}]^T[(\bm{X}^{(i)})^{\odot_{i\neq j}}]+\mu\bm{I}\Big)^{-1},
\end{aligned}
\end{equation}
\begin{equation}\label{eq_trpca_y}
\bm{Y}^{(j)}_t=\Phi_{\lambda_x/\mu}\left(\bm{X}^{(j)}_t-\bm{Z}^{(j)}/\mu\right),
\end{equation}
\begin{equation}\label{eq_trpca_z}
\bm{Z}^{(j)}\longleftarrow \bm{Z}^{(j)}+\mu(\bm{Y}^{(j)}_t-\bm{X}^{(j)}_t),
\end{equation}
\begin{equation}\label{eq_soft_thresh}
\bm{\mathcal{E}}_t=\Psi_{\lambda_e}\left( \bm{\mathcal{D}}-\bm{\mathcal{I}}\times_{1}\bm{X}_t^{(1)}\times_{2}\bm{X}_t^{(2)}\ldots\times_{d}\bm{X}_t^{(d)}\right).
\end{equation}
$\Psi$ is the element-wise soft-thresholding operator \cite{parikh2014proximal}  defined by
$$\Psi_{\lambda_e}(v)=\textup{sign}(v)\max(0,\vert v\vert-{\lambda_e}).$$

\begin{algorithm}[h]
\caption{RTPCA-ENR ($p=1/d$) solved by ADMM}
\label{alg_RTPCA}
\begin{algorithmic}[1]
\Require
$\bm{\mathcal{D}}$, $k$, $\lambda_x$, $\lambda_e$, $\mu$, $t_\mathrm{max}$.
\State Initialize $\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}$ with \eqref{eq_x_ini}; for $j=1,\ldots, d$, let $\bm{Y}^{(j)}_0=\bm{X}^{(j)}_0$ and $\bm{Z}^{(j)}=\bm{0}$; let $t=0$ and $\bm{\mathcal{E}}=\bm{0}$.
\Repeat
\State $t\leftarrow t+1$.
	  \For{$j=1,2,\ldots,d$}
	\State Compute $\bm{X}_t^{(j)}$ by \eqref{eq_trpca_x}.
	\State Compute $\bm{Y}_t^{(j)}$ by \eqref{eq_trpca_y}.
	\State Update $\bm{Z}^{(j)}$ by \eqref{eq_trpca_z}.
      \EndFor
\State Compute $\bm{\mathcal{E}}_{t}$ by \eqref{eq_soft_thresh}.
\Until{converged or $t=t_\mathrm{max}$}
\Ensure $\hat{\bm{\mathcal{T}}}=\bm{\mathcal{I}}\times_{1}\bm{X}_t^{(1)}\times_{2}\bm{X}_t^{(2)}\ldots\times_{d}\bm{X}_t^{(d)}$.
\end{algorithmic}
\end{algorithm}

In \eqref{problem_RTPCA_CVEN}, when $p\notin\lbrace 1/d, 2/d\rbrace$, the optimization becomes more difficult because all $d$ groups of the regularizers are nonconvex and nonsmooth. Thanks to Theorem \ref{theorem_sp_asym}, especially its (b), we can obtain arbitrarily sharp Schatten-$p$ quasi-norm regularization by using only one group of nonconvex and nonsmooth regularizers on the component vectors. The corresponding problem is
\begin{equation}\label{problem_RTPCA_CVEN_asym}
\begin{aligned}
\mathop{\text{minimize}}_{\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]},\ \bm{\mathcal{E}}}\ &\dfrac{1}{2}\left\Vert \bm{\mathcal{D}}-\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)}-\bm{\mathcal{E}}\right\Vert_F^2\\
&+\lambda_x\sum_{i=1}^k\Big(\dfrac{1}{q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\dfrac{1}{2}\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert^2\Big)+\lambda_e\Vert\bm{\mathcal{E}}\Vert_1.
\end{aligned}
\end{equation}
We propose to solve \eqref{problem_RTPCA_CVEN_asym} by ADMM with iteratively reweighted update (Algorithm \ref{alg_q_IRLS}) embedded, which is shown in Algorithm \ref{alg_RTPCA_asym}. In the algorithm, for $j=2,\ldots,d$, the subproblem of $\bm{X}_t^{(j)}$ has a closed-form solution, which makes the algorithm more efficient than Algorithm \ref{alg_RTPCA}, especially when $d$ is large.

\begin{algorithm}[h]
\caption{RTPCA-ENR when $p\notin\lbrace 1/d, 2/d\rbrace$}
\label{alg_RTPCA_asym}
\begin{algorithmic}[1]
\Require
$\bm{\mathcal{D}}$, $k$, $q$, $\lambda_x$, $\lambda_e$, $\mu$, $t_\mathrm{max}$.
\State Initialize $\lbrace\bm{x}_{i}^{(j)}\rbrace_{i\in[k]}^{j\in[d]}$ with \eqref{eq_x_ini};  let $t=0$ and $\bm{\mathcal{E}}=\bm{0}$; let $\bm{Y}^{(1)}_0=\bm{X}^{(1)}_0$ and $\bm{Z}^{(1)}=\bm{0}$.
\Repeat
\State $t\leftarrow t+1$.
	\State Compute $\bm{X}^{(1)}_t$ by \eqref{eq_trpca_x}.
	\State Compute $\bm{Y}^{(1)}_t$ by Algorithm \ref{alg_q_IRLS}.
	\State Update $\bm{Z}^{(1)}$ by \eqref{eq_trpca_z}.
	  \For{$j=2,3,\ldots,d$}
	\State Compute $\bm{X}^{(j)}_t$ by $\bm{X}^{(j)}_t=(\bm{D}_{(j)}-\bm{E}_{(j)})[(\bm{X}^{(i)})^{\odot_{i\neq j}}]\Big([(\bm{X}^{(i)})^{\odot_{i\neq j}}]^T[(\bm{X}^{(i)})^{\odot_{i\neq j}}]+\lambda_x\bm{I}\Big)^{-1}$.
      \EndFor
\State Compute $\bm{\mathcal{E}}_{t}$ by \eqref{eq_soft_thresh}.
\Until{converged or $t=t_\mathrm{max}$}
\Ensure $\hat{\bm{\mathcal{T}}}=\bm{\mathcal{I}}\times_{1}\bm{X}_t^{(1)}\times_{2}\bm{X}_t^{(2)}\ldots\times_{d}\bm{X}_t^{(d)}$.
\end{algorithmic}
\end{algorithm}

\section{Generalization Error Bound of LRTC}\label{sec:bound}
\subsection{Bound for the decomposition model}\label{sec:bound_CP}
Generalization error (or excess risk) measures the prediction capacity on new data of the machine learning model obtained from the training data. For low-rank matrix or tensor completion, the training data are the observed entries, and the generalization error is for the missing entries \cite{srebro2005rank,wangxu2012}. Usually, we expect that the generalization error is not too larger than the training error. So a good matrix or tensor completion method should have a tight generalization error bound. Foster and Risteski \cite{foster2019sum} showed that, in noisy completion of a rank-$r$ $n\times n\times n$ orthogonal tensor using sum-of-squares relaxation of tensor nuclear norm, the error bound is $\tilde{O}(r^2n^{3/2}/\vert\Omega\vert)$. In this paper, we show tensor completion error bounds for CP decomposition without using the orthogonal assumption.

First, we define the following tensor set:
\begin{equation*}
\begin{aligned}
\mathcal{S}_{dnk}=\lbrace  &\bm{\mathcal{X}}\in\mathbb{R}^{n\times \ldots\times n}: \bm{\mathcal{X}}=\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)};\\
&\ \Vert\bm{x}_i^{(j)}\Vert\leq\delta_{ij},\ 1\leq i\leq k,\ 1\leq j\leq d\rbrace.
\end{aligned}
\end{equation*}
Without loss of generality, we let the size of each mode be the same.
We consider the observation model \eqref{eq_LRTC_ob}, i.e. observing $\vert\Omega\vert$ entries of $\bm{\mathcal{D}}=\bm{\mathcal{T}}+\bm{\mathcal{N}}$ uniformly at random. The following theorem shows that deviation of the sampling operator $\mathcal{P}_\Omega$ from its expectation is controlled by $k/\vert\Omega\vert$. Since the rank of $\bm{\mathcal{X}}$ is upper-bounded by $k$, the bound in Theorem \ref{theorem_GE0} can be tighter when the rank of $\bm{\mathcal{X}}$ is smaller.
\begin{theorem}\label{theorem_GE0}
Suppose $\bm{\mathcal{X}}\in\mathcal{S}_{dnk}$ and $\max\lbrace{\Vert \bm{\mathcal{D}}\Vert_{\infty},\Vert \bm{\mathcal{X}}\Vert_{\infty}\rbrace}\leq\gamma$. Then with probability at least $1-2n^{-d}$, there exists a constant $c$ such that
\begin{align*}
&\sup_{{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\left\vert\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{D}}-\bm{\mathcal{X}}\Vert_F-\dfrac{1}{\sqrt{\vert\Omega\vert}}\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bm{\mathcal{X}})\Vert_F\right\vert\\
\leq& c\gamma\left(\dfrac{kdn\log \left(\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)\right)}{\vert\Omega\vert}\right)^{1/4}.
\end{align*}
\end{theorem}
For noisy tensor completion, we have the following generalization error bound.
\begin{theorem}\label{theorem_GE1}
Suppose $\bm{\mathcal{X}}\in\mathcal{S}_{dnk}$ and $\max\lbrace{\Vert \bm{\mathcal{D}}\Vert_{\infty},\Vert \bm{\mathcal{X}}\Vert_{\infty}\rbrace}\leq\gamma$. Then with probability at least $1-2n^{-d}$, there exists a constant $c$ such that
\begin{align*}
\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{T}}-\bm{\mathcal{X}}\Vert_F\leq &\dfrac{1}{\sqrt{\vert\Omega\vert}}\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bm{\mathcal{X}})\Vert_F+\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{N}}\Vert_F\\
&+c\gamma\left(\dfrac{kdn\log \left(\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)\right)}{\vert\Omega\vert}\right)^{1/4}.
\end{align*}
\end{theorem}
Note that in the theorem, $\bm{\mathcal{T}}$ is not necessarily in $\mathcal{S}_{dnk}$ and in fact can be any tensor. Hence, the error bound is useful even when $\bm{\mathcal{T}}$ is approximately low-rank (not exactly low-rank). According to Theorem \ref{theorem_GE1}, smaller $k$, larger $\vert\Omega\vert$, and smaller Euclidean norms of the component vectors provide a tighter error bound for tensor completion. Therefore, if the initial rank is large and there is no regularization, the tensor completion error can be large in practice. Notice that $kdn$ is exactly the maximum number of  parameters required to determine $\bm{\mathcal{X}}$ uniquely. When $k$ is the rank of $\bm{\mathcal{X}}$, $kdn$ is the number of degrees of freedom in $\bm{\mathcal{X}}$. To reduce the error of tensor completion, we need to reduce $\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bm{\mathcal{X}})\Vert_F$, $k$, and the Euclidean norms of the component vectors, which is exactly what our LRTC-ENR aims at.

In Theorem \ref{theorem_GE1}, when $\bm{\mathcal{N}}=\bm{0}$, $d=3$, and $k=r$, the error rate becomes $\tilde{O}((rn)^{1/4}/\vert\Omega\vert^{1/4})$. It cannot be directly compared with the error rate $\tilde{O}(r^2n^{3/2}/\vert\Omega\vert)$ given by Foster and Risteski \cite{foster2019sum} because the latter is based on the squared error. Using the inequality $\vert a-b\vert\leq\sqrt{\vert a^2-b^2\vert}$, we convert their bound to $\tilde{O}(rn^{3/4}/\vert\Omega\vert^{1/2})$, which can then be compared with our rate. When $\vert\Omega\vert<r^3n^2$ (much larger than the number of parameters in the CP decomposition of the tensor), our bound is tighter than theirs. Suppose $\vert\Omega\vert=rn^2$, our rate is $\tilde{O}(n^{-1/4})$ while their rate is $\tilde{O}(r^{1/2}n^{-1/4})$. Moreover, in contrast to Foster and Risteski \cite{foster2019sum}, we do not need the orthogonal assumption in CP decomposition, which is not easy to hold in practice.

\subsection{Bound for the Schatten-$p$ quasi-norm model}
The error bounds presented in Section \ref{sec:bound_CP} are not related to the tensor nuclear norm or Schatten-$p$ quasi-norm. In addition,  the tensor completion error bounds proposed in \cite{barak2016noisy} and \cite{foster2019sum} are for the sum-of-squares hierarchy of tensor nuclear norm. We aim to provide a generalization bound for Schatten-$p$ quasi-norm.
Considering Problem \eqref{problem_LRTC_Sp}, we expect that the tensor completion error when $p<1$ is less than that when $p=1$. In practical applications, we want to find a $p$ to make the tensor completion error as small as possible.

The following theorem provides an upper bound of the average squared errors for the observed entries in $\bm{\mathcal{T}}$.

\begin{theorem}\label{theorem_errorbound_sp}
Suppose $\hat{\bm{\mathcal{T}}}$ is the solution\footnote{The choice of $\lambda$ can be found in the proof.} of \eqref{problem_LRTC_Sp}. Let $c_{\kappa}=(2\kappa-1)(2\kappa)\kappa^{-1/(2\kappa-1)}$ where $\kappa=(2-p)/(2-2p)$. Let $\tilde{n}=\max\lbrace n,k \rbrace$. Then: \\
(a) with probability at least $1-2n^{-d}$, we have
\begin{align*}
&\dfrac{1}{\vert \Omega\vert}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber\\
\leq &\dfrac{8\sqrt{2}\sigma}{\vert \Omega\vert}\sqrt{dn\log(5d)+d\log(n)}\Vert {\bm{\mathcal{T}}}\Vert_\ast;
\end{align*}
(b) for any $0<p<1$, there exists some bounded constant $C=C(p,\vert\Omega\vert,\sigma)>0$ independent of $\tilde{n}$ such that with probability at least $1-C\exp(-d \tilde{n})$, we have
\begin{align*}
&\dfrac{1}{\vert \Omega\vert}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber\\
\leq &16c_{\kappa}(C^2/p)^{1-p/2}(d\tilde{n}/\vert\Omega\vert)^{1-p/2}\Vert\bm{\mathcal{T}}\Vert_{S_p}^p.
\end{align*}
Additionally, $C<p$ provided that $\vert\Omega\vert$ is sufficiently large and $p$ is sufficiently far from $0$ and $1$.
\end{theorem}
Notice that in Theorem \ref{theorem_errorbound_sp}(b), when $p\rightarrow 1$, with high probability, the right-hand side of the inequality approaches $c'\sigma\vert\Omega\vert^{-1}\sqrt{d\tilde{n}}\Vert\bm{\mathcal{T}}\Vert_\ast$ with some constant $c'$, which is comparable to Theorem \ref{theorem_errorbound_sp}(a). In Theorem \ref{theorem_errorbound_sp}(b): 1) $(C^2/p)^{1-p/2}$ decreases with descending $p$ provided that $p$ is not too small and $\vert\Omega\vert$ is large enough; 2)$(d\tilde{n}/\vert\Omega\vert)^{1-p/2}$ decreases with descending $p$ provided that  $\vert\Omega\vert>d\tilde{n}$; 3) $c_{\kappa}$ increases slowly when $p$ decreases but stays far from $0$; when $p$ decreases, $\Vert\bm{\mathcal{T}}\Vert_{S_p}^p$ increases if the singular values are relatively small and decreases otherwise. Putting these pieces together, we deduce that a smaller $p$ but not too small leads to a tighter generalization error bound provided that $\vert\Omega\vert$ is sufficiently large. Though Theorem \ref{theorem_errorbound_sp} relies on the optimal solution of \eqref{problem_LRTC_Sp}, it reveals the effect of $p$ on the tensor recovery error theoretically.

It is worth mentioning that in Theorem \ref{theorem_errorbound_sp}, the bound is for $\mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})$ and hence measures the recovery error for the entries in $\Omega$, while the error for the entries not in $\Omega$ remains unknown. We expect that $\tfrac{1}{n^d}\Vert\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F^2\simeq\tfrac{1}{\Omega}\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2$. The following theorem shows that the expected recovery error is small if $\vert\Omega\vert$ is large, $k$ is small, or/and $p$ (to some extent) is small.
\begin{theorem}\label{theorem_bound_overall}
Suppose $\hat{\bm{\mathcal{T}}}\in\mathcal{S}_{dnk}$ is the solution of \eqref{problem_LRTC_Sp}. Follow the same notations used in Theorem \ref{theorem_errorbound_sp}.
Suppose $\max\lbrace{\Vert \bm{\mathcal{T}}\Vert_{\infty},\Vert \hat{\bm{\mathcal{T}}}\Vert_{\infty}\rbrace}\leq\gamma$. Then with probability at least $1-4n^{-d}$, for any $0<p\leq1$,  there exists a constant $c$ and $C$ (stated in Theorem \ref{theorem_errorbound_sp}) such that
\begin{align*}
&\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F\\
\leq&4\beta_1\sigma^{1/2}\left((dn\log(5d)+d\log(n))/\vert\Omega\vert^2\right)^{1/4}(\Vert {\bm{\mathcal{T}}}\Vert_\ast)^{1/2}\\
&+4\beta_2c_{\kappa}^{1/2}(C^2/p)^{1/2-p/4}(d\tilde{n}/\vert\Omega\vert)^{1/2-p/4}(\Vert\bm{\mathcal{T}}\Vert_{S_p}^p)^{1/2}\\
&+c\gamma\left(\dfrac{kdn\log \left(\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)\right)}{\vert\Omega\vert}\right)^{1/4},
\end{align*}
where $\beta_1=1$ and $\beta_2=0$ when $p=1$ and $\beta_1=0$ and $\beta_2=1$ when $0<p<1$.
\end{theorem}

The right-hand-side of the bound in Theorem \ref{theorem_bound_overall} includes three parts. The first two parts involve tensor nuclear norm and Schatten-$p$ quasi-norm while the last part is related to the decomposition size of the tensor. When $\vert \Omega \vert$ increases, the three parts become smaller and hence the tensor completion error can be smaller. In our LRTC-ENR, we try to use a small $p$ and the algorithms explicitly reduces $k$ and $\delta_{ij}$. Therefore, our LRTC-ENR is able to provide high accuracy especially when a small $p$ is used.

\section{Recovery guarantee for TRPCA}\label{sec:bound_TRPCA}
The following theorem shows the upper bound of the recovery error for $\bm{\mathcal{T}}$ and $\bm{\mathcal{E}}$ by solving the Schatten-$p$ (quasi) norm regularized TRPCA problem \eqref{problem_RTPCA_Sp}.
\begin{theorem}\label{theorem_TRPCA}
Suppose the $\bm{\mathcal{E}}$ in \eqref{eq.TRPCA_model} has $s_e$ nonzero elements drawn from $\mathcal{N}(0,\varsigma^2)$. Let $c_N=2\Vert\bm{\mathcal{N}}\Vert_\infty$ and $c_T=\Vert\bm{\mathcal{T}}\Vert_\infty$. Let $c_{\kappa}=(2\kappa-1)(2\kappa)\kappa^{-1/(2\kappa-1)}$ where $\kappa=(2-p)/(2-2p)$. Let $\tilde{n}=\max\lbrace n,k \rbrace$. Suppose $\hat{\bm{\mathcal{T}}}$ and $\hat{\bm{\mathcal{E}}}$ are the solution\footnote{The choice of $\lambda_x$ and $\lambda_e$ can be found in the proof.} of \eqref{problem_RTPCA_Sp}. Then: \\
(a) with probability at least $1-4n^{-d}$, we have
\begin{align}
&\dfrac{1}{2}\Vert \bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F^2+\dfrac{1}{2}\Vert {\bm{\mathcal{E}}}-\hat{\bm{\mathcal{E}}}\Vert_F^2 \nonumber \\
\leq &4\sqrt{2}(\sigma+\varsigma)\sqrt{dn\log(5d)+d\log(n)}\Vert\bm{\mathcal{T}}\Vert_\ast\nonumber\\
&+(c_T+c_N)\Vert\bm{\mathcal{E}}\Vert_1+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty\Vert\hat{\bm{\mathcal{E}}}\Vert_1;\nonumber
\end{align}
(b) for any $0<p<1$, there exist some bounded constants $C_1=C_1(p,n^d,\sigma)>0$, $C_2=C_2(p,s_e,\varsigma)>0$ independent of $\tilde{n}$ such that with probability at least $1-2\max(C_1,C_2)\exp(-\tilde{n}/\max(C_1,C_2)^2)$, we have
\begin{align}
&\dfrac{1}{2}\Vert \bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F^2+\dfrac{1}{2}\Vert {\bm{\mathcal{E}}}-\hat{\bm{\mathcal{E}}}\Vert_F^2 \nonumber \\
\leq &8c_{\kappa}\tilde{n}\big((\tfrac{C_1^2}{p})^{1-\frac{p}{2}}(\tfrac{n^d}{\tilde{n}})^{\frac{p}{2}}+(\tfrac{C_2^2}{p})^{1-\frac{p}{2}}(\tfrac{s_e}{\tilde{n}})^{\frac{p}{2}}\big)\Vert {\bm{\mathcal{T}}}\Vert_{S_p}^p \nonumber \\
&+(c_T+c_N)\Vert\bm{\mathcal{E}}\Vert_1+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty\Vert\hat{\bm{\mathcal{E}}}\Vert_1.\nonumber
\end{align}
Additionally, $C_1,C_2<p$ provided that $p$ stays far from $0$ and $1$.
\end{theorem}

In Theorem \ref{theorem_TRPCA}, the lower bounds of $\lambda_x$ and $\lambda_e$ are in the proof (Appendix \ref{appendix_the_TRPCA}). In Theorem \ref{theorem_TRPCA} (a) and (b), the bounds become tighter when $\sigma$ is smaller. In real applications, usually, $\bm{\mathcal{N}}$ is small Gaussian noise and $\bm{\mathcal{E}}$ is sufficiently sparse but $\varsigma$ can be relatively large. Given $\bm{\mathcal{T}}$ and $\bm{\mathcal{E}}$,  the bounds are dominated by $p$, $\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty$, and $\Vert\hat{\bm{\mathcal{E}}}\Vert_1$, where $\Vert\hat{\bm{\mathcal{E}}}\Vert_1$ has appeared in the objective function. We may add one more constraint $\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty\leq C_T$ to ensure the bounds tight, where $C_T$ is some pre-defined constant.

In Theorem \ref{theorem_TRPCA} (b), $(\tfrac{n^d}{\tilde{n}})^{\frac{p}{2}}$ and $(\tfrac{s_e}{\tilde{n}})^{\frac{p}{2}}$ decrease quickly when $p$ becomes smaller because $n^d$ and $s_e$ are much larger than $\tilde{n}$. Using the definition of $c_\kappa$, $C_1$, and $C_2$, we conclude that a smaller $p$ but not too small provides tighter error bound of TRPCA. Thus, we may use a relatively small $p$ in \eqref{eq.TRPCA_model} or \eqref{problem_RTPCA_CVEN} to achieve high recovery accuracy for $\bm{\mathcal{T}}$.

\section{Experiments of LRTC}\label{sec_exp_LRTC}

We compare our tensor completion method LRTC-ENR with HaLRTC \cite{liu2012tensor}, TenALS \cite{{jain2014provable}}, TMac \cite{xu2013parallel}, BCPF \cite{zhao2015bayesian}, Rprecon \cite{kasai2016low}, KBR-TC \cite{8000407}, TRLRF \cite{yuan2019tensor}. In TenALS, TMac, BCPF, Rprencon, TRLRF, and our LRTC-ENR, we need to determine the factorization size (or the initial rank in other word) beforehand. These methods, compared to HaLRTC and KBR-TC, may be applicable to big tensors efficiently provided that the factorization sizes are sufficiently small. Particularly, in TMac, BCPF, and our LRTC-ENR, the rank is adjusted adaptively.
The tensor operation and computation are based on the MATLAB Tensor Tool Box  \cite{TTB_Software}.
\subsection{Synthetic data}
\subsubsection{Data generating model}\label{data_gen_syn}
We generate synthetic tensors by $\bm{\mathcal{D}}=\bm{\mathcal{T}}+\bm{\mathcal{N}}$,
where $\bm{\mathcal{T}}=\sum_{i=1}^r\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\circ\bm{x}_i^{(3)}$.
The entries of $\bm{x}_i^{(j)}\in\mathbb{R}^{n}$ ($i\in[r]$, $j\in[3]$) are drawn from $\mathcal{N}(0,1)$. The entries of the noise tensor $\bm{\mathcal{N}}$ are drawn from $\mathcal{N}(0,(0.1\sigma_{\bm{\mathcal{T}}})^2)$, where $\sigma_{\bm{\mathcal{T}}}$ denotes the standard deviation of the entries of $\bm{\mathcal{T}}\in\mathbb{R}^{n\times n\times n}$. We set $n=50$  and randomly remove a fraction (which we call missing rate) of the entries to test the performance of tensor completion, measured by
$$\text{Relative recovery error}=\Vert \mathcal{P}_{\bar{\Omega}}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F/\Vert \mathcal{P}_{\bar{\Omega}}(\bm{\mathcal{T}})\Vert_F,$$
where $\bar{\Omega}$ denotes the locations of unknown entries and $\hat{\bm{\mathcal{T}}}$ denotes the recovered tensor.

\subsubsection{LRTC with different regularizers}\label{sec_compare_p_syn}
We compare the performance of LRTC with Schatten-$p$ (quasi) norm regularization \eqref{problem_LRTC_Sp_min} and Euclidean norm regularization (solved by BCDE and LBFGS) in Figure \ref{fig_p}, in which we set $r=10$ and considered different noise levels, missing rates, and $p$ values. In all methods, we set $k=2r=20$, $t_\mathrm{max}=500$, and use grid search (within $[0.01,500]$) to determine $\lambda$. We have the following results.
\begin{itemize}
\item\textit{Regularization is helpful}
\begin{itemize}
\item In Figure \ref{fig_p} (a-d), the recovery error when there is no regularization (i.e. $\lambda=0$) is higher than those of the regularized methods.
\end{itemize}
\item\textit{Smaller $p$ leads to lower recovery error?}
\begin{itemize}
\item In almost all cases, $p<1$ outperforms $p=1$, verifying the superiority of Schatten-$p$ quasi-norm over nuclear norm in LRTC (both LRTC-ENR and LRTC-Schatten-$p$).
\item In LRTC-ENR solved by BCDE and LBFGS, smaller $p$ provides lower recovery error provided that $p$ is not too small (e.g. $p\geq 1/6$), which matches Theorem \ref{theorem_errorbound_sp}.
\end{itemize}
\item\textit{LRTC-ENR v.s. LRTC-Schatten-$p$}
\begin{itemize}
\item Compared to Euclidean norm minimization, direct Schatten-$p$ norm minimization has higher recovery error and time cost because the optimization problem is more difficult to solve.
\end{itemize}
\item\textit{BCDE v.s. LBFGS}
\begin{itemize}
\item When $p$ is too small, the recovery error of LRTC-ENR solved by BCDE does not change but the recovery error of LRTC-ENR solved by LBFGS increases. One possible reason is that LBFGS is not effective in handling nonsmooth objective function and is often stuck in bad local minima or even saddle points, especially when $p$ is very small. Notice that when $p<1/3$, LRTC-ENR has at least one nonsmooth nonconvex subproblem, which brings difficulty to LBFGS.
\item LBFGS is more efficient than BCDE because the latter is in the manner alternating update scheme. In addition, LRTC-ENR solved by BCDE requires to determine $\hat{L}_{j,t}$, for which the estimation given by \eqref{eq_estimateL} may not be a very good choice especially when the missing rate is high.
\end{itemize}
\item\textit{Symmetric regularizer v.s. asymmetric regularizer}
\begin{itemize}
\item For LRTC-ENR, the asymmetric regularization slightly outperforms the symmetric regularization in a few cases because the former has only one nonsmooth subproblem.
\end{itemize}
\end{itemize}
Based on the above results, we suggest using BCDE to solve LRTC-ENR on small datasets and using LBFGS to solve LRTC-ENR on large datasets. In Schatten-$p$ quasi norm, we suggest using $p=1/3$ or $p=1/6$.

\begin{figure}[h!]
\includegraphics[width=8cm,trim={20 10 40 20},clip]{fig_syn_p}
\caption{Performance (average of 20 repeated trials) of Schatten-$p$ (quasi) norm with different $p$ in the case of different noise level $v$ for the noise distribution $\mathcal{N}(0,(v\sigma_{\bm{\mathcal{T}}})^2)$; ``NoReg'' means LRTC without regularization \cite{acar2011scalable}; the time costs (500 iterations) of the five methods are 8.6s, 3.0s, 2.1s, 2.8s, and 2.3s respectively).}\label{fig_p}
\end{figure}

\subsubsection{LRTC-ENR v.s. seven baseline methods}
Figure \ref{fig_syn_all_r10} shows the recovery of LRTC-ENR ($p=1/3$, solved by LBFGS) and seven baselines on the synthetic data with $r=10$. In all methods except HaLRTC and KBR-TC, one has to determine the initial rank beforehand. We have set the initial rank to $2r$ because in practice it is difficult to know the true rank. Though $r$ is very small compared to the size of the tensor, the recovery errors of HaLRTC and Rprecon increase quickly when the missing rate is high. TMac, BCPF, and LRTC-ENR are able to adjust rank adaptively and the latter two perform much better than TMac.

\begin{figure}[h!]
\centering
\includegraphics[width=7cm,trim={0 0 20 20},clip]{fig_syn_ALL_r5}
\caption{Low-rank tensor completion on synthetic data ($r=10$): LRTC-ENR ($p=1/3$) v.s. seven baseline methods}\label{fig_syn_all_r10}
\end{figure}

In Figure \ref{fig_syn_all_r10}, when the missing rate is $0.9$, we observed $12,500$ entries, much more than the minimum number (1500) of freedom parameters required to determine the tensor uniquely. It means the tensor completion problem is too easy. Hence we increase $r$ to 50 and report the recovery error and computational time of KBR-TC, TenALS, BCPF, TMac, and LRTC-ENR in Figure \ref{fig_syn_all_r50}. We see that the recovery error of BCPF and LRTC-ENR are lower than those of other methods. When the missing rate is sufficiently high, LRTC-ENR has smaller recovery error than BCPF. Moreover, LRTC-ENR (200 iterations) is at least 15 times faster than BCPF and 5 times faster than KBR-TC.

\begin{figure}[h!]
\centering
\includegraphics[width=8cm,trim={35 0 50 0},clip]{fig_syn_ALL_r50}
\caption{Low-rank tensor completion on synthetic data: LRTC-ENR ($p=1/3$) v.s. baseline methods, $r=50$; (a) recovery error with different missing rates; (b) average time cost in all cases of missing rates.}\label{fig_syn_all_r50}
\end{figure}

\begin{figure*}[h!]
\centering
\includegraphics[width=14cm,trim={15 120 50 30},clip]{fig_MSI_im1_99}
\caption{Inpainting performance on the first image (band 16) of the MSI dataset when the missing rate is 0.99.}\label{fig_real_msi_1_99}
\end{figure*}

\begin{figure*}[h!]
\centering
\includegraphics[width=14cm,trim={15 120 50 30},clip]{fig_MSI_im1_97}
\caption{Inpainting performance on the first image (band 16)  of the MSI dataset when the missing rate is 0.97.}\label{fig_real_msi_1_97}
\end{figure*}

\subsection{Multi-spectral image inpainting}
We use the Columbia multi-spectral image database \cite{CAVE_0293} to show the effectiveness of our tensor completion method in image inpainting problem. The dataset consists of the multi-spectral images of 32 real-world scenes of a variety of real-world materials and objects. The images have a spatial resolution of $512\times 512$ and the number of spectral bands as $31$. Thus the data of each image is a tensor of size $512\times 512\times 31$. In our experiments, we pre-scale all images of every band to $[0,1]$.  Since TRLRF \cite{yuan2019tensor} and Rprecon \cite{kasai2016low} have extremely high computational costs on these tensors, we omit their implementations here. In \cite{8000407}, it has been shown that BCPF \cite{zhao2015bayesian} significantly outperformed HaLRTC \cite{liu2012tensor} as well as a few other baselines on this dataset, so we will not compare HaLRTC. The initial rank of TMac \cite{xu2013parallel} and our LRTC-ENR ($\lambda=0.1$) are set to $100$. When the initial rank is $100$, it takes TenALS \cite{jain2014provable} and BCPF \cite{zhao2015bayesian} more than $1500$ seconds on each image, thus we set the initial rank to $50$. Since the tensor rank is very small compared to the size and the images are noiseless, we here consider highly incomplete images of missing rate no less than 0.95.

As an example, Figure \ref{fig_real_msi_p} shows the recovery errors of LRTC-ENR with different $p$ and the compared methods on the first image  of the dataset when the missing rate is 0.95, 0.97, and 0.99. It can be found that when the missing rate is 0.95, KBR-TC outperformed our LRTC-ENR. When the missing rate is 0.99, our LRTC-ENR outperformed all other methods. In Figure \ref{fig_real_msi_p} (a), for LRTC-ENR, $p=1/6$ is the best choice, while in Figure \ref{fig_real_msi_p} (b) and (c), $p=1/3$ has least recovery error. These results are consistent with Theorem \ref{theorem_errorbound_sp}.

\begin{figure}[h!]
\includegraphics[width=8.5cm,trim={50 0 50 0},clip]{fig_MSI_p}
\caption{Recovery error on the first image of the MSI dataset.}\label{fig_real_msi_p}
\end{figure}

The average results (PSNR) on the 32 images are reported in Table \ref{tab_CAVE}. Obviously, KBR-TC and LRTC-ENR ($p=1/3$) outperformed other methods significantly. When the missing rate is 0.95, KBR-TC outperformed LRTC-ENR. When the missing rate is 0.97 or 0.99, LRTC-ENR outperformed KBR-TC, and the improvement is significant according to the $p$-value of paired t-test. In addition, the time cost of LRTC-ENR is only $15\%$ of that of KBR-TC. Figure \ref{fig_real_msi_1_97} and Figure  \ref{fig_real_msi_1_99} present the recovered images given by the compared methods on the first image of the dataset. Visually,  the recovery performance of our LRTC-ENR ($p=1/3$) is much better than those of other methods.

\begin{table}[h!]
\caption{Average recovery performance (PSNR) on the Columbia MSI dataset of 32 images (the $p$-value is from the paired t-test between KBR-TC and LRTC-ENR)}
\begin{tabular}{l|ccc||c}
\hline
Missing rate& \multicolumn{1}{c}{0.95}  &  \multicolumn{1}{c}{0.97} & \multicolumn{1}{c||}{0.99} & Time cost \\ \hline
TenALS \cite{jain2014provable} &30.1$\pm$4.4 &25.9$\pm$3.22 &13.3$\pm$4.9 & 500s\\
BCPF \cite{zhao2015bayesian}	&30.4$\pm$5.8 &26.8$\pm$5.6  &20.2$\pm$3.3 & 450s \\
TMac \cite{xu2013parallel}	&33.3$\pm$5.4 &30.6$\pm$4.9 &17.9$\pm$3.1 & 190s\\
KBR-TC \cite{8000407} &\textbf{37.1}$\pm$4.9 &31.4$\pm$4.9  &20.5$\pm$2.8 & 730s \\
LRTC-ENR &35.4$\pm$4.4 &\textbf{32.4}$\pm$5.1 & \textbf{26.1}$\pm$5.4 & \textbf{110}s\\ \hline \hline
$p$-value {\scriptsize (t-test)} &5$\times 10^{-11}$ &5$\times 10^{-4}$ & 2$\times 10^{-8}$& 0 \\ \hline
\end{tabular}\label{tab_CAVE}
\vspace{-10pt}
\end{table}

\begin{figure}
\includegraphics[width=8cm,trim={40 0 50 5},clip]{fig_TRPCA_syn}
\caption{Performance of TRPCA on synthetic data (average of 20 repeated trials): (a) KBR-RPCA v.s. TRPCA-ENR (initial rank=2r); (b)  relative recovery error of TRPCA-ENR with different initial rank (noise density=0.5).}\label{fig_TRPCA_syn}
\end{figure}

\section{Experiments of Tensor Robust  PCA}\label{sec_exp_TRPCA}
\subsection{Synthetic data}\label{sec_trpca_syn}
We generate synthetic tensors by $\bm{\mathcal{D}}=\bm{\mathcal{T}}+\bm{\mathcal{N}}+\bm{\mathcal{E}}$.
The low-rank tensor $\bm{\mathcal{T}}$ is given by $\bm{\mathcal{T}}=\sum_{i=1}^rw_i\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\circ\bm{x}_i^{(3)}$, where $w_i=i/r$ and the entries of $\bm{x}_i^{(j)}\in\mathbb{R}^{n}$ ($i\in[r]$, $j\in[3]$) are drawn from $\mathcal{N}(0,1)$. $\bm{\mathcal{N}}$ is a dense noise tensor drawn from $\mathcal{N}(0,(0.1\sigma_{\bm{\mathcal{T}}})^2)$ and $\bm{\mathcal{E}}$ is a sparse noise tensor drawn from $\mathcal{N}(0,\sigma_{\bm{\mathcal{T}}}^2)$.  We set $n=50$, $r=25$ and evaluate the TRPCA performance in recovering $\bm{\mathcal{T}}$ in the case of different sparsity (noise density) of  $\bm{\mathcal{E}}$. The evaluation metric is the relative recovery error defined by
$$\text{Relative recovery error}=\Vert \bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F/\Vert \bm{\mathcal{T}}\Vert_F.$$
We compare our TRPCA-ENR with KBR-PCA \cite{8000407}. The TRPCA (denoted by T-TRPCA) proposed by \cite{lu2019tensor} and the robust tensor decomposition method OITNN-L proposed by \cite{wang2020robust} are under the assumption of $t$-product and hence not suitable in our setting of synthetic data. We will compare them on real data. The hyper-parameters of KBR-RPCA and our TRPCA-ENR\footnote{Throughout this paper, we set $\mu=10$ and $t_\mathrm{max}=500$ for Algorithms \ref{alg_RTPCA} and  \ref{alg_RTPCA_asym}.} are sufficiently tuned to provide their best performance. Figure \ref{fig_TRPCA_syn} (a) shows the recovery errors of KBR-RPCA and TRPCA-ENR (with $p=2/3$, $1/3$, and $1/6$) when the noise density increased from 0.1 to 0.7. We see that TRPCA-ENR outperformed KBR-RPCA consistently. In TRPCA-ENR, $p =1/6$ is better than $p=2/3$ and $1/3$ but the difference is not obvious. In Figure \ref{fig_TRPCA_syn} (b), we show the performance of different $p$'s. It can be found that a smaller $p$ (but not too small, e.g. larger than $10^{-2}$) yields less recovery error. This result is consistent with Theorem \ref{theorem_TRPCA}.

\begin{figure}
\centering
\includegraphics[width=8cm,trim={55 10 90 0},clip]{fig_9images}
\caption{Nine color images for denoising}\label{fig_9images}
\vspace{-10pt}
\end{figure}

\vspace{-10pt}
\subsection{Color image denoising}
We compare our TRPCA-ENR with KBR-PCA \cite{8000407}, T-TRPCA \cite{lu2019tensor}, and OITNN-L \cite{wang2020robust} in the task of color image denoising. Shown in Figure \ref{fig_9images}, nine color images of size $256\times 256\times 3$ used in \cite{wang2020robust} are considered in this study.
For each image, we randomly set $10\%$ of the tensor elements to random values in $[0, 1]$. We tuned all hyper-parameters carefully to provided the best denoising performance of the methods. As a result, in KBR-RPCA, we set $\lambda=2$ or $3$. In OITNN-L, we set $\lambda_L=4$ or $5$ and $\lambda_S=0.16$. In T-TRPCA,  we set $\lambda=1/\sqrt{768}$ or $2/\sqrt{768}$. In our TRPCA-ENR, we set the initial rank to a relative large value 200 because the approximate rank of the image tensors is not too small and our methods can adjust the rank adaptively; we set $\lambda_e=0.02$ and consider the cases of $p=2/3$, $1/3$, $1/6$, and $1/10$, for which $\lambda_x=0.1$, $0.1$, $0.02$, and $0.015$.

\begin{table*}[h]
\centering
\caption{PSNR of color image denoising}
\begin{tabular}{lccccccccc}
\hline

& Image 1 & Image 2		& Image 3		& Image 4		& Image5		& Image 6		& Image 7  & Image 8		& Image 9  \\ \hline
KBR-RPCA \cite{8000407} & 30.25$\pm$0.06	&	33.36$\pm$0.22&31.93$\pm$0.12 &31.87$\pm$0.08	&35.08$\pm$0.12	&31.36$\pm$0.09	&30.08$\pm$0.10 &32.41$\pm$0.05&34.28$\pm$0.12 \\
T-TRPCA \cite{lu2019tensor} & 31.16$\pm$0.05&32.75$\pm$0.06 &32.87$\pm$0.11 & 32.52$\pm$0.06&33.91$\pm$0.10 &28.96$\pm$0.06 & 29.36$\pm$0.09&31.49$\pm$0.11&33.67$\pm$0.14 \\
OITNN-O \cite{wang2020robust} & 28.46$\pm$0.04&29.99$\pm$0.03 &29.16$\pm$0.05 & 28.13$\pm$0.04&29.86$\pm$0.03 &27.39$\pm$0.05 & 27.46$\pm$0.05&28.61$\pm$0.07&30.06$\pm$0.06 \\
TRPCA-ENR ($p$=$\tfrac{2}{3}$)  & 34.23$\pm$0.12&35.44$\pm$0.09& 34.82$\pm$0.18	&34.23$\pm$0.14	&37.14$\pm$0.06	&35.09$\pm$0.12	& 33.67$\pm$0.10&35.66$\pm$0.16&36.12$\pm$0.05 \\
TRPCA-ENR ($p$=$\tfrac{1}{3}$)	 &34.68$\pm$0.18&36.88$\pm$0.16&37.52$\pm$0.12	&34.53$\pm$0.31	&37.78$\pm$0.04	&33.21$\pm$0.54	& 34.10$\pm$0.23&36.09$\pm$0.14&37.95$\pm$0.16\\

TRPCA-ENR ($p$=$\tfrac{1}{6}$)	 &34.96$\pm$0.40&36.04$\pm$1.08&37.57$\pm$0.33	&34.46$\pm$1.21&38.14$\pm$0.11&31.93$\pm$0.98	& 34.44$\pm$0.12&36.49$\pm$0.08&37.96$\pm$0.10\\

TRPCA-ENR ($p$=$\tfrac{1}{10}$)	 & \textbf{35.45}$\pm$0.12&\textbf{37.14}$\pm$0.10&\textbf{37.76}$\pm$0.14	&\textbf{35.38}$\pm$0.46	&\textbf{38.28}$\pm$0.13	&\textbf{35.18}$\pm$0.35	& \textbf{34.57}$\pm$0.14&\textbf{36.71}$\pm$0.06&\textbf{38.01}$\pm$0.16\\  \hline
\end{tabular}\label{tab_TRPCA_image}
\end{table*}

\begin{figure*}[h]
\centering
\includegraphics[width=15cm,trim={100 110 80 30},clip]{fig_TRPCA_im1}
\caption{Denoising performance on Image 1.}\label{fig_TRPCA_im1}
\end{figure*}

The average PSNR with standard deviation of 20 repeated trials are reported in Tables \ref{tab_TRPCA_image}. We see that the PSNRs of our TRPCA-ENR methods are much lower than those of other methods in all cases. In addition, TRPCA-ENR with smaller $p$ has higher PSNR and TRPCA-ENR with $p=1/10$ performs the best. Figure \ref{fig_TRPCA_im1} visualizes the denoising performance on Image 1, showing that TRPCA-ENR methods are better than other methods. These results provided evidence that tensor Schatten-$p$ quasi-norm  minimization is able to provide higher recovery accuracy in TRPCA when a relatively smaller $p$ is used.

\section{Conclusion}\label{sec:conclusion}
In this work, we focused on efficient and accurate low-rank tensor recovery and proposed a new class of nonconvex rank regularizers based on the Euclidean norms of component vectors in CP decomposition. We  have shown that these regularizers are monotonic transformations of tensor nuclear norm and Schatten-$p$ quasi-norms. Thus, the regularizers enable us to implicitly and efficiently minimize Schattern-$p$ quasi-norms serving as arbitrarily sharp rank proxies provided that $p$ is sufficiently small. Although based on CP decomposition, the regularizers are able to determine the rank adaptively and hence are not sensitive to the initial rank. We applied the regularizers to LRTC and TRPCA and provided efficient optimization algorithms. We proved that Schattern-$p$ quasi-norm with a smaller $p$ but not too small provides tighter error bounds for LRTC and TRPCA.

Numerical results on synthetic data and natural images demonstrated the advantages of our methods in LRTC and TRPCA against state-of-the-art methods. For instance, our LRTC-ENR is more accurate and efficient than BCPF \cite{zhao2015bayesian} and KBR-TC \cite{8000407} in multi-spectral image inpainting. The numerical results also showed that a smaller $p$ is potentially able to result in higher recovery accuracy in LRTC and TRPCA. All code used in this paper are available at xxxx.

Future work may focus on the following two points. First,  good initializations for the noncovex LRTC and TRPCA methods may improve the convergence speed and recovery accuracy. But unfortunately, we currently only have a random initialization for the tensor case, while the popular spectral initialization in matrix case cannot be used. Second, in our experiments, we observed that even with random initialization, LRTC-ENR and TRPCA-ENR often achieved optimal solutions. It is possible to extend the theoretical analysis on the global optimality of low-rank matrix factorization in \cite{haeffele2019structured} to our low-rank tensor decomposition models: if the number of the factors (i.e. $k$) is large enough, under certain mild conditions, any local minimizer for the factors is a global minimizer, although the optimization problem is nonconvex or/and nonsmooth.

\appendices
\section{Proof of Theorem \ref{theorem_sp}}
\begin{proof}
Let $\bm{x}_i^{(j)}=\alpha_{ij}\bar{\bm{x}}_i^{(j)}$, where $\Vert \bar{\bm{x}}_i^{(j)}\Vert=1$. Then $\prod_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert=\prod_{j=1}^d\alpha_{ij}=\vert\lambda_i\vert$. We have
\begin{equation*}
\begin{aligned}
\dfrac{1}{d}\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^q&\geq\sum_{i=1}^k\left(\prod_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^q\right)^{1/d}\\
&=\sum_{i=1}^k\prod_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^{q/d}=\sum_{i=1}^k\vert\lambda_i\vert^{q/d},
\end{aligned}
\end{equation*}
in which the inequality holds according to the AM-GM inequality. When $\alpha_{i1}=\alpha_{i2}=\ldots\alpha_{id}=\vert\lambda_i\vert^{1/d}$, the equality of AM-GM inequality is true. Replacing $q$ with $pd$, we have
$$\dfrac{1}{d}\sum_{i=1}^k\sum_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^{pd}\geq\sum_{i=1}^k\vert\lambda_i\vert^{p}\geq \Vert \bm{\mathcal{X}}\Vert_{S_p}^p.$$
\end{proof}

\section{Proof of Theorem \ref{theorem_sp_asym}}

\begin{proof}
(1) We have
\begin{align*}
\begin{aligned}
&\sum_{i=1}^k\Big(\dfrac{1}{q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert\Big)
=\sum_{i=1}^k\Big(\sum_{l=1}^{1/q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert\Big)\\
\geq&\sum_{i=1}^k(\dfrac{1}{q}+d-1)\left(\prod_{l=1}^{1/q}\Vert\bm{x}_i^{(1)}\Vert^{q}\prod_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert\right)^{\tfrac{q}{1+qd-q}}\\
=&\sum_{i=1}^k(\dfrac{1}{q}+d-1)\left(\prod_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert\right)^{\tfrac{q}{1+qd-q}}\\
=&\dfrac{1+qd-q}{q}\sum_{i=1}^k\vert\lambda_i\vert^{\tfrac{q}{1+qd-q}}
\geq\dfrac{1+qd-q}{q}\Vert \bm{\mathcal{X}}\Vert_{S_{q/(1+qd-q)}}^{q/(1+qd-q)}.
\end{aligned}
\end{align*}

(2) We have
\begin{align*}
\begin{aligned}
&\sum_{i=1}^k\Big(\dfrac{2}{q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert^2\Big)
=\sum_{i=1}^k\Big(\sum_{l=1}^{2/q}\Vert\bm{x}_i^{(1)}\Vert^{q}+\sum_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert^2\Big)\\
\geq&\sum_{i=1}^k(\dfrac{2}{q}+d-1)\left(\prod_{l=1}^{2/q}\Vert\bm{x}_i^{(1)}\Vert^{q}\prod_{j=2}^d\Vert\bm{x}_i^{(j)}\Vert^2\right)^{\tfrac{q}{2+qd-q}}\\
=&\sum_{i=1}^k(\dfrac{2}{q}+d-1)\left(\prod_{j=1}^d\Vert\bm{x}_i^{(j)}\Vert^2\right)^{\tfrac{q}{2+qd-q}}\\
=&\dfrac{2+qd-q}{q}\sum_{i=1}^k\vert\lambda_i\vert^{\tfrac{2q}{2+qd-q}}
\geq\dfrac{2+qd-q}{q}\Vert \bm{\mathcal{X}}\Vert_{S_{2q/(2+qd-q)}}^{2q/(2+qd-q)}.
\end{aligned}
\end{align*}
The equality in the second row of formula (1) (also (2)) holds when all terms in the parentheses are equal, for each $i$.
\end{proof}

\section{Proof of the last two in Tabel\ref{tab_reg3}}\label{proof_S25asym}
\begin{proof}
Recall the definition of $\lambda_1,\ldots,\lambda_k$.
We have
\begin{align*}
\begin{aligned}
&\dfrac{16^{1/5}}{5}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^{2}+\Vert\bm{x}_i^{(2)}\Vert+\Vert\bm{x}_i^{(3)}\Vert\right)\\
=&\dfrac{16^{1/5}}{5}\sum_{i=1}^k\bigg(\Vert\bm{x}_i^{(1)}\Vert^{2}+\dfrac{1}{2}\Vert\bm{x}_i^{(2)}\Vert+\dfrac{1}{2}\Vert\bm{x}_i^{(2)}\Vert\\
&\qquad \qquad +\dfrac{1}{2}\Vert\bm{x}_i^{(3)}\Vert+\dfrac{1}{2}\Vert\bm{x}_i^{(3)}\Vert\bigg)\\
\geq&\dfrac{16^{1/5}}{5}\sum_{i=1}^k5\left(\dfrac{1}{16}\Vert\bm{x}_i^{(1)}\Vert^{2}\Vert\bm{x}_i^{(2)}\Vert\Vert\bm{x}_i^{(2)}\Vert\Vert\bm{x}_i^{(3)}\Vert\Vert\bm{x}_i^{(3)}\Vert\right)^{1/5}\\
=&\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^{2}\Vert\bm{x}_i^{(2)}\Vert^{2}\Vert\bm{x}_i^{(3)}\Vert^2\right)^{1/5}\\
=&\sum_{i=1}^k\vert\lambda_i\vert^{2/5}
\geq\Vert \bm{\mathcal{X}}\Vert_{S_{2/5}}^{2/5}.
\end{aligned}
\end{align*}

Similarly, we have
\begin{align*}
\begin{aligned}
&\dfrac{81^{1/7}}{7}\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert^{3}+\Vert\bm{x}_i^{(2)}\Vert+\Vert\bm{x}_i^{(3)}\Vert\right)\\
=&\dfrac{81^{1/7}}{7}\sum_{i=1}^k\bigg(\Vert\bm{x}_i^{(1)}\Vert^{3}+\sum_{j=1}^3\Big(\dfrac{1}{3}\Vert\bm{x}_i^{(2)}\Vert+\dfrac{1}{3}\Vert\bm{x}_i^{(3)}\Vert\Big)\bigg)\\
\geq&\dfrac{81^{1/7}}{7}\sum_{i=1}^k7\left(\dfrac{1}{3^6}\Vert\bm{x}_i^{(1)}\Vert^3\Vert\bm{x}_i^{(2)}\Vert^3\Vert\bm{x}_i^{(3)}\Vert^3\right)^{1/7}\\
=&\sum_{i=1}^k\left(\Vert\bm{x}_i^{(1)}\Vert\Vert\bm{x}_i^{(2)}\Vert\Vert\bm{x}_i^{(3)}\Vert\right)^{2/7}\\
=&\sum_{i=1}^k\vert\lambda_i\vert^{3/7}
\geq\Vert \bm{\mathcal{X}}\Vert_{S_{3/7}}^{3/7}.
\end{aligned}
\end{align*}
\end{proof}

\section{Lipschitz constant estimation for $\bm{\nabla}\mathcal{L}$ in \eqref{eq_estimateL}}\label{appendix_Lips}
For convenience, we denote a perturbed copy of $\hat{\bm{X}}_{t-1}^{(j)}$ by $\hat{\bm{Z}}_{t-1}^{(j)}$. We have
\begin{align*}
\begin{aligned}
&\Vert\bm{\nabla}_{\hat{\bm{X}}_{t-1}^{(j)}}\mathcal{L}-\bm{\nabla}_{\hat{\bm{Z}}_{t-1}^{(j)}}\mathcal{L}\Vert_F\\
=&\Vert\left(\bm{M}_{(j)}\ast\big(\bm{D}_{(j)}-\hat{\bm{X}}_{t-1}^{(j)}[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]^\top\big)\right)\big(-[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\big)\\
&-\left(\bm{M}_{(j)}\ast\big(\bm{D}_{(j)}-\hat{\bm{Z}}_{t-1}^{(j)}[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]^\top\big)\right)\big(-[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\big)\Vert_F\\
\leq&\Vert[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\Vert_2\Vert\bm{M}_{(j)}\ast\big(\hat{\bm{X}}_{t-1}^{(j)}[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]^\top\\
&-\hat{\bm{Z}}_{t-1}^{(j)}[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]^\top\big)\Vert_F\\
\leq& \Vert[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\Vert_2\varrho\sqrt{\dfrac{\vert\Omega\vert}{\prod_{i=1}^dn_i}}\Vert\big(\hat{\bm{X}}_{t-1}^{(j)}-\hat{\bm{Z}}_{t-1}^{(j)}\big)[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]^\top\Vert_F\\
\leq &\varrho\sqrt{\dfrac{\vert\Omega\vert}{\prod_{i=1}^dn_i}}\Vert[(\bm{X}^{(i)})_t^{\odot_{i\neq j}}]\Vert_2^2\Vert\hat{\bm{X}}_{t-1}^{(j)}-\hat{\bm{Z}}_{t-1}^{(j)}\Vert_F.
\end{aligned}
\end{align*}
Here $0<\varrho\leq \sqrt{\dfrac{\prod_{i=1}^dn_i}{\vert\Omega\vert}}$ is some suitable constant.

\section{Proof of Theorem \ref{theorem_GE0}}
Before proof, we give the following two lemmas.
\begin{lemma}[Covering number of tensor]\label{lemma_cov_num}
Let $\mathcal{S}_{dnk}=\lbrace  \bm{\mathcal{X}}\in\mathbb{R}^{n\times \ldots\times n}: \bm{\mathcal{X}}=\sum_{i=1}^k\bm{x}_i^{(1)}\circ\bm{x}_i^{(2)}\ldots\circ\bm{x}_i^{(d)},\ k\in\mathbb{N};\ \Vert\bm{x}_i^{(j)}\Vert\leq\delta_{ij},\ 1\leq i\leq k,\ 1\leq j\leq d\rbrace$. Then the covering numbers of $\mathcal{S}_{dnk}$ with respect to the Frobenius norm satisfy
$$\mathcal{N}(\mathcal{S}_{dnk},\Vert\cdot\Vert_F,\epsilon)\leq\left(\dfrac{3d}{\epsilon}\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)^{1/2}\right)^{dnk}.$$
\end{lemma}

\begin{lemma}[(Hoeffding inequality for sampling
without replacement]\label{lemma_samp_no_rep}
Let $X_1,X_2,\ldots,Xs$ be a set of samples taken without replacement from a distribution $\lbrace x_1,x_2,\ldots,x_N\rbrace$ of mean $u$ and variance $\sigma^2$. Denote $a=\min_i x_i$ and $b=\max_i x_i$. Then
$$P\left[\Big\vert\dfrac{1}{s}\sum_{i=1}^sX_i-u\Big\vert\geq t\right]\leq2\exp\left(-\dfrac{2st^2}{(1-(s-1)/N)(b-a)^2}\right).$$
\end{lemma}

For convenience, we define
\begin{align*}
\hat{h}(\bm{\mathcal{X}})&=\dfrac{1}{\vert\Omega\vert}\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bm{\mathcal{X}})\Vert_F^2,\\
h(\bm{\mathcal{X}})&=\dfrac{1}{n^d}\Vert\bm{\mathcal{D}}-\bm{\mathcal{X}}\Vert_F^2.
\end{align*}
Suppose $\max\lbrace{\Vert \bm{\mathcal{D}}\Vert_{\infty},\Vert \bm{\mathcal{X}}\Vert_{\infty}\rbrace}\leq\gamma$.
According to Lemma \ref{lemma_samp_no_rep}, we have
$$P\left[\vert\hat{h}-h\vert\geq t\right]\leq2\exp\left(-\dfrac{2\vert\Omega\vert t^2}{(1-(\vert\Omega \vert-1)/n^d)\varsigma^2}\right),$$
where $\varsigma=4\gamma^2$.
Using union bound for all $\bar{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}$ yields
\begin{align*}
&P\left[\sup_{\bar{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\vert\hat{h}(\bar{\bm{\mathcal{X}}})-h(\bar{\bm{\mathcal{X}}})\vert\geq t\right]\\\leq&2\vert {\mathcal{S}}_{dnk}\vert\exp\left(-\dfrac{2\vert\Omega\vert t^2}{(1-(\vert\Omega \vert-1)/n^d)\varsigma^2}\right).
\end{align*}
Or equivalently, with probability at least $1-2n^{-d}$,
\begin{align*}
&\sup_{\bar{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\vert\hat{h}(\bar{\bm{\mathcal{X}}})-h(\bar{\bm{\mathcal{X}}})\vert\\
\leq&\sqrt{\dfrac{\varsigma^2\log\left(\vert {\mathcal{S}}_{dnk}\vert n^d\right)}{2}\left(\dfrac{1}{\vert\Omega\vert}-\dfrac{1}{n^d}+\dfrac{1}{n^d\vert\Omega \vert}\right)}.
\end{align*}
Using Lemma \ref{lemma_cov_num}, we have
\begin{align*}
&g(\Omega)\triangleq\sup_{\bar{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\vert\hat{h}(\bar{\bm{\mathcal{X}}})-h(\bar{\bm{\mathcal{X}}})\vert\\
\leq&\sqrt{\dfrac{\varsigma^2}{2}\left(d\log n+kdn\log \theta\right)\left(\dfrac{1}{\vert\Omega\vert}-\dfrac{1}{n^d}+\dfrac{1}{n^d\vert\Omega \vert}\right)},
\end{align*}
where $\theta=\dfrac{3d}{\epsilon}\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)^{1/2}$.

Since $\vert\sqrt{u}-\sqrt{v}\vert\leq\sqrt{\vert u-v\vert}$ holds for any non-negative $u$ and $v$, we have
$$\sup_{\bar{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\left\vert\sqrt{\hat{h}(\bar{\bm{\mathcal{X}}})}-\sqrt{h(\bar{\bm{\mathcal{X}}})}\right\vert\leq \sqrt{g(\Omega)}.$$
Recall that $\epsilon\geq\Vert \bm{\mathcal{X}}-\bar{\bm{\mathcal{X}}}\Vert_F\geq\Vert \mathcal{P}(\bm{\mathcal{X}}-\bar{\bm{\mathcal{X}}})\Vert_F$, we have
\begin{align*}
&\left\vert\sqrt{{h}({\bm{\mathcal{X}}})}-\sqrt{{h}(\bar{\bm{\mathcal{X}}})}\right\vert\\
=&\dfrac{1}{\sqrt{n^d}}\Big\vert\Vert\bm{\mathcal{D}}-\bm{\mathcal{X}}\Vert_F-\Vert\bm{\mathcal{D}}-\bar{\bm{\mathcal{X}}}\Vert_F \Big\vert\leq\dfrac{\epsilon}{\sqrt{n^d}}
\end{align*}
and
\begin{align*}
&\left\vert\sqrt{\hat{h}({\bm{\mathcal{X}}})}-\sqrt{\hat{h}(\bar{\bm{\mathcal{X}}})}\right\vert\\
=&\dfrac{1}{\sqrt{\vert\Omega\vert}}\Big\vert\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bm{\mathcal{X}})\Vert_F-\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bar{\bm{\mathcal{X}}})\Vert_F \Big\vert\leq\dfrac{\epsilon}{\sqrt{\vert\Omega\vert}}.
\end{align*}
It follows that
\begin{align*}
&\sup_{{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\left\vert\sqrt{\hat{h}({\bm{\mathcal{X}}})}-\sqrt{h({\bm{\mathcal{X}}})}\right\vert\\
\leq &\sup_{{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\left\vert\sqrt{\hat{h}({\bm{\mathcal{X}}})}-\sqrt{\hat{h}(\bar{\bm{\mathcal{X}}})}\right\vert\\
&\qquad\ +\left\vert\sqrt{\hat{h}(\bar{\bm{\mathcal{X}}})}-\sqrt{{h}(\bar{\bm{\mathcal{X}}})}\right\vert+\left\vert\sqrt{{h}(\bar{\bm{\mathcal{X}}})}-\sqrt{h({\bm{\mathcal{X}}})}\right\vert\\
\leq&\dfrac{\epsilon}{\sqrt{\vert\Omega\vert}}+\sqrt{g(\Omega)}+\dfrac{\epsilon}{\sqrt{n^d}}.
\end{align*}
Using the definition of $g(\Omega)$ and letting $\epsilon=3d$, we have
\begin{align*}
&\sup_{{\bm{\mathcal{X}}}\in{\mathcal{S}}_{dnk}}\left\vert\sqrt{\hat{h}({\bm{\mathcal{X}}})}-\sqrt{h({\bm{\mathcal{X}}})}\right\vert\\
\leq&\dfrac{2\epsilon}{\sqrt{\vert\Omega\vert}}+\left(\dfrac{\varsigma^2}{2}\left(d\log n+kdn\log \theta\right)\left(\dfrac{1}{\vert\Omega\vert}-\dfrac{1}{n^d}+\dfrac{1}{n^d\vert\Omega \vert}\right)\right)^{1/4}\\
\leq&\dfrac{6d}{\sqrt{\vert\Omega\vert}}+\left(\dfrac{2\varsigma^2kdn\log \left(\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)^{1/2}\right)}{\vert\Omega\vert}\right)^{1/4}\\
\leq& c\gamma\left(\dfrac{kdn\log \left(\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)\right)}{\vert\Omega\vert}\right)^{1/4},
\end{align*}
where $c$ is some numerical constant.

\section{Proof of Theorem \ref{theorem_GE1}}
\begin{proof}
Denote $\Delta(\Omega)=c\gamma\left(\dfrac{kdn\log \left(\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)\right)}{\vert\Omega\vert}\right)^{1/4}$. According to the definition of $\bm{\mathcal{D}}$ and Theorem \ref{theorem_GE0}, we have
\begin{align*}
&\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{T}}-\bm{\mathcal{X}}\Vert_F=\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{D}}-\bm{\mathcal{N}}-\bm{\mathcal{X}}\Vert_F\\
\leq&\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{D}}-\bm{\mathcal{X}}\Vert_F+\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{N}}\Vert_F\\
\leq&\dfrac{1}{\sqrt{\vert\Omega\vert}}\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\bm{\mathcal{X}})\Vert_F+\Delta(\Omega)+\dfrac{1}{\sqrt{n^d}}\Vert\bm{\mathcal{N}}\Vert_F.
\end{align*}
\end{proof}

\section{Proof of Theorem \ref{theorem_errorbound_sp}}
We first give the following two lemmas.
\begin{lemma}\label{lemma_Nspectral}
Suppose the entries of $\bm{\mathcal{N}}$ are drawn from $\mathcal{N}(0,\sigma^2)$ independently. Then the following inequality holds with probability at least $1-2n^{-d}$
\begin{align*}
\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{N}})\Vert_2\leq 2\sqrt{2}\sigma\sqrt{dn\log(5d)+d\log(n)}.
\end{align*}
\end{lemma}
This is a special case of Corollary 2 in \cite{tomioka2014spectral}.

\begin{lemma}\label{lemma_pless1}
Suppose $0<p<1$ and let $\tilde{n}=\max\lbrace n,k \rbrace$. Let $c_{\kappa}=(2\kappa-1)(2\kappa)\kappa^{-1/(2\kappa-1)}$ where $\kappa=(2-p)/(2-2p)$. Then for any fixed $\nu>0$, $\vartheta\geq C^2$, and $\tau=c_{\kappa}(\vartheta/p)^{1-p/2}(\tilde{n}/\vert\Omega\vert)^{1-p/2}$, we have
\begin{align*}
&\left\vert \dfrac{1}{\vert\Omega\vert}\langle \mathcal{P}_{\Omega}(\bm{\mathcal{N}}),\hat{\bm{\mathcal{T}}}-\bm{\mathcal{X}}\rangle\right\vert\\
\leq&\dfrac{\nu}{2\vert\Omega\vert}\Vert \mathcal{P}_{\Omega}(\hat{\bm{\mathcal{T}}}-\bm{\mathcal{X}})\Vert_F^2+\tau\nu^{p-1}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{X}}\Vert_{S_p}^p
\end{align*}
with probability at least $1-C\exp(-\vartheta \tilde{n}/C^2)$ for some constant $C=C(p,\vert\Omega\vert,\sigma)>0$ independent of $\tilde{n}$, satisfying $\sup_{0<p\leq q}C(p,\vert\Omega\vert,\sigma)<\infty$ for all $q<1$, and $C<p$ provided that $\vert\Omega\vert$ is sufficiently large and $p$ is sufficiently far from $0$.
\end{lemma}

Now we prove Theorem \ref{theorem_sp}.
By the definition of $\hat{\bm{\mathcal{T}}}$, we have
\begin{align*}
&\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{D}}-\hat{\bm{\mathcal{T}}})\Vert_F^2+\lambda\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p\\
\leq &\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}\left(\bm{\mathcal{D}}-\bm{\mathcal{T}}\right)\Vert_F^2+\lambda\Vert \bm{\mathcal{T}}\Vert_{S_p}^p.
\end{align*}
It follows that
\begin{align}\label{Eq_ineq_N}
&\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber\\
\leq &\left\langle \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}, \mathcal{P}_{\Omega}\left(\bm{\mathcal{N}} \right)\right\rangle+\lambda\left(\Vert \bm{\mathcal{T}}\Vert_{S_p}^p-\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p\right).
\end{align}

When $p=1$, we use the trace duality to bound the inner product in \eqref{Eq_ineq_N}:
\begin{align*}
&\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber\\
\leq &\Vert\hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_\ast\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{N}})\Vert_2+\lambda\left(\Vert \bm{\mathcal{T}}\Vert_\ast-\Vert \hat{\bm{\mathcal{T}}}\Vert_\ast\right)\\
\leq &\left(\Vert \bm{\mathcal{T}}\Vert_\ast+\Vert \hat{\bm{\mathcal{T}}}\Vert_\ast\right)\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{N}})\Vert_2+\lambda\left(\Vert \bm{\mathcal{T}}\Vert_\ast-\Vert \hat{\bm{\mathcal{T}}}\Vert_\ast\right).
\end{align*}
Using Lemma \ref{lemma_Nspectral} and letting $\lambda=\Vert\mathcal{P}_{\Omega}(\bm{\mathcal{N}})\Vert_2$, we obtain
\begin{align*}
&\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber\\
\leq &4\sqrt{2}\sigma\sqrt{dn\log(5d)+d\log(n)}\Vert {\bm{\mathcal{T}}}\Vert_\ast
\end{align*}
with probability at least $1-2n^{-d}$.

When $0<p<1$, we use Lemma \ref{lemma_pless1} to bound the inner product in \eqref{Eq_ineq_N}:
\begin{align*}
&\dfrac{1}{2}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber\\
\leq &\dfrac{\nu}{2}\Vert \mathcal{P}_{\Omega}(\hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}})\Vert_F^2+\tau\vert\Omega\vert\nu^{p-1}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_{S_p}^p\\
&+\lambda\left(\Vert \bm{\mathcal{T}}\Vert_{S_p}^p-\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p\right).
\end{align*}
Let $\nu=1/2$ and $\lambda=\tau\vert\Omega\vert\nu^{p-1}$, we have
\begin{align*}
\dfrac{1}{4}\Vert \mathcal{P}_{\Omega}(\bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}})\Vert_F^2 \nonumber
\leq 4\tau\vert\Omega\vert\Vert\bm{\mathcal{T}}\Vert_{S_p}^p,
\end{align*}
where $\tau=c_{\kappa}(\vartheta/p)^{1-p/2}(\tilde{n}/\vert\Omega\vert)^{1-p/2}$. Let $\vartheta=dC^2$. Then $\tau=c_{\kappa}(C^2/p)^{1-p/2}(d\tilde{n}/\vert\Omega\vert)^{1-p/2}$.

\section{Proof of Theorem \ref{theorem_bound_overall} }
\begin{proof}
The theorem can be obtained by combining Theorem \ref{theorem_GE0} and Theorem \ref{theorem_errorbound_sp}.
\end{proof}

\section{Proof of Theorem \ref{theorem_TRPCA}}\label{appendix_the_TRPCA}
\begin{proof}
Let $\hat{\bm{\mathcal{T}}}$ (for $\bm{\mathcal{X}}$) and  $\hat{\bm{\mathcal{E}}}$ be the optimal solution of \eqref{problem_RTPCA_Sp}. We have
\begin{align}
&\dfrac{1}{2}\Vert \bm{\mathcal{D}}-\hat{\bm{\mathcal{T}}}-\hat{\bm{\mathcal{E}}}\Vert_F^2+\lambda_x\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p+\lambda_e\Vert\hat{\bm{\mathcal{E}}}\Vert_1 \nonumber \\
\leq & \dfrac{1}{2}\Vert \bm{\mathcal{D}}-{\bm{\mathcal{T}}}-{\bm{\mathcal{E}}}\Vert_F^2+\lambda_x\Vert {\bm{\mathcal{T}}}\Vert_{S_p}^p+\lambda_e\Vert{\bm{\mathcal{E}}}\Vert_1 \nonumber,
\end{align}
where $\bm{\mathcal{D}}={\bm{\mathcal{T}}}+\bm{\mathcal{N}}+{\bm{\mathcal{E}}}$. It follows that
\begin{align}\label{eq.delta_bound1}
\Delta:=&\dfrac{1}{2}\Vert \bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F^2+\dfrac{1}{2}\Vert {\bm{\mathcal{E}}}-\hat{\bm{\mathcal{E}}}\Vert_F^2 \nonumber \\
\leq &\lambda_x(\Vert {\bm{\mathcal{T}}}\Vert_{S_p}^p-\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p)+\lambda_e(\Vert{\bm{\mathcal{E}}}\Vert_1-\Vert\hat{\bm{\mathcal{E}}}\Vert_1) \nonumber \nonumber \\
& + \left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{N}}\right\rangle+\left\langle\hat{\bm{\mathcal{E}}} -\bm{\mathcal{E}},
\bm{\mathcal{N}}\right\rangle \nonumber \\
& + \left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{E}}-\hat{\bm{\mathcal{E}}}\right\rangle. \nonumber\\
\leq &\lambda_x(\Vert {\bm{\mathcal{T}}}\Vert_{S_p}^p-\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p)+\lambda_e(\Vert{\bm{\mathcal{E}}}\Vert_1-\Vert\hat{\bm{\mathcal{E}}}\Vert_1) \nonumber \nonumber \\
& + \left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{N}}\right\rangle+\Vert\hat{\bm{\mathcal{E}}} -\bm{\mathcal{E}}\Vert_1\Vert\bm{\mathcal{N}}\Vert_\infty \nonumber \\
& + \left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{E}}\right\rangle+(\Vert\bm{\mathcal{T}}\Vert_\infty+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty)\Vert\hat{\bm{\mathcal{E}}}\Vert_1.
\end{align}
When $p=1$, using Lemma \ref{lemma_Nspectral}, we have
\begin{align}
&\left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{N}}\right\rangle+\left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{E}}\right\rangle \nonumber\\
\leq &\Vert\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}}\Vert_\ast(\Vert\bm{\mathcal{N}}\Vert_2+\Vert\bm{\mathcal{E}}\Vert_2)\nonumber\\
\leq &2\sqrt{2}(\sigma+\varsigma)\sqrt{dn\log(5d)+d\log(n)}\Vert\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}}\Vert_\ast
\end{align}
with probability at least $1-4n^{-d}$. Then we arrive at
\begin{align}
\Delta
\leq &\lambda_x(\Vert {\bm{\mathcal{T}}}\Vert_\ast-\Vert \hat{\bm{\mathcal{T}}}\Vert_\ast)+\lambda_e(\Vert{\bm{\mathcal{E}}}\Vert_1-\Vert\hat{\bm{\mathcal{E}}}\Vert_1) \nonumber \nonumber \\
& +(\Vert\hat{\bm{\mathcal{E}}}\Vert_1+\Vert\bm{\mathcal{E}}\Vert_1)\Vert\bm{\mathcal{N}}\Vert_\infty \nonumber+(\Vert\bm{\mathcal{T}}\Vert_\infty+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty)\Vert\hat{\bm{\mathcal{E}}}\Vert_1. \nonumber\\
&+2\sqrt{2}(\sigma+\varsigma)\sqrt{dn\log(5d)+d\log(n)}(\Vert\hat{\bm{\mathcal{T}}}\Vert_\ast+\Vert\bm{\mathcal{T}}\Vert_\ast).
\end{align}
Let $\lambda_x=\bar{\lambda}_x+2\sqrt{2}(\sigma+\varsigma)\sqrt{dn\log(5d)+d\log(n)}$ and $\lambda_e=\bar{\lambda}_e+\Vert\bm{\mathcal{N}}\Vert_\infty+\Vert\bm{\mathcal{T}}\Vert_\infty$, where $\bar{\lambda}_x\geq 0$ and  $\bar{\lambda}_e\geq 0$. We have
\begin{align}
\Delta
\leq &\Big(\bar{\lambda}_x+4\sqrt{2}(\sigma+\varsigma)\sqrt{dn\log(5d)+d\log(n)}\Big)\Vert\bm{\mathcal{T}}\Vert_\ast\nonumber\\
&+(\bar{\lambda}_e+2\Vert\bm{\mathcal{N}}\Vert_\infty+\Vert\bm{\mathcal{T}}\Vert_\infty)\Vert\bm{\mathcal{E}}\Vert_1+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty\Vert\hat{\bm{\mathcal{E}}}\Vert_1.
\end{align}

When $0<p<1$, let $\tilde{n}=\max\lbrace n,k \rbrace$, $c_{\kappa}=(2\kappa-1)(2\kappa)\kappa^{-1/(2\kappa-1)}$, where $\kappa=(2-p)/(2-2p)$. Using Lemma \ref{lemma_pless1}, for any fixed $\nu_1>0$, $\vartheta_1\geq C_1^2$, and $\tau_1=c_{\kappa}(\vartheta_1/p)^{1-p/2}(\tilde{n}/n^d)^{1-p/2}$, we have
\begin{align}\label{eq.TTN_bound}
&\left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{N}}\right\rangle \nonumber\\
\leq &\dfrac{\nu_1}{2}\Vert\hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_F^2+n^d\tau_1\nu_1^{p-1}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_{S_p}^p
\end{align}
with probability at least $1-C_1\exp(-\vartheta_1 \tilde{n}/C_1^2)$ for some constant $C_1=C_1(p,n^d,\sigma)>0$ independent of $\tilde{n}$, satisfying $\sup_{0<p\leq q}C_1(p,n^d,\sigma)<\infty$ for all $q<1$, and $C_1<p$ provided that $p$ is sufficiently far from $0$.

Similarly, for any fixed $\nu_2>0$, $\vartheta_2\geq C_2^2$, and $\tau_2=c_{\kappa}(\vartheta_2/p)^{1-p/2}(\tilde{n}/s_e)^{1-p/2}$, we have
\begin{align}\label{eq.TTE_bound}
&\left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{E}}\right\rangle \nonumber\\
\leq &\dfrac{\nu_2}{2}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_F^2+s_e\tau_2\nu_2^{p-1}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_{S_p}^p
\end{align}
with probability at least $1-C_2\exp(-\vartheta \tilde{n}/C_2^2)$ for some constant $C_2=C_2(p,s_e,\varsigma_2)>0$ independent of $\tilde{n}$, satisfying $\sup_{0<p\leq q}C_2(p,s_e,\sigma)<\infty$ for all $q<1$., and $C_2<p$ provided that $s_e$ is sufficiently large and $p$ is sufficiently far from $0$.

Now letting $\nu_1=\nu_2=1/4$, we get
\begin{align}\label{eq.TTNE_bound}
&\left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{N}}\right\rangle+\left\langle\hat{\bm{\mathcal{T}}} -\bm{\mathcal{T}},
\bm{\mathcal{E}}\right\rangle \nonumber\\
\leq &\dfrac{1}{4}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_F^2+4c_{\kappa}(C_1^2/p)^{1-p/2}\tilde{n}(n^d/\tilde{n})^{p/2}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_{S_p}^p\nonumber\\
&+4c_{\kappa}(C_2^2/p)^{1-p/2}\tilde{n}(s_e/\tilde{n})^{p/2}\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_{S_p}^p
\end{align}
with probability at least $1-2\max(C_1,C_2)\exp(-\tilde{n}/\max(C_1,C_2)^2)$.
Substituting \eqref{eq.TTNE_bound} into \eqref{eq.delta_bound1} , we have
\begin{align}
&\dfrac{1}{2}\Vert \bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F^2+\dfrac{1}{2}\Vert {\bm{\mathcal{E}}}-\hat{\bm{\mathcal{E}}}\Vert_F^2 \nonumber \\
\leq &\lambda_x(\Vert {\bm{\mathcal{T}}}\Vert_{S_p}^p-\Vert \hat{\bm{\mathcal{T}}}\Vert_{S_p}^p)+\lambda_e(\Vert{\bm{\mathcal{E}}}\Vert_1-\Vert\hat{\bm{\mathcal{E}}}\Vert_1) \nonumber \nonumber \\
& +(\Vert\hat{\bm{\mathcal{E}}}\Vert_1+\Vert\bm{\mathcal{E}}\Vert_1)\Vert\bm{\mathcal{N}}\Vert_\infty \nonumber+(\Vert\bm{\mathcal{T}}\Vert_\infty+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty)\Vert\hat{\bm{\mathcal{E}}}\Vert_1\nonumber\\
& + 4c_{\kappa}\tilde{n}\Big((\tfrac{C_1^2}{p})^{1-\frac{p}{2}}(\tfrac{n^d}{\tilde{n}})^{\frac{p}{2}}+(\tfrac{C_2^2}{p})^{1-\frac{p}{2}}(\tfrac{s_e}{\tilde{n}})^{p/2}\Big)\Vert \hat{\bm{\mathcal{T}}}-\bm{\mathcal{T}}\Vert_{S_p}^p.
\end{align}

Let $\lambda_x=\bar{\lambda}_x+4c_{\kappa}\tilde{n}\Big((\tfrac{C_1^2}{p})^{1-\frac{p}{2}}(\tfrac{n^d}{\tilde{n}})^{\frac{p}{2}}+(\tfrac{C_2^2}{p})^{1-\frac{p}{2}}(\tfrac{s_e}{\tilde{n}})^{p/2}\Big)$ and $\lambda_e=\bar{\lambda}_e+\Vert\bm{\mathcal{N}}\Vert_\infty+\Vert\bm{\mathcal{T}}\Vert_\infty$, where $\bar{\lambda}_x\geq 0$ and  $\bar{\lambda}_e\geq 0$.
We have
\begin{align}
&\dfrac{1}{2}\Vert \bm{\mathcal{T}}-\hat{\bm{\mathcal{T}}}\Vert_F^2+\dfrac{1}{2}\Vert {\bm{\mathcal{E}}}-\hat{\bm{\mathcal{E}}}\Vert_F^2 \nonumber \\
\leq &\Big(\bar{\lambda}_x+8c_{\kappa}\tilde{n}\big((\tfrac{C_1^2}{p})^{1-\frac{p}{2}}(\tfrac{n^d}{\tilde{n}})^{\frac{p}{2}}+(\tfrac{C_2^2}{p})^{1-\frac{p}{2}}(\tfrac{s_e}{\tilde{n}})^{p/2}\big)\Big)\Vert {\bm{\mathcal{T}}}\Vert_{S_p}^p \nonumber \\
&+(\bar{\lambda}_e+2\Vert\bm{\mathcal{N}}\Vert_\infty+\Vert\bm{\mathcal{T}}\Vert_\infty)\Vert\bm{\mathcal{E}}\Vert_1+\Vert\hat{\bm{\mathcal{T}}}\Vert_\infty\Vert\hat{\bm{\mathcal{E}}}\Vert_1.
\end{align}
Letting $\bar{\lambda}_x=\bar{\lambda}_e=0$, we finished the proof.
\end{proof}

\section{Proof of Lemma \ref{lemma_cov_num}}
\begin{proof}
We rewrite $\bm{\mathcal{X}}$ as
$$\bm{\mathcal{X}}=\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\bm{X}^{(2)}\ldots\times_{d}\bm{X}^{(d)}.$$
Let $S_{nk}=\lbrace \bm{X}\in\mathbb{R}^{n\times k}: \Vert\bm{X}\Vert_F\leq\alpha\rbrace$. Then there exists an $\epsilon$-net $\bar{S}_{nk}$ obeying
$$\mathcal{N}(S_{nk},\Vert\cdot\Vert_F,\epsilon)\leq\left(\dfrac{3\alpha}{\epsilon}\right)^{nk}$$
such that  $\Vert\bm{X}-\bar{\bm{X}}\Vert_F\leq\epsilon$.
Replacing $\epsilon$ with $\epsilon/\beta$, we have
$\mathcal{N}(S_{nk},\Vert\cdot\Vert_F,\epsilon/\beta)\leq\left(\dfrac{3\alpha\beta}{\epsilon}\right)^{nk}$.
Let $\alpha\beta=d\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)^{1/2}$, then
$\Vert\bm{X}^{(j)}-\bar{\bm{X}}^{(j)}\Vert_F\leq\epsilon/\beta_j$,
where
$\beta_j=d\prod_{l\in[d],l\neq j}\left(\sum_{i=1}^k\delta_{il}^2\right)^{1/2}$.
We have
\begin{align*}
&\Vert\bm{\mathcal{X}}-\bar{\bm{\mathcal{X}}}\Vert_F\\
=&\Vert\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\ldots\times_{d}\bm{X}^{(d)}-\bm{\mathcal{I}}\times_{1}\bar{\bm{X}}^{(1)}\times_{2}\ldots\times_{d}\bar{\bm{X}}^{(d)}\Vert_F\\
=&\Vert\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\ldots\times_{d}\bm{X}^{(d)}\pm\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\ldots\times_{d}\bar{\bm{X}}^{(d)}\\
&\pm\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\ldots\times_{d-1}\bar{\bm{X}}^{(d-1)}\times_{d}\bar{\bm{X}}^{(d)}\\
&\pm\ldots\pm\bm{\mathcal{I}}\times_{1}\bar{\bm{X}}^{(1)}\times_{2}\ldots\times_{d-1}\bar{\bm{X}}^{(d-1)}\times_{d}\bar{\bm{X}}^{(d)}\\
&-\bm{\mathcal{I}}\times_{1}\bar{\bm{X}}^{(1)}\times_{2}\ldots\times_{d}\bar{\bm{X}}^{(d)}\Vert_F\\
\leq&\Vert\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\ldots\times_{d}(\bm{X}^{(d)}-\bar{\bm{X}}^{(d)})\Vert_F\\
&+\Vert\bm{\mathcal{I}}\times_{1}\bm{X}^{(1)}\times_{2}\ldots\times_{d-1}(\bm{X}^{(d-1)}-\bar{\bm{X}}^{(d-1)})\times_{d}\bar{\bm{X}}^{(d)}\Vert_F\\
&+\ldots+\Vert\bm{\mathcal{I}}\times_{1}(\bm{X}^{(1)}-\bar{\bm{X}}^{(d-1)})\times_{2}\bar{\bm{X}}^{(2)}\ldots\times_{d}\bar{\bm{X}}^{(d)}\Vert_F\\
\leq& \Vert\bm{\mathcal{I}}\Vert_2\Vert\bm{X}^{(1)}\Vert_F\ldots\Vert\bm{X}^{(d-1)}\Vert_F\Vert\bm{X}^{(d)}-\bar{\bm{X}}^{(d)}\Vert_F\\
&+\Vert\bm{\mathcal{I}}\Vert_2\Vert\bm{X}^{(1)}\Vert_F\ldots\Vert\bm{X}^{(d-1)}-\bar{\bm{X}}^{(d-1)}\Vert_F\Vert\bar{\bm{X}}^{(d)}\Vert_F\\
&+\ldots+\Vert\bm{\mathcal{I}}\Vert_2\Vert\bm{X}^{(1)}-\bar{\bm{X}}^{(1)}\Vert_F\Vert\bar{\bm{X}}^{(2)}\Vert_F\ldots\Vert\bar{\bm{X}}^{(d-1)}\Vert_F\Vert\bar{\bm{X}}^{(d)}\Vert_F\\
\leq&\Vert\bm{X}^{(1)}\Vert_F\ldots\Vert\bm{X}^{(d-1)}\Vert_F\beta_d
+\Vert\bm{X}^{(1)}\Vert_F\ldots\beta_{d-1}\Vert_F\Vert\bar{\bm{X}}^{(d)}\Vert_F\\
&+\ldots+\beta_1\Vert\bar{\bm{X}}^{(2)}\Vert_F\ldots\Vert\bar{\bm{X}}^{(d-1)}\Vert_F\Vert\bar{\bm{X}}^{(d)}\Vert_F\\
\leq&\epsilon
\end{align*}
Thus $\bar{\mathcal{S}}_{dnk}$ is an $\epsilon$-cover of $\mathcal{S}_{dnk}$. Therefore, we have
\begin{align*}
&\mathcal{N}(\mathcal{S}_{dnk},\Vert\cdot\Vert_F,\epsilon)\\
\leq&\prod_{j=1}^d\left(\dfrac{3\alpha_j\beta_j}{\epsilon}\right)^{nk}\\
=&\prod_{j=1}^d\left(\dfrac{3d\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)^{1/2}}{\epsilon}\right)^{nk}\\
=&\left(\dfrac{3d\prod_{j=1}^d\left(\sum_{i=1}^k\delta_{ij}^2\right)^{1/2}}{\epsilon}\right)^{dnk}.
\end{align*}
\end{proof}

\section{Proof of Lemma \ref{lemma_pless1}}
\begin{proof}
The lemma is a generalization of Lemma 5 in \cite{rohde2011estimation}. The main difference is that in tensor, the rank can be larger than the side lengths. So the entropy number of tensor Schatten-$p$ quasi-norm ball should be changed accordingly. Specifically, replace the $M$ of Proposition 1 and Corollary 7  in \cite{rohde2011estimation} with $\max(k,n)$, and
replace $\mathcal{G}_{M,p}$ (page 37 of \cite{rohde2011estimation}) with $$
\mathcal{G}_{n,p}=\left\lbrace \bm{\mathcal{A}}\in\mathbb{R}^{n\times n\ldots\times n}: \left(\tfrac{\max(k,n)}{p}\right)^{\tfrac{2-p}{2p}}\bm{\mathcal{A}}\in\mathcal{B}(S_p^n)\right\rbrace.$$
Then following the same procedures in \cite{rohde2011estimation}, we accomplish the proof.

Note that according to Claim 1 of \cite{rohde2011estimation} and using the corresponding notations, we have
$$T\geq c(q)c_0^{p/(2(2-p))}p\alpha_0(p)\dfrac{2-p}{2-2p}\triangleq \tilde{c}(q),$$
where $\alpha_0(p)\sim O(1/p)$ and $c_0=1/N$. Hence there exists a $c(q)$ such that $\tilde{c}(q)<p$ when $p$ stays away from $0$ and $1$ and $N$ is large enough. Similar result can be expanded to Claim 2 of \cite{rohde2011estimation}. Therefore, in our Lemma \ref{lemma_pless1}, there exists $C<p$ provided that $p$ is far from $0$ and $1$ and $\vert\Omega\vert$ is large enough.

\end{proof}

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
\else

  \section*{Acknowledgment}
\fi

The authors would like to thank...

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}