Editing Conjugate gradient method (section)

===Numerical example===

Consider the linear system '''Ax''' = '''b''' given by

:<math>\mathbf{A} \mathbf{x}= \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix},</math>

we will perform two steps of the conjugate gradient method beginning with the initial guess

:<math>\mathbf{x}_0 = \begin{bmatrix} 2 \\ 1 \end{bmatrix}</math>

in order to find an approximate solution to the system.

====Solution====
For reference, the exact solution is 
:<math> \mathbf{x} = \begin{bmatrix} \frac{1}{11} \\\\ \frac{7}{11} \end{bmatrix} \approx \begin{bmatrix} 0.0909 \\\\ 0.6364 \end{bmatrix}</math>
Our first step is to calculate the residual vector '''r'''<sub>0</sub> associated with '''x'''<sub>0</sub>.  This residual is computed from the formula '''r'''<sub>0</sub> = '''b''' - '''Ax'''<sub>0</sub>, and in our case is equal to

:<math>\mathbf{r}_0 = \begin{bmatrix} 1 \\ 2 \end{bmatrix} - 
\begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix}
\begin{bmatrix} 2 \\ 1 \end{bmatrix} = 
\begin{bmatrix}-8 \\ -3 \end{bmatrix} = \mathbf{p}_0.</math>

Since this is the first iteration, we will use the residual vector '''r'''<sub>0</sub> as our initial search direction '''p'''<sub>0</sub>; the method of selecting '''p'''<sub>''k''</sub> will change in further iterations.

We now compute the scalar {{math|''α''<sub>0</sub>}} using the relationship

:<math> \alpha_0 = \frac{\mathbf{r}_0^\mathsf{T} \mathbf{r}_0}{\mathbf{p}_0^\mathsf{T} \mathbf{A p}_0} = \frac{\begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix}}{  \begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix}  } =\frac{73}{331}\approx0.2205</math>

We can now compute '''x'''<sub>1</sub> using the formula

:<math>\mathbf{x}_1 = \mathbf{x}_0 + \alpha_0\mathbf{p}_0 = \begin{bmatrix} 2 \\ 1 \end{bmatrix} + \frac{73}{331} \begin{bmatrix} -8 \\ -3 \end{bmatrix} \approx \begin{bmatrix} 0.2356 \\ 0.3384 \end{bmatrix}.</math>

This result completes the first iteration, the result being an "improved" approximate solution to the system, '''x'''<sub>1</sub>. We may now move on and compute the next residual vector '''r'''<sub>1</sub> using the formula

:<math>\mathbf{r}_1 = \mathbf{r}_0 - \alpha_0 \mathbf{A} \mathbf{p}_0 = \begin{bmatrix} -8 \\ -3 \end{bmatrix} - \frac{73}{331} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix} \approx \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}.</math>

Our next step in the process is to compute the scalar {{math|''β''<sub>0</sub>}} that will eventually be used to determine the next search direction '''p'''<sub>1</sub>.

:<math>\beta_0 = \frac{\mathbf{r}_1^\mathsf{T} \mathbf{r}_1}{\mathbf{r}_0^\mathsf{T} \mathbf{r}_0} \approx \frac{\begin{bmatrix} -0.2810 & 0.7492 \end{bmatrix} \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}}{\begin{bmatrix} -8 & -3 \end{bmatrix} \begin{bmatrix} -8 \\ -3 \end{bmatrix}} = 0.0088.</math>

Now, using this scalar {{math|''β''<sub>0</sub>}}, we can compute the next search direction '''p'''<sub>1</sub> using the relationship

:<math>\mathbf{p}_1 = \mathbf{r}_1 + \beta_0 \mathbf{p}_0 \approx \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix} + 0.0088 \begin{bmatrix} -8 \\ -3 \end{bmatrix} = \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix}.</math>

We now compute the scalar {{math|''α''<sub>1</sub>}} using our newly acquired '''p'''<sub>1</sub> using the same method as that used for {{math|''α''<sub>0</sub>}}.

:<math> \alpha_1 = \frac{\mathbf{r}_1^\mathsf{T} \mathbf{r}_1}{\mathbf{p}_1^\mathsf{T} \mathbf{A p}_1} \approx \frac{\begin{bmatrix} -0.2810 & 0.7492 \end{bmatrix} \begin{bmatrix} -0.2810 \\ 0.7492 \end{bmatrix}}{  \begin{bmatrix} -0.3511 & 0.7229 \end{bmatrix} \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix}  } = 0.4122.</math>

Finally, we find '''x'''<sub>2</sub> using the same method as that used to find '''x'''<sub>1</sub>.

:<math>\mathbf{x}_2 = \mathbf{x}_1 + \alpha_1 \mathbf{p}_1 \approx \begin{bmatrix} 0.2356 \\ 0.3384 \end{bmatrix} + 0.4122 \begin{bmatrix} -0.3511 \\ 0.7229 \end{bmatrix} = \begin{bmatrix} 0.0909 \\ 0.6364 \end{bmatrix}.</math>

The result, '''x'''<sub>2</sub>, is a "better" approximation to the system's solution than '''x'''<sub>1</sub> and '''x'''<sub>0</sub>.  If exact arithmetic were to be used in this example instead of limited-precision, then the exact solution would theoretically have been reached after ''n'' = 2 iterations (''n'' being the order of the system).