如何快速生成满足线性不等式组的单位超立方体随机点？

问题背景 #

变量维度n≈10，需满足约60个线性不等式（每个不等式格式为[c0,c1,...,cn]，对应c0 + c1x1 + ... + cnxn ≤ 0），可行域包含于单位超立方体0≤xi≤1。原随机采样后全量验证的方法效率极低，需优化。

优化思路与实现 #

1. 基础优化：向量运算加速验证 #

原验证函数用列表循环计算，改用numpy向量运算可大幅提升单样本验证速度，同时提前预处理不等式为numpy数组，避免重复转换。

import numpy as np
import random

# 预处理不等式为numpy数组，每个行对应一个不等式的系数[c0,c1,...,cn]
def preprocess_inequalities(inequalities):
    return np.array(inequalities, dtype=np.float64)

def satisfies_inequalities_np(point, inequalities_np):
    # point为numpy数组，计算每个不等式的左边值，判断是否全部≤0
    lhs = inequalities_np[:, 0] + inequalities_np[:, 1:] @ point
    return np.all(lhs <= 1e-8)  # 加入微小容差处理浮点误差

def random_point_in_unit_cube_np(n):
    return np.random.rand(n)

2. 拒绝采样优化：分层验证提前终止 #

将不等式随机排序，验证时一旦发现不满足的约束立即终止，减少不必要的计算。结合向量运算，进一步提升效率。

def random_point_satisfying_inequalities_fast(n, inequalities):
    inequalities_np = preprocess_inequalities(inequalities)
    
    while True:
        point = random_point_in_unit_cube_np(n)
        # 随机打乱不等式顺序，避免每次验证相同约束序列
        shuffled_ineqs = inequalities_np[np.random.permutation(len(inequalities_np))]
        valid = True
        for ineq in shuffled_ineqs:
            lhs = ineq[0] + ineq[1:] @ point
            if lhs > 1e-8:
                valid = False
                break
        if valid:
            return point.tolist()

3. MCMC采样（Metropolis-Hastings算法） #

当可行域占单位超立方体比例极低时，拒绝采样仍效率低下。此时可采用MCMC方法：从一个可行点出发，在可行域内随机游走生成候选点，接受满足约束的点，快速遍历可行域。

def find_initial_feasible_point(n, inequalities):
    # 用线性规划找任意可行点：目标函数取0，约束为所有不等式和0≤xi≤1
    inequalities_np = preprocess_inequalities(inequalities)
    # 构造线性规划的约束矩阵A和向量b：A@x ≤ b
    A = inequalities_np[:, 1:]
    b = -inequalities_np[:, 0]
    # 加入单位超立方体约束：0≤xi≤1
    A = np.vstack([A, -np.eye(n), np.eye(n)])
    b = np.hstack([b, np.zeros(n), np.ones(n)])
    
    from scipy.optimize import linprog
    res = linprog(np.zeros(n), A_ub=A, b_ub=b, bounds=(0,1), method='highs')
    if res.success:
        return res.x
    else:
        raise ValueError("No feasible point exists")

def mcmc_sample_feasible_point(n, inequalities, num_steps=1000, step_size=0.1):
    inequalities_np = preprocess_inequalities(inequalities)
    current_point = find_initial_feasible_point(n, inequalities)
    
    for _ in range(num_steps):
        # 在当前点附近生成候选点
        candidate = current_point + step_size * (np.random.rand(n) - 0.5)
        # 约束候选点在单位超立方体内
        candidate = np.clip(candidate, 0.0, 1.0)
        # 验证候选点是否满足所有不等式
        if satisfies_inequalities_np(candidate, inequalities_np):
            current_point = candidate
    
    return current_point.tolist()

4. 凸多面体顶点凸组合采样 #

线性不等式组的可行域是凸多面体，可行点可表示为顶点的凸组合（非负系数和为1）。通过随机采样顶点，再生成凸组合可得到近似均匀分布的点。

def sample_vertex(n, inequalities):
    inequalities_np = preprocess_inequalities(inequalities)
    # 加入单位超立方体的约束
    cube_constraints = []
    for i in range(n):
        cube_constraints.append([0.0] + [-1.0 if j==i else 0.0 for j in range(n)])  # xi≥0
        cube_constraints.append([-1.0] + [1.0 if j==i else 0.0 for j in range(n)]) # xi≤1
    all_constraints = np.vstack([inequalities_np, np.array(cube_constraints)])
    
    while True:
        # 随机选n个线性无关的约束求解交点
        selected_idx = np.random.choice(len(all_constraints), n, replace=False)
        selected = all_constraints[selected_idx]
        A = selected[:, 1:]
        b = -selected[:, 0]
        try:
            vertex = np.linalg.solve(A, b)
        except np.linalg.LinAlgError:
            continue  # 线性相关，重新选择
        # 验证顶点是否满足所有约束
        if np.all(vertex >= -1e-8) and np.all(vertex <= 1.0 + 1e-8):
            if satisfies_inequalities_np(vertex, inequalities_np):
                return vertex

def convex_combination_sample(n, inequalities, num_vertices=20):
    # 采样多个顶点
    vertices = [sample_vertex(n, inequalities) for _ in range(num_vertices)]
    # 生成凸组合系数（非负，和为1）
    weights = np.random.dirichlet(np.ones(num_vertices))
    # 计算凸组合
    point = np.sum(weights[:, np.newaxis] * vertices, axis=0)
    # 修正浮点误差，确保在单位超立方体内
    point = np.clip(point, 0.0, 1.0)
    return point.tolist()

方案选择建议 #

• 若可行域占单位超立方体比例≥1%：优先使用分层验证的拒绝采样，实现简单且速度快
• 若可行域占比极低：使用MCMC采样，能快速收敛到可行域内
• 若需要严格均匀分布：使用凸组合采样，但顶点采样可能耗时，适合对分布均匀性要求高的场景

内容的提问来源于stack exchange，提问作者Michael Zheng