Simple, Efficient, and Neural Algorithms for Sparse Coding

Introduction

Sparse coding or dictionary learning consists of learning to express (i.e., code) a set of input vectors, say image patches, as linear combinations of a small number of vectors chosen from a large dictionary. It is a basic task in many fields. In signal processing, a wide variety of signals turn out to be sparse in an appropriately chosen basis (see references in Mallat (1998)). In neuroscience, sparse representations are believed to improve energy efficiency of the brain by allowing most neurons to be inactive at any given time. In machine learning, imposing sparsity as a constraint on the representation is a useful way to avoid over-fitting. Additionally, methods for sparse coding can be thought of as a tool for feature extraction and are the basis for a number of important tasks in image processing such as segmentation, retrieval, de-noising and super-resolution (see references in Elad (2010)), as well as a building block for some deep learning architectures Ranzato et al. (2007). It is also a basic problem in linear algebra itself since it involves finding a better basis.

The notion was introduced by neuroscientists Olshausen and Field (1997a) who formalized it as follows: Given a dataset $y^{(1)},y^{(2)},\ldots,y^{(p)}\in\Re^{n}$, our goal is to find a set of basis vectors $A_{1},A_{2},\ldots,A_{m}\in\Re^{n}$ and sparse coefficient vectors $x^{(1)},x^{(2)},\ldots,x^{(p)}\in\Re^{m}$ that minimize the reconstruction error

where $A$ is the $n\times m$ coding matrix whose $j$th column is $A_{j}$ and $S(\cdot)$ is a nonlinear penalty function that is used to encourage sparsity. This function is nonconvex because both $A$ and the $x^{(i)}$’s are unknown. Their paper as well as subsequent work chooses $m$ to be larger than $n$ (so-called overcomplete case) because this allows greater flexibility in adapting the representation to the data. We remark that sparse coding should not be confused with the related — and usually easier — problem of finding the sparse representations of the $y^{(i)}$’s given the coding matrix $A$, variously called compressed sensing or sparse recovery Candes et al. (2006); Candes and Tao (2005).

Olshausen and Field also gave a local search/gradient descent heuristic for trying to minimize the nonconvex energy function (1). They gave experimental evidence that it produces coding matrices for image patches that resemble known features (such as Gabor filters) in $V1$ portion of the visual cortex. A related paper of the same authors Olshausen and Field (1997b) (and also Lewicki and Sejnowski (2000)) places sparse coding in a more familiar generative model setting whereby the data points $y^{(i)}$’s are assumed to be probabilistically generated according to a model $y^{(i)}=A^{*}\cdot{x^{*(i)}}+\mbox{noise}$ where ${x^{*(1)}},{x^{*(2)}},\ldots,{x^{*(p)}}$ are samples from some appropriate distribution and $A^{*}$ is an unknown code. Then one can define the maximum likelihood estimate, and this leads to a different and usually more complicated energy function — and associated heuristics — compared to (1).

Surprisingly, maximum likelihood-based approaches seem unnecessary in practice and local search/gradient descent on the energy function (1) with hard constraints works well, as do related algorithms such as MOD Aharon et al. (2006) and $k$-SVD Engan et al. (1999). In fact these methods are so effective that sparse coding is considered in practice to be a solved problem, even though it has no polynomial time algorithm per se.

Recently, there has been rapid progress on designing polynomial time algorithms for sparse coding with provable guarantees (the relevant papers are discussed below). All of these adopt the generative model viewpoint sketched above. But the surprising success of the simple descent heuristics has remained largely unexplained. Empirically, these heuristics far out perform — in running time, sample complexity, and solution quality — the new algorithms, and this (startling) observation was in fact the starting point for the current work.

Of course, the famous example of simplex vs ellipsoid for linear programming reminds us that it can be much more challenging to analyze the behavior of an empirically successful algorithm than it is to to design a new polynomial time algorithm from scratch! But for sparse coding the simple intuitive heuristics are important for another reason beyond just their algorithmic efficiency: they appear to be implementable in neural architectures. (Roughly speaking, this means that the algorithm stores the code matrix $A$ as synapse weights in a neural network and updates the entries using differences in potentials of the synapse’s endpoints.) Since neural computation — and also deep learning — have proven to be difficult to analyze in general, analyzing sparse coding thoroughly seems to be a natural first step for theory. Our algorithm is a close relative of the Olshausen-Field algorithm and thus inherits its neural implementability; see Appendix E for further discussion.

Here we present a rigorous analysis of the simple energy minimization heuristic, and as a side benefit this yields bounds on running time and sample complexity for sparse coding that are better (in some cases, dramatically so) than the algorithms in recent papers. This adds to the recent literature on analyzing alternating minimization Jain et al. (2013); Hardt (2013); Netrapalli et al. (2013, 2014) but these work in a setting where there is a convex program that is known to work too, and in our setting, the only known convex program runs in time exponential in a natural parameter Barak et al. (2014).

1 Recent Work

A common thread in recent work on sparse coding is to assume a generative model; the precise details vary, but each has the property that given enough samples the solution is essentially unique. Spielman et al. (2012) gave an algorithm that succeeds when $A^{*}$ has full column rank (in particular $m\leq n$) which works up to sparsity roughly $\sqrt{n}$. However this algorithm is not applicable in the more prevalent overcomplete setting. Arora et al. (2014) and Agarwal et al. (2013, 2014) independently gave algorithms in the overcomplete case assuming that $A^{*}$ is $\mu$-incoherent (which we define in the next section). The former gave an algorithm that works up to sparsity $n^{1/2-\gamma}/\mu$ for any $\gamma>0$ but the running time is $n^{\Theta(1/\gamma)}$; Agarwal et al. (2013, 2014) gave an algorithm that works up to sparsity either $n^{1/4}/\mu$ or $n^{1/6}/\mu$ depending on the particular assumptions on the model. These works also analyze alternating minimization but assume that it starts from an estimate $A$ that is column-wise $1/\mbox{poly}(n)$-close to $A^{*}$, in which case the objective function is essentially convex.

Barak et al. (2014) gave a new approach based on the sum-of-squares hierarchy that works for sparsity up to $n^{1-\gamma}$ for any $\gamma>0$. But in order to output an estimate that is column-wise $\epsilon$-close to $A^{*}$ the running time of the algorithm is $n^{1/\epsilon^{O(1)}}$. In most applications, one needs to set (say) $\epsilon=1/k$ in order to get a useful estimate. However in this case their algorithm runs in exponential time. The sample complexity of the above algorithms is also rather large, and is at least $\Omega(m^{2})$ if not much larger. Here we will give simple and more efficient algorithms based on alternating minimization whose column-wise error decreases geometrically, and that work for sparsity up to $n^{1/2}/\mu\log n$. We remark that even empirically, alternating minimization does not appear to work much beyond this bound.

2 Model, Notation and Results

We will work with the following family of generative models (similar to those in earlier papers)The casual reader should just think of $x^{*}$ as being drawn from some distribution that has independent coordinates. Even in this simpler setting —which has polynomial time algorithms using Independent Component Analysis—we do not know of any rigorous analysis of heuristics like Olshausen-Field. The earlier papers were only interested in polynomial-time algorithms, so did not wish to assume independence.:

Each sample is generated as $y=A^{*}x^{*}+\mbox{noise}$ where $A^{*}$ is a ground truth dictionary and $x^{*}$ is drawn from an unknown distribution $\mathcal{D}$ where

the support $S=\mbox{supp}(x^{*})$ is of size at most $k$, $\mathop{\bf Pr\/}[i\in S]=\Theta(k/m)$ and $\mathop{\bf Pr\/}[i,j\in S]=\Theta(k^{2}/m^{2})$

the distribution is normalized so that $\mathop{\bf E\/}[x^{*}_{i}|x^{*}_{j}\neq 0]=0$; $\mathop{\bf E\/}[{x^{*}_{i}}^{2}|x^{*}_{i}\neq 0]=1$ and when $x^{*}_{i}\neq 0$, $|x^{*}_{i}|\geq C$ for some constant $C\leq 1$ and

the non-zero entries are pairwise independent and subgaussian, conditioned on the support.

The noise is Gaussian and independent across coordinates.

Such models are natural since the original motivation behind sparse coding was to discover a code whose representations have the property that the coordinates are almost independent. We can relax most of the requirements above, at the expense of further restricting the sparsity, but will not detail such tradeoffs.

The rest of the paper ignores the iid noise: it has little effect on our basic steps like computing inner products of samples or taking singular vectors, and easily tolerated so long as it stays smaller than the “signal.”

We assume $A^{*}$ is an incoherent dictionary, since these are widespread in signal processing Elad (2010) and statistics Donoho and Huo (1999), and include various families of wavelets, Gabor filters as well as randomly generated dictionaries.

An $n\times m$ matrix $A$ whose columns are unit vectors is $\mu$-incoherent if for all $i\neq j$ we have $\langle A_{i},A_{j}\rangle\leq\mu/\sqrt{n}.$

We also require that $\|A^{*}\|=O(\sqrt{m/n})$. However this can be relaxed within polylogarithmic factors by tightening the bound on the sparsity by the same factor. Throughout this paper we will say that $A^{s}$ is $(\delta,\kappa)$-near to $A^{*}$ if after a permutation and sign flips its columns are within distance $\delta$ and we have $\|A^{s}-A^{*}\|\leq\kappa\|A^{*}\|$. See also Definition 3.1. We will use this notion to measure the progress of our algorithms. Moreover we will use $g(n)=O^{*}(f(n))$ to signify that $g(n)$ is upper bounded by $Cf(n)$ for some small enough constant $C$. Finally, throughout this paper we will assume that $k\leq O^{*}(\sqrt{n}/\mu\log n)$ and $m=O(n)$. Again, $m$ can be allowed to be higher by lowering the sparsity. We assume all these conditions in our main theorems.

In Section 2 we give a general framework for analyzing alternating minimization. Instead of thinking of the algorithm as trying to minimize a known non-convex function, we view it as trying to minimize an unknown convex function. Various update rules are shown to provide good approximations to the gradient of the unknown function. See Lemma 3.4, Lemma B.1 and Lemma B.7 for examples. We then leverage our framework to analyze existing heuristics and to design new ones also with provable guarantees. In Section 3, we prove:

There is a neurally plausible algorithm which when initialized with an estimate $A^{0}$ that is $(\delta,2)$-near to $A^{*}$ for $\delta=O^{*}(1/\log n)$, converges at a geometric rate to $A^{*}$ until the column-wise error is $O(\sqrt{k/n})$. Furthermore the running time is $O(mnp)$ and the sample complexity is $p=\widetilde{O}(mk)$ for each step.

Additionally we give a neural architecture implementing our algorithm in Appendix E. To the best of our knowledge, this is the first neurally plausible algorithm for sparse coding with provable convergence.

Having set up our general framework and analysis technique we can use it on other variants of alternating minimization. Section 3.1.2 gives a new update rule whose bias (i.e., error) is negligible:

There is an algorithm which when initialized with an estimate $A^{0}$ that is $(\delta,2)$-near to $A^{*}$ for $\delta=O^{*}(1/\log n)$, converges at a geometric rate to $A^{*}$ until the column-wise error is $O(n^{-\omega(1)})$. Furthermore each step runs in time $O(mnp)$ and the sample complexity $p$ is polynomial.

This algorithm is based on a modification where we carefully project out components along the column currently being updated. We complement the above theorems by revisiting the Olshausen-Field rule and analyzing a variant of it in Section 3.1.1 (Theorem 3.6). However its analysis is more complex because we need to bound some quadratic error terms. It uses convex programming.

What remains is to give a method to initialize these iterative algorithms. We give a new approach based on pair-wise reweighting and we prove that it returns an estimate $A^{0}$ that is $(\delta,2)$-near to $A^{*}$ for $\delta=O^{*}(1/\log n)$ with high probability. As an additional benefit, this algorithm can be used even in settings where $m$ is not known and this could help solve another problem in practice — that of model selection. In Section 5 we prove:

There is an algorithm which returns an estimate $A^{0}$ that is $(\delta,2)$-near to $A^{*}$ for $\delta=O^{*}(1/\log n)$. Furthermore the running time is $\widetilde{O}(mn^{2}p)$ and the sample complexity $p=\widetilde{O}(mk)$.

This algorithm also admits a neural implementation, which is sketched in Appendix E. The proof currently requires a projection step that increases the run time though we suspect it is not needed.

We remark that these algorithms work up to sparsity $O^{*}(\sqrt{n}/\mu\log n)$ which is within a logarithmic factor of the information theoretic threshold for sparse recovery on incoherent dictionaries Donoho and Huo (1999); Gribonval and Nielsen (2003). All previous known algorithms that approach Arora et al. (2014) or surpass this sparsity Barak et al. (2014) run in time exponential in some natural parameter. Moreover, our algorithms are simple to describe and implement, and involve only basic operations. We believe that our framework will have applications beyond sparse coding, and could be used to show that simple, iterative algorithms can be powerful in other contexts as well by suggesting new ways to analyze them.

Our Framework, and an Overview

Here we describe our framework for analyzing alternating minimization. The generic scheme we will be interested in is given in Algorithm 1 and it alternates between updating the estimates $A$ and $X$. It is a heuristic for minimizing the non-convex function in (1) where the penalty function is a hard constraint. The crucial step is if we fix $X$ and compute the gradient of (1) with respect to $A$, we get:

We then take a step in the opposite direction to update $A$. Here and throughout the paper $\eta$ is the learning rate, and needs to be set appropriately. The challenge in analyzing this general algorithm is to identify a suitable “measure of progress”— called a Lyapunov function in dynamical systems and control theory — and show that it improves at each step (with high probability over the samples). We will measure the progress of our algorithms by the maximum column-wise difference between $A$ and $A^{*}$.

In the next subsection, we identify sufficient conditions that guarantee progress. They are inspired by proofs in convex optimization. We view Algorithm 1 as trying to minimize an unknown convex function, specifically $f(A)=\mathcal{E}(A,X^{*})$, which is strictly convex and hence has a unique optimum that can be reached via gradient descent. This function is unknown since the algorithm does not know $X^{*}$. The analysis will show that the direction of movement is correlated with $A^{*}-A^{s}$, which in turn is the gradient of the above function. An independent paper of Balakrishnan et al. (2014) proposes a similar framework for analysing EM algorithms for hidden variable models. The difference is that their condition is really about the geometry of the objective function, though ours is about the property of the direction of movement. Therefore we have the flexibility to choose different decoding procedures. This flexibility allows us to have a closed form of $X^{s}$ and obtain a useful functional form of $g^{s}$. The setup is reminiscent of stochastic gradient descent, which moves in a direction whose expectation is the gradient of a known convex function. By contrast, here the function $f()$ is unknown, and furthermore the expectation of $g^{s}$ is not the true gradient and has bias. Due to the bias, we will only be able to prove that our algorithms reach an approximate optimum up to some error whose magnitude is determined by the bias. We can make the bias negligible using more complicated algorithms.

Consider a general iterative algorithm that is trying to get to a desired solution $z^{*}$ (in our case $z^{*}=A^{*}_{i}$ for some $i$). At step $s$ it starts with a guess $z^{s}$, computes some direction $g^{s}$, and updates its estimate as: $z^{s+1}=z^{s}-\eta g^{s}$. The natural progress measure is $\|z^{*}-z^{s}\|^{2}$, and below we will identify a sufficient condition for it to decrease in each step:

A vector $g^{s}$ is $(\alpha,\beta,\epsilon_{s})$-correlated with $z^{*}$ if

Remark: The traditional analysis of convex optimization corresponds to the setting where $z^{*}$ is the global optimum of some convex function $f$, and $\epsilon_{s}=0$. Specifically, if $f(\cdot)$ is $2\alpha$-strongly convex and $1/(2\beta)$-smooth, then $g^{s}=\nabla f(z^{s})$ $(\alpha,\beta,0)$-correlated with $z^{*}$. Also we will refer to $\epsilon_{s}$ as the bias.

Suppose $g^{s}$ satisfies Definition 2.1 for $s=1,2,\dots,T$, and $\eta$ satisfies $0<\eta\leq 2\beta$ and $\epsilon=\max_{s=1}^{T}\epsilon_{s}$. Then for any $s=1,\dots,T$,

In particular, the update rule above converges to $z^{*}$ geometrically with systematic error $\epsilon/\alpha$ in the sense that

Furthermore, if $\epsilon_{s}<\frac{\alpha}{2}\|z^{s}-z^{*}\|^{2}$ for $s=1,\dots,T$, then

The proof closely follows existing proofs in convex optimization:

Then solving this recurrence we have $\|z^{s+1}-z^{*}\|^{2}\leq(1-2\alpha\eta)^{s+1}R^{2}+\frac{\epsilon}{\alpha}$ where $R=\|z^{0}-z^{*}\|$. And furthermore if $\epsilon_{s}<\frac{\alpha}{2}\|z^{s}-z^{*}\|^{2}$ we have instead

and this yields the second part of the theorem too.

In fact, we can extend the analysis above to obtain identical results for the case of constrained optimization. Suppose we are interested in optimizing a convex function $f(z)$ over a convex set $\mathcal{B}$. The standard approach is to take a step in the direction of the gradient (or $g^{s}$ in our case) and then project into $\mathcal{B}$ after each iteration, namely, replace $z^{s+1}$ by $\mbox{Proj}_{\mathcal{B}}\,z^{s+1}$ which is the closest point in $\mathcal{B}$ to $z^{s+1}$ in Euclidean distance. It is well-known that if $z^{*}\in\mathcal{B}$, then $\|\mbox{Proj}_{\mathcal{B}}\;z-z^{*}\|\leq\|z-z^{*}\|$. Therefore we obtain the following as an immediate corollary to the above analysis:

Suppose $g^{s}$ satisfies Definition 2.1 for $s=1,2,\dots,T$ and set $0<\eta\leq 2\beta$ and $\epsilon=\max_{s=1}^{T}\epsilon_{s}$. Further suppose that $z^{*}$ lies in a convex set $\mathcal{B}$. Then the update rule $z^{s+1}=\mbox{Proj}_{\mathcal{B}}(z^{s}-\eta g^{s})$ satisfies that for any $s=1,\dots,T$,

In particular, $z^{s}$ converges to $z^{*}$ geometrically with systematic error $\epsilon/\alpha$. Additionally if $\epsilon_{s}<\frac{\alpha}{2}\|z^{s}-z^{*}\|^{2}$ for $s=1,\dots,T$, then

What remains is to derive a functional form for various update rules and show that these rules move in a direction $g^{s}$ that approximately points in the direction of the desired solution $z^{*}$ (under the assumption that our data is generated from a stochastic model that meets certain conditions).

An Overview of Applying Our Framework

Our framework clarifies that any improvement step meeting Definition 2.1 will also converge to an approximate optimum, which enables us to engineer other update rules that turn out to be easier to analyze. Indeed we first analyze a simpler update rule with $g^{s}=\mathop{\bf E\/}[(y-A^{s}x)\mbox{sgn}(x^{T})]$ in Section 3. Here $\mbox{sgn}(\cdot)$ is the coordinate-wise sign function. We then return to the Olshausen-Field update rule and analyze a variant of it in Section 3.1.1 using approximate projected gradient descent. Finally, we design a new update rule in Section 3.1.2 where we carefully project out components along the column currently being updated. This has the effect of replacing one error term with another and results in an update rule with negligible bias. The main steps in showing that these update rules fit into our framework are given in Lemma 3.4, Lemma B.1 and Lemma B.7.

The rest of the analysis can be described as follows: If the signs and support of $x$ are recovered correctly, then alternating minimization makes progress in each step. In fact this holds each for much larger values of $k$ than we consider; as high as $n/(\log n)^{O(1)}$. (However, the explicit decoding rule fails for $k>\sqrt{n}/\mu\log n$.) Thus it only remains to properly initialize $A^{0}$ so that it is close enough to $A^{*}$ to let the above decoding rule succeed. In Section 5 we give a new initialization procedure based on pair-wise reweighting that we prove works with high probability. This section may be of independent interest, since this algorithm can be used even in settings where $m$ is not known and could help solve another problem in practice — that of model selection. See Lemma 5.2.

A Neurally Plausible Algorithm with Provable Guarantees

Here we will design and analyze a neurally plausible algorithm for sparse coding which is given in Algorithm 2, and we give a neural architecture implementing our algorithm in Appendix E. The fact that such a simple algorithm provably works sheds new light on how sparse coding might be accomplished in nature. Here and throughout this paper we will work with the following measure of closeness:

$A$ is $\delta$-close to $A^{*}$ if there is a permutation $\pi:[m]\rightarrow[m]$ and a choice of signs $\sigma:[m]\rightarrow\{\pm 1\}$ such that \|\sigma(i)A_{\pi(i)}-A^{*}_{i}\|\leq\delta\mbox{ for alli} We say $A$ is $(\delta,\kappa)$-near to $A^{*}$ if in addition $\|A-A^{*}\|\leq\kappa\|A^{*}\|$ too.

This is a natural measure to use, since we can only hope to learn the columns of $A^{*}$ up to relabeling and sign-flips. In our analysis, we will assume throughout that $\pi(\cdot)$ is the identity permutation and $\sigma(\cdot)\equiv+1$ because our family of generative models is invariant under this relabeling and it will simplify our notation.

Let $\mbox{sgn}(\cdot)$ denote the coordinate-wise sign function and recall that $\eta$ is the learning rate, which we will soon set. Also we fix both $\delta,\delta_{0}=O^{*}(1/\log n)$. We will also assume that in each iteration, our algorithm is given a fresh set of $p$ samples. Our main theorem is:

Suppose that $A^{0}$ is $(2\delta,2)$-near to $A^{*}$ and that $\eta=\Theta(m/k)$. Then if each update step in Algorithm 2 uses $p=\widetilde{\Omega}(mk)$ fresh samples, we have

for some $0<\tau<1/2$ and for any $s=1,2,...,T$. In particular it converges to $A^{*}$ geometrically, until the column-wise error is $O(\sqrt{k/n})$.

Our strategy is to prove that $\widehat{g}^{s}$ is $(\alpha,\beta,\epsilon)$-correlated (see Definition 2.1) with the desired solution $A^{*}$, and then to prove that $\|A\|$ never gets too large. We will first prove that if $A$ is somewhat close to $A^{*}$ then the estimate $x$ for the representation almost always has the correct support. Here and elsewhere in the paper, we use “very high probability” to mean that an event happens with probability at least $1-1/n^{\omega(1)}$.

Suppose that $A^{s}$ is $\delta$-close to $A^{*}$. Then with very high probability over the choice of the random sample $y=A^{*}x^{*}$:

We prove a more general version of this lemma (Lemma A.1) in Appendix A; it is an ingredient in analyzing all of the update rules we consider in this paper. However this is just one step on the way towards proving that $\widehat{g}^{s}$ is correlated with the true solution.

The next step in our proof is to use the properties of the generative model to derive a new formula for $\widehat{g}^{s}$ that is more amenable to analysis. We define $g^{s}$ to be the expectation of $\widehat{g}^{s}$

where $x:=\textrm{threshold}_{C/2}((A^{s})^{T}y)$ is the decoding of $y$. Let $q_{i}=\mathop{\bf Pr\/}[x^{*}_{i}\neq 0]$ and $q_{i,j}=\mathop{\bf Pr\/}[x^{*}_{i}x^{*}_{j}\neq 0]$, and define $p_{i}=\mathop{\bf E\/}[x^{*}_{i}\mbox{sgn}(x^{*}_{i})|x^{*}_{i}\neq 0]$.

Here and in the rest of the paper, we will let $\gamma$ denote any vector whose norm is negligible (i.e. smaller than $1/n^{C}$ for any large constant $C>1$). This will simplify our calculations. Also let $A^{*}_{-i}$ denote the matrix obtained from deleting the $i$th column of $A^{*}$. The following lemma is the main step in our analysis.

Suppose that $A^{s}$ is $(2\delta,2)$-near to $A^{*}$. Then the update step in Algorithm 2 takes the form $\mathop{\mathbf{E}}\displaylimits[A_{i}^{s+1}]=A_{i}^{s}-\eta g_{i}^{s}$ where $g_{i}^{s}=p_{i}q_{i}\left(\lambda^{s}_{i}A_{i}^{s}-A^{*}_{i}+\epsilon^{s}_{i}\pm\gamma\right)$, and $\lambda^{s}_{i}=\langle A_{i}^{s},A^{*}_{i}\rangle$ and

Moreover the norm of $\epsilon^{s}_{i}$ can be bounded as $\|\epsilon^{s}_{i}\|\leq O(k/n)$.

Note that $p_{i}q_{i}$ is a scaling constant and $\lambda_{i}\approx 1$; hence from the above formula we should expect that $g_{i}^{s}$ is well-correlated with $A_{i}^{s}-A^{*}_{i}$.

Since $A^{s}$ is $(2\delta,2)$-near to $A^{*}$, $A^{s}$ is $2\delta$-close to $A^{*}$. We can now invoke Lemma 3.3 and conclude that with high probability, $\mbox{sgn}(x^{*})=\mbox{sgn}(x)$. Let $\mathcal{F}_{x^{*}}$ be the event that $\mbox{sgn}(x^{*})=\mbox{sgn}(x)$, and let $\mathbf{1}_{\mathcal{F}_{x^{*}}}$ be the indicator function of this event.

The key is that this allows us to essentially replace $\mbox{sgn}(x)$ with $\mbox{sgn}(x^{*})$. Moreover, let $S=\mbox{supp}(x^{*})$. Note that when $\mathcal{F}_{x^{*}}$ happens $S$ is also the support of $x$. Recall that according to the decoding rule (where we have replaced $A^{s}$ by $B$ for notational simplicity) $x=\textrm{threshold}_{C/2}(B^{T}y)$. Therefore, $x_{S}=(B^{T}y)_{S}=B_{S}^{T}y=B_{S}^{T}A^{*}x^{*}$. Using the fact that the support of $x$ is $S$ again, we have $Bx=B_{S}^{T}B_{S}A^{*}x^{*}$. Plugging it into equation (3):

where again we have used the fact that $\mathcal{F}_{x^{*}}$ happens with very high probability. Now we rewrite the expectation above using subconditioning where we first choose the support $S$ of $x^{*}$, and then we choose the nonzero values $x^{*}_{S}$.

where we use the fact that $\mathop{\bf E\/}[x^{*}_{i}\cdot\mbox{sgn}(x^{*}_{i})|S]=p_{i}$. Let $R=S-\{i\}$. Using the fact that $B_{S}B_{S}^{T}=B_{i}B_{i}^{T}+B_{R}B_{R}^{T}$, we can split the quantity above into two parts,

where $\mbox{diag}(q_{i,j})$ is a $m\times m$ diagonal matrix whose $(j,j)$-th entry is equal to $q_{i,j}$, and $B_{-i}$ is the matrix obtained by zeroing out the $i$th column of $B$. Here we used the fact that $\mathop{\bf Pr\/}[i\in S]=q_{i}$ and $\mathop{\bf Pr\/}[i,j\in S]=q_{ij}$.

Now we set $B=A^{s}$, and rearranging the terms, we have $g_{i}^{s}=p_{i}q_{i}\left(\langle A_{i}^{s},A^{*}_{i}\rangle A_{i}^{s}-A^{*}_{i}+\epsilon^{s}_{i}\pm\gamma\right)$ where \epsilon_{i}^{s}=\Big{(}A^{s}_{-i}\mbox{diag}(q_{i,j})\left(A_{-i}^{s}\right)^{T}\Big{)}A^{*}_{i}/q_{i}, which can be bounded as follows

where the last step used the fact that $\frac{\max_{i\neq j}q_{i,j}}{\min q_{i}}\leq O(k/m)$, which is an assumption of our generative model.

We will be able to reuse much of this analysis in Lemma B.1 and Lemma B.7 because we have derived to for a general decoding matrix $B$.

In Section 4 we complete the analysis of Algorithm 2 in the infinite sample setting. In particular, in Section 4.1, we prove that if $A^{s}$ is $(2\delta,2)$-near to $A^{*}$ then $g_{i}^{s}$ is indeed $(\alpha,\beta,\epsilon)$-correlated with $A_{i}$ (Lemma 4.3). Finally we prove that if $A^{s}$ is $(2\delta,2)$-near to $A^{*}$ then $\|A^{s+1}-A^{*}\|\leq 2\|A^{*}\|$ (Lemma 4.6). These lemmas together with Theorem 2.2 imply Theorem 4.1, the simplified version of Theorem 3.2 where the number of samples $p$ is assumed to be infinite (i.e. we have access to the true expectation $g^{s}$). In Appendix D we prove the sample complexity bounds we need and this completes the proof of Theorem 3.2.

Here we apply our framework to design and analyze further variants of alternating minimization.

In this subsection we analyze a variant of the Olshausen-Field update rule. However there are quadratic error terms that arise in the expressions we derive for $g^{s}$ and bounding them is more challenging. We will also need to make (slightly) stronger assumptions on the distributional model that for distinct $i_{1},i_{2},i_{3}$ we have $q_{i_{1},i_{2},i_{3}}=O(k^{3}/m^{3})$ where $q_{i_{1},i_{2},i_{3}}=\mathop{\bf Pr\/}[i_{1},i_{2},i_{3}\in S]$.

Suppose that $A^{0}$ is $(2\delta,2)$-near to $A^{*}$ and that $\eta=\Theta(m/k)$. There is a variant of Olshausen-Field (given in Algorithm 4 in Appendix B.1) for which at each step $s$ we have

for some $0<\tau<1/2$ and for any $s=1,2,...,T$. In particular it converges to $A^{*}$ geometrically until the error in Frobenius norm is $O(\sqrt{m}k/n)$.

We defer the proof of the main theorem to Appendix B.1. Currently it uses a projection step (using convex programming) that may not be needed but the proof requires it.

1.2 Removing the Systemic Error

In this subsection, we design and analyze a new update rule that converges geometrically until the column-wise error is $n^{-\omega(1)}$. The basic idea is to engineer a new decoding matrix that projects out the components along the column currently being updated. This has the effect of replacing a certain error term in Lemma 3.4 with another term that goes to zero as $A$ gets closer to $A^{*}$ (the earlier rules we have analyzed do not have this property).

We will use $B^{(s,i)}$ to denote the decoding matrix used when updating the $i$th column in the $s$th step. Then we set B^{(s,i)}_{i}=A_{i}\mbox{ and }B^{(s,i)}_{j}=\mbox{Proj}_{A_{i}^{\perp}}A_{j}\mbox{ forj\neq i.} Note that $B^{(s,i)}_{-i}$ (i.e. $B^{(s,i)}$ with the $i$th column removed) is now orthogonal to $A_{i}$. We will rely on this fact when we bound the error. We defer the proof of the main theorem to Appendix B.2.

Suppose that $A^{0}$ is $(2\delta,2)$-near to $A^{*}$ and that $\eta=\Theta(m/k)$. There is an algorithm (given in Algorithm 5 given in Appendix B.2) for which at each step $s$, we have

for some $0<\tau<1/2$ and for any $s=1,2,...,T$. In particular it converges to $A^{*}$ geometrically until the column-wise error is $n^{-\omega(1)}$.

Analysis of the Neural Algorithm

In Lemma 3.4 we gave a new (and more useful) expression that describes the update direction under the assumptions of our generative model. Here we will make crucial use of Lemma 3.4 in order to prove that $g^{s}_{i}$ is $(\alpha,\beta,\epsilon)$-correlated with $A_{i}$ (Lemma 4.2). Moreover we use Lemma 3.4 again to show that $\|A^{s+1}-A^{*}\|\leq 2\|A^{*}\|$ (Lemma 4.6). Together, these auxiliary lemmas imply that the column-wise error decreases in the next step and moreover the errors across columns are uncorrelated.

We assume that each iteration of Algorithm 2 takes infinite number of samples, and prove the corresponding simplified version of Theorem 3.2. The proof of this Theorem highlights the essential ideas of behind the proof of the Theorem 3.2, which can be found at Appendix D.

Suppose that $A^{0}$ is $(2\delta,2)$-near to $A^{*}$ and that $\eta=\Theta(m/k)$. Then if each update step in Algorithm 2 uses infinite number of samples at each iteration, we have

for some $0<\tau<1/2$ and for any $s=1,2,...,T$. In particular it converges to $A^{*}$ geometrically until the column-wise error is $O(k/n)$.

The proof is deferred to the end of this section.

In Lemma 3.4 we showed that $g^{s}_{i}=p_{i}q_{i}(\lambda_{i}A_{i}^{s}-A^{*}_{i}+\epsilon^{s}_{i}+\gamma)$ where $\lambda_{i}=\langle A_{i},A^{*}_{i}\rangle$. Here we will prove that $g_{i}^{s}$ is $(\alpha,\beta,\epsilon)$-correlated with $A^{*}_{i}$. Recall that we fixed $\delta=O^{*}(1/\log n)$. The main intuition is that $g^{s}_{i}$ is mostly equal to $p_{i}q_{i}(A_{i}^{s}-A^{*}_{i})$ with a small error term.

If a vector $g^{s}_{i}$ is equal to $4\alpha(A_{i}^{s}-A^{*}_{i})+v$ where $\|v\|\leq\alpha\|A_{i}^{s}-A^{*}_{i}\|+\zeta$, then $g^{s}_{i}$ is $(\alpha,1/100\alpha,\zeta^{2}/\alpha)$-correlated with $A^{*}_{i}$, more specifically,

In particular, $g_{i}^{s}$ is $(\alpha,\beta,\epsilon)$-correlated with $A^{*}_{i}$, where $\alpha=\Omega(k/m)$, $\beta\geq\Omega(m/k)$ and $\epsilon=O(k^{3}/mn^{2})$. We can now apply Theorem 2.2 and conclude that the column-wise error gets smaller in the next step:

If $A^{s}$ is $(2\delta,2)$-near to $A^{*}$ and $\eta\leq\min_{i}(p_{i}q_{i}(1-\delta))=O(m/k)$, then $g_{i}^{s}=p_{i}q_{i}(\lambda_{i}A_{i}^{s}-A^{*}_{i}+\epsilon^{s}_{i}+\gamma)$ is $(\Omega(k/m),\Omega(m/k),O(k^{3}/mn^{2}))$-correlated with $A^{*}_{i}$, and further

Using these bounds, we will show $\langle g_{i},A_{i}-A^{*}_{i}\rangle-\alpha\|A_{i}-A^{*}_{i}\|^{2}-\frac{1}{100\alpha}\|g_{i}\|^{2}+\zeta^{2}/\alpha\geq 0$. Indeed we have

We use the form in Lemma 3.4, $g^{s}_{i}=p_{i}q_{i}(\lambda_{i}A_{i}^{s}-A^{*}_{i}+\epsilon^{s}_{i}+\gamma)$ where $\lambda_{i}=\langle A_{i},A^{*}_{i}\rangle$. We can write $g^{s}_{i}=p_{i}q_{i}(A_{i}^{s}-A^{*}_{i})+p_{i}q_{i}((1-\lambda_{i})A_{i}^{s}+\epsilon^{s}_{i}+\gamma)$, so when applying Lemma 4.2 we can use $4\alpha=p_{i}q_{i}=\Theta(k/m)$ and $v=p_{i}q_{i}((1-\lambda_{i})A_{i}^{s}+\epsilon^{s}_{i}+\gamma)$. The norm of $v$ can be bounded in two terms, the first term $p_{i}q_{i}(1-\lambda_{i})A_{i}^{s}$ has norm $p_{i}q_{i}(1-\lambda_{i})$ which is smaller than $p_{i}q_{i}\|A_{i}^{s}-A^{*}_{i}\|$, and the second term has norm bounded by $\zeta=O(k^{2}/mn)$.

By Lemma 4.2 we know the vector $g^{i}_{s}$ is $(\Omega(k/m),\Omega(m/k),O(k^{3}/mn^{2}))$-correlated with $A^{s}$. Then by Theorem 2.2 we have the last part of the corollary.

2 Maintaining Nearness

Suppose that $A^{s}$ is $(2\delta,2)$-near to $A^{*}$. Then $\|A^{s+1}-A^{*}\|\leq 2\|A^{*}\|$ in Algorithm 2.

As in the proof of the previous lemma, we will make crucial use of Lemma 3.4. Substituting and rearranging terms we have:

Our first goal is to write this equation in a more convenient form. In particular let $U$ and $V$ be matrices such that $U_{i}=p_{i}q_{i}(1-\lambda_{i}^{s})A_{i}^{s}$ and V_{i}=p_{i}\Big{(}A^{s}_{-i}\mbox{diag}(q_{i,j})\left(A_{-i}^{s}\right)^{T}\Big{)}A^{*}_{i}. Then we can re-write the above equation as:

where $\mbox{diag}(1-\eta p_{i}q_{i})$ is the $m\times m$ diagonal matrix whose entries along the diagonal are $1-\eta p_{i}q_{i}$.

We will bound the spectral norm of $A^{s+1}-A$ by bounding the spectral norm of each of the matrices of right hand side. The first two terms are straightforward to bound:

where the last inequality uses the assumption that $p_{i}=\Theta(1)$ and $q_{i}\leq O(k/m)$, and the assumption that $\|A^{s}-A^{*}\|\leq 2\|A^{*}\|$.

From the definition of $U$ it follows that $U=A^{s}\mbox{diag}(p_{i}q_{i}(1-\lambda_{i}^{s}))$, and therefore

where we have used the fact that $\lambda_{i}^{s}\geq 1-\delta$ and $\delta=o(1)$, and $\|A^{s}\|\leq\|A^{s}-A^{*}\|+\|A^{*}\|=O(\|A^{*}\|)$.

What remains is to bound the third term, and let us first introduce an auxiliary matrix $Q$ which we define as follows: $Q_{ii}=0$ and $Q_{i,j}=q_{i,j}\langle A^{s}_{i},A^{*}_{i}\rangle$ for $i\neq j$. It is easy to verify that the following claim:

The $i$th column of $A^{s}Q$ is equal to \Big{(}A^{s}_{-i}\mbox{diag}(q_{i,j})\left(A_{-i}^{s}\right)^{T}\Big{)}A^{*}_{i}

Therefore we can write $V=A^{s}Q\mbox{diag}(p_{i})$. We will bound the spectral norm of $Q$ by bounding its Frobenus norm instead. Then from the definition of $A$, we have that:

Moreover since ${A^{*}}^{T}A^{s}$ is an $m\times m$ matrix, its Frobenius norm can be at most a $\sqrt{m}$ factor larger than its spectral norm. Hence we have

where the last inequality uses the fact that $k=O(\sqrt{n}/\log n)$ and $\|A^{s}\|\leq O(\|A^{*}\|)$.

Therefore, putting the pieces together we have:

and this completes the proof of the lemma.

3 Proof of Theorem 4.1

We prove by induction on $s$. Our induction hypothesis is that the theorem is true at each step $s$ and $A^{s}$ is $(2\delta,2)$-near to $A^{*}$. The hypothesis is trivially true for $s=0$. Now assuming the inductive hypothesis is true. Recall that Corollary 4.3 of Section 4.1 says that if $A^{s}$ is $(2\delta,2)$-near to $A^{*}$, which is guaranteed by the inductive hypothesis, by then $g_{i}^{s}$ is indeed $(\Omega(k/m),\Omega(m/k),O(k^{3}/mn^{2}))$-correlated with $A^{*}_{i}$. Invoking our framework of analysis (Theorem 2.2), we have that

Therefore it also follows that $A^{s+1}$ is $2\delta$-close to $A^{*}$. Then we invoke Lemma 4.6 to prove $A^{s+1}$ has not too large spectral norm $\|A^{s+1}-A^{*}\|\leq 2\|A^{*}\|$, which completes the induction.

Initialization

There is a large gap between theory and practice in terms of how to initialize alternating minimization. The usual approach is to set $A$ randomly or to populate its columns with samples $y^{(i)}$. These often work but we do not know how to analyse them. Here we give a novel method for initialization which we show succeeds with very high probability. Our algorithm works by pairwise reweighting. Let $u=A^{*}\alpha$ and $v=A^{*}\alpha^{\prime}$ be two samples from our model whose supports are $U$ and $V$ respectively. The main idea is that if we reweight fresh samples $y$ with a factor $\langle y,u\rangle\langle y,v\rangle$ and compute

then the top singular vectors will correspond to columns $A^{*}_{j}$ where $j\in U\cap V$. (This is reminiscent of ideas in recent papers on dictionary learning, but more sample efficient.)

Throughout this section we will assume that the algorithm is given two sets of samples of size $p_{1}$ and $p_{2}$ respectively. Let $p=p_{1}+p_{2}$. We use the first set of samples for the pairs $u,v$ that are used in reweighting and we use the second set to compute $\widehat{M}_{u,v}$ (that is, the same set of $p_{2}$ samples is used for each $u,v$ throughout the execution of the algorithm). Our main theorem is:

Suppose that Algorithm 3 is given $p_{1}=\widetilde{\Omega}(m)$ and $p_{2}=\widetilde{\Omega}(mk)$ fresh samples and moreover $(a)$ $A^{*}$ is $\mu$-incoherent with $\mu=O^{*}(\frac{\sqrt{n}}{k\log^{3}n})$, $(b)$ $m=O(n)$ and $(c)$ $\|A^{*}\|\leq O(\sqrt{\frac{m}{n}})$. Then with high probability $A$ is $(\delta,2)$-near to $A^{*}$ where $\delta=O^{*}(1/\log n)$.

We will defer the proof of this theorem to Appendix C. The key idea is the following main lemma. We will invoke this lemma several times in order to analyze Algorithm 3 to verify whether or not the supports of $u$ and $v$ share a common element, and again to show that if they do we can approximately recover the corresponding column of $A^{*}$ from the top singular vector of $M_{u,v}$.

Suppose $u=A^{*}\alpha$ and $v=A^{*}\alpha^{\prime}$ are two random samples with supports $U$, $V$ respectively. Let $\beta={A^{*}}^{T}u$ and $\beta^{\prime}={A^{*}}^{T}v$. Let $y=A^{*}x^{*}$ be random sample that is independent of $u$,$v$, then

where $q_{i}=\mathop{\bf Pr\/}[i\in S]$, $c_{i}=\mathop{\bf E\/}[{x^{*}_{i}}^{4}|i\neq S]$, and the error terms are:

Moreover the error terms $E_{1}+E_{2}+E_{3}$ has spectral norm bounded by $O^{*}(k/m\log n)$, $|\beta_{i}|\geq\Omega(1)$ for all $i\in\mbox{{\em supp}}(\alpha)$ and $|\beta^{\prime}_{i}|\geq\Omega(1)$ for all $i\in\mbox{{\em supp}}(\alpha^{\prime})$.

We will prove this lemma in three parts. First we compute the expectation and show it has the desired form. Recall that $y=A^{*}x^{*}$, and so:

where the second-to-last line follows because the entries in $x^{*}_{S}$ are independent and have mean zero, and the only non-zero terms come from ${x^{*}_{i}}^{4}$ (whose expectation is $c_{i}$) and ${x^{*}_{i}}^{2}{x^{*}_{j}}^{2}$ (whose expectation is one). Equation (4) now follows by rearranging the terms in the last line. What remains is to bound the spectral norm of $E_{1},E_{2}$ and $E_{3}$.

Next we establish some useful properties of $\beta$ and $\beta^{\prime}$:

With high probability it holds that $(a)$ for each $i$ we have $|\beta_{i}-\alpha_{i}|\leq\frac{\mu k\log m}{\sqrt{n}}$ and $(b)$ $\|\beta\|\leq O(\sqrt{mk/n})$.

In particular since the difference between $\beta_{i}$ an $\alpha_{i}$ is $o(1)$ for our setting of parameters, we conclude that if $\alpha_{i}\neq 0$ then $C-o(1)\leq|\beta_{i}|\leq O(\log m)$ and if $\alpha_{i}=0$ then $|\beta_{i}|\leq\frac{\mu k\log m}{\sqrt{n}}$.

Recall that $U$ is the support of $\alpha$ and let $R=U\backslash\{i\}$. Then:

and since $A^{*}$ is incoherent we have that $\|{A^{*}_{i}}^{T}A^{*}_{R}\|\leq\mu\sqrt{k/n}$. Moreover the entries in $\alpha_{R}$ are independent and subgaussian random variables, and so with high probability $|\langle{A^{*}_{i}}^{T}A^{*}_{R},\alpha_{R}\rangle|\leq\frac{\mu k\log m}{\sqrt{n}}$ and this implies the first part of the claim.

For the second part, we can bound $\|\beta\|\leq\|A^{*}\|\|A^{*}_{U}\|\|\alpha\|$. Since $\alpha$ is a $k$-sparse vector with independent and subgaussian entries, with high probability $\|\alpha\|\leq O(\sqrt{k})$ in which case it follows that $\|\beta\|\leq O(\sqrt{mk/n})$.

Now we are ready to bound the error terms.

With high probability each of the error terms $E_{1},E_{2}$ and $E_{3}$ in (4) has spectral norm bounded at most $O^{*}(k/m\log m)$.

Let $S=[m]\backslash(U\cap V)$, then $E_{1}=A^{*}_{S}D_{1}{A^{*}_{S}}^{T}$ where $D_{1}$ is a diagonal matrix whose entries are $q_{i}c_{i}\beta_{i}\beta^{\prime}_{i}$. We first bound $\|D_{1}\|$. To this end, we can invoke the first part of Claim 2 to conclude that $|\beta_{i}\beta^{\prime}_{i}|\leq\frac{\mu^{2}k^{2}\log^{2}m}{n}$. Also $q_{i}c_{i}=\Theta(k/m)$ and so

Finally $\|A_{S}\|\leq\|A\|\leq O(1)$ (where we have used the assumption that $m=O(n)$), and this yields the desired bound on $\|E_{1}\|$.

The second term $E_{2}$ is a sum of positive semidefinite matrices and we will make crucial use of this fact below:

Here the first inequality follows by using bounds on $q_{i,j}$ and then completing the square. The second inequality uses Cauchy-Schwartz. We can now invoke the second part of Claim 2 and conclude that $\|E_{2}\|\leq O(k^{2}/m^{2})\|\beta\|\|\beta^{\prime}\|\|A^{*}\|^{2}\leq O^{*}(k/m\log m)$ (where we have used the assumption that $m=O(n)$).

For the third error term $E_{3}$, by symmetry we need only consider terms of the form $q_{i,j}\beta_{i}\beta_{j}^{\prime}A^{*}_{i}{A^{*}_{j}}^{T}$. We can collect these terms and write them as $A^{*}Q{A^{*}}^{T}$, where $Q_{i,j}=0$ if $i=j$ and $Q_{i,j}=q_{i,j}\beta_{i}\beta_{j}^{\prime}$ if $i\neq j$. First, we bound the Frobenius norm of $Q$:

Finally $\|E_{3}\|\leq 2\|A^{*}\|^{2}\|Q\|\leq O(m/n\cdot k^{2}/m^{2})\|\beta\|\|\beta^{\prime}\|\leq O^{*}(k/m\log m)$, and this completes the proof of the claim.

The proof of the main lemma is now complete.

Conclusions

Going beyond $\sqrt{n}$ sparsity requires new ideas as alternating minimization appears to break down. Mysterious properties of alternating minimization are also left to explore, such as why a random initialization works. Are these heuristics information theoretically optimal in terms of their sample complexity? Finally, can we analyse energy minimization in other contexts as well?

References

Appendix A Threshold Decoding

Here we show that a simple thresholding method recovers the support of each sample with high probability (over the randomness of $x^{*}$). This corresponds to the fact that sparse recovery for incoherent dictionaries is much easier when the non-zero coefficients do not take on a wide range of values; in particular, one does not need iterative pursuit algorithms in this case. As usual let $y=A^{*}x^{*}$ be a sample from the model, and let $S$ be the support of $x^{*}$. Moreover suppose that $A^{*}$ is $\mu$-incoherent and let $A$ be column-wise $\delta$-close to $A$. Then

If $\frac{\mu}{\sqrt{n}}\leq\frac{1}{2k}$ and $k=\Omega^{*}(\log m)$ and $\delta=O^{*}(1/\sqrt{\log m})$, then with high probability (over the choice of $x^{*}$) we have $S=\{i~{}\mbox{:}~{}|\langle A_{i},y\rangle|>C/2\}$. Also for all $i\in S$ $\mbox{sgn}(\langle A_{i},y\rangle)=\mbox{sgn}(x^{*}_{i})$.