multiphysicsL9.tex

%
% Copyright © 2014 Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
% for template copy, run:
%
% ~/bin/ct multiphysicsL1  multiphysicsLectureN tl1
%
%\input{../blogpost.tex}
%\renewcommand{\basename}{multiphysicsL9}
%\renewcommand{\dirname}{notes/ece1254/}
%\newcommand{\keywords}{ECE1254H}
%\input{../peeter_prologue_print2.tex}
%
%\usepackage{kbordermatrix}
%\usepackage{easybmat}
%
%\beginArtNoToc
%\generatetitle{ECE1254H Modeling of Multiphysics Systems.  Lecture 9: Conjugate gradient methods.  Taught by Prof.\ Piero Triverio}
%\chapter{Conjugate gradient methods}
\label{chap:multiphysicsL9}
%
%\section{Disclaimer}
%
%Peeter's lecture notes from class.  These may be incoherent and rough.
%
\section{Conjugate gradient convergence.}
\index{conjugate gradient!convergence}
%
For \( k \ll n \), convergence orders are given in \cref{tab:multiphysicsL9:20}.  A matrix norm, similar to \( \Norm{ \Bx } = \sqrt{ \Bx^\T \Bx } \), is defined
%
\captionedTable{Convergence.}{tab:multiphysicsL9:20}{
\begin{tabular}{|l|l|}
\hline
       & Full     \\ \hline
Direct & \(O(n^3)\)   \\ \hline
C.G.   & \(O(k n^2)\) \\ \hline
\end{tabular}
}
%
\begin{equation}\label{eqn:multiphysicsL9:40}
\Norm{\Bx}_M \equiv \sqrt{ \Bx^\T M \Bx }.
\end{equation}
%
Note that this norm is real valued for CG which only applies to positive definite matrices (or it will not converge), so this norm is real valued.
%
... lots on slides...
%
\begin{equation}\label{eqn:multiphysicsL9:60}
K(M) = \frac
{\sigma_{\mathrm{max}}}
{\sigma_{\mathrm{min}}}
\end{equation}
%
...
%
Some fast ways to estimate the conditioning number are required.
%
\sectionAndIndex{Gershgorin circle theorem.}
%
\maketheorem{Gershgorin circle theorem.}{thm:multiphysicsL9:1}{
Given \( M \), for any eigenvalue of \( M \) there is an \( i \in [1, n] \) such that
%
\begin{equation*}
\Abs{ \lambda - M_{i i} } \le \sum_{j \ne i} \Abs{ M_{i j} }
\end{equation*}
}
%
Consider this in the complex plane for row \( i \)
%
\begin{equation}\label{eqn:multiphysicsL9:80}
\begin{bmatrix}
M_{i 1} &
M_{i 2} &  \cdots M_{i i} & \cdots M_{i n}
\end{bmatrix}
\end{equation}
%
This inequality covers a circular region in the complex plane as illustrated in \cref{fig:lecture9:lecture9Fig1} for a two eigenvalue system.
%
\imageFigure{../figures/ece1254-multiphysics/lecture9Fig1}{Gershgorin circles.}{fig:lecture9:lecture9Fig1}{0.3}
%
These are called Gershgorin circles.
%
\makeexample{Leaky bar.}{example:multiphysicsL9:1}{
For the leaky bar of \cref{fig:lecture9:lecture9Fig2}, the matrix is
%
%\imageFigure{../figures/ece1254-multiphysics/lecture9Fig2}{Leaky bar}{fig:lecture9:lecture9Fig2}{0.3}
\imageFigure{../figures/ece1254-multiphysics/leaky-bar.pdf}{Leaky bar.}{fig:lecture9:lecture9Fig2}{0.3}
%
\begin{equation}\label{eqn:multiphysicsL9:100}
M =
\begin{bmatrix}
2  & -1 &        &    &    \\
-1 & 3  & -1     &    &    \\
   & -1 & 3      & -1 &    \\
   &    &        & \ddots &    \\
   &    &        &    & -1     \\
   &    &        & -1 & 2    \\
\end{bmatrix}
\end{equation}
%
The Gershgorin circles are \cref{fig:lecture9:lecture9Fig3}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture9Fig3}{Gershgorin circles for leaky bar.}{fig:lecture9:lecture9Fig3}{0.3}
%
This puts a bound on the eigenvalues
%
\begin{equation}\label{eqn:multiphysicsL9:120}
1 \le \lambda(M) \le 5,
\end{equation}
%
so that
%
\begin{equation}\label{eqn:multiphysicsL9:140}
K(M) =
\frac{\lambda_{\mathrm{max}}}
{\lambda_{\mathrm{max}}}
\le 5.
\end{equation}
%
On slides: example with smaller leakage to ground.
On slides: example with no leakage to ground.
%
These had, progressively larger and larger (possibly indefinite for the latter) conditioning number estimates.
%
The latter had the form of
%
\begin{equation}\label{eqn:multiphysicsL9:160}
M =
\begin{bmatrix}
     1  & -1  &  0  &  0 \\
    -1  &  2  & -1  &  0 \\
     0  & -1  &  2  &  1 \\
     0  &  0  & -1  &  1
\end{bmatrix}
\end{equation}
%
The exact eigenvalues for this system happens to be
\begin{equation}\label{eqn:multiphysicsL9:180}
\lambda \in
\{
 3.1069
 0.2833,
 1.3049 \pm 0.7545i
\}
\end{equation}
%
so the exact conditioning number is \( 3.1/0.28 \approx 11 \).
%
Compare this to the estimates, which are
%
\begin{equation}\label{eqn:multiphysicsL9:320}
\begin{aligned}
\Abs{ \lambda_1 - 1 } &\le 1 \\
\Abs{ \lambda_2 - 2 } &\le 2 \\
\Abs{ \lambda_3 - 2 } &\le 2 \\
\Abs{ \lambda_4 - 1 } &\le 1 \\
\end{aligned}
\end{equation}
%
These are two circles at \( z = 1 \) of radius 1, and two circles at \( z = 2 \) of radius 2, as plotted in \cref{fig:GershgorinL9:GershgorinL9Fig9}.
%
\imageFigure{../figures/ece1254-multiphysics/GershgorinL9Fig9}{Gershgorin circles for 4 eigenvalue system.}{fig:GershgorinL9:GershgorinL9Fig9}{0.3}
}
%
\section{Preconditioning.}
\index{preconditioning}
%
Goal is to take
%
\begin{equation}\label{eqn:multiphysicsL9:280}
M \Bx = \Bb
\end{equation}
%
and introduce an easy to invert matrix \( P \) to change the problem to
%
\begin{equation}\label{eqn:multiphysicsL9:200}
P^{-1} M \Bx = P^{-1} \Bb.
\end{equation}
%
This system has the same solution, but allows for choosing \( P \) to maximize the convergence speed.
%
\section{Symmetric preconditioning.}
\index{symmetric preconditioning}
%
Because the conjugate gradient methods requires a symmetric matrix, it is desirable to pick a preconditioning method that preserves the symmetric (and positive definite) nature of the matrix.  This is possible by splitting \( P \) into square root factors
%
\begin{equation}\label{eqn:multiphysicsL9:220}
P = P^{1/2} P^{1/2},
\end{equation}
%
and apply to \( M \Bx = \Bb \) as
%
\begin{equation}\label{eqn:multiphysicsL9:240}
P^{-1/2} M
\mathLabelBox{P^{-1/2} P^{1/2}}{\(= I\)} \Bx = P^{-1/2} \Bb.
\end{equation}
%
Now introduce a change of variables \( \By = P^{1/2} \Bx \), transforming the system to solve into
%
\begin{equation}\label{eqn:multiphysicsL9:260}
P^{-1/2} M P^{-1/2} \By = \Bb'.
\end{equation}
%
Some options
\begin{itemize}
\item Diagonal preconditioner: \( \BP = \text{diag}\{M\} \)
\item Incomplete LU or Cholesky factorization.  Cheap, approximate decomposition where a preconditioner \( M \simeq L U = P \) is picked.  An incomplete LU factorization would be easy to invert since lower or upper triangular matrices are easy to invert.  In Matlab the \( \text{ilu}() \) function can be used to do an incomplete LU factorization.
\item ... (many preconditioner are available).
\end{itemize}
%
For a symmetric positive definite matrix \( M \), an LU decomposition of the form \( M = L L^\T \), is called the \textAndIndex{Cholesky factorization}.
%
As an example consider the matrix
%
\begin{equation}\label{eqn:multiphysicsL9:340}
M = \left[
\begin{BMAT}(b){ccccc}{ccccc}
\begin{matrix} 2 & \\ & 2 \end{matrix} & & & \text{(small)} &\\
& \ddots & & & \\
& & \begin{matrix} 6 & \\ & 6 \end{matrix} & & \\
& & & \ddots & \\
& \text{(small)} & & & \begin{matrix} 50 & \\ & 50 \end{matrix} \\
\end{BMAT}
\right],
\end{equation}
%
for which a diagonal preconditioner can be used
%
\begin{equation}\label{eqn:multiphysicsL9:360}
P = \left[
\begin{BMAT}(b){ccccc}{ccccc}
\begin{matrix} 1/2 & \\ & 1/2 \end{matrix} & & & &\\
& \ddots & & & \\
& & \begin{matrix} 1/6 & \\ & 1/6 \end{matrix} & & \\
& & & \ddots & \\
& & & & \begin{matrix} 1/50 & \\ & 1/50 \end{matrix} \\
\end{BMAT}
\right].
\end{equation}
%
The preconditioned matrix will now have the form
%
\begin{equation}\label{eqn:multiphysicsL9:380}
M' = \left[
\begin{BMAT}(b){ccccc}{ccccc}
\begin{matrix} 1 & \\ & 1 \end{matrix} & & & \text{(small)} &\\
& \ddots & & & \\
& & \begin{matrix} 1 & \\ & 1 \end{matrix} & & \\
& & & \ddots & \\
& \text{(small)} & & & \begin{matrix} 1 & \\ & 1 \end{matrix} \\
\end{BMAT}
\right],
\end{equation}
%
so that the Gershgorin circles can all be found within a small radius of unity as sketched in \cref{fig:lecture9:lecture9Fig8}.
%
\imageFigure{../figures/ece1254-multiphysics/lecture9Fig8}{Gershgorin circles after preconditioning.}{fig:lecture9:lecture9Fig8}{0.3}
%
\section{Preconditioned conjugate gradient.}
\index{conjugate gradient!preconditioned}
%
It is possible to avoid inverting the preconditioner by requiring that the LU decomposition of \( P \) be easily computed.  Then
%
\begin{equation}\label{eqn:multiphysicsL9:300}
P \Bz^k = \Br^k
\end{equation}
%
can be solved by successive forward and back substitution.
%
More on slides...
%
%\EndArticle
%\EndNoBibArticle