[R] not positive definite D matrix in quadprog

[Molins, Jordi]

Hello to everybody,

I have a quadratic programming problem that I am trying to solve by various
methods. One of them is to use the quadprog package in R.

When I check positive definiteness of the D matrix, I get that one of the
eigenvalues is negative of order 10^(-8). All the others are positive. When
I set this particular eigenvalue to 0.0 and I recheck the eigenvalues in R,
the last eigenvalue is positive of order 10^(-12). I try to use solve.QP,
but I get an error message that matrix D in quadratic function is not
positive definite. For reference, a fully R session is listed below.

Is 10^(-12) too close to 0? i.e., does R consider that with an eigenvalue of
order +10^(-12) the matrix is not positive definite but positive
semidefinite?

In general, has somebody know a way (in R or outside R, maybe in c++) to
solve quadratic programming with  positive semidefinite matrices? 

In particular, my problem is not so hard: given y an n x 1 matrix, and beta
an n x m matrix,  I want to find an  m x 1 matrix x s. t. sum(y - beta *
x)^2 is minimum. The particularity is that I want to impose restrictions on
x: all x components should be between 0 and 1, and there are also
constraints of the type A x = b, where A and b have the necessary dimensions
to ensure consistency.

I have tried with some other packages, and they do not give a correct
solution when the system increases in size (e.g., 24 variables and 9
constraint equations) ... some idea?

thanks!

Jordi

_____________________________


The problem:
 library(MASS)
 library(quadprog)

 D <-
matrix(c(439.5883658,438.8445615,438.1007572,2430.285506,2426.162884,2422.04
0262,44.21800696,44.14261394,
 
438.8445615,438.1020157,437.3594699,2426.173348,2422.057702,2417.942056,44.1
43188,44.06792255,
 
438.1007572,437.3594699,449.6727418,2445.212326,2542.83573,2643.780669,50.19
455336,52.04059805,
 
2430.285506,2426.173348,2445.212326,13491.19467,13614.55046,13737.90625,253.
4897678,255.745654,
 
2426.162884,2422.057702,2542.83573,13614.55046,14687.86142,15923.99043,313.8
180838,336.4239658,
 
2422.040262,2417.942056,2643.780669,13737.90625,15923.99043,19107.7405,410.9
729841,472.5104919,
 
44.21800696,44.143188,50.19455336,253.4897678,313.8180838,410.9729841,9.5462
51262,11.57677661,
 
44.14261394,44.06792255,52.04059805,255.745654,336.4239658,472.5104919,11.57
677661,14.51245153),8,8)

 D.vectors <- eigen(D,only.values=F)$vectors
 D.values <- eigen(D,only.values=F)$values

#the last value is negative
 D.values
[1]  4.609489e+04  2.458166e+03  8.232288e+01  1.961199e+00  5.976441e-01
[6]  2.810968e-01  1.253157e-09 -2.685763e-08

 D.quad <- matrix(0,8,8)
 diag(D.quad) <- D.values

#checking that the eigenvalue decomposition works fine
 D.vectors%*%D.quad%*%ginv(D.vectors)

 D.quad[8,8]
[1] -2.685763e-08
 D.quad[8,8] <- 0.0

#checking; nothing changes too much
D.vectors%*%D.quad%*%ginv(D.vectors)

#now all eigenvalues are positive:
 D.values.new <-
eigen(D.vectors%*%D.quad%*%ginv(D.vectors),only.values=F)$values
 D.values.new
[1] 4.609489e+04 2.458166e+03 8.232288e+01 1.961199e+00 5.976441e-01
[6] 2.810968e-01 1.253140e-09 1.428534e-12

Dmat <- D.vectors%*%D.quad%*%ginv(D.vectors)

dvec <-
-c(-2910.533769,-2905.609008,-3012.223863,-16274.97455,-17222.46423,-18380.6
391,-357.8878464,-379.6371849)

#this ensures that coefficients are positive:
 Amat <- matrix(0,8,8)
 diag(Amat) <- 1
 bvec <- rep(0,8)

#it says D is not positive definite ...
 solve.QP(Dmat,dvec,Amat,bvec=bvec)
Error in solve.QP(Dmat, dvec, Amat, bvec = bvec) : 
        matrix D in quadratic function is not positive definite!


--------------------------------------------------------------------------------
The information contained herein is confidential and is inte...{{dropped}}