bp_simplex.m

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

function c = BP_Simplex( x, F, c0 )

% BP_Simplex.m -- Basis Pursuit
%----------------------------------------------------------------------
%
%       Min || c ||_1  subject to \Phi * c = x
%		
%    via a dense implementation of the Simplex method.
%
% Usage
%    c = BP_Simplex( x, F[, c0] );
%
% Inputs
%    x   1-d signal; column vector
%    F   Dictionary matrix; each column corresponds to an atom
%    c0  Initialization for c; default = 0
%
% Outputs
%    c   Coefficients of BP representation; column vector
%
% Description
%    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
%    BP_Interior, BP_Interior2, lssol
%----------------------------------------------------------------------

%----------------------------------------------------------------------
% 18 Mar 1997: Original version of BP_Simplex.m.
%              Shaobing Chen, Dept of Statistics, Stanford University.
% 14 Apr 2001: (M. Saunders) Minor editing.
%----------------------------------------------------------------------

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

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

% Set up the parameters for lssol.
  l   = [zeros(m2,1);         xx];
  u   = [ ones(m2,1) * 1e+20; xx];
  tp  = 'lp';
  F   = [F (-F)];
  c   = ones(m2,1);
  msg = 0;
  kx  = [1:m2]';

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

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

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