Created
March 8, 2015 18:29
-
-
Save rwalk/666bde369571568b2fb9 to your computer and use it in GitHub Desktop.
Quadratic Programming: Compare three methods for solving a randomly generated QP in R
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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