Let's say I have the following quadratic programming problem where the numbers are pretty big:
> require(quadprog)
> Q <- diag(c(34890781627, 34890781627, 34890781627, 34890781627, 34890781627))
> c <- c(133013236723, 29459621018, 31362634710, 24032348447, 23332381578)
> ( A <- t(kronecker(diag(1,5),as.matrix(c(1,-1)))) )
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 -1 0 0 0 0 0 0 0 0
[2,] 0 0 1 -1 0 0 0 0 0 0
[3,] 0 0 0 0 1 -1 0 0 0 0
[4,] 0 0 0 0 0 0 1 -1 0 0
[5,] 0 0 0 0 0 0 0 0 1 -1
> ( e <- rep(c(0,1),5) )
[1] 0 1 0 1 0 1 0 1 0 1
The constraint is Ax >= -e: I am constraining all the variables to be within [0,1]. I guess because the numbers involved here are too huge, it causes some numerical issues. If I run
> solve.QP(Q,c,A,bvec=-e)$solution
I will get Error in solve.QP(Q, c, A, bvec = -e) : constraints are inconsistent, no solution!
When I try to divide everything by a large number to shrink these crazy numbers then the function is able to generate some outputs, but NOT CORRECT:
> 1E4*solve.QP(1E-8*Q,1E-4*c,A,bvec=-1E-4*e)$solution
[1] 1 1 1 1 1
I think the correct answer is
> w <- (1/diag(Q))*c
> w[w>1] <- 1
> w[w<0] <- 0
> w
[1] 1.0000000 0.8443382 0.8988803 0.6887879 0.6687263
I don't really need to QP here because Q matrix is diagonal, but I just want to use this as an illustration and get some advice on handling the numerical issues related to solve.QP(). Let's imagine that I have a dense Q matrix, what should I do in order to get the correct result?
Thanks in advance!
Sorry guys, I just realized that I made a mistake when trying to shrink matrix Q and vector c.
min x^T Q x /2 - c^T x is equivalent to min a * (x^T Q x /2 - c^T x) with any positive a. So all I need to do is multiply both Q and c by the same small number, like the following:
> solve.QP(1E-8*Q,1E-8*c,A,bvec=-e)$solution
[1] 1.0000000 0.8443382 0.8988803 0.6887879 0.6687263
And this time it gives the correct result. Previously I was thinking of transforming x = ay, but my logic is incorrect. The correct way is the following:
min x^T Q x /2 - c^T x is equivalent to min y^T (a^2*Q) y /2 - a*c^T y with any positive a. And the constraint will also change accordingly to Ay >= -e/a. So you can also do
1E-4*solve.QP(1E-8*Q,1E-4*c,A,bvec=-1E4*e)$solution
[1] 1.0000000 0.8443382 0.8988803 0.6887879 0.6687263
In general, when should we consider this numerical issue and what is the optimal shrink size a? Or is there ways other than shrinking that can better solve this problem?