Note-Optimization

Use of this document

This is a study note for using \(stats\), \(lpSolve\), \(quadprog\), \(optimx\) package for Optimization. For more details on the study material see:

Prerequisites

# Optimization in R: https://www.is.uni-freiburg.de/resources/computational-economics/5_OptimizationR.pdf

library(stats) # General purpose
library(lpSolve) # Linear Programming
library(quadprog) # Quadratic Programming
library(optimx) # Non-Linear Programming

options(digits=6)
library(dplyr)

1. Introduciton

Optimizationusesarigorousmathematicalmodeltodeterminethe most efficient solution to a described problem.

Common R packages for optimization

Problem type	Package::function
General purpose (1-dim.)	stats::optimize()
General purpose (n-dim.)	stats::optim()
Linear Programming	lpSolve::lp()
Quadratic Programming	quadprog::solve.QP()
Non-Linear Programming	optimx::optimx()
General interface	ROI::ROI_solve()

1.1 Common Method

• Linear Programming
• Quadratic Programming
• Non-Linear Programming

1.2 Elements

• Objective function
• Constraints

2. General purpose

stats::optim(par, fun, ...) is for n-dimensional general purpose optimization.

• Argument par sets the initial values of the search
• fn specifies a function to optimize over
• Optional, named argument method chooses an algorithm

# Define objective function and initial value
objfun <- function(x) 2*(x[1]-1)^2 + 5*(x[2]-3)^2 + 10
init_x <- c(1, 1)
# Call optimization routine
op_result <- stats::optim(par = init_x, fn = objfun)
# optimizatized Objective 
op_result$par

## [1] 1.00017 3.00023

# Optimal input values
op_result$value

## [1] 10

3. Linear Programming

Linear Programming in Dim(x) = 2

Linear Programming Mathematical specification

lpSolve::lp(direction="min", objective.in, const.mat, const.dir, const.rhs) is for n-dimensional Linear Programming optimization.

e.g. maximize the total profit

# objective function  
const.mat <- matrix(c(20, 12, 1/15, 1/15), nrow=2, byrow=TRUE)
# initial value
objective.in <- c(25, 20)
# constraints
const.dir <- c("<=", "<=")
const.rhs <- c(1800, 8)
# Call optimization routine
op_result <- lpSolve::lp(direction="max",  objective.in = objective.in, 
                         const.mat = const.mat, const.rhs = const.rhs,
                         const.dir = const.dir)
# Optimal input values
op_result$solution

## [1] 45 75

# optimizatized Objective 
op_result$objval

## [1] 2625

4. Quadratic Programming

Quadratic Programming in Dim(x) = 1

Quadratic Programming Mathematical specification

quadprog::solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE) is for n-dimensional Quadratic Programming optimization.

e.g.,portfolio optimization, circus tent problem, demand response

# objective function  
Dmat <- matrix(0,3,3); diag(Dmat) <- 1
dvec <- c(0,5,0)
# initial value
## N.A.
# constraints
Amat <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
bvec <- c(-8,2,0)
# Call optimization routine
op_result <- quadprog::solve.QP(Dmat,dvec,Amat,bvec=bvec)
# Optimal input values
op_result$solution

## [1] 0.47619 1.04762 2.09524

# optimizatized Objective 
op_result$value

## [1] -2.38095

5. Non-Linear Programming

Overview: Non-Linear Optimization

optimx::optimx(par, fn, gr=Null, Hess=Null, lower=inf, upper=inf, method='', itnmax=Null, ...) provides collection of for Non-Linear Multi-Dimensional Programming methods.

5.1 Method

• Gradient based: Gradient descent methods (method = "CG")
• Hessian based (method = "BFGS", method = "L-BFGS-B"): Newton and quasi-Newton methods
• Non-gradient based (method = "Nelder-Mead" as default): Golden section search, Nelder-Mead, etc. ### 5.1.1 Golden Section Search

Basic steps of Golden Section Search

5.1.2 Nelder-Mead

Basic steps of Nelder-Mead

5.1.3 Conjugate Gradients

Basic steps of Conjugate Gradients

5.1.4 Newton-Raphson

Basic steps of Newton-Raphson

5.1.5 BFGS and L-BFGS-B

Broyden-Fletcher-Goldfarb-Shanno(BFGS) algorithm builds on the idea of Newton’s method to take gradient information into account. Gradient information comes from an approximation of the Hessian matrix. No guaranteed conversion; expecially problematic if Taylor expansion does not fit well.

L-BFGS-Bstandsforlimited-memory-BFGS-box:

• Extension of BFGS
• Memory efficient implementation
• Additionally handles box constraints

5.1 1-dimensional Non-Linear optimization.

stats::optimize(f, interval, tol, maximum = FALSE) is for 1-dimensional Non-Linear optimization using Golden Section Search.

# ifferentiable function
f <- function(x)(x - 1/3)^2 
op_result <- stats::optimize(f, interval = c(0, 1), tol = 0.0001)
# Non-differentiable function
f <- function(x) return(abs(x-2) + 2*abs(x-1))
op_result <- optimize(f, interval = c(0, 3), tol = 0.0001)
# Optimal input values
op_result$minimum

## [1] 1.00001

# optimizatized Objective 
op_result$objective

## [1] 1.00001

5.2 n-dimensional Non-Linear optimization.

# objective function: Himmelblau’s function
fn <- function(para){ # Vector of the parameters 
  matrix.A <- matrix(para, ncol=2)
  x <- matrix.A[,1]
  y <- matrix.A[,2]
  f.x <- (x^2+y-11)^2+(x+y^2-7)^2
  return(f.x) 
}

xy <- as.matrix(expand.grid(seq(-5,5,length = 101), seq(-5,5,length = 101)))
colnames(xy) <- c("x", "y")
df <- data.frame(fnxy = fn(xy), xy)
library(lattice)
wireframe(fnxy ~ x*y, data = df, shade = TRUE, drape=FALSE,
          scales = list(arrows = FALSE), screen = list(z=-240, x=-70, y=0))

# initial value
par <- c(1,1)
# Call optimization routine: compare search methods
optimx(par, fn, control=list(all.methods=TRUE, save.failures=TRUE, trace=0)) %>%
  `rownames<-`(rownames(.)) %>% 
  mutate(method = rownames(.)) %>% 
  select(method, everything()) %>% 
  mutate_if(is.numeric, signif, digits = 6) %>% 
  DT::datatable(.)

# Call optimization routine: apply specific method
optimx(par, fn, method = c("CG")) %>% mutate_if(is.numeric, signif, digits = 6) %>% DT::datatable(.)