Search code examples
rmathematical-optimization

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


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!


Solution

  • 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.