mp5: done 1-3

mp5/A_eig.tex

@@ -0,0 +1,126 @@
+% This file was created by matlab2tikz.
+%
+\definecolor{mycolor1}{rgb}{0.00000,0.44700,0.74100}%
+%
+\begin{tikzpicture}
+
+\begin{axis}[%
+width=6.028in,
+height=4.754in,
+at={(1.011in,0.642in)},
+scale only axis,
+xmin=0,
+xmax=100,
+ymode=log,
+ymin=0.001,
+ymax=10000,
+yminorticks=true,
+axis background/.style={fill=white},
+title style={font=\bfseries},
+title={Eigenvalues of A}
+]
+\addplot [color=mycolor1, forget plot]
+  table[row sep=crcr]{%
+1	0.00183340367988759\\
+2	0.00453790786113484\\
+3	0.0107896856046935\\
+4	0.0211235851960644\\
+5	0.0305491923444261\\
+6	0.0424763662116699\\
+7	0.0652345604159706\\
+8	0.102757346585257\\
+9	0.153538421986227\\
+10	0.211241943436195\\
+11	0.252211715498981\\
+12	0.289989070850436\\
+13	0.351131851500607\\
+14	0.435706363023172\\
+15	0.475701188846017\\
+16	0.536773942313087\\
+17	0.570786145601508\\
+18	0.651735983378915\\
+19	0.727246558118914\\
+20	0.787465786548366\\
+21	0.876077658577838\\
+22	0.959634982712415\\
+23	1.11272341205416\\
+24	1.16964510997899\\
+25	1.19317301472187\\
+26	1.34288169306462\\
+27	1.48612799031325\\
+28	1.56948387062071\\
+29	1.78341752583741\\
+30	1.93501326653969\\
+31	1.94944936154279\\
+32	2.0834385818203\\
+33	2.25830482514161\\
+34	2.3104130395112\\
+35	2.44838431259137\\
+36	2.78327385790993\\
+37	2.90201594529841\\
+38	3.13205713047726\\
+39	3.25341826203778\\
+40	3.35818273116376\\
+41	3.70494882713741\\
+42	3.85238780788471\\
+43	3.9846975992715\\
+44	4.09652512284048\\
+45	4.26425186014381\\
+46	4.42615565949843\\
+47	4.64205273607693\\
+48	4.83516700026364\\
+49	4.89609731947842\\
+50	5.14541243767704\\
+51	5.47177991705298\\
+52	5.83426853647913\\
+53	6.24919735017955\\
+54	6.51069828517495\\
+55	6.76926716800134\\
+56	6.96357007137143\\
+57	7.22399026534903\\
+58	7.45959355109047\\
+59	7.6215883910995\\
+60	8.13211718148974\\
+61	8.2348761782867\\
+62	9.02938692612482\\
+63	9.26003297236477\\
+64	9.59443602091006\\
+65	9.78230272154912\\
+66	9.9397509124274\\
+67	10.8200374064158\\
+68	11.2239205229841\\
+69	11.5014660930885\\
+70	11.8095643578441\\
+71	11.9722770578147\\
+72	12.1797060122863\\
+73	12.3294752069053\\
+74	13.5048869251576\\
+75	13.6918740486823\\
+76	14.0726278232322\\
+77	14.6109572787373\\
+78	14.8717535190965\\
+79	15.5900717137855\\
+80	16.4268311711295\\
+81	16.8632629397018\\
+82	17.1804683132309\\
+83	17.5656109813257\\
+84	18.4405918720105\\
+85	18.9256543387976\\
+86	19.925816210275\\
+87	20.5420086066435\\
+88	20.8431836934263\\
+89	21.2272436238988\\
+90	22.028625968725\\
+91	22.7106432954137\\
+92	23.6099490290742\\
+93	24.9559178747221\\
+94	26.0355909133376\\
+95	26.4089094809298\\
+96	27.6083096450885\\
+97	28.4262392990311\\
+98	29.2220595414911\\
+99	30.9554105903231\\
+100	2511.69406680491\\
+};
+\end{axis}
+\end{tikzpicture}%

mp5/Makefile

@@ -0,0 +1,15 @@
+filename=template
+
+pdf:
+	pdflatex ${filename}
+	#bibtex ${filename}
+	pdflatex ${filename}
+	pdflatex ${filename}
+	make clean
+
+read:
+	evince ${filename}.pdf &
+
+
+clean:
+	rm -f *.out *.log *.bbl *.blg *.aux ${filename}.log  ${filename}.ps ${filename}.aux ${filename}.out ${filename}.dvi ${filename}.bbl ${filename}.blg

mp5/Project_5_Maggioni_Claudio.tex

@@ -0,0 +1,126 @@
+\documentclass[unicode,11pt,a4paper,oneside,numbers=endperiod,openany]{scrartcl}
+
+\usepackage{graphicx}
+\usepackage{subcaption}
+\usepackage{amsmath}
+\input{assignment.sty}
+\usepackage{pgfplots}
+\pgfplotsset{compat=newest}
+\usetikzlibrary{plotmarks}
+\usetikzlibrary{arrows.meta}
+\usepgfplotslibrary{patchplots}
+\usepackage{grffile}
+\usepackage{amsmath}
+\usepackage{subcaption}
+\usepgfplotslibrary{external}
+\tikzexternalize
+
+\begin{document}
+
+\setassignment
+\setduedate{Wednesday, December 02, 2020, 11:59 PM}
+
+\serieheader{Numerical Computing}{2020}{Student: Claudio Maggioni}{Discussed with: --}{Solution for Project 5}{}
+\newline
+
+\assignmentpolicy
+
+The purpose of this assignment is to gain insight on the theoretical and
+numerical properties of the Conjugate Gradient method. Here we use this method
+in an image processing application with the goal of deblur an image given the
+exact (noise-free) blurred image and the original transformation matrix. Note
+that the ``noise-free" simplification is essential for us to solve this problem
+in the scope of this assignment.
+
+\section{General Questions [10 points]}
+\subsection{What is the size of the matrix $A$?}
+$A$ is an $n^2$ by $n^2$ matrix, where $n$ is the width and height in pixels
+of the image to transform.
+
+\subsection{How many diagonals bands does $A$ have?}
+$A$ as $d^2$ diagonal bands, where $d$ is strictly an order of magnitude below
+$n$.
+
+\subsection{What is the length of the vectorized blur image $b$?}
+$b$ is the row-vectorized form of the image pixel matrix, and thus has
+dimensions 1 by $n^2$.
+
+\section{Properties of A [10 points]}
+\subsection{ If $A$ is not symmetric, how would this affect $\tilde{A}$?}
+If $A$ were not symmetric, then $\tilde{A}$ would not be positive definite since
+by definition $\tilde{A} = A A^T$, thus not satifying the assumptions taken when
+solving the system.
+
+\subsection{ Explain why solving $Ax = b$ for $x$ is equivalent to
+minimizing $x^T A x - b^T x$ over $x$.}
+
+First, we can say that:
+
+\[f(x) = \frac12 x^T A x - b^T x = \frac12\langle Ax,x \rangle -
+\langle b,x \rangle\]
+
+Then by taking the derivative of $f(x)$ w.r.t. $x$ we have (assuming $A$ is
+\textit{spd}, which it is):
+
+\[f'(x) = \frac12 A^T x + \frac12Ax - b = Ax - b\]
+
+which for $f'(x) = 0$ will be equivalent to solving $Ax = b$. By taking the
+second derivative we have:
+
+\[f''(x) = A > 0\]
+
+since $A$ is positive definite. Therefore, we can say that the absolute minima
+of $f(x)$ is the solution for $Ax = b$.
+
+
+\section{Conjugate Gradient [40 points]}
+\subsection{ Write a function for the conjugate gradient solver
+\texttt{[x,rvec]=myCG(A,b,x0,max\_itr,tol)}, where \texttt{x}
+and \texttt{rvec} are, respectively, the solution value and a
+vector containing the residual at every iteration.}
+
+The implementation can be found in file \texttt{myCG.m} in the source directory.
+The test code for the function \texttt{myCG} can be found in the \texttt{test.m} file.
+
+\subsection{ In order to validate your implementation, solve the system
+defined by \texttt{A\_test.mat} and \texttt{b\_test.mat}. Plot
+the convergence (residual vs iteration).}
+
+The plot of the squared residual 2-norms over all iterations can be found in Figure
+\ref{fig:plot1}.
+
+\begin{figure}[h]
+	\input{test_semilogy}
+	\caption{Semilog plot of the plot of the squared residual 2-norms over all iterations}\label{fig:plot1}
+\end{figure}
+
+\subsection{ Plot the eigenvalues of \texttt{A\_test.mat} and comment on the
+condition number and convergence rate.}
+
+The eigenvalues of A can be found in figure
+\ref{fig:plot2}. The condition number for matrix $A$ according to \texttt{rcond(...)} is $\approx 3.2720 \cdot 10^7$, which is very low without sitting in the denormalized range (i.e. $< \text{eps}$) and thus very good for the Conjugate Gradient algorithm.
+This well conditioning is also reflected in the eigenvalue plot, which shows a not so
+drastic increase of the first eigenvalues ordered in increasing order.
+
+\begin{figure}[h]
+	\input{A_eig}
+	\caption{Semilog plot of the eigenvalues of A}\label{fig:plot2}
+\end{figure}
+
+\section{Debluring problem [40 points]}
+\subsection{ Solve the debluring problem for the blurred image matrix
+\texttt{B.mat} and transformation matrix \texttt{A.mat} using
+your routine \texttt{myCG} and Matlab's preconditioned conjugate
+gradient \texttt{pcg}. As a preconditioner, use \texttt{ichol}
+to get the incomplete Cholesky factors and set routine type to
+\texttt{nofill} with $\alpha=0.01$ for the diagonal shift (see
+Matlab documentation). Solve the system with both solvers using
+$max\_iter=200$ $tol= 10^{-6}$. Plot the convergence (residual
+vs iteration) of each solver and display the original and final
+deblurred image.}
+\subsection{ When would \texttt{pcg} be worth the added computational cost?
+What about if you are debluring lots of images with the same
+blur operator?}
+
+
+\end{document}

mp5/Project_5_Maggioni_Claudio/code_template.m

@@ -0,0 +1,29 @@
+close all;
+clear; clc;
+
+%% Load Default Img Data
+load('blur_data/B.mat');
+B=double(B);
+
+% Show Image
+figure
+im_l=min(min(B));
+im_u=max(max(B));
+imshow(B,[im_l,im_u])
+title('Blured Image')
+
+% Vectorize the image (row by row)
+b=B';
+b=b(:);
+
+
+
+%% Validate Test values
+load('test_data/A_test.mat');
+load('test_data/x_test_exact.mat');
+load('test_data/b_test.mat');
+
+%res=||x^*-A^{-1}b||
+res=x_test_exact-inv(A_test)*b_test
+norm(res)
+%(Now do it with your CG and Matlab's PCG routine!!!)

mp5/Project_5_Maggioni_Claudio/myCG.m

@@ -0,0 +1,26 @@
+function [x,rvec] = myCG(A, b, x0, max_itr, tol)
+    x = x0;
+    rvec = zeros(1,max_itr+1);
+    r = b - A * x0;
+    rvec(1) = norm(r, 2);
+    d = r;
+    delta_old = dot(r, r);
+    
+    for i = 1:max_itr
+        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;
+        rvec(i + 1) = delta_new; 
+        
+        if delta_new < tol
+            rvec = rvec(1:i+1);
+            break
+        end
+    end
+end
+

mp5/Project_5_Maggioni_Claudio/test.m

@@ -0,0 +1,18 @@
+addpath /Users/maggicl/Git/matlab2tikz/src/;
+
+load('test_data/A_test.mat');
+load('test_data/b_test.mat');
+load('test_data/x_test_exact.mat');
+
+[xa, rveca] = myCG(A_test, b_test, b_test * 0, 1000, cutOff);
+
+figure;
+semilogy(rveca);
+title("Residual vector squared 2-norm (log) over iterations");
+matlab2tikz('showInfo', false, '../test_semilogy.tex');
+
+v = eig(A_test);
+figure;
+semilogy(v);
+title("Eigenvalues of A");
+matlab2tikz('showInfo', false, '../A_eig.tex');

mp5/assignment.sty

@@ -0,0 +1,95 @@
+\usepackage{ifthen}
+\usepackage[utf8]{inputenc}
+\usepackage{graphics}
+\usepackage{graphicx}
+\usepackage{hyperref}
+
+\pagestyle{plain}
+\voffset -5mm
+\oddsidemargin  0mm
+\evensidemargin -11mm
+\marginparwidth 2cm
+\marginparsep 0pt
+\topmargin 0mm
+\headheight 0pt
+\headsep 0pt
+\topskip 0pt        
+\textheight 255mm
+\textwidth 165mm
+
+\newcommand{\duedate} {}
+\newcommand{\setduedate}[1]{%
+\renewcommand\duedate {Due date:~ #1}}
+\newcommand\isassignment {false}
+\newcommand{\setassignment}{\renewcommand\isassignment {true}}
+\newcommand{\ifassignment}[1]{\ifthenelse{\boolean{\isassignment}}{#1}{}}
+\newcommand{\ifnotassignment}[1]{\ifthenelse{\boolean{\isassignment}}{}{#1}}
+
+\newcommand{\assignmentpolicy}{
+\begin{table}[h]
+\begin{center}
+\scalebox{0.8} {%
+\begin{tabular}{|p{0.02cm}p{16cm}|}
+\hline
+&\\
+\multicolumn{2}{|c|}{\Large\textbf{Numerical Computing 2020 ---  Submission Instructions}}\\
+\multicolumn{2}{|c|}{\large\textbf{(Please, notice that following instructions are mandatory: }}\\
+\multicolumn{2}{|c|}{\large\textbf{submissions that don't comply with, won't be considered)}}\\
+&\\
+\textbullet & Assignments must be submitted to \href{https://www.icorsi.ch/course/view.php?id=10018}{iCorsi} (i.e. in electronic format).\\
+\textbullet & Provide both executable package and sources (e.g. C/C++ files, Matlab). 
+If you are using libraries, please add them in the file. Sources must be organized in directories called:\\
+\multicolumn{2}{|c|}{\textit{Project\_number\_lastname\_firstname}}\\
+& and  the  file must be called:\\
+\multicolumn{2}{|c|}{\textit{project\_number\_lastname\_firstname.zip}}\\
+\multicolumn{2}{|c|}{\textit{project\_number\_lastname\_firstname.pdf}}\\
+\textbullet &  The TAs will grade your project by reviewing your project write-up, and looking at the implementation 
+                 you attempted, and benchmarking your code's performance.\\
+
+\textbullet & You are allowed to discuss all questions with anyone you like; however: (i) your submission must list anyone you discussed problems with and (ii) you must write up your submission independently.\\
+\hline
+\end{tabular}
+}
+\end{center}
+\end{table}
+}
+\newcommand{\punkte}[1]{\hspace{1ex}\emph{\mdseries\hfill(#1~\ifcase#1{Points}\or{Points}\else{Points}\fi)}}
+
+
+\newcommand\serieheader[6]{
+\thispagestyle{empty}%
+\begin{flushleft}
+\includegraphics[width=0.4\textwidth]{usi_inf}
+\end{flushleft}
+  \noindent%
+  {\large\ignorespaces{\textbf{#1}}\hspace{\fill}\ignorespaces{ \textbf{#2}}}\\ \\%
+  {\large\ignorespaces #3 \hspace{\fill}\ignorespaces #4}\\
+  \noindent%
+  \bigskip
+  \hrule\par\bigskip\noindent%
+  \bigskip {\ignorespaces {\Large{\textbf{#5}}}
+  \hspace{\fill}\ignorespaces \large \ifthenelse{\boolean{\isassignment}}{\duedate}{#6}}
+  \hrule\par\bigskip\noindent%  \linebreak
+ }
+
+\makeatletter
+\def\enumerateMod{\ifnum \@enumdepth >3 \@toodeep\else
+      \advance\@enumdepth \@ne
+      \edef\@enumctr{enum\romannumeral\the\@enumdepth}\list
+      {\csname label\@enumctr\endcsname}{\usecounter
+        {\@enumctr}%%%? the following differs from "enumerate"
+	\topsep0pt%
+	\partopsep0pt%
+	\itemsep0pt%
+	\def\makelabel##1{\hss\llap{##1}}}\fi}
+\let\endenumerateMod =\endlist
+\makeatother
+
+
+
+
+\usepackage{textcomp}
+
+
+
+

mp5/test_semilogy.tex

@@ -0,0 +1,205 @@
+% This file was created by matlab2tikz.
+%
+\definecolor{mycolor1}{rgb}{0.00000,0.44700,0.74100}%
+%
+\begin{tikzpicture}
+
+\begin{axis}[%
+width=6.028in,
+height=4.754in,
+at={(1.011in,0.642in)},
+scale only axis,
+xmin=0,
+xmax=180,
+ymode=log,
+ymin=1e-12,
+ymax=1000000,
+yminorticks=true,
+axis background/.style={fill=white},
+title style={font=\bfseries},
+title={Residual vector squared 2-norm (log) over iterations}
+]
+\addplot [color=mycolor1, forget plot]
+  table[row sep=crcr]{%
+1	11909.9916990842\\
+2	1215.89010103615\\
+3	119.561506273628\\
+4	25.7035634544525\\
+5	10.5611870908368\\
+6	4.9270658032822\\
+7	2.13962279549106\\
+8	123.990459337534\\
+9	1.23020970291512\\
+10	0.726292742104297\\
+11	0.437655877749282\\
+12	0.260873707832738\\
+13	0.17072420219991\\
+14	0.101340305340657\\
+15	0.189833128810975\\
+16	0.131825297283515\\
+17	0.0572076044770027\\
+18	0.039668263144492\\
+19	0.0319287317681079\\
+20	0.0195740744402844\\
+21	0.0149039165591494\\
+22	0.0123032827360602\\
+23	0.712044976614692\\
+24	0.00995641116934014\\
+25	0.00893611054884436\\
+26	0.00630098046101872\\
+27	0.00529335587330463\\
+28	0.00589902032947632\\
+29	0.00588515046013994\\
+30	0.406172316998164\\
+31	0.0040950751378787\\
+32	0.00373404399719963\\
+33	0.0035997277309428\\
+34	0.00299194637307639\\
+35	0.00177328859853697\\
+36	0.00129725288116854\\
+37	0.049535388437175\\
+38	0.0011091440769831\\
+39	0.00118263501318397\\
+40	0.000966559105586302\\
+41	0.000718169185559145\\
+42	0.000724441141532492\\
+43	0.000768917036778476\\
+44	0.00269020183225918\\
+45	0.000570706808899611\\
+46	0.000831736774154204\\
+47	0.000649300990333462\\
+48	0.000310381139613716\\
+49	0.000235403925330503\\
+50	0.000340608947493281\\
+51	0.001318689210481\\
+52	0.000205980081464339\\
+53	0.000279777304655239\\
+54	0.000109664496511463\\
+55	0.000141033277838556\\
+56	0.000112654872491888\\
+57	0.000165756718289282\\
+58	0.00084800673419117\\
+59	0.000130276460617259\\
+60	0.000203034526747493\\
+61	0.000176255235921362\\
+62	9.4197961691087e-05\\
+63	8.20523976696223e-05\\
+64	0.00370668271425266\\
+65	9.93383601592671e-05\\
+66	0.000153079242536576\\
+67	9.14044749794555e-05\\
+68	5.72233807063173e-05\\
+69	6.34922495601101e-05\\
+70	7.27636976386079e-05\\
+71	0.00225428493244681\\
+72	0.000113799191067991\\
+73	5.35440047147865e-05\\
+74	8.41355153412546e-05\\
+75	8.48518650076496e-05\\
+76	6.22215998223378e-05\\
+77	0.000817207214704047\\
+78	0.000111352409122667\\
+79	4.21172150194178e-05\\
+80	2.87734653228315e-05\\
+81	4.10050601226006e-05\\
+82	3.24394526517843e-05\\
+83	4.16294986387063e-05\\
+84	0.00214569803932431\\
+85	8.55723713819309e-05\\
+86	3.32557658527889e-05\\
+87	2.52240728364592e-05\\
+88	6.2097593050376e-05\\
+89	2.18185972524218e-05\\
+90	0.000295303841112173\\
+91	4.87144579957892e-05\\
+92	2.8215056148989e-05\\
+93	2.35361887050575e-05\\
+94	2.75041794828811e-05\\
+95	5.06188272864469e-05\\
+96	2.42061490399217e-05\\
+97	0.00027580878447057\\
+98	2.74664478340053e-05\\
+99	3.50733084924626e-05\\
+100	3.1778479688316e-05\\
+101	2.92589989035591e-05\\
+102	0.000129072395418128\\
+103	0.0047343379748716\\
+104	4.09472473509763e-05\\
+105	4.10312409025682e-05\\
+106	4.76401831045035e-05\\
+107	4.27117411441813e-05\\
+108	0.000157541717376013\\
+109	0.000168844397829638\\
+110	7.4065189988145e-05\\
+111	0.00022062018323934\\
+112	7.88776820700891e-05\\
+113	0.000116991745030062\\
+114	9.06776322106616e-05\\
+115	4.3360143026987e-05\\
+116	0.00329403333122333\\
+117	4.76599036875217e-05\\
+118	4.10940764030938e-05\\
+119	3.58177666192589e-05\\
+120	7.89625389146185e-06\\
+121	1.61751062518811e-05\\
+122	0.000441952575307103\\
+123	2.4016614151962e-05\\
+124	2.73115201713094e-05\\
+125	3.22102681494731e-05\\
+126	2.08196860675751e-05\\
+127	3.39310585133428e-05\\
+128	2.35747264136633e-05\\
+129	0.000318952227696151\\
+130	9.37353712503354e-05\\
+131	4.2670279781618e-05\\
+132	8.46751510649503e-05\\
+133	2.41046768091548e-05\\
+134	0.000118220396031136\\
+135	0.00984156577501441\\
+136	7.3484752201193e-05\\
+137	4.51598504488994e-05\\
+138	4.51515273101713e-05\\
+139	1.16508810374033e-05\\
+140	1.03229247762991e-05\\
+141	0.000629284559580057\\
+142	2.10831113801501e-06\\
+143	8.64626377047462e-06\\
+144	1.90060956957351e-06\\
+145	2.9904838823846e-06\\
+146	1.1698091976711e-05\\
+147	0.00157994086647156\\
+148	3.90577499999355e-06\\
+149	3.63123271261456e-05\\
+150	1.13464425483349e-05\\
+151	6.23630758043014e-06\\
+152	5.33727229652419e-05\\
+153	0.00766964489351792\\
+154	4.40019782836122e-05\\
+155	4.91424078393586e-06\\
+156	5.23015596899946e-06\\
+157	1.59109137910169e-05\\
+158	2.44921817021047e-06\\
+159	8.92395259274082e-05\\
+160	1.65115695649721e-06\\
+161	2.09460062061473e-06\\
+162	5.05915189772462e-07\\
+163	1.32841240174248e-07\\
+164	1.28360359603615e-06\\
+165	1.43607180361307e-05\\
+166	5.04094734598852e-06\\
+167	7.41197749055341e-06\\
+168	5.04948498688418e-07\\
+169	1.54446803267932e-08\\
+170	1.81444161747323e-07\\
+171	8.78831793802036e-08\\
+172	1.53613965535421e-06\\
+173	2.81415884074516e-06\\
+174	4.43530815074645e-07\\
+175	2.7224336876613e-08\\
+176	1.43297134609307e-06\\
+177	1.62527199539552e-09\\
+178	2.50972491576067e-10\\
+179	1.74075823612926e-11\\
+};
+\end{axis}
+\end{tikzpicture}%