Monday Morning Algorithm 12: Boosted Subspace Pursuit using Re-weigthed Lp for Compressive Sensing

[Warning: Before you read this entry, read this one first ]
Today's algorithm is going to be a composition of two algorithms already implemented. When I implemented the algorithm of Wei Dai and Olgica Milenkovic on Subspace Pursuit (Monday Morning Algorithm 10), my implementation and my simple example would work most of the time but not all the time. I also noticed that the final solution was close to the expected solution. Then I remembered what Emmanuel Candes said in his lecture on re-weighted L1 (specifically, the last one entitled Applications, experiments and open problems) : find a solution close enough and then use re-weighted L1 to refine the solution further. Since the re-weighted Lp solution technique of Rick Chartrand and Wotao Yin has also been implemented (Monday Morning Algorithm 2), I rewrote a line into it so that it would take the solution of the subspace Pursuit as the initial guess of its solution (instead of the L2 solution as is done in the original algorithm). By running SPtrial of the initial subspace pursuit implementation one may get this result:


The solution is wrong but not that wrong. By using this solution as a guess into the re-weighted Lp algorithm, one eventually obtains:

function x = lp_re1(A,y,p,x0)
%
% lp-Reweighted Least Squares Problem Solver
% Algorithm implementing an approach suggested in
% "Iteratively Reweighted Algorithms for Compressive Sensing"
% by Rick Chartrand and Wotao Yin
% Algorithm implemented as featured in:
% http://math.lanl.gov/Research/Publications/Docs/chartrand-2008-iteratively.pdf
%
%
% lp_re solves problems of the following form:
% minimize ||x||_p
%
% subject to A*x = y
%
% where A and y are problem data and x is variable (described below).
%
% CALLING SEQUENCES
% [x] = lp_re1(A,y,p,x0);
%
% INPUT
% A : mxn matrix; input data. columns correspond to features.
%
% m : number of measurements ( or rows) of A
% n : dimension of ambient space (or columns) of A
%
% y : vector of dimension m; result.
%
% p : norm order; p = 1 is l1 norm
%
% x0 : guess for x, vector of size n
%
% OUTPUT
% x : vector of vector n; initial unknown
%
% USAGE EXAMPLES
% x = lp_re(A,y,p,x0);
%
% AUTHOR Igor Carron
% UPDATE ver. 0.1, Feb 7 2008
%
% Creative Commons Licence
% http://en.wikipedia.org/wiki/Creative_Commons
%
epsilon = 1
% u_0 is a guess given to the subroutine
u_0 = x0;
u_old = u_0;
j=0;
while epsilon > 10^(-8)
j;
epsilon;
j = j + 1
w = (u_old.^(2) + epsilon).^(p/2-1);
v = 1./w;
Q_n = diag(v,0);
tu = inv(A*Q_n*A');
u_new = Q_n * A' * tu * y;
if lt(norm(u_new - u_old,2),epsilon^(1/2)/100)
epsilon = epsilon /10;
end
u_old = u_new;
end
x=u_new;

This is obviously a toy example that can be started by running SP1trial. In larger cases, one needs to pay attention to the issue Joel Tropp mentioned earlier (fast multiply) as re_lp1 is real slow otherwise. The code and example can be found on the Monday Morning Algorithm page.


Credit:NASA/JPL/GSFC/SwRI/SSI, Enceladus Heat map as seen by Cassini 19 days ago.