Signal Processing Course: Sparse Regularization of Inverse Problems

Slide 1

Slide 1 text

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

Slide 2

Slide 2 text

Overview • Inverse Problems Regularization • Sparse Synthesis Regularization • Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution

Slide 3

Slide 3 text

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

Slide 4

Slide 4 text

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

Slide 5

Slide 5 text

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

Slide 6

Slide 6 text

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

Slide 7

Slide 7 text

Kf = (p k ) 1 k K Inverse Problem in Medical Imaging

Slide 8

Slide 8 text

Magnetic resonance imaging (MRI): Kf = (p k ) 1 k K Kf = ( ˆ f( )) Inverse Problem in Medical Imaging ˆ f

Slide 9

Slide 9 text

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

Slide 10

Slide 10 text

Noisy measurements: y = Kf0 + w. Prior model: J : RQ R assigns a score to images. Inverse Problem Regularization

Slide 11

Slide 11 text

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

Slide 12

Slide 12 text

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

Slide 13

Slide 13 text

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

Slide 14

Slide 14 text

J(f) = || f(x)||2dx Smooth and Cartoon Priors | f|2

Slide 15

Slide 15 text

J(f) = || f(x)||2dx J(f) = || f(x)||dx J(f) = R length(Ct )dt Smooth and Cartoon Priors | f|2 | f|

Slide 16

Slide 16 text

Inpainting Example Input y = Kf0 + w Sobolev Total variation

Slide 17

Slide 17 text

Overview • Inverse Problems Regularization • Sparse Synthesis Regularization • Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution

Slide 18

Slide 18 text

Q N Dictionary = ( m ) m RQ N , N Q. Redundant Dictionaries

Slide 19

Slide 19 text

m = ei ·, m frequency Fourier: Q N Dictionary = ( m ) m RQ N , N Q. Redundant Dictionaries

Slide 20

Slide 20 text

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

Slide 21

Slide 21 text

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

Slide 22

Slide 22 text

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

Slide 23

Slide 23 text

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

Slide 24

Slide 24 text

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

Slide 25

Slide 25 text

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

Slide 26

Slide 26 text

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

Slide 27

Slide 27 text

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

Slide 28

Slide 28 text

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

Slide 29

Slide 29 text

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

Slide 30

Slide 30 text

L1 Regularization coe cients x0 RN

Slide 31

Slide 31 text

L1 Regularization coe cients image x0 RN f0 = x0 RQ

Slide 32

Slide 32 text

L1 Regularization observations w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP

Slide 33

Slide 33 text

L1 Regularization observations = K ⇥ ⇥ RP N w coe cients image K x0 RN f0 = x0 RQ y = Kf0 + w RP

Slide 34

Slide 34 text

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

Slide 35

Slide 35 text

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

Slide 36

Slide 36 text

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

Slide 37

Slide 37 text

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

Slide 38

Slide 38 text

Regularization Data ﬁdelity y = x0 + w Noisy measurements: x argmin x RQ 1 2 ||y x||2 + ||x||1 Noisy Sparse Regularization

Slide 39

Slide 39 text

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

Slide 40

Slide 40 text

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

Slide 41

Slide 41 text

Overview • Inverse Problems Regularization • Sparse Synthesis Regularization • Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution

Slide 42

Slide 42 text

Image De-blurring Original f0 y = h f0 + w

Slide 43

Slide 43 text

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

Slide 44

Slide 44 text

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

Slide 45

Slide 45 text

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

Slide 46

Slide 46 text

K y = Kf0 + w Measures: Inpainting Problem (Kf)(x) = 0 if x , f(x) if x / .

Slide 47

Slide 47 text

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

Slide 48

Slide 48 text

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

Slide 49

Slide 49 text

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.

Slide 50

Slide 50 text

Examples of Decompositions

Slide 51

Slide 51 text

Cartoon+Texture Separation

Slide 52

Slide 52 text

Overview • Inverse Problems Regularization • Sparse Synthesis Regularization • Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution

Slide 53

Slide 53 text

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

Slide 54

Slide 54 text

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

Slide 55

Slide 55 text

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(·)

Slide 56

Slide 56 text

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(·)

Slide 57

Slide 57 text

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(·)

Slide 58

Slide 58 text

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(·)

Slide 59

Slide 59 text

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(·)

Slide 60

Slide 60 text

Overview • Inverse Problems Regularization • Sparse Synthesis Regularization • Examples: Sparse Wavelet Regularizations • Iterative Soft Thresholding • Sparse Seismic Deconvolution

Slide 61

Slide 61 text

Seismic Imaging

Slide 62

Slide 62 text

1D propagation convolution P 1D Idealization Initial condition: “wavelet” f f = h f = band pass ﬁlter h h(x) ˆ h( ) y = f0 h + w

Slide 63

Slide 63 text

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+

Slide 64

Slide 64 text

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

Slide 65

Slide 65 text

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

Slide 66

Slide 66 text

oracle, minimize ||f0 f || Numerical Example y = f h + w f SNR(f0, f ) f0 Choosing optimal :

Slide 67

Slide 67 text

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