$30 off During Our Annual Pro Sale. View Details »

Off-the-Grid Sparse Super-Resolution

Gabriel Peyré
February 19, 2017

Off-the-Grid Sparse Super-Resolution

Gabriel Peyré

February 19, 2017
Tweet

More Decks by Gabriel Peyré

Other Decks in Research

Transcript

  1. Gabriel Peyré
    Off-the-Grid
    Sparse Super-resolution
    Joint work with
    Vincent Duval & Quentin Denoyelle
    É C O L E N O R M A L E
    S U P É R I E U R E
    Outline
    • What is the Lasso
    • Lasso with an orthogonal design
    • From projected gradient to proximal grad
    • Optimality conditions and subgradients (L
    • Coordinate descent algorithm
    … with some demos
    www.numerical-tours.com

    View Slide

  2. Sparse Super-resolution
    Neural spikes (1D)
    Astrophysics (2D)
    Seismic imaging (1.5D)
    A
    B
    Fig. 2. Experimental validation of the m
    of the system designed to create a slope
    of the maximum intensity projection of i
    TIRF illumination. (Scale bar: 5 μm.) The ev
    illumination angle θ of two selected bead
    sponding fitting theoretical model (con
    depth (respectively 10 and 89 nm). (D) Dep
    [Boulanger et al. 2014]
    Single-molecule
    fluorescence (3-D)
    Recover pointwise sources
    from noisy low-resolution
    observations.

    View Slide

  3. Sparse Super-resolution
    Neural spikes (1D)
    Astrophysics (2D)
    Seismic imaging (1.5D)
    A
    B
    Fig. 2. Experimental validation of the m
    of the system designed to create a slope
    of the maximum intensity projection of i
    TIRF illumination. (Scale bar: 5 μm.) The ev
    illumination angle θ of two selected bead
    sponding fitting theoretical model (con
    depth (respectively 10 and 89 nm). (D) Dep
    [Boulanger et al. 2014]
    Single-molecule
    fluorescence (3-D)
    Recover pointwise sources
    from noisy low-resolution
    observations.
    Practice: Scalable algorithms?
    Theory: Rayleigh limit?

    View Slide

  4. Overview
    • Sparse Spikes Super-resolution
    • Robust Support Recovery
    • Asymptotic Positive Measure Recovery
    • Off-the-Grid Optimization Algorithms
    • Application: Laplace Inversion for TIRF Imaging

    View Slide

  5. Discrete measure:
    ma,x =
    PN
    i=1
    ai xi
    , a 2 RN , x 2 TN
    Super-resolution of Measures
    ma,x
    m
    Radon measure m on T =

    (R/Z)d
    Rd
    .

    View Slide

  6. Discrete measure:
    ma,x =
    PN
    i=1
    ai xi
    , a 2 RN , x 2 TN
    Super-resolution of Measures
    ma,x
    m
    Radon measure m on T =

    (R/Z)d
    Rd
    .
    y = (m) + w
    Linear measurements: '(x) 2 H
    ' continuous.
    (m) def.
    =
    Z
    T
    '(x)dm(x) 2 H

    View Slide

  7. Discrete measure:
    ma,x =
    PN
    i=1
    ai xi
    , a 2 RN , x 2 TN
    Super-resolution of Measures
    ma,x
    m
    Radon measure m on T =

    (R/Z)d
    Rd
    .
    '(x) = ˜
    '(· x)
    Deconvolution:
    ! Signal-dependent recovery criteria.
    y
    = 2/fc
    y
    = 0.5/fc
    '(0)
    y = (m) + w
    Linear measurements: '(x) 2 H
    ' continuous.
    (m) def.
    =
    Z
    T
    '(x)dm(x) 2 H

    View Slide

  8. Discrete measure:
    ma,x =
    PN
    i=1
    ai xi
    , a 2 RN , x 2 TN
    Super-resolution of Measures
    ma,x
    m
    Radon measure m on T =

    (R/Z)d
    Rd
    .
    '(x) = ˜
    '(· x)
    Deconvolution:
    ! Signal-dependent recovery criteria.
    y
    = 2/fc
    y
    = 0.5/fc
    '(0)
    Fourier: '(x) = (ei`x)fc
    `= fc
    2 C2fc+1
    y = (m) + w
    Linear measurements: '(x) 2 H
    ' continuous.
    (m) def.
    =
    Z
    T
    '(x)dm(x) 2 H

    View Slide

  9. Discrete measure:
    ma,x =
    PN
    i=1
    ai xi
    , a 2 RN , x 2 TN
    Super-resolution of Measures
    ma,x
    m
    Radon measure m on T =

    (R/Z)d
    Rd
    .
    '(x) = ˜
    '(· x)
    Deconvolution:
    ! Signal-dependent recovery criteria.
    y
    = 2/fc
    y
    = 0.5/fc
    '(0)
    Fourier: '(x) = (ei`x)fc
    `= fc
    2 C2fc+1
    Laplace: '(x) = e x·
    2 H
    def.
    = L2(R+)
    y = (m) + w
    Linear measurements: '(x) 2 H
    ' continuous.
    (m) def.
    =
    Z
    T
    '(x)dm(x) 2 H

    View Slide

  10. Grid-free Sparse Recovery
    Grid-free regularization: total variation of measures:
    |m|(T) = sup
    R
    ⌘dm : ⌘ 2 C(T), ||⌘||1
    6 1
    |ma,x
    |(T) = ||a||`1
    ma,x
    |m|(T) =
    R
    |f| = ||f||L1
    dm(x) = f(x)dx

    View Slide

  11. Grid-free Sparse Recovery
    Grid-free regularization: total variation of measures:
    |m|(T) = sup
    R
    ⌘dm : ⌘ 2 C(T), ||⌘||1
    6 1
    |ma,x
    |(T) = ||a||`1
    ma,x
    |m|(T) =
    R
    |f| = ||f||L1
    dm(x) = f(x)dx
    min
    m
    1
    2
    || (m) y||2 + |m|(T) (P (y))
    Sparse recovery:
    (P0(y))
    min
    m
    {|m|(T) ; m = y}
    ! 0+

    View Slide

  12. Grid-free Sparse Recovery
    Grid-free regularization: total variation of measures:
    |m|(T) = sup
    R
    ⌘dm : ⌘ 2 C(T), ||⌘||1
    6 1
    |ma,x
    |(T) = ||a||`1
    ma,x
    |m|(T) =
    R
    |f| = ||f||L1
    dm(x) = f(x)dx
    [Fischer Jerome, 1974]
    If dim(Im( )) < +1, 9(a, x) 2 RN ⇥ TN
    with N 6 dim(Im( ))
    such that ma,x is a solution to P (y).
    Proposition:
    min
    m
    1
    2
    || (m) y||2 + |m|(T) (P (y))
    Sparse recovery:
    (P0(y))
    min
    m
    {|m|(T) ; m = y}
    ! 0+

    View Slide

  13. Grid-free Sparse Recovery
    Grid-free regularization: total variation of measures:
    |m|(T) = sup
    R
    ⌘dm : ⌘ 2 C(T), ||⌘||1
    6 1
    |ma,x
    |(T) = ||a||`1
    ma,x
    |m|(T) =
    R
    |f| = ||f||L1
    dm(x) = f(x)dx
    [Fischer Jerome, 1974]
    If dim(Im( )) < +1, 9(a, x) 2 RN ⇥ TN
    with N 6 dim(Im( ))
    such that ma,x is a solution to P (y).
    Proposition:
    min
    m
    1
    2
    || (m) y||2 + |m|(T) (P (y))
    Sparse recovery:
    (P0(y))
    min
    m
    {|m|(T) ; m = y}
    ! 0+
    Other approaches:
    Greedy (MP/OMP/etc.)
    Prony (MUSIC/FRI/etc.)
    similar to Frank-Wolfe
    better/less general.

    View Slide

  14. (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) ?
    = mini6=j
    |xi xj
    |

    View Slide

  15. (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).
    = mini6=j
    |xi xj
    |

    View Slide

  16. (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).
    = mini6=j
    |xi xj
    |
    ! [Cand`
    es, Fernandez-G. 2012]
    ! [Fernandez-G.][de Castro 2012]
    Support approximation:

    View Slide

  17. (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).
    = mini6=j
    |xi xj
    |
    ! [Cand`
    es, Fernandez-G. 2012]
    ! [Fernandez-G.][de Castro 2012]
    Support approximation:
    General kernels?
    Support recovery? No separation?

    View Slide

  18. Overview
    • Sparse Spikes Super-resolution
    • Robust Support Recovery
    • Asymptotic Positive Measure Recovery
    • Off-the-Grid Optimization Algorithms
    • Application: Laplace Inversion for TIRF Imaging

    View Slide

  19. = 1/fc
    = 0.6/fc
    Limit Certificate
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    Proposition:


    ⌘ def.
    =
    1 ⇤(y m )
    m solves (P (y)) , ⌘ 2 @|m |(T)
    , |⌘ | 6 1 and ⌘ (xi) = sign(ai)

    View Slide

  20. = 1/fc
    = 0.6/fc
    Limit Certificate
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    ⌘0
    def.
    = argmin
    ⌘= ⇤p
    ||p|| s.t.

    8 i, ⌘(xi) = sign(ai),
    ||⌘||1
    6 1.
    Proposition:


    ⌘ def.
    =
    1 ⇤(y m )
    m solves (P (y)) , ⌘ 2 @|m |(T)
    , |⌘ | 6 1 and ⌘ (xi) = sign(ai)
    Theorem:
    m ! m0 = ma,x
    If ( , ||w||/ ) ! 0, then ⌘ ! ⌘0.

    View Slide

  21. −1
    1 η
    0
    η
    V
    ⌘0
    6= ⌘V
    −1
    1 η
    0
    η
    V
    ⌘0 = ⌘V
    Vanishing Derivative Pre-certificate
    Input measure: m0 = ma,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.

    View Slide

  22. −1
    1 η
    0
    η
    V
    ⌘0
    6= ⌘V
    −1
    1 η
    0
    η
    V
    ⌘0 = ⌘V
    Vanishing Derivative Pre-certificate
    Input measure: m0 = ma,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.
    ⌘V =
    P
    i
    ai
    h'(xi), '(·)i + bi
    h'0(xi), '(·)i
    ! 2N unknown (a, b), 2N equations

    View Slide

  23. −1
    1 η
    0
    η
    V
    ⌘0
    6= ⌘V
    −1
    1 η
    0
    η
    V
    ⌘0 = ⌘V
    Vanishing Derivative Pre-certificate
    Input measure: m0 = ma,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.
    Theorem: ⌘V
    2 ND(m0) =) ⌘V = ⌘0
    ()
    Non-degenerate certificate: ⌘ 2 ND(ma,x) :
    8 t /
    2 {x1, . . . , xN
    }, |⌘(t)| < 1 and 8 i, ⌘00(xi) 6= 0
    ⌘V =
    P
    i
    ai
    h'(xi), '(·)i + bi
    h'0(xi), '(·)i
    ! 2N unknown (a, b), 2N equations

    View Slide

  24. Support Stability Theorem
    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 = ma,x, then
    m =
    PN
    i=1 a?
    i x?
    i
    where ||(x, a) (x?
    , a?)|| = O(||w||)
    ||w||
    Stable
    x

    View Slide

  25. Support Stability Theorem
    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 = ma,x, then
    m =
    PN
    i=1 a?
    i x?
    i
    where ||(x, a) (x?
    , a?)|| = O(||w||)
    ||w||
    Stable
    x
    ||w||
    Unstable
    x

    View Slide

  26. When is Non-degenerate ?
    ⌘V
    ⌘V
    ⌘V
    ⌘V
    Input measure:
    '(0)
    . . .
    ˆ
    '(0) = 1[ fc,fc]
    m0 = ma, x, ! 0
    = 1/fc
    [Cand`
    es, F. Granda] > 1.3

    View Slide

  27. When is Non-degenerate ?
    ⌘V
    ⌘V
    ⌘V
    ⌘V
    Input measure:
    '(0)
    . . .
    ˆ
    '(0) = 1[ fc,fc]
    ⌘V
    ⌘V
    ⌘V
    '(0)
    m0 = ma, x, ! 0
    = 1/fc
    [Cand`
    es, F. Granda] > 1.3

    View Slide

  28. When is Non-degenerate ?
    ⌘V
    ⌘V
    ⌘V
    ⌘V
    Input measure:
    '(0)
    . . .
    ˆ
    '(0) = 1[ fc,fc]
    ⌘V
    ⌘V
    ⌘V
    '(0)
    m0 = ma, x, ! 0
    = 1/fc
    [Cand`
    es, F. Granda] > 1.3
    Valid for:
    Theorem: [Tang, Recht, 2013][Denoyelle 2017]
    9C, ( > C ) =) (⌘V is non degenerate)
    '(x) = e (x ·)2/ 2
    '(x) = 1
    2+(x ·)2

    View Slide

  29. Overview
    • Sparse Spikes Super-resolution
    • Robust Support Recovery
    • Asymptotic Positive Measure Recovery
    • Off-the-Grid Optimization Algorithms
    • Application: Laplace Inversion for TIRF Imaging

    View Slide

  30. Super-resolution for Positive Measures
    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
    Input measure: m0 = ma,x where a 2 RN
    +
    .
    '(x) = (ei`x)fc
    `= fc
    for N 6 fc and ⇢ small enough, ⌘S
    2 ND(m0).
    -1
    1
    ⌘S
    -1
    1 ⌘S

    View Slide

  31. Super-resolution for Positive Measures
    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
    Input measure: m0 = ma,x where a 2 RN
    +
    .
    '(x) = (ei`x)fc
    `= fc
    for N 6 fc and ⇢ small enough, ⌘S
    2 ND(m0).
    ! Extends to sampled Gaussian [Schiebinger et al 2015]
    -1
    1
    ⌘S
    -1
    1 ⌘S

    View Slide

  32. Super-resolution for Positive Measures
    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
    Input measure: m0 = ma,x where a 2 RN
    +
    .
    [Morgenshtern, Cand`
    es, 2015] discrete `1 robustness.
    [Demanet, Nguyen, 2015] discrete `0 robustness.
    ! behavior as 8 i, xi
    ! 0 ?
    '(x) = (ei`x)fc
    `= fc
    for N 6 fc and ⇢ small enough, ⌘S
    2 ND(m0).
    ! Extends to sampled Gaussian [Schiebinger et al 2015]
    -1
    1
    ⌘S
    -1
    1 ⌘S

    View Slide

  33. Super-resolution for Positive Measures
    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
    Input measure: m0 = ma,x where a 2 RN
    +
    .
    [Morgenshtern, Cand`
    es, 2015] discrete `1 robustness.
    [Demanet, Nguyen, 2015] discrete `0 robustness.
    ! behavior as 8 i, xi
    ! 0 ?
    '(x) = (ei`x)fc
    `= fc
    for N 6 fc and ⇢ small enough, ⌘S
    2 ND(m0).
    ! Extends to sampled Gaussian [Schiebinger et al 2015]
    -1
    1
    ⌘S
    -1
    1 ⌘S

    View Slide

  34. Asymptotic of Vanishing Certificate
    1
    ⌘V
    Vanishing Derivative pre-certificate:
    ⌘V
    def.
    = argmin
    ⌘= ⇤p
    ||p||
    m0 = ma, x where ! 0
    s.t. 8 i,

    ⌘( xi) = 1,
    ⌘0( xi) = 0.
    Valid only in 1-D, i.e. T = R or T = R/Z.

    View Slide

  35. Asymptotic of Vanishing Certificate
    1
    ⌘V
    1
    1
    1
    ⌘V
    ⌘V
    ⌘W
    Vanishing Derivative pre-certificate:
    ⌘V
    def.
    = argmin
    ⌘= ⇤p
    ||p||
    m0 = ma, x where ! 0
    s.t. 8 i,

    ⌘( xi) = 1,
    ⌘0( xi) = 0.
    Valid only in 1-D, i.e. T = R or T = R/Z.

    View Slide

  36. 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 = ma, x where ! 0
    s.t. 8 i,

    ⌘( xi) = 1,
    ⌘0( xi) = 0.
    Valid only in 1-D, i.e. T = R or T = R/Z.

    View Slide

  37. 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

    View Slide

  38. 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)

    View Slide

  39. Asymptotic Robustness
    ||w
    ||
    N = 2
    ||w
    ||
    N = 1
    Theorem:
    the solution of P (y) for y = (m0) + w is
    for w
    ,
    w
    2N 1 , 2N 1
    = O(1)
    PN
    i=1
    a?
    i x?
    i
    If ⌘W
    2 NDN , letting m0 = ma, x, then
    where ||(x, a) (x
    ?
    , a
    ?)|| = O

    ||w|| +
    2N 1

    [Denoyelle, D., P. 2015]
    ! signal/noise ⇠ 1/t2N 1 for super-resolution.

    View Slide

  40. Asymptotic Robustness
    ||w
    ||
    N = 2
    ||w
    ||
    N = 1
    Theorem:
    the solution of P (y) for y = (m0) + w is
    for w
    ,
    w
    2N 1 , 2N 1
    = O(1)
    PN
    i=1
    a?
    i x?
    i
    If ⌘W
    2 NDN , letting m0 = ma, x, then
    where ||(x, a) (x
    ?
    , a
    ?)|| = O

    ||w|| +
    2N 1

    [Denoyelle, D., P. 2015]
    ! signal/noise ⇠ 1/t2N 1 for super-resolution.
    ! Extends to clusters:
    ⇠ t
    ⇠ t
    ⇠ t

    View Slide

  41. Asymptotic Robustness
    ||w
    ||
    N = 2
    ||w
    ||
    N = 1
    Theorem:
    the solution of P (y) for y = (m0) + w is
    for w
    ,
    w
    2N 1 , 2N 1
    = O(1)
    PN
    i=1
    a?
    i x?
    i
    If ⌘W
    2 NDN , letting m0 = ma, x, then
    where ||(x, a) (x
    ?
    , a
    ?)|| = O

    ||w|| +
    2N 1

    [Denoyelle, D., P. 2015]
    ! signal/noise ⇠ 1/t2N 1 for super-resolution.
    ! Extends to clusters:
    ⇠ t
    ⇠ t
    ⇠ t
    [Poon, Peyr´
    e 2017]
    signal/noise ⇠ 1/t4
    ! Extends in dimension > 2 for N = 2
    N = 2 N = 3 N = 4 N = 5
    N = 2 N = 3 N = 4 N =
    N = 2 N = 2 N = 2 N =
    N = 1 N = 2 N = 2 N =
    Gaussian MEG-EEG

    View Slide

  42. When is Non-degenerate ?
    ⌘W
    Proposition: one has ⌘(2N)
    W
    (0) < 0.
    ! “locally” non-degenerate.
    (for convolutions)

    View Slide

  43. When is Non-degenerate ?
    ⌘W
    Proposition: one has ⌘(2N)
    W
    (0) < 0.
    ! “locally” non-degenerate.
    ˆ
    '
    ⌘W
    ⌘W
    ⌘W
    N = 2
    N = 3
    N = 4
    '(0)
    (for convolutions)

    View Slide

  44. Gaussian Deconvolution
    Gaussian convolution:
    Proposition: ⌘W (x) = e x2
    4 2
    N 1
    X
    k=0
    (x/2 )2k
    k!
    In particular, ⌘W is non-degenerate.
    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
    ⌘W
    ⌘W ⌘W
    ⌘W
    '(x) = e |x ·|2
    2 2 (m) def.
    =
    Z
    '(x)dm(x)
    '(0)

    View Slide

  45. Overview
    • Sparse Spikes Super-resolution
    • Robust Support Recovery
    • Asymptotic Positive Measure Recovery
    • Off-the-Grid Optimization Algorithms
    • Application: Laplace Inversion for TIRF Imaging

    View Slide

  46. Algorithms
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    D (y)
    = sup
    || ⇤p||1
    61
    hp, yi
    2
    ||p||2
    Primal Dual
    ! 1-dimensional ! 1-many constraints

    View Slide

  47. Algorithms
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    D (y)
    = sup
    || ⇤p||1
    61
    hp, yi
    2
    ||p||2
    Primal Dual
    ! 1-dimensional ! 1-many constraints
    Algorithms:
    “-”: artifacts, slow.
    ! Lasso/Basis-Pursuit: discretize m. [Chen, Donoho, Saunders, 99]
    [Tibshirani, 96]

    View Slide

  48. Algorithms
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    D (y)
    = sup
    || ⇤p||1
    61
    hp, yi
    2
    ||p||2
    Primal Dual
    ! 1-dimensional ! 1-many constraints
    Algorithms:
    “-”: only works for Fourier.
    “-”: artifacts, slow.
    ! SDP-represent D . [Cand`
    es, Fernandez-G. 2012]
    ! Lasso/Basis-Pursuit: discretize m. [Chen, Donoho, Saunders, 99]
    [Tibshirani, 96]

    View Slide

  49. Algorithms
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    D (y)
    = sup
    || ⇤p||1
    61
    hp, yi
    2
    ||p||2
    Primal Dual
    ! 1-dimensional ! 1-many constraints
    Algorithms:
    “-”: only works for Fourier.
    “-”: artifacts, slow.
    ! SDP-represent D . [Cand`
    es, Fernandez-G. 2012]
    [Bredies,Pikkarainen 2010]
    ! Frank-Wolfe on P .
    ! Lasso/Basis-Pursuit: discretize m. [Chen, Donoho, Saunders, 99]
    [Tibshirani, 96]

    View Slide

  50. Algorithms
    min
    m
    |m|(T) +
    1
    2
    || m y||2
    P (y)
    D (y)
    = sup
    || ⇤p||1
    61
    hp, yi
    2
    ||p||2
    Primal Dual
    ! 1-dimensional ! 1-many constraints
    Algorithms:
    “-”: only works for Fourier.
    “-”: artifacts, slow.
    ! SDP-represent D . [Cand`
    es, Fernandez-G. 2012]
    [Bredies,Pikkarainen 2010]
    ! Frank-Wolfe on P .
    ! Lasso/Basis-Pursuit: discretize m. [Chen, Donoho, Saunders, 99]
    [Tibshirani, 96]
    Competitors: Prony’s methods (MUSIC, ESPRIT, FRI).
    “+”: always works when w = 0, less sensitive to sign.
    “-”: only for specific ' (e.g. Fourier), non trivial in 2-D.

    View Slide

  51. y
    zk
    Iterative Soft Thresholding Algorithm
    Computation grid z = (zk)K
    k=1
    .
    Basis-pursuit / Lasso:
    ¯a def.
    =
    P
    k
    ak'(zk)
    ! Force m = ma,z in (P (y)).
    min
    a2RK
    1
    2
    ||y ¯a||2 + ||a||1

    View Slide

  52. y
    zk
    Iterative Soft Thresholding Algorithm
    Computation grid z = (zk)K
    k=1
    .
    Basis-pursuit / Lasso:
    Forward-Backward [Lions, Mercier 1979] (a.k.a. ISTA)
    ¯a def.
    =
    P
    k
    ak'(zk)
    ! Force m = ma,z in (P (y)).
    min
    a2RK
    1
    2
    ||y ¯a||2 + ||a||1
    S (a)
    a
    a(`+1) def.
    = S⌧ (a(`) ⌧ ¯⇤(¯a(`) y)) ⌧ < 2
    ||¯||2

    View Slide

  53. y
    zk
    Iterative Soft Thresholding Algorithm
    Computation grid z = (zk)K
    k=1
    .
    Basis-pursuit / Lasso:
    m0
    ⇤y
    Forward-Backward [Lions, Mercier 1979] (a.k.a. ISTA)
    ¯a def.
    =
    P
    k
    ak'(zk)
    ! Force m = ma,z in (P (y)).
    ! Approximate Diracs by density ! post-processing.
    min
    a2RK
    1
    2
    ||y ¯a||2 + ||a||1
    S (a)
    a
    a(`+1) def.
    = S⌧ (a(`) ⌧ ¯⇤(¯a(`) y)) ⌧ < 2
    ||¯||2
    ! Slow for large K, ¯ ill-posed.

    View Slide

  54. Frank-Wolfe Based Methods
    Initialize: (a(0), x(0)) def.
    = (;, ;)
    min
    m
    1
    2
    || m y||2 + |m|(T) , min
    a,x
    1
    2
    || ma,x y||2 + |a|1
    ` 0

    View Slide

  55. Frank-Wolfe Based Methods
    Initialize: (a(0), x(0)) def.
    = (;, ;)
    min
    m
    1
    2
    || m y||2 + |m|(T) , min
    a,x
    1
    2
    || ma,x y||2 + |a|1
    ⌘(0)
    x?
    Grid generation:
    where ⌘(`) def.
    =
    1 ⇤(y ma(`),x(`)
    )
    x? def.
    = argmaxx
    |⌘(`)(x)|
    ` 0

    View Slide

  56. Frank-Wolfe Based Methods
    Initialize: (a(0), x(0)) def.
    = (;, ;)
    min
    m
    1
    2
    || m y||2 + |m|(T) , min
    a,x
    1
    2
    || ma,x y||2 + |a|1
    Grid deformation: convex non-convex
    (a(`+1), x(`+1)) def.
    = argmin
    a,x
    1
    2
    || ma,x y|| + ||a||1
    BFGS initialized at ([a(`), 0], [x(`), x?])
    ⌘(0)
    x?
    Grid generation:
    where ⌘(`) def.
    =
    1 ⇤(y ma(`),x(`)
    )
    x? def.
    = argmaxx
    |⌘(`)(x)|
    ` 0

    View Slide

  57. Frank-Wolfe Based Methods
    Initialize: (a(0), x(0)) def.
    = (;, ;)
    min
    m
    1
    2
    || m y||2 + |m|(T) , min
    a,x
    1
    2
    || ma,x y||2 + |a|1
    Grid deformation: convex non-convex
    (a(`+1), x(`+1)) def.
    = argmin
    a,x
    1
    2
    || ma,x y|| + ||a||1
    BFGS initialized at ([a(`), 0], [x(`), x?])
    ⌘(3)
    ⌘(0)
    x?
    Grid generation:
    where ⌘(`) def.
    =
    1 ⇤(y ma(`),x(`)
    )
    x? def.
    = argmaxx
    |⌘(`)(x)| ⌘(1) x?
    ⌘(2)
    x?
    ` 0

    View Slide

  58. Frank-Wolfe Based Methods
    Initialize: (a(0), x(0)) def.
    = (;, ;)
    ! Without moving x?: Frank-Wolfe, convergent.
    ! Non-convex update: still convergent, surprisingly e cient.
    min
    m
    1
    2
    || m y||2 + |m|(T) , min
    a,x
    1
    2
    || ma,x y||2 + |a|1
    Grid deformation: convex non-convex
    (a(`+1), x(`+1)) def.
    = argmin
    a,x
    1
    2
    || ma,x y|| + ||a||1
    BFGS initialized at ([a(`), 0], [x(`), x?])
    ⌘(3)
    ⌘(0)
    x?
    Grid generation:
    where ⌘(`) def.
    =
    1 ⇤(y ma(`),x(`)
    )
    x? def.
    = argmaxx
    |⌘(`)(x)| ⌘(1) x?
    ⌘(2)
    x?
    ` 0

    View Slide

  59. Overview
    • Sparse Spikes Super-resolution
    • Robust Support Recovery
    • Asymptotic Positive Measure Recovery
    • Off-the-Grid Optimization Algorithms
    • Application: Laplace Inversion for TIRF Imaging

    View Slide

  60. Laplace Transform Inversion
    Laplace transform:
    x = 2
    x = 20
    t
    (m1)
    (m2)
    x
    m1
    t
    m2
    x
    [with E. Soubies]
    (m) def.
    =
    Z
    '(x)dm(x)
    '(x) = e x·
    '(x)

    View Slide

  61. Laplace Transform Inversion
    Laplace transform:
    x = 2
    x = 20
    t
    (m1)
    (m2)
    x
    m1
    t
    m2
    x
    [with E. Soubies]
    (m) def.
    =
    Z
    '(x)dm(x)
    '(x) = e x·
    '(x)
    Total internal reflection fluorescence microscopy (TIRFM)
    [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
    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
    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 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
    C
    D
    perimental validation of the multiangle TIRF model. (A) Schema
    em designed to create a slope of fluorescent beads. (B) Overlay
    mum intensity projection of image stack acquired with WF and
    ✓(t)
    y(t)
    ! multiple angles ✓(t).
    light
    depth x
    cell
    y(t) = m(t)
    ✓(t)

    View Slide

  62. Laplace Transform Inversion
    Laplace transform:
    x = 2
    x = 20
    t
    (m1)
    (m2)
    x
    m1
    t
    m2
    x
    [with E. Soubies]
    N = 1
    N = 3
    ⌘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
    (m) def.
    =
    Z
    '(x)dm(x)
    '(x) = e x·
    '(x)
    Total internal reflection fluorescence microscopy (TIRFM)
    [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
    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
    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 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
    C
    D
    perimental validation of the multiangle TIRF model. (A) Schema
    em designed to create a slope of fluorescent beads. (B) Overlay
    mum intensity projection of image stack acquired with WF and
    ✓(t)
    y(t)
    ! multiple angles ✓(t).
    light
    depth x
    cell
    y(t) = m(t)
    ✓(t)

    View Slide

  63. Spikes Detection Benchmark
    R =
    TP
    TP + FN
    P =
    TP
    TP + FP
    ||w||
    FW+BFGS ISTA ISTA+Postprocessing
    R
    P
    R
    P

    View Slide

  64. Palm-TIRF Hybridisation
    x x
    y
    z
    x x
    y
    z
    Shot #1
    ...
    Shot #2
    BLASSO
    ...
    Kernel: '(x, y, z) = k(||(x, y) ·||)e z·
    fuse
    Gaussian
    Laplace
    Gaussian
    Laplace
    (joint work with Emmanuel Soubi`
    es)

    View Slide

  65. Conclusion
    varying the azimuth φ during the ex
    modeled by the following expression:
    gðθÞ =
    Z2π
    0
    Z∞
    0
    Z∞
    −∞
    Iðz; α; φÞρ

    θ
    Ω=
    A
    B
    C
    D
    Fig. 2. Experimental validation of the mult
    of the system designed to create a slope of
    of the maximum intensity projection of imag
    TIRF illumination. (Scale bar: 5 μm.) The evolu
    illumination angle θ of two selected beads a
    sponding fitting theoretical model (continu
    depth (respectively 10 and 89 nm). (D) Depth o
    fitting the theoretical TIRF model (in red) an
    estimated by fitting a Gaussian model in the
    [Boulanger et al. 2014]
    Super-resolution should be o↵-the-grid!
    Theory:
    `2 errors meaningless.
    6= compressed sensing.
    Super-resolution , signal/noise vs t

    View Slide

  66. Conclusion
    varying the azimuth φ during the ex
    modeled by the following expression:
    gðθÞ =
    Z2π
    0
    Z∞
    0
    Z∞
    −∞
    Iðz; α; φÞρ

    θ
    Ω=
    A
    B
    C
    D
    Fig. 2. Experimental validation of the mult
    of the system designed to create a slope of
    of the maximum intensity projection of imag
    TIRF illumination. (Scale bar: 5 μm.) The evolu
    illumination angle θ of two selected beads a
    sponding fitting theoretical model (continu
    depth (respectively 10 and 89 nm). (D) Depth o
    fitting the theoretical TIRF model (in red) an
    estimated by fitting a Gaussian model in the
    [Boulanger et al. 2014]
    Super-resolution should be o↵-the-grid!
    Theory:
    `2 errors meaningless.
    6= compressed sensing.
    Practice:
    Adaptive grid refinement.
    Non-convex step crucial.
    Surprisingly e cient.
    Super-resolution , signal/noise vs t

    View Slide

  67. Conclusion
    Open problem: other regularizations (e.g. piecewise constant) ?
    see [Chambolle, Duval, Peyr´
    e, Poon 2016] for TV denoising.
    varying the azimuth φ during the ex
    modeled by the following expression:
    gðθÞ =
    Z2π
    0
    Z∞
    0
    Z∞
    −∞
    Iðz; α; φÞρ

    θ
    Ω=
    A
    B
    C
    D
    Fig. 2. Experimental validation of the mult
    of the system designed to create a slope of
    of the maximum intensity projection of imag
    TIRF illumination. (Scale bar: 5 μm.) The evolu
    illumination angle θ of two selected beads a
    sponding fitting theoretical model (continu
    depth (respectively 10 and 89 nm). (D) Depth o
    fitting the theoretical TIRF model (in red) an
    estimated by fitting a Gaussian model in the
    [Boulanger et al. 2014]
    Super-resolution should be o↵-the-grid!
    Theory:
    `2 errors meaningless.
    6= compressed sensing.
    Practice:
    Adaptive grid refinement.
    Non-convex step crucial.
    Surprisingly e cient.
    Super-resolution , signal/noise vs t

    View Slide