random vector generation which meets a certain condition

Question

Suppose G is a matrix of (100 , 20 ) and h is vector of ( 100).

How do I randomly generate a vector x of size 100 that meets the condition G x <= h ?

Clearly, I can do something like the following. But, that is inefficient if I want to generate thousands of these points. Thanks.

import numpy as np

# Define G and h
G = np.random.rand(100, 20)
h = np.random.rand(100)

def is_feasible(x):
  return np.all(np.dot(G, x) <= h)

while True:
  # Generate random vector x
  x = np.random.rand(20)

  # Check if x is feasible
  if is_feasible(x):
    break

Answer 1

I believe you could try following. You have unequalities

G x <= h

x_i < 1 for any i

x_i >= 0 for any i

You find intersection convex polygon satisfying all above mentioned conditions. And from that polygon for each dimension you get

max x_i and min x_i, obviously those values being smaller than 1 and larger or equal to 0.

Then you sample your point in the hypercube defined by min/max value, and only then you apply matrix check

Pseudocode

while True:

for i in range(0,100):

x_i = min x_i + (max x_i - min x_i) U01

if G x <= h: break

Answer 2

This code generates a vector randomly checks for feasibility. If it fails subtracts random small amount from each value of the vector and tries again. That way I was able to generate 1000 non-zero vectors under 3 seconds:

import numpy as np


def is_feasible(x):
    return np.all(np.dot(G, x) <= h)


def generate(n=1000):
    out = []
    for _ in range(n):
        x = np.random.rand(20)

        while True:
            to_substract = np.random.uniform(0.01, 0.00001, size=20)
            x = np.clip(x - to_substract, 0, 1)

            # discard all-zero vector
            if np.all(x == 0):
                x = np.random.rand(20)
                continue

            if is_feasible(x):
                out.append(x)
                break

    return out


# Define G and h
G = np.random.rand(100, 20)
h = np.random.rand(100)

# generate and show first 20 vectors with most elements non-zero
out = sorted(generate(), key=lambda a: -np.sum(a > 0))
for row in out[:20]:
    print(*map(lambda v: "{:.5f}".format(v) if v != 0 else "       ", row), sep="|")

Prints (for example):

       |       |0.01784|0.00423|       |       |       |       |       |0.01945|       |       |       |       |       |       |       |0.00355|       |       
       |       |       |0.01594|0.01037|       |       |       |       |       |       |       |0.00797|       |       |0.01140|       |       |       |       
0.01059|       |       |       |0.02336|       |       |       |       |       |       |       |       |       |       |       |       |0.01277|       |       
       |       |0.03377|       |       |0.00754|       |       |       |       |0.02260|       |       |       |       |       |       |       |       |       
       |       |0.04551|       |       |       |       |0.01515|       |       |       |       |       |0.00466|       |       |       |       |       |       
       |0.02042|       |       |       |       |       |       |       |       |       |       |       |       |       |       |0.00385|0.00374|       |       
0.02001|       |0.00690|       |       |       |       |       |       |       |       |       |       |       |       |0.01926|       |       |       |       
       |       |       |0.01995|       |       |       |       |       |0.02459|0.00761|       |       |       |       |       |       |       |       |       
       |       |       |       |       |       |       |       |       |0.01517|       |       |0.00088|       |       |       |0.01188|       |       |       
       |       |       |       |       |       |       |       |       |0.01883|       |0.00382|       |0.00260|       |       |       |       |       |       
       |       |       |       |       |       |       |       |0.02289|0.00556|       |       |       |       |       |       |       |0.00122|       |       
       |       |0.00943|       |       |       |       |       |       |       |       |0.00766|       |       |0.01971|       |       |       |       |       
0.01891|       |0.00308|       |       |       |       |       |       |       |       |0.00916|       |       |       |       |       |       |       |       
       |       |       |       |       |       |       |0.00428|       |       |       |       |       |0.00848|       |       |0.00597|       |       |       
       |       |       |0.03437|       |       |       |       |0.01687|       |       |       |       |       |       |0.00748|       |       |       |       
       |       |       |0.02499|       |       |       |       |       |0.01509|       |       |0.00726|       |       |       |       |       |       |       
       |       |       |       |       |       |       |       |       |       |0.01138|       |       |       |       |       |0.00951|       |       |0.01423
       |       |0.00362|       |       |       |       |0.02474|       |       |       |       |       |       |       |       |       |0.00115|       |       
       |       |       |0.00695|       |       |0.01630|       |       |       |       |       |       |       |       |0.00665|       |       |       |       
       |       |       |0.00485|       |0.00809|0.00404|       |       |       |       |       |       |       |       |       |       |       |       |       

real    0m2,585s
user    0m0,007s
sys     0m0,008s

Answer 3

Another possible approach is to use linear programming, as @pjs suggests in a comment, where the solutions of the linear program are the randomly generated samples. The objective function had to be perturbed, to get different solutions at each run.

from scipy.optimize import linprog

G = np.random.rand(100, 20)
h = np.random.rand(100)

# Perturbed objective function
c = np.random.normal(0, 0.01, 20)

# Using linear programming
z = linprog(c, A_ub=G, b_ub=h, method='highs')

if z.success:
    x = z.x
    print(x)