Title: UNDERSTANDING 4 -BASED DICTIONARY LEARNING: INTERPRETATION, STABILITY, AND ROBUSTNESS

Abstract: Recently, the 4 -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 4 -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.

1. INTRODUCTION:

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 = [y 1 , y 2 , . . . , y p ] ∈ R n×p that contains p n-dimensional samples, one aims to find a linear transformation (i.e. a dictionary) D ∈ R n×m and an associated maximally sparse representation X = [x 1 , x 2 , . . . , x p ] ∈ R m×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 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 learning 1 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 = D o X o ∈ R n×p satisfies that D o ∈ O(n; R) is orthonormal and X o ∈ R n×p is

Bernoulli-Gaussian sparse, 2 then maximizing the 4 -norm 3 of AY over O(n; R):

EQUATION

is able to find the ground truth dictionary D o 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:

EQUATION

where U t V * t are from the singular value decomposition:

U t Σ t V * t = SVD[(A t Y ) •3 Y * ].

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 1 4 AY 4 4 = ∇ A Λ, AA * -I for some Lagrangian multipliers Λ, implies:

EQUATION

)

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:

EQUATION

Hence the critical point A can be viewed as a "fixed point" of the map:

P O(n;R) ((•)Y ) •3 Y * 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 D o by solving the program:

EQUATION

where St(n, k; R) is the Stiefel manifold. 5 The orthogonal group O(n; R) and the unit sphere S n-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 Journée et al. (2010) in the context of sparse PCA, which has provided convergence guarantees for this class of programs.

1 Complete 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) . 3 We abuse the notation a bit, by denoting • 4 4 as the sum of element-wise 4 th power of all entries of a vector and matrix, that is, ∀a ∈ R n , a 4 4 = n i=1 a 4 i and ∀A ∈ R n×m , A 4 4 = i,j a 4 i,j . 4 For any symmetric matrix S ∈ R n×n and an orthogonal matrix A ∈ O(n; R), the projection of SA onto the orthogonal group is A: P O(n;R) [SA] = A, one may see Absil & Malick (2012) for details. 5 For any

1 ≤ k ≤ n, St(n, k; R) . = {W ∈ R n×k : W * W = I k }.

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 (Hyvä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.

2 SDL VERSUS PCA AND ICA 

EQUATION

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 4 4 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 . = max a∈S n-1 a * Y 4 .

Relation with the Statistical Interpretation of PCA. View each column y j , j ∈ [p] of data matrix Y ∈ R n×p as an n dimensional random vector that is i.i.d. drawn from a distribution of random variable y and let Y c denote the centered Y :

Y c . = Y [I -1 p 11 * ], where 1 ∈ R p is a vector of all 1's. Then, finding the top k principal components of Y c : max W ∈St(n,k;R) 1 2 W * Y c 2 2

is to find k uncorrelated projections of y ∈ R n that has the top k sample variance (i.e. 2 nd order moment) (Jolliffe, 2011; Helwig, 2017) . Similar to PCA, the 4 -norm maximization of centered data matrix Y c : max W ∈St(n,k;R)

1 4 W * Y c 4

4 can be viewed as finding k uncorrelated projections of y that have the top k sample 4 th 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) (Hyvärinen & Oja, 1997; 2000) : In order to identify one component of a given random vector y ∈ R n , ICA aims to find a unit vector (a direction) w ∈ S n-1 that maximizes or minimizes the kurtosis of w * y, defined as:

kurt(w * y) = E (w * y) 4 -3 w 4 2 .

(8) Kurtosis is widely used for evaluating the nonnormality of a random variable; see DeCarlo (1997); Hyvä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 4 th order moment helps understand such statistics of real datasets (Hyvärinen et al., 2009) and even their topology (Carlsson, 2009) . One may also find that the 4maximization based dictionary learning formulation is similar to maximizing kurtosis (equation 8) with spherical constraint w ∈ S n-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.

2.2. FIXED-POINT STYLE ALGORITHMS:

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 Journé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 ∈ R n×p , finding the top singular value of Y is equivalent to solving the following optimization problem:

EQUATION

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 = P S n-1 [∇ w ϕ(w)]. The classic power-iteration method

EQUATION

)

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 S n-1 ,

EQUATION

Hyvärinen & Oja (1997) has proposed the following fixed-point type iteration:

EQUATION

)

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 4 4 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:

EQUATION

where

U t Σ t V * t = SVD[Y (Y * W t ) •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 S n-1 and St(n, k; R), respectively. Table 1 summarizes these striking similarities.

Objectives Constraint Sets Algorithms:

Power Iteration ϕ(w) .

= 1 2 w * Y 2 2 w ∈ S n-1 wt+1 = P S n-1 [∇wϕ(wt)] FastICA ψ(w) . = 1 4 kurt[w * y] w ∈ S n-1 wt+1 = P S n-1 [∇wψ(wt)] MSP φ(W ) . = 1 4 W * Y 4 4 W ∈ St(n, k; R) Wt+1 = P St(n,k;R) [∇ W φ(Wt)]

Table 1 : Similarities among fixed-point algorithms for Power Iteration, FastICA, and MSP.

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 (Candès et al., 2011; Xu et al., 2010; 2012) . In this section, we want to examine the stability and robustness of the 4maximization for different types of imperfect measurement models: small noise, outliers, and sparse corruptions of large magnitude.

3.1. DIFFERENT MODELS FOR IMPERFECT MEASUREMENTS:

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 ∈ R n×p is produced by the product of a ground truth orthogonal dictionary D o and a Bernoulli-Gaussian matrix

EQUATION

Now let us assume we only observe different types of imperfect measurements of Y :

Noisy Measurements: Y N := Y +G, where G ∈ R n×p is matrix that satisfies g i,j ∼ iid N (0, η 2 )

and η > 0 controls the variance of the noise.

Measurements with Outliers:

Y O := [Y , G ],

where Y O contains extra columns (G ∈ R n×τ 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: Y C := Y +σB •S, where σ > 0 controls the scale of corrupting entries,7 B ∈ R n×p is a Bernoulli matrix so b i,j ∼ iid Ber(β) with β ∈ (0, 1) controlling the ratio of the sparse corruptions, and S ∈ R n×p has entries s i,j drawn i.i.d. from a Rademacher distribution:

EQUATION

)

3.2. STATISTICAL ANALYSIS AND JUSTIFICATION:

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 4 4 , subject to W ∈ O(n; R) (16) satisfies W * D o ∈ SP(n).

EQUATION

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 X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Gaussian matrix G ∈ R n×p , g i,j ∼ iid N (0, η 2 ) independent of X o , let Y N = Y + G denote

the input with noise, then:

EQUATION

)

when p = Ω (1 + η 2 ) 4 n 2 ln n/δ 2 . Proof See Appendix A.2. Proposition 3.3 (Expectation of Objective with Outliers) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthog- onal matrix W ∈ O(n; R) and any random Gaussian matrix G ∈ R n×τ p , g i,j ∼ iid N (0, 1) independent of X o , let Y O = [Y , G

] denote the data with outliers G . Then the expectation of

1 np W * Y O 4

4 satisfies:

EQUATION

where C θ,τ is a constant depends on θ, τ .

Proof See Appendix A.3

Theorem 3.4 (Concentration of Objective with Outliers) ∀θ ∈ (0, 1), let

X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Gaussian matrix G ∈ R n×τ p , g i,j ∼ iid N (0, 1) indepen- dent of X o , let Y O = [Y , G

] denote the input with outlier G , then:

EQUATION

)

when p = Ω(τ 2 n 2 ln n/δ 2 ).

Proof See Appendix A.4.

In the above results, Proposition 3.1 and 3.3 reveal that both normalized 

1 np E W * Y N 4 4 , 1 np E W * Y O

EQUATION

4 concentrate around their expectation with high probability. Therefore, the 4 -norm maximization formulation W * Y 4 4 for dictionary learning is insensitive to dense Gaussian noise and robust to Gaussian outliers.9 

EQUATION

where C θ,σ,β is a constant depending on θ, σ and β.

Proof See Appendix A.5. 

EQUATION

)

when p = Ω σ 8 βn 2 ln n/δ 2 .

Proof See Appendix A.6. As shown in Figure 1 , this is indeed the case whenever: (a) the sparsity level θ of ground truth signal X o , 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 intrinsic 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. 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. (Zhai et al., 2019b) with prior complete dictionary learning algorithms: Subgradient method (Bai et al., 2018) and Riemannian Trust Region (Sun et al., 2015) methods in different models with fixed ground truth sparsity θ = 0.3. Note that SG only learns a unit vector each time and does not guarantee orthogonality; we therefore project the dictionary learned to O(n; R) for fair comparison.

4. SIMULATIONS AND EXPERIMENTS:



4.2. QUALITATIVE EVALUATION: EXPERIMENTS ON REAL IMAGES AND PATCHES:

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 (Le-Cun et al., 1998) 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 Y N , and the learned top bases are visualized in Figure 5 . We observed a similar effect when using an 8 × 8 patch size. 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 R 192 ) 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 Y N 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 patches from the CIFAR-10 data-set (Krizhevsky et al., 2009) . A noisy data matrix Y N , 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. 

5. CONCLUSION AND DISCUSSIONS:

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.

A PROOF OF SECTION 3 

EQUATION

) Proof Let W * D o = 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 (X o ). We abuse the notation a bit let G = W * G in the following calculation, since W * G is also a Gaussian matrix independent of X o and it will not affect the final result.

EQUATION

EQUATION

which completes the proof.

A.2 PROOF OF THEOREM 3.2

Claim A.2 (Concentration of Objective with Small Noise) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Gaussian matrix G ∈ R n×p , g i,j ∼ iid N (0, η 2 ) independent of X o , let Y N = Y + G denote

the input with noise, then:

EQUATION

when p = Ω (1 + η 2 ) 4 n 2 ln n/δ 2 .

Proof According to Proposition 3.1, we know that

EQUATION

Moreover, since Y N = [y 1 + g 1 , y 2 + g 2 , . . . , y p + g p ], whose columns are independent. Let

EQUATION

note that we abuse the notation a bit by assuming

D * o G = G, since D * o G

is also a Gaussian matrix independent of X o and it will not affect the final result. Define function

EQUATION

Next, we will check assumption 1, 2, and 3 then apply Lemma B.4.

Assumption 1: µ(n, p). Since z i,j = x i,j +g i,j , we know that with probability θ, z ∼ N (0, 1+η2 ) and with probability 1 -θ, z ∼ N (0, 1). Hence,

EQUATION

Specifically, we set B = ln p, which yields

EQUATION

)

Assumption 2: Lipschitz Constant L f . By Proposition 3.1, ∀W ∈ O(n; R), we know that

EQUATION

EQUATION

Hence, we know that

EQUATION

which yields

EQUATION

Assumption 2: Lipschitz Constant Lf . Since f z (W ) = W * z 4 4 , we have

W * 1 z 4 4 -W * 2 z 4 4 = i (W * 1 z) 2 i -(W * 2 z) 2 i (W * 1 z) 2 i + (W * 2 z) 2 i ≤ i (W * 1 z) 2 i -(W * 2 z) 2 i 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

For Γ 1 , we have

EQUATION

)

For Γ 2 , we have

EQUATION

Hence, we know that

EQUATION

which yields Lf = 4n 2 B 4 . (41)

EQUATION

Assumption 3: Upper Bound R 2 . Assume the support for x is S, so we know that know that ∀i ∈ [n],

EQUATION

where v ∼ N (0, I) is a Gaussian vector. Let P S : R n → R n be a projection that satisfies ∀q ∈ R n

EQUATION

Moreover, ∀W ∈ O(n; R), let w i denote the i th column vector of W , hence

EQUATION

)

Notice that v ∼ N (0, I), therefore ∀i ∈ [n], we have 

EQUATION

E S 4 l=1 (1 + η 2 1 k l ∈S )w 2 i,k1 w 2 i,k2 w 2 j,k3 w 2 j,k4 ≤105(1 + η 2 ) 4 n i=1 n j=1 k1,k2.k3,k4

[w 2 i,k1 w 2 i,k2 w 2 j,k3 w 2 j,k4 ] = 105(1 + η 2 ) 4 n 2 .

(47)

Therefore, we can conclude that 

EQUATION

EQUATION

EQUATION

Moreover,

• Γ 1 :

EQUATION

• Γ 2 :

EQUATION

• Γ 3 : 

EQUATION

EQUATION

when p = Ω σ 8 βn 2 ln n/δ 2 .

Proof According to Proposition 3.5, we know that

EQUATION

where C = D * o (B • S). Note that the columns of Z are independent, we can define function

EQUATION

Next, we will check assumption 1, 2, and 3 then apply Lemma B.4.

Assumption 1: µ(n, p). Since Z = X o + σC, we have

EQUATION

Moreover, since max i,j

EQUATION

we have {max i,j

EQUATION

Thus we know that

EQUATION

Specifically, we set B = p 1 4 , which yields

EQUATION

Moreover, we have:

EQUATION

We now provide upper bound for E ( P S w i , v r ) and E ( m i , 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

EQUATION

Now we discuss these cases separately:

EQUATION

• With probability c 2 β 3 (c 2 ≤ 1), only three among k 1 , k 2 , k 3 , k 4 are in S , in this case, we have k1,k2,k3,k4

EQUATION

• With probability c 3 β 2 (c 3 ≤ 1), only two among k 1 , k 2 , k 3 , k 4 are in S , in this case, we have k1,k2,k3,k4

EQUATION

EQUATION

Hence, combine equation 93, equation 94, equation 95, and equation 96, we know that

EQUATION

for a constant c > 1. Similarly, one can conclude that

EQUATION

Moreover, since v ∼ N (0, I), we know that

EQUATION

Hence, with the similar technique applied above, we have

EQUATION

for a constant c > 1. Similarly, we have

EQUATION

) for a constant c > 1. Therefore, combine equation 91, we have

EQUATION

for a constant c > 1. Hence, combine equation 90, we have

EQUATION

for a constant c > 1. Therefore, we can conclude that R 2 = cσ 8 βn 2 , (104) for a constant c > 1.

Applying Lemma B.4 for Concentration. Now we apply Lemma B.4 with 1.

EQUATION

EQUATION

EQUATION

for a constant c > 1, we have

EQUATION

for a constant C > 96, when p = Ω σ 8 βn 2 ln n/δ 2 .

B RELATED LEMMAS AND INEQUALITIES:

Lemma B.1 (Two-sided Bernstein's Inequality) Given p random variables x 1 , x 2 , . .

EQUATION

)

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 ∈ R n×p , where x i,j ∼ iid BG(θ) and let • ∞ denote the maximum element (in absolute value) of a matrix, then

EQUATION

)

Proof A Bernoulli Gaussian variable x i,j , ∀i ∈ [n], j ∈ [p] satisfies x i,j = b i,j • g i,j , where b i,j ∼ iid Ber(θ), g i,j ∼ iid N (0, 1) and therefore

EQUATION

)

By union bound, we have: Assume that:

P max i,j |x i,j | ≥ t ≤

1. P (max i,j |z i,j | > B) < µ(n, p), where µ(n, p) → 0 as p increase, B depends on p. (121)

2.. Ef:

EQUATION

Tail bound within each B(W l , ). Then, ∀l ∈ [S ], we apply point-wise control to a given point W l ∈ S. Later we will provide a uniform concentration bound over O(n; R). By triangle inequality, we have

EQUATION

∀W ∈ B(W l , ), by assumption 1 and 3, we have (126)

EQUATION

:

In case τ p is not an integer, we round τ p to the closest integer.

In our context, σ = 1 is already corruption of large magnitude, since the variance of the sparse signal is 1.

SP(n) is the signed permutation group, a group of orthogonal matrices that only contain 0, ±1.

In 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 ofBryc (2012).

Bases with opposite intensity are considered as the same base, since they only differ by a sign.

Γ2.(37)

A similar result can be found in Lemma 4.5 ofRecht et al. (2010).

Figures and Tables:

Note here Y denotes different imperfect measurements: noisy (Y N ), with outliers (Y O ), and with sparse corruptions (Y C ). Below we provide the expectation and concentration results for W * Y 4 4 over the data distribution. We show that E W * Y quantity that indicates a "distance" from W * D o to SP(n). As shown in Lemma 2.3 and Lemma 2.4 in Zhai et al. (2019b), the only global maximizers of W * D o 4 4 are signed permutation matrices, and W * D o converges to a signed permutation matrix as W * D o 4 4 reaches its global maximum. Proposition 3.1 (Expectation of Objective with Small Noise) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Gaussian matrix

are only determined by W * D o 4 4 . Moreover, as shown in Theorem 3.2 and 3.4, when p is large enough (

Proposition 3.5 (Expectation of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Bernoulli matrix B ∈ R n×p , b i,j ∼ iid Ber(β) independent of X o , let Y C = Y + σB • S denote the data with sparse corruptions, and S ∈ R n×p is defined in equation 15. Then the expectation of 1 np W * Y C 4 4 satisfies:

Theorem 3.6 (Concentration of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Bernoulli matrix B ∈ R n×p , b i,j ∼ iid Ber(β) independent of X o , let Y C = Y + σB • S denote the input with sparse corruptions, and S ∈ R n×p is defined in equation 15, then:

Figure 1: Comparison between y = 3x(1 -x) and y = σ 4 |x(1 -3x)| when x ∈ [0, 1] with different σ.

Figure 2: Normalized W * D o 4 4 /n of the MSP algorithm for dictionary learning, using imperfect measurements Y N , Y O , Y C , respectively.

into a single matrix Y ∈ R 784×60,000 , then create imperfect measurements based on models specified in Section 3: Y N (MNIST with noise), Y O (MNIST with outliers), Y C (MNIST with sparse corruptions). We run the MSP algorithm 3 with Y , Y N , Y C , Y O and compare the learned bases. Figure 4(a), (c), (e), and (g) show examples of Y , Y N , Y O , and Y C , and Figure 4(b), (d), (f), and (h) show top 10 bases learned from Y , Y N , Y C , Y O , 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.

Figure 4: Left: Examples of MNIST and its different imperfect measurements. Right: Learned bases from MNIST and its different imperfect measurements using the MSP algorithm 3. 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 ofZhai 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 Figure5). We then run the MSP algorithm with 100 iterations on both Y and a noisy version Y N , and the learned top bases are visualized in Figure5.

Figure 5: The top 12 bases learned from all 16 × 16 patches of Barbara, both with (b) and without (a) noise.The noisy image is produced by adding Gaussian noise to the clean image, resulting in a signal-to-noise ratio (SNR) of 5.87. We observed a similar effect when using an 8 × 8 patch size. Analogously, we apply the same scheme to a 256×256 color image, "Duck" (see Figure6), converting each 8×8×3 patch into a column vector (in R 192 ) 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 Y N in Figure6.

Figure 6: The top 12 bases learned from all 8 × 8 × 3 color patches of the clean and noisy image, respectively. Here, the SNR of the noisy image is 6.56.

of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Bernoulli matrix B ∈ R n×p , b i,j ∼ iid Ber(β) independent of X o , let Y C = Y + σB • S denote the data with sparse corruptions, and S ∈ R n×p is defined in equation 15. Then the expectation of 1 np W * Y C 4 4 satisfies:

|x i,j | ≥ t ≤ 2npθ exp --Net Covering of Stiefel Manifolds) 11 There is a covering -net S for Stiefel manifold M = {W ∈ R n×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 ∈ R n×p be a random matrix whose columns z 1 , z 2 , . . . , z p are i.i.d. drawn from a distribution P. ∀z ∼ P, let f z (•) : O(n; R) → R denote a function that maps O(n; R) to R. ∀B > 0, let Z denote the truncation of Z: zi,j = z i,j if |z i,j | ≤ B,

∀W ∈ O(n; R), f z (W ) ≤ R 1 , E[f 2 z (W )] ≤ R 2 . zj (W ) -Ef zj (W ) ≥ δUniform bound over O(n; R). ∀ > 0, lemma B.3 shows there exists an -nets S that covers O(n; R):S = {W 1 , W 2 , . . . , W |S | }, O(n; |S | satisfies |S | ≤ (6/ ) n 2. Together with Lipschitz assumption 2, we know thatsup W ∈B(W l , ) zj (W l ) -Ef z (W l ) + sup W ∈B(W l , ) |Ef z (W l ) -Ef z (zj (W l ) -Ef z (W l ) + (L f + Lf ) .

zj (W l ) -Ef z (W l ) ≥ nδ -(L f + Lf ) zj (W l ) -Ef z (W l )

W )P(max i,j |z i,j | > B) ≤ R 2 µ(n, p).(123)Let ρ be the lower bound of p, such that ∀p > ρ, we have R 2 µ(n, p) < nδ 4 . Hence, when p > ρ, we have1 p p j=1 f zj (W l ) -Ef zj (W l ) ≤ 1 p p j=1 f zj (W l ) -Ef zj (W l ) + nδ4 inequality is achieved by Bernstein's inequality (Lemma B.1), along with assumption 3 and Ef 2 z (W ) ≤ Ef 2 z (W ) ≤ R 2 . Combine equation 125 and equation 121, we have exp n 2 ln 12(L f + Lf ) nδ P   sup W ∈B(W l ; )1 p p j=1 f zj (W l ) -Ef z (W l ) ≥ nδ 2   < exp -pn 2 δ 2 32R 2 + 8R 1 nδ/3+ n 2 ln 12(L f + Lf ) nδ + ln 2 .

Alg. Error Time Error Time Error Time Error Time Error Time Error Time Error Time

Table content: <table><tr><td>Clean</td><td colspan="2">Noise</td><td>Outlier</td><td/><td colspan="2">Corruption</td></tr><tr><td>0.2</td><td/><td>0.4</td><td>0.2</td><td>0.4</td><td>0.2</td><td>0.4</td></tr><tr><td colspan="8">n p 25 10k MSP 0.34% 1.14s 0.45% 1.11s 0.99% 1.12s 1.11% 1.29s 1.82% 1.45s 1.27% 1.05s 2.85% 1.03s</td></tr><tr><td colspan="8">SG 0.00% 8.00m 5.54% 19.3m 25.7% 26.4m 0.00% 9.62m 0.00% 11.1m 87.0% 16.6m 88.4% 34.2m</td></tr><tr><td colspan="8">RTR 0.00% 3.28m 0.03% 17.0m 1.34% 20.2m 0.00% 4.63m 0.00% 5.87m 3.20% 18.0m 33.9% 16.3m</td></tr><tr><td colspan="8">50 20k MSP 0.34% 4.82s 0.47% 4.65s 1.02% 5.44s 1.15% 5.31s 2.01% 6.45s 1.33% 4.80s 3.04% 4.41s</td></tr><tr><td>SG 0.00% 1.34h N/A</td><td>&gt;2h</td><td>N/A</td><td colspan="2">&gt;2h 3.42% 1.40h 0.00% 1.81h N/A</td><td>&gt;2h</td><td>N/A</td><td>&gt;2h</td></tr><tr><td colspan="8">RTR 0.00% 23.0m 0.04% 1.57h 1.65% 1.38h 0.00% 30.7m 0.00% 41.8m 2.17% 1.25h 82.57% 1.29h</td></tr></table>

Comparison of the MSP algorithm

Table content: <table/>

Table 3 reports the statistics.

Table content: <table><tr><td colspan="5">Minimum Lower Quartile Median Upper Quartile Maximum</td></tr><tr><td>Barbara 0.3048</td><td>0.8471</td><td>0.9941</td><td>0.9993</td><td>1.0000</td></tr><tr><td>Duck 0.2510</td><td>0.9782</td><td>0.9891</td><td>0.9971</td><td>1.0000</td></tr><tr><td>CIFAR-10 0.5147</td><td>0.7203</td><td>0.9892</td><td>0.9998</td><td>1.0000</td></tr></table>

Statistics about the inner products between the top 20 noisy bases and their corresponding closest top-20 clean bases.

Table content: <table/>

Claim A.1 (Expectation of Objective with Small Noise) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Gaussian matrix G ∈ R n×p , g i,j ∼ iid N (0, η 2 ) independent of X o , let Y N = Y +G denote the data with noise. Then the expectation of 1 np W * Y N

Table content: <table><tr><td>A.1 PROOF OF PROPOSITION 3.1</td></tr><tr><td>4 4 satisfies:</td></tr></table>

• Γ 4 , Γ 5 : E Xo,B,S Γ 4 = 0, E Xo,B,S Γ 5 = 0. (69) Substitute EΓ 1 , EΓ 2 , EΓ 3 , EΓ 4 , EΓ 5 back to equation 65, yields Claim A.6 (Concentration of Objective with Sparse Corruptions) ∀θ ∈ (0, 1), let X o ∈ R n×p , x i,j ∼ iid BG(θ). Let D o ∈ O(n; R) be an orthogonal matrix and assume Y = D o X o . For any orthogonal matrix W ∈ O(n; R) and any random Bernoulli matrix B ∈ R n×p , b i,j ∼ iid Ber(β) independent of X o , let Y C = Y + σB • S denote the input with sparse corruptions, and S ∈ R n×p is defined in equation 15, then:

Table content: <table><tr><td>1 np</td><td>E Xo,B,S W Y C</td><td>4 4 = 3θ(1-θ)</td><td>W  *  D o n</td><td>4 4</td><td>+σ 4 β(1-3β)</td><td>W n</td><td>4 4</td><td>+3θ 2 +6σ 2 θβ+3σ 4 β 2 . (70)</td></tr><tr><td colspan="3">A.6 PROOF OF THEOREM 3.6</td><td/><td/><td/><td/><td/><td/></tr></table>

Peer review (iclr2020):
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. Overall this paper makes significant contributions by extending the work in the paper referenced above to noisy dictionary learning settings and I would vote to accept based on these results. The connections between Complete Dictionary Learning, PCA and ICA are interesting, but the algorithmic analogies seem superficial in my opinion. There are a lot of algorithms which follow a projected/proximal gradient descent scheme. If there are any deeper connections between the specific algorithms discussed, they should be spelled out more clearly. One point of clarification that I would like to raise is the similarity between the kurtosis and l4 objectives. This paper could be strengthened by delineating the conditions under which one would learn an ICA basis vs a Complete Dictionary. It seems to me that the only difference is in the generative model, and that maximizing the same objective under different data conditions could return an ICA basis or a Complete Dictionary. The robustness theory and experiments on synthetic data are reasonable and demonstrate that complete dictionary learning is robust to the different noise conditions. I would like to how this technique compares to other complete dictionary learning algorithms (ER-SpUD, Complete dictionary learning over the sphere - Sun, Qu, Wright 2015) and whether the l4 objective is unique in providing this robustness. Another central claim of Zhai et. al. 2019 seems to be that l4 maximization is able to recover the entire dictionary at once, vs other algorithms that recover the dictionary one column at a time. To test this, I would like to see runtime evaluations and comparisons to other algorithms. While the claim of recovering the entire dictionary is true, it seems to me that requiring an SVD at each iteration would be very expensive. I am not completely convinced that the approach of estimating the entire dictionary would indeed be faster. To summarize, I believe this paper would be a good addition to the literature on l4 maximization algorithms for dictionary learning. I am willing to adjust my score based on responses to the above concerns.