Square Root Principal Component Pursuit: Tuning-Free Noisy Robust Matrix Recovery
Year: 2021

Abstract:
We propose a new framework – Square Root Principal Component Pursuit – for low-rank matrix recovery from observations corrupted with noise and outliers. Inspired by the square root Lasso, this new formulation does not require prior knowledge of the noise level. We show that a single, universal choice of the regularization parameter suffices to achieve reconstruction error proportional to the (a priori unknown) noise level. In comparison, previous formulations such as stable PCP rely on noise-dependent parameters to achieve similar performance, and are therefore challenging to deploy in applications where the noise level is unknown. We validate the effectiveness of our new method through experiments on simulated and real datasets. Our simulations corroborate the claim that a universal choice of the regularization parameter yields near optimal performance across a range of noise levels, indicating that the proposed method outperforms the (somewhat loose) bound proved here.

heading: 1 Introduction
text: The problem of recovering a low-rank matrix from unreliable observations arises in a wide range of engineering applications, including collaborative filtering [1], latent semantic indexing [2], image and video analysis [3, 4, 5] and so on. This problem can be formalized in terms of the following observation model: given an observationD which is a superposition
D = L0 low-rank + S0 sparse + Z0 noise
. (1.1)
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
of an unknown low-rank matrix L0, sparse corruptions S0 and dense noise Z0, our goal is to accurately estimate both L0 and S0. This model has been intensely studied, leading to algorithmic theory for methods based on both convex and nonconvex optimization [6, 7, 8]. One virtue of the convex approach is that in the noise-free setting (Z0 = 0), it is possible to exactly recover a broad range of low-rank and sparse pairs (L0,S0), with a universal choice of regularization parameters, which does not depend on either the rank or sparsity. This makes it possible to deploy this method in a “hands-free” manner, provided the dataset of interest indeed has low-rank and sparse structure. In the presence of noise, however, the situation becomes more complicated: all efficient, guaranteed estimators require knowledge of the noise level (or the rank and sparsity) [8, 9, 10]. This is problematic, since in most applications the noise level is not known ahead of time. In standard convex formulations, the appropriate regularization parameter depends on the noise standard deviation, leaving the user with a painful and time-consuming task of tuning these parameters on a per-dataset basis. Motivated by this issue, we revisit this classical matrix recovery problem. The main contribution of this paper is the proposal and analysis of a new formulation for robust matrix recovery, which stably recovers L0 and S0 without requiring prior knowledge of the rank, sparsity, or noise level. In particular, our approach admits a single, universal choice of regularization parameters, which under standard hypotheses on L0 and S0, yields an estimation error proportional to the noise standard deviation σ. To our knowledge, our method and analysis are the first to achieve this. Our approach is based on a combination of two natural ideas. For matrix recovery, we draw on the stable principal component pursuit [9], a natural convex relaxation, which minimizes a combination of the nuclear norm of L, the `1 norm of S and the squared Frobenius norm ‖Z‖2F of the noise. This is a principled approach to handling both the structured components L0,S0 and the noise: ‖Z‖2F can bemotivated naturally from the negative log-likelihood of the gaussian distribution. Moreover, under mild assumptions on the rank and singular vectors of L0 and the sparsity pattern of S0, the reconstruction error of stable PCP is O(‖Z0‖F ) [9]. On the other hand, optimally balancing these terms requires knowledge of the standard deviation σ of the true noise distribution. To address this issue, we draw inspiration from the square root Lasso [11]. The square root Lasso is a sparse estimator which achieves minimax optimal estimation with a universal choice of parameters, which does not depend on the noise level. The core idea is very simple: instead of penalizing the squared Frobenius norm ‖Z‖2F of the noise, one penalizes its square root, ‖Z‖F . We call the resulting formulation square root principal component pursuit ( √ PCP ). Our new formulation has the benefit that with a noise-independent universal choice of regularization parameters, essentially the same level of reconstruction error can be achieved. This makes √ PCP a more practical approach to low-rank recovery in unknown noise. Due to the square root term, the objective function is no longer smooth or differentiable, and so we cannot apply algorithms such as the proximal gradient method1. Nevertheless, our new formulation remains convex and separable, i.e. the objective is the sum of functions of different variables, making Alternating Direction Method of Multipliers (ADMM) a suitable solver [12]. We test our new formulation with ADMM on both simulated data and real data in image processing. The experimental results show the effectiveness of our proposed new formulation in recovering the low-rank and the sparse matrix, and suggest that √ PCP has better performance than anticipated by our loose upper bound of the reconstruction error.

heading: 1.1 Notations and Assumptions
text: We use ‖A‖, ‖A‖F , and ‖A‖∗ to denote the spectral, Frobenius, and nuclear norm of the matrixA , andA∗ its (conjugate) transpose. For convenience, we letX0 = (L0,S0) be the concatenation ofL0 andS0. We assume thatL0 is a low-rankmatrix of rank rwhose compact
1which requires the objective to be the sum of a smooth and a non-smooth function
SVD is L0 = UΣV ∗, U ∈ Rn1×r, V ∈ Rn2×r, and without loss of generality, n1 ≥ n2. Let T = {UQ∗ + RV ∗|Q ∈ Rn2×r,R ∈ Rn1×r} denote the tangent space of rank r matrices at L0. In addition, we assume that S0 is sparse with support in Ω. Since it is impossible to disentangle L0 and S0 if the low-rank matrix L0 is sparse, or if the sparse matrix S0 is low-rank, we make the following two assumptions:
Assumption 1.1 The low-rank matrix L0 satisfies the incoherence property with parameter ν2, i.e.
max i ‖U∗ei‖2 ≤
νr n1 , max i ‖V ∗ei‖2 ≤ νr n2 , ‖UV ∗‖∞ ≤
√ νr
n1n2 .
Assumption 1.2 The support Ω is chosen uniformly among all sets of cardinalitym, and the signs of supports are random, i.e. P [(S0)i,j > 0|(i, j) ∈ Ω] = P [(S0)i,j < 0|(i, j) ∈ Ω] = 0.5.
These assumptions follow [6]; indeed, our proof makes use of a dual certificate constructed for noiseless low-rank and sparse recovery in that paper.

heading: 1.2 Problem Formulation and Main Results
text: Inspired by square root Lasso, we propose to solve the robust noisy matrix recovery problem through the following optimization problem:
√ PCP : min
L,S ‖L‖∗ + λ‖S‖1 + µ ‖L + S −D‖F . (1.2)
The parameter λ that balances the low-rank and the sparse regularizers is studied in [6], where it is shown thatλ = 1/√n1 gives exact recoverywhenZ0 = 0 and theµ ‖L + S −D‖F penalty term in (1.2) is replaced with the constraint L + S = D. In this work, we build on this result and focus on the parameter µ. Our main result is that under the aforementioned (standard) hypotheses on L0 and S0, using a single, universal choice µ = √ n2/2, √ PCP recovers L0 and S0, with an estimation error that is proportional to the norm of the noise:
Theorem 1.1 Under Assumptions 1.1 and 1.2, provided that r,m satisfies
r ≤ ρrn2ν−1(log n1)−2, m ≤ ρsn1n2, (1.3)
where ρr ≤ 1/10, ρs are some positive constants. Then there is a numerical constant c such that with probability at least 1− cn−101 , the √ PCP problem (1.2) with λ = 1/√n1 and µ = √ n2/2 produces
a solution X̂ = (L̂, Ŝ) such that
‖X̂ −X0‖F ≤ 560 √ n1n2‖Z0‖F . (1.4)
Why is it possible to achieve accurate estimation with a single choice of µ? We draw intuition from a connection to the stable principal component pursuit formulations studied in [9]. This work studies both constrained and unconstrained formulations:
StablePCPc : min L,S ‖L‖∗ + λ‖S‖1 s.t. ‖L + S −D‖F ≤ δ. (1.5)
StablePCPu : min L,S ‖L‖∗ + λ‖S‖1 +
µ̄ 2 ‖L + S −D‖2F . (1.6)
These formulations are equivalent, and equivalent to √ PCP in the following sense: for each problem instance, there is a calibration of parameters δ ↔ µ̄ ↔ µ such that √ PCP , StablePCPu and StablePCPc have exactly the same set of optimal solutions. However, (1.5)- (1.6) require that the parameters δ and µ̄ be determined on an instance-by-instance basis, based on the noise level. Choosing these parameters correctly is essential: in (1.5), δ should be chosen to be larger than ‖Z0‖F . For square (n1 = n2 = n) matrices, in a stochastic setting in which (Z0)ij are iid N (0, σ2), following [9], µ̄ can be chosen as 12σ√n . This follows from
2ν ≥ 1 since ‖U‖2F = r, and somaxi ‖U∗ei‖2 ≥ rn1 .
the fact that in this setting n−1/2‖Z0‖ → 2σ almost surely; setting µ̄ in this fashion ensures that the singular value shrinkage induced by the nuclear norm regularizer ‖ · ‖∗ is greater than the largest singular value of Z0. In contrast to this σ-dependent penalty parameter, fixing S = S0, √ PCP formulation (1.2) requires that 0 ∈ ∂(‖L̂‖∗+µ‖L̂−Z0−L0‖F ), which translates into−µ L̂−Z0−L0‖L̂−Z0−L0‖F ∈ ∂‖L̂‖∗. With the hope that L̂ ≈ L0, and by the subdifferential formulation3, we have µ ‖Z0‖‖Z0‖F ≈ 1. The concentration stated above then gives an intuitive choice for µ ≈ nσ
2σ √ n , or µ = c0
√ n for
some c0 > 0. The magic of √ PCP is that by using the Frobenius norm instead of its square, the objective function becomes homogeneous, i.e. the gradient of the penalty term at the ground truthX0 becomes σ independent, making a universal penalty parameter possible.

heading: 1.3 Relationship to the Literature
text: The problem of low-rank matrix recovery from gross sparse corruption can be considered a form of robust PCA, [6, 13], and has been studied extensively in the literature. Algorithmic theory has been developed for both convex [6, 13, 14, 15, 16], and nonconvex optimization methods [7, 8, 17, 18, 19]. While many of the aforementioned works pertain to noiseless data, a line of work has studied extensions to noisy data. [9] studied the problem of robust matrix recovery with bounded noise under the incoherence assumption, and proved a bound on the recovery error, with linear dependence on the noise level but suboptimal dependence on the matrix size. [20] studied the problem with a weaker assumption about spikiness using decomposable regularizers and restricted strong convexity (RSC), and obtained essentially optimal bounds on the reconstruction error when the noise level is large. These weaker assumptions are not sufficient to ensure exact recovery, and so when the noise standard deviation σ is small, this approach does not yield a reconstruction error proportional to the noise level. [10] formulated robust PCA as a semidefinite programming problem, which requires strong assumptions about square matrices and positive semidefiniteness of the lowrank matrix. Some other works [21, 22] further assumed partial observation of the matrix, and also derived tighter bounds on the recovery error. The recent work of [8] achieves optimal error bounds for both large and small noise, using a novel analysis that leverages an auxiliary nonconvex program. Taken together, these results give efficient and provably effective methods, whose statistical performance is nearly optimal. Compared to e.g., [8, 20], the stability guarantees provided by our theory are worse by a dimension-dependent factor. Nevertheless, all of the above works regarding robust PCA with noise, the optimization involves parameters that must be set based on the noise level distribution, and therefore challenging in actual applications. On the other hand, there has been existingwork in the literature on structured signal recovery without needing to know the noise level. Our proposed square root PCP is directly inspired by the square root Lasso [11], which proposed the idea of replacing the squared loss with its “square root” version. This allows for a choice of the parameter independent of noise level, while maintaining near-oracle performance. Later works have extended this idea to other scenarios, such as group lasso [23], SLOPE variable selection [24], elastic net [25] and matrix completion [26], etc. Notably, the work of [26] studied the matrix completion where one aims to recover a low-rank matrix from noisy linear observations, aka matrix completion. Compared with that paper, this work aims to solve a different problem where the observation also contains a sparse outlier matrix. To the best of our knowledge, this paper is the first to propose a provable algorithm for Robust PCA with noisy observation that does not require knowledge of the noise level beforehand. On the algorithmic side, interior point method and first order method are used to solve the square root Lasso in [11], while later works apply ADMM to the problem [27][28]. In our problem, the objective function can be transformed into a separable form, making ADMM a reasonable choice. We note that for large σ the error bound established in this paper is suboptimal compared with [11] and [26]. The problem of square root lasso enjoys benign properties (lower bounds on restricted eigenvalues) which we do not have in square root PCP. The paper [26] on matrix completion makes a spikiness assumption, and proves that a square root
3The subdifferential of a norm satisfies ∂‖x‖ = {z | 〈x, z〉 = ‖x‖, ‖z‖∗ ≤ 1}.
lasso-inspired formulation achieves essentially optimal estimation when the noise is large. As with robust matrix recovery, the spikiness assumption is not strong enough to imply exact recovery in the noiseless case. Compared to these works, the principal differences in this paper are (i) the problem formulation: we consider robust PCA with sparse errors, (ii) the analysis, which proceeds down different lines, and (iii) that our bounds are linear in the noise level, for both large and small noise. However, in contrast to [26], our analysis does not yield minimax optimal estimation errors; it is worse by a dimension-dependent factor. Improving this dependence is an important direction for future work.

heading: 2 Analysis
text: The proof of the main Theorem 1.1 is different from the standard approach in [11] due to a lack of the Restricted Strong Convexity property for the map (L,S) → ‖L + S −D‖F . Instead, our approach has three key ingredients:
• The result from StablePCPc (Theorem 2.1) shows a recovery error ‖(L̂, Ŝ)− (L0,S0)‖F which depends linearly on the parameter δ. • The intimate connection between √ PCP formulation (1.2) and StablePCPc formulation
(1.5) can help translate the above solution property to √ PCP (Lemma 2.2).
• The powerful dual certificate construction proposed in [6] (restated in Lemma 2.3) can be used as an approximate subgradient to bound the regularizer at X̂ .
The proof of the main theorem has two steps. First, it uses the optimality condition and the subgradient to provide an upper and an lower bound for the regularization terms at X̂ . Second, the result in Theorem 2.1 is translated into the √ PCP setting, and together with the bounds obtained above, we get the desired result. The proof is given in the supplementary material, and below we provide three ingredients. First, we state the main theorem for StablePCPc problem:
Theorem 2.1 (Theorem 2 in [9]) Under Assumptions 1.1 and 1.2, assuming further that r ≤ ρ′rn2ν −1(log n1) −2 andm ≤ ρ′sn1n2 where ρ′r, ρ′s are some positive constants, there is a numerical constant c′ such that with probability at least 1− c′n−101 , for any Z0 with ‖Z0‖F ≤ δ, the solution X̂ = (L̂, Ŝ) to the StablePCPc problem 1.5 with λ = 1/ √ n1 satisfies
‖X̂ −X0‖F ≤ √ 320n1n2 + 4 · δ. (2.1)
Note that choosing δ = ‖Z0‖F allows a reconstruction error that is O( √ n1n2‖Z0‖F ). In the
case when Z0 = 0, StablePCPc recovers the matrices exactly: X̂ = X0. This is in agreement with the result in [6]. The next lemma connects the two formulations √ PCP and StablePCPc and the proof is provided in the supplementary material:
Lemma 2.2 Consider the √ PCP problem parameterized by µ and denote the result as
L̂root(µ), Ŝroot(µ), as well as the StablePCPc formulation parameterized by δ and denote the result as L̂stable(δ), Ŝstable(δ). Define δ(µ) = ‖D − L̂root(µ)− Ŝroot(µ)‖F , then
L̂stable(δ(µ)), Ŝstable(δ(µ)) = L̂root(µ), Ŝroot(µ). (2.2)
Lastly, we show an adapted dual certificate construction:
Lemma 2.3 (Adapted from [6]) Under Assumptions 1.1 and 1.2, assume that r,m satisfies
r ≤ ρrn2ν−1(log n1)−2, m ≤ ρsn1n2, (2.3)
where ρr ≤ 1/10, ρs are some positive constants. Then there is a numerical constant c such that with probability at least 1− cn−101 , there exists W ,F ,H such that
UV ∗ + W = λ(sign(S0) + F + PΩH), (2.4)
whereW ∈ T⊥, ‖W ‖ ≤ 12 , PΩF = 0, ‖F ‖∞ ≤ 1 2 , and ‖PΩH‖F ≤ 1 260 √ 2 .
The dual construction in [6] satisfies ‖PΩH‖F ≤ 14 . However, the proof for Lemma 2.8(b) indicates that ‖PΩH‖F ≤ √ r n21 ≤ √ r/n2 n1.51 and we only need to make sure that √ r/n2 n1.51 ≤ 1 260 √ 2 . This is a very mild condition, especially when it comes to the high dimensional real data (such as video). If we require that r ≤ n2/10, then problems of reasonably large dimension suffice, say, n1 ≥ 120. And in the extreme case, we can set ρr ≤ 1/(260 √ 2)2.
3 Solving √ PCP with ADMM
Different from [9] where StablePCPu (1.6) is solved viaAccelerated Proximal Gradientmethod, we solve √ PCP (and StablePCPu) viaADMM-splitting since the objective is not differentiable. To avoid multi-block ADMM which is not guaranteed to converge [29], we define variables X∗1 = (L ∗ 1,S ∗ 1 ,Z ∗), X∗2 = (L ∗ 2,S ∗ 2), and reformulate problem (1.2) as:
min X1,X2 f(X1) := ‖L1‖∗ + λ‖S1‖1 + µ ‖Z‖F (3.1)
s.t. X1 + [−I 0 0 −I I I ] X2 = [ 0 0 D ] .
The problem (3.1) can be separated into 2 blocks nicely (X1 and X2), which guarantees convergence of ADMM (under additional mild conditions)[12]. Define dual variables Y ∗ = (Y ∗1 ,Y ∗2 ,Y ∗3 ), the Lagrangian can be written as
Lρ(X1,X2,Y ) = ‖L1‖∗ + λ‖S1‖1 + µ ‖Z‖F + 〈L1 −L2,Y1〉+ ρ
2 ‖L1 −L2‖2F + 〈S1 − S2,Y2〉
+ ρ
2 ‖S1 − S2‖2F + 〈L2 + S2 + Z −D,Y3〉+
ρ 2 ‖L2 + S2 + Z −D‖2F .
We present the update rules4 as well as the stopping criteria adapted from [12] in Algorithm 1 and helper() function in the supplementary material. The stopping criteria takes into account the primal and the dual feasibility conditions, and the algorithm stops when the tolerances set using an absolute and relative criterion are reached.
If we modify the update of Z in Algorithm 1 to Z ← ( D −L2 − S2 − 1ρY3 ) /(1 + µ/ρ), we get ADMM for StablePCPu (1.6). Algorithm 1 Algorithm for √ PCP
Input: D ∈ Rn1×n2 , λ, µ. Output: L,S ∈ Rn1×n2 .
# Tolerance levels, max iterations abs ← 10−6, rel ← 10−6, N ← 5000 # Initialization L1,L2,S1,S2,Z,Y1,Y2,Y3 ← 0n1×n2 ρ← 0.1 for i = 1, i ≤ N, i+ + do
# Save old values temporarily (L′2,S ′ 2)← (L2,S2) # ADMM updates L1 ← prox 1
ρ‖·‖∗ ( L2 − 1ρY1 ) S1 ← proxλ
ρ ‖·‖1
( S2 − 1ρY2 )
Z ← proxµ ρ ‖·‖F ( D −L2 − S2 − 1ρY3 ) L2 ← (D−Z+2L1−S1+ 1ρ (2Y1−Y2−Y3)) 3 S2 ← (D−Z+2S1−L1+ 1ρ (2Y2−Y1−Y3)) 3 Y1 ← Y1 + ρ(L1 −L2) Y2 ← Y2 + ρ(S1 − S2) Y3 ← Y3 + ρ(L2 + S2 + Z −D) # Update ρ and check convergence ρ, ifConverge← helper() if ifConverge then
break end if
end for (L,S)← ((L1 + L2)/2, (S1 + S2)/2) return L,S
4Recall that proxγ‖·‖∗ (Z) = ∑ imax(λi − γ, 0)uiv ∗ i , where Z = ∑ i λiuiv ∗ i is the SVD,
[proxγ‖·‖1 (Z)]i,j = max(|Zi,j | − γ, 0) · sign(Zi,j), and proxγ‖·‖F (Z) = max(‖Z‖F − γ, 0) Z ‖Z‖F .

heading: 4 Experiments
text: To show the effectiveness of our new formulation, we test √ PCP on simulated data as well as real-world video datasets. The experiments suggest that our error bound in Theorem 1.1 has a correct dependency on the noise level of Z0, but loses a factor of n (the dimension of the problem). In addition, the solutions produced by √ PCP with our proposed noiseindependent µ and StablePCPu with the noise-dependent µ often look very similar to each other. Moreover, experiments on real-world datasets with natural noise also show the denoising effect of √ PCP , making √ PCP a practical approach with good performance in this robust noisy low-rank matrix recovery setting. Additional experiments of √ PCP on simulated data with varying µ also suggest that √ n2/2 can provide performance (recovery error) close to the optimal µ, justifying our proposed choice of µ = √ n2/2.

heading: 4.1 Simulations with Varying Noise Levels and Dimension
text: In this set of experiments, we are interested in how our error bound in Theorem 1.1 compare with the actual reconstruction error. We simulate (L0,S0,Z0) with varying noise levels of Z0 and problem dimension n1, n2. To simulate L0 ∈ Rn1×n2 of rank r, we generate U ∈ Rn1×r,V ∈ Rn2×r as the unnormalized singular vectors such that U ,V are entrywise i.i.d. N (0, 1/n1) andN (0, 1/n2) respectively and let L0 = UV ∗. For S0, we let P [(i, j) ∈ Ω] = ρS and for (i, j) in support Ω, (S0)(i,j) ∈ {0.05,−0.05} with equal probability. For the noise Z0, we generate it as entrywise i.i.d. N (0, σ2). In addition, in the experiments we take n1 = n2 = n, so we choose λ = 1/ √ n, µstable = 1/(2σ √ n) (the noise level σ is known), and µroot = √ n/2. Theoretical analysis in Theorem
1.1 and 2.1 shows that with these parameters, ‖X̂ −X0‖F = O(n‖Z0‖F ).
To test the dependency of the error on σ, we take n = 200, r = 10, and so ‖L0‖2F ≈ r = 10. For the outlier S0, we take ρS = 0.1, so ‖S0‖2F ≈ 0.052n2ρS = 10. For the noise Z0, we take σ ∈ {0, 0.001, . . . , 0.015}5, so ‖Z0‖2F ≈ σ2n2 ∈ [0, 9]. For each σ in the given set, we randomly generate 20 ground truth (L0,S0,Z0) triplets and run √ PCP and StablePCPu on
them. We use the root-mean-squared (RMS) error defined as ( 120 ∑20 k=1 ‖L̂(k) − L0‖2F )1/2
and ( 120 ∑20 k=1 ‖Ŝ(k)−S0‖2F )1/2 for evaluation. In Figure 1(a) we show the RMS error over 20 trials for the low-rank and the sparse. It is clear from the plot that ‖L̂−L0‖F and ‖Ŝ−S0‖F are O(σ) for both √ PCP and StablePCPu, which confirms that the reconstruction error is linear in the noise level σ. We also notice that the recovery error in Figures 1(a) and 1(c) is linear in the noise level for small σ, but exhibits a sublinear behavior for larger σ. This behavior reflects a general phenomenon in recovery/denoising using structured models (sparse, low-rank, etc.): the minimax noise sensitivity η = supσ>0 1σE[‖x̂− x0‖] is obtained as σ → 0. This means that for small σ, we expect a linear trend with slope η, while for larger σ, the dependence can be sublinear. This behavior has a general geometric explanation. For simplicity we sketch how this plays out in a simpler norm denoising problem, in which the target is to recover a structured signal x0, and we observe y = x0 + σz. For simplicity, assume that we know that ‖x0‖1 ≤ τ , and solve min‖x‖1≤τ ‖x − y‖2. For small σ, the estimation error x̂ − x0 is simply the projection of the noise σz onto the descent cone of the norm ball {‖x‖1 ≤ τ} at x0; its size is linear in σ. For larger σ, there is additional denoising due to the fact that the L1 ball is smaller than the descent cone at x0 — this leads to the behavior observed here. To test the dependency of the error on the problem dimension, we vary n ∈ {200, 300, . . . , 1000} and take r = 0.1n. We keep the setting for S0, and take σ = 0.01 as the noise level for Z0. Figure 1(b) shows the RMS error. Note that for fixed σ, ‖Z0‖F ∼ nσ, so the results in Theorem 1.1 and Theorem 2.1 bound the reconstruction error as O(n2).
5When σ = 0 and µstable = +∞, StablePCPu is equivalent to StablePCPc with δ = 0.
However, the analysis provides only a loose error bound. As can be seen from this set of experiment, the error is closer to O(n). We provide experiments with different distributions of the noise in the appendix.

heading: 4.2 Real Data with Added Noise: Surveillance Video
text: Many imaging datasets can be modeled as the sum of a low-rank matrix L0, a sparse outlier S0, and noise Z0. For instance, video data often consists of an almost fixed background which can be seen as low-rank, and a foreground (such as people) that only occupies a small fraction of the image pixels for a short amount of time, which can be considered as sparse. Thus, videos can naturally fit into our robust PCA framework. In this set of experiments, we use √ PCP and StablePCPu to separate the background and the foreground for surveillance video data. We assume that the original video is noiseless, and manually add noise Z0 that is entrywise i.i.d. N (0, σ2) to test the dependency of the reconstruction error on the noise level σ. We use the “hall dataset” in [30], a 200-frame video of a hall that has people walking around. Each frame has resolution 144× 176, and is flattened as one column of the noiseless observation matrixD, so we have n1 = 144× 176 and n2 = 200. Each pixel is represented by a number in [0, 255], and the mean value among all pixels is 150.3295, with standard deviation 45.7438, and median 155.0000. For the added noise, we choose σ ∈ {0, 30, 60, 90, 120}, and denote the recovered matrices as X̂(σ)root/stable. In addition, we let X0 = 1 2 (X̂ (0) root + X̂ (0) stable) be the ground truth, and evaluate the error using ‖L̂(σ)root/stable −L0‖F and ‖Ŝ (σ) root/stable − S0‖F . We take λ = 1/ √ n1,
µroot = √ n2/2 and µstable = 1σ(√n1+√n2) following the same intuition as in Section 1.
6 We run the experiments on a laptop with 2.3 GHz Dual-Core Intel Core i5, and set the maximal iteration of ADMM to be 5000. All of these experiments on real datasets end within 1 hour. For full details, please see the supplementary material. In Figure 1(c), we show the reconstruction error with varying noise levels. It can be seen that the error is indeed linear in σ, as predicted by our analysis. In Figures 2, we present the first frame (i.e. the first column) of the original video (with noise σ = 0, 30), and the √ PCP recovered low-rank and sparse matrices. Although the added noise blurs the videos, our√ PCP is still stable and successfully decompose the background and the foreground.

heading: 4.3 Real Data with Natural Noise: Low Light Video
text: Low light videos are known to have very large observation noise due to limited photon counts. In this experiment, we apply our √ PCP to the Dark RawVideo (DRV) dataset in [32] (under MIT License) for foreground background separation and denoising. This dataset of 6Recall that E[‖Z0‖] ≤ σ( √ n1 + √ n2) for rectangular matrices, e.g. from [31]
RGB videos, approximately 110 frames each, 3672× 5496 in resolution, was collected at low light settings, so the signal-to-noise ratio (SNR) is extremely low (negative if measured in dB)[32]. For the experiments, we choose videoM0001 (basketball player), M0004 (toy windmill), and M0009 (billiard table) from DRV. As preprocessing, we convert the RGB videos to grayscale using rgb2gray() in Matlab, and crop and downsample each frame to reduce data size. The final resolution is 322× 440 for M0001, 294× 440 for M0004, and 306× 458 for M0009. We apply √ PCP with λ = 1/√n1, and µ = √ n2/2 to these 3 videos, and present the results for frame 30 in Figure 3. The denoising effect of √ PCP can be seen by comparing D with
L̂+Ŝ. In addition, L̂ recovers the background prettywell, Ŝ captures themoving foreground but is still mixed with noise, which we believe is due to the extremely low SNR.

heading: 4.4 Real Data with Natural Noise: Optical Coherence Tomography
text: In medical imaging, Optical Coherence Tomography can be used for micro-scale resolution, quick scanning of biological phenomenon [33]. These scans of the same scene over time, called time-lapse B-scan, are often noisy, but fit into our low rank/sparse model. In this experiment, we apply √ PCP to the time-lapse B-scans (250 frames of resolution 300× 150) of human trachea samples containing motile cilia (demo dataset of [33] under CC0 License). We present the recovered frame 50 and 100 in Figure 4. As expected, L̂ captures the static background, and Ŝ captures the motion of cilia.

heading: 4.5 Optimal Choice of µ
text: Our main result Theorem 1.1 suggests a tuning-free µ = √ n2/2. Here, we investigate experimentally if this choice of µ is optimal. We vary the problem dimensions n1 and n2, the rank-dimension ratio ρL := r/n (n = n1 = n2), and the noise standard deviation σ. For each choice of parameters (n1, n2, ρL, σ), we generate 10 pairs of (L0,S0,Z0) using the same method as in Section 4.1, run √ PCP with λ = 1/√n1 and µ = cµ0, where µ0 = √ n2 and c is a varying coefficient (c = 1/ √
2 ≈ 0.71 corresponds to our proposed value). Settings of all parameters are included in the supplementary material.
We use the 10-average of ‖(L̂, Ŝ) − (L0,S0)‖F as the evaluation metric. In Figure 5, we show the heatmaps of this metric relative to the optimal µ = cµ0 among all tested c, i.e. ηrel(µ) =
‖(L̂(µ),Ŝ(µ))−(L0,S0)‖F minµ′=cµ0 ‖(L̂(µ ′),Ŝ(µ′))−(L0,S0)‖F , so the optimal µ has value 1 in each row of the
heatmaps. From Figure 5, we see that varying n1, n2 has little effect on the optimal choice
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
coefficient
300
400
500
600
700
800
900
1000
n 1
1.375
1.413
1.427
1.434
1.438
1.435
1.432
1.425
1.191
1.22
1.23
1.233
1.236
1.236
1.236
1.235
1.069
1.094
1.103
1.106
1.107
1.107
1.107
1.106
1.005
1.021
1.028
1.03
1.032
1.032
1.032
1.031
1
1
1
1
1
1
1
1
1.066
1.037
1.024
1.018
1.014
1.012
1.011
1.01
1.15
1.144
1.11
1.091
1.08
1.073
1.069
1.067
1.15
1.177
1.191
1.204
1.213
1.196
1.185
1.178
1.15
1.177
1.191
1.204
1.214
1.222
1.229
1.239
1.15
1.177
1.191
1.204
1.214
1.222
1.229
1.239
1.633
1.689
1.704
1.702
1.69
1.674
1.657
1.638
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
(b) vary n1
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
coefficient
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
L
1.496
1.537
1.275
1.373
1.467
1.526
1.532
1.519
1.499
1.479
1.456
1.436
1.126
1.19
1.257
1.307
1.338
1.345
1.344
1.34
1.328
1.318
1.035
1.07
1.114
1.146
1.167
1.182
1.197
1.201
1.201
1.198
1
1.005
1.028
1.045
1.056
1.063
1.077
1.084
1.089
1.092
1.03
1
1
1
1
1.001
1.008
1.012
1.017
1.019
1.137
1.065
1.038
1.02
1.008
1
1
1
1
1
1.266
1.149
1.097
1.062
1.038
1.023
1.019
1.015
1.013
1.012
1.266
1.149
1.097
1.062
1.038
1.023
1.019
1.015
1.013
1.012
1.266
1.149
1.097
1.062
1.038
1.023
1.019
1.015
1.013
1.012
1.266
1.149
1.097
1.062
1.038
1.023
1.019
1.015
1.013
1.012
1.628
1.727
1.749
1.721
1.679
1.638
1.604
1.568
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
(c) vary ρL
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
coefficient
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
1.631
1.368
1.231
1.149
1.097
1.054
1.03
1.009
1
1.531
1.374
1.226
1.138
1.083
1.05
1.023
1.01
1
1.001
1.251
1.189
1.117
1.064
1.03
1.013
1
1
1.002
1.016
1.096
1.068
1.039
1.015
1
1
1.003
1.018
1.035
1.065
1.016
1.004
1
1
1.006
1.025
1.047
1.081
1.117
1.166
1
1
1.013
1.036
1.064
1.106
1.15
1.206
1.266
1.34
1.055
1.066
1.095
1.141
1.193
1.26
1.331
1.413
1.5
1.602
1.124
1.151
1.198
1.265
1.34
1.43
1.524
1.632
1.742
1.124
1.151
1.198
1.265
1.34
1.43
1.524
1.632
1.742
1.124
1.151
1.198
1.265
1.34
1.43
1.524
1.632
1.742
1.124
1.151
1.198
1.265
1.34
1.43
1.524
1.632
1.742
2.128
1.873 1.873 1.873 1.873
1
1.2
1.4
1.6
1.8
2
(d) vary σ
of µ, which is approximately between 0.7√n2 and 0.75 √ n2, close to our
√ n2/2. However,
decreasing ρL or increasing σ suggests a smaller value of optimalµ. Thismakes sense because with higher level of noise or smaller rank (thus smaller norm ‖L0‖F ), the SNR is smaller, so we should put smaller penalty on ‖L+S−D‖F . Nevertheless, in all these settings, choosing µ = √ n2/2 still gives satisfying results, as the recovery errors for µ = 0.7 √ n2 are close to the optimal performance: the error ratios are below 1.2. From these results, we believe that while µ = √ n2/2 may not be the optimal choice with respect to the recovery error, it can provide performance close to the optimal, and is therefore a very effective choice.

heading: 5 Conclusion
text: In this work, we propose √ PCP , a convex optimization approach for noisy robust low-rank matrix recovery. The benefit of our approach as compared to previous methods such as stable PCP is that it enables tuning-free recovery of low-rank matrices: theoretical analysis and simulations show that a single universal penalty parameter yields stable recovery at any noise standard deviation. Real video data experiments show suggest that many real life models fit into this low-rank plus sparse setting, and √ PCP (as well as stable PCP) does a good job in denosing and recovering the patterns of interest. The presented experiments suggest the potential for both positive and negative societal impacts: visual surveillance can be abused, leading to significant negative impacts; at the same time, the denoising and foreground/background separation ability of √ PCP can help improve the quality of noisy data in biomedical and scientific research (e.g. medical imaging), and people’s life (e.g. low light video).

heading: Acknowledgement
text: Wewould like to thankMarianthi-Anna Kioumourtzoglou, Jeff Goldsmith, Elizabeth Gibson, Rachel Tao and Lawrence Chillrud for many helpful discussions around matrix modeling of environmental data and the need for tuning-free solutions. This work was partially funded by the National Institute of Environmental Health Sciences (NIEHS) grant R01 ES028805. We also thank Christine Hendon for helpful pointers regarding optical coherence tomography data. Jingkai Yan also gratefully acknowledges support from the Wei Family Foundation.

Review summaries:
Low-rank matrix recovery in the presence of sparse corruptions and random noise, under incoherence assumptions on the low-rank and on the sparse components, have been considered. Following an existing line of work, a tuning-free optimization programs has been proposed for demixing and denoising. The analysis, utilizing existing techniques, establishes a bound on the reconstruction error. The bound has been observed to be loose by a dimension-dependent factor (according to the experiments and compared to other estimators for this problem). The estimator is a convex optimization problem which authors solve via ADMM.  Numerical experiments support the usefulness of the proposed approach in a set of synthetic and real data scenarios. 

This work combines stable Principal Component Pursuit and Square Root Lasso and proposes Square Root Principal Component Pursuit for low-rank + sparse recovery problem. Compared to stable PCP, Square Root PCP has the advantage that the choice of the regularization parameter is noise-independent, so it is easier to be deployed in applications where the noise level is unknown. Theoretical analysis and experiments confirm the effectiveness of the proposed method.