Signal Processing Course: Sparse Regularization of Inverse Problems

Signal Processing Course: Sparse Regularization of Inverse Problems

E34ded36efe4b7abb12510d4e525fee8?s=128

Gabriel Peyré

January 01, 2012
Tweet

Transcript

  1. Sparse Regularization Of Inverse Problems Gabriel Peyré www.numerical-tours.com

  2. Overview • Inverse Problems Regularization • Sparse Synthesis Regularization •

    Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution
  3. Forward model: Observations Operator Noise (Unknown) Input : RQ RP

    Inverse Problems y = K f0 + w RP
  4. Forward model: Observations Operator Noise (Unknown) Input : RQ RP

    Denoising: K = Id Q, P = Q. Inverse Problems y = K f0 + w RP
  5. Forward model: Observations Operator Noise (Unknown) Input (Kf)(x) = 0

    if x , f(x) if x / . K : RQ RP Denoising: K = Id Q, P = Q. Inpainting: set of missing pixels, P = Q | |. Inverse Problems y = K f0 + w RP
  6. Forward model: Observations Operator Noise (Unknown) Input (Kf)(x) = 0

    if x , f(x) if x / . K : RQ RP Denoising: K = Id Q, P = Q. Inpainting: set of missing pixels, P = Q | |. Super-resolution: Kf = (f k) , P = Q/ . Inverse Problems K y = K f0 + w RP
  7. Kf = (p k ) 1 k K Inverse Problem

    in Medical Imaging
  8. Magnetic resonance imaging (MRI): Kf = (p k ) 1

    k K Kf = ( ˆ f( )) Inverse Problem in Medical Imaging ˆ f
  9. Magnetic resonance imaging (MRI): Other examples: MEG, EEG, . .

    . Kf = (p k ) 1 k K Kf = ( ˆ f( )) Inverse Problem in Medical Imaging ˆ f
  10. Noisy measurements: y = Kf0 + w. Prior model: J

    : RQ R assigns a score to images. Inverse Problem Regularization
  11. Noisy measurements: y = Kf0 + w. f argmin f

    RQ 1 2 ||y Kf||2 + J(f) Prior model: J : RQ R assigns a score to images. Inverse Problem Regularization Data fidelity Regularity
  12. Noisy measurements: y = Kf0 + w. Choice of :

    tradeo ||w|| Regularity of f0 J(f0 ) Noise level f argmin f RQ 1 2 ||y Kf||2 + J(f) Prior model: J : RQ R assigns a score to images. Inverse Problem Regularization Data fidelity Regularity
  13. Noisy measurements: y = Kf0 + w. No noise: 0+,

    minimize Choice of : tradeo ||w|| Regularity of f0 J(f0 ) Noise level f argmin f RQ 1 2 ||y Kf||2 + J(f) Prior model: J : RQ R assigns a score to images. f argmin f RQ,Kf=y J(f) Inverse Problem Regularization Data fidelity Regularity
  14. J(f) = || f(x)||2dx Smooth and Cartoon Priors | f|2

  15. J(f) = || f(x)||2dx J(f) = || f(x)||dx J(f) =

    R length(Ct )dt Smooth and Cartoon Priors | f|2 | f|
  16. Inpainting Example Input y = Kf0 + w Sobolev Total

    variation
  17. Overview • Inverse Problems Regularization • Sparse Synthesis Regularization •

    Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution
  18. Q N Dictionary = ( m ) m RQ N

    , N Q. Redundant Dictionaries
  19. m = ei ·, m frequency Fourier: Q N Dictionary

    = ( m ) m RQ N , N Q. Redundant Dictionaries
  20. Wavelets: m = (2 jR x n) m = (j,

    , n) scale orientation position m = ei ·, m frequency Fourier: Q N Dictionary = ( m ) m RQ N , N Q. Redundant Dictionaries = 2 = 1
  21. Wavelets: m = (2 jR x n) m = (j,

    , n) scale orientation position m = ei ·, m frequency Fourier: DCT, Curvelets, bandlets, . . . Q N Dictionary = ( m ) m RQ N , N Q. Redundant Dictionaries = 2 = 1
  22. Synthesis: f = m xm m = x. x f

    Image f = x Coe cients x Wavelets: m = (2 jR x n) m = (j, , n) scale orientation position m = ei ·, m frequency Fourier: DCT, Curvelets, bandlets, . . . Q N Dictionary = ( m ) m RQ N , N Q. Redundant Dictionaries = = 2 = 1
  23. Ideal sparsity: for most m, xm = 0. J0 (x)

    = # {m \ xm = 0} Sparse Priors Image f0 Coe cients x
  24. Ideal sparsity: for most m, xm = 0. J0 (x)

    = # {m \ xm = 0} Sparse approximation: f = x where argmin x 2RN || f0 x ||2 + T 2 J0( x ) Sparse Priors Image f0 Coe cients x
  25. Ideal sparsity: for most m, xm = 0. J0 (x)

    = # {m \ xm = 0} Sparse approximation: f = x where Orthogonal : = = Id N xm = f0, m if | f0, m | > T, 0 otherwise. ST f = ST (f0 ) argmin x 2RN || f0 x ||2 + T 2 J0( x ) Sparse Priors Image f0 Coe cients x
  26. Ideal sparsity: for most m, xm = 0. J0 (x)

    = # {m \ xm = 0} Sparse approximation: f = x where Orthogonal : = = Id N xm = f0, m if | f0, m | > T, 0 otherwise. ST Non-orthogonal : NP-hard. f = ST (f0 ) argmin x 2RN || f0 x ||2 + T 2 J0( x ) Sparse Priors Image f0 Coe cients x
  27. Image with 2 pixels: q = 0 J0 (x) =

    # {m \ xm = 0} J0 (x) = 0 null image. J0 (x) = 1 sparse image. J0 (x) = 2 non-sparse image. x2 Convex Relaxation: L1 Prior x1
  28. Image with 2 pixels: q = 0 q = 1

    q = 2 q = 3/2 q = 1/2 Jq (x) = m |xm |q J0 (x) = # {m \ xm = 0} J0 (x) = 0 null image. J0 (x) = 1 sparse image. J0 (x) = 2 non-sparse image. x2 Convex Relaxation: L1 Prior q priors: (convex for q 1) x1
  29. Image with 2 pixels: q = 0 q = 1

    q = 2 q = 3/2 q = 1/2 Jq (x) = m |xm |q J0 (x) = # {m \ xm = 0} J0 (x) = 0 null image. J0 (x) = 1 sparse image. J0 (x) = 2 non-sparse image. x2 J1 (x) = m |xm | Convex Relaxation: L1 Prior Sparse 1 prior: q priors: (convex for q 1) x1
  30. L1 Regularization coe cients x0 RN

  31. L1 Regularization coe cients image x0 RN f0 = x0

    RQ
  32. L1 Regularization observations w coe cients image K x0 RN

    f0 = x0 RQ y = Kf0 + w RP
  33. L1 Regularization observations = K ⇥ ⇥ RP N w

    coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP
  34. Fidelity Regularization min x RN 1 2 ||y x||2 +

    ||x||1 L1 Regularization Sparse recovery: f = x where x solves observations = K ⇥ ⇥ RP N w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP
  35. x argmin x=y m |xm | x x = y

    Noiseless Sparse Regularization Noiseless measurements: y = x0
  36. x argmin x=y m |xm | x x = y

    x argmin x=y m |xm |2 Noiseless Sparse Regularization x x = y Noiseless measurements: y = x0
  37. x argmin x=y m |xm | x x = y

    x argmin x=y m |xm |2 Noiseless Sparse Regularization Convex linear program. Interior points, cf. [Chen, Donoho, Saunders] “basis pursuit”. Douglas-Rachford splitting, see [Combettes, Pesquet]. x x = y Noiseless measurements: y = x0
  38. Regularization Data fidelity y = x0 + w Noisy measurements:

    x argmin x RQ 1 2 ||y x||2 + ||x||1 Noisy Sparse Regularization
  39. Regularization Data fidelity Equivalence || x = y|| y =

    x0 + w Noisy measurements: x argmin x RQ 1 2 ||y x||2 + ||x||1 x argmin || x y|| ||x||1 Noisy Sparse Regularization x
  40. Iterative soft thresholding Forward-backward splitting Algorithms: Regularization Data fidelity Equivalence

    || x = y|| y = x0 + w Noisy measurements: x argmin x RQ 1 2 ||y x||2 + ||x||1 x argmin || x y|| ||x||1 Noisy Sparse Regularization Nesterov multi-steps schemes. see [Daubechies et al], [Pesquet et al], etc x
  41. Overview • Inverse Problems Regularization • Sparse Synthesis Regularization •

    Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution
  42. Image De-blurring Original f0 y = h f0 + w

  43. f = argmin f RN ||f ⇥ h y||2 +

    ||⇥f||2 ˆ f (⇥) = ˆ h(⇥) |ˆ h(⇥)|2 + |⇥|2 ˆ y(⇥) Sobolev regularization: Image De-blurring Original f0 y = h f0 + w Sobolev SNR=22.7dB
  44. f = argmin f RN ||f ⇥ h y||2 +

    ||⇥f||2 ˆ f (⇥) = ˆ h(⇥) |ˆ h(⇥)|2 + |⇥|2 ˆ y(⇥) Sobolev regularization: = translation invariant wavelets. x argmin x 1 2 ||h ( x) y||2 + ||x||1 f = x where Sparsity Image De-blurring Original f0 y = h f0 + w Sparsity regularization: SNR=24.7dB Sobolev SNR=22.7dB
  45. SNR=23.7dB SNR=24.7dB L2 Sobolev Sparsity Invariant Comparison of Regularizations SNR

    SNR SNR opt opt opt Sparsity regularization Sobolev regularization L2 regularization SNR=21.7dB SNR=22.7dB
  46. K y = Kf0 + w Measures: Inpainting Problem (Kf)(x)

    = 0 if x , f(x) if x / .
  47. Image Separation Model: f = f1 + f2 + w,

    (f1, f2 ) components, w noise.
  48. Image Separation Model: f = f1 + f2 + w,

    (f1, f2 ) components, w noise.
  49. Union dictionary: (x1 , x2 ) argmin x=(x1,x2) RN 1

    2 ||f x||2 + ||x||1 = [ 1, 2 ] RQ (N1+N2) Image Separation Recovered component: fi = ixi . Model: f = f1 + f2 + w, (f1, f2 ) components, w noise.
  50. Examples of Decompositions

  51. Cartoon+Texture Separation

  52. Overview • Inverse Problems Regularization • Sparse Synthesis Regularization •

    Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution
  53. Orthogonal-basis: ⇤ = IdN , x = ⇤ f .

    Regularization-based denoising: x? = argmin x 2RN 1 2 || x y ||2 + J ( x ) Denoising: y = x0 + w 2 RN , K = Id. Sparse regularization: J ( x ) = P m | xm |q (where |a|0 = (a)) Sparse Regularization Denoising
  54. Orthogonal-basis: ⇤ = IdN , x = ⇤ f .

    Regularization-based denoising: x? = argmin x 2RN 1 2 || x y ||2 + J ( x ) Denoising: y = x0 + w 2 RN , K = Id. Sparse regularization: J ( x ) = P m | xm |q x ? m = S q T ( xm) (where |a|0 = (a)) Sparse Regularization Denoising
  55. Sparse regularization: x? 2 argmin x 2RN E ( x

    ) = 1 2 || y x ||2 + || x ||1 E ( x, ˜ x ) = E ( x ) 1 2 || ( x ˜ x )||2 + 1 2 ⌧ || x ˜ x ||2 E (· , ˜ x ) x ˜ x S ⌧ (u) ⌧ < 1/|| ⇤ || Surrogate Functionals Surrogate functional: E(·)
  56. Sparse regularization: x? 2 argmin x 2RN E ( x

    ) = 1 2 || y x ||2 + || x ||1 E ( x, ˜ x ) = E ( x ) 1 2 || ( x ˜ x )||2 + 1 2 ⌧ || x ˜ x ||2 Proof: E ( x, ˜ x ) / 1 2 || u x ||2 + || x ||1+ cst. argmin x E ( x, ˜ x ) = S ⌧ ( u ) Theorem: where u = x ⌧ ⇤( x ˜ x ) E (· , ˜ x ) x ˜ x S ⌧ (u) ⌧ < 1/|| ⇤ || Surrogate Functionals Surrogate functional: E(·)
  57. Algorithm: x ( ` +1) = argmin x E (

    x, x ( ` )) x (`+1) = S 1 ⌧ ( u (`)) u (`) = x (`) ⌧ ⇤( x (`) ⌧y ) Initialize x (0), set ` = 0. x (0) x (1) x (2) Iterative Thresholding x E(·)
  58. Algorithm: x ( ` +1) = argmin x E (

    x, x ( ` )) x (`+1) = S 1 ⌧ ( u (`)) u (`) = x (`) ⌧ ⇤( x (`) ⌧y ) Remark: x (`) 7! u (`) is a gradient descent of || x y ||2 . S 1 `⌧ is the proximal step of associated to || x ||1. Initialize x (0), set ` = 0. x (0) x (1) x (2) Iterative Thresholding x E(·)
  59. Theorem: if ⌧ < 2 / || ⇤ ||, then

    x (`) ! x ?. Algorithm: x ( ` +1) = argmin x E ( x, x ( ` )) x (`+1) = S 1 ⌧ ( u (`)) u (`) = x (`) ⌧ ⇤( x (`) ⌧y ) Remark: x (`) 7! u (`) is a gradient descent of || x y ||2 . S 1 `⌧ is the proximal step of associated to || x ||1. Initialize x (0), set ` = 0. x (0) x (1) x (2) Iterative Thresholding x E(·)
  60. Overview • Inverse Problems Regularization • Sparse Synthesis Regularization •

    Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution
  61. Seismic Imaging

  62. 1D propagation convolution P 1D Idealization Initial condition: “wavelet” f

    f = h f = band pass filter h h(x) ˆ h( ) y = f0 h + w
  63. Pseudo-inverse: Stabilization: ˆ h+(⇥) = 0 if |ˆ h(⇥)| Pseudo

    Inverse ˆ f+( ) = ˆ y( ) h( ) = f+ = h+ ⇥ y = f0 + h+ ⇥ w where ˆ h+( ) = ˆ h( ) 1 ˆ h( ) 1/ˆ h( ) y = h f0 + w f+
  64. f with small ||f||0 Sparse Spikes Deconvolution ⇥m [x] =

    [x m] Sparsity basis: Diracs y = f h + w f = argmin f RN 1 2 ||f ⇥ h y||2 + m |f[m]|.
  65. f with small ||f||0 Sparse Spikes Deconvolution ⇥m [x] =

    [x m] Sparsity basis: Diracs y = f h + w Algorithm: • Inversion: ˜ f(k) = f(k) h ⇥ (h ⇥ f(k) y). f(k+1)[m] = S1 ⇥ ( ˜ f(k)[m]) < 2/|| || = 2/max |ˆ h(⇥)|2 f = argmin f RN 1 2 ||f ⇥ h y||2 + m |f[m]|.
  66. oracle, minimize ||f0 f || Numerical Example y = f

    h + w f SNR(f0, f ) f0 Choosing optimal :
  67. log 10 (E(f(k))/E(f ) 1) log 10 (||f(k) f ||/||f0

    ||) k Not strictly convex = no convergence speed. Convergence Study k Energy: E(f) = 1 2 ||h ⇥ f y||2 + m |f[m]|. Sparse deconvolution: f = argmin f RN E(f).
  68. Conclusion

  69. Conclusion

  70. Conclusion