conjugateGradientAlgorithm.tex

%
% Copyright © 2014  Ahmed Dorrah, Peeter Joot.  All Rights Reserved.
% Licenced as described in the file LICENSE under the root directory of this GIT repository.
%
\paragraph{CG without preconditioning}
\index{CG!see {conjugate gradient}}
The \textAndIndex{conjugate gradient} algorithm presented in the slides (without preconditioning) was
\makealgorithm{Conjugate gradient (CG).}{alg:conjugateGradientAlgorithm:401}{
\STATE \( \Bd^{(0)} = \Br^{(0)} \)
\REPEAT
\STATE \( \alpha_k = \frac{ \lr{ \Bd^{(k)} }^\T \Br^{(k)} }{ \lr{ \Bd^{(k)} }^\T M \Bd^{(k)} } \)
\STATE \( \Bx^{(k+1)} = \Bx^{(k)} + \alpha_k \Bd^{(k)} \)
\STATE \( \Br^{(k+1)} = \Br^{(k)} - \alpha_k M \Bd^{(k)} \)
\STATE \( \beta_k = \frac{ \lr{ M \Bd^{(k)} }^\T \Br^{(k+1)} }{ \lr{ M \Bd^{(k)} }^\T \Bd^{(k)} } \)
\STATE \( \Bd^{(k+1)} = \Br^{(k+1)} - \beta_k \Bd^{(k)} \)
\UNTIL{converged}
}
%
The repeated calculations are undesirable for actually coding this algorithm.  First introduce a temporary for the matrix product
%
\makealgorithm{CG with temporaries.}{alg:conjugateGradientAlgorithm:402}{
\STATE \( \Bd^{(0)} = \Br^{(0)} \)
\REPEAT
\STATE \( \Bq = M \Bd^{(k)} \)
\STATE \( \alpha_k = \frac{ \lr{ \Bd^{(k)} }^\T \Br^{(k)} }{ \lr{ \Bd^{(k)} }^\T \Bq } \)
\STATE \( \Bx^{(k+1)} = \Bx^{(k)} + \alpha_k \Bd^{(k)} \)
\STATE \( \Br^{(k+1)} = \Br^{(k)} - \alpha_k \Bq \)
\STATE \( \beta_k = \frac{ \Bq^\T \Br^{(k+1)} }{ \Bq^\T \Bd^{(k)} } \)
\STATE \( \Bd^{(k+1)} = \Br^{(k+1)} - \beta_k \Bd^{(k)} \)
\UNTIL{converged}
}
%
This has a lot more computation than the algorithm specified in \citep{shewchuk1994introduction} \S B.2.  It looks like the orthogonality properties can be used to recast the \( \Bd^{(k)} \cdot \Br^{(k)} \) products in terms of \( \Br^{(k)} \)
%
\begin{equation}\label{eqn:conjugateGradientAlgorithm:400}
\lr{ \Bd^{(k)} }^\T \Br^{(k)}
=
\Br^{(k)} \cdot \lr{ \Br^{(k)} + \beta_{k-1} \Bd^{(k-1) }
}.
\end{equation}
%
Since the new residual is orthogonal to all the previous search directions \( \Br^{(k)} \cdot \Bd^{(k-1)} = 0 \), the \( \beta_{k-1} \) term is killed leaving just \( \Br^{(k)} \cdot \Br^{(k)} \).
%
The numerator of \( \beta_k \) can be tackled by noting that the transformed direction vector \( \Bq \) is a scaled difference of residuals.  Taking dot products
%
\begin{equation}\label{eqn:conjugateGradientAlgorithm:420}
\begin{aligned}
\Bq \cdot \Br^{(k+1)}
&=
\inv{\alpha_k} \lr{ \Br^{(k)} - \Br^{(k+1)} } \cdot \Br^{(k+1)}
\\ &=
\inv{\alpha_k} \lr{ \cancel{\Bd^{(k)}} - \beta_{k-1} \cancel{\Bd^{(k-1)}} - \Br^{(k+1)} } \cdot \Br^{(k+1)}
\\ &=
-\inv{\alpha_k} \Br^{(k+1)} \cdot \Br^{(k+1)}.
\end{aligned}
\end{equation}
%
This gives
%
\begin{subequations}
\begin{equation}\label{eqn:conjugateGradientAlgorithm:440}
\alpha_k = \frac{ \lr{ \Br^{(k)} }^\T \Br^{(k)} }{ \lr{ \Bd^{(k)} }^\T \Bq }
\end{equation}
\begin{equation}\label{eqn:conjugateGradientAlgorithm:460}
\beta_k = -\frac{ \lr{\Br^{(k+1)}}^\T \Br^{(k+1)} }{ \alpha_k \Bq^\T \Bd^{(k)} },
\end{equation}
\end{subequations}
%
A final elimination of \( \alpha_k \) from \cref{eqn:conjugateGradientAlgorithm:460} gives
%
\begin{equation}\label{eqn:conjugateGradientAlgorithm:480}
\beta_k =
-\frac{ \lr{\Br^{(k+1)}}^\T \Br^{(k+1)} }{
\lr{ \Br^{(k)} }^\T \Br^{(k)}
}.
\end{equation}
%
All the pieces put together yield
%
\makealgorithm{Optimized conjugate gradient.}{alg:conjugateGradientAlgorithm:403}{
\STATE \( \Bd^{(0)} = \Br^{(0)} \)
\REPEAT
\STATE \( \Bq = M \Bd^{(k)} \)
\STATE \( \alpha_k = \frac{ \lr{ \Br^{(k)} }^\T \Br^{(k)} }{ \lr{ \Bd^{(k)} }^\T \Bq } \)
\STATE \( \Bx^{(k+1)} = \Bx^{(k)} + \alpha_k \Bd^{(k)} \)
\STATE \( \Br^{(k+1)} = \Br^{(k)} - \alpha_k \Bq \)
\STATE \( \beta_k =
-\frac{ \lr{\Br^{(k+1)}}^\T \Br^{(k+1)} }{
\lr{ \Br^{(k)} }^\T \Br^{(k)}
}
\)
\STATE \( \Bd^{(k+1)} = \Br^{(k+1)} - \beta_k \Bd^{(k)} \)
\UNTIL{converged}
}
%
This is now consistent with eqns 45-49 of \citep{shewchuk1994introduction}, with the exception that the sign of the \( \beta_k \) term has been flipped.
%
Since all the previous state does not have to be tracked the indexes can be dropped after introducing a couple temporary variables for the squared residuals
%
\makealgorithm{Optimized conjugate gradient, w/ temps.}{alg:conjugateGradientAlgorithm:404}{
\STATE \( \Bx = \Bx^{(0)} \)
\STATE \( \Bq = M \Bx \)
\STATE \( \Br = \Bb - \Bq\)
\STATE \( \Bd = \Br \)
\REPEAT
\STATE \( \Bq = M \Bd \)
\STATE \( \delta = \Br^\T \Br \)
\STATE \( \alpha = \frac{ \delta }{ \Bd^\T \Bq } \)
\STATE \( \Bx = \Bx + \alpha \Bd \)
\STATE \( \Br = \Br - \alpha \Bq \)
\STATE \( \delta' = \Br^\T \Br \)
\STATE \( \beta = \frac{ \delta' }{ \delta } \)
\STATE \( \delta = \delta' \)
\STATE \( \Bd = \Br + \beta \Bd \)
\UNTIL{converged}
}
%
This is coded in \nbcite{ps2a:conjugateGradientPainlessB2.m}{ps2a/conjugateGradientPainlessB2.m}.  That implementation allows for a preconditioner, but applies the preconditioner in a dump and inefficient way.
%
\paragraph{CG with preconditioning}
\index{conjugate gradient!preconditioning}
%
%%%This is now in good shape to consider application of a preconditioner.  Solution of the problem \( N \By = \Bc \), given the variables
%%%
%%%\begin{equation}\label{eqn:conjugateGradientAlgorithm:n}
%%%\begin{aligned}
%%%\Bc &= P^{-1/2} \Bb \\
%%%\By^{(0)} &= P^{1/2} \Bx^{(0)} \\
%%%\By &= P^{1/2} \Bx \\
%%%N &= P^{-1/2} M P^{-1/2}
%%%\end{aligned},
%%%\end{equation}
%%%
%%%has the form \( P^{-1/2} M P^{-1/2} P^{1/2} \Bx = P^{-1/2} \Bb \), which is equivalent to the original problem \( M \Bx = \Bb \).  Changing to the new variables in the algorithm gives
%%%
%%%\begin{algorithmic}
%%%\STATE \( \By = \By^{(0)} \)
%%%\STATE \( \Bq = N \By \)
%%%\STATE \( \Br = \Bc - \Bq\)
%%%\STATE \( \Bd = \Br \)
%%%\REPEAT
%%%\STATE \( \Bq = N \Bd \)
%%%\STATE \( \delta = \Br^\T \Br \)
%%%\STATE \( \alpha = \frac{ \delta }{ \Bd^\T \Bq } \)
%%%\STATE \( \By = \By + \alpha \Bd \)
%%%\STATE \( \Br = \Br - \alpha \Bq \)
%%%\STATE \( \delta' = \Br^\T \Br \)
%%%\STATE \( \beta = \frac{ \delta' }{ \delta } \)
%%%\STATE \( \delta = \delta' \)
%%%\STATE \( \Bd = \Br + \beta \Bd \)
%%%\UNTIL{converged}
%%%\RETURN \( P^{-1/2} \By \)
%%%\end{algorithmic}
%%%
%%%Now insert the original variables
%%%
%%%\begin{algorithmic}
%%%\STATE \( \By = P^{1/2} \Bx^{(0)} \)
%%%\STATE \( \Bq = P^{-1/2} M P^{-1/2} \By \)
%%%\STATE \( \Br = P^{-1/2} \Bb - \Bq\)
%%%\STATE \( \Bd = \Br \)
%%%\REPEAT
%%%\STATE \( \Bq = P^{-1/2} M P^{-1/2} \Bd \)
%%%\STATE \( \delta = \Br^\T \Br \)
%%%\STATE \( \alpha = \frac{ \delta }{ \Bd^\T \Bq } \)
%%%\STATE \( \By = \By + \alpha \Bd \)
%%%\STATE \( \Br = \Br - \alpha \Bq \)
%%%\STATE \( \delta' = \Br^\T \Br \)
%%%\STATE \( \beta = \frac{ \delta' }{ \delta } \)
%%%\STATE \( \delta = \delta' \)
%%%\STATE \( \Bd = \Br + \beta \Bd \)
%%%\UNTIL{converged}
%%%\RETURN \( P^{-1/2} \By \).
%%%\end{algorithmic}
%%%
%%%Let \( \Bs = P^{1/2} \Br, \Bz = P^{1/2} \By, \Bf = P^{1/2} \Bd, \Bt = P^{1/2} \Bq \)
%%%
%%%%\begin{algorithmic}
%%%%\STATE \( P^{1/2} \By = P \Bx^{(0)} \)
%%%%\STATE \( P^{1/2} \Bq = M P^{-1} P^{1/2} \By \)
%%%%\STATE \( P^{1/2} \Br = \Bb - P^{1/2} \Bq\)
%%%%\STATE \( P^{1/2} \Bd = P^{1/2} \Br \)
%%%%\REPEAT
%%%%\STATE \( P^{1/2} \Bq = M P^{-1} P^{1/2} \Bd \)
%%%%\STATE \( \delta = \Br^\T \Br \)
%%%%\STATE \( \alpha = \frac{ \delta }{ \Bd^\T \Bq } \)
%%%%\STATE \( P^{1/2}\By = P^{1/2}\By + \alpha P^{1/2}\Bd \)
%%%%\STATE \( P^{1/2}\Br = P^{1/2}\Br - \alpha P^{1/2}\Bq \)
%%%%\STATE \( \delta' = \Br^\T \Br \)
%%%%\STATE \( \beta = \frac{ \delta' }{ \delta } \)
%%%%\STATE \( \delta = \delta' \)
%%%%\STATE \( P^{1/2}\Bd = P^{1/2}\Br + \beta P^{1/2}\Bd \)
%%%%\UNTIL{converged}
%%%%\RETURN \( P^{-1} P^{1/2} \By \)
%%%%\end{algorithmic}
%%%
%%%\begin{algorithmic}
%%%%\STATE \( \Bz = P \Bx^{(0)} \)
%%%%\STATE \( \Bt = M P^{-1} \Bz \)
%%%\STATE \( \Bt = M \Bx^{(0)} \)
%%%\STATE \( \Bs = \Bb - \Bt\)
%%%\STATE \( \Bf = \Bs \)
%%%\REPEAT
%%%\STATE \( \Bt = M P^{-1} \Bf \)
%%%\STATE \( \delta = \Br^\T \Br \)
%%%\STATE \( \alpha = \frac{ \delta }{ \Bd^\T \Bq } \)
%%%\STATE \( \Bz = \Bz + \alpha \Bf \)
%%%\STATE \( \Bs = \Bs - \alpha \Bt \)
%%%\STATE \( \delta' = \Br^\T \Br \)
%%%\STATE \( \beta = \frac{ \delta' }{ \delta } \)
%%%\STATE \( \delta = \delta' \)
%%%\STATE \( \Bf = \Bs + \beta \Bf \)
%%%\UNTIL{converged}
%%%\RETURN \( P^{-1} \Bz \)
%%%\end{algorithmic}
%%%
%%%With symmetric \( P \) the dot products simplify \( \Br^\T \Br = \Bs^\T P^{-1} \Bs \), and \( \Bd^\T \Bq = \Bf^\T P^{-1} \Bt \), for
%%%
%%%\begin{algorithmic}
%%%\STATE \( \Bt = M \Bx^{(0)} \)
%%%\STATE \( \Bs = \Bb - \Bt\)
%%%\STATE \( \Bf = \Bs \)
%%%\REPEAT
%%%\STATE \( \Bt = M P^{-1} \Bf \)
%%%\STATE \( \delta = \Bs^\T P^{-1} \Bs \)
%%%\STATE \( \alpha = \frac{ \delta }{ \Bf^\T P^{-1} \Bt } \)
%%%\STATE \( \Bz = \Bz + \alpha \Bf \)
%%%\STATE \( \Bs = \Bs - \alpha \Bt \)
%%%\STATE \( \delta' = \Bs^\T P^{-1} \Bs \)
%%%\STATE \( \beta = \frac{ \delta' }{ \delta } \)
%%%\STATE \( \delta = \delta' \)
%%%\STATE \( \Bf = \Bs + \beta \Bf \)
%%%\UNTIL{converged}
%%%\RETURN \( P^{-1} \Bz \)
%%%\end{algorithmic}
%%%
%
An optimized preconditioned CG algorithm can be found in \citep{quarteroni2007numerical} \S 4.3.5.  Given \( \Bx^{(0)} \), that algorithm is
%
\makealgorithm{CG with preconditioning.}{alg:conjugateGradientAlgorithm:405}{
\STATE \( \Br^{(0)} = \Bb - A \Bx^{(0)} \)
\STATE \( \Bz^{(0)} = P^{-1} \Br^{(0)} \)
\STATE \( \Bp^{(0)} = \Bz^{(0)} \)
\REPEAT
\STATE \( \alpha_k = \frac{ \lr{ \Bz^{(k)}}^\T \Br^{(k)} }{ \lr{ \Bp^{(k)} }^\T A \Bp^{(k)} } \)
\STATE \( \Bx^{(k+1)} = \Bx^{(k)} + \alpha_k \Bp^{(k)} \)
\STATE \( \Br^{(k+1)} = \Br^{(k)} - \alpha_k A \Bp^{(k)} \)
\STATE \( P \Bz^{(k+1)} = \Br^{(k+1)} \)
\STATE \( \beta_k = \frac{ \lr{ \Bz^{(k+1)}}^\T \Br^{(k+1)} }{ \lr{ \Bz^{(k)} }^\T \Br^{(k)}} \)
\STATE \( \Bp^{(k+1)} = \Bz^{(k+1)} + \beta_k \Bp^{(k)} \)
\UNTIL{converged}
}
%
To adapt this to code, drop the indexes and introduce some temporaries
%
\makealgorithm{CG with preconditioning and temporaries.}{alg:conjugateGradientAlgorithm:406}{
\STATE \( \Br = \Bb - A \Bx \)
\STATE \( P \Bz = \Br \)
\STATE \( \Bp = \Bz \)
\STATE \( \delta = \Bz^\T \Br \)
\REPEAT
\STATE \( \Bq = A \Bp \)
\STATE \( \alpha = \delta /\lr{ \Bp^\T \Bq } \)
\STATE \( \Bx = \Bx + \alpha \Bp \)
\STATE \( \Br = \Br - \alpha \Bq \)
\STATE \( P \Bz = \Br \)
\STATE \( \delta' = \Bz^\T \Br \)
\STATE \( \beta = \delta' / \delta  \)
\STATE \( \delta = \delta' \)
\STATE \( \Bp = \Bz + \beta \Bp \)
\UNTIL{converged}
}
%
This is coded in \nbcite{ps2a:cgQuarteroniPrecond.m}{ps2a/cgQuarteroniPrecond.m}