function rosenbrock_extrust

% nearly-exact solution of the trust-region problem

x0=[-1.2; 1];
xs=[1;1];


% Trust-region/dogleg algorithm
Delta = 1.0;
xk    = x0;
df0   = df(x0);

Nmax  = 1000;

for i = 1:Nmax,
    g =  df(xk);
    B = df2(xk);
    %
    printf('%3d: f=%15.6e, |df|=%15.6e, |x-x*| = %15.6e Delta = %15.6e\n', ...
           i, f(xk), norm(g), norm(xk-xs), Delta);
    %
    if norm(g) <= 1.0E-12*norm(df0)
        % converged!
        break;
    end
    %
    p = solve_trust(g,B,Delta);
    %
    %  check the actual reduction
    rho = (f(xk) - f(xk + p)) / (-g'*p - p'*B*p/2);
    % update the trust-region radius
    if rho < 0.25
        Delta = 0.25*Delta;
    elseif rho > 0.75 && Delta-norm(p) >= -1.0E-06*Delta
        Delta = min(2*Delta, 2.0);
    end
    % accept the step, if the ratio is OK
    if rho > 1.0E-03
        xk = xk+p;
    end
end

function pC = cauchy(g,B,Delta)
if g'*B*g <= 0
    pC = -Delta/norm(g)*g;
else
    alphaS = g'*g/(g'*B*g);
    if norm(alphaS*g) <= Delta
        pC = -alphaS*g;
    else
        pC = -Delta/norm(g)*g;
    end
end

function p = solve_trust(g,B,Delta)
lambda1 = min(eig(B));
if lambda1>=0
    pB = -B\g;
    if norm(pB) <= Delta
        p = pB;
        return;
    end
end
lammin  = max(-lambda1, 0) + 1.0E-08;
lambda  = lammin+1;
for i = 1:10,
    R = chol(B+lambda*eye(size(B)));
    p = -R\(R'\g);
    if abs(norm(p) - Delta) <= 1.0E-04*Delta
        return;
    end
    q = R'\p;
    lambda = max(lammin, ...
                 lambda + (norm(p)/norm(q))^2*(norm(p)-Delta)/Delta);
end
% did not converge - return e.g. Cauchy point
p = cauchy(g,B,Delta);

function r=f(x)
r = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;

function r=df(x)
r=[-400*(x(2)-x(1)^2)*x(1) - 2*(1-x(1)); 200*(x(2)-x(1)^2)];

function r=df2(x)
r=[[400*(3*x(1)^2-x(2))+2, -400*x(1)];  [-400*x(1), 200]];