Revision f8badb04eb41adf950890050d90e6eb8963cba8f authored by Michael Friedlander on 12 June 2010, 22:27:58 UTC, committed by Michael Friedlander on 12 June 2010, 22:27:58 UTC
1 parent 00ab835
spgdemo.m
function spgdemo(interactive)
%DEMO Demonstrates the use of the SPGL1 solver
%
% See also SPGL1.
% demo.m
% $Id: spgdemo.m 1079 2008-08-20 21:34:15Z ewout78 $
%
% ----------------------------------------------------------------------
% This file is part of SPGL1 (Spectral Projected Gradient for L1).
%
% Copyright (C) 2007 Ewout van den Berg and Michael P. Friedlander,
% Department of Computer Science, University of British Columbia, Canada.
% All rights reserved. E-mail: <{ewout78,mpf}@cs.ubc.ca>.
%
% SPGL1 is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation; either version 2.1 of the
% License, or (at your option) any later version.
%
% SPGL1 is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
% Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public
% License along with SPGL1; if not, write to the Free Software
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
% USA
% ----------------------------------------------------------------------
if nargin < 1 || isempty(interactive), interactive = true; end
% Initialize random number generators
rand('state',0);
randn('state',0);
% Create random m-by-n encoding matrix and sparse vector
m = 50; n = 128; k = 14;
[A,Rtmp] = qr(randn(n,m),0);
A = A';
p = randperm(n); p = p(1:k);
x0 = zeros(n,1); x0(p) = randn(k,1);
% -----------------------------------------------------------
% Solve the underdetermined LASSO problem for ||x||_1 <= pi:
%
% minimize ||Ax-b||_2 subject to ||x||_1 <= 3.14159...
%
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve the underdetermined LASSO problem for \n');
fprintf('%% \n');
fprintf('%% minimize ||Ax-b||_2 subject to ||x||_1 <= 3.14159...\n');
fprintf('%% \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Set up vector b, and run solver
b = A * x0;
tau = pi;
x = spg_lasso(A, b, tau);
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf(['nonzeros(x) = %i, ' ...
'||x||_1 = %12.6e, ' ...
'||x||_1 - pi = %13.6e\n'], ...
length(find(abs(x)>1e-5)), norm(x,1), norm(x,1)-pi);
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Solve the basis pursuit (BP) problem:
%
% minimize ||x||_1 subject to Ax = b
%
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve the basis pursuit (BP) problem:\n');
fprintf('%% \n');
fprintf('%% minimize ||x||_1 subject to Ax = b \n');
fprintf('%% \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Set up vector b, and run solver
b = A * x0; % Signal
opts = spgSetParms('verbosity',1);
x = spg_bp(A, b, opts);
figure(1); subplot(2,4,1);
plot(1:n,x,'b', 1:n,x0,'ro');
legend('Recovered coefficients','Original coefficients');
title('(a) Basis Pursuit');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(a).\n');
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Solve the basis pursuit denoise (BPDN) problem:
%
% minimize ||x||_1 subject to ||Ax - b||_2 <= 0.1
%
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve the basis pursuit denoise (BPDN) problem: \n');
fprintf('%% \n');
fprintf('%% minimize ||x||_1 subject to ||Ax - b||_2 <= 0.1\n');
fprintf('%% \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Set up vector b, and run solver
b = A * x0 + randn(m,1) * 0.075;
sigma = 0.10; % Desired ||Ax - b||_2
opts = spgSetParms('verbosity',1);
x = spg_bpdn(A, b, sigma, opts);
figure(1); subplot(2,4,2);
plot(1:n,x,'b', 1:n,x0,'ro');
legend('Recovered coefficients','Original coefficients');
title('(b) Basis Pursuit Denoise');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(b).\n');
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Solve the basis pursuit (BP) problem in COMPLEX variables:
%
% minimize ||z||_1 subject to Az = b
%
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve the basis pursuit (BP) problem in COMPLEX variables:\n');
fprintf('%% \n');
fprintf('%% minimize ||z||_1 subject to Az = b \n');
fprintf('%% \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Create partial Fourier operator with rows idx
idx = randperm(n); idx = idx(1:m);
opA = @(x,mode) partialFourier(idx,n,x,mode);
% Create sparse coefficients and b = 'A' * z0;
z0 = zeros(n,1);
z0(p) = randn(k,1) + sqrt(-1) * randn(k,1);
b = opA(z0,1);
opts = spgSetParms('verbosity',1);
z = spg_bp(opA,b,opts);
figure(1); subplot(2,4,3);
plot(1:n,real(z),'b+',1:n,real(z0),'bo', ...
1:n,imag(z),'r+',1:n,imag(z0),'ro');
legend('Recovered (real)', 'Original (real)', ...
'Recovered (imag)', 'Original (imag)');
title('(c) Complex Basis Pursuit');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(c).\n');
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Sample the Pareto frontier at 100 points:
%
% phi(tau) = minimize ||Ax-b||_2 subject to ||x|| <= tau
%
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Sample the Pareto frontier at 100 points:\n');
fprintf('%% \n');
fprintf('%% phi(tau) = minimize ||Ax-b||_2 subject to ||x|| <= tau\n');
fprintf('%% \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
fprintf('\nComputing sample');
% Set up vector b, and run solver
b = A*x0;
x = zeros(n,1);
tau = linspace(0,1.05 * norm(x0,1),100);
phi = zeros(size(tau));
opts = spgSetParms('iterations',1000,'verbosity',0);
for i=1:length(tau)
[x,r] = spgl1(A,b,tau(i),[],x,opts);
phi(i) = norm(r);
if ~mod(i,10), fprintf('...%i',i); end
end
fprintf('\n');
figure(1); subplot(2,4,4);
plot(tau,phi);
title('(d) Pareto frontier');
xlabel('||x||_1'); ylabel('||Ax-b||_2');
fprintf('\n');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(d).\n');
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Solve
%
% minimize ||y||_1 subject to AW^{-1}y = b
%
% and the weighted basis pursuit (BP) problem:
%
% minimize ||Wx||_1 subject to Ax = b
%
% followed by setting y = Wx.
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve \n');
fprintf('%% \n');
fprintf('%% (1) minimize ||y||_1 subject to AW^{-1}y = b \n');
fprintf('%% \n');
fprintf('%% and the weighted basis pursuit (BP) problem: \n');
fprintf('%% \n');
fprintf('%% (2) minimize ||Wx||_1 subject to Ax = b \n');
fprintf('%% \n');
fprintf('%% followed by setting y = Wx. \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Sparsify vector x0 a bit more to get exact recovery
k = 9;
x0 = zeros(n,1); x0(p(1:k)) = randn(k,1);
% Set up weights w and vector b
w = rand(n,1) + 0.1; % Weights
b = A * (x0 ./ w); % Signal
% Run solver for both variants
opts = spgSetParms('iterations',1000,'verbosity',1);
AW = A * spdiags(1./w,0,n,n);
x = spg_bp(AW, b, opts);
x1 = x; % Reconstruct solution, no weighting
opts = spgSetParms('iterations',1000,'verbosity',1,'weights',w);
x = spg_bp(A, b, opts);
x2 = x .* w; % Reconstructed solution, with weighting
figure(1); subplot(2,4,5);
plot(1:n,x1,'m*',1:n,x2,'b', 1:n,x0,'ro');
legend('Coefficients (1)','Coefficients (2)','Original coefficients');
title('(e) Weighted Basis Pursuit');
fprintf('\n');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(e).\n');
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Solve the multiple measurement vector (MMV) problem
%
% minimize ||Y||_1,2 subject to AW^{-1}Y = B
%
% and the weighted MMV problem (weights on the rows of X):
%
% minimize ||WX||_1,2 subject to AX = B
%
% followed by setting Y = WX.
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve the multiple measurement vector (MMV) problem \n');
fprintf('%% \n');
fprintf('%% (1) minimize ||Y||_1,2 subject to AW^{-1}Y = B \n');
fprintf('%% \n');
fprintf('%% and the weighted MMV problem (weights on the rows of X): \n');
fprintf('%% \n');
fprintf('%% (2) minimize ||WX||_1,2 subject to AX = B \n');
fprintf('%% \n');
fprintf('%% followed by setting Y = WX. \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Initialize random number generator
randn('state',0); rand('state',0);
% Create problem
m = 100; n = 150; k = 12; l = 6;
A = randn(m,n);
p = randperm(n); p = p(1:k);
X0= zeros(n,l); X0(p,:) = randn(k,l);
weights = 3 * rand(n,1) + 0.1;
W = spdiags(1./weights,0,n,n);
B = A*W*X0;
% Solve unweighted version
opts = spgSetParms('verbosity',1);
x = spg_mmv(A*W,B,0,opts);
x1 = x;
% Solve weighted version
opts = spgSetParms('verbosity',1,'weights',weights);
x = spg_mmv(A,B,0,opts);
x2 = spdiags(weights,0,n,n) * x;
% Plot results
figure(1); subplot(2,4,6);
plot(x1(:,1),'b-'); hold on;
plot(x2(:,1),'b.');
plot(X0,'ro');
plot(x1(:,2:end),'-');
plot(x2(:,2:end),'b.');
legend('Coefficients (1)','Coefficients (2)','Original coefficients');
title('(f) Weighted Basis Pursuit with Multiple Measurement Vectors');
fprintf('\n');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(f).\n');
fprintf([repmat('-',1,80), '\n']);
fprintf('\n\n');
% -----------------------------------------------------------
% Solve the group-sparse Basis Pursuit problem
%
% minimize sum_i ||y(group == i)||_2
% subject to AW^{-1}y = b,
%
% with W(i,i) = w(group(i)), and the weighted group-sparse
% problem
%
% minimize sum_i w(i)*||x(group == i)||_2
% subject to Ax = b,
%
% followed by setting y = Wx.
% -----------------------------------------------------------
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('%% Solve the group-sparse Basis Pursuit problem \n');
fprintf('%% \n');
fprintf('%% (1) minimize sum_i ||y(group == i)||_2 \n');
fprintf('%% subject to AW^{-1}y = b, \n');
fprintf('%% \n');
fprintf('%% with W(i,i) = w(group(i)), and the weighted group-sparse\n');
fprintf('%% problem \n');
fprintf('%% \n');
fprintf('%% (2) minimize sum_i w(i)*||x(group == i)||_2 \n');
fprintf('%% subject to Ax = b, \n');
fprintf('%% \n');
fprintf('%% followed by setting y = Wx. \n');
fprintf(['%% ', repmat('-',1,78), '\n']);
fprintf('\nPress <return> to continue ... \n');
if interactive, pause; end
% Initialize random number generator
randn('state',0); rand('state',2); % 2
% Set problem size and number of groups
m = 100; n = 150; nGroups = 25; groups = [];
% Generate groups with desired number of unique groups
while (length(unique(groups)) ~= nGroups)
groups = sort(ceil(rand(n,1) * nGroups)); % Sort for display purpose
end
% Determine weight for each group
weights = 3*rand(nGroups,1) + 0.1;
W = spdiags(1./weights(groups),0,n,n);
% Create sparse vector x0 and observation vector b
p = randperm(nGroups); p = p(1:3);
idx = ismember(groups,p);
x0 = zeros(n,1); x0(idx) = randn(sum(idx),1);
b = A*W*x0;
% Solve unweighted version
opts = spgSetParms('verbosity',1);
x = spg_group(A*W,b,groups,0,opts);
x1 = x;
% Solve weighted version
opts = spgSetParms('verbosity',1,'weights',weights);
x = spg_group(A,b,groups,0,opts);
x2 = spdiags(weights(groups),0,n,n) * x;
% Plot results
figure(1); subplot(2,4,7);
plot(x1); hold on;
plot(x2,'b+');
plot(x0,'ro'); hold off;
legend('Coefficients (1)','Coefficients (2)','Original coefficients');
title('(g) Weighted Group-sparse Basis Pursuit');
fprintf('\n');
fprintf([repmat('-',1,35), ' Solution ', repmat('-',1,35), '\n']);
fprintf('See figure 1(g).\n');
fprintf([repmat('-',1,80), '\n']);
end % function demo
function y = partialFourier(idx,n,x,mode)
if mode==1
% y = P(idx) * FFT(x)
z = fft(x) / sqrt(n);
y = z(idx);
else
z = zeros(n,1);
z(idx) = x;
y = ifft(z) * sqrt(n);
end
end % function partialFourier
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...