Quadratic Programming: Compare three methods for solving a randomly generated QP in R

qp_experiment_random.R

#  This gist compares three methods for solving a randomly generated quadratic program in R.
#  author: R. Walker  (r_walker@zoho.com)
#  LICENSE: MIT
library(quadprog)
library(kernlab)
library(ipoptr)
library(ggplot2)
library(reshape2)
#########################################################
# Random QP Generation (in the style of quadprog)
#########################################################
randomQP <- function(N, d){
  # Generate matrices for a random QP
  # 
  # This method returns Dmat an SPD matrix and dvec a 
  # postive random vector defining the quadratic form
  #   -d^T x + 1/2 x^T D x
  # We also return Amat and bvec which give the constraint
  # relation
  #    Ax >= b
  # bvec is a ranod positive random vector of length N/d and Amat is
  # a binary matrix of dimension (N/d) x N.  Amat is defined
  # so that (Ax)_i is the sum of the i through (i+d) components of x.
  #
  # NOTE: differs with quadprog input format by using A rather than A^T
  # in the constraint.
  if(N%%d) stop("N must be divisible by d")
  
  # build an SPD matrix
  G <- matrix(rnorm(N^2, 0, 1), N, N)
  M <- lower.tri(matrix(1, N, N), diag=FALSE)
  G[!M] <- 0
  G <- t(G)%*%G 
  G <- G/sqrt(sum(G^2)) + diag(1,N,N)
  
  # QP data
  qpdata <- NULL
  qpdata$Dmat <- G
  qpdata$dvec <- 1+matrix(runif(N), N, 1)   # dvec with values near 0 seems to degrade stability in ipop?
  qpdata$Amat <- rbind(kronecker(diag(N/d), matrix(1, 1, d)))
  qpdata$bvec <- rbind(matrix(runif(N), N/d, 1))
  return(qpdata)
}

# Plot a sparse matrix -- similar to MATLAB spy
spy <- function(A){
  M <- ifelse(A==0,0,1)
  image(t(M)[ncol(M):1,], col=c("white","blue"), axes=FALSE, pch = 3, lty=0)
  title("Sparse Matrix Plot")
  box()
}

######################################################
# QP solver staging
######################################################

# quadprog
quadprogStage <- function(Dmat, dvec, Amat, bvec){
  f <- function(){
    # add posiitivy constraint
    N <- length(dvec)
    Amat <- rbind(Amat, diag(1,N,N))
    bvec <- rbind(bvec, matrix(0,N,1))
    
    # note, quadprog constraint set requires a transpose
    sol <- solve.QP(Dmat, dvec, t(Amat), bvec)
    return(sol)
  }
  return(f)
}

# ipop
ipopStage <- function(Dmat, dvec, Amat, bvec){ 
  # 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
  ub <- 1e8
  n <- length(dvec)
  cc <- -dvec     # note the sign difference relative to quadprog
  H <- Dmat
  b <- bvec
  r <- matrix(ub, nrow=length(bvec), ncol=1)
  A <- Amat  
  
  # box constraint over x
  l <- matrix(0, nrow=n, ncol=1)
  u <- matrix(ub, nrow=n, ncol=1)
  
  f <- function(){
    sol <- ipop(cc, H, A, b, l, u, r, margin=1e-3)
    return(sol)
  }
  return(f)
}

# ipoptr
ipoptrStage <- function(Dmat, dvec, Amat, bvec, ub=1e10){
  # 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)
  N <- length(dvec)
  
  # 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(Amat%*%x)
  eval_jac_g <- function(x) return(j_vals)
  eval_jac_g_structure <- j_mask
  constraint_lb <- bvec
  constraint_ub <- rep( ub, n)
  lb <- rep(0, 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,
                  lb = lb,
                  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)
}

#######################################################
# Experiment setup and outputs
#######################################################
profilePlot <- function(times){
  df <- melt(times[,1: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","yellow")) +
    theme(legend.position="bottom", legend.title=element_blank())
}

trial <- function(N){
  qpdata <- randomQP(N, 4)
  ipop_solve <- ipopStage(qpdata$Dmat, qpdata$dvec, qpdata$Amat, qpdata$bvec)
  quadprog_solve <- quadprogStage(qpdata$Dmat, qpdata$dvec, qpdata$Amat, qpdata$bvec)
  ipoptr_solve <- ipoptrStage(qpdata$Dmat, qpdata$dvec, qpdata$Amat, qpdata$bvec)
  quadprog.time <- system.time(sol1 <- quadprog_solve())[3]
  ipop.time <- system.time(sol2 <- ipop_solve())[3]
  ipoptr.time <- system.time(sol3 <- ipoptr_solve())[3]
  ipop.diff <- sqrt(sum(sol1$solution - primal(sol2))^2)
  ipoptr.diff <- sqrt(sum(sol1$solution - sol3$solution)^2)
  experiment <- c(N, quadprog.time, ipop.time, ipoptr.time, ipop.diff, ipoptr.diff)
  names(experiment) <- c("problem.size", "quadprog.time", "ipop.time", "ipoptr.time", "ipop.diff", "ipoptr.diff")
  return(experiment)
}

######################################################
# Run the experiment
######################################################
problemSize <- seq(from = 500, to = 7500, by = 1000)
experiment <- as.data.frame(t(sapply(problemSize, trial)))
profilePlot(experiment)