# 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:

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
• You want to find a random point in a polytope. There are known algorithms for this. – Yuval Filmus Oct 28 '18 at 1:16
• Thank you for the pointers. Requirement 1 = cube. Requirement 2 = parallelepiped. Intersection = convex ploytope. Which kind of polytope should I google? – R zu Oct 28 '18 at 14:57
• Found solution for convex polytope with barycentric coordinate. Thanks! – R zu Oct 28 '18 at 15:09
• 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. – R zu Oct 28 '18 at 15:33
• There are methods that work for any polytope. – Yuval Filmus Oct 28 '18 at 15:45

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()

• A good observation. Usually, the simplest algorithm to implement is rejection sampling and its performance is often acceptable. – Discrete lizard Apr 26 '19 at 17:13