给定100×20矩阵G与100维向量h，如何高效生成满足Gx≤h的随机向量x？

问题分析 #

你当前的随机生成后校验方法，本质是在整个[0,1]^20空间盲目采样，当可行域占比极低时，效率会指数级下降，完全不适合生成数千个样本的场景。下面是几种更高效的解决方案：

方法1：Hit-and-Run 采样（凸可行域首选，生成均匀分布样本） #

Hit-and-Run是针对凸集的高效采样算法，核心逻辑是：先找到一个可行初始点，然后在可行域内随机生成方向，沿着方向找到可行边界，再在边界区间内随机选取新点，重复此过程快速遍历可行域，生成均匀分布的样本。

实现代码 #

import numpy as np
from scipy.optimize import linprog

def hit_and_run_sample(G, h, num_samples=1000, init_x=None):
    n = G.shape[1]
    samples = []
    
    # 生成初始可行点（若未提供）
    if init_x is None:
        # 通过线性规划找可行点：最小化0*x，约束Gx ≤ h，x∈[0,1]（可根据实际需求调整边界）
        res = linprog(np.zeros(n), A_ub=G, b_ub=h, bounds=[(0,1)]*n)
        if not res.success:
            raise ValueError("可行域为空，请检查G和h的取值")
        init_x = res.x
    
    current_x = init_x
    samples.append(current_x.copy())
    
    for _ in range(num_samples-1):
        # 生成随机单位方向向量
        direction = np.random.randn(n)
        direction /= np.linalg.norm(direction)
        
        # 计算当前点沿方向的可行区间[t_min, t_max]
        denominators = G @ direction
        numerators = h - G @ current_x
        
        t_max_list = []
        t_min_list = []
        
        # 处理Gx ≤ h的约束
        for d, num in zip(denominators, numerators):
            if d > 1e-9:
                t_max_list.append(num / d)
            elif d < -1e-9:
                t_min_list.append(num / d)
            else:
                if num < -1e-9:
                    raise ValueError("当前点不可行，算法执行出错")
        
        # 处理x∈[0,1]^n的边界约束
        for i in range(n):
            dir_i = direction[i]
            if dir_i > 1e-9:
                t_max_list.append((1 - current_x[i]) / dir_i)
                t_min_list.append(-current_x[i] / dir_i)
            elif dir_i < -1e-9:
                t_max_list.append(-current_x[i] / dir_i)
                t_min_list.append((1 - current_x[i]) / dir_i)
        
        t_max = min(t_max_list) if t_max_list else np.inf
        t_min = max(t_min_list) if t_min_list else -np.inf
        
        # 在可行区间内随机取t，更新当前点
        t = np.random.uniform(t_min, t_max)
        current_x += t * direction
        samples.append(current_x.copy())
    
    return np.array(samples)

# 测试示例
G = np.random.rand(100, 20)
# 构造确保可行域非空的h（可根据实际需求调整）
h = G @ np.ones(20)*0.6 + np.random.rand(100)*0.2

samples = hit_and_run_sample(G, h, num_samples=1000)
# 验证所有样本满足约束（浮点精度允许微小误差）
assert np.all(G @ samples.T <= h[:, np.newaxis] + 1e-6)

方法2：投影法（实现简单，适合可行域占比不极低的场景） #

先在扩大的随机空间生成点，再将其投影到可行域内。投影问题可转化为二次规划：最小化||x - x0||²，约束Gx ≤ h且x∈[0,1]^n，实现成本低。

实现代码 #

import numpy as np
from scipy.optimize import minimize

def project_to_feasible(x0, G, h):
    n = G.shape[1]
    # 目标函数：最小化投影点与原始点的距离平方
    def objective(x):
        return np.sum((x - x0)**2)
    
    # 约束条件
    constraints = [{'type': 'ineq', 'fun': lambda x: h - G @ x}]
    bounds = [(0, 1) for _ in range(n)]
    
    res = minimize(objective, x0, method='SLSQP', constraints=constraints, bounds=bounds)
    if res.success:
        return res.x
    else:
        # 若投影失败，重新生成初始点重试
        return project_to_feasible(np.random.rand(n), G, h)

def generate_samples_via_projection(G, h, num_samples=1000):
    samples = []
    for _ in range(num_samples):
        # 扩大随机生成范围，提升投影到可行域的概率
        x0 = np.random.rand(G.shape[1]) * 1.5
        x_feasible = project_to_feasible(x0, G, h)
        samples.append(x_feasible)
    return np.array(samples)

# 测试示例
samples_proj = generate_samples_via_projection(G, h, num_samples=1000)
assert np.all(G @ samples_proj.T <= h[:, np.newaxis] + 1e-6)

方法选择建议 #

• 若需要均匀分布的样本，优先选Hit-and-Run算法，它是凸集采样的标准方案，效率远高于随机校验，样本分布更均匀。
• 若追求实现简单，投影法更易理解和编写，但样本分布均匀性不如前者，且当可行域极小时，投影计算成本会上升。
• 所有方法的前提是可行域非空，可先用线性规划（如示例中的linprog）验证是否存在可行解。

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