From 3f8034123ef44b795f8eae1a00de020efc7bcc37 Mon Sep 17 00:00:00 2001
From: "Claudio Maggioni (maggicl)" <maggicl@kolabnow.ch>
Date: Mon, 3 May 2021 16:32:30 +0200
Subject: [PATCH] midterm: hw2 matlab mostly done

---
 Claudio_Maggioni_midterm/cauchy.m       | 12 +++++
 Claudio_Maggioni_midterm/dogleg.m       | 20 +++++++
 Claudio_Maggioni_midterm/trust_region.m | 71 +++++++++++++++++++++++++
 3 files changed, 103 insertions(+)
 create mode 100644 Claudio_Maggioni_midterm/cauchy.m
 create mode 100644 Claudio_Maggioni_midterm/dogleg.m
 create mode 100644 Claudio_Maggioni_midterm/trust_region.m

diff --git a/Claudio_Maggioni_midterm/cauchy.m b/Claudio_Maggioni_midterm/cauchy.m
new file mode 100644
index 0000000..3a336a5
--- /dev/null
+++ b/Claudio_Maggioni_midterm/cauchy.m
@@ -0,0 +1,12 @@
+function pk = cauchy(B, g, deltak)
+   gbg = (g' * B * g);
+   
+   if gbg <= 0
+       tau = 1;
+   else
+       tau = min(norm(g, 2)^3 / (deltak * gbg), 1); 
+   end
+   
+   pk = -tau * deltak / norm(g, 2) * g;
+end
+
diff --git a/Claudio_Maggioni_midterm/dogleg.m b/Claudio_Maggioni_midterm/dogleg.m
new file mode 100644
index 0000000..35df6bf
--- /dev/null
+++ b/Claudio_Maggioni_midterm/dogleg.m
@@ -0,0 +1,20 @@
+function pk = dogleg(B, g, deltak) 
+    pnewton = - (B \ g);
+    
+    if norm(pnewton, 2) <= deltak
+        pk = pnewton;
+    else
+        pu = - dot(g, g) / (g' * B * g) * g;
+        
+        if norm(pu, 2) > deltak
+            pk = cauchy(B, g, deltak);
+        else
+            syms taux
+            eqn = norm(pu + taux * (pnewton - pu))^2 == deltak^2;
+            tau = solve(eqn, taux);
+            tau = double(max(tau));
+            pk = pu + tau * (pnewton - pu);
+        end
+    end
+end
+
diff --git a/Claudio_Maggioni_midterm/trust_region.m b/Claudio_Maggioni_midterm/trust_region.m
new file mode 100644
index 0000000..09f22c4
--- /dev/null
+++ b/Claudio_Maggioni_midterm/trust_region.m
@@ -0,0 +1,71 @@
+
+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