From b29cdc34e352f393337ef00da0b39dee29aea0d0 Mon Sep 17 00:00:00 2001
From: Claudio Maggioni <maggicl@usi.ch>
Date: Mon, 19 Apr 2021 11:23:32 +0200
Subject: [PATCH] hw3: started 1.

---
 Claudio_Maggioni_3/GD.m           | 37 ++++++++++++++++++++++++++++++
 Claudio_Maggioni_3/Newton.m       | 38 +++++++++++++++++++++++++++++++
 Claudio_Maggioni_3/at.m           |  4 ++++
 Claudio_Maggioni_3/backtracking.m |  8 +++++++
 4 files changed, 87 insertions(+)
 create mode 100644 Claudio_Maggioni_3/GD.m
 create mode 100644 Claudio_Maggioni_3/Newton.m
 create mode 100644 Claudio_Maggioni_3/at.m
 create mode 100644 Claudio_Maggioni_3/backtracking.m

diff --git a/Claudio_Maggioni_3/GD.m b/Claudio_Maggioni_3/GD.m
new file mode 100644
index 0000000..6d8a937
--- /dev/null
+++ b/Claudio_Maggioni_3/GD.m
@@ -0,0 +1,37 @@
+syms x y;
+f = (1 - x)^2 + 100 * (y - x^2)^2;
+
+[x, xs, gs] = gd(f, [0;0], 500, 1e-6, true);
+display(x);
+display(xs);
+display(gs);
+plot(xs(1, :), xs(2, :), '-o');
+
+function [xk, xs, gs] = gd(f, x0, max_itr, tol, bt)
+    syms x y;
+    xk = x0;
+    gf = gradient(f);
+    hf = hessian(f, [x, y]);
+    tnorm = norm(at(gf, xk), 2);
+    xs = zeros(size(x0, 1), max_itr + 1);
+    xs(:, 1) = x0;
+    gs = zeros(1, max_itr + 1);
+    gs(1) = tnorm;
+    k = 1;
+    
+    while tnorm > tol && k <= max_itr 
+        p = -at(gf, xk);
+        H = at(hf, xk);
+        if bt
+            alpha = backtracking(f, xk, 1, p, 0.01, 1e-4)
+        else
+            alpha = dot(-H * xk, p) / dot(H * p, p)
+        end
+        xk = xk + alpha * p
+        k = k + 1;
+        xs(:, k) = xk;
+        gs(k) = tnorm;
+    end
+    xs = xs(:, 1:k);
+    gs = gs(:, 1:k);
+end
\ No newline at end of file
diff --git a/Claudio_Maggioni_3/Newton.m b/Claudio_Maggioni_3/Newton.m
new file mode 100644
index 0000000..2e59bd9
--- /dev/null
+++ b/Claudio_Maggioni_3/Newton.m
@@ -0,0 +1,38 @@
+
+
+syms x y;
+f = (1 - x)^2 + 100 * (y - x^2)^2;
+
+[x, xs, gs] = newton(f, [0;0], 50000, 1e-6);
+display(x);
+display(xs);
+display(gs);
+plot(xs(1, :), xs(2, :), '-o');
+
+function y = at(f, xk)
+    syms x y;
+    y = double(subs(f, [x, y], [xk(1), xk(2)]));
+end
+
+function [xk, xs, gs] = newton(f, x0, max_itr, tol)
+    syms x y;
+    xk = x0;
+    gf = gradient(f);
+    hf = hessian(f, [x, y]);
+    tnorm = norm(at(gf, xk), 2);
+    xs = zeros(size(x0, 1), max_itr + 1);
+    xs(:, 1) = x0;
+    gs = zeros(1, max_itr + 1);
+    gs(1) = tnorm;
+    k = 1;
+    
+    while tnorm > tol && k <= max_itr + 1
+        xk = xk - (at(hf, xk) \ at(gf, xk));
+        tnorm = norm(at(gf, xk), 2);
+        k = k + 1;
+        xs(:, k) = xk;
+        gs(k) = tnorm;
+    end
+    xs = xs(:, 1:k);
+    gs = gs(:, 1:k);
+end
\ No newline at end of file
diff --git a/Claudio_Maggioni_3/at.m b/Claudio_Maggioni_3/at.m
new file mode 100644
index 0000000..230a7aa
--- /dev/null
+++ b/Claudio_Maggioni_3/at.m
@@ -0,0 +1,4 @@
+function y = at(f, xk)
+    syms x y;
+    y = double(subs(f, [x, y], [xk(1), xk(2)]));
+end
diff --git a/Claudio_Maggioni_3/backtracking.m b/Claudio_Maggioni_3/backtracking.m
new file mode 100644
index 0000000..adecdf7
--- /dev/null
+++ b/Claudio_Maggioni_3/backtracking.m
@@ -0,0 +1,8 @@
+function alpha = backtracking(f, xk, alpha_bar, pk, rho, c)
+    alpha = alpha_bar;
+    gf = gradient(f);
+    while at(f, xk + alpha * pk) > at(f, xk) - c * alpha * dot(at(gf, xk), pk)
+        alpha = alpha * rho
+    end 
+end
+      
\ No newline at end of file