lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updating week39

Browse Source
This commit is contained in:
mhjensen
2020-09-24 06:54:17 +02:00
parent 9694755694
commit 6df59f6e6d
78 changed files with 11412 additions and 13885 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/src/week39/week39.do.txt
+371 -500
View File
@@ -386,6 +386,372 @@ o A norm is any function that satisfy the following properties
Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).
!split
===== Standard steepest descent =====
Before we proceed, we would like to discuss the approach called the
_standard Steepest descent_, which again leads to us having to be able
to compute a matrix. It belongs to the class of Conjugate Gradient methods (CG).
"The success of the CG method":"https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf"
for finding solutions of non-linear problems is based on the theory
of conjugate gradients for linear systems of equations. It belongs to
the class of iterative methods for solving problems from linear
algebra of the type
!bt
\begin{equation*}
\hat{A}\hat{x} = \hat{b}.
\end{equation*}
!et
In the iterative process we end up with a problem like
!bt
\begin{equation*}
\hat{r}= \hat{b}-\hat{A}\hat{x},
\end{equation*}
!et
where $\hat{r}$ is the so-called residual or error in the iterative process.
When we have found the exact solution, $\hat{r}=0$.
!split
===== Gradient method =====
The residual is zero when we reach the minimum of the quadratic equation
!bt
\begin{equation*}
P(\hat{x})=\frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T\hat{b},
\end{equation*}
!et
with the constraint that the matrix $\hat{A}$ is positive definite and
symmetric. This defines also the Hessian and we want it to be positive definite.
!split
===== Steepest descent method =====
We denote the initial guess for $\hat{x}$ as $\hat{x}_0$.
We can assume without loss of generality that
!bt
\begin{equation*}
\hat{x}_0=0,
\end{equation*}
!et
or consider the system
!bt
\begin{equation*}
\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
\end{equation*}
!et
instead.
!split
===== Steepest descent method =====
!bblock
One can show that the solution $\hat{x}$ is also the unique minimizer of the quadratic form
!bt
\begin{equation*}
f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
\end{equation*}
!et
This suggests taking the first basis vector $\hat{r}_1$ (see below for definition)
to be the gradient of $f$ at $\hat{x}=\hat{x}_0$,
which equals
!bt
\begin{equation*}
\hat{A}\hat{x}_0-\hat{b},
\end{equation*}
!et
and
$\hat{x}_0=0$ it is equal $-\hat{b}$.
!eblock
!split
===== Final expressions =====
!bblock
We can compute the residual iteratively as
!bt
\begin{equation*}
\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
\end{equation*}
!et
which equals
!bt
\begin{equation*}
\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{r}_k),
\end{equation*}
!et
or
!bt
\begin{equation*}
(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{r}_k,
\end{equation*}
!et
which gives
!bt
\[
\alpha_k = \frac{\hat{r}_k^T\hat{r}_k}{\hat{r}_k^T\hat{A}\hat{r}_k}
\]
!et
leading to the iterative scheme
!bt
\begin{equation*}
\hat{x}_{k+1}=\hat{x}_k-\alpha_k\hat{r}_{k},
\end{equation*}
!et
!eblock
!split
===== Steepest descent example =====
!bc pycod
import numpy as np
import numpy.linalg as la
import scipy.optimize as sopt
import matplotlib.pyplot as pt
from mpl_toolkits.mplot3d import axes3d
def f(x):
return 0.5*x[0]**2 + 2.5*x[1]**2
def df(x):
return np.array([x[0], 5*x[1]])
fig = pt.figure()
ax = fig.gca(projection="3d")
xmesh, ymesh = np.mgrid[-2:2:50j,-2:2:50j]
fmesh = f(np.array([xmesh, ymesh]))
ax.plot_surface(xmesh, ymesh, fmesh)
!ec
And then as countor plot
!bc pycod
pt.axis("equal")
pt.contour(xmesh, ymesh, fmesh)
guesses = [np.array([2, 2./5])]
!ec
Find guesses
!bc pycod
x = guesses[-1]
s = -df(x)
!ec
Run it!
!bc pycod
def f1d(alpha):
return f(x + alpha*s)
alpha_opt = sopt.golden(f1d)
next_guess = x + alpha_opt * s
guesses.append(next_guess)
print(next_guess)
!ec
What happened?
!bc pycod
pt.axis("equal")
pt.contour(xmesh, ymesh, fmesh, 50)
it_array = np.array(guesses)
pt.plot(it_array.T[0], it_array.T[1], "x-")
!ec
!split
===== Conjugate gradient method =====
!bblock
In the CG method we define so-called conjugate directions and two vectors
$\hat{s}$ and $\hat{t}$
are said to be
conjugate if
!bt
\begin{equation*}
\hat{s}^T\hat{A}\hat{t}= 0.
\end{equation*}
!et
The philosophy of the CG method is to perform searches in various conjugate directions
of our vectors $\hat{x}_i$ obeying the above criterion, namely
!bt
\begin{equation*}
\hat{x}_i^T\hat{A}\hat{x}_j= 0.
\end{equation*}
!et
Two vectors are conjugate if they are orthogonal with respect to
this inner product. Being conjugate is a symmetric relation: if $\hat{s}$ is conjugate to $\hat{t}$, then $\hat{t}$ is conjugate to $\hat{s}$.
!eblock
!split
===== Conjugate gradient method =====
!bblock
An example is given by the eigenvectors of the matrix
!bt
\begin{equation*}
\hat{v}_i^T\hat{A}\hat{v}_j= \lambda\hat{v}_i^T\hat{v}_j,
\end{equation*}
!et
which is zero unless $i=j$.
!eblock
!split
===== Conjugate gradient method =====
!bblock
Assume now that we have a symmetric positive-definite matrix $\hat{A}$ of size
$n\times n$. At each iteration $i+1$ we obtain the conjugate direction of a vector
!bt
\begin{equation*}
\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
\end{equation*}
!et
We assume that $\hat{p}_{i}$ is a sequence of $n$ mutually conjugate directions.
Then the $\hat{p}_{i}$ form a basis of $R^n$ and we can expand the solution
$ \hat{A}\hat{x} = \hat{b}$ in this basis, namely
!bt
\begin{equation*}
\hat{x} = \sum^{n}_{i=1} \alpha_i \hat{p}_i.
\end{equation*}
!et
!eblock
!split
===== Conjugate gradient method =====
!bblock
The coefficients are given by
!bt
\begin{equation*}
\mathbf{A}\mathbf{x} = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.
\end{equation*}
!et
Multiplying with $\hat{p}_k^T$ from the left gives
!bt
\begin{equation*}
\hat{p}_k^T \hat{A}\hat{x} = \sum^{n}_{i=1} \alpha_i\hat{p}_k^T \hat{A}\hat{p}_i= \hat{p}_k^T \hat{b},
\end{equation*}
!et
and we can define the coefficients $\alpha_k$ as
!bt
\begin{equation*}
\alpha_k = \frac{\hat{p}_k^T \hat{b}}{\hat{p}_k^T \hat{A} \hat{p}_k}
\end{equation*}
!et
!eblock
!split
===== Conjugate gradient method and iterations =====
!bblock
If we choose the conjugate vectors $\hat{p}_k$ carefully,
then we may not need all of them to obtain a good approximation to the solution
$\hat{x}$.
We want to regard the conjugate gradient method as an iterative method.
This will us to solve systems where $n$ is so large that the direct
method would take too much time.
We denote the initial guess for $\hat{x}$ as $\hat{x}_0$.
We can assume without loss of generality that
!bt
\begin{equation*}
\hat{x}_0=0,
\end{equation*}
!et
or consider the system
!bt
\begin{equation*}
\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
\end{equation*}
!et
instead.
!eblock
!split
===== Conjugate gradient method =====
!bblock
One can show that the solution $\hat{x}$ is also the unique minimizer of the quadratic form
!bt
\begin{equation*}
f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
\end{equation*}
!et
This suggests taking the first basis vector $\hat{p}_1$
to be the gradient of $f$ at $\hat{x}=\hat{x}_0$,
which equals
!bt
\begin{equation*}
\hat{A}\hat{x}_0-\hat{b},
\end{equation*}
!et
and
$\hat{x}_0=0$ it is equal $-\hat{b}$.
The other vectors in the basis will be conjugate to the gradient,
hence the name conjugate gradient method.
!eblock
!split
===== Conjugate gradient method =====
!bblock
Let $\hat{r}_k$ be the residual at the $k$-th step:
!bt
\begin{equation*}
\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
\end{equation*}
!et
Note that $\hat{r}_k$ is the negative gradient of $f$ at
$\hat{x}=\hat{x}_k$,
so the gradient descent method would be to move in the direction $\hat{r}_k$.
Here, we insist that the directions $\hat{p}_k$ are conjugate to each other,
so we take the direction closest to the gradient $\hat{r}_k$
under the conjugacy constraint.
This gives the following expression
!bt
\begin{equation*}
\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
\end{equation*}
!et
!eblock
!split
===== Conjugate gradient method =====
!bblock
We can also compute the residual iteratively as
!bt
\begin{equation*}
\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
\end{equation*}
!et
which equals
!bt
\begin{equation*}
\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
\end{equation*}
!et
or
!bt
\begin{equation*}
(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
\end{equation*}
!et
which gives
!bt
\begin{equation*}
\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
\end{equation*}
!et
!eblock
!split
@@ -429,9 +795,9 @@ It is convenient to write $\mathbf{\hat{y}} = X\beta$ where $X \in \mathbb{R}^{1
!bt
\[
X \equiv \begin{bmatrix}
1 & x_1 \\
\vdots & \vdots \\
1 & x_{100} & \\
1 &; x_1 \\
\vdots &; \vdots \\
1 &; x_{100} &; \\
\end{bmatrix}.
\]
!et
@@ -462,8 +828,8 @@ The Hessian matrix of $C(\beta)$ is given by
!bt
\[
\hat{H} \equiv \begin{bmatrix}
\frac{\partial^2 C(\beta)}{\partial \beta_0^2} & \frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} \\
\frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} & \frac{\partial^2 C(\beta)}{\partial \beta_1^2} & \\
\frac{\partial^2 C(\beta)}{\partial \beta_0^2} &; \frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} \\
\frac{\partial^2 C(\beta)}{\partial \beta_0 \partial \beta_1} &; \frac{\partial^2 C(\beta)}{\partial \beta_1^2} &; \\
\end{bmatrix} = 2X^T X.
\]
!et
@@ -1497,498 +1863,3 @@ a /=b
!ec
!split
===== Standard steepest descent =====
Before we proceed, we would like to discuss the approach called the
_standard Steepest descent_, which again leads to us having to be able
to compute a matrix. It belongs to the class of Conjugate Gradient methods (CG).
"The success of the CG method":"https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf"
for finding solutions of non-linear problems is based on the theory
of conjugate gradients for linear systems of equations. It belongs to
the class of iterative methods for solving problems from linear
algebra of the type
!bt
\begin{equation*}
\hat{A}\hat{x} = \hat{b}.
\end{equation*}
!et
In the iterative process we end up with a problem like
!bt
\begin{equation*}
\hat{r}= \hat{b}-\hat{A}\hat{x},
\end{equation*}
!et
where $\hat{r}$ is the so-called residual or error in the iterative process.
When we have found the exact solution, $\hat{r}=0$.
!split
===== Gradient method =====
The residual is zero when we reach the minimum of the quadratic equation
!bt
\begin{equation*}
P(\hat{x})=\frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T\hat{b},
\end{equation*}
!et
with the constraint that the matrix $\hat{A}$ is positive definite and
symmetric. This defines also the Hessian and we want it to be positive definite.
!split
===== Steepest descent method =====
We denote the initial guess for $\hat{x}$ as $\hat{x}_0$.
We can assume without loss of generality that
!bt
\begin{equation*}
\hat{x}_0=0,
\end{equation*}
!et
or consider the system
!bt
\begin{equation*}
\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
\end{equation*}
!et
instead.
!split
===== Steepest descent method =====
!bblock
One can show that the solution $\hat{x}$ is also the unique minimizer of the quadratic form
!bt
\begin{equation*}
f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
\end{equation*}
!et
This suggests taking the first basis vector $\hat{r}_1$ (see below for definition)
to be the gradient of $f$ at $\hat{x}=\hat{x}_0$,
which equals
!bt
\begin{equation*}
\hat{A}\hat{x}_0-\hat{b},
\end{equation*}
!et
and
$\hat{x}_0=0$ it is equal $-\hat{b}$.
!eblock
!split
===== Final expressions =====
!bblock
We can compute the residual iteratively as
!bt
\begin{equation*}
\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
\end{equation*}
!et
which equals
!bt
\begin{equation*}
\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{r}_k),
\end{equation*}
!et
or
!bt
\begin{equation*}
(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{r}_k,
\end{equation*}
!et
which gives
!bt
\[
\alpha_k = \frac{\hat{r}_k^T\hat{r}_k}{\hat{r}_k^T\hat{A}\hat{r}_k}
\]
!et
leading to the iterative scheme
!bt
\begin{equation*}
\hat{x}_{k+1}=\hat{x}_k-\alpha_k\hat{r}_{k},
\end{equation*}
!et
!eblock
!split
===== Code examples for steepest descent =====
!split
===== Simple codes for steepest descent and conjugate gradient using a $2\times 2$ matrix, in c++, Python code to come =====
!bblock
!bc cppcod
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "vectormatrixclass.h"
using namespace std;
// Main function begins here
int main(int argc, char * argv[]){
int dim = 2;
Vector x(dim),xsd(dim), b(dim),x0(dim);
Matrix A(dim,dim);
// Set our initial guess
x0(0) = x0(1) = 0;
// Set the matrix
A(0,0) = 3; A(1,0) = 2; A(0,1) = 2; A(1,1) = 6;
b(0) = 2; b(1) = -8;
cout << "The Matrix A that we are using: " << endl;
A.Print();
cout << endl;
xsd = SteepestDescent(A,b,x0);
cout << "The approximate solution using Steepest Descent is: " << endl;
xsd.Print();
cout << endl;
}
!ec
!eblock
!split
===== The routine for the steepest descent method =====
!bblock
!bc cppcod
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
int IterMax, i;
int dim = x0.Dimension();
const double tolerance = 1.0e-14;
Vector x(dim),f(dim),z(dim);
double c,alpha,d;
IterMax = 30;
x = x0;
r = A*x-b;
i = 0;
while (i <= IterMax){
z = A*r;
c = dot(r,r);
alpha = c/dot(r,z);
x = x - alpha*r;
r = A*x-b;
if(sqrt(dot(r,r)) < tolerance) break;
i++;
}
return x;
}
!ec
!eblock
!split
===== Steepest descent example =====
!bc pycod
import numpy as np
import numpy.linalg as la
import scipy.optimize as sopt
import matplotlib.pyplot as pt
from mpl_toolkits.mplot3d import axes3d
def f(x):
return 0.5*x[0]**2 + 2.5*x[1]**2
def df(x):
return np.array([x[0], 5*x[1]])
fig = pt.figure()
ax = fig.gca(projection="3d")
xmesh, ymesh = np.mgrid[-2:2:50j,-2:2:50j]
fmesh = f(np.array([xmesh, ymesh]))
ax.plot_surface(xmesh, ymesh, fmesh)
!ec
And then as countor plot
!bc pycod
pt.axis("equal")
pt.contour(xmesh, ymesh, fmesh)
guesses = [np.array([2, 2./5])]
!ec
Find guesses
!bc pycod
x = guesses[-1]
s = -df(x)
!ec
Run it!
!bc pycod
def f1d(alpha):
return f(x + alpha*s)
alpha_opt = sopt.golden(f1d)
next_guess = x + alpha_opt * s
guesses.append(next_guess)
print(next_guess)
!ec
What happened?
!bc pycod
pt.axis("equal")
pt.contour(xmesh, ymesh, fmesh, 50)
it_array = np.array(guesses)
pt.plot(it_array.T[0], it_array.T[1], "x-")
!ec
!split
===== Conjugate gradient method =====
!bblock
In the CG method we define so-called conjugate directions and two vectors
$\hat{s}$ and $\hat{t}$
are said to be
conjugate if
!bt
\begin{equation*}
\hat{s}^T\hat{A}\hat{t}= 0.
\end{equation*}
!et
The philosophy of the CG method is to perform searches in various conjugate directions
of our vectors $\hat{x}_i$ obeying the above criterion, namely
!bt
\begin{equation*}
\hat{x}_i^T\hat{A}\hat{x}_j= 0.
\end{equation*}
!et
Two vectors are conjugate if they are orthogonal with respect to
this inner product. Being conjugate is a symmetric relation: if $\hat{s}$ is conjugate to $\hat{t}$, then $\hat{t}$ is conjugate to $\hat{s}$.
!eblock
!split
===== Conjugate gradient method =====
!bblock
An example is given by the eigenvectors of the matrix
!bt
\begin{equation*}
\hat{v}_i^T\hat{A}\hat{v}_j= \lambda\hat{v}_i^T\hat{v}_j,
\end{equation*}
!et
which is zero unless $i=j$.
!eblock
!split
===== Conjugate gradient method =====
!bblock
Assume now that we have a symmetric positive-definite matrix $\hat{A}$ of size
$n\times n$. At each iteration $i+1$ we obtain the conjugate direction of a vector
!bt
\begin{equation*}
\hat{x}_{i+1}=\hat{x}_{i}+\alpha_i\hat{p}_{i}.
\end{equation*}
!et
We assume that $\hat{p}_{i}$ is a sequence of $n$ mutually conjugate directions.
Then the $\hat{p}_{i}$ form a basis of $R^n$ and we can expand the solution
$ \hat{A}\hat{x} = \hat{b}$ in this basis, namely
!bt
\begin{equation*}
\hat{x} = \sum^{n}_{i=1} \alpha_i \hat{p}_i.
\end{equation*}
!et
!eblock
!split
===== Conjugate gradient method =====
!bblock
The coefficients are given by
!bt
\begin{equation*}
\mathbf{A}\mathbf{x} = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.
\end{equation*}
!et
Multiplying with $\hat{p}_k^T$ from the left gives
!bt
\begin{equation*}
\hat{p}_k^T \hat{A}\hat{x} = \sum^{n}_{i=1} \alpha_i\hat{p}_k^T \hat{A}\hat{p}_i= \hat{p}_k^T \hat{b},
\end{equation*}
!et
and we can define the coefficients $\alpha_k$ as
!bt
\begin{equation*}
\alpha_k = \frac{\hat{p}_k^T \hat{b}}{\hat{p}_k^T \hat{A} \hat{p}_k}
\end{equation*}
!et
!eblock
!split
===== Conjugate gradient method and iterations =====
!bblock
If we choose the conjugate vectors $\hat{p}_k$ carefully,
then we may not need all of them to obtain a good approximation to the solution
$\hat{x}$.
We want to regard the conjugate gradient method as an iterative method.
This will us to solve systems where $n$ is so large that the direct
method would take too much time.
We denote the initial guess for $\hat{x}$ as $\hat{x}_0$.
We can assume without loss of generality that
!bt
\begin{equation*}
\hat{x}_0=0,
\end{equation*}
!et
or consider the system
!bt
\begin{equation*}
\hat{A}\hat{z} = \hat{b}-\hat{A}\hat{x}_0,
\end{equation*}
!et
instead.
!eblock
!split
===== Conjugate gradient method =====
!bblock
One can show that the solution $\hat{x}$ is also the unique minimizer of the quadratic form
!bt
\begin{equation*}
f(\hat{x}) = \frac{1}{2}\hat{x}^T\hat{A}\hat{x} - \hat{x}^T \hat{x} , \quad \hat{x}\in\mathbf{R}^n.
\end{equation*}
!et
This suggests taking the first basis vector $\hat{p}_1$
to be the gradient of $f$ at $\hat{x}=\hat{x}_0$,
which equals
!bt
\begin{equation*}
\hat{A}\hat{x}_0-\hat{b},
\end{equation*}
!et
and
$\hat{x}_0=0$ it is equal $-\hat{b}$.
The other vectors in the basis will be conjugate to the gradient,
hence the name conjugate gradient method.
!eblock
!split
===== Conjugate gradient method =====
!bblock
Let $\hat{r}_k$ be the residual at the $k$-th step:
!bt
\begin{equation*}
\hat{r}_k=\hat{b}-\hat{A}\hat{x}_k.
\end{equation*}
!et
Note that $\hat{r}_k$ is the negative gradient of $f$ at
$\hat{x}=\hat{x}_k$,
so the gradient descent method would be to move in the direction $\hat{r}_k$.
Here, we insist that the directions $\hat{p}_k$ are conjugate to each other,
so we take the direction closest to the gradient $\hat{r}_k$
under the conjugacy constraint.
This gives the following expression
!bt
\begin{equation*}
\hat{p}_{k+1}=\hat{r}_k-\frac{\hat{p}_k^T \hat{A}\hat{r}_k}{\hat{p}_k^T\hat{A}\hat{p}_k} \hat{p}_k.
\end{equation*}
!et
!eblock
!split
===== Conjugate gradient method =====
!bblock
We can also compute the residual iteratively as
!bt
\begin{equation*}
\hat{r}_{k+1}=\hat{b}-\hat{A}\hat{x}_{k+1},
\end{equation*}
!et
which equals
!bt
\begin{equation*}
\hat{b}-\hat{A}(\hat{x}_k+\alpha_k\hat{p}_k),
\end{equation*}
!et
or
!bt
\begin{equation*}
(\hat{b}-\hat{A}\hat{x}_k)-\alpha_k\hat{A}\hat{p}_k,
\end{equation*}
!et
which gives
!bt
\begin{equation*}
\hat{r}_{k+1}=\hat{r}_k-\hat{A}\hat{p}_{k},
\end{equation*}
!et
!eblock
!split
===== Simple implementation of the Conjugate gradient algorithm =====
!bblock
!bc cppcod
Vector ConjugateGradient(Matrix A, Vector b, Vector x0){
int dim = x0.Dimension();
const double tolerance = 1.0e-14;
Vector x(dim),r(dim),v(dim),z(dim);
double c,t,d;
x = x0;
r = b - A*x;
v = r;
c = dot(r,r);
int i = 0; IterMax = dim;
while(i <= IterMax){
z = A*v;
t = c/dot(v,z);
x = x + t*v;
r = r - t*z;
d = dot(r,r);
if(sqrt(d) < tolerance)
break;
v = r + (d/c)*v;
c = d; i++;
}
return x;
}
!ec
!eblock
!split
===== Broyden–Fletcher–Goldfarb–Shanno algorithm =====
!bblock
The optimization problem is to minimize $f(\mathbf {x} )$ where $\mathbf {x}$ is a vector in $R^{n}$, and $f$ is a differentiable scalar function. There are no constraints on the values that $\mathbf {x}$ can take.
The algorithm begins at an initial estimate for the optimal value $\mathbf {x}_{0}$ and proceeds iteratively to get a better estimate at each stage.
The search direction $p_k$ at stage $k$ is given by the solution of the analogue of the Newton equation
!bt
\[
B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x}_{k}),
\]
!et
where $B_{k}$ is an approximation to the Hessian matrix, which is
updated iteratively at each stage, and $\nabla f(\mathbf {x} _{k})$
is the gradient of the function
evaluated at $x_k$.
A line search in the direction $p_k$ is then used to
find the next point $x_{k+1}$ by minimising
!bt
\[
f(\mathbf {x}_{k}+\alpha \mathbf {p}_{k}),
\]
!et
over the scalar $\alpha > 0$.
!eblock