# mpf/spgl1

### Subversion checkout URL

You can clone with HTTPS or Subversion.

Fetching contributors…

Cannot retrieve contributors at this time

439 lines (366 sloc) 16.195 kb
 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 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 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 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 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 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 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 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 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
Something went wrong with that request. Please try again.