Quadratic programming demo: Circus Tent

bigTent.R

#  This gist uses rgl and quadprog to draw a circus tent
#  by 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(rgl)
library(quadprog)

# below are only need if you want to run profiling
library(kernlab)  
library(reshape2)
library(ggplot2)
######################################################
# 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)
}

# 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()
}

######################################################
# 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)
}

plotEnsemble<-function(tentbase){
  x <- tentbase$x
  y <- tentbase$y
  z <- tentbase$z
  open3d()
  rgl.surface(x,y, z, color='blue',alpha=1, smooth=TRUE)
  rgl.surface(x,y, matrix(1/2,length(x), length(y)), color='red', alpha=.4)
  rgl.bg(color="white")
}

plotTentSolution <- function(tentbase, tent){
  x <- tentbase$x
  y <- tentbase$y
  z <- tentbase$z
  
  open3d()
  rgl.surface(x, y , z, color='blue',alpha=1, smooth=TRUE)
  rgl.surface(x, y , tent, color='red', alpha=.4)
  rgl.bg(color="white")
}

######################################################
# QP solver staging
######################################################
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)
}

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)
}

#######################################################
# Solver time profiling
#######################################################

profile <- function(Nbase, itvec){
  ipop_times <- c()
  quadprog_times <- c()
  errors <- 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
    qmat <-NULL
    
    # 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
    ipmat <- NULL
    errors <- c(errors, sqrt(sum((sol1$solution-primal(sol2))^2)))
  }
  
  times <- data.frame(problem.size=dims, 
                   quadprog.time=quadprog_times,
                   ipop.time = ipop_times,
                   difference = errors)  
  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 examples
######################################################
# INSTRUCTIONS: 
# Set nFold = 1 to see original MATLAB demo
# Set nFold >=2 to build more complex tents using symmetries
# 
# Caution: running with nBase >= 2**5 and nFold>=2 will need at least 750 MB
# of available memory.
tentbase <-makeTentEnsemble(nBase=2**5, nFold=1)
plotEnsemble(tentbase)
qpmat <-quadprogStage(tentbase)
spy(qpmat$Dmat)
sol  <- solve.QP(qpmat$Dmat, qpmat$dvec, qpmat$Amat, qpmat$bvec)
tent <- matrix(sol$solution, nrow(tentbase$z), ncol(tentbase$z), byrow=FALSE)  # reshape before passing to plotter
plotTentSolution(tentbase, tent)

# PROFILING [can be memory intensive]
# time the solvers for the QPs
times <- profile(2**3, 1:4)
profilePlot(times)