Skip to content

Instantly share code, notes, and snippets.

@rwalk
Created March 8, 2015 23:22
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save rwalk/8f3a38064da2876acb8a to your computer and use it in GitHub Desktop.
Save rwalk/8f3a38064da2876acb8a to your computer and use it in GitHub Desktop.
Quadratic programming: Use matrix inputs to solve a QP with ipoptr
ipoptr_qp <- function(Dmat, dvec, Amat, bvec, ub=100){
# Solve the quadratic program
#
# min -d^T x + 1/2 x^T D x
# s.t. A%*%x>= b
#
# with ipoptr.
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)
# initialize with the global unconstrained minimum, as done in quadprog package
# NOTE: This will only work if lb <= x0 <= ub. If this is not the case,
# use x0 = lb can be used instead.
x0 <- solve(Dmat, dvec)
# call the solver
sol <- 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(sol)
}
@vate01
Copy link

vate01 commented Dec 30, 2016

Hello, thanks for this code. I have been trying to reproduce the solve.QP example from quadprog R library using ipoptr_qp, but I have not been success. These are the lines I have tested:

>source('quadprog_ipoptr_translate.R')
>Dmat       <- matrix(0,3,3)
>diag(Dmat) <- 1
>dvec       <- c(0,5,0)
>Amat       <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
>bvec       <- c(-8,2,0)
# solve.QP(Dmat,dvec,Amat,bvec=bvec)
> ipoptr_qp(Dmat,dvec,Amat,bvec=bvec)
Error in do.call(c, sapply(1:length(h_mask), function(i) Dmat[i, h_mask[[i]]])) : 
  second argument must be a list

Do you have any hint to solve this?

Thanks in advance!

@adam-bec
Copy link

adam-bec commented Oct 24, 2017

Hi, Thanks for the implementation, but I have an error running it:

> error in ipoptr_qp(Dmat, dvec, t(Amat), bvec) : Function "make.sparse" not found

any idea?
Thanks

@nfultz
Copy link

nfultz commented Mar 16, 2020

@vate01 - be careful with Amat, it's transposed in solve.QP but since yours is square it's producing different numbers.

@adam-bec you need to load the ipoptr package first.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment