At the moment I'm trying to do a minimization (optimization) problem in R, where I have a vector X1 that I want to approximate through a weighted average of a matrix X2 and a vector of weights w. That means I want to minimize
wg <- function(w)
{
t(X1 - X2 %*% w) %*% (X1 - X2 %*% w)
}
the constraints on the weights are w[i]>= 0 and sum(w) = 1 .
At the moment I'm using the DEoptim package to do my optimization, but I feel like it can't deal well with corner solutions.
I'm replicating a method that was used in an economics paper and in that paper almost all of the weights turned out to be zero. I expected a similar result in my case ( I want to model Arizona through a weighted average of the other states), especially due to the heterogeneity in the economic situation.
At the moment I feel like it's more of a problem with the DEoptim package than with my methodology and I don't really trust the results. Which other package can I use, preferrably ones that are stronger in looking for corner solutions?
my DEoptim is set up as follows:
controlDE <- list(reltol=.0000000001,steptol=150, itermax = 5000,trace = 250)
#control parameters
outDEoptim <- DEoptim(fn = wg, lower = rep(0, N), upper = rep(1, N),
control = controlDE)
Any help would be much appreciated!
A stochastic solver such as DEoptim will by nature have difficulties finding optimal solutions on lower dimensional subsets such as the one defined by sum(w) = 1.
There is a first, not quite correct way of doing this by reducing the problem to (n-1) dimensions by setting w <- c(w, 1-sum(w)). The last component might get less than 0, but normally it won't. Now apply DEoptim or optim:
set.seed(1357); m <- 4; n <- 5
X2 <- round(matrix(runif(20), m, n), 2)
X1 <- X2 %*% c(0, 1, 0, 0, 0) # solution c(0,1,0,0,0)
wg <- function(w) { # length(w) == 4
w <- c(w, 1 - sum(w))
t(X1 - X2 %*% w) %*% (X1 - X2 %*% w) # sum((X1 - X2 %*% w)^2)
}
w0 <- rep(1/n, n-1) # initial point (1/n, ..., 1/n)
optim(w0, wg, lower = rep(0, n), upper = rep(1, n),
method = "L-BFGS-B", control = list(factr = 1e-8))
## $par
## [1] 0 1 0 0 # wmin = c(0,1,0,0,0)
Or you apply one of the solvers in R that can handle equality constraints, for example Rdonlp2 (on R-Forge), auglag in package alabama, or slsqp in package nloptr. But I feel this would be overshooting.