Randomly construct linear combination that is within bounds over two basis

Question

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

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.clf()
plt.imshow(heatmap.T, extent=extent, origin='lower')
plt.plot(v_x, v_y, color='r')
plt.show()