Robust Uncertainty Principles: Exact Signal ... - Semantic Scholar
1
Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information
` California Institute of Technology Emmanuel Candes,
SIAM Conference on Imaging Science, Salt Lake City, Utah, May 2004
Collaborators: Justin Romberg (Caltech), Terence Tao (UCLA)
2
Incomplete Fourier Information Observe Fourier samples fˆ(ω) on a domain Ω.
22 radial lines, ≈ 8% coverage
3
Classical Reconstruction Backprojection: essentially reconstruct g ∗ with fˆ(ω) ω ∈ Ω g ˆ∗ (ω) = 0 ω 6∈ Ω Original Phantom (Logan−Shepp)
Naive Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
g∗
200
250
4
Interpolation? A Row of the Fourier Matrix 25
20
15
10
5
0
−5
−10
−15
−20
−25
original
0
50
100
150
g∗
200
250
300
5
Total Variation Reconstruction Reconstruct g ∗ with min kgkT V g
s.t. g ˆ(ω) = fˆ(ω), ω ∈ Ω
Original Phantom (Logan−Shepp)
Reconstruction: min BV + nonnegativity constraint
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
250
g ∗ = original — perfect reconstruction!
6
Sparse Spike Train Sparse sequence of NT spikes
Observe NΩ Fourier coefficients 5
4
3
2
1
0
−1
−2
−3
−4
−5
0
20
40
60
80
100
120
140
7
Interpolation?
8
`1 Reconstruction Reconstruct by solving min g
X
|gt |
s.t. g ˆ(ω) = fˆ(ω), ω ∈ Ω
t
For NT ∼ NΩ /2, we recover f perfectly.
original
recovered from 30 Fourier samples
9
Extension to TV kgkT V
=
X
|gi+1 − gi | = `1 -norm of finite differences
i
Given frequency observations on Ω, using min kgkT V
s.t. g ˆ(ω) = fˆ(ω), ω ∈ Ω
we can perfectly reconstruct signals with a small number of jumps.
10
Reconstructed perfectly from 30 Fourier samples
11
Model Problem • Signal made out of T spikes • Observed at only |Ω| frequency locations • Extensions –
Piecewise constant signal
–
Spikes in higher-dimensions; 2D, 3D, etc.
–
Piecewise constant images
–
Many others
12
Sharp Uncertainty Principles • Signal is sparse in time, only |T | spikes • Solve combinatorial optimization problem (P0 )
min kgk`0 := #{t, g(t) 6= 0}, g
g ˆ|Ω = fˆ|Ω
Theorem 1 N (sample size) is prime (i) Assume that |T | ≤ |Ω|/2, then (P0 ) reconstructs exactly. (ii) Assume that |T | > |Ω|/2, then (P0 ) fails at exactly reconstructing f ; ∃f1 , f2 with kf1 k`0 + kf2 k`0 = |Ω| + 1 and fˆ1 (ω) = fˆ2 (ω),
∀ω ∈ Ω
13
`1 Relaxation? Solve convex optimization problem (LP for real-valued signals) X (P1 ) min kgk`1 := |g(t)|, g ˆ|Ω = fˆ|Ω g
t
• Example: Dirac’s comb √ – N equispaced spikes (N perfect square). Invariant through Fourier transform fˆ = f √ – Can find |Ω| = N − N with fˆ(ω) = 0, ∀ω ∈ Ω. –
–
Can’t reconstruct
• More dramatic examples exist • But all these examples are very special
14
Dirac’s Comb f(t)
f(ω) N
N
t
f
ω
fˆ
15
Main Result Theorem 2 Suppose |T | ≤ α(M ) ·
|Ω| log N
Then min-`1 reconstructs exactly with prob. greater than 1 − O(N −M ). (n.b. one can choose α(M ) ∼ [29.6(M + 1)]−1 .
Extensions • |T |, number of jump discontinuities (TV reconstruction) • |T |, number of 2D, 3D spikes. • |T |, number of 2D jump discontinuities (2D TV reconstruction)
16
Heuristics: Robust Uncertainty Principles f unique minimizer of (P1 ) iff X X |f (t) + h(t)| > |f (t)|, t
ˆ |Ω = 0 ∀h, h
t
Triangle inequality X X X X X |f (t) + h(t)| = |f (t) + h(t)| + |ht | ≥ |f (t)| − |h(t)| + |ht | Tc
T
Tc
T
Sufficient condition X T
|h(t)| ≤
X Tc
|h(t)|
⇔
X T
|h(t)| ≤
1 2
khk`1
ˆ |Ω = 0, it is impossible to Conclusion: f unique minimizer if for all h, s.t. h ‘concentrate’ h on T
17
Connections: • Donoho & Stark (88) • Donoho & Huo (01) • Gribonval & Nielsen (03) • Tropp (03) and (04) • Donoho & Elad (03)
18
Dual Viewpoint • Convex problem has a dual • Dual polynomial P (t) =
X
Pˆ (ω)eiωt
ω∈Ω
–
P (t) = sgn(f )(t), ∀t ∈ T
–
|P (t)| < 1, ∀t ∈ T c
–
Pˆ supported on set Ω of visible frequencies
Theorem 3 (i) If FT →Ω and there exits a dual polynomial, then the (P1 ) minimizer (P1 ) is unique and is equal to f . (ii) Conversely, if f is the unique minimizer of (P1 ), then there exists a dual polynomial.
19
Dual Polynomial ^ P(ω)
P(t)
ω
t
Space
Frequency
20
Construction of the Dual Polynomial P (t) =
X ω∈Ω
• P interpolates sgn(f ) on T • P has minimum energy
Pˆ (ω)eiωt
21
Auxilary matrices Hf (t) := −
iω(t−t0 )
X
X
ω∈Ω
t0 ∈E:t0 6=t
e
f (t0 ).
Restriction: • ι∗ is the restriction map, ι∗ f := f |T • ι is the obvious embedding obtained by extending by zero outside of T • Identity ι∗ ι is simply the identity operator on T .
P := (ι −
1 |Ω|
1
∗
H)(ι ι −
|Ω|
ι∗ H)−1 ι∗ sgn(f ).
• Frequency support. P has Fourier transform supported in Ω • Spatial interpolation. P obeys ∗
∗
ι P = (ι ι −
1 |Ω|
∗
∗
ι T )(ι ι −
and so P agrees with sgn(f ) on T .
1 |Ω|
ι∗ T )−1 ι∗ sgn(f ) = ι∗ sgn(f ),
22
Hard Things P := (ι −
• (ι∗ ι −
1 ∗ ι H) |Ω|
1 |Ω|
invertible
• |P (t)| < 1, t ∈ /T
∗
H)(ι ι −
1 |Ω|
ι∗ H)−1 ι∗ sgn(f ).
23
Invertibility
(ι∗ ι−
1 |Ω|
ι∗ H) = IT −
0
Fact: |H0 (t, t )| ∼
p
2
1 |Ω|
H0 ,
0 H0 (t, t0 ) = − P
t = t0 0
iω(t−t ) e . ω∈Ω
|Ω|
kH0 k ≤
∗ Tr(H0 H0 )
=
X
|H0 (t, t0 )|2 ∼ |T |2 · |Ω|
t,t0
Want kH0 k ≤ |Ω|, and therefore 2
2
|T | · |Ω| = O(|Ω| )
⇔
|T | = O(
p
|Ω|)
t 6= t0
24
Key Estimates • Want to show largest eigenvalue of H0 (self-adjoint) is less than Ω. • Take large powers of random matrices 2n
Tr(H0
2n ) = λ2n 1 + . . . + λT
• Key estimate: develop bounds on E[Tr(H02n )] • Key intermediate result: p p kH0 k ≤ γ log |T | |T | |Ω| with large-probability • A lot of combinatorics!
25
Numerical Results • Signal length N = 1024 • Randomly place Nt spikes, observe Nw random frequencies • Measure % recovered perfectly • red = always recovered, blue = never recovered
Nw
Nt /Nw
26
Other Phantoms, I Original Phantom
Classical Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = classical reconstruction
250
27
Original Phantom
Total Variation Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = TV reconstruction = Exact!
250
28
Other Phantoms, II Original Phantom
Classical Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = classical reconstruction
250
29
Original Phantom
Total Variation Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = TV reconstruction = Exact!
250
30
Scanlines A Scanline of the Original Phantom
Classical (Black) and TV (Red) Reconstructions
2
2
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
50
100
150
original
200
250
300
0
0
50
100
150
200
250
g ∗ = classical reconstruction
300
31
Summary • Exact reconstruction • Tied to new uncertainty principles • Stability • Robustness • Optimality • Many extensions: e.g. arbitrary synthesis/measurement pairs Contact: [email protected]
Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information
` California Institute of Technology Emmanuel Candes,
SIAM Conference on Imaging Science, Salt Lake City, Utah, May 2004
Collaborators: Justin Romberg (Caltech), Terence Tao (UCLA)
2
Incomplete Fourier Information Observe Fourier samples fˆ(ω) on a domain Ω.
22 radial lines, ≈ 8% coverage
3
Classical Reconstruction Backprojection: essentially reconstruct g ∗ with fˆ(ω) ω ∈ Ω g ˆ∗ (ω) = 0 ω 6∈ Ω Original Phantom (Logan−Shepp)
Naive Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
g∗
200
250
4
Interpolation? A Row of the Fourier Matrix 25
20
15
10
5
0
−5
−10
−15
−20
−25
original
0
50
100
150
g∗
200
250
300
5
Total Variation Reconstruction Reconstruct g ∗ with min kgkT V g
s.t. g ˆ(ω) = fˆ(ω), ω ∈ Ω
Original Phantom (Logan−Shepp)
Reconstruction: min BV + nonnegativity constraint
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
250
g ∗ = original — perfect reconstruction!
6
Sparse Spike Train Sparse sequence of NT spikes
Observe NΩ Fourier coefficients 5
4
3
2
1
0
−1
−2
−3
−4
−5
0
20
40
60
80
100
120
140
7
Interpolation?
8
`1 Reconstruction Reconstruct by solving min g
X
|gt |
s.t. g ˆ(ω) = fˆ(ω), ω ∈ Ω
t
For NT ∼ NΩ /2, we recover f perfectly.
original
recovered from 30 Fourier samples
9
Extension to TV kgkT V
=
X
|gi+1 − gi | = `1 -norm of finite differences
i
Given frequency observations on Ω, using min kgkT V
s.t. g ˆ(ω) = fˆ(ω), ω ∈ Ω
we can perfectly reconstruct signals with a small number of jumps.
10
Reconstructed perfectly from 30 Fourier samples
11
Model Problem • Signal made out of T spikes • Observed at only |Ω| frequency locations • Extensions –
Piecewise constant signal
–
Spikes in higher-dimensions; 2D, 3D, etc.
–
Piecewise constant images
–
Many others
12
Sharp Uncertainty Principles • Signal is sparse in time, only |T | spikes • Solve combinatorial optimization problem (P0 )
min kgk`0 := #{t, g(t) 6= 0}, g
g ˆ|Ω = fˆ|Ω
Theorem 1 N (sample size) is prime (i) Assume that |T | ≤ |Ω|/2, then (P0 ) reconstructs exactly. (ii) Assume that |T | > |Ω|/2, then (P0 ) fails at exactly reconstructing f ; ∃f1 , f2 with kf1 k`0 + kf2 k`0 = |Ω| + 1 and fˆ1 (ω) = fˆ2 (ω),
∀ω ∈ Ω
13
`1 Relaxation? Solve convex optimization problem (LP for real-valued signals) X (P1 ) min kgk`1 := |g(t)|, g ˆ|Ω = fˆ|Ω g
t
• Example: Dirac’s comb √ – N equispaced spikes (N perfect square). Invariant through Fourier transform fˆ = f √ – Can find |Ω| = N − N with fˆ(ω) = 0, ∀ω ∈ Ω. –
–
Can’t reconstruct
• More dramatic examples exist • But all these examples are very special
14
Dirac’s Comb f(t)
f(ω) N
N
t
f
ω
fˆ
15
Main Result Theorem 2 Suppose |T | ≤ α(M ) ·
|Ω| log N
Then min-`1 reconstructs exactly with prob. greater than 1 − O(N −M ). (n.b. one can choose α(M ) ∼ [29.6(M + 1)]−1 .
Extensions • |T |, number of jump discontinuities (TV reconstruction) • |T |, number of 2D, 3D spikes. • |T |, number of 2D jump discontinuities (2D TV reconstruction)
16
Heuristics: Robust Uncertainty Principles f unique minimizer of (P1 ) iff X X |f (t) + h(t)| > |f (t)|, t
ˆ |Ω = 0 ∀h, h
t
Triangle inequality X X X X X |f (t) + h(t)| = |f (t) + h(t)| + |ht | ≥ |f (t)| − |h(t)| + |ht | Tc
T
Tc
T
Sufficient condition X T
|h(t)| ≤
X Tc
|h(t)|
⇔
X T
|h(t)| ≤
1 2
khk`1
ˆ |Ω = 0, it is impossible to Conclusion: f unique minimizer if for all h, s.t. h ‘concentrate’ h on T
17
Connections: • Donoho & Stark (88) • Donoho & Huo (01) • Gribonval & Nielsen (03) • Tropp (03) and (04) • Donoho & Elad (03)
18
Dual Viewpoint • Convex problem has a dual • Dual polynomial P (t) =
X
Pˆ (ω)eiωt
ω∈Ω
–
P (t) = sgn(f )(t), ∀t ∈ T
–
|P (t)| < 1, ∀t ∈ T c
–
Pˆ supported on set Ω of visible frequencies
Theorem 3 (i) If FT →Ω and there exits a dual polynomial, then the (P1 ) minimizer (P1 ) is unique and is equal to f . (ii) Conversely, if f is the unique minimizer of (P1 ), then there exists a dual polynomial.
19
Dual Polynomial ^ P(ω)
P(t)
ω
t
Space
Frequency
20
Construction of the Dual Polynomial P (t) =
X ω∈Ω
• P interpolates sgn(f ) on T • P has minimum energy
Pˆ (ω)eiωt
21
Auxilary matrices Hf (t) := −
iω(t−t0 )
X
X
ω∈Ω
t0 ∈E:t0 6=t
e
f (t0 ).
Restriction: • ι∗ is the restriction map, ι∗ f := f |T • ι is the obvious embedding obtained by extending by zero outside of T • Identity ι∗ ι is simply the identity operator on T .
P := (ι −
1 |Ω|
1
∗
H)(ι ι −
|Ω|
ι∗ H)−1 ι∗ sgn(f ).
• Frequency support. P has Fourier transform supported in Ω • Spatial interpolation. P obeys ∗
∗
ι P = (ι ι −
1 |Ω|
∗
∗
ι T )(ι ι −
and so P agrees with sgn(f ) on T .
1 |Ω|
ι∗ T )−1 ι∗ sgn(f ) = ι∗ sgn(f ),
22
Hard Things P := (ι −
• (ι∗ ι −
1 ∗ ι H) |Ω|
1 |Ω|
invertible
• |P (t)| < 1, t ∈ /T
∗
H)(ι ι −
1 |Ω|
ι∗ H)−1 ι∗ sgn(f ).
23
Invertibility
(ι∗ ι−
1 |Ω|
ι∗ H) = IT −
0
Fact: |H0 (t, t )| ∼
p
2
1 |Ω|
H0 ,
0 H0 (t, t0 ) = − P
t = t0 0
iω(t−t ) e . ω∈Ω
|Ω|
kH0 k ≤
∗ Tr(H0 H0 )
=
X
|H0 (t, t0 )|2 ∼ |T |2 · |Ω|
t,t0
Want kH0 k ≤ |Ω|, and therefore 2
2
|T | · |Ω| = O(|Ω| )
⇔
|T | = O(
p
|Ω|)
t 6= t0
24
Key Estimates • Want to show largest eigenvalue of H0 (self-adjoint) is less than Ω. • Take large powers of random matrices 2n
Tr(H0
2n ) = λ2n 1 + . . . + λT
• Key estimate: develop bounds on E[Tr(H02n )] • Key intermediate result: p p kH0 k ≤ γ log |T | |T | |Ω| with large-probability • A lot of combinatorics!
25
Numerical Results • Signal length N = 1024 • Randomly place Nt spikes, observe Nw random frequencies • Measure % recovered perfectly • red = always recovered, blue = never recovered
Nw
Nt /Nw
26
Other Phantoms, I Original Phantom
Classical Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = classical reconstruction
250
27
Original Phantom
Total Variation Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = TV reconstruction = Exact!
250
28
Other Phantoms, II Original Phantom
Classical Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = classical reconstruction
250
29
Original Phantom
Total Variation Reconstruction
50
50
100
100
150
150
200
200
250
250 50
100
150
original
200
250
50
100
150
200
g ∗ = TV reconstruction = Exact!
250
30
Scanlines A Scanline of the Original Phantom
Classical (Black) and TV (Red) Reconstructions
2
2
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
50
100
150
original
200
250
300
0
0
50
100
150
200
250
g ∗ = classical reconstruction
300
31
Summary • Exact reconstruction • Tied to new uncertainty principles • Stability • Robustness • Optimality • Many extensions: e.g. arbitrary synthesis/measurement pairs Contact: [email protected]
