cleaning up splines
This commit is contained in:
mhjensen
+250
-257
@@ -615,7 +615,256 @@ pt.plot(it_array.T[0], it_array.T[1], "x-")
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Conjugate gradient =====
|
||||
===== 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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1456,262 +1705,6 @@ plt.show()
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Momentum based methods =====
|
||||
|
||||
|
||||
|
||||
!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
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Using gradient descent methods, limitations =====
|
||||
|
||||
Reference in New Issue
Block a user