Creating random portfolios with weights $x_i$ can be thought of as sampling from the surface of a simplex given by $$Ex = f$$ and $$Ax \le b$$ Where $E$ and $A$ are constraint matrices for equality and inequality constraints, and $f$ and $b$ are solutions for some portfolio you want to match on, $x_0$. Adjusting the previous equations to reflect this gives us: $$Ex = f = Ex_0$$ and $$Ax \le b = Ax_0.$$
The equality constraints are relatively easy to satisfy, the solution is simply $x = x_0 + Z r$, where $Zr$ is in the null-space of $E$, and this can be the basis of a Monte Carlo random-walk. However there is no straightforward, mathematical way to link the inequality constraints with the equality constraints. The solution seems to be to check if your walk takes you over the boundary, and then reflect back over it, however as the number of $x$'s gets very large, the number of faces to reflect over grows with $n$ and this starts to pose a problem with computation time. Is there a better way to do this algorithm, or is this problem just hard?
Edit: Here is some sample R code with the random walk and reflecting over the boundaries, only handing the inequality constraint that all $x$'s must be positive:
require(MASS)
getWeights <- function(Emat, x0, n, verbose = FALSE) {
Z = Null(t(Emat))
ret = matrix(0, nrow = length(x0), ncol = n + 1)
## Would it be better to use apply here?
nc = ncol(Z)
mn = mean(x0)
ret[, 1] = x0 + Z %*% rnorm(nc, 0, mn)/sqrt(nc)
k = 0
if(verbose) cat("Created Vectors: ")
if(verbose) cat(paste(k))
for (i in 2:(n + 1)) {
ret[, i] = ret[, i - 1] + Z %*% rnorm(nc, 0, mn)/sqrt(nc)
m = k + 1;
while(any(ret[, i] < 0)) {
reflection = rep(0, ncol(Emat))
reflection[which(ret[, i] < 0)] = ret[, i][which(ret[, i] < 0)]
for (j in 1:ncol(Z)) {
ret[, i] = ret[, i] - 2* Z[, j] * (reflection %*% Z[, j])/sqrt(Z[,
j] %*% Z[, j])
}
##for(i in 1:nchar(paste(k))) cat("\b")
##if(verbose) cat(paste(m))
##k = m
##m = k + 1
}
if(verbose) for(i in 1:nchar(paste(k))) cat("\b")
if(verbose) cat(paste(m))
k = m
}
ret = ret[, 2:(n + 1)]
if(verbose) cat("\n")
return(ret)
}
Emat = matrix(1, ncol = 1000, nrow = 1)
x0 = rep(1/1000, 1000)
w = getWeights(Emat, x0, 1000, TRUE)
As you can see this code simply goes too slow. (Do you think it would be best to try and implement this code faster, using C and multiple cores, or would it be more worthwhile changing the algorithm?)