Exact Recovery of Sparse-Used Dictionaries

1 Introduction

In the Sparsely-Used Dictionary Learning Problem [M. Aharon and Bruckstein, 2006; Georgiev et al., 2005; Vainsencher et al., 2011], one is given a matrix Y ∈ Rn×p and asked to find a pair of matrices A ∈ Rn×m and X ∈ Rm×p so that ‖Y -AX‖ is small and so that X is sparse - X has only a few nonzero elements. We examine solutions to this problem in which A is a basis, so m = n, and without the presence of noise, in which case we insist Y = AX. Variants of this problem arise in different contexts in machine learning, signal processing, and even computational neuroscience [Olshausen and Field, 1996; Engan et al., 1999; Aharon et al., 2006; Mairal et al., 2009].

In this work [Spielman et al., 2012], we contribute to the understanding of the dictionary learning problem in the case when A is square and nonsingular. We prove that O(n log n) samples are sufficient to uniquely determine the decomposition with high probability, under the assumption X is generated by a Bernoulli-Gaussian or Bernoulli-Rademacher process.

Our argument for uniqueness suggests a new, efficient dictionary learning algorithm, which we call Exact Recovery of Sparsely-Used Dictionaries (ER-SpUD). This algorithm solves a sequence of linear programs with varying constraints.

∗The paper on which this extended abstract is based was the recipient of the best paper award of the 2012 Conference On Learning Theory [Spielman et al., 2012].

We prove that under the aforementioned assumptions, the algorithm exactly recovers A and X with high probability. This result holds when the expected number of nonzero elements in each column of X is at most O(√n) and the number of samples p is at least Ω(n2 log2 n). To the best of our knowledge, this result is the first to demonstrate an efficient algorithm for dictionary learning with provable guarantees.

Moreover, we prove that this result is tight to within a log factor: for the Bernoulli-Gaussian case, when the expected number of nonzeros in each column is Ω(√n logn), algorithms of this style fail with high probability.

Our algorithm is related to previous proposals by Zibulevsky and Pearlmutter [Zibulevsky and Pearlmutter, 2001] (for source separation) and Gottlieb and Neylon [Gottlieb and Neylon, 2010] (for dictionary learning), but involves several new techniques that seem to be important for obtaining provable correct recovery - in particular, the use of sample vectors in the constraints. Other related recent proposals include [Plumbley, 2007; Jaillet et al., 2010].

2 Notation

We write ‖v‖p for the standard ℓp norm of a vector v, and we write ‖M‖p for the induced operator norm on a matrix M. ‖v‖0 denotes the number of non-zero entries in v. We denote the Hadamard (point-wise) product by ⊙. [n] denotes the first n positive integers, {1, 2, . . . , n}. For a set of indices I, we let PI denote the projection matrix onto the subspace of vectors supported on indices I, zeroing out the other coordinates. For a matrix X and a set of indices J, we let XJ (XJ ) denote the submatrix containing just the rows (columns) indexed by J. We write the standard basis vector that is non-zero in coordinate i as ei. For a matrix X we let row(X) denote the span of its rows. For a set S, |S| is its cardinality.

3 The Probabilistic Models

We analyze the dictionary learning problem under the assumption that A is an arbitrary nonsingular n-by-n matrix, but that X is a random sparse n-by-p matrix with i.i.d. entries. In the Bernoulli(θ)-Gaussian model, the entries of X are independent random variables, each of which has the form Xi,j = ςτ , where ς ∼ N(0, 1) is a standard Gaussian, and τ is 1 with probability θ and 0 with probability 1 - θ, independent of ς.

We also consider a Bernoulli(θ)-Rademacher model, in which the non-zero entries are chosen uniformly in ±1.

4 When is the Factorization Unique?

At first glance, it seems the number of samples p required to identify A could be quite large. For example, Aharon et. al. view the given data matrix X as having sparse columns, each with at most k nonzero entries. If the given samples yj = Axj lie on an arrangement of (nk) k-dimensional subspaces range(AI), corresponding to possible support sets I, A is identifiable.

On the other hand, the most immediate lower bound on the number of samples required comes from the simple fact that to recover A we need to see at least one linear combination involving each of its columns. The "coupon collection" phenomenon tells us that p = Ω( 1 θ logn) samples are required for this to occur with constant probability, where θ is the probability that an element Xij is nonzero. When θ is as small as O(1/n), this means p must be at least proportional to n log n. Our next result shows that, in fact, this lower bound is tight - the problem becomes well-posed once we have observed cn log n samples.

Theorem 4.1 (Uniqueness) Under the Bernoulli(θ)- Gaussian and Bernoulli(θ)-Rademacher models, if 1/n ≤ θ ≤ 1/C and p > Cn logn, then with probability at least 1-exp{-c′p}, for any alternative factorization Y = A′X′ such that maxi ‖eTi X′‖0 ≤ maxi ‖eTi X‖0, we have A′ = AΠΛ and X′ = Λ-1ΠT X, for some permutation matrix Π and nonsingular diagonal matrix Λ, for some absolute constants C and c′ .

5 Exact Recovery

Theorem 4.1 suggests that we can recover X by looking for sparse vectors in the row space of Y . Any vector in this space can be generated by taking a linear combination wT Y of the rows of Y (here, wT denotes the vector transpose). We arrive at the optimization problem

Theorem 4.1 implies that any solution to this problem must satisfy wT Y = λeTj X for some j ∈ [n], λ ̸= 0. Unfortunately, both the objective and constraint are nonconvex. We therefore replace the ℓ0 norm with its convex envelope, the ℓ1 norm, and prevent w from being the zero vector by constraining it to lie in an affine hyperplane {rT w = 1}. This gives a linear programming problem of the form

Or equivalently under the change of variables z = AT w, b = A-1r,

ER-SpUD(SC): Exact Recovery of Sparsely-Used Dictionaries using single columns of Y as constraint vectors. For j = 1 . . . p Solve minw ‖wT Y ‖1s.t. (Y ej)T w = 1, and set sj =wT Y .

5.1 The Algorithms

ER-SpUD(DC): Exact Recovery of Sparsely-Used Dictionaries using the sum of two columns of Y as constraint vectors. 1. Randomly pair the columns of Y into p/2 groups gj = {Y ej1 , Y ej2 }. 2. For j = 1 . . . p/2 Let rj = Y ej1 + Y ej2 , where gj = {Y ej1 , Y ej2 }. Solveminw ‖wT Y ‖1 subject to rTj w = 1, and set sj =wT Y .

Our algorithms are divided into two stages. In the first stage, we collect many potential rows of X by solving problems of the form (1). In the simpler Algorithm ER-SpUD(SC) ("single column"), we do this by using each column of Y as the constraint vector r in the optimization. In the slightly better Algorithm ER-SpUD(DC) ("double column"), we pair up all the columns of Y and then substitue the sum of each pair for r. In the second stage, we use a greedy algorithm (Algorithm Greedy) to select a subset of n of the rows produced. In particular, we choose a linearly independent subset among those with the fewest non-zero elements. From the proof of the uniqueness of the decomposition, we know with high probability that the rows of X are the sparsest n vectors in row(Y ). Moreover, for p ≥ Ω(n logn), Theorems 6.1 and 6.2, along with the coupon collection phenomenon, tell us that a scaled multiple of each of the rows of X is returned by the first phase of our algorithm, with high probability.

Greedy: A Greedy Algorithm to Reconstruct X and A. 1. REQUIRE: S = {s1, . . . , sT } ⊂ Rp . 2. For i = 1 . . . n REPEAT l ← argminsl∈S ‖sl‖0, breaking ties arbitrarily xi = sl S = S\{sl} UNTIL rank([x1, . . . , xi])= i 3. Set X = [x1, . . . , xn]T , and A = Y Y T (XY T )-1 .

6 Main Theoretical Results

The intuitive explanations in the previous section can be made rigorous. In particular, under our random models, we can prove that when the number of samples is reasonably large compared to the dimension, (say p ∼ n2 log2 n), with high probability in X the algorithm will succeed. We conjecture it is possible to decrease the dependency on p to O(n log n).

Theorem 6.1 (Correct recovery (single-column))

Suppose X is Bernoulli(θ)-Gaussian. Then provided p > c1n2 log2 n, and

with probability at least 1 - cfp-10 , the Algorithm ERSpUD(SC) recovers all n rows of X. That is, all n rows of X are included in the p potential vectors wT1 Y , . . . , wTp Y . Above, c1, α and cf are positive numerical constants.

The upper bound of α/√n log n on θ has two sources: an upper bound of α/√n is imposed by the requirement that b be sparse. An additional factor of log n comes from the need for a gap between |b|(1) and |b|(2) of the k i.i.d. Gaussian random variables. On the other hand, using the sum of two columns of Y as r can save the factor of log n in the requirement on θ since the "collision" of non-zero entries in the two columns of X creates a larger gap between |b|(1) and |b|(2). More importantly, the resulting algorithm is less dependent on the magnitudes of the nonzero elements in X. The algorithm using a single column exploited the fact that there exists a reasonable gap between |b|(1) and |b|(2), whereas the twocolumn variant ER-SpUD(DC) succeeds even if the nonzeros all have the same magnitude.

Theorem 6.2 (Correct recovery (two-column)) Suppose

X is Bernoulli(θ)-Gaussian or Bernoulli(θ)-Rademacher. For some α > 0 and for all n larger than some n0, and p > c1n2 log2 n, if the probability of non-zero entries θ satisfies

Then with overwhelming probability, the Algorithm ERSpUD(DC) recovers all n rows of X. That is, all n rows of X are included in the p/2 potential vectors wT1 Y , . . . , wTp/2Y .

Hence, as we choose p to grow faster than n2 log2 n, the algorithm will succeed with probability approaching one. That the algorithm succeeds is interesting, perhaps even unexpected. There is potentially a great deal of symmetry in the problem - all of the rows of X might have similar ℓ1 -norm. The vectors r break this symmetry, preferring one particular solution at each step, at least in the regime where X is sparse. To be precise, the expected number of nonzero entries in each column must be bounded by √n log n.

It is natural to wonder whether this is an artifact of the analysis, or whether such a bound is necessary. We can prove that for Algorithm ER-SpUD(DC), the sparsity demands in Theorem 6.2 cannot be improved by more than a factor of √log n. Consider the optimization problem (2). One can show that for each i, ‖eTi X‖1 ≈ θp. Hence, if we set z = ej⋆ /bj⋆ , where j⋆ is the index of the largest element of b in magnitude, then

If we consider the alternative solution v = sign(b)/‖b‖1, a calculation shows that

Hence, if θ > c√log n/n for sufficiently large c, the second solution will have smaller objective function. These calculations are carried through rigorously in the full version, giving:

Theorem 6.3 For any fixed β and sufficiently large n, and p ≥ C(β)n, the following occurs. If the probability of nonzeros θ satisfies

then the probability (in X) that solving the optimization problem (1) with r = Y ei or r = Y ei +Y ej recovers one of the rows of X is at most n-c(β) , where c(β) > 0.

This implies that the result in Theorem 6.1 is nearly the best possible for this algorithm, at least in terms of its demands on θ.

7 Simulations

In this section we systematically evaluate our algorithm, and compare it with the state-of-the-art dictionary learning algorithms, including K-SVD [Aharon et al., 2006], online dictionary learning [Mairal et al., 2009], SIV [Gottlieb and Neylon, 2010], and the relative Newton method for source separation [Zibulevsky, 2003]. The first two methods are not limited to square dictionaries, while the final two methods, like ours, exploit properties of the square case. The method of [Zibulevsky, 2003] is similar in provenance to the incremental nonconvex approach of [Zibulevsky and Pearlmutter, 2001], but seeks to recover all of the rows of X simultaneously, by seeking a local minimum of a larger nonconvex problem. We found in the experiments that a slight variant of the greedy ER-SPUD algorithm, we call the ER-SPUD(proj)[?], works even better than the greedy scheme. And thus we also add its result to the comparison list. As our emphasis in this paper is mostly on correctness of the solution, we modify the default settings of these packages to obtain more accurate results (and hence a fairer comparison). For K-SVD, we use high accuracy mode, and switch the number of iterations from 10 to 30. Similarly, for relative Newton, we allow 1,000 iterations. For online dictionary learning, we allow 1,000. We observed diminishing returns beyond these numbers. Since K-SVD and online dictionary learning tend to get stuck at local optimum, for each trial we restart K-SVD and Online learning algorithm 5 times with randomized initializations and report the best performance. We measure accuracy in terms of the relative error, after permutation-scale ambiguity has been removed:

Phase transition graph. In our experiments we have chosen A to be a an n-by-n matrix of independent Gaussian random variables. The coefficient matrix X is n-by-p, where p = 5n loge n. Each column of X has k randomly chosen nonzero entries. In our experiments we have varied n between 10 and 60 and k between 1 and 10. Figure 1 shows the results for each method, with the average relative error reported in greyscale. White means zero error and black is 1. The best performing algorithm is ER-SpUD with iterative projections, which solves almost all the cases except when n = 10 and k ≥ 6. For the other algorithm, When n is small, the relative Newton method appears to be able to handle a denser X, while as n grows large, the greedy ER-SpUD is more precise. In fact, empirically the phase transition between success and failure for ER-SpUD is quite sharp - problems below the boundary are solved to high numerical accuracy, while beyond the boundary the algorithm breaks down. In contrast, both online dictionary learning and relative Newton exhibit neither the same accuracy, nor the same sharp transition to failure - even in the black region of the graph, they still return solutions that are not completely wrong. The breakdown boundary of K-SVD is clear compared to online learning and relative Newton. As an active set algorithm, when it reaches a correct solution, the numerical accuracy is quite high. However, in our simulations we observe that both K-SVD and online learning may be trapped into a local optimum even for relatively sparse problems.

Acknowledgments

This material is based in part upon work supported by the National Science Foundation under Grant No. 0915487. JW also acknowledges support from Columbia University.