[R] Constraint Optimization with constrOptim

Ravi Varadhan ravi.varadhan at jhu.edu
Tue Sep 18 18:25:05 CEST 2012


Hi Raimund,



You can use spg() in the BB package. This requires that you have a projection rule for projecting an infeasible point onto the simplex defined by:  0 <= x <= 1, sum(x) = 1.



I show you here how to project onto the unit simplex.  So, this is how your problem can be solved:



require(BB)

 

proj.simplex <- function(y, c=1) {

#####

# project an n-dim vector y to the simplex S_n

# S_n = { x | x \in R^n, 0 <= x <= c, sum(x) = c}

# Copyright:  Ravi Varadhan, Johns Hopkins University

# August 8, 2012

#####

n <- length(y)

sy <- sort(y, decreasing=TRUE)

csy <- cumsum(sy)

rho <- max(which(sy > (csy - c)/(1:n)))

theta <- (csy[rho] - c) / rho

return(pmax(0, y - theta))

}

#Rendite - data
dataset <- matrix(c(0.019120, 0.033820, -0.053180, 0.036220, -0.021480, -0.029380,
-0.012180, -0.076980, -0.060380,
0.038901, -0.032107, -0.034153, -0.031944, 0.006494, -0.021062,
-0.058497, 0.050473, 0.026086,
0.004916, -0.015996, 0.003968, 0.030721, 0.034774, -0.004972,
0.019459, 0.000196, -0.001849), nrow=3, ncol=9, byrow=TRUE)

#Computation of the variance-covariance matrix
covarianceMatrix <- ((t(dataset)) %*% dataset) / 3



fkt <- function(x, covMat) {t(x) %*% covMat %*% x}
startingEstimates = rep(1/9, 9)



ans <- spg(startingEstimates, fn=fkt, project=proj.simplex, covMat = covarianceMatrix)



ans

sum(ans$par)



ans <- spg(runif(9), fn=fkt, project=proj.simplex, covMat = covarianceMatrix) # non-uniqueness of the optimal solution



ans

sum(ans$par)



# You see that there is an infinite number of solutiions. Because your covariance matrix is rank-deficient, it only has rank 3.



# You get a unique solution when the covariance Matrix is positive-definite

aMat <- matrix(rnorm(81), 9, 9)

covMat <- aMat %*% t(aMat)



ans <- spg(runif(9), fn=fkt, project=proj.simplex, covMat=covMat) # non-uniqueness of the



ans

sum(ans$par)


Hope this is helpful,

Ravi




More information about the R-help mailing list