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:

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

  • $\begingroup$ Can you be more precise on what you want to do ? Your constraints make me think you want to optimize something but it not very clear. $\endgroup$ Jun 20, 2013 at 13:15
  • $\begingroup$ I can see how you would think it is an optimization problem because the constraints are the same as an LP problem, however, I just want to sample from all possible portfolios given that those constraints are true. $\endgroup$
    – Mike Flynn
    Jun 20, 2013 at 16:00
  • $\begingroup$ There is somerhing wrong here: the text states you require $Ax < b$, but the code checks for values $<0$, which would imply $Ax \ge 0$. In this case there is a faster way to solve your problem. $\endgroup$
    – user1157
    Feb 12, 2014 at 22:03

1 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.


  • 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.