hw2: Done 1.1-5

Claudio_Maggioni_2/ex1.m

@@ -1,24 +1,52 @@
 %% Homework 2 - Optimization Methods
 % Author: Claudio Maggioni
 %
+% Note: exercises are not in the right order due to matlab constraints of
+% functions inside of scripts.
+%
 % Sources:
 
 clear
 clc
 close all
 
-[A, b] = build_poisson(4);
-disp(A)
-disp(b)
+%% 1.4 - Solution for 1D Poisson for N=1000 using CG
+
+n = 1000;
+[A, b] = build_poisson(n);
+
+A(2, 1) = 0;
+A(n-1, n) = 0;
+
+[x, ys, gnorms] = CGSolve(A, b, zeros(n,1), 1000, 1e-8);
+display(x);
+
+%% 1.5 - Plots for the 1D Poisson solution
+
+plot(ys);
+sgtitle("Objective function values per iteration");
+
+figure;
+semilogy(gnorms);
+sgtitle("Log of gradient norm per iteration");
+
+
+
+
+%% 1.1 - Build the Poisson matrix A and vector b (check this)
+
+% Answer is a energy function does not exist. Since A is not symmetric
+% (even if it is pd), the minimizer used for the c.g. method
+% (i.e. (1/2)x^TAx - b^x) won't work
+% since x^TAx might be negative and thus the minimizer does not point to
+% the solution of Ax=B necessairly
 
-% 1.1
 function [A,b] = build_poisson(n)
     A = diag(2 * ones(1,n));
     
     A(1,1) = 1;
     A(n,n) = 1;
     
-    
     for i = 2:n-1
         A(i, i+1) = -1;
         A(i, i-1) = -1;
@@ -31,14 +59,6 @@ function [A,b] = build_poisson(n)
     b(n) = 0;
 end
 
-%% 1.1 (check this)
-
-% Answer is a energy function does not exist. Since A is not symmetric
-% (even if it is pd), the minimizer used for the c.g. method
-% (i.e. (1/2)x^TAx - b^x) won't work
-% since x^TAx might be negative and thus the minimizer does not point to
-% the solution of Ax=B necessairly
-
 %% 1.2
 
 % we already enforce x(1) = x(n) = 0, since b(1) = b(n) = 0 and thus
@@ -48,3 +68,33 @@ end
 % The objective is therefore \phi(x) = (1/2)x^T\overline{A}x - b^x with a and b
 % defined above, gradient is = \overline{A}x - b, hessian is = \overline{A} 
 
+%% 1.3 - Implementation of Conjugate Gradient
+
+function [x, ys, gnorms] = CGSolve(A, b, x, max_itr, tol)
+    ys = zeros(max_itr, 1);
+    gnorms = zeros(max_itr, 1);
+    
+    r = b - A * x;
+    d = r;
+    delta_old = dot(r, r);
+    
+    for i = 1:max_itr
+        ys(i) = 1/2 * dot(A*x, x) - dot(b, x);
+        gnorms(i) = sqrt(delta_old);
+        
+        s = A * d;
+        alpha = delta_old / dot(d, s);
+        x = x + alpha * d;
+        r = r - alpha * s;
+        delta_new = dot(r, r);
+        beta = delta_new / delta_old;
+        d = r + beta * d;
+        delta_old = delta_new;
+        
+        if delta_new / norm(b) < tol
+            ys = [ys(1:i); 1/2 * dot(A*x, x) - dot(b, x)];
+            gnorms = [gnorms(1:i); sqrt(delta_old)];
+            break
+        end
+    end
+end