Single variable functions

1. Find \(f'(x)\)
2. Set \(f'(x)=0\) and solve for \(x\). Call all \(x_0\) such that \(f'(x_0)=0\) critical values
3. Find \(f''(x)\). Evaluate at each \(x_0\)
   • If \(f''(x) > 0\), concave up, and therefore a local minimum
   • If \(f''(x) < 0\), concave down, and therefore a local maximum
   • If it’s the global maximum/minimum, it will produce the largest/smallest value for \(f(x)\)
   • On a closed range along the domain, check the endpoints as well

data_frame(x = c(-4, 4)) %>%
  ggplot(aes(x)) +
  stat_function(fun = function(x) -x^2, size = .5) +
  geom_vline(xintercept = c(-3, 3), linetype = 2) +
  labs(title = expression(f(x) == -x^2),
       x = expression(x),
       y = expression(f(x)))

data_frame(x = c(-5, 5)) %>%
  ggplot(aes(x)) +
  stat_function(fun = function(x) x^3, size = .5) +
  geom_vline(xintercept = c(-3, 3), linetype = 2) +
  labs(title = expression(f(x) == x^3),
       x = expression(x),
       y = expression(f(x)))

Multivariable functions

First derivative test: Gradient

Suppose \(f:X \rightarrow \Re^{n}\) with \(X \subset \Re^{1}\) is a differentiable function. Define the gradient vector of \(f\) at \(\boldsymbol{x}_{0}\), \(\nabla f(\boldsymbol{x}_{0})\) as

\[ \begin{eqnarray} \nabla f (\boldsymbol{x}_{0}) & = & \left(\frac{\partial f (\boldsymbol{x}_{0}) }{\partial x_{1} }, \frac{\partial f (\boldsymbol{x}_{0}) }{\partial x_{2} }, \frac{\partial f (\boldsymbol{x}_{0}) }{\partial x_{3} }, \ldots, \frac{\partial f (\boldsymbol{x}_{0}) }{\partial x_{n} } \right) \end{eqnarray} \]

It is the first partial derivatives for each variable \(x_n\) stored in a vector. So if \(\boldsymbol{a} \in X\) is a local extremum, then,

\[ \begin{eqnarray} \nabla f(\boldsymbol{a}) & = & \boldsymbol{0} \\ & = & (0, 0, \ldots, 0) \end{eqnarray} \]

That is, the root(s) of the gradient are where \(f(\boldsymbol{a})\) equals \(\boldsymbol{0}\) in \(n\)-dimensional space.

Second derivative test: Hessian

Suppose \(f:X \rightarrow \Re^{1}\) , \(X \subset \Re^{n}\), with \(f\) a twice differentiable function. We will define the Hessian matrix as the matrix of second derivatives at \(\boldsymbol{x}^{*} \in X\),

\[ \begin{eqnarray} \boldsymbol{H}(f)(\boldsymbol{x}^{*} ) & = & \begin{pmatrix} \frac{\partial^{2} f }{\partial x_{1} \partial x_{1} } (\boldsymbol{x}^{*} ) & \frac{\partial^{2} f }{\partial x_{1} \partial x_{2} } (\boldsymbol{x}^{*} ) & \ldots & \frac{\partial^{2} f }{\partial x_{1} \partial x_{n} } (\boldsymbol{x}^{*} ) \\ \frac{\partial^{2} f }{\partial x_{2} \partial x_{1} } (\boldsymbol{x}^{*} ) & \frac{\partial^{2} f }{\partial x_{2} \partial x_{2} } (\boldsymbol{x}^{*} ) & \ldots & \frac{\partial^{2} f }{\partial x_{2} \partial x_{n} } (\boldsymbol{x}^{*} ) \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} f }{\partial x_{n} \partial x_{1} } (\boldsymbol{x}^{*} ) & \frac{\partial^{2} f }{\partial x_{n} \partial x_{2} } (\boldsymbol{x}^{*} ) & \ldots & \frac{\partial^{2} f }{\partial x_{n} \partial x_{n} } (\boldsymbol{x}^{*} ) \\ \end{pmatrix} \nonumber \end{eqnarray} \]

• Requires differentiating on the entire gradient with respect to each \(x_n\)

Basic procedure summarized

1. Calculate gradient (first derivative)
2. Set equal to zero, solve system of equations
3. Calculate Hessian (second derivative)
4. Assess Hessian at critical values
5. Boundary values? (if relevant)

Example: Maximum likelihood estimation for a normal distribution

Suppose that we draw an independent and identically distributed random sample of \(n\) observations from a normal distribution,

\[ \begin{eqnarray} Y_{i} & \sim & \text{Normal}(\mu, \sigma^2) \\ \boldsymbol{Y} & = & (Y_{1}, Y_{2}, \ldots, Y_{n} ) \end{eqnarray} \]

The likelihood function defines the plausibility of a statistical model’s parameter value, given specific observed data. It is proportional to the following function:

\[ \begin{eqnarray} L(\mu, \sigma^2 | \boldsymbol{Y} ) & \propto & \prod_{i=1}^{n} f(Y_{i}|\mu, \sigma^2) \\ &\propto & \prod_{i=1}^{N} \frac{\exp[ - \frac{ (Y_{i} - \mu)^2 }{2\sigma^2} ]}{\sqrt{2 \pi \sigma^2}} \\ & \propto & \frac{\exp[ -\sum_{i=1}^{n} \frac{(Y_{i} - \mu)^2}{2\sigma^2} ]}{ (2\pi)^{n/2} \sigma^{2n/2} } \end{eqnarray} \]

Taking the logarithm, we have

\[ \begin{eqnarray} \log L(\mu, \sigma^2|\boldsymbol{Y} ) & = & -\sum_{i=1}^{n} \frac{(Y_{i} - \mu)^2}{2\sigma^2} - \frac{n}{2} log(2 \pi) - \frac{n}{2} \log (\sigma^2) \end{eqnarray} \]

In R, let’s draw 10,000 realizations from \(Y_{i} \sim \text{Normal}(0.25, 100)\).

set.seed(1234)

y <- rnorm(n = 10000, mean = 0.25, sd = 10)

From this observed data, we want to obtain the optimal values for \(\mu, \sigma^2\). To employ a grid search, we evaluate the function \(\log L(\mu, \sigma^2|\boldsymbol{Y})\) iteratively with different combinations of potential values for \(\mu\) and \(\sigma^2\). Each time we evaluate this likelihood function we store the resulting log-likelihood value. This optimization problem requires maximization (maximizing the log-likelihood), so whichever combination of parameter values produces the largest log-likelihood value is the solution.

# define the log-likelihood function for a normal distribution
log.like <- function(mu, sigma.2, y){
    part1 <- - (1 / (2 * sigma.2)) * sum((y - mu)^2) 
    part2 <- - (length(y)/2) * log(sigma.2)         
    out <- part1 + part2
    return(out)
    }

# define search space
(search_space <- expand.grid(mu = seq(-2, 2, by = .01),
            sigma2 = seq(8^2, 12^2, by = .1)) %>%
  as_tibble)

## # A tibble: 321,201 x 2
##       mu sigma2
##    <dbl>  <dbl>
##  1 -2        64
##  2 -1.99     64
##  3 -1.98     64
##  4 -1.97     64
##  5 -1.96     64
##  6 -1.95     64
##  7 -1.94     64
##  8 -1.93     64
##  9 -1.92     64
## 10 -1.91     64
## # ... with 321,191 more rows

# evaluate parameters
system.time((search_space <- search_space %>%
  mutate(logLik = map2_dbl(mu, sigma2, log.like, y = y))))

##    user  system elapsed 
##    11.9     3.1    15.4

# which parameter values maximize the likelihood function
filter(search_space, logLik == max(logLik))

## # A tibble: 1 x 3
##      mu sigma2  logLik
##   <dbl>  <dbl>   <dbl>
## 1  0.31   97.5 -27900.

# draw a heatmap of the log-likelihood values
ggplot(search_space, aes(mu, sigma2, fill = logLik)) +
  geom_raster() +
  geom_contour(aes(z = logLik)) +
  scale_fill_continuous(type = "viridis") +
  labs(x = expression(mu),
       y = expression(sigma^2))

## Warning: Removed 1 rows containing non-finite values (stat_contour).

Here, the grid search provides a reasonable approximation of \(\hat{\mu}, \hat{\sigma}^2\).

mean(y)

## [1] 0.311

var(y)

## [1] 97.5

However, the drawbacks are clear. This exhaustive search took over 18 seconds on a high-powered laptop for a relatively simple likelihood function and only two parameters. Grid searches suffer from the curse of dimensionality – as the number of parameters in the search space increases, the sample space increases rapidly. Grid searches can be parallelized to improve performance, but still require the researcher to define the domain of the sample space manually.

Description of algorithm

We want to find the root value \(x_1\) based on a starting point of \(x_0\)

\[ \begin{eqnarray} 0 & \cong & f(x_0) + \frac{f^{'}(x_0)}{1} (x_1 - x_0) \end{eqnarray} \]

where the quality of the approximation increases with better and better guesses of \(x_0\). With a little rearrangement:

\[ \begin{eqnarray} x_1 & \cong & x_0 - \frac{f(x_0)}{f'(x_0)} \end{eqnarray} \]

we calculate candidate values for \(x_1\). The procedure is algorithmic and iterative because our initial guess for \(x_1\) will not be optimal. However, we can use the same procedure multiple times substituting the new value for \(x_1\) into the function as \(f_0\) and updating \(x_1\). The generalized form is:

\[ x_{n+1} \cong x_n - \frac{f(x_n)}{f'(x_n)} \]

Repeat this step sufficiently until \(f(x_{j+1})\) is sufficiently close to zero, then stop.

Newton-Raphson is not a foolproof method in that it can fail to converge to a solution. In fact, there are no “perfect” computational methods that will always converge on a solution. Particularly, the assumption \(f''(x)\) exists and is continuous near the actual root \(r\) must be made.

Notice for optimization problems, we want to find the roots of the first derivative of the initial function \(f(x)\), not the roots of \(f(x)\) directly.

Implementation in R

The unitroot() 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:

\[y = x^3 + 2x - 5\] This function looks like this:

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

# draw a plot of the function
range <- data_frame(x = c(-5, 5))

f0 <- range %>%
  ggplot(aes(x)) +
  stat_function(fun = func2, size = .5, linetype = 2) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  labs(x = expression(x),
       y = expression(f(x)))
f0

It looks like the function equals 0 when \(x\) is approximately 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, interval = c(2, 3))

## $root
## [1] 2.09
## 
## $f.root
## [1] -0.000115
## 
## $iter
## [1] 5
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.1e-05

It took 5 iterations for Newton-Raphson to converge on this solution. 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)\).

# f - the function to optimize
# a - lower bound for the search
# b - upper bound for the search
# tol - tolerance (stopping criteria for the algorithm)
# n - maximum number of iterations to attempt. will not exceed even if
#     tolerance is not achieved

newton_raphson <- function(f, a, b, tol = 1e-5, n_iter = 1000) {
  require(numDeriv) # Package for computing f'(x)
  
  x0 <- a # Set start value to supplied lower bound
  k <- vector("numeric", length = n_iter) # 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_iter) {
    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 <- x1
      res <- list('root approximation' = root.approx, 'iterations' = k[1:i])
      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:

newton_raphson(f = func2, a = 2, b = 3)

## Loading required package: numDeriv

## $`root approximation`
## [1] 2.09
## 
## $iterations
## [1] 2.10 2.09 2.09 2.09

Implementation in R

The following code is an example of implementing gradient descent algorithm to find the minimum of the function \(f{(x)} = x^{4} - 3x^{3} + 2\).

# define the function
obj_func <- function(x) {
  x^4 - 3 * x^3 + 2
}

# draw a plot of the function
range <- data_frame(x = c(-1, 3.25))

f0 <- range %>%
  ggplot(aes(x)) +
  stat_function(fun = obj_func, size = .5, linetype = 2) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  labs(x = expression(x),
       y = expression(f(x)))
f0

Note that we are looking for \(f\)’s minimum by solving its derivative being equal to zero.

\[\nabla f(x) = 4x^3 - 9x^2 = 0\]

# define the gradient
gradient <- function(x) {
  (4 * x^3) - (9 * x^2)
}

f0 <- range %>%
  ggplot(aes(x)) +
  stat_function(fun = gradient, size = .5, linetype = 2) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  labs(x = expression(x),
       y = expression(f(x)))
f0

And the \(x\) can be updated with gradient descent method every iteration in the form of

\[x_{(k+1)} = x_{(k)} - \alpha \nabla f\bigl(x_{(k)}\bigr)\]

where \(k = 1, 2, \ldots,\) maximum iteration, and \(\alpha\) is the step size.

# func - the function to optimize
# grad - gradient of the function to optimize
# stepsize - size of each step
# tol - tolerance (stopping criteria for the algorithm)
# iter - maximum number of iterations to attempt. will not exceed even if
#        tolerance is not achieved

grad_desc <- function(func, grad, stepsize = 0.003, tol = 1e-5, iter = 500){
  # randomly initialize a value to x
  set.seed(100)
  x <- floor(runif(1)*10)
  
  # create a vector to contain all xs for all steps
  x.All <- vector("numeric", iter)
  
  # gradient descent method to find the minimum
  for(i in 1:iter){
    x1 <- x - stepsize * grad(x)
    x.All[[i]] <- x1
    
    # Once the difference between x0 and x1 becomes sufficiently small, output the results.
    if (abs(x - x1) < tol) {
      root.approx <- x
      res <- list('root approximation' = root.approx, 'iterations' = x.All[1:i])
      return(res)
    }
    
    # If gradient descent has not yet reached convergence set x as x0 and continue
    x <- x1
  }
  
  print('Too many iterations in method')
}

grad_desc(func = obj_func, grad = gradient)

## $`root approximation`
## [1] 2.25
## 
## $iterations
##   [1] 2.92 2.85 2.79 2.74 2.70 2.66 2.62 2.59 2.56 2.54 2.52 2.50 2.48 2.46
##  [15] 2.45 2.43 2.42 2.41 2.40 2.39 2.38 2.37 2.36 2.35 2.35 2.34 2.33 2.33
##  [29] 2.32 2.32 2.31 2.31 2.31 2.30 2.30 2.30 2.29 2.29 2.29 2.29 2.28 2.28
##  [43] 2.28 2.28 2.28 2.27 2.27 2.27 2.27 2.27 2.27 2.27 2.27 2.26 2.26 2.26
##  [57] 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26 2.26
##  [71] 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25
##  [85] 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25
##  [99] 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25
## [113] 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25 2.25
## [127] 2.25

Example: one-dimensional utility

Let’s try to figure out what policy a politician preferred given their utility function. Here’s the utility function:

politician.utility <- function(policy.content){
        politician.support <- -(policy.content - 1)^2 + 8
        return(politician.support)
}

And here is how this utility function looked like:

range <- data_frame(x = c(-2, 4))

range %>%
  ggplot(aes(x)) +
  stat_function(fun = politician.utility, size = .5, linetype = 2) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  labs(title = "Legislator's Utility Function",
       x = "Policy Ideology",
       y = "Legislator Utility From Policy")

Where is the peak of this function? Let’s use optimize() to find out.

First, we pass our function politician.utility() to optimize(). When using optimizers you will need to write a function that returns an output you would like to maximize or minimize. Second, we need to provide an interval over which optimize() can search for a maximum value. Third, note that we need to specify that we want to find the maximum of this function. The default setting in optimize() is to find the minimum of a function.

optimize(f = politician.utility, interval = c(-2, 4), maximum = TRUE)

## $maximum
## [1] 1
## 
## $objective
## [1] 8

• maximum tells us value for \(x\) at which the function is maximized
• objective tells us the value of the function \(f(x)\) evaluated at this point

Example: two-dimensional utility

When politicians decide whether to support a bill there are often more than a single dimension at play. For example, an issue might touch on both economic and social issues. Suppose a politician has different utility functions over two aspects of a bill. What combination of economic and social policy content might they prefer from the bill?

Here is a new function describing politican support for a policy proposal:

politician.utility.2d <- function(params){
  # split parameters
  economic.content <- params[1]
  social.content <- params[2]
  
  # calculate utility
  politician.support <- (-(economic.content - 1)^2 + 8 ) + (-(social.content + 2)^2 + 8)
  return(politician.support)
}

# generate data values
expand.grid(economic.substance = seq(from=-8,to=8, by=.2),
            social.substance = seq(from=-8, to=8, by=.2)) %>%
  as_tibble %>%
  mutate(utility = map2_dbl(economic.substance, social.substance,
                            ~ politician.utility.2d(c(.x, .y)))) %>%
  ggplot(aes(economic.substance, social.substance, fill = utility)) +
  geom_raster() +
  geom_contour(aes(z = utility)) +
  scale_fill_continuous(type = "viridis") +
  labs(title = "Utility Based on Policy Content on Two Dimensions",
       x = "Economic Policy Content",
       y = "Social Policy Content")

Note the function receives only one argument (params) which will contain two elements. The function then assigns the first element of params to the first dimension of the legislator’s utility function, and the second one to the social dimension. This is because optim() will optimize only the first argument of our function – regardless of how many elements (i.e., variables) it contains.

Now we can run optim(). The function takes at least three arguments. The first argument is par, which specifies the starting points at which optim() will start to assess optima. The second argument is fun – the function we want to optimize. As the function optimize() minimizes by default, in order to find the maximum of our function we need to add the argument “control=list(fnscale=-1)”. An alternative would be to have our politician.utility.2d() function return the negative of politician support. This would accomplish the same thing without needing to change around the optim parameters.

Let’s optimize our function now:

optim(par = c(-1, 0), fn = politician.utility.2d, control = list(fnscale = -1))

## $par
## [1]  1 -2
## 
## $value
## [1] 16
## 
## $counts
## function gradient 
##       61       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

optim() has several outputs. The most important ones are the two elements contained in the par output. They show the optimized values for each parameter.