Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Quadratic Programming: Compare three methods for solving a randomly generated QP in 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.