\documentclass[11pt,a4paper]{article}
\usepackage{graphicx}
\usepackage{amsmath,amssymb}
\usepackage{xcolor}
\usepackage{blkarray}
\usepackage{graphicx}
\usepackage{xeCJK}
\usepackage{hyperref}
\usepackage{subcaption}
\newcommand{\overbar}[1]{\mkern 1.5mu\overline{\mkern-1.5mu#1\mkern-1.5mu}\mkern 1.5mu}
\usepackage[OT2,T1]{fontenc}
\usepackage[utf8]{inputenc}
\DeclareSymbolFont{cyrletters}{OT2}{wncyr}{m}{n}
\DeclareMathSymbol{\Sha}{\mathalpha}{cyrletters}{"58}
\usepackage{scrextend}
\usepackage{lipsum}
\usepackage{listings}
\usepackage{multicol}
\usepackage{mathtools}
 \usepackage{nccmath}
\usepackage[ruled,vlined]{algorithm2e}
\setlength{\columnseprule}{1pt}
\def\columnseprulecolor{\color{blue}}
\setlength{\columnsep}{2mm}
\setlength\columnwidth{4cm}

\usepackage[scale=0.77]{geometry}


\begin{document}

\title{Conjugate Gradient Method: Steepest Descent with Orthogonal Directions under Change of Coordinate}
\author{J.R}
\date{\today}

\maketitle

In this note, I want to guide you through the motivation and derivation of the bizarre Conjugate Gradient Method (CG Method) 
and show you that a method as intricate as CG Method is nothing more than a steepest descent with a change of coordinate.

\tableofcontents

\newpage
\section{Linear Equation and Quadratic Form}
Let's begin with a linear equation 
\begin{align*}
	Ax = b \tag{1.1}
\end{align*}
where $A$ is symmetric and positive definite. In most situation, we can solve this linear equation by factoring $A$ in to a product of a 
lower-triangle matrix and an upper-triangle matrix and doing backward/forward substitution. However, when $A$ is huge, for example, 
1 million by 1 million, the factoring seems to be intractable. Thus we have to come up with other methods for solving a large-scale linear 
system.

The quadratic form of (1.1) is 
\begin{align*}
	f(x) = \frac{1}{2} x^TAx - b^Tx + c	\tag{1.2}
\end{align*}
We call (1.2) the quadratic form of (1.1) because the solution to (1.2) $x^* = A^{-1}b$ is the solution to (1.1). In other words, if we can solve
optimization problem (1.2), we can solve the linear equation (1.1). Let's introduce some definition which will be used across this note.
\begin{align*}
	x^* &= A^{-1}b = \text{argmin}_{x} f(x)	\\
	e_{(i)} &= x_{(i)} - x^*\\
	r_{(i)} &= -Ae_{(i)} = b - Ax_{(i)} = - \nabla f(x_{(i)})
\end{align*}
$e^{(i)}$ is called the error at $x_{(i)}$ and $r_{(i)}$ is called the  residual at $x_{(i)}$.

And in this note, we will have a running example where 
$A = 2 \cdot \begin{psmallmatrix}0.416 & 0.168\\ 0.168 & 0.416\end{psmallmatrix}$,
$b= \begin{psmallmatrix}5\\ 2.5\end{psmallmatrix}$ 
and $c=6.25$, so that 
$f(x) = (x-u)^TP(x-u)$, where $P$ has eigenvalue decomposition $P=$
$\begin{psmallmatrix}2/\sqrt{5} & -1\sqrt{5}\\ 1\sqrt{5} & 2\sqrt{5}\end{psmallmatrix}$
$\begin{psmallmatrix}1/4 & 0\\ 0 & 1/25\end{psmallmatrix}$
$\begin{psmallmatrix}2\sqrt{5} & -1\sqrt{5}\\ 1\sqrt{5} & 2\sqrt{5}\end{psmallmatrix}^T$
and $u=[2,1]$.
The contour lines are showed below, which are 
ellipsoids tilted and centered at $[2,1]$.

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.35]{eg_contour.png}
\end{center}
\caption{Contour lines of running example}
\end{figure}

\newpage
\section{Steepest Descent}
(1.2) is a convex optimization problem, we can solve it with many methods like newton method. However, in second-order methods, we need
to invert $A$, which is too expansive so that we want to avoid. 
So we may want to use some first-order method. The first method coming to mind is gradient
descent, but choosing step size in gradient descent can be tricky. For a problem has a simple quadratic form,
we can easily do line search, which leads us to steepest descent. At each step in steepest descent, we 
evaluate the gradient, and go in the direction of negative gradient, which is $r_{(i)}$, until we reach a 
minimum along the direction of $r_{(i)}$.  The step size is calculated by setting the gradient to 0
\begin{align*}
	\frac{d}{d \alpha} f(x_{(i)} + \alpha_{(i)} r_{(i)}) = 0	\tag{2.1} \\
	r_{(i)}^T [ b - A(x_{(i)} + \alpha_{(i)} r_{(i)}) ] = 0		\tag{2.2} \\	
	\alpha_{(i)} = \frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}	\tag{2.3}
\end{align*}

Steepest descent on the running example with different starting point is shown in the following picture, 
stopping criterion is when $||r_{(i)}||_2 \leq 1e-6$.

\begin{figure}[h!]

			\hspace{-15mm}
			\begin{subfigure}[b]{0.55\textwidth}
			\includegraphics[scale=0.35]{steepest_descent_1.png}
			\caption{}
			\end{subfigure}
			\hspace{-8mm}
			\begin{subfigure}[b]{0.55\textwidth}
			\includegraphics[scale=0.35]{steepest_descent_2.png}
			\caption{}
			\end{subfigure}
		\caption{Steepest descent on running example}
\end{figure}

Notice that search directions are orthogonal to the previous one. We can easily get this result from (2.1).
\begin{align*}	
	  &	\frac{d}{d \alpha} f(x_{(i)} + \alpha_{(i)} r_{(i)})	\\
	=& 	r_{(i)}^T \nabla f(x_{(i)} + \alpha_{(i)} r_{(i)}) \\
	=&	-r_{(i)}^Tr_{(i+1)}	\tag{2.4}
\end{align*}
when we set $\frac{d}{d \alpha} f(x_{(i)} + \alpha_{(i)} r_{(i)})$ to $0$, $r_{(i)}^Tr_{(i+1)}=0$ naturally follows.

One problem with this method is that we step in the same direction for many times. For example, 
in figure 2(a), we 
first take a step to the right(roughly speaking), then take a step down, and right, and down, .... Wouldn't it be
nice if we can take a full step to the right and never go in that direction again. To be more concrete, after 
take a step in the direction of $d_{(i)}$ with step size $\alpha_{(i)}$, the error at $x^{(i)} + \alpha_{(i)}d_{(i)}$
is orthogonal to previous searching direction $d_{(i)}$ so that we don't need to go in direction of
$d_{(i)}$ ever again, as illustrated in the following picture. 

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.25]{orthogonal_directions.png}
\end{center}
\caption{ }
\end{figure}

In figure 3, $\hat{x}$ is the minimizer subject to direction $d_{(i)}$. We can see that, to make $e_{(i+1)}$ 
orthogonal to the previous searching direction, we have to overshoot or undershoot compare to $\hat{x}$
depending on the shape of the ellipsoid. 
We can calculate the step size $\alpha_{(i)}$ by
\begin{align*}
	d_{(i)}^T e^{(i+1)} &= 0	\tag{2.5} \\
	d_{(i)}^T (e^{(i)} + \alpha_{(i)} d_{(i)}) &= 0	\\
	\alpha_{(i)} &= - \frac{d_{(i)}^Te^{(i)}}{d_{(i)}^Td_{(i)}} \tag{2.6}
\end{align*}

That means we cannot calculate $\alpha_{(i)}$ without knowing $e_{(i)}$, and the "orthogonal search 
directions" idea is infeasible. The problem with "orthogonal search directions" idea is that we have to 
"overshoot" or "undershoot". With that in mind, can we come up with something just like the 
"orthogonal search directions" but in which we don't have to overshoot? Fortunately, yes. And that is called
"Conjugate Directions".

\newpage
\section{Conjugate Directions}
In the previous section, we found that if we want $e_{(i+1)}$ to be orthogonal to previous search direction
$d_{(i)}$, we have to "overshoot" or "undershoot". 

The reason why we have to overshoot or undershoot
is the asymmetry of ellipsoid. When we have a unit circle or unit ball, which is homogeneous in all direction, 
the resulting error after minimizing along a 
direction is guaranteed to be orthogonal 
to the search direction, as illustrated in the following picture. 

\begin{figure}[h!]
			\hspace{-15mm}
			\begin{subfigure}[b]{0.5\textwidth}
			\includegraphics[scale=0.25]{circle_1.png}
			\caption{}
			\end{subfigure}
			\begin{subfigure}[b]{0.5\textwidth}
			\includegraphics[scale=0.25]{circle_2.png}
			\caption{}
			\end{subfigure}
		\caption{Orthogonal directions on unit circle}
\end{figure}
We can see from the picture that when minimize a quadratic function with $A=I$ along a direction, resulting 
error is orthogonal to error. Mathematical proof can be easily given.

Now the problem is, we have a quadratic function $f(x) = \frac{1}{2}x^TAx - b^Tx + c$, whose contour lines
are ellipsoids, not unit balls. How can we make it a 
unit ball? The answer is, not surprisingly, yes and by changing coordinate. 

Let $u = A^{1/2}x$, we have
\begin{align*}
	f(x) = \tilde{f}(u) 
		&= \frac{1}{2} u^Tu - (A^{-1/2}b)^Tu + c	\\
		&= \frac{1}{2}\bigg[ (u - A^{-1/2}b)^T (u - A^{-1/2}b) \bigg] - \frac{b^TA^{-1}b}{2} + c \tag{3.1}
\end{align*}
by dropping constant terms, we can see the contour lines of $\tilde{f}$ are unit circles centered at $A^{-1/2}b$, 
as illustrated in the following picture. 
\newpage
\begin{figure}[h!]
\hspace{-10mm}
\includegraphics[scale=0.47]{transform.png}
\caption{Change of coordinate $u = A^{1/2}x$}
\end{figure}

We can easily apply orthogonal directions as search direction in the new coordinate system, and when we 
find $u^* = A^{-1/2}b$, we can get $x^* = A^{-1/2}u^*$. (Actually we do not have to do this, but this is the idea.)

Thus, instead of orthogonal directions in the old system of $x$ where the ellipsoids lie, we now need a group
of direction which are orthogonal to each other in the new coordinate system.

Suppose $d_{(i)}$ and $d_{(j)}$ are two vectors in old coordinate system, under change of coordinate, we have
$\tilde{d}_{(i)} = A^{1/2}d_{(i)}$ and $\tilde{d}_{(j)} = A^{1/2}d_{(j)}$. Since $d_{(i)}$ and $d_{(j)}$
are orthogonal under after of coordinate, we have
\begin{align*}
	(A^{1/2}d_{(i)})^TA^{1/2}d_{(j)} &= 0	\\
	d_{(i)}^TAd_{(j)} &= 0	\tag{3.2}
\end{align*}

If $d_{(i)}$ and $d_{(i)}$ satisfies (3.2), we say $d_{(i)}$ and $d_{(i)}$ are $A$-orthogonal, or in more 
fancy words, $A$-conjugate. Roughly speaking, in this context, conjugate
means "orthogonal under linear transformation". The following pictures shows some conjugate vectors
when $A = \begin{psmallmatrix} 1/4 & 0 \\ 0 & 1\end{psmallmatrix}$.

\begin{figure}[h!]
\includegraphics[scale=0.22]{conjugates.png}
\caption{Conjugate vectors}
\end{figure}


\section{Steepest Descent and Conjugate Directions}
When applied to quadratic programming, steepest descent, combined with conjugate directions, has 
very interesting property and this can be seen with a change of coordinate. This section will be divided into 
two columns, the left column is what we do in the original coordinate system, the right column is what happens
in the new coordinate system $u = A^{1/2}x$.  

% two columns page 23
\begin{multicols}{2}
\begin{align*}
	f(x) = \frac{1}{2}x^TAx - b^Tx + c
\end{align*}
$x^*$ solves $Ax = b$.
\begin{align*}
	e_{(i)} &= x_{(i)} - x^*	\\
	\nabla f(x) &= Ax - b	\\
	r_{(i)} &= A(x^* - x_{(i)}) = -Ae_{(i)}	\\
		& = -Ax_{(i)} + b = -\nabla f(x_{(i)})
\end{align*}


\columnbreak
\begin{align*}
	\tilde{f}(x) = \frac{1}{2}u^Tu - (A^{-1/2}b)^Tu + c
\end{align*}
$u^*$ solves $Iu^* = A^{-1/2}b$
\begin{align*}
	\tilde{e}_{(i)} &= u_{(i)} - u^*	\\
	\nabla \tilde{f}(u) &= u - A^{-1/2}b	\\
	\tilde{r} &= I(u^* - u_{(i)}) = -I \tilde{e}_{(i)} \\
		&= A^{-1/2}b - u_{(i)} = - \nabla \tilde{f}(u_{(i)})
\end{align*}
\end{multicols}


When we search in direction of $d_{(i)}$ in the original coordinate system, we are searching in direction 
$A^{1/2}d_{(i)}$ in the new coordinate system. And we step in the searching direction with a step size 
$\alpha_{(i)}$ such that the resulting point minimize the function along this direction. The step size can be 
calculated by setting the gradient to 0.

\begin{multicols}{2}
\begin{align*}
	x_{(i+1)} = x_{(i)} + \alpha_{(i)} d_{(i)}	\\
	\frac{d}{d\alpha_{(i)}} f(x_{(i)} + \alpha_{(i)}d_{(i)}) &= 0 \\
	\nabla f(x_{(i)} + \alpha_{(i)}d_{(i)})^T d_{(i)} &= 0\\
	e_{(i+1)}^T A d^{(i)} &= 0
\end{align*}



\columnbreak
\begin{align*}
	u_{(i+1)} = u_{(i)} + \alpha_{(i)} A^{1/2}d_{(i)}	\\
	\frac{d}{d\alpha_{(i)}} \tilde{f}(u_{(i)} + \alpha_{(i)} A^{1/2}d_{(i)}) &= 0\\
	\nabla \tilde{f}(u_{(i)} + \alpha_{(i)}A^{1/2}d_{(i)})^T (A^{1/2}d_{(i)}) &= 0\\
	\tilde{e}_{(i+1)}^TA^{1/2}d_{(i)} &= 0
\end{align*}
\end{multicols}

Let's translate the above math into English. When we go from $x_{(i)}$ in the direction $d_{(i)}$, 
we are going from $u_{(i)}$ in the direction $A^{1/2}d_{(i)}$ in the new coordinate system. 
If $x^{(i)} + \alpha_{(i)}d_{(i)}$ minimizes $f$ in the direction of $d_{(i)}$, 
then $u_{(i)} + \alpha_{(i)}A^{1/2}d_{(i)}$ minimizes $\tilde{f}$ in the direction of $A^{1/2}d_{(i)}$ 
(since it's a linear transformation). However, in the new coordinate system,
because the contour lines of $\tilde{f}$ are circles/balls, the resulting
error $\tilde{e}_{(i+1)} = A^{1/2}e_{(i+1)}$ is guaranteed to be orthogonal to the previous 
searching direction $A^{1/2}d_{(i)}$. Thus we can know that $e_{(i+1)}^T A d^{(i)} = 0$. 

In other words, in the original coordinate system, when we minimize the function along a direction $d_{(i)}$, 
the resulting error $e_{(i+1)}$ and searching direction $d_{(i)}$ are $A$-orthogonal, or say, conjugate, as 
shown in the following picture.

\newpage
\begin{figure}[h!]
\includegraphics[scale=0.3]{conjugate_steepest_descent_transform.png}
\caption{Steepest descent under linear transformation}
\end{figure}

\section{Almost Conjugate Gradient Method}
In the previous section, we see that when we apply steepest descent to a quadratic function, in the new
coordinate system, the resulting error is orthogonal to the searching direction. That means, in the new
coordinate system, after minimizing along a direction, we don't need to go into that direction again, or
say, the later searching directions are gonna be orthogonal to the current searching direction. 

Thus, if we have $n$ orthogonal vectors that spans the space of new coordinate system, where $n$ is the 
dimension of the space, then we can
iteratively minimize over each of these directions. 
In other word, we need a group of vectors $d_{(i)}$ 
in the original coordinate system that are $A$-orthogonal. 

Instead of having $n$ $A$-orthogonal vectors at the beginning, we can gradually build up a group of them by
Gram-Schmidt conjugation.  Let me brutally show you the "almost conjugate gradient method" and then 
explain it.

\begin{align*}
	d_{(0)} &= r_{(0)}	\tag{5.1}\\
	d_{(i)} &= r_{(i)} - \sum_{j=0}^{i-1} \frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}d_{(j)}	\tag{5.2}\\
	\alpha_{(i)} &= \frac{r_{(i)}^Td_{(i)}}{d_{(i)}^TAd_{(i)}} \tag{5.3}\\
	x_{(i+1)} &= x^{(i)} + \alpha_{(i)} d_{(i)} \tag{5.4}\\
\end{align*}

At first step, we perform steepest descent along the direction of negative gradient $r_{(0)}$.
At each following step at $x_{(i)}$, we evaluate the negative gradient $r_{(i)}$ and build $d_{(i)}$ which is 
$A$-orthogonal to all previous searching directions. First, we take $r_{(i)}$ and previous searching directions 
$d_{(j)} (j<i)$ to the new coordinate system, we get 
\begin{align*}
	\tilde{r}_{(i)} &= A^{1/2}r_{(i)}	\tag{5.5}\\
	\tilde{d}_{(j)} &= A^{1/2}d_{(j)}	~~~~~(j<i) \tag{5.6}
\end{align*}
Then we project $\tilde{r}_{(i)}$ on each of $\tilde{d}_{(j)}$, the projection of $\tilde{r}_{(i)}$ on $\tilde{d}_{(j)}$ is
\begin{align*}
	P_{\tilde{d}_{(j)}}(\tilde{r}_{(i)})	
	&=	\tilde{d}_{(j)} \bigg[ \tilde{d}_{(j)}^T\tilde{d}_{(j)} \bigg]^{-1}\tilde{d}_{(j)}^T \tilde{r}_{(i)} \\
	&= A^{1/2}d_{(j)} \bigg[(A^{1/2}d_{(j)})^T (A^{1/2}d_{(j)}) \bigg]^{-1} (A^{1/2}d_{(j)})^T A^{1/2}r_{(i)}	\\
	&=\frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}A^{1/2}d_{(j)}	\tag{5.7}
\end{align*}
then we subtract $\tilde{r}_{(i)}$ by it's projection on all $\tilde{d}_{(j)}$ such that $j<i$, 
then the remaining part is surly 
orthogonal to all previous $\tilde{d}_{(j)}$ (in the new coordinate system)
\begin{align*}
	\tilde{d}_{(i)} &= \tilde{r}_{(i)} - \sum_{j=0}^{i-1}\frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}A^{1/2}d_{(j)}	\\
			&= A^{1/2}r_{(i)}- \sum_{j=0}^{i-1}\frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}A^{1/2}d_{(j)} \tag{5.8}
\end{align*}
We have $\tilde{d}_{(i)}$ orthogonal to all previous $\tilde{d}_{(j)}$ in the new coordinate system. The next 
thing to do is to bring $\tilde{d}_{(i)}$ to the original coordinate system
\begin{align*}
	d_{(i)} &= A^{-1/2}\tilde{d}_{(i)} \\
		&= r_{(i)}- \sum_{j=0}^{i-1}\frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}d_{(j)}
\end{align*}
thus we get (5.2). 

Since $\tilde{d}_{(i)}$ is orthogonal to all $\tilde{d}_{(j)}$ in new coordinate system, $d_{(i)}$ is $A$-orthogonal
to. all ${d}_{(j)}$ in the original system.

Now we have the searching direction $d_{(i)}$, what we need is the step size for steepest descent 
$\alpha_{(i)}$, which we can calculate by setting the gradient to $0$.
\begin{align*}
	\frac{d}{d\alpha_{(i)}} f(x_{(i)} + \alpha_{(i)} d_{(i)}) &= 0	\\
	\nabla f(x_{(i)} + \alpha_{(i)} d_{(i)})^Td_{(i)} &= 0	\\
	\bigg[A(x_{(i)} + \alpha_{(i)} d_{(i)}) - b\bigg]^Td_{(i)} &= 0	\\
	-r_{(i)}^Td_{(i)} + \alpha_{(i)}d_{(i)}^TAd_{(i)} &= 0 \\
	\alpha_{(i)} &= \frac{r_{(i)}^Td_{(i)}}{d_{(i)}^TAd_{(i)}}
\end{align*}
thus we get (5.3).

Figure 8 and Figure 9 show progress of the "almost conjugate gradient" with different starting points.

\begin{figure}[h!]
\hspace{-8mm}
\includegraphics[scale=0.25]{almost_cg1.png}
\caption{Almost conjugate gradient method}
\end{figure}

\begin{figure}[h!]
\hspace{-8mm}
\includegraphics[scale=0.25]{almost_cg2.png}
\caption{Almost conjugate gradient method}
\end{figure}


As you can see, applying steepest descent with conjugate directions in the original coordinate system is 
equivalent to applying steepest descent with orthogonal directions in the new coordinate system, and the reason
we can do this in the new coordinate system is that the contour lines in the new coordinate system are circles and
resulting error of steepest descent is guaranteed to be orthogonal to previous searching direction.

\newpage
\section{Simplifing "Almost CG Method"}
If you have understand the "almost conjugate gradient method", you now get the core idea of the conjugate 
gradient method. Now we are one step away from the real conjugate gradient method, and that step is to 
reduce the complexity of  the "almost conjugate gradient method". One thing you may notice that when we build
$d_{(i)}$ according to (5.2), we must use all previous $d_{(j)}$ so that we must keep all $d_{(j)}$ in memory and
that can be a large storage. However, by exploring the properties of spaces spanned by all $d_{(j)}$ and by 
$r_{(i)}$, we can simplify the process of building $d_{(i)}$. 


Suppose $D_i = \text{span}\{d_{(0)}, d_{(1)}, ..., d_{(i)}\}$ is the space spanned by first $i$ searching directions 
and $\text{span}\{r_{(0)}, r_{(1)}, ..., r_{(i)}\}$ is the space spanned by first $i$ residuals. Now we can build 
the fact that $\text{span}\{r_{(0)}, r_{(1)}, ..., r_{(i)}\}$ = $D_i$ using mathematical induction.

\begin{align*}
	& \text{if}~~\text{span}\{r_{0}, ..., r_{(i)}\} = \text{span}\{d_{(0)}, d_{(1)}, ..., d_{(i)}\} = D_i	\\
	& \because d_{(i+1)} = r_{(i+1)} - \sum_{j=0}^{i} \frac{d_{(j)}^TAr_{(i+1)}}{d_{(j)}^TAd_{(j)}}d_{(j)} \\
	& \therefore \text{span}\{d_{(0)}, ..., d_{(i+1)}\} = \text{span}\{d_{(0)},...,d_{(i)}, r_{(i+1)}\}	\\
	&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ = \text{span}\{D_i, r_{(i+1)}\} \\
	&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ = \text{span}\{r_{(0)},...,r_{(i)}, r_{(i+1)}\} \\
	& \because d_{(0)} = r_{(0)} \\
	& \therefore \text{span}\{d_{(0)}\} = \text{span}\{r_{(0)}\} \\
	& \therefore \text{span}\{d_{(0)}, d_{(1)}, ..., d_{(i)}\} = \text{span}\{r_{0}, ..., r_{(i)}\} ~\forall ~ i
\end{align*}

Based on this conclusion, we can get more equations which can help us to simplify the 
"almost conjugate gradient method". 

First, since all search directions in the transformed coordinate system
$\tilde{d}_{(k)}$ are orthogonal to each other, the error in transformed coordinate system
can be written as $\tilde{e}_{(j)} = \sum_{k=j}^{n-1} \delta_{(k)}\tilde{d}_{(k)}$. Taking this to the original
coordinate system, we can write $e_{(j)}=\sum_{k=j}^{n-1} \delta_{(k)}d_{(k)}$. Multiply both sides by
$d_{(i)}^TA$ where $i<j$, we have
\begin{align*}
	d_{(i)}^TAe_{(j)} &= \sum_{k=j}^{n-1} \delta_{(k)} d_{(i)}^TAd_{(k)} = 0\tag{6.1}
\end{align*}
since $d_{(i)}$'s are $A$-orthogonal, $d_{(i)}^TAd_{(k)}$ unless $i=k$. Thus we have
\begin{align*}
	&d_{(i)}^TAe_{(j)} = 0~~ (i<j) \\
	&d_{(i)}^T r_{(j)} = 0 ~~(i<j) \tag{6.2}	\\
	&\therefore r_{(j)} \perp \text{span}\{d_{(0)}, d_{(1)}, d_{(i)} | i<j\} \\
	& \because \text{span}\{d_{(0)}, d_{(1)}, ..., d_{(i)}\} = \text{span}\{r_{0}, ..., r_{(i)}\} ~\forall ~ i \\
	& \therefore r_{(j)} \perp \text{span}\{r_{(0)}, r_{(1)}, r_{(i)} | i<j\} \\
	& \therefore r_{(i)} \perp r_{(j)} \forall ~~i \neq j	\tag{6.3}
\end{align*}

Once we have $r_{(i)} \perp r_{(j)} \forall ~~i \neq j$, we multiply both sides of (5.2) by $r_{(i)}^T$, we have
\begin{align*}
	r_{(i)}^Td_{(i)} &= r_{(i)}^Tr_{(i)} - \sum_{j=0}^{i-1} \frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}r_{(i)}^Td_{(j)} \\
	r_{(i)}^Td_{(i)} &= r_{(i)}^Tr_{(i)} \tag{6.4}
\end{align*}


\section{Finally, Conjugate Gradient Method}
In the previous section, by exploiting the spaces spanned by $d_{(i)}$'s and by $r_{(i)}$'s, we get two identities
\begin{itemize}
	\item[1.] $r_{(i)} \perp r_{(j)} \forall ~~i \neq j$
	\item[2.] $r_{(i)}^Td_{(i)} = r_{(i)}^Tr_{(i)}$
\end{itemize}
Now we use these two identities to simplify the "almost conjugate gradient method" and get the real conjugate 
gradient method. First, by the definition of the residual, we have
\begin{align*}
	r_{(j+1)} &=  b - Ax_{(j+1)} \\	
			&= b - A(x_{(j)} + \alpha_{(j)}d_{(j)}) \\
			&= r_{(j)} - \alpha_{(j)}Ad_{(j)}	\tag{7.1}
\end{align*}
and multiply both sides of (7.1) by $r_{(i)}^T$
\begin{align*}
	r_{(i)}^Tr_{(j+1)}	&= r_{(i)}^Tr_{(j)} - \alpha_{(j)}r_{(i)}^TAd_{(j)}	\\
	r_{(i)}^TAd_{(j)} 	&= \frac{r_{(i)}^Tr_{(j)} - r_{(i)}^Tr_{(j+1)}}{\alpha_{(j)}} \\
					&= 
					\begin{cases}
						\frac{r_{(i)}^Tr_{(i)}}{\alpha_{(i)}} & i=j \\
						-\frac{r_{(i)}^Tr_{(i)}}{\alpha_{(i-1)}} & i = j+1\\
						0	& \text{otherwise} \\
					\end{cases} \tag{7.2}
\end{align*}
From (7.2) we can see that most term in the summation in (5.2) appears to be 0. By plugging (7.2) into (5.2),
we get
\begin{align*}
	d_{(i)} &= r_{(i)} - \sum_{j=0}^{i-1} \frac{d_{(j)}^TAr_{(i)}}{d_{(j)}^TAd_{(j)}}d_{(j)} \\
		&= r_{(i)} - \frac{d_{(i-1)}^TAr_{(i)}}{d_{(i-1)}^TAd_{(i-1)}}d_{(i-1)} 	\\
		&= r_{(i)} + \frac{r_{(i)}^Tr_{(i)}}{\alpha_{(i-1)} d_{(i-1)}^TAd_{(i-1)}}d_{(i-1)} \tag{7.3}
\end{align*}
then plug (5.3) into (7.3), we get
\begin{align*}
	d_{(i)} &= r_{(i)} + \frac{r_{(i)}^Tr_{(i)}}{\alpha_{(i-1)} d_{(i-1)}^TAd_{(i-1)}}d_{(i-1)}	\\
			&= \frac{r_{(i)}^Tr_{(i)}}{r_{(i-1)}^Td_{(i-1)}}d_{(i-1)}	\\
			&= \frac{r_{(i)}^Tr_{(i)}}{r_{(i-1)}^Tr_{(i-1)}}d_{(i-1)}	\tag{7.4}
\end{align*}

Thus we have a new update rule for $d_{(i)}$ and $\alpha_{(i)}$, and this lead us to the final conjugate gradient 
method
\begin{align*}
	d_{(0)} 	&= r_{(0)}	\\
	d_{(i)} 	&= \frac{r_{(i)}^Tr_{(i)}}{r_{(i-1)}^Tr_{(i-1)}}d_{(i-1)} 	\tag{7.5}\\
	\alpha_{(i)} &= \frac{r_{(i)}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}} 	\tag{7.6}\\
	x_{(i+1)}	&= x_{(i)} + \alpha_{(i)} d_{(i)}	\tag{7.7}
\end{align*}
Notice that this is basicly the same thing as "almost conjugate gradient method" except that it save some 
computation by exploiting the innate property of the relationship between $d_{(i)}$'s and $r_{(i)}$'s.
Figure 10 and Figure 11 shows the progress of conjugate gradient method on the running example with different
starting points. They are almost the same as Figure 8 and Figure 9 because these two algorithms are doing 
the same thing, but with different complexities. 

\begin{figure}[h!]
\hspace{-8mm}
\includegraphics[scale=0.25]{cg1_transform.png}
\caption{Conjugate gradient method}
\end{figure}

\newpage
\begin{figure}[h!]
\hspace{-8mm}
\includegraphics[scale=0.25]{cg2_transform.png}
\caption{Conjugate gradient method}
\end{figure}

Figure 12 shows the residual norm and error norm when apply conjugate gradient method to a 2000x2000 
system.
\begin{figure}[h!]
\hspace{-4mm}
\includegraphics[scale=0.7]{cg_2000.png}
\caption{Conjugate gradient method on 2000x2000 system}
\end{figure}

\section{Conclusion}
This note is an introduction to the intricate conjugate gradient descent. Hopefully I have made the core idea 
of cg method clear to you. However, there are a lot of things about cg method that's not included in this note like
convergence analysis of cg method, krylov space and spectral analysis of krylov space
(this is available in EE364b). 

As you can see, the awfully complicated conjugate gradient method is just steepest descent with orthogonal
searching direction under a change of coordinate. Actually, not just in cg method,  the technique
of changing coordinate applies everywhere in math and engineering. Some examples are listed 
\begin{itemize}
	\item Gradient descent + changing coordinate (according to ellipsoid induced by local hessian) = 
	Newton's method
	\item Cutting-plane method + changing coordinate (according to ignorance ellipsoid) = Ellipsoid method
	\item Fourier transform of radial function+ changing coordinate (cartesian to polar) = Radon Transform
\end{itemize}

The philosophy here is, many sophisticated things are just combination of many simple but clever things. 
The hard part is what simples things to use and how to assemble these simple things. The problem 
itself, in some sense, gives you hints about what technique to use(changing coordinates, Fourier analysis, etc) 
with its properties
(asymmetric, radial, periodic, etc. In the example of cg method, we change coordinate to turn
an "inhomogeneous" ellipsoid into a "homogeneous" circle/ball). 
In other words, the solution to a problem is right there in the problem itself, however, identifying those 
properties and correctly applying those techniques to the very specific problem is where being smart 
makes the difference.

\end{document}



















