diff --git a/Claudio_Maggioni_2/ex1.m b/Claudio_Maggioni_2/ex1.m index 8a9d5b5..6703782 100644 --- a/Claudio_Maggioni_2/ex1.m +++ b/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 \ No newline at end of file