Rapid, Robust, and Reliable Blind Deconvolution via Nonconvex Optimization

Introduction

Suppose we are given the convolution of two signals, $\bm{y}=\bm{f}\ast\bm{g}$. When, under which conditions, and how can we reconstruct $\bm{f}$ and $\bm{g}$ from the knowledge of $\bm{y}$ if both $\bm{f}$ and $\bm{g}$ are unknown? This challenging problem, commonly referred to as blind deconvolution problem, arises in many areas of science and technology, including astronomy, medical imaging, optics, and communications engineering, see e.g. . Indeed, the quest for finding a fast and reliable algorithm for blind deconvolution has confounded researchers for many decades.

It is clear that without any additional assumptions, the blind deconvolution problem is ill-posed. One common and useful assumption is to stipulate that $\bm{f}$ and $\bm{g}$ belong to known subspaces . This assumption is reasonable in various applications, provides flexibility and at the same time lends itself to mathematical rigor. We also adopt this subspace assumption in our algorithmic framework (see Section 2 for details). But even with this assumption, blind deconvolution is a very difficult non-convex optimization problem that suffers from an overabundance of local minima, making its numerical solution rather challenging.

In this paper, we present a numerically efficient blind deconvolution algorithm that converges geometrically to the optimal solution. Our regularized gradient descent algorithm comes with rigorous mathematical convergence guarantees. The number of measurements required for the algorithm to succeed is only slightly larger than the information theoretic minimum. Moreover, our algorithm is also robust against noise. To the best of our knowledge, the proposed algorithm is the first blind deconvolution algorithm that is numerically efficient, robust against noise, and comes with rigorous recovery guarantees under certain subspace conditions.

Since blind deconvolution problems are ubiquitous in science and engineering, it is not surprising that there is extensive literature on this topic. It is beyond the scope of this paper to review the existing literature; instead we briefly discuss those results that are closest to our approach. We gladly acknowledge that those papers that are closest to ours, namely , are also the ones that greatly influenced our research in this project.

In the inspiring article , Ahmed, Recht, and Romberg develop a convex optimization framework for blind deconvolution. The formal setup of our blind deconvolution problem follows essentially their setting. Using the meanwhile well-known lifting trick, transforms the blind deconvolution problem into the problem of recovering a rank-one matrix from an underdetermined system of linear equations. By replacing the rank condition by a nuclear norm condition, the computationally infeasible rank minimization problem turns into a convenient convex problem. The authors provide explicit conditions under which the resulting semidefinite program is guaranteed to have the same solution as the original problem. In fact, the number of required measurements is not too far from the theoretical minimum. The only drawback of this otherwise very appealing convex optimization approach is that the computational complexity of solving the semidefinite program is rather high for large-scale data and/or for applications where computation time is of the essence. Overcoming this drawback was one of the main motivations for our paper. While does suggest a fast matrix-factorization based algorithm to solve the semidefinite program, the convergence of this algorithm to the optimal solution is not established in that paper. The theoretical number of measurements required in for the semidefinite program to succeed is essentially comparable to that for our non-convex algorithm to succeed. The advantage of the proposed non-convex algorithm is of course that it is dramatically faster. Furthermore, numerical simulations indicate that the empirically observed number of measurements for our non-convex approach is actually smaller than for the convex approach.

The philosophy underlying the method presented in our paper is strongly motivated by the non-convex optimization algorithm for phase retrieval proposed in , see also . In the pioneering paper the authors use a two-step approach: (i) Construct in a numerically efficient manner a good initial guess; (ii) Based on this initial guess, show that simple gradient descent will converge to the true solution. Our paper follows a similar two-step scheme. At first glance one would assume that many of the proof techniques from should carry over to the blind deconvolution problem. Alas, we quickly found out that despite some general similarities between the two problems, phase retrieval and blind deconvolution are indeed surprisingly different. At the end, we mainly adopted some of the general “proof principles” from (for instance we also have a notion of local regularity condition - although it deviates significantly from the one in ), but the actual proofs are quite different. For instance, in and convergence of the gradient descent algorithm is shown by directly proving that the distance between the true solution and the iterates decreases. The key conditions (a local regularity condition and a local smoothness condition) are tailored to this aim. For the blind deconvolution problem we needed to go a different route. We first show that the objective function decreases during the iterations and then use a certain local restricted isometry property to transfer this decrease to the iterates to establish convergence to the solution.

We also gladly acknowledge being influenced by the papers by Montanari and coauthors on matrix completion via non-convex methods. While the setup and analyzed problems are quite different from ours, their approach informed our strategy in various ways. In , the authors propose an algorithm which comprises a two-step procedure. First, an initial guess is computed via a spectral method, and then a nonconvex problem is formulated and solved via an iterative method. The authors prove convergence to a low-rank solution, but do not establish a rate of convergence. As mentioned before, we also employ a two-step strategy. Moreover, our approach to prove stability of the proposed algorithm draws from ideas in .

We also benefitted tremendously from . In that paper, Sun and Luo devise a non-convex algorithm for low-rank matrix completion and provide theoretical guarantees for convergence to the correct solution. We got the idea of adding a penalty term to the objective function from (as well as from ). Indeed, the particular structure of our penalty term closely resembles that in . In addition, the introduction of the various neighborhood regions around the true solution, that are used to eventually characterize a “basin of attraction”, stems from . These correspondences may not be too surprising, given the connections between low-rank matrix completion and blind deconvolution. Yet, like also discussed in the previous paragraph, despite some obvious similarities between the two problems, it turned out that many steps and tools in our proofs differ significantly from those in . Also this should not come as a surprise, since the measurement matrices and setup differ significantlyAnyone who has experience in the proofs behind compressive sensing and matrix completion is well aware of the substantial challenges one can already face by “simply” changing the sensing matrix from, say, a Gaussian random matrix to one with less randomness.. Moreover, unlike , our paper also provides robustness guarantees for the case of noisy data. Indeed, it seems plausible that some of our techniques to establish robustness against noise are applicable to the analysis of various recent matrix completion algorithms, such as e.g. .

We briefly discuss other interesting papers that are to some extent related to our work. proposes a projected gradient descent algorithm based on matrix factorizations and provide a convergence analysis to recover sparse signals from subsampled convolution. However, this projection step can be hard to implement, which does impact the efficiency and practical use of this method. As suggested in , one can avoid this expensive projection step by resorting to a heuristic approximate projection, but then the global convergence is not fully guaranteed. On the other hand, the papers consider identifiability issue of blind deconvolution problem with both $\bm{f}$ and $\bm{g}$ in random linear subspaces and achieve nearly optimal result of sampling complexity in terms of information theoretic limits. Very recently, improved the result from by using techniques from algebraic geometry.

The past few years have also witnessed an increasing number of excellent works other than blind deconvolution but related to nonconvex optimization . The paper analyzes the problem of recovering a low-rank positive semidefinite matrix from linear measurements via a gradient descent algorithm. The authors assume that the measurement matrix fulfills the standard and convenient restricted isometry property, a condition that is not suitable for the blind deconvolution problem (besides the fact that the positive semidefinite assumption is not satisfied in our setting). In , Chen and Wainwright study various the solution of low-rank estimation problems by projected gradient descent. The very recent paper investigates matrix completion for rectangular matrices. By “lifting”, they convert the unknown matrix into a positive semidefinite one and apply matrix factorization combined with gradient descent to reconstruct the unknown entries of the matrix. considers an interesting blind calibration problem with a special type of measurement matrix via nonconvex optimization. Besides some general similarities, there is little overlap of the aforementioned papers with our framework. Finally, a convex optimization approach to blind deconvolution and self-calibration that extends the work of can be found in , while also covers the joint blind deconvolution-blind demixing problem.

2 Organization of our paper

This paper is organized as follows. We introduce some notation used throughout the paper in the remainder of this section. The model setup and problem formulation are presented in Section 2. Section 3 describes the proposed algorithm and our main theoretical result establishing the convergence of our algorithm. Numerical simulations can be found in Section 4. Section 5 is devoted to the proof of the main theorem. Since the proof is quite involved, we have split this section into several subsections. Some auxiliary results are collected in the Appendix.

3 Notation

We introduce notation which will be used throughout the paper. Matrices and vectors are denoted in boldface such as $\bm{Z}$ and $\bm{z}$. The individual entries of a matrix or a vector are denoted in normal font such as $Z_{ij}$ or $z_{i}.$ For any matrix $\bm{Z}$, $\|\bm{Z}\|_{*}$ denotes its nuclear norm, i.e., the sum of its singular values; $\|\bm{Z}\|$ denotes its operator norm, i.e., the largest singular value, and $\|\bm{Z}\|_{F}$ denotes its the Frobenius norm, i.e., $\|\bm{Z}\|_{F}=\sqrt{\sum_{ij}|Z_{ij}|^{2}}$. For any vector $\bm{z}$, $\|\bm{z}\|$ denotes its Euclidean norm. For both matrices and vectors, $\bm{Z}^{\top}$ and $\bm{z}^{\top}$ stand for the transpose of $\bm{Z}$ and $\bm{z}$ respectively while $\bm{Z}^{*}$ and $\bm{z}^{*}$ denote their complex conjugate transpose. We equip the matrix space $\hbox{\msbm{C}}^{K\times N}$ with the inner product defined as $\left\langle\bm{U},\bm{V}\right\rangle:=\text{Tr}(\bm{U}^{*}\bm{V}).$ A special case is the inner product of two vectors, i.e., $\left\langle\bm{u},\bm{v}\right\rangle=\text{Tr}(\bm{u}^{*}\bm{v})=\bm{u}^{*}\bm{v}.$ For a given vector $\bm{v}$, $\operatorname{diag}(\bm{v})$ represents the diagonal matrix whose diagonal entries are given by the vector $\bm{v}$. For any $z\in\hbox{\msbm{R}}$, denote $z_{+}$ as $z_{+}=\frac{z+|z|}{2}.$ $C$ is an absolute constant and $C_{\gamma}$ is a constant which depends linearly on $\gamma$, but on no other parameters.

Problem setup

We consider the blind deconvolution model

where $\bm{y}$ is given, but $\bm{f}$ and $\bm{g}$ are unknown. Here “$\ast$” denotes circular convolutionAs usual, ordinary convolution can be well approximated by circulant convolution, as long as the function $\bm{f}$ decays sufficiently fast .. We will usually consider $\bm{f}$ as the “blurring function” and $\bm{g}$ as the signal of interest. It is clear that without any further assumption it is impossible to recover $\bm{f}$ and $\bm{g}$ from $\bm{y}$. We want to impose conditions on $\bm{f}$ and $\bm{g}$ that are realistic, flexible, and not tied to one particular application (such as, say, image deblurring). At the same time, these conditions should be concrete enough to lend themselves to meaningful mathematical analysis.

A natural setup that fits these demands is to assume that $\bm{f}$ and $\bm{g}$ belong to known linear subspaces. Concerning the blurring function, it is reasonable in many applications to assume that $\bm{f}$ is either compactly supported or that $\bm{f}$ decays sufficiently fast so that it can be well approximated by a compactly supported function. Therefore, we assume that $\bm{f}\in\hbox{\msbm{C}}^{L}$ satisfies

For our theoretical analysis as well as for numerical purposes, it is much more convenient to express (2.1) in the Fourier domain, see also . To that end, let $\bm{F}$ be the $L\times L$ unitary Discrete Fourier Transform (DFT) matrix and let the $L\times K$ matrix $\bm{B}$ be given by the first $K$ columns of $\bm{F}$ (then $\bm{B}^{*}\bm{B}=\bm{I}_{K}$ ). By applying the scaled DFT matrix $\sqrt{L}\bm{F}$ to both sides of (2.1) we get

which follows from the property of circular convolution and Discrete Fourier Transform. Here, $\operatorname{diag}(\sqrt{L}\bm{F}\bm{f})(\sqrt{L}\bm{F}\bm{g})=(\sqrt{L}\bm{F}\bm{f})\odot(\sqrt{L}\bm{F}\bm{g})$ where $\odot$ denotes pointwise product. By definition of $\bm{f}$ in (2.2), we have

where $\overline{\bm{A}}=\bm{F}\bm{C}$ (we use $\overline{\bm{A}}$ instead of $\bm{A}$ simply because it gives rise to a more convenient notation later, see e.g. (2.8)).

For the remainder of the paper, instead of considering the original blind deconvolution problem (2.1), we focus on its mathematically equivalent version (2.5), where $\bm{h}\in\hbox{\msbm{C}}^{K\times 1}$, $\bm{x}\in\hbox{\msbm{C}}^{N\times 1}$, $\bm{y}\in\hbox{\msbm{C}}^{L\times 1}$, $\bm{B}\in\hbox{\msbm{C}}^{L\times K}$ and $\bm{A}\in\hbox{\msbm{C}}^{L\times N}$. As mentioned before, $\bm{h}_{0}$ and $\bm{x}_{0}$ are the ground truth. Our goal is to recover $\bm{h}_{0}$ and $\bm{x}_{0}$ when $\bm{B}$, $\bm{A}$ and $\bm{y}$ are given. It is clear that if $(\bm{h}_{0},\bm{x}_{0})$ is a solution to the blind deconvolution problem. then so is $(\alpha\bm{h}_{0},\alpha^{-1}\bm{x}_{0})$ for any $\alpha\neq 0$. Thus, all we can hope for in the absence of any further information, is to recover a solution from the equivalence class $(\alpha\bm{h}_{0},\alpha^{-1}\bm{x}_{0}),\alpha\neq 0$. Hence, we can as well assume that $\|\bm{h}_{0}\|=\|\bm{x}_{0}\|:=\sqrt{d}_{0}$.

As already mentioned, we choose $\bm{B}$ to be the “low-frequency” discrete Fourier matrix, i.e., the first $K$ columns of an $L\times L$ unitary DFT (Discrete Fourier Transform) matrix. Moreover, we choose $\bm{A}$ to be an $L\times N$ complex Gaussian random matrix, i.e.,

We define the matrix-valued linear operator $\mathcal{A}(\cdot)$ via

where $\bm{b}_{l}$ denotes the $l$-th column of $\bm{B}^{*}$ and $\bm{a}_{l}$ is the $l$-th column of $\bm{A}^{*}.$ Immediately, we have $\sum_{l=1}^{L}\bm{b}_{l}\bm{b}_{l}^{*}=\bm{B}^{*}\bm{B}=\bm{I}_{K}$, $\|\bm{b}_{l}\|=\frac{K}{L}$ and $\operatorname{E}(\bm{a}_{l}\bm{a}_{l}^{*})=\bm{I}_{N}$ for all $1\leq l\leq L$. This is essentially a smart and popular trick called “lifting” , which is able to convert a class of nonlinear models into linear models at the costs of increasing the dimension of the solution space.

It seems natural and tempting to recover $(\bm{h}_{0},\bm{x}_{0})$ obeying (2.5) by solving the following optimization problem

which is a special case of $F(\bm{h},\bm{x})$ when $\bm{e}=\bm{0}.$ Furthermore, we define $\delta(\bm{h},\bm{x})$, an important quantity throughout our discussion, via

When there is no danger of ambiguity, we will often denote $\delta(\bm{h},\bm{x})$ simply by $\delta$. But let us remember that $\delta(\bm{h},\bm{x})$ is always a function of $(\bm{h},\bm{x})$ and measures the relative approximation error of $(\bm{h},\bm{x})$.

Obviously, minimizing (2.8) becomes a nonlinear least square problem, i.e., one wants to find a pair of vectors $(\bm{h},\bm{x})$ or a rank-1 matrix $\bm{h}\bm{x}^{*}$ which fits the measurement equation in (2.8) best. Solving (2.7) is a challenging optimization problem since it is highly nonconvex and most of the available algorithms, such as alternating minimization and gradient descent, may suffer from getting easily trapped in some local minima. Another possibility is to consider a convex relaxation of (2.7) at the cost of having to solve an expensive semidefinite program. In the next section we will describe how to avoid this dilemma and design an efficient gradient descent algorithm that, under reasonable conditions, will always converge to the true solution.

Algorithm and main result

In this section we introduce our algorithm as well as our main theorems which establish convergence of the proposed algorithm to the true solution. As mentioned above, in a nutshell our algorithm consists of two parts: First we use a carefully chosen initial guess, and second we use a variation of gradient descent, starting at the initial guess to converge to the true solution. One of the most critical aspects is of course that we must avoid getting stuck in local minimum or saddle point. Hence, we need to ensure that our iterates are inside some properly chosen basin of attraction of the true solution. The appropriate characterization of such a basin of attraction requires some diligence, a task that will occupy us in the next subsection. We will then proceed to introducing our algorithm and analyzing its convergence.

The road toward designing a proper basin of attraction is basically paved by three observations, described below. These observations prompt us to introduce three neighborhoods (inspired by ), whose intersection will form the desired basin of attraction of the solution.

Observation 1 - Nonuniqueness of the solution: As pointed out earlier, if $(\bm{h},\bm{x})$ is a solution to (2.5), then so is $(\alpha\bm{h},\alpha^{-1}\bm{x})$ for any $\alpha\neq 0$. Thus, without any prior information about $\|\bm{h}\|$ and/or $\|\bm{x}\|$, it is clear that we can only recover the true solution up to such an unknown constant $\alpha$. Fortunately, this suffices for most applications. From the viewpoint of numerical stability however, we do want to avoid, while $\|\bm{h}\|\|\bm{x}\|$ remains bounded, that $\|\bm{h}\|\to 0$ and $\|\bm{x}\|\to\infty$ (or vice versa). To that end we introduce the following neighborhood:

(Recall that $d_{0}=\|\bm{h}_{0}\|\|\bm{x}_{0}\|$.)

Observation 2 - Incoherence: Our numerical simulations indicate that the number of measurements required for solving the blind deconvolution problem with the proposed algorithm does depend (among others) on how much $\bm{h}_{0}$ is correlated with the rows of the matrix $\bm{B}$ — the smaller the correlation the better. A similar effect has been observed in blind deconvolution via convex programming . We quantify this property by defining the incoherence between the rows of $\bm{B}$ and $\bm{h}_{0}$ via

It is easy to see that $1\leq\mu_{h}^{2}\leq K$ and both lower and upper bounds are tight; i.e., $\mu^{2}_{h}=K$ if $\bm{h}_{0}$ is parallel to one of $\{\bm{b}_{l}\}_{l=1}^{L}$ and $\mu^{2}_{h}=1$ if $\bm{h}_{0}$ is a 1-sparse vector of length $K$. Note that in our setup, we assume that $\bm{A}$ is a random matrix and $\bm{x}_{0}$ is fixed, thus with high probability, $\bm{x}_{0}$ is already sufficiently incoherent with the rows of $\bm{A}$ and thus we only need to worry about the incoherence between $\bm{B}$ and $\bm{h}_{0}$.

It should not come as a complete surprise that the incoherence between $\bm{h}_{0}$ and the rows of $\bm{B}$ is important. The reader may recall that in matrix completion the left and right singular vectors of the solution cannot be “too aligned” with those of the measurement matrices. A similar philosophy seems to apply here. Being able to control the incoherence of the solution is instrumental in deriving rigorous convergence guarantees of our algorithm. For that reason, we introduce the neighborhood

Observation 3 - Initial guess: It is clear that due to the non-convexity of the objective function, we need a carefully chosen initial guess. We quantify the distance to the true solution via the following neighborhood

where $\varepsilon$ is a predetermined parameter in $(0,\frac{1}{15}]$.

It is evident that the true solution $(\bm{h}_{0},\bm{x}_{0})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}$. Note that $(\bm{h},\bm{x})\in\mathcal{N}_{d_{0}}\bigcap\mathcal{N}_{\varepsilon}$ implies $\|\bm{h}\bm{x}^{*}\|\geq(1-\varepsilon)d_{0}$ and $\frac{1}{\|\bm{h}\|}\leq\frac{\|\bm{x}\|}{(1-\varepsilon)d_{0}}\leq\frac{2}{(1-\varepsilon)\sqrt{d_{0}}}$. Therefore, for any element $(\bm{h},\bm{x})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon},$ its incoherence can be well controlled by

2 Objective function and key ideas of the algorithm

Our approach consists of two parts: We first construct an initial guess that is inside the “basin of attraction” $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$. We then apply a carefully regularized gradient descent algorithm that will ensure that all the iterates remain inside $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$.

Due to the difficulties of directly projecting onto $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}$ (the neigbourhood $\mathcal{N}_{\varepsilon}$ is easier to manage) we add instead a regularizer $G(\bm{h},\bm{x})$ to the objective function $F(\bm{h},\bm{x})$ to enforce that the iterates remain inside $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}$. While the idea of adding a penalty function to control incoherence is proposed in different forms to solve matrix completion problems, see e.g., , our version is mainly inspired by .

Hence, we aim to minimize the following regularized objective function to solve the blind deconvolution problem:

where $F(\bm{h},\bm{x})$ is defined in (2.8) and $G(\bm{h},\bm{x})$, the penalty function, is of the form

where $G_{0}(z)=\max\{z-1,0\}^{2}$ and $\rho\geq d^{2}+2\|\bm{e}\|^{2}$. Here we assume $\frac{9}{10}d_{0}\leq d\leq\frac{11}{10}d_{0}$ and $\mu\geq\mu_{h}$.

The idea behind this, at first glance complicated, penalty function is quite simple. The first two terms in (3.6) enforce the projection of $(\bm{h},\bm{x})$ onto $\mathcal{N}_{d_{0}}$ while the last term is related to $\mathcal{N}_{\mu}$; it will be shown later that any $(\bm{h},\bm{x})\in\frac{1}{\sqrt{3}}\mathcal{N}_{d_{0}}\bigcap\frac{1}{\sqrt{3}}\mathcal{N}_{\mu}$ gives $G(\bm{h},\bm{x})=0$ if $\frac{9}{10}d_{0}\leq d\leq\frac{11}{10}d_{0}$. Since $G_{0}(z)$ is a truncated quadratic function, it is obvious that $G_{0}^{\prime}(z)=2\sqrt{G_{0}(z)}$ and $G(\bm{h},\bm{x})$ is a continuously differentiable function. Those two properties play a crucial role in proving geometric convergence of our algorithm presented later.

3 Wirtinger derivative of the objective function and algorithm

In particular, we denote $\nabla\widetilde{F}_{\bm{h}}:=\frac{\partial\widetilde{F}}{\partial\bar{\bm{h}}}$ and $\nabla\widetilde{F}_{\bm{x}}:=\frac{\partial\widetilde{F}}{\partial\bar{\bm{x}}}$.

We also introduce the adjoint operator of $\mathcal{A}:\hbox{\msbm{C}}^{L}\rightarrow\hbox{\msbm{C}}^{K\times N}$, given by

Both $\nabla\widetilde{F}_{\bm{h}}$ and $\nabla\widetilde{F}_{\bm{x}}$ can now be expressed as

Our algorithm consists of two steps: initialization and gradient descent with constant stepsize. The initialization is achieved via a spectral method followed by projection. The idea behind spectral method is that

and hence one can hope that the leading singular value and vectors of $\mathcal{A}^{*}(\bm{y})$ can be a good approximation of $d_{0}$ and $(\bm{h}_{0},\bm{x}_{0})$ respectively. The projection step ensures $\bm{u}_{0}\in\mathcal{N}_{\mu}$, which the spectral method alone might not guarantee. We will address the implementation and computational complexity issue in Section 4.

4 Main results

Our main finding is that with a diligently chosen initial guess $(\bm{u}_{0},\bm{v}_{0})$, simply running gradient descent to minimize the regularized non-convex objective function $\widetilde{F}(\bm{h},\bm{x})$ will not only guarantee linear convergence of the sequence $(\bm{u}_{t},\bm{v}_{t})$ to the global minimum $(\bm{h}_{0},\bm{x}_{0})$ in the noiseless case, but also provide robust recovery in the presence of noise. The results are summarized in the following two theorems.

The initialization obtained via Algorithm 1 satisfies

holds with probability at least $1-L^{-\gamma}$ if the number of measurements satisfies

Here $\varepsilon$ is any predetermined constant in $(0,\frac{1}{15}]$, and $C_{\gamma}$ is a constant only linearly depending on $\gamma$ with $\gamma\geq 1$.

The proof of Theorem 3.1 is given in Section 5.5. While the initial guess is carefully chosen, it is in general not of sufficient accuracy to already be used as good approximation to the true solution. The following theorem establishes that as long as the initial guess lies inside the basin of attraction of the true solution, regularized gradient descent will indeed converge to this solution (or to a solution nearby in case of noisy data).

Algorithm 2 will create a sequence $(\bm{u}_{t},\bm{v}_{t})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$ which converges geometrically to $(\bm{h}_{0},\bm{x}_{0})$ in the sense that with probability at least $1-4L^{-\gamma}-\frac{1}{\gamma}\exp(-(K+N))$, there holds

where $d_{t}:=\|\bm{u}_{t}\|\|\bm{v}_{t}\|$, $\omega>0$, $\eta$ is the fixed stepsize and $\angle(\bm{u}_{t},\bm{h}_{0})$ is the angle between $\bm{u}_{t}$ and $\bm{h}_{0}$. Here

While the setup in (2.5) assumes that $\bm{B}$ is a matrix consisting of the first $K$ columns of the DFT matrix, this is actually not necessary for Theorem 3.2. As the proof will show, the only conditions on $\bm{B}$ are that $\bm{B}^{\ast}\bm{B}=\bm{I}_{K}$ and that the norm of the $l$-th row of $\bm{B}$ satisfies $\|\bm{b}_{l}\|^{2}\leq C\frac{K}{L}$ for some numerical constant $C$.

The minimum number of measurements required for our method to succeed is roughly comparable to that of the convex approach proposed in (up to log-factors). Thus there is no price to be paid for trading a slow, convex-optimization based approach with a fast non-convex based approach. Indeed, numerical experiments indicate that the non-convex approach even requires a smaller number of measurements compared to the convex approach, see Section 4.

The convergence rate of our algorithm is completely determined by $\eta\omega$. Here, the regularity constant $\omega=\mathcal{O}(d_{0})$ is specified in (5.7) and $\eta\leq\frac{1}{C_{L}}$ where $C_{L}=\mathcal{O}(d_{0}(N\log L+\frac{\rho L}{d_{0}^{2}\mu^{2}}))$. The attentive reader may have noted that $C_{L}$ depends essentially linearly on $\frac{\rho L}{\mu^{2}}$, which actually reflects a tradeoff between sampling complexity (or statistical estimation quality) and computation time. Note that if $L$ gets larger, the number of constraints is also increasing and hence leads to a larger $C_{L}$. However, this issue can be solved by choosing parameters smartly. Theorem 3.2 tells us that $\mu^{2}$ should be roughly between $\mu_{h}^{2}$ and $\frac{L}{(K+N)\log^{2}L}$. Therefore, by choosing $\mu^{2}=\mathcal{O}(\frac{L}{(K+N)\log^{2}L})$ and $\rho\approx d^{2}+2\|\bm{e}\|^{2}$, $C_{L}$ is optimized and

which is shown in details in Section 5.4.

Relations (3.15) and (3.16) are basically equivalent to the following:

which says that $(\bm{u}_{t},\bm{v}_{t})$ converges to an element of the equivalence class associated with the true solution $(\bm{h}_{0},\bm{x}_{0})$ (up to a deviation governed by the amount of additive noise).

The matrix $\mathcal{A}^{*}(\bm{e})=\sum_{l=1}^{L}e_{l}\bm{b}_{l}\bm{a}_{l}^{*}$, as a sum of $L$ rank-1 random matrices, has nice concentration of measure properties under the assumption of Theorem 3.2. Asymptotically, $\|\mathcal{A}^{*}(\bm{e})\|$ converges to with rate $\mathcal{O}(L^{-1/2})$, which will be justified in Lemma 5.20 of Section 5.5 (see also ). Note that

If one lets $L\rightarrow\infty$, then $\|\bm{e}\|^{2}\sim\frac{\sigma^{2}d_{0}^{2}}{2L}\chi^{2}_{2L}$ will converge almost surely to $\sigma^{2}d_{0}^{2}$ under the Law of Large Numbers and the cross term $\operatorname{Re}(\left\langle\bm{h}\bm{x}^{*}-\bm{h}_{0}\bm{x}_{0}^{*},\mathcal{A}^{*}(\bm{e})\right\rangle)$ will converge to . In other words, asymptotically,

for all fixed $(\bm{h},\bm{x})$. This implies that if the number of measurements is large, then $F(\bm{h},\bm{x})$ behaves “almost like” $F_{0}(\bm{h},\bm{x})=\|\mathcal{A}(\bm{h}\bm{x}^{*}-\bm{h}_{0}\bm{x}_{0}^{*})\|^{2}$, the noiseless version of $F(\bm{h},\bm{x})$. This provides the key insight into analyzing the robustness of our algorithm, which is reflected in the so-called “Robustness Condition” in (5.2). Moreover, the asymptotic property of $\mathcal{A}^{*}(\bm{e})$ is also seen in our main result (3.18). Suppose $L$ is becoming larger and larger, the effect of noise diminishes and heuristically, we might just rewrite our result as $\|\bm{u}_{t}\bm{v}_{t}^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}\leq\frac{2}{3}(1-\eta\omega)^{t/2}\varepsilon d_{0}+o_{p}(1)$, which is consistent with the result without noise.

Numerical simulations

We present empirical evaluation of our proposed gradient descent algorithm (Algorithm 2) using simulated data as well as examples from blind deconvolution problems appearing in communications and in image processing.

We first investigate how many measurements are necessary in order for an algorithm to reliably recover two signals from their convolution. We compare Algorithm 2, a gradient descent algorithm for the sum of the loss function $F(\bm{h},\bm{x})$ and the regularization term $G(\bm{h},\bm{x})$, with the gradient descent algorithm only applied to $F(\bm{h},\bm{x})$ and the nuclear norm minimization proposed in . These three tested algorithms are abbreviated as regGrad, Grad and NNM respectively. To make fair comparisons, both regGrad and Grad are initialized with the normalized leading singular vectors of $\mathcal{A}^{*}(\bm{y})$, which are computed by running the power method for $50$ iterations. Though we do not further compute the projection of $\hat{\bm{h}}_{0}$ for regGrad as stated in the third step of Algorithm 1, we emphasize that the projection can be computed efficiently as it is a linear programming on $K$-dimensional vectors. A careful reader may notice that in addition to the computational cost for the loss function $F(\bm{h},\bm{x})$, regGrad also requires to evaluate $G(\bm{h},\bm{x})$ and its gradient in each iteration. When $B$ consists of the first $K$ columns of a unitary DFT matrix, we can evaluate $\left\{\bm{b}_{l}^{*}\bm{h}\right\}_{l=1}^{L}$ and the gradient of $\sum_{l=1}^{L}G_{0}\left(\frac{L|\bm{b}_{l}^{*}\bm{h}|^{2}}{8d\mu^{2}}\right)$ using FFT. Thus the additional per iteration computational cost for the regularization term $G(\bm{h},\bm{x})$ is only $O(L\log L)$ flops. The stepsizes in both regGrad and Grad are selected adaptively in each iteration via backtracking. As suggested by the theory, the choices for $\rho$ and $\mu$ in regGrad are $\rho=d^{2}/100$ and $\mu=6\sqrt{L/(K+N)}/\log L$.

We conduct tests on random Gaussian signals $\bm{h}\in\hbox{\msbm{C}}^{K\times 1}$ and $\bm{x}\in\hbox{\msbm{C}}^{N\times 1}$ with $K=N=50$. The matrix $\bm{B}\in\hbox{\msbm{C}}^{L\times K}$ is the first $K$ columns of a unitary $L\times L$ DFT matrix, while $\bm{A}\in\hbox{\msbm{C}}^{L\times N}$ is either a Gaussian random matrix or a partial Hadamard matrix with randomly selected $N$ columns and then multiplied by a random sign matrix from the left. When $\bm{A}$ is a Gaussian random matrix, $L$ takes $16$ equal spaced values from $K+N$ to $4(K+N)$. When $\bm{A}$ is a partial Hadamard matrix, we only test $L=2^{s}$ with $6\leq s\leq 10$ being integers. For each given triple $(K,N,L)$, fifty random tests are conducted. We consider an algorithm to have successfully recovered $(\bm{h}_{0},\bm{x}_{0})$ if it returns a matrix $\hat{\bm{X}}$ which satisfies

We present the probability of successful recovery plotted against the number of measurements in Figure 1.

𝐾𝑁L/(K+N) and vertical axis probability of successful recovery out of $50$ random tests. It can be observed that regGrad and Grad have similar performance, and both of them require a significantly smaller number of measurements than NNM to achieve successful recovery of high probability.

2 Number of measurements vs incoherence

Theorem 3.2 indicates that the number of measurements $L$ required for Algorithm 2 to achieve successful recovery scales linearly with $\mu_{h}^{2}$. We conduct numerical experiments to investigate the dependence of $L$ on $\mu_{h}^{2}$ empirically. The tests are conducted with $\mu_{h}^{2}$ taking on $10$ values $\mu_{h}^{2}\in\{3,6,\cdots,30\}$. For each fixed $\mu_{h}^{2}$, we choose $\bm{h}_{0}$ to be a vector whose first $\mu_{h}^{2}$ entries are $1$ and the others are so that its incoherence is equal to $\mu_{h}^{2}$ when $\bm{B}$ is low frequency Fourier matrix. Then the tests are repeated for random Gaussian matrices $\bm{A}$ and random Gaussian vectors $\bm{x}_{0}$. The empirical probability of successful recovery on the $(\mu_{h}^{2},L)$ plane is presented in Figure 2, which suggests that $L$ does scale linearly with $\mu_{h}^{2}$.

While regGrad and Grad have similar performances in the simulation when $\bm{h}_{0}$ is a random Gaussian signal (see Figure 1), we investigate their performances on a fixed $\bm{h}_{0}$ with a large incoherence. The tests are conducted for $K=N=200$, $\bm{x}_{0}\in\hbox{\msbm{C}}^{N\times 1}$ being a random Gaussian signal, $\bm{A}\in\hbox{\msbm{C}}^{L\times N}$ being a random Gaussian matrix, and $\bm{B}\in\hbox{\msbm{C}}^{L\times K}$ being a low frequency Fourier matrix. The signal $\bm{h}_{0}$ with $\mu_{h}^{2}=100$ is formed in the same way as in Section 4.2; that is, the first $100$ entries of $\bm{h}_{0}$ are one and the other entries are zero. The number of measurements $L$ varies from $3(K+N)$ to $8(K+N)$. For each $L$, $100$ random tests are conducted. Figure 3 shows the probability of successful recovery for regGrad and Grad. It can be observed that the successful recovery probability of regGrad is at least $10\%$ larger than that of Grad when $L\geq 6(K+N)$.

4 Robustness to additive noise

We explore the robustness of Algorithm 2 when the measurements are contaminated by additive noise. The tests are conducted with $K=N=100$, $L\in\{500,1000\}$ when $\bm{A}$ is a random Gaussian matrix, and $L\in\{512,1024\}$ when $\bm{A}$ is a partial Hadamard matrix. Tests with additive noise have the measurement vector $\bm{y}$ corrupted by the vector

where $\bm{w}\in\hbox{\msbm{C}}^{L\times 1}$ is standard Gaussian random vector, and $\sigma$ takes nine different values from $10^{-4}$ to $1$. For each $\sigma$, fifty random tests are conducted. The average reconstruction error in dB plotted against the signal to noise ratio (SNR) is presented in Fig. 4. First the plots clearly show the desirable linear scaling between the noise levels and the relative reconstruction errors. Moreover, as desired, the relative reconstruction error decreases linearly on a $\log$-$\log$ scale as the number of measurements $L$ increases.

5 An example from communications

In order to demonstrate the effectiveness of Algorithm 2 for real world applications, we first test the algorithm on a blind deconvolution problem arising in communications. Indeed, blind deconvolution problems and their efficient numerical solution are expected to play an increasingly important role in connection with the emerging Internet-of-Things . Assume we want to transmit a signal from one place to another over a so-called time-invariant multi-path communication channel, but the receiver has no information about the channel, except its delay spread (i.e., the support of the impulse response). In many communication settings it is reasonable to assume that the receiver has information about the signal encoding matrix—in other words, we know the subspace $\bm{A}$ to which $\bm{x}_{0}$ belongs to. Translating this communications jargon into mathematical terminology, this simply means that we are dealing with a blind deconvolution problem of the form (2.5).

The plots of successful recovery probability are presented in Figure 5, which shows that Algorithm 2 can successfully recover the real channel $\bm{h}_{0}$ with high probability if $L\gtrsim 2.5(L+N)$ when $\bm{A}$ is a random Gaussian matrix and if $L\gtrsim 4.5(L+N)$ when $\bm{A}$ is a partial Hadamard matrix. Thus our theory seems a bit pessimistic. It is gratifying to see that very little additional measurements are required compared to the number of unknowns in order to recover the transmitted signal.

6 An example from image processing

Next, we test Algorithm 2 on an image deblurring problem, inspired by . The observed image (Figure 6c) is a convolution of a $512\times 512$ MRI image (Figure 6a) with a motion blurring kernel (Figure 6b). Since the MRI image is approximately sparse in the Haar wavelet basis, we can assume it belongs to a low dimensional subspace; that is, $\bm{g}=\bm{C}\bm{x}_{0}$, where $\bm{g}\in\hbox{\msbm{C}}^{L}$ with $L=262,144$ denotes the MRI image reshaped into a vector, $\bm{C}\in\hbox{\msbm{C}}^{L\times N}$ represents the wavelet subspace and $\bm{x}_{0}\in\hbox{\msbm{C}}^{N}$ is the vector of wavelet coefficients. The blurring kernel $\bm{f}\in\hbox{\msbm{C}}^{L\times 1}$ is supported on a low frequency region. Therefore $\widehat{\bm{f}}=\bm{B}\bm{h}_{0}$, where $\bm{B}\in\hbox{\msbm{C}}^{L\times K}$ is a reshaped $2$D low frequency Fourier matrix and $\bm{h}_{0}\in\hbox{\msbm{C}}^{K}$ is a short vector.

Figure 6d shows the initial guess for Algorithm 2 in the image domain, which is obtained by running the power method for fifty iterations. While this initial guess is clearly not a good approximation to the true solution, it suffices as a starting point for gradient descent. In the first experiment, we take $\bm{C}$ to be the wavelet subspace corresponding to the $N=20000$ largest Haar wavelet coefficients of the original MRI image, and we also assume the locations of the $K=65$ nonzero entries of the kernel are known. Figure 6e shows the reconstructed image in this ideal setting. It can be observed that the recovered image is visually indistinguishable from the ground truth MRI image. In the second experiment, we test a more realistic setting, where both the support of the MRI image is the wavelet domain and the support of the kernel are not known. We take the Haar wavelet transform of the blurred image (Figure 6c) and select $\bm{C}$ to be the wavelet subspace corresponding to the $N=35000$ largest wavelet coefficients. We do not assume the exact support of the kernel is known, but assume that its support is contained in a small box region. The reconstructed image in this setting is shown in Figure 6f. Despite not knowing the subspaces exactly, Algorithm 2 is still able to return a reasonable reconstruction.

Yet, this second experiment also demonstrates that there is clearly room for improvement in the case when the subspaces are unknown. One natural idea to improve upon the result depicted in Figure 6f is to include an additional total-variation penalty in the reconstruction algorithm. We leave this line of work for future research.

Proof of the main theorem

This section is devoted to the proof of Theorems 3.1 and 3.2. Since proving Theorem 3.2 is a bit more involved, we briefly describe the architecture of its proof. In Subsection 5.1 we state four key conditions: The Local Regularity Condition will allow us to show that the objective function decreases; the Local Restricted Isometry Property enables us to transfer the decrease in the objective function to a decrease of the error between the iterates and the true solution; the Local Smoothness Condition yields the actual rate of convergence, and finally, the Robustness Condition establishes robustness of the proposed algorithm against additive noise. Armed with these conditions, we will show how the three regions defined in Section 3 characterize the convergence neighborhood of the solution, i.e., if the initial guess is inside this neighborhood, the sequence generated via gradient descent will always stay inside this neighborhood as well. In Subsections 5.2–5.5 we justify the aforementioned four conditions, show that they are valid under the assumptions stated in Theorem 3.2, and conclude with a proof of Theorem 3.1.

The following local Restricted Isometry Property (RIP) for $\mathcal{A}$ holds uniformly for all $(\bm{h},\bm{x})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}:$

if $L\geq C_{\gamma}(\frac{\sigma^{2}}{\varepsilon^{2}}+\frac{\sigma}{\varepsilon})\max\{K,N\}\log L$.

This condition follows directly from (3.17). It is quite essential when we analyze the behavior of Algorithm 2 under Gaussian noise. With those two conditions above in hand, the lower and upper bounds of $F(\bm{h},\bm{x})$ are well approximated over $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$ by two quadratic functions of $\delta$, where $\delta$ is defined in (2.10). A similar approach towards noisy matrix completion problem can be found in . For any $(\bm{h},\bm{x})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$, applying Condition 5.1 to (3.19) leads to

because $\|\cdot\|$ and $\|\cdot\|_{*}$ is a pair of dual norm and $\operatorname{\text{rank}}(\bm{h}\bm{x}^{*}-\bm{h}_{0}\bm{x}_{0}^{*})\leq 2$. Moreover, with the Condition 5.2, (5.3) and (5.4) yield the followings:

The third condition is about the regularity condition of $\widetilde{F}(\bm{h},\bm{x})$, which is the key to establishing linear convergence later. The proof will be given in Lemma 5.18.

Let $\widetilde{F}(\bm{h},\bm{x})$ be as defined in (3.5) and $\nabla\widetilde{F}(\bm{h},\bm{x}):=(\nabla\widetilde{F}_{\bm{h}},\nabla\widetilde{F}_{\bm{x}})\in\hbox{\msbm{C}}^{K+N}$. Then there exists a regularity constant $\omega=\frac{d_{0}}{5000}>0$ such that

for all $(\bm{h},\bm{x})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$ where $c=\|\bm{e}\|^{2}+a\|\mathcal{A}^{*}(\bm{e})\|^{2}$ with $a=1700$. In particular, in the noiseless case, i.e., $\bm{e}=\bm{0}$, we have

Besides the three regions defined in (3.1) to (3.3), we define another region $\mathcal{N}_{\widetilde{F}}$ via

for proof technical purposes. $\mathcal{N}_{\widetilde{F}}$ is actually the sublevel set of the nonconvex function $\widetilde{F}$.

Finally we introduce the last condition called Local smoothness condition and its corresponding quantity $C_{L}$ which characterizes the choice of stepsize $\eta$ and the rate of linear convergence.

Denote $\bm{z}:=(\bm{h},\bm{x})$. There exists a constant $C_{L}$ such that

for all $\{(\bm{z},\Delta\bm{z})|\bm{z}+t\Delta\bm{z}\in\mathcal{N}_{\varepsilon}\bigcap\mathcal{N}_{\widetilde{F}},\forall 0\leq t\leq 1\}$, i.e., the whole segment connecting $\bm{z}$ and $\bm{z}+\Delta\bm{z}$, which is parametrized by $t$, belongs to the nonconvex set $\mathcal{N}_{\varepsilon}\bigcap\mathcal{N}_{\widetilde{F}}.$

The upper bound of $C_{L}$, which scales with $\mathcal{O}(d_{0}(1+\sigma^{2})(K+N)\log^{2}L)$, will be given in Section 5.4. We will show later in Lemma 5.8 that the stepsize $\eta$ is chosen to be smaller than $\frac{1}{C_{L}}.$ Hence $\eta=\mathcal{O}((d_{0}(1+\sigma^{2})(K+N)\log^{2}L)^{-1}).$

There holds $\mathcal{N}_{\widetilde{F}}\subset\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}$; under Condition 5.1 and 5.2, we have $\mathcal{N}_{\widetilde{F}}\cap\mathcal{N}_{\varepsilon}\subset\mathcal{N}_{\frac{9}{10}\varepsilon}$.

If $(\bm{h},\bm{x})\notin\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}$, by the definition of $G$ in (3.6), at least one component in $G$ exceeds $\rho G_{0}\left(\frac{2d_{0}}{d}\right)$. We have

where $\rho\geq d^{2}+2\|\bm{e}\|^{2}$ and $0.9d_{0}\leq d\leq 1.1d_{0}.$ This implies $(\bm{h},\bm{x})\notin\mathcal{N}_{\widetilde{F}}$ and hence $\mathcal{N}_{\widetilde{F}}\subset\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}$.

For any $(\bm{h},\bm{x})\in\mathcal{N}_{\widetilde{F}}\cap\mathcal{N}_{\varepsilon}$, we have $(\bm{h},\bm{x})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$ now. By (5.6),

Therefore, $(\bm{h},\bm{x})\in\mathcal{N}_{\frac{9}{10}\varepsilon}$ and $\mathcal{N}_{\widetilde{F}}\cap\mathcal{N}_{\varepsilon}\subset\mathcal{N}_{\frac{9}{10}\varepsilon}$.

This lemma implies that the intersection of $\mathcal{N}_{\widetilde{F}}$ and the boundary of $\mathcal{N}_{\varepsilon}$ is empty. One might believe this suggests that $\mathcal{N}_{\widetilde{F}}\subset\mathcal{N}_{\varepsilon}$. This may not be true. A more reasonable interpretation is that $\mathcal{N}_{\widetilde{F}}$ consists of several disconnected regions due to the non-convexity of $\widetilde{F}(\bm{h},\bm{x})$, and one or several of them are contained in $\mathcal{N}_{\varepsilon}$.

Denote $\bm{z}_{1}=(\bm{h}_{1},\bm{x}_{1})$ and $\bm{z}_{2}=(\bm{h}_{2},\bm{x}_{2})$. Let $\bm{z}(\lambda):=(1-\lambda)\bm{z}_{1}+\lambda\bm{z}_{2}$. If $\bm{z}_{1}\in\mathcal{N}_{\varepsilon}$ and $\bm{z}(\lambda)\in\mathcal{N}_{\widetilde{F}}$ for all $\lambda\in$, we have $\bm{z}_{2}\in\mathcal{N}_{\varepsilon}$.

Let us prove the claim by contradiction. If $\bm{z}_{2}\notin\mathcal{N}_{\varepsilon}$, since $\bm{z}_{1}\in\mathcal{N}_{\varepsilon}$, there exists $\bm{z}(\lambda_{0}):=(\bm{h}(\lambda_{0}),\bm{x}(\lambda_{0}))\in\mathcal{N}_{\varepsilon}$ for some $\lambda_{0}\in$, such that $\|\bm{h}(\lambda_{0})\bm{x}(\lambda_{0})^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}=\varepsilon d_{0}$. However, since $\bm{z}(\lambda_{0})\in\mathcal{N}_{\widetilde{F}}$, by Lemma 5.5, we have $\|\bm{h}(\lambda_{0})\bm{x}(\lambda_{0})^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}\leq\frac{9}{10}\varepsilon d_{0}$. This leads to a contradiction.

Lemma 5.6 tells us that if one line segment is completely inside $\mathcal{N}_{\widetilde{F}}$ with one end point in $\mathcal{N}_{\varepsilon}$, then this whole line segment lies in $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}.$

Let the stepsize $\eta\leq\frac{1}{C_{L}}$, $\bm{z}_{t}:=(\bm{u}_{t},\bm{v}_{t})\in\hbox{\msbm{C}}^{K+N}$ and $C_{L}$ be the constant defined in (5.9). Then, as long as $\bm{z}_{t}\in\mathcal{N}_{\varepsilon}\cap\mathcal{N}_{\widetilde{F}}$, we have $\bm{z}_{t+1}\in\mathcal{N}_{\varepsilon}\cap\mathcal{N}_{\widetilde{F}}$ and

It suffices to prove (5.10). If $\nabla\widetilde{F}(\bm{z}_{t})=\bm{0}$, then $\bm{z}_{t+1}=\bm{z}_{t}$, which implies (5.10) directly. So we only consider the case when $\nabla\widetilde{F}(\bm{z}_{t})\neq\bm{0}$. Define the function

since $\varphi(\lambda)$ is a real-valued function with complex variables (See (6.1) for details). By the definition of derivatives, we know there exists $\eta_{0}>0$, such that $\varphi(\lambda)<\varphi(0)$ for all $0<\lambda\leq\eta_{0}$. Now we will first prove that $\varphi(\lambda)\leq\varphi(0)$ for all $0\leq\lambda\leq\eta$ by contradiction. Assume there exists some $\eta_{1}\in(\eta_{0},\eta]$ such that $\varphi(\eta_{1})>\varphi(0)$. Then there exists $\eta_{2}\in(\eta_{0},\eta_{1})$, such that $\varphi(\eta_{2})=\varphi(0)$ and $\varphi(\lambda)<\varphi(0)$ for all $0<\lambda<\eta_{2}$, since $\varphi(\lambda)$ is a continuous function. This implies

since $\widetilde{F}(\bm{z}_{t}-\lambda\nabla\widetilde{F}(\bm{z}_{t}))\leq\widetilde{F}(\bm{z}_{t})$ for $0\leq\lambda\leq\eta_{2}.$ By Lemma 5.6 and the assumption $\bm{z}_{t}\in\mathcal{N}_{\varepsilon}$, we have

Then, by using the modified descent lemma (Lemma 6.1),

where the final inequality is due to $\eta_{2}/\eta<1$, $\eta_{2}>\eta_{0}\geq 0$ and $\nabla\widetilde{F}(\bm{z}_{t})\neq\bm{0}$. This contradicts $\widetilde{F}(\bm{z}_{t}-\eta_{2}\nabla\widetilde{F}(\bm{z}_{t}))=\varphi(\eta_{2})=\varphi(0)=\widetilde{F}(\bm{z}_{t})$.

Therefore, there holds $\varphi(\lambda)\leq\varphi(0)$ for all $0\leq\lambda\leq\eta$. Similarly, we can prove

which implies $\bm{z}_{t+1}=\bm{z}_{t}-\eta\nabla\widetilde{F}(\bm{z}_{t})\in\mathcal{N}_{\varepsilon}\cap\mathcal{N}_{\widetilde{F}}$. Again, by using Lemma 6.1 we can prove

where the final inequality is due to $\eta\leq\frac{1}{C_{L}}$.

We conclude this subsection by proving Theorem 3.2 under the Local regularity condition, the Local RIP condition, the Robustness condition, and the Local smoothness condition. The next subsections are devoted to justifying these conditions and showing that they hold under the assumptions of Theorem 3.2.

[of Theorem 3.2] Suppose that the initial guess $\bm{z}_{0}:=(\bm{u}_{0},\bm{v}_{0})\in\frac{1}{\sqrt{3}}\mathcal{N}_{d_{0}}\bigcap\frac{1}{\sqrt{3}}\mathcal{N}_{\mu}\bigcap\mathcal{N}_{\frac{2}{5}\varepsilon}$, we have $G(\bm{u}_{0},\bm{v}_{0})=0$. This holds, because

where $\|\bm{u}_{0}\|\leq\frac{2\sqrt{d_{0}}}{\sqrt{3}}$, $\sqrt{L}\|\bm{B}\bm{u}_{0}\|_{\infty}\leq\frac{4\sqrt{d_{0}}\mu}{\sqrt{3}}$ and $\frac{9}{10}d_{0}\leq d\leq\frac{11}{10}d_{0}.$ Therefore $G_{0}\left(\frac{\|\bm{u}_{0}\|^{2}}{2d}\right)=G_{0}\left(\frac{\|\bm{v}_{0}\|^{2}}{2d}\right)=G_{0}\left(\frac{L|\bm{b}_{l}^{*}\bm{u}_{0}|^{2}}{8d\mu^{2}}\right)=0$ for all $1\leq l\leq L$ and $G(\bm{u}_{0},\bm{v}_{0})=0.$ Since $(\bm{u}_{0},\bm{v}_{0})\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$, (5.5) combined with $\delta(\bm{z}_{0}):=\frac{\|\bm{u}_{0}\bm{v}_{0}^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}}{d_{0}}\leq\frac{2\varepsilon}{5}$ imply that

and hence $\bm{z}_{0}=(\bm{u}_{0},\bm{v}_{0})\in\mathcal{N}_{\varepsilon}\bigcap\mathcal{N}_{\widetilde{F}}.$ Denote $\bm{z}_{t}:=(\bm{u}_{t},\bm{v}_{t}).$ Combining Lemma 5.8 by choosing $\eta\leq\frac{1}{C_{L}}$ with Condition 5.3, we have

with $c=\|\bm{e}\|^{2}+a\|\mathcal{A}^{*}(\bm{e})\|^{2}$, $a=1700$ and $\bm{z}_{t}\in\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$ for all $t\geq 0.$ Obviously, the inequality above implies

and by monotonicity of $z_{+}=\frac{z+|z|}{2}$, there holds

where $\widetilde{F}(\bm{z}_{0})\leq\frac{1}{3}\varepsilon^{2}d_{0}^{2}+\|\bm{e}\|^{2}$ and hence $\left[\widetilde{F}(\bm{z}_{0})-c\right]_{+}\leq\left[\frac{1}{3}\varepsilon^{2}d_{0}^{2}-a\|\mathcal{A}^{*}(\bm{e})\|^{2}\right]_{+}\leq\frac{1}{3}\varepsilon^{2}d_{0}^{2}.$ Now we can conclude that $\left[\widetilde{F}(\bm{z}_{t})-c\right]_{+}$ converges to geometrically. Note that over $\mathcal{N}_{d_{0}}\cap\mathcal{N}_{\mu}\cap\mathcal{N}_{\varepsilon}$,

where $\delta(\bm{z}_{t}):=\frac{\|\bm{u}_{t}\bm{v}_{t}^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}}{d_{0}}$, $F_{0}$ is defined in (2.9) and $G(\bm{z}_{t})\geq 0$. There holds

Solving the inequality above for $\delta(\bm{z}_{t})$, we have

Let $d_{t}:=\|\bm{u}_{t}\|\|\bm{v}_{t}\|$, $t\geq 1.$ By (5.11) and triangle inequality, we immediately conclude that

Now we derive the upper bound for $\sin\angle(\bm{u}_{t},\bm{h}_{0})$ and $\sin\angle(\bm{v}_{t},\bm{x}_{0}).$ Due to symmetry, it suffices to consider $\sin\angle(\bm{u}_{t},\bm{h}_{0})$. The bound follows from standard linear algebra arguments:

where the second equality uses $\left(\bm{I}-\frac{\bm{h}_{0}\bm{h}_{0}^{*}}{d_{0}}\right)\bm{h}_{0}=\bm{0}.$

2 Supporting lemmata

This subsection introduces several lemmata, especially Lemma 5.9, 5.12 and 5.13, which are central for justifying Conditions 5.1 and 5.3. After that, we will prove the Local RIP Condition in Lemma 5.14 based on those three lemmata. We start with defining a linear space $T$, which contains $\bm{h}_{0}\bm{x}_{0}^{*}$, via

Denote $\mathcal{P}_{T}$ to be the projection operator from $\hbox{\msbm{C}}^{K\times N}$ onto $T$.

For any $\bm{h}$ and $\bm{x}$, there are unique orthogonal decompositions

Recall that $\|\bm{h}_{0}\|=\|\bm{x}_{0}\|=\sqrt{d_{0}}$. If $\delta:=\frac{\|\bm{h}\bm{x}^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}}{d_{0}}<1$, we have the following useful bounds

This lemma is actually a simple version of singular value/vector perturbation. It says that if $\frac{\|\bm{h}\bm{x}^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}}{d_{0}}$ is of $\mathcal{O}(\delta)$, then the individual vectors $(\bm{h},\bm{x})$ are also close to $(\bm{h}_{0},\bm{x}_{0})$, with the error of order $\mathcal{O}(\delta).$

The equality (5.13) implies that $\|\alpha_{1}\bm{h}_{0}\|\leq\|\bm{h}\|$, so there holds $|\alpha_{1}|\leq\frac{\|\bm{h}\|}{\|\bm{h}_{0}\|}$. Since $\|\bm{h}\bm{x}^{*}-\bm{h}_{0}\bm{x}_{0}^{*}\|_{F}=\delta d_{0}$, by (5.14), we have

In the following, we introduce and prove a series of local and global properties of $\mathcal{A}$:

with probability at least $1-L^{-\gamma}.$

Let $\mathcal{A}$ be the operator defined in (2.6), then on an event $E_{1}$ with probability at least $1-L^{-\gamma}$, $\mathcal{A}$ restricted on $T$ is well-conditioned, i.e.,

where $\mathcal{P}_{T}$ is the projection operator from $\hbox{\msbm{C}}^{K\times N}$ onto $T$, provided $L\geq C_{\gamma}\max\{K,\mu_{h}^{2}N\}\log^{2}(L)$.

Now we introduce a property of $\mathcal{A}$ when restricted on rank-one matrices.

On an event $E_{2}$ with probability at least $1-L^{-\gamma}-\frac{1}{\gamma}\exp(-(K+N))$, we have

uniformly for any $\bm{u}$ and $\bm{v}$, provided $L\geq C_{\gamma}(K+N)\log L$.

Due to the homogeneity, without loss of generality we can assume $\|\bm{u}\|=\|\bm{v}\|=1$. Define

It suffices to prove that $f(\bm{u},\bm{v})\leq\frac{4}{3}$ uniformly for all $(\bm{u},\bm{v})\in\mathcal{S}^{K-1}\times\mathcal{S}^{N-1}$ with high probability, where $\mathcal{S}^{K-1}$ is the unit sphere in $\hbox{\msbm{C}}^{K}.$ For fixed $(\bm{u},\bm{v})\in\mathcal{S}^{K-1}\times\mathcal{S}^{N-1}$, notice that

is the sum of subexponential variables with expectation $\hbox{\msbm{E}}\|\mathcal{A}(\bm{u}\bm{v}^{*})\|^{2}=\sum\limits_{l=1}^{L}|\bm{b}_{l}^{*}\bm{u}|^{2}=1$. For any generalized $\chi_{n}^{2}$ variable $Y\sim\sum_{i=1}^{n}c_{i}\xi_{i}^{2}$ satisfies

where $\{\xi_{i}\}$ are i.i.d. $\chi^{2}_{1}$ random variables and $\bm{c}=(c_{1},\cdots,c_{n})^{T}\in\hbox{\msbm{R}}^{n}$. Here we set $|\bm{a}_{l}^{*}\bm{v}|^{2}=\frac{1}{2}\xi_{2l-1}^{2}+\frac{1}{2}\xi_{2l}^{2}$, $c_{2l-1}=c_{2l}=\frac{|\bm{b}_{l}^{*}\bm{u}|^{2}}{2}$ and $n=2L$. Therefore,

That is, $f(\bm{u},\bm{v})\leq 1$ with probability at least $1-\exp\left(-2(K+N)(\log L)\right)$. We define $\mathcal{K}$ and $\mathcal{N}$ as $\varepsilon_{0}$-nets of $\mathcal{S}^{K-1}$ and $\mathcal{S}^{N-1}$, respectively. Then, $|\mathcal{K}|\leq(1+\frac{2}{\varepsilon_{0}})^{2K}$ and $|\mathcal{N}|\leq(1+\frac{2}{\varepsilon_{0}})^{2N}$ follow from the covering numbers of the sphere (Lemma 5.2 in ).

By taking the union bounds over $\mathcal{K}\times\mathcal{N},$ we have $f(\bm{u},\bm{v})\leq 1$ holds uniformly for all $(\bm{u},\bm{v})\in\mathcal{K}\times\mathcal{N}$ with probability at least

Our goal is to show that $f(\bm{u},\bm{v})\leq\frac{4}{3}$ uniformly for all $(\bm{u},\bm{v})\in\mathcal{S}^{K-1}\times\mathcal{S}^{N-1}$ with the same probability. For any $(\bm{u},\bm{v})\in\mathcal{S}^{K-1}\times\mathcal{S}^{N-1}$, we can find its closest $(\bm{u}_{0},\bm{v}_{0})\in\mathcal{K}\times\mathcal{N}$ satisfying $\|\bm{u}-\bm{u}_{0}\|\leq\varepsilon_{0}$ and $\|\bm{v}-\bm{v}_{0}\|\leq\varepsilon_{0}$. By Lemma 5.11, with probability at least $1-L^{-\gamma}$, we have $\|\mathcal{A}\|\leq\sqrt{(N+\gamma)\log L}$. Then straightforward calculation gives