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}}
