On the Global Geometry of Sphere-Constrained Sparse Blind Deconvolution

TL;DR

This paper explores the geometric landscape of sphere-constrained sparse blind deconvolution, revealing that local minima are near true solutions, and proposes a two-stage algorithm for effective recovery in microscopy and image deblurring.

Contribution

It characterizes the global geometry of the nonconvex sphere-constrained blind deconvolution problem and introduces a novel two-stage algorithm inspired by this geometry.

Findings

Local minima are close to shifted ground truths under certain conditions.

The proposed algorithm successfully recovers ground truth in microscopy and deblurring tasks.

Insights extend to convolutional dictionary learning with multiple signals.

Abstract

Blind deconvolution is the problem of recovering a convolutional kernel $\boldsymbol a_0$ and an activation signal $\boldsymbol x_0$ from their convolution $\boldsymbol y = \boldsymbol a_0 \circledast \boldsymbol x_0$. This problem is ill-posed without further constraints or priors. This paper studies the situation where the nonzero entries in the activation signal are sparsely and randomly populated. We normalize the convolution kernel to have unit Frobenius norm and cast the sparse blind deconvolution problem as a nonconvex optimization problem over the sphere. With this spherical constraint, every spurious local minimum turns out to be close to some signed shift truncation of the ground truth, under certain hypotheses. This benign property motivates an effective two stage algorithm that recovers the ground truth from the partial information offered by a suboptimal local minimum. This geometry-inspired algorithm recovers the ground truth for certain microscopy problems, also exhibits promising performance in the more challenging image deblurring problem. Our insights into the global geometry and the two stage algorithm extend to the convolutional dictionary learning problem, where a superposition of multiple convolution signals is observed.

Equations

$$\mathbf{y}({\mathbf{a}}_0,{\mathbf{x}}_0)={\mathbf{a}}_0⊛{\mathbf{x}}_0=\widetilde{{\mathbf{a}}_0}⊛{\mathbf{x}}_0\in {\mathbb{R}}^m.$$

$$\mathbf{y}({\mathbf{a}}_0,{\mathbf{x}}_0)={\mathbf{C}}_{\widetilde{{\mathbf{a}}_0}}{\mathbf{x}}_0={\mathbf{C}}_{{\mathbf{x}}_0}\widetilde{{\mathbf{a}}_0}.$$

$$s_{\tau }\left[\mathbf{v}\right](i)=\mathbf{v}([i-\tau {]}_m),\forall i\in [0,\cdots ,m-1].$$

$$\mathbf{y}=\left(\alpha s_{\tau }[\widetilde{{\mathbf{a}}_0}]\right)⊛\left({\alpha }^{-1}{s}_{-\tau }[{\mathbf{x}}_0]\right).$$

$$\left({\iota }^{\ast }{s}_{\tau }[\widetilde{{\mathbf{a}}_0}]\right)⊛s_{-\tau }[{\mathbf{x}}_0]\approx \mathbf{y},$$

$$\underset{\mathbf{a}\in {\mathbb{S}}^{k-1},\mathbf{x}}{\min }\psi (\mathbf{a},\mathbf{x})\doteq \frac{1}{2}\Vert \mathbf{y}-\mathbf{a}⊛\mathbf{x}{\Vert }_2^2+\lambda r(\mathbf{x}).$$

$$\phi (\mathbf{a})\doteq \underset{\mathbf{x}}{\min }\psi (\mathbf{a},\mathbf{x}),$$

$${\mathbf{x}}^{\ast }(\mathbf{a})=\arg \underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{y}-\mathbf{a}⊛\mathbf{x}\Vert }_2^2+\lambda {\Vert \mathbf{x}\Vert }_1.$$

$$I=\mathrm{supp}\left({\mathbf{x}}^{\ast }\right),\sigma =\mathrm{sign}({\mathbf{x}}^{\ast }),$$

$${\mathbb{S}}^{k-1}=\bigcup _{\sigma }{\mathcal{R}}_{\sigma },{\mathcal{R}}_{\sigma }=\left\{\mathbf{a}\mid \mathrm{sign}({\mathbf{x}}^{\ast }(\mathbf{a}))=\sigma \right\}.$$

$${\mathbf{x}}_I^{\ast }(\mathbf{a})={\left({\mathbf{C}}_{\mathbf{a}}^{\ast }{\mathbf{C}}_{\mathbf{a}}\right)}_I^{-1}{\left({\mathbf{C}}_{\mathbf{a}}^{\ast }\mathbf{y}-\lambda \sigma \right)}_I.$$

$${\phi }_{\sigma }(\mathbf{a})=-\frac{1}{2}{\left({\mathbf{C}}_{\mathbf{a}}^{\ast }\mathbf{y}-\lambda \sigma \right)}_I^{\ast }{\left({\mathbf{C}}_{\mathbf{a}}^{\ast }{\mathbf{C}}_{\mathbf{a}}\right)}_I^{-1}{\left({\mathbf{C}}_{\mathbf{a}}^{\ast }\mathbf{y}-\lambda \sigma \right)}_I+\frac{1}{2}{\Vert \mathbf{y}\Vert }_2^2.$$

$$\text{minimize}\hat{\phi }(\mathbf{a})\text{subject to}\mathbf{a}\in {\mathbb{S}}^{k-1},$$

$$\hat{\phi }(\mathbf{a})\doteq \underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \widetilde{{\mathbf{a}}_0}\Vert }_2^2+\frac{1}{2}{\Vert \mathbf{x}\Vert }_2^2-⟨\tilde{\mathbf{a}}⊛\mathbf{x},\widetilde{{\mathbf{a}}_0}⟩+\lambda {\Vert \mathbf{x}\Vert }_1.$$

$${\mathbf{x}}^{\ast }(\mathbf{a})={\mathrm{SOFT}}_{\lambda }\left[{\mathbf{C}}_{\mathbf{a}}^{\ast }\widetilde{{\mathbf{a}}_0}\right]={\mathrm{SOFT}}_{\lambda }\left[{\check{\mathbf{C}}}_{{\mathbf{a}}_0}^{\ast }\iota \mathbf{a}\right],$$

$$\check{\mathbf{C}}_{\mathbf{a}_{0}}=\left[\begin{array}[]{c|c|c|c}s_{0}\left[\widetilde{\mathbf{a}_{0}}\right]&s_{-1}\left[\widetilde{\mathbf{a}_{0}}\right]&\dots&s_{-(m-1)}\left[\widetilde{\mathbf{a}_{0}}\right]\end{array}\right].$$

$${\hat{\phi }}_{\sigma }(\mathbf{a})=-\frac{1}{2}{\Vert {\left({\mathbf{C}}_{\mathbf{a}}^{\ast }\widetilde{{\mathbf{a}}_0}-\lambda \sigma \right)}_I\Vert }_2^2+\frac{1}{2}{\Vert \widetilde{{\mathbf{a}}_0}\Vert }_2^2.$$

$$\mathcal{I}=\left\{\mathrm{supp}\left({\mathrm{SOFT}}_{\lambda }\left[{\check{\mathbf{C}}}_{{\mathbf{a}}_0}^{\ast }\iota \mathbf{a}\right]\right)\mid \mathbf{a}\in {\mathbb{S}}^{k-1}\right\}.\text{vspace}-.05in$$

$${\mathbf{W}}_I=\left[\frac{{\iota }^{\ast }{s}_{-i_1}\left[\widetilde{{\mathbf{a}}_0}\right]}{{\Vert {\iota }^{\ast }{s}_{-i_1}\left[\widetilde{{\mathbf{a}}_0}\right]\Vert }_F}\left|\frac{{\iota }^{\ast }{s}_{-i_2}\left[\widetilde{{\mathbf{a}}_0}\right]}{{\Vert {\iota }^{\ast }{s}_{-i_2}\left[\widetilde{{\mathbf{a}}_0}\right]\Vert }_F}\right|\dots |\frac{{\iota }^{\ast }{s}_{-i_{|I|}}\left[\widetilde{{\mathbf{a}}_0}\right]}{{\Vert {\iota }^{\ast }{s}_{-i_{|I|}}\left[\widetilde{{\mathbf{a}}_0}\right]\Vert }_F}\right]$$

$$\Vert {\mathbf{W}}_I^{\ast }{\mathbf{W}}_I-\mathbf{I}\Vert <\frac{{\lambda }^2}{6},$$

$$\bar{\mathbf{a}}=\pm \frac{{\iota }^{\ast }{s}_{\tau }\left[\widetilde{{\mathbf{a}}_0}\right]}{{\Vert {\iota }^{\ast }{s}_{\tau }\left[\widetilde{{\mathbf{a}}_0}\right]\Vert }_F}=\pm {\mathcal{P}}_{\mathbb{S}}\left[{\iota }^{\ast }{s}_{-\tau }\left[\widetilde{{\mathbf{a}}_0}\right]\right]$$

$$h_{\mu }(\mathbf{x})=\sum _{|x_i|\le \mu }\left(\frac{x_i^2}{2\mu }+\frac{\mu }{2}\right)+\sum _{|x_i|>\mu }\left|x_i\right|$$

$${\mathcal{R}}_{h,0}:=\{\mathbf{a}\in {\mathbb{S}}^{k-1}:{\Vert {\mathbf{x}}^{\ast }(\mathbf{a})\Vert }_{\infty }\le \mu \}.$$

$${\phi }_{h_{\mu }}(\mathbf{a})=\frac{1}{2}\Vert \mathbf{y}-\mathbf{a}⊛\mathbf{x}{\Vert }_2^2+\frac{\lambda }{2\mu }\Vert \mathbf{x}{\Vert }_2^2+\frac{\mu n}{2},$$

$${\mathbf{x}}^{\ast }(\mathbf{a})={\left({\mathbf{C}}_{\mathbf{a}}^{\ast }{\mathbf{C}}_{\mathbf{a}}+\frac{\lambda }{\mu }\mathbf{I}\right)}^{-1}{\mathbf{C}}_{\mathbf{a}}^{\ast }\mathbf{y}\approx \frac{\mu }{\lambda }{\mathbf{C}}_{\mathbf{a}}^{\ast }\mathbf{y}.$$

$${\phi }_{h_{\mu }}(\mathbf{a})\approx -\frac{\mu }{2\lambda }{\Vert \mathbf{y}⊛\mathbf{a}\Vert }_2^2+\frac{1}{2}{\Vert \mathbf{y}\Vert }_2^2+\frac{\mu n}{2}.$$

$$\left|⟨{\mathbf{a}}_0,\iota s_{\tau \ne 0}\left[\widetilde{{\mathbf{a}}_0}\right]⟩\right|<{\lambda }_{rel}^2-\left(2+1/{\lambda }_{rel}^2\right)\sqrt{1-{\lambda }_{rel}^2}$$

$${\nabla }_x\mathbf{y}={\mathbf{a}}_0⊛{\nabla }_x{\mathbf{x}}_0,{\nabla }_y\mathbf{y}={\mathbf{a}}_0⊛{\nabla }_y{\mathbf{x}}_0.$$

$$\underset{\mathbf{a}\in {\mathbb{S}}_{+}^{k-1},{\mathbf{x}}_1,{\mathbf{x}}_2}{\min }$$

$$\mathbf{y}={\sum }_{n=1}^N{\mathbf{a}}_{0n}⊛{\mathbf{x}}_{0n}.$$

Full text

1 Introduction

Blind deconvolution aims to recover two unknown signals: a kernel $\mathbf{a}_{0}$ and some underlying signal $\mathbf{x}_{0}$ from their convolution $\mathbf{y}=\mathbf{a}_{0}\circledast\mathbf{x}_{0}$. Blind deconvolution is ill-posed in general: there are infinitely many pairs of signals rendering the same convolution. To render the problem well-posed, one may exploit prior knowledge about the structure of $\mathbf{a}_{0}$ and $\mathbf{x}_{0}$. For example, the underlying signal $\mathbf{x}_{0}$ is sparse in many engineering and scientific applications:

Microscopy data analysis: In the crystal lattice of nanoscale materials, there exist randomly and sparsely distributed “defects”, whose locations and signatures encode crucial information about the electronic structure of the material. Accurate recovery of such information can facilitate investigation of the detailed structure of materials [CLC*+*17].

Neural spike sorting: Neurons communicate by firing brief voltage spikes, whose characteristics reflect important features of the neuron. These spikes occur randomly and sparsely in time. Neurophysiologists are interested in assigning stereotyped spikes to putative cells, as well as in knowing their respective spike times [ETS11, Lew98].

Image deblurring: Motion blur can be modeled as the convolution of a latent sharp image and a kernel capturing the motion of the camera, usually assumed to be invariant across the image [FSH*+*06]. The inverse process of recovering the original sharp image from a blurry image has been widely studied [CL09, KH96]. Many well-performing approaches leverage the observation that sharp natural images typically have (approximately) sparse gradients [CW98, LWDF11, PF14].

All of these applications lead to instances of the sparse blind deconvolution problem. The dominant algorithmic approach to sparse blind deconvolution involves nonconvex optimization111In signal processing, a number of elegant convex relaxations of the problem have been developed [ARR12, Chi16, LLSW18]. However, these approaches typically require stronger prior information (subspace constraint rather than sparsity) or exhibit suboptimal scalings.. Nonconvex formulations for deconvolution can be derived via several probabilistic formalisms (ML/MAP, VB, ect.), or simply from heuristics. For example, in image deblurring, the kernel $\mathbf{a}$ can be modeled as residing on a simplex [LWDF11, GTZ*+*16, KTF11, LCM14]. This is natural from a modeling prospective222Since entries of $\mathbf{a}$ roughly represent the fraction of the camera exposure time at a given location., but problematic for optimization: natural formulations of deconvolution over the simplex admit trivial global minimizers (corresponding to spiky convolutional kernels $\mathbf{a}=\mathbf{\delta}$) [PF14, BVG13], which provide no information about the ground truth. Practical remedies for this problem include exploiting additional data priors [GTZ*+*16, LCM14, XWHH16] or careful initialization via edge restoration or multi-scale refinement [KTF11, XJ10], to avoid the trivial spiky global minima.

In contrast, motivated by a careful comparison of MAP and VB approaches, [WZ13, ZWZ13] propose to instead constrain $\mathbf{a}$ to have unit Frobenius norm – i.e., to reside on a high-dimensional sphere.333[WZ13] contains a wealth of additional ideas about the role of sparsity-promoting priors in obtaining good local minima, and on the probabilistic underpinnings of deconvolution. Our experiments support the viewpoint that the key insight in [WZ13] is the role of the spherical constraint in avoiding bad minimizers. This choice is arguably more appropriate for certain scientific applications – such as microscopy – in which the kernel $\mathbf{a}$ can have negative entries. For image deblurring, $\mathbf{a}$ can be assumed to be nonnegative, and the sphere constraint seems less natural from a modeling perspective.

In this paper, we study the geometry of sphere-constrained sparse blind deconvolution. Our goal is to understand when simple algorithms based on nonconvex optimization can exactly recover the convolutional kernel $\mathbf{a}$ and the sparse signal $\mathbf{x}$. This goal is motivated by the applications described above – in particular, microscopy data analysis – in which there is a strong, physical sparsity prior and a clear, physical notion of the ground truth. We develop our theory and algorithms under the assumption that $\mathbf{a}$ is a short kernel, and that $\mathbf{x}$ is sparsely and randomly supported. We demonstrate through a theoretical analysis of certain (idealized) cases and many numerical experiments that when these assumptions are satisfied, the proposed algorithm correctly recovers $\mathbf{a}$, and hence $\mathbf{x}$. These results stem from a striking geometric property of sphere-constrained sparse blind deconvolution: although the problem is still nonconvex, every local minimizer $\mathbf{\bar{a}}$ is very close to a signed shift-truncation of the ground truth kernel $\mathbf{a}_{0}$. This observation provides a geometric explanation of how the sphere constraint can facilitate sparse blind deconvolution.

The remainder of this paper is organized as follows. Section 2 discusses the intrinsic symmetries associated with the convolutional operator and their implication on the geometry of sphere-constrained sparse blind deconvolution. Section 3 introduces the optimization-based two stage algorithm and some related technical details. Section 4 discusses two other important extensions in image deblurring and convolutional dictionary learning. Section 5 gives experimental corroboration of our theory, and shows promising results on microscopy data analysis, image deblurring, and convolutional dictionary learning. Section 6 discusses directions for future work.

For simplicity, we assume that the convolutional signals are one dimensional in both our problem formulation and technical proofs; all of our results extend naturally to higher-dimensional signals. Throughout this paper, vectors $\mathbf{v}\in\mathbb{R}^{k}$ are indexed as $\mathbf{v}=[v_{0},v_{1},\cdots,v_{k-1}]$, and $[\cdot]_{m}$ denotes the modulo-$m$ operation. We use $\left\|\cdot\right\|$ to denote the operator norm, and $\left\|\cdot\right\|_{p}$ to denote the entry wise $\ell^{p}$ norm. A projection onto the Frobenius sphere is denoted with $\mathcal{P}_{\mathbb{S}}\left[\cdot\right]=\frac{\cdot}{\left\|\cdot\right\|_{2}}$, and a projection onto subset $I$ is denoted with $(\cdot)_{I}$.

2 Symmetry and Global Geometry

Without loss of generality, we assume the observation data $\mathbf{y}$ is generated via a circular convolution $\circledast$ of the ground truth $\mathbf{a}_{0}\in\mathbb{R}^{k}$ and $\mathbf{x}_{0}\in\mathbb{R}^{m}$:

[TABLE]

Here, $\widetilde{\mathbf{a}_{0}}\in\mathbb{R}^{m}$ denote the zero padded $m$-length version of $\mathbf{a}_{0}$, which can be expressed as $\widetilde{\mathbf{a}_{0}}=\mathbf{\iota}\mathbf{a}_{0}$ with $\mathbf{\iota}:\mathbb{R}^{k}\to\mathbb{R}^{m}$ be a zero padding operator. Its adjoint $\mathbf{\iota}^{*}:\mathbb{R}^{m}\to\mathbb{R}^{k}$ acts as a projection to lower dimensional space by keeping the first $k$ components. Equivalently, we can write

[TABLE]

Here, $\mathbf{C}_{\mathbf{v}}\in\mathbb{R}^{m\times m}$ is the circulant matrix generated from vector $\mathbf{v}$, whose $j$-th column is a cyclic shift $s_{j-1}[\mathbf{v}]$ of the vector $\mathbf{v}$:

[TABLE]

2.1 Symmetries and Symmetry Breaking

The SBD problem exhibits a scaled-shift symmetry, which derives from the symmetries of the convolution operator. Namely, given a pair $(\mathbf{a}_{0},\mathbf{x}_{0})$ satisfying $\mathbf{y}=\mathbf{a}_{0}\circledast\mathbf{x}_{0}$, for any nonzero scalar $\alpha$ and integer $\tau$

[TABLE]

Note that a scaled shift $\alpha^{-1}s_{-\tau}[\mathbf{x}_{0}]$ of a sparse signal $\mathbf{x}_{0}$ remains sparse, and that a scaled shift $\alpha s_{\tau}[\widetilde{\mathbf{a}_{0}}]$ of a length-$k$ kernel $\mathbf{a}_{0}$ still has $k$-nonzero entries. So, these symmetries are intrinsic to the SBD problem, as formulated here. We can only hope to recover $(\mathbf{a}_{0},\mathbf{x}_{0})$ up to this symmetry.

The presence of nontrivial symmetries is a hallmark of bilinear problems arising in practice – see, e.g., [SQW17a, SQW17b] for examples from dictionary learning and generalized phase retrieval. Symmetries render straightforward approaches to convexify the problem ineffective.444Given any set of points where the convex objective function achieves equal values, the function value will be no larger at any convex combination of them. They also raise challenges for nonconvex optimization: equivalent symmetric solutions correspond to multiple disconnected global optima. This creates a very complicated objective landscape, which could potentially also contain spurious local optimizers. Certain highly symmetric nonconvex problems arising in signal processing do not exhibit spurious minimizers [SQW17a, SQW17b], however, proving this can be challenging.

2.1.1 Symmetry breaking.

We employ a weak symmetry breaking mechanism by constraining $\mathbf{a}\in\mathbb{S}^{k-1}$ 555This is motivated in part by [SQW17a], which demonstrates that a certain formulation of the dictionary learning problem over the sphere has no spurious local minimizers, even for relatively dense target representations. The “simplex constrained” analogue of that work, which optimizes over hyperplanes, requires the target solution to be much sparser [SWW12].: we reduce the scale ambiguity to a sign ambiguity, by constraining $\mathbf{a}$ to have unit Frobenius norm; we mitigate the shift ambiguity by constraining $\mathbf{a}$ to be supported on the first $k$ entries.

In general, $s_{\tau}[\widetilde{\mathbf{a}_{0}}]$ are not be supported on the first $k$ entries, hence constraining $\mathbf{a}$ to be supported on the first $k$ entries removes the shift symmetry. However, effects of such shift symmetry still persist. Since the restriction $\mathbf{\iota}^{*}s_{\tau}[\widetilde{\mathbf{a}_{0}}]$ to the first $k$ entries can be convolved with the sparse signal $s_{-\tau}[\mathbf{x}_{0}]$ to produce a near approximation to $\mathbf{y}$:

[TABLE]

especially when the shift $\left|\tau\right|$ is small. We will see that (i) these symmetric solutions $\mathbf{\iota}^{*}s_{\tau}[\widetilde{\mathbf{a}_{0}}]$ persist as local minima of a natural optimization formulation of the SBD problem, but that (ii) under conditions, these are the only local minima.

2.2 Global Geometry on the Sphere

We study the following objective function, which can be viewed as balancing sparsity of $\mathbf{x}$ with fidelity to the observation $\mathbf{y}$: 666Similar formulation can be found in lot of sparse representation problems [MBP14].

[TABLE]

When $\mathbf{x}_{0}$ is long and random, it is more convenient to study this function through its “marginalization”

[TABLE]

which is defined over the sphere $\mathbb{S}^{k-1}$.

In Figure 1, we plot the function value of $\varphi(\mathbf{a})$ on the sphere $\mathbf{a}\in\mathbb{S}^{2}$: red and blue imply larger and smaller objective value respectively and there are several local minima. For this highly nonconvex function, the ground truth $\mathbf{a}_{0}$ achieves the global minimum, while other local minima $\mathbf{\bar{a}}$ are very close to certain signed shift truncations of the ground truth. Figure 1 (right) exhibits an example of a local minimum in a higher-dimensional problem.

2.2.1 Analysis under Restricted Settings

Demonstrating that this observation holds in general is challenging: for most reasonable choices of the regularizer $r$, there is no closed form expression for the objective $\varphi(\mathbf{a})$. We develop an analysis under several simplifying assumptions. Throughout, we let $r(\mathbf{x})$ be the $\ell^{1}$ norm, although similar conclusions hold for other sparsifying regularizers. With this choice, we can simplify the objective $\varphi$ by dividing the sphere $\mathbb{S}^{k-1}$ via the sign-support pattern of the minimizing $\mathbf{x}^{*}(\mathbf{a})$:

[TABLE]

Let $I$ and $\mathbf{\sigma}$ denote the support and sign of $\mathbf{x}^{*}$

[TABLE]

then the whole sphere can be divided via the sign support pattern

[TABLE]

On each $\mathcal{R}_{\mathbf{\sigma}}$ where the sign support pattern $\mathbf{\sigma}$ remains the same, the stationarity condition for minimizer $\mathbf{x}^{*}(\mathbf{a})$ implies

[TABLE]

Plugging above expression back to the original objective function $\varphi(\mathbf{a})$ yields

[TABLE]

Although the objective function $\varphi(\mathbf{a})$ can be substantially simplified by removing the $\mathbf{x}$ variable in this way, it still maintains a complicated dependence on $\mathbf{a}$. To obtain some preliminary insights, we make two simplifications for easier calculation while preserving important characteristics of the geometry of $\varphi$:

Simplification I: $\mathbf{x}_{0}=\mathbf{\delta}$. We maximally simplify the underlying sparse signal as a single spike $\mathbf{\delta}$ and the observation will be $\mathbf{y}=\mathbf{a}_{0}\circledast\mathbf{x}_{0}=\widetilde{\mathbf{a}_{0}}$. This case itself is trivial, but its function geometry is a basic but important case to be understood. This simple case also yields intuitions that carry over to the less trivial situation in which $\mathbf{x}_{0}$ is a long random vector.

Simplification II: $\mathbf{C}^{*}_{\mathbf{a}}\mathbf{C}_{\mathbf{a}}\to\mathbf{I}$. For a random $\mathbf{a}\in\mathbb{S}^{k-1}$, its expectation satisfies $\mathbb{E}\left[\mathbf{C}^{*}_{\mathbf{a}}\mathbf{C}_{\mathbf{a}}\right]=\mathbf{I}$. Here, we simply use the identity matrix to replace any $\mathbf{C}^{*}_{\mathbf{a}}\mathbf{C}_{\mathbf{a}}$ and therefore reduce the complexity of Equation 2.12.

With these two simplifications, the original objective problem can be replaced with the following:

[TABLE]

where

[TABLE]

In this case, the minimizing $\mathbf{x}^{*}(\mathbf{a})$ has a simple closed form solution:

[TABLE]

here, $\mathrm{SOFT}_{\lambda}\left[u\right]=\operatorname{sign}(u)\max\left\{|u|-\lambda,0\right\}$ is the entry-wise soft-thresholding operator and $\check{\mathbf{C}}_{\mathbf{a}_{0}}\in\mathbb{R}^{m\times m}$ is the reversed circulant matrix for $\mathbf{a}_{0}$ defined via

[TABLE]

On a constant sign support pattern $\mathbf{\sigma}$, $\widehat{\varphi}(\mathbf{a})$ can be written into a simpler quadratic form:

[TABLE]

For this surrogate $\widehat{\varphi}(\mathbf{a})$, we can show that if $\lambda$ is sufficiently large compared to the magnitude of $\mathbf{x}_{0}$, every strict local minimizer is a signed shift truncation of the ground truth:

Theorem 2.1**.**

Define the set of possible supports of minimizer $\mathbf{x}^{*}$ with

[TABLE]

For each nonempty support $I=\left\{i_{1}<i_{2}<\cdots<i_{|I|}\right\}$, let

[TABLE]

Suppose that $\lambda<1$ and that for every nonempty $I\in\mathcal{I}$,

[TABLE]

then every local minimum $\bar{\mathbf{a}}$ of $\widehat{\varphi}$ over $\mathbb{S}^{k-1}$ satisfies either $\bar{\mathbf{a}}\in\mathcal{R}_{\mathbf{0}}$ (in which case $\bar{\mathbf{a}}$ is also a global maximum), or

[TABLE]

with $\mathbf{x}^{\star}(\bar{\mathbf{a}})=\pm\mathrm{SOFT}_{\lambda}\left[\left\|\mathbf{\iota}^{*}s_{\tau}\left[\widetilde{\mathbf{a}_{0}}\right]\right\|_{F}\right]s_{-\tau}\left[\mathbf{x}_{0}\right]$ for some shift $\tau$.

Proof Please refer to the supplement.

This theorem says that the only local minima in this idealized case are signed shift truncations of the ground truth $\mathbf{a}_{0}$, with certain choice of $\lambda$. Moreover, on those local minima $\bar{\mathbf{a}}$, the minimizing sparse $\mathbf{x}^{\star}(\bar{\mathbf{a}})$ correspond to the soft thresholded, oppositely shifted ground truth $\mathbf{x}_{0}$. The quantity $\left\|\mathbf{W}_{I}^{*}\mathbf{W}_{I}-\mathbf{I}\right\|$ measures the orthogonality of different shifts of $\mathbf{a}_{0}$. In particular, if $\mathbf{a}_{0}$ is benign enough in the sense that any two different shifts of $\mathbf{a}_{0}$ are uncorrelated, or $\left\|\mathbf{W}_{I}^{*}\mathbf{W}_{I}-\mathbf{I}\right\|\to 0$ for any $I$, then any nonzero $\lambda$ guarantees the desired geometry. On the other hand, for a fixed $\mathbf{a}_{0}$, both $\left|\mathcal{I}\right|$ and $\left|I\right|$ becomes smaller as $\lambda$ increases, therefore the constraint $\left\|\mathbf{W}_{I}^{*}\mathbf{W}_{I}-\mathbf{I}\right\|<\lambda^{2}/6$ is more likely to be satisfied.

A similar result holds when $\mathbf{x}_{0}$ is separated enough that copies of the kernel do not overlap. It also holds if $\mathbf{x}_{0}$ is a long, sufficiently sparse random vector. For example, if the entries of $\mathbf{x}_{0}$ satisfy a Bernoulli-Gaussian distribution $\mathbf{x}_{0}(i)=\Omega(i)\mathbf{v}(i)$, with $\Omega(i)\sim\mathrm{Ber}(\theta)$ and $\mathbf{v}(i)\sim\mathcal{N}(0,1)$, and the probability $\theta$ diminishes sufficiently quickly with $k$. We conjecture that this phenomenon holds much more broadly. In particular, determining how slowly $\theta$ can diminish with $k$ is an open problem.

2.3 Global Geometry on the Simplex

In the application to image deblurring, the blur kernel is always positive and sums to $1$, which naturally leads to a simplex-constrained optimization problem. However, the optimization landscape changes drastically when the convolutional kernel $\mathbf{a}$ is constrained to live on the $\ell^{1}$ norm ball. The objective value of the same objective function over the $\ell^{1}$ ball is shown in Figure 2.

There is a significant difference induced by these two constraints: the trivial spike kernel ($\mathbf{a}=\mathbf{\delta}$) becomes the global minima and other meaningful solutions become local minima with the $\ell^{1}$ norm constraint [BVG13], while the spherical constraint always renders local minima close to some signed shift truncation of the ground truth. This important empirical knowledge of the structure of the local minima enables us to infer the ground truth from any local minimizer.

3 A Two-Stage Algorithm

Inspired by the geometric property that every local minimum of the simplified problem $\widehat{\varphi}$ is a signed shift-truncation of the ground truth $\mathbf{a}_{0}$, we present a two stage algorithm for reliable recovery of the ground truth $\mathbf{a}_{0}$ in this section. In the first stage, the algorithm recovers some signed shift truncation of the ground truth, and the following stage infers the ground truth from this partial recovery.

3.1 Stage I: Find the Signed Shift Truncation

Theorem 2.1 suggests that $\lambda$ needs to be relatively large to guarantee that all local minimizers are signed shift truncations of the ground truth. However, this is not sufficient to guarantee the success of an optimization algorithm due to the non-differentiability of the $\ell^{1}$ regularizer at $\mathbf{x}^{\star}=\mathbf{0}$. Because of this non-differentiability, when $\lambda$ is too large, there is a nonzero measure set of $\mathbf{a}$ where $\left\|\mathbf{C}^{*}_{\mathbf{a}}\mathbf{y}\right\|_{\infty}\leq\lambda$ and therefore $\mathbf{x}^{\star}(\mathbf{a})=\mathbf{0}$. These $\mathbf{a}$ are not correlated with any signed shift truncation of $\mathbf{a}_{0}$ and are the global maxima of $\varphi$. The objective function is constant over this region, so there is no way to escape using only local information.

One way to cope with this flat global maxima region is to replace the nondifferentiable $\ell^{1}$ sparsity penalty with a differentiable one. A natural choice is the huber-$\mu$ function, which can be seen as an $\ell^{1}$ penalty but with a rounded bottom for $\left|x_{i}\right|\leq\mu$:

[TABLE]

As we choose $\mu\ll\lambda$, the $\mu$-huber function closely approximates the $\ell^{1}$ norm, which still maintaining the effect of “smoothing” the flat region. The flat region for the $\ell^{1}$ penalty objective occurs when $\mathbf{x}^{*}(\mathbf{a})=\mathbf{0}$, correspondingly we define a region as $\mathcal{R}_{h,0}$ with small $\mathbf{x}^{*}(\mathbf{a})$ such that

[TABLE]

Within the $\mathcal{R}_{h,0}$ region, the original objective function can be rewritten into a simpler form:

[TABLE]

thus the optimality condition for $\mathbf{x}^{*}$ implies

[TABLE]

Plugging $\mathbf{x}^{*}(\mathbf{a})$ back to (2.6) and ignoring the higher order term $O(\frac{\mu^{2}}{\lambda^{2}})$ yields

[TABLE]

In this case, minimization of the objective function $\varphi_{h_{\mu}}(\mathbf{a})$ within the region $\mathcal{R}_{h,0}$ is equivalent to finding the maximum eigenvalue of the matrix $\mathbf{\iota}^{*}\mathbf{C}_{\mathbf{y}}^{*}\mathbf{C}_{\mathbf{y}}\mathbf{\iota}$ with the corresponding leading eigenvectors $e_{1}(\mathbf{\iota}^{*}\mathbf{C}_{\mathbf{y}}^{*}\mathbf{C}_{\mathbf{y}}\mathbf{\iota})$ achieving the local minima. However, these points can be excluded from $\mathcal{R}_{h,0}$ by setting $\lambda<\min_{\mathbf{v}\in e_{1}(\mathbf{\iota}^{*}\mathbf{C}^{*}_{\mathbf{y}}\mathbf{C}_{\mathbf{y}}\mathbf{\iota})}\left\|\mathbf{C}^{*}_{\mathbf{y}}\mathbf{\iota}\mathbf{v}\right\|_{\infty}$.777A computationally easier upper bound would be $\sqrt{\frac{\lambda_{1}(\mathbf{\iota}^{*}\mathbf{C}^{*}_{\mathbf{y}}\mathbf{C}_{\mathbf{y}}\mathbf{\iota})}{k}}$. Note that in some scenario, there exists a local minima appearing in either region $\mathcal{R}_{h,0}$ or $\mathcal{R}_{h,0}^{c}$ regardless of how we set $\lambda$. Such extreme case happens when the ground truth convolutional kernel is only supported on a small consecutive portion of its full size, hence a tight estimate of the kernel size is preferred.

With above modifications, the original flat local maxima region $\mathcal{R}_{h,0}$ becomes concave and always have a direction of negative curvature for the algorithm to escape $\mathcal{R}_{h,0}$. Hence, the first stage of the algorithm can find a signed shift-truncation of the ground truth $\bar{\mathbf{a}}=\pm\frac{\mathbf{\iota}^{*}s_{\tau}\left[\widetilde{\mathbf{a}_{0}}\right]}{\left\|\mathbf{\iota}^{*}s_{\tau}\left[\widetilde{\mathbf{a}_{0}}\right]\right\|_{F}}$ as desired.

3.2 Stage II: Infer the Ground Truth

The second stage of the algorithm aims to recover the ground truth from its signed shift truncation $\bar{\mathbf{a}}$. To recover the truncated part, we first put $\bar{\mathbf{a}}$ in a higher dimensional sphere by zero padding (Figure 3), and then recover the ground truth $\mathbf{a}_{0}$ on this higher dimensional sphere. Intuitively, as $\bar{\mathbf{a}}$ still captures a considerable portion of the ground truth $\mathbf{a}_{0}$ (the zero padded $\bar{\mathbf{a}}$ is close to the shifted $\mathbf{a}_{0}$ in a higher dimensional space), the zero padded $\bar{\mathbf{a}}$ serves as a good initialization. This intuition is made rigorous in the following lemma:

Lemma 3.1**.**

Let $\lambda_{rel}=\lambda/\left\|\mathbf{x}_{0}\right\|_{\infty}$, suppose the ground truth $\mathbf{a}_{0}$ satisfies

[TABLE]

for any nonzero shift $\tau$, and $\mathbf{x}_{0}$ is separated enough such that any two nonzero components are at least $2k$ entries away from each other. If initialized at some $\mathbf{a}\in\mathbb{S}^{k-1}$ that $\left|\left\langle\mathbf{a},\mathbf{a}_{0}\right\rangle\right|>\lambda_{rel}$, a small-stepping projected gradient method minimizing $\varphi(\mathbf{a})$ recovers the signed ground truth $\pm\mathbf{a}_{0}$.

Proof Please refer to the supplement.

This lemma says when the initial point $\mathbf{a}$ is close enough to the ground truth, the gradient always points to $\mathbf{a}_{0}$ as long as $\left|\left\langle\mathbf{a}_{0},\mathbf{\iota}s_{\tau\neq 0}\left[\widetilde{\mathbf{a}_{0}}\right]\right\rangle\right|$ is sufficiently small. Theorem 2.1 suggests that the first stage of the algorithm finds one local minimum $\mathbf{\bar{a}}$ that $\left|\left\langle\mathbf{\bar{a}},\mathbf{a}_{0}\right\rangle\right|\geq\lambda_{rel}$. Hence, the second stage of the algorithm, which minimizes the same objective function but on a higher dimensional sphere, recovers the ground truth up to sign shift ambiguity as desired.

To ensure accurate recovery, it is important to take the effect of $\lambda$ on the function geometry into consideration. A larger $\lambda$ encourages a sparser $\mathbf{x}$ and induces a simpler and smoother function landscape, which effectively eliminates undesirable local minima that are not close to any signed shift truncations, as shown in Figure 4. On the other hand, a smaller $\lambda$ emphasizes more on the accurate recovery of the signal, therefore the global minima of (2.6) will be closer to the ground truth when $\lambda$ decreases.

This geometric effect induced by $\lambda$ suggests a continuation method in the second stage of the algorithm. We start with a relatively big $\lambda$ for smoother function geometry, which encourages the algorithm to converge to one meaningful local minimum close to some signed shift truncation of the ground truth. Then run the same algorithm with decreasing sequence of $\lambda$ to produce a finer approximation of the ground truth. The overall algorithm is described in Algorithm 1.

We need to note that solving $\mathbf{a}=\arg\min\varphi_{\lambda}(\mathbf{a})$ in Algorithm 1 involves iteration between (i) finding the marginalization over $\mathbf{x}^{*}(\mathbf{a})$ step, and (ii) updating $\mathbf{a}$ based on the gradient/Hessian of $\varphi_{\lambda}(\mathbf{a})$. This could be very computationally consuming, a more efficient variant would be to optimize over the cross space of $\mathbf{a}$ and $\mathbf{x}$ together. The corresponding algorithm can be easily adapted to fit into the same general framework. The only things we want to emphasize in the proposed algorithm are the dimension lifting of the sphere and the continuation of $\lambda$.

4 Further Extensions

In this section, we extend our algorithm to handle two other deconvolution problems of practical interests: image deblurring and convolutional dictionary learning. The proposed two stage algorithm can be modified and applied to these more complicated applications.

4.1 Image Deblurring

Image deblurring aims to recover a sharp natural image from its blurred observation due to unknown photographic processes such as camera shake or defocus. Although the natural images are not necessarily sparse, it is widely acknowledged that their gradients are approximately sparse. Let $\mathbf{y}=\mathbf{a}_{0}\circledast\mathbf{x}_{0}$ denote the observed blurry image, which is the convolution of the original sharp image $\mathbf{x}_{0}$ and the blurring kernel $\mathbf{a}_{0}$. Because of the linearity of the convolution operator, the gradient of the blurred image equals the convolution of the kernel and gradient of the original sharp image, which is usually sparse as desired

[TABLE]

Here, $\nabla_{x}$ and $\nabla_{y}$ denote derivatives in the $x$ and $y$ directions. In this application, $\nabla_{x}\mathbf{x}_{0}$ and $\nabla_{y}\mathbf{x}_{0}$ are the underlying sparse signals, and the blind image deblurring problem can be cast as solving:

[TABLE]

Here, $\mathbb{S}_{+}^{k-1}$ denotes the intersection of the unit sphere and the positive orthant. In this application, the non-negativity of the blurring kernel removes the sign ambiguity. We observe in experiments that local minimizers are all near some shift truncation of the ground truth kernel. The same two stage algorithm can therefore be applied to infer the ground truth.

4.2 Convolutional Dictionary Learning

Convolutional dictionary learning (CDL) is an important problem in machine learning for images, speeches, as well as scientific problems like microscopy data analysis and neural spike sorting. The observation signal $\mathbf{y}$ is the superposition of convolutions of $N$ pairs of kernels $\mathbf{a}_{0n}$ and corresponding coefficients $\mathbf{x}_{0n}$:

[TABLE]

Blind deconvolution can be seen as a special case of CDL with $N=1$. If the coefficients $\mathbf{x}_{0n}$ are sparse, a natural way to extend our knowledge of SBD would be to assume all $N$ convolutional kernels having unit Frobenius norm and cast it as minimizing following objective function over the product of $N$ spheres:

[TABLE]

We anticipate that all the local minima are near signed shift truncations of the ground truth, provided the target kernels $\mathbf{a}_{0n}$ are sufficiently diverse. The modified two stage algorithm still manages to capture the partial information offered by local minima and hence recovers the ground truth. Experimental results are provided in Section 5.4 to corroborate this claim.

5 Experiments

In this section, we investigate the performance of our algorithm on both synthetic and real data. We first report a systematic investigation, performed in [CLC*+*17], of the performance of our algorithm on synthetic data, which are designed to mimic properties of the microscopy data analysis problem. In Sections 5.2-5.4, we present experiment results showing how our method performs on real data from microscopy and image deblurring.

5.1 Evaluation on Synthetic Data

Noise-free data: we generate the noise-free observation signal of size $m=256\times 256$ through circular convolution between a kernel of size $k$ and a random underlying activation signal with a Bernoulli distribution with sparsity $\theta$, i.e. $x_{i}\overset{\text{i.i.d.}}{\sim}\text{Ber}(\theta)$, or $x\sim\text{Ber}(\theta)$. We plot the kernel recovery error for varying kernel size $k$ and sparsity level $\theta$ in the left of Figure 5 [CLC*+*17]. Each point on the diagram is the average of 20 independent measurements. The algorithm performs excellently in the blue regions, but begins to fail in the red regions, where either the kernel size is large or the underlying activation signal is dense. The region where typical STM measurements are performed are bounded below by the white dashed line, where the proposed algorithm achieves satisfying performance.

Noisy data: we generate convolutional signals by convolving fixed kernel of dimension $k$, $k/m=0.14$ with the random activation map $x\sim\text{Ber}(\theta)$ of dimension $m$, and applying additive Gaussian noise. We test the performance of our algorithm for varying sparsity $\theta$ and noise power. The result is shown in Figure 5 (right): the algorithm achieves noise-robust recovery when the sparsity constraint is satisfied.

5.2 Microscopy Data Analysis

We apply our algorithm on experimental microscopy data obtained from a NaFeCoAs sample. Our results shown in Figure 6 indicate that the proposed algorithm manages to recover the missing details of the ripples in the Fourier domain of the defect, which encode the physical scattering processes of electrons at work.

5.3 Image Deblurring

We test our algorithm on the image deblurring dataset from [LWDF11], solving (4.1) to recover the convolutional kernel. To clearly separate the inaccuracy of the algorithm and the universal blurring kernel model, all the experiments are done on three kinds of blurred images: (i) synthetic blurred images generated by the convolution of sharp images and blurry kernels; (ii) noisy blurred images generated by adding Gaussian noise to the clean synthetic blurred images (SNR=100); and (iii) real blurry images taken with camera shakes [LWDF11].

We compare with algorithms by Zhang et al.[ZWZ13], Krishnan et al.[KTF11], Sun et al.[SCWH13], and Liu et al.[LCM14].888We use the default parameters for these algorithms. It’s possible that better performance could be obtained by tuning the parameters more carefully. In our algorithm, we fix the $\lambda$’s to be $0.1,0.01,0.001,0.001$ for all the instances.

Because of the shift ambiguity, we evaluate the accuracy of the recovered blurring kernel considering all possible shifts. The kernel recovery error is defined as $\min_{\tau}\left\|\mathbf{\iota}^{*}s_{\tau}\left[\widetilde{\mathbf{a}}\right]/\left\|\mathbf{a}\right\|_{1}-\mathbf{a}_{0}/\left\|\mathbf{a}_{0}\right\|_{1}\right\|_{F}$, and the cumulative distribution is shown in Figure 7.

We use the same non-blind deblurring algorithm from [KF09], with the same parameter. We consider the blurred image using the ground truth kernel to be the bench mark, and evaluate the quality of the deblurred image by calculating the Frobenius norm of its difference to such bench mark. Results are shown are in Figure 8.

Our algorithm achieves better convolutional kernel recovery for all three types of images, but its improvement on deblurred image is less obvious, especially for real images. This could be due to (i) the convolutional kernel in this dataset is not strictly uniform across the image, and (ii) the non-blind deconvolution algorithm exploits the heavy-tailed distribution of a natural image’s gradient and becomes less sensitive to the accuracy of the recovered convolutional kernel.

5.4 Convolutional Dictionary Learning

We show results of recovering multiple convolutional kernels on both synthetic data (Figure 9) and real STM data (Figure 10). In the synthetic data, the three convolutional kernels are of size $16\times 16$ and their corresponding activation signals are generated through a Bernoulli model of sparsity $0.005$. Results of both stages of the algorithm are shown in Figure 9: the first stage returns kernels close to some shift truncations of the ground truth, and the second stage recovers the ground truth on a higher dimensional space.

We repeat this experiment with microscopy data obtained from a NaFeAs sample. The algorithm manages to differentiate the two convolutional kernels (defect patterns), as shown in Figure 10. For this material, the kernel orientations depend on the history of the material (stress, temperature, etc.), and using convolutional dictionary learning can be used to automatically detect these features.

6 Generalizations: Matching Loss and Constraints

In the Lasso-like objective function studied in this paper, the approximation error of $\mathbf{y}-\mathbf{a}\circledast\mathbf{x}$ are measured in the squared Frobenius norm (or entry-wise $\ell^{2}$ norm), which is usually adopted to penalize gaussian addictive noise in the observation. To obtain reliable solutions, a spherical constraint on the short kernel turned out to be crucial. Moreover, several theoretical papers [ZKW18, KZLW18] studying variations of the Lasso-like objective function also adopt the unit Frobenius constraint, and show that the claimed geometry continues to hold under much more general conditions. However, noise could be more complicated and does not always satisfy the gaussian model in real applications. If given prior knowledge of the noise, other penalties for the loss term $\mathbf{y}-\mathbf{a}\circledast\mathbf{x}$ could be preferable. For example, $\|\mathbf{y}-\mathbf{a}\circledast\mathbf{x}\|_{p}^{p}$ with $p<2$ performs better to penalize heavy tailed noise, and $\|\mathbf{y}-\mathbf{a}\circledast\mathbf{x}\|_{p}^{p}$ with $p>2$ performs better control on the magnitude of the noise.

In this section, we briefly discuss the following more general formulation for the short-and-sparse deconvolution problem

[TABLE]

Similarly, we write $\varphi_{p}(\mathbf{a})\doteq\min_{\mathbf{x}}\psi_{p}(\mathbf{a},\mathbf{x})$. In this case, we use entrywise $\ell^{p}$ norm to measure the approximation error of $\mathbf{a}\circledast\mathbf{x}$, and assume the optimization constraint for $\mathbf{a}$ to be a unit $\ell^{q}$ ($2<q<\infty$) sphere. Comprehensive study of the function landscape of above general formulation is even more challenging. However, we can demonstrate that once $p=q$, there exist local solutions that share the same geometric property (a local optimum is close to some scaled shifted truncation) through local analysis. As in the objective optimization problem (2.6), the inherent shift ambiguity leads to these structured local minimizers. Because of the randomness in $\mathbf{x}_{0}$, the derived results hold for more general $\mathbf{x}_{0}$.

Lemma 6.1**.**

Suppose $\mathbf{y}=\mathbf{a}_{0}\circledast\mathbf{x}_{0}$ with $\mathbf{x}_{0}=\mathbf{e}_{0}$, and $p=q\geq 2$. Then for any shift $\tau$, positive scalar $\varepsilon$ and $\lambda$ such that at every $\mathbf{a}\in\mathbb{S}_{q}\cap\mathbb{B}\left(\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}},\varepsilon\right)$ the solution $\mathbf{x}^{*}(\mathbf{a})\doteq\arg\min_{\mathbf{x}}\psi_{p}(\mathbf{a},\mathbf{x})$ is (i) unique and (ii) is supported on the $\tau$-th entry, the point $\bar{\mathbf{a}}\doteq\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}}$ is a strict local minimizer of the cost $\varphi_{p}(\mathbf{a})$ over the manifold $\left\|\mathbf{a}\right\|_{p}=1$.

Proof Please refer to the supplement.

This lemma says that as long as the constraint matches the loss, i.e., we choose $q=p\geq 2$, a scaled shift truncation $\bar{\mathbf{a}}=\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]/\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}$ achieves the local minimum. In contrast, when the constraint does not match the loss, i.e., $q\neq p$, this property is not satisfied: the Riemannian gradient vanishes at a stationary point $\mathbf{a}$, which satisfies

[TABLE]

and

[TABLE]

with $\alpha^{\prime}$ denoting another non-zero scalar of arbitrary value. This point is not a shift truncation of the ground truth.999When $p\neq q$, with $\lambda$ decreasing, (6.2) can be rewritten as $\alpha\mathbf{a}=\mathbf{\iota}^{*}s_{-\tau}\left[\widetilde{\mathbf{a}_{0}}\right]+\alpha^{\prime}\mathbf{a}^{\circ\frac{q-1}{p-1}}$. With the stationary condition for $\mathbf{x}^{*}(\mathbf{a})$, we have $\alpha\mathbf{a}\to\mathbf{\iota}^{*}s_{-\tau}\left[\widetilde{\mathbf{a}_{0}}\right]$ as $\lambda\to 0$. Algorithmically, this helps to explain why practical algorithms solving (6.1) with $p\neq q$ can sometimes recover the kernel, where a careful continuation in $\lambda$ is always necessary.

Hence, when the constraint matches the loss, there is good local geometry. In this paper, we have empirically demonstrated that when $p=q=2$, there is good global geometry, under appropriate conditions. We can potentially leverage this property to find global solutions to the general $\ell^{p}$ ($p\geq 2$) problem: one first minimizes $\ell^{2}$ to obtain a point near the ground truth and then locally minimizes $\varphi_{p}$ over an $\ell^{p}$ ball to obtain an estimate that uses the statistical characteristics of the noise.

Although our analysis only pertains to $q\geq 2$, we believe that there is good local geometry even when $q<2$. In this situation, the constraint set is no longer twice differentiable and the shift truncations are always on the non-smooth region of the constraint set. As the $p$ decrease, the normal cone at the nonsmooth points gets wider. For example, if $q\geq 2$ the normal cone is the normal direction. Therefore, as $q$ decreases, the normal cone is larger and hence more like to contain the gradient direction, then this shift truncation is likely to be a stationary point.

7 Discussions

This work studies the global geometry of a nonconvex optimization problem for SBD when the kernel is assumed to have unit Frobenius norm. In this setting, we find that all the local minima are benign, in the sense that they are close to some signed shift truncation of the ground truth. With this insight, we propose a two stage algorithm that recovers the ground truth by exploiting the information hidden in local minima.

This problem reveals the challenges faced when analyzing the SBD problem via a geometrical approach. For problems enjoying stronger symmetry properties [SQW17a, SQW17b], similar approaches yield a global understanding of the function geometry and recovery guarantees. We expect that the weak symmetries in SBD make a major contribution to the difficulties encountered for this problem.

There are lots of additional further directions could be of great interests for both theory and application: Our empirical results show that the our characterization of local minima carries through to the convolutional dictionary learning problem, which can also be efficiently solved by slight adaptation of the proposed algorithm. However, the theory part is open, it would be interesting to know how many kinds of kernels, or what kinds of kernels are recoverable, probably by some measures of incoherence, which is a common assumption in dictionary learning problem.

Two other imperfections we encounter in scientific measurement are resolution limit and measurement error, which inspire us to consider (i) if it’s possible to integrate blind deconvolution and super-resolution process together; (ii) if we can come up with a robust blind deconvolution algorithm to automatically rule out noisy entries.

Acknowledgement

The authors gratefully acknowledge support from NSF 1343282, NSF CCF 1527809, and NSF IIS 1546411.

Appendix A Notations

This supplementary material contains a complete proof of Theorem 2.1 and Lemma 3.1 of the main paper. Section 1 discusses some technical notations and geometric observations, Section 2 and 3 presents a complete proof for the theorem and lemma respectively, and Section 4 contains two auxiliary lemmas used in the main proof.

Before commencing our main proof, we introduce some notations, recap the simplified objective function $\widehat{\varphi}$ and a few basic observations on its geometry. In this note, we use $\mathbf{C}_{\mathbf{a}}$ or $\mathbf{C}_{\widetilde{\mathbf{a}}}$ exchangeably to denote the $m\times m$ circulant matrix generated by a $k$-length short convolutional kernel $\mathbf{a}$ without ambiguity.

Reversed Circulant Matrix Let $\check{\mathbf{C}}_{\mathbf{a}_{0}}$ be the reversed circulant matrix for $\mathbf{a}_{0}$ such that

[TABLE]

Here, $s_{\tau}\left[\widetilde{\mathbf{a}_{0}}\right]$ denotes a cyclic shift as defined in Equation (3) of the main paper. Since the $i$-th entry of $\mathbf{C}_{\mathbf{a}}^{*}\widetilde{\mathbf{a}_{0}}$ satisfies

[TABLE]

therefore following equation always holds

[TABLE]

Projection onto $I$ Let $\mathbf{e}_{0},\dots,\mathbf{e}_{m-1}$ denote the standard basis vectors. For index set $I=\left\{i_{1},\dots,i_{|I|}\right\}\neq\emptyset$, let

[TABLE]

and then

[TABLE]

denote the projection onto the coordinates indexed by $I$ while setting other entries zero.

Partial Signed Support Let $\mathbf{u}$, $\mathbf{v}$ be two vectors of the same dimension. If $\mathrm{supp}\left(\mathbf{u}\right)\subseteq\mathrm{supp}\left(\mathbf{v}\right)$ and $\mathbf{u}(i)\mathbf{v}(i)\geq 0$ for all $i\in\mathrm{supp}\left(\mathbf{v}\right)$, then $\mathbf{u}$ attains the partial signed support of $\mathbf{v}$, denoted as $\mathbf{u}\unlhd\mathbf{v}$.

Piecewise Quadratic Function Let $\mathbf{x}^{*}(\mathbf{a})$ denote the minimizer for simplified objective function

[TABLE]

with sign $\mathbf{\sigma}$ and support $I$ defined as

[TABLE]

By stationary condition for $\mathbf{x}^{*}(\mathbf{a})$, we obtain

[TABLE]

where $\mathrm{SOFT}_{\lambda}\left[u\right]=\operatorname{sign}(u)\max\left\{|u|-\lambda,0\right\}$ is the entry-wise soft-thresholding operator.

For each sign pattern $\mathbf{\sigma}=\mathrm{supp}\left(\mathbf{x}^{*}\right)\in\left\{-1,0,1\right\}^{m}$, there exists corresponding region on the sphere such that

[TABLE]

On the relative interior of each $R_{\mathbf{\sigma}}$, the function $\widehat{\varphi}$ has a simple expression:

[TABLE]

Therefore, the objective function is piecewise quadratic and can be rewritten as

[TABLE]

with

[TABLE]

With above notations clarified, we are now ready to present a proof for the main theorem.

Appendix B Proof of Theorem 2.1

Proof On the relative interior of each $R_{\mathbf{\sigma}}$, the simplified objective function

[TABLE]

has Euclidean derivative and Hessian

[TABLE]

As we assume $\mathbf{a}$ to have unit Frobenius norm, or to live on a sphere, the more natural Riemannian gradient and Hessian are defined as

[TABLE]

Here, $\mathbf{P}_{\mathbf{a}^{\perp}}=\mathbf{I}-\mathbf{a}\mathbf{a}^{*}$ denotes projection onto the tangent space over the sphere at $\mathbf{a}$. As in the Euclidean space, a stationary point on the sphere needs to satisfy $\mathrm{grad}\left[\widehat{\varphi}_{\mathbf{\sigma}}\right]\left(\mathbf{a}\right)=\mathbf{0}$. At a stationary point $\bar{\mathbf{a}}$, if $\mathrm{Hess}\left[\widehat{\varphi}_{\mathbf{\sigma}}\right]\left(\bar{\mathbf{a}}\right)$ is positive semidefinite, the function is convex and $\bar{\mathbf{a}}$ is a local minimum; if $\mathrm{Hess}\left[\widehat{\varphi}_{\mathbf{\sigma}}\right]\left(\bar{\mathbf{a}}\right)$ has a negative eigenvalue, then there exists a direction alone which the objective value decreases and hence $\bar{\mathbf{a}}$ is a saddle point [AMS07].

Let $I=\left\{i_{1}<i_{2}<\dots<i_{|I|}\right\}$ and define

[TABLE]

and

[TABLE]

Here, columns of $\mathbf{U}$ have unit $\ell^{2}$ norm. Then we have

[TABLE]

As $I$ is defined via soft thresholding, we have $\left|[\check{\mathbf{C}}_{\mathbf{a}_{0}}^{*}\mathbf{\iota}\mathbf{a}]_{i}\right|>\lambda$ holds for every $i\in I$. Hence,

[TABLE]

holds at any $\mathbf{a}\in\mathrm{cl}\left(R_{\mathbf{\sigma}}\right)\setminus R_{\mathbf{0}}$.

Stationary point and implications

Consider any stationary point $\bar{\mathbf{a}}\in\mathrm{cl}\left(R_{\mathbf{\sigma}}\right)\setminus R_{\mathbf{0}}$ of $\widehat{\varphi}$. By continuity of the gradient of $\widehat{\varphi}$ (proved in Lemma B.1), $\bar{\mathbf{a}}$ is also a stationary point of $\widehat{\varphi}_{\mathbf{\sigma}}$. By definition, $\mathrm{grad}\left[\widehat{\varphi}_{\mathbf{\sigma}}\right]\left(\bar{\mathbf{a}}\right)=\mathbf{0}$, which implies that

[TABLE]

Note that since $\bar{\mathbf{a}}^{*}\mathbf{M}_{\mathbf{\sigma}}\bar{\mathbf{a}}+\mathbf{b}_{\mathbf{\sigma}}^{*}\bar{\mathbf{a}}\neq 0$ and $\mathbf{b}_{\mathbf{\sigma}}\in\mathrm{range}(\mathbf{M}_{\mathbf{\sigma}})$, this implies that $\bar{\mathbf{a}}\in\mathrm{range}(\mathbf{M}_{\mathbf{\sigma}})$.

Let $\gamma=-\left(\mathbf{a}^{*}\mathbf{M}_{\mathbf{\sigma}}\mathbf{a}+\mathbf{b}_{\mathbf{\sigma}}^{*}\mathbf{a}\right)>0$, then the condition for a stationary point $\bar{\mathbf{a}}$ becomes

[TABLE]

Let $\mathbf{\alpha}=\mathbf{U}^{*}\bar{\mathbf{a}}$, and note that for each $j$, $\alpha_{j}>0$ and $\alpha_{j}\eta_{j}>\lambda$. In terms of $\mathbf{U}$, the stationarity condition becomes

[TABLE]

Since the diagonal elements of $\mathbf{U}^{*}\mathbf{U}$ are all ones, and hence can be written as

[TABLE]

We have

[TABLE]

As $\mathbf{\alpha}\succ\mathbf{0}$ and $\mathrm{diag}(\mathbf{\eta})\mathbf{\alpha}\succ\lambda\cdot\mathbf{1}$111111Here, $\succ$ denotes element-wise inequality between vectors., together with an auxiliary Lemma E.2 proved in Section E, we have

[TABLE]

whence

[TABLE]

Suppose that $\mathbf{\Delta}$ is small enough that the right hand side of (B.24) is bounded by $\lambda^{2}/2$, i.e.,

[TABLE]

Plugging back into the stationary condition \operatorname{diag}(\mathbf{\eta})^{2}\mathbf{\alpha}-\gamma\mathbf{\alpha}=\lambda\mathbf{\eta}-\mathbf{\Delta}\Bigl{\{}\operatorname{diag}(\mathbf{\eta})^{2}\mathbf{\alpha}-\lambda\mathbf{\eta}\Bigr{\}} gives

[TABLE]

Since $\alpha_{i}<1$ and $\eta_{i}>\lambda$ for all $i$, which implies that

[TABLE]

where $\eta_{\min}$ is the smallest of the $\eta_{i}$.

Negative curvature in Hessian Recall the Riemannian Hessian on the sphere is defined as

[TABLE]

Below argument shows that this Riemmanian Hessian has negative eigenvalues. Let $\tilde{\mathbf{U}}$ be an orthonormal matrix generated via

[TABLE]

Whenever $\left\|\mathbf{\Delta}\right\|_{\ell^{2}\to\ell^{2}}<1/2$ holds, Lemma E.2 guarantees $\left\|\mathbf{U}-\tilde{\mathbf{U}}\right\|_{\ell^{2}\to\ell^{2}}<3\left\|\mathbf{\Delta}\right\|_{\ell^{2}\to\ell^{2}}$. Under this condition, we can lower bound the smallest nonzero eigenvalue of $\mathbf{U}\operatorname{diag}(\mathbf{\eta})^{2}\mathbf{U}^{*}$, as

[TABLE]

Since $\lambda<\eta_{\min}\leq\eta_{\max}\leq 1$, additionally if $\left\|\mathbf{\Delta}\right\|_{\ell^{2}\to\ell^{2}}\leq\frac{\lambda}{6}$, we have

[TABLE]

Together with (B.25) and (B.27), we can obtain

[TABLE]

or

[TABLE]

Thus, whenever the following conditions are satisfied

[TABLE]

we have $\lambda_{\max}(\mathbf{M}_{\mathbf{\sigma}})<-\gamma$ as desired.

Above calculations imply that for every $\mathbf{\xi}\in\mathrm{range}\left(\mathbf{M}_{\mathbf{\sigma}}\right)\subseteq\mathbb{R}^{k}$,

[TABLE]

Since

[TABLE]

for $\mathbf{\xi}\in\bar{\mathbf{a}}^{\perp}\cap\mathrm{range}\left(\mathbf{M}_{\mathbf{\sigma}}\right)$,

[TABLE]

Hence, on the relative interior $\mathrm{relint}\left(R_{\mathbf{\sigma}}\right)$, $\widehat{\varphi}\equiv\widehat{\varphi}_{\mathbf{\sigma}}$ obtains, and so this implies that for $\left\|\mathbf{\sigma}\right\|_{0}>1$, there are no local minima in $\mathrm{relint}\left(R_{\mathbf{\sigma}}\right)$.

Relative boundaries We first note that if $\left\|\mathbf{\sigma}\right\|_{0}=1$ and $I=\left\{i\right\}$, either $R_{\mathbf{\sigma}}$ is empty when $\left\|\mathbf{\iota}^{*}s_{-i}\left[\widetilde{\mathbf{a}_{0}}\right]\right\|_{2}\leq\lambda$, or it contains an open ball around $\mathrm{range}\left(\mathbf{M}_{\mathbf{\sigma}}\right)\bigcap\mathbb{S}^{k-1}=\pm\frac{\mathbf{\iota}^{*}s_{-i}\left[\widetilde{\mathbf{a}_{0}}\right]}{\left\|\mathbf{\iota}^{*}s_{-i}\left[\widetilde{\mathbf{a}_{0}}\right]\right\|_{2}}$. Hence, if $\bar{\mathbf{a}}\in\mathrm{relbdy}\left(R_{\mathbf{\sigma}}\right)$ is a stationary point and $\mathbf{\sigma}\neq\mathbf{0}$, we necessarily have $\left\|\mathbf{\sigma}\right\|_{0}\geq 2$.

Since $\bar{\mathbf{a}}$ is on the boundary of $R_{\mathbf{\sigma}}$, it is also in $\mathrm{relbdy}\left(\mathrm{cl}\left(R_{\mathbf{\sigma}^{\prime}}\right)\right)$ for some $\mathbf{\sigma}^{\prime}\neq\mathbf{\sigma}$. Let

[TABLE]

Suppose that for every $\mathbf{\sigma}^{\prime}\in\Xi$, $\mathbf{\sigma}\unlhd\mathbf{\sigma}^{\prime}$. Hence, $\mathrm{range}\left(\mathbf{M}_{\mathbf{\sigma}}\right)\subseteq\mathrm{range}\left(\mathbf{M}_{\mathbf{\sigma}^{\prime}}\right)$ for every $\mathbf{\sigma}^{\prime}\in\Xi$ and

[TABLE]

By continuity of the gradients, $\bar{\mathbf{a}}$ is a stationary point for every $\widehat{\varphi}_{\mathbf{\sigma}^{\prime}}$ such that $\mathbf{\sigma}^{\prime}\in\Xi$. If we choose an arbitrary nonzero $\mathbf{\xi}\in\mathrm{range}(\mathbf{M}_{\mathbf{\sigma}})$, we have that for every $\mathbf{\sigma}^{\prime}\in\Xi$,

[TABLE]

There exists a neighborhood $N$ of $\bar{\mathbf{a}}$ for which, at every $\mathbf{a}\in N\cap R_{\mathbf{\sigma}^{\prime}}$, $\widehat{\varphi}(\mathbf{a})=\widehat{\varphi}_{\mathbf{\sigma}^{\prime}}(\mathbf{a})\leq\widehat{\varphi}(\bar{\mathbf{a}})$ for some $\mathbf{\sigma}^{\prime}\in\Xi$. Hence, $\bar{\mathbf{a}}$ is not a local minimum of $\widehat{\varphi}$.

Local minima If $\left\|\mathbf{\sigma}\right\|_{0}=1$ and $I=\left\{i\right\}$, then the simplified objective function is

[TABLE]

The minimizer appears at the boundary for $\left\langle\sigma_{i}\mathbf{\iota}^{*}s_{-i}\left[\widetilde{\mathbf{a}_{0}}\right],\mathbf{a}\right\rangle$, namely $\left\|\mathbf{\iota}^{*}s_{-i}\left[\widetilde{\mathbf{a}_{0}}\right]\right\|_{2}$121212The other boundary point is $\left\langle\sigma_{i}\mathbf{\iota}^{*}s_{-i}\left[\widetilde{\mathbf{a}_{0}}\right],\mathbf{a}\right\rangle=\lambda$, which achieves a smaller objective value. obtained by

[TABLE]

It can be easily verified that

[TABLE]

Global maxima If $\left\|\mathbf{\sigma}\right\|_{0}=0$, then the objective remains constant and achieves the global maximum.

Lemma B.1** (Continuity of the Gradient of $\widehat{\varphi}$).**

$\nabla\widehat{\varphi}$* is a continuous function of $\mathbf{a}$.*

Proof Recall that for a given $\mathbf{\sigma}$, the gradient

[TABLE]

This is a continuous function within the relative interior of $R_{\mathbf{\sigma}}$. Next, we show this function is continuous at the relative boundary of $R_{\mathbf{\sigma}}$. Let $\mathbf{a}^{\prime}\in\mathrm{relbdy}\left(R_{\mathbf{\sigma}}\right)$, and $\mathbf{\sigma}^{\prime}=\operatorname{sign}(\mathbf{a}^{\prime})$, $I=\mathrm{supp}\left((\right)\mathbf{\sigma}^{\prime})$ are the corresponding sign and support. Without loss of generality, we assume $\mathbf{\sigma}^{\prime}\unlhd\mathbf{\sigma}$, denote $\mathbf{a}=\mathbf{a}^{\prime}+\varepsilon\mathbf{\delta}$ ($\left\|\mathbf{\delta}\right\|_{2}=1$) and $I_{\mathbf{\delta}}=I\setminus I^{\prime}$, then

[TABLE]

Since $\left\|\mathbf{P}_{I_{\delta}}\left(\check{\mathbf{C}}_{\mathbf{a}_{0}}^{*}\mathbf{\iota}\mathbf{a}-\lambda\mathbf{\sigma}\right)\right\|_{\infty}=\varepsilon\left\|\mathbf{P}_{I_{\delta}}\check{\mathbf{C}}_{\mathbf{a}_{0}}^{*}\mathbf{\iota}\mathbf{\delta}\right\|_{\infty}$, we have $\left\|\nabla\widehat{\varphi}_{\mathbf{\sigma}}(\mathbf{a})-\nabla\widehat{\varphi}_{\mathbf{\sigma}^{\prime}}(\mathbf{a}^{\prime})\right\|_{\infty}\leq\mathcal{O}(\varepsilon)$.

Appendix C Proof of Lemma 3.1

Lemma C.1**.**

Let $\lambda_{rel}=\lambda/\left\|\mathbf{x}_{0}\right\|_{\infty}$, suppose the ground truth $\mathbf{a}_{0}$ satisfies

[TABLE]

for any nonzero shift $\tau$, and $\mathbf{x}_{0}$ is separated enough such that any two nonzero components are at least $2k$ entries away from each other. If initialized at some $\mathbf{a}\in\mathbb{S}^{k-1}$ that $\left|\left\langle\mathbf{a},\mathbf{a}_{0}\right\rangle\right|>\lambda/\left\|\mathbf{x}_{0}\right\|_{\infty}$, a small-stepping projected gradient method minimizing $\varphi(\mathbf{a})$ recovers the signed ground truth $\pm\mathbf{a}_{0}$.

Proof Without loss of generality, we are going to assume $\left\|\mathbf{x}_{0}\right\|_{\infty}=1$ for simplicity. Given that

[TABLE]

and $\mathbf{a}=\left\langle\mathbf{a},\mathbf{a}_{0}\right\rangle\mathbf{a}_{0}+\mathbf{\delta}$ with $\left\|\mathbf{\delta}\right\|_{2}\leq\sqrt{1-\lambda^{2}}$, therefore

[TABLE]

Moreover, as $\mathbf{x}_{0}$ is sufficiently separated, we have

[TABLE]

Hence, there exists a unique nonzero minimizer satisfying $\mathrm{supp}\left(\mathbf{x}^{\star}\right)\subset\mathrm{supp}\left(\mathbf{x}_{0}\right)$, and the optimality condition for $\mathbf{x}^{\star}$ implies

[TABLE]

In this case, we can calculate the Euclidean gradient

[TABLE]

and the Riemannian gradient

[TABLE]

It is easy to check that at any point along the geodesic curve between $\mathbf{a}_{0}$ and $\mathbf{a}$, support recovery of $\mathbf{x}^{\star}$ is achieved. A small-stepping gradient descent algorithm moves towards the signed ground truth $\pm\mathbf{a}_{0}$, as desired.

Appendix D Proof of Lemma 6.1

Lemma D.1**.**

Suppose $\mathbf{y}=\mathbf{a}_{0}\circledast\mathbf{x}_{0}$ with $\mathbf{x}_{0}=\mathbf{e}_{0}$, and $p=q\geq 2$. Then for any shift $\tau$, positive scalar $\varepsilon$ and $\lambda$ such that at every $\mathbf{a}\in\mathbb{S}_{q}\cap\mathbb{B}\left(\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}},\varepsilon\right)$ the solution $\mathbf{x}^{*}(\mathbf{a})\doteq\arg\min_{\mathbf{x}}\psi_{p}(\mathbf{a},\mathbf{x})$ is (i) unique and (ii) supported on the $\tau$-th entry, the point $\bar{\mathbf{a}}\doteq\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}}$ is a strict local minimizer of the cost $\varphi_{p}(\mathbf{a})$ over the manifold $\left\|\mathbf{a}\right\|_{q}=1$.

Proof Write $\mathbf{x}^{*}(\mathbf{a})=\alpha\mathbf{e}_{\tau}$ at point $\mathbf{a}\in\mathbb{S}_{p}\cap\mathbb{B}\left(\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}},\varepsilon\right)$, and suppose without loss of generality that $\alpha>0$. Then using the KKT conditions, we have

[TABLE]

Magnitude of $\mathbf{x}^{*}$

At point $\bar{\mathbf{a}}=\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}}$, there exists closed form solution

[TABLE]

As the shifted support recovery of $\mathbf{x}^{*}(\mathbf{a})$ holds at point $\mathbf{a}\in\mathbb{B}(\bar{\mathbf{a}},\varepsilon)$, with implicit function theorem, we can obtain that $\alpha$ is a differentiable function of $\mathbf{a}$, and

[TABLE]

Euclidean Gradient and Hessian

The Euclidean gradient of $\varphi_{p}$ at $\bar{\mathbf{a}}$ is

[TABLE]

And the Euclidean Hessian of $\varphi_{p}$ at $\bar{\mathbf{a}}$ is

[TABLE]

Riemannian Gradient and Hessian

We call the optimization domain manifold $\mathcal{M}_{q}$ where $\left\|\mathbf{a}\right\|_{q}=1$ with $2\leq q<\infty$, the Riemannian gradient at $\mathbf{a}$ can be written as

[TABLE]

with

[TABLE]

Therefore, if $p=q$, $\bar{\mathbf{a}}$ is a stationary point with vanishing Riemannian gradient

[TABLE]

The Riemannian Hessian $\mathrm{Hess}[\varphi_{p}](\mathbf{a}):T_{\mathbf{a}}\mathcal{M}_{q}\to T_{\mathbf{a}}\mathcal{M}_{q}$ is defined via the formula

[TABLE]

where in the second line we have used Proposition 5.3.2 of the [AMS09]. Here, $\tilde{\nabla}$ is the Riemannian connection on $\mathcal{M}_{q}$, $\nabla$ is the standard Euclidean connection (directional derivative) on $\mathbb{R}^{k}$ and $\mathrm{grad}[\varphi_{p}]$ is any smooth extension of the Riemannian gradient to $\mathbb{R}^{k}$. Plugging in and using the definition of Euclidean gradient $\nabla$, we have

[TABLE]

For simplicity, we write

[TABLE]

and

[TABLE]

Notice that the vector-valued function $\mathbf{\xi}(\mathbf{a})$ is differentiable away from $\mathbf{a}=\mathbf{0}$, and that

[TABLE]

Since that $\mathbb{P}_{\mathbf{a}^{\perp}}\mathbf{\xi}(\mathbf{a})=\mathbf{0}$, we have

[TABLE]

Calculus gives

[TABLE]

hence

[TABLE]

At last, we can obtain the expression

[TABLE]

for $\mathbf{\delta}\in T_{\mathbf{a}}\mathcal{M}_{q}$. Noting that for any such $\mathbf{\delta}$, $\mathbb{P}_{\mathbf{a}^{\perp}}\mathbf{\delta}=\mathbf{\delta}$, we can identify the Riemannian Hessian with the matrix

[TABLE]

Local Convexity

Plugging in $\bar{\mathbf{a}}=\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}}$ with $p=q$, we have

[TABLE]

Hence, the function is strictly convex at point $\bar{\mathbf{a}}=\frac{\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]}{\left\|\mathbf{\iota}^{*}s_{-\tau}\left[\mathbf{a}_{0}\right]\right\|_{q}}$ along direction $\mathbf{\delta}\in T_{\mathbf{a}}\mathcal{M}_{q}$ as long as whose support are contained in the support of $\mathbf{a}$. We further consider a perturbation $\mathbf{\delta}_{c}$ such that $\mathrm{supp}\left(\mathbf{\delta}_{c}\right)\cap\mathrm{supp}\left(\bar{\mathbf{a}}\right)=\emptyset$, and $\left\|\mathbf{\delta}_{c}\right\|_{2}\leq\varepsilon$, then the objective value at $\frac{\bar{\mathbf{a}}+\mathbf{\delta}_{c}}{\left\|\bar{\mathbf{a}}+\mathbf{\delta}_{c}\right\|_{q}}$ satisfies

[TABLE]

This completes the proof.

Appendix E Auxiliary Lemmas

Lemma E.1**.**

*(Lemma B.2 of [SQW17a])

Suppose that $\mathbf{A}\succ\mathbf{0}$ is a positive definite matrix. For any symmetic matrix $\mathbf{\Delta}$ with $\left\|\mathbf{\Delta}\right\|_{\ell^{2}\to\ell^{2}}\leq\sigma_{\min}(\mathbf{A})/2$,*

[TABLE]

Lemma E.2**.**

Let $\mathbf{U}$ be a matrix such that $\mathbf{U}^{*}\mathbf{U}=\mathbf{I}+\mathbf{\Delta}$, with $\left\|\mathbf{\Delta}\right\|_{\ell^{2}\to\ell^{2}}<1/2$. Then $\mathbf{U}^{*}\mathbf{U}$ is invertible,

[TABLE]

and

[TABLE]

Proof Upper bound for the first quantity can be derived

[TABLE]

Applying Lemma E.1 for the second term

[TABLE]

Hence, we can obtain the claim by using $2\sqrt{3/2}<3$ to simplify the constant.