Gabriel Peyré
January 01, 2012
1.4k

# Signal Processing Course: Sparse Regularization of Inverse Problems

January 01, 2012

## Transcript

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 ﬁdelity 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 ﬁdelity 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 ﬁdelity Regularity

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

R length(Ct )dt Smooth and Cartoon Priors | f|2 | f|

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

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 ﬁdelity y = x0 + w Noisy measurements:

x argmin x RQ 1 2 ||y x||2 + ||x||1 Noisy Sparse Regularization
39. ### Regularization Data ﬁdelity 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 ﬁdelity 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

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.

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

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

f = h f = band pass ﬁlter 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).