Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Claire Boyer

C3bc10b8a72ed3c3bfd843793b8a9868?s=47 S³ Seminar
January 20, 2017

Claire Boyer

(LSTA, UPMC)

https://s3-seminar.github.io/seminars/claire-boyer

Title — Adapting to unknown noise level in super-resolution

Abstract — We study sparse spikes deconvolution over the space of complex-valued measures when the input measure is a finite sum of Dirac masses. We introduce a new procedure to handle the spike deconvolution when the noise level is unknown. Prediction and localization results will be presented for this approach. An insight on the probabilistic tools used in the proofs could be briefly given as well.

C3bc10b8a72ed3c3bfd843793b8a9868?s=128

S³ Seminar

January 20, 2017
Tweet

Transcript

  1. Spike deconvolution with unknown noise level Claire Boyer, Yohann de

    Castro, Joseph Salmon, UPMC, LSTA, France January 20, 2017
  2. 1 Introduction 2 The proposed approach Results Some words about

    the proof 3 Numerical aspects 4 Numerical experiments
  3. Introduction Motivations Point sources in applications Astronomy Microscopy (fluorescent molecules)

    Spectroscopy 2 / 37 Spike deconvolution
  4. Introduction Motivations Point sources in applications Astronomy Microscopy (fluorescent molecules)

    Spectroscopy Off-the-grid Point sources do not live on a cartesian grid. 2 / 37 Spike deconvolution
  5. Introduction What is spike deconvolution? Setting Aim: reconstruct a discrete

    measure µ0 = s i=1 a0 i δt0 i supp(µ0) = {t0 1 , . . . , t0 s } ⊂ T (a0 i )1≤i≤s ∈ Cs 3 / 37 Spike deconvolution
  6. Introduction What is spike deconvolution? Setting Aim: reconstruct a discrete

    measure µ0 = s i=1 a0 i δt0 i supp(µ0) = {t0 1 , . . . , t0 s } ⊂ T (a0 i )1≤i≤s ∈ Cs From linear measurements y = Φµ0 + ε Φ is a filter ε is a complex Gaussian vector, ε = ε(1) + iε(2) with ε(i) ∼ N (0, σ2 0 Id) 3 / 37 Spike deconvolution
  7. Introduction What is spike deconvolution? Setting Aim: reconstruct a discrete

    measure µ0 = s i=1 a0 i δt0 i supp(µ0) = {t0 1 , . . . , t0 s } ⊂ T (a0 i )1≤i≤s ∈ Cs From linear measurements y = Φµ0 + ε Φ is a filter ε is a complex Gaussian vector, ε = ε(1) + iε(2) with ε(i) ∼ N (0, σ2 0 Id) In this work Low-pass filter : Φ = Fn Fn(µ) := (ck (µ))|k|≤fc , ck (µ) := T exp(−2πıkt)µ(dt) 3 / 37 Spike deconvolution
  8. Introduction How to reconstruct µ0? 1-deconvolution: Beurling lasso (Blasso) min

    µ 1 2 y − Fn (µ) 2 2 + λ µ TV µ TV = sup | j µ(Bj )| over all finite partitions (Bj ) of [0, 1] µ TV = s i=1 |ai | when the measure is discrete 4 / 37 Spike deconvolution
  9. Introduction Beurling minimal extrapolation Beurling minimal extrapolation (1938) µ ∈

    argmin µ µ TV s.t. T Φdµ = y with Φ = (ϕ1, . . . , ϕn ) Infinite-dimensional variable µ Finitely many constraints 5 / 37 Spike deconvolution
  10. Introduction Beurling minimal extrapolation Beurling minimal extrapolation (1938) µ ∈

    argmin µ µ TV s.t. T Φdµ = y with Φ = (ϕ1, . . . , ϕn ) Infinite-dimensional variable µ Finitely many constraints Fenchel dual program ˆ u ∈ arg max y, u s.t. n k=0 uk ϕk ∞ ≤ 1 Finite-dimensional variable Infinitely many constraints Strong duality 5 / 37 Spike deconvolution
  11. Introduction How to reconstruct µ0? 1-deconvolution: Beurling lasso (Blasso) min

    µ 1 2 y − Fn (µ) 2 2 + λ µ TV µ TV = sup | j µ(Bj )| over all finite partitions (Bj ) of [0, 1] µ TV = s i=1 |ai | when the measure is discrete Algorithms proximal-based [Bredies, Pikkarainen 2012] root-finding [Candés, Fernandez-Granda 2012] 6 / 37 Spike deconvolution
  12. Introduction How to reconstruct µ0? 1-deconvolution: Beurling lasso (Blasso) min

    µ 1 2 y − Fn (µ) 2 2 + λ µ TV µ TV = sup | j µ(Bj )| over all finite partitions (Bj ) of [0, 1] µ TV = s i=1 |ai | when the measure is discrete Algorithms proximal-based [Bredies, Pikkarainen 2012] root-finding [Candés, Fernandez-Granda 2012] 6 / 37 Spike deconvolution
  13. Introduction How to reconstruct µ0? 1-deconvolution: Beurling lasso (Blasso) min

    µ 1 2 y − Fn (µ) 2 2 + λ µ TV µ TV = sup | j µ(Bj )| over all finite partitions (Bj ) of [0, 1] µ TV = s i=1 |ai | when the measure is discrete Algorithms proximal-based [Bredies, Pikkarainen 2012] root-finding [Candés, Fernandez-Granda 2012] Others based on Prony’s methods : MUSIC, ESPRIT, FRI 6 / 37 Spike deconvolution
  14. Introduction How to reconstruct µ0? 1-deconvolution: Beurling lasso (Blasso) min

    µ 1 2 y − Fn (µ) 2 2 + λ µ TV µ TV = sup | j µ(Bj )| over all finite partitions (Bj ) of [0, 1] µ TV = s i=1 |ai | when the measure is discrete Algorithms proximal-based [Bredies, Pikkarainen 2012] root-finding [Candés, Fernandez-Granda 2012] Others based on Prony’s methods : MUSIC, ESPRIT, FRI Link with 1-regularization ? 6 / 37 Spike deconvolution
  15. Introduction Link to compressed sensing Spike deconvolution Compressed sensing Setting

    ∞-dim finite-dim (or easy ∞-dim) Object of interest a discrete measure a sparse signal (mainly) µ = s i=1 ai δti x ∈ Cn 1illustrations from Carlos Fernandez-Granda 7 / 37 Spike deconvolution
  16. Introduction Link to compressed sensing Spike deconvolution Compressed sensing Setting

    ∞-dim finite-dim (or easy ∞-dim) Object of interest a discrete measure a sparse signal (mainly) µ = s i=1 ai δti x ∈ Cn Measurements low-pass filter random 1 spectrum extrapolation spectrum interpolation 1illustrations from Carlos Fernandez-Granda 7 / 37 Spike deconvolution
  17. Introduction Link to compressed sensing Spike deconvolution Compressed sensing Setting

    ∞-dim finite-dim (or easy ∞-dim) Object of interest a discrete measure a sparse signal (mainly) µ = s i=1 ai δti x ∈ Cn Measurements low-pass filter random 1 spectrum extrapolation spectrum interpolation Resolution minµ 1 2 y − Fn (µ) 2 2 + λ µ TV minx 1 2 y − Frd (x) 2 2 + λ x 1 Beurling Lasso estimator Lasso Key feature Minimum separation Measurement incoherence + degree of sparsity 1illustrations from Carlos Fernandez-Granda 7 / 37 Spike deconvolution
  18. Introduction Link to compressed sensing Spike deconvolution Compressed sensing Setting

    ∞-dim finite-dim (or easy ∞-dim) Object of interest a discrete measure a sparse signal (mainly) µ = s i=1 ai δti x ∈ Cn Measurements low-pass filter random 1 spectrum extrapolation spectrum interpolation Resolution minµ 1 2 y − Fn (µ) 2 2 + λ µ TV minx 1 2 y − Frd (x) 2 2 + λ x 1 Beurling Lasso estimator Lasso Key feature Minimum separation Measurement incoherence + degree of sparsity Algorithm Root-finding (...) Easy convex programming 1illustrations from Carlos Fernandez-Granda 7 / 37 Spike deconvolution
  19. Introduction Contributions Handling unknown noise level 8 / 37 Spike

    deconvolution
  20. Introduction Contributions Handling unknown noise level Assessing the noise level

    using the Rice method for a non-Gaussian process 8 / 37 Spike deconvolution
  21. Introduction Contributions Handling unknown noise level Assessing the noise level

    using the Rice method for a non-Gaussian process Prediction & strong localization accuracy 8 / 37 Spike deconvolution
  22. Introduction Contributions Handling unknown noise level Assessing the noise level

    using the Rice method for a non-Gaussian process Prediction & strong localization accuracy Numerics : the root-finding can always be done (open question of [1, Eq. (4.5)]) E. J. Candès and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. Communications on Pure and Applied Mathematics, 67(6):906–956, 2014. 8 / 37 Spike deconvolution
  23. Introduction Contributions Handling unknown noise level Assessing the noise level

    using the Rice method for a non-Gaussian process Prediction & strong localization accuracy Numerics : the root-finding can always be done (open question of [1, Eq. (4.5)]) Claire Boyer, Yohann de Castro, and Joseph Salmon. Adapting to unknown noise level in sparse deconvolution. Information and Inference, abs/1606.04760, 2016. E. J. Candès and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. Communications on Pure and Applied Mathematics, 67(6):906–956, 2014. 8 / 37 Spike deconvolution
  24. 1 Introduction 2 The proposed approach Results Some words about

    the proof 3 Numerical aspects 4 Numerical experiments
  25. The proposed approach The proposed method Data y = Fnµ0

    + ε y ∈ Cn with n = 2fc + 1 ε is a complex Gaussian vector, ε = ε(1) + iε(2) with ε(j) ∼ N(0, σ2 0 Id) 9 / 37 Spike deconvolution
  26. The proposed approach The proposed method Data y = Fnµ0

    + ε y ∈ Cn with n = 2fc + 1 ε is a complex Gaussian vector, ε = ε(1) + iε(2) with ε(j) ∼ N(0, σ2 0 Id) Handling unknown noise level : CBLasso (ˆ µ, ˆ σ) ∈ argmin (µ,σ)∈E∗×R++ 1 2nσ y − Fn (µ) 2 2 + σ 2 + λ µ TV , Owen 07, Antoniadis 10, Belloni Chernozhukov Wang 11, Sun Zhang 12, van de Geer 15 : square-root lasso and scaled-lasso. 9 / 37 Spike deconvolution
  27. The proposed approach Where is the square-root? Why square-root? The

    CBLasso reads ([2] ∼ scaled lasso) (ˆ µ, ˆ σ) ∈ argmin (µ,σ)∈E∗×R++ 1 2nσ y − Fn (µ) 2 2 + σ 2 + λ µ TV . When the solution is reached for ˆ σ > 0, ([1] ∼ square-root lasso) ˆ σ = y − Fn (ˆ µ) 2/ √ n ˆ µ ∈ argmin µ∈E∗ y − Fn (µ) 2/ √ n + λ µ TV A. Belloni, V. Chernozhukov, and L. Wang. Square-root Lasso: Pivotal recovery of sparse signals via conic programming. Biometrika, 98(4):791–806, 2011. T. Sun and C.-H. Zhang. Scaled sparse linear regression. Biometrika, 99(4):879–898, 2012. 10 / 37 Spike deconvolution
  28. The proposed approach Compatibility limits Sufficient conditions for oracle inequalities

    RIP =⇒ REC =⇒ Compatibility Compatibility condition : C(L, S) > 0 C(L, S) := inf |S| Fn(ν) 2 2 /n s.t. supp(ν) = S, νS TV = 1, νSc TV ≤ L . Consider ν as a difference of two measures: ˆ µ and µ0 11 / 37 Spike deconvolution
  29. The proposed approach Compatibility limits Sufficient conditions for oracle inequalities

    RIP =⇒ REC =⇒ Compatibility Compatibility condition : C(L, S) > 0 C(L, S) := inf |S| Fn(ν) 2 2 /n s.t. supp(ν) = S, νS TV = 1, νSc TV ≤ L . Consider ν as a difference of two measures: ˆ µ and µ0 Compatibility does not hold νε = δ−ε + δε @@@@@ Compatibility =⇒ $$ REC =⇒ ¨ ¨ RIP highly coherent designs: close Dirac masses share almost the same Fourier coefficients Non-uniform approach Measure-dependent reconstruction 11 / 37 Spike deconvolution
  30. The proposed approach – Results Standard assumptions Assumption (Sampling rate

    condition). λSNR ≤ √ 17 − 4 2 0.0616 with SNR := µ0 TV E[ ε 2 2 ]/n = µ0 TV √ 2σ0 . λ ≥ 2 2 log n/n =⇒ n/log n ≥ C SNR2 12 / 37 Spike deconvolution
  31. The proposed approach – Results Standard assumptions Assumption (Sampling rate

    condition). λSNR ≤ √ 17 − 4 2 0.0616 with SNR := µ0 TV E[ ε 2 2 ]/n = µ0 TV √ 2σ0 . λ ≥ 2 2 log n/n =⇒ n/log n ≥ C SNR2 Assumption (No-overfitting). ˆ σ > 0 12 / 37 Spike deconvolution
  32. The proposed approach – Results Standard assumptions Assumption (Sampling rate

    condition). λSNR ≤ √ 17 − 4 2 0.0616 with SNR := µ0 TV E[ ε 2 2 ]/n = µ0 TV √ 2σ0 . λ ≥ 2 2 log n/n =⇒ n/log n ≥ C SNR2 Assumption (No-overfitting). ˆ σ > 0 Assumption (Separation condition). For supp(µ0) = {t0 1 , . . . , t0 s0 }, ∀i = j, d(t0 i , t0 j ) ≥ 1.26 fc 12 / 37 Spike deconvolution
  33. The proposed approach – Results Prediction result Theorem (B., de

    Castro, Salmon, 2016). Let C > 2 √ 2. Set C > 0, that may depend on C. Assume the sampling rate condition, the separation condition. The estimator ˆ µ solution to CBLasso with a choice λ ≥ C log n/n satisfies 1 n Fn (ˆ µ − µ0) 2 2 ≤ C s0 λ2 σ0, with high probability. 13 / 37 Spike deconvolution
  34. The proposed approach – Results Prediction result Theorem (B., de

    Castro, Salmon, 2016). Let C > 2 √ 2. Set C > 0, that may depend on C. Assume the sampling rate condition, the separation condition. The estimator ˆ µ solution to CBLasso with a choice λ ≥ C log n/n satisfies 1 n Fn (ˆ µ − µ0) 2 2 ≤ C s0 λ2 σ2 0 = O s0σ2 0 log n n , with high probability. "fast rate of convergence" of [1] G. Tang, B. N. Bhaskar, and B. Recht. Near minimax line spectral estimation. Information Theory, IEEE Transactions on, 61(1):499–512, 2015. 13 / 37 Spike deconvolution
  35. The proposed approach – Results Localization results Near regions ∀j

    ∈ [s0 ], Nj := t ∈ [0, 1]; d(t, t0 j ) ≤ c1 fc , with 0 < c1 < 1.26/2. Far region F := [0, 1] \ j∈[s0] Nj 14 / 37 Spike deconvolution
  36. The proposed approach – Results Localization results Near regions ∀j

    ∈ [s0 ], Nj := t ∈ [0, 1]; d(t, t0 j ) ≤ c1 fc , with 0 < c1 < 1.26/2. Far region F := [0, 1] \ j∈[s0] Nj Theorem (B., de Castro, Salmon, 2016). Let C > 2 √ 2. Let C > 0 may depend on C. The estimator ˆ µ, solution to CBLasso with a choice λ ≥ C log n/n, satisfies 1 ∀j ∈ [s0], a0 j − {k: ˆ tk ∈Nj } ˆ ak ≤ C σ0s0 log n/n , 14 / 37 Spike deconvolution
  37. The proposed approach – Results Localization results Near regions ∀j

    ∈ [s0 ], Nj := t ∈ [0, 1]; d(t, t0 j ) ≤ c1 fc , with 0 < c1 < 1.26/2. Far region F := [0, 1] \ j∈[s0] Nj Theorem (B., de Castro, Salmon, 2016). Let C > 2 √ 2. Let C > 0 may depend on C. The estimator ˆ µ, solution to CBLasso with a choice λ ≥ C log n/n, satisfies 1 ∀j ∈ [s0], a0 j − {k: ˆ tk ∈Nj } ˆ ak ≤ C σ0s0 log n/n , 2 ∀j ∈ [s0], {k: ˆ tk ∈Nj } |ˆ ak |d2(t0 j , ˆ tk ) ≤ C σ0s0 log n/n/n2 , 14 / 37 Spike deconvolution
  38. The proposed approach – Results Localization results Near regions ∀j

    ∈ [s0 ], Nj := t ∈ [0, 1]; d(t, t0 j ) ≤ c1 fc , with 0 < c1 < 1.26/2. Far region F := [0, 1] \ j∈[s0] Nj Theorem (B., de Castro, Salmon, 2016). Let C > 2 √ 2. Let C > 0 may depend on C. The estimator ˆ µ, solution to CBLasso with a choice λ ≥ C log n/n, satisfies 1 ∀j ∈ [s0], a0 j − {k: ˆ tk ∈Nj } ˆ ak ≤ C σ0s0 log n/n , 2 ∀j ∈ [s0], {k: ˆ tk ∈Nj } |ˆ ak |d2(t0 j , ˆ tk ) ≤ C σ0s0 log n/n/n2 , 3 {k: ˆ tk ∈F} |ˆ ak | ≤ C σ0s0 log n/n , with high probability. 14 / 37 Spike deconvolution
  39. The proposed approach – Results Localization results Corollary (B., de

    Castro, Salmon, 2016). For any t0 j in the support of µ0 such that a0 j > C σ0 s0λ, there exists an element ˆ tk in the support of ˆ µ such that d(t0 j , ˆ tk ) ≤ C σ0 s0λ |a0 j | − C σ0 s0λ 1 n , with high probability. independent of the other spikes magnitude 15 / 37 Spike deconvolution
  40. The proposed approach – Results Noise level estimation Proposition (B.,

    de Castro, Salmon, 2016). Under the sampling rate assumption, if n > −8 log(α), √ nˆ σ ε 2 − 1 ≤ 1 2 , with probability larger than 1 − α 2 √ 2 n + 2 √ 3+3 3 . 16 / 37 Spike deconvolution
  41. The proposed approach – Some words about the proof Sketch

    of proof KKT condition 1 n F∗ n (y − Fn (ˆ µ)) = ˆ σλˆ p ˆ p ∞ ≤ 1 R T ˆ p(t)ˆ µ(dt) = ˆ µ TV (ˆ p subgradient of · TV at ˆ µ) 17 / 37 Spike deconvolution
  42. The proposed approach – Some words about the proof Sketch

    of proof KKT condition 1 n F∗ n (y − Fn (ˆ µ)) = ˆ σλˆ p ˆ p ∞ ≤ 1 R T ˆ p(t)ˆ µ(dt) = ˆ µ TV (ˆ p subgradient of · TV at ˆ µ) 17 / 37 Spike deconvolution
  43. The proposed approach – Some words about the proof Sketch

    of proof KKT condition 1 n F∗ n (y − Fn (ˆ µ)) = ˆ σλˆ p ˆ p ∞ ≤ 1 R T ˆ p(t)ˆ µ(dt) = ˆ µ TV (ˆ p subgradient of · TV at ˆ µ) Dual certificate Construction of q(t) = fc k=−fc ck e2iπkt obeying to q(tj ) = v0 j := a0 j |a0 j | ∀tj ∈ supp µ0 |q(t)| < 1 ∀t / ∈ supp µ0 17 / 37 Spike deconvolution
  44. The proposed approach – Some words about the proof Sketch

    of proof KKT condition 1 n F∗ n (y − Fn (ˆ µ)) = ˆ σλˆ p ˆ p ∞ ≤ 1 R T ˆ p(t)ˆ µ(dt) = ˆ µ TV (ˆ p subgradient of · TV at ˆ µ) Dual certificate Construction of q(t) = fc k=−fc ck e2iπkt obeying to q(tj ) = v0 j := a0 j |a0 j | ∀tj ∈ supp µ0 q (tj ) = 0 |q(t)| < 1 ∀t / ∈ supp µ0 17 / 37 Spike deconvolution
  45. The proposed approach – Some words about the proof Sketch

    of proof Spike localization Proofs are based on the work of [1, 2, 3] amended by Rice method Prediction adapt the proof of [3] to our setting Rice method for a non-Gaussian process J.-M. Azaïs, Y. De Castro, and F. Gamboa. Spike detection from inaccurate samplings. Applied and Computational Harmonic Analysis, 38(2):177–195, 2015. C. Fernandez-Granda. Support detection in super-resolution. arXiv preprint arXiv:1302.3921, 2013. G. Tang, B. N. Bhaskar, and B. Recht. Near minimax line spectral estimation. Information Theory, IEEE Transactions on, 61(1):499–512, 2015. 18 / 37 Spike deconvolution
  46. 1 Introduction 2 The proposed approach Results Some words about

    the proof 3 Numerical aspects 4 Numerical experiments
  47. Numerical aspects Primal-dual Proposition. Denoting Dn = c ∈ Cn

    : F∗ n (c) ∞ ≤ 1, nλ2 c 2 ≤ 1 , the dual formulation of the CBLasso reads ˆ c ∈ arg max c∈Dn λ y, c . (1) Then, we have the link-equation between primal and dual solutions y = nλˆ σˆ c + Fn (ˆ µ) . (2) as well as a link between the coefficient and the polynomial F∗ n (ˆ c) = ˆ p . (3) 19 / 37 Spike deconvolution
  48. Numerical aspects Primal-dual Proposition. Denoting Dn = c ∈ Cn

    : F∗ n (c) ∞ ≤ 1, nλ2 c 2 ≤ 1 , the dual formulation of the CBLasso reads ˆ c ∈ arg max c∈Dn λ y, c . (1) Then, we have the link-equation between primal and dual solutions y = nλˆ σˆ c + Fn (ˆ µ) . (2) as well as a link between the coefficient and the polynomial F∗ n (ˆ c) = ˆ p . (3) If the dual polynomial ˆ p associated is not constant the support of ˆ µ is finite, the support of ˆ µ is included in the set of its derivative roots, i.e. where the polynomial saturates at 1. the proof based on a dual certificate construction is possible. 19 / 37 Spike deconvolution
  49. Numerical aspects Primal-dual : comparison with Blasso Dual of the

    Blasso ˆ c ∈ arg max F∗ n (c) ∞≤1 y, c . Dual of the CBLasso ˆ c ∈ arg max F∗ n (c) ∞≤1 nλ2 c 2≤1 λ y, c . 20 / 37 Spike deconvolution
  50. Numerical aspects Primal-dual : comparison with Blasso Dual of the

    Blasso ˆ c ∈ arg max F∗ n (c) ∞≤1 y, c . Open question When the dual polynomial is non- constant ? Dual of the CBLasso ˆ c ∈ arg max F∗ n (c) ∞≤1 nλ2 c 2≤1 λ y, c . For the CBLasso We answer this. 20 / 37 Spike deconvolution
  51. Numerical aspects Showing that the dual polynomial is non-constant |ˆ

    p|2 is of constant modulus 1 ⇒ ˆ p = vϕk with v ∈ C and ϕk (·) = exp(2πık·) for some −fc ≤ k ≤ fc . if |v| < 1, using Holder’s inequality on R T ˆ p(t)ˆ µ(dt) = ˆ µ TV leads to ˆ µ = 0. if |v| = 1, we also have ˆ c ∈ Dn, so ˆ c 2 ≤ 1/( √ nλ), leading to |v| ≤ 1/( √ nλ). Since λ > 2 √ logn/ √ n ⇒ |v| < 1, which contradicts |v| = 1. One can then conclude that a dual polynomial of constant modulus never occurs in the CBLasso setup 21 / 37 Spike deconvolution
  52. Numerical aspects Discussion on the value of λ 22 /

    37 Spike deconvolution
  53. Numerical aspects SDP formulation of the CBLasso One can represent

    the dual feasible set Dn as an SDP condition. The dual problem can be cast as follows max c∈Cn λ y, c such that ∃Q ∈ Cn×n Q c c∗ 1 0 n−j i=1 Qi,i+j = 1 if j = 0 0 j = 1, . . . , n − 1. 23 / 37 Spike deconvolution
  54. Numerical aspects SDP formulation of the CBLasso One can represent

    the dual feasible set Dn as an SDP condition. The dual problem can be cast as follows max c∈Cn λ y, c such that ∃Q ∈ Cn×n Q c c∗ 1 0 n−j i=1 Qi,i+j = 1 if j = 0 0 j = 1, . . . , n − 1. In λ √ nc λ √ nc∗ 1 0 23 / 37 Spike deconvolution
  55. Numerical aspects Proposed algorithm Algorithm Given the data y ∈

    Cn 1 solve the dual problem to find the coefficients ˆ c of the dual polynomial ˆ p (cvx Matlab toolbox); 24 / 37 Spike deconvolution
  56. Numerical aspects Proposed algorithm Algorithm Given the data y ∈

    Cn 1 solve the dual problem to find the coefficients ˆ c of the dual polynomial ˆ p (cvx Matlab toolbox); 2 identify supp(ˆ µ) = {ˆ t1, . . . , ˆ tˆ s } using the roots of 1 − |ˆ p|2; 24 / 37 Spike deconvolution
  57. Numerical aspects Proposed algorithm Algorithm Given the data y ∈

    Cn 1 solve the dual problem to find the coefficients ˆ c of the dual polynomial ˆ p (cvx Matlab toolbox); 2 identify supp(ˆ µ) = {ˆ t1, . . . , ˆ tˆ s } using the roots of 1 − |ˆ p|2; 3 fix the design matrix X ∈ Cn׈ s , defined by Xk,j = ϕk (ˆ tj ); 24 / 37 Spike deconvolution
  58. Numerical aspects Proposed algorithm Algorithm Given the data y ∈

    Cn 1 solve the dual problem to find the coefficients ˆ c of the dual polynomial ˆ p (cvx Matlab toolbox); 2 identify supp(ˆ µ) = {ˆ t1, . . . , ˆ tˆ s } using the roots of 1 − |ˆ p|2; 3 fix the design matrix X ∈ Cn׈ s , defined by Xk,j = ϕk (ˆ tj ); 4 solve now the following finite-dimensional problem (ˆ a, ˆ σ) ∈ argmin (a,σ)∈Cˆ s ×R++ 1 2nσ y − Xa 2 2 + σ 2 + λ a 1 , 24 / 37 Spike deconvolution
  59. Numerical aspects Proposed algorithm Algorithm Given the data y ∈

    Cn 1 solve the dual problem to find the coefficients ˆ c of the dual polynomial ˆ p (cvx Matlab toolbox); 2 identify supp(ˆ µ) = {ˆ t1, . . . , ˆ tˆ s } using the roots of 1 − |ˆ p|2; 3 fix the design matrix X ∈ Cn׈ s , defined by Xk,j = ϕk (ˆ tj ); 4 solve now the following finite-dimensional problem (ˆ a, ˆ σ) ∈ argmin (a,σ)∈Cˆ s ×R++ 1 2nσ y − Xa 2 2 + σ 2 + λ a 1 , as follows: for an initial value of ˆ σ, until some stopping criterion, alternate the following steps solve the previous problem for a fixed ˆ σ to compute ˆ a, ˆ a = X+y − λˆ σ(X∗X)−1 sign(X∗ ˆ c) update ˆ σ = y − X ˆ a 2/ √ n using the new value of ˆ a, 24 / 37 Spike deconvolution
  60. 1 Introduction 2 The proposed approach Results Some words about

    the proof 3 Numerical aspects 4 Numerical experiments
  61. Numerical experiments Source code Code available at https://github.com/claireBoyer/CBLasso http://www.lsta.upmc.fr/boyer/ 25

    / 37 Spike deconvolution
  62. Numerical experiments Measure reconstruction 0 0.2 0.4 0.6 0.8 1

    −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Original BLasso CBLasso Reconstruction of a discrete measure. The original measure µ0 is composed of 3 spikes (in black). The reconstructed measure ˆ µ using our proposed CBLasso (in blue). In comparison, we plot the reconstructed measure using the BLasso, (in red). 26 / 37 Spike deconvolution
  63. Numerical experiments Noise level estimation ε = ε(1) + iε(2)

    with ε(j) ∼ N(0, σ0 Id) with σ0 = 1/ √ 2. 27 / 37 Spike deconvolution
  64. Numerical experiments Noise level estimation ε = ε(1) + iε(2)

    with ε(j) ∼ N(0, σ0 Id) with σ0 = 1/ √ 2. 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 BLasso CBLasso Boxplot on ˆ σ for 100 CBLasso consistent estimations of √ 2σ0 = 1. We compare our method to ˆ σBLasso = y − Fn (ˆ µBLasso) 2 √ n − ˆ sBLasso proposed in [1] where ˆ µBLasso is the reconstructed measure supported on ˆ sBLasso spikes via BLasso. S. Reid, R. Tibshirani, and J. Friedman. A study of error variance estimation in lasso regression. arXiv preprint arXiv:1311.5274, 2014. 27 / 37 Spike deconvolution
  65. Numerical experiments Bias Noise level estimation √ nˆ σ ε

    2 − 1 ≤ 1 2 , with high probability. Bias ˆ σ √ 2σ0 × E g 2 √ 2n = √ 2σ0 × Γ(n + 1/2) √ nΓ(n) → √ 2σ0 , showing that ˆ σ/ √ 2 is a consistent estimator of σ0 . 28 / 37 Spike deconvolution
  66. Numerical experiments Conclusion The CBLasso new approach to handle unknown

    noise level in spike detection theoretical contributions : prediction for CBLasso localization for CBLasso closing the question of constant polynomial in this setting numerical method is also proposed to tackle the CBLasso 29 / 37 Spike deconvolution
  67. Numerical experiments Conclusion The CBLasso new approach to handle unknown

    noise level in spike detection theoretical contributions : prediction for CBLasso localization for CBLasso closing the question of constant polynomial in this setting numerical method is also proposed to tackle the CBLasso Claire Boyer, Yohann de Castro, and Joseph Salmon. Adapting to unknown noise level in sparse deconvolution. Information and Inference, abs/1606.04760, 2016. 29 / 37 Spike deconvolution
  68. Numerical experiments Thank you for your attention

  69. Discussion on λ Discussion on λ λmax(y) = F∗ n

    (y) ∞/ √ n y 2 (for Blasso F∗ n (y) ∞/n).
  70. Proposition. Defining λmin (y) = 1/( ˆ c(BME) 2 √

    n) and the problem ˆ µ(BME) ∈ argmin Fn(µ)=y µ TV , (4) and its dual formulation ˆ c(BME) ∈ arg max c∈Cn y, c s.t. F∗ n (c) ∞ ≤ 1 . (5) the following statements are equivalent (i) λ ∈]0, λmin(y)], (ii) ˆ c = ˆ c(BME), (iii) ˆ σ = 0 (overfitting). Remark that λmin (y) = 1/( ˆ c(BME) 2 √ n) > 1/ √ n.
  71. Proof tricks Rice method One needs to upper bound the

    probabilities of the following events { F∗ n (ε) ∞ ≤ λ} and F∗ n (ε) ∞ √ n ε 2 ≤ R , where F∗ n (ε)(t) = fc k=−fc εk exp(2πıkt) , for all t ∈ [0, 1]. The first event can be handled with a Rice formula for stationary Gaussian processes as in [1]. Due to the presence of the denominator, the second event cannot be described by a Gaussian process and its analysis is a bit more involved.
  72. Proof tricks Rice method : notation We start with some

    notation. Let z(1) = (z(1) −fc , . . . , z(1) 0 , . . . , z(1) fc ), z(2) = (z(2) −fc , . . . , z(2) 0 , . . . , z(2) fc ) , be i.i.d Nn (0, Idn ) random vectors. Set, for any t ∈ [0, 1], X(t) = z(1) 0 + fc k=1 (z(1) k + z(1) −k ) cos(2πkt) + fc k=1 (z(2) −k − z(2) k ) sin(2πkt) , Y (t) = z(2) 0 + fc k=1 (z(2) k + z(2) −k ) cos(2πkt) + fc k=1 (z(1) k − z(1) −k ) sin(2πkt) , Z(t) = X(t) + ıY (t) . Then, note that σ0 Z ∞ d = F∗ n (ε) ∞ and sup t∈[0,1] |Z(t)| √ n( z(1) 2 2 + z(2) 2 2 )1 2 d = F∗ n (ε) ∞ √ n ε 2 , where σ0 > 0 is the (unknown) standard deviation of the noise ε.
  73. Proof tricks For the Gaussian process Lemma 1. For a

    complex valued centered Gaussian random vector ε as defined in (??), it holds ∀u > 0, P { F∗ n (ε) ∞ > u} ≤ n exp − u2 4nσ2 0 , where σ0 > 0 denotes the noise level. Observe that X(t) and Y (t) are two independent stationary Gaussian processes with the same auto-covariance function Σ given by ∀t ∈ [0, 1], Σ(t) = 1 + 2 fc k=1 cos(2πkt) =: Dfc (t) , where Dfc denotes the Dirichlet kernel. Set σn 2 = var X(t) = Dfc (0) = n . (6)
  74. Proof tricks For the Gaussian process We use the following

    inequalities, for any u > 0, P{ Z ∞ > u} ≤ P{ X ∞ > u/ √ 2}+P{ Y ∞ > u/ √ 2} = 2P{ X ∞ > u/ √ 2} , (7) and (by symmetry of the process X) P{ X ∞ > u/ √ 2} ≤ 2P{ sup t∈[0,1] X(t) > u/ √ 2} . (8) To give bounds to the right hand side of (8), we use the Rice method (see Azais & Wschebor [page 93]) P{ sup t∈[0,1] X(t) > u/ √ 2} = P{∀t ∈ [0, 1]; X(t) > u/ √ 2} + P{U u/ √ 2 > 0} , ≤ 1 √ 2π +∞ u/( √ 2σn) exp − v2 2 d + E(U u/ √ 2 ) , where Uv is the number of up-crossings of the level v by the process X(t) on the interval [0, 1]. By the Rice formula E(U u/ √ 2 ) = 1 2π var X (t)) 1 σn exp − u2 4σn 2 , where var X (t) = −Σ (0) = 2(2π)2 fc k=1 k2 = 4π2 3 fc (fc + 1)n . (9)
  75. Proof tricks For the Gaussian process A Chernoff argument provides

    for any w > 0, +∞ w exp(−v2/2)d/ √ 2π ≤ exp(−w2/2), which yields P sup t∈[0,1] X(t) > u √ 2 ≤ exp − u2 4n + fc (fc + 1) 3 exp − u2 4n ≤ n 2 exp − u2 4n . The result follows with (8).