spgdemo.m

Spectral Projected L1 solver

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']);