R function solve.QP error "constraints are inconsistent, no solution!"

Question

I am trying to run an optimzation by using the solve.QP function (from the quadprog package) with the following parameters

R = matrix( c( 2.231113e-05,-4.816095e-05,-5.115287e-05, 0,2.989584e-05,4.212173e-06,0,0, 5.504990e-05), ncol=3, byrow=T)
b = c(-1,0,rep(0,ncol(R)))
C = cbind(rep(1,ncol(R)), diag(ncol(R)))
C = cbind(-rep(1,ncol(R)),C)
d = as.matrix(c(57621264,78057622,171342351),ncol=1)
H = solve.QP(Dmat = R, factorized = TRUE, dvec = d, Amat = C, bvec = b)

But the error I receive is

Error in solve.QP(Dmat = R, factorized = TRUE, dvec = d, Amat = C, bvec = b) : 
constraints are inconsistent, no solution!

However, when I am using a different matrix for R

R2 = matrix( c( 0.05365071,-0.06364421,-0.04102565, 0, 0.08423283,-0.04048879,0,0,0.09659707), ncol=3, byrow=T)

the solve.QP call

H = solve.QP(Dmat = R2, factorized = TRUE, dvec = d, Amat = C, bvec = b)

does not cause any problems. My question is now why in the former case issues are arising.

Any help is greatly appreciated!

Answer 1

The elements in your matrix dvec are quite huge. Your matrix R is the inverse of upper triangular matrix in the cholesky factorization of Dmat, i.e. Dmat = M^T M where M = R^{-1}. So in unfactorized form:

M <- solve(R) 
Dmat <- t(M)%*%M

Dmat unfactorized is huge and roughly on the same scale as dvec:

#  [,1]       [,2]       [,3]
#[1,] 2008893283 3236243201 1619059085
#[2,] 3236243201 6332319710 2522625866
#[3,] 1619059085 2522625866 1641403882

Thus your issue is likely due to some kind of overflow error. To work around this, you can scale Dmat (unfactorized) and dvec:

sc <- norm(Dmat,"2")
solve.QP(Dmat = Dmat/sc, dvec=d/sc, Amat=C, bvec=b, meq=0, factorized=FALSE )

The solution is

# $solution
# [1] -1.220832e-17  0.000000e+00  1.043877e-01

which matches nicely with the solution the solution found by constrOptim in another answer.

Scaling Dmat and dvec is fine to do because the constrained minimizer of -d^T b + 1/2 b^T D b is the same as the constrained minimizer of sc*(-d^T b + 1/2 b^T D b) for any constant sc.

EDIT: To solve the problem in factorized form you could try something like the following scaling:

nn2 = sqrt(norm(d,"2"))
H = solve.QP(Dmat = R*nn2, dvec = d/(nn2^2), Amat = C, bvec = b, factorized=TRUE)
#$solution
#[1] 0.0000000 0.0000000 0.1043877

Answer 2

Maybe two problems. First you haven't input the parameter meq which tells the solver how many of your constraint equations are equality constraints vs. how many are inequality so it defaults to meq=0 which makes them all equality constraints so you've overdetermined your solution. Looking at your problem, I might guess that at least the last three constraints are inequality ones; ie components of solution vector should all be > 0. Next, the second equation seems a bit confusing. If it is an inequality constraint, it is redundant with the last three. If it is an equality constraint, it conflicts with the first one. Maybe it shouldn't be there or should be changed in some way.

------------ edit ---------------------------------------------

I responded to your first post a bit too quickly and now understand that all your constraints are inequality constraints so you can use the default for meq. It still seems to me that the second constraint is redundant with the remaining ones but that doesn't seem to be causing any problems so its not important for now. Your edit also helped me better understand your problem and I agree that your example with the R matrix should be solvable with the given constraints. It could be that the size of the matrix elements in R may be causing solution problems for solve.QP so I tried a more general nonlinear constrained optimizer, constrOptim. This gives solutions for both your R and R2 matrices with the R2 solution very close to the solve.QP solution for R2.

R2 = matrix( c( 0.05365071,-0.06364421,-0.04102565, 0, 0.08423283,-0.04048879,0,0,0.09659707), ncol=3, byrow=T)
R = matrix( c( 2.231113e-05,-4.816095e-05,-5.115287e-05, 0,2.989584e-05,4.212173e-06,0,0, 5.504990e-05), ncol=3, byrow=T)
d = as.matrix(c(57621264,78057622,171342351),ncol=1)
start <- rep(1/(ncol(R)+1), ncol(R))
min_fn <- function(b, dvec, Dmat)  -t(dvec)%*%b +t(b)%*%Dmat%*%b/2
grad_min_fn <- function(b, dvec, Dmat) -dvec + Dmat%*%b
b = c(-1., 0, rep(0,ncol(R)))
C = cbind(rep(1,ncol(R)), diag(ncol(R)))
C = cbind(-rep(1,ncol(R)),C)
D <- t(solve(R))%*%solve(R)
constrOptim(theta=start, f=min_fn, grad=grad_min_fn, ui=t(C), ci=b,     control=list(reltol=10*.Machine$double.eps), 
        dvec=d, Dmat=D )

solve.QP solution for R2
$solution
[1] -1.025463e-10  0.000000e+00  1.000000e+00

constrOptim solution for R2
$par
[1] 2.479855e-15 1.178762e-14 1.000000e+00

as well as giving the solution for R

 $par
[1] 9.272547e-17 5.958225e-14 1.040137e-01

constrOptim gives more information about its path to the solution than solve.QP so might be more helpful for numerically sensitive problems.
I don't know if constrOptim will work for you as a replacement for solve.QP but at least it may help you investigate properties of solutions to your problem.