lecture10.tex

%Template
%Copyright (C) 2019  Patrick Diehl
%
%This program is free software: you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation, either version 3 of the License, or
%(at your option) any later version.

%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.

%You should have received a copy of the GNU General Public License
%along with this program.  If not, see <http://www.gnu.org/licenses/>.
\providecommand\classoption{12pt}
\documentclass[\classoption]{beamer}

\input{./template/packages}

\input{template/variables.tex}

\pgfplotsset{compat=1.13}

% frame slide
\title{\coursename}
\subtitle{Lecture 9: Solvers, Conjugate gradient method, and BlazeIterative}

\author{\tiny Patrick Diehl \orcid{0000-0003-3922-8419}}
%\institute {
%    \href{}{\tt \scriptsize \today}
%}
\date {
 \tiny \url{\courseurl}
\vspace{2cm}
\doclicenseThis  
  
}

\begin{document} {
    \setbeamertemplate{footline}{}
    \frame {
        \titlepage
    }
}

\frame{

\tableofcontents

}

\AtBeginSection[]{
  \begin{frame}
  \vfill
  \centering
  \begin{beamercolorbox}[sep=8pt,center,shadow=true,rounded=true]{title}
    \usebeamerfont{title}\insertsectionhead\par%
  \end{beamercolorbox}
  \vfill
  \end{frame}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Reminder}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}{Lecture 9}
\begin{block}{What you should know from last lecture}
\begin{itemize}
\item Vectors and matrices
\item How to use Blaze for matrix and vector operations
\item How to compile a program using a external library
\end{itemize}
\end{block}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Solving linear equation systems}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}{Illustration of the linear system}

\begin{center}
\begin{tikzpicture}
\begin{axis}[grid=both,ymin=-6,ymax=4,xmax=8,xmin=-4,
               minor tick num=1,axis lines = middle,xlabel=$x_1$,ylabel=$x_2$]
\addplot [domain=-4:6, samples=101]{-1.5*x+1} node[pos=0.425] (endofplotsquare) {};
\node [right] at (endofplotsquare) {$3x_1+2x_2=2$};
\addplot [domain=-4:6, samples=101]{-(1/3)*x-(8/6)} node[pos=0.675] (endofplotsquare2) {};
\node [right] at (endofplotsquare2) {$2x_1+6x_2=-8$};
\end{axis}
\end{tikzpicture}
\end{center}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conjugate gradient method}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}{Conjugate gradient method}

\begin{block}{Properties:}
\begin{itemize}
\item Most popular iterative method for solving large systems of linear equations
\item Developed by Hestenes and Stiefel in 1952~\cite{hestenes1952methods}
\item Solves linear equation systems $\mathbf{A} \mathbf{x} = \mathbf{b}$
\item Each iteration does one matrix-vector multiplication and some computation of inner products
\end{itemize}
\end{block}

\begin{block}{Matrix}
\begin{itemize}
\item Symmetry $\mathbf{A}^T = \mathbf{A}$
\item Positive-definite $\mathbf{x}^T \mathbf{A} \mathbf{x} > 0, \forall \mathbf{x}>0$
\end{itemize}
\end{block}

More details about iterative methods~\cite{briggs2000multigrid}.
\end{frame}


\begin{frame}{The quadratic form}
Let us define the problem as a matrix:\\
\begin{center}
$\mathbf{A} \mathbf{x}=\mathbf{b}$
\end{center}
with 
\begin{center}
$\mathbf{A} = \begin{pmatrix}
 3 & 2 \\ 2 & 6 
\end{pmatrix}, \mathbf{x}=\begin{pmatrix}
x_1 \\ x_2  
\end{pmatrix}, \text{ and }  \mathbf{b}=\begin{pmatrix}
2 \\ -8  
\end{pmatrix}\text{.}$
\end{center}
\vspace{0.25cm}


Instead of solving $\mathbf{A}\mathbf{x}=\mathbf{b}$, the quadratic form, which is a function of $\mathbf{x}$ can be 

\begin{center}
$ f(\mathbf{x}) = \frac{1}{2} \mathbf{x}^T A \mathbf{x} - b^T \mathbf{x} + c$
\end{center}
can be minimized to find the solution $\mathbf{x}$.
\vspace{0.25cm}

\end{frame}

\begin{frame}{Plot of the quadratic form $f(\mathbf{x})$}

\begin{center}
\begin{tikzpicture}
\begin{axis}[ymin=-6,ymax=4,xmax=8,xmin=-4,zmin=0,zmax=160,xlabel=$x_1$,ylabel=$x_2$,zlabel={$f(\mathbf{x})$}]
\addplot3[
    surf,
    colormap/cool,
    %shader=interp,
    %shader=flat,
    samples=60] 
    {x*(1.5*x+y)-2*x+y*(x+3*y)+8*y};
\end{axis}
\end{tikzpicture}
\end{center}
\vspace{-0.5cm}
Finding the minimal point of $\mathbf{x}$ corresponds to the solution of $\mathbf{A} \mathbf{x}=\mathbf{b}$.
\end{frame}



\begin{frame}{Contour plot of the quadratic form $f(\mathbf{x})$}

\begin{center}
\begin{tikzpicture}
\begin{axis}[ymin=-6,ymax=4,xmax=6,xmin=-4,xlabel=$x_1$,ylabel=$x_2$]
\addplot gnuplot[raw gnuplot,thick,mark=none]
    {
        unset surface;
        set cntrparam levels auto 10;
        set contour;
        set yrange [-6:4];
        splot [x=-4:6] x*(1.5*x+y)-2*x+y*(x+3*y)+8*y;
    };
\addplot[mark=*] coordinates {(2,-2)};
\end{axis}
\end{tikzpicture}
\end{center}
\vspace{-0.5cm}
\end{frame}

\begin{frame}{Gradient of the quadratic form}

\begin{block}{Definition of the gradient:}
\vspace{0.25cm}
\begin{center}
$f'(\mathbf{x})=\begin{pmatrix}
\frac{\partial }{\partial x_1} f(\mathbf{x}) \\
\frac{\partial }{\partial x_2} f(\mathbf{x}) \\
\vdots \\
\frac{\partial }{\partial x_n} f(\mathbf{x})
\end{pmatrix}$
\end{center}
\end{block}
Applying a little bit of maths:

\begin{center}
$f'(\mathbf{x}) = \frac{1}{2} \mathbf{A}^T \mathbf{x} + \frac{1}{2} \mathbf{A} \mathbf{x}-\mathbf{b}$
\end{center}
and for a symmetric matrix $\mathbf{A}$, we get
\begin{center}
$f'(\mathbf{x})= \mathbf{A}\mathbf{x}-\mathbf{b}$
\end{center}
\end{frame}


\begin{frame}[fragile]{Gradient field}

\vspace{-0.25cm}
\begin{center}
\begin{tikzpicture}[scale=0.5]
\foreach \x in {-4,...,6} {
\foreach \y in {-6,...,4} {

\pgfmathsetmacro\mx{0.01*(3*\x+2*\y-2)};
\pgfmathsetmacro\my{0.01*(2*\x+6*\y+8)};

%\def\l{{sqrt(\u+\v)}};
%\node at (\x,\y) {\ml};
\draw[->,azure] (\x,\y) -- (\x+\mx,\y+\my);
}
}

\node at (2,-2)  {$\mathbf{x}$};

\draw[->,thin,gray] (0,-7) -- (0,5);
\node[above] at (0,5) {\small $x_1$};
\draw[->,thin,gray] (-5,0) -- (7,0);
\node[above] at (7,0) {\small $x_2$};
\end{tikzpicture}
\end{center}
\vspace{-0.25cm}
Since the gradient at the solution $\mathbf{x}$ is zero, we can set $f'(\mathbf{x})$ to zero to minimize $f(\mathbf{x})$.

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The method of the steepest decent}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}{The method of the steepest decent}

\begin{itemize}
\item We chose an random point $\mathbf{x}_0$
\item and slide down to the bottom of the quadratic form $f(\mathbf{x})$
\item by taking a series of steps $\mathbf{x}_1,\mathbf{x}_2,\ldots$
\item Each step we go to the direction which $f$ decreases most which is the opposite of $f'(\mathbf{x}_i)$ which is
\begin{center}
$-f'(\mathbf{x}_i)= \mathbf{b}-\mathbf{A}\mathbf{x}_i$
\end{center}
\end{itemize}

\end{frame}


\begin{frame}{The method of the steepest decent}

\begin{block}{Error}
\begin{center}
$\mathbf{e}_i = \mathbf{x}_i - \mathbf{x}$
\end{center}
Defines how far way we are from the exact solution at iteration $i$.
\end{block}

\begin{block}{Residual}
\begin{center}
$\mathbf{r}_i = \mathbf{b} - \mathbf{A}\mathbf{x}_i= -f'(\mathbf{x}_i)$
\end{center}
Defines how far away we are from the correct value for $\mathbf{b}$ in iteration $i$. 
\end{block}

\end{frame}


\begin{frame}{Visualization of the residual and error}
\begin{center}
\begin{tikzpicture}
\begin{axis}[ymin=-6,ymax=4,xmax=6,xmin=-4,xlabel=$x_1$,ylabel=$x_2$]
\addplot gnuplot[raw gnuplot,thick,mark=none]
    {
        unset surface;
        set cntrparam levels auto 10;
        set contour;
        set yrange [-6:4];
        splot [x=-4:6] x*(1.5*x+y)-2*x+y*(x+3*y)+8*y;
    };
\addplot[mark=*] coordinates {(2,-2)};
\node[below] at (2,-2) {\textcolor{black}{$\mathbf{x}$}};
\addplot[mark=*] coordinates {(-2,-2)};
\node[below] at (-2,-2) {$\mathbf{x}_0$};
\draw (2,-2) -- (-2,-2);
\node[below] at (-0,-2) {$\mathbf{e}_0$};
\draw (-2,-2) -- (8,12);
\node[left] at (-0,1) {$\mathbf{r}_0$};
\end{axis}
\end{tikzpicture}
\end{center}
\begin{center}
How far to go along the residual vector?
\end{center}
\end{frame}


\begin{frame}{Line search}

\begin{itemize}
\item We look at a starting point $\mathbf{x}_0=[-2,-2]^T$
\item from this point, we go along the direction of the steepest decent
\begin{center}
$\mathbf{x}_1 = \mathbf{x}_0 + \alpha \mathbf{r}_0$
\end{center}
\end{itemize}
How large to chose $\alpha$?
\end{frame}


\begin{frame}{Two surfaces}
\begin{center}
\begin{tikzpicture}
\begin{axis}[ymin=-6,ymax=4,xmax=8,xmin=-4,zmin=0,zmax=160,xlabel=$x_1$,ylabel=$x_2$,zlabel=$x_3$,view={53}{338}]
\addplot3[
    surf,
    colormap/cool,
    %shader=interp,
    %shader=flat,
    samples=60] 
    {x*(1.5*x+y)-2*x+y*(x+3*y)+8*y};
    \addplot3[
    domain=-4:10,
    samples = 60,
    samples y=0,
]
({-2+12*x},
{-2+8*x},
{x});
\addplot3[
    surf,
    colormap/cool,
    %shader=interp,
    %shader=flat,
    samples=60] 
    {x*(1.5*x+y)-2*x+y*(x+3*y)+8*y};
    \addplot3[
    domain=-4:10,
    samples = 60,
    samples y=0,
]
({-2+12*x},
{-2+8*x},
{25+x});
\end{axis}
\end{tikzpicture}
\end{center}
We need to find the point on the intersection of the two surfaces which minimizes $f$.
\end{frame}


\begin{frame}[fragile]{Parabola by the intersection of the two surfaces }

\begin{center}
\begin{tikzpicture}
\begin{axis}[xmin=-0.5,xmax=1,xlabel=$\alpha$,ylabel=$f(\mathbf{x}_0 + \alpha \mathbf{r}_0)$,axis lines = middle, grid= both]
\addplot[samples=100]{(-2+12*x)*(1.5*(-2+12*x)+(-2+8*x))-2*(-2+12*x)+(-2+8*x)*((-2+12*x)+3*(-2+8*x))+8*(-2+8*x)};
\end{axis}
\end{tikzpicture}
\end{center}
The minimum of this function is as $\frac{d}{d\alpha} f(\mathbf{x}_0 + \alpha \mathbf{r}_0) = 0$. 

\end{frame}

\begin{frame}{How to determine $\alpha$?}

Applying the chain rule: $\frac{d}{d\alpha} f(\mathbf{x}_0 + \alpha \mathbf{r}_0)= f'(\mathbf{x}_0 + \alpha \mathbf{r}_0)^T \mathbf{r}_0$. This expression is zero, if the two vectors are orthogonal. 

\begin{align*}
\mathbf{r}^T_1 \mathbf{r}_0 &= 0 \\
(\mathbb{b}-\mathbf{A} \mathbf{x}_1)^T \mathbf{r}_0 & = 0 \\
(\mathbb{b}-\mathbf{A} (\mathbf{x}_0 + \alpha \mathbf{r}_0))^T \mathbf{r}_0 & = 0 \\
(\mathbf{b}- \mathbf{A} \mathbf{x}_0)^T \mathbf{r}_0 - \alpha(\mathbf{A}\mathbf{r}_0)^T \mathbf{r}_0 &= 0 \\
(\mathbf{b}- \mathbf{A} \mathbf{x}_0)^T \mathbf{r}_0  &= \alpha(\mathbf{A}\mathbf{r}_0)^T \mathbf{r}_0 \\
\mathbf{r}^T_0 \mathbf{r}_0 &= \alpha \mathbf{r}^T_0 (\mathbf{A}\mathbf{r}_0) \\
\alpha &= \frac{\mathbf{r}^T_0 \mathbf{r}_0}{\mathbf{r}^T_0 \mathbf{A}\mathbf{r}_0}
\end{align*}

\end{frame}

\begin{frame}{Visualization of gradient of the previous step}
\begin{center}
\begin{tikzpicture}
\begin{axis}[ymin=-6,ymax=4,xmax=6,xmin=-4,xlabel=$x_1$,ylabel=$x_2$]
\addplot gnuplot[raw gnuplot,thick,mark=none]
    {
        unset surface;
        set cntrparam levels auto 10;
        set contour;
        set yrange [-6:4];
        splot [x=-4:6] x*(1.5*x+y)-2*x+y*(x+3*y)+8*y;
    };
\addplot[mark=*] coordinates {(2,-2)};
\node[below] at (2,-2) {$\mathbf{x}$};
\addplot[mark=*] coordinates {(-2,-2)};
\node[below] at (-2,-2) {$\mathbf{x}_0$};
\addplot[mark=*] coordinates {(0.08,-0.61333)};
\draw (-2,-2) -- (0.08,-0.61333);
\node[below] at (0.08,-0.61333) {$\mathbf{x}_1$};
\draw[->] (0.08,-0.61333) -- (0.1*-2.98666667, 0.1*4.48); 
%\node[left] at (-0,1) {$\mathbf{r}_0$};
\end{axis}
\end{tikzpicture}
\end{center}
\begin{center}
The gradient at $\mathbf{x}_1$ is orthogonal to $\mathbf{x}_{0}$.
\end{center}
\end{frame}

\begin{frame}{Algorithm}

\begin{enumerate}
\item $\mathbf{r}_0 = \mathbf{b} - \mathbf{A} \mathbf{x}_0$
\item If $\vert \mathbf{r}_0 \vert < \epsilon$ return $\mathbf{x}_0$
\item $\mathbf{r}_i = \mathbf{b} - \mathbf{A} \mathbf{x}_i$
\item $\alpha_i = \frac{\mathbf{r}^T_0 \mathbf{r}_0}{\mathbf{r}^T_0 \mathbf{A}\mathbf{r}_0}$
\item $\mathbf{x}_{i+1} = \mathbf{x}_i \alpha_i \mathbf{r}_i$
\item If $\vert \mathbf{r}_i \vert < \epsilon$ return  $\mathbf{x}_i$ 
\item Go to (3)
\end{enumerate}


\end{frame}


\begin{frame}{Visualization of the line search}
\begin{center}
\begin{tikzpicture}
\begin{axis}[ymin=-6,ymax=4,xmax=6,xmin=-4,xlabel=$x_1$,ylabel=$x_2$]
\addplot gnuplot[raw gnuplot,thick,mark=none]
    {
        unset surface;
        set cntrparam levels auto 10;
        set contour;
        set yrange [-6:4];
        splot [x=-4:6] x*(1.5*x+y)-2*x+y*(x+3*y)+8*y;
    };
\addplot[mark=*] coordinates {(2,-2)};
\node[above] at (2,-2) {$\mathbf{x}$};
\addplot[mark=*] coordinates {(-2,-2)};
\node[below] at (-2,-2) {$\mathbf{x}_0$};
\addplot[mark=*] coordinates {(0.08,-0.61333)};
\draw (-2,-2) -- (0.08,-0.61333);
\node[below] at (0.08,-0.61333) {$\mathbf{x}_1$};
\addplot[mark=*] coordinates {(1.00444444, -2.)};
\draw (0.08,-0.61333) -- (1.00444444, -2.);
\node[below] at (1.00444444, -2.) {$\mathbf{x}_2$};
\addplot[mark=*] coordinates {(1.52213333, -1.65487407)};
\draw (1.00444444, -2.) -- (1.52213333, -1.65487407);
\node[above] at (1.52213333, -1.65487407) {$\mathbf{x}_3$};
\addplot[mark=*] coordinates {(1.75221728, -2.)};
\draw (1.52213333, -1.65487407) -- (1.75221728, -2.);
\node[below] at (1.75221728, -2.) {$\mathbf{x}_4$};
\draw (1.75221728, -2.) -- (1.8810643 , -1.91410199);
\draw (1.8810643 , -1.91410199) -- (1.93832964, -2.);
\end{axis}
\end{tikzpicture}
\end{center}
\begin{center}
The solution after five steps $\mathbf{x}_5=[1.93832964, -2.]^T$.
\end{center}
\end{frame}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Blaze Iterative}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{frame}[fragile]{About BlazeIterative\footnote{\tiny\url{https://github.com/STEllAR-GROUP/BlazeIterative}}}
This is a set of iterative linear system solvers intended for use with the Blaze library, a high-performance C++ linear algebra library. The API is currently based on a tag-dispatch system to choose a particular algorithm.

\begin{block}{Usage}
\begin{lstlisting}[language=bash]
#Install
tar -xvf blaze_iterative.gz
cd blaze_iterative
cp -r ./blaze_iterative /home/patrick/

#Compile
g++ -I/home/diehlpk/blaze 
    -I/home/patrick/blaze_iterative BlazeTest.cpp
\end{lstlisting}
\end{block}


\end{frame}

\begin{frame}[fragile]{Conjugate gradient example}
\begin{lstlisting}
#include "BlazeIterative.hpp"

using namespace blaze;
using namespace blaze::iterative;

std::size_t N = 10;
DynamicMatrix<double,false> A(N,N, 0.0);
DynamicVector<double> b(N, 0.0);
DynamicVector<double> x1(N, 0.);

//Initialize the matrix

// Solve the system
ConjugateGradientTag tag;
auto x2 = solve(A,b,tag);

\end{lstlisting}

\end{frame}

\begin{frame}{Available algorithms}

\begin{block}{Solvers}
\begin{itemize}
\item Conjugate Gradient
\item Preconditioned CG
\item BiCGSTAB
\item Generalized minimal residual method (GMRES),
\end{itemize}
\end{block}


\begin{block}{Eigenvalues}
\begin{itemize}
\item Lanczos
\end{itemize}
\end{block}
More details about solvers~\cite{barrett1994templates}.
\end{frame}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Summary}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}{Summary}
\begin{block}{After this lecture, you should know}
\begin{itemize}
\item Linear equation systems
\item Conjugate gradient method
\item BlazeIterative
\end{itemize}
\end{block}
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgment}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}{Acknowledgment}
\begin{itemize}
\item The very nice example for the introduction of the conjugate gradient method was adapted from:\\
\vspace{0.25cm}
Shewchuk, Jonathan Richard. "An introduction to the conjugate gradient method without the agonizing pain." (1994).

\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{References}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}[t, allowframebreaks]
\frametitle{References}
\bibliographystyle{plain}
\bibliography{bib}
\end{frame}


\end{document}