Sparsity and optimality of splines: Deterministic vs. statistical justification Michael Unser Biomedical Imaging Group EPFL, Lausanne, Switzerland Mathematics and Image Analysis (MIA), 18-20 January 2016, Paris, France. OUTLINE ■ Sparsity and signal reconstruction ■ Inverse problems in bio-imaging ■ Compressed sensing: towards a continuous-domain formulation ■ Deterministic formulation ■ Splines and operators ■ New optimality results for generalized TV ■ Statistical formulation ■ Sparse stochastic processes ■ Derivation of MAP estimators 2
Inverse problems in bio-imaging y = Hs + n Linear forward model noise H n Integral operator s Problem: recover s from noisy measurements y Reconstruction as an optimization problem s rec = arg min k y � Hs k 2 + λ k Ls k p p = 1 , 2 , 2 p | {z } | {z } data consistency regularization − log Prob( s ) : prior likelihood 3 Inverse problems in imaging: Current status Higher reconstruction quality : Sparsity-promoting schemes almost sys- tematically outperform the classical linear reconstruction methods in MRI, x-ray tomography, deconvolution microscopy, etc... (Lustig et al. 2007) Increased complexity : Resolution of linear inverse problems using ` 1 regularization requires more sophisticated algorithms (iterative and non- linear); efficient solutions (FISTA, ADMM) have emerged during the past decade. (Chambolle 2004; Figueiredo 2004; Beck-Teboule 2009; Boyd 2011) The paradigm is supported by the theory of compressed sensing (Candes-Romberg-Tao; Donoho, 2006) Outstanding research issues Beyond ` 1 and TV: Connection with statistical modeling & learning Beyond matrix algebra: Continuous-domain formulation Guarantees of optimality (either deterministic or statistical) 4
Sparsity and continuous-domain modeling Compressed sensing (CS) Generalized sampling and infinite-dimensional CS (Adcock-Hansen, 2011) Xampling: CS of analog signals (Eldar, 2011) Splines and approximation theory L 1 splines (Fisher-Jerome, 1975) Locally-adaptive regression splines (Mammen-van de Geer, 1997) Generalized TV (Steidl et al. 2005; Bredies et al. 2010) Statistical modeling (Unser et al. 2011-2014) Sparse stochastic processes 5 The spline connection Photo courtesy of Carl De Boor 6
Splines are intrinsically sparse L {·} : admissible differential operator δ ( · − x 0 ) : Dirac impulse shifted by x 0 ∈ R d Definition The function s : R d → R is a (non-uniform) L -spline with knots ( x k ) K k =1 if K X L { s } = a k δ ( · − x k ) = w δ : spline’s innovation k =1 L = d Spline theory: (Schultz-Varga, 1967) d x a k x k x k +1 FIR signal processing: Innovation variables (2 K ) (Vetterli et al., 2002) Location of singularities (knots) : { x k } K k =1 Strength of singularities (linear weights): { a k } K k =1 7 Splines and operators Definition A linear operator L : X → Y , where X ⊃ S ( R d ) and Y are appropriate sub- spaces of S 0 ( R d ) , is called spline-admissible if 1. it is linear shift-invariant (LSI); 2. its null space N L = { p ∈ X : L { p } = 0 } is finite-dimensional of size N 0 ; 3. there exists a function ρ L : R d → R of slow growth (the Green’s function of L ) such that L { ρ L } = δ . Structure of null space: N L = span { p n } N 0 n =1 Admits some basis p = ( p 1 , · · · , p N 0 ) Only includes elements of the form x m e j h ω 0 , x i with | m | = P d i =1 m i ≤ n 0 8
Spline synthesis: example L = D = d N D = span { p 1 } , p 1 ( x ) = 1 d x ρ D ( x ) = + ( x ) : Heaviside function X w δ ( x ) = a k δ ( x − x k ) k x x 1 X s ( x ) = b 1 p 1 ( x ) + + ( x − x k ) a k k a 1 b 1 x 9 Spline synthesis: generalization L : spline admissible operator (LSI) ρ L ( x ) : Green’s function of L N L = span { p n } N 0 n =1 X w δ ( x ) = a k δ ( x − x k ) Spline’s innovation: k N 0 X X s ( x ) = a k ρ L ( x − x k ) + b n p n ( x ) ⇒ n =1 k Requires specification of boundary conditions a 1 x 10
Principled operator-based approach Biorthogonal basis of N L = span { p n } N 0 n =1 φ = ( φ 1 , · · · , φ N 0 ) such that h φ m , p n i = δ m,n N 0 X p = h φ n , p i p n for all p 2 N L n =1 Operator-based spline synthesis h s, φ n i = b n , n = 1 , · · · , N 0 Boundary conditions: X L { s } = w δ = a k δ ( · � x k ) Spline’s innovation: k N 0 X s ( x ) = L − 1 φ { w δ } ( x ) + b n p n ( x ) n =1 Existence of L − 1 φ as a stable right-inverse of L ? (see Theorem 1 ) LL − 1 φ w = w h φ , L − 1 φ w i = 0 11 Beyond splines ... Photo courtesy of Carl De Boor 12
From Dirac impulses to Borel measures S ( R d ) : Schwartz’s space of smooth and rapidly decaying test functions on R d S 0 ( R d ) : Schwartz’s space of tempered distributions Space of real-valued, countably additive Borel measures on R d � 0 = M ( R d ) = C 0 ( R d ) w 2 S 0 ( R d ) : k w k TV = � � sup h w, ϕ i < 1 , ϕ 2 S ( R d ): k ϕ k ∞ =1 R where w : ϕ 7! h w, ϕ i = R d ϕ ( r )d w ( r ) Equivalent definition of “total variation” norm k w k TV = h w, ϕ i sup ϕ 2 C 0 ( R d ): k ϕ k ∞ =1 Basic inclusions δ ( · � x 0 ) 2 M ( R d ) with k δ ( · � x 0 ) k TV = 1 for any x 0 2 R d k f k TV = k f k L 1 ( R d ) for all f 2 L 1 ( R d ) L 1 ( R d ) ✓ M ( R d ) ) 13 Generalized Beppo-Levi spaces L : spline-admissible operator Generalized “total variation” semi-norm (gTV) gTV( f ) = k L { f } k TV Generalized Beppo-Levi spaces f : R d ! R � k L f k TV < 1 M L ( R d ) = � � f 2 M L ( R d ) , L { f } 2 M ( R d ) Classical Beppo-Levi spaces: ( M ( R d ) , L) → ( L p ( R ) , D n ) (Deny-Lions, 1954) Inclusion of non-uniform L -splines N 0 X X X s = a k ρ L ( · � x k ) + b n p n ) L { s } = a k δ ( · � x k ) n =1 k k X gTV( s ) = k L { s } k TV = | a k | = k a k ` 1 k 14
New optimality result: universality of splines L : spline-admissible operator � M L ( R ) = f : gTV( f ) = k L { f } k TV = sup h L { f } , ϕ i < 1 k ϕ k ∞ 1 Generalized total variation : gTV( f ) = k L { f } k L 1 when L { f } 2 L 1 ( R d ) Linear measurement operator M L ( R ) ! R M : f 7! z = ν ( f ) , z m = h f, ν m i Theorem [U.-Fageot-Ward, 2015]: The generic linear-inverse problem k y � ν ( f ) k 2 � � min 2 + λ k L { f } k TV f ∈ M L ( R ) K N 0 X X admits a global solution of the form f ( x ) = a k ρ L ( x � x k ) + b n p n ( x ) n =1 k =1 with K < M , which is a non-uniform L -spline with knots ( x k ) K k =1 . 15 Optimality result for Dirac measures F : linear continuous map M ( R d ) ! R M C : convex compact subset of R M Generic constrained TV minimization problem V = arg w ∈ M ( R d ) : F ( w ) ∈ C k w k TV min Generalized Fisher-Jerome theorem The solution set V is a convex, weak ⇤ -compact subset of M ( R d ) with extremal points of the form K X w δ = a k δ ( · � x k ) k =1 with K M and x k 2 R d . Jerome-Fisher, 1975: Compact domain & scalar intervals 16
Existence of stable right-inverse operator L ∞ ,n 0 ( R d ) = { f : R d ! R : sup x ∈ R d ( | f ( x ) | (1 + k x k ) n 0 ) < + 1 } Theorem 1 [U.-Fageot-Ward, preprint] Let L be spline-admissible operator with a N 0 -dimensional null space N L ✓ L ∞ , − n 0 ( R d ) such that p = P N 0 n =1 h p, φ n i p n for all p 2 N L . Then, there exists a unique and sta- ble operator L − 1 φ : M ( R d ) ! L ∞ , − n 0 ( R d ) such that, for all w 2 M ( R d ) , • Right-inverse property: LL − 1 φ w = w , • Boundary conditions: h φ , L − 1 φ w i = 0 with φ = ( φ 1 , · · · , φ N 0 ) . Its generalized impulse response g φ ( x , y ) = L − 1 φ { δ ( · � y ) } ( x ) is given by N 0 X g φ ( x , y ) = ρ L ( x � y ) � p n ( x ) q n ( y ) n =1 with ρ L such that L { ρ L } = δ and q n ( y ) = h φ n , ρ L ( · � y ) i . 17 Characterization of generalized Beppo-Levi spaces Regularization operator L : M L ( R d ) ! M ( R d ) f 2 M L ( R d ) , gTV( f ) = k L { f } k TV < 1 Theorem 2 [U.-Fageot-Ward, preprint] Let L be a spline-admissible operator that admits a stable right-inverse L � 1 of the φ form specified by Theorem 1. Then, any f 2 M L ( R d ) has a unique representation as f = L � 1 φ w + p, where w = L { f } 2 M ( R d ) and p = P N 0 � 0 . n =1 h φ n , f i p n 2 N L with φ n 2 � M L ( R d ) Moreover, M L ( R d ) ✓ L 1 , � n 0 ( R d ) and is a Banach space equipped with the norm k f k M L , φ = k L f k TV + kh f, φ ik 2 . M L ( R d ) = M L , φ ( R d ) � N L Generalized Beppo-Levi space: M L , φ ( R d ) = � f 2 M L ( R d ) : h φ , f i = 0 � p 2 M L ( R d ) : L { p } = 0 N L = 18
Link with sparse stochastic processes EDEE Course 19 Random spline: archetype of sparse signal non-uniform spline of degree 0 0 2 4 6 8 10 X D s ( t ) = a n δ ( t − t n ) = w ( t ) n Random weights { a n } i.i.d. and random knots { t n } (Poisson with rate λ ) Anti-derivative operators Z t Shift-invariant solution: D − 1 ϕ ( t ) = ( + ∗ ϕ )( t ) = ϕ ( τ )d τ −∞ Z t Scale-invariant solution: D − 1 (see Theorem 1 with φ 1 = δ ) φ 1 ϕ ( t ) = ϕ ( τ )d τ 0 20
Recommend
More recommend
