Exact Support Recovery for Sparse Spikes Deconvolution

Initial talk at SIAM IS14, updated for the Workshop "IFIP TC7.4 Workshop on Inverse Problems and Imaging", updated for the conference SPARS 2015.

Transcript

1. Gabriel Peyré www.numerical-tours.com Exact Support Recovery for Sparse Spikes Deconvolution Joint work with Vincent Duval & Quentin Denoyelle VISI N
2. Sparse Deconvolution Neural spikes (1D) y m0 = ? ginal Signal m0 t 1 x2 x3 ⇤ Low-pass filter ' t 0.5 0.5 + Noise w 0 1 + y = ' ? m0 + w m0 is “sparse” ' w
3. Sparse Deconvolution Neural spikes (1D) Seismic imaging (1.5D) y m0 = ? ginal Signal m0 t 1 x2 x3 ⇤ Low-pass filter ' t 0.5 0.5 + Noise w 0 1 + y = ' ? m0 + w m0 is “sparse” ' w
4. Sparse Deconvolution Neural spikes (1D) Astrophysics (2D) Seismic imaging (1.5D) y m0 = ? ginal Signal m0 t 1 x2 x3 ⇤ Low-pass filter ' t 0.5 0.5 + Noise w 0 1 + y = ' ? m0 + w m0 is “sparse” ' w
5. Sparse Deconvolution Neural spikes (1D) Astrophysics (2D) Seismic imaging (1.5D) y m0 = ? ginal Signal m0 t 1 x2 x3 ⇤ Low-pass filter ' t 0.5 0.5 + Noise w 0 1 + y = ' ? m0 + w m0 is “sparse” Presented results n D problems extend to ' w
6. Overview • Sparse Spikes Super-resolution • Robust Support Recovery • Asymptotic Positive Measure Recovery
7. Radon measure m on T = ( R/Z ) d . Deconvolution of Measures m
8. Discrete measure: ma,x = P N i =1 ai xi , a 2 RN , x 2 TN Radon measure m on T = ( R/Z ) d . Deconvolution of Measures m a,x m
9. Discrete measure: ma,x = P N i =1 ai xi , a 2 RN , x 2 TN Radon measure m on T = ( R/Z ) d . Deconvolution of Measures m a,x m ' 2 C2(T ⇥ T) y = (m) + w ( m ) = Z T ' ( x, ·)d m ( x ) Linear measurements:
10. y Discrete measure: ma,x = P N i =1 ai xi , a 2 RN , x 2 TN Radon measure m on T = ( R/Z ) d . Deconvolution of Measures m a,x ' m ' 2 C2(T ⇥ T) y = (m) + w ( m ) = Z T ' ( x, ·)d m ( x ) Linear measurements: Example: 1-D ( d = 1) convolution ' ( x, t ) = ' ( t x )
11. y Discrete measure: ma,x = P N i =1 ai xi , a 2 RN , x 2 TN Minimum separation: = mini6=j | xi xj | ! Signal-dependent recovery criteria. Radon measure m on T = ( R/Z ) d . Deconvolution of Measures m a,x ' = 2/fc y = 0.5/fc m ' 2 C2(T ⇥ T) y = (m) + w ( m ) = Z T ' ( x, ·)d m ( x ) Linear measurements: Example: 1-D ( d = 1) convolution ' ( x, t ) = ' ( t x )
12. Discrete `1 regularization: Computation grid z = ( zk) K k=1. Sparse l1 Deconvolution y zk
13. Discrete `1 regularization: Computation grid z = ( zk) K k=1. Basis-pursuit / Lasso: Sparse l1 Deconvolution y zk ¯ : a 2 RK 7! (ma,z) = X k ak'(zk, ·) 2 Im( ) min a2RN 1 2 ||y ¯a||2 + ||a||1
14. Discrete `1 regularization: Computation grid z = ( zk) K k=1. Basis-pursuit / Lasso: Sparse l1 Deconvolution y zk ¯ : a 2 RK 7! (ma,z) = X k ak'(zk, ·) 2 Im( ) min a2RN 1 2 ||y ¯a||2 + ||a||1 Why `1? q = 0 a1 a2 “`0 ball”
15. Discrete `1 regularization: Computation grid z = ( zk) K k=1. Basis-pursuit / Lasso: Sparse l1 Deconvolution y zk ¯ : a 2 RK 7! (ma,z) = X k ak'(zk, ·) 2 Im( ) min a2RN 1 2 ||y ¯a||2 + ||a||1 q = 1 q = 2 q = 3/2 q = 1/2 Why `1? `q ball a 2 RK ; P k |ak |q 6 1 q = 0 a1 a2 “`0 ball” !
16. Discrete `1 regularization: Computation grid z = ( zk) K k=1. Basis-pursuit / Lasso: Sparse l1 Deconvolution y zk ¯ : a 2 RK 7! (ma,z) = X k ak'(zk, ·) 2 Im( ) min a2RN 1 2 ||y ¯a||2 + ||a||1 q = 1 q = 2 q = 3/2 q = 1/2 Why `1? `q ball a 2 RK ; P k |ak |q 6 1 q = 0 a1 a2 sparse “`0 ball” !
17. Discrete `1 regularization: Computation grid z = ( zk) K k=1. Basis-pursuit / Lasso: Sparse l1 Deconvolution y zk ¯ : a 2 RK 7! (ma,z) = X k ak'(zk, ·) 2 Im( ) min a2RN 1 2 ||y ¯a||2 + ||a||1 q = 1 q = 2 q = 3/2 q = 1/2 Why `1? `q ball a 2 RK ; P k |ak |q 6 1 q = 0 a1 a2 sparse convex “`0 ball” !
18. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1
19. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x
20. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x
21. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x
22. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x
23. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: (P0(y)) min m {|m|(T) ; m = y} ! 0+ |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x
24. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: (P0(y)) min m {|m|(T) ; m = y} ! 0+ |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x Proposition: If dim(Im( )) < +1, 9( a, x ) 2 RN ⇥ TN with N 6 dim(Im( )) such that m a,x is a solution to P ( y ).
25. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 ! Algorithms: [Bredies, Pikkarainen, 2010] (proximal-based) [Cand` es, Fernandez-G. 2012] (root finding) min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: (P0(y)) min m {|m|(T) ; m = y} ! 0+ |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x Proposition: If dim(Im( )) < +1, 9( a, x ) 2 RN ⇥ TN with N 6 dim(Im( )) such that m a,x is a solution to P ( y ).
26. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 ! Algorithms: [Bredies, Pikkarainen, 2010] (proximal-based) [Cand` es, Fernandez-G. 2012] (root finding) min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: (P0(y)) min m {|m|(T) ; m = y} ! 0+ Competitors: Prony’s methods (MUSIC, ESPRIT, FRI). |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x Proposition: If dim(Im( )) < +1, 9( a, x ) 2 RN ⇥ TN with N 6 dim(Im( )) such that m a,x is a solution to P ( y ).
27. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 ! Algorithms: [Bredies, Pikkarainen, 2010] (proximal-based) [Cand` es, Fernandez-G. 2012] (root finding) min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: (P0(y)) min m {|m|(T) ; m = y} ! 0+ Competitors: Prony’s methods (MUSIC, ESPRIT, FRI). “+”: always works when w = 0, less sensitive to sign. |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x Proposition: If dim(Im( )) < +1, 9( a, x ) 2 RN ⇥ TN with N 6 dim(Im( )) such that m a,x is a solution to P ( y ).
28. Grid-free Sparse Recovery Grid-free regularization: total variation of measures: |m|(T) = sup R ⌘dm : ⌘ 2 C(T), ||⌘||1 6 1 ! Algorithms: [Bredies, Pikkarainen, 2010] (proximal-based) [Cand` es, Fernandez-G. 2012] (root finding) min m 1 2 || (m) y||2 + |m|(T) (P (y)) Sparse recovery: (P0(y)) min m {|m|(T) ; m = y} ! 0+ Competitors: Prony’s methods (MUSIC, ESPRIT, FRI). “+”: always works when w = 0, less sensitive to sign. “-”: only for convolution operator, '(x, t) = '(x t) |m a,x |(T) = ||a|| ` 1 m a,x |m|(T) = R |f| = ||f||L1 d m ( x ) = f ( x )d x Proposition: If dim(Im( )) < +1, 9( a, x ) 2 RN ⇥ TN with N 6 dim(Im( )) such that m a,x is a solution to P ( y ).
29. Overview • Sparse Spikes Super-resolution • Robust Support Recovery • Asymptotic Positive Measure Recovery
30. (P0(y)) Robustness and Support-stability = 0.55/fc = 0.45/fc = 0.1/fc = 0.3/fc min m {|m|(T) ; m = y} Low-pass filter supp( ˆ ' ) = [ fc, fc]. When is m0 solution of P0( m0) ?
31. (P0(y)) Robustness and Support-stability = 0.55/fc = 0.45/fc = 0.1/fc = 0.3/fc min m {|m|(T) ; m = y} Low-pass filter supp( ˆ ' ) = [ fc, fc]. When is m0 solution of P0( m0) ? Theorem: [Cand` es, Fernandez G.] > 1.26 fc ) m0 solves P0( m0).
32. (P0(y)) Robustness and Support-stability = 0.55/fc = 0.45/fc = 0.1/fc = 0.3/fc min m {|m|(T) ; m = y} Low-pass filter supp( ˆ ' ) = [ fc, fc]. are solutions of P ( m0 + w )? How close to m0 When is m0 solution of P0( m0) ? Theorem: [Cand` es, Fernandez G.] > 1.26 fc ) m0 solves P0( m0).
33. ! [Cand` es, Fernandez-G. 2012] (P0(y)) Robustness and Support-stability = 0.55/fc = 0.45/fc = 0.1/fc = 0.3/fc min m {|m|(T) ; m = y} Low-pass filter supp( ˆ ' ) = [ fc, fc]. are solutions of P ( m0 + w )? ! [Fernandez-G.][de Castro 2012] Weighted L2 error: Support localization: How close to m0 When is m0 solution of P0( m0) ? Theorem: [Cand` es, Fernandez G.] > 1.26 fc ) m0 solves P0( m0).
34. ! [Cand` es, Fernandez-G. 2012] (P0(y)) Robustness and Support-stability = 0.55/fc = 0.45/fc = 0.1/fc = 0.3/fc min m {|m|(T) ; m = y} Low-pass filter supp( ˆ ' ) = [ fc, fc]. are solutions of P ( m0 + w )? ! [Fernandez-G.][de Castro 2012] General kernels? Exact support recovery? Open problems: Weighted L2 error: Support localization: How close to m0 When is m0 solution of P0( m0) ? Theorem: [Cand` es, Fernandez G.] > 1.26 fc ) m0 solves P0( m0).
35. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9
36. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9 = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i
37. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9 m 2 @◆||·||1 61 (⌘) = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i
38. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9 m 2 @◆||·||1 61 (⌘) = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i
39. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9 Ideal low-pass filter: ! ⌘ = ⇤p trigonometric polynomial. m 2 @◆||·||1 61 (⌘) ⌘ = ⇤ m y = ⇤p = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i
40. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9 Ideal low-pass filter: ! ⌘ = ⇤p trigonometric polynomial. m 2 @◆||·||1 61 (⌘) ⌘ = ⇤ m y = ⇤p = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i D (y) = sup || ⇤p||1 61 hp, yi 2 ||p||2
41. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) 9 Ideal low-pass filter: ! ⌘ = ⇤p trigonometric polynomial. m 2 @◆||·||1 61 (⌘) , ⌘ 2 @|m|(T) ⌘ = ⇤ m y = ⇤p = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i D (y) = sup || ⇤p||1 61 hp, yi 2 ||p||2
42. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) ⇢ ⌘ ; 8 t, | ⌘ ( t )| 6 1 8 i, ⌘ ( xi) = sign( ai) @|m a,x |(T) = 9 Ideal low-pass filter: ! ⌘ = ⇤p trigonometric polynomial. m 2 @◆||·||1 61 (⌘) , ⌘ 2 @|m|(T) ⌘ = ⇤ m y = ⇤p = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i D (y) = sup || ⇤p||1 61 hp, yi 2 ||p||2
43. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) ⇢ ⌘ ; 8 t, | ⌘ ( t )| 6 1 8 i, ⌘ ( xi) = sign( ai) @|m a,x |(T) = 9 ! Interpolates spikes location and sign. Ideal low-pass filter: ! ⌘ = ⇤p trigonometric polynomial. m 2 @◆||·||1 61 (⌘) , ⌘ 2 @|m|(T) ⌘ = ⇤ m y = ⇤p = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i D (y) = sup || ⇤p||1 61 hp, yi 2 ||p||2
44. From Primal to Dual min m |m|(T) + 1 2 || m y||2 P (y) ⇢ ⌘ ; 8 t, | ⌘ ( t )| 6 1 8 i, ⌘ ( xi) = sign( ai) @|m a,x |(T) = 9 ! Interpolates spikes location and sign. ! |⌘ ( t ) |2 = 1: polynomial equation of supp( m ). Ideal low-pass filter: ! ⌘ = ⇤p trigonometric polynomial. m 2 @◆||·||1 61 (⌘) , ⌘ 2 @|m|(T) ⌘ = ⇤ m y = ⇤p = min m h sup ||⌘||1 61 h⌘, mi + 1 2 || m y||2 i = sup ||⌘||1 61 h min m h⌘, mi + 1 2 || m y||2 i D (y) = sup || ⇤p||1 61 hp, yi 2 ||p||2
45. Asymptotic Dual and Certificate min m |m|(T) + 1 2 || m y||2 P (y) D (y) p def. = argmax || ⇤p||1 61 hp, yi 2 ||p||2
46. Asymptotic Dual and Certificate min m |m|(T) + 1 2 || m y||2 P (y) D (y) P0(y) ! 0+ m0 2 argmin m=y |m|(T) p def. = argmax || ⇤p||1 61 hp, yi 2 ||p||2
47. Asymptotic Dual and Certificate min m |m|(T) + 1 2 || m y||2 P (y) D (y) P0(y) ! 0+ m0 2 argmin m=y |m|(T) D0(y) D0( y ) = argmax || ⇤p||1 61 hp, yi p def. = argmax || ⇤p||1 61 hp, yi 2 ||p||2
48. Asymptotic Dual and Certificate min m |m|(T) + 1 2 || m y||2 P (y) D (y) P0(y) ! 0+ m0 2 argmin m=y |m|(T) p0 def. = argmax p2D0(y) 1 2 ||p||2 ! 0+ D0(y) D0( y ) = argmax || ⇤p||1 61 hp, yi p def. = argmax || ⇤p||1 61 hp, yi 2 ||p||2
49. Asymptotic Dual and Certificate min m |m|(T) + 1 2 || m y||2 P (y) D (y) P0(y) ! 0+ m0 2 argmin m=y |m|(T) p0 def. = argmax p2D0(y) 1 2 ||p||2 ! 0+ D0(y) D0( y ) = argmax || ⇤p||1 61 hp, yi p def. = argmax || ⇤p||1 61 hp, yi 2 ||p||2 D0(y) = {p ; ⇤p 2 @|m0 |(T)} Lemma:
50. Asymptotic Dual and Certificate min m |m|(T) + 1 2 || m y||2 P (y) D (y) P0(y) ! 0+ m0 2 argmin m=y |m|(T) p0 def. = argmax p2D0(y) 1 2 ||p||2 ! 0+ D0(y) D0( y ) = argmax || ⇤p||1 61 hp, yi = 1/fc = 0.6/fc ⌘0 ⌘0 p def. = argmax || ⇤p||1 61 hp, yi 2 ||p||2 D0(y) = {p ; ⇤p 2 @|m0 |(T)} Lemma: Definition: for any m0 solution of P0( y ), ⌘0 = ⇤p0 = argmin ⌘= ⇤p2@|m0 |(T) ||p||
51. −1 1 η 0 η V Vanishing Derivative Pre-certificate 9⌘0 () m0 solves P0( m0) Input measure: m0 = m a,x . ⌘0 def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , || ⌘ ||1 6 1 .
52. −1 1 η 0 η V ⌘0 = ⌘V Vanishing Derivative Pre-certificate 9⌘0 () m0 solves P0( m0) Input measure: m0 = m a,x . ⌘0 def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , || ⌘ ||1 6 1 . ⌘V def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , 8 i, ⌘ 0( xi) = 0 .
53. −1 1 η 0 η V ⌘0 = ⌘V Vanishing Derivative Pre-certificate 9⌘0 () m0 solves P0( m0) Input measure: m0 = m a,x . ⌘0 def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , || ⌘ ||1 6 1 . ⌘V def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , 8 i, ⌘ 0( xi) = 0 . where Ax ( b ) = P i b 1 i ' ( xi, ·) + b 2 i ' 0( xi, ·) Proposition: ⌘ V = ⇤A+ x (sign(a); 0)
54. −1 1 η 0 η V ⌘0 = ⌘V Vanishing Derivative Pre-certificate 9⌘0 () m0 solves P0( m0) Input measure: m0 = m a,x . ⌘0 def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , || ⌘ ||1 6 1 . ⌘V def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , 8 i, ⌘ 0( xi) = 0 . () Non-degenerate certificate: ⌘ 2 ND(m a,x ) : 8 t / 2 { x1, . . . , xN } , | ⌘ ( t )| < 1 and 8 i, ⌘ 00( xi) 6= 0 where Ax ( b ) = P i b 1 i ' ( xi, ·) + b 2 i ' 0( xi, ·) Proposition: ⌘ V = ⇤A+ x (sign(a); 0)
55. −1 1 η 0 η V ⌘0 6= ⌘V −1 1 η 0 η V ⌘0 = ⌘V Vanishing Derivative Pre-certificate 9⌘0 () m0 solves P0( m0) Theorem: ⌘V 2 ND(m0) =) ⌘V = ⌘0 Input measure: m0 = m a,x . ⌘0 def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , || ⌘ ||1 6 1 . ⌘V def. = argmin ⌘= ⇤p || p || s.t. ⇢ 8 i, ⌘ ( xi) = sign( ai) , 8 i, ⌘ 0( xi) = 0 . () Non-degenerate certificate: ⌘ 2 ND(m a,x ) : 8 t / 2 { x1, . . . , xN } , | ⌘ ( t )| < 1 and 8 i, ⌘ 00( xi) 6= 0 where Ax ( b ) = P i b 1 i ' ( xi, ·) + b 2 i ' 0( xi, ·) Proposition: ⌘ V = ⇤A+ x (sign(a); 0)
56. Support Stability Theorem ⌘ = ⇤p !0 ! ⌘0 = ⇤p0 supp(m ) ⇢ {|⌘ | = 1}
57. Support Stability Theorem ⌘ = ⇤p !0 ! ⌘0 = ⇤p0 supp(m ) ⇢ {|⌘ | = 1} If ⌘0 2 ND(m0) then supp(m ) ! supp(m0) ⌘0 ⌘ x1 x2 x ? 2 x ? 1
58. Support Stability Theorem x ? i max Noiseless w = 0. min ⇠ ||w|| x ? i max ||w|| ⌘ = ⇤p !0 ! ⌘0 = ⇤p0 supp(m ) ⇢ {|⌘ | = 1} If ⌘0 2 ND(m0) then supp(m ) ! supp(m0) ⌘0 ⌘ x1 x2 x ? 2 x ? 1 Theorem: the solution of P ( y ) for y = ( m0) + w is for ( ||w||/ , ) = O (1), [Duval, Peyr´ e 2014] If ⌘ V 2 ND( m0) for m0 = m a,x, then m = P N i =1 a? i x ? i where ||(x, a) (x? , a?)|| = O(||w||)
59. When is Non-degenerate ? ⌘V ' ' ⌘V ⌘V ⌘V ⌘V ⌘V ⌘V Input measure: m0 = m a, x , ! 0
60. When is Non-degenerate ? ⌘V ' ' ⌘V ⌘V ⌘V ⌘V ⌘V ⌘V Input measure: Theorem: [Tang, Recht, 2013] Valid for: ' ( x ) = e x 2 / 2 ' ( x ) = (1 + ( x/ )2) 1 . . . 9C, ( > C ) = ) ( ⌘V is non degenerate) m0 = m a, x , ! 0
61. Overview • Sparse Spikes Super-resolution • Robust Support Recovery • Asymptotic Positive Measure Recovery
62. Recovery of Positive Measures m = ( R e 2i⇡ktdm(t))fc k= fc Theorem: let and [de Castro et al. 2011] ! m0 is recovered when there is no noise. ⌘S( t ) = 1 ⇢ QN i=1 sin( ⇡ ( t xi))2 for N 6 fc and ⇢ small enough, ⌘S 2 ¯ D ( m0). -1 1 -1 1 -1 1 -1 1 ⌘S ⌘S ⌘S ⌘S Input measure: m0 = m a,x where a 2 RN + .
63. Recovery of Positive Measures m = ( R e 2i⇡ktdm(t))fc k= fc Theorem: let and [de Castro et al. 2011] ! m0 is recovered when there is no noise. ⌘S( t ) = 1 ⇢ QN i=1 sin( ⇡ ( t xi))2 for N 6 fc and ⇢ small enough, ⌘S 2 ¯ D ( m0). -1 1 -1 1 -1 1 -1 1 ⌘S ⌘S ⌘S ⌘S [Morgenshtern, Cand` es, 2015] discrete `1 robustness. [Demanet, Nguyen, 2015] discrete `0 robustness. Input measure: m0 = m a,x where a 2 RN + . ! behavior as 8 i, xi ! 0 ?
64. Recovery of Positive Measures m = ( R e 2i⇡ktdm(t))fc k= fc Theorem: let and [de Castro et al. 2011] ! m0 is recovered when there is no noise. ! noise robustness of support recovery ? ⌘S( t ) = 1 ⇢ QN i=1 sin( ⇡ ( t xi))2 for N 6 fc and ⇢ small enough, ⌘S 2 ¯ D ( m0). -1 1 -1 1 -1 1 -1 1 ⌘S ⌘S ⌘S ⌘S [Morgenshtern, Cand` es, 2015] discrete `1 robustness. [Demanet, Nguyen, 2015] discrete `0 robustness. Input measure: m0 = m a,x where a 2 RN + . ! behavior as 8 i, xi ! 0 ?
65. Comparison of Certificates -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 ⌘S ⌘V
66. Asymptotic of Vanishing Certificate 1 ⌘V Vanishing Derivative pre-certificate: ⌘V def. = argmin ⌘= ⇤p ||p|| m0 = m a, x where ! 0 s.t. 8 i, ⇢ ⌘ ( xi) = 1 , ⌘ 0( xi) = 0 .
67. Asymptotic of Vanishing Certificate 1 ⌘V 1 1 1 ⌘V ⌘V ⌘W Vanishing Derivative pre-certificate: ⌘V def. = argmin ⌘= ⇤p ||p|| m0 = m a, x where ! 0 s.t. 8 i, ⇢ ⌘ ( xi) = 1 , ⌘ 0( xi) = 0 .
68. Asymptotic of Vanishing Certificate 1 ⌘V 1 1 1 ⌘V ⌘V ⌘W s.t. ⇢ ⌘(0) = 1, ⌘0(0) = . . . = ⌘(2N 1)(0) = 0. Asymptotic pre-certificate: ⌘W def. = argmin ⌘= ⇤p ||p|| Vanishing Derivative pre-certificate: ⌘V def. = argmin ⌘= ⇤p ||p|| ! 0 m0 = m a, x where ! 0 s.t. 8 i, ⇢ ⌘ ( xi) = 1 , ⌘ 0( xi) = 0 .
69. Asymptotic Certificate 1 1 1 1 ⌘V = ⌘W ⌘W ⌘W ⌘W N = 1 N = 2 N = 3 N = 4 (2N 1) -Non degenerate: () ⌘W (2N)(0) 6= 0 ⇢ 8 t 6= 0, |⌘W (t)| < 1 ⌘W 2 NDN
70. Asymptotic Certificate 1 1 1 1 ⌘V = ⌘W ⌘W ⌘W ⌘W N = 1 N = 2 N = 3 N = 4 (2N 1) -Non degenerate: () ⌘W (2N)(0) 6= 0 ⇢ 8 t 6= 0, |⌘W (t)| < 1 ⌘W 2 NDN Lemma: ! ⌘W govern stability as ! 0. If ⌘W 2 NDN , 9 0 > 0, 8 < 0, ⌘ V 2 ND(m x,a )
71. Asymptotic Robustness Theorem: the solution of P ( y ) for y = ( m0) + w is for w , w 2N 1 , 2N 1 = O (1) P N i =1 a? i x ? i ||w || N = 2 ||w || N = 1 If ⌘ W 2 ND N , letting m0 = m a, x , then where ||(x, a) (x?, a?)|| = O ✓ ||w|| + 2N 1 ◆ [Denoyelle, D., P. 2015]
72. Asymptotic Robustness Theorem: the solution of P ( y ) for y = ( m0) + w is x ? i 0 ↵ < 2N 1 Noise: w = w0. for w , w 2N 1 , 2N 1 = O (1) P N i =1 a? i x ? i Regularization: = 0 ↵ 0 max ↵ = 2N 1 x ? i x0 ||w || N = 2 ||w || N = 1 If ⌘ W 2 ND N , letting m0 = m a, x , then where ||(x, a) (x?, a?)|| = O ✓ ||w|| + 2N 1 ◆ y = m a, x + w [Denoyelle, D., P. 2015]
73. When is Non-degenerate ? ⌘W Proposition: one has ⌘(2N) W (0) < 0. ! “locally” non-degenerate.
74. When is Non-degenerate ? ⌘W Proposition: one has ⌘(2N) W (0) < 0. ! “locally” non-degenerate. ' ˆ ' ⌘W ⌘W ⌘W N = 2 N = 3 N = 4
75. Gaussian Deconvolution Gaussian convolution: ' ( x, t ) = e | x t |2 2 2 Proposition: ⌘W ( x ) = e x 2 4 2 N 1 X k=0 ( x/ 2 )2k k ! In particular, ⌘W is non-degenerate. ( m ) def. = Z ' ( x, ·)d m ( x ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.2 0 0.2 0.4 0.6 0.8 1 1 1 1 1 ! Gaussian deconvolution is support-stable. N = 1 N = 2 N = 3 N = 4 '(0, ·) ⌘W ⌘W ⌘W ⌘W
76. Laplace Transform Inversion ( m ) def. = Z ' ( x, ·)d m ( x ) ' ( x, t ) = e xt Laplace transform: ' ( x, ·) x = 2 x = 20 t [with E. Soubies]
77. Laplace Transform Inversion ( m ) def. = Z ' ( x, ·)d m ( x ) ' ( x, t ) = e xt Laplace transform: ' ( x, ·) x = 2 x = 20 t (m1) (m2) x m1 t m2 x [with E. Soubies]
78. Laplace Transform Inversion ( m ) def. = Z ' ( x, ·)d m ( x ) ' ( x, t ) = e xt Laplace transform: ' ( x, ·) x = 2 x = 20 Total internal reflection fluorescence microscopy (TIRFM) t (m1) (m2) x m1 t m2 x [with E. Soubies] [Boulanger et al. 2014] varying the azimuth φ during the exposure time and can be modeled by the following expression: gðθÞ = Z2π 0 Z∞ 0 Z∞ −∞ Iðz; α; φÞρ  θ − α Ω=cos θ  f À z Á dαdzdφ; where fðzÞ is the density of fluorophores in the medium con- volved by the emission point spread function and ρð · Þ represents the laser beam profile of divergence Ω. The function Iðz; α; φÞ slope of the glass slide recovered (Fig. 2D), the latter falling within the confidence interval deducted from the accuracy of the mea- surement of the different characteristic dimensions of the sample. Finally, from the dispersion of the estimated depth around the average slope (Fig. 2D), we can conclude that the localization precision obtained with this approach is higher than the corre- sponding precision given by estimating the location of the beads in the WF image stack as already mentioned (17). Estimating the 3D density of fluorophores convolved by the emission point spread function then would simply boil down to inverting the linear system. Some care has to be taken when inverting such system, as the inverse problem is at best badly con- ditioned. Nevertheless, constraints can be imposed to the solution such as positivity, and, in the case of time-lapse acquisitions, a multiframe regularization can be used in addition to the spatial and temporal regularization smoothness to solve the reconstruction problem. Moreover, to be effective, such a positivity constraint requires a correct knowledge of the background level. As a conse- quence, for each multiangle image stack, a background image is obtained by driving the beam out of the objective. Given that several convex constraints have to be satisfied at the same time, we propose to rely on a flavor of the PPXA algorithm (26) to estimate the tridimensional density of fluorophores (Fig. S4). More detailed information on how noise, object depth, and the required number of angles can be taken into account is discussed in SI Imaging Model and Reconstruction and Fig. S5. Finally, to take into account the variations of the medium index, we select an effective index within a predefined range by minimizing the reconstruction error at each pixel under a spatial smoothness constraint (Fig. S6). It is worth noting that the computation time for the reconstruction on 10 planes from a stack 512 × 512 images corresponding to 21 in- cidence angles ranges from 1 to 5 min depending on the number of iterations. Imaging in Vitro and in Vivo Actin Assembly. The proposed multi- angle TIRF image reconstruction approach was then tested on complex samples such as actin network architectures for which spatial resolution and dynamics remain an issue. We first chal- lenged the spatial organization of actin nucleation geometry using an in vitro assay based on micropatterning method (27). A B C D Fig. 2. Experimental validation of the multiangle TIRF model. (A) Schema of the system designed to create a slope of fluorescent beads. (B) Overlay of the maximum intensity projection of image stack acquired with WF and TIRF illumination. (Scale bar: 5 μm.) The evolution of the intensity versus the illumination angle θ of two selected beads are plotted in C with the corre- sponding fitting theoretical model (continuous line) for their estimated depth (respectively 10 and 89 nm). (D) Depth of all of the beads estimated by fitting the theoretical TIRF model (in red) and the depth of the same beads estimated by fitting a Gaussian model in the WF image stack (in green). BIOPHYSICS AND COMPUTATIONAL BIOLOGY slope of the glass slide recovered (Fig. 2D), the confidence interval deducted from the surement of the different characteristic dim Finally, from the dispersion of the estima average slope (Fig. 2D), we can conclude precision obtained with this approach is h sponding precision given by estimating the l the WF image stack as already mentioned Estimating the 3D density of fluoroph emission point spread function then woul inverting the linear system. Some care h inverting such system, as the inverse proble ditioned. Nevertheless, constraints can be i such as positivity, and, in the case of tim multiframe regularization can be used in add temporal regularization smoothness to so problem. Moreover, to be effective, such requires a correct knowledge of the backgro quence, for each multiangle image stack, obtained by driving the beam out of the several convex constraints have to be satisfie propose to rely on a flavor of the PPXA alg A B C D Fig. 2. Experimental validation of the multiangle TIRF model. (A) Schema of the system designed to create a slope of fluorescent beads. (B) Overlay of the maximum intensity projection of image stack acquired with WF and ✓(t) y(t) light depth x cell y(t) = m(t) ✓(t) ! multiple angles ✓(t).
79. Laplace Transform Inversion ( m ) def. = Z ' ( x, ·)d m ( x ) ' ( x, t ) = e xt Laplace transform: ' ( x, ·) x = 2 x = 20 Total internal reflection fluorescence microscopy (TIRFM) t (m1) (m2) x m1 t m2 x [with E. Soubies] [Boulanger et al. 2014] varying the azimuth φ during the exposure time and can be modeled by the following expression: gðθÞ = Z2π 0 Z∞ 0 Z∞ −∞ Iðz; α; φÞρ  θ − α Ω=cos θ  f À z Á dαdzdφ; where fðzÞ is the density of fluorophores in the medium con- volved by the emission point spread function and ρð · Þ represents the laser beam profile of divergence Ω. The function Iðz; α; φÞ slope of the glass slide recovered (Fig. 2D), the latter falling within the confidence interval deducted from the accuracy of the mea- surement of the different characteristic dimensions of the sample. Finally, from the dispersion of the estimated depth around the average slope (Fig. 2D), we can conclude that the localization precision obtained with this approach is higher than the corre- sponding precision given by estimating the location of the beads in the WF image stack as already mentioned (17). Estimating the 3D density of fluorophores convolved by the emission point spread function then would simply boil down to inverting the linear system. Some care has to be taken when inverting such system, as the inverse problem is at best badly con- ditioned. Nevertheless, constraints can be imposed to the solution such as positivity, and, in the case of time-lapse acquisitions, a multiframe regularization can be used in addition to the spatial and temporal regularization smoothness to solve the reconstruction problem. Moreover, to be effective, such a positivity constraint requires a correct knowledge of the background level. As a conse- quence, for each multiangle image stack, a background image is obtained by driving the beam out of the objective. Given that several convex constraints have to be satisfied at the same time, we propose to rely on a flavor of the PPXA algorithm (26) to estimate the tridimensional density of fluorophores (Fig. S4). More detailed information on how noise, object depth, and the required number of angles can be taken into account is discussed in SI Imaging Model and Reconstruction and Fig. S5. Finally, to take into account the variations of the medium index, we select an effective index within a predefined range by minimizing the reconstruction error at each pixel under a spatial smoothness constraint (Fig. S6). It is worth noting that the computation time for the reconstruction on 10 planes from a stack 512 × 512 images corresponding to 21 in- cidence angles ranges from 1 to 5 min depending on the number of iterations. Imaging in Vitro and in Vivo Actin Assembly. The proposed multi- angle TIRF image reconstruction approach was then tested on complex samples such as actin network architectures for which spatial resolution and dynamics remain an issue. We first chal- lenged the spatial organization of actin nucleation geometry using an in vitro assay based on micropatterning method (27). A B C D Fig. 2. Experimental validation of the multiangle TIRF model. (A) Schema of the system designed to create a slope of fluorescent beads. (B) Overlay of the maximum intensity projection of image stack acquired with WF and TIRF illumination. (Scale bar: 5 μm.) The evolution of the intensity versus the illumination angle θ of two selected beads are plotted in C with the corre- sponding fitting theoretical model (continuous line) for their estimated depth (respectively 10 and 89 nm). (D) Depth of all of the beads estimated by fitting the theoretical TIRF model (in red) and the depth of the same beads estimated by fitting a Gaussian model in the WF image stack (in green). BIOPHYSICS AND COMPUTATIONAL BIOLOGY slope of the glass slide recovered (Fig. 2D), the confidence interval deducted from the surement of the different characteristic dim Finally, from the dispersion of the estima average slope (Fig. 2D), we can conclude precision obtained with this approach is h sponding precision given by estimating the l the WF image stack as already mentioned Estimating the 3D density of fluoroph emission point spread function then woul inverting the linear system. Some care h inverting such system, as the inverse proble ditioned. Nevertheless, constraints can be i such as positivity, and, in the case of tim multiframe regularization can be used in add temporal regularization smoothness to so problem. Moreover, to be effective, such requires a correct knowledge of the backgro quence, for each multiangle image stack, obtained by driving the beam out of the several convex constraints have to be satisfie propose to rely on a flavor of the PPXA alg A B C D Fig. 2. Experimental validation of the multiangle TIRF model. (A) Schema of the system designed to create a slope of fluorescent beads. (B) Overlay of the maximum intensity projection of image stack acquired with WF and ✓(t) y(t) light depth x cell y(t) = m(t) ✓(t) ! multiple angles ✓(t). N = 1 N = 2 N = 3 ⌘W ⌘W ⌘W ¯ x = 2 ¯ x = 20 Non-translation-invariant operator ¯ x x1 x2 ! ⌘W depends on ¯ x!
80. Laplace Transform Inversion ( m ) def. = Z ' ( x, ·)d m ( x ) ' ( x, t ) = e xt Laplace transform: ' ( x, ·) x = 2 x = 20 Total internal reflection fluorescence microscopy (TIRFM) t (m1) (m2) x m1 t m2 x [with E. Soubies] [Boulanger et al. 2014] varying the azimuth φ during the exposure time and can be modeled by the following expression: gðθÞ = Z2π 0 Z∞ 0 Z∞ −∞ Iðz; α; φÞρ  θ − α Ω=cos θ  f À z Á dαdzdφ; where fðzÞ is the density of fluorophores in the medium con- volved by the emission point spread function and ρð · Þ represents the laser beam profile of divergence Ω. The function Iðz; α; φÞ slope of the glass slide recovered (Fig. 2D), the latter falling within the confidence interval deducted from the accuracy of the mea- surement of the different characteristic dimensions of the sample. Finally, from the dispersion of the estimated depth around the average slope (Fig. 2D), we can conclude that the localization precision obtained with this approach is higher than the corre- sponding precision given by estimating the location of the beads in the WF image stack as already mentioned (17). Estimating the 3D density of fluorophores convolved by the emission point spread function then would simply boil down to inverting the linear system. Some care has to be taken when inverting such system, as the inverse problem is at best badly con- ditioned. Nevertheless, constraints can be imposed to the solution such as positivity, and, in the case of time-lapse acquisitions, a multiframe regularization can be used in addition to the spatial and temporal regularization smoothness to solve the reconstruction problem. Moreover, to be effective, such a positivity constraint requires a correct knowledge of the background level. As a conse- quence, for each multiangle image stack, a background image is obtained by driving the beam out of the objective. Given that several convex constraints have to be satisfied at the same time, we propose to rely on a flavor of the PPXA algorithm (26) to estimate the tridimensional density of fluorophores (Fig. S4). More detailed information on how noise, object depth, and the required number of angles can be taken into account is discussed in SI Imaging Model and Reconstruction and Fig. S5. Finally, to take into account the variations of the medium index, we select an effective index within a predefined range by minimizing the reconstruction error at each pixel under a spatial smoothness constraint (Fig. S6). It is worth noting that the computation time for the reconstruction on 10 planes from a stack 512 × 512 images corresponding to 21 in- cidence angles ranges from 1 to 5 min depending on the number of iterations. Imaging in Vitro and in Vivo Actin Assembly. The proposed multi- angle TIRF image reconstruction approach was then tested on complex samples such as actin network architectures for which spatial resolution and dynamics remain an issue. We first chal- lenged the spatial organization of actin nucleation geometry using an in vitro assay based on micropatterning method (27). A B C D Fig. 2. Experimental validation of the multiangle TIRF model. (A) Schema of the system designed to create a slope of fluorescent beads. (B) Overlay of the maximum intensity projection of image stack acquired with WF and TIRF illumination. (Scale bar: 5 μm.) The evolution of the intensity versus the illumination angle θ of two selected beads are plotted in C with the corre- sponding fitting theoretical model (continuous line) for their estimated depth (respectively 10 and 89 nm). (D) Depth of all of the beads estimated by fitting the theoretical TIRF model (in red) and the depth of the same beads estimated by fitting a Gaussian model in the WF image stack (in green). BIOPHYSICS AND COMPUTATIONAL BIOLOGY slope of the glass slide recovered (Fig. 2D), the confidence interval deducted from the surement of the different characteristic dim Finally, from the dispersion of the estima average slope (Fig. 2D), we can conclude precision obtained with this approach is h sponding precision given by estimating the l the WF image stack as already mentioned Estimating the 3D density of fluoroph emission point spread function then woul inverting the linear system. Some care h inverting such system, as the inverse proble ditioned. Nevertheless, constraints can be i such as positivity, and, in the case of tim multiframe regularization can be used in add temporal regularization smoothness to so problem. Moreover, to be effective, such requires a correct knowledge of the backgro quence, for each multiangle image stack, obtained by driving the beam out of the several convex constraints have to be satisfie propose to rely on a flavor of the PPXA alg A B C D Fig. 2. Experimental validation of the multiangle TIRF model. (A) Schema of the system designed to create a slope of fluorescent beads. (B) Overlay of the maximum intensity projection of image stack acquired with WF and ✓(t) y(t) light depth x cell y(t) = m(t) ✓(t) ! multiple angles ✓(t). N = 1 N = 2 N = 3 ⌘W ⌘W ⌘W ¯ x = 2 ¯ x = 20 Non-translation-invariant operator ¯ x x1 x2 ! ⌘W depends on ¯ x! Proposition: In particular, ⌘W is non-degenerate. ⌘W ( x ) = 1 ✓ x ¯ x x + ¯ x ◆2N
81. Deconvolution of measures: ! L2 errors are not well-suited. Weak-* convergence. Optimal transport distance. Exact support estimation. ... Conclusion
82. Deconvolution of measures: ! L2 errors are not well-suited. Weak-* convergence. Optimal transport distance. Exact support estimation. ... Conclusion Low-noise behavior: ! dictated by ⌘0. ! checkable via ⌘V . ! asymptotic via ⌘W .
83. Lasso on discrete grids: Deconvolution of measures: ! L2 errors are not well-suited. Weak-* convergence. Optimal transport distance. Exact support estimation. ... similar ⌘0-analysis applies. ! Relate discrete and continuous recoveries. Conclusion Low-noise behavior: ! dictated by ⌘0. ! checkable via ⌘V . ! asymptotic via ⌘W .
84. Lasso on discrete grids: Deconvolution of measures: ! L2 errors are not well-suited. Weak-* convergence. Optimal transport distance. Exact support estimation. ... similar ⌘0-analysis applies. ! Relate discrete and continuous recoveries. Open problem: other regularizations (e.g. piecewise constant) ? Conclusion Low-noise behavior: ! dictated by ⌘0. ! checkable via ⌘V . ! asymptotic via ⌘W . see [Chambolle, Duval, Peyr´ e, Poon 2016] for TV denoising.