Positive Constrained Minimization

eric_positive_armijo.m

function x = eric_positive_armijo(x0,f,g,h)
% Calculate a Newton-Step type minimization using the given 
% f function, its gradient g, and its hessian (second order)
% h.

scalar_mode = false;
if isscalar(x0)
	scalar_mode = true;
end

% Construct new barrier gradients
fb = @(x_,b_) f(x_) + fbarrier(x_,b_);
gb = @(x_,b_) g(x_) + gbarrier(x_,b_);
hb = @(x_,b_) h(x_) + hbarrier(x_,b_);

ns_maxit = 100;		% Maximum newton iterations
ls_maxit = 100;		% Maximum line search iterations
conv_tol = 1e-10;
barrier  = 10;
barrierdec = 0.5;
a = 1;
c = 0.5;
tau = 0.5;
x = x0;
for ns_iter = 1:ns_maxit
	G = gb(x,barrier);
	H = hb(x,barrier);
	xprev = x;

	if scalar_mode
		step = -G./H;
	else
		step = -H\G;
	end

	% Find the value of scaling term using armijo method
	fx  = fb(proj(x),barrier);
	fxp = fb(proj(x+(a.*step)),barrier);
	M   = step'*G;
	
	t = -c.*M;
    ls_iter = 1;
	while ((fx - fxp) < a*t) && (ls_iter < ls_maxit)
		a = tau*a;	          % Scale a down
		fxp = fb(proj(x+(a.*step)),barrier); % Recalculate the new score
		ls_iter = ls_iter + 1;
	end


	x = proj(x + a*step);
    fprintf('[%d] x : %0.2e | a : %0.2e (%d iter)\n',ns_iter,x,a,ls_iter);

	barrier = barrier.*barrierdec;
    
    if scalar_mode
        if (1 - x/xprev) < conv_tol
            break;
        end
    else        
        if mean((x-xprev).^2) < conv_tol
            break;
        end
    end
end

function x = proj(x)
	% Put everything back into positive territory
	x(x<=0) = 1e-18;


function f = fbarrier(x,b)
	f = -b.*log(x);

function g = gbarrier(x,b)
	g = -b./x;

function h = hbarrier(x,b)
	h = b./x.^2;