Signal Processing Course: Sparse Regularization of Inverse Problems

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