bpdn_primal.m

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

function [xrec, coef] = BPDN_Primal(y, NameOfDict, par1, par2, par3, BPAccuracy, CGAccuracy)
% BPDN_Primal -- Basis Pursuit DeNoising:
%
%		Min .5 * || y - \Phi * c ||_2^2 + lambda * || c ||_1
%
%		via a Primal Logrithmic Barrier Interior Point method.
%  Usage
%	[xrec, coef] = BPDN_Primal(y, NameOfDict, par1, par2, par3
%			[, BPAccuracy, CGAccuracy])
%  Inputs
%       y               1-d signal to be de-noised; column vector
%	NameOfDict	string; name of the dictionary
%	par1,par2,par3	parameters of the dictionary
%	BPAccuracy	bound on the barrier parameter, default = 1e-3
%	CGAccuracy	accuracy of the CG solver, default = 1e-3
%
%	Use 'help dictionary' for dictionary objects: NameOfDict,par1,par2,par3
%  Outputs
%       xrec            cleaned 1-d signal; column vector
%       coef            coeff of cleaned signal; column vector
%  Description
%	1. Assumes noise level == 1
%	2. BPDN_Interior stops when the barrier parameter \mu gets below
%	   the thrshold defined by the parameter BPAccuracy.
%	3. In order to achieve more accurate optimility, one should set
%	   both BPAccuracy and CGAccuracy smaller. We found the default
%	   values give satisfying approximate optimum
%  Reference
%	S. Chen
%	``Basis Pursuit''
%	Ph.D. thesis, Department of Statistics, Stanford Univ., 1995
%
%	M.A. Saunders.
%	Cholesky-based methods for sparse least squares:
%		the benefits of regularization.
%	Technical Report SOL 95-1, Stanford University, July 1995.
%  See Also
%	BPDN_Simplex
%

%*** Size of the problem
y = y(:);
n = length(y);
[m L] = SizeOfDict(n, NameOfDict, par1, par2, par3);

%*** Setting the controlling parameters
if nargin < 7,
	CGAccuracy = 1e-3;
end
if nargin < 6,
	BPAccuracy = 1e-3;
end
MaxCGIter = log2(m)*m;
lambda = (2 * log(m)) ^ .5;

%*** Initializing x, mu
fprintf('\nBPDN-Primal:\n\n');
disp('Initializing BPDN-Primal by MOF ...')
x = MOF(y, NameOfDict, par1, par2, par3);
x = [ x .* (x > 0) ; (-x) .* (x < 0) ];
x = max(x, 1);
mu = 1e-1;
itn = 0;
converged = 0;
zeros2m = zeros(2*m, 1);

%*** BPDN Iterations
disp(' ');
disp('BPDN-Interior Iterations ...')
disp('       Itn        mu            Step            obj       CGStep')
r = y - FastSS(x, n, NameOfDict, par1, par2, par3);
obj = lambda * sum(x) + r' * r / 2;
fprintf('Step %4g: %14.7e %14.7e %14.7e %4g\n', 0, mu, 0, obj, 0);


while ~converged,
	r = y - FastSS(x, n, NameOfDict, par1, par2, par3);
	obj = lambda * sum(x) + r' * r / 2; 

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%% (DA'AD + mu I) z = mu e + D A' r - Dc
	%%%Solve the equation by Conjugate Gradient
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	r = mu + x .* FastAA(r, NameOfDict, par1, par2, par3) - lambda * x;
	thresh = CGAccuracy * norm(r); thresh = thresh^2;

	z = zeros2m;
        rho_1 = norm(r) ^ 2;
        for k = 1:MaxCGIter,
                %fprintf('      CG Step: %3g:   %g\n', k, rho_1');
                if rho_1  < thresh,
                        break;
                end
                if k == 1,
                        p = r;
                else
                        beta = rho_1 / rho_0;
                        p = r + beta * p;
                end
		w = x .* FastAA(FastSS(x.*p, n, NameOfDict, par1, par2, par3), ...
				NameOfDict, par1, par2, par3) + mu * p;
                alpha = rho_1 / (p' * w);
                r = r - alpha * w;
		rho_0 = rho_1;
                rho_1 = norm(r) ^2;
                z = z + alpha * p;
        end

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if k == MaxCGIter,
		disp('Too Many CG Iterations');
		break;
	end

	dx = x .* z;
	
	negdx = dx < -1e-14;
	if sum(negdx) > 0,
		R = min(x(negdx) ./ abs(dx(negdx)));
	else
		disp('Optimality achieved!');
		break;
	end
	xold = x;
	x = x + .99 * R * dx;
	mu = (1 - min(R, .99)) * mu;
	itn = itn + 1;
	converged = mu < BPAccuracy;
	Step = .99 * R * norm(dx) / norm(xold);
	fprintf('Step %4g: %14.7e %14.7e %14.7e %4g\n', itn, mu, Step, obj, k);
end

coef = x(1:m) - x((m+1):(2*m));
xrec = FastSS(x, n, NameOfDict, par1, par2, par3);
fprintf('Obj = %14.7e\n', obj);