Randomly construct linear combination that is within bounds over two basis Given rank deficient matrix $A$, I want to randomly construct vectors $\vec{x}$ such that:


$0 \le x_{j} \le 1$
$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


Generate $\vec{y}$ with entries between 0 and 1. 
Generate $\vec{b}$ with entries between 0 and 1. 
Minimize $(\vec{y} - \vec{x})^{2}$ subjected to constraint $A\vec{x} = \vec{b}$. 
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 &amp; 0\\
I &amp; 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
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()