Given rank deficient matrix $A$, I want to randomly construct vectors $\vec{x}$ such that:

  1. $0 \le x_{j} \le 1$
  2. $0 \le b_{j} \le 1$ where $\vec{b} = A\vec{x}$

Matrix $A$ is about 10 x 15. I want to generate $10^{5}$ $\vec{x}$ as an initial guess in each iteration of a solver.

Below is my attempts.

Least square and rejection

  1. Generate $\vec{y}$ with entries between 0 and 1.
  2. Generate $\vec{b}$ with entries between 0 and 1.
  3. Minimize $(\vec{y} - \vec{x})^{2}$ subjected to constraint $A\vec{x} = \vec{b}$.
  4. Reject $\vec{x}$ that are out of bounds.

Problem: the rejection rate is very high. Very impractical even for a 10 x 15 matrix $A$.

Step 2 is equivalent to solving

$\begin{bmatrix} A & 0\\ I & A^{T} \end{bmatrix} \begin{bmatrix} \vec{x}\\ \vec{\lambda}\\ \end{bmatrix} = \begin{bmatrix} \vec{b}\\ \vec{y}\\ \end{bmatrix} $

Constrained Optimization

  • Generate random $\vec{y}$ within bounds.
  • Objective: $(\vec{x} - \vec{y})^2$.
  • Constraint: $A\vec{x} = b$
  • Bounds: between 0 and 1.

Problem: Optimization is too slow and often gives result out of bounds. Impractical.

Other ideas

  • Null space
  • Simplex, random planes, and projection
  • Hit and run sampling
  • $\begingroup$ You want to find a random point in a polytope. There are known algorithms for this. $\endgroup$ Oct 28, 2018 at 1:16
  • $\begingroup$ Thank you for the pointers. Requirement 1 = cube. Requirement 2 = parallelepiped. Intersection = convex ploytope. Which kind of polytope should I google? $\endgroup$
    – R zu
    Oct 28, 2018 at 14:57
  • $\begingroup$ Found solution for convex polytope with barycentric coordinate. Thanks! $\endgroup$
    – R zu
    Oct 28, 2018 at 15:09
  • $\begingroup$ But the method with barycentric coordinate is not uniform sampling. For each $A$ and $\vec{b}$, the polytope that bounds $\vec{x}$ is different. Since I only need one point for each polytope, hit-and-run sampling might need quite a few iterations to spread to all space within the polytope. $\endgroup$
    – R zu
    Oct 28, 2018 at 15:33
  • $\begingroup$ There are methods that work for any polytope. $\endgroup$ Oct 28, 2018 at 15:45

1 Answer 1


Using barycentric coordinate and keeping sum of coordinate equal to 1 do not sample a convex polytope uniformly.


import numpy as np
import matplotlib.pyplot as plt

N = 1000000  # Number of samples
V = 5 # Generate points within a V-gon

# generate columns that each sum to 1
s = np.exp(-np.random.uniform(0, 1, size=(V, N)))
s /= s.sum(axis=0)

# generate vertex of regular polygon
r = 2
theta = np.linspace(0, 2*np.pi, V+1)[1:]
v_x = r*np.cos(theta)
v_y = r*np.sin(theta)

# generate random points within the polygon
# by linear combination of coordinates of
# vetex where the coefficients sum to 1.
x = (v_x[:, np.newaxis] * s).sum(axis=0)
y = (v_y[:, np.newaxis] * s).sum(axis=0)

# plot density
heatmap, xedges, yedges = np.histogram2d(x, y, bins=50)
extent = [-r, r, -r, r]

plt.imshow(heatmap.T, extent=extent, origin='lower')
plt.plot(v_x, v_y, color='r')
  • 1
    $\begingroup$ A good observation. Usually, the simplest algorithm to implement is rejection sampling and its performance is often acceptable. $\endgroup$
    – Discrete lizard
    Apr 26, 2019 at 17:13

Your Answer

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

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