Rémi Gribonval INRIA Rennes - Bretagne Atlantique, France · Inverse problems and sparse models...

Inverse problems and sparse models (4/5)

Rémi Gribonval INRIA Rennes - Bretagne Atlantique, France

[email protected]

2012R. GRIBONVAL - SPARSE METHODS - SISEA

Exercice: Matlab code for MP

• Full clean code would include some checking (column normalization, dimension checking, etc.)

function [x res] = mp(b,A,k)% compute k-sparse approximation to b with matrix A using Matching pursuit[m,N] = size(A); x = zeros(N,1); res = b; for i=1:k

% compute correlation between residual and columns of Acorr = A’*res; % find position n (and value c) of the maximally correlated column [c n] = max(abs(corr)); % NB: modern Matlab notation allows [~ n] = max(abs(corr)) % update the representation x(n) = x(n) + corr(n); % update the residualres = res - corr(n)*A(:,n)

end

102


Overview of greedy algorithms

103

Matching Pursuit OMP Stagewise OMPSelection

Update

MP & OMP: Mallat & Zhang 1993StOMP: Donoho & al 2006 (similar to MCA, Bobin & al 2006)

A = [A1, . . .AN ]

�i := arg maxn

|ATnri�1| �i := {n | |AT

nri�1| > �i}

⇥i = ⇥i�1 � �i

xi = xi�1 + A+�i

ri�1

⇥i = ⇥i�1 � �i

xi = A+�i

bri = b�A�ixi

b = Axi + ri

ri = ri�1 �A�iA+�i

ri�1


Pseudo-inversion

• Over-determined case: m>k • Under-determined case m<k

104

B

B+B+

Bm

k

m

B+ = (BTB)�1BT B+ = BT (BBT )�1

B+B = Idk BB+ = Idm


Exercice: Matlab code for OMP

• Full clean code would include some checkingfunction [x res] = omp(b,A,k)% compute k-sparse approximation to b with matrix A using Matching pursuit[m,N] = size(A); x = zeros(N,1); res = b; support = []; % empty supportfor i=1:k

% compute correlation between residual and columns of Acorr = A’*res; % find position n (and value c) of the maximally correlated column [c n] = max(abs(corr)); % NB: modern Matlab notation allows [~ n] = max(abs(corr)) % extend the supportsupport(end+1) = n;% update the representation x(support) = pinv(A(:,support))*b; % update the residualres = res - A(:,support)*x(support);

end

105


Overview of the course

• Session 1: Introduction

• Session 2: Complexity & Feasibility✓ Difficulty of ideal sparse approximation✓ Greedy algorithms

• Session 3: Convex Pursuit Algorithms

• Session 4: Iterative Thresholding; Recovery Guarantees• Session 5: Compressive Sensing; Summary

106


Hard-thresholding (p=0)

• Solution of

107

H�(c)

c

minx

12(c� x)2 + � · |x|0

�2�

�⇥

2�


Soft-thresholding (p=1)

• Solution of

�

��

108

S�(c)

c

minx

12(c� x)2 + � · |x|


Matlab code ?

• Soft thresholding• @softthresh(c,lambda)(sign(c).*max(abs(c)-lambda,0))

• x = softthresh(c,lambda);

• Hard-thresholding• @hardthresh(c,lambda)(c.*(abs(c)>=sqrt(2*lambda)))

• x = hardthresh(c,lambda);

•

109


Iterative thresholding• Definition: proximity operator

• Goal = compute

• Iterative algorithm:✓ gradient descent on fidelity term

✓ thresholding

110

arg minx

12⇥Ax� b⇥2

2 + �⇥x⇥pp

x(i+1/2) := x(i) + �(i)AT (b�Ax(i))

�p�(c) = arg min

x

12(x� c)2 + �|x|p

x(i+1) := �p�(i)(x(i+1/2))


Iterative Thresholding

• Theorem : [Daubechies, de Mol, Defrise 2004, Combettes & Pesquet 2008]

✓ consider the iterates defined by the thresholding function, with

✓ assume that and✓ then, the iterates converge strongly to a limit

✓ the limit is a global minimum of

✓ if p>1, or if A is invertible, is the unique minimum

111

x(i+1) = f(x(i))

x�⇥x, ⇤Ax⇤2

2 � c⇤x⇤22 � < 2/c

x�

⇤x(i) � x�⇤2 ⇥i�⇥ 012⇥Ax� b⇥2

2 + �⇥x⇥pp

x�

p � 1

f(x) = �p�⇥(x + �AT (b�Ax))


Iterative Thresholding: convex penalties

• Strong convergence to global minimum

• Accelerated convergence: ✓ Nesterov schemes✓ see e.g. Beck & Teboulle 2009;

• Many variants of iterative thresholding ✓ depends on properties of penalty terms

✦ smoothness ✦ strong convexity✦ etc.

112


Iterative Thresholding: nonconvex penalties

• Example: Iterative Hard Thresholding for L0✓ keep components above threshold✓ or rather keep k largest components

✦ [IHT: Blumensath & Davies 2009]

• More generally, with nonconvex cost functions✓ Possible ‘spurious’ local minima✓ Convergence: fixed point, under certain assumptions ✓ Limit = global min: under certain assumptions (RIP)

• Pruning strategies:✓ ex: keep 2k components, project, keep k components

✦ ex: CoSAMP [Needell &Tropp 2008], ALPS [Cevher 2011], ...

113


Code for Iterative Thresholding?

• Proximal operator (or prox)

• Prox of the absolute value = soft-thresholding@prox(c,lambda)(sign(c).*max(abs(c)-lambda,0))

• Iterative thresholding with general proxfunction xhat = iterate_thresh(b,A,prox,step,niter)

xhat = 0;for i=1:niter

xhat =prox(xhat+ step * A’*(b-A*xhat))end

114

prox

f

(c) := arg min

x

⇢1

2

kx� ck22 + f(x)

�


Principleiterative decomposition• select new components• update residual

Tuning quality/sparsity regularization parameter

stopping criterion(nb of iterations, error level, ...)

Variants• choice of sparsity measure p• optimization algorithm • initialization

•selection criterion (weak, stagewise ...)•update strategy (orthogonal ...)

Iterative greedy algorithmsGlobal optimization

Summary

115

�

ri = b�Axi

⇥ri⇥ � �

minx

12⇥Ax� b⇥2

2 + �⇥x⇥pp

⇥xi⇥0 � k

NOVEMBER 2011R. GRIBONVAL - SPARSE METHODS - SISEA 2011

Two geometric viewpoints

• High-dim coeff domain • Low-dim signal domain

116

{x s.t.b = Ax}

Find closest subspacethrough correlations AT b

Find sparsest representationthrough (convex) optimization

b

2012 - R. GRIBONVAL - SPARSE METHODS - SISEA

Sparse recovery:Provably good algorithms?


• Unknown image with N pixels

• Partially observed image: ✓ m < N observed pixels

• Measurement matrix

Example: Inpainting Problem

118

b y

yy 2 RN

b[�p] = y[�p], �p 2 Observed

b = My


• Unknown image with N pixels

• Sparse Model in wavelet domain✦ wavelets coefficients are sparse

✦ sparse representation of unknown image

• Measurement matrix

Example: Inpainting Problem

119

yy 2 RN

b = My

y ⇡ �x

x ⇡ �T y

b ⇡M�x

A := M�


Signal space ~ RN

Set of signals of interest

Observation space ~ RM M<<N

Linear projection

Nonlinear Approximation =

Sparse recovery

Inverse problems

120

Courtesy: M. Davies, U. Edinburgh


Monte-Carlo simulations

121

�x0�0

P (x� = x0) x�p = arg min

Ax=Ax0�x�p

p=1

x0 b := Ax0P (x0) draw ground truth direct model

inverse problem

Typical observation

k1(A)


Exercice for next time

• Implement in Matlab / Scilab:✓ Matching Pursuit, Orthonormal Matching Pursuit✓ L1 minimization (with CVX)

• Generate test problems✓ Gaussian A + normalization of columns

• Measures of performance & computation time

• Curves of success as function of sparsity k

122

Rémi Gribonval INRIA Rennes - Bretagne Atlantique, France · Inverse problems and sparse models...

Documents

Transcript of Rémi Gribonval INRIA Rennes - Bretagne Atlantique, France · Inverse problems and sparse models...