diff --git a/Claudio_Maggioni_2/Claudio_Maggioni_2.pdf b/Claudio_Maggioni_2/Claudio_Maggioni_2.pdf new file mode 100644 index 0000000..4a61d6c Binary files /dev/null and b/Claudio_Maggioni_2/Claudio_Maggioni_2.pdf differ diff --git a/Claudio_Maggioni_2/Claudio_Maggioni_2.tex b/Claudio_Maggioni_2/Claudio_Maggioni_2.tex new file mode 100755 index 0000000..3c0bba5 --- /dev/null +++ b/Claudio_Maggioni_2/Claudio_Maggioni_2.tex @@ -0,0 +1,122 @@ +\documentclass{scrartcl} +\usepackage[utf8]{inputenc} +\usepackage{graphicx} +\usepackage{subcaption} +\usepackage{amsmath} +\usepackage{pgfplots} +\pgfplotsset{compat=newest} +\usetikzlibrary{plotmarks} +\usetikzlibrary{arrows.meta} +\usepgfplotslibrary{patchplots} +\usepackage{grffile} +\usepackage{amsmath} +\usepackage{subcaption} +\usepgfplotslibrary{external} +\tikzexternalize +\usepackage[margin=2.5cm]{geometry} + +% To compile: +% sed -i 's#title style={font=\\bfseries#title style={yshift=1ex, font=\\tiny\\bfseries#' *.tex +% luatex -enable-write18 -shellescape main.tex + +\pgfplotsset{every x tick label/.append style={font=\tiny, yshift=0.5ex}} +\pgfplotsset{every title/.append style={font=\tiny, align=center}} +\pgfplotsset{every y tick label/.append style={font=\tiny, xshift=0.5ex}} +\pgfplotsset{every z tick label/.append style={font=\tiny, xshift=0.5ex}} + +\setlength{\parindent}{0cm} +\setlength{\parskip}{0.5\baselineskip} + +\title{Optimization methods -- Homework 2} +\author{Claudio Maggioni} + +\begin{document} + +\maketitle + +\section{Exercise 1} + +\subsection{Implement the matrix $A$ and the vector $b$, for the moment, without taking into consideration the +boundary conditions. As you can see, the matrix $A$ is not symmetric. Does an energy function of +the problem exist? Consider $N = 4$ and show your answer, explaining why it can or cannot exist.} + +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. $\frac12 x^T A x - b^T x$ won't work +since $x^T A x$ might be negative and thus the minimizer does not point to +the solution of $Ax = b$ necessairly + +\subsection{Once the new matrix has been derived, write the energy function related to the new problem +and the corresponding gradient and Hessian.} + +we already enforce x(1) = x(n) = 0, since b(1) = b(n) = 0 and thus +A(1, :) * x = b(0) = 0 and same for n can be solved only for x(1) = x(n) += 0size(A, 1) + +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}$ + +\subsection{Write the Conjugate Gradient algorithm in the pdf and implement it Matlab code in a function +called \texttt{CGSolve}.} + +See page 112 (133 for pdf) for the algorithm implementation + +The solution of this task can be found in Section 1.3 of the script \texttt{main.m}. + +\subsection{Solve the Poisson problem.} + +The solution of this task can be found in Section 1.4 of the script \texttt{main.m}. + +\subsection{Plot the value of energy function and the norm of the gradient (here, +use semilogy) as functions of the iterations.} + +The solution of this task can be found in Section 1.5 of the script \texttt{main.m}. + +\subsection{Finally, explain why the Conjugate Gradient method is a Krylov subspace method.} + +Because theorem 5.3 holds, which itself holds mainly because of this (5.10, page 106 [127]): + +\[r_{k+1} = r_k + a_k * A * p_k\] + +\section{Exercise 2} + +Consider the linear system $Ax = b$, where the matrix $A$ is constructed in three different ways: + +\begin{itemize} + \item $A =$ diag([1:10]) + \item $A =$ diag(ones(1,10)) + \item $A =$ diag([1, 1, 1, 3, 4, 5, 5, 5, 10, 10]) + \item $A =$ diag([1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]) +\end{itemize} + +\subsection{How many distinct eigenvalues has each matrix?} + +Each matrix has a distinct number of eigenvalues equal to the number of distinct +elements on its diagonal. So, in order, each A has respectively 10, 1, 5, and 10 distinct eigenvalues. + +\subsection{Construct a right-hand side $b=$rand(10,1) and apply the +Conjugate Gradient method to solve the system for each $A$.} + +The solution of this task can be found in section 2.2 of the \texttt{main.m} MATLAB script. + +\subsection{Compute the logarithm energy norm of the error for each matrix +and plot it with respect to the number of iteration.} + +The solution of this task can be found in section 2.3 of the \texttt{main.m} MATLAB script. + +\subsection{Comment on the convergence of the method for the different matrices. What can you say observing +the number of iterations obtained and the number of clusters of the eigenvalues of the related +matrix?} + +The method converges quickly for each matrix. The fastest convergence surely happens for $A2$, which is +the identity matrix and therefore makes the $Ax = b$ problem trivial. + +For all the other matrices, we observe the energy norm of the error decreasing exponentially as the iterations +increase, eventually reaching $0$ for the cases where the method converges exactly (namely on matrices $A1$ and $A3$). + +Other than for the fourth matrix, the number of iterations is exactly equal +to the number of distinct eigenvalues for the matrix. That exception on the fourth matrix is simply due to the +tolerance termination condition holding true for an earlier iteration, i.e. we terminate early since we find an +approximation of $x$ with residual norm below $10^{-8}$. + +\end{document} diff --git a/Claudio_Maggioni_2/ex1.m b/Claudio_Maggioni_2/ex1.m deleted file mode 100644 index 6703782..0000000 --- a/Claudio_Maggioni_2/ex1.m +++ /dev/null @@ -1,100 +0,0 @@ -%% 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 - -%% 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 - -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; - end - - h = 1 / (n - 1); - - b = h^2 * ones(n, 1); - b(1) = 0; - b(n) = 0; -end - -%% 1.2 - -% we already enforce x(1) = x(n) = 0, since b(1) = b(n) = 0 and thus -% A(1, :) * x = b(0) = 0 and same for n can be solved only for x(1) = x(n) -% = 0 -% -% 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 diff --git a/Claudio_Maggioni_2/main.m b/Claudio_Maggioni_2/main.m new file mode 100644 index 0000000..2bff263 --- /dev/null +++ b/Claudio_Maggioni_2/main.m @@ -0,0 +1,132 @@ +%% Homework 2 - Optimization Methods +% Author: Claudio Maggioni +% +% Note: exercises are not in the right order due to matlab constraints of +% functions inside of scripts. + +clear +clc +close all +plots = 1; + +%% 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), n, 1e-8); +display(x); + +%% 1.5 - Plots for the 1D Poisson solution + +if plots + plot(0:(size(ys,1)-1), ys); + sgtitle("Objective function values per iteration"); + axis([-1 500 -inf +inf]); + + figure; + semilogy(0:(size(gnorms,1)-1), gnorms); + sgtitle("Log of gradient norm per iteration"); + axis([-1 500 -inf +inf]); +end +%% 2.1 - Matrix definitions + +A1 = diag([1:10]); +A2 = diag(ones(1,10)); +A3 = diag([1, 1, 1, 3, 4, 5, 5, 5, 10, 10]); +A4 = diag([1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]); + +%% 2.2 - Application of CG for each AX/b couple + +rng(0); % random seed fixed due to reproducibility purposes +b = rand(10,1); +display(b); + +n = 10; +[x1, ys1, gnorms1, xs1] = CGSolve(A1, b, zeros(n,1), n, 1e-8); +[x2, ys2, gnorms2, xs2] = CGSolve(A2, b, zeros(n,1), n, 1e-8); +[x3, ys3, gnorms3, xs3] = CGSolve(A3, b, zeros(n,1), n, 1e-8); +[x4, ys4, gnorms4, xs4] = CGSolve(A4, b, zeros(n,1), n, 1e-8); + +%% 2.3 - Logarithm energy norm of the error computation + +if plots + enl_plot(x1, xs1, A1); + sgtitle("Log energy norm of the error per iter. (matrix A1)"); + enl_plot(x2, xs2, A2); + sgtitle("Log energy norm of the error per iter. (matrix A2)"); + enl_plot(x3, xs3, A3) + sgtitle("Log energy norm of the error per iter. (matrix A3)"); + enl_plot(x4, xs4, A4); + sgtitle("Log energy norm of the error per iter. (matrix A4)"); +end + +function enl_plot(xsol, xs, A) + enls = zeros(size(xs, 2), 1); + for i = 1:size(xs, 2) + x = xs(:, i); + enls(i) = log((xsol - x)' * A * (xsol - x)); + end + + figure; + plot(0:(size(xs, 2)-1), enls, '-k.'); + axis([-1 11 -35 +2]); +end + +%% 1.1 - Build the Poisson matrix A and vector b + +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; + end + + h = 1 / (n - 1); + + b = h^2 * ones(n, 1); + b(1) = 0; + b(n) = 0; +end + +%% 1.3 - Implementation of Conjugate Gradient + +function [x, ys, gnorms, xs] = CGSolve(A, b, x, max_itr, tol) + ys = zeros(max_itr + 1, 1); + gnorms = zeros(max_itr + 1, 1); + xs = zeros(size(x, 1), max_itr + 1); + + r = A * x - b; + p = -r; + + gnorms(1) = norm(p, 2); + xs(:, 1) = x; + ys(1) = 1/2 * dot(A*x, x) - dot(b, x); + + for i = 1:(max_itr+1) + + alpha = -r' * p / (p' * A * p); + x = x + alpha * p; + r = A * x - b; + beta = (r' * A * p) / (p' * A * p); + p = -r + beta * p; + + gnorms(i+1) = norm(p, 2); + xs(:, i+1) = x; + ys(i+1) = 1/2 * dot(A*x, x) - dot(b, x); + + if gnorms(i+1) < tol + ys = ys(1:(i+1)); + gnorms = gnorms(1:(i+1)); + xs = xs(:, 1:(i+1)); + break + end + end +end