OM/Claudio_Maggioni_midterm/trust_region.m

68 lines
1.9 KiB
Matlab

function [xk, xs, gnorms] = trust_region(f, delta_hat, delta0, eta, ...
x0, tol, max_n)
xs = zeros(2, max_n);
gnorms = zeros(max_n);
xk = x0;
deltak = delta0;
fl = vecLambda(matlabFunction(f));
gl = vecLambda(matlabFunction(gradient(f)));
hl = vecLambda(matlabFunction(hessian(f)));
xs(:, 1) = x0;
for i = 2:max_n
fk = fl(xk);
B = hl(xk);
g = gl(xk);
gnorms(i - 1) = norm(g);
if gnorms(i - 1) < tol
gnorms = gnorms(1:i-1);
xs = xs(:, 1:i-1);
break
end
pk = dogleg(B, g, deltak);
rho_k = (fl(xk) - fl(xk + pk)) / ...
(qf(B, g, [0;0], fk) - qf(B, g, pk, fk));
if rho_k < 1/4
deltak = 1/4 * deltak;
% When comparing the method's execution with some classmates, we
% found some numerical instability in the comparison between the
% step's norm and the trust region radius. To sum up, using
% different Matlab versions (2020b vs 2019a) we obtained different
% results on that comparison (true vs. false) on seemingly
% identical values. It seems there is some difference w.r.t.
% comparisons in the unnormalized double range, and therefore we
% both approximated that equality with the subtraction you will find
% below
elseif rho_k > 3/4 && (norm(pk, 2) - deltak) < eps
deltak = min(2 * deltak, delta_hat);
end
% otherwhise do not change delta
if rho_k > eta
xk = xk + pk;
end
% otherwise do not change xk
xs(:, i) = xk;
end
end
% Convert lambda to accept vector parameters
function vl = vecLambda(fl)
vl = @(x) fl(x(1), x(2));
end
% Compute quadratic form
function y = qf(B, g, p, fk)
y = fk + 1/2 * p' * B * p + dot(g, p);
end