Compressed Sampling
This numerical tour explores compressed/copressive sensing/sampling reconstruction with both greedy matching pursuit and convex L1 regularization.
Contents
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab commands you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Compressed Sensing Acquisition
Compressed sensing acquisition corresponds to a random projection of a signal x on a few linear vectors. For the recovery of x to be possible, it is assumed to be sparsely represented in an orthogonal basis. Up to a change of basis, we suppose that the vector x itself is sparse.
We begin by generating a sparse signal with s randomized values.
% dimension n = 512; % sparsity s = 17; % location of the diracs set_rand_seeds(123456,123456); sel = randperm(n); sel = sel(1:s); x = zeros(n,1); x(sel)=1; % % randomization of the signs and values x = x.*sign(randn(n,1)).*(1-.5*rand(n,1)); % save for future use x0 = x;
We create a random Gaussian matrix and normalize its columns. This defines the acquisition operator U. The columns are thus random point on the unit sphere of R^p.
% number of measures p = 100; % Gaussian matrix U = randn(p,n); % normalization U = U ./ repmat( sqrt(sum(U.^2)), [p 1] );
We perform random measurements without noise.
y = U*x;
Compressed sensing recovery corresponds to solving the inverse problem y=U*x, which is ill posed because x is higher dimensional than x. The reconstruction can be perform with L1 minimization, which regularizes the problems by using the sparsity of the solution.
min_{x1} norm(x1,1) subject to U*x1=y.
Using the homotopy algorithm, we compute the solution path for a large number of different regularization parameter lambda values, and retrieve the last step, which corresponds to the minimum L1 norm solution.
options.niter = Inf; [X1,lambda_list,sparsity_list] = perform_homotopy(U,y,options); xbp = X1(:,end);
We display the solution. You can see that the location of the Dirac is perfectly recovered.
clf; subplot(2,1,1); plot_sparse_diracs(x); set_graphic_sizes([], 20); title('Original Signal'); subplot(2,1,2); plot_sparse_diracs(xbp); set_graphic_sizes([], 20); title('Recovered by L1 minimization');
Exercice 1: (the solution is exo1.m) For sparsity values s in [14,34], measure the ratio of vector that are perfectly recovered. To do this, generate many input signal (something like 100) and perform CS recovery on these signals to see wether they are recovered or not.
exo1;
Compressed Sensing Recovery with Orthogonal Matching Pursuit
Instead of using L1 minimization, it is possible to use matching pursuit using the columns of U as atoms of a redundant dictionary.
x = x0; y = U*x; s = sum(x~=0);
We perform s steps of orthogonal pursuit. This also recovers perfectly the signal.
options.nbr_max_atoms = s; xomp = perform_omp(U,y,options); clf; subplot(2,1,1); plot_sparse_diracs(x); set_graphic_sizes([], 20); title('Original Signal'); subplot(2,1,2); plot_sparse_diracs(xomp); set_graphic_sizes([], 20); title('Recovered by OMP');
Exercice 2: (the solution is exo2.m) Solve the compressed sensing reconstruction using s steps of orthogonal matching pursuits.
exo2;
Compressed Sensing on Compressible Signals
Exactely s-sparse signals are of little practical use because natural signals are compressible rather than exactely sparse. We study here the relative performance of non-linear approximation and compressive sampling approximation on signals with polynomially decaying coefficients.
A typical exemple of signal that is compressible in the Dirac basis is a signal with fast decaying coefficients. We flip here the sign, although it is not relevant (because the sensing matrix is random).
n = 512; alpha = 1.5; x = (1:n)'.^(-alpha); x(2:2:n) = -x(2:2:n);
Random signal are generated from this template by randomizing the coefficients.
clf; subplot(2,1,1); plot_sparse_diracs( x(randperm(n)) ); subplot(2,1,2); plot_sparse_diracs( x(randperm(n)) );
The measure of sparsity here is alpha, which should be >1 for compressible signal, and in (1/2,1) for not so compressible signal.
The non-linear approximation error with M coefficients of such a signal decays like norm(f-f_M)^2=M^{-2*alpha+1}, so that the slope in log/log plot is -2*alpha+1.
xs = sort(abs(x)); if xs(1)>xs(n) xs = reverse(xs); end err_nl = reverse(cumsum( xs.^2 )); err_nl = err_nl/err_nl(1); % display clf; plot( log10(1:n), log10(err_nl) ); axis('tight'); set_label('log10(M)', 'log10(|f-f_M|^2)');
Exercice 3: (the solution is exo3.m) Compute the alpha slope for the wavelet coefficients of a natural image (for instance 'lena'). You should perform a linear fit of the slope of ordered coefficients for M for instance in the range [.005 .1]*n^2 (where n^2 is the number of pixels). Use the function polyfit to compute the linear regression.
exo3;
The compressed sensing recovery from a compressible signal from noiseless measurements can be performed by pure L1 minimization. It is however better to use a relaxed L1 minimzation
min_{x1} 1/2*norm(y-U*x1)^2 + lambda*norm(x1,1)
for a well chosen lambda.
We compute the full regularization path for all the lambda.
% measures y = U*x; % L1 resolution options.niter = Inf; [X1,lambda_list,sparsity_list] = perform_homotopy(U,y,options);
We can display the L2 compressed sensing approximation error as a function of the iteration step. You can see that it is not strictly decaying. A slight regularization somehow decreases the error.
m = size(X1,2); err = sum( (X1 - repmat(x, [1 m])).^2 ); clf; plot(err / norm(x)); axis([1 m 0 .2]); set_label('#step', 'Error');
The number of samples is p. We can compute the number of q of sample that generates the same non-linear approximation error. The ratio p/q>1 gives the compressed sensing sub-optimality constant. The closer to 1, the better.
% the best possible compressed sensing error [e,i] = compute_min(err(:), 1); xcs = X1(:,i); % compute the equivalent number of NL terms [tmp,q] = compute_min(abs(e-err_nl(:)), 1); disp(strcat(['Compressed sensing sub-optimality factor: p/q=' num2str(p/q,2)]));
Compressed sensing sub-optimality factor: p/q=9.1