Year: 2020

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.

heading: 1 INTRODUCTION
text: 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 a0 P Rp0 (the motif) with a longer sparse signal x0 P Rm (m " p0) which encodes the locations of the motifs in the sample1:
y “ ιa0 f x0. (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 a0 and x0 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 a0 (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
˚YL and QQ contributed equally to this work. The full version of this work can be found at https: //arxiv.org/abs/1908.10959.
1For simplicity, (1) uses cyclic convolution; algorithms are results also apply to linear convolution with minor modifications. Here ι denotes the zero padding operator.
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 pa0,x0q 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 ι : Rp Ñ Rm. Projection of a vector v P Rp onto the sphere is denoted by PSp´1pvq .“ v{ }v}2, and Pzpvq .“ v ´ xv, zy z denotes projection onto the tangent space of z P Sp´1. The Riemannian gradient of a function f : Sp´1 Ñ R at point z on the sphere is given by grad fpzq .“ Pzp∇fpzqq.
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/.

heading: 2 SYMMETRY AND GEOMETRY IN SASD
text: 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 a0 and x0, which exhibits the trivial ambiguity2
y “ ιa0 f x0 “ pαs` rιa0sq f ` 1 αs´` rx0s ˘ ,
for any nonzero scalar α and cyclic shift s` r¨s. These scale and shift symmetries create several acceptable candidates for a0 and x0, and in the absence of further information we only expect to recover a0 and x0 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 a0, its difficulty also depends highly on the “similarity” of the a0 to its shifts. Here we capture this notion using the shift-coherence of a0,
µpa0q .“ max
` “0 |xιa0, s` rιa0sy| P r0, 1s . (2)
Intuitively, the shifts of a0 become closer together as µpa0q increases (Figure 10), making objective landscapes for optimization less favorable for recovering any specific shift of a0.

heading: 2.1 LANDSCAPE GEOMETRY UNDER SHIFT-INCOHERENCE
text: 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 af x, plus a `1-norm sparsity penalty on x,
min aPSp´1,xPRm
”
ΨBLpa,xq .“ 12 }y ´ ιaf x} 2 2 ` λ }x}1
ı
. (3)
Later in this section, we will see that the kernel length p should be set slightly larger than p0.
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
2We therefore assume w.l.o.g. that }a0}2 “ 1 in this paper.
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
ΨABLpa,xq » ΨBLpa,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. µpa0q » 0 (Figures 1d and 1e). Under its marginalization4
ϕABLpaq .“ minxPRm ΨABLpa,xq, (4)
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 ą p0, which ensures that we can find an a » α1s`1ra0s ` α2s`2ra0s P Sp´1 that lies near the span of two shifts of a0. If α1 » ˘1 (or α2 » 0) then, under suitable conditions on a0 and x0, Kuo et al. (2019) asserts that a lies in a strongly convex region of ϕABL, containing a single minimizer near s`1ra0s (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 a0 — 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 a0 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 a0 and x0. 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 a0, but the sphere induces positive curvature near individual shifts to create strong convexity. These two properties combine to ensure recoverability of a0. 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 a0 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
3As 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.
4Minimizing ϕABL, this is equivalent to minimizing ΨABL as x can be recovered via convex optimization.
itself a sparse superposition of shifts. Therefore, one initialization strategy is to randomly choose a length-p0 window ryi .“ ryi yi`1 . . . yi`p0´1s T from the observation and set
ap0q .“ PSp´1 ` r 0p0´1 ; ryi ; 0p0´1 s ˘ . (5)
This brings ap0q suitably close to the sum of a few shifts of a0 (Appendix A.2); any truncation effects are absorbed by padding the ends of ryi, which also sets the length for a to be p “ 3p0 ´ 2. Implications for practical computation. The (regionally) benign optimization landscape of ϕABL guarantees that efficient recovery is possible for SaSD when a0 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 a0 is smooth and coherent (i.e. µpa0q 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 a0 are empirically similar when a0 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 ap0q P Sp´1, x P Rm; observation y P Rm; penalty λ ě 0; momentum α P r0, 1q. Output: papkq,xpkqq, a local minimizer of ΨBL.
Initialize ap1q “ ap0q, xp1q “ xp0q. for k “ 1, 2, . . . until converged do
Update x with accelerated proximal gradient step:
wpkq Ð xpkq ` α ¨ ` xpkq ´ xpk´1q ˘
xpk`1q Ð softλtk “ wpkq ´ tk ¨∇xψλ ` apkq,wpkq ˘‰ ,
where softλpvq .“ signpvq dmaxp|v ´ λ| ,0q denotes the soft-thresholding operator. Update a with accelerated Riemannian gradient step:
zpkq Ð PSp´1 ` apkq ` αxapkq,apk´1qy ¨ Papk´1q ` apkq ˘˘
apk`1q Ð PSp´1 ` zpkq ´ τk ¨ grada ψλ ` zpkq,xpk`1q ˘˘
. end for

heading: 3 DESIGNING A PRACTICAL SASD ALGORITHM
text: 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 a0 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 a0, 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 fpzq with stepsize τ with the term w,
wpkq Ð zpkq ` α ¨ pzpkq ´ zpk´1qq (6) zpk`1q Ð wpkq ´ τ ¨∇fpwpkqq. (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 Rm, motif size p0; momentum α P r0, 1q; initial λp1q final λ‹, penalty decrease η P p0, 1q; precision factor δ P p0, 1q.
Homotopy continuation. It is also possible to improve optimization by modifying the objective ΨBL directly through the sparsity penalty λ. Variations of this idea appear in both Zhang et al. (2017) and Kuo et al. (2019), and can also help to mitigate the effects of large shift-coherence.
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 ϕBLpaq
.“ minx ΨBLpa,xq closer to signed-shifts of a0 by emphasizing reconstruction quality. When µpa0q is large, however, ϕBL becomes ill-conditioned as λ Ñ 0 due to the poor spectral conditioning of a0, 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, p0, 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 pâp1q, x̂p1qq
5This is because the circulant matrix Ca0 is ill-conditioned. 6Setting α “ 0 removes momentum and reverts to standard gradient descent. 7The stepsizes tk and τk are obtained by backtracking (Nocedal & Wright, 2006; Pock & Sabach, 2016) to
ensure sufficient decrease for ΨBL ` apkq,wpkq ˘ ´ΨBL ` apkq,xpk`1q ˘ , and vice versa.
is obtained by solving (3) with iADM using a large choice for λp1q. This estimate is refined via a solution path pâpnq, x̂pnq;λpnqq (
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.

heading: 4 EXPERIMENTS

heading: 4.1 SYNTHETIC EXPERIMENTS
text: Here we perform SaSD in simulations on both coherent and incoherent settings. Coherent kernels are discretized from the Gaussian window function a0 “ gp0,0.5, where gp,σ .“ PSp´1 `“ exp `
´ p2i´p´1q2 σ2pp´1q2 ˘‰p i“1 ˘ . Incoherent kernels a0 „ UnifpSp0´1q are sampled uniformly on the sphere.
Recovery performance. We test recovery probability for varying kernel lengths p0 and sparsity rates θ. To ensure the problem size is sufficiently large, we set m “ 100p0. For each p0 and θ, we randomly generate8 x „i.i.d. BRpθq for both coherent and incoherent a0. We solve ten trials of (3) on clean observation data a0 f x0 using iADM with λ “ 10
´2 ? p0θ . The probability of recovering a signed
8BRpθq denotes the Bernoulli-Rademacher distribution, which has values ˘1 w.p. θ{2 and zero w.p. 1´ θ.
shift of a0 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´1k`2 ; see Pock & Sabach (2016)) and without (α“0). This is done by minimizing ΨBL with initialization (5), using clean observations with p0 “ 102, m “ 104, and θ “ p´3{40 for coherent and incoherent a0. We also apply homotopy (Algorithm 2) with λp1q “ max`|xs`rap0qs,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 apkq with xpkq 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 a0 instead of the unit sphere.
For each method, we deconvolve signals with p0 “ 50,m “ 100p0, and θ “ p´3{40 for both coherent and incoherent a0. For both iADM, iADM with homotopy, and iPALM we set α “ 0.3. For homotopy, we set λp1q “ max`|xs`rap0qs,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 a0 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.

heading: 4.2 IMAGING APPLICATIONS
text: 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 a0 and the spike train x0. Whilst a0 and x0 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 a0, a shiftcoherent kernel that is challenging for deconvolution, see e.g. Friedrich et al. (2017). We set x0 „i.i.d. Bernoullipp´4{50 q P R10 4
with additive noise n „i.i.d. N p0, 5 ¨ 10´2q. Figures 5a and 5b demonstrate accurate recovery of a0 and x0 in this synthetic setting. Next, we test our method on real data10; Figures 5c and 5d demonstrate recovery of spike locations. Although iADM provides
9Other superresolution methods for microscopy include photoactivated localization microscopy (PALM) (Betzig et al., 2006), and fluorescence photoactivation localization microscopy (fPALM) (Hess et al., 2006).
10Obtained at http://spikefinder.codeneuro.org.
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
Yt lomon
STORM frame
“ ιA0 lomon
point spread function
f X0,t lomon
sparse point sources
` Nt lomon
noise
, (8)
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 PSFA0 and the point source mapsX0,t simultaneously. We apply iADM with reweighting (Appendix B) on frames of size 128ˆ 128 from the video sequence “Tubulin”; each pixel is of 100nm2 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
11Data can be accessed at http://bigwww.epfl.ch/smlm/datasets/index.html. 12Here 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.
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 Yt “ ιA1 fX1,t ` ιA2 fX2,t, where each summand consists of neurons or dendrites exclusively. By extending Algorithm 2 to recover each of the kernels Ak and maps Xk, 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.

heading: 5 DISCUSSION
text: 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), Bayesian approaches (Babacan et al., 2008; Wipf & Zhang, 2014), and hierarchical SaS models (Chen et al., 2013).

heading: ACKNOWLEDGMENTS
text: This work was funded by NSF 1343282, NSF CCF 1527809, and NSF IIS 1546411. QQ also acknowledges supports from Microsoft PhD fellowship and the Moore-Sloan fellowship. We would like to thank Gongguo Tang, Shuyang Ling, Carlos Fernandez-Granda, Ruoxi Sun, and Liam Paninski for fruitful discussions.

heading: A APPROXIMATE BILINEAR LASSO OBJECTIVE
text: Recall from Section 2.2 of the main text that SaSD can be formulated as the Bilinear Lasso problem
min aPSp´1,xPRm
”
ΨBLpa,xq .“ 12 }y ´ ιaf x} 2 2 ` λ }x}1
ı
. (9)
Unfortunately, this objective is challenging for analysis. A major culprit is that its marginalization
ϕBLpaq .“ min
x
!
1 2 }y ´ ιaf x} 2 2 ` λ }x}1
)
, (10)
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
min aPSp´1,xPRm
”
ΨABLpa,xq .“ 12 }x} 2 2 ´ xιaf x,yy ` }y} 2 2 ` λ }x}1
ı
. (11)
We refer to (11) as the Approximate Bilinear Lasso formulation, and it is quite easy to see that ΨABLpa,xq « ΨBLpa,xq when }af x}2 « }x}2, i.e. if a is shift-incoherent, or µpaq « 0. The marginalized objective function ϕBLpaq .“ minx ΨDQpa,xq now has the closed form expression
ϕABLpaq .“ ´12 }softλ rqaf ys} 2 2 . (12)
Here soft denotes the elementwise soft-thresholding operator softtpxiq “ signpxiq ¨maxp|xi| ´ t, 0q, and qa denotes the adjoint kernel of a, i.e. the kernel s.t. xιaf u,vy “ xu, qaf vy @u,v P Rm.
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 a0. This language is made precise in the form of the subsphere
SI .“
ř `PI α`s` rιa0s : α` P R ( Ş Sp´1, (13)
spanned by a small set of cyclic shifts of ιa0. 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θp0, where θ is the probability that any individual entry of x0 is nonzero. Suppose we have a « ř
`PI α`s` rιa0s for some appropriate index set I. Note that if µsa0 « 0, then µsa « 0, @a P SI . Now let αp1q and αp2q be the first and second largest coordinates of the shifts participating in a, and let sp1qra0s and sp2qra0s 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 sp1qra0s.
• If ˇ ˇ
ˇ αp2q αp1q
ˇ
ˇ ˇ « 1, then a is near a saddle-point, with negative curvature pointing towards sp1qra0s
and sp2qra0s. If ˇ ˇ ˇ αp3q αp2q ˇ ˇ ˇ « 0 , i.e. sp1qra0s and sp2qra0s are the only two participating shifts, then ϕABL is also characterized by positive curvature in all orthogonal directions.
• Otherwise, x´gradϕABLpaq, z ´ ay takes on a large positive value, for either u “ sp1qra0s or u “ sp2qra0s, i.e. the negative Riemannian gradient is large and points towards one of the participating shifts.
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 SI must correspond to signed-shifts of a0, this guarantees that the Approximate Bilinear Lasso formulation can be efficiently solved to recover a0 (and subsequently x0) for incoherent a0, as long as a is initialized near some appropriate subsphere and the sparsity coherence tradeoff p0θ Æ pµspa0qq´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 xp0q “ 0, so that our initialization is sparse. For the optimization variable a P Rn, recall from Section 2.2 in the main text that it is desirable to obtain an initialization a0 which is close to the intersection of Sp´1 and a subsphere SI spanned by a few shifts of a0. When x0 is sparse, our measurement y is a linear combination of a few shifts of a0. Therefore, an arbitrary consecutive p0-length window ryi .“ ryi yi`1 . . . yi`p0´1s T of the data y should be not far away from such a subspace SI . As illustrated in Figure 9, one step of the generalized power method (Kuo et al., 2019)
rap0q .“ PSp´1 ` r 0p´1 ; ryi ; 0p´1 s ˘
(14)
ap0q “ PSp´1 ´ ´∇ϕABL ´ rap0q ¯¯
(15)
produces a refined initialization that is very close to a subspace SI spanned by a few shifts of a0 with |I| « θp0. However, (15) is a relatively complicated for a simple idea. In practice, we find that the simple initialization ap0q “ rap0q from (14) works suitably well for solving SaSD with (9).
A.3 COMPARISON TO THE BILINEAR LASSO
Although it is easy to see that ΨABLpaq and ΨBLpaq 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 µpa0q 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 a0. 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 a0 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 a0 become indistinguishable as µpa0q Ñ 1, meaning that one only needs to ensure that a lands in the “global valley” to achieve good estimates of a0 up to symmetry.

heading: B REWEIGHTED SPARSE PENALIZATION
text: When a0 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),
min aPSp´1,xPRm
ΨwBLpa,xq .“ 12 }y ´ af x} 2 2 ` λ }w d x}1 (16)
where w P Rm` and d denotes the Hadamard product. Here we will set the weightsw to be roughly inverse to the magnitude of the true signal x0, i.e.,
wi “ 1
|x0,i| ` ε . (17)
Algorithm 3 Reweighted Bilinear Lasso
Input: Initializations âp0q, x̂p0q, penalty λ ą 0 Output: Local minimizers âpjq, x̂pjq of Ψw pjq
BL . Initialize wp1q “ 1m, j Ð 1. while not converged do
Using the initialization ` âpj´1q, x̂pj´1q ˘ and weight wpjq, solve (16) — e.g. with iADM — to obtain solution ` âpjq, x̂pjq ˘
; Set ε with (19) and update the weights as
wpj`1q “ 1ˇ ˇx̂pjq ˇ ˇ` ε . (18)
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 !
|x|prn{ logpm{nqsq , 10 ´3
)
. (19)
Starting with the initial weights wp0q “ 1m, 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|xi| ` ε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 xi. 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.

heading: C EXTENSION FOR CONVOLUTIONAL DICTIONARY LEARNING
text: 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
the sum of N convolutions between short kernels ta0,kuNk“1 and sparse maps tx0,ku N k“1,
y “ N ÿ
k“1 ιa0,k f x0,k, a0,k P Rp0 , x0,k P Rm, p1 ď k ď Nq. (20)
The natural extension of SaSD, then, is to recover ta0,kuNk“1 and tx0,ku N k“1 up to signed, shift, and permutation ambiguities, leading to the SaS convolutional dictionary learning (SaS-CDL) problem. The SaSD problem can be seen as a special case of SaS-CDL with N “ 1. Based on the Bilinear Lasso formulation in (9) for solving SaSD, we constrain all kernels a0,k over the sphere, and consider the following nonconvex objective:
min takuNk“1, txku N k“1
1
2
›
›
›
›
›
y ´ N ÿ
k“1 ak f xk
›
›
›
›
›
2
2
` λ N ÿ
k“1 }xk}1 , s.t. ak P S
p´1 p1 ď k ď Nq. (21)
Similar to the idea of solving the Bilinear Lasso in (9), we optimize (21) via iADM, by taking alternating descent steps on takuNk“1 and txku N k“1 with the other variable fixed.

heading: D SUPER-RESOLUTION WITH STORM IMAGING
text: For point source localization in STORM frames, recall that we use the SaS model from Section 4.2.2,
Yt lomon
STORM frame
“ ιA0 lomon
point spread function
f X0,t lomon
sparse point sources
` Nt lomon
noise
. (22)
We then apply our SaSD method to recover both A0 and X0,t from Yt. We show our recovery of X0,t as well as the super-resolved image using all available frames in Figure 6 of the main text. Since the main objective of STORM imaging is to recover the point sources, we have deferred the recovered PSFA0 to Figure 13 here.

Review summaries:
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.