SPGL1: A Brief Tour
Michael P. Friedlander, June 2008
SPGL1 is a solver for one-norm regularized least-squares. It is suitable for problems that are large-scale and in the complex domain. This document gives a brief tour of SPGL1's main features.
The examples below are based on randomly generated data. We'll begin by creating a random m-by-n encoding matrix and a sparse solution vector:
rand('twister',0); randn('state',0); m = 50; n = 128; k = 14; % No. of rows (m), columns (n), and nonzeros (k) [A,Rtmp] = qr(randn(n,m),0); % Random encoding matrix with orthogonal rows A = A'; % ... A is m-by-n p = randperm(n); p = p(1:k); % Location of k nonzeros in x x0 = zeros(n,1); x0(p) = randn(k,1); % The k-sparse solution b = A*x0; % The right-hand side corresponding to x0
The basis pursuit (BP) problem
minimize ||x||_1 subject to Ax = b
often finds the sparsest solution to the linear system of equations Ax = b. The call to SPGL1 is straightforward:
opts = spgSetParms('verbosity',0); % Turn off the SPGL1 log output x = spg_bp(A, b, opts);
Let's verify that the true sparse solution was recovered correctly:
plot(x0,'r*'); hold on stem(x ,'b '); hold off legend('Original coefficients','Recovered coefficients'); title('Basis Pursuit');
As a final check, the residual should be nearly zero:
norm(b - A*x)
ans = 8.3531e-05
If there was noise inherent in the signal that we measured, it would be better to instead solve the basis pursuit denoise (BPDN) problem
minimize ||x||_1 subject to ||Ax - b||_2 <= sigma,
which only approximately fits the equation Ax=b. In this experiment, we'll set sigma = 0.1:
b = A * x0 + randn(m,1) * 0.075; sigma = 0.1; % Desired ||Ax - b||_2 opts = spgSetParms('verbosity',0); x = spg_bpdn(A, b, sigma, opts);
We don't expect the recovery to be exact this time because of the noise:
plot(x0,'r*'); hold on stem(x ,'b '); hold off legend('Original coefficients','Recovered coefficients'); title('Basis Pursuit Denoise');
The corresponding residual should have a norm (nearly) equal to sigma (0.1):
norm(b - A*x)
ans = 0.0999
The Lasso problem puts a constraint directly on the one-norm of the solution:
minimize ||Ax - b||_2 subject to ||x||_1 <= tau
Here, we'll recover a least-squares solution with a one-norm equal to pi -- an easily recognizable number! We'll turn on SPGL1's output so that we can see it's progress:
opts = spgSetParms('verbosity',1); tau = pi; x = spg_lasso(A, b, tau, opts);
================================================================================ SPGL1 v. 1017 (Mon, 16 Jun 2008) ================================================================================ No. rows : 50 No. columns : 128 Initial tau : 3.14e+00 Two-norm of b : 2.61e+00 Optimality tol : 1.00e-04 Target one-norm of x : 3.14e+00 Basis pursuit tol : 1.00e-06 Maximum iterations : 500 Iter Objective Relative Gap gNorm stepG nnzX nnzG 0 2.6116134e+00 8.8203906e-01 9.57e-01 0.0 0 0 1 1.7351041e+00 1.5596782e-01 5.16e-01 0.0 10 1 2 1.7035044e+00 3.5470990e-02 4.43e-01 0.0 6 1 3 1.6985890e+00 1.0356556e-02 4.16e-01 0.0 6 1 4 1.6984830e+00 1.3137969e-03 4.09e-01 0.0 6 3 5 1.6984820e+00 3.5582204e-04 4.09e-01 0.0 6 6 6 1.6984819e+00 2.5425067e-05 4.09e-01 0.0 6 6 EXIT -- Optimal solution found Products with A : 7 Total time (secs) : 0.0 Products with A' : 7 Project time (secs) : 0.0 Newton iterations : 0 Mat-vec time (secs) : 0.0 Line search its : 0 Subspace iterations : 0
The one-norm of the solution:
ans = 3.1416
Many signal-processing problems involve complex-valued coefficients. In this example, we create a solution that is sparse in the Fourier domain. We simulate missing data by randomly restricting the rows of the Fourier operator. The resulting BP problem is
minimize ||z||_1 subject to R*F*z = R*b,
where R is a restriction operator and F is a Fourier operator. Keep in mind that z is a complex variable, and that
||z||_1 = sum_j ( sqrt( x_j^2 + y_j^2 ) ), where z_j = x_j + i*y_j.
We use the partial Fourier operator
z = partialFourier(idx,n,x,mode);
the index set idx
idx = randperm(n); idx = idx(1:m);
determines which rows of the Fourier operator are kept. We create an operator that conforms to SPGL1's requirements:
opA = @(x,mode) partialFourier(idx,n,x,mode); % This is now "A"
Create sparse coefficients and the right-hand side b:
z0 = zeros(n,1); z0(p) = randn(k,1) + sqrt(-1) * randn(k,1); Rb = opA(z0,1); % Rb = R*A*z0
Fire up SPG_BP in the usual way. Note that in this, case, the first argument is a function:
opts = spgSetParms('verbosity',0); z = spg_bp(opA,Rb,opts);
Plot the recovered and original coefficients:
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('Complex Basis Pursuit');
The Pareto curve traces the optimal tradeoff between the one-norm of the solution and the two-norm of the residual. We take advantage here of SPGL1's warm-start capability and sample the Pareto curve at a bunch of points.
We do this by solving a sequence of LASSO problems for increasing values of tau. Larger values of tau will yield smaller residuals. In fact, if phi(tau) is the norm of the residual corresponding to a given tau, then phi is a continuously differentiable function [BergFriedlander08].
Because here we wish to warm-start SPGL1, we need to use its expert interface, rather than the simpler interface available via SPG_LASSO.
b = A*x0; tau = linspace(0,1.01*norm(x0,1),50); % Values of tau at which to sample phi = zeros(size(tau)); % Preallocate the vector of residuals. x = zeros(n,1); % Initialize the solution. opts = spgSetParms('verbosity',0); for i=1:length(tau) [x,r] = spgl1(A,b,tau(i),,x,opts); % Re-solve LASSO(tau) with a different tau. phi(i) = norm(r,2); % SPGL1's 2nd output arg is the residual. end
Plot the curve.
plot(tau,phi,'b.'); title('Pareto frontier'); xlabel('||x||_1'); ylabel('||Ax-b||_2');
[BergFriedlander] E. van den Berg and M. P. Friedlander, "Probing the Pareto frontier for basis pursuit solutions", January 2008 (revised May 2008). To appear in SIAM Journal on Scientific Computing.
% $Id: spgexamples.m 1077 2008-08-20 06:15:16Z ewout78 $