带约束的拉丁超立方（Latin Hypercube）采样优化问询

问题背景 #

我想用拉丁超立方采样生成约500个包含36个变量的样本，确保参数空间的良好覆盖，同时要求每个样本满足协方差矩阵正定的约束。目前我基于scipy.stats.qmc实现了代码，但遇到两个问题：

1. 生成符合要求的样本耗时极长；
2. 经过约束过滤后，最终样本可能失去拉丁超立方采样原本的空间覆盖特性，违背了使用LHS的初衷。

我的实现代码如下：

from scipy.stats import qmc
import numpy as np

def unpack(params: np.ndarray) -> tuple[
    float,       # p
    np.ndarray,  # variable means
    np.ndarray,  # standard errors
    np.ndarray,  # correlation coefficients
    np.ndarray,  # covariance
]:
    (p,), means, dev_diag, X_triu = np.split(params, (1, 8, 15))
    dev = np.diag(dev_diag)
    n = dev_diag.size
    X = np.zeros((n, n))
    X[np.triu_indices(n=n, k=1)] = X_triu
    X += X.T + np.eye(n)
    cov = dev @ X @ dev
    return p, means, dev, X, cov

def positive_definite(params: np.ndarray) -> np.ndarray:
    p, means, dev, X, cov = unpack(params)
    return np.real(np.linalg.eigvals(cov))

def init_pop_qml(bounds, popsize,multipler=10):
    limits = np.array(bounds, dtype='float').T
    parameter_count = len(bounds)
    num_population_members = popsize * parameter_count

    sampler = qmc.LatinHypercube(d=parameter_count,seed=123456)
    sample = sampler.random(n=num_population_members*multipler)
    l_bounds = limits[0]
    u_bounds = limits[1]
    samples = qmc.scale(sample, l_bounds, u_bounds)

    checks = np.zeros((num_population_members*multipler, 7))
    for row in range(num_population_members*multipler):
        check = positive_definite(samples[row, :])
        checks[row, :] = check

    population = samples[(checks > 0).all(axis=1), :]

    while population.shape[0] < num_population_members:
        sample2 = sampler.random(n=num_population_members*multipler)
        samples2 = qmc.scale(sample2, l_bounds, u_bounds)
        checks2 = np.zeros((num_population_members*multipler, 7))
        for row in range(num_population_members*multipler):
            check2 = positive_definite(samples2[row, :])
            checks2[row,:] = check2
        population2 = samples2[(checks2 > 0).all(axis=1), :]
        population = np.vstack((population, population2))
        print(population.shape)

    return population[0:num_population_members,:]

bounds = (((0.0, 1.0),) * 8 + ((0.0, 1.0),) + ((0.0, 0.5),) * 6 + ((-1.0,1.0),)*21)
init_pop=init_pop_qml(bounds, 15, 10)
print(init_pop)

优化建议 #

针对你遇到的两个问题，我给出几个实用的优化方向：

1. 大幅提升约束检查的效率

当前代码里最耗时的部分是逐行循环计算协方差矩阵的特征值，完全可以用更高效的方式替代：

• 用Cholesky分解判断正定：正定矩阵能成功完成Cholesky分解，反之则会报错。这个操作比计算所有特征值快得多，我们可以利用这一点改写约束检查函数：

def is_positive_definite(params: np.ndarray) -> bool:
    _, _, dev, X, cov = unpack(params)
    try:
        np.linalg.cholesky(cov)
        return True
    except np.linalg.LinAlgError:
        return False

• 进一步用矢量化操作替代循环：可以用np.apply_along_axis批量处理样本矩阵，避免手动逐行循环的开销。

2. 从“先采样后过滤”转向“带约束采样”

单纯的过滤会破坏LHS的空间覆盖性，我们可以在采样阶段就考虑约束：

• 优化接受-拒绝逻辑：不要一次性生成大量样本，改为每次生成小批量样本，检查约束后保留合格的，直到满足所需数量。这样既减少内存占用，也能避免无效计算。比如修改后的采样函数：

def init_pop_qml_optimized(bounds, popsize, batch_size=50):
    limits = np.array(bounds, dtype='float').T
    l_bounds, u_bounds = limits
    parameter_count = len(bounds)
    num_needed = popsize * parameter_count
    sampler = qmc.LatinHypercube(d=parameter_count, seed=123456)
    
    population = []
    while len(population) < num_needed:
        # 生成小批量样本
        sample = sampler.random(n=batch_size)
        scaled_samples = qmc.scale(sample, l_bounds, u_bounds)
        # 筛选合格样本
        valid = [s for s in scaled_samples if is_positive_definite(s)]
        population.extend(valid)
        print(f"Collected {len(population)}/{num_needed} samples")
    
    return np.array(population[:num_needed])

• 尝试条件拉丁超立方采样（CLHS）：这种方法本身就支持在采样时纳入约束，能更好地保持LHS的空间覆盖特性。你可以手动实现CLHS的核心逻辑：每次选择一个未采样的参数组合，确保它满足约束，同时与已采样样本的空间分布尽可能均匀。

3. 简化约束验证的前置逻辑

你的协方差矩阵是dev @ X @ dev，其中dev是对角矩阵（元素在0-0.5之间，都是正数），X是对称的相关系数矩阵（对角线为1）。根据线性代数知识，只要X是正定矩阵，那么协方差矩阵必然是正定的。所以我们可以直接检查X的正定性，省去构建协方差矩阵的步骤，进一步减少计算量：

def is_correlation_matrix_positive_definite(params: np.ndarray) -> bool:
    _, _, dev_diag, X_triu = np.split(params, (1, 8, 15))
    n = dev_diag.size
    X = np.zeros((n, n))
    X[np.triu_indices(n=n, k=1)] = X_triu
    X += X.T + np.eye(n)
    try:
        np.linalg.cholesky(X)
        return True
    except np.linalg.LinAlgError:
        return False

4. 调整采样参数

• 不用固定multipler=10，可以根据约束的严格程度动态调整批量大小：如果约束较严，批量小一点，循环次数多一点，但每次处理更快；如果约束宽松，可以适当增大批量。
• 考虑使用qmc.LatinHypercube的centered=True参数，让采样点更均匀分布，后续过滤后保留的覆盖性也会更好。

备注：内容来源于stack exchange，提问作者jasmine