Sparse stochastic processes and biomedical image reconstruction - - PDF document
Sparse stochastic processes and biomedical image reconstruction Michael Unser Biomedical Imaging Group EPFL, Lausanne, Switzerland Joint work with P. Tafti, Q. Sun, A. Amini, M. Guerquin-Kern, E. Bostan, etc. Tutorial presentation, FRUMAN,
Tutorial presentation, FRUMAN, University Aix-Marseille, Feb. 4 2013.
2
noise
2
data consistency
regularization
3
Formal linear solution:
s
: Whitening filter Quadratic regularization (Tikhonov)
Statistical formulation under Gaussian hypothesis Wiener (LMMSE) solution = Gauss MMSE = Gauss MAP
2
data consistency
regularization
Signal covariance: Cs = E{s · sT }
2
Data Log likelihood
s
2
Gaussian prior likelihood
4
2
data consistency
regularization
Wavelet expansion: s = Wv (typically, sparse) Wavelet-domain sparsity-constraint:
with
Iterated shrinkage-thresholding algorithm (ISTA, FISTA, FWISTA)
Iterative reweighted least squares (IRLS) or FISTA
5
1
2
3
6
Constructed by Paul Lévy (circa 1930)
Example: compound Poisson process (piecewise-constant with random jumps)
Generalization of Brownian motion Independent increments: i.i.d. infinitely divisible (from Gaussian to heavy tailed)
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Compound Poisson Gaussian: Brownian motion
SαS (Cauchy)
7
φ(x)
Poisson; H = 0.50
s(x)
D = d dx
) hs, ψHaar(· y0)i = hs, Dφ(· y0)i = hD∗s, φ(· y0)i Compound Poisson process = piecewise-constant signal Wavelet as a smoothed derivative
sparse innovations (train of Dirac impulses) D∗ = − d
dx (adjoint)
“Sparse derivative” property:
Ds(t) = P
n anδ(x − xn) with xn jump locations
8
DCT→ Karhunen-Lo` eve transform
9
DCT→ Karhunen-Lo` eve transform
10
11
! The real world is continuous
! Input signal ! Imaging physics
! Principled formulation
! Stochastic differential equations (rather than reverse engineering) ! Invariance to coordinate transformations ! Specification of optimal estimators (MAP, MMSE)
! The power of continuous mathematics
! Full backward compatibility with Gaussian theory, link with TV ! Integral operators, characteristic form ! Derivation of joint PDF in any transformed domain
! Operational definition of “sparsity” based on existence considerations:
12
■ Gaussian (Wiener) vs. sparse (Lévy) signals ✔ ■ The spline connection
■
L-splines and signals with finite rate of innovation
■
Green functions as elementary building blocks
■ Sparse stochastic processes
■
Generalized innovation model
■
Gelfand’s theory of generalized stochastic processes
■
Statistical characterization of sparse stochastic processes
■ Implications of innovation model
■
Link with regularization
■
Wavelet representation of sparse processes
■
Determination of transform-domain statistics
■ Sparse processes and signal reconstruction
■
MAP estimator
■
MRI examples
13
Photo courtesy of Carl De Boor
14
(Vetterli et al., 2002)
Location of singularities (knots) : {xn}n=1,...,N Strength of singularities (linear weights): {an}n=1,...,N
Definition The function s(x) is a (non-uniform) L-spline with knots {xn} iff.
N
n=1
15
δ(x)
δ(x) L−1{·} (+ null-space component?)
k∈Z
X
k∈Z
akδ(x − xk)
ρL(x) s(x) = pL(x) + X
k∈Z
akρL(x − xk)
Green function = Impulse response
Translation invariance Linearity
16
δ(x) δ(x − x0) ρ(x)
a[k]δ(x − k) ρ(x − x0)
s(x) =
a[k]ρ(x − k)
d dx = D ⇒ L−1: integrator
17
18
Random jumps with rate λ (Poisson point process)
d dx
k
Jump size distribution:
0.0 0.2 0.4 0.6 0.8 1.0
19
Gaussian white noise
d dx
Non-stationary Self-similar: “1/ω” spectral decay Independent increments = defining property of L´ evy processes
u[k] = s(k) − s(k − 1):
i.i.d. infinitely divisible
0.0 0.2 0.4 0.6 0.8 1.0
w(x) s(x)
White noise (Gaussian, Poisson or Lévy) Generalized stochastic process
Shaping filter
(appropriate boundary conditions)
Whitening operator
What is white noise ?
The problem: Continuous-domain white noise does not have a pointwise interpretation. Standard stochastic calculus. Statisticians circumvent the difficulty by working with ran- dom measures (dW(x) = w(x)dx) and stochastic integrals; i.e, s(x) =
R
R ρ(x, x0)dW(x0)
where ρ(x, x0) is a suitable kernel. Innovation model. The white noise interpretation is more appealing for engineers (cf. Papoulis), but it requires a proper distributional formulation and operator calculus.
21
White noise Characterization of continuous-domain white noise
Mixing operator Whitening operator
s = L−1w w
Characterization of generalized stochastic process Specification of inverse operator Characterization of transform-domain statistics
Multi-scale wavelet analysis
ψi = L∗φi
Functional analysis solution of SDE Very easy ! (after solving 1. & 2.) Easy when: Higher mathematics: generalized functions (Schwartz) measures on topological vector spaces
Gelfand’s theory of generalized stochastic processes Infinite divisibility (Lévy-Khintchine formula)
1 2 4 3
22
White noise Whitening operator
L−1 L
s = L−1w w
Difficulty 2: w is an infinite-dimensional random entity; its “pdf” can be formally specified by a measure Pw(E) where E ⊆ S0 White noise property: independence at every point + stationarity
5 10 15 20Difficulty 1: w 6= w(x) is too rough to have a pointwise interpretation
) X = hw, ϕi for any ϕ 2 S
Infinite-dimensional counterpart of characteristic function (Gelfand, 1955) Characteristic functional:
d Pw(ϕ) = E{ejhw,ϕi} = Z
S0 ejhs,ϕiPw(ds),
for any ϕ ∈ S
5 10 15 20
23
Characteristic functional:
c Ps(ϕ) = G(ϕ) = e 1
2 kϕk2 L2 = exp
✓Z
R
f
◆
Lévy exponent:
2ω2
Discrete white Gaussian noise
X = (X1, . . . , XN) with Xn i.i.d standardized Gaussian
Continuous-domain white Gaussian noise
Infinite-dimensional entity w with generic observations Xn = hw, ϕni
Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê Ê 5 10 15 20
ˆ pXn(ω) = E{ejωhw,ϕni} = E{ejhw,ωϕni} = c Ps(ωϕn) = e 1
2 ω2kϕnk2 L2
Characteristic function:
ˆ pX(ω) = g(ω) = exp PN
n=1 f(ωn)
2 kωk2
24
(based on Schoenberg’s correspondence theorem + Minlos-Bochner theorem)
Functional characterization (Gelfand-Vilenkin)
The characteristic form
d Pw(ϕ) = E{ejhw,ϕi} = exp ✓Z
Rd f
◆
defines a white noise w over S0(Rd)
⇔ f(ω) is first-order conditionally-positive definite
evy exponent)
25
Definition
Example:
f(ω) = −|ω|α, 0 < α ≤ 2
Schoenberg’s correspondence theorem The function eτf(ω) is positive-definite for any τ > 0 if and only if f(ω) is conditionally positive-definite of order one ; i.e.,
N
X
m=1 N
X
n=1
f(ωm − ωn)ξmξn ≥ 0
under the condition PN
m=1 ξm = 0 for every possible choice of ω1, . . . , ωN ∈ R, ξ1, . . . , ξN ∈
C and N ∈ Z+. ⇔ ˆ pX(ω) = ef(ω) is the characteristic function of an infinitely divisible random variable
A continuous, complex-valued function f : R → C such that f(0) = 0 is a valid Lévy exponent if and only if ˆ
pXτ (ω) = eτf(ω) is a valid characteristic function for any τ > 0.
Continuous-domain white noise is highly singular; its points values are undefined A given brand uniquely specified by pid(x) = F−1{ef(ω)}(x)
d Pw
⇥ pid(x) = pXid(x)
with Xid = ⇤rect(· k), w⌅ (i.i.d.)
f(ω) = − 1
2|ω|2
⇔ pid(x): normalized Gaussian f(ω) = −|ω|α with α ∈ (0, 2] ⇔ pid(x): Symmetric-α-stable (SαS)
Also allowed: compound Poisson, Beta, Student, Cauchy, etc. (typically heavy tailed)
27
4 2 2 4 0.1 0.2 0.3 0.4 4 2 2 4 0.1 0.2 0.3 0.4 0.5 4 2 2 4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 4 2 2 4 0.05 0.10 0.15 0.20 0.25 0.30(a) Gaussian (b) Laplace (c) Compound Poisson (d) Cauchy (stable)
2000 5 5 2000 5 5 2000 5 5 2000 5 5
Sparser
f(ω) = λ R
R(ejxω − 1)p(x)dx
f(ω) = log ⇣
1 1+ω2
⌘ pid(x)
Complete mathematical characterization:
d Pw(ϕ) = exp ✓Z
Rd f
◆ f(ω) = −
1 2σ2
0 |ω|2
f(ω) = −s|ω|
28
Definition: A random variable X with generic pdf pid(x) is infinitely divisible (id) iff., for any N ∈ Z+, there exist i.i.d. random variables X1, . . . , XN such that X has the same distribution as X1 + · · · + XN.
L´ evy-Khinchine theorem
pid(x) is an infinitely divisible (id) PDF iff. its characteristic function can be written as ˆ pid(ω) = Z
R
pid(x)ejωxdx = ef(ω)
with L´ evy exponent
f(ω) = jb1ω − b2ω2 2 + Z
R\{0}
{|a|<1}(a)
where b1 ∈ R and b2 ∈ R+ are some constants, and where V is some (positive) Borel measure such that
R
R min(a2, 1)V (da) < ∞.
29
⇒ L−1wδ is a L-spline with random knots
Impulsive Poisson noise
wδ(x) = X
k∈Z
akδ(x − xk) xk: random point locations in Rd with Poisson density λ ak: i.i.d. random variables with amplitude pdf pA(a)
Theorem [U.-Tafti, IEEE-SP 2011] The characteristic form of impulsive Poisson noise is
c Pwδ(ϕ) = E{ejhwδ,ϕi} = exp ✓Z
Rd fPoisson
◆
with L´ evy exponent
fPoisson(ω) = λ Z
R
(ejaω − 1)pA(a)da = λ(ˆ pA(ω) − 1).
30
White noise Whitening operator
L−1 L
s = L−1w w
s = L−1w
⌅ϕ, s⇧ = ⌅ϕ, L−1w⇧ = ⌅L−1∗ϕ | {z }, w⇧
⇒ c Ps(ϕ) = E{ejhs,ϕi} = d Pw(L1⇤ϕ) = exp ✓Z
Rd f
◆
Mathematical conditions on L−1∗ ?
31
Implication for innovation model
Find an acceptable (Lp stable) inverse operator: T = L−1∗ Theorem [U.-Tafti-Sun, preprint] Let f is a valid L´ evy exponent and T is a linear operator acting on ϕ 2 S(Rd) such that any one of the conditions below is met:
kTϕkLp CkϕkLp for all ϕ 2 S(Rd).
Then, c
Ps(ϕ) = E{ejhs,ϕi} = exp R
R f
generalized stochastic process over S0(Rd).
32
Ds = w s = D−1
0 w
⌅ϕ, s⇧ = ⌅D−1∗ ϕ, w⇧ (by Parseval) L2-stable anti-derivative: D−1∗ ϕ(x) = Z
R
ˆ ϕ(ω) − ˆ ϕ(0) −jω ejωx dω 2π
(Blu-U., IEEE-SP 2007)
Dγs = w
Stabilization ⇔ non-stationary behavior
Characteristic form of fractional Brownian motion
c Ps(ϕ) = exp −1 2 Z
R
ϕ(ω) − ˆ ϕ(0) |ω|γ
2π ! (unstable SDE !)
Characteristic form of Brownian motion (a.k.a. Wiener process)
c Ps(ϕ) = exp ✓ 1 2kD−1∗ ϕk2
L2
◆ = exp 1 2 Z
R
ϕ(ω) ˆ ϕ(0) jω
2π !
33
fBm; H = 0.50 fBm; H = 0.75 fBm; H = 1.25 fBm; H = 1.50 Poisson; H = 0.50 Poisson; H = 0.75 Poisson; H = 1.25 Poisson; H = 1.50
H=.5 H=.75 H=1.25 H=1.5
F
2
Fractional Brownian motion (Mandelbrot, 1968) (U.-Tafti, IEEE-SP 2010)
34
F
35
γ 2
F
36
(Sun-U., 2012) (U.-Tafti, 2011)
L ρ = L1δ L1ϕ L1⇤ϕ Dγ xγ1
+
Γ(γ) Z
R
ejωx 1 (jω)γ ˆ ϕ(ω)dω 2π Z
R
ejωx ˆ ϕ(ω) ˆ ϕ(0) (jω)γ dω 2π ∂γ
τ
|x|γ1 Γ(γ)
γ / ⇥ N Z
R
ejωx 1 (jω)
γ 2 τ(jω) γ 2 +τ ˆ
ϕ(ω)dω 2π Z
R
ejωx ˆ ϕ(ω) ˆ ϕ(0) (jω)
γ 2 τ(jω) γ 2 +τ
dω 2π (∆)
γ 2
Cγ⇤x⇤γd, γ d / ⇥ 2N Z
Rd
ejhx,ωi 1 ⇤ω⇤γ ˆ ϕ(ω) dω (2π)d Z
Rd ejhx,ωi ˆ
ϕ(ω) ˆ ϕ(0) ⇤ω⇤γ dω (2π)d 0 < γ < 1 + d/2
Theorem (Generalized Riesz potentials) Unique left-inverse of L⇤ = (∆)
γ 2 that is Lp-stable and scale-invariant:
Iγ
pϕ(x) =
Z
Rdejhx,ωi ˆ ϕ(ω)P⇥γd+ d
p ⇤ |k|=0
ˆ ϕ(k)(0) ωk
k!
kωkγ dω (2π)d = L1⇤ϕ(x)
⇧ϕ ⌅ S(Rd), ⌃Iγ
pϕ⌃Lp(Rd) < C · ⌃ϕ⌃Lp(Rd)
for γ ⌅ R+\Z+ and 1 ⇥ p ⇥ +⇤.
37
H=.5 H=.75 H=1.25 H=1.75
H+1 2 s(x) = w(x)
(U.-Tafti, IEEE-SP 2010)
38
39
■ Optimized analysis tools = B-splines ■ Decoupling sparse ■ Wavelet analysis ■ Link with regularization ■ Signal reconstruction algorithm (MAP)
40
1 3
d Pw(ϕ) = exp ✓Z
Rd f
◆ c Ps(ϕ) = d Pw(L−1∗ϕ)
White noise Whitening operator
L−1 L
s = L−1w w
2
Generic test function ϕ ∈ S plays the role of index variable
Regularization operator vs. wavelet analysis 4
Solution of SDE
White “noise” signature:
pid(x) = F−1{ef(ω)}(x)
41
Ls = w s = L−1w
Lds(x) = X
k∈Zd
d[k]s(x − k)
βL(x) = LdL−1δ(x) = X
k∈Zd
d[k]ρL(x − k)
Green function ρL(x) such that LρL = δ
⇒ βL should be well-defined
Quality of discrete approximation:
Lds(x) = Ld L−1L | {z }
Id
s(x) = βL ∗ Ls(x)
42
Ds(x) = w(x) s(x) = Z x w(y)dy
Dds(x) = s(x) − s(x − 1)
SDE for Lévy process B-spline of minimal support: β(0)(x) ∈ Lp(R) for p > 0
β(0)(x) = 1+(x) − 1+(x − 1) = rect(x − 1 2)
Green function ρD(x) = 1+(x) (unit step)
43
u = Lds = LdL−1w = βL ∗ w hu, ϕi = hβL ⇤ w, ϕi = hw, β∨
L ⇤ ϕi
with
β∨
L(x) = βL(x)
L ∗ ϕ)
= exp ✓Z
Rd f
L ∗ ϕ(x)
◆ f(ω): Lévy exponent of innovation process
44
u(x) = Lds(x) ⇔ u = Ls
(matrix notation)
u = Lds is stationary and infinitely divisible
Characteristic function of X = hu, ϕi with ϕ = δ:
ˆ pX(ω) = E{ejωX} = d Pu(ωδ) = exp ✓Z
Rd f
◆
Quality of decoupling depends upon support of B-spline βL(x)
E{ejhu,ϕi} = d Pu(ϕ) = exp ✓Z
Rd f
L ∗ ϕ(x)
◆
45
46
φi: smoothing kernel at wavelet scale i vi(x) = ⇥ψi(· x), s⇤ = ⇥L∗φi(· x), L−1w⇤ = ⇥φi(· x), w⇤
Wavelet coefficients vi are stationary with characteristic function d
Pw(ωφi)
Quality of decoupling depends upon support of wavelet/smoothing kernel φi
47
48
1 2 kg Hsk2 2 + σ2 P n ΦX([Ls]n)
pX(x) = F−1{d Pw(ωβ∨
L)}(x)
ΦX(x) = − log pX(x)
u = Ls
(matrix notation)
Ls = w s = L−1w
49
10 20 30 5 10 15 20 25 30
Student potentials: r = 2, 4, 8, 32 (fixed variance)
Sparser
1 2 kg Hsk2 2 + σ2 P n ΦX([Ls]n)
Gaussian: pX(x) =
1 √ 2πσ0 e−x2/(2σ2
0)
⇒ ΦX(x) =
1 2σ2
0 x2
Laplace: pX(x) = λ
2 e−λ|x|
⇒ ΦX(x) = λ|x|
Student: pX(x) =
1 B
2
1 x2 + 1 ◆r+ 1
2
⇒ ΦX(x) =
2
50
FWISTA (Guerquin-Kern TMI 2011), IRWL1 (Candès)
Auxiliary innovation variable: u = Ls
s∗ = argmin
s
1 2 kg Hsk2
2 + σ2 X n
ΦX([u]n) ! s.t. u = Ls AL/ADM (Ramani-Fessler TMI 2011)
Innovation variable: u = Ls
LA(s, u, α) = 1 2 ⇥g Hs⇥2
2+λ
X
n
ΦX([u]n) + αT (Ls u) + µ 2 ⇥Ls u⇥2
2
!
51
αk+1 = αk − µ
Linear inverse problem: Nonlinear denoising:
sk+1 ← arg min
s∈RN LA(s, uk, αk)
uk+1 ← arg min
u∈RN LA(sk+1, u, αk)
sk+1 =
−1 HT y + zk+1
with
zk+1 = LT µuk+1 − α
2 4
1 2 3
Cauchy prior with increasing s
LA(s, u, α) = 1 2 kg Hsk2
2+λ
X
n
ΦX([u]n) + αT (Ls u) + µ 2 kLs uk2
2
!
Sequential minimization Proximal operator taylored to stochastic model
proxΦX(y; λ) = arg min
u
1 2|y − u|2 + λΦX(u)
Original SL Phantom Fourier Sampling Pattern 12 Angles Student prior (log)
Laplace prior (TV)
Original Phantom (Guerquin-Kern TMI 2012) Gaussian prior (Tikhonov) SER =17.69 dB Laplace prior (TV) SER = 21.37 dB Student prior SER = 27.22 dB
53
Real T2 Brain Image MR Angiography Image k-space sampling pattern 40 radial lines
Gaussian Estimator Laplace Estimator Student’s Estimator T2 brain Image 8.71 16.08 15.79 MR Angiography Image 6.31 14.48 14.97
Reconstruction results in dB
54
Astrocytes cells bovine pulmonary artery cells human embryonic stem cells Gaussian Estimator Laplace Estimator Student’s Estimator Astrocytes cells 12.18 10.48 10.52 Pulmonary cells 16.90 19.04 18.34 Stem cells 15.81 20.19 20.50
Deconvolution results in dB
Disk shaped PSF (7x7)
56
! Unifying continuous-domain innovation model
! Backward compatibility with classical Gaussian theory ! Operator-based formulation: Lévy-driven SDEs or SPDEs ! Gaussian vs. sparse (generalized Poisson, student, SαS) ! Focus on unstable SDEs ⇒ non-stationary, self-similar processes
! Regularization
! Central role of B-spline ! Sparsification via “operator-like” behavior
! Theoretical framework for sparse signal recovery
! New statistically-founded sparsity priors ! Analytical determination of PDF in any transformed domain ! Derivation of optimal estimators (MAP, MMSE) ! Guide for the development of novel algorithms
58
Theory of generalized stochastic processes Algorithms and imaging applications
aberlin, K.P . Pruessmann, M. Unser, “A Fast Wavelet-Based Reconstruc- tion Method for Magnetic Resonance Imaging,” IEEE Trans. Medical Imaging, vol. 30, no. 9, pp. 1649-1660, September 2011.
. Tafti, An introduction to sparse stochastic processes, e-book at http://www.sparseprocesses.org.
. Tafti, and Q. Sun, “A unified formulation of Gaussian vs. sparse stochastic pro- cesses—Part I: Continuous-domain theory,” preprint, available at http://arxiv.org/abs/1108.6150.
.D. Tafti, ”Stochastic models for sparse and piecewise-smooth signals”, IEEE Trans. Signal Processing, vol. 59, no. 3, pp. 989-1006, March 2011.
vances in Computational Mathematics, vol. 36, no. 3, pp. 399-441, April 2012. P .D. Tafti, D. Van De Ville, M. Unser, “Invariances, Laplacian-Like Wavelet Bases, and the Whitening
1-59
! Preprints and demos: http://bigwww.epfl.ch/