Update on steepest descent
This commit is contained in:
mhjensen
@@ -15,6 +15,203 @@ the cost function. Ideally we would be able to solve for $\beta$
|
||||
analytically, however this is not possible in general and we must use
|
||||
some approximative/numerical method to compute the minimum.
|
||||
|
||||
|
||||
!split
|
||||
===== Revisiting our Logistic Regression case =====
|
||||
|
||||
In our discussion on Logistic Regression we defined we studied first the
|
||||
case of
|
||||
two classes, with $y_i$ either
|
||||
$0$ or $1$. Furthermore we assumed also that we have only two
|
||||
parameters $\beta$ in our fitting of the Sigmoid function, that is we
|
||||
defined probabilities
|
||||
|
||||
!bt
|
||||
\begin{align*}
|
||||
p(y_i=1|x_i,\hat{\beta}) &= \frac{\exp{(\beta_0+\beta_1x_i)}}{1+\exp{(\beta_0+\beta_1x_i)}},\nonumber\\
|
||||
p(y_i=0|x_i,\hat{\beta}) &= 1 - p(y_i=1|x_i,\hat{\beta}),
|
||||
\end{align*}
|
||||
!et
|
||||
where $\hat{\beta}$ are the weights we wish to extract from data, in our case $\beta_0$ and $\beta_1$.
|
||||
|
||||
!split
|
||||
===== The equations to solve =====
|
||||
|
||||
Our compact equations used a definition of a vector $\hat{y}$ with $n$
|
||||
elements $y_i$, an $n\times p$ matrix $\hat{X}$ which contains the
|
||||
$x_i$ values and a vector $\hat{p}$ of fitted probabilities
|
||||
$p(y_i\vert x_i,\hat{\beta})$. We rewrote in a more compact form
|
||||
the first derivative of the cost function as
|
||||
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}} = -\hat{X}^T\left(\hat{y}-\hat{p}\right).
|
||||
\]
|
||||
!et
|
||||
|
||||
If we in addition define a diagonal matrix $\hat{W}$ with elements
|
||||
$p(y_i\vert x_i,\hat{\beta})(1-p(y_i\vert x_i,\hat{\beta})$, we can obtain a compact expression of the second derivative as
|
||||
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T} = \hat{X}^T\hat{W}\hat{X}.
|
||||
\]
|
||||
!et
|
||||
This defines what we call the Hessian.
|
||||
|
||||
!split
|
||||
===== Solving using Newton-Raphson's method =====
|
||||
|
||||
If we can set up these equations, Newton-Raphson's iterative method is the nomrally the method of choice. It requires however that we setting the matrices that define the first and second derivatives.
|
||||
|
||||
Our iterative scheme is then given by
|
||||
|
||||
!bt
|
||||
\[
|
||||
\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\frac{\partial^2 \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}\partial \hat{\beta}^T}\right)^{-1}\times \left(\frac{\partial \mathcal{C}(\hat{\beta})}{\partial \hat{\beta}}\right)_{\hat{\beta}^{\mathrm{old}}},
|
||||
\]
|
||||
!et
|
||||
or in matrix form as
|
||||
|
||||
!bt
|
||||
\[
|
||||
\hat{\beta}^{\mathrm{new}} = \hat{\beta}^{\mathrm{old}}-\left(\hat{X}^T\hat{W}\hat{X} \right)^{-1}\times \left(-\hat{X}^T(\hat{y}-\hat{p}) \right)_{\hat{\beta}^{\mathrm{old}}}.
|
||||
\]
|
||||
!et
|
||||
The right-hand side is computed with the old values of $\beta$.
|
||||
|
||||
If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement.
|
||||
|
||||
|
||||
!split
|
||||
===== Brief reminder on Newton-Raphson's method =====
|
||||
|
||||
Let us quicly remind ourselves how we derive the above method.
|
||||
|
||||
Perhaps the most celebrated of all one-dimensional root-finding
|
||||
routines is Newton's method, also called the Newton-Raphson
|
||||
method. This method is distinguished from the previously discussed
|
||||
methods by the fact that it requires the evaluation of both the
|
||||
function $f$ and its derivative $f'$ at arbitrary points. In this
|
||||
sense, it is taylored to cases with e.g., transcendental equations.
|
||||
If you can only calculate the derivative
|
||||
numerically and/or your function is not of the smooth type, we
|
||||
discourage the use of this method.
|
||||
|
||||
!split
|
||||
===== The equations =====
|
||||
|
||||
The Newton-Raphson formula consists geometrically of extending the
|
||||
tangent line at a current point until it crosses zero, then setting
|
||||
the next guess to the abscissa of that zero-crossing. The mathematics
|
||||
behind this method is rather simple. Employing a Taylor expansion for
|
||||
$x$ sufficiently close to the solution $s$, we have
|
||||
|
||||
|
||||
!bt
|
||||
\[
|
||||
f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
|
||||
\label{eq:taylornr}
|
||||
\]
|
||||
!et
|
||||
|
||||
For small enough values of the function and for well-behaved
|
||||
functions, the terms beyond linear are unimportant, hence we obtain
|
||||
|
||||
|
||||
!bt
|
||||
\[
|
||||
f(x)+(s-x)f'(x)\approx 0,
|
||||
\]
|
||||
!et
|
||||
yielding
|
||||
!bt
|
||||
\[
|
||||
s\approx x-\frac{f(x)}{f'(x)}.
|
||||
\]
|
||||
!et
|
||||
|
||||
Having in mind an iterative procedure, it is natural to start iterating with
|
||||
!bt
|
||||
\[
|
||||
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
|
||||
\]
|
||||
!et
|
||||
|
||||
!split
|
||||
===== Simple geometric interpretation =====
|
||||
|
||||
The above is Newton-Raphson's method. It has a simple geometric
|
||||
interpretation, namely $x_{n+1}$ is the point where the tangent from
|
||||
$(x_n,f(x_n))$ crosses the $x-$axis. Close to the solution,
|
||||
Newton-Raphson converges fast to the desired result. However, if we
|
||||
are far from a root, where the higher-order terms in the series are
|
||||
important, the Newton-Raphson formula can give grossly inaccurate
|
||||
results. For instance, the initial guess for the root might be so far
|
||||
from the true root as to let the search interval include a local
|
||||
maximum or minimum of the function. If an iteration places a trial
|
||||
guess near such a local extremum, so that the first derivative nearly
|
||||
vanishes, then Newton-Raphson may fail totally
|
||||
|
||||
|
||||
!split
|
||||
===== Extending to more than one variable =====
|
||||
|
||||
Newton's method can be generalized to systems of several non-linear equations
|
||||
and variables. Consider the case with two equations
|
||||
!bt
|
||||
\[
|
||||
\begin{array}{cc} f_1(x_1,x_2) &=0\\
|
||||
f_2(x_1,x_2) &=0\end{array},
|
||||
\]
|
||||
!et
|
||||
which we Taylor expand to obtain
|
||||
|
||||
!bt
|
||||
\[
|
||||
\begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1
|
||||
\partial f_1/\partial x_1+h_2
|
||||
\partial f_1/\partial x_2+\dots\\
|
||||
0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1
|
||||
\partial f_2/\partial x_1+h_2
|
||||
\partial f_2/\partial x_2+\dots
|
||||
\end{array}.
|
||||
\]
|
||||
!et
|
||||
Defining the Jacobian matrix ${\bf \hat{J}}$ we have
|
||||
!bt
|
||||
\[
|
||||
{\bf \hat{J}}=\left( \begin{array}{cc}
|
||||
\partial f_1/\partial x_1 & \partial f_1/\partial x_2 \\
|
||||
\partial f_2/\partial x_1 &\partial f_2/\partial x_2
|
||||
\end{array} \right),
|
||||
\]
|
||||
!et
|
||||
we can rephrase Newton's method as
|
||||
!bt
|
||||
\[
|
||||
\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
|
||||
\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
|
||||
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
|
||||
\]
|
||||
!et
|
||||
where we have defined
|
||||
!bt
|
||||
\[
|
||||
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
|
||||
-{\bf \hat{J}}^{-1}
|
||||
\left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
|
||||
\]
|
||||
!et
|
||||
We need thus to compute the inverse of the Jacobian matrix and it
|
||||
is to understand that difficulties may
|
||||
arise in case ${\bf \hat{J}}$ is nearly singular.
|
||||
|
||||
It is rather straightforward to extend the above scheme to systems of
|
||||
more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Steepest descent =====
|
||||
|
||||
@@ -178,6 +375,118 @@ 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 mention the approach called the _standard Steepest descent_, which again leads to us having to be able to compute a matrix.
|
||||
|
||||
"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
|
||||
===== Conjugate 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. 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.
|
||||
|
||||
More details will be added here soon.
|
||||
|
||||
!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;
|
||||
x = ConjugateGradient(A,b,x0);
|
||||
xsd = SteepestDescent(A,b,x0);
|
||||
cout << "The approximate solution using Conjugate Gradient is: " << endl;
|
||||
x.Print();
|
||||
cout << endl;
|
||||
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;
|
||||
f = A*x-b;
|
||||
i = 0;
|
||||
while (i <= IterMax){
|
||||
z = A*f;
|
||||
c = dot(f,f);
|
||||
alpha = c/dot(f,z);
|
||||
x = x - alpha*f;
|
||||
f = A*x-b;
|
||||
if(sqrt(dot(f,f)) < tolerance) break;
|
||||
i++;
|
||||
}
|
||||
return x;
|
||||
}
|
||||
!ec
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Revisiting our first homework =====
|
||||
|
||||
@@ -839,3 +1148,4 @@ which gives
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user