A note on the sample complexity of the Er-SpUD algorithm by Spielman, Wang and Wright for exact recovery of sparsely used dictionaries

Introduction

Learning sparsely-used dictionaries has recently attracted considerable attention in connection to applications in machine learning, signal processing or computational neuroscience. In particular, two important fields of applications are dictionary learning and blind source separation . We do not discuss these applications and refer the Reader to the aforesaid articles for details.

Among many approaches to this problem a particularly successful one has been presented by Spielman, Wang and Wright , who considered the noiseless-invertible case:

Consider an invertible $n\times n$ matrix $A$ and a random $n\times p$ sparse matrix $X$. Denote $Y=AX$. The objective is to reconstruct $A$ and $X$ (up to scaling and permutation of columns of $A$ and rows of $X$) based on the observable data $Y$.

The Authors of provide an algorithm which with high probability successfully recovers the matrices $A$ and $X$ up to rescaling and permutation of the columns of $A$ and rows of $X$, provided that $X$ is a sparse random matrix satisfying the following probabilistic assumptions.

$\chi_{ij},R_{ij}$ are independent random variables,

$R_{ij}$ are i.i.d., with mean zero and satisfy

Following we will say that matrices satisfying the above assumptions follow the Bernoulli-Subgaussian model with parameter $\theta$.

We remark that the constant 1/10 above is of no importance and has been chosen following and .

The approach of Spielman, Wang and Wright consists of two steps. At the first step (given by the Er-SpUD algorithm we describe below) one gathers $p/2$ candidates for the rows of $X$. The second, greedy step (Greedy algorithm, also described below) selects from the candidates the set of $n$ sparsest vectors, which form a matrix of rank $n$.

ER-SpUD(DC): Exact Recovery of Sparsely-Used Dictionaries using the sum of two columns of Y as constraint vectors.

Randomly pair the columns of $Y$ into $p/2$ groups $g_{j}=\{Ye_{j_{1}},Ye_{j_{2}}\}$.

For $j=1,\ldots,p/2$ Let $r_{j}=Ye_{j_{1}}+Ye_{j_{2}}$, where $g_{j}=\{Ye_{j_{1}},Ye_{j_{2}}\}$. Solve $\min_{w}\|w^{T}Y\|_{1}$ subject to $r_{j}^{T}w=1$, and set $s_{j}=w^{T}Y$.

Above we use the convention that if $r_{j}=0$ (which happens with nonzero probability), and as a consequence the minimization problem has no solution, then we skip the corresponding step of the algorithm.

The second stage, described below, is run on the set $S$ of vectors $s_{i}$ returned at the first stage (for notational simplicity we relabel them if $r_{j}=0$ for some $j$). We use the standard notation that $\|x\|_{0}$ denotes the number of nonzero coordinates of a vector $x$.

Greedy: A Greedy Algorithm to Reconstruct X and A.

For $i=1,\ldots,n$ REPEAT $l\leftarrow{\rm argmin\;}_{s_{l}\in S}\|s_{l}\|_{0}$, breaking ties arbitrarily $x_{i}=s_{l}$, $S=S\setminus\{s_{l}\}$ UNTIL $rank([x_{1},\ldots,x_{i}])=i$

Set $X=[x_{1},\ldots,x_{n}]^{T}$ and $A=YY^{T}(XY^{T})^{-1}$.

In it was proved that there exist positive constants $C,\alpha$, such that if

In particular it follows that the Greedy algorithm will extract from the set $\{s_{1},\ldots,s_{T}\}$ multiples of all $n$ rows of $X$ (note that all $s_{j}$’s are in the row space of $Y$ and thus also in the row space of $X$). Since, as one can prove, $X$ is with high probability of rank $n$, one easily proves that one can recover $A$ by the formula used in the 3rd step of the algorithm. We remark that in Luh and Vu obtained the same results concerning sparsity of linear combinations of rows of $X$ without the assumptions about the symmetry of the variables $R_{ij}$.

Note also that for $\theta$ of the order $n^{-1}$, $p=Cn\log n$ is necessary for uniqueness of the solution in the sense described above, otherwise with significant probability some of the rows of $X$ may be zero, which means that some columns of $A$ do not influence the matrix $Y$.

In it was also proved that if $p>Cn\log n$, $\theta>C^{\prime}\sqrt{\frac{\log n}{n}}$, then with high probability the ER-SpUD algorithm does not recover any of the rows of $X$.

Spielman, Wang and Wright have conjectured that their algorithm works with high probability provided that $p>Cn\log n$ (which, as mentioned above is required for well-posedness of the problem).

Recently, Luh and Vu have proved that the algorithm works for $p>Cn\log^{4}n$, which differs from the conjectured number of samples just by a polylogarithmic factor.

In this note we will consider a modified version of the algorithm with a slightly different first stage. Namely, instead of using only $p/2$ pairs of columns of $Y$, we will use all $\binom{p}{2}$ pairs. For fixed $p$ it clearly increases the time complexity of the algorithm (which however remains polynomial), but the advantage of this modification is the possibility of proving that it requires only $p=Cn\log n$ to recover $X$ and $A$ with high probability, which as explained above is optimal. More specifically, we will consider the following algorithm.

Modified ER-SpUD(DC): Exact Recovery of Sparsely-Used Dictionaries using the sum of two columns of Y as constraint vectors.

For $i=1,\ldots,p-1$ aa For $j=i+1,\ldots,p$ aaaaLet $r_{ij}=Ye_{i}+Ye_{j}$ aaaaSolve $\min_{w}\|w^{T}Y\|_{1}$ subject to $r_{ij}^{T}w=1$, and set $s_{ij}=w^{T}Y$.

There exist absolute constants $C,\alpha\in(0,\infty)$ such that if

and $X$ follows the Bernoulli-Subgaussian model with parameter $\theta$, then for $p\geq Cn\log n$, with probability at least $1-1/p$ the modified ER-SpUD algorithm successfully recovers all the rows of $X$, i.e. multiples of all the rows of $X$ are present among the vectors $s_{ij}$ returned by the algorithm.

Very recently in , Sun, Qing and Wright proposed an algorithm with polynomial sample complexity, which recovers well conditioned dictionaries under the assumption that the variables $R_{ij}$ are i.i.d. standard Gaussian and $\theta\leq 1/2$, thus allowing for the first time for a linear number of nonzero entries per column of the matrix $X$. Their novel approach is based on non-convex optimization. The sample complexity of the algorithms in is however higher then for the Er-SpUD algorithm; as mentioned by the Authors, numerical simulations suggest that it is at least $p=\Omega(n^{2}\log n)$ even in the case of orthogonal matrix $A$. The Authors of conjecture that algorithms with sample complexity $p=O(n\log n)$ should be possible also for large $\theta$.

Proof of Theorem 1.1

We will follow the general approach presented in and . The main new part of the argument is an improved bound on the sample complexity for empirical approximation of first moments of arbitrary marginals of the columns of the matrix $X$, given in Proposition 2.1 below. So as not to reproduce technical and lengthy parts of the original proof, we organize this section as follows. First, we present the crucial Proposition 2.1 and provide a brief discussion of its mathematical content. Next, we present an overview of the main steps in the proof scheme of . For parts of the proof not related to Proposition 2.1 or to the modification of the algorithm considered here, we only indicate the relevant statements from , while for the part involving the use of Proposition 2.1 and for the conclusion of the proof we provide the full argument. Proposition 2.1 is proved in Section 3.

Define the random vectors $Z_{1},\ldots,Z_{p}$ with the equality $Z_{i}(j)=U_{i}(j)\chi_{i}(j)$ for $1\leq i\leq p$, $1\leq j\leq n$ and consider the random variable

Then, for some universal constant $C$ and every $q\geq\max(2,\log n)$,

Let us also remark that in the above proposition we do not require independence between components of the random vectors $U_{i}$ or $\chi_{i}$ for fixed $i$, but just independence between the random vectors $U_{i},\chi_{i},i=1,\ldots,p$.

As announced, we will now present an outline of the proof of Theorem 1.1, indicating which steps differ from the original argument in .

Recall that $r_{ij}$ are sums of two columns of the matrix $Y$. At the first step of the proof, instead of looking at the original optimization problem

one performs a change of variables $z=A^{T}w$, $b_{ij}=A^{-1}r_{ij}$, arriving at the optimization problem

Note that one cannot solve (6) since it involves the unknown matrices $X$ and $A$. The goal of the subsequent steps is to prove that with probability separated from zero the solution $z_{\ast}$ of (6) is a multiple of one of the basis vectors $e_{1},\ldots,e_{n}$, say $z_{\ast}=\lambda e_{k}$. This means that $w_{\ast}^{T}Y=z_{\ast}^{T}X=\lambda e_{k}^{T}X$, i.e. (5) recovers the $k$-th row of $X$ up to scaling.

Step 2. The solution $z_{\ast}$ satisfies $\operatorname{{\rm supp}}(z_{\ast})\subseteq\operatorname{{\rm supp}}(b_{ij})$.

At this step we prove the following lemma, which is a counterpart of Lemma 11 in . It is weaker in that we do not consider arbitrary vectors $b_{ij}$, but only sums of two distinct columns of $X$ (which is enough for the application in the proof of Theorem 1.1). On the other hand it works already for $p>Cn\log n$ and not for $p>Cn^{2}\log n$ as the original lemma from .

For $1\leq i<j\leq p$, define $b_{ij}=Xe_{i}+Xe_{j}$, $I_{ij}=(\operatorname{{\rm supp}}Xe_{i})\cup(\operatorname{{\rm supp}}Xe_{j})$. There exist numerical constants $C,\alpha>0$ such that if $2/n\leq\theta\leq\alpha/\sqrt{n}$ and $p>Cn\log n$, then with probability at least $1-p^{-2}$ the random matrix $X$ has the following property:

(P1) For every $1\leq i<j\leq p$ either $|I_{ij}|\in\{0\}\cup(1/(8\theta),n]$ or every solution $z_{\ast}$ to the optimization problem (6) satisfies $\operatorname{{\rm supp}}z_{\ast}\subseteq I_{ij}$.

To prove the above lemma, one first shows a counterpart of Lemma 16 in .

Let $\varepsilon_{1},\ldots,\varepsilon_{n}$ be a sequence of i.i.d. Rademacher variables, independent of $Z$. By standard symmetrization inequalities (see e.g. Lemma 6.3. in ),

The next lemma is an improvement of Lemma 17 in , which is crucial for obtaining Lemma 2.2.

Note first that by increasing the set $S$, we increase $\|v^{T}X_{S}\|_{1}$, so without loss of generality we can assume that $|S|=\lfloor p/4\rfloor$. Apply Proposition 2.1 with the vectors $U_{j}=(R_{1j},\ldots,R_{nj})$ and $\chi_{j}=(\chi_{1j},\ldots,\chi_{nj})$ and $q=8\log p$. Note that our integrability assumptions on $R_{ij}$ imply (1) with $M$ being a universal constant. Therefore, for some absolute constant $C$ and $p\geq Cn\log n$, with probability at least $1-p^{-8}$ we have

where we used that for $C$ sufficiently large, $p/\log p\geq n\geq 1/\theta$.

In particular this means that (using the notation of Proposition 2.1)

Now, by Lemma 2.3 and the assumed bound on the cardinality of $S$, we get

for $p>C^{\prime}n\log n$, where $C^{\prime}$ is another absolute constant. ∎

We are now in position to prove Lemma 2.2.

We will show that for each $1\leq i<j\leq p$ the probability that $0<|I_{ij}|\leq 1/(8\theta)$ and there exists a solution to (6) not supported on $I_{ij}$ is bounded from above by $1/p^{4}$. By the union bound over all $i<j$, this implies the lemma.

Fix $i,j$ and let $S=\{l\in[p]\colon\exists_{k\in I_{ij}}X_{kl}\neq 0\}$. Denote by $\mathcal{F}_{1}$ the $\sigma$-field generated by $Xe_{i}$ and $Xe_{j}$. Then $\mathcal{A}=\{0<|I_{ij}|\leq 1/(8\theta)\}\in\mathcal{F}_{1}$. By independence, for each $k\notin\{i,j\}$, on the event $\mathcal{A}$,

where the second inequality holds if $\alpha$ is sufficiently small.

Thus, by independence of columns of $X$ and Hoeffding’s inequality,

where we used the fact that $z_{1}^{T}Xe_{i}=z_{1}^{T}Xe_{j}=0$.

Denote by $\mathcal{F}_{2}$ the $\sigma$-field generated by $Xe_{i},Xe_{j}$ and the rows of $X$ labeled by $I_{ij}$ (note that $I_{ij}$ is itself random, but this will not be a problem in what follows). The random set $S$ is measurable with respect to $\mathcal{F}_{2}$. Moreover, due to independence and identical distribution of the entries of $X$, conditionally on $\mathcal{F}_{2}$ the matrix $X^{\prime}$ still follows the Bernoulli-Subgaussian model with parameter $\theta$. Therefore, by Lemma 2.4, if $C$ is large enough, then on $\{|S^{\prime}|\leq p/4\}$ we have

Note that by the definition of $z_{0}$, we have $b_{ij}^{T}z_{0}=b_{ij}^{T}z=1$, therefore $z_{0}$ is a feasible candidate for the solution of the optimization problem (6). Thus, we have $\|z_{1}^{T}X^{\prime}\|_{1}-2\|z_{1}^{T}X^{\prime}_{S^{\prime}}\|_{1}\leq 0$ and as a consequence, on the event $\{|S^{\prime}|\leq p/4\}$,

Thus, denoting \mathcal{B}=\{\textrm{for some solutionz_{\ast}to \eqref{eq:changed-problem},}\;z_{1}\neq 0\;\textrm{and}\;0<|I_{ij}^{c}|<1/(8\theta)\}, we get by (7) and (8),

for $p>Cn\log n$ with a sufficiently large absolute constant $C$. ∎

Step 3. With high probability $z_{\ast}=\lambda e_{k}$ for $k=\operatorname{{\rm argmax}}_{1\leq l\leq n}|b_{ij}(l)|$.

At this step one proves the following lemma (Lemma 12 in ). Since no changes with respect to the original argument are required (we do not use Proposition 2.1 here), we do not reproduce the proof and refer the Reader to for details. We remark that although the lemma is formulated in for symmetric variables, the symmetry assumption is not used in its proof.

Below, by $|b|^{\downarrow}_{1}\geq|b|^{\downarrow}_{2}\geq\ldots\geq|b|^{\downarrow}_{n}$, we denote the nonincreasing rearrangement of the sequence $|b_{1}|,\ldots,|b_{n}|$, while for $J\subseteq[n]$, $X^{J}$ denotes the matrix obtained from $X$ by selecting the rows indexed by the set $J$.

with probability at least $1-4p^{-10}$, the random matrix $X$ has the following property.

is unique, 1-sparse, and is supported on the index of the largest entry of $b$.

Set $s=12\theta n+1$. Our first goal is to prove that with probability at least $1-1/p^{2}$, for all $k\in[n]$, there exist $i,j\in[p]$, $i\neq j$ such that the vector $b=Xe_{i}+Xe_{j}$ satisfies the assumptions of Lemma 2.5, $|b|^{\downarrow}_{1}=|b_{k}|$ and $I_{ij}:=(\operatorname{{\rm supp}}Xe_{i})\cup(\operatorname{{\rm supp}}Xe_{j})$ satisfies $0<|I_{ij}|\leq 1/(8\theta)$, which will allow us to take advantage of Lemma 2.2.

We will assume that $p\geq 2Cn\log n$ for some numerical constant $C$ to be fixed later on. For $k\in[n]$, consider the events

We will first show that for all $k\in[n]$,

Let us start with the proof of (10). Set $\mathcal{B}_{ki}=\{|\{r\in[n]\setminus\{k\}\colon\chi_{rk}=1\}|\leq(s-1)/2\}$. By independence we have

where we used the inequality $p/\log p\geq 16c_{1}^{-1}n$ for $p\geq Cn\log n$. We have thus established (10).

Let us now pass to (11). Denote by $\mathcal{F}_{1}$ the $\sigma$-field generated by $\chi_{ki},R_{ki},k\in[n],1\leq i\leq\lfloor p/2\rfloor$.

For $\omega\in\mathcal{A}_{k}$ define $i_{\min}(\omega)=\min\{1\leq i\leq\lfloor p/2\rfloor\colon\omega\in\mathcal{E}_{ki}\}$. Note that on $\mathcal{A}_{k}$,

Similarly as in the argument leading to (10), for fixed $j$, using the independence of the variables $\chi_{lm},R_{lm}$ we obtain

Now recall that $\theta\leq\frac{\alpha}{n}$ for some universal constant $\alpha$. If $\alpha$ is small enough then $1-\theta\geq e^{-2\theta}$ and

Since $\frac{2\theta(n-1)}{s-1}\leq\frac{1}{6}$, this implies that

for some positive universal constant $c_{2}$. Since the events \mathcal{E}_{kj}\cap\Big{\{}\{l\in[n]\colon\chi_{lj_{\min}}=\chi_{lk}=1\}=\{k\}\Big{\}}, $\lfloor p/2\rfloor<k\leq p$ are conditionally independent, given $\mathcal{F}_{1}$, we obtain that on $\mathcal{A}_{k}$,

provided $C$ is a sufficiently large universal constant. Now, using (10), we get

Taking the union bound over $k\in[n]$, we get

the largest entry of $b$ (in absolute value) equals $b_{k}\geq 2q>0$ whereas the remaining entries do not exceed $q$,

In particular, by property P1 we obtain that any solution $z_{\ast}$ to the problem (6) satisfies $\operatorname{{\rm supp}}z_{\ast}\subseteq I_{ij}$. Therefore for some (any) $J\supseteq I_{ij}$ with $|J|=s$, we obtain (identifying vectors supported on $J$ with their restrictions to $J$), that $z_{\ast}$ is in fact a solution to the restricted problem (9) with $b=b_{ij}$, which by property P2 implies that $z_{\ast}=\lambda e_{k}$ for some $\lambda\neq 0$.

According to the discussion at the beginning of Step 1, this means that the solution $w_{\ast}$ to (5) satisfies $w_{\ast}^{T}Y=\lambda e_{k}^{T}X$, i.e. the algorithm, when analyzing the vector $b_{ij}$, will add a multiple of the $k$-th row of $X$ to the collection $S$.

Proof of Proposition 2.1

The first tool we will need is the classical Bernstein’s inequality (see e.g. Lemma 2.2.11 in ).

Another (also quite standard) tool we will rely on is the contraction principle for empirical processes due to Talagrand (see Theorem 4.12. in ).

Let $\varepsilon_{1},\ldots,\varepsilon_{p}$ be i.i.d. Rademacher variables, independent of the sequences $(U_{i})$, $(\chi_{i})$. By the symmetrization inequality (see e.g. Lemma 6.3. in ) we have

Now, since the function $t\mapsto|t|$ is a contraction, an application of Lemma 3.2, conditionally on $Z_{i}$, gives

Now, for every $i,j$ and every integer $k\geq 2$ we have

with $v=4\theta M^{2}$. Thus by the moment version (12) of Bernstein’s inequality for some universal constant $C$ we get

which, when combined with (3), yields for $q\geq\log n$,

The first part of the proposition follows by adjusting the constant $C$. The tail bound is a direct consequence of the Chebyshev inequality for the $q$-th moment. ∎

References