Complete Dictionary Recovery over the Sphere I: Overview and the Geometric Picture

I Introduction

Traditionally, concise signal representations have relied heavily on explicit analytic bases constructed in nonlinear approximation and harmonic analysis. This constructive approach has proved highly successful; the numerous theoretical advances in these fields (see, e.g., for summary of relevant results) provide ever more powerful representations, ranging from the classic Fourier basis to modern multidimensional, multidirectional, multiresolution bases, including wavelets, curvelets, ridgelets, and so on. However, two challenges confront practitioners in adapting these results to new domains: which function class best describes signals at hand, and consequently which representation is most appropriate. These challenges are coupled, as function classes with known “good” analytic bases are rare. As Donoho et al put it, “…in effect, uncovering the optimal codebook structure of naturally occurring data involves more challenging empirical questions than any that have ever been solved in empirical work in the mathematical sciences.”

Around 1996, neuroscientists Olshausen and Field discovered that sparse coding, the principle of encoding a signal with few atoms from a learned dictionary, reproduces important properties of the receptive fields of the simple cells that perform early visual processing . The discovery has spurred a flurry of algorithmic developments and successful applications for DL in the past two decades, spanning classical image processing, visual recognition, compressive signal acquisition, and also recent deep architectures for signal classification (see, e.g., for review of this development).

The learning approach is particularly relevant to modern signal processing and machine learning, which deal with data of huge volume and great variety (e.g., images, audios, graphs, texts, genome sequences, time series, etc). The proliferation of problems and data seems to preclude analytically deriving optimal representations for each new class of data in a timely manner. On the other hand, as datasets grow, learning dictionaries directly from data looks increasingly attractive and promising. When armed with sufficiently many data samples of one signal class, by solving the model DL problem, one would expect to obtain a dictionary that allows sparse representation for the whole class. This hope has been borne out in a number of successful examples and theories .

In contrast to the above empirical successes, theoretical study of dictionary learning is still developing. For applications in which dictionary learning is to be applied in a “hands-free” manner, it is desirable to have efficient algorithms which are guaranteed to perform correctly, when the input data admit a sparse model. There have been several important recent results in this direction, which we will review in Section I-E, after our sketching main results. Nevertheless, obtaining algorithms that provably succeed under broad and realistic conditions remains an important research challenge.

To understand where the difficulties arise, we can consider a model formulation, in which we attempt to obtain the dictionary $\mathbf{A}$ and coefficients $\mathbf{X}$ which best trade-off sparsity and fidelity to the observed data:

Here, $\left\|\mathbf{X}\right\|_{1}\doteq\sum_{i,j}\left|X_{ij}\right|$ promotes sparsity of the coefficients, $\lambda\geq 0$ trades off the level of coefficient sparsity and quality of approximation, and $\mathcal{A}$ imposes desired structures on the dictionary.

This formulation is nonconvex: the admissible set $\mathcal{A}$ is typically nonconvex (e.g., orthogonal group, matrices with normalized columns)For example, in nonlinear approximation and harmonic analysis, orthonormal basis or (tight-)frames are preferred; to fix the scale ambiguity discussed in the text, a common practice is to require that $\mathbf{A}$ to be column-normalized. , while the most daunting nonconvexity comes from the bilinear mapping: $\left(\mathbf{A},\mathbf{X}\right)\mapsto\mathbf{A}\mathbf{X}$. Because $\left(\mathbf{A},\mathbf{X}\right)$ and $\left(\mathbf{A}\mathbf{\Pi}\mathbf{\Sigma},\mathbf{\Sigma}^{-1}\mathbf{\Pi}^{*}\mathbf{X}\right)$ result in the same objective value for the conceptual formulation (I.1), where $\mathbf{\Pi}$ is any permutation matrix, and $\mathbf{\Sigma}$ any diagonal matrix with diagonal entries in $\{\pm 1\}$, and $\left(\cdot\right)^{*}$ denotes matrix transpose. Thus, we should expect the problem to have combinatorially many global minimizers. These global minimizers are generally isolated, likely jeopardizing natural convex relaxation (see similar discussions in, e.g., and ).Simple convex relaxations normally replace the objective function with a convex surrogate, and the constraint set with its convex hull. When there are multiple isolated global minimizers for the original nonconvex problem, any point in the convex hull of these global minimizers are necessarily feasible for the relaxed version, and such points tend to produce smaller or equal values than that of the original global minimizers by the relaxed objective function, due to convexity. This implies such relaxations are bound to be loose. Semidefinite programming (SDP) lifting may be one useful general strategy to convexify bilinear inverse problems, see, e.g., . However, for problems with general nonlinear constraints, it is unclear whether the lifting always yields tight relaxation; consider, e.g., and the identification issue in blind deconvolution . This contrasts sharply with problems in sparse recovery and compressed sensing, in which simple convex relaxations are often provably effective . Is there any hope to obtain global solutions to the DL problem?

I-B An Intriguing Numerical Experiment with Real Images

We provide empirical evidence in support of a positive answer to the above question. Specifically, we learn orthogonal bases (orthobases) for real images patches. Orthobases are of interest because typical hand-designed dictionaries such as discrete cosine (DCT) and wavelet bases are orthogonal, and orthobases seem competitive in performance for applications such as image denoising, as compared to overcomplete dictionaries See Section I-C for more detailed discussions of this point. also gave motivations and algorithms for learning (union of) orthobases as dictionaries. .

We divide a given greyscale image into $8\times 8$ non-overlapping patches, which are converted into $64$-dimensional vectors and stacked column-wise into a data matrix $\mathbf{Y}$. Specializing (I.1) to this setting, we obtain the optimization problem:

where $O_{n}$ is the set of order $n$ orthogonal matrices, i.e., order-$n$ orthogonal group. To derive a concrete algorithm for (I.2), one can deploy the alternating direction method (ADM)This method is also called alternating minimization or (block) coordinate descent method. see, e.g., for classic results and for several interesting recent developments. , i.e., alternately minimizing the objective function with respect to (w.r.t.) one variable while fixing the other. The iteration sequence actually takes very simple form: for $k=1,2,3,\dots$,

where $\mathcal{S}_{\lambda}\left[\cdot\right]$ denotes the well-known soft-thresholding operator acting elementwise on matrices, i.e., $\mathcal{S}_{\lambda}\left[x\right]\doteq\operatorname{sign}\left(x\right)\max\left(\left|x\right|-\lambda,0\right)$ for any scalar $x$.

Fig. 1 shows what we obtained using the simple ADM algorithm, with independent and randomized initializations:

The algorithm seems to always produce the same optimal value, regardless of the initialization.

This observation is consistent with the possibility that the heuristic ADM algorithm may always converge to a global minimizer! Technically, the convergence to global solutions is surprising because even convergence of ADM to critical points is not guaranteed in general, see, e.g., and references therein. Equally surprising is that the phenomenon has been observed on real imagesActually the same phenomenon is also observed for simulated data when the coefficient matrix obeys the Bernoulli-Gaussian model, which is defined later. The result on real images supports that previously claimed empirical successes over two decades may be non-incidental. . One may imagine only random data typically have “favorable” structures; in fact, almost all existing theories for DL pertain only to random data .

I-C Dictionary Recovery and Our Results

To define a reasonably simple and structured problem, we make the following assumptions:

The target dictionary $\mathbf{A}_{0}$ is complete, i.e., square and invertible ($m=n$). In particular, this class includes orthogonal dictionaries. Admittedly overcomplete dictionaries tend to be more powerful for modeling and to allow sparser representations. Nevertheless, most classic hand-designed dictionaries in common use are orthogonal. Orthobases are competitive in performance for certain tasks such as image denoising , and admit faster algorithms for learning and encoding. Empirically, there is no systematic evidence supporting that overcomplete dictionaries are strictly necessary for good performance in all published applications (though argues for the necessity from a neuroscience perspective). Some of the ideas and tools developed here for complete dictionaries may also apply to certain classes of structured overcomplete dictionaries, such as tight frames. See Section III for relevant discussion.

In this paper, we provide a nonconvex formulation for the DR problem, and characterize the geometric structure of the formulation that allows development of efficient algorithms for optimization. In the companion paper , we derive an efficient algorithm taking advantage of the structure, and describe a complete algorithmic pipeline for efficient recovery. Together, we prove the following result:

Obviously, even if $\mathbf{X}_{0}$ is known, one needs $p\geq n$ to make the identification problem well posed. Under our particular probabilistic model, a simple coupon collection argument implies that one needs $p\geq\Omega\left(\tfrac{1}{\theta}\log n\right)$ to ensure all atoms in $\mathbf{A}_{0}$ are observed with high probability (w.h.p.). Ensuring that an efficient algorithm exists may demand more. Our result implies when $p$ is polynomial in $n$, $1/\theta$ and $\kappa(\mathbf{A}_{0})$, recovery with an efficient algorithm is possible.

The parameter $\theta$ controls the sparsity level of $\mathbf{X}_{0}$. Intuitively, the recovery problem is easy for small $\theta$ and becomes harder for large $\theta$.Indeed, when $\theta$ is small enough such that columns of $\mathbf{X}_{0}$ are predominately $1$-sparse, one directly observes scaled versions of the atoms (i.e., columns of $\mathbf{X}_{0}$); when $\mathbf{X}_{0}$ is fully dense corresponding to $\theta=1$, recovery is never possible as one can easily find another complete $\mathbf{A}_{0}^{\prime}$ and fully dense $\mathbf{X}_{0}^{\prime}$ such that $\mathbf{Y}=\mathbf{A}_{0}^{\prime}\mathbf{X}_{0}^{\prime}$ with $\mathbf{A}_{0}^{\prime}$ not equivalent to $\mathbf{A}_{0}$. It is perhaps surprising that an efficient algorithm can succeed up to constant $\theta$, i.e., linear sparsity in $\mathbf{X}_{0}$. Compared to the case when $\mathbf{A}_{0}$ is known, there is only at most a constant gap in the sparsity level one can deal with.

For DL, our result gives the first efficient algorithm that provably recovers complete $\mathbf{A}_{0}$ and $\mathbf{X}_{0}$ when $\mathbf{X}_{0}$ has $O(n)$ nonzeros per column under appropriate probability model. Section I-E provides detailed comparison of our result with other recent recovery results for complete and overcomplete dictionaries.

I-D Main Ingredients and Innovations

In this section we describe three main ingredients that we use to obtain the stated result.

The constraint means at least one coordinate of $\mathbf{q}^{*}\mathbf{Y}$ has unit magnitudeThe sign ambiguity is tolerable here. . Thus, (I.4) reduces to a sequence of convex (linear) programs. has shown that (see also ) solving (I.4) recovers $\left(\mathbf{A}_{0},\mathbf{X}_{0}\right)$ for very sparse $\mathbf{X}_{0}$, but the idea provably breaks down when $\theta$ is slightly above $O(1/\sqrt{n})$, or equivalently when each column of $\mathbf{X}_{0}$ has more than $O\left(\sqrt{n}\right)$ nonzeros.

Inspired by our previous image experiment, we work with a nonconvex alternativeA similar formulation has been proposed in in the context of blind source separation; see also . :

which is infinitely differentiable and $\mu$ controls the smoothing level.In fact, there is nothing special about this choice and we believe that any valid smooth (twice continuously differentiable) approximation to $\left|\cdot\right|$ would work and yield qualitatively similar results. We also have some preliminary results showing the latter geometric picture remains the same for certain nonsmooth functions, such as a modified version of the Huber function, though the analysis involves handling a different set of technical subtleties. The algorithm also needs additional modifications.

I-D2 A Glimpse into High-dimensional Function Landscape

Two challenges stand out when implementing this idea. For geometry, one has to show similar structure exists for general complete $\mathbf{A}_{0}$, in high dimensions ($n\geq 3$), when the number of observations $p$ is finite (vs. the expectation in the experiment). For algorithms, we need to be able to take advantage of this structure without knowing $\mathbf{A}_{0}$ ahead of time. In Section I-D3, we describe a Riemannian trust region method which addresses the latter challenge.

To study the function on this exemplar region, we again invoke the projection trick described above, this time onto the equatorial section $\mathbf{e}_{n}^{\perp}$. This can be formally captured by the reparameterization mapping:

It can be verified the exemplar we chose to work with is strictly contained in this setIndeed, if $\left\langle\mathbf{q},\mathbf{e}_{n}\right\rangle\geq\left|\left\langle\mathbf{q},\mathbf{e}_{i}\right\rangle\right|$ for all $i\neq n$, $1-\left\|\mathbf{w}\right\|^{2}=q_{n}^{2}\geq 1/n$, implying $\left\|\mathbf{w}\right\|^{2}\leq\tfrac{n-1}{n}<\tfrac{4n-1}{4n}$. The reason we have defined an open set instead of a closed (compact) one is to avoid potential trivial local minimizers located on the boundary. We study behavior of $g$ over this slightly larger set $\Gamma$, instead of just the projection of the chosen symmetric section, to conveniently deal with the boundary effect: if we choose to work with just projection of the chosen symmetric section, there would be considerable technical subtleties at the boundaries when we call the union argument to cover the whole sphere. . This is illustrated for the case $n=3$ in Fig. 4 (right).

Our analysis characterizes the properties of $g\left(\mathbf{w};\mathbf{X}_{0}\right)$ by studying three quantities

respectively over three consecutive regions moving away from the origin, corresponding to the three regions in Fig. 3 (right). In particular, through typical expectation-concentration style arguments, we show that there exists a positive constant $c$ such that

over the respective regions w.h.p., confirming our low-dimensional observations described above. In particular, the favorable structure we observed for $n=3$ persists in high dimensions, w.h.p., even when $p$ is large yet finite, for the case $\mathbf{A}_{0}$ is orthogonal. Moreover, the local minimizer of $g\left(\mathbf{w};\mathbf{X}_{0}\right)$ over $\Gamma$ is very close to $\mathbf{0}$, within a distance of $O\left(\mu\right)$When $p\to\infty$, the local minimizer is exactly $\mathbf{0}$; deviation from $\mathbf{0}$ that we described is due to finite-sample perturbation. The deviation distance depends both the $h_{\mu}(\cdot)$ and $p$; see Theorem II.1 for example. .

For general complete dictionaries $\mathbf{A}_{0}$, we hope that the function $f$ retains the nice geometric structure discussed above. We can ensure this by “preconditioning” $\mathbf{Y}$ such that the output looks as if being generated from a certain orthogonal matrix, possibly plus a small perturbation. We can then argue that the perturbation does not significantly affect qualitative properties of the objective landscape. Write

where $\mathtt{SVD}(\mathbf{A}_{0})=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{*}$. It is easy to see $\mathbf{U}\mathbf{V}^{*}$ is an orthogonal matrix. Hence the preconditioning scheme we have introduced is technically sound.

Our analysis shows that $\overline{\mathbf{Y}}$ can be written as

where $\mathbf{\Xi}$ is a matrix with a small magnitude. Simple perturbation argument shows that the constant $c$ in (I.9) is at most shrunk to $c/2$ for all $\mathbf{w}$ when $p$ is sufficiently large. Thus, the qualitative aspects of the geometry have not been changed by the perturbation.

I-D3 A Second-order Algorithm on Manifold: Riemannian Trust-Region Method

We do not know $\mathbf{A}_{0}$ ahead of time, so our algorithm needs to take advantage of the structure described above without knowledge of $\mathbf{A}_{0}$. Intuitively, this seems possible as the descent direction in the $\mathbf{w}$ space appears to also be a local descent direction for $f$ over the sphere. Another issue is that although the optimization problem has no spurious local minimizers, it does have many saddle points with indefinite Hessian, which we call ridable saddles See and . (Fig. 3). We can use second-order information to guarantee to escape from such saddle points. In the companion paper , we derive an algorithm based on the Riemannian trust region method (TRM) for this purpose. There are other algorithmic possibilities; see, e.g., .

We provide here only the basic intuition why a local minimizer can be retrieved by the second-order trust-region method. Consider an unconstrained optimization problem

Typical (second-order) TRM proceeds by successively forming a second-order approximation to $\phi$ at the current iterate,

where $\mathbf{Q}(\mathbf{x}^{(r-1)})$ is a proxy for the Hessian matrix $\nabla^{2}\phi(\mathbf{x}^{(r-1)})$, which encodes the second-order geometry. The next movement direction is determined by seeking a minimum of $\widehat{\phi}(\mathbf{\delta};\mathbf{x}^{(r-1)})$ over a small region, normally a norm ball $\|\mathbf{\delta}\|_{p}\leq\Delta$, called the trust region, inducing the well-studied trust-region subproblem that can efficiently solved:

where $\Delta$ is called the trust-region radius that controls how far the movement can be made. If we take $\mathbf{Q}(\mathbf{x}^{(r-1)})=\nabla^{2}\phi(\mathbf{x}^{(r-1)})$ for all $r$, then whenever the gradient is nonvanishing or the Hessian is indefinite, we expect to decrease the objective function by a concrete amount provided $\|\mathbf{\delta}\|$ is sufficiently small. Since the domain is compact, the iterate sequence ultimately moves into the strongly convex region, where the trust-region algorithm behaves like a typical Newton algorithm. All these are generalized to our objective over the sphere and made rigorous in the companion paper .

I-E Prior Arts and Connections

It is far too ambitious to include here a comprehensive review of the exciting developments of DL algorithms and applications after the pioneer work . We refer the reader to Chapter 12 - 15 of the book and the survey paper for summaries of relevant developments in image analysis and visual recognition. In the following, we focus on reviewing recent developments on the theoretical side of dictionary learning, and draw connections to problems and techniques that are relevant to the current work.

The theoretical study of DL in the recovery setting started only very recently. was the first to provide an algorithmic procedure to correctly extract the generating dictionary. The algorithm requires exponentially many samples and has exponential running time; see also . Subsequent work studied when the target dictionary is a local optimizer of natural recovery criteria. These meticulous analyses show that polynomially many samples are sufficient to ensure local correctness under natural assumptions. However, these results do not imply that one can design efficient algorithms to obtain the desired local optimizer and hence the dictionary.

initiated the on-going research effort to provide efficient algorithms that globally solve DR. They showed that one can recover a complete dictionary $\mathbf{A}_{0}$ from $\mathbf{Y}=\mathbf{A}_{0}\mathbf{X}_{0}$ by solving a certain sequence of linear programs, when $\mathbf{X}_{0}$ is a sparse random matrix (under the Bernoulli-Subgaussian model) with $O(\sqrt{n})$ nonzeros per column (and the method provably breaks down when $\mathbf{X}_{0}$ contains slightly more than $\Omega(\sqrt{n})$ nonzeros per column). and gave efficient algorithms that provably recover overcomplete ($m\geq n$), incoherent dictionaries, based on a combination of {clustering or spectral initialization} and local refinement. These algorithms again succeed when $\mathbf{X}_{0}$ has $\widetilde{O}(\sqrt{n})$ The $\widetilde{O}$ suppresses some logarithm factors. nonzeros per column. Recent work provided the first polynomial-time algorithm that provably recovers most “nice” overcomplete dictionaries when $\mathbf{X}_{0}$ has $O(n^{1-\delta})$ nonzeros per column for any constant $\delta\in(0,1)$. However, the proposed algorithm runs in super-polynomial (quasipolynomial) time when the sparsity level goes up to $O(n)$. Similarly, also proposed a super-polynomial time algorithm that guarantees recovery with (almost) $O\left(n\right)$ nonzeros per column. Detailed models for those methods dealing with overcomplete dictionaries are differ from one another; nevertheless, they all assume each column of $\mathbf{X}_{0}$ has bounded sparsity levels, and the nonzero coefficients have certain sub-Gaussian magnitudesThus, one may anticipate that the performances of those methods do not change much qualitatively, if the BG model for the coefficients had been assumed. . By comparison, we give the first polynomial-time algorithm that provably recovers complete dictionary $\mathbf{A}_{0}$ when $\mathbf{X}_{0}$ has $O\left(n\right)$ nonzeros per column, under the BG model.

Aside from efficient recovery, other theoretical work on DL includes results on identifiability , generalization bounds , and noise stability .

We have followed and cast the core problem as finding the sparsest vectors in a given linear subspace, which is also of independent interest. Under a planted sparse model… where one sparse vector embedded in an otherwise random subspace., showed that solving a sequence of linear programs similar to can recover sparse vectors with sparsity up to $O\left(p/\sqrt{n}\right)$, sublinear in the vector dimension. improved the recovery limit to $O\left(p\right)$ by solving a nonconvex sphere-constrained problem similar to (I.5)The only difference is that they chose to work with the Huber function as a proxy of the $\left\|\cdot\right\|_{1}$ function. via an ADM algorithm. The idea of seeking rows of $\mathbf{X}_{0}$ sequentially by solving the above core problem sees precursors in for blind source separation, and for matrix sparsification. also proposed a nonconvex optimization similar to (I.5) here and that employed in .

For other nonconvex optimization problems of recovery of structured signalsThis is a body of recent work studying nonconvex recovery up to statistical precision, including, e.g., . , including low-rank matrix completion/recovery , phase retreival , tensor recovery , mixed regression , structured element pursuit , and recovery of simultaneously structured signals , numerical linear algebra and optimization , the initialization plus local refinement strategy adopted in theoretical DL is also crucial: nearness to the target solution enables exploiting the local property of the optimizing objective to ensure that the local refinement succeeds.The powerful framework to establish local convergence of ADM algorithms to critical points applies to DL/DR also, see, e.g., . However, these results do not guarantee to produce global optima. By comparison, we provide a complete characterization of the global geometry, which admits efficient algorithms without any special initialization.

DL can also be considered in the general framework of matrix factorization problems, which encompass the classic principal component analysis (PCA), ICA, and clustering, and more recent problems such as nonnegative matrix factorization (NMF), multi-layer neural nets (deep learning architectures). Most of these problems are NP-hard. Identifying tractable cases of practical interest and providing provable efficient algorithms are subject of on-going research endeavors; see, e.g., recent progresses on NMF , and learning deep neural nets .

ICA factors a data matrix $\mathbf{Y}$ as $\mathbf{Y}=\mathbf{A}\mathbf{X}$ such that $\mathbf{A}$ is square and rows of $\mathbf{X}$ achieve maximal statistical independence . In theoretical study of the recovery problem, it is often assumed that rows of $\mathbf{X}_{0}$ are (weakly) independent (see, e.g., ). Our i.i.d. probability model on $\mathbf{X}_{0}$ implies rows of $\mathbf{X}_{0}$ are independent, aligning our problem perfectly with the ICA problem. More interestingly, the $\log\cosh$ objective we analyze here was proposed as a general-purpose contrast function in ICA that has not been thoroughly analyzed . Algorithms and analysis with another popular contrast function, the fourth-order cumulants, however, indeed overlap with ours considerably Nevertheless, the objective functions are apparently different. Moreover, we have provided a complete geometric characterization of the objective, in contrast to . We believe the geometric characterization could not only provide insight to the algorithm, but also help improve the algorithm in terms of stability and also finding all components. . While this interesting connection potentially helps port our analysis to ICA, it is a fundamental question to ask what is playing a more vital role for DR, sparsity or independence.

Fig. 5 helps shed some light in this direction, where we again plot the asymptotic objective landscape with the natural reparameterization as in Section I-D2. From the left and central panels, it is evident that even without independence, $\mathbf{X}_{0}$ with sparse columns induces the familiar geometric structures we saw in Fig. 3; such structures are broken when the sparsity level becomes large. We believe all our later analyses can be generalized to the correlated cases we experimented with. On the other hand, from the right panelWe have not showed the results on the BG model here, as it seems the structure persists even when $\theta$ approaches $1$. We suspect the “phase transition” of the landscape occurs at different points for different distributions and Gaussian is the outlying case where the transition occurs at $1$. , it seems that with independence, the function landscape undergoes a transition, as sparsity level grows: target solution goes from minimizers of the objective to the maximizers of the objective. Without adequate knowledge of the true sparsity, it is unclear whether one would like to minimize or maximize the objective.For solving the ICA problem, this suggests the $\log\cosh$ contrast function, that works well empirically , may not work for all distributions (rotation-invariant Gaussian excluded of course), at least when one does not process the data (say perform certain whitening or scaling). This suggests that sparsity, instead of independence, makes our current algorithm for DR work.

Besides ICA discussed above, it turns out that a handful of other practical problems arising in signal processing and machine learning induce the “no spurious minimizers, all saddles are second-order” structure under natural setting, including the eigenvalue problem, generalized phase retrieval , orthogonal tensor decomposition , low-rank matrix recovery/completion , noisy phase synchronization and community detection , linear neural nets learning . gave a review of these problems, and discussed how the methodology developed in this and the companion paper can be generalized to solve those problems.

I-F Notations, and Reproducible Research

The codes to reproduce all the figures and experimental results are available online:

II The High-dimensional Function Landscape

In particular, we focus our attention to the smaller set

Suppose $\mathbf{A}_{0}=\mathbf{I}$ and hence $\mathbf{Y}=\mathbf{A}_{0}\mathbf{X}_{0}=\mathbf{X}_{0}$. There exist positive constants $c_{\star}$ and $C$, such that for any $\theta\in(0,1/2)$ and $\mu<c_{a}\min\left\{\theta n^{-1},n^{-5/4}\right\}$, whenever

the following hold simultaneously with probability at least $1-c_{b}p^{-6}$:

and the function $g(\mathbf{w};\mathbf{X}_{0})$ has exactly one local minimizer $\mathbf{w}_{\star}$ over the open set $\Gamma\doteq\left\{\mathbf{w}:\left\|\mathbf{w}\right\|<\sqrt{\tfrac{4n-1}{4n}}\right\}$, which satisfies

Here $c_{a}$ through $c_{c}$ are all positive constants.

Here $c_{a}$ to $c_{c}$ are positive constants.

By Theorem II.1, over $\mathbf{q}\left(\Gamma\right)$, $\mathbf{q}\left(\mathbf{w}_{\star}\right)$ is the unique local minimizer. Suppose not. Then there exist $\mathbf{q}^{\prime}\in\mathbf{q}\left(\Gamma\right)$ with $\mathbf{q}^{\prime}\neq\mathbf{q}\left(\mathbf{w}_{\star}\right)$ and $\varepsilon>0$, such that $f\left(\mathbf{q}^{\prime};\mathbf{X}_{0}\right)\leq f\left(\mathbf{q};\mathbf{X}_{0}\right)$ for all $\mathbf{q}\in\mathbf{q}\left(\Gamma\right)$ satisfying $\left\|\mathbf{q}^{\prime}-\mathbf{q}\right\|<\varepsilon$. Since the mapping $\mathbf{w}\mapsto\mathbf{q}\left(\mathbf{w}\right)$ is $2\sqrt{n}$-Lipschitz (Lemma IV.8), $g\left(\mathbf{w}\left(\mathbf{q}^{\prime}\right);\mathbf{X}_{0}\right)\leq g\left(\mathbf{w}\left(\mathbf{q}\right);\mathbf{X}_{0}\right)$ for all $\mathbf{w}\in\Gamma$ satisfying $\left\|\mathbf{w}\left(\mathbf{q}^{\prime}\right)-\mathbf{w}\left(\mathbf{q}\right)\right\|<\varepsilon/\left(2\sqrt{n}\right)$, implying $\mathbf{w}\left(\mathbf{q}^{\prime}\right)$ is a local minimizer different from $\mathbf{w}_{\star}$, a contradiction. Let $\left\|\mathbf{w}_{\star}-\mathbf{0}\right\|=\eta$. Straightforward calculation shows

Repeating the argument $2n$ times in the vicinity of other signed basis vectors $\pm\mathbf{e}_{i}$ gives $2n$ local minimizers of $f$. Indeed, the $2n$ symmetric sections cover the sphere with certain overlaps. We claim that none of the $2n$ local minimizers lies in the overlapped regions. This is due to the nearness of these local minimizers to standard basis vectors. To see this, w.l.o.g., suppose $\mathbf{q}$, which is the local minimizer next to $\mathbf{e}_{n}$, is in the overlapped region determined by $e_{n}$ and $e_{i}$ for some $i\neq n$. This implies that

by the definition of our symmetric sections. On the other hand, we know

Thus, so long as $1-\eta^{2}\geq\tfrac{4n-1}{4n}$, or $\eta\leq 1/(2\sqrt{n})$, a contradiction arises. Since $\eta\in O(\mu)$ and $\mu\leq O(n^{-5/4})$ by our assumption, our claim is confirmed. There are no extra local minimizers, as any extra local minimizer must be contained in at least one of the $2n$ symmetric sections, making two different local minimizers in one section, contradicting the uniqueness result we obtained above. ∎

Though the $2n$ isolated local minimizers may have different objective values, they are equally good in the sense each of them helps produce a close approximation to a certain row of $\mathbf{X}_{0}$. As discussed in Section I-D2, for cases $\mathbf{A}_{0}$ is an orthobasis other than $\mathbf{I}$, the landscape of $f\left(\mathbf{q};\mathbf{Y}\right)$ is simply a rotated version of the one we characterized above.

Suppose $\mathbf{A}_{0}$ is complete with its condition number $\kappa\left(\mathbf{A}_{0}\right)$. There exist positive constants $c_{\star}$ (particularly, the same constant as in Theorem II.1) and $C$, such that for any $\theta\in(0,1/2)$ and $\mu<c_{a}\min\left\{\theta n^{-1},n^{-5/4}\right\}$, when

and $\overline{\mathbf{Y}}\doteq\sqrt{p\theta}\left(\mathbf{Y}\mathbf{Y}^{*}\right)^{-1/2}\mathbf{Y}$, $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{*}=\mathtt{SVD}\left(\mathbf{A}_{0}\right)$, the following hold simultaneously with probability at least $1-c_{b}p^{-6}$:

and the function $g(\mathbf{w};\mathbf{V}\mathbf{U}^{*}\overline{\mathbf{Y}})$ has exactly one local minimizer $\mathbf{w}_{\star}$ over the open set $\Gamma\doteq\left\{\mathbf{w}:\left\|\mathbf{w}\right\|<\sqrt{\tfrac{4n-1}{4n}}\right\}$, which satisfies

Here $c_{a},a_{b}$ are both positive constants.

Here $c_{a},c_{b}$ are both positive constants.

II-B Useful Technical Lemmas and Proof Ideas for Orthogonal Dictionaries

Proving Theorem II.1 is conceptually straightforward: one shows that the expectation of each quantity of interest has the claimed property, and then proves that each quantity concentrates uniformly about its expectation. The detailed calculations are nontrivial.

The next three propositions show that in the expected function landscape, we see successively strongly convex region, large gradient region, and negative directional curvature region when moving away from zero, as depicted in Fig. 3 and sketched in Section I-D2.

For any $\theta\in(0,1/2)$, if $\mu\leq c\min\left\{\theta n^{-1},n^{-5/4}\right\}$, it holds for all $\mathbf{w}$ with $1/\left(20\sqrt{5}\right)\leq\left\|\mathbf{w}\right\|\leq\sqrt{(4n-1)/(4n)}$ that

For any $\theta\in(0,1/2)$, if $\mu\leq 9/50$, it holds for all $\mathbf{w}$ with $\mu/(4\sqrt{2})\leq\left\|\mathbf{w}\right\|\leq 1/(20\sqrt{5})$ that

For any $\theta\in(0,1/2)$, if $\mu\leq 1/(20\sqrt{n})$, it holds for all $\mathbf{w}$ with $\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$ that

To prove that the above hold qualitatively for finite $p$, i.e., the function $g\left(\mathbf{w};\mathbf{X}_{0}\right)$, we will need first prove that for a fixed $\mathbf{w}$ each of the quantity of interest concentrates about their expectation w.h.p., and the function is nice enough (Lipschitz) such that we can extend the results to all $\mathbf{w}$ via a discretization argument. The next three propositions provide the desired pointwise concentration results.

For every $\mathbf{w}\in\Gamma$, it holds that for any $t>0$,

Suppose $0<\mu\leq 1/\sqrt{n}$. For every $\mathbf{w}\in\Gamma$, it holds that for any $t>0$,

Suppose $0<\mu\leq 1/\sqrt{n}$. For every $\mathbf{w}\in\Gamma\cap\left\{\mathbf{w}:\left\|\mathbf{w}\right\|\leq 1/4\right\}$, it holds that for any $t>0$,

The next three propositions provide the desired Lipschitz results.

Fix any $r_{\fgecap}\in\left(0,1\right)$. Over the set $\Gamma\cap\left\{\mathbf{w}:\left\|\mathbf{w}\right\|\geq r_{\fgecap}\right\}$, $\mathbf{w}^{*}\nabla^{2}g(\mathbf{w};\mathbf{X}_{0})\mathbf{w}/\left\|\mathbf{w}\right\|^{2}$ is $L_{\fgecap}$-Lipschitz with

Fix any $r_{g}\in\left(0,1\right)$. Over the set $\Gamma\cap\left\{\mathbf{w}:\left\|\mathbf{w}\right\|\geq r_{g}\right\}$, $\mathbf{w}^{*}\nabla g(\mathbf{w};\mathbf{X}_{0})/\left\|\mathbf{w}\right\|$ is $L_{g}$-Lipschitz with

Fix any $r_{\fgecup}\in(0,1/2)$. Over the set $\Gamma\cap\left\{\mathbf{w}:\left\|\mathbf{w}\right\|\leq r_{\fgecup}\right\}$, $\nabla^{2}g(\mathbf{w};\mathbf{X}_{0})$ is $L_{\fgecup}$-Lipschitz with

Integrating the above pieces, Section IV-B provides a complete proof of Theorem II.1.

II-C Extending to Complete Dictionaries

As hinted in Section I-D2, instead of proving things from scratch, we build on the results we have obtained for orthogonal dictionaries. In particular, we will work with the preconditioned data matrix

and show that the function landscape $f\left(\mathbf{q};\overline{\mathbf{Y}}\right)$ looks qualitatively like that of orthogonal dictionaries (up to a global rotation), provided that $p$ is large enough.

The next lemma shows $\overline{\mathbf{Y}}$ can be treated as being generated from an orthobasis with the same BG coefficients, plus small noise.

for a certain $\mathbf{\Xi}$ obeying $\left\|\mathbf{\Xi}\right\|\leq 20\kappa^{4}\left(\mathbf{A}\right)\sqrt{\frac{\theta n\log p}{p}}$, with probability at least $1-p^{-8}$. Here $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{*}=\mathtt{SVD}\left(\mathbf{A}_{0}\right)$, and $C>0$ is a constant.

Notice that $\mathbf{U}\mathbf{V}^{*}$ above is orthogonal, and that landscape of $f(\mathbf{q};\overline{\mathbf{Y}})$ is simply a rotated version of that of $f(\mathbf{q};\mathbf{V}\mathbf{U}^{*}\overline{\mathbf{Y}})$, or using the notation in the above lemma, that of $f(\mathbf{q};\mathbf{X}_{0}+\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}\mathbf{X}_{0})=f(\mathbf{q};\mathbf{X}_{0}+\widetilde{\mathbf{\Xi}}\mathbf{X}_{0})$ with $\widetilde{\mathbf{\Xi}}\doteq\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}$. So similar to the orthogonal case, it is enough to consider this “canonical” case, and its “canonical” reparametrization:

The following lemma provides quantitative comparison between the gradient and Hessian of $g\left(\mathbf{w};\mathbf{X}_{0}+\widetilde{\mathbf{\Xi}}\mathbf{X}_{0}\right)$ and that of $g\left(\mathbf{w};\mathbf{X}_{0}\right)$.

with probability at least $1-\theta\left(np\right)^{-7}-\exp\left(-0.3\theta np\right)$. Here $C_{a},C_{b}$ are positive constants.

Combining the above two lemmas, it is easy to see when $p$ is large enough, $\|\widetilde{\mathbf{\Xi}}\|=\left\|\mathbf{\Xi}\right\|$ is then small enough (Lemma II.14), and hence changes to the gradient and Hessian caused by the perturbation are small. This gives the results presented in Theorem II.3; see Section IV-C for a detailed proof. In particular, for the $p$ chosen in Theorem II.3, it holds that

for a certain constant $c$ which can be made arbitrarily small by making the constant $C$ in $p$ large.

III Discussion

It is possible to extend the current analysis to other dictionary settings. Our geometric structures (and algorithms) allow plug-and-play noise analysis. Nevertheless, we believe a more stable way of dealing with noise is to directly extract the whole dictionary, i.e., to consider geometry and optimization (and perturbation) over the orthogonal group. This will require additional nontrivial technical work, but likely feasible thanks to the relatively complete knowledge of the orthogonal group . A substantial leap forward would be to extend the methodology to recovery of structured overcomplete dictionaries, such as tight frames. Though there is no natural elimination of one variable, one can consider the marginalization of the objective function w.r.t. the coefficients and work with implicit functions. This recent work on overcomplete DR has used a similar idea. The marginalization taken there is near to the global optimum of one variable, where the function is well-behaved. Studying the global properties of the marginalization may introduce additional challenges. For the coefficient model, as we alluded to in Section I-E, our analysis and results likely can be carried through to coefficients with statistical dependence and physical constraints.

The connection to ICA we discussed in Section I-E suggests our geometric characterization and algorithms can be modified for the ICA problem. This likely will provide new theoretical insights and computational schemes to ICA. In the surge of theoretical understanding of nonconvex heuristics , the initialization plus local refinement strategy mostly differs from practice, whereby random initializations seem to work well, and the analytic techniques developed in that line are mostly fragmented and highly specialized. The analytic and algorithmic framework we developed here holds promise to providing a coherent account of these problems, see . In particular, we have intentionally separated the geometric characterization and algorithm development, hoping to making both parts modular. It is interesting to see how far we can streamline the geometric characterization. Moreover, the separation allows development of more provable and practical algorithms, say in the direction of .

IV Proofs of Technical Results

The proof involves some delicate analysis, particularly polynomial approximation of the function $f\left(t\right)=1/\left(1+t\right)^{2}$ over $t\in\left[0,1\right]$. This is naturally induced by the $1-\tanh^{2}\left(\cdot\right)$ function. The next lemma characterizes one polynomial approximation of $f\left(t\right)$.

In particular, one can choose $b_{k}=(-1)^{k}(k+1)\beta^{k}$ with $\beta=1-1/\sqrt{T}<1$ such that

Moreover, such sequence satisfies $0<\sum_{k=0}^{\infty}\frac{b_{k}}{(1+k)^{3}}<\sum_{k=0}^{\infty}\frac{\left|b_{k}\right|}{(1+k)^{3}}<2$.

Let $X\sim\mathcal{N}\left(0,\sigma_{X}^{2}\right)$ and $Y\sim\mathcal{N}\left(0,\sigma_{Y}^{2}\right)$ be independent. We have

For $X+Y\geq 0$, let $Z=\exp\left(-2(X+Y)/\mu\right)\in$, then

First fix any $T>1$. By Lemma IV.1, we choose the polynomial $p_{\beta}\left(Z\right)=\frac{1}{\left(1+\beta Z\right)^{2}}$ with $\beta=1-1/\sqrt{T}$ to upper bound $f\left(Z\right)=\frac{1}{\left(1+Z\right)^{2}}$. So we have

where $b_{k}=(-1)^{k}(k+1)\beta^{k}$, and the exchange of infinite summation and expectation above is justified due to

and the dominated convergence theorem (see, e.g., theorem 2.24 and 2.25 of ). By Lemma B.1, we have

where we have applied Type I upper and lower bounds for $\Phi^{c}\left(\cdot\right)$ to even $k$ and odd $k$ respectively and rearranged the terms to obtain the last line. Using the following estimates (see Lemma IV.1)

Since the above holds for any $T>1$, we obtain the claimed result by letting $T\to\infty$ such that $\beta\to 1$. ∎

Let $X\sim\mathcal{N}\left(0,\sigma_{X}^{2}\right)$ and $Y\sim\mathcal{N}\left(0,\sigma_{Y}^{2}\right)$ be independent. We have

For $X+Y>0$, let $z=\exp\left(-2(X+Y)/\mu\right)\in$, then

First fix any $T>1$. By Lemma IV.1, we choose the polynomial $p_{\beta}\left(Z\right)=\frac{1}{\left(1+\beta Z\right)^{2}}$ with $\beta=1-1/\sqrt{T}$ to upper bound $f\left(Z\right)=\frac{1}{\left(1+Z\right)^{2}}$. So we have

where $b_{k}=(-1)^{k}(k+1)\beta^{k}$, and exchange of infinite summation and expectation is justified, due to

and the dominated convergence theorem (see, e.g., theorem 2.24 and 2.25 of ). By Lemma B.1, we have

where we have applied Type I lower and upper bounds for $\Phi^{c}\left(\cdot\right)$ to even $k$ and odd $k$ respectively and rearranged the terms to obtain the last line. Using the following estimates (see Lemma IV.1)

For the second term, by Lemma IV.1, we have

Since the above bound holds for any $T>1$, we obtain the claimed result by letting $T\to\infty$ such that $\beta\to 1$ and $1/\sqrt{T}\to 0$. ∎

Let $X\sim\mathcal{N}\left(0,\sigma_{X}^{2}\right)$ and $Y\sim\mathcal{N}\left(0,\sigma_{Y}^{2}\right)$ be independent. We have

Similar to the proof of the above lemma, for $X+Y>0$, let $Z\doteq\exp\left(-2\frac{X+Y}{\mu}\right)$ and $f\left(Z\right)\doteq\frac{1}{\left(1+Z\right)^{2}}$. First fix any $T>1$. We will use $4zp_{\beta}\left(Z\right)=\frac{4Z}{\left(1+\beta Z\right)^{2}}$ to approximate the $1-\tanh^{2}\left(\frac{X+Y}{\mu}\right)=4Zf\left(Z\right)$ function from above, where again $\beta=1-1/\sqrt{T}$. So we obtain

where we have applied Type I upper and lower bounds for $\Phi^{c}\left(\cdot\right)$ to odd $k$ and even $k$ respectively and rearranged the terms to obtain the last line. Using the following estimates (see Lemma IV.1)

For the second term, by Lemma B.1 and Lemma IV.1, we have

where we have also used Type I upper bound for $\Phi^{c}\left(\cdot\right)$. Combining the above estimates, we get

Since the above holds for any $T>1$, we obtain the claimed result by letting $T\to\infty$, such that $\beta\to 1$ and $1/\sqrt{T}\to 0$. ∎

(of Proposition II.5) For any $i\in[n-1]$, we have

The above holds for any pair of $i,j\in[n-1]$, so it follows that

An upper bound for $(\mathcal{A})$. When $x_{n}$ is not in support set of $\mathbf{x}$, the term reduces to

where to obtain the last line we used that $t\mapsto\exp(-2t/\mu)t^{2}$ for $t>0$ is maximized at $\mu$.

When $x_{n}$ is in the support set, we expand the square term inside the expectation and obtain

where conditioned on each support set $\mathcal{J}$, we let $X\doteq q_{n}\left(\mathbf{w}\right)v_{n}\sim\mathcal{N}\left(0,q_{n}^{2}\left(\mathbf{w}\right)\right)$ and $Y\doteq\mathbf{w}_{\mathcal{J}}^{*}\overline{\mathbf{v}}\sim\mathcal{N}\left(0,\left\|\mathbf{w}_{\mathcal{J}}\right\|^{2}\right)$. An upper bound for the above is obtained by calling the estimates in Lemma IV.2 and Lemma IV.3:

where we have used $\mu<q_{n}\left(\mathbf{w}\right)\leq\left\|\mathbf{q}_{\mathcal{I}}\right\|$ and $\left\|\mathbf{w}_{\mathcal{J}}\right\|\leq\left\|\mathbf{q}_{\mathcal{I}}\right\|$ and $\left\|\mathbf{w}\right\|\leq 1$ and $\theta\in\left(0,1/2\right)$ to simplify the intermediate quantities to obtain the last line.

A lower bound for $(\mathcal{B})$. Similarly, we obtain

Collecting the above estimates, we obtain

where to obtain the last line we have invoked the association inequality in Lemma A.2, as both $\left\|\mathbf{w}_{\mathcal{J}^{c}}\right\|^{2}$ and $1/\left\|\mathbf{q}_{\mathcal{I}}\right\|^{3}$ both coordinatewise nonincreasing w.r.t. the index set. Substituting the upper bound for $\mu$ into (IV.4) and noting $q_{n}\left(\mathbf{w}\right)\geq 1/(2\sqrt{n})$ (implied by the assumption $\left\|\mathbf{w}\right\|\leq\sqrt{(4n-1)/(4n)}$), we obtain the claimed result. ∎

IV-A2 Proof of Proposition II.6

By similar consideration as proof of the above proposition, the following is justified:

A lower bound for $(\mathcal{A})$. We have

where $X\doteq q_{n}\left(\mathbf{w}\right)v_{n}\sim\mathcal{N}\left(0,q_{n}^{2}\left(\mathbf{w}\right)\right)$ and $Y\doteq\mathbf{w}^{*}_{\mathcal{J}}\overline{\mathbf{v}}\sim\mathcal{N}\left(0,\left\|\mathbf{w}_{\mathcal{J}}\right\|^{2}\right)$. Now by Lemma A.2 we obtain

as $\tanh\left(\frac{X+Y}{\mu}\right)$ and $X$ are both coordinatewise nondecreasing function of $X$ and $Y$. Using $\tanh\left(z\right)\geq\left(1-\exp\left(-2z\right)\right)/2$ for all $z\geq 0$ and integral results in Lemma B.1, we obtain

where at the second line we have used the assumption that $\left\|\mathbf{w}\right\|\geq\mu/(4\sqrt{2})$ and also the fact that $\sqrt{1+x^{2}}\geq x+\frac{1}{10x}$ for $x\geq 1/(4\sqrt{2})$.

An upper bound for $(\mathcal{B})$. We have

as $\tanh\left(\cdot\right)$ is bounded by one in magnitude

Plugging the results of (IV.6) and (IV.7) into (IV.5) and noticing that $q_{n}\left(\mathbf{w}\right)^{2}+\left\|\mathbf{w}\right\|^{2}=1$ we obtain

where we have used $\frac{2\left\|\mathbf{w}\right\|}{\sqrt{1-\left\|\mathbf{w}\right\|^{2}}}\leq\frac{1}{10}\left(1-\theta\right)$ when $\left\|\mathbf{w}\right\|\leq 1/(20\sqrt{5})$ and $\theta\leq 1/2$, completing the proof. ∎

IV-A3 Proof of Proposition II.7

By consideration similar to proof of Proposition II.5, we can exchange the Hessian and expectation, i.e.,

We are interested in the expected Hessian matrix

in the region that $0\leq\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$.

When $\mathbf{w}=\mathbf{0}$, by Lemma B.1, we have

Simple calculation based on Lemma B.1 shows

Invoking the assumptions $\mu\leq 1/(20\sqrt{n})\leq 1/20$ and $\theta<1/2$, we obtain

When $0<\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$, we aim to derive a semidefinite lower bound for

We will first provide bounds for $(\mathcal{C})$ and $(\mathcal{D})$, which are relatively simple. Then we will bound $(\mathcal{A})$ and $(\mathcal{B})$, which are slightly more tricky.

An upper bound for $(\mathcal{C})$. We have

where to obtain the final bound we have invoked the assumptions: $\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$, $\mu\leq 1/(20\sqrt{n})$, and $\theta\leq 1/2$.

A lower bound for $(\mathcal{D})$. We directly drop the first expectation which is positive, and derive an upper for the second expectation term as:

where we have again used $\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$, $\mu\leq 1/(20\sqrt{n})$, and $q_{n}\left(\mathbf{w}\right)\geq 1/(2\sqrt{n})$ to obtain the final bound.

An upper bound for $(\mathcal{B})$. Similar to the way we bound $(\mathcal{D})$,

A lower bound for $(\mathcal{A})$. First note that

Thus, we set out to lower bound the expectation as

To see the above implication, suppose the latter claimed holds, then for any $\mathbf{z}$ with unit norm,

Now for any fixed support set $\mathcal{S}\subset[k]$, $\mathbf{z}=\mathcal{P}_{\widetilde{\mathbf{w}}_{\mathcal{S}}}\mathbf{z}+\left(\mathbf{I}-\mathcal{P}_{\widetilde{\mathbf{w}}_{\mathcal{S}}}\right)\mathbf{z}$. So we have

Using expectation result from Lemma B.1, the above lower bound is further bounded as:

where to obtain the last line we have used $\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$. On the other hand, we similarly obtain

So we can take $\beta=\frac{1}{\sqrt{2\pi}}\left(2-\frac{3}{4}\sqrt{2}\right)<1$.

Putting together the above estimates for the case $\mathbf{w}\neq\mathbf{0}$, we obtain

Hence for all $\mathbf{w}$, we can take the $\frac{1}{5\sqrt{2\pi}}\frac{\theta}{\mu}$ as the lower bound, completing the proof. ∎

IV-A4 Proof of Pointwise Concentration Results

We first establish a useful comparison lemma between random i.i.d. Bernoulli random vectors random i.i.d. normal random vectors.

Now, we are ready to prove Proposition II.8 to Proposition II.10 as follows.

then $\mathbf{w}^{*}\nabla g(\mathbf{w})/\left\|\mathbf{w}\right\|=\frac{1}{p}\sum_{k=1}^{p}X_{k}$. For each $X_{k},k\in[p]$, from (IV.1), we know that

as the magnitude of $\tanh\left(\cdot\right)$ is bounded by one. Because

invoking Lemma IV.5, we obtain for every integer $m\geq 2$ that

then $\mathbf{w}^{*}\nabla^{2}g(\mathbf{w})\mathbf{w}/\left\|\mathbf{w}\right\|^{2}=\frac{1}{p}\sum_{k=1}^{p}Y_{k}$. For each $Y_{k}$ ($k\in[p]$), from (IV.2), we know that

Then by similar argument as in proof to Proposition II.8, we have for all integers $m\geq 2$ that

provided that $\mu\leq 1/\sqrt{n}$, as desired. ∎

(of Proposition II.10) Let $\mathbf{Z}_{k}=\nabla^{2}_{\mathbf{w}}h_{\mu}\left(\mathbf{q}(\mathbf{w})^{*}(\mathbf{x}_{0})_{k}\right)$, then $\nabla^{2}_{\mathbf{w}}g\left(\mathbf{w}\right)=\frac{1}{p}\sum_{k=1}^{p}\mathbf{Z}_{k}$. From (IV.2), we know that

where we have used the fact that $\left\|\mathbf{w}\right\|^{2}/q_{n}^{2}(\mathbf{w})=\left\|\mathbf{w}\right\|^{2}/(1-\left\|\mathbf{w}\right\|^{2})\leq 1$ for $\left\|\mathbf{w}\right\|_{2}\leq 1/4$ and Lemma IV.5 to obtain the last line. By Lemma A.6, we obtain

where we have simplified the final result using $\mu\leq 1/\sqrt{n}$. ∎

IV-A5 Proof of Lipschitz Results

We need the following lemmas to prove the Lipschitz results.

Suppose that $\varphi_{1}:U\to V$ is an $L$-Lipschitz map from a normed space $U$ to a normed space $V$, and that $\varphi_{2}:V\to W$ is an $L^{\prime}$-Lipschitz map from $V$ to a normed space $W$. Then the composition $\varphi_{2}\circ\varphi_{1}:U\to W$ is $LL^{\prime}$-Lipschitz.

For every $\mathbf{w},\mathbf{w}^{\prime}\in\Gamma$, and every fixed $\mathbf{x}$, we have

where we have used the fact $q_{n}\left(\mathbf{w}\right)\geq 1/(2\sqrt{n})$ to get the final result. Hence the mapping $\mathbf{w}\mapsto\mathbf{q}(\mathbf{w})$ is $2\sqrt{n}$-Lipschitz over $\Gamma$. Moreover it is easy to see $\mathbf{q}\mapsto\mathbf{q}^{*}\mathbf{x}$ is $\left\|\mathbf{x}\right\|_{2}$-Lipschitz. By Lemma A.1 and the composition rule in Lemma IV.6, we obtain the desired claims. ∎

For any fixed $\mathbf{x}$, consider the function

defined over $\mathbf{w}\in\Gamma$. Then, for all $\mathbf{w},\mathbf{w}^{\prime}$ in $\Gamma$ such that $\left\|\mathbf{w}\right\|\geq r$ and $\left\|\mathbf{w}^{\prime}\right\|\geq r$ for any constant $r\in\left(0,1\right)$, it holds that

where we have used the assumption that $q_{n}\left(\mathbf{w}\right)\geq 1/(2\sqrt{n})$ to simplify the final result. The claim about $\left|t_{\mathbf{x}}^{2}\left(\mathbf{w}\right)\right|$ follows immediately. Now

where we have used the assumption that $\left\|\mathbf{w}\right\|\geq r$ to simplify the result. Noticing that $t\mapsto t/\sqrt{1-t^{2}}$ is continuous over $\left[a,b\right]$ and differentiable over $\left(a,b\right)$ for any $0<a<b<1$, by mean value theorem,

where we have again used the assumption that $q_{n}\left(\mathbf{w}\right)\geq 1/(2\sqrt{n})$ to simplify the last result. Collecting the above estimates, we obtain

leading to the claimed result once we substitute estimates of the involved quantities. ∎

For any fixed $\mathbf{x}$, consider the function

defined over $\mathbf{w}\in\Gamma$. Then, for all $\mathbf{w},\mathbf{w}^{\prime}\in\Gamma$ such that $\left\|\mathbf{w}\right\|<r$ and $\left\|\mathbf{w}^{\prime}\right\|<r$ with any constant $r\in\left(0,1/2\right)$, it holds that

where we have applied the estimate for $\left|q_{n}\left(\mathbf{w}\right)-q_{n}\left(\mathbf{w}^{\prime}\right)\right|$ as established in Lemma IV.8 and also used $\left\|\mathbf{w}\right\|\leq 1/2$ and $\left\|\mathbf{w}^{\prime}\right\|\leq 1/2$ to simplify the above result. Further noticing $t\mapsto t^{2}/\left(1-t^{2}\right)^{3/2}$ is differentiable over $t\in\left(0,1\right)$, we apply the mean value theorem and obtain

Combining the above estimates gives the claimed result. ∎

For any fixed $\mathbf{x}$, consider the function

defined over $\mathbf{w}\in\Gamma$. Then, for all $\mathbf{w},\mathbf{w}^{\prime}\in\Gamma$ such that $\left\|\mathbf{w}\right\|\leq r$ and $\left\|\mathbf{w}^{\prime}\right\|\leq r$ for any constant $r\in\left(0,1/2\right)$, it holds that

We have $\left\|\mathbf{w}\right\|^{2}/q_{n}^{2}\left(\mathbf{w}\right)\leq 1/3$ when $\left\|\mathbf{w}\right\|\leq r<1/2$, hence it holds that

Now, we are ready to prove all the Lipschitz propositions.

Then, $\mathbf{w}^{*}\nabla^{2}g(\mathbf{w};\mathbf{X}_{0})\mathbf{w}/\left\|\mathbf{w}\right\|^{2}=\frac{1}{p}\sum_{k=1}^{p}F_{k}(\mathbf{w})$. Noticing that $\ddot{h}_{\mu}\left(\mathbf{q}(\mathbf{w})^{*}(\mathbf{x}_{0})_{k}\right)$ is bounded by $1/\mu$ and $\dot{h}_{\mu}\left(\mathbf{q}(\mathbf{w})^{*}(\mathbf{x}_{0})_{k}\right)$ is bounded by $1$, both in magnitude. Applying Lemma IV.7, Lemma IV.8 and Lemma IV.9, we can see $F_{k}(\mathbf{w})$ is $L_{\fgecap}^{k}$-Lipschitz with

Thus, $\frac{1}{\left\|\mathbf{w}\right\|_{2}}\mathbf{w}^{*}\nabla^{2}g(\mathbf{w};\mathbf{X}_{0})\mathbf{w}$ is $L_{\fgecap}$-Lipschitz with

where $\dot{h}_{\mu}(t)=\tanh(t/\mu)$ is bounded by one in magnitude, and $t_{(\mathbf{x}_{0})_{k}}(\mathbf{w})$ and $t_{(\mathbf{x}_{0})_{k}^{\prime}}(\mathbf{w})$ is defined as in Lemma IV.9. By Lemma IV.7, Lemma IV.8 and Lemma IV.9, we know that $\dot{h}_{\mu}\left(\mathbf{q}(\mathbf{w})^{*}(\mathbf{x}_{0})_{k}\right)t_{(\mathbf{x}_{0})_{k}}\left(\mathbf{w}\right)$ is $L_{k}$-Lipschitz with constant

with $\mathbf{\zeta}_{k}(\mathbf{w})=\overline{\mathbf{x}_{0}}_{k}-\frac{x_{0k}\left(n\right)}{q_{n}(\mathbf{w})}\mathbf{w}$ and $\mathbf{\Phi}_{k}(\mathbf{w})=\frac{x_{0k}\left(n\right)}{q_{n}(\mathbf{w})}\mathbf{I}+\frac{x_{0k}(n)}{q_{n}^{3}(\mathbf{w})}\mathbf{w}\mathbf{w}^{*}$. Then, $\nabla^{2}g(\mathbf{w})=\frac{1}{p}\sum_{k=1}^{p}\mathbf{F}_{k}(\mathbf{w})$. Using Lemma IV.7, Lemma IV.8, Lemma IV.10 and Lemma IV.11, and the facts that $\ddot{h}_{\mu}(t)$ is bounded by $1/\mu$ and that $\ddot{h}_{\mu}(t)$ is bounded by $1$ in magnitude, we can see $\mathbf{F}_{k}(\mathbf{w})$ is $L_{\fgecup}^{k}$-Lipschitz continuous with

IV-B Proofs of Theorem II.1

Before proving Theorem II.1, we record one useful lemma.

For convenience, we define three regions for the range of $\mathbf{w}$:

Lipschitz by Proposition II.13. Set $\varepsilon=\frac{c_{1}\theta}{3\mu L_{1}}$, so

On $\mathcal{E}_{1}\cap\mathcal{E}_{\infty}$,

and so on $\mathcal{E}_{1}\cap\mathcal{E}_{\infty}$, (II.5) holds for any constant $c_{\star}\leq c_{1}/3$. Setting $t=c_{1}\theta/3\mu$ in Proposition II.10, we obtain that for any fixed $\mathbf{w}$,

Similarly, for the gradient quantity, for $\mathbf{w}\in R_{2}$, Proposition II.6 shows that

Moreover, on $\mathcal{E}_{\infty}$, $\mathbf{w}^{*}\nabla g(\mathbf{w};\mathbf{X}_{0})/\left\|\mathbf{w}\right\|$ is

Lipschitz by Proposition II.12. For any $\varepsilon<\frac{1}{20\sqrt{5}}$, the set $R_{2}$ has an $\varepsilon$-net $N_{2}$ of size at most $\left(\frac{3}{20\varepsilon\sqrt{5}}\right)^{n}$. Set $\varepsilon=\frac{c_{6}\theta}{3L_{2}}$, so

On $\mathcal{E}_{2}\cap\mathcal{E}_{\infty}$,

and so on $\mathcal{E}_{2}\cap\mathcal{E}_{\infty}$, (II.6) holds for any constant $c_{\star}\leq c_{6}/3$. Setting $t=c_{6}\theta/3$ in Proposition II.8, we obtain that for any fixed $\mathbf{w}\in R_{2}$,

Finally, for any $\mathbf{w}\in R_{3}$, Proposition II.5 shows that

On $\mathcal{E}_{\infty}$, $\mathbf{w}^{*}\nabla^{2}g(\mathbf{w};\mathbf{X}_{0})\mathbf{w}/\left\|\mathbf{w}\right\|^{2}$ is

Lipschitz by Proposition II.11. As above, for any $\varepsilon\leq\sqrt{\frac{4n-1}{4n}}$, $R_{3}$ has an $\varepsilon$-net $N_{3}$ of size at most $(3/\varepsilon)^{n}$. Set $\varepsilon=c_{9}\theta/3L_{3}$. Then

On $\mathcal{E}_{3}\cap\mathcal{E}_{\infty}$,

and (II.7) holds with any constant $c_{\star}<c_{9}/3$. Setting $t=c_{9}\theta/3$ in Proposition II.9 and taking a union bound, we obtain

Let $\mathcal{E}_{g}$ be the event that the bounds (II.5)-(II.7) hold. On $\mathcal{E}_{g}$, the function $g$ is $\frac{c_{\star}\theta}{\mu}$-strongly convex over $R_{1}=\left\{\mathbf{w}:\left\|\mathbf{w}\right\|\leq\mu/\left(4\sqrt{2}\right)\right\}$. This implies that $f$ has at most one local minimum on $R_{1}$. It also implies that for any $\mathbf{w}\in R_{1}$,

So, if $g(\mathbf{w};\mathbf{X}_{0})\leq g(\mathbf{0};\mathbf{X}_{0})$, we necessarily have

Then $g(\mathbf{w};\mathbf{X}_{0})\leq g(\mathbf{0};\mathbf{X}_{0})$ implies that $\left\|\mathbf{w}\right\|\leq\mu/16$. By Wierstrass’s theorem, $g(\mathbf{w};\mathbf{X}_{0})$ has at least one minimizer $\mathbf{w}_{\star}$ over the compact set $S=\left\{\mathbf{w}:\left\|\mathbf{w}\right\|\leq\mu/10\right\}$. By the above reasoning, $\left\|\mathbf{w}_{\star}\right\|\leq\mu/16$, and hence $\mathbf{w}_{\star}$ does not lie on the boundary of $S$. This implies that $\mathbf{w}_{\star}$ is a local minimizer of $g$. Moreover, as above,

We now use the vector Bernstein inequality to show that with our choice of $p$, (IV.11) is satisfied w.h.p. Notice that

and $\dot{h}_{\mu}$ is bounded by one in magnitude, so for any integer $m\geq 2$,

where we have applied the moment estimate for the $\chi\left(n\right)$ distribution shown in Lemma A.7. Applying the vector Bernstein inequality in Corollary A.10 with $R=\sqrt{n}$ and $\sigma^{2}=2n$, we obtain

for all $t>0$. Using this inequality, it is not difficult to show that there exist constants $C_{13},C_{14}>0$ such that when $p\geq C_{13}n\log n$, with probability at least $1-4np^{-10}$,

When $\frac{p}{\log p}\geq\frac{C_{14}n}{\theta^{2}}$, for appropriately large $C_{14}$, (IV.12) implies (IV.11). Summing up failure probabilities completes the proof. ∎

IV-C Proofs for Section II-C and Theorem II.3

(of Lemma II.14) By the generative model,

On the other hand, by Lemma B.3, when $p\geq C_{1}n^{2}\log n$, $\left\|\tfrac{1}{p\theta}\mathbf{X}_{0}\mathbf{X}_{0}^{*}-\mathbf{I}\right\|\leq 10\sqrt{\tfrac{\theta n\log p}{p}}$ with probability at least $1-p^{-8}$. Thus, when $p\geq C_{2}\kappa^{4}\left(\mathbf{A}_{0}\right)\theta n^{2}\log(n\theta\kappa\left(\mathbf{A}_{0}\right))$,

where $L_{\dot{h}_{\mu}}$ denotes the Lipschitz constant for $\dot{h}_{\mu}\left(\cdot\right)$. Similarly, suppose $\left\|\widetilde{\mathbf{\Xi}}\right\|\leq 1/(2n)$, and also notice that

where $L_{\ddot{h}_{\mu}}$ denotes the Lipschitz constant for $\ddot{h}_{\mu}\left(\cdot\right)$. Since

and by Lemma IV.12, $\left\|\mathbf{X}\right\|_{\infty}\leq 4\sqrt{\log\left(np\right)}$ with probability at least $1-\theta\left(np\right)^{-7}-\exp\left(-0.3\theta np\right)$, we obtain

(of Theorem II.3) Here $c_{\star}$ is as defined in Theorem II.1. By Lemma II.14, when

the magnitude of the perturbation is bounded as

where $C_{2}$ can be made arbitrarily small by making $C_{1}$ large. Combining this result with Lemma II.15, we obtain that for all $\mathbf{w}\in\Gamma$,

with probability at least $1-p^{-8}-\theta\left(np\right)^{-7}-\exp\left(-0.3\theta np\right)$. In view of (II.13) in Theorem II.1, we have