Compare 3 QP solvers on a demo problem in R

qp_solver_compare.R

#  This gist compares three methods for solving a quadratic program in R.
#  The original problem is from the MathWorks MATLAB demo at:
#  http://www.mathworks.com/help/optim/examples/large-scale-bound-constrained-quadratic-programming.html
#
#  More information on this problem at:
#  http://quantitate.blogspot.com/2014/04/the-circus-tent-problem-with-rs-quadprog.html
#
#  author: R. Walker  (r_walker@zoho.com)
#  LICENSE: MIT
library(quadprog)
library(kernlab)
library(ipoptr)

######################################################
# Helper functions
######################################################
Laplacian2D <- function(lx,ly, nx,ny){
  hx <- lx/(nx-1)
  hy <- ly/(ny-1)
  tr_x <- c(2,-1,rep(0,nx-2))
  tr_y <- c(2,-1,rep(0,ny-2))
  Tx <- toeplitz(tr_x)/(hx^2)
  Ty <- toeplitz(tr_y)/(hy^2)
  Ix <- diag(nx)
  Iy <- diag(ny)
  L <- kronecker(Tx,Iy) + kronecker(Ix,Ty)
  return(L)
}

######################################################
# Tent base builder and plotters
######################################################
# make a tent pole ensamble; recommeneded to use power of 2 for N to maintain symmetry
makeTentEnsemble <- function(nBase=2**5, nFold=1){
  n <- floor(nBase/2)       # number of rows and number of columns per quandrant
  m <- floor(n/2)       # midpoint row and column index in the first quadrant
  q2 = matrix(0,n,n)
  q2[m:(m+1), m:(m+1)] = .3
  q2[(n-1):n, (n-1):n] = 1/2
  q1 <- q2[,n:1]
  top  <- cbind(q2,q1)
  z <- rbind(top,top[n:1,])
  
  # use nFold to stack copies of z in an nFold by nFold grid
  # rep and stack columns then rep and stack the result down the row
  if(nFold>1) z<-apply(apply(z,2, rep, nFold), 1, rep, nFold)  
  
  # xy axis builder
  dd <- dim(z)
  x <- (1:dd[1])/dd[1]
  y <- (1:dd[2])/dd[2]
  
  # pack result in a list
  res = list(x=x, y=y, z=z)
  return(res)
}

######################################################
# QP solver staging
######################################################
# quadprog
quadprogStage <- function(tentbase){ 
  # System matrices
  # Note: quadprog's quadratic form is -d^T x + 1/2 x^T D x with t(A)%*%x>= b
  theta <- 1  # material coefficient on first order term in energy function; may want to adjust
  x <- tentbase$x
  y <- tentbase$y
  z <- tentbase$z
  Ny  <- nrow(z)
  Nx  <- ncol(z)
  N   <- Nx*Ny
  hx  <- 1/(Nx-1)
  hy  <- 1/(Ny-1)
  res <- list()
  res$Dmat <- hx*hy*Laplacian2D(1,1,Nx, Ny)                               
  res$dvec <- -theta*(hx*hy)*rep(1,N)       
  res$Amat <- t(diag(N))
  res$bvec <- matrix(z,N,1,byrow=FALSE)     # lower bound
  return(res)
}

# ipop
ipopStage <- function(tentbase){ 
  # solve the QP with kernlab solver
  # ipop solves the quadratic programming problem :
  # \min(c'*x + 1/2 * x' * H * x)
  # subject to: 
  # b <= A * x <= b + r
  # l <= x <= u
  n <- length(tentbase$z)
  res <- list()
  qpmat <- quadprogStage(tentbase)
  res$c <- -qpmat$dvec     # note the sign difference relative to quadprog
  res$H <- qpmat$Dmat
  res$b <- qpmat$bvec
  res$r <-  matrix(1, nrow=n, ncol=1)
  
  # box constraint over x  (note redundant with constraint on Ax)
  res$A <- diag(1,nrow=n)
  res$l <- matrix(0, nrow=length(tentbase$z), ncol=1)
  res$u <- matrix(1, nrow=length(tentbase$z), ncol=1)
  return(res)
}

#
# ipoptr
#

ipoptr_qp_stage <- function(Dmat, dvec, Amat, bvec, ub=100){
  # Stage the quadratic program
  #
  # min -d^T x + 1/2 x^T D x
  #   s.t. t(A)%*%x>= b
  #
  # for ipoptr given the matrix inputs.
  
  n <- length(bvec)
  # Jacobian structural components
  j_vals <- unlist(Amat[Amat!=0])
  j_mask <- make.sparse(ifelse(Amat!=0, 1, 0))
  
  # Hessian structural components
  h_mask <- make.sparse(ifelse(upper.tri(Dmat, diag=TRUE) & Dmat!=0, 1, 0))
  h_vals <- do.call(c, sapply(1:length(h_mask), function(i) Dmat[i,h_mask[[i]]]))
  
  # build the ipoptr inputs
  eval_f <- function(x) return(-t(dvec)%*%x + 0.5*t(x)%*%Dmat%*%x)
  eval_grad_f <- function(x) return(-dvec + Dmat%*%x)
  eval_h <- function(x, obj_factor, hessian_lambda) return(obj_factor*h_vals)
  eval_h_structure <- h_mask
  eval_g <- function(x) return(t(Amat)%*%x)
  eval_jac_g <- function(x) return(j_vals)
  eval_jac_g_structure <- j_mask
  constraint_lb <- bvec
  constraint_ub <- rep( ub, n)
  
  # wrapper function that will solves the ipoptr problem when called
  f <- function() {
    # initialize with the global unconstrained minimum.
    x0 <- solve(Dmat, dvec)
    
    # call the solver
    res <- ipoptr(x0 = x0, 
                  eval_f = eval_f,
                  eval_grad_f = eval_grad_f, 
                  eval_g = eval_g, 
                  eval_jac_g = eval_jac_g,
                  eval_jac_g_structure = eval_jac_g_structure,
                  constraint_lb = constraint_lb,
                  constraint_ub = constraint_ub,
                  eval_h = eval_h,
                  eval_h_structure = eval_h_structure)
    return(res)
  }
  
  return(f)
}

ipoptr_qp_tent_build <- function(tentbase){
  qpmat <- quadprogStage(tentbase)
  solveme <- ipoptr_qp_stage(qpmat$Dmat, qpmat$dvec, qpmat$Amat, qpmat$bvec)
  return(solveme)  
}

#######################################################
# Solver time profiling
#######################################################
profile <- function(Nbase, itvec){
  ipop_times <- c()
  quadprog_times <- c()
  ipoptr_times <- c()
  errors <- c()
  errors2 <- c()   
  dims <- c()
  for(i in itvec){
    tentbase <-makeTentEnsemble(Nbase, i)
    qpmat <-quadprogStage(tentbase)
    
    # time quadprog solver
    N<-dim(tentbase$z)[1]*dim(tentbase$z)[1]
    print(sprintf("Solving QP with N=%d using quadprog", N))
    time <- system.time(sol1 <- solve.QP(qpmat$Dmat, qpmat$dvec, qpmat$Amat, qpmat$bvec))
    dims <- c(dims, length(sol1$solution))
    quadprog_times <- c(quadprog_times, time[3])  # elapsed time
    
    # time ipop solver
    ipmat <-ipopStage(tentbase)
    print(sprintf("Solving QP with N=%d using ipop", N))
    time <- system.time(sol2 <- ipop(ipmat$c, ipmat$H, ipmat$A, ipmat$b, ipmat$l, ipmat$u, ipmat$r))
    ipop_times <- c(ipop_times, time[3])  # elapsed time
    errors <- c(errors, sqrt(sum((sol1$solution-primal(sol2))^2)))
    
    # time ipoptr solver
    solveme <-ipoptr_qp_tent_build(tentbase)
    print(sprintf("Solving QP with N=%d using ipoptr", N))
    time <- system.time(sol3 <- solveme())
    ipoptr_times <- c(ipoptr_times, time[3])  # elapsed time
    errors2 <- c(errors2, sqrt(sum((sol1$solution-sol3$solution)^2)))
  }
  
  times <- data.frame(problem.size=dims, 
                      quadprog.time=quadprog_times,
                      ipop.time = ipop_times,
                      ipoptr.time = ipoptr_times,
                      difference = errors,
                      difference2 = errors2)  
  return(times)
}

profilePlot <- function(times){
  df <- melt(times[,-4], id = "problem.size", value.name = "time", variable.name="method")
  ggplot(df, aes(x=problem.size, y=time, color=method)) + 
    ggtitle("Solve time vs problem size for tent problem")+
    geom_line(size=1.5) + geom_point(size=5) + scale_colour_manual(values=c("red","blue")) +
    theme(legend.position="bottom", legend.title=element_blank())
}

######################################################
# Run the experiment
######################################################
# Caution: memory and CPU intensive!
times <- profile(2**3, 1:10)