*5.1. Synthetic Data*

We first verify the correctness of Theorem 1. Specifically, we check whether the following two statements indicated in Theorem 1 hold in experiments on synthetic data sets:


**Signal Generation.** With a given tubal rank *r*0, we first generate the underlying tensor L0 <sup>∈</sup> <sup>R</sup>*n*1×*n*2×*n*<sup>3</sup> by L0 <sup>=</sup> <sup>A</sup> <sup>∗</sup> B/*n*3, where tensors <sup>A</sup> <sup>∈</sup> <sup>R</sup>*n*1×*r*0×*n*<sup>3</sup> and <sup>B</sup> <sup>∈</sup> <sup>R</sup>*r*0×*n*2×*n*<sup>3</sup> are generated with *i.i.d.* standard Gaussian elements. Then, the sparse corruption tensor S0 is generated by choosing its support uniformly at random. The non-zero elements of S0 will be *i.i.d.* sampled from a certain distribution that will be specified afterwards. Furthermore, the noise tensor E0 is generated with entries sampled *i.i.d.* from N (0, *<sup>σ</sup>*2) with *<sup>σ</sup>* = *<sup>c</sup>*L0F/ <sup>√</sup>*n*1*n*2*n*3, where we set constant *<sup>c</sup>* is to control the signal noise ratio. Finally, the observed tensor M is formed by M = L0 + S0 + E0.

### 5.1.1. Exact Recovery in the Noiseless Setting

We first check *Statement (I)*, i.e., exact recovery in the noiseless setting. Specifically, we will show that Algorithm 1 and Algorithm 2 can exactly recover the underlying tensor L0 and the sparse corruption S0. We first test the recovery performance of different tensor sizes by setting *<sup>n</sup>* = *<sup>n</sup>*<sup>1</sup> = *<sup>n</sup>*<sup>2</sup> ∈ {100, 160, 200} and *<sup>n</sup>*<sup>3</sup> = 20, with (*r*tubal(L0), S00)=(0.05*n*, 0.05*n*2*n*3). The non-zero elements of tensor S0 is sampled from *i.i.d.* symmetric Bernoulli distribution, i.e., the possibility of being 1 or −1 are 1/2. The results are shown in Table 1. It can be seen that both Algorithm 1 and Algorithm 2 can obtain relative standard error (RSE) smaller than 1e − 5 by which we can say that L0 and S0 are exact recovered. We can also see that Algorithm 2 runs much faster than Algorithm 1.

**Table 1.** Performance of Algorithm 1 and Algorithm 2 in both accuracy and speed for different tensor sizes when the gross corruption. Outliers from symmetric Bernoulli, observation tensor <sup>M</sup> <sup>∈</sup> <sup>R</sup>*n*×*n*×*n*<sup>3</sup> , *<sup>n</sup>*<sup>3</sup> <sup>=</sup> 30, *<sup>r</sup>*tubal(L0) = 0.05*n*, S0<sup>1</sup> <sup>=</sup> 0.05*n*2*n*3, noise level *<sup>c</sup>* <sup>=</sup> 0, *<sup>r</sup>* <sup>=</sup> max 2*r*tubal(L0), 15 .


We then test whether the recovery performance can be affected by the distribution of the corruptions. This is done by choosing the non-zeros elements of S0 from *i.i.d.* standard Gaussian distribution. The experimental results are reported in Table 2. We can find that both Algorithm 1 and Algorithm 2 can exactly recover the true L0 and S0 and Algorithm 2 runs much faster than Algorithm 1.

We also conduct STPCP by Algorithm 1 and Algorithm 2 with missing entries. After generating L0, S0 and E0, we get the observation by Model (27). We choose the support of B uniformly at random with possibility 0.8 and then set elements in the chosen support to be 1. Thus, %20 of the entries are missing. The corrupted observation *M* is then formed by M = B (L0 + S0 + E0). We show the recover results in Table 3. We can see that the underlying low-rank tensor L0 can be exactly recovered and the observed part of the corruption tensor B S0 can also be exactly recovered (Please note that it is impossible to recover the unobserved entries of a sparse tensor S0 [52]).

**Table 2.** Performance of Algorithm 1 and Algorithm 2 in both accuracy and speed for different tensor sizes when the gross corruption. Outliers from standard Gaussian distribution, observation tensor <sup>M</sup> <sup>∈</sup> <sup>R</sup>*n*×*n*×*n*<sup>3</sup> , *<sup>n</sup>*<sup>3</sup> <sup>=</sup> 30, *<sup>r</sup>*tubal(L0) = 0.05*n*, S0<sup>1</sup> <sup>=</sup> 0.05*n*2*n*3, noise level *<sup>c</sup>* <sup>=</sup> 0, *<sup>r</sup>* <sup>=</sup> max 2*r*tubal(L0), 15 .


**Table 3.** Performance of Algorithm 1 and Algorithm 2 in both accuracy and speed for different tensor sizes when the gross corruption. Outliers from symmetric Bernoulli, observation tensor <sup>M</sup> <sup>∈</sup> <sup>R</sup>*n*×*n*×*n*<sup>3</sup> , *<sup>n</sup>*<sup>3</sup> <sup>=</sup> 30, *<sup>r</sup>*tubal(L0) = 0.05*n*, S0<sup>1</sup> <sup>=</sup> 0.05*n*2*n*3, noise level *<sup>c</sup>* <sup>=</sup> 0, *<sup>r</sup>* <sup>=</sup> max 2*r*tubal(L0), 15 , with %20 random missing entries.


### 5.1.2. Linear Scaling of Errors with the Noise Level

We then verify *Statement (II)* that the estimation errors have linear scale behavior with respect to the noise level. The estimation errors are measured using the mean-squared-error (MSE):

$$\text{MSE}(\underline{\text{L}}) = \frac{||\underline{\text{L}} - \underline{\text{L}}\_{0}||\_{\text{F}}^{2}}{n\_{1}n\_{2}n\_{3}}, \quad \text{MSE}(\underline{\text{S}}) = \frac{||\underline{\text{S}} - \underline{\text{S}}\_{0}||\_{\text{F}}^{2}}{n\_{1}n\_{2}n\_{3}}.$$

for the low rank component L0 and the sparse component S0, respectively. We test tensors of 3 different size by choosing *n* ∈ {60, 80, 100} and *n*<sup>3</sup> = 20. The tubal rank *r*tubal(L0) of L0 and sparsity *s* of S0 are set as (*r*tubal(L0),*s*)=(5, 0.1*n*2*n*3). We vary the signal noise ratio *c* = 0.03 : 0.03 : 0.6 which is in proportional of the noise level *δ*. We run the proposed Algorithm 1, test 50 trials, and report the averaged MSEs. The MSEs of Lˆ and S0 versus *c*<sup>2</sup> are shown in sub-figures (a) and (b) in Figure 4. We can see that the estimation error has linear scaling behavior along with the noise level as Theorem 1 indicates. Since the results for *n* = 80 and *n* = 100 are quite similar to the case of *n* = 60, they are simply omitted.

**Figure 4.** The MSEs of Lˆ and S0 versus *<sup>c</sup>*<sup>2</sup> for tensors of size 60 <sup>×</sup> <sup>60</sup> <sup>×</sup> 20 where the tubal rank *<sup>r</sup>*tubal(L0) of L0 and sparsity *s* of S0 are set as (*r*tubal(L0),*s*)=(5, 0.1*n*2*n*3). (**a**): MSE of Lˆ vs *c*2. (**b**): MSE of Sˆ vs *c*2.

### *5.2. Real Data Sets*

In this section, experiments on real data sets (color images and videos) are carried out to evaluate the effectiveness and efficiency of the proposed Algorithms 1 and 2. Besides noises and sparse corruptions, we also consider missing values which is more challenging. The proposed algorithms are compared with the following typical models:

(I). NN-I: tensor recovery based on matrix nuclear norms of frontal slices formulated as follows:

$$\min\_{\mathbf{L}, \mathbf{S}} \frac{1}{2} ||\underline{\mathbf{B}} \odot (\underline{\mathbf{M}} - \underline{\mathbf{L}} - \underline{\mathbf{E}})||\_{\mathbf{F}} + \gamma \sum\_{k=1}^{n\_{\mathbf{S}}} (||\mathbf{L}\_{\{k\}}||\_{\*} + \lambda ||\mathbf{S}^{(k)}||\_{1}).\tag{53}$$

This model will be used for image restoration in Section 5.2.1. Please note that Model (53) is equivalent to parallel matrix recovery on each frontal slice.

(II). NN-II: tensor recovery based on matrix nuclear norm formulated as follows:

$$\min\_{\mathbf{L}\in\underline{\mathbf{S}}} \frac{1}{2} \|\underline{\mathbf{B}} \odot (\underline{\mathbf{M}} - \underline{\mathbf{L}} - \underline{\mathbf{E}})\|\_{\mathbf{F}} + \gamma \|\mathbf{L}\|\_{\*} + \gamma \lambda \|\underline{\mathbf{S}}\|\_{1\prime} \tag{54}$$

where **<sup>L</sup>** = [**l**1,**l**2, ··· ,**l***n*<sup>3</sup> ] <sup>∈</sup> <sup>R</sup>*n*1*n*2×*n*<sup>3</sup> with **<sup>l</sup>***<sup>k</sup>* :<sup>=</sup> vec(**L**(*k*)) <sup>∈</sup> <sup>R</sup>*n*1*n*<sup>2</sup> defined as the vectorization [40] of frontal slices **L**(*k*) , for all *k* = 1, 2, ··· , *n*3. This model will be used for video restoration in Section 5.2.2.

(III). SNN: tensor recovery based on SNN formulated as follows:

$$\min\_{\mathsf{LLS}} \frac{1}{2} \|\mathsf{B} \odot (\mathsf{M} - \mathsf{L} - \mathsf{E})\|\_{\mathsf{F}} + \gamma \sum\_{i=1}^{3} \alpha\_{m} \|\mathsf{L}\_{(i)}\|\_{\*} + \gamma \|\mathsf{S}\|\_{1\prime} \tag{55}$$

where **<sup>L</sup>**(*i*) <sup>∈</sup> <sup>R</sup>*ni*×∏*j*=*<sup>i</sup> nj* is the mode-*<sup>i</sup>* matriculation of tensor L <sup>∈</sup> <sup>R</sup>*n*1×*n*2×*n*<sup>3</sup> , for all *<sup>i</sup>* <sup>=</sup> 1, 2, 3.

We solve the above Model (53)–(55) using ADMM implemented by ourselves in Matlab. The effectiveness of the algorithms is measured by Peaks Signal Noise Ratio (PSNR):

$$\text{PSNR} := 10 \log\_{10} \left( \frac{n\_1 n\_2 n\_3 ||\underline{\mathbf{L}}\_0||\_\infty^2}{||\hat{\underline{\mathbf{L}}} - \underline{\mathbf{L}}\_0||\_\mathbb{F}^2} \right).$$

Please note that a larger PSNR value indicates higher quality of Lˆ.

### 5.2.1. Color Images

Color images are the most commonly used 3-way tensors. We test the twenty 256-by-256-by-3 color images which have been used in [37], and carry out robust tensor recovery with missing entries (see Figure 5). Following [37], for a color image L0 <sup>∈</sup> <sup>R</sup>*n*×*n*×3, we choose its support uniformly at random with ratio *ρ*<sup>s</sup> and fill in the values with *i.i.d.* symmetric Bernoulli variables to generate S0. The noise tensor E0 is generated with *i.i.d.* zero-mean Gaussian entries whose standard deviation is given by *σ* = 0.05L0F/ √ 3*n*2. Then, we form the binary observation mask B by choosing its support uniformly at random with ratio *ρ*obs. Finally, the partially observed corruption M = B (L0 + S0 + E0) are formed.

**Figure 5.** The 20 color images used.

We consider two scenarios by setting (*ρ*obs, *ρ*s) ∈ {(0.9, 0.1),(0.8, 0.2)}. For NN (Model (53)), we set the regularization parameters *<sup>λ</sup>* <sup>=</sup> 1/√*nρ*obs (suggested by [46]), and set parameter *<sup>γ</sup>* <sup>=</sup> E0sp where E0sp is estimated as 6.5*σ* 3*ρ*obs*<sup>n</sup>* log(6*n*) (suggested by [5]). For SNN, the parameters are chosen to satisfy *<sup>γ</sup>* = 0.05, *<sup>α</sup>*<sup>1</sup> = *<sup>α</sup>*<sup>2</sup> = 3*nρ*obs, *<sup>α</sup>*<sup>3</sup> = 0.013*nρ*obs. For Algorithm 1 and Algorithm 2, we set *γ* = 0.3E0sp, and *λ* = 1/ 3*nρ*obs. The initialized rank *<sup>r</sup>* in Algorithm 2 is set as 60. In each setting, we test each color image for 10 times and report the averaged PSNR and time. For quantitative comparison, we show the PSNR values and running times in Figures 6 and 7 for settings of (*ρ*obs, *ρ*s)=(0.9, 0.1) and (*ρ*obs, *ρ*s)=(0.8, 0.2), respectively. Several visual examples are shown in Figure 8 for qualitative comparison for the setting of (*ρ*obs, *ρ*s)=(0.8, 0.2). We can see from Figures 6–8 that the proposed Algorithm 1 has the highest recovery quality and the proposed Algorithm 2 has the second highest quality but the fastest running time.

**Figure 6.** The quantitative comparison in PSNR and time on color images. First, 10% entries of each image is corrupted by *i.i.d.* symmetric Bernoulli variable, then polluted by Gaussian noise of noise level *c* = 0.05, and finally 10% of the corrupted entries are missing uniformly at random. (**a**): the PSNR values of each algorithm; (**b**): the running time of each algorithm.

**Figure 7.** The quantitative comparison in PSNR and time on color images. First, 20% entries of each image is corrupted by *i.i.d.* symmetric Bernoulli variable, then polluted by Gaussian noise of noise level *c* = 0.05, and finally 20% of the corrupted entries are missing uniformly at random. (**a**): the PSNR values of each algorithm; (**b**): the running time of each algorithm.

**Figure 8.** The visual results for image recovery of different algorithms. First, 20% entries of each image is corrupted by *i.i.d.* symmetric Bernoulli variable, then polluted by Gaussian noise of noise level *c* = 0.05, and finally 20% of the corrupted entries are missing uniformly at random. (**a**): the original image; (**b**): the corrupted image; (**c**) image recovered by Algorithm 1; (**d**) image recovered by Algorithm 2; (**e**) image recovered by the matrix nuclear norm (NN)-based Model (53); (**f**) image recovered by the SNN-based Model (55).

### 5.2.2. Videos

In this subsection, video restoration experiments are conducted on four broadly used YUV videos (They can be downloaded from https://sites.google.com/site/subudhibadri/fewhelpfuldownloads: Akiyo\_qcif, Scilent\_qcif, Carphone\_qcif, and Claire\_qcif.) Owing to computational limitation, we simply use the first 32 frames of the Y components of all the videos which results in four 144-by-176-by-30 tensors. For a 3-way data tensor L0 <sup>∈</sup> <sup>R</sup>*n*1×*n*2×*n*<sup>3</sup> , To generate corruption S0, the support is chosen uniformly at random with ratio *ρ*<sup>s</sup> and then elements in the support are filled in with *i.i.d.* symmetric Bernoulli variables. The noise tensor E0 is also generated with *i.i.d.* zero-mean Gaussian entries whose standard deviation is given by *σ* = 0.05L0F/ <sup>√</sup>*n*1*n*2*n*3. Then, the binary observation mask B is formed thorough choosing its support uniformly at random with ratio *ρ*obs. Finally, the partially observed corruption M = B (L0 + S0 + E0) are formed.

We also consider two scenarios by setting (*ρ*obs, *ρ*s) ∈ {(0.9, 0.1),(0.8, 0.2)}. NN-II Model (54) is used in video restoration. For NN-II, we set the regularization parameters *<sup>λ</sup>* <sup>=</sup> 1/√*n*1*n*2*ρ*obs (suggested by [46]), and set parameter *γ* = E0sp where E0sp is estimated as 6.5*σ ρ*obs*n*1*n*<sup>3</sup> log((*n*<sup>1</sup> + *<sup>n</sup>*2)*n*3) (suggested by [5]). For SNN, the parameters are chosen to satisfy *<sup>γ</sup>* <sup>=</sup> 0.05, *<sup>α</sup>*<sup>1</sup> <sup>=</sup> *<sup>α</sup>*<sup>2</sup> <sup>=</sup> <sup>√</sup>*n*1*n*3*ρ*obs, *<sup>α</sup>*<sup>3</sup> <sup>=</sup> 5 <sup>√</sup>*n*1*n*3*ρ*obs. For Algorithm <sup>1</sup> and Algorithm 2, we set *<sup>γ</sup>* <sup>=</sup> 0.3E0sp, and *<sup>λ</sup>* <sup>=</sup> 1/ max{*n*1, *<sup>n</sup>*2}*n*3*ρ*obs after careful parameter tuning. The initialized rank *r* in Algorithm 2 is set as 60. In each setting, we test each video for 10 times and report the averaged PSNR and time. For quantitative comparison, we show the PSNR values and running times in Table 4. It can be seen that Algorithm 1 has the highest recovery quality and the proposed Algorithm 2 has the second highest quality but the fastest running time.


**Table 4.** PSNR values and running time (in seconds) of different algorithms on video data. First, *ρ*s*n*1*n*2*n*<sup>3</sup> entries of each image is corrupted by *i.i.d.* symmetric Bernoulli variable, then polluted by Gaussian noise of noise level *c* = 0.05, and finally (1 − *ρ*obs)*n*1*n*2*n*<sup>3</sup> of the corrupted entries are missing uniformly at random. The items with **highest PSNR values** are highlighted with bold face, and the items with shortest running time are highlighted with underline.

#### **6. Conclusions**

This paper studied the problem of stable tensor principal component pursuit which aims to recover a tensor from noises and sparse corruptions. We proposed a constrained tubal nuclear norm-based model and established upper bounds on the estimation error. In contrast to prior work [37], our theory can guarantee exact recovery in the noiseless setting. We also designed two algorithms, the first ADMM algorithm can be accelerated by the second Algorithm which adopts a factorization strategy. We validated the correctness of our theory by simulations on synthetic data, and evaluated the effectiveness and efficiency of the proposed algorithms via experiments on color images and videos.

For future directions, it is a natural and interesting extension to consider recovery of 4-way tensors [35] with arbitrary linear transformation [53,54]. It is also interesting to use tensor factorization-based methods [55,56] for STPCP. Another challenging future direction is developing tools to verify whether the unknown tensor satisfies the tensor incoherence condition from its incomplete or corrupted observations.

For extensions of the proposed approach to higher-way tensors, we produce the following two ideas:


**Author Contributions:** Conceptualization, W.F. Data curation, D.W. and R.Z. Formal analysis, W.F. Methodology, W.F., D.W and R.Z. Software, D.W. and R.Z. Writing, original draft, W.F., D.W. and R.Z.

**Funding:** This research was funded by the Key Projects of Natural Science Research in Universities in Anhui Province under grant number KJ2019A0994.

**Acknowledgments:** We sincerely thank Andong Wang who shared the codes of [37] and gave us some suggestions of the proof.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A. Proofs of Lemmas and Theorems**
