FRANALASSO - Frame LASSO regression

Usage

tc = franalasso(F,f,lambda)
tc = franalasso(F,f,lambda,C,tol,maxit)
[tc,relres,iter,frec,cd] = franalasso(...)

Input parameters

Output parameters

Description

franalasso(F,f,lambda) solves the LASSO (or basis pursuit denoising) regression problem for a general frame: minimize a functional of the synthesis coefficients defined as the sum of half the \(l^2\) norm of the approximation error and the \(l^1\) norm of the coefficient sequence, with a penalization coefficient lambda such that

The solution is obtained via an iterative procedure, called Landweber iteration, involving iterative soft thresholdings.

The following flags determining an algorithm to be used are recognized at the end of the argument list:

'ista'

The basic (Iterative Soft Thresholding) algorithm given F and f and parameters C\(>0\), lambda \(>0\) acts as follows:

Initialize c
repeat until stopping criterion is met
  c <- soft(c + F*(f - Fc)/C,lambda/C)
end

'fista'

The fast version of the previous. This is the default option.

Initialize c(0),z,tau(0)=1,n=1
repeat until stopping criterion is met
   c(n) <- soft(z + F*(f - Fz)/C,lambda/C)
   tau(n) <- (1+sqrt(1+4*tau(n-1)^2))/2
   z   <- c(n) + (c(n)-c(n-1))(tau(n-1)-1)/tau(n)
   n   <- n + 1
end

[tc,relres,iter] = franalasso(...) returns the residuals relres in a vector and the number of iteration steps done iter.

[tc,relres,iter,frec,cd] = franalasso(...) returns the reconstructed signal from the coefficients, frec and coefficients cd obtained by analysing using the canonical dual system. Note that this requires additional computations.

The relationship between the output coefficients and frec is given by

frec = frsyn(F,tc);

The function takes the following optional parameters at the end of the line of input arguments:

'C',cval

Landweber iteration parameter: must be larger than square of upper frame bound. Default value is the upper frame bound.

'tol',tol

Stopping criterion: minimum relative difference between norms in two consecutive iterations. Default value is 1e-2.

'maxit',maxit

Stopping criterion: maximal number of iterations to do. Default value is 100.

'print'

Display the progress.

'quiet'

Don't print anything, this is the default.

'printstep',p

If 'print' is specified, then print every p'th iteration. Default value is 10;

'debias'

Performs debiasing (improves approximation accuracy) of the solution using the conjugate gradient algorithm. This amounts to doing best-fit with respect to a subdictionary selected by ISTA.

'pcgmaxit',pcgmaxit

Maximum allowed number of iterations for pcg. Only used together with the 'debias' flag. The default value is 100.

'pcgtol',pcgtol

Tolerance of pcg. Only used together with the 'debias' flag. The default value is 1e-6.

The parameters C, itermax and tol may also be specified on the command line in that order: franalasso(F,x,lambda,C,tol,maxit).

Note: If you do not specify C, it will be obtained as the upper framebound. Depending on the structure of the frame, this can be an expensive operation.

Examples:

The following example shows how franalasso produces a sparse representation of a test signal greasy:

f = greasy;
% Gabor frame with redundancy 8
F = frame('dgtreal','gauss',64,512);
% Choosing lambda (weight of the sparse regularization param.)
lambda = 0.1;
% Solve the basis pursuit problem
[c1,~,~,frec1,cd] = franalasso(F,f,lambda);
% Solve again with debiasing enabled
[c2,~,~,frec2,cd] = franalasso(F,f,lambda,'debias');
% Plot sparse coefficients
figure(1); plotframe(F,c1,'dynrange',50);
% Plot debiased coefficients
figure(2); plotframe(F,c2,'dynrange',50);
% Plot coefficients obtained by applying an analysis operator of a
% dual Gabor system to f
figure(3); plotframe(F,cd,'dynrange',50);

% Normalized MSE of the approximation
% frec is obtained by applying the synthesis operator of frame F
% to sparse coefficients c.
fprintf('Normalized MSE:                  %.2f dB\n',20*log10(norm(f-frec1)/norm(f)));
fprintf('Normalized MSE (with debiasing): %.2f dB\n',20*log10(norm(f-frec2)/norm(f)));

% Compare decay of coefficients sorted by absolute values
% (compressibility of coefficients)
figure(3);
semilogx([sort(abs(c1),'descend')/max(abs(c1)),...
sort(abs(c2),'descend')/max(abs(c2)),...
sort(abs(cd),'descend')/max(abs(cd))]);
legend({'sparsified coefficients','with debiasing','dual system coefficients'});