72 lines
2 KiB
Mathematica
72 lines
2 KiB
Mathematica
|
|
||
|
|
syms x y
|
||
|
|
f1 = (y - 4 * x^2)^2 + (1 - x)^2;
|
||
|
|
[x, xs, gnorms] = trust_reg(f1, 2, 1, 0.2, [0;0], 1e-8, 1000);
|
||
|
|
|
||
|
|
% 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
|
||
|
|
|
||
|
|
function [xk, xs, gnorms] = trust_reg(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
|