Lab5: General Optimization

1 General Optimization

In mathematics, nonlinear programming (NLP) is the process of solving an optimization problem where some of the constraints or the objective function are nonlinear. An optimization problem is one of calculation of the extrema (maxima, minima or stationary points) of an objective function over a set of unknown real variables and conditional to the satisfaction of a system of equalities and inequalities, collectively termed constraints. It is the sub-field of mathematical optimization that deals with problems that are not linear.

1.1 Two-dimensional function

The blue region (shown in the diagram below) is the feasible region. The tangency of the line with the feasible region represents the solution. The line is the best achievable contour line (area with a given value of the objective function).

2-dimensional example

This simple problem can be defined by the constraints:

\[ \begin{eqnarray} x_1 & \ge & 0 \\ x_2 & \ge & 0 \\ x_1^2 + x_2^2 & \ge & 1 \\ x_1^2 + x_2^2 & \le & 2 \end{eqnarray} \] with an objective function to be maximized

\[ f(x) = x_1 + x_2 , \text{where } x = (x_1, x_2). \]

1.2 Three-dimensional function

The tangency (see diagram bellow) of the top surface with the constrained space in the center represents the solution.

3-dimensional example

This simple problem can be defined by the constraints:

\[ \begin{eqnarray} x_1^2 − x_2^2 + x_3^2 \le 2 \\ x_1^2 + x_2^2 + x_3^2 \le 10 \\ \end{eqnarray} \]

with an objective function to be maximized

\[f(x) = x_1x_2 + x_2x_3, \text{where } x = (x_1, x_2, x_3).\]

2 Steepest Descent (SD)

The method of steepest descent works on functions which have a single derivative. It is used most often in problems involving more than 1 variable. The essential idea of steepest descent is that the function decreases most quickly in the direction of the negative gradient. Let’s assume we have the following function:

\[f(x)=f(x_1,x_2,\cdots, x_n)\]

The objective is to find the maximum or minimum value (according to our purpose).

2.1 SD Algorithm

• The method starts at an initial guess \(x\).
• The next guess is made by moving in the direction of the negative gradient. The location of the minimum along this line can then be found by using a one-dimensional search algorithm such as golden section search.
• The nth update is then

\[x_n=x_{n−1}−\alpha f'(x_{n−1})\]

where \(\alpha\) is chosen to minimize the one-dimensional function:

\[g(\alpha)=f(x_{n−1}−\alpha f'(x_{n−1}))\]

In order to use golden section, we need to assume that \(\alpha\) is in an interval. So in this case, we take the interval to be \([0,h]\) where \(h\) is a value that we need to choose.

2.2 SD Algorithm in R

steepestdescent <- function(f, fprime, start, h,
 tol=1e-7, maxiter=100) {
 x <- start
 g <- function(alpha) { f(x - alpha*fpx) }
 niter <- 0
 while(niter < maxiter & sum(abs(fprime(x))) > tol) {
 fpx <- fprime(x)
 alpha <- golden(g, 0, h)
 x <- x - alpha*fpx
 niter <- niter + 1
 }
 if (niter == maxiter) {
 print("Warning: Maximum number of iterations reached")
 }
 c("Minimizer" = x)
 }

2.3 SD Algorithm in Visualization

You can see the visualization process of Steepest Descent in the following graphic, or click here to see another example of Steepest descent method, specifically for a quadratic function. more

3-dimensional example

2.4 Two-dimensional Example

A simple example of a function of 2 variables to be minimized is

\[f(x)=f(x_1,x_2)={(2−x_1)^2\over 2x_2^2} + {(3−x_1)^2 \over 2x_2^2}+log(x_2)\]

Note that \(x_2\) should be positive, so we might need to protect against negative values of \(x_2\).

• First : We need functions for both the function and the gradient

f1 <- function(x) {
 (2-x[1])^2/(2*x[2]^2) +(3-x[1])^2/(2*x[2]^2) + log(x[2])
 }
 f1prime <- function(x) {
 c(-(2-x[1])/x[2]^2 - (3-x[1])/x[2]^2,-(2-x[1])^2/x[2]^3 -(3-x[1])^2/x[2]^3 + 1/x[2])
 }

• Second: The location of the minimum along this line can then be found by using a one-dimensional search algorithm such as golden section search.

golden <- function (f, a, b, tol = 0.0000001)
 {
 ratio <- 2 / (sqrt(5) + 1)
 x1 <- b - ratio * (b - a)
 x2 <- a + ratio * (b - a)
 f1 <- f(x1)
 f2 <- f(x2)

 while(abs(b - a) > tol) {
 if (f2 > f1) {
 b <- x2
 x2 <- x1
 f2 <- f1
 x1 <- b - ratio * (b - a)
 f1 <- f(x1)
 } else {
 a <- x1
 x1 <- x2
 f1 <- f2
 x2 <- a + ratio * (b - a)
 f2 <- f(x2)
 }
 }
 return((a + b) / 2)
 }

• Third: Let’s try a starting value of x=(.1,.1).

steepestdescent(f1, f1prime, start=c(.1,.1), h=.1)

## [1] "Warning: Maximum number of iterations reached"

## Minimizer1 Minimizer2 
##  2.4992967  0.7123675

We haven’t converged yet.

• Fourth: One possibility is to run the procedure again, using the most recent result as our starting guess.

steepestdescent(f1, f1prime,start=c(2.4992967, 0.7123675), h=.1)

## Minimizer1 Minimizer2 
##  2.5000000  0.7071068

steepestdescent(f1, f1prime, start=c(.1, .1), h=.1,maxiter=200)

## Minimizer1 Minimizer2 
##  2.5000000  0.7071068

Done. The value has been converged.

2.5 Multivariate Normal

One can use steepest descent to compute the maximum likelihood estimate of the mean in a multivariate Normal density, given a sample of data. However, when the data are highly correlated, as they are in the simulated example below, the log-likelihood surface can be come difficult to optimize. In such cases, a very narrow ridge develops in the log-likelihood that can be difficult for the steepest descent algorithm to navigate.

In the example below, we actually compute the negative log-likelihood because the algorithm is designed to minimize functions.

set.seed(2020-09-30)
mu <- c(1, 2)
S <- rbind(c(1, .9), c(.9, 1))
x <- MASS::mvrnorm(500, mu, S)
nloglike <- function(mu1, mu2) {
        dmv <- mvtnorm::dmvnorm(x, c(mu1, mu2), S, log = TRUE)
        -sum(dmv)
}
nloglike <- Vectorize(nloglike, c("mu1", "mu2"))
nx <- 40
ny <- 40
xg <- seq(-5, 5, len = nx)
yg <- seq(-5, 6, len = ny)
g <- expand.grid(xg, yg)
nLL <- nloglike(g[, 1], g[, 2])
z <- matrix(nLL, nx, ny)
par(mar = c(4.5, 4.5, 1, 1))
contour(xg, yg, z, nlevels = 40, xlab = expression(mu[1]), 
        ylab = expression(mu[2]))
abline(h = 0, v = 0, lty = 2)

Note that in the figure above the surface is highly stretched and that the minimum (1,2) lies in the middle of a narrow valley. For the steepest descent algorithm we will start at the point (−5,−2) and track the path of the algorithm.

library(dplyr, warn.conflicts = FALSE)
norm <- function(x) x / sqrt(sum(x^2))
Sinv <- solve(S)                          ## I know I said not to do this!
step1 <- function(mu, alpha = 1) {
        D <- sweep(x, 2, mu, "-")
        score <- colSums(D) %>% norm
        mu + alpha * drop(Sinv %*% score)
}
steep <- function(mu, n = 10, ...) {
        results <- vector("list", length = n)
        for(i in seq_len(n)) {
                results[[i]] <- step1(mu, ...)
                mu <- results[[i]]
        }
        results
}
m <- do.call("rbind", steep(c(-5, -2), 8))
m <- rbind(c(-5, -2), m)

par(mar = c(4.5, 4.5, 1, 1))
contour(xg, yg, z, nlevels = 40, xlab = expression(mu[1]), 
        ylab = expression(mu[2]))
abline(h = 0, v = 0, lty = 2)
points(m, pch = 20, type = "b")

We can see that the path of the algorthm is rather winding as it traverses the narrow valley. Now, we have fixed the step-length in this case, which is probably not optimal. However, one can still see that the algorithm has some difficulty navigating the surface because the direction of steepest descent does not take one directly towards the minimum ever.

3 The Newton Direction

Newton’s method is an iterative method for finding the roots of a differentiable function F, which are solutions to the equation \(F(x) = 0.\) In optimization, Newton’s method is applied to the derivative \(f\) ′ of a twice-differentiable function \(f\) to find the roots of the derivative (solutions to \(f ′(x) = 0\), also known as the stationary points of \(f\). These solutions may be minima, maxima, or saddle points. The central problem of optimization is minimization of functions. Let us first consider the case of univariate functions, i.e., functions of a single real variable. We will later consider the more general and more practically useful multivariate case.

Given a twice differentiable function \({\displaystyle f:\mathbb {R} \to \mathbb {R} }{\displaystyle f:\mathbb {R} \to \mathbb {R} },\) we seek to solve the optimization problem

\({\displaystyle \min _{x\in \mathbb {R} }f(x).}{\displaystyle \min _{x\in \mathbb {R} }f(x).}\)

Newton’s method attempts to solve this problem by constructing a sequence \({\displaystyle \{x_{k}\}}\) from an initial guess (starting point) \({\displaystyle x_{0}\in \mathbb {R} }\) that converges towards a minimizer \(x_*\) of \(f\) by using a sequence of second-order Taylor approximations of \({\displaystyle f}\) around the iterates. The second-order Taylor expansion of f around \(x_k\) is

\[{\displaystyle f(x_{k}+t)\approx f(x_{k})+f'(x_{k})t+{\frac {1}{2}}f''(x_{k})t^{2}.}{\displaystyle f(x_{k}+t)\approx f(x_{k})+f'(x_{k})t+{\frac {1}{2}}f''(x_{k})t^{2}.}\] The next iterate \({\displaystyle x_{k+1}}\) is defined so as to minimize this quadratic approximation in \({\displaystyle t}\), and setting \({\displaystyle x_{k+1}=x_{k}+t}\). If the second derivative is positive, the quadratic approximation is a convex function of \({\displaystyle t}\), and its minimum can be found by setting the derivative to zero. Since

\[{\displaystyle \displaystyle 0={\frac {\rm {d}}{{\rm {d}}t}}\left(f(x_{k})+f'(x_{k})t+{\frac {1}{2}}f''(x_{k})t^{2}\right)=f'(x_{k})+f''(x_{k})t,}\] the minimum is achieved for

\[{\displaystyle t=-{\frac {f'(x_{k})}{f''(x_{k})}}.}\] Putting everything together, Newton’s method performs the iteration

\[{\displaystyle x_{k+1}=x_{k}+t=x_{k}-{\frac {f'(x_{k})}{f''(x_{k})}}.}\]

3.1 Newton in Visualization

A comparison of gradient descent (green) and Newton’s method (red) for minimizing a function (with small step sizes). Newton’s method uses curvature information to take a more direct route.

knitr::include_graphics("newton.png")

3.2 Geometric interpretation

The geometric interpretation of Newton’s method is that at each iteration, it amounts to the fitting of a paraboloid to the surface of \({\displaystyle f(x)}\) at the trial value \(x_{k}\), having the same slopes and curvature as the surface at that point, and then proceeding to the maximum or minimum of that paraboloid (in higher dimensions, this may also be a saddle point).

3.3 Higher dimensions

The above iterative scheme can be generalized to \(d\) dimensions by replacing the derivative with the gradient (different authors use different notation for the gradient, including \({\displaystyle f'(x)=\nabla f(x)=g_{f}(x)\in \mathbb {R} ^{d}}\), and the reciprocal of the second derivative with the inverse of the Hessian matrix (different authors use different notation for the Hessian, including \({\displaystyle f''(x)=\nabla ^{2}f(x)=H_{f}(x)\in \mathbb {R} ^{d\times d}}\). One thus obtains the iterative scheme \({\displaystyle x_{k+1}=x_{k}-[f''(x_{k})]^{-1}f'(x_{k}),k\geq 0.}\) Often Newton’s method is modified to include a small step size \({\displaystyle 0<\gamma \leq 1}\) instead of \({\displaystyle \gamma =1}\) \[{\displaystyle x_{k+1}=x_{k}-\gamma [f''(x_{k})]^{-1}f'(x_{k}).}\] This is often done to ensure that the Wolfe conditions are satisfied at each step of the method. For step sizes other than 1, the method is often referred to as the relaxed or damped Newton’s method.

3.4 Convergence

If \(f\) is a strongly convex function with Lipschitz Hessian, then provided that \(x_{0}\) is close enough to \({\displaystyle x_{*}=\arg \min f(x)}\), the sequence \({\displaystyle x_{0},x_{1},x_{2},\dots }\) generated by Newton’s method will converge to the (necessarily unique) minimizer \(x_{*}\) of \(f\) quadratically fast That is, \[{\displaystyle \|x_{k+1}-x_{*}\|\leq {\frac {1}{2}}\|x_{k}-x_{*}\|^{2},\qquad \forall k\geq 0.}\]

3.5 Computing the Newton direction

Finding the inverse of the Hessian in high dimensions to compute the Newton direction \(h=-(f''(x_{k}))^{-1}f'(x_{k})\) can be an expensive operation. In such cases, instead of directly inverting the Hessian, it’s better to calculate the vector \(h\) as the solution to the system of linear equations \[{\displaystyle [f''(x_{k})]h=-f'(x_{k})}\] which may be solved by various factorizations or approximately (but to great accuracy) using iterative methods. Many of these methods are only applicable to certain types of equations, for example the Cholesky factorization and conjugate gradient will only work if \(f''(x_{k})\) is a positive definite matrix. While this may seem like a limitation, it’s often a useful indicator of something gone wrong; for example if a minimization problem is being approached and \(f''(x_{k})\) is not positive definite, then the iterations are converging to a saddle point and not a minimum. On the other hand, if a constrained optimization is done (for example, with Lagrange multipliers), the problem may become one of saddle point finding, in which case the Hessian will be symmetric indefinite and the solution of \(x_{k+1}\) will need to be done with a method that will work for such, such as the \(LDL\) variant of Cholesky factorization or the conjugate residual method. There also exist various quasi-Newton methods, where an approximation for the Hessian (or its inverse directly) is built up from changes in the gradient. If the Hessian is close to a non-invertible matrix, the inverted Hessian can be numerically unstable and the solution may diverge. In this case, certain workarounds have been tried in the past, which have varied success with certain problems. One can, for example, modify the Hessian by adding a correction matrix \(B_k\) so as to make \(f''(x_{k})\)+\(B_{k}\) positive definite. One approach is to diagonalize the Hessian and choose \(B_{k}\) so that \(f''(x_{k})\)+\(B_k\) has the same eigenvectors as the Hessian, but with each negative eigenvalue replaced by \({\epsilon >0}\) An approach exploited in the Levenberg–Marquardt algorithm (which uses an approximate Hessian) is to add a scaled identity matrix to the Hessian, \({\mu I}\), with the scale adjusted at every iteration as needed. For large \({\mu }\) and small Hessian, the iterations will behave like gradient descent with step size \({1/\mu }\) . This results in slower but more reliable convergence where the Hessian doesn’t provide useful information.

3.6 Examples of Newton

Before proceeding to an implementation of the Newton-Raphson method in R, it is worth working through some examples to get an understanding of the definitions and equations above. The NR method can be used to approximate square roots such as \(√10.\) The square root of 10 is about three, so we can use that as a good starting value. It often helps to plot the function to see where the roots occur. The function is first rearranged to be an expression of \(f(x)\) before plotting. \[x=√10\] \[x−√10=0\] \[f(x)=x2−10\] Plotting the function yields the graph below.

func <- function(x) {x^2 - 10}
curve(func, col = 'red', lwd = 2, from = 0, n = 100, xlim=c(0, 5), ylab='f(x)')
abline(a=0, b=0, lty = 5)

It can be seen from the graph the function crosses the x-axis at \(√10\) on an interval [3, 4]. The derivative of the function is \(2x\). And so on until the xn estimates are within a particular level of tolerance. This example converges in three iterations. Thus 3.162278 is the estimated root of the function. This result can be verified by simply taking the square root of 10.

sqrt(10)

## [1] 3.162278

3.7 Newton Method in R

The \(uniroot\) function in R provides an implementation of Newton-Raphson for finding the root of an equation. The function is only capable of finding one root in the given interval. The \(rootSolve\) package features the uniroot.all function which extends the uniroot routine to detect multiple roots should they exist.

The equation for which we wish to find the root is: \[f(x)=x^3−2x−5\] First, write a function to express the equation above.

func2 <- function(x) {x^3 - 2* x - 5}

Plot the function to visualize how the equation behaves and where any roots may be located.

curve(func2, xlim=c(-5,5), col='blue', lwd=2, lty=2, ylab='f(x)')
abline(h=0)
abline(v=0)

It looks like the function equals 0 when y is about 2. To find the root of the equation, use the uniroot function with a starting value of 2 and upper bound of 3.

uniroot(func2, c(2,3))

## $root
## [1] 2.094541
## 
## $f.root
## [1] -0.0001147009
## 
## $iter
## [1] 5
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

As suspected, the root of the function is very close to 2 at 2.0945412. The iterations can be written as follows:

\[x_{n+1}=x_n-{\frac {f(x_{n})}{f'(x_{n})}}=x_n-{\frac {f(x^3_{n}-2x-5)}{3x^2_n-2}}\] \[x_1={\frac {(2)^3-2(2)-5}{3(2)^2-2}}=1.281239\] \[x_2={\frac {(1.281239)^3-2(1.281239)-5}{3(1.281239)^2-2}}=3.147821\] \[x_3={\frac {(3.147821)^3-2(3.147821)-5}{3(3.147821)^2-2}}=2.430257\] \[x_4={\frac {(2.430257)^3-2(2.430257)-5}{3(2.430257)^2-2}}=2.144418\] \[x_5={\frac {(2.144418)^3-2(2.144418)-5}{3(2.144418)^2-2}}=2.095897\] \[x_5={\frac {(2.095897)^3-2(2.095897)-5}{3(2.095897)^2-2}}=2.094552\] The manual calculations equal the output of the function to four decimals. The minuscule difference between results is likely the cause of roundoff and significance errors in the manual calculation. A quick and dirty function to perform the method in R can be implemented to further verify our understanding of the Newton-Raphson method. The numDeriv package is used to compute the derivative \(f′(x)\), though this could also be done by taking the limit with a sufficiently small \(h\).

newton.raphson <- function(f, a, b, tol = 1e-5, n = 1000) {
  require(numDeriv) # Package for computing f'(x)
  
  x0 <- a # Set start value to supplied lower bound
  k <- n # Initialize for iteration results
  
  # Check the upper and lower bounds to see if approximations result in 0
  fa <- f(a)
  if (fa == 0.0) {
    return(a)
  }
  
  fb <- f(b)
  if (fb == 0.0) {
    return(b)
  }

  for (i in 1:n) {
    dx <- genD(func = f, x = x0)$D[1] # First-order derivative f'(x0)
    x1 <- x0 - (f(x0) / dx) # Calculate next value x1
    k[i] <- x1 # Store x1
    # Once the difference between x0 and x1 becomes sufficiently small, output the results.
    if (abs(x1 - x0) < tol) {
      root.approx <- tail(k, n=1)
      res <- list('root approximation' = root.approx, 'iterations' = k)
      return(res)
    }
    # If Newton-Raphson has not yet reached convergence set x1 as x0 and continue
    x0 <- x1
  }
  print('Too many iterations in method')
}

Computing the root of the above equation with the newton.raphson function yields:

library (numDeriv)
newton.raphson(func2, 2, 3)

## $`root approximation`
## [1] 2.094551
## 
## $iterations
## [1] 2.100000 2.094568 2.094551 2.094551

The custom Newton-Raphson function is slightly off from the uniroot result, however; they are equal up to 4 decimals and is adequate for our purposes.

#Summary The Newton-Raphson method is a powerful and relatively straightforward method for finding the roots of an equation. It has many advantages but suffers from several drawbacks such as a readily calculated derivative, inconsistencies when \(f′(x)=0\) and so on. As explored in the post, it often makes sense first to plot the function in question to visualize how it behaves before attempting to locate the \(root(s)\) as bad starting values can be detrimental to results. Other numerical methods for estimating roots of equations such as the Bisection, Secant and Brent’s methods will be examined in future posts.