2
$\begingroup$

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
|cite|improve this question
$\endgroup$
  • $\begingroup$ You want to find a random point in a polytope. There are known algorithms for this. $\endgroup$ – Yuval Filmus Oct 28 '18 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 '18 at 14:57
  • $\begingroup$ Found solution for convex polytope with barycentric coordinate. Thanks! $\endgroup$ – R zu Oct 28 '18 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 '18 at 15:33
  • $\begingroup$ There are methods that work for any polytope. $\endgroup$ – Yuval Filmus Oct 28 '18 at 15:45
1
$\begingroup$

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

Barycentric_coordinate_sampling

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()
|cite|improve this answer
$\endgroup$
  • 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 at 17:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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