mhjensen
542 lines
14 KiB
Plaintext
542 lines
14 KiB
Plaintext
TITLE: Data Analysis and Machine Learning Lectures: Cubic Splines and Gradient Methods
|
|
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
|
|
DATE: today
|
|
|
|
|
|
!split
|
|
===== Cubic Splines =====
|
|
!bblock
|
|
Cubic spline interpolation is among one of the most used
|
|
methods for interpolating between data points where the arguments
|
|
are organized as ascending series. In the library program we supply
|
|
such a function, based on the so-called cubic spline method to be
|
|
described below.
|
|
|
|
A spline function consists of polynomial pieces defined on
|
|
subintervals. The different subintervals are connected via
|
|
various continuity relations.
|
|
|
|
Assume we have at our disposal $n+1$ points $x_0, x_1, \dots x_n$
|
|
arranged so that $x_0 < x_1 < x_2 < \dots x_{n-1} < x_n$ (such points are called
|
|
knots). A spline function $s$ of degree $k$ with $n+1$ knots is defined
|
|
as follows
|
|
* On every subinterval $[x_{i-1},x_i)$ *s* is a polynomial of degree $\le k$.
|
|
* $s$ has $k-1$ continuous derivatives in the whole interval $[x_0,x_n]$.
|
|
!eblock
|
|
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
As an example, consider a spline function of degree $k=1$ defined as follows
|
|
!bt
|
|
\[
|
|
s(x)=\begin{bmatrix} s_0(x)=a_0x+b_0 & x\in [x_0, x_1) \\
|
|
s_1(x)=a_1x+b_1 & x\in [x_1, x_2) \\
|
|
\dots & \dots \\
|
|
s_{n-1}(x)=a_{n-1}x+b_{n-1} & x\in
|
|
[x_{n-1}, x_n] \end{bmatrix}.
|
|
\]
|
|
!et
|
|
In this case the polynomial consists of series of straight lines
|
|
connected to each other at every endpoint. The number of continuous
|
|
derivatives is then $k-1=0$, as expected when we deal with straight lines.
|
|
Such a polynomial is quite easy to construct given
|
|
$n+1$ points $x_0, x_1, \dots x_n$ and their corresponding
|
|
function values.
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
The most commonly used spline function is the one with $k=3$, the so-called
|
|
cubic spline function.
|
|
Assume that we have in adddition to the $n+1$ knots a series of
|
|
functions values $y_0=f(x_0), y_1=f(x_1), \dots y_n=f(x_n)$.
|
|
By definition, the polynomials $s_{i-1}$ and $s_i$
|
|
are thence supposed to interpolate the same point $i$, that is
|
|
!bt
|
|
\[
|
|
s_{i-1}(x_i)= y_i = s_i(x_i),
|
|
\]
|
|
!et
|
|
with $1 \le i \le n-1$. In total we have $n$ polynomials of the
|
|
type
|
|
!bt
|
|
\[
|
|
s_i(x)=a_{i0}+a_{i1}x+a_{i2}x^2+a_{i2}x^3,
|
|
\]
|
|
!et
|
|
yielding $4n$ coefficients to determine.
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
Every subinterval provides in addition the $2n$ conditions
|
|
!bt
|
|
\[
|
|
y_i = s(x_i),
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
s(x_{i+1})= y_{i+1},
|
|
\]
|
|
!et
|
|
to be fulfilled. If we also assume that $s'$ and $s''$ are continuous,
|
|
then
|
|
!bt
|
|
\[
|
|
s'_{i-1}(x_i)= s'_i(x_i),
|
|
\]
|
|
!et
|
|
yields $n-1$ conditions. Similarly,
|
|
!bt
|
|
\[
|
|
s''_{i-1}(x_i)= s''_i(x_i),
|
|
\]
|
|
!et
|
|
results in additional $n-1$ conditions. In total we have $4n$ coefficients
|
|
and $4n-2$ equations to determine them, leaving us with $2$ degrees of
|
|
freedom to be determined.
|
|
!eblock
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
Using the last equation we define two values for the second derivative, namely
|
|
!bt
|
|
\[
|
|
s''_{i}(x_i)= f_i,
|
|
\]
|
|
!et
|
|
and
|
|
!bt
|
|
\[
|
|
s''_{i}(x_{i+1})= f_{i+1},
|
|
\]
|
|
!et
|
|
and setting up a straight line between $f_i$ and $f_{i+1}$ we have
|
|
!bt
|
|
\[
|
|
s_i''(x) = \frac{f_i}{x_{i+1}-x_i}(x_{i+1}-x)+
|
|
\frac{f_{i+1}}{x_{i+1}-x_i}(x-x_i),
|
|
\]
|
|
!et
|
|
and integrating twice one obtains
|
|
!bt
|
|
\[
|
|
s_i(x) = \frac{f_i}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3+
|
|
\frac{f_{i+1}}{6(x_{i+1}-x_i)}(x-x_i)^3
|
|
+c(x-x_i)+d(x_{i+1}-x).
|
|
\]
|
|
!et
|
|
!eblock
|
|
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
Using the conditions $s_i(x_i)=y_i$ and $s_i(x_{i+1})=y_{i+1}$
|
|
we can in turn determine the constants $c$ and $d$ resulting in
|
|
!bt
|
|
\begin{align}
|
|
s_i(x) =&\frac{f_i}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3+
|
|
\frac{f_{i+1}}{6(x_{i+1}-x_i)}(x-x_i)^3 \nonumber \\
|
|
+&(\frac{y_{i+1}}{x_{i+1}-x_i}-\frac{f_{i+1}(x_{i+1}-x_i)}{6})
|
|
(x-x_i)+
|
|
(\frac{y_{i}}{x_{i+1}-x_i}-\frac{f_{i}(x_{i+1}-x_i)}{6})
|
|
(x_{i+1}-x).
|
|
\end{align}
|
|
!et
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
How to determine the values of the second
|
|
derivatives $f_{i}$ and $f_{i+1}$? We use the continuity assumption
|
|
of the first derivatives
|
|
!bt
|
|
\[
|
|
s'_{i-1}(x_i)= s'_i(x_i),
|
|
\]
|
|
!et
|
|
and set $x=x_i$. Defining $h_i=x_{i+1}-x_i$ we obtain finally
|
|
the following expression
|
|
!bt
|
|
\[
|
|
h_{i-1}f_{i-1}+2(h_{i}+h_{i-1})f_i+h_if_{i+1}=
|
|
\frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_{i}-y_{i-1}),
|
|
\]
|
|
!et
|
|
and introducing the shorthands $u_i=2(h_{i}+h_{i-1})$,
|
|
$v_i=\frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_{i}-y_{i-1})$,
|
|
we can reformulate the problem as a set of linear equations to be
|
|
solved through e.g., Gaussian elemination
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
Gaussian elimination
|
|
!bt
|
|
\[
|
|
\begin{bmatrix} u_1 & h_1 &0 &\dots & & & & \\
|
|
h_1 & u_2 & h_2 &0 &\dots & & & \\
|
|
0 & h_2 & u_3 & h_3 &0 &\dots & & \\
|
|
\dots& & \dots &\dots &\dots &\dots &\dots & \\
|
|
&\dots & & &0 &h_{n-3} &u_{n-2} &h_{n-2} \\
|
|
& && & &0 &h_{n-2} &u_{n-1} \end{bmatrix}
|
|
\begin{bmatrix} f_1 \\
|
|
f_2 \\
|
|
f_3\\
|
|
\dots \\
|
|
f_{n-2} \\
|
|
f_{n-1} \end{bmatrix} =
|
|
\begin{bmatrix} v_1 \\
|
|
v_2 \\
|
|
v_3\\
|
|
\dots \\
|
|
v_{n-2}\\
|
|
v_{n-1} \end{bmatrix}.
|
|
\]
|
|
!et
|
|
Note that this is a set of tridiagonal equations and can be solved
|
|
through only $O(n)$ operations.
|
|
!eblock
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
The functions supplied in the program library are *spline* and *splint*.
|
|
In order to use cubic spline interpolation you need first to call
|
|
!bc cppcod
|
|
spline(double x[], double y[], int n, double yp1, double yp2, double y2[])
|
|
!ec
|
|
This function takes as
|
|
input $x[0,..,n - 1]$ and $y[0,..,n - 1]$ containing a tabulation
|
|
$y_i = f(x_i)$ with $x_0 < x_1 < .. < x_{n - 1}$
|
|
together with the
|
|
first derivatives of $f(x)$ at $x_0$ and $x_{n-1}$, respectively. Then the
|
|
function returns $y2[0,..,n-1]$ which contains the second derivatives of
|
|
$f(x_i)$ at each point $x_i$. $n$ is the number of points.
|
|
This function provides the cubic spline interpolation for all subintervals
|
|
and is called only once.
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Splines =====
|
|
!bblock
|
|
Thereafter, if you wish to make various interpolations, you need to call the function
|
|
!bc cppcod
|
|
splint(double x[], double y[], double y2a[], int n, double x, double *y)
|
|
!ec
|
|
which takes as input
|
|
the tabulated values $x[0,..,n - 1]$ and $y[0,..,n - 1]$ and the output
|
|
y2a[0,..,n - 1] from *spline*. It returns the value $y$ corresponding
|
|
to the point $x$.
|
|
!eblock
|
|
|
|
|
|
|
|
!split
|
|
===== Conjugate gradient (CG) method =====
|
|
!bblock
|
|
The success of the CG method 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$.
|
|
!eblock
|
|
|
|
|
|
!split
|
|
===== Conjugate gradient method =====
|
|
!bblock
|
|
|
|
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.
|
|
If we search for a minimum of the quantum mechanical variance, then the matrix
|
|
$\hat{A}$, which is called the Hessian, is given by the second-derivative of the function we want to minimize. This quantity is always positive definite. In our case this corresponds normally to the second derivative of the energy.
|
|
!eblock
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
===== Conjugate gradient method, Newton's method first =====
|
|
!bblock
|
|
We seek the minimum of the energy or the variance as function of various variational parameters.
|
|
In our case we have thus a function $f$ whose minimum we are seeking.
|
|
In Newton's method we set $\nabla f = 0$ and we can thus compute the next iteration point
|
|
!bt
|
|
\begin{equation*}
|
|
\hat{x}-\hat{x}_i=\hat{A}^{-1}\nabla f(\hat{x}_i).
|
|
\end{equation*}
|
|
!et
|
|
Subtracting this equation from that of $\hat{x}_{i+1}$ we have
|
|
!bt
|
|
\begin{equation*}
|
|
\hat{x}_{i+1}-\hat{x}_i=\hat{A}^{-1}(\nabla f(\hat{x}_{i+1})-\nabla f(\hat{x}_i)).
|
|
\end{equation*}
|
|
!et
|
|
!eblock
|
|
|
|
!split
|
|
===== Simple example and demonstration =====
|
|
!bblock
|
|
The function $f$ can be either the energy or the variance. If we choose the energy then we have
|
|
!bt
|
|
\begin{equation*}
|
|
\hat{\alpha}_{i+1}-\hat{\alpha}_i=\hat{A}^{-1}(\nabla E(\hat{\alpha}_{i+1})-\nabla E(\hat{\alpha}_i)).
|
|
\end{equation*}
|
|
!et
|
|
In the simple harmonic oscillator model, the gradient and the Hessian $\hat{A}$ are
|
|
!bt
|
|
\begin{equation*}
|
|
\frac{d\langle E_L[\alpha]\rangle}{d\alpha} = \alpha-\frac{1}{4\alpha^3}
|
|
\end{equation*}
|
|
!et
|
|
and a second derivative which is always positive (meaning that we find a minimum)
|
|
!bt
|
|
\begin{equation*}
|
|
\hat{A}= \frac{d^2\langle E_L[\alpha]\rangle}{d\alpha^2} = 1+\frac{3}{4\alpha^4}
|
|
\end{equation*}
|
|
!et
|
|
!eblock
|
|
|
|
!split
|
|
===== Simple example and demonstration =====
|
|
!bblock
|
|
We get then
|
|
!bt
|
|
\begin{equation*}
|
|
\alpha_{i+1}=\frac{4}{3}\alpha_i-\frac{\alpha_i^4}{3\alpha_{i+1}^3},
|
|
\end{equation*}
|
|
!et
|
|
which can be rewritten as
|
|
!bt
|
|
\begin{equation*}
|
|
\alpha_{i+1}^4-\frac{4}{3}\alpha_i\alpha_{i+1}^4+\frac{1}{3}\alpha_i^4.
|
|
\end{equation*}
|
|
!et
|
|
!eblock
|
|
|
|
!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
|
|
|
|
|