Table of Contents

Introduction: Optimization and non linear methods¶

This notebook explores introductory level algorithms and tools that can be used for non linear optimization. It starts with root finding algorithms in 1 dimensions using a simple example and then moves onto optimization methods (minimum finding) and multidimensional cases. For each example, a graph of the function and a visualization of the steps taken by the optimization algorithm are shown.

Of interest are the worked examples done with "pen and paper" as well as the direct implementations using python and numpy only. A study of the optimization of the well known Rosenbrock function is done using both the Gradient Descent and Newton Rapshon approaches.

Sources¶

The following is a collection of sources of interest which were used to build this notebook:

• Duke university: Optimization and Non-linear Methods
• Sebastien Ruder's blog
• Nonlinear Optimization course at the University of Jyväskylä
• Numerical python tours
• Louis Tiao's blog: visualizing-and-animating-optimization-algorithms-with-matplotlib

Libraries¶

Finding roots of non linear functions (1 - dimension)¶

Example used throughout this section:¶

Consider the following polynomial function which is continuous, differentiable and well behaved. It is an interesting example to use when testing root finding algorithms because it exhibits several basins of attractions in 1D. A basin of attraction is a set of initial values leading to the same solution (i.e. root in this case), see here for some details. In this example, different starting points will lead to different solutions, or in some cases no solutions at all when the algorithm "blows up"

The next section will investigate different methods of finding the roots of:

$f(x) = x^3 - 2x^2 - 11x + 12$

$f'(x) = 3x^2 - 4x - 11$

Bisection Method¶

The bisection method is one of the simplest methods for finding zeroes of a non-linear function. It is guaranteed to find a root - but it can be slow. The main idea comes from the intermediate value theorem: If $f(a)$ and $f(b)$ have different signs and $f$ is continous, then f must have a zero between a and b. We evaluate the function at the midpoint, $c=\frac{1}{2}(a+b).$ $f(c)$ is either zero, has the same sign as f(a) or the same sign as $f(b)$. Suppose $f(c)$ has the same sign as f(a) (as pictured below). We then repeat the process on the interval $[c,b]$.

Newton Raphson's method¶

Newton's method, also known as Newton-Raphson's method, is a very famous and widely used method for solving nonlinear algebraic equations. Compared to the other methods we will consider, it is generally the fastest one (usually by far). It does not guarantee that an existing solution will be found, however.

Newton's method is normally formulated with an iteration index n:

$x_{n+1} = x_n - \large \frac{f(x_n)}{f'(x_n)}$

To avoid storing intermediate values in memory, one can overwrite the previous value and work with a single variable:

$x := x - \frac{f(x)}{f'(x)}$

In this example:

$f(x) = x^3 - 2x^2 - 11x + 12$

$f'(x) = 3x^2 - 4x - 11$

Newton Raphson's using Scipy¶

Similar results are obtained using Scipy's built in optimizer

Secant method¶

When finding the derivative $f′(x)$ in Newton's method is problematic, or when function evaluations take too long; we may adjust the method slightly. Instead of using tangent lines to the graph we may use secants. The approach is referred to as the secant method.

The idea of the secant method is to think as in Newton's method, but instead of using $f′(x_n)$, we approximate this derivative by a finite difference or the secant, i.e., the slope of the straight line that goes through the two most recent approximations $x_n$ and $x_n{−1}$. This slope reads

$f'(x_n) \approx \large \frac{f(x_n) - f(x_{n-1})}{x_n - x_{n-1}}$

Inserting this expression for $f'(x_n)$ into Newton's method simply gives us the secant method:

$x_{n+1} = x_n - f(x_n) \times \large \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$

As with the Bisection Search - the secant method requires two starting points

Finding extrema of a function¶

In this section, we are concerned in optimizing a function rather than finding a root, or solution. As a result the algorithms will try to find local / global minimum and maximum points. A well known function used to compare optimization algorithms is the Rosenbrock function, described next.

Introducing the Rosenbrock function¶

In mathematical optimization, the Rosenbrock function is a non-convex function, introduced by Howard H. Rosenbrock in 1960, which is used as a performance test problem for optimization algorithms. It is also known as Rosenbrock's valley or Rosenbrock's banana function.

The global minimum is inside a long, narrow, parabolic shaped flat valley. To find the valley is trivial. To converge to the global minimum, however, is difficult.

The function is defined by:

$f(x,y) = (a - x)^2 + b(y - x^2)^2$ and has a global minimum at $(x,y) = (a,a^2)$. Usually we choose parameters $(a,b) = (1, 100)$ giving:

$f(x,y) = (1 - x)^2 + 100(y - x^2)^2$ with minimum at $(1,1)$

Gradient

$\nabla f = [-400xy + 400x^3 + 2x -2, 200y - 200x^2]$

Hessian

$H = \begin{bmatrix}-400y + 1200x^2 + 2 & -400x \\ -400x & 200 \end{bmatrix} $

Gradient descent method¶

Gradient descent (or steepest descent) is a first-order iterative optimization algorithm for finding the minimum of a function. To find a local minimum of a function using gradient descent, one takes steps proportional to the negative of the gradient (or approximate gradient) of the function at the current point. If instead one takes steps proportional to the positive of the gradient, one approaches a local maximum of that function; the procedure is then known as gradient ascent.

The approach is based on the observation that if a multivariable function is defined and differentiable in a neighborhood of a point, then the function decreases fastest in the direction of the negative gradient. It follows that if:

$\mathbf{x_{n+1} = x_n} - \gamma \nabla \mathbf{F(x)}$

For $\gamma$ small enough, then $F(x_n) \geq F(x_{n+1})$

The value of $\gamma$ is called the stepsize and is allowed to change at every iteration. Some more advanced implementations find the optimal value of $gamma$ at each step by using line search methods...

Comments

• The gradient is everywhere perpendicular to the contour lines.
• After each line minimization the new gradient is always orthogonal to the previous step direction (true of any line minimization.)
• Consequently, the iterates tend to zig-zag down the valley in a very inefficient manner

Gradient Descent implementation¶

Plotting the Rosenbrock function and the Gradient Descent iteration steps¶

Comments¶

As expected, the algorithm finds the "valley" - however it struggles to reach the global minimum and gets stuck in a zigzag behaviour. Even after several thousand iterations, the minimum is not found. Trying different initial values leads to vastly different results, many of which are not close to the global minimum, and in many cases leading to infinite values (i.e. algorithm does not converge)

The behaviour is also very dependent on the initial choice of the $\gamma$ value.

Gradient descent on a simpler function (quadratic)¶

Consider a simpler function to optimize.

$f(x,y) = .01x^2 + .1y^2$

As this function is convex and significantly easier to optimize, we would expect the Gradient Descent algorithm to give better results than for the Rosenbrock function.

Improving the Gradient descent with line search (to be completed)¶

We can improve the algorithm by replacing the constant $\gamma$ value with a value that changes at each step.

To be completed

Newton's method¶

Similarly to the root finding algorithm, Newton's method for finding extrema takes the form:

$x_{n+1} = x_n - \frac{f'(x_n)}{f''(x_n)}$

Extending to multivariate case, we replace the derivative with the gradient $\nabla f(x_n)$ and the reciprocal of the second derivative with the inverse of the Hessian matrix:

$x_{n+1} = x_n - [\mathbf{H} \ f(x_n)]^{-1} \nabla f(x_n)$

Finding the inverse of the Hessian in high dimensions can be computationally expensive. In such cases, instead of inverting the Hessian it is better to calculate the vector $\delta x = x_{n+1} - x_n$ as the solution to the linear system of equations:

$[\mathbf{H} \ f(x_n)] \nabla \mathbf{x_n} = - \nabla f(x_n)$

Note that since Newton's method converges quadratically, any (smooth) quadratic function will be optimized in one step, as is shown in the trivial example below.

Example 1) Simple case

$f(x_1, x_2) = 0.5 x_1^2 + 2.5x_2^2$ where $\mathbf{x} = \begin{bmatrix}x_1 \\ x_2 \end{bmatrix}$

$\nabla \ f(x_1,x_2) = \begin{bmatrix}x_1 \\ 5x_2 \end{bmatrix} $ and $H = J(\nabla^T) = \begin{bmatrix} 1 & 0 \\ 0 & 5 \end{bmatrix} $

$H^{-1} = \begin{bmatrix}1 & 0 \\ 0 & 1/5 \end{bmatrix}$

$\mathbf{x_{n+1} = x_n - H^{-1}_n} \nabla f(\mathbf{x_n})$

$\begin{bmatrix}x_1\\ x_2 \end{bmatrix}_{n+1} := \begin{bmatrix}x_1\\ x_2 \end{bmatrix}_{n} - \begin{bmatrix}1 & 0 \\ 0 & 1/5 \end{bmatrix} \begin{bmatrix}x_1 \\ 5x_2 \end{bmatrix}_n = \begin{bmatrix}x_1 \\ x_2 \end{bmatrix}_n - \begin{bmatrix}x_1 \\ x_2 \end{bmatrix}_n = \begin{bmatrix}0 \\ 0 \end{bmatrix}$

Plotting the iterations on top of the surface and contour plot¶

Comments¶

The Newton Raphson method converges extremely rapidly to the global solution, despite the significant challenge posed by the Rosenbrock function. Irrespective of the starting point (within the scope of the graph above) the minimum is found within 4 - 6 iterations. Compare this to the several thousand iterations required by the gradient descent approach.

Newton Raphson optimization using Scipy¶

Comparing Newton and Gradient Descent in presence of a single saddle point¶

$f_2(x,y) = .01x^2 - .1y^2$

$f'_2(x,y) = .02x - .2y$

Comparing Newton and Gradient Descent in presence of a multiple saddle points¶

Using the himmerlblau function¶

$f(x,y)=(x^{2}+y-11)^{2}+(x+y^{2}-7)^{2}.\quad $

$\nabla f(x,y) = [2 (-7 + x + y^2 + 2 x (-11 + x^2 + y)), 2 (-11 + x^2 + y + 2 y (-7 + x + y^2))] $

$H f(x,y) = \begin{bmatrix} 4 (x^2 + y - 11) + 8 x^2 + 2 & 4 x + 4 y \\ 4 x + 4 y & 4 (x + y^2 - 7) + 8 y^2 + 2) \end{bmatrix}$

Drawbacks of Newton's method¶

A combination of two reasons:

• Newton method attracts to saddle points;
• saddle points are common in machine learning, or in fact any multivariable optimization.

Look at the function $$f=x^2-y^2$$

If you apply multivariate Newton method, you get the following. $$\mathbf{x}_{n+1} = \mathbf{x}_n - [\mathbf{H}f(\mathbf{x}_n)]^{-1} \nabla f(\mathbf{x}_n)$$

Let's get the Hessian: $$\mathbf{H}= \begin{bmatrix} \dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1\,\partial x_n} \\[2.2ex] \dfrac{\partial^2 f}{\partial x_2\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2\,\partial x_n} \\[2.2ex] \vdots & \vdots & \ddots & \vdots \\[2.2ex] \dfrac{\partial^2 f}{\partial x_n\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_n\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2} \end{bmatrix}.$$

$$\mathbf{H}= \begin{bmatrix} 2 & 0 \\[2.2ex] 0 & -2 \end{bmatrix}$$

Invert it: $$[\mathbf{H} f]^{-1}= \begin{bmatrix} 1/2 & 0 \\[2.2ex] 0 & -1/2 \end{bmatrix}$$

Get the gradient: $$\nabla f=\begin{bmatrix} 2x \\[2.2ex] -2y \end{bmatrix}$$

Get the final equation: $$\mathbf{\begin{bmatrix} x \\[2.2ex] y \end{bmatrix}}_{n+1} = \begin{bmatrix} x \\[2.2ex] y \end{bmatrix}_n -\begin{bmatrix} 1/2 & 0 \\[2.2ex] 0 & -1/2 \end{bmatrix} \begin{bmatrix} 2x_n \\[2.2ex] -2y_n \end{bmatrix}= \mathbf{\begin{bmatrix} x \\[2.2ex] y \end{bmatrix}}_n - \begin{bmatrix} x \[2.2ex] y

\end{bmatrix}_n¶

\begin{bmatrix} 0 \\[2.2ex] 0 \end{bmatrix} $$

So, you see how the Newton method led you to the saddle point at $x=0,y=0$.

In contrast, the gradient descent method will not lead to the saddle point. The gradient is zero at the saddle point, but a tiny step out would pull the optimization away as you can see from the gradient above - its gradient on y-variable is negative.

Source: https://stats.stackexchange.com/questions/253632/why-is-newtons-method-not-widely-used-in-machine-learning/253655

A comparison between gradient descent and Newton's method¶

You asked two questions: Why don't more people use Newton's method, and why do so many people use stochastic gradient descent? These questions have different answers, because there are many algorithms that lessen the computational burden of Newton's method but often work better than SGD.

First: Newton's Method takes a long time per iteration and is memory-intensive. As jwimberley points out, Newton's Method requires computing the second derivative, $H$, which is $O(N^2)$, where $N$ is the number of features, while computing the gradient, $g$, is only $O(N)$. But the next step is $H^{-1} g$, which is $O(N^3)$ to compute. So while computing the Hessian is expensive, inverting it or solving least squares is often even worse. (If you have sparse features, the asymptotics look better, but other methods also perform better, so sparsity doesn't make Newton relatively more appealing.)

Second, many methods, not just gradient descent, are used more often than Newton; they are often knockoffs of Newton's method, in the sense that they approximate a Newton step at a lower computational cost per step but take more iterations to converge. Some examples:

• Because of the expense of inverting the Hessian, ``quasi-Newton" methods like BFGS approximate the inverse Hessian, $H^{-1}$, by looking at how the gradient has changed over the last few steps.

• BFGS is still very memory-intensive in high-dimensional settings because it requires storing the entire $O(N^2)$ approximate inverse Hessian. Limited memory BFGS (L-BFGS) calculates the next step direction as the approximate inverse Hessian times the gradient, but it only requires storing the last several gradient updates; it doesn't explicitly store the approximate inverse Hessian.

• When you don't want to deal with approximating second derivatives at all, gradient descent is appealing because it only uses only first-order information. Gradient descent is implicitly approximating the inverse Hessian as the learning rate times the identity matrix. I, personally, rarely use gradient descent: L-BFGS is just as easy to implement, since it only requires specifying the objective function and gradient; it has a better inverse Hessian approximation than gradient descent; and because gradient descent requires tuning the learning rate.

• Sometimes you have a very large number of observations (data points), but you could learn almost as well from a smaller number of observations. When that is the case, you can use "batch methods", like stochastic gradient descent, that cycle through using subsets of the observations.

More advanced discussion¶

PRO

Indeed the ratio of the number of saddle points to local minima increases exponentially with the dimensionality N.

While gradient descent dynamics are repelled away from a saddle point to lower error by following directions of negative curvature, ...the Newton method does not treat saddle points appropriately; as argued below, saddle-points instead become attractive under the Newton dynamics.

https://arxiv.org/pdf/1406.2572.pdf

CONS

More people should be using Newton's method in machine learning*. I say this as someone with a background in numerical optimization, who has dabbled in machine learning over the past couple of years.

The drawbacks in answers here (and even in the literature) are not an issue if you use Newton's method correctly. Moreover, the drawbacks that do matter also slow down gradient descent the same amount or more, but through less obvious mechanisms.

• Using linesearch with the Wolfe conditions or using or trust regions prevents convergence to saddle points. A proper gradient descent implementation should be doing this too. The paper referenced in Cam.Davidson.Pilon's answer points out problems with "Newton's method" in the presence of saddle points, but the fix they advocate is also a Newton method.

• Using Newton's method does not require constructing the whole (dense) Hessian; you can apply the inverse of the Hessian to a vector with iterative methods that only use matrix-vector products (e.g., Krylov methods like conjugate gradient). See, for example, the CG-Steihaug trust region method.

• You can compute Hessian matrix-vector products efficiently by solving two higher order adjoint equations of the same form as the adjoint equation that is already used to compute the gradient (e.g., the work of two backpropagation steps in neural network training).

• Ill conditioning slows the convergence of iterative linear solvers, but it also slows gradient descent equally or worse. Using Newton's method instead of gradient descent shifts the difficulty from the nonlinear optimization stage (where not much can be done to improve the situation) to the linear algebra stage (where we can attack it with the entire arsenal of numerical linear algebra preconditioning techniques).

• Also, the computation shifts from "many many cheap steps" to "a few costly steps", opening up more opportunities for parallelism at the sub-step (linear algebra) level.

For background information about these concepts, I recommend the book "Numerical Optimization" by Nocedal and Wright.

*Of course, Newton's method will not help you with L1 or other similar compressed sensing/sparsity promoting penalty functions, since they lack the required smoothness.

Appendix - code snipets and to do¶

Rosenbrock function returning vector result¶

Finding non linear solutions - Multivariate¶

Example:

$f_2: (x,y) = x^4 + 2x^2y^2 + y^4$ where $x,y \in \mathbb{R}$

$f_3: (x,y) = x^2 + y^2$ where $x,y \in \mathbb{R}$

Defining and plotting the functions¶

Newton method - Multivariate¶

Generalizing the 1D approach by:

$x_{n+1} = x_n - J^{-1} \ f(x_n)$

• $x$ and $f(x)$ are now vectors
• $J^{-1}$ is the inverse Jacobian matrix.
• In general the Jacobian is not square and cannot be inverted, so we use the generalized inverse: $(J^TJ)^{-1}J^T$

$x_{n+1} = x_n - (J^TJ)^{-1}J^T \ f(x_n)$

Example for f3

$f_3: (x,y) = x^2 + y^2$ where $x,y \in \mathbb{R}$

$f'_3: (x,y) = 2x + 2y$

$ \frac{d \ f_3}{dx} = 2x $ and $ \frac{d \ f_3}{dy} = 2y $

$J = \left( 2x , 2y \right)$