Is creating constrained random portfolios a hard problem?

Question

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?)

Answer 1

The approach of reflecting is expensive, since the $d$-simplex has $d$ maximal faces, all of which have to be checked for intersection at each step. Additionally, if the random walk moves into a corner, the number of moves which have to be discarded can become very high. Depending on the configuration of the constraints this could well be your best solution.

If I understood your question correctly, the question boils down to sampling from the intersection of the boundary of one simplex ($E$) and the volume of another ($A$) simultaneously.

Suggestions:

• Rejection: if the admissible regions for the two conditions are not too different, it could be efficient to create points inside $A$, test for inclusion in $E$ and reject the point if not.
• Restart: if the random walk creates too many discarded moves, start from a new location using the rejection method above.
• Projection: Create an interior point for $E$ and check if the closest point on $A$ is still inside $E$.
• Adaptive step-size: You can calculate the distance from the boundary of $A$ and $E$ and choose half the step size, of course this will add a bias to the generation proces, this could be a problem.

Algorithm for one simplex: There is a nice trick for the unit simplex by Donald Rubin, which is discussed on the computer science site, to sample from the unit simplex. The result can be scaled afterwards to match a general simplex:

• Create $n-1$ uniformly distributed random values between 0 and 1 and sort them.
  [Example: 0.1, 0.3, 0.5, 0.55]
• Add 0 and 1 to the list.
  [Example: 0, 0.1, 0.3, 0.5, 0.55, 1]
• Use the differences between the values, they will sum to 1.
  [Example: 0.1, 0.2, 0.2, 0.05, 0.45]

Edit: Strictly speaking it is not a hard problem, since it can be solved in polynomial time; it would be though if you were looking for integer weights.