A O(n^2) procedure for uniform random sampling of points in R^n contained in the hypercube and an arbitrary hyperplane

In this note we consider a procedure for the uniform random sampling of points in Rn that satisfy two constraints:

1. Each point lies on a given hyperplane. Equivalently, each point satisfies a given linear equation.
2. Each point lies in the interior of the hypercube at the origin. Equivalently, the extension of each point along any dimension is between -1 and 1.

We first note a procedure for uniform random sampling of points in Rn that satisfy only the first constraint of being contained by a hyperplane:

Procedure 1

• Input: An array of N reals representing coefficients of a linear equation describing a hyperplane.
• Output: An array of N reals representing a point in Rn that satisfies the input equation.

1. For each of the first N-1 dimensions uniformly sample a real number and place it in the output array.
2. For the remaining N’th dimension solve the input equation with the prior N-1 variables assigned using the samples from Step 1.

We take it as self-evident this uniformly samples points on the input hyperplane.

Points generated in this way are not guaranteed to lie within a hypercube. The most straightforward way to modify Procedure 1 to handle this would be to ensure that in step 1 of Procedure 1 we are sampling numbers within [-1, 1]. However in step 2 it is not guaranteed that the last variable, determined by solving the hyperplane equation, will be contained within [-1, 1].

We will describe a solution to this problem that gives a quadratic computational cost per sample. It is not known (by this author) if this is optimal.

The procedure will be similar to Procedure 1 but will look ahead at each sample to determine if the range to sample from should be smaller than [-1, 1].

Let’s describe the solution procedure by example application to an input. We arbitrarily assume that we are in R4 and are given an input coefficient (C_i) array of [4, -2, 1, -3]. We want to sample the first element of our point called v_0. Consider that in the linear equation for any value of v_0 if I increase the value (e.g. by adding a small delta) the remaining three terms in the equation must collectively decrease by the same amount (the delta times C_0) to maintain the equality.

The key observation is that if we are at the point (v_0, 1, -1, 1) it is impossible for the remaining three terms to collectively decrease without breaking the box constraints. For instance the term (C_1)(v_1) = (-2)(1) cannot become less than -2 without v_1 becoming greater than 1.

We name this observation positive (or negative) saturated terms. A term with coefficient 5 and variable assignment 1 (5)(1) is positive saturated; it cannot increase anymore without the variable breaking the box constraint. (-5)(-1) is similarly positive saturated (to increase it the -1 must become less than -1). Examples of negative saturated terms (terms that cannot decrease anymore) are (5)(-1) and (-5)(1).

A breakpoint for a particular variable v_i is a point where all of the other terms in the sum are all positively saturated or all negatively saturated. If this is the case, and only if this is the case, it follows that a movement along that vi in the positive (negative) direction will force the point off of the hyperplane. These are precisely the points we need to avoid to keep our uniform samples legal.

There are always two breakpoints for v_i with a right and left direction determined directly by the sign of the coefficients. One breakpoint can be constructed as the non-v_i variables assigned all 1’s with signs matching their hyperplane coefficients, the other breakpoint is the all 1’s with signs opposite their hyperplane coefficients. Before sampling each variable we compute its two breakpoints and update the sample range if a breakpoint lies within [-1, 1] for the dimension we are sampling. Once a variable has been sampled it is assigned that value in the equation and the process iterates.

Procedure 2:

1. Iterate over the dimensions d in arbitrary order. For the current d:
1. If it is the last d solve the input equation and return. By the procedure this is guaranteed to lie in [-1, 1].
2. Else,
1. Substitute already sampled d’s (if any) into the input equation.
2. Determine left and right breakpoints based on the sign of the input coefficients.
3. Solve each breakpoint for v_k. If v_k is contained within [-1, 1] update the legal sample range correspondingly.
4.  Sample v_k uniformly from within legal range.

We claim this procedure uniformly samples points in the intersection of a hypercube and hyperplane in R^n, and runs in time quadratic in n.

Sample python code of a successful implementation:

def GetRange(C, i, rhs): 

  #Get rightward breakpoint

  bp1 = [(1 if C[i]<0 else -1)*math.copysign(1, j) for j in C] 

  #Solve v_0 for rightward breakpoint

  v_r = (rhs - sum([C[j]*bp1[j] for j in range(i+1,len(C))])) / float(C[i]) 

  #Get leftward breakpoint

  bp2 = [-1*j for j in bp1]

  #Solve v_0 for leftward breakpoint

  v_l = (rhs - sum([C[j]*bp2[j] for j in range(i+1,len(C))])) / float(C[i])    

  return max(-1, v_l), min(1, v_r)

def Sample(C, rhs):

  sample = []

  for i in range(len(C)):

    if i == len(C)-1:

      sample.append(rhs / float(C[i]))

      return sample

    else:

      lr, rr = GetRange(C, i, rhs)

      v_i = random.uniform(lr, rr)

      sample.append(v_i)

      rhs = rhs - v_i*C[i]