Year: 2020

Abstract:
Recently, the `-norm maximization has been proposed to solve the sparse dictionary learning (SDL) problem. The simple MSP (matching, stretching, and projection) algorithm proposed by Zhai et al. (2019a) has shown to be surprisingly efficient and effective. This paper aims to better understand this algorithm from its strong geometric and statistical connections with the classic PCA and ICA, as well as their associated fixed-point style algorithms. Such connections provide a unified way of viewing problems that pursue principal, independent, or sparse components of high-dimensional data. Our studies reveal additional good properties of `-maximization: not only is the MSP algorithm for sparse coding insensitive to small noise, but it is also robust to outliers and resilient to sparse corruptions. We provide statistical justification for such inherently nice properties. To corroborate the theoretical analysis, we also provide extensive and compelling experimental evidence with both synthetic data and real images.

heading: 1 INTRODUCTION
text: The explosion of massive amounts of high-dimensional data has become the modern-day norm for a large number of scientific and engineering disciplines and hence presents a daunting challenge for both computation and learning. Rising to this challenge, sparse dictionary learning (SDL) provides a potent framework in representation learning that exploits the blessing of dimensionality: real data tends to lie in or near some low-dimensional subspaces or manifolds, even though the ambient dimension is often extremely large (e.g. the number of raw pixels in an image). More specifically, SDL (Olshausen & Field (1997); Mairal et al. (2008; 2012; 2014); Spielman et al. (2012); Sun et al. (2015); Bai et al. (2018); Qu et al. (2019)) concerns the problem of learning a compact, sparse representation from raw, unlabelled data: given a data matrix Y = [y1,y2, . . . ,yp] ∈ Rn×p that contains p n-dimensional samples, one aims to find a linear transformation (i.e. a dictionary) D ∈ Rn×m and an associated maximally sparse representation X = [x1,x2, . . . ,xp] ∈ Rm×p that satisfies Y = DX. (1) As the data matrix Y can represent a variety of signals (e.g. images, audios, languages, and genetics etc) in practical applications, SDL provides a versatile structure-seeking formulation that has found widespread applications in computational neuroscience, image processing, computer vision, and machine learning at large (Olshausen & Field, 1996; 1997; Argyriou et al., 2008; Ranzato et al., 2007; Elad & Aharon, 2006; Wright et al., 2008; Yang et al., 2010; Zhang et al., 2013; Mairal et al., 2014; Zhang et al., 2014; 2019).
Related Works. Motivated by this practical significance, there has been a growing surge of interest recently (e.g. Rambhatla et al. (2019); Bai et al. (2018); Gilboa et al. (2018); Nguyen et al. (2018); Chatterji & Bartlett (2017); Mensch et al. (2016)) that aims to tackle SDL. In attempts to recover the sparse signals X , these existing work adopt an `0- or `1-penalty function to promote the underlying sparsity and give various optimization algorithms for the resulting objectives (some of those are heuristics while a few others have theoretical convergence guarantees). Although these penalty
∗Work done when interning at ByteDance. YZ would like to thank Chong You at UC Berkeley for meaningful discussions. Correspondence to: Yuexiang Zhai, ysz@berkeley.edu or Yi Ma, yima@eecs.berkeley.edu. Codes are available at https://github.com/hermish/ZMZM-ICLR-2020.
functions are indeed sparsity-promoting, the resulting optimization problems must be solved one row at a time, hence resulting as many optimization problems as the ambient dimension n. Consequently, `0- or `1-based objectives result only in local methods (i.e. cannot yield the entire solution at once) and hence entail prohibitive computational burden. Another prominent approach in SDL is Sumof-Squares (SOS), proposed by and articulated in a series of recent work (Barak et al. (2015); Ma et al. (2016); Schramm & Steurer (2017)). The key idea there is to utilize the properties of higher order SOS polynomials to correctly recover one column of the dictionary at a time, for which there are m columns in total. However, the computational complexity of these recovery methods are quasi-polynomial, thus again resulting in large computational expense.
Very recently, in the complete dictionary learning1 setting, a novel global approach has been suggested in Zhai et al. (2019a;b) that presents a formulation which can efficiently recover the sparse signal matrix X once and for all. In particular, Zhai et al. (2019b) has shown that if the generative model for Y = DoXo ∈ Rn×p satisfies that Do ∈ O(n;R) is orthonormal and Xo ∈ Rn×p is Bernoulli-Gaussian sparse,2 then maximizing the `4-norm3 of AY over O(n;R):
max A
1 4 ‖AY ‖44 subject to A ∈ O(n;R) (or AA ∗ = I), (2)
is able to find the ground truth dictionary Do up to an arbitrary signed permutation. Moreover, Zhai et al. (2019b) has proposed the simple “Matching, Stretching, and Projection” (MSP) algorithm, which is shown to be experimentally efficient and effective, for solving the program in equation 2:
MSP: At+1 = PO(n;R) [ (AtY ) ◦3Y ∗ ] = UtV ∗ t , (3)
where UtV ∗t are from the singular value decomposition: UtΣtV ∗ t = SVD[(AtY ) ◦3Y ∗].
In this paper, we here give an alternative (arguably simpler and more revealing) derivation of the MSP algorithm (3). Consider the Lagrangian formulation of the constrained optimization problem given in equation 2. The necessary condition of critical points,∇A 14 ‖AY ‖ 4 4 = ∇A 〈Λ,AA∗ − I〉 for some Lagrangian multipliers Λ, implies:
(AY )◦3Y ∗ = (Λ + Λ∗)A. (4)
As the optimization is over the orthogonal group O(n;R), restricting the condition in equation 4 onto the orthogonal group yields a necessary condition for any critical point A:4
PO(n;R) [ (AY )◦3Y ∗ ] = A. (5)
Hence the critical point A can be viewed as a “fixed point” of the map: PO(n;R) [ ((·)Y )◦3Y ∗ ] from O(n;R) to itself. The MSP algorithm in equation 3 is to find the fixed point(s) of this map.
Notice that the orthonormal constraint A ∈ O(n;R) in equation 2 can be viewed as enforcing the orthogonality of n unit vectors simultaneously. So, more flexibly and generally, one may choose to compute any k, for 1 ≤ k ≤ n, leading orthonormal bases of Do by solving the program:
max W
1 4 ‖W ∗Y ‖44 subject to W ∈ St(n, k;R) ⊂ R n×k, (6)
where St(n, k;R) is the Stiefel manifold.5 The orthogonal group O(n;R) and the unit sphere Sn−1 can be viewed as two special cases of the Stiefel manifold St(n, k;R), with k = n and k = 1, respectively. In some specific tasks such as dictionary learning and blind deconvolution, optimization over the unit sphere has been widely practiced, such as in Sun et al. (2015); Bai et al. (2018); Zhang et al. (2018); Kuo et al. (2019). The more general setting of maximizing a convex function over any compact set also has been studied by Journée et al. (2010) in the context of sparse PCA, which has provided convergence guarantees for this class of programs.
1Complete dictionary learning requires the learned dictionary D in equation 1 to be square and invertible. 2 Each entry xi,j of X can be represented as the product of a Bernoulli variable and a normal Gaussian variable: xi,j = Ωi,jVi,j , where Ωi,j ∼iid Ber(θ) and Vi,j ∼iid N (0, 1), similar for vectors or scalars. This is the standard setting adopted in Spielman et al. (2012); Sun et al. (2015); Bai et al. (2018).
3We abuse the notation a bit, by denoting ‖·‖44 as the sum of element-wise 4 th power of all entries of a vector and matrix, that is, ∀a ∈ Rn, ‖a‖44 = ∑n i=1 a 4 i and ∀A ∈ Rn×m, ‖A‖44 = ∑ i,j a 4 i,j .
4For any symmetric matrix S ∈ Rn×n and an orthogonal matrix A ∈ O(n;R), the projection of SA onto the orthogonal group is A: PO(n;R) [SA] = A, one may see Absil & Malick (2012) for details.
5For any 1 ≤ k ≤ n, St(n, k;R) .= {W ∈ Rn×k : W ∗W = Ik}.
Our Contributions. Our contributions are twofold. First, by taking a suitable analytical angle, we reveal novel geometric and statistical connections between PCA, ICA and the `4-norm maximization based SDL. We then show that algorithm-wise, the fixed-point type MSP algorithm for `4-norm maximization has the same nature as the classic power-iteration method for PCA (Jolliffe, 2011) and the FastICA algorithm for ICA (Hyvärinen & Oja, 1997). This interpretation gives a unified view for problems that pursue principal, independent, or sparse components from high-dimensional data and enriches our understanding of low-dimensional structure recovery frameworks, classical and new, at both formulation and algorithmic fronts.
Second, and more importantly from a practical perspective, we examine how MSP performs under a variety of more realistic conditions, when the measurements Y could be contaminated with noise, outliers, or sparse corruptions. We show that, similar to PCA, `4-norm maximization and the MSP algorithm are inherently stable to small noise. Somewhat surprisingly though, unlike PCA, the MSP algorithm is further robust to outliers and resilient to sparse gross errors! We provide characterizations of these desirable properties of MSP. The claims are further corroborated with extensive experiments on both synthetic data and real images. Taken as a whole, our results contribute to the broad landscape of dictionary learning by affirming that `4-maximization based SDL and the corresponding global algorithm MSP provide a valuable toolkit to the existing literature.

heading: 2 SDL VERSUS PCA AND ICA

heading: 2.1 PURSUIT OF PRINCIPAL, INDEPENDENT, OR SPARSE COMPONENTS
text: Relation with the Geometric Interpretation of PCA. For a data matrix Y ∈ Rn×p, Principal Component Analysis (PCA), aiming to find the top (top k) left singular vector (vectors) of Y ,
max W
1 2 ‖W ∗Y ‖2F subject to W ∈ St(n, k;R), (7)
can be understood as finding a direction (a k-dimensional subspace) in row(Y ) in which Y has the largest `2-norm (Frobenius norm). For instance, finding the direction with the largest `2-norm over the unit sphere can be viewed as calculating the spectral norm (or the largest singular value), of matrix Y . In comparison, we may view equation 6:
max W
1 4 ‖W ∗Y ‖44 subject to W ∈ St(n, k;R)
as finding a direction, or a k-dimensional subspace, in row(Y ) where the projection of Y has the largest `4-norm. For instance, finding the direction with the largest `4-norm over the unit sphere (equation 6) can be viewed as calculating the induced ‖·‖2,4 norm of matrix Y : ‖Y ‖2,4 . = maxa∈Sn−1 ‖a∗Y ‖4 .
Relation with the Statistical Interpretation of PCA. View each column yj , j ∈ [p] of data matrix Y ∈ Rn×p as an n dimensional random vector that is i.i.d. drawn from a distribution of random variable y and let Yc denote the centered Y : Yc . = Y [I − 1p11
∗], where 1 ∈ Rp is a vector of all 1’s. Then, finding the top k principal components of Yc: maxW∈St(n,k;R) 12 ‖W
∗Yc‖22 is to find k uncorrelated projections of y ∈ Rn that has the top k sample variance (i.e. 2nd order moment) (Jolliffe, 2011; Helwig, 2017). Similar to PCA, the `4-norm maximization of centered data matrix Yc: maxW∈St(n,k;R) 14 ‖W
∗Yc‖44 can be viewed as finding k uncorrelated projections of y that have the top k sample 4th order moment, whose statistical meaning is better revealed below.
Relation with ICA and Nonnormality. The `4-norm maximization over the Stiefel manifold is strongly related to finding the maximal or minimal kurtosis in Independent Component Analysis (ICA) (Hyvärinen & Oja, 1997; 2000): In order to identify one component of a given random vector y ∈ Rn, ICA aims to find a unit vector (a direction) w ∈ Sn−1 that maximizes or minimizes the kurtosis of w∗y, defined as:
kurt(w∗y) = E (w∗y)4 − 3 ‖w‖42 . (8) Kurtosis is widely used for evaluating the nonnormality of a random variable; see DeCarlo (1997); Hyvärinen & Oja (1997; 2000). According to Huber (1985), the nonnormality of data carries “abnormal” hence interesting information in real data for many applications (e.g. Lee et al. (2003); Cain
et al. (2017)). Thus, extracting the 4th order moment helps understand such statistics of real datasets (Hyvärinen et al., 2009) and even their topology (Carlsson, 2009). One may also find that the `4- maximization based dictionary learning formulation is similar to maximizing kurtosis (equation 8) with spherical constraint w ∈ Sn−1, in fact, these two formulations are exactly the same if one only wants to find one column of the dictionary. Intuitively, such coincidence occurs due to the fact that maximizing `4-norm or kurtosis have the same effect— they both promote the “spikiness” (Zhang et al. (2018); Li & Bresler (2018)) (or “peak” (DeCarlo, 1997)) of a distribution. More rigorous analysis regarding the similarity between `4-based dictionary learning and ICA is beyond the scope of this paper and we leave them for future research.

heading: 2.2 FIXED-POINT STYLE ALGORITHMS
text: In optimization, the `4-norm maximization in equation 6 over the Stiefel manifold St(k, n;R) is a special type of nonconvex optimization problem – convex maximization over a compact set. Although Journée et al. (2010); Zhai et al. (2019b) have shown that the MSP algorithm is guaranteed to find critical points, the experiments in Zhai et al. (2019b) suggest that the MSP algorithm finds global maxima of the `4-norm efficiently and effectively. For better understanding, in this section we illustrate some striking similarities between the MSP algorithm and the “power-iteration” type algorithms for solving PCA as well as ICA.
Fixed-point Perspective of Power Iteration. For a general data matrix Y ∈ Rn×p, finding the top singular value of Y is equivalent to solving the following optimization problem:
max w
ϕ(w) . =
1 2 ‖w∗Y ‖22 subject to w ∈ S n−1. (9)
For this constrained optimization, the Lagrangian multiplier method gives the necessary condition: ∇wϕ(w) = Y Y ∗w = λw, similar to equation 4. If we restrict this condition onto the sphere, we obtain the fixed point condition w = PSn−1 [∇wϕ(w)]. The classic power-iteration method
wt+1 = PSn−1 [∇wϕ(wt)] = Y Y ∗wt ‖Y Y ∗wt‖2 , (10)
is precisely computing this fixed point, which is arguably the most efficient and widely used algorithm to solve equation 9, for PCA (or computing SVD of Y ).
Fixed-point Perspective of FastICA. In order to maximize (or minimize) the kurtosis over Sn−1,
max w
ψ(w) . =
1 4 kurt[w∗y] = 1 4 E [w∗y]4 − 3 4 ‖w‖42 subject to w ∈ S n−1, (11)
Hyvärinen & Oja (1997) has proposed the following fixed-point type iteration: wt+1 = PSn−1 [∇wψ(wt)] = E [ y (y∗wt) 3]− 3 ‖wt‖22 wt∥∥E [y (y∗wt)3]− 3 ‖wt‖22 wt∥∥2 , (12) which enjoys cubic (at least quadratic) rate of convergence, under the ICA model assumption.
Fixed-point Perspective of MSP. For the `4-norm maximization program:
max W
φ(W ) . =
1 4 ‖W ∗Y ‖44 subject to W ∈ St(n, k;R),
through a similar derivation to that in Section 1, one can show that the MSP iteration in equation 5 for the orthogonal group generalizes to the Stiefel manifold case as:
Wt+1 = PSt(n,k;R) [∇Wφ(Wt)] = UtV ∗t , (13)
where UtΣtV ∗t = SVD[Y (Y ∗Wt) ◦3]. The above iteration has the same nature as the power iteration in equation 10 and equation 12, since they all solve a fixed-point type problem, by projecting gradient of the objective function ∇ϕ(·),∇φ(·),∇ψ(·) onto the constraint manifold Sn−1 and St(n, k;R), respectively. Table 1 summarizes these striking similarities.
3 STABILITY AND ROBUSTNESS OF `4-NORM MAXIMIZATION
Even though the MSP algorithm for `4-norm maximization is similar to power-iteration for PCA, in real applications, PCA often requires modifications to improve its robustness (Candès et al., 2011; Xu et al., 2010; 2012). In this section, we want to examine the stability and robustness of the `4- maximization for different types of imperfect measurement models: small noise, outliers, and sparse corruptions of large magnitude.

heading: 3.1 DIFFERENT MODELS FOR IMPERFECT MEASUREMENTS
text: We adopt the same Bernoulli-Gaussian model as in prior works (Spielman et al., 2012; Sun et al., 2015; Bai et al., 2018; Zhai et al., 2019b) to test the stability and robustness of the `4-maximization framework. Assume our clean observation matrix Y ∈ Rn×p is produced by the product of a ground truth orthogonal dictionary Do and a Bernoulli-Gaussian matrix Xo ∈ Rn×p:
Y = DoXo, Do ∈ O(n;R), {Xo}i,j ∼iid BG(θ). (14)
Now let us assume we only observe different types of imperfect measurements of Y :
Noisy Measurements: YN := Y +G, where G ∈ Rn×p is matrix that satisfies gi,j ∼iid N (0, η2) and η > 0 controls the variance of the noise.
Measurements with Outliers: YO := [Y ,G′], where YO contains extra columns (G′ ∈ Rn×τp)6 generated from an independent Gaussian process g′i,j ∼iid N (0, 1). Here, τ controls the portion of the outliers w.r.t. the clean data size p.
Measurements with Sparse Corruptions: YC := Y +σB◦S, where σ > 0 controls the scale of corrupting entries,7 B ∈ Rn×p is a Bernoulli matrix so bi,j ∼iid Ber(β) with β ∈ (0, 1) controlling the ratio of the sparse corruptions, and S ∈ Rn×p has entries si,j drawn i.i.d. from a Rademacher distribution:
si,j = { 1 with probability 1/2 −1 with probability 1/2 . (15)

heading: 3.2 STATISTICAL ANALYSIS AND JUSTIFICATION
text: The analysis for the stability and robustness of the `4-norm maximization follows similar statistical analysis techniques in Zhai et al. (2019b) to show that the global maximum of
W? ∈ arg max W E ‖W ∗Y ‖44 , subject to W ∈ O(n;R) (16)
satisfies W ∗?Do ∈ SP(n).8 Note here Y denotes different imperfect measurements: noisy (YN ), with outliers (YO), and with sparse corruptions (YC). Below we provide the expectation and concentration results for ‖W ∗Y ‖44 over the data distribution. We show that E ‖W ∗Y ‖ 4 4 is largely determined by ‖WDo‖44, a quantity that indicates a “distance” from W ∗Do to SP(n). As shown in Lemma 2.3 and Lemma 2.4 in Zhai et al. (2019b), the only global maximizers of ‖W ∗Do‖44 are signed permutation matrices, and W ∗Do converges to a signed permutation matrix as ‖W ∗Do‖44 reaches its global maximum.
6In case τp is not an integer, we round τp to the closest integer. 7In our context, σ = 1 is already corruption of large magnitude, since the variance of the sparse signal is 1. 8SP(n) is the signed permutation group, a group of orthogonal matrices that only contain 0,±1.
Proposition 3.1 (Expectation of Objective with Small Noise) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G ∈ Rn×p, gi,j ∼iid N (0, η2) independent of Xo, let YN = Y + G denote the data with noise. Then the expectation of 1 np ‖W
∗YN‖44 satisfies: 1
np EXo,G ‖W ∗YN‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + Cθ,η, (17)
where Cθ,η is a constant which depends on θ and η.
Proof See Appendix A.1.
Theorem 3.2 (Concentration of Objective with Small Noise) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G ∈ Rn×p, gi,j ∼iid N (0, η2) independent of Xo, let YN = Y + G denote the input with noise, then:
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YN‖44 − E ‖W ∗YN‖44∣∣∣ ≥ δ ) < 1
p , (18)
when p = Ω ( (1 + η2)4n2 lnn/δ2 ) .
Proof See Appendix A.2.
Proposition 3.3 (Expectation of Objective with Outliers) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G′ ∈ Rn×τp, g′i,j ∼iid N (0, 1) independent of Xo, let YO = [Y ,G′] denote the data with outliers G′. Then the expectation of 1 np ‖W ∗YO‖44 satisfies:
1
np EXo,G′ ‖W ∗YO‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + Cθ,τ , (19)
where Cθ,τ is a constant depends on θ, τ .
Proof See Appendix A.3
Theorem 3.4 (Concentration of Objective with Outliers) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G′ ∈ Rn×τp, g′i,j ∼iid N (0, 1) independent of Xo, let YO = [Y ,G′] denote the input with outlier G′, then:
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YO‖44 − E ‖W ∗YO‖44∣∣∣ ≥ δ ) < 1
p , (20)
when p = Ω(τ2n2 lnn/δ2).
Proof See Appendix A.4.
In the above results, Proposition 3.1 and 3.3 reveal that both normalized 1npE ‖W ∗YN‖44,
1 npE ‖W ∗YO‖44 are only determined by ‖W ∗Do‖ 4 4. Moreover, as shown in Theorem 3.2 and 3.4, when p is large enough (p = Ω((1 + η2)4n2 lnn),Ω(τ2n2 lnn) respectively), both 1np ‖W ∗YN‖44,
1 np ‖W ∗YO‖44 concentrate around their expectation with high probability. Therefore, the `4-norm maximization formulation ‖W ∗Y ‖44 for dictionary learning is insensitive to dense Gaussian noise and robust to Gaussian outliers.9
9In the outlier case, Proposition 3.3 and Theorem 3.4 can be generalized to any rotation invariant distributions, e.g., uniform distribution on the sphere. More details regarding rotation invariant distributions can be found in Chapter 4 of Bryc (2012).
Proposition 3.5 (Expectation of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Bernoulli matrix B ∈ Rn×p, bi,j ∼iid Ber(β) independent of Xo, let YC = Y + σB ◦ S denote the data with sparse corruptions, and S ∈ Rn×p is defined in equation 15. Then the expectation of 1np ‖W
∗YC‖44 satisfies: 1
np EXo,B,S ‖W ∗YC‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + σ4β(1− 3β) ‖W ‖44 n + Cθ,σ,β , (21)
where Cθ,σ,β is a constant depending on θ, σ and β.
Proof See Appendix A.5.
Theorem 3.6 (Concentration of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Bernoulli matrix B ∈ Rn×p, bi,j ∼iid Ber(β) independent of Xo, let YC = Y + σB ◦ S denote the input with sparse corruptions, and S ∈ Rn×p is defined in equation 15, then:
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YC‖44 − E ‖W ∗YC‖44∣∣∣ ≥ δ ) < 1
p , (22)
when p = Ω ( σ8βn2 lnn/δ2 ) .
Proof See Appendix A.6.
Unlike the cases with noise and outliers, Proposition 3.5 and Theorem 3.6 indicate that 1npE ‖W
∗YC‖44 depends on both ‖W ∗Do‖44 and ‖W ‖ 4 4; when p = Ω(σ
8βn2 lnn), the objective 1 np ‖W
∗YC‖44 concentrates around this expectation with high probability. Nevertheless, when the magnitude of σ4β(1− 3β) is significantly smaller than 3θ(1 − θ), the landscape of the objective ‖W ∗YC‖44 would largely be determined by ‖W ∗Do‖ 4 4 only. As shown in Figure 1, this is indeed the case whenever: (a) the sparsity level θ of ground truth signal Xo, is “reasonably” small (neither diminishing to 0 nor larger than 0.5);(b) β, the sparsity level of the corruption, is small (smaller than 0.5); (c) σ, the magnitude of the sparse errors, is not significantly larger than the in-
trinsic variance of the sparse signal (the intrinsic variance of the sparse signal Bernoulli-Gaussian model is 1).
Therefore, besides the insensitivity to small noise and robustness to outliers, the `4-maximization based dictionary learning also shows resilience to sparse corruptions under reasonable conditions.

heading: 4 SIMULATIONS AND EXPERIMENTS

heading: 4.1 QUANTITATIVE EVALUATION: SIMULATIONS ON SYNTHETIC DATA
text: Single Trial of MSP. In this simulation, we run the MSP algorithm from equation 3, using the imperfect measurements Y of different models (YN ,YO,YC). As shown in Figure 2, the normalized value of ‖W ∗Do‖44 /n reaches global maximum with all types of inputs when varying the level of noise, outliers, and sparse corruptions. Moreover, as the scale of imperfect measurements increase, Figure 2 shows that (a) the iterations for convergence increases and (b) the final objective value ‖W ∗Do‖44 decreases almost negligibly. This numerical experiment suggests that the MSP algorithm is able to identify the ground truth orthogonal transformation Do despite different types of imperfect measurements.
Phase Transition. Next, we conduct extensive simulations to study the relation between recovery accuracy and sample size p. We run the experiments by increasing the sample size p w.r.t. the scale of imperfect measurements η, τ, β, respectively. As shown in Figure 3, the MSP algorithm (3)
demonstrates a clear phase transition behavior w.r.t. noise, outliers, and sparse corruptions. Such phenomena suggest that the algorithm is inherently stable and robust to certain amounts of noise, outliers, and sparse corruptions. The results corroborate our concentration results in Theorem 3.2, 3.4 and 3.6—a larger sample size p increases the accuracy and robustness of the MSP algorithm (3) for all types of nuisances.
Comparison with Prior Arts. We also compare the MSP algorithm (Zhai et al., 2019b) with previous complete dictionary learning algorithms: Subgradient (SG) (Bai et al., 2018) and Riemannian Trust Region (RTR) methods (Sun et al., 2015). While SG and RTR demonstrate slightly better accuracy in some cases (e.g. small Gaussian noise and Gaussian outliers), both algorithms appear unstable to sparse corruptions. Meanwhile, the MSP algorithm is stable to all imperfect measurement models mentioned in this paper and runs significantly faster than the others.

heading: 4.2 QUALITATIVE EVALUATION: EXPERIMENTS ON REAL IMAGES AND PATCHES
text: Besides simulations, we also conduct extensive experiments to verify the stability and robustness of the MSP algorithm with real imagery data, at both image level and patch level. Throughout these experiments, rather than visualize all bases, we routinely show the top bases learned—heuristically, the top bases are those with the largest coefficients (here, in terms of `1-norm). Image Level. At the image level, we first vectorize all 60,000 images in the MNIST dataset (LeCun et al., 1998) into a single matrix Y ∈ R784×60,000, then create imperfect measurements based on models specified in Section 3: YN (MNIST with noise), YO (MNIST with outliers), YC (MNIST with sparse corruptions). We run the MSP algorithm 3 with Y ,YN ,YC ,YO and compare the learned bases. Figure 4(a), (c), (e), and (g) show examples of Y , YN , YO, and YC , and Figure 4(b), (d), (f), and (h) show top 10 bases learned from Y ,YN ,YC ,YO, respectively. Despite that we use different types of imperfect measurements of MNIST, the top bases learned from MSP algorithm 3 are very much the same.10 This result corroborates our analysis: the `4-maximization and the MSP algorithm is inherently insensitive to noise, robust to outliers, and resilient to sparse corruptions.
Patch Level. A classic application of dictionary learning involves learning sparse representations of image patches (Elad & Aharon, 2006; Mairal et al., 2007). In this section, we extend the experiments of Zhai et al. (2019b) to learn patches from gray scale and color images. First, we construct a data matrix Y by vectorizing each 8× 8 patch from the 512× 512 gray scale image, “Barbara” (see Figure 5). We then run the MSP algorithm with 100 iterations on both Y and a noisy version YN , and the learned top bases are visualized in Figure 5.
Analogously, we apply the same scheme to a 256×256 color image, “Duck” (see Figure 6), converting each 8×8×3 patch into a column vector (in R192) of Y . Notice this forces the algorithm to learn bases involving all three channels simultaneously, rather than one at a time. After running the MSP algorithm for 100 iterations, we visualize the top bases learned from both Y and corresponding YN in Figure 6.
We next consider the problem of learning a “global dictionary” (Mairal et al., 2007) for patches from many different images. To construct our data matrix Y , we randomly sample 100,000 8 × 8 × 3
10Bases with opposite intensity are considered as the same base, since they only differ by a sign.
patches from the CIFAR-10 data-set (Krizhevsky et al., 2009). A noisy data matrix YN , is then generated by adding Gaussian noise. Again, we apply the MSP algorithm with 200 iterations to learn 192 bases and visualize the results in Figure 7. We leave the experiments of CIFAR-10 with outliers and sparse corruptions in the Appendix due to limited space. In each of these experiments, the top bases in the learned dictionary remain relatively unchanged with the addition of noise. To quantify this similarity, we take the top bases from the noisy dictionary and find the closest top clean base for each. If the bases are nearly identical, then the inner product of each of these pairs should be close to 1. Table 3 reports the statistics.

heading: 5 CONCLUSION AND DISCUSSIONS
text: In this paper, we find the `4-norm maximization based dictionary learning and corresponding MSP algorithm introduced by Zhai et al. (2019b) have strong geometric and statistical connections to classic data analysis methods PCA and ICA. These connections seem to be the reason why they all admit similar simple and efficient algorithms. Empirically, we have observed that `4-norm maximization is surprisingly insensitive to noise, robust to outliers, and resilient to sparse corruptions. Moreover, such empirical observations corroborate our concentration analysis—larger data samples (p) improve the stability and robustness of the `4-maximization based dictionary learning.
From experiments on real images, we observed that top bases learned are rather stable but tail bases can be less stable (see Figure 10 in Appendix C). We conjecture this phenomenon happens due to the fact that real images generally do not follow the uniformly sparse Bernoulli Gaussian model (of equation 14). Generalizing dictionary learning to non-uniformly sparsely generated data would be a good topic for future study. Finally, note that Bai et al. (2018) has established subgradient descent can find global optimal solutions for dictionary learning, an interesting theoretical guarantee since dictionary learning problems do not satisfy the standard structural (i.e. generalized convexity) assumptions that guarantee global convergence (Ben-Tal & Nemirovski, 2001; Zhou et al., 2017a; Nesterov, 2013; Zhou et al., 2017b; Bubeck et al., 2015; Zhou et al., 2017c; 2018; Ma et al., 2018; Chi et al., 2019). It would be desirable, although highly challenging, to establish similar global convergence guarantees for our MSP algorithm. We leave that for future work.

heading: A PROOF OF SECTION 3

heading: A.1 PROOF OF PROPOSITION 3.1
text: Claim A.1 (Expectation of Objective with Small Noise) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G ∈ Rn×p, gi,j ∼iid N (0, η2) independent of Xo, let YN = Y +G denote the data with noise. Then the expectation of 1np ‖W ∗YN‖44 satisfies:
1
np EXo,G ‖W ∗YN‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + 3θ2 + 6θη2 + 3η4. (23)
Proof Let W ∗Do = M ∈ O(n;R), notice that the orthogonal transformation (W ∗G) of a Gaussian matrix (G) is still a Gaussian matrix and satisfies {W ∗G}i,j ∼ N (0, 1), and it is independent of Y (Xo). We abuse the notation a bit let G = W ∗G in the following calculation, since W ∗G is also a Gaussian matrix independent of Xo and it will not affect the final result.
EXo,G ‖W ∗YN‖ 4 4 = EXo,G ‖MXo + G‖ 4 4
= p∑ j=1 n∑ i=1 EXo,G ( n∑ k=1 mi,kxk,j + gi,j )4
= p∑ j=1 n∑ i=1 EXo,G ( n∑ k=1 mi,kxk,j )4 + 6EXo,G n2i,j ( n∑ k=1 mi,kxk,j )2+ EGg4i,j 
= p∑ j=1 n∑ i=1 EXo,G ( n∑ k=1 mi,kxk,j )4 + 6η2EXo ( n∑ k=1 mi,kxk,j )2+ 3npη4 =
p∑ j=1 n∑ i=1 EXo,G ( n∑ k=1 mi,kxk,j )4+ 6npθη2+3npη4 =EXo ‖MXo‖ 4 4 + 6npθη 2 + 3npη4 =3pθ(1− θ) ‖M‖44 + 3npθ 2 + 6npθη2 + 3npη4,
(24)
therefore, 1
np EXo,G ‖W ∗YN‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + 3θ2 + 6θη2 + 3η4, (25)
which completes the proof.

heading: A.2 PROOF OF THEOREM 3.2
text: Claim A.2 (Concentration of Objective with Small Noise) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G ∈ Rn×p, gi,j ∼iid N (0, η2) independent of Xo, let YN = Y + G denote the input with noise, then:
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YN‖44 − E ‖W ∗YN‖44∣∣∣ ≥ δ ) < 1
p , (26)
when p = Ω ( (1 + η2)4n2 lnn/δ2 ) .
Proof According to Proposition 3.1, we know that
1
np EXo,G ‖W ∗YN‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + 3θ2 + 6θη2 + 3η4. (27)
Moreover, since YN = [y1 + g1,y2 + g2, . . . ,yp + gp], whose columns are independent. Let
Z = [z1, z2, . . . ,zp] = D ∗ oYN = [x1 + g1,x2 + g2, . . . ,xp + gp] = Xo + G, (28)
note that we abuse the notation a bit by assuming D∗oG = G, since D ∗ oG is also a Gaussian matrix independent of Xo and it will not affect the final result. Define function fz(W ) : O(n;R) 7→ R as
fz(W ) = ‖W ∗z‖44 . (29) Next, we will check assumption 1, 2, and 3 then apply Lemma B.4.
Assumption 1: µ(n, p). Since zi,j = xi,j+gi,j , we know that with probability θ, z ∼ N (0, 1+η2) and with probability 1− θ, z ∼ N (0, 1). Hence,
P (
max i,j |zi,j | > B
) ≤ 2npθ exp ( − B 2
2 + 2η2
) + 2np(1− θ) exp ( −B 2
2
) = µ(n, p). (30)
Specifically, we set B = ln p, which yields
µ(n, p) = 2npθ exp ( − (ln p) 2
2 + 2η2
) + 2np(1− θ) exp ( − (ln p) 2
2
) . (31)
Assumption 2: Lipschitz Constant Lf . By Proposition 3.1, ∀W ∈ O(n;R), we know that
Efz(W ) = 3θ(1− θ) ‖W ∗Do‖44 + 3nθ 2 + 6nθη2 + 3nη4. (32)
Hence, ∀W1,W2 ∈ O(n;R), we have∣∣∣E ‖W ∗1 z‖44 − E ‖W ∗1 z‖44∣∣∣ = 3θ(1− θ) ∣∣∣‖W ∗1 Do‖44 − ‖W ∗2 Do‖44∣∣∣ , (33) moreover,∣∣∣∣∣∣ ∑ i,j [ (W ∗1 Do) 4 i,j − (W ∗2 Do)4i,j
]∣∣∣∣∣∣ =
∣∣∣∣∣∣ ∑ i,j {[ (W ∗1 Do) 2 i,j − (W ∗2 Do)2i,j ] [ (W ∗1 Do) 2 i,j + (W ∗ 2 Do) 2 i,j ]}∣∣∣∣∣∣ ≤
∑ i,j [ (W ∗1 Do) 2 i,j − (W ∗2 Do)2i,j ]2 1/2∑ i,j [ (W ∗1 Do) 2 i,j + (W ∗ 2 Do) 2 i,j ]2 1/2
= ∥∥(W ∗1 Do)◦2 − (W ∗1 Do)◦2∥∥F ∥∥(W ∗1 Do)◦2 + (W ∗1 Do)◦2∥∥F ≤‖W ∗1 Do −W ∗2 Do‖F ‖W ∗ 1 Do + W ∗ 2 Do‖F ( ‖W ∗1 Do‖ 2 F + ‖W ∗ 2 Do‖ 2 F
) ≤4n2 ‖W ∗1 Do −W ∗2 Do‖2 = 4n 2 ‖W1 −W2‖2 .
(34)
Hence, we know that ∣∣∣E ‖W ∗1 z‖44 − E ‖W ∗1 z‖44∣∣∣ ≤ 12n2θ(1− θ), (35) which yields
Lf = 12n 2θ(1− θ). (36)
Assumption 2: Lipschitz Constant L̄f . Since fz̄(W ) = ‖W ∗z̄‖44, we have∣∣∣‖W ∗1 z̄‖44 − ‖W ∗2 z̄‖44∣∣∣ = ∣∣∣∣∣∑ i {[ (W ∗1 z̄) 2 i − (W ∗2 z̄)2i ] [ (W ∗1 z̄) 2 i + (W ∗ 2 z̄) 2 i ]}∣∣∣∣∣ ≤
{∑ i [ (W ∗1 z̄) 2 i − (W ∗2 z̄)2i ]2}1/2{∑ i [ (W ∗1 z̄) 2 i + (W ∗ 2 z̄) 2 i ]2}1/2 = ∥∥(W ∗1 z̄)◦2 − (W ∗2 z̄)◦2∥∥2︸ ︷︷ ︸
Γ1
∥∥(W ∗1 z̄)◦2 + (W ∗2 z̄)◦2∥∥2︸ ︷︷ ︸ Γ2 .
(37)
For Γ1, we have Γ1 = ∥∥(W ∗1 z̄)◦2 − (W ∗2 z̄)◦2∥∥2 = ‖(W ∗1 z̄)− (W ∗2 z̄) ◦ (W ∗1 z̄) + (W ∗2 z̄)‖2 ≤‖W1 −W2‖2 ‖W1 + W2‖2︸ ︷︷ ︸
≤2
‖z̄‖22︸︷︷︸ ≤nB2 ≤ 2nB2 ‖W1 −W2‖2 . (38)
For Γ2, we have Γ2 = ∥∥(W ∗1 z̄)◦2 + (W ∗2 z̄)◦2∥∥2 ≤ ∥∥(W ∗1 z̄)◦2∥∥2 + ∥∥(W ∗2 z̄)◦2∥∥2 ≤‖W ∗1 z̄‖ 2 2 + ‖W ∗ 2 z̄‖ 2 2 ≤ 2nB 2. (39) Hence, we know that ∣∣∣‖W ∗1 z̄‖44 − ‖W ∗2 z̄‖44∣∣∣ ≤ 4n2B4 ‖W1 −W2‖2 , (40) which yields
L̄f = 4n 2B4. (41)
Assumption 3: Upper Bound R1. Notice that ∀W ∈ O(n;R),
fz̄(W ) = ‖W ∗z̄‖44 = ‖W ∗z̄‖42 ∥∥∥∥ W ∗z̄‖W ∗z̄‖2 ∥∥∥∥4 4 ≤ ‖W ∗z̄‖42 = ‖z̄‖ 4 2 ≤ n 2B4 = R1. (42)
Assumption 3: Upper Bound R2. Assume the support for x is S , so we know that know that ∀i ∈ [n],
zi =
{√ 1 + η2vi, if i ∈ S,
vi, otherwise, (43)
where v ∼ N (0, I) is a Gaussian vector. Let PS : Rn 7→ Rn be a projection that satisfies ∀q ∈ Rn( PSq ) i = {√ 1 + η2qi, if i ∈ S,
qi, otherwise. (44)
Moreover, ∀W ∈ O(n;R), let wi denote the ith column vector of W , hence
E[f2z(W )] = E [( ‖W ∗z‖44 )2] = E  n∑ i=1 n∑ j=1 〈wi, z〉4 〈wj , z〉4 
≤ES n∑ i=1 n∑ j=1 E [ 〈PSwi,v〉4 〈PSwj ,v〉4 ] ≤ ES n∑ i=1 n∑ j=1 [ E 〈PSwi,v〉8 E 〈PSwj ,v〉8 ] 1 2 . (45)
Notice that v ∼ N (0, I), therefore ∀i ∈ [n], we have
E 〈PSwi,v〉8 = 105 ‖PSwi‖82 , (46)
hence
ES n∑ i=1 n∑ j=1 [ E 〈PSwi,v〉8 E 〈PSwj ,v〉8 ] 1 2 = 105ES n∑ i=1 n∑ j=1 ‖PSwi‖42 ‖PSwi‖ 4 2
=105 n∑ i=1 n∑ j=1 ∑ k1,k2.k3,k4 ES [ 4∏ l=1 (1 + η21kl∈S)w 2 i,k1w 2 i,k2w 2 j,k3w 2 j,k4 ]
≤105(1 + η2)4 n∑ i=1 n∑ j=1 ∑ k1,k2.k3,k4 [w2i,k1w 2 i,k2w 2 j,k3w 2 j,k4 ] = 105(1 + η 2)4n2.
(47)
Therefore, we can conclude that
E[f2z(W )] ≤ 105(1 + η2)4 = R2. (48)
Applying Lemma B.4 for Concentration. Now we apply Lemma B.4 with
1.
B = ln p, µ(n, p) = 2npθ exp ( − (ln p) 2
4
) + 2np(1− θ) exp ( − (ln p) 2
2
) , (49)
2. Lf = 12n 2θ(1− θ), L̄f = 4n2(ln p)4, (50)
3. R1 = n 2(ln p)4, R2 = 105(1 + η 2)4n2, (51)
we have
P  sup W∈O(n;R) 1 np ∣∣∣∣∣∣ p∑ j=1 [ fzj (W )− Efzj (W ) ]∣∣∣∣∣∣ ≥ δ 
=P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YN‖44 − E ‖W ∗YN‖44∣∣∣ ≥ δ )
< exp [ − pn 2δ2
32R2 + 8R1nδ/3 + n2 ln
( 12(Lf + L̄f )
nδ
) + ln 2 ] + µ(n, p)
< exp [ − 3pδ 2
C(1 + η2)4 + 8n(ln p)4δ + n2 ln
( 60n(ln p)4
δ
) + ln 2 ] + 2npθ exp ( − (ln p) 2
4
) + 2np(1− θ) exp ( − (ln p) 2
2
) < 1
p ,
(52)
for a constant C > 104, when p = Ω ( (1 + η2)4n2 lnn/δ2 ) .

heading: A.3 PROOF OF PROPOSITION 3.3
text: Claim A.3 (Expectation of Objective with Outliers) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G′ ∈ Rn×τp, g′i,j ∼iid N (0, 1) independent of Xo, let YO = [Y ,G′] denote the data with outliers G′. Then the expectation of 1 np ‖W ∗YO‖44 satisfies:
1
np EXo,G′ ‖W ∗YO‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + 3θ2 + 3τ. (53)
Proof Notice that
EXo,G′ ‖W ∗YO‖ 4 4 = EXo ‖W ∗Y ‖44 + EG′ ‖W ∗G′‖44 , (54)
and EXo ‖W ∗Y ‖ 4 4 = EXo ‖W ∗DoXo‖44 = 3pθ(1− θ) ‖W ∗Do‖44 + 3npθ 2. (55) Moreover, the orthogonal rotation (W ∗G′) of a standard Gaussian matrix G′ is also a standard Gaussian matrix, therefore,
EG′ ‖W ∗G′‖ 4 4 = 3τnp. (56)
Hence, 1
np EXo,G′ ‖W ∗YO‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + 3θ2 + 3τ, (57)
which completes the proof.

heading: A.4 PROOF OF THEOREM 3.4
text: Claim A.4 (Concentration of Objective with Outliers) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Gaussian matrix G′ ∈ Rn×τp, g′i,j ∼iid N (0, 1) independent of Xo, let YO = [Y ,G′] denote the input with outlier G′, then:
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YO‖44 − E ‖W ∗YO‖44∣∣∣ ≥ δ ) < 1
p , (58)
when p = Ω(τ2n2 lnn/δ2).
Proof Notice that
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YO‖44 − E ‖W ∗YO‖44∣∣∣ ≥ δ )
≤P ( sup
W∈O(n;R)
1
np ∣∣∣‖W ∗Y ‖44 − E ‖W ∗Y ‖44∣∣∣ ≥ δ2 )
︸ ︷︷ ︸ Γ1
+ P ( sup
W∈O(n;R)
1
np ∣∣∣‖W ∗G′‖44 − E ‖W ∗G′‖44∣∣∣ ≥ δ2 )
︸ ︷︷ ︸ Γ2
.
(59)
Apply Lemma B.5 (substitute δ with δ2 ) to Γ1, we know that
Γ1 =P ( sup
W∈O(n;R)
1
np ∣∣∣‖W ∗Y ‖44 − E ‖W ∗Y ‖44∣∣∣ ≥ δ2 )
< exp ( − 3pδ 2
c1θ + 16n(ln p)4δ + n2 ln (120np(ln p)4 δ
))
+ exp ( − pδ 2
c2θ + n2 ln (120np(ln p)4 δ )) + 2npθ exp ( − (ln p) 2 2 ) < 1 2p ,
(60)
for two constants c1 ≥ 4 × 104, c2 > 13, 440, and the last inequality holds when p = Ω(θn2 lnn/δ2). Moreover, apply Lemma B.5 again (let θ = 1,Y = G′ ∈ Rn×τp and substitute δ with δ2τ ), we have
Γ2 =P ( sup
W∈O(n;R)
1
np ∣∣∣‖W ∗G′‖44 − E ‖W ∗G′‖44∣∣∣ ≥ δ2 )
=P ( sup
W∈O(n;R)
1
τnp ∣∣∣‖W ∗G′‖44 − E ‖W ∗G′‖44∣∣∣ ≥ δ2τ )
< exp ( − 3p(δ/τ) 2
c1 + 16n(ln p)4δ/τ + n2 ln (120np(ln p)4 δ/τ
))
+ exp ( − p(δ/τ) 2
c2 + n2 ln (120np(ln p)4 δ/τ )) + 2np exp ( − (ln p) 2 2 ) < 1 2p ,
(61)
for two constants c1 ≥ 4 × 104, c2 > 13, 440, and the last inequality holds when p = Ω(τ2n2 lnn/δ2). Combine θ ≤ 1 and τ can be larger than 1, we conclude that
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YO‖44 − E ‖W ∗YO‖44∣∣∣ ≥ δ ) < Γ1 + Γ2 < 1
p (62)
when p = Ω ( τ2n2 lnn/δ2 ) .

heading: A.5 PROOF OF PROPOSITION 3.5
text: Claim A.5 (Expectation of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Bernoulli matrix B ∈ Rn×p, bi,j ∼iid Ber(β) independent of Xo, let YC = Y + σB ◦S denote the data with sparse corruptions, and S ∈ Rn×p is defined in equation 15. Then the expectation of 1np ‖W ∗YC‖44 satisfies:
1
np EXo,B,S ‖W ∗YC‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + σ4β(1− 3β) ‖W ‖44 n + 3θ2 + 6σ2θβ + 3σ4β2.
(63)
Proof Let W ∗Do = M ∈ O(n;R), notice that
‖W ∗YC‖44 = ‖MXo + σW ∗(B ◦ S)‖44 , (64)
hence
EXo,B,S ‖W ∗YC‖ 4 4 = EXo,B,S p∑ j=1 n∑ i=1 ( n∑ k=1 mi,kxk,j + σ n∑ k=1 wk,ibk,jsk,j )4
= p∑ j=1 n∑ i=1 EXo,B,S [( n∑ k=1 mi,kxk,j )4 ︸ ︷︷ ︸
Γ1
+6 ( n∑ k=1 mi,kxk,j )2( σ n∑ k=1 wk,ibk,jsk,j )2 ︸ ︷︷ ︸
Γ2
+ ( σ
n∑ k=1 wk,ibk,jsk,j )4 ︸ ︷︷ ︸
Γ3
+4 ( n∑ k=1 mi,kxk,j )3( σ n∑ k=1 wk,ibk,jsk,j ) ︸ ︷︷ ︸
Γ4
+ 4 ( n∑ k=1 mi,kxk,j )( σ n∑ k=1 wk,ibk,jsk,j )3 ︸ ︷︷ ︸
Γ5
] .
(65)
Moreover,
• Γ1 :
EXo,B,SΓ1 = EXoΓ1 = 3θ n∑ k=1 m4i,k + 6θ 2 ∑ 1≤k1<k2≤n m2i,k1m 2 i,k2
= 3θ(1− θ) n∑ k=1 m4i,k + 3θ 2,
(66)
• Γ2 :
EXo,B,SΓ2 = [ EXo ( n∑ k=1 mi,kxk,j )2 ][ EB,S ( σ n∑ k=1 wk,ibk,jsk,j )2 ]
= θ ( n∑ k=1 m2i,k ) σ2β ( n∑ k=1 w2k,i ) = σ2θβ,
(67)
• Γ3 :
EXo,B,SΓ3 =EB,SΓ3 = σ4β n∑ k=1 w4k,i + 6σ 4β2 ∑ 1≤k1<k2≤n w2k1,iw 2 k2,i
=σ4β(1− 3β) n∑ k=1 w4k,i + 3σ 4β2,
(68)
• Γ4,Γ5 : EXo,B,SΓ4 = 0, EXo,B,SΓ5 = 0. (69)
Substitute EΓ1,EΓ2,EΓ3,EΓ4,EΓ5 back to equation 65, yields
1
np EXo,B,S ‖WYC‖ 4 4 = 3θ(1−θ) ‖W ∗Do‖44 n +σ4β(1−3β) ‖W ‖44 n +3θ2+6σ2θβ+3σ4β2. (70)

heading: A.6 PROOF OF THEOREM 3.6
text: Claim A.6 (Concentration of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let Xo ∈ Rn×p, xi,j ∼iid BG(θ). Let Do ∈ O(n;R) be an orthogonal matrix and assume Y = DoXo. For any orthogonal matrix W ∈ O(n;R) and any random Bernoulli matrix B ∈ Rn×p, bi,j ∼iid Ber(β) independent of Xo, let YC = Y +σB ◦S denote the input with sparse corruptions, and S ∈ Rn×p is defined in equation 15, then:
P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YC‖44 − E ‖W ∗YC‖44∣∣∣ ≥ δ ) < 1
p , (71)
when p = Ω ( σ8βn2 lnn/δ2 ) .
Proof According to Proposition 3.5, we know that
1
np EXo,B,S ‖W ∗YC‖ 4 4 = 3θ(1− θ) ‖W ∗Do‖44 n + σ4β(1− 3β) ‖W ‖44 n + 3θ2 + 6σ2θβ + 3σ4β2. (72) Moreover, since YC = [y1 + σb1 ◦ s1,y2 + σb2 ◦ s2, . . . ,yp + σbp ◦ sp], whose columns are independent. Let
Z = [z1, z2, . . . ,zp] = D ∗ oYC = Xo + σC, (73)
where C = D∗o(B ◦ S). Note that the columns of Z are independent, we can define function fz(W ) : O(n;R) 7→ R as fz(W ) = ‖W ∗z‖44 . (74) Next, we will check assumption 1, 2, and 3 then apply Lemma B.4.
Assumption 1: µ(n, p). Since Z = Xo + σC, we have
{max i,j |zi,j | > B} = {max i,j |xi,j + σci,j | > B} ⊆ {max i,j |xi,j | > B − σmax i,j |ci,j |}. (75)
Moreover, since max i,j |ci,j | ≤ max i ‖di‖1 ≤ √ n, (76)
we have {max i,j |zi,j | > B} ⊆ {max i,j |xi,j | > B − σ √ n}. (77)
Thus we know that P (
max i,j |zi,j | > B
) ≤ P ( max i,j |xi,j | > B − σ √ n ) ≤ 2npθ exp ( − (B − σ √ n) 2 2 ) . (78)
Specifically, we set B = p 1 4 , which yields
µ(n, p) = 2npθ exp − ( p 1 4 − σ √ n )2
2
 . (79)
Assumption 2: Lipschitz Constant Lf . By Proposition 3.5, ∀W ∈ O(n;R), we know that
Efz(W ) = 3θ(1− θ) ‖W ∗Do‖44 + σ 4β(1− 3β) ‖W ‖44 + 3nθ 2 + 6nσ2θβ + 3nσ4β2. (80)
Hence ∀W1,W2 ∈ O(n;R), we have∣∣∣E ‖W ∗1 z‖44 − E ‖W ∗wz‖44∣∣∣ ≤3θ(1− θ)
∣∣∣‖W ∗1 Do‖44 − ‖W ∗1 Do‖44∣∣∣+ σ4β(1− 3β) ∣∣∣‖W1‖44 − ‖W2‖44∣∣∣ , (81) from the derivation of equation 34, we have∣∣∣‖W ∗1 Do‖44 − ‖W ∗2 Do‖44∣∣∣ ≤ 4n2 ‖W1 −W2‖2 , (82) and therefore, we have∣∣∣E ‖W ∗1 z‖44 − E ‖W ∗wz‖44∣∣∣ ≤ [12n2θ(1− θ) + 4n2σ4β(1− 3β)] ‖W1 −W2‖2 , (83) which yields
Lf = 12n 2θ(1− θ) + 4n2σ4β(1− 3β). (84)
Assumption 2: Lipschitz Constant L̄f . Since fz̄(W ) = ‖W ∗z̄‖44, apply the same derivation as equation 37, we have ∣∣∣‖W ∗1 z̄‖44 − ‖W ∗2 z̄‖44∣∣∣ ≤ 4n2B4 ‖W1 −W2‖2 , (85) which yields
L̄f = 4n 2B4. (86)
Assumption 3: Upper Bound R1. Notice that ∀W ∈ O(n;R),
fz̄(W ) = ‖W ∗z̄‖44 = ‖W ∗z̄‖44 = ‖W ∗z̄‖42 ∥∥∥∥ W ∗z̄‖W ∗z̄‖2 ∥∥∥∥4 4 ≤ ‖W ∗z̄‖42 = ‖z̄‖ 4 2 ≤ n 2B4 = R1.
(87)
Assumption 3: Upper Bound R2. Assume the support for x is S, so we know that ∀i ∈ [n],
zi = { vi + ci, if i ∈ S, ci, otherwise,
(88)
where v ∼ N (0, I) is a Gaussian vector and c = D∗o(b ◦ s) is the sparse corruption vector after orthogonal rotation D∗o . Let PS : Rn 7→ Rn be a projection onto the support S, that is, ∀q ∈ Rn:(
PSq ) i = { qi, if i ∈ S, 0, otherwise.
(89)
Let M = DoW , and h = b ◦ s, we have Ef2z(W ) = E ( ‖W ∗z‖84 ) =E
n∑ i=1 n∑ j=1 [ (〈wi,x〉+ σ 〈mi,h〉)4 (〈wj ,x〉+ σ 〈mj ,h〉)4 ] ≤
n∑ i=1 n∑ j=1 [ E (〈wi,x〉+ σ 〈mi,h〉)8 E (〈wj ,x〉+ σ 〈mj ,h〉)8 ] 1 2
= n∑ i=1 n∑ j=1 [ E (〈PSwi,v〉+ σ 〈mi,h〉)8 E (〈PSwj ,v〉+ σ 〈mj ,h〉)8 ] 1 2 .
(90)
Moreover, we have:
E (〈PSwi,v〉+ σ 〈mi,h〉)8 =E ( 〈PSwi,v〉8 ) + 28σ2E ( 〈PSwi,v〉6 〈mi,h〉2 ) + 70σ4E ( 〈PSwi,v〉4 〈mi,h〉4 ) + 28σ6E ( 〈PSwi,v〉6 〈mi,h〉2 ) + σ8E ( 〈mi,h〉8 ) . (91)
We now provide upper bound for E (〈PSwi,v〉r) and E (〈mi,h〉r) for r = 2, 4, 6, 8 respectively. Since Eh = 0 and h = b ◦ s. Let S ′ denote the support of h, we have
E ( 〈mi,h〉8 ) = E 8∏ l=1
( n∑
k1=1
mkl,ihi ) =
∑ k1,k2,k3,k4 ES′ ( m2k1,i1k1∈S′m 2 k2,i1k2∈S′m 2 k3,i1k3∈S′m 2 k4,i1k4∈S′ ) .
(92)
Now we discuss these cases separately:
• With probability c1β4(c1 ≤ 1), all k1, k2, k3, k4 are in S, in this case, we have∑ k1,k2,k3,k4 ES′ ( m2k1,i1k1∈S′m 2 k2,i1k2∈S′m 2 k3,i1k3∈S′m 2 k4,i1k4∈S′ ) =
∑ k1,k2,k3,k4 ( m2k1,im 2 k2,im 2 k3,im 2 k4,i ) = 1.
(93)
• With probability c2β3(c2 ≤ 1), only three among k1, k2, k3, k4 are in S ′, in this case, we have ∑
k1,k2,k3,k4
ES′ ( m2k1,i1k1∈S′m 2 k2,i1k2∈S′m 2 k3,i1k3∈S′m 2 k4,i1k4∈S′ ) =
∑ k1,k2,k3 ( m2k1,im 2 k2,im 4 k3,i ) = ‖mi‖44 ≤ 1.
(94)
• With probability c3β2(c3 ≤ 1), only two among k1, k2, k3, k4 are in S ′, in this case, we have ∑
k1,k2,k3,k4
ES′ ( m2k1,i1k1∈S′m 2 k2,i1k2∈S′m 2 k3,i1k3∈S′m 2 k4,i1k4∈S′ ) = ∑ k1,k2 ( m4k1,im 4 k2,i ) + ∑ k1,k2 ( m2k1,im 6 k2,i ) ≤ ‖mi‖44 + ‖mi‖ 6 6 ≤ 2.
(95)
• With probability c4β(c4 ≤ 1), only one among k1, k2, k3, k4 is in S ′, in this case, we have∑ k1,k2,k3,k4 ES′ ( m2k1,i1k1∈S′m 2 k2,i1k2∈S′m 2 k3,i1k3∈S′m 2 k4,i1k4∈S′ ) = ∑ k1 m8k1,i = ‖mi‖ 8 8 ≤ 1. (96)
Hence, combine equation 93, equation 94, equation 95, and equation 96, we know that E ( 〈mi,h〉8 ) ≤ cβ, (97)
for a constant c > 1. Similarly, one can conclude that
E (〈mi,h〉r) ≤ cβ, ∀r = 2, 4, 6. (98)
Moreover, since v ∼ N (0, I), we know that
E (〈PSwi,v〉r) = ‖PSwi‖r2 , ∀r = 2, 4, 6, 8. (99)
Hence, with the similar technique applied above, we have E ( 〈PSwi,v〉8 ) = 105E ‖PSwi‖82
=105 ∑
k1,k2,k3,k4
ES′ ( w2k1,i1k1∈S′w 2 k2,i1k2∈S′w 2 k3,i1k3∈S′w 2 k4,i1k4∈S′ ) ≤ 105cθ, (100)
for a constant c > 1. Similarly, we have E (〈PSwi,v〉r) = ‖PSwi‖r2 ≤ c(r − 1)!!θ, ∀r = 2, 4, 6, (101)
for a constant c > 1. Therefore, combine equation 91, we have
E (〈PSwi,v〉+ σ 〈mi,h〉)8 =E ( 〈PSwi,v〉8 ) + 28σ2E ( 〈PSwi,v〉6 〈mi,h〉2 ) + 70σ4E ( 〈PSwi,v〉4 〈mi,h〉4 ) + 28σ6E ( 〈PSwi,v〉6 〈mi,h〉2 ) + σ8E ( 〈mi,h〉8
) ≤cθ + 28cσ2θβ + 70cσ4θβ + 28cσ6θβ + cσ8β ≤ cσ8β,
(102)
for a constant c > 1. Hence, combine equation 90, we have Ef2z(W ) = E ( ‖W ∗z‖84 ) =
n∑ i=1 n∑ j=1 [ E (〈PSwi,v〉+ σ 〈mi,h〉)8 E (〈PSwj ,v〉+ σ 〈mj ,h〉)8 ] 1 2 ≤ cσ8βn2,
(103)
for a constant c > 1. Therefore, we can conclude that R2 = cσ
8βn2, (104) for a constant c > 1.
Applying Lemma B.4 for Concentration. Now we apply Lemma B.4 with
1.
B = p 1 4 , µ(n, p) = 2npθ exp − ( p 1 4 − σ √ n )2
2  , (105) 2.
Lf = 12n 2θ(1− θ) + 4n2σ4β(1− 3β), L̄f = 4n2p, (106)
3. R1 = n 2p, R2 = cσ 8βn2, (107)
for a constant c > 1, we have
P  sup W∈O(n;R) 1 np ∣∣∣∣∣∣ p∑ j=1 [ fzj (W )− Efzj (W ) ]∣∣∣∣∣∣ ≥ δ 
=P ( sup
W∈O(n;R)
1
np
∣∣∣‖W ∗YC‖44 − E ‖W ∗YC‖44∣∣∣ ≥ δ )
< exp [ − pn 2δ2
32R2 + 8R1nδ/3 + n2 ln
( 12(Lf + L̄f )
nδ
) + ln 2 ] + µ(n, p)
< exp [ − 3pδ 2
Cσ8β + 8npδ + n2 ln
( 60np
δ
) + ln 2 ]
+ 2npθ exp − ( p 1 4 − σ √ n )2
2
 < 1 p ,
(108)
for a constant C > 96, when p = Ω ( σ8βn2 lnn/δ2 ) .

heading: B RELATED LEMMAS AND INEQUALITIES
text: Lemma B.1 (Two-sided Bernstein’s Inequality) Given p random variables x1, x2, . . . xp, if ∀i ∈ [p], |xi| ≤ b almost surely, then
P
( 1
p ∣∣∣∣∣ p∑ i=1 [xi − E[xi]] ∣∣∣∣∣ ≥ t ) ≤ 2 exp ( − pt 2 2 p ∑p i=1 E[x2i ] + 2bt/3 ) . (109)
Proof See Proposition 2.14 in Wainwright (2019), one can easily generalize it to two-sided case.
Lemma B.2 (Entry-wise Truncation of a Bernoulli Gaussian Matrix) Let X ∈ Rn×p, where xi,j ∼iid BG(θ) and let ‖·‖∞ denote the maximum element (in absolute value) of a matrix, then
P (
max i,j |xi,j | ≥ t
) ≤ 2npθ exp ( − t 2
2
) . (110)
Proof A Bernoulli Gaussian variable xi,j ,∀i ∈ [n], j ∈ [p] satisfies xi,j = bi,j · gi,j , where bi,j ∼iid Ber(θ), gi,j ∼iid N (0, 1) and therefore
P ( |xi,j | ≥ t ) = θ · P ( |gi,j | ≥ t ) ≤ 2θ exp ( − t 2
2
) . (111)
By union bound, we have:
P (
max i,j |xi,j | ≥ t
) ≤ n∑ i=1 p∑ j=1 P ( |xi,j | ≥ t ) ≤ 2npθ exp
( − t 2
2
) . (112)
Lemma B.3 ( −Net Covering of Stiefel Manifolds) 11 There is a covering −net S for Stiefel manifoldM = {W ∈ Rn×r|W ∗W = I}, (n ≥ r) in operator norm
∀W ∈M, ∃W ′ ∈ S subject to ‖W −W ′‖2 ≤ , (113) of size |S | ≤ ( 6 )nr .
Proof See Lemma D.4 in Zhai et al. (2019b).
Lemma B.4 (Uniform Concentration bound over O(n;R)) Let Z ∈ Rn×p be a random matrix whose columns z1, z2, . . . ,zp are i.i.d. drawn from a distributionP . ∀z ∼ P , let fz(·) : O(n;R) 7→ R denote a function that maps O(n;R) to R. ∀B > 0, let Z̄ denote the truncation of Z:
z̄i,j = { zi,j if |zi,j | ≤ B, 0 otherwise.
(114)

heading: Assume that:
text: 1. P (maxi,j |zi,j | > B) < µ(n, p), where µ(n, p)→ 0 as p increase, B depends on p.
2. Efz(·) is Lf -Lipschitz and fz̄(·) is L̄f -Lipschitz, that is, ∀W1,W2 ∈ O(n;R):
|Efz(W1)− Efz(W2)| ≤ Lf ‖W1 −W2‖2 , |fz̄(W1)− fz̄(W2)| ≤ L̄f ‖W1 −W2‖2 .
(115)
3. ∀W ∈ O(n;R), fz̄(W ) ≤ R1,E[f2z(W )] ≤ R2. 11A similar result can be found in Lemma 4.5 of Recht et al. (2010).

heading: Then:
text: P  sup W∈O(n;R) 1 np ∣∣∣∣∣∣ p∑ j=1 [ fzj (W )− Efzj (W ) ]∣∣∣∣∣∣ ≥ δ 
< exp [ − pn 2δ2
32R2 + 8R1nδ/3 + n2 ln
( 12(Lf + L̄f )
nδ
) + ln 2 ] + µ(n, p),
(116)
when p > ρ, where ρ depends on n.
Proof By assumption 1, we have
P  sup W∈O(n;R) 1 np ∣∣∣∣∣∣ p∑ j=1 [ fzj (W )− Efzj (W ) ]∣∣∣∣∣∣ ≥ δ 
≤P(Z 6= Z̄) + P  sup W∈O(n;R) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≥ nδ 
≤µ(n, p) + P  sup W∈O(n;R) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≥ nδ  .
(117)
Uniform bound over O(n;R). ∀ > 0, lemma B.3 shows there exists an -nets S that covers O(n;R):
S = {W1,W2, . . . ,W|S |}, O(n;R) ⊂ |S |⋃ l=1 B(Wl, ), (118)
and |S | satisfies |S | ≤ (6/ )n 2 . Together with Lipschitz assumption 2, we know that
sup W∈B(Wl, ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≤ sup
W∈B(Wl, ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz(Wl) ∣∣∣∣∣∣+ supW∈B(Wl, ) |Efz(Wl)− Efz(W )| + sup
W∈B(W , ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− 1 p p∑ j=1 fz̄j (W ) ∣∣∣∣∣∣ ≤ sup
W∈B(Wl, ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz(Wl) ∣∣∣∣∣∣+ (Lf + L̄f ) .
(119)
Hence, let
= nδ
2(Lf + L̄f ) , (120)
we have
P  sup W∈O(n;R) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≥ nδ 
≤ |S |∑ l=1 P  sup W∈B(Wl; ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≥ nδ 
<
( 6 )n2 P  sup W∈B(Wl; ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz(Wl) ∣∣∣∣∣∣ ≥ nδ − (Lf + L̄f ) 
= exp [ n2 ln ( 12(Lf + L̄f )
nδ
)] P  sup W∈B(Wl; ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz(Wl) ∣∣∣∣∣∣ ≥ nδ2  .
(121)
Tail bound within each B(Wl, ). Then, ∀l ∈ [S ], we apply point-wise control to a given point Wl ∈ S. Later we will provide a uniform concentration bound over O(n;R). By triangle inequality, we have∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efzj (Wl) ∣∣∣∣∣∣ ≤ ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz̄(Wl) ∣∣∣∣∣∣+ ∣∣Efz̄j (Wl)− Efz(Wl)∣∣ . (122)
∀W ∈ B(Wl, ), by assumption 1 and 3, we have∣∣Efz̄j (W )− Efzj (W )∣∣ = ∣∣E [fz(W ) · 1Z 6=Z̄]∣∣ ≤√E2fz(W )√∣∣E1Z 6=Z̄∣∣ = √ E2fz(W )P(max
i,j |zi,j | > B) ≤
√ R2µ(n, p).
(123)
Let ρ be the lower bound of p, such that ∀p > ρ, we have √ R2µ(n, p) < nδ 4 . Hence, when p > ρ, we have ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efzj (Wl) ∣∣∣∣∣∣ ≤ ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz̄j (Wl)
∣∣∣∣∣∣+ nδ4 , (124) which implies
P ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≥ nδ2  ≤ P ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz̄(W ) ∣∣∣∣∣∣ ≥ nδ4 
≤ 2 exp ( − pn 2δ2
32R2 + 8R1nδ/3
) ,
(125)
where the last inequality is achieved by Bernstein’s inequality (Lemma B.1), along with assumption 3 and Ef2z̄(W ) ≤ Ef2z(W ) ≤ R2. Combine equation 125 and equation 121, we have
P  sup W∈O(n;R) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (W )− Efz(W ) ∣∣∣∣∣∣ ≥ nδ 
< exp [ n2 ln ( 12(Lf + L̄f )
nδ
)] P  sup W∈B(Wl; ) ∣∣∣∣∣∣1p p∑ j=1 fz̄j (Wl)− Efz(Wl) ∣∣∣∣∣∣ ≥ nδ2 
< exp [ − pn 2δ2
32R2 + 8R1nδ/3 + n2 ln
( 12(Lf + L̄f )
nδ
) + ln 2 ] .
(126)
Summary. Therefore, we conclude that
P  sup W∈O(n;R) 1 np ∣∣∣∣∣∣ p∑ j=1 [ fzj (W )− Efzj (W ) ]∣∣∣∣∣∣ ≥ δ 
< exp [ − pn 2δ2
32R2 + 8R1nδ/3 + n2 ln
( 12(Lf + L̄f )
nδ
) + ln 2 ] + µ(n, p),
(127)
when p > ρ.
Lemma B.5 (Concentration Bound of the Clean Objective over O(n;R)) ∀θ ∈ (0, 1], if X ∈ Rn×p, xi,j ∼iid BG(θ), for any δ > 0, the following inequality holds
P ( sup
W∈O(n;R)
1
np
∣∣∣‖WX‖44 − E ‖WX‖44∣∣∣ ≥ δ )
< exp ( − 3pδ 2
c1θ + 8n(ln p)4δ + n2 ln (60np(ln p)4 δ
))
+ exp ( − pδ 2
c2θ + n2 ln (60np(ln p)4 δ )) + 2npθ exp ( − (ln p) 2 2 ) ,
(128)
for some constants c1 > 104, c2 > 3360. Moreover
exp ( − 3pδ 2
c1θ + 8n(ln p)4δ + n2 ln (60np(ln p)4 δ
))
+ exp ( − pδ 2
c2θ + n2 ln (60np(ln p)4 δ )) + 2npθ exp ( − (ln p) 2 2 ) ≤ 1 p ,
(129)
when p = Ω(θn2 lnn/δ2).
Proof See Lemma 2.2 in Zhai et al. (2019b), note that the sparsity condition θ ∈ (0, 1) of the original Lemma in Zhai et al. (2019b) can be easily generalized to θ = 1.

heading: C ADDITIONAL EXPERIMENTAL RESULTS

Review summaries:
This paper presents results on Dictionary Learning through l4 maximization. The authors base this paper heavily off of the formulation and algorithm in Zhai et. al. (2019) "Complete dictionary learning via l4-norm maximization over the orthogonal group". The paper draws connections between complete dictionary learning, PCA, and ICA by pointing out similarities between the objectives functions that are optimized as well as the algorithms used. The paper further presents results on dictionary learning in the presence of different types of noise (AWGN, sparse corruptions, outliers) and show that the l4 objective is robust to different types of noise. Finally the authors apply different types of noise to synthetic and real images and show that the dictionaries that they learn are robust to the noise applied.