oneprojector.m

Spectral Projected L1 solver

function [x,itn] = oneProjector(b,d,tau)
% ONEPROJECTOR  Projects b onto the weighted one-norm ball of radius tau
%
%    [X,ITN] = ONEPROJECTOR(B,D,TAU) returns the orthogonal
%    projection of the vector b onto the weighted one-norm ball or
%    radius tau. That is, it solves
%
%            minimize  ||b-x||_2  st  ||Dx||_1 <= tau.
%               x
%
%    If D is empty or omitted all weights are set to one. The
%    return value ITN given the number of elements of B that were
%    thresholded. 
%
% See also spgl1.

%   oneProjector.m
%   $Id $
%
%   ----------------------------------------------------------------------
%   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
%   ----------------------------------------------------------------------

% Check arguments
if nargin < 2
  error('The oneProjector function requires at least two parameters');
end
if nargin < 3
  tau = d;
  d   = [];
end

% Check weight vector
if isempty(d), d = 1; end;

if ~isscalar(d) && (length(b) ~= length(d))
  error('Vectors b and d must have the same length');
end

if isscalar(d) && (d == 0)
  x = b; itn = 0; return;
end


% Get sign of b and set to absolute values
s = sign(b); b = abs(b);

% Perform the projection
if isscalar(d)
  [x,itn] = oneProjectorMex(b,tau/d);
else
  d   = abs(d);
  idx = find(d > eps); % Get index of all non-zero entries of d
  x   = b;             % Ensure x_i = b_i for all i not in index set idx
  [x(idx),itn] = oneProjectorMex(b(idx),d(idx),tau);
end

% Restore signs in x
x = x.*s;