Short and Sparse Deconvolution --- A Geometric Approach
Authors: Yenson Lau
Qing Qu
Han-Wen Kuo
Pengcheng Zhou
Yuqian Zhang
John Wright
Venue:  (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 \citep{zhang2017global,kuo2019geometry}, which characterize the optimization landscape of a particular nonconvex formulation of SaSD. This is used to derive a provable algorithm that exactly solves certain non-practical instances of the SaSD problem. We leverage the key ideas from this theory (sphere constraints, data-driven 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.

Published as a conference paper at ICLR 2020
SHORT AND SPARSE DECONVOLUTION A GEOMET- RIC APPROACH
Yenson Lau Electrical Engineering Columbia University y.lau@columbia.co
Qing Qu*
Han-wen Kuo
Center for Data Science New York University qq213@ @ nyu.edu
Electrical Engineering
Columbia University hk2673 coluubia.edu
Pengeheng Zhou
Yugian Zhang
John Wright
Department of Statistics Columbia University
Electrical & Computer Engineering
Electrical Engineering
Rutgers University
Columbia University jw2966 weto6//@oolumbb...e. edu
zhoupe201 @ @@gmail.com
yqz.zhang @ rutgers.edu edu
ABSTRACT
Short-and-sparse deconvolution (SaSD) is the problem of extracting localized, recurring motifs in signals with spatial 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 natu- ral 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 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, data- driven initialization) to develop a practical algorithm, which performs well on data arising from range of application posed by the ill-conditioning of real SaSD problems, and suggest heuristics (accel- eration, continuation, reweighting) mitigate them. Experiments demonstrate the performance and generality of the proposed method.
We highlight key additional challenges
INTRODUCTION
Many signals arising in science and engineering can be modeled as superpositions of basic, recurring motifs, which encode critical information about physical process of interest. Signals of this type can be modeled as the convolution of a zero-paddes zero-padded short kernel ag E IPPo (the motif) with longer sparse signal æg E I (m Po) which encodes the locations of the motifs in the sample':
y Lag
(1)
We term this short-and-sparse (SaS) model. Since often only y is observed, shoor--add-pprse deconvolution (SaSD) is the problem of recovering both ag and Dg from y. Variants of SaSD arise in areas such as microscopy (Cheung et al., 2018), astronomy (Briers al., 2013), and neuroscience (Song et al., 2018). 8). 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 ag (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
$ YI. and QQ contributed equally this work. The full version of this work can be found https: / a aarxii..a v.
www.ii.org/bbb....oo. wwwg/bg./////ooooooo 908. 10959. For simplicity, (1) uses cyclic convolution; algorithms are results also apply to linear convolution with minor modifications. Here denotes the zero padding operator.

---

Published as a conference paper at ICLR 2020
formulations for SaSD and show that the landscape is largely driven by the problem symmetries of SaSD. They derive provable methods for idcalized problem instances, which exactly recover (ao, zo) up to trivial ambiguities. While inspiring, these methods not practical and perform poorly real problem instances. Where the emphasis of Zhang et al. (2017) and Kuo et al. (2019) is on theoretical guarantees, here focus 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. critical issue in moving deconvolution problems: we
from theory to practice is the poor conditioning of show how address this with combination of ideas from sparse optimization, such as momentum, continuation, and reweighting. The end result is general purpose method, which we demonstrate on data from neural spike sorting, calcium imaging and fluorescence microscopy.
The zero-padding operator is denoted by L: RP - R Projection of vector v R2 onto the sphere is denoted by Pa.-1 (v) v/ |12. and P= (= e (v, x> z denotes projection onto the tangent space of z SP-1 The Riemannian gradient of a function f S R at point on the sphere is given by grad f(z) = Px (V)())).
Notation.
Reproducible research.
The code for implementations of our algorithms can be found online:
https: http:/ttuu...mmmiiggegtot/ htttpp//gggttaao.cccc////ggg///p htw.:pprrreeeecoooollo....ooo.o.
For more details of our work on SaSD, we refer interested readers to our project website
https:://eeconvll github.
SYMMETRY AND GEOMETRY IN SASD
In this section, we begin by describing two intrinsic properties for SaSD. Later, 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 admits multiple equivalent solutions. This is purely due to the cyclic convolution between dg and æp, which exhibits the trivial ambiguity²
y
10g To
(ase [lap]) (>
for any nonzero scalar a and cyclic shift 5t I-]. These scale and shift symmetries create several acceptable candidates for ao and Fo, and in the absence of further information we only expect to recover ag and Vy 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 ag, its difficulty also depends highly on the "similarity" of the an to its shifts. Here we capture this notion using the shift-coherence of do.
pl (n)
max Kla0, se [<ao]>| E [0, 1]. 240
Intuitively, the shifts of ag become closer together as pl(a) increases (Figure 10), making objective landscapes for optimization less favorable for recovering any specific shift of ag.
2.1 LANDSCAPE GEOMETRY UNDER SHIFT-INCOIERRENC
A natural approach to solving SaSD is formulate it 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
æ. plus -norm sparsity penalty x.
min
WEL (a, x) y ca c 2| + x
(3)
Later in this section, will see that the kernel length p should be set slightly larger than Po-
The Bilinear Lasso is a nonconvex optimization problem, as the shift symmetries of SaSD create dis- crete local minimizers in the objective landscape. The regularization created by problem symmetries
We therefore assume w.l.o.g. that |aol2 = in this paper.
2

---

Published as a conference paper at ICLR 2020
(a) Near one shift
(b) Two shifts
(c) Multiple shifts
FABL
(e) PUL
Figure I: Geometry of sPArI near superpositions shifts of as (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 μta (aa)) O, the Bilinear Lasso (a)((a min, VaL(a, (a. x) and ABL (aAL (a) are empirically similar the span three shifts.
in nonconvex inverse problems are fairly general phenomenon (Sun 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) objective* AABL which satisfies
Y ABL x) VBL x),
when pla) O.
This non-practical objective serves as valid simplification of the Bilinear Lasso for analysis when the true kernel is itself incoherent, i.e. u(an) = 0 (Figures Id and le). Under its marginalization'
(a) = min WaBL(a.C (a.x),
(4)
certain crucial properties regarding its curvature be characterized for generic choices of æc The reason we choose to partial minimize æ instead of a is because (i) the problem (4) is convex w.r.t. a and (ii) the dimension of the subspace of significantly smaller than that of æ (i.e., p & m), which is the place that the measure concentrates.
Curvature in the span of a few shifts. Suppose we set > Po. which ensures that we can find an = 01 [ao] + [25t,[(a)
gP-1 that lies near the span of two shifts of ag- If a1 = #1 (or 02 > 0) then, under suitable conditions on ag and æn, Kuo et al. (2019) asserts that a lies in strongly convex region of PABL, containing single minimizer near se, [a.] (Figure la); the converse is also true. A saddle-point exists nearby when 0₁ = 022 is balanced, characterized by large negative curvature along the two shifts and positive curvature orthogonal directions (Figure lb). Interpolating between these two cases, large negative gradients point towards individual shifts.
The behavior of PABI. between two shifts of ag - strong convexity near single shifts, and saddle- points near balanced points extends to regions of the sphere spanned by several shifts (Figure Ic); we elaborate on this further in Appendix A.I. This regional landscape guarantees that ag can be efficiently recovered up signed shift using methods for first and second-order descent, as soon 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 is enforced to break the scaling symmetry between and Zy. Choosing the ln-nomm however, has surprisingly strong implications for optimization. The ABI objectivee for example, is piecewise concave whenever is sufficiently far away from any shift of ag but the sphere induces positive curvature near individual shifts to create strong convexity. These two properties combine to ensure recoverability of ag. In contrast, enforcing f--ormm constraints often leads to spurious minimizers for deconvolution problems (Levin et al., 2011; Benichoux et al., 2013; Zhang et al., 2017).
Initializing a few shifts. The landscape of PABL. makes single shifts of as easy locate if initialized near a span of a few shifts. Fortunately, this a relatively simple matter in SaSD, as is
3 "s the intention here apply some key intuition from ABL ohjective towards the Bilinear Lasso itself, we intentionally omit the concrete form of V VABL(a) Readers may refer to Appendix A for more details. Minimizing PABL this equivalent minimizing ADL as a be recovered via convex optimization.
3

---

Published as a conference paper at ICLR 2020
itself sparse superposition of shifts. Therefore, initialization strategy is to randomly choose length-po window D. [t] Vi+1--. ... po from the observation and set
a a Ps Ppp-1 y: Pp-
This brings (1) suitably close to the sum of few shifts of ao (Appendix A.2); any truncation effects are absorbed by padding the ends of y/s which also sets the length for a to be 3pu 2.
Implications for practical computation. The (regionally) benign optimization landscape of PARI guarantees that efficient recovery possible for SaSD when ag is incoherent, Applications of sparse deconvolution, however, are often motivated by sharpening resolution tasks (Huang et al., 2009; Candès & Srmaandz-Grradda 2014; Campisi & Egiazarian, 2016) where the motif ag is smooth and coherent (i.e. μ(ao) is large). The ABL objective is 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 Ie, empirically similar when du is incoherent. Although this breaks down in the coherent case, as we illustrate in Appendix A.3, the symmetry breaking properties of yPn remain present. This allows us to apply the geometric intuition discussed here to create an optimization method that, with the help of number of computational heuristics, performs well in for SaSD in general problem instances.
that -iwddimensional subspheres spanned by shifts of a,
Algorithm Inertial Alternating Descent Method (iADM)
Input: Initializations (0) ga-1, E R", observation é R""; penalty À > 0; momentum e [0, 1).
Output: (a (a) (x())), local minimizer of VBL.
Initialize a (0) (x) (D)
1,2, until converged
Update æ with accelerated proximal gradient step:
w (k) - x (4) + a (k) x (1)
(k+1) - soft w (k) - 4. Vaúi(a w ))
where (v) = sign(v) O max(/v - Al 0) denotes the soft wot-threesholdng operator. Update a with accelerated Riemannian gradient step:
(k) - (a + (ath) 175 Park 1) (a
(k+1) Per- (2 To grad, via (x(2(1) ()+
end for
(a) Gradient descent
(b) CD with momentum (C,.)
wm
GG
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.
3 DESIGNING A PRACTICAL SASD ALGORITHM
Several algorithms for SaSD-ty SaSD-type problems have been developed for specific applications, such as image deblurring (Levin et al., 201 I; Briers et al., 2013; Campisi & Egiazarian, 2016), neuroscience (Rey et al., 2015; Friedrich al., 2017; Song et al., 2018), and image super-resolution (Baker Kanade, 2002; Shtengel et al. 2009; Yang al., 2010), (Wipf& & Zhang, 2014; Ling & Strohmer, 2017; Walk et al., 2017).
augmented with additional structure
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, minimize the Bilinear Lasso (3) instead of the ABL (4) to more accurately account for interactions between shifts of an in highly shift-coherent settings. Furthermore, we also address the negative effects of large coherence using number of heuristics, leading to efficient algorithm for SaSD.
4

---

Published as a conference paper at ICLR 2020
Momentum acceleration. In shift-coherent settings, the Hessian of BL becomes ill-conditioned" near shifts of ag 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(z) with stepsize 7 with the term w,
(k) z + (2(1) - (k-11)
(k+1) ((k)
Here, a controls the momentum added 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 VnL (Alg (Algorithm 1), which modifies iPALM (Pock & Sabach, 2016) to perform updates on via retraction on the sphere (Absil et al., 2009)".
Algorithm 2 SaS-BD BD with homotopy continuation
Input: Observation y RT.. motif size po: momentum a [0, 1); initial ( final X penalty decrease 7 e (0, 1): precision factor 8 (0, 1).
Output: Solution path {(a(m) (,)); X ( for SaSD.
Ser number of iterations - log(AA 1) logn].
Initialize ( RP u (5), (0) 0 R.
for n = 1,..., N do Minimize V (n) precision (2)m) with Algorithm 1:
(a (n)) iADM( (111)) a(n (n11)) 11),(0),))
Update X (m+ < m))(m))
end for
(a) x 10-1
(b) 5
(c) A = 5 10
Figure 3: Bilinear-lasso objective ęP) the sphere 3P-1, "or 3 and varying A; brighter colors indicate higher values. The function landscape yox flattens as sparse penalty X decreases from left right.
Homotopy continuation. It also possible to improve optimization by modifying the objective VBL directly through the sparsity penalty X. Variations of this idea appear in both Zhang et al. (2017) and Kuo 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 X encourage sparser solutions for æ. Conversely, smaller choices of place local minimizers of the marginal objective (a) - min (a, ()) closer to signed-shifts of ag by emphasizing reconstruction quality. When (aa)) is large, however, PBL becomes ill-conditioned conditioning of ag. leading to severe flatness near local minimizers and the creation spurious local minimizers when noise is present (Figure 3). Conversely, larger values of limit a to a small set of support patterns and simplify the landscape ofy/L. at the expense of precision.
A 0 due the poor spectral
is therefore important both for fast convergence and accurate recovery for be chosen appro- priately. When problem parameters such noise level, po. or 0 - are not known a priori, homotopy continuation method (Hale et al., 2008; Wright et al., 2009; Xiao & Zhang, 2013) can be used obtain a range of solutions for SaSD. Using initialization (5), a rough estimate (a(1), (1)
*This is because the circulant matrix Cap is ill-conditioned.
'Setting 0 removes momentum and reverts standard gradient descent. 7 'The stepsizes 4 and nu are obtained by backtracking (Nocedal & Wright, 2006; Pock & Sabach, 2016) to ensure sufficient decrease for V BL (a (8)) (n)(a!)) .n (k++)))) and vice versa.
5

---

Published as a conference paper at ICLR 2020
is obtained by solving (3) with iADM using a large choice for ( This estimate is refined via solution path {(((m) , { by gradually decreasing X By ensuring that a remains sparse along the solution path, the objective VEL enjoys restricted strong convexity w.n.t both a and æ 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 discretized from the Gaussian window function Og- 9 where gp.r = Par-1 ([ [[exp(
(2i-p-1)* 1)²
1
Incoherent kernels ag- Unif({Po-1) 1) are sampled uniformly on the sphere.
(a) Incoherent dg (b) Coherent ag
(c) Incoherent doo
(d) Coherent dg
2004
APM
Sonelay 4DM
Figure 4: Synthetic experiments for Bilinear Lasso. Success probability (a, b): Zo midd (R( 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 au, and vice versa. Success with fixed sparsity rate is more likely when au is incoherent. Algorithmic convergence (c, d): iterate convergence for iADM with an = (k- - 1)/(k + 1) vs. convergence rate, and momentum improves convergence when aa is coherent.
0 (ADM); with and without homotopy. Homotopy significantly improves
(a) Simulated kernel recovery
(c) Real calcium signal vs. reconstruction
0.3
Groundrrsh
teay data
Planosonnc
0.2
Reureee
Flurrescnnee
nace dracted AABM
mecerntrettel tounighigg
0.
-0.1
0.8
80
100
120
140
160
180
Times [s]
Times
(b) Spike train estimates (simulated)
(d) Spike train estimates (real data)
30
cemereaneel of JADM
Floveacone
recovered IXDM
recoverdd www.gging
Floxeecuns
soaeeed of iADM- fuing
0.5
10
12
14
16
18
Tinks [5]
Times [s]
Figure 5: Deconvolution for calcium imaging using Algorithm 2 with iADM and with reweighting (Ap-
pendix B). Simulated data: (a) recovered AR2 kernel: (b) estimate of spike train. Real data: (c) reconstructed calcium signal (d) estimate of spike train. Reweighting improves estimation quality in each case.
Recovery performance. We test recovery probability for varying kernel lengths Pu and sparsity rates 0. To ensure the problem size is sufficiently large, set 100pp. For each po and 6 randomly generate* x iidd. (0) for both coherent and incoherent ag. We solve ten trials of (3) on clean observation data ag By using iADM with X 10 The probability of recovering signed via0
(0) denotes the Bernoulli- Bernoulli-Raeemache Rademacher distribution, which has values +1 w.p. 0/2 and zero w.p. 8.

---

Published as a conference paper at ICLR 2020
shift of dg 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, test the performance of Algorithm I with momentum k-1, see Pock & Sabach (2016)) and without (ax - 0). This is done by minimizing Vn. with (ax K-22
initialization (5), using clean observations with po = 102, m = 10 and 0 = Po 3//4 for coherent and incoherent ag. We also apply homotopy (Algorithm 2) with A = max [r[[a[)], أ - see Xiao & Zhang (2013), X v'px 7) = 0.8, and 8 0.1. The final solve of (3) uses precision regardless of method. Figures 4c and 4d show the comparison results on coherent problem settings.
10
Comparison to existing methods. Finally, we compare iADM, and iADM with homotopy, against a number of existing methods for minimizing PBL. The first is alternating minimization (Kuo et al., 2019), which at each iteration & minimizes a (k) with (()) 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 ag instead of the unit sphere.
For each method, we deconvolve signals with Po = 50, m 100po, and Pu and incoherent ag- For both iADM, iADM with homotopy, and iPALM we set
3/44 for both coherent 0.3. For
homotopy. we set (1) = max {sæ [a 11 y>l. X 0.3 and 5 = 0.5. Furthermore we set 17 = 0.5 v
or 17 0.8 and for ADMM, we set the slack parameter p 0.7 p 0.5 for incoherent and coherent dg 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.
(a) Incoherent do
(b) Coherent Go
(c) Calcium AR2
losu (min ((0,, e/A))
-0.6
JAUIE
IPALM
00.2
SADSE AADSlaowwy APALL
3.5
ADMB
-0.5
-1.2
0.5
3.0
-1.4
FFT operations
ogerations
PPT onerations
Figure 6: Algorithmic comparison. (a) Convergence of various methods minimizing DaL with incoherent do
over FFT operations used (for computing convolutions). The y-axis denotes the log of the angle between a (k) and the nearest shift of au. and each marker denotes five iterations. (b) Convergence for coherent ap, and (c) with AR2 kernel for modeling caleium signals.
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 mieroscopy (STORM) (Rust et al., 2006: Huang et al., 2008; 2010). superresolution technique for in vivo microscopy".
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 convolution between the transient ao and the spike train æo- Whilst ag and 2g 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 ao. a shift- coherent kernel that is challenging for deconvolution, see c.g. Friedrich et al. (2017). We set Ip Lid. Bernoulli(p) 4/53 E IR 10 with additive noise n Lidd (0,5 10-2). Figures 5a and 5b demonstrate accurate recovery of ag and To 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
OOtee superresolution methods for microscopy include photoactivated localization microscopy (PALM) (Betzig et al., 2006), and fluorescence photoactivation localization microscopy (fPALM) (Hess et al., 2006). OOtainee at http:/ h/tp:///iiiii finefrinerrrr codeneuro. org.
7

---

Published as a conference paper at ICLR 2020
decent performance, in the presence of large noise estimation quality 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-ar methods.
(a) Frame 100, time- = 4s
(b) Frame 200, time 8s
(c) Original
(d) Resolved
Figure 7: SaSD for STORM imaging, b) Individual frames (left) and predicted point process map using
SaSD (right). (c, d) shows original microscopy and super-resolved image obtained by our method.
(a) Calcium image Y
(c) Reconstruction A + X (k 1,2)
(d) Predicted activation maps X
(h) Estimated kernels A
Figure 8: Classification of calcium images. (a) Original calcium image; (b) respective kernel estimates; (c) reconstructed images with the (left) neuron and (right) dendrite kernels; (d) respective occurence map estimates.
Super-resolution for fluorescence mieroscopy. 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 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 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 super-resolution microscopy image with nanoscale resolution. For each image frame, the localization task can be formulated via the SaS model
LAo
X 0,4
N noise
STORM frame
spread Function
sparse point sources
where [)] denotes 2D convolution, Here we will solve this task on the single-molecule localization microscopy (SMLM) benchmarking dataser via SaSD, recovering both the PSF Ao and the point source maps Xo,r simultaneously. We apply iADM with reweighting (Appendix B) on frames size 128 x 128 from the video sequence "Tubulin"" each pixel is of 100mm² resolution the fluorescence wavelength is 690mm, 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 produce higher resolution microscopy images.
Localization in calcium images. Our methods easily extended 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
can be accessed at http: h/tp:///ttt //bi gwww. epf1 ch/ sm m/datasetee indx. htm1. Here solve SaSD the same 128 x 128 grid. In practice, the localization problem is solved on finer grid, that the resulting resolution reach 20 30
8

---

Published as a conference paper at ICLR 2020
background noise Pnevmatikakis al. (2016); Giovannucci et al. (2019). We consider frames video obtained via the two-photon calcium microscopy dataset from the Allen Institute for Brain Science 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, A.1 + X1, + 1.2g W X2,t where each summand consists of neurons dendrites exclusively. By extending Algorithm 2 to recover each of the kernels A and maps Xx., we can solve this convolutional dictionary learning (SaS-CDL; see Appendix C) problem which allows 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 denoising or analysis tool for calcium imaging videos provides a very promising direction for future research.
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 heuristies, 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 general bascline for studying and developing extensions to SaSD, such SaS-CDL (Bristow & Lucey, 2014; Chun & Fessler, 2017; Garcia- Carcia-Crddonaaaa & Wohlberg, 2018), Bayesian approaches (Babacan al., 2008; Wipf & Zhang, 2014), and hierarchical SaS models (Chen et al., 2013).
ACKNOWLEDGMENTS
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-Grande Ruoxi Sun, and Liam Paninski for fruitful discussions.
REFERENCES
Pierre-Antoine, Absil, Robert Mahoney, and Rodolphe Sepulchre. Optimization Algorithmms Matrix Manifolds. Princeton University Press, 2009.
Alekh Agarwal, Sahand Negahban, and Martin J Wainwright. Fast global convergence rates of gradient methods for high- digh-dimensional statistical recovery. In Advances in Neural Information Processing Systems, pp. 37-45, 2010.
S Derin Babacan, Rafael Molina, and Aggelos K Katsaggelos. Variational bayesian blind decon- volution using total variation prior. IEEE Transactions on Image Processing, 18(1):12-26. 2008.
Simon Baker and Takeo Kanade. Limits on super-resolution and how to break them. IEEE Transac- lions on Pattern Analysis and Machine Intelligence, 24(9):1167-1183, 2002.
Amir Beck and Mare Teboulle. A fast iterative shrinkage-thressholdnng algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183-202, 2009.
Alexis Benichoux, Emmanuel Vincent, and Rémi Gribonval. A fundamental pitfall blind de- convolution with sparse and shift-invariant priors. In ICASSP-38th International Conference Acoustics, Speech, and Signal Processing- 20113, 2013.
Eric Betzig, George H Patterson, Rachid Sougrat, O Wolf Lindwasser. Scott Olenych, Juan S Bonifacino, Michael W Davidson, Jennifer Lippincott-Schwar and Harald F Hess. Imaging intracellular fluorescent proteins nanometer resolution. Science,
2006.
Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine learning, 3(1):1-122, 201
Jonathan Eckstein, et al. Distributed optimization
Ohttaine at http://observatoory brai inninn-ma p--ap--or www.org/vv.. /vv sual coding/.

---

Published as a conference paper at ICLR 2020
David Briers, Donald D Duncan, Evan R Hirst, Sean J Kirkpatrick, Mareus Larsson, Wiendelt Steen- bergen, Tomas Stromberg, and Oliver Thompson. Laser speckle contrast imaging: theoretical and practical limitations. Journal of biomedical optics, 18(6):060118, 2013.
Hilton Bristow and Simon Lucey. Optimization methods for convolutional sparse coding. arXiv preprint 1406vv:400.222 2407, 2014.
Patrizio Campisi and Karen Egiazarian. Blind image deconvolution: theory and applications. CRC press, 2016.
Emmanuel Candès and Carlos Grandad Towards mathematical theory of super- resolution. Communications on pure and applied Mathematics, 67(6))90-9956, 2014.
Emmanuel J Candes, Michael B Wakin, and Stephen P Boyd. Enhancing sparsity by reweighted & minimization. Journal of Fourier analysis and applications, 14(5-6):877-905, 2008.
Bo Chen, Gungor Polatkan, Guillermo Sapiro. David Blei, David Dunson, and Lawrence Carin. Deep learning with hierarchical convolutional factor analysis. IEEE transactions on pattern analysis and machine intelligence, 35(8):1887-190) 2013.
Sky Cheung, John Shin, Yenson Lau, Zhengyu Chen, Ju Sun, Yuqian Zhang, John Wright, and Abhay N Pasupathy. Dictionary learning in fourier transform scanning tunneling spectroscopy. arXiv preprint arXiv: 10752, 2018.
Il Yong Chun and Jeffrey A Fessler. Convolutional dictionary learning: Acceleration and convergence.
IEEE Transactions on Image Processing, 27(4):1697-172,, 2017.
Chaitanya Ekanadham, Daniel Tranchina, and Eero P Simoncelli. A blind sparse deconvolution method for neural spike identification. In Advances in Neural Information Processing Systems, pp. 1440-1448, 2011.
Johannes Friedrich, Pengeheng Zhou, and Liam Paninski. Fast online deconvolution of calcium imaging data. PLoS Computational Biology, 1305)2e00 2017.
Cristina
and Brendt Wohlberg, Convolutional dictionary learning: A comparative review and new algorithms. IEEE Transactions on Computational Imaging, 4(3):366-38,, 2018.
Andrea Giovannucci, Johannes Friedrich, Pat Jeremie Kalfon, Brandon L. Brown, Sue Ann Koay, Jiannis Taxidis, Farzaneh Najafi, Jeffrey L. Gauthier, Pengcheng Zhou, al. Caiman an open source tool for scalable calcium imaging data analysis. Elife, 8:e38173, 2019.
Christine Grienberger and Arthur Konnerth. Imaging calcium in neurons. Neuron, 73(5):862-885. 2012.
Elaine T Hale, Wotao Yin, and Yin Zhang. Fixed-point continuation for \ell_ 1-minmmizatiinn: Methodology and convergence. SIAM Journal on Optimization, 19(3): 100-)11301110 2008.
Samuel T Hess, Thanu PK Girirajan, and Michael D Mason. Ultra-high resolution imaging by fluorescence photoactivation localization microscopy. Biophysical journal, 9((11):4258-427 2006.
Bo Huang, Wenqin Wang, Mark Bates, and Xiaowei Zhuang. Three-dimennional super-resolution imaging by stochastic optical reconstruction microscopy. Science, 319(566))8100-888 2008.
Bo Huang, Mark Bates, and Xiaowei Zhuang. Super-resolution fluorescence microscopy. Annual Review of Biochemistry, 78:993-1016, 2009.
Bo Huang, Hazen Babcock, and Xiaowei Zhuang. Breaking the diffraction barrier: super-resolution imaging of cells. Cell. 143(7):0447-155; 2010.
Chi Jin, Rong Ge, Praneeth Netrapalli, Sham M Kakade, and Michael I Jordan. How to escape saddle points efficiently. Proceedings of the 34th International Conference on Machine Learning, pp. 1724-1732, 2017.
10

---

Published as a conference paper at ICLR 2020
Chi Jin, Praneeth Netrapalli, and Michael I Jordan. Accelerated gradient descent escapes saddle points faster than gradient descent. In Conference On Learning Theory, pp. 1042-1085, 2018.
Han-Wen Kuo, Yuqian Zhang, Yenson Lau, and John Wright. Geometry and symmetry in short-and- sparse deconvolution. In International Conference on Machine Learning (ICML), June 2019.
Anat Levin, Yair Weiss, Fredo Durand, and William Freeman. Understanding blind deconvolution algorithms. IEEE Transactions on Pattern Analysis and Machine Intelligence, 331223)23333))) 2011.
Shuyang Ling and Thomas Strohmer. Blind deconvolution meets blind demixing: Algorithms and performance bounds. IEEE Transactions on Information Theory, 63(7):497-44500 2017.
Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2013.
Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science & Business Media, 2006.
Eftychios A Pnevmatikakis, Daniel Soudry, Yuanjun Gao, Timothy A Machado, Josh Merel, David Pfau, Thomas Reardon, Yu Mu, Clay Lacefield, Weijian Yang, et al. Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron, 89(2):285 220)2999999, 2016.
Thomas Pock and Shoham Sabach. Inertial proximal alternating linearized minimization (ipalm) for nonconvex and nonsmooth problems. SIAM Journal on Imaging Sciences, 94):1756-1777, 2016.
Boris T Polyak. Some methods of speeding up the convergence of iteration methods. USSR
Computational Mathematics and Mathematical Physics, 4(5):1-17, 1964.
Hernan Gonzalo Rey, Carlos Pedreira, and Rodrigo Quian Quiroga. Past, present and future of spike sorting techniques. Brain Research Bulletin, 119:106-117, 2015.
Michael J Rust, Mark Bates, and Xiaowei Zhuang. Sub-diffraction-lini imaging by stochastic optical reconstruction microscopy (storm). Nature Methods, 3(10):793, 2006.
Gleb Shtengel, James A Galbraith, Catherine G Galbraith, Jennifer Sippincott-Schwart Jennifer M Gillette, Suliana Manley, Rachid Sougrat, Clare M Waterman, Pakorn Kanchanawong, Michael W Davidson, al. Interferometric fluorescent super-resolution microscopy resolves 3d cellular ultrastructure. Proceedings of the National Academy of Sciences, 106(9):3125-3133 2009.
Andrew H Song, Francisco Flores, and Demba Ba. Spike sorting by convolutional dictionary learning. arXiv preprint 18Xiv:1806.00 01999, 2018
Ju Sun, Qing Qu, and John Wright. When are nonconvex problems not scary? 15.0 06096, 2015.
arXiv preprint
Philipp Walk, Peter Jung, Götz Pfander, and Babak Hassibi. Blind deconvolution with additional 04499, 2017.
autocorrelations via convex programs. arXiv preprint
David Wipf and Haichao Zhang. Revisiting bayesian blind deconvolution. The Journal of Machine Learning Research, 15(1)3559-3633, 2014.
Stephen J Wright, Robert D Nowak, and Mário AT Figueiredo. Sparse reconstruction by separable approximation. JEEE Transactions on Signal Processing, 57(7):2479-244 24933 2009.
Lin Xiao and Tong Zhang. A proximal-gradient homotopy method for the sparse least-squares problem. SIAM Journal Optimization, 1362):1062-1011, 2013. 2013.
Jianchao Yang, John Wright, Thomas S Huang, and Yi Ma. Image super-resolution via sparse representation. IEEE Transactions Image Processing, 19(11):2861 19(11)2861-2873, 2010.
Florence Yellin, Benjamin D Haeffele, and René Vidal. Blood cell detection and counting in holographic lens-free imaging by convolutional sparse dictionary learning and coding. In IEEE 14th International Symposium Biomedical Imaging, pp. 650-653. IEEE, 2017.
11

---

Published as a conference paper at ICLR 2020
Yuqian Zhang, Yenson Lau, Han-Wen Wen Kuo, Sky Cheung. Abhay Pasupathy, and John Wright. On the global geometry of ctbhere-oonsraineet sparse blind deconvolution. In Computer Vision and Patter 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 multispectral convolutional sparse coding. In Proceedings of the IEFE Conference Computer Vision and Pattern Recognition, pp. 308 3081-3088, 2014.
A APPROXIMATE BILINEAR LASSO OBJECTIVE
Recall from Section 2.2 of the main text that SaSD can be formulated the Bilinear Lasso problem
min mm
(VBL(, x:) b aa + 201 X
|
Unfortunately, this objective is challenging for analysis. A major culprit is that its marginalization
(a) = mmin ly
x ||22|||
(10)
generally does not admit closed form solutions due the convolution with a in the squared error term. This motivates Kuo al. (2019) to study the nonconvex formulation
min
MaBL x) = 1 (va @ x. + 1B12 A |+||]]
(11)
dESP- mER
We refer to (11) as the Approximate Bilinear Lasso formulation, and it is quite easy to see that V VABI(a, x) & (a, when la @ æ= | i.e. if a is shift-incoherent, or u(a) s The marginalized objective function PaL (a) = min (a x)) now has the closed form expression
(a) = - |soft, à y)1
(12)
Here soft denotes the elementwise soft-thresholding operator soft, sign(z;) and à denotes the adjoint kemnel of a, i.e. the kernel s.t. (La @ u. u) = (u, a VE
A.1 LANDSCAPE GEOMETRY
The rest of Section 2.2 discusses the regional characterization of PABL in the span of small number of shifts from ag. This language is made precise in the form of the subsphere
ST = EIET 048; [1ao] ag E R
SP
(13)
spanned by a small set of cyclic shifts of cag- Although we will not discuss the explicit distance function here, the characterization by Kuo et al. (2019) holds whenever is close enough to such a subsphere with [Z] & 40po. where 9 is the probability that any individual entry of Zo is nonzero. Suppose have a Eiez asse [caa] for some appropriate index I. Note that if 14a0 x then ,,a & O, Va E ST. Now let ((()) and (((2) be the first and second largest coordinates of the shifts participating in a, and let ((1)[a)) and $t2)[as) be the corresponding shifts. Then
0, then is in strongly convex region of PABL containing single local minimizer corresponding s(y][a].
If
P122 CiE
If (12)) o(1)
R 1, then a is near maddle--oint, with negative curvature pointing towards $(1) [ao]
0. i.e. sca][ao] and s(2) [ao] are the only two participating shifts, then PBBL is also characterized by positive curvature in all orthogonal directions.
and 8(2) [an]. If (:
(2:2)
Otherwise, (-)rad(AALL())) - a) takes on a large positive value, for either 8(1) [an] or = S(2) [q.] i.e. the negative Riemannian gradient large and points towards one the participating shifts.
12

---

Published as a conference paper at ICLR 2020
Data y
Kernel ao
Sparse wo
Truncation a (-1)
Initialization (0)
048; [48;[a0] + 4288 [cg]
Figure 9: Data-driven initialization for using piece of the observed data to generate good initial point (0) Top: data ao @ zo a superposition of shifts of the true kernel at. Bottom: a length-po window contains pieces of just a few shifts. Bottom-center: one step of the generalized power method approximately fills in the missing pieces, yiclding initialization that is close lincar combination of shifts of as (right).
This is an example of a ridable saddle property (Jin al., 2017) that allows many first and second- order methods to locate local minimizers. Since all local minimizers of PABL near Sy must correspond to signed-shifts
ag. this guarantees that the Approximate Bilinear Lasso formulation can be efficiently solved to recover ag (and subsequently ai)) for incoherent do, long is initialized near some appropriate subsphere and the sparsity coherence tradeoff paf E (H,(aa))) (1o2))//² is satisfied. We note that this is poor tradeoff rate, which reflects that the Approximate Bilinear Lasso formulation non-pren con-practical and cannot handle SaSD problems involving kernels with high shift-coherence.
A.2
INITIALIZATION
For the SaS-BD problem, usually initialize x: by (0) 0, that our initialization is sparse. For the optimization variable a R", recall from Section 2.2 in the main text that it is desirable to obtain an initialization a" which is close the intersection of SP-1 and subsphere Sz spanned by a few shifts of ag. When 2p is sparse, our measurement is linear combination of a few shifts of ag- Therefore, data y should be not far away from such a subspace Sz. As illustrated Figure 9, one step of the generalized power method (Kuo et al., 2019)
arbitrary consecutive Po-length window y: = [3/ V++1... ... -1-ppo1] of the
Pay-1 (I 0p-1 V: Op-1 1)
(14)
(0) Par-1 (-V (a())
(15)
produces refined initialization that very close to subspace Sn spanned by few shifts of as with [Z] & Opo. However, (15) is relatively complicated for simple idea. In practice, we find that the simple initialization a (0) (0)) from (14) works suitably well for solving SaSD with (9).
A.3 COMPARISON TO THE BILINEAR LASSO
Although is easy to see that V VBLL(a and VBL (a) are similar as long as plaa & is also clear that these two quantities can be very different when μ((a) is large. This is especially significant when (((a)) 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 ag Compared to the incoherent case, şBL also takes on small values in regions between adjacent shifts, creating "global valley" on the subsphere. Theoretically, this makes it difficult to ensure exact recovery of up to symmetry when ag is coherent, and the objective function becomes much more complicated. This is not significant issue in terms of practical computation, however, since adjacent shifts of ag become indistinguishable as u(an) - 1, meaning that one only needs to ensure that lands in the "global valley" to achieve good estimates of as up symmetry.
13

---

Published as a conference paper at ICLR 2020
(a) SPABL, incoherent
(b) Pat., incoherent
(c) PBL., coherent
Figure 10: Low-dimensional subspheres spanned by shifts of do. Subfigures (a.b) present the optimization
landscapes (PABL (a) and you (a), for 3P [pan (ao, si [a[o [aoo]88(aa0) ) with higher values heing brighter. red dots denote shifts of ag. Subfigure (c) shows landscape FaL when ag is coherent, which significantly departs from the landscapes of (a,b), but still retains symmetry breaking curvature.
B REWEIGHTED SPARSE PENALIZATION
When do is shift-coherent, minimization of the objective VBL with respect to becomes sensitive to perturbations, creating "smudging" effects on the recovered map æ. 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 weighted variant of the original Bilinear Lasso problem from (9).
min m
uP (a, x))
y a A 2
w O a
(16)
where w E and denotes the Hadamard product. Here will set the weights w be roughly inverse to the magnitude of the true signal Zo. i.e.,
(17)
2o.el +
Algorithm Reweighted Bilinear Lasso
Input: Initializations ("), (0) penalty À > 0
Output: Local minimizers a ((i) of U Initialize (1) mm j - 1.
while not converged do
Using the initialization (a (00-11,,(i-11) 20-11) and weight w(), solve (16) e.g. with iADM to obtain solution (a(s), .(:)); Set e with (19) and update the weights as
(
(18)
[b +
Update --+t end while
In addition to choosing À 0, here trades off between sparsification strength (small e) and algorithmic stability (large e). Let [x](s) denote the i-th largest entry of [æ].. For experiments in the main text, we set
g = max { | (m/ log(m/n))) 10
(19)
Starting with the initial weights (0) = 1m, Algorithm successively solves (16), updating the weights using (17) at each outer loop iteration j. As j this method becomes equivalent to replacing the noom in (9) with the nonconvex penalty Σ, log [ + E) (Candes et al., 2008).
We easily adopt our iADM algorithm to solve this subproblem, by taking the proximal gradient on z with a different penalty A; for each entry z;. Figure 11, as well as calcium imaging experiments in Section 4.2, Figure of the main text, demonstrate improved estimation result of this method.
14

---

Published as a conference paper at ICLR 2020
(a) True map æb
(c) Noisy y, & only
(e) Noisy y, reweighted
(b) True motif as
(d) Noisy a. & only
(f) Noisy a, reweighted
Figure 11: Recovery of ma with fi-reweighting, (a, b) Truth signals. (c) Solving ming VaL (a, a)) with noisy data and coherent an leads to low-quality estimates of æ: (d) performance suffers further when a is noisy estimate of ap. (e, f) Reweighted & minimization alleviates this issue significantly.
C EXTENSION FOR CONVOLUTIONAL DICTIONARY LEARNING
1
+
40,1
20,1
4g22
Figure 12: Convolutional dictionary learning. Simultaneous recovery for multiple unknown kernels
{ao, k and sparse activation maps (0o, from y = Σn ao.k c Lo
1=
(a) PSF in 2D
(b) PSF in 3D
Figure 13: Estimated PSF for STORM imaging. The left hand side shows the estimated 8 x 8 PSF 2D, the right hand side visualizes the PSF in 3D.
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 Murcia-Cardona & Wohlberg (2018)), which have broad applications in mieroscopy data analysis (Yellin et al., 2017; Zhou et al., 2014; Cheung et al., 2018) and neural spike sorting (Ekanadham al., 201 1; Rey al., 2015; Song al., 2018). As illustrated in Figure 12, the observation is
15

---

Published as a conference paper at ICLR 2020
the sum of N convolutions between short kernels (au,e and sparse maps ( N
N
Y
LOook
0o,k E
IPpo 20,6
Rm
(1 N).
(20)
The natural extension of SaSD, then, is recover (ao and {an up 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 Based on the Bilinear Lasso formulation in (9) for solving SaSD, we constrain all kernels do,k over the sphere, and consider the following nonconvex objective:
1 (aulé-1. (ar l 2
min
y
Σ R-1
+
s.t. ag E (IEREN). (21)
Similar to the idea of solving the Bilinear Lasso in (9), we optimize (21) via iADM, by taking alternating descent steps (ax) and { with the other variable fixed.
D SUPER-RESOLUTTIN WITH STORM IMAGING
For point source localization in STORM frames, recall that the SaS model from Section 4.2.2,
LA
X
N
(22)
STORM frame
point spread function
sparse point SOUICES
noise
We then apply our SaSD method to recover both Ao and Xo from Y. We show our recovery of Xo, well the super-resolved image using all available frames Figure of the main text. Since the main objective of STORM imaging is to recover the point sources, we have deferred the recovered PSF Ag to Figure 13 here.
16


Reviews:
[{"review_id": "rJlh3raLqH", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper analyzes some of the optimization challenges presented by a particular formulation of \"short-and-sparse deconvolution\" and proposes a new general purpose algorithm to try to alleviate them. The paper is reasonably clearly written and the experimental results seem impressive (as a non-specialist in this area). However the experimental investigation on real data does not compare to a baseline, which seems like an essential part of investigating the proposed method. If this were corrected I would recommend an accept.\n\nI would like to make it clear that I have little experience in this area and am not familiar with relevant previous literature, so I was only able to roughly judge novelty of the ideas and how the experimental results might compare with existing approaches.\n\n\nMajor comments:\n\nFor the experimental investigations on real data there was no baseline presented. Are there no task-specific algorithms that have previously been developed for deconvolution of calcium signals, fluoresence microscopy and calcium imaging? At the very least it would be helpful to compare to one of the other general purpose methods used for figure 6.\n\nIt seemed like a lot of the paper is taken up with a recap of the results presented in Kuo et al. (2019), and it didn't always seem clear exactly what the relevance of these results were to the present paper. In particular it seemed strange to me to devote so much space in section 2 to the ABL objective given that it is not used (as far as I can tell) for the proposed algorithm.\n\n\nMinor comments:\n\nThe abstract says \"We leverage... sphere constraints, data-driven initialization\". Is the effect of these actually investigated experimentally?\n\nIn the abstract, \"This is used to derive a provable algorithm\" makes it sounds to me like this will be done in the current paper.\n\nIn the abstract, a reference for the \"due to the spectral decay of the kernel a_0\" claim would be helpful.\n\nIn the notation section, the definition of the Riemannian gradient seems a little sloppy mathematically: strictly f needs to be defined on $\\reals^p$ (or an open subset of $\\reals^p$ containing $S^{p - 1}$) in order for the right side of the gradient expression to be defined.\n\nAt the start of section 2, it would be helpful to give a one-sentence description of the unifying theme of the section.\n\nIt might be helpful to explicitly state that \"coherence\" is related to the strength of temporal / spatial correlations to aid developing the reader's intuition for the meaning of $\\mu(a)$.\n\nThroughout the paper the problems that multiple equivalent local optima (due to shift symmetry) cause come up repeatedly. What happens if this symmetry is removed straightforwardly by adding a term to the cost function to encourage a particular value of $a$ to be the largest?\n\nIn section 2.1, for \"It analyzes an ABL...\", it's not completely clear what \"It\" refers to (I presume Kuo et. al (2019) from context).\n\nMarginalization in my experience refers to a \"partial summing\" operation, whereas just above (4) it is used to refer to a \"partial maximization\" operation. This seems non-standard to me, but is this usage standard in this field?\n\nI didn't understand the relevance of the sentence \"Under its marginalization... smaller dimension p << m.\" to the present paper. The sentence also seemed vague and difficult to understand if you were not already familiar with this result. It should also have a reference to justify this claim.\n\nIt wasn't clear to me whether figure 1 were schematic \"rough intuition\" diagrams intended merely to be suggestive, or provably showing the type of behavior that happens in all cases. Also, what are the axes, generic \"parameter space\", I guess? I also didn't follow why (a) appears projected on to a plane while (b) and (c) appear projected on to a sphere(?)\n\nIn section 3, under \"momentum acceleration\", it would be helpful to justify the claim that \"In shift-coherent settings, the Hessian... ill-conditioned...\".\n\nIn figure 4, the axis labels are much too small to read. In figure 5 (b), the red poorer results obscure the green better results.\n\nIn figure 4 (d), is it fair to say that the main convergence speed improvement for homotopy-iADM is in going faster from \"quite near\" the optimum to \"really near\" it, rather than getting \"quite near\" it in the first place? If so, isn't the latter more often what's relevant in practical applications?\n\nIn figure 5, it's a little confusing to switch color meanings between (a) and (b).\n\nReweighting seems to have a dramatic positive effect in figure 5. Is it worth investigating its effect in the other experiments as well, for example in figure 4?\n\nIn figure 6 (b), is there any reason not to compare against standard black box optimizers like vanilla SGD, ADAM, etc (possibly with projection to satisfy the sphere constraint if necessary)? Would they perform very badly?\n\nTypo \"Whilst $a_0$ nor $x_0$\" should be \"Whilst $a_0$ and $x_0$\".\n\nWhat does the square-boxed convolution operator in (8) mean?\n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HklZ3uqptH", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "* Summary\n\nThe 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.\n\nI 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.\nOn 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.\n\n* Detailed Comments\n\nWhile 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?\n\nSomehow 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?\n\nIn Fig 1 an explanation of the axis and projections would be needed.\nWhile 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.\n\nWhat are similarities / differences to related reconstruction problems such as non-negative matrix factorization with sparseness constraints?\n\n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: 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).

The paper extends the work by Kuo et al. (2019) by providing a practical algorithm for solving those inverse problems. A focus of the paper is to solve the bilinear lasso instead of the approximate bilinear lasso, because this approximation is poor for coherent problems. Having read the rebuttal and the paper, I believe the authors addressed the issues raised by Reviewer #2 in a sufficient way.

small things:
- it would be good to define $\iota$ (zero-padding operator) in (1)
- it would be good to define $p, p_0$ just below (3). They seem to be appearing out of the blue without any direct relation to anything mentioned prior in section 2.
- it would be good to cite some older/historic references for various optimization methods , e.g. [1] below. 


[1] Richter & deCarlo 
Continuation methods: Theory and applications
IEEE Transactions on Systems, Man, and Cybernetics, 1983
https://ieeexplore.ieee.org/abstract/document/6313131
Decision: Accept (Poster)