Which library can efficiently deal with corner solutions in R?

Question

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!

Answer 1

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.