bpdn_simplex.m

来自「％ Atomizer Main Directory, Version .802 」· M 代码 · 共 82 行

function [xrec, coef] = BPDN_Simplex( y, F, c0 )

% BPDN_Simplex.m -- Basis Pursuit DeNoising
%----------------------------------------------------------------------
%
%       min .5 * || y - F * c ||_2^2 + lambda * || c ||_1
%
%    via a dense implementation of the Simplex method.
%
% Usage
%       [xrec, coef] = BPDN_Simplex( y, F[, c0] );
%
% Inputs
%    y     1-d signal to be de-noised; column vector
%    F     Dictionary matrix with each column corresponding to an atom
%    c0    Initialization for coef; default = 0
%
% Outputs
%    xrec  Cleaned signal; column vector
%	c  Coefficients of cleaned signal; column vector
%
% Description
%    1. Assumes noise level == 1.
%    2. BP_Simplex uses the LSSOL package, a simplex code for linear
%       programming, quadratic programming and constrained linear
%       least squares.
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% Reference
%    P. E. Gill, S. J. Hammarling, W. Murray, M. A. Saunders
%    and M. H. Wright,
%    User's guide for LSSOL (Version 1.0):
%    A Fortran package for constrained linear least-squares
%    and convex quadratic programming,
%    Report SOL 86-1, Systems Optimization Laboratory,
%    Department of Operations Research, Stanford University, Jan 1986.
%
%  See also
%    BPDN_Interior, BPDN_Interior2, lssol
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 09 Jan 1996: Original version of BPDN_Simplex.m.
%              Shaobing Chen, Dept of Statistics, Stanford University.
% 14 Apr 2001: (M. Saunders) Minor editing.
%----------------------------------------------------------------------

  y      = y(:);
  [n, m] = size(F);
  m2     = 2*m;

% Initialization.
  if nargin < 3
     x0 = zeros(m2,1);
  else
     x0 = [c0 .* (c0 >= 0); - c0 .* (c0 < 0)];
  end

% Set up the parameters for lssol.
  thresh = (2 * log(m)) ^ .5;
  c      = ones(m2, 1) * thresh;
  b      = y;
  F      = [F (-F)];
  l      = zeros(m2,1);
  u      =  ones(m2,1) * 1e+20;
  tp     = 'ls2';
  msg    = 0; 
  kx     = [1:m2]';

% Run lssol.
  fprintf('\nBPDN_Simplex:\n');
  time   = cputime;

  [coef, obj, lambda, istate, iter, inform, A, b, kx] = ...
      lssol(F, b, c, [], x0, l, u, tp, msg, kx);

  coef  = coef(1:m) - coef((m+1):m2);
  xrec  = F(:, 1:m) * coef;
  fprintf('\nlssol time =%10.1f    sum(c) =%12.5e\n', time, sum(abs(coef)));