Title: SHORT AND SPARSE DECONVOLUTION -A GEOMET-RIC APPROACH

Abstract: Short-and-sparse deconvolution (SaSD) is the problem of extracting localized, recurring motifs in signals with spatial or temporal structure. Variants of this problem arise in applications such as image deblurring, microscopy, neural spike sorting, and more. The problem is challenging in both theory and practice, as natural optimization formulations are nonconvex. Moreover, practical deconvolution problems involve smooth motifs (kernels) whose spectra decay rapidly, resulting in poor conditioning and numerical challenges. This paper is motivated by recent theoretical advances (Zhang et al., 2017; Kuo et al., 2019) , which characterize the optimization landscape of a particular nonconvex formulation of SaSD and give a provable algorithm which exactly solves certain non-practical instances of the SaSD problem. We leverage the key ideas from this theory (sphere constraints, datadriven initialization) to develop a practical algorithm, which performs well on data arising from a range of application areas. We highlight key additional challenges posed by the ill-conditioning of real SaSD problems, and suggest heuristics (acceleration, continuation, reweighting) to mitigate them. Experiments demonstrate the performance and generality of the proposed method.
˚YL and QQ contributed equally to this work. The full version of this work can be found at https: //arxiv.org/abs/1908.10959.
1 For simplicity, (1) uses cyclic convolution; algorithms are results also apply to linear convolution with minor modifications. Here ι denotes the zero padding operator.

1. INTRODUCTION:

Many signals arising in science and engineering can be modeled as superpositions of basic, recurring motifs, which encode critical information about a physical process of interest. Signals of this type can be modeled as the convolution of a zero-padded short kernel a 0 P R p0 (the motif) with a longer sparse signal x 0 P R m (m " p 0 ) which encodes the locations of the motifs in the sample 1 : y " ιa 0 f x 0 .

(1)

We term this a short-and-sparse (SaS) model. Since often only y is observed, short-and-sparse deconvolution (SaSD) is the problem of recovering both a 0 and x 0 from y. Variants of SaSD arise in areas such as microscopy (Cheung et al., 2018) , astronomy (Briers et al., 2013) , and neuroscience (Song et al., 2018) . SaSD is a challenging inverse problem in both theory and practice. Natural formulations are nonconvex, and very little algorithmic theory was available. Moreover, practical instances are often ill-conditioned, due to the spectral decay of the kernel a 0 (Cheung et al., 2018) .

This paper is motivated by recent theoretical advances in nonconvex optimization and, in particular, on the geometry of SaSD. Zhang et al. (2017) and Kuo et al. (2019) study particular optimization formulations for SaSD and show that the landscape is largely driven by the problem symmetries of SaSD. They derive provable methods for idealized problem instances, which exactly recover pa 0 , x 0 q up to trivial ambiguities. While inspiring, these methods are not practical and perform poorly on real problem instances. Where the emphasis of Zhang et al. (2017) and Kuo et al. (2019) is on theoretical guarantees, here we focus on practical computation. We show how to combine ideas from this theory with heuristics that better address the properties of practical deconvolution problems, to build a novel method that performs well on data arising in a range of application areas. A critical issue in moving from theory to practice is the poor conditioning of naturally-occurring deconvolution problems: we show how to address this with a combination of ideas from sparse optimization, such as momentum, continuation, and reweighting. The end result is a general purpose method, which we demonstrate on data from neural spike sorting, calcium imaging and fluorescence microscopy.

Notation. The zero-padding operator is denoted by ι : R p Ñ R m . Projection of a vector v P R p onto the sphere is denoted by P S p´1 pvq . " v{ }v} 2 , and P z pvq . " v ´xv, zy z denotes projection onto the tangent space of z P S p´1 . The Riemannian gradient of a function f : S p´1 Ñ R at point z on the sphere is given by grad f pzq . " P z p∇f pzqq.

Reproducible research. The code for implementations of our algorithms can be found online:

https://github.com/qingqu06/sparse_deconvolution.

For more details of our work on SaSD, we refer interested readers to our project website https://deconvlab.github.io/.

2. SYMMETRY AND GEOMETRY IN SASD:

In this section, we begin by describing two intrinsic properties for SaSD. Later, we show how these play an important role in the geometry of optimization and the design of efficient methods.

An important observation of the SaSD problem is that it admits multiple equivalent solutions. This is purely due to the cyclic convolution between a 0 and x 0 , which exhibits the trivial ambiguity2 

y " ιa 0 f x 0 " pαs rιa 0 sq f `1 α s ´ rx 0 s ˘,

for any nonzero scalar α and cyclic shift s r¨s. These scale and shift symmetries create several acceptable candidates for a 0 and x 0 , and in the absence of further information we only expect to recover a 0 and x 0 up to symmetry. Furthermore, they largely drive the behavior of certain nonconvex optimization problems formulated for SaSD. Since the success of SaSD requires distinguishing between overlapping copies of a 0 , its difficulty also depends highly on the "similarity" of the a 0 to its shifts. Here we capture this notion using the shift-coherence of a 0 ,

EQUATION

Intuitively, the shifts of a 0 become closer together as µpa 0 q increases (Figure 10 ), making objective landscapes for optimization less favorable for recovering any specific shift of a 0 .

2.1. LANDSCAPE GEOMETRY UNDER SHIFT-INCOHERENCE:

A natural approach to solving SaSD is to formulate it as a suitable optimization problem. In this paper we will focus on the Bilinear Lasso (BL) problem, which minimizes the squared error between the observation y and its reconstruction a f x, plus a 1 -norm sparsity penalty on x,

EQUATION

Later in this section, we will see that the kernel length p should be set slightly larger than p 0 .

The Bilinear Lasso is a nonconvex optimization problem, as the shift symmetries of SaSD create discrete local minimizers in the objective landscape. The regularization created by problem symmetries in nonconvex inverse problems are a fairly general phenomenon (Sun et al., 2015) and, as Kuo et al. (2019) shows, its influence in SaSD extends beyond the neighborhoods of these local minimizers. Kuo et al. analyzed an Approximate Bilinear Lasso (ABL) objective3 Ψ ABL , which satisfies

Ψ ABL pa, xq » Ψ BL pa, xq, when µpaq » 0.

This non-practical objective serves as a valid simplification of the Bilinear Lasso for analysis when the true kernel is itself incoherent, i.e. µpa 0 q » 0 (Figures 1d and 1e ). Under its marginalization4 

EQUATION

certain crucial properties regarding its curvature can be characterized for generic choices of x. The reason we choose to partial minimize x instead of a is because (i) the problem (4) is convex w.r.t. x, and (ii) the dimension of the subspace of a is significantly smaller than that of x (i.e., p ! m), which is the place that the measure concentrates.

Curvature in the span of a few shifts. Suppose we set p ą p 0 , which ensures that we can find an a » α 1 s 1 ra 0 s `α2 s 2 ra 0 s P S p´1 that lies near the span of two shifts of a 0 . If α 1 » ˘1 (or α 2 » 0) then, under suitable conditions on a 0 and x 0 , Kuo et al. (2019) asserts that a lies in a strongly convex region of ϕ ABL , containing a single minimizer near s 1 ra 0 s (Figure 1a ); the converse is also true. A saddle-point exists nearby when α 1 » α 2 is balanced, characterized by large negative curvature along the two shifts and positive curvature in orthogonal directions (Figure 1b ). Interpolating between these two cases, large negative gradients point towards individual shifts.

The behavior of ϕ ABL between two shifts of a 0 -strong convexity near single shifts, and saddlepoints near balanced points -extends to regions of the sphere spanned by several shifts (Figure 1c ); we elaborate on this further in Appendix A.1. This regional landscape guarantees that a 0 can be efficiently recovered up to a signed shift using methods for first and second-order descent, as soon as a can be brought sufficiently close to the span of a few shifts.

Optimization over the sphere. For both the Bilinear Lasso and ABL, a unit-norm constraint on a is enforced to break the scaling symmetry between a 0 and x 0 . Choosing the 2 -norm, however, has surprisingly strong implications for optimization. The ABL objective, for example, is piecewise concave whenever a is sufficiently far away from any shift of a 0 , but the sphere induces positive curvature near individual shifts to create strong convexity. These two properties combine to ensure recoverability of a 0 . In contrast, enforcing 1 -norm constraints often leads to spurious minimizers for deconvolution problems (Levin et al., 2011; Benichoux et al., 2013; Zhang et al., 2017) .

Initializing near a few shifts. The landscape of ϕ ABL makes single shifts of a 0 easy to locate if a is initialized near a span of a few shifts. Fortunately, this is a relatively simple matter in SaSD, as y is itself a sparse superposition of shifts. Therefore, one initialization strategy is to randomly choose a length-p 0 window r y i . " ry i y i`1 . . . y i`p0´1 s T from the observation and set a p0q .

" P S p´1 `r 0 p0´1 ; r y i ; 0 p0´1 s ˘.

(5)

This brings a p0q suitably close to the sum of a few shifts of a 0 (Appendix A.2); any truncation effects are absorbed by padding the ends of r y i , which also sets the length for a to be p " 3p 0 ´2.

Implications for practical computation. The (regionally) benign optimization landscape of ϕ ABL guarantees that efficient recovery is possible for SaSD when a 0 is incoherent. Applications of sparse deconvolution, however, are often motivated by sharpening or resolution tasks (Huang et al., 2009; Candès & Fernandez-Granda, 2014; Campisi & Egiazarian, 2016) where the motif a 0 is smooth and coherent (i.e. µpa 0 q is large). The ABL objective is a poor approximation of the Bilinear Lasso in such cases and fails to yield practical algorithms, so we should optimize the Bilinear Lasso directly.

From Figures 1d and 1e , we can see that low-dimensional subspheres spanned by shifts of a 0 are empirically similar when a 0 is incoherent. Although this breaks down in the coherent case, as we illustrate in Appendix A.3, the symmetry breaking properties of ϕ BL remain present. This allows us to apply the geometric intuition discussed here to create an optimization method that, with the help of a number of computational heuristics, performs well in for SaSD even in general problem instances.

Algorithm 1 Inertial Alternating Descent Method (iADM)

Input: Initializations a p0q P S p´1 , x P R m ; observation y P R m ; penalty λ ě 0; momentum α P r0, 1q.

Output: pa pkq , x pkq q, a local minimizer of ΨBL.

Initialize a p1q " a p0q , x p1q " x p0q . for k " 1, 2, . . . until converged do Update x with accelerated proximal gradient step:

w pkq Ð x pkq `α ¨`x pkq ´xpk´1q xpk`1q Ð soft λt k " w pkq ´tk ¨∇xψλ `apkq , w pkq ˘‰,

where soft λ pvq .

" signpvq d maxp|v ´λ| , 0q denotes the soft-thresholding operator. Update a with accelerated Riemannian gradient step:

z pkq Ð P S p´1

`apkq `α xa pkq ,a pk´1q y ¨Pa pk´1q `apkq ˘ȃ 

3. DESIGNING A PRACTICAL SASD ALGORITHM:

Several algorithms for SaSD-type problems have been developed for specific applications, such as image deblurring (Levin et al., 2011; Briers et al., 2013; Campisi & Egiazarian, 2016) , neuroscience (Rey et al., 2015; Friedrich et al., 2017; Song et al., 2018) , and image super-resolution (Baker & Kanade, 2002; Shtengel et al., 2009; Yang et al., 2010) , or are augmented with additional structure (Wipf & Zhang, 2014; Ling & Strohmer, 2017; Walk et al., 2017) .

Here, we instead leverage the theory from Section 2 to build an algorithm for general practical settings. In addition to applying an appropriate initialization scheme (5) and optimizing on the sphere, we minimize the Bilinear Lasso (3) instead of the ABL (4) to more accurately account for interactions between shifts of a 0 in highly shift-coherent settings. Furthermore, we also address the negative effects of large coherence using a number of heuristics, leading to an efficient algorithm for SaSD.

Momentum acceleration. In shift-coherent settings, the Hessian of Ψ BL becomes ill-conditioned5 near shifts of a 0 , a situation known to cause slow convergence for first-order methods (Nesterov, 2013) . A remedy is to add momentum (Polyak, 1964; Beck & Teboulle, 2009) to first-order iterations, for instance, by augmenting gradient descent on some smooth f pzq with stepsize τ with the term w, w pkq Ð z pkq `α ¨pz pkq ´zpk´1q q (6) z pk`1q Ð w pkq ´τ ¨∇f pw pkq q.

(7)

Here, α controls the momentum added6 . As illustrated in Figure 2 , this additional term improves convergence by reducing oscillations of the iterates for ill-conditioned problems. Momentum has been shown to improve convergence for nonconvex and nonsmooth problems (Pock & Sabach, 2016; Jin et al., 2018) . Here we provide an inertial alternating descent method (iADM) for finding local minimizers of Ψ BL (Algorithm 1), which modifies iPALM (Pock & Sabach, 2016) to perform updates on a via retraction on the sphere (Absil et al., 2009) 7 .

Algorithm 2 SaS-BD with homotopy continuation

Input: Observation y P R m , motif size p0; momentum α P r0, 1q; initial λ p1q final λ ‹ , penalty decrease η P p0, 1q; precision factor δ P p0, 1q. Output: Solution path pâ pnq , xpnq ; λ pnq q ( for SaSD. Set number of iterations N Ð X logpλ ‹ {λ p1q q { log η \ . Initialize âp0q P R 3p 0 ´2 using ( 5 When solving (3) in the noise-free case, it is clear that larger choices of λ encourage sparser solutions for x. Conversely, smaller choices of λ place local minimizers of the marginal objective ϕ BL paq . " min x Ψ BL pa, xq closer to signed-shifts of a 0 by emphasizing reconstruction quality. When µpa 0 q is large, however, ϕ BL becomes ill-conditioned as λ Ñ 0 due to the poor spectral conditioning of a 0 , leading to severe flatness near local minimizers and the creation spurious local minimizers when noise is present (Figure 3 ). Conversely, larger values of λ limit x to a small set of support patterns and simplify the landscape ofϕ BL , at the expense of precision.

It is therefore important both for fast convergence and accurate recovery for λ to be chosen appropriately. When problem parameters -such as noise level, p 0 , or θ -are not known a priori, a homotopy continuation method (Hale et al., 2008; Wright et al., 2009; Xiao & Zhang, 2013 ) can be used to obtain a range of solutions for SaSD. Using initialization (5), a rough estimate pâ p1q , xp1q q is obtained by solving (3) with iADM using a large choice for λ p1q . This estimate is refined via a solution path pâ pnq , xpnq ; λ pnq q ( by gradually decreasing λ pnq . By ensuring that x remains sparse along the solution path, the objective Ψ BL enjoys restricted strong convexity w.r.t. both a and x throughout optimization (Agarwal et al., 2010) . As a result, homotopy achieves linear convergence for SaSD where sublinear convergence is expected otherwise (Figures 4c and 4d ). We provide a complete algorithm for SaSD combining Bilinear Lasso and homotopy continuation in Algorithm 2.

4. EXPERIMENTS:



4.1. SYNTHETIC EXPERIMENTS:

Here we perform SaSD in simulations on both coherent and incoherent settings. Coherent kernels are discretized from the Gaussian window function a 0 " g p0,0.5 , where g p,σ

.

" P S p´1 `" exp `ṕ2i´p´1q 2 σ 2 pp´1q 2 ˘‰p i"1 ˘.

Incoherent kernels a 0 " UnifpS p0´1 q are sampled uniformly on the sphere. Recovery performance. We test recovery probability for varying kernel lengths p 0 and sparsity rates θ. To ensure the problem size is sufficiently large, we set m " 100p 0 . For each p 0 and θ, we randomly generate 8 x " i.i.d. BRpθq for both coherent and incoherent a 0 . We solve ten trials of (3) on clean observation data a 0 f x 0 using iADM with λ "10 ´2 ? p0θ . The probability of recovering a signed shift of a 0 is shown in Figure 4 . Recovery is likely when sparsity is low compared to the kernel length. The coherent problem setting has a smaller success region compared to the incoherent setting.

Momentum and homotopy. Next, we test the performance of Algorithm 1 with momentum (α k " k´1 k`2 ; see Pock & Sabach (2016) ) and without (α " 0). This is done by minimizing Ψ BL with initialization (5), using clean observations with p 0 " 10 2 , m " 10 4 , and θ " p ´3{4 0 for coherent and incoherent a 0 . We also apply homotopy (Algorithm 2) with λ p1q " max |xs ra p0q s, yy| -see Xiao & Zhang (2013) , λ ‹ " 0.3 ? p0λ , η " 0.8, and δ " 0.1. The final solve of (3) uses precision ε ‹ " 10 ´6, regardless of method. Figures 4c and 4d show the comparison results on coherent problem settings.

Comparison to existing methods. Finally, we compare iADM, and iADM with homotopy, against a number of existing methods for minimizing ϕ BL . The first is alternating minimization (Kuo et al., 2019) , which at each iteration k minimizes a pkq with x pkq fixed using accelerated (Riemannian) gradient descent with backtracking, and vice versa. The next method is the popular alternating direction method of multipliers (Boyd et al., 2011) . Finally, we compare against iPALM (Pock & Sabach, 2016) with backtracking, using the unit ball constraint on a 0 instead of the unit sphere.

For each method, we deconvolve signals with p 0 " 50, m " 100p 0 , and θ " p ´3{4 0 for both coherent and incoherent a 0 . For both iADM, iADM with homotopy, and iPALM we set α " 0.3. For homotopy, we set λ p1q " max |xs ra p0q s, yy|, λ ‹ " 0.3 ? p0λ , and δ " 0.5. Furthermore we set η " 0.5 or η " 0.8 and for ADMM, we set the slack parameter to ρ " 0.7 or ρ " 0.5 for incoherent and coherent a 0 respectively. From Figure 6 , we can see that ADMM performs better than iADM in the incoherent case, but becomes less reliable in the coherent case. In both cases, iADM with homotopy is the best performer. Finally, we observe roughly equal performance between iPALM and iADM. 

4.2. IMAGING APPLICATIONS:

Here we demonstrate the performance and generality of the proposed method. We begin with calcium fluorescence imaging, a popular modality for studying spiking activity in large neuronal populations (Grienberger & Konnerth, 2012) , followed by stochastic optical reconstruction microscopy (STORM) (Rust et al., 2006; Huang et al., 2008; 2010) , a superresolution technique for in vivo microscopy9 .

Sparse deconvolution of calcium signals. Neural spike trains created by action potentials, each inducing a transient response in the calcium concentration of the surrounding environment. The aggregate signal can be modeled as a convolution between the transient a 0 and the spike train x 0 . Whilst a 0 and x 0 both encode valuable information, neither are perfectly known ahead of time.

Here, we first test our method on synthetic data generated using an AR2 model for a 0 , a shiftcoherent kernel that is challenging for deconvolution, see e.g. Friedrich et al. (2017) . We set x 0 " i.i.d. Bernoullipp

´4{5 0 q P R 10 4 with additive noise n " i.i.d. N p0, 5 ¨10 ´2q. Figures 5a and 5b demonstrate accurate recovery of a 0 and x 0 in this synthetic setting. Next, we test our method on real data 10 ; Figures 5c and 5d demonstrate recovery of spike locations. Although iADM provides decent performance, in the presence of large noise estimation quality can be improved by stronger sparsification methods, such as the reweighting technique by Candes et al. ( 2008), which we elaborate on in Appendix B. Additionally, Figure 6c shows that the proposed method converges to higher precision in comparison with state-of-the-art methods. Super-resolution for fluorescence microscopy. Fluorescence microscopy is often spatially limited by the diffraction of light; its wavelength (several hundred nanometers) is often larger than typical molecular length-scales in cells, preventing a detailed characterization of subcellular structures. The STORM technique overcomes this resolution limit by using photoswitchable fluorescent probes to multiplex the image into multiple frames, each containing a subset of the molecules present (Figure 7 ). If the location of these molecules can be precisely determined for each frame, synthesizing all deconvolved frames will produce a super-resolution microscopy image with nanoscale resolution.

For each image frame, the localization task can be formulated via the SaS model

EQUATION

where f denotes 2D convolution. Here we will solve this task on the single-molecule localization microscopy (SMLM) benchmarking dataset11 via SaSD, recovering both the PSF A 0 and the point source maps X 0,t simultaneously. We apply iADM with reweighting (Appendix B) on frames of size 128 ˆ128 from the video sequence "Tubulin"; each pixel is of 100nm 2 resolution12 , the fluorescence wavelength is 690nm, and the framerate is f " 25Hz. Figure 7 shows examples of recovered activation maps, and the aggregated super-resolution image from all 500 frames, accurately predicting the PSF (see Appendix D) and the activation map for each video frame to produce higher resolution microscopy images.

Localization in calcium images. Our methods are easily extended to handle superpositions of multiple SaS signals. In calcium imaging, this can potentially be used to track the neurons in video sequences, a challenging task due to (non-) rigid motion, overlapping sources, and irregular background noise Pnevmatikakis et al. ( 2016); Giovannucci et al. (2019) . We consider frames video obtained via the two-photon calcium microscopy dataset from the Allen Institute for Brain Science13 , shown in Figure 8 . Each frame contains the cross section of several neurons and dendrites, which have distinct sizes. We model this as the SaS signal Y t " ιA 1 f X 1,t `ιA 2 f X 2,t , where each summand consists of neurons or dendrites exclusively. By extending Algorithm 2 to recover each of the kernels A k and maps X k , we can solve this convolutional dictionary learning (SaS-CDL; see Appendix C) problem which allows us to separate the dendritic and neuronal components from this image for localization of firing activity, etc. As a result, the application of SaS-CDL as a denoising or analysis tool for calcium imaging videos provides a very promising direction for future research.

5. DISCUSSION:

Many nonconvex inverse problems, such as SaSD, are strongly regulated by their problem symmetries.

Understanding this regularity and when or how it breaks down is important for developing effective algorithms. We illustrate this by combining geometric intuition with practical heuristics, motivated by common challenges in real deconvolution, to produce an efficient and general purpose method that performs well on data arising from a range of application areas. Our approach, therefore, can serve as a general baseline for studying and developing extensions to SaSD, such as SaS-CDL (Bristow & Lucey, 2014; Chun & Fessler, 2017; Garcia-Cardona & Wohlberg, 2018) 

A APPROXIMATE BILINEAR LASSO OBJECTIVE:

Recall from Section 2.2 of the main text that SaSD can be formulated as the Bilinear Lasso problem

EQUATION

Unfortunately, this objective is challenging for analysis. A major culprit is that its marginalization

EQUATION

generally does not admit closed form solutions due the convolution with a in the squared error term. This motivates Kuo et al. (2019) to study the nonconvex formulation

EQUATION

We refer to (11) as the Approximate Bilinear Lasso formulation, and it is quite easy to see that Here soft denotes the elementwise soft-thresholding operator soft t px i q " signpx i q ¨maxp|x i | ´t, 0q, and q a denotes the adjoint kernel of a, i.e. the kernel s.t. xιa f u, vy " xu, q a f vy @u, v P R m .

Ψ ABL pa, xq « Ψ BL pa, xq when }a f x} 2 « }x} 2 , i.e. if

A.1 LANDSCAPE GEOMETRY

The rest of Section 2.2 discusses the regional characterization of ϕ ABL in the span of a small number of shifts from a 0 . This language is made precise in the form of the subsphere

EQUATION

spanned by a small set of cyclic shifts of ιa 0 . Although we will not discuss the explicit distance function here, the characterization by Kuo et al. (2019) holds whenever a is close enough to such a subsphere with |I| ď 4θp 0 , where θ is the probability that any individual entry of x 0 is nonzero. Suppose we have a « ř PI α s rιa 0 s for some appropriate index set I. Note that if µ s a 0 « 0, then µ s a « 0, @a P S I . Now let α p1q and α p2q be the first and second largest coordinates of the shifts participating in a, and let s p1q ra 0 s and s p2q ra 0 s be the corresponding shifts. Then

• If ˇˇα p2q

α p1q ˇˇ« 0, then a is in a strongly convex region of ϕ ABL , containing a single local minimizer corresponding to s p1q ra 0 s.

• If ˇˇα p2q:

α p1q ˇˇ« 1, then a is near a saddle-point, with negative curvature pointing towards s p1q ra 0 s and s p2q ra 0 s. If ˇˇα p3q α p2q ˇˇ« 0 , i.e. s p1q ra 0 s and s p2q ra 0 s are the only two participating shifts, then ϕ ABL is also characterized by positive curvature in all orthogonal directions.

• Otherwise, x´gradϕ ABL paq, z ´ay takes on a large positive value, for either u " s p1q ra 0 s or u " s p2q ra 0 s, i.e. the negative Riemannian gradient is large and points towards one of the participating shifts.

Data y Kernel a 0 Sparse x 0

Truncation a p´1q Initialization a p0q " « α i s i ra 0 s `αj s j ra 0 s

Figure 9 : Data-driven initialization for a: using a piece of the observed data y to generate a good initial point a p0q . Top: data y " a0 f x0 is a superposition of shifts of the true kernel a0. Bottom: a length-p0 window contains pieces of just a few shifts. Bottom-center: one step of the generalized power method approximately fills in the missing pieces, yielding an initialization that is close to a linear combination of shifts of a0 (right).

This is an example of a ridable saddle property (Jin et al., 2017) that allows many first and secondorder methods to locate local minimizers. Since all local minimizers of ϕ ABL near S I must correspond to signed-shifts of a 0 , this guarantees that the Approximate Bilinear Lasso formulation can be efficiently solved to recover a 0 (and subsequently x 0 ) for incoherent a 0 , as long as a is initialized near some appropriate subsphere and the sparsity coherence tradeoff p 0 θ AE pµ s pa 0 qq ´1{2 is satisfied. We note that this is a poor tradeoff rate, which reflects that the Approximate Bilinear Lasso formulation is non-practical and cannot handle SaSD problems involving kernels with high shift-coherence.

A.2 DATA-DRIVEN INITIALIZATION:

For the SaS-BD problem, we usually initialize x by x p0q " 0, so that our initialization is sparse. For the optimization variable a P R n , recall from Section 2.2 in the main text that it is desirable to obtain an initialization a 0 which is close to the intersection of S p´1 and a subsphere S I spanned by a few shifts of a 0 . When x 0 is sparse, our measurement y is a linear combination of a few shifts of a 0 . Therefore, an arbitrary consecutive p 0 -length window r y i . " ry i y i`1 . . . y i`p0´1 s T of the data y should be not far away from such a subspace S I . As illustrated in Figure 9 , one step of the generalized power method (Kuo et al., 2019) r a p0q .

" P S p´1 `r 0 p´1 ; r y i ; 0 p´1 s ˘(14) a p0q " P S p´1 ´´∇ϕ ABL ´r a p0q ¯¯(15) produces a refined initialization that is very close to a subspace S I spanned by a few shifts of a 0 with |I| « θp 0 . However, ( 15) is a relatively complicated for a simple idea. In practice, we find that the simple initialization a p0q " r a p0q from ( 14) works suitably well for solving SaSD with (9).

A.3 COMPARISON TO THE BILINEAR LASSO:

Although it is easy to see that Ψ ABL paq and Ψ BL paq are similar as long as µpaq « 0, it is also clear that these two quantities can be very different when µpaq is large. This is especially significant when µpa 0 q is itself large, as the desired solutions for a are then also coherent.

From Figure 10 , we can see that these changes are reflected in the low-dimensional subspheres (13) spanned by adjacent shifts of a 0 . Compared to the incoherent case, ϕ BL also takes on small values in regions between adjacent shifts, creating a "global valley" on the subsphere. Theoretically, this makes it difficult to ensure exact recovery of up to symmetry when a 0 is coherent, and the objective function becomes much more complicated. This is not a significant issue in terms of practical computation, however, since adjacent shifts of a 0 become indistinguishable as µpa 0 q Ñ 1, meaning that one only needs to ensure that a lands in the "global valley" to achieve good estimates of a 0 up to symmetry. 

B REWEIGHTED SPARSE PENALIZATION:

When a 0 is shift-coherent, minimization of the objective Ψ BL with respect to x becomes sensitive to perturbations, creating "smudging" effects on the recovered map x. These resolution issues can be remedied with stronger concave regularizers. A simple way of facilitating this with the Bilinear Lasso is to use a reweighting technique (Candes et al., 2008) . The basic idea is to adaptively adjust the penalty by considering a weighted variant of the original Bilinear Lasso problem from ( 9 

where w P R m `and d denotes the Hadamard product. Here we will set the weights w to be roughly inverse to the magnitude of the true signal x 0 , i.e.,

EQUATION

Algorithm 3 Reweighted Bilinear Lasso Input: Initializations âp0q , xp0q , penalty λ ą 0 Output: Local minimizers âpjq , xpjq of Ψ w pjq BL . Initialize w p1q " 1 m , j Ð 1. while not converged do Using the initialization `â pj´1q , xpj´1q ˘and weight w pjq , solve ( 16) -e.g. with iADM -to obtain solution `â pjq , xpjq ˘;

Set ε with (19) and update the weights as

EQUATION

)

Update Ð `1. end while

In addition to choosing λ ą 0, here ε ą 0 trades off between sparsification strength (small ε) and algorithmic stability (large ε). Let |x| piq denote the i-th largest entry of |x|. For experiments in the main text, we set ε " max

EQUATION

)

Starting with the initial weights w p0q " 1 m , Algorithm 3 successively solves ( 16), updating the weights using (17) at each outer loop iteration j. As j Ñ 8, this method becomes equivalent to replacing the 1 -norm in (9) with the nonconvex penalty ř i logp|x i | `εq (Candes et al., 2008) . We can easily adopt our iADM algorithm to solve this subproblem, by taking the proximal gradient on x with a different penalty λ i for each entry x i . Figure 11 , as well as calcium imaging experiments in Section 4.2, Figure 5 of the main text, demonstrate improved estimation as a result of this method. The optimization methods we introduced for SaSD here can be naturally extended for sparse blind deconvolution problems with multiple kernels/motifs (a.k.a. convolutional dictionary learning; see Garcia-Cardona & Wohlberg (2018) ), which have broad applications in microscopy data analysis (Yellin et al., 2017; Zhou et al., 2014; Cheung et al., 2018) and neural spike sorting (Ekanadham et al., 2011; Rey et al., 2015; Song et al., 2018) . As illustrated in Figure 12 , the new observation y is

:

We therefore assume w.l.o.g. that }a0} 2 " 1 in this paper.

As the intention here is apply some key intuition from the ABL objective towards the Bilinear Lasso itself, we intentionally omit the concrete form of ΨABLpaq. Readers may refer to Appendix A for more details.

Minimizing ϕABL, this is equivalent to minimizing ΨABL as x can be recovered via convex optimization.

This is because the circulant matrix Ca 0 is ill-conditioned.

Setting α " 0 removes momentum and reverts to standard gradient descent.

The stepsizes t k and τ k are obtained by backtracking(Nocedal & Wright, 2006;Pock & Sabach, 2016) to ensure sufficient decrease for ΨBL `apkq , w pkq ˘´ΨBL `apkq , x pk`1q ˘, and vice versa.

BRpθq denotes the Bernoulli-Rademacher distribution, which has values ˘1 w.p. θ{2 and zero w.p. 1 ´θ.

Other superresolution methods for microscopy include photoactivated localization microscopy (PALM)(Betzig et al., 2006), and fluorescence photoactivation localization microscopy (fPALM)(Hess et al., 2006).

Obtained at http://spikefinder.codeneuro.org.

Data can be accessed at http://bigwww.epfl.ch/smlm/datasets/index.html.

Here we solve SaSD on the same 128 ˆ128 grid. In practice, the localization problem is solved on a finer grid, so that the resulting resolution can reach 20 ´30 nm.

Obtained at http://observatory.brain-map.org/visualcoding/.

Figures and Tables:

Figure 1: Geometry of ϕ ABL near superpositions of shifts of a0 (Kuo et al., 2019). (a) Regions near single shifts are strongly convex. (b) Regions between two shifts contain a saddle-point, with negative curvature towards each shift and positive curvature orthogonally. (c) The span of three shifts. For each figure, the top shows the function value in height, and the bottom shows function value over the sphere. (d,e) When µspa0q « 0, the Bilinear Lasso ϕBLpaq . " minx ΨBLpa, xq and ABL ϕABLpaq are empirically similar in the span of three shifts.

Figure 2: Momentum acceleration. a) Iterates of gradient descent oscillate on ill-conditioned functions; each marker denotes one iteration. b) Momentum dampens oscillation and speeds up convergence.

Figure 3: Bilinear-lasso objective ϕ λ on the sphere S p´1 , for p " 3 and varying λ; brighter colors indicate higher values. The function landscape of ϕ λ flattens as sparse penalty λ decreases from left to right.

Figure 4: Synthetic experiments for Bilinear Lasso. Success probability (a, b): x0 "i.i.d. BRpθq, the success probability of SaS-BD by solving (3), shown by increasing brightness, is large when the sparsity rate θ is sufficiently small compared to the length of a0, and vice versa. Success with a fixed sparsity rate is more likely when a0 is incoherent. Algorithmic convergence (c, d): iterate convergence for iADM with α k " pk ´1q{pk `1q vs. α k " 0 (ADM); with and without homotopy. Homotopy significantly improves convergence rate, and momentum improves convergence when a0 is coherent.

Figure 6: Algorithmic comparison. (a) Convergence of various methods minimizing ΨBL with incoherent a0 over FFT operations used (for computing convolutions). The y-axis denotes the log of the angle between a pkq and the nearest shift of a0, and each marker denotes five iterations. (b) Convergence for coherent a0, and (c) with an AR2 kernel for modeling calcium signals.

Figure 7: SaSD for STORM imaging. (a, b) Individual frames (left) and predicted point process map using SaSD (right). (c, d) shows the original microscopy and the super-resolved image obtained by our method.

Figure 10: Low-dimensional subspheres spanned by shifts of a0. Subfigures (a,b) present the optimization landscapes of ϕABLpaq and ϕBLpaq, for a P S p´1 Ş spanta0, s1ra0s, s2ra0u, with higher values being brighter. The red dots denote the shifts of a0. Subfigure (c) shows the landscape ϕBL when a0 is coherent, which significantly departs from the landscapes of (a,b), but still retains symmetry breaking curvature.

True map x0 (b) True motif a0 (c) Noisy y, 1 only (d) Noisy a, 1 only (e) Noisy y, reweighted (f) Noisy a, reweighted

Figure 11: Recovery of x0 with 1-reweighting. (a, b) Truth signals. (c) Solving minx ΨBLpa, xq with noisy data and coherent a0 leads to low-quality estimates of x; (d) performance suffers further when a is a noisy estimate of a0. (e, f) Reweighted 1 minimization alleviates this issue significantly.

, Bayesian approaches(Babacan et al., 2008;Wipf & Zhang, 2014), and hierarchical SaS models(Chen et al., 2013). Yuqian Zhang, Yenson Lau, Han-Wen Kuo, Sky Cheung, Abhay Pasupathy, and John Wright. On the global geometry of sphere-constrained sparse blind deconvolution. In Computer Vision and Pattern Recognition (CVPR), 2017 IEEE Conference on, pp. 4381-4389. IEEE, 2017. Yin Zhou, Hang Chang, Kenneth Barner, Paul Spellman, and Bahram Parvin. Classification of histology sections via multispectral convolutional sparse coding. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 3081-3088, 2014.

Table content: <table/>

a is shift-incoherent, or µpaq « 0. The marginalized objective function ϕ BL paq . " min x Ψ DQ pa, xq now has the closed form expression

Table content: <table><tr><td>ϕ ABL paq</td><td>. " ´1 2 }soft λ rq a f ys} 2 2 .</td><td>(12)</td></tr></table>

Peer review (iclr2020):
* Summary The work considers sparse and short blind deconvolution problem, which is to inverse a convolution of a sparse source (such as spikes at cell locations in microscopy) with a short (of limited spatial size) kernel or point spread function, not known in advance. This is posed as a bilinear lasso optimization problem. The work applies a non-linear optimization method with some practical improvements (such as data-driven initialization, momentum, homotopy continuation) and applies it to 3 real imagine problems. I think the work can be described as bridging from the known problem formulation, available theoretical understanding of its properties and available selection of optimization methods, to an implementation that can be applied in practical cases. The practical improvements made are not specifically novel, but result in a well-fit optimization method. On the practical end, it is shown that the method can be applied in multiple cases but it is not demonstrated to give practical improvements over any alternative reconstruction approaches. The emphasis is more on showcasing possible further extensions (many deferred to appendices). I view this work as not very strong but a valid contribution. * Detailed Comments While the ideas from a provable algorithm by Kuo et al. are used here (sphere constraints, data-driven initialization), can the authors show in their setting, this leads as well to optimal recovery guarantees for sufficiently simple problems? Somehow the discussion advocates spherical constraint, because different shifts of the kernel become different local minima. But at the same time, this increases the non-linearity of the problem, and thus makes it more difficult to solve. Although, these multiple local optima provide equivalent solutions, I find this somewhat counter-intuitive. Cannot the shift ambiguity be resolved in a convex manner, e.g. by fixing the mean value of the kernel? In Fig 1 an explanation of the axis and projections would be needed. While in the introduction, the setting m>>p0 is assumed, in the experiments n0 is used to denote kernel width and some experiments actually work in the setting of comparable values such as n0=50, m=100. What are similarities / differences to related reconstruction problems such as non-negative matrix factorization with sparseness constraints?