Beyond Compressed Sensing: The Effectiveness of Convex Programming - - PowerPoint PPT Presentation

Beyond Compressed Sensing: The Effectiveness of Convex Programming in the Information and Physical Sciences

Emmanuel Cand` es EUSIPCO 2015, Nice, September 2015

Three stories about signal recovery from missing information

Today I want to tell you three stories from my life. That’s it. No big deal. Just three stories Steve Jobs

Three stories about signal recovery from missing information

Today I want to tell you three stories from my life. That’s it. No big deal. Just three stories Steve Jobs (1) Missing phase (phase retrieval) (2) Missing and/or corrupted entries in data matrix (robust PCA) (3) Missing high-frequency spectrum (super-resolution) Makes signal/data recovery difficult

This lecture

Convex programming usually (but not always) returns the right answer!

Story # 1: Phase Retrieval

Collaborators: Y. Eldar, X. Li, T. Strohmer, V. Voroninski

X-ray crystallography

Method for determining atomic structure within a crystal principle typical setup 10 Nobel Prizes in X-ray crystallography, and counting...

Importance

principle Franklin’s photograph

Missing phase problem

Detectors record intensities of diffracted rays → phaseless data only! Fraunhofer diffraction − → intensity of electrical field |ˆ x(f1, f2)|2 =

• x(t1, t2)e−i2π(f1t1+f2t2) dt1dt2
• 2

Missing phase problem

Detectors record intensities of diffracted rays → phaseless data only! Fraunhofer diffraction − → intensity of electrical field |ˆ x(f1, f2)|2 =

• x(t1, t2)e−i2π(f1t1+f2t2) dt1dt2
• 2

Phase retrieval problem (inversion)

How can we recover the phase (or equivalently signal x(t1, t2)) from |ˆ x(f1, f2)|?

About the importance of phase...

About the importance of phase...

DFT DFT

About the importance of phase...

DFT DFT keep magnitude swap phase

About the importance of phase...

DFT DFT keep magnitude swap phase

X-ray imaging: now and then

R¨

• ntgen (1895)

Dierolf (2010)

Ultrashort X-ray pulses

Imaging single large protein complexes

Discrete mathematical model

Phaseless measurements about x0 ∈ Cn bk = |ak, x0|2 k ∈ {1, . . . , m} = [m] Phase retrieval is feasibility problem find x subject to |ak, x|2 = bk k ∈ [m] Solving quadratic equations is NP hard in general

Discrete mathematical model

Phaseless measurements about x0 ∈ Cn bk = |ak, x0|2 k ∈ {1, . . . , m} = [m] Phase retrieval is feasibility problem find x subject to |ak, x|2 = bk k ∈ [m] Solving quadratic equations is NP hard in general 1985 Nobel Prize for Hauptman & Karle: use very specific prior knowledge

Discrete mathematical model

Phaseless measurements about x0 ∈ Cn bk = |ak, x0|2 k ∈ {1, . . . , m} = [m] Phase retrieval is feasibility problem find x subject to |ak, x|2 = bk k ∈ [m] Solving quadratic equations is NP hard in general 1985 Nobel Prize for Hauptman & Karle: use very specific prior knowledge Standard approach: Gerchberg–Saxton (or Fienup) iterative algorithm Sometimes works well Sometimes does not

Quadratic equations: geometric view I

Quadratic equations: geometric view I

Quadratic equations: geometric view I

Quadratic equations: geometric view I

Quadratic equations: geometric view I

Quadratic equations: geometric view I

Quadratic equations: geometric view II

PhaseLift

Lifting: X = xx∗ bk = |ak, x|2 = a∗

kxx∗ak = a∗ kXak

Turns quadratic measurements into linear measurements b := A(X) about xx∗

PhaseLift

Lifting: X = xx∗ bk = |ak, x|2 = a∗

kxx∗ak = a∗ kXak

Turns quadratic measurements into linear measurements b := A(X) about xx∗

Phase retrieval: equivalent formulation

find X

• s. t.

A(X) = b X 0, rank(X) = 1 ⇐ ⇒ min rank(X)

• s. t.

A(X) = b X 0 Combinatorially hard

PhaseLift

Lifting: X = xx∗ bk = |ak, x|2 = a∗

kxx∗ak = a∗ kXak

Turns quadratic measurements into linear measurements b := A(X) about xx∗

PhaseLift: tractable semidefinite relaxation

minimize Tr(X) subject to A(X) = b X 0 This is a semidefinite program (SDP) Trace is convex proxy for rank

Semidefinite programming (SDP)

Special class of convex optimization problems Relatively natural extension of linear programming (LP) ‘Efficient’ numerical solvers (interior point methods)

LP (std. form): x ∈ Rn

minimize c, x subject to aT

k x = bk k = 1, . . .

x ≥ 0

SDP (std. form): X ∈ Rn×n

minimize C, X subject to Ak, X = bk k = 1, . . . X 0 Standard inner product: C, X = Tr(C∗X)

This is not really new...

Relaxation of quadratically constrained QP’s Shor (87) [Lower bounds on nonconvex quadratic optimization problems] Goemans and Williamson (95) [MAX-CUT] Ben-Tal and Nemirovskii (01) [Monograph] ... Similar approach for array imaging: Chai, Moscoso, Papanicolaou (11)

From overdetermined to highly underdetermined

Quadratic equations b = A(xx∗) bk = |ak, x|2 k ∈ [m] Lift b = A(X) minimize Tr(X) subject to A(X) = b X 0

Have we made things worse?

• verdetermined (m > n)

→ highly underdetermined (m ≪ n2)

Exact phase retrieval via SDP

Quadratic equations b = A(xx∗) bk = |ak, x|2 k ∈ [m] Simplest model: ak independently and uniformly sampled on unit sphere

• f Cn if x ∈ Cn (complex-valued problem)
• f Rn if x ∈ Rn (real-valued problem)

Exact phase retrieval via SDP

Quadratic equations b = A(xx∗) bk = |ak, x|2 k ∈ [m] Simplest model: ak independently and uniformly sampled on unit sphere

• f Cn if x ∈ Cn (complex-valued problem)
• f Rn if x ∈ Rn (real-valued problem)

Theorem (C. and Li (’12); C., Strohmer and Voroninski (’11))

Assume m n. With prob. 1 − O(e−γm), for all x ∈ Cn, only point in feasible set {X : A(X) = b and X 0} is xx∗

Exact phase retrieval via SDP

Quadratic equations b = A(xx∗) bk = |ak, x|2 k ∈ [m] Simplest model: ak independently and uniformly sampled on unit sphere

• f Cn if x ∈ Cn (complex-valued problem)
• f Rn if x ∈ Rn (real-valued problem)

Theorem (C. and Li (’12); C., Strohmer and Voroninski (’11))

Assume m n. With prob. 1 − O(e−γm), for all x ∈ Cn, only point in feasible set {X : A(X) = b and X 0} is xx∗ Injectivity if m ≥ 4n − 2 (Balan, Bodmann, Casazza, Edidin ’09)

Extensions to physical setups

Random masking + diffraction Similar theory: C. , Li and Soltanolkotabi (’13)

Numerical results: noiseless recovery

(a) Smooth signal (real part) (b) Random signal (real part)

Figure: Recovery (with reweighting) of n-dimensional complex signal (2n unknowns) from 4n quadratic measurements (random binary masks)

How is this possible?

How can feasible set {X 0} ∩ {A(X) = b} have a unique point? Intersection of

• x

y y z

• 0 with affine space

Correct representation

Rank-1 matrices are on the boundary (extreme rays) of PSD cone

My mental representation

My mental representation

My mental representation

My mental representation

My mental representation

Numerical results: noiseless 2D images

50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250

• riginal image

3 Gaussian masks

50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250

8 binary masks error with 8 binary masks

Courtesy

• S. Marchesini (LBL)

Numerical results: noisy 2D data

Low SNR (20dB) High SNR (60 dB)

Reconstructions from noisy data using 32 Gaussian random masks Similar output with 16 masks

With noise

Figure: Relative MSE vs. signal-to-noise ratio on dB scale (binary masks)

Story #2: Robust Principal Component Analysis

Collaborators: X. Li, Y. Ma, J. Wright

The separation problem (Chandrasekahran et al.)

M = L + S M: data matrix (observed) L: low-rank (unobserved) S: sparse (unobserved)

The separation problem (Chandrasekahran et al.)

M = L + S M: data matrix (observed) L: low-rank (unobserved) S: sparse (unobserved)

Problem: can we recover L and S accurately?

Again, missing information

Motivation: robust principal component analysis (RPCA)

PCA sensitive to outliers: breaks down with one (badly) corrupted data point      x11 x12 . . . x1n x21 x22 . . . x2n . . . . . . . . . . . . xd1 xd2 . . . xdn      ❆

Motivation: robust principal component analysis (RPCA)

PCA sensitive to outliers: breaks down with one (badly) corrupted data point      x11 x12 . . . x1n x21 x22 . . . x2n . . . . . . . . . . . . xd1 xd2 . . . xdn      = ⇒      x11 x12 . . . x1n x21 x22 . . . x2n . . . . . . . . . . . . xd1 ❆ . . . xdn     

Robust PCA

Data increasingly high dimensional Gross errors frequently occur in many applications

Image processing Web data analysis Bioinformatics ... Occlusions Malicious tampering Sensor failures ...

     ❆ x12 . . . x1n x21 x22 . . . ❆ . . . . . . . . . . . . xd1 ❆ . . . xdn      Important to make PCA robust

Gross errors

Movies Users         × × ❆ ❆ × ❆ × × × ❆ ×         Observe corrupted entries Yij = Lij + Sij (i, j) ∈ Ωobs L low-rank matrix S entries that have been tampered with (impulsive noise)

Problem

Recover L from missing and corrupted samples

When does separation make sense?

M = L + S Sparse component cannot be low rank: S =      ∗ · · · ∗ · · · . . . . . . . . . . . . . . . . . . ∗ · · ·      Sparsity pattern will be assumed (uniform) random

When does separation make sense?

M = L + S Sparse component cannot be low rank: S =      ∗ · · · ∗ · · · . . . . . . . . . . . . . . . . . . ∗ · · ·      Sparsity pattern will be assumed (uniform) random Low-rank component cannot be sparse: L =      ∗ ∗ ∗ ∗ · · · ∗ ∗ · · · . . . . . . . . . . . . . . . . . . · · ·     

Sparse component cannot be low-rank

L =      x1 x2 · · · xn−1 xn x1 x2 · · · xn−1 xn . . . . . . . . . . . . . . . x1 x2 · · · xn−1 xn     

• 1x∗

⇒ L + S =      ❆ x2 · · · xn−1 xn ❆ x2 · · · xn−1 xn . . . . . . . . . . . . . . . ❆ x2 · · · xn−1 xn     

Low-rank component cannot be sparse

L =          x1 x2 x3 x4 · · · xn−1 xn x1 x2 x3 x4 · · · xn−1 xn · · · · · · . . . . . . . . . . . . . . . . . . · · ·          Incoherent condition [C. and Recht (’08)]: column and row spaces not aligned with coordinate axes (cannot have small subsets of rows and/or columns that are singular)

Low-rank component cannot be sparse

M =           x1 x2 ❆ x4 · · · xn−1 xn x1 x2 ❆ x4 · · · xn−1 ❆ ❆ ❆ · · · ❆ · · · . . . . . . . . . . . . . . . . . . ❆ · · · ❆           Incoherent condition [C. and Recht (’08)]: column and row spaces not aligned with coordinate axes (cannot have small subsets of rows and/or columns that are singular)

Demixing by convex programming

M = L + S L unknown (rank unknown) S unknown (# of entries = 0, locations, magnitudes all unknown)

Demixing by convex programming

M = L + S L unknown (rank unknown) S unknown (# of entries = 0, locations, magnitudes all unknown)

Recovery via SDP

minimize ˆ L∗ + λ ˆ S1 subject to ˆ L + ˆ S = M See also Chandrasekaran, Sanghavi, Parrilo, Willsky (’09) nuclear norm: L∗ =

i σi(L) (sum of sing. values)

ℓ1 norm: S1 =

ij |Sij| (sum of abs. values)

Exact recovery via SDP

min ˆ L∗ + λ ˆ S1

• s. t.

ˆ L + ˆ S = M ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆

Exact recovery via SDP

min ˆ L∗ + λ ˆ S1

• s. t.

ˆ L + ˆ S = M

Theorem

L is n × n of rank(L) ≤ ρrn (log n)−2 and incoherent S is n × n, random sparsity pattern of cardinality at most ρsn2 Then with probability 1 − O(n−10), SDP with λ = 1/√n is exact: ˆ L = L, ˆ S = S Same conclusion for rectangular matrices with λ = 1/ √ max dim ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆ ❆

Exact recovery via SDP

min ˆ L∗ + λ ˆ S1

• s. t.

ˆ L + ˆ S = M

Theorem

L is n × n of rank(L) ≤ ρrn (log n)−2 and incoherent S is n × n, random sparsity pattern of cardinality at most ρsn2 Then with probability 1 − O(n−10), SDP with λ = 1/√n is exact: ˆ L = L, ˆ S = S Same conclusion for rectangular matrices with λ = 1/ √ max dim No tuning parameter! Whatever the magnitudes of L and S          × ❆ ❆ ❆ × ❆ ❆ ❆ × × ❆ ❆ × ❆ ❆ × ❆ ❆ ❆ ❆ × ❆ ❆ × × ❆ ❆ ❆ ❆ ❆ ❆ ❆ × × ❆ ❆         

Phase transitions in probability of success

(a) RPCA, Random Signs (b) RPCA, Coherent Signs (c) Matrix Completion

L = XY T is a product of independent n × r i.i.d. N(0, 1/n) matrices

Missing and corrupted

RPCA

min ˆ L∗ + λ ˆ S1

• s. t.

ˆ Lij + ˆ Sij = Lij + Sij (i, j) ∈ Ωobs

        × ❆ ? ? × ? ? ? × ❆ ? ? × ? ? × ? ? ? ? × ? ? ❆ × ? ❆ ? ? ? ? ? × ❆ ? ?        

Same theorem: with high prob. λ = 1 √

• frac. observed × max dim

= ⇒ ˆ L = L!

Video surveillance

Sequence of 200 video frames (144 × 172 pixels) with a static background Problem: detect any activity in the foreground

… …

RPCA

L + S background subtraction

L + S background subtraction

From GoDec

L + S reconstruction of MR angiography

L + S L S automatic and improved background suppression

Joint with R. Otazo and D. Sodickson

Free-breathing MRI of the liver

NUFFT Standard L + S Motion-Guided L + S 12.8 fold acceleration min L∗ + λS1

• s. t.

A(L + S) = y

Joint with R. Otazo and D. Sodickson

Free-breathing MRI of the liver

NUFFT Standard L + S Motion-Guided L + S Temporal blurring

Joint with R. Otazo and D. Sodickson

Free-breathing MRI of the kidneys

NUFFT Standard L + S Motion-Guided L + S 12.8 fold acceleration min L∗ + λS1

• s. t.

A(L + S) = y

Joint with R. Otazo and D. Sodickson

Free-breathing MRI of the kidneys

NUFFT Standard L + S Motion-Guided L + S

Joint with R. Otazo and D. Sodickson

Story #3: Super-resolution

Collaborator: C. Fernandez-Granda

Limits of resolution

In any optical imaging system, diffraction imposes fundamental limit on resolution

The physical phenomenon called diffraction is of the utmost importance in the theory of optical imaging systems (Joseph Goodman)

Interested in usual bandlimited imaging systems

Pupil Airy disk Cross section

Rayleigh resolution limit

Lord Rayleigh

The super-resolution problem

⇐

• bjective

data Retrieve fine scale information from low-pass data

⇐

Equivalent description: extrapolate spectrum Fundamental problem Radar Microscopy Spectroscopy Medical imaging Astronomy Geophysics ...

Single molecule imaging

Microscope receives light from fluorescent molecules

Problem

Resolution is much coarser than size of individual molecules (low-pass data) Can we ‘beat’ the diffraction limit and super-resolve those molecules? Higher molecule density − → faster imaging

Mathematical model

Signal x =

j ajδτj

aj ∈ C, τj ∈ [0, 1] Data y = Fnx: n = 2flo + 1 low-frequency coefficients (Nyquist sampling) y(k) = 1 e−i2πktx(dt) =

• j

aje−i2πkτj k ∈ Z, |k| ≤ flo Resolution limit: (λlo/2 is Rayleigh distance) 1/flo = λlo

Mathematical model

Signal x =

j ajδτj

aj ∈ C, τj ∈ [0, 1] Data y = Fnx: n = 2flo + 1 low-frequency coefficients (Nyquist sampling) y(k) = 1 e−i2πktx(dt) =

• j

aje−i2πkτj k ∈ Z, |k| ≤ flo Resolution limit: (λlo/2 is Rayleigh distance) 1/flo = λlo

Question

Can we resolve the signal beyond this limit?

Can you find the spikes?

Low-frequency data about spike train

Can you find the spikes?

Low-frequency data about spike train

Recovery by minimum total variation

Recovery by cvx prog.

min ˜ xTV subject to Fn ˜ x = y xTV =

• |x(dt)| is continuous analog of discrete ℓ1 norm xℓ1 =

t |xt|

x =

• j

ajδτj = ⇒ xTV =

• j

|aj| Work on ℓ1 minimization: Logan, Donoho, Stark, Tropp, Elad, C. , Tao...

Recovery by convex programming

y(k) = 1 e−i2πktx(dt) |k| ≤ flo

Recovery by convex programming

y(k) = 1 e−i2πktx(dt) |k| ≤ flo

Theorem (C. and Fernandez Granda (2012))

If spikes are separated by at least 1.86 /flo := 1.86 λlo then min TV solution is exact! (Current state of the art 1.4λlo)

Recovery by convex programming

y(k) = 1 e−i2πktx(dt) |k| ≤ flo

Theorem (C. and Fernandez Granda (2012))

If spikes are separated by at least 1.86 /flo := 1.86 λlo then min TV solution is exact! (Current state of the art 1.4λlo) Infinite precision! (Whatever the amplitudes)

Recovery by convex programming

y(k) = 1 e−i2πktx(dt) |k| ≤ flo

Theorem (C. and Fernandez Granda (2012))

If spikes are separated by at least 1.86 /flo := 1.86 λlo then min TV solution is exact! (Current state of the art 1.4λlo) Infinite precision! (Whatever the amplitudes) Cannot go below λlo

Recovery by convex programming

y(k) = 1 e−i2πktx(dt) |k| ≤ flo

Theorem (C. and Fernandez Granda (2012))

If spikes are separated by at least 1.86 /flo := 1.86 λlo then min TV solution is exact! (Current state of the art 1.4λlo) Infinite precision! (Whatever the amplitudes) Cannot go below λlo Can recover (2λlo)−1 = flo/2 = n/4 spikes from n low-freq. samples

Recovery by convex programming

y(k) = 1 e−i2πktx(dt) |k| ≤ flo

Theorem (C. and Fernandez Granda (2012))

If spikes are separated by at least 1.86 /flo := 1.86 λlo then min TV solution is exact! (Current state of the art 1.4λlo) Infinite precision! (Whatever the amplitudes) Cannot go below λlo Can recover (2λlo)−1 = flo/2 = n/4 spikes from n low-freq. samples Essentially same result in higher dimensions

Formulation as a finite-dimensional problem

Primal problem min xTV s. t. Fnx = y Infinite-dimensional variable x Finitely many constraints Dual problem max Rey, c s. t. F∗

nc∞ ≤ 1

Finite-dimensional variable c Infinitely many constraints (F∗

n c)(t) =

• |k|≤flo

ckei2πkt

Formulation as a finite-dimensional problem

Primal problem min xTV s. t. Fnx = y Infinite-dimensional variable x Finitely many constraints Dual problem max Rey, c s. t. F∗

nc∞ ≤ 1

Finite-dimensional variable c Infinitely many constraints (F∗

n c)(t) =

• |k|≤flo

ckei2πkt

Semidefinite representability

|(F∗

n c)(t)| ≤ 1 for all t ∈ [0, 1] equivalent to

(1) there is Q Hermitian s. t.

• Q

c c∗ 1

• (2) Tr(Q) = 1

(3) sums along superdiagonals vanish: n−j

i=1 Qi,i+j = 0 for 1 ≤ j ≤ n − 1

Dual solution

c: coeffs. of low-pass trig. poly.

k ckei2πkt interpolating sign of primal solution

Super-resolution via semidefinite programming

Super-resolution via semidefinite programming

Super-resolution via semidefinite programming

• 1. Solve semidefinite program to obtain dual solution

Super-resolution via semidefinite programming

• 2. Locate points at which corresponding polynomial has unit magnitude

Super-resolution via semidefinite programming

Signal Estimate

• 3. Estimate amplitudes via least squares

Support-location accuracy

fc 25 50 75 100 Average error 6.66 10−9 1.70 10−9 5.58 10−10 2.96 10−10 Maximum error 1.83 10−7 8.14 10−8 2.55 10−8 2.31 10−8 For each fc, 100 random signals with |T| = fc/4 and ∆(T) ≥ 2/fc

Example

Minimum separation : 1.5 λc

Example

SNR 20 dB

Noisy Noiseless

Example

SNR 20 dB

Signal Estimate

Example

SNR 15 dB

Noisy Noiseless

Example

SNR 15 dB

Signal Estimate

Example

SNR 5 dB

Noisy Noiseless

Example

SNR 5 dB

Signal Estimate

Noise is amenable to convex modeling

Phase retrieval minimize b − A(X)1 subject to X 0 Robust PCA minimize L∗ + λS1 + γY − L − S2

F

Super-resolution minimize xTV subject to y − Fnx2

2 ≤ σ2

Stability in all cases, sometimes near the information theoretic limit

Apologies: things I have not talked

Noise Algorithms Applications Avalanche of related works

A small sample of papers I have greatly enjoyed

Phase retrieval

Netrapalli, Jain, Sanghavi, Phase retrieval using alternating minimization (’13) Waldspurger, d’Aspremont, Mallat, Phase recovery, MaxCut and complex semidefinite programming (’12)

Robust PCA

Gross, Recovering low-rank matrices from few coefficients in any basis (’09) Chandrasekaran, Parrilo and Willsky, Latent variable graphical model selection via convex optimization (’11) Hsu, Kakade and Zhang, Robust matrix decomposition with outliers (’11)

Super-resolution

Kahane, Analyse et synth` ese harmoniques (’11) Slepian, Prolate spheroidal wave functions, Fourier analysis, and uncertainty. V - The discrete case (’78)

Summary

Three important problems with missing data

Phase retrieval Matrix completion/RPCA Super-resolution

Three simple and model-free recovery procedures via convex programming Three near-perfect solutions — Three theorems Three fundamental modeling/computational tools Matrices Probability (Convex) Optimization

Numerical results: noiseless recovery

(a) Smooth signal (real part) (b) Random signal (real part)

Figure: Recovery (with reweighting) of n-dimensional complex signal (2n unknowns) from 4n quadratic measurements (random binary masks)

About separation: sparsity is not enough!

CS: sparse signals are ‘away’ from null space of sampling operator Super-res: this is not the case Signal Spectrum

10

−10

10

−5

10

About separation: sparsity is not enough!

CS: sparse signals are ‘away’ from null space of sampling operator Super-res: this is not the case Signal Spectrum

10

−10

10

−5

10 10

−10

10

−5

10

Analysis via prolate spheroidal functions

David Slepian If distance between spikes less than λlo/2 (Rayleigh), problem hopelessly ill posed