Complete Dictionary Recovery over the Sphere II: Recovery by Riemannian Trust-region Method

I Introduction

Recently, there is a surge of research studying nonconvex formulations and provable algorithms for a number of central problems in signal processing and machine learning, including, e.g., low-rank matrix completion/recovery , phase retreival , tensor recovery , mixed regression , structured element pursuit , blind deconvolution , noisy phase synchronization and community detection , deep learning , numerical linear algebra and optimization . The research efforts are fruitful in producing more practical and scalable algorithms and even significantly better performance guarantees than known convex methods.

where $\widehat{\mathbf{Y}}$ is a proxy of $\mathbf{Y}$ (i.e., after appropriate processing), $\widehat{\mathbf{y}}_{k}$ is the $k$-th column of $\widehat{\mathbf{Y}}$, and $h_{\mu}(z)\doteq\mu\log\cosh(z/\mu)$ is a (convex) smooth approximation to the absolute-value function. The spherical constraint renders the problem nonconvex.

Despite the apparent nonconvexity, our prior analysis in has showed that all local minimizers of (I.1) are qualitatively equally good, because each of them produces a close approximation to certain row of $\mathbf{X}_{0}$ (Corollary II.4 in ). So the central issue is how to escape from saddle points. Fortunately, our previous results (Theorem II.3 in ) imply that all saddle points under consideration are ridable, i.e., the associated Hessians have both strictly positive and strictly negative values (see the recapitulation in Section II-B). Particularly, eigenvectors of the negative eigenvalues are direction of negative curvature, which intuitively serve as directions of local descent.

Second-order methods can naturally exploit the curvature information to escape from ridable saddle points. To gain some intuition, consider an unconstrained optimization problem

The second-order Taylor expansion of $\phi$ at a saddle point $\mathbf{x}_{0}$ is

Thus, minimizing $\widehat{\phi}(\mathbf{\delta};\mathbf{x}_{0})$ returns a direction $\mathbf{\delta}_{\star}$ that tends to decrease the objective $\phi$, provided local approximation of $\widehat{\phi}$ to $\phi$ is reasonably accurate. Based on this intuition, we derive a (second-order) Riemannian trust-region algorithm that exploits the second-order information to escape from saddle points and provably returns a local minimizer to (I.1), from arbitrary initializations. We provide rigorous guarantees for recovering a local minimizer in Section II.

Obtaining a local minimizer only helps approximate one row of $\mathbf{X}_{0}$. To recover the row, we derive a simple linear programming rounding procedure that provably works. To recover all rows of $\mathbf{X}_{0}$, one repeats the above process based on a carefully designed deflation process. The whole algorithmic pipeline and the related recovery guarantees are provided in Section III. Particularly, we show that when $p$ is reasonably large, with high probability (w.h.p.), our pipeline efficiently recovers $\mathbf{A}_{0}$ and $\mathbf{X}_{0}$, even when each column of $\mathbf{X}_{0}$ contains $O(n)$ nonzeros.

In Section II.E of the companion paper , we provide detailed comparisons of our results with prior theoretical results on DR; we conclude that this is the first algorithmic framework that guarantees efficient recovery of complete dictionaries when the coefficients have up to constant fraction of nonzeros. We also draw methodological connections to work on understanding nonconvex heuristics, and other nonconvex problems with similar geometric structures. Here we focus on drawing detailed connections to the optimization literature.

Trust-region method (TRM) has a rich history dating back to 40’s; see the monograph for accounts of the history and developments. The main motivation for early developments was to address limitations of the classic Newton’s method (see, e.g., Section 3 of ). The limitations include the technical subtleties to establish local and global convergence results. Moreover, when the Hessian is singular or indefinite, the movement direction is either not well-defined, or does not improve the objective function. initialized the line of work that addresses the limitations. Particularly, proposed using local second-order Taylor approximation as model function in the trust-region framework for unconstrained optimization. They showed that under mild conditions, the trust-region iterate sequence has a limit point that is critical and has positive semidefinite Hessian; see also Section 6.5-6.6 of . Upon inspecting the relevant proofs, it seems not hard to strengthen the results to sequence convergence to local minimizers, under a ridable saddle condition as ours, for unconstrained optimization.

Research activities to port theories and algorithms of optimization in Euclidean space to Riemannian manifolds are best summarized by three monographs: . developed Newton and conjugate-gradient methods for the Stiefel manifolds, of which the sphere is a special case; presents a complete set of first- and second-order Riemannian algorithms and convergence analyses; see also the excellent associated optimization software toolbox . Among these, trust-region method was first ported to the Riemannian setting in , with emphasis on efficient implementation which only approximately solves the trust-region subproblem according to the Cauchy point scheme. The Cauchy point definition adopted there was the usual form based on the gradient, not strong enough to ensure the algorithm escape from ridable saddle points even if the true Hessian is in use in local approximation. In comparison, in this work we assume that the trust-region subproblem is exactly solved, such that ridable saddles (the only possible saddles for our problem) are properly skipped. By this, we obtain the strong guarantee that the iterate sequence converges to a local minimizer, in contrast to the weak global convergence (gradient sequence converging to zero) or local convergence (sequence converging to a local minimizer within a small radius) established in . To the best of our knowledge, our convergence result is first of its kind for a specific problem on sphere. After our initial submission, has recently established worst-case iteration complexity of Riemannian TRM to converge to second-order critical points (i.e., critical points with positive semidefinite Hessians), echoing the results in the Euclidean case . Their results are under mild Lipschitz-type assumptions and allow inexact subproblem solvers, and hence are very practical and general. However, on our particular problem, their result is considerably pessimistic, compared to our convergence result obtained from a specialized analysis.

Solving the trust-region subproblem exactly is expensive. Practically, often a reasonable approximate solution with controlled quality is adequate to guarantee convergence. In this regard, the truncated conjugate gradient (tCG) solver with a good initial search direction is commonly employed in practice (see, e.g., Section 7.5 in ). To ensure ridable saddle points be properly escaped from, the eigenpoint idea (see, e.g., Section 6.6 of ) is particularly relevant; see also Algorithm 3 and Lemma 10 in .

The benign function landscape we characterized in the first paper allows any reasonable iterative method that is capable of escaping from ridable saddles to find a local minimizer, with possibly different performance guarantees. The trust-region method we focus on here, and the curviliear search method are second-order methods that guarantee global optimization from arbitrary initializations. Typical first-order methods such as the vanilla gradient descent can only guarantee convergence to a critical point. Nonetheless, for our particular function, noisy/stochastic gradient method guarantees to find a local minimizer from an arbitrary initialization with high probability .

I-B Notations, and Reproducible Research

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

II Finding One Local Minimizer via the Riemannian Trust-Region Method

We are interested to seek a local minimizer of (I.1). The presence of saddle points have motivated us to develop a second-order Riemannian trust-region algorithm over the sphere; the existence of descent directions at nonoptimal points drives the trust-region iteration sequence towards one of the minimizers asymptotically. We will prove that under our modeling assumptions, this algorithm with an arbitrary initialization efficiently produces an accurate approximationBy “accurate” we mean one can achieve an arbitrary numerical accuracy $\varepsilon>0$ with a reasonable amount of time. Here the running time of the algorithm is on the order of $\log\log(1/\varepsilon)$ in the target accuracy $\varepsilon$, and polynomial in other problem parameters. to one of the minimizers. Throughout the exposition, basic knowledge of Riemannian geometry is assumed. We will try to keep the technical requirement minimal possible; the reader can consult the excellent monograph for relevant background and details.

defines a smooth curve on the sphere that satisfies $\gamma(0)=\mathbf{q}$ and $\dot{\gamma}(0)=\mathbf{\delta}$. Geometrically, $\gamma(t)$ is a segment of the great circle that passes $\mathbf{q}$ and has $\mathbf{\delta}$ as its tangent vector at $\mathbf{q}$. The exponential map for $\mathbf{\delta}$ is defined as

It is a canonical way of pulling $\mathbf{\delta}$ to the sphere.

Thus, the above quadratic approximation can be rewritten compactly as

If $\operatorname{Hess}f(\mathbf{q};\widehat{\mathbf{Y}})$ is positive semidefinite and has “full rank” $n-1$ (hence “nondegenerate”Note that the $n\times n$ matrix $\operatorname{Hess}f(\mathbf{q};\widehat{\mathbf{Y}})$ has rank at most $n-1$, as the nonzero $\mathbf{q}$ obviously is in its null space. When $\operatorname{Hess}f(\mathbf{q};\widehat{\mathbf{Y}})$ has rank $n-1$, it has no null direction in the tangent space. Thus, in this case it acts on the tangent space like a full-rank matrix. ), the unique solution $\mathbf{\delta}_{\star}$ is

II-B The Geometric Results from [3]

Geometrically, this corresponds to projection of the function $f$ above the equatorial section $\mathbf{e}_{n}^{\perp}$ onto $\mathbf{e}_{n}^{\perp}$ (see Fig. 2 (right) for illustration). In particular, we focus our attention to the smaller set of the ball:

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.

Recall that the reason we just need to characterize the geometry for the case $\mathbf{A}_{0}=\mathbf{I}$ is that for other orthogonal $\mathbf{A}_{0}$, the function landscape is simply a rotated version of that of $\mathbf{A}_{0}=\mathbf{I}$.

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},c_{b}$ are both positive constants.

It is very easy to verify the following fact (see proof of Lemma II.13 on page VI-G):

II-C The Riemannian Trust-Region Algorithm over the Sphere

Solution to (II.21) can then be recovered as $\widehat{\mathbf{\delta}}=\mathbf{U}\widehat{\mathbf{\xi}}$.

The problem (II.23) is an instance of the classic trust region subproblem, i.e., minimizing a quadratic function subject to a single quadratic constraint. Albeit potentially nonconvex, this notable subproblem can be solved in polynomial time by several numerical methods . Approximate solution of the subproblem suffices to guarantee convergence in theory, and lessens the storage and computational burden in practice. We will deploy the approximate version in simulations. For simplicity, however, our subsequent analysis assumes the subproblem is solved exactly. We next briefly describe how one can deploy the semidefinite programming (SDP) approach to solve the subproblem exactly. This choice is due to the well-known effectiveness and robustness of the SDP approach on this problem. We introduce

where $\mathbf{A}=\mathbf{U}^{*}\operatorname{Hess}f(\mathbf{q}^{(r-1)};\widehat{\mathbf{Y}})\mathbf{U}$ and $\mathbf{b}=\mathbf{U}^{*}\operatorname{grad}\nabla f(\mathbf{q}^{(r-1)};\widehat{\mathbf{Y}})$. The resulting SDP to solve is

where $\mathbf{E}_{n+1}=\mathbf{e}_{n+1}\mathbf{e}_{n+1}^{*}$. Once the problem (II.27) is solved to its optimum $\mathbf{\Theta}_{\star}$, one can provably recover the minimizer $\mathbf{\xi}_{\star}$ of (II.23) by computing the SVD of $\mathbf{\Theta}_{\star}=\widetilde{\mathbf{U}}\mathbf{\Sigma}\widetilde{\mathbf{V}}^{*}$, and extract as a subvector the first $n-1$ coordinates of the principal eigenvector $\widetilde{\mathbf{u}}_{1}$ (see Appendix B of ).

II-D Main Convergence Results

Using general convergence results on Riemannian TRM (see, e.g., Chapter 7 of ), it is not difficult to prove that the gradient sequence $\operatorname{grad}f(q^{(r)};\widehat{\mathbf{Y}})$ produced by TRM converges to zero (i.e., global convergence), or the sequence converges (at quadratic rate) to a local minimizer if the initialization is already close a local minimizer (i.e., local convergence). In this section, we show that under our probabilistic assumptions, these results can be substantially strengthened. In particular, the algorithm is guaranteed to produce an accurate approximation to a local minimizer of the objective function, in a number of iterations that is polynomial in the problem size, from arbitrary initializations. The arguments in the companion paper showed that w.h.p. every local minimizer of $f$ produces a close approximation to a row of $\mathbf{X}_{0}$. Taken together, this implies that the algorithm efficiently produces a close approximation to one row of $\mathbf{X}_{0}$.

Thorough the analysis, we assume the trust-region subproblem is exactly solved and the step size parameter $\Delta$ is fixed. Our next two theorems summarize the convergence results for orthogonal and complete dictionaries, respectively.

Suppose the dictionary $\mathbf{A}_{0}$ is orthogonal. There exists a positive constant $C$, such that for all $\theta\in\left(0,1/2\right)$ and $\mu<c_{a}\min\left\{\theta n^{-1},n^{-5/4}\right\}$, whenever

with probability at least $1-c_{b}p^{-6},$ the Riemannian trust-region algorithm with input data matrix $\widehat{\mathbf{Y}}=\mathbf{Y}$, any initialization $\mathbf{q}^{(0)}$ on the sphere, and a step size satisfying

iterations. Here $c_{\star}$ is as defined in Theorem II.1, and $c_{a}$ through $c_{f}$ are all positive constants.

Suppose the dictionary $\mathbf{A}_{0}$ is complete with condition number $\kappa\left(\mathbf{A}_{0}\right)$. There exists a positive constant $C$, such that for all $\theta\in\left(0,1/2\right)$, and $\mu<c_{a}\min\left\{\theta n^{-1},n^{-5/4}\right\}$, whenever

with probability at least $1-c_{b}p^{-6},$ the Riemannian trust-region algorithm with input data matrix $\overline{\mathbf{Y}}\doteq\sqrt{p\theta}\left(\mathbf{Y}\mathbf{Y}^{*}\right)^{-1/2}\mathbf{Y}$ where $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{*}=\mathtt{SVD}\left(\mathbf{A}_{0}\right)$, any initialization $\mathbf{q}^{(0)}$ on the sphere and a step size satisfying

iterations. Here $c_{\star}$ is as in Theorem II.1, and $c_{a}$ through $c_{f}$ are all positive constants.

Our convergence result shows that for any target accuracy $\varepsilon>0$ the algorithm terminates within polynomially many steps. Specifically, the first summand in (II.29) or (II.31) is the number of steps the sequence takes to enter the strongly convex region and be “reasonably” close to a local minimizer. All subsequent trust-region subproblems are then unconstrained (proved below) – the constraint is inactive at optimal point, and hence the steps behave like Newton steps. The second summand reflects the typical quadratic local convergence of the Newton steps.

Our estimate of the number of steps is pessimistic: the running time is a relatively high-degree polynomial in $p$ and $n$. We will discuss practical implementation details that help speed up in Section IV. Our goal in stating the above results is not to provide a tight analysis, but to prove that the Riemannian TRM algorithm finds a local minimizer in polynomial time. For nonconvex problems, this is not entirely trivial – results of show that in general it is NP-hard to find a local minimizer of a nonconvex function.

II-E Sketch of Proof for Orthogonal Dictionaries

Note that for any orthogonal $\mathbf{A}_{0}$, $f\left(\mathbf{q};\mathbf{A}_{0}\mathbf{X}_{0}\right)=f\left(\mathbf{A}_{0}^{*}\mathbf{q};\mathbf{X}_{0}\right)$. In words, this is the established fact that the function landscape of $f(\mathbf{q};\mathbf{A}_{0}\mathbf{X}_{0})$ is a rotated version of that of $f(\mathbf{q};\mathbf{X}_{0})$. Thus, any local minimizer $\mathbf{q}_{\star}$ of $f(\mathbf{q};\mathbf{X}_{0})$ is rotated to $\mathbf{A}_{0}\mathbf{q}_{\star}$, a local minimizer of $f(\mathbf{q};\mathbf{A}_{0}\mathbf{X}_{0})$. Also if our algorithm generates iteration sequence $\mathbf{q}_{0},\mathbf{q}_{1},\mathbf{q}_{2},\dots$ for $f(\mathbf{q};\mathbf{X}_{0})$ upon initialization $\mathbf{q}_{0}$, it will generate the iteration sequence $\mathbf{A}_{0}\mathbf{q}_{0},\mathbf{A}_{0}\mathbf{q}_{1},\mathbf{A}_{0}\mathbf{q}_{2},\dots$ for $f\left(\mathbf{q};\mathbf{A}_{0}\mathbf{X}_{0}\right)$. So w.l.o.g. it is adequate that we prove the convergence results for the case $\mathbf{A}_{0}=\mathbf{I}$. So in this section (Section II-E), we write $f(\mathbf{q})$ to mean $f(\mathbf{q};\mathbf{X}_{0})$.

We partition the sphere into three regions, for which we label as $R_{\mathtt{I}}$, $R_{\mathtt{II}}$, $R_{\mathtt{III}}$, corresponding to the strongly convex, nonzero gradient, and negative curvature regions, respectively (see Theorem II.1). That is, $R_{\mathtt{I}}$ consists of a union of $2n$ spherical caps of radius $\mu/(4\sqrt{2})$, each centered around a signed standard basis vector $\pm\mathbf{e}_{i}$. $R_{\mathtt{II}}$ consist of the set difference of a union of $2n$ spherical caps of radius $1/(20\sqrt{5})$, centered around the standard basis vectors $\pm\mathbf{e}_{i}$, and $R_{\mathtt{I}}$. Finally, $R_{\mathtt{III}}$ covers the rest of the sphere. We say a trust-region step takes an $R_{\mathtt{I}}$ step if the current iterate is in $R_{\mathtt{I}}$; similarly for $R_{\mathtt{II}}$ and $R_{\mathtt{III}}$ steps. Since we use the geometric structures derived in Theorem II.1 and Corollary II.2 in , the conditions

At step $r$ of the algorithm, suppose $\mathbf{\delta}^{(r)}$ is the minimizer of the trust-region subproblem (II.21). We call the step “constrained” if $\left\|\mathbf{\delta}^{(r)}\right\|=\Delta$ (the minimizer lies on the boundary and hence the constraint is active), and call it “unconstrained” if $\|\mathbf{\delta}^{(r)}\|<\Delta$ (the minimizer lies in the relative interior and hence the constraint is not in force). Thus, in the unconstrained case the optimality condition is (II.6).

The next lemma provides some estimates about $\nabla f$ and $\nabla^{2}f$ that are useful in various contexts.

We have the following estimates about $\nabla f$ and $\nabla^{2}f$:

Our next lemma says if the trust-region step size $\Delta$ is small enough, one Riemannian trust-region step reduces the objective value by a certain amount when there is any descent direction.

where $\eta_{f}\doteq M_{\nabla}+2M_{\nabla^{2}}+L_{\nabla}+L_{\nabla^{2}}$ and $M_{\nabla}$, $M_{\nabla^{2}}$, $L_{\nabla}$, $L_{\nabla^{2}}$ are the quantities defined in Lemma II.5.

To show decrease in objective value for $R_{\mathtt{II}}$ and $R_{\mathtt{III}}$, now it is enough to exhibit a descent direction for each point in these regions. The next two lemmas help us almost accomplish the goal. For convenience again we choose to state the results for the “canonical” section that is in the vicinity of $\mathbf{e}_{n}$ and the projection map $\mathbf{q}\left(\mathbf{w}\right)=[\mathbf{w};(1-\left\|\mathbf{w}\right\|^{2})^{1/2}]$, with the idea that similar statements hold for other symmetric sections.

One can take $\beta_{g}=\beta_{\fgecap}=c_{\star}\theta$ as shown in Theorem II.1, and take the Lipschitz results in Proposition B.4 and Proposition B.3 (note that $\left\|\mathbf{X}_{0}\right\|_{\infty}\leq 4\log^{1/2}(np)$ w.h.p. by Lemma B.6), repeat the argument for other $2n-1$ symmetric regions, and conclude that w.h.p. the objective value decreases by at least a constant amount. The next proposition summarizes the results.

Assume (II.32). In regions $R_{\mathtt{II}}$ and $R_{\mathtt{III}}$, each trust-region step reduces the objective value by at least

where $c_{a}$ to $c_{c}$ are positive constants, and $c_{\star}$ is as defined in Theorem II.1.

We only consider the symmetric section in the vicinity of $\mathbf{e}_{n}$ and the claims carry on to others by symmetry. If the current iterate $\mathbf{q}^{(r)}$ is in the region $R_{\mathtt{II}}$, by Theorem II.1, w.h.p., we have $\mathbf{w}^{*}g\left(\mathbf{w}\right)/\left\|\mathbf{w}\right\|\geq c_{\star}\theta$ for the constant $c_{\star}$. By Proposition B.4 and Lemma B.6, w.h.p., $\mathbf{w}^{*}g\left(\mathbf{w}\right)/\left\|\mathbf{w}\right\|$ is $C_{2}n^{2}\log\left(np\right)/\mu$-Lipschitz. Therefore, By Lemma II.6 and Lemma II.7, a trust-region step decreases the objective value by at least

Similarly, if $\mathbf{q}^{(r)}$ is in the region $R_{\mathtt{III}}$, by Proposition B.3, Theorem II.1 and Lemma B.6, w.h.p., $\mathbf{w}^{*}\nabla^{2}g\left(\mathbf{w}\right)\mathbf{w}/\left\|\mathbf{w}\right\|^{2}$ is $C_{3}n^{3}\log^{3/2}\left(np\right)/\mu^{2}$-Lipschitz and upper bounded by $-c_{\star}\theta$. By Lemma II.6 and Lemma II.8, a trust-region step decreases the objective value by at least

It can be easily verified that when $\Delta$ obeys (II.33), (II.34) holds. ∎

The analysis for $R_{\mathtt{I}}$ is slightly trickier. In this region, near each local minimizer, the objective function is strongly convex. So we still expect each trust-region step decreases the objective value. On the other hand, it is very unlikely that we can provide a universal lower bound for the amount of decrease - as the iteration sequence approaches a local minimizer, the movement is expected to be diminishing. Nevertheless, close to the minimizer the trust-region algorithm takes “unconstrained” steps. For constrained $R_{\mathtt{I}}$ steps, we will again show reduction in objective value by at least a fixed amount; for unconstrained step, we will show the distance between the iterate and the nearest local minimizer drops down rapidly.

The next lemma concerns the function value reduction for constrained $R_{\mathtt{I}}$ steps.

where $\eta_{f}$ is defined the same as Lemma II.6.

The next lemma provides an estimate of $m_{H}$. Again we will only state the result for the “canonical” section with the “canonical” $\mathbf{q}(\mathbf{w})$ mapping.

There exists a positive constant $C$, such that for all $\theta\in\left(0,1/2\right)$ and $\mu<\theta/10$, whenever $p\geq Cn^{3}\log\frac{n}{\theta\mu}/(\mu\theta^{2})$, it holds with probability at least $1-cp^{-7}$ that for all $\mathbf{q}$ with $\left\|\mathbf{w}\left(\mathbf{q}\right)\right\|\leq\mu/(4\sqrt{2})$,

Here $c_{\star}$ is as in Theorem II.1 and Theorem II.2, and $c>0$ is another constant.

We know that $\left\|\mathbf{X}_{0}\right\|_{\infty}\leq 4\log^{1/2}(np)$ w.h.p., and hence by the definition of Riemannian Hessian and Lemma II.5,

Combining this estimate and Lemma II.11 and Lemma II.6, we obtain a concrete lower bound for the reduction of objective value for each constrained $R_{\mathtt{I}}$ step.

Assume (II.32). Each constrained $R_{\mathtt{I}}$ trust-region step (i.e., $\left\|\mathbf{\delta}\right\|=\Delta$) reduces the objective value by at least

Here $c_{\star}$ is as in Theorem II.1 and Theorem II.2, and $c,c^{\prime}$ are positive constants.

We only consider the symmetric section in the vicinity of $\mathbf{e}_{n}$ and the claims carry on to others by symmetry. We have that w.h.p.

Combining these estimates with Lemma II.6 and Lemma II.10, one trust-region step will find next iterate $\mathbf{q}^{(r+1)}$ that decreases the objective value by at least

Finally, by the condition on $\Delta$ in (II.36) and the assumed conditions (II.32), we obtain

By the proof strategy for $R_{\mathtt{I}}$ we sketched before Lemma II.10, we expect the iteration sequence ultimately always takes unconstrained steps when it moves very close to a local minimizer. We will show that the following is true: when $\Delta$ is small enough, once the iteration sequence starts to take unconstrained $R_{\mathtt{I}}$ step, it will take consecutive unconstrained $R_{\mathtt{I}}$ steps afterwards. It takes two steps to show this: (1) upon an unconstrained $R_{\mathtt{I}}$ step, the next iterate will stay in $R_{\mathtt{I}}$. It is obvious we can make $\Delta\in O(1)$ to ensure the next iterate stays in $R_{\mathtt{I}}\cup R_{\mathtt{II}}$. To strengthen the result, we use the gradient information. From Theorem II.1, we expect the magnitudes of the gradients in $R_{\mathtt{II}}$ to be lower bounded; on the other hand, in $R_{\mathtt{I}}$ where points are near local minimizers, continuity argument implies that the magnitudes of gradients should be upper bounded. We will show that when $\Delta$ is small enough, there is a gap between these two bounds, implying the next iterate stays in $R_{\mathtt{I}}$; (2) when $\Delta$ is small enough, the step is in fact unconstrained. Again we will only state the result for the “canonical” section with the “canonical” $\mathbf{q}(\mathbf{w})$ mapping. The next lemma exhibits an absolute lower bound for magnitudes of gradients in $R_{\mathtt{II}}$.

For all $\mathbf{q}$ satisfying $\mu/(4\sqrt{2})\leq\left\|\mathbf{w}\left(\mathbf{q}\right)\right\|\leq 1/(20\sqrt{5})$, it holds that

Assuming (II.32), Theorem II.1 gives that w.h.p. $\mathbf{w}^{*}\nabla g(\mathbf{w})/\left\|\mathbf{w}\right\|\geq c_{\star}\theta$. Thus, w.h.p, $\left\|\operatorname{grad}f(\mathbf{q})\right\|\geq 9c_{\star}\theta/10$ for all $\mathbf{q}\in R_{\mathtt{II}}$. The next lemma compares the magnitudes of gradients before and after taking one unconstrained $R_{\mathtt{I}}$ step. This is crucial to providing upper bound for magnitude of gradient for the next iterate, and also to establishing the ultimate (quadratic) sequence convergence.

where $L_{H}\doteq 5n^{3/2}/({2\mu^{2}})\left\|\mathbf{X}_{0}\right\|_{\infty}^{3}+9n/\mu\left\|\mathbf{X}_{0}\right\|_{\infty}^{2}+9\sqrt{n}\left\|\mathbf{X}_{0}\right\|_{\infty}$.

We can now bound the Riemannian gradient of the next iterate as

Obviously, one can make the upper bound small by tuning down $\Delta$. Combining the above lower bound for $\left\|\operatorname{grad}f(\mathbf{q})\right\|$ for $\mathbf{q}\in R_{\mathtt{II}}$, one can conclude that when $\Delta$ is small, the next iterate $\mathbf{q}^{(r+1)}$ stays in $R_{\mathtt{I}}$. Another application of the optimality condition (II.6) gives conditions on $\Delta$ that guarantees the next trust-region step is also unconstrained. Detailed argument can be found in proof of the following proposition.

Assume (II.32). W.h.p, once the trust-region algorithm takes an unconstrained $R_{\mathtt{I}}$ step (i.e., $\left\|\mathbf{\delta}\right\|<\Delta$), it always takes unconstrained $R_{\mathtt{I}}$ steps, provided that

Here $c_{\star}$ is as in Theorem II.1 and Theorem II.2, and $c>0$ is another constant.

We only consider the symmetric section in the vicinity of $\mathbf{e}_{n}$ and the claims carry on to others by symmetry. Suppose that step $k$ is an unconstrained $R_{\mathtt{I}}$ step. Then

Thus, if $\Delta\leq\tfrac{1}{20\sqrt{5}}-\tfrac{\mu}{4\sqrt{2}}$, $\mathbf{q}^{(r+1)}$ will be in $R_{\mathtt{I}}\cup R_{\mathtt{II}}$. Next, we show that if $\Delta$ is sufficiently small, $\mathbf{q}^{(r+1)}$ will be indeed in $R_{\mathtt{I}}$. By Lemma II.14,

as the step is unconstrained. On the other hand, by Theorem II.1 and Lemma II.13, w.h.p.

we have $\mathbf{q}^{(r+1)}\in R_{\mathtt{I}}$.

We next show that when $\Delta$ is small enough, the next step is also unconstrained. Straight forward calculations give

in words, the minimizer to the trust-region subproblem for the next step lies in the relative interior of the trust region - the constraint is inactive. By Lemma II.14 and Lemma B.6, we have

w.h.p.. Combining this and our previous estimates of $m_{H}$, $M_{H}$, we conclude whenever

w.h.p., our next trust-region step is also an unconstrained $R_{\mathtt{I}}$ step. Simplifying the above bound completes the proof. ∎

Finally, we want to show that ultimate unconstrained $R_{\mathtt{I}}$ iterates actually converges to one nearby local minimizer rapidly. Lemma II.14 has established the gradient is diminishing. The next lemma shows the magnitude of gradient serves as a good proxy for distance to the local minimizer.

To see this relates the magnitude of gradient to the distance away from the nearby local minimizer, w.l.o.g., one can assume $\tau=1$ and consider the point $\mathbf{q}=\exp_{\mathbf{q}_{\star}}(\mathbf{\delta})$. Then

where at the last inequality above we have used Lemma II.16. Hence, combining this observation with Lemma II.14, we can derive the asymptotic sequence convergence rate as follows.

Assume (II.32) and the conditions in Lemma II.15. Let $\mathbf{q}^{(r_{0})}\in R_{\mathtt{I}}$ and the $r_{0}$-th step the first unconstrained $R_{\mathtt{I}}$ step and $\mathbf{q}_{\star}$ be the unique local minimizer of $f$ over one connected component of $R_{\mathtt{I}}$ that contains $\mathbf{q}^{(r_{0})}$. Then w.h.p., for any positive integer $r^{\prime}\geq 1$,

Here $c_{\star}$ is as in Theorem II.1 and Theorem II.2, and $c$, $c^{\prime}$ are both positive constants.

By the geometric characterization in Theorem II.1 and corollary II.2 in , $f$ has $2n$ separated local minimizers, each located in $R_{\mathtt{I}}$ and within distance $\sqrt{2}\mu/16$ of one of the $2n$ signed basis vectors $\{\pm\mathbf{e}_{i}\}_{i\in[n]}$. Moreover, it is obvious when $\mu\leq 1$, $R_{\mathtt{I}}$ consists of $2n$ disjoint connected components. We only consider the symmetric component in the vicinity of $\mathbf{e}_{n}$ and the claims carry on to others by symmetry.

Suppose that $r_{0}$ is the index of the first unconstrained iterate in region $R_{\mathtt{I}}$, i.e., $\mathbf{q}^{(r_{0})}\in R_{\mathtt{I}}$. By Lemma II.14, for any integer $r^{\prime}\geq 1$, we have

where $L_{H}$ is as defined in Lemma II.14, $m_{H}$ as the strong convexity parameter for $R_{\mathtt{I}}$ defined above.

Now suppose $\mathbf{q}_{\star}$ is the unique local minimizer of $f$, lies in the same $R_{\mathtt{I}}$ component that $q^{(r_{0})}$ is located. Let $\gamma_{r^{\prime}}(t)=\exp_{\mathbf{q}_{\star}}\left(t\mathbf{\delta}\right)$ to be the unique geodesic that connects $\mathbf{q}_{\star}$ and $\mathbf{q}^{(r_{0}+r^{\prime})}$ with $\gamma_{r^{\prime}}(0)=\mathbf{q}_{\star}$ and $\gamma_{r^{\prime}}(1)=\mathbf{q}^{(r_{0}+r^{\prime})}$. We have

where at the second line we have repeatedly applied Lemma II.16.

By the optimality condition (II.6) and the fact that $\left\|\mathbf{\delta}^{(r_{0})}\right\|<\Delta$, we have

we can combine the above results and obtain

Based on the previous estimates for $m_{H}$, $M_{H}$ and $L_{H}$, we obtain that w.h.p.,

Moreover, by (II.46), w.h.p., it is sufficient to have the trust region size

Now we are ready to piece together the above technical proposition to prove Theorem II.3.

(of Theorem II.3) Assuming (II.32) and in addition that

it can be verified that the conditions of all the above propositions are satisfied.

By the preceding four propositions, a step will either be $R_{\mathtt{III}}$, $R_{\mathtt{II}}$, or constrained $R_{\mathtt{I}}$ step that decreases the objective value by at least a certain fixed amount (we call this Type A), or be an unconstrained $R_{\mathtt{I}}$ step (Type B), such that all future steps are unconstrained $R_{\mathtt{I}}$ and the sequence converges to a local minimizer quadratically. Hence, regardless the initialization, the whole iteration sequence consists of consecutive Type A steps, followed by consecutive Type B steps. Depending on the initialization, either the Type A phase or the Type B phase can be absent. In any case, from $\mathbf{q}^{(0)}$ it takes at most (note $f(\mathbf{q})\geq 0$ always holds)

steps for the iterate sequence to start take consecutive unconstrained $R_{\mathtt{I}}$ steps, or to already terminate. In case the iterate sequence continues to take consecutive unconstrained $R_{\mathtt{I}}$ steps, Proposition II.17 implies that it takes at most

steps to obtain an $\varepsilon$-near solution to the $\mathbf{q}_{\star}$ that is contained in the connected subset of $R_{\mathtt{I}}$ that the sequence entered.

Thus, the number of iterations to obtain an $\varepsilon$-near solution to $\mathbf{q}_{\star}$ can be grossly bounded by

Finally, the claimed failure probability comes from a simple union bound with careful bookkeeping. ∎

II-F Extending to Convergence for Complete Dictionaries

Recall that in this case we consider the preconditioned input

Note that for any complete $\mathbf{A}_{0}$ with condition number $\kappa\left(\mathbf{A}_{0}\right)$, from Lemma B.7 we know when $p$ is large enough, w.h.p. one can write the preconditioned $\overline{\mathbf{Y}}$ as

for a certain $\mathbf{\Xi}$ with small magnitude, and $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{*}=\mathtt{SVD}\left(\mathbf{A}_{0}\right)$. Particularly, when $p$ is chosen by Theorem II.2, the perturbation is bounded as

for a certain constant $c$ which can be made arbitrarily small by making the constant $C$ in $p$ large. Since $\mathbf{U}\mathbf{V}^{*}$ is orthogonal,

In words, the function landscape of $f(\mathbf{q};\mathbf{U}\mathbf{V}^{*}\mathbf{X}_{0}+\mathbf{\Xi}\mathbf{X}_{0})$ is a rotated version of that of $f(\mathbf{q};\mathbf{X}_{0}+\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}\mathbf{X}_{0})$. Thus, any local minimizer $\mathbf{q}_{\star}$ of $f(\mathbf{q};\mathbf{X}_{0}+\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}\mathbf{X}_{0})$ is rotated to $\mathbf{U}\mathbf{V}^{*}\mathbf{q}_{\star}$, one minimizer of $f(\mathbf{q};\mathbf{U}\mathbf{V}^{*}\mathbf{X}_{0}+\mathbf{\Xi}\mathbf{X}_{0})$. Also if our algorithm generates iteration sequence $\mathbf{q}_{0},\mathbf{q}_{1},\mathbf{q}_{2},\dots$ for $f(\mathbf{q};\mathbf{X}_{0}+\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}\mathbf{X}_{0})$ upon initialization $\mathbf{q}_{0}$, it will generate the iteration sequence $\mathbf{U}\mathbf{V}^{*}\mathbf{q}_{0}$, $\mathbf{U}\mathbf{V}^{*}\mathbf{q}_{1}$, $\mathbf{U}\mathbf{V}^{*}\mathbf{q}_{2},\dots$ for $f\left(\mathbf{q};\mathbf{U}\mathbf{V}^{*}\mathbf{X}_{0}+\mathbf{\Xi}\mathbf{X}_{0}\right)$. So w.l.o.g. it is adequate that we prove the convergence results for the case $f(\mathbf{q};\mathbf{X}_{0}+\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}\mathbf{X}_{0})$, corresponding to $\bm{A}_{0}=\mathbf{I}$ with perturbation $\widetilde{\mathbf{\Xi}}\doteq\mathbf{V}\mathbf{U}^{*}\mathbf{\Xi}$. So in this section (Section II-F), we write $f(\mathbf{q};\widetilde{\mathbf{X}_{0}})$ to mean $f(\mathbf{q};\mathbf{X}_{0}+\widetilde{\mathbf{\Xi}}\mathbf{X}_{0})$.

the geometric structure of the landscape is qualitatively unchanged from the orthogonal case, and the parameter $c_{\star}$ constant can be replaced with $c_{\star}/2$. Particularly, for this choice of $p$, Lemma B.7 implies

for a constant $c$ that can be made arbitrarily small by setting the constant $C$ in $p$ sufficiently large. The whole proof is quite similar to that of orthogonal case in the last section. We will only sketch the major changes below. To distinguish with the corresponding quantities in the last section, we use $\widetilde{\cdot}$ to denote the corresponding perturbed quantities here.

where by (II.52) we have used $\|\widetilde{\mathbf{\Xi}}\|\leq 1/(2\sqrt{n})$ to simplify the above result. So we obtain

Lemma II.7 and Lemma II.8 are generic and nothing changes.

Proposition II.9: We have now $\mathbf{w}^{*}\mathbf{g}(\mathbf{w};\widetilde{\mathbf{X}_{0}})/\left\|\mathbf{w}\right\|\geq c_{\star}\theta/2$ by Theorem II.2, w.h.p. $\mathbf{w}^{*}\nabla g(\mathbf{w};\widetilde{\mathbf{X}_{0}})/\left\|\mathbf{w}\right\|$ is $C_{1}n^{2}\log(np)/\mu$-Lipschitz by Proposition B.4, and $\left\|\mathbf{X}_{0}+\widetilde{\mathbf{\Xi}}\mathbf{X}_{0}\right\|_{\infty}\leq 3\left\|\mathbf{X}_{0}\right\|_{\infty}/2$ as shown above. Similarly, $\mathbf{w}^{*}\mathbf{g}(\mathbf{w};\widetilde{\mathbf{X}_{0}})/\left\|\mathbf{w}\right\|\leq-c_{\star}\theta/2$ by Theorem II.2, and $\mathbf{w}^{*}\nabla^{2}g(\mathbf{w};\widetilde{\mathbf{X}_{0}})\mathbf{w}/\left\|\mathbf{w}\right\|^{2}$ is $C_{2}n^{3}\log^{3/2}(np)/\mu^{2}$-Lipschitz. Moreover, $\widetilde{\eta}_{f}\leq 4\eta_{f}$ as shown above. Since there are only multiplicative constant changes to the various quantities, we conclude

Lemma II.10: $\eta_{f}$ is changed to $\widetilde{\eta}_{f}$ with $\widetilde{\eta}_{f}\leq 4\eta_{f}$ as shown above.

where $L_{\ddot{h}}$ is the Lipschitz constant for the function $\ddot{h}_{\mu}\left(\cdot\right)$ and we have used the fact that $\|\widetilde{\mathbf{\Xi}}\|\leq 1$. Similarly, by II.2,

where $L_{\dot{h}}$ is the Lipschitz constant for the function $\dot{h}_{\mu}\left(\cdot\right)$. Since $L_{\ddot{h}}\leq 2/\mu^{2}$ and $L_{\dot{h}}\leq 1/\mu$, and $\left\|\mathbf{X}_{0}\right\|_{\infty}\leq 4\sqrt{\log(np)}$ w.h.p. (Lemma B.6). By (II.52), w.h.p. we have

provided the constant $C$ in (II.51) for $p$ is large enough. Thus, by (II.5) and the above estimates we have

Proposition II.12: From the estimate of $M_{H}$ above Proposition II.12 and the last point, we have

Also since $\widetilde{\eta}_{f}\leq 4\eta_{f}$ in Lemma II.6 and Lemma II.10, there are only multiplicative constant change to the various quantities. We conclude that

Lemma II.13 is generic and nothing changes.

Lemma II.14: $\widetilde{L}_{H}\leq 27L_{H}/8$.

Proposition II.15: All the quantities involved in determining $\Delta$, $m_{H}$, $M_{H}$, and $L_{H}$, $\beta_{\operatorname{grad}}$ are modified by at most constant multiplicative factors and changed to their respective tilde version, so we conclude that the TRM algorithm always takes unconstrained $R_{\mathtt{I}}$ step after taking one, provided that

Lemma II.16:is generic and nothing changes.

Proposition II.17: Again $m_{H}$, $M_{H}$, $L_{H}$ are changed to $\widetilde{m_{H}}$, $\widetilde{M_{H}}$, and $\widetilde{L_{H}}$, respectively, differing by at most constant multiplicative factors. So we conclude for any integer $k^{\prime}\geq 1$,

The final proof to Theorem II.2 is almost identical to that of Theorem II.1, except that $d_{\mathtt{I}}$, $d_{\mathtt{II}}$, and $d_{\mathtt{III}}$ are changed to $\widetilde{d_{\mathtt{I}}}$, $\widetilde{d_{\mathtt{II}}}$, and $\widetilde{d_{\mathtt{III}}}$ as defined above, respectively. The final iteration complexity to each an $\varepsilon$-near solution is hence

Hence overall the qualitative behavior of the algorithm is not changed, as compared to that for the orthogonal case.

III Complete Algorithm Pipeline and Main Results

Recovering one row of $\mathbf{X}_{0}$ by rounding. To obtain the target solution $\mathbf{q}_{\star}$ and hence recover (up to scale) one row of $\mathbf{X}_{0}$, we solve the following linear program:

with $\mathbf{r}=\widehat{\mathbf{q}}$. We show in Lemma III.2 (resp. Lemma III.4) that when $\left\langle\widehat{\mathbf{q}},\mathbf{q}_{\star}\right\rangle$ is sufficiently large, implied by $\mu$ being sufficiently small, w.h.p. the minimizer of (III.1) is exactly $\mathbf{q}_{\star}$, and hence one row of $\mathbf{X}_{0}$ is recovered by $\mathbf{q}_{\star}^{*}\widehat{\mathbf{Y}}$.

Reconstructing the dictionary $\mathbf{A}_{0}$. By solving the linear system $\mathbf{Y}=\mathbf{A}\mathbf{X}_{0}$, one can obtain the dictionary $\mathbf{A}_{0}=\mathbf{Y}\mathbf{X}_{0}^{*}\left(\mathbf{X}_{0}\mathbf{X}_{0}^{*}\right)^{-1}$.

Assume the dictionary $\mathbf{A}_{0}$ is orthogonal and we take $\widehat{\mathbf{Y}}=\mathbf{Y}$. Suppose $\theta\in\left(0,1/3\right)$, $\mu_{\star}<c_{a}\min\left\{\theta n^{-1},n^{-5/4}\right\}$, and $p\geq Cn^{3}\log\frac{n}{\mu_{\star}\theta}/\left(\mu_{\star}^{2}\theta^{2}\right)$. The above algorithmic pipeline with parameter setting

recovers the dictionary $\mathbf{A}_{0}$ and $\mathbf{X}_{0}$ in polynomial time, with failure probability bounded by $c_{c}p^{-6}$. Here $c_{\star}$ is as defined in Theorem II.1, and $c_{a}$ through $c_{c}$, and $C$ are all positive constants.

Towards a proof of the above theorem, it remains to be shown the correctness of the rounding and deflation procedures.

The following lemma shows w.h.p. the rounding will return the desired $\mathbf{q}_{\star}$, provided the estimated $\widehat{\mathbf{q}}$ is already near to it.

For any $\theta\in\left(0,1/3\right)$, whenever $p\geq Cn^{2}\log(n/\theta)/\theta$, with probability at least $1-cp^{-6},$ the rounding procedure (III.1) returns $\mathbf{q}_{\star}$ for any input vector $\mathbf{r}$ that satisfies

Since $\left\langle\widehat{\mathbf{q}},\mathbf{q}_{\star}\right\rangle=1-\|\widehat{\mathbf{q}}-\mathbf{q}_{\star}\|^{2}/2$, and $\left\|\widehat{\mathbf{q}}-\mathbf{q}_{\star}\right\|\in O(\mu)$, it is sufficient when $\mu$ is smaller than some small constant.

Proof sketch of deflation.

Thus, by Lemma III.2, one successfully recovers $\mathbf{U}\mathbf{z}_{\star}$ from $\mathbf{U}\widehat{\mathbf{z}}$ w.h.p. when $\mu_{\star}$ is smaller than a constant. The overall failure probability can be obtained via a simple union bound and simplification of the exponential tails with inverse polynomials in $p$.

III-B Recovering Complete Dictionaries

By working with the preconditioned data samples $\widehat{\mathbf{Y}}=\overline{\mathbf{Y}}\doteq\sqrt{\theta p}\left(\mathbf{Y}\mathbf{Y}^{*}\right)^{-1/2}\mathbf{Y}$,In practice, the parameter $\theta$ might not be know beforehand. However, because it only scales the problem, it does not affect the overall qualitative aspect of results. we can use the same procedure as described above to recover complete dictionaries.

Assume the dictionary $\mathbf{A}_{0}$ is complete with a condition number $\kappa\left(\mathbf{A}_{0}\right)$ and we take $\widehat{\mathbf{Y}}=\overline{\mathbf{Y}}$. Suppose $\theta\in\left(0,1/3\right)$, $\mu_{\star}<c_{a}\min\left\{\theta n^{-1},n^{-5/4}\right\}$, and $p\geq\frac{C}{c_{\star}^{2}\theta^{2}}\max\left\{\frac{n^{4}}{\mu^{4}},\frac{n^{5}}{\mu^{2}}\right\}\kappa^{8}\left(\mathbf{A}_{0}\right)\log^{4}\left(\frac{\kappa\left(\mathbf{A}_{0}\right)n}{\mu\theta}\right)$. The algorithmic pipeline with parameter setting

recovers the dictionary $\mathbf{A}_{0}$ and $\mathbf{X}_{0}$ in polynomial time, with failure probability bounded by $c_{b}p^{-6}$. Here $c_{\star}$ is as defined in Theorem II.1, and $c_{a},c_{b}$ are both positive constants.

Similar to the orthogonal case, we need to show the correctness of the rounding and deflation procedures so that the theorem above holds.

The result of the LP rounding is only slightly different from that of the orthogonal case in Lemma III.2, so is the proof.

For any $\theta\in\left(0,1/3\right)$, whenever

with probability at least $1-cp^{-6},$ the rounding procedure (III.1) returns $\mathbf{q}_{\star}$ for any input vector $\mathbf{r}$ that satisfies

Proof sketch of deflation.

By Lemma B.7 and (II.50), we can write $\overline{\mathbf{Y}}=(\mathbf{Q}+\mathbf{\Xi})\mathbf{X}_{0}$ for some orthogonal matrix $\mathbf{Q}$ and small perturbation $\mathbf{\Xi}$ with $\left\|\mathbf{\Xi}\right\|\leq\delta<1/10$ for some large $p$ as usual. Similar to the orthogonal case, we have

where $\mathbf{z}_{\star}$ with $\left\|\mathbf{z}_{\star}\right\|=1$ is the exact solution. More specifically, Corollary II.4 in implies that

Next, we show that $\widehat{\mathbf{z}}$ is also very near to the exact solution $\mathbf{z}_{\star}$. Indeed, the identity (III.8) suggests

where $\mathbf{W}^{\dagger}=(\mathbf{W}^{*}\mathbf{W})^{-1}\mathbf{W}^{*}$ denotes the pseudo inverse of a matrix $\mathbf{W}$ with full column rank. Hence, by (III.9) we can bound the distance between $\widehat{\mathbf{z}}_{\star}$ and $\mathbf{z}_{\star}$ by

By Lemma B.1, when $p\geq\Omega(n^{2}\log n)$, w.h.p.,

Hence, combined with Lemma III.5, we obtain

which implies that $\left\|\widehat{\mathbf{z}}_{\star}-\mathbf{z}_{\star}\right\|\leq 28\left\|\mathbf{\Xi}\right\|$. Thus, combining the results above, we obtain

Lemma B.7, and in particular (II.50), for our choice of $p$ as in Theorem II.2, $\left\|\mathbf{\Xi}\right\|\leq c\mu_{\star}^{2}n^{-3/2}$, where $c$ can be made smaller by making the constant in $p$ larger. For $\mu_{\star}$ sufficiently small, we conclude that

In words, the TRM algorithm returns a $\widehat{\mathbf{z}}$ such that $\mathbf{U}\widehat{\mathbf{z}}$ is very near to one of the unit vectors $\left\{\mathbf{q}_{\star}^{i}\right\}_{i=1}^{n}$, such that $(\mathbf{q}_{\star}^{i})^{*}\overline{\mathbf{Y}}=\alpha\mathbf{e}_{i}^{*}\mathbf{X}_{0}$ for some $\alpha\neq 0$. For $\mu_{\star}$ smaller than a fixed constant, one will have

and hence by Lemma III.4, the LP rounding exactly returns the optimal solution $\mathbf{q}_{\star}^{i}$ upon the input $\mathbf{U}\widehat{\mathbf{z}}$.

The proof sketch above explains why the recursive TRM plus rounding works. The overall failure probability can be obtained via a simple union bound and simplifications of the exponential tails with inverse polynomials in $p$.

IV Simulations

Fixing a small step size and solving the trust-region subproblem exactly eases the analysis, but also renders the TRM algorithm impractical. In practice, the trust-region subproblem is never exactly solved, and the trust-region step size is adjusted to the local geometry, say by backtracking. It is possible to modify our algorithmic analysis to account for inexact subproblem solvers and adaptive step size; for sake of brevity, we do not pursue it here. Recent theoretical results on the practical version include .

Here we describe a practical implementation based on the Manopt toolbox Available online: http://www.manopt.org. . Manopt is a user-friendly Matlab toolbox that implements several sophisticated solvers for tackling optimization problems over Riemannian manifolds. The most developed solver is based on the TRM. This solver uses the truncated conjugate gradient (tCG; see, e.g., Section 7.5.4 of ) method to (approximately) solve the trust-region subproblem (vs. the exact solver in our analysis). It also dynamically adjusts the step size using backtracking. However, the original implementation (Manopt 2.0) is not adequate for our purposes. Their tCG solver uses the gradient as the initial search direction, which does not ensure that the TRM solver can escape from saddle points . We modify the tCG solver, such that when the current gradient is small and there is a negative curvature direction (i.e., the current point is near a saddle point or a local maximizer of $f(\mathbf{q})$), the tCG solver explicitly uses the negative curvature direction…adjusted in sign to ensure positive correlation with the gradient – if it does not vanish. as the initial search direction. This modification ensures the TRM solver always escape from saddle points/local maximizers with negative directional curvature. Hence, the modified TRM algorithm based on Manopt is expected to have the same qualitative behavior as the idealized version we analyzed above, with better scalability. We will perform our numerical simulations using the modified TRM algorithm whenever necessary. Algorithm 3 together with Lemmas 9 and 10 and the surrounding discussion in the very recent work provides a detailed description of this practical version.

IV-B Simulated Data

To corroborate our theory, we experiment with dictionary recovery on simulated data.The code is available online: https://github.com/sunju/dl_focm For simplicity, we focus on recovering orthogonal dictionaries and we declare success once a single row of the coefficient matrix is recovered.

One trial is determined to be a success once $\mathtt{RE}\leq\mu$, with the idea that this indicates $\widehat{\mathbf{q}}$ is already very near the target and the target can likely be recovered via the LP rounding we described (which we do not implement here).

We consider two settings: (1) fix $p=5n^{2}\log n$ and vary the dimension $n$ and sparsity $k$; (2) fix the sparsity level as $\lceil 0.2\cdot n\rceil$ and vary the dimension $n$ and number of samples $p$. For each pair of $(k,n)$ for (1), and each pair of $(p,n)$ for (2), we repeat the simulations independently for $T=5$ times.

Fig. 3 shows the phase transition for the two settings. It seems that our TRM algorithm can work well into the linear region whenever $p\in O(n^{2}\log n)$ (Fig. 3-Left), but $p$ should have order greater than $\Omega(n)$ (Fig. 3-Right). The sample complexity from our theory is significantly suboptimal compared to this.

IV-C Image Data Again

Our algorithmic framework has been derived based on the BG model on the coefficients. Real data may not admit sparse representations w.r.t. complete dictionaries, or even so, the coefficients may not obey the BG model. In this experiment, we explore how our algorithm performs in learning complete dictionaries for image patches, emulating our motivational experiment in the companion paper (Section I.B). Thanks to research on image compression, we know patches of natural images tend to admit sparse representation, even w.r.t. simple orthogonal bases, such as Fourier basis or wavelets.

We take the three images that we used in the motivational experiment. For each image, we divide it into $8\times 8$ non-overlapping patches, vectorize the patches, and then stack the vectorized patches into a data matrix $\mathbf{Y}$. $\mathbf{Y}$ is preconditioned as

and the resulting $\overline{\mathbf{Y}}$ is fed to the dictionary learning pipeline described in Section III. The smoothing parameter $\mu$ is fixed to $10^{-2}$. Fig. 4 contains the learned dictionaries: the dictionaries generally contain localized, directional features that resemble subset of wavelets and generalizations. These are very reasonable representing elements for natural images. Thus, the BG coefficient model may be a sensible, simple model for natural images.

Another piece of strong evidence in support of the above claim is as follows. For each image, we repeat the learning pipeline for one hundred times, with independent initializations across the runs. Let $\widehat{\mathbf{A}}$ be the final learned dictionary for each run, we plot the value of $\|\widehat{\mathbf{A}}^{-1}\mathbf{Y}\|_{1}$ across the one hundred independent runs. Strikingly, the values are virtually the same, with a relative difference of $10^{-3}$! This is predicted by our theory, under the BG model. If the model is unreasonable for natural images, the preconditioning, benign function landscape, LP rounding, and the deflation process that hinge on this model would have completely fallen down.

For this image experiment, $n=64$ and $p=4096$. A single run of the learning pipeline, including solving $64$ instances of the optimization over the sphere (with varying dimensions) and solving $64$ instances of the LP rounding (using CVX), lasts about $20$ minutes on a mid-range modern laptop. So with careful implementation we discussed above, the learning pipeline is actually not far from practical.

V Discussion

Our experiments seem to suggest the necessary complexity level lies between $\Omega(n)$ and $O(n^{2}\log n)$ even for the orthogonal case. While it is interesting to determine the true complexity requirement for the TRM, there could be other efficient algorithms that demand less. For example, simulations in seem to suggest $O(n\log n)$ samples suffice to guarantee efficient recovery. The simulations run an alternating direction algorithm fed with problem-specific initializations, nevertheless.

Our analysis is based on exact trust-region subproblem solver and fixed step size. The convergence result for the practical version from , based on approximate solver and adaptive step size, is general, but pessimistic. It seems not difficult to adapt their analysis according to our objective geometry, and obtain a tight, practical convergence result.

Our motivating experiment on real images in introduction of our companion paper remains mysterious. If we were to believe that real image data are “nice” and our objective there does not have spurious local minima either, it is surprising ADM would escape all other critical points – this is not predicted by classic or modern theories. One reasonable place to start is to look at how gradient descent algorithms with generic initializations can escape from ridable saddle points (at least with high probability). The recent work has showed that randomly perturbing each iterate can help gradient algorithm to achieve this with high probability.

VI Proof of Convergence for the Trust-Region Algorithm

Using the fact $\tanh\left(\cdot\right)$ and $1-\tanh^{2}\left(\cdot\right)$ are bounded by one in magnitude, by (II.2) and (II.3) we have

where at the last line we have used the fact the mapping $\mathbf{q}\mapsto\mathbf{q}^{*}(\mathbf{x}_{0})_{k}/\mu$ is $\left\|(\mathbf{x}_{0})_{k}\right\|/\mu$ Lipschitz, and $x\mapsto\tanh\left(x\right)$ is $1$-Lipschitz, and the composition rule of Lipschitz functions (i.e., Lemma V.5 of ). Similar argument yields the final bound. ∎

VI-B Proof of Lemma II.6

as claimed. Next we establish the first result. Let $\mathbf{\delta}_{0}=\frac{\mathbf{\delta}}{\left\|\mathbf{\delta}\right\|}$, and $t=\left\|\mathbf{\delta}\right\|$. Consider the composite function

We next develop a bound on $\left|\ddot{\zeta}(t)-\ddot{\zeta}(0)\right|$. Using the triangle inequality, we can casually bound this difference as

where in the final line we have used the fact $1-\cos x=2\sin^{2}\left(x/2\right)$ and that $\sin x\leq x$ for $x\in\left[0,1\right]$, and $M_{\nabla}$, $M_{\nabla^{2}}$, $L_{\nabla}$ and $L_{\nabla^{2}}$ are the quantities defined in Lemma II.5. By the integral form of Taylor’s theorem in Lemma A.7 and the result above, we have

with $t=\left\|\mathbf{\delta}\right\|$ we obtain the desired result. ∎

VI-C Proof of Lemma II.7

By the integral form of Taylor’s theorem in Lemma A.7, for any $t\in\left[0,\frac{3\Delta}{2\pi\sqrt{n}}\right]$, we have

Minimizing this function over $t\in\left[0,\frac{3\Delta}{2\pi\sqrt{n}}\right]$, we obtain that there exists a $\mathbf{w}^{\prime}\in\mathcal{B}\left(\mathbf{w},\frac{3\Delta}{2\pi\sqrt{n}}\right)$ such that

which means that $\sin\left(\left\|\mathbf{\delta}\right\|/2\right)\leq 3\Delta/\left(2\pi\right)$. Because $\sin x\geq\tfrac{3}{\pi}x$ over $x\in\left[0,\pi/6\right]$, it implies that $\left\|\mathbf{\delta}\right\|\leq\Delta$. Since $g(\mathbf{w})=f(\mathbf{q}(\mathbf{w}))$, by summarizing all the results, we conclude that there exists a $\mathbf{\delta}$ with $\left\|\mathbf{\delta}\right\|\leq\Delta$, such that

VI-D Proof of Lemma II.8

Minimizing this function over $t\in\left[0,\frac{3\Delta}{2\pi\sqrt{n}}\right]$, we obtain

and there exists a $\mathbf{w}^{\prime}=\mathbf{w}-t_{\star}\sigma\frac{\mathbf{w}}{\left\|\mathbf{w}\right\|}$ such that

VI-E Proof of Lemma II.10

For any $t\in\left[0,\frac{\Delta}{\left\|\operatorname{grad}f\left(\mathbf{q}^{(r)}\right)\right\|}\right]$, it holds that $\left\|t\;\operatorname{grad}f\left(\mathbf{q}^{(r)}\right)\right\|\leq\Delta$, and the quadratic approximation

Taking $t_{0}=\min\left\{\frac{\Delta}{\left\|\operatorname{grad}f\left(\mathbf{q}^{(r)}\right)\right\|},\frac{1}{M_{H}}\right\}$, we obtain

which means that $\left\|\operatorname{grad}f\left(\mathbf{q}^{(r)}\right)\right\|\geq m_{H}\Delta$. Substituting this into (VI.1), we obtain

By the key comparison result established in proof of Lemma II.6, we have

VI-F Proof of Lemma II.11

It takes certain delicate work to prove Lemma II.11. Basically to use discretization argument, the degree of continuity of the Hessian is needed. The tricky part is that for continuity, we need to compare the Hessian operators at different points, while these Hessian operators are only defined on the respective tangent planes. This is the place where parallel translation comes into play. The next two lemmas compute spectral bounds for the forward and inverse parallel translation operators.

For $\tau\in$ and $\left\|\mathbf{\delta}\right\|\leq 1/2$, we have

where we have used the fact $\sin\left(t\right)\leq t$ and $1-\cos x=2\sin^{2}\left(x/2\right)$. Moreover, $\mathcal{P}_{\gamma}^{0\leftarrow\tau}$ is in the form of $\left(\mathbf{I}+\mathbf{u}\mathbf{v}^{*}\right)^{-1}$ for some vectors $\mathbf{u}$ and $\mathbf{v}$. By the Sherman-Morrison matrix inverse formula, i.e., $\left(\mathbf{I}+\mathbf{u}\mathbf{v}^{*}\right)^{-1}=\mathbf{I}-\mathbf{u}\mathbf{v}^{*}/\left(1+\mathbf{v}^{*}\mathbf{u}\right)$ (justified as $\left\|\left(\cos(\tau\left\|\mathbf{\delta}\right\|)-1\right)\frac{\mathbf{\delta}\mathbf{\delta}^{*}}{\left\|\mathbf{\delta}\right\|^{2}}-\mathbf{q}\sin\left(\tau\left\|\mathbf{\delta}\right\|\right)\frac{\mathbf{\delta}^{*}}{\left\|\mathbf{\delta}\right\|}\right\|\leq 5\tau\left\|\mathbf{\delta}\right\|/4\leq 5/8<1$ as shown above), we have

The next lemma establishes the “local-Lipschitz” property of the Riemannian Hessian.

where $L_{H}=5n^{3/2}\left\|\mathbf{X}_{0}\right\|_{\infty}^{3}/(2\mu^{2})+\frac{9}{\mu}n\left\|\mathbf{X}_{0}\right\|_{\infty}^{2}+9\sqrt{n}\left\|\mathbf{X}_{0}\right\|_{\infty}$.

First of all, by (II.5) and using the fact that the operator norm of a projection operator is unitary bounded, we have

By the estimates in Lemma II.5, we obtain

where at the last line we have used the following estimates:

By Lemma II.5 and substituting the estimate in (VI.5), we obtain the claimed result. ∎

Expectation of the operator. By definition of the Riemannian Hessian in (II.5), expressions of $\nabla^{2}f$ and $\nabla f$ in (II.2) and (II.3), and exchange of differential and expectation operators, we obtain

where we have used $\|\mathbf{w}\|\leq\mu/(4\sqrt{2})$ and $\mu\leq 1/10$ to obtain the last lower bound. Combining the above with the fact that $\|\mathbf{z}\|=1$, we obtain

where we have simplified the expression using $\theta\leq 1/2$. To bound the second term,

where at the last inequality we have applied Gaussian tail upper bound of Type II in Lemma A.2. Since $\left\|\mathbf{w}_{\mathcal{J}}\right\|^{2}+q_{n}^{2}\geq q_{n}^{2}=1-\left\|\mathbf{w}\right\|^{2}\geq 1-\mu^{2}/32\geq 31/32$ for $\left\|\mathbf{w}\right\|\leq\mu/(4\sqrt{2})$ and $\mu\leq 1$, we obtain

Collecting the above estimates, we obtain

where we have used the fact $\mu\leq\theta/10$ to obtain the final lower bound.

Concentration. Next we perform concentration analysis. For any $\mathbf{q}$, we can write

where we have used Lemma B.8 to obtain the last inequality. By Lemma A.4, we obtain

where at the first inequality we used the fact $\left|\tanh(\cdot)\right|\leq 1$, at the second we invoked Lemma B.8, and at the third we invoked Lemma A.3. Taking $R_{Z}=\sigma^{2}_{Z}=1$, by Lemma A.5, we obtain

for any $t>0$. Gathering (VI.9) and (VI.10), we obtain that for any $t>0$,

Uniformizing the bound. Now we are ready to pull above results together for a discretization argument. For any $\varepsilon\in(0,\mu/(4\sqrt{2}))$, there is an $\varepsilon$-net $N_{\varepsilon}$ of size at most $(3\mu/(4\sqrt{2}\varepsilon))^{n}$ that covers the region $\left\{\mathbf{q}:\left\|\mathbf{w}(\mathbf{q})\right\|\leq\mu/(4\sqrt{2})\right\}$. By Lemma VI.2, the function $\operatorname{Hess}f(\mathbf{q})$ is locally Lipschitz within each normal ball of radius

with Lipschitz constant $L_{H}$ (as defined in Lemma VI.2). Note that $\varepsilon<\mu/(4\sqrt{2})<1/(4\sqrt{2})<1/\sqrt{5}$ for $\mu<1$, so any choice of $\varepsilon\in(0,\mu/(4\sqrt{2}))$ makes the Lipschitz constant $L_{H}$ valid within each $\varepsilon$-ball centered around one element of the $\varepsilon$-net. Let

Set $\varepsilon=\frac{\theta}{12\sqrt{2\pi}\mu L_{H}}<\mu/(4\sqrt{2})$, so

Let $\mathcal{E}_{H}$ denote the event that

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

So on $\mathcal{E}_{\infty}\cap\mathcal{E}_{H}$, we have

for any $c_{\sharp}\leq 1/(12\sqrt{2\pi})$. We take $c_{\sharp}=c_{\star}$ for simplicity. Setting $t=\frac{\theta}{12\sqrt{2\pi}\mu}$ in (2), we obtain that for any fixed $\mathbf{q}$ in this region,