拉丁超立方体抽样（LHS）
拉丁超立方体抽样（LHS）的理论基础
LHS的基本原理
老师，拉丁超立方体抽样比普通的随机抽样有什么优势呢？只听名字就特别复杂…
名字看似复杂，但原理其实很简单。将各输入变量的值域均等分成N份，然后从每一份中恰好选一个点。普通随机抽样可能会出现某个区域点太多、某个区域点太少的情况——我们称之为聚集效应——但LHS通过分层结构来防止这种现象。
仅仅均等分割就有效果吗？
效果很明显。相比蒙特卡罗法，同样的样本数可以使方差降低30～50%。换句话说，如果蒙特卡罗需要100次模拟来达到某个精度，LHS可能用50～60次就够了。在CAE中，一次模拟往往需要几小时，所以样本数的减少意味着巨大的成本节省。
实际应用中什么时候用LHS呢？
主要有三个场景：
- 不确定性传播（UQ）: 如材料物性的变异如何影响应力分布
- 灵敏度分析: 找出哪个输入参数对输出的影响最大
- 试验设计（DOE）的初始点布置: 为Kriging等代理模型准备训练数据点
比如在汽车碰撞分析中，如果板厚、屈服应力和摩擦系数有变异，用LHS生成50个样本点可以达到蒙特卡罗100个点的精度。如果每次碰撞分析需要4小时，那就省了200小时的计算时间。
分层抽样的数理基础
怎样用数学语言描述"均等分割"这件事呢？
对于$d$维输入变量$\mathbf{X} = (X_1, X_2, \ldots, X_d)$和$N$个样本，设第$j$个变量的累积分布函数为$F_j$。LHS的核心是将各变量的周边分布等概率地分成$N$层。
第$j$个变量的第$k$层（$k = 1, 2, \ldots, N$）定义为：
从各层抽样时，任意样本点落在第$k$层的概率为：
具体地，通过均匀随机数$U_{j,k} \sim \text{Uniform}(0,1)$，第$k$个样本的第$j$个变量值为：
其中$\pi_j$是$\{1, 2, \ldots, N\}$的随机置换，对每个变量独立生成。这个置换使得各变量的层组合随机化。
也就是说，先建立"均等的架子"，再在每个架子内随机选点？
完全正确。二维情况下，想象一个$N \times N$的国际象棋棋盘，每行、每列恰好放一个点——这就是"拉丁方阵"的概念。"拉丁"这个名字就来自18世纪数学家欧拉研究的拉丁方阵，其中每个符号在每一行和每一列中恰好出现一次。
方差降低效应的证明
为什么LHS的方差会更小呢？能不能用数学证明？
当然可以。设输出$Y = g(\mathbf{X})$，估计期望值的推估量$\bar{Y}_N = \frac{1}{N}\sum_{i=1}^{N} g(\mathbf{x}^{(i)})$的方差。McKay, Beckman & Conover (1979) 的经典结果表明：
第一项$\sigma^2 / N$就是标准蒙特卡罗的方差。第二项是LHS的方差降低量，当$g$对各输入变量表现出单调性时（加法成分大）这一项就很大。因此：
等号仅当$g$对所有输入变量都完全非加法（纯交互作用）时才成立。在实际CAE问题中，加法成分几乎总是存在的，所以LHS总是优于蒙特卡罗。
也就是说，输出变化越"规律"，LHS的优势越明显？
正是这样。结构分析中应力对厚度的响应基本是线性的，所以LHS效果特别好。相比之下，乱流的瞬时值可能有复杂的交互作用，LHS的优势会略微减弱。但理论上讲，LHS永远不会比蒙特卡罗更差，所以完全可以作为默认选择。
空间填充准则（Space-Filling Criteria）
即使用了LHS，有时生成的样本点好像还是有点聚集…是这样吗？
观察很敏锐。基础LHS保证了周边分布的分层，但不保证多维空间中的均匀性。比如二维时可能出现"点在右上和左下聚集，中间区域稀疏"的现象。这就需要用空间填充准则来优化点的位置。
主要有三种准则：
1. 最大最小距离准则
最大化所有样本点对之间的最小距离：
其中$\| \cdot \|_p$是$L_p$范数（通常取$p=2$的欧几里得距离）。直观上就是"让点之间尽可能分散"。
2. 最大最小距离的$\phi_p$准则
Morris & Mitchell (1995) 提出的连续松弛版本：
当$p$增大时收敛到最大最小准则。实践中$p = 15 \sim 50$效果很好，而且函数可微分，可用梯度方法优化。
3. 差异度准则
衡量点集与均匀分布的偏离程度：
这是星差异度$D_N^{*}$，值越小说明分布越均匀。通过Koksma-Hlawka不等式，数值积分误差的上界与差异度直接关联。
最大最小距离和差异度，实务中用哪个比较好？
用Kriging或RBF做代理模型时，最大最小距离是首选。因为预测精度很大程度上取决于预测点到最近学习点的距离。而如果目标是数值积分（估计期望值），差异度准则在理论上更严格。拿不准的话就用最大最小距离，不会出大问题。
相关性控制（Correlation Control）
用LHS生成的各变量之间会不会意外产生相关性？
很好的问题。因为每个变量独立用随机置换，偶然的伪相关确实会出现。当$N$较小时相关性可能很明显。为了防止这种情况，可以用Iman & Conover (1982)的相关性控制方法。
算法的核心步骤：
- 设定目标相关矩阵$\mathbf{C}_{\text{target}}$（独立则为单位矩阵$\mathbf{I}$）
- 计算LHS的秩矩阵$\mathbf{R}$的现有秩相关矩阵$\mathbf{C}_R$
- 对目标相关矩阵做Cholesky分解$\mathbf{C}_{\text{target}} = \mathbf{L}\mathbf{L}^T$
- 用修正矩阵$\mathbf{Q} = \mathbf{L} \cdot \mathbf{C}_R^{-1/2}$重新排列秩
修正后的秩相关矩阵会接近目标值：
输入变量本身有相关性的时候（比如弹性模量和泊松比相关），也能处理吗？
完全可以。材料物性之间的真实相关性可以通过Iman-Conover法反映在样本中，这对实务应用特别重要。结合Nataf变换或Copula方法，还能处理更复杂的相关结构。
算法与实现
基本LHS生成算法
怎样用程序实现LHS的生成？有什么具体步骤吗？
非常直白。对于$N$个样本、$d$个变量：
- 对每个变量$j = 1, \ldots, d$，生成$\{1, 2, \ldots, N\}$的随机置换$\pi_j$
- 对每个样本$i$和变量$j$，生成均匀随机数$u_{ij} \sim U(0,1)$
- 计算单位超立方体上的坐标：$s_{ij} = \frac{\pi_j(i) - 1 + u_{ij}}{N}$
- 用逆CDF变换映射到目标分布：$x_{ij} = F_j^{-1}(s_{ij})$
第3步是LHS的关键：$\pi_j(i) - 1$确定所在的层，$u_{ij}$在层内随机选位置。
计算量是多少？
基础LHS的复杂度是$O(Nd)$，与蒙特卡罗相同。优化LHS（最大最小距离准则）因为有优化循环，复杂度为$O(N^2 d \cdot T)$（$T$为迭代次数）。不过样本生成成本相对于一次CAE分析完全可以忽略。
优化LHS
怎样优化最大最小距离？全局穷举肯定不行吧？
组合爆炸是$(N!)^d$，绝对不可行。实用的方法有：
- ESE法（增强随机演化）: 随机交换两行，保留改进的配置。Jin et al. (2005) 的方法实现简单
- 模拟退火: Morris & Mitchell (1995) 提出，能避免局部最优
- 遗传算法: 大规模问题有效，但计算成本较高
- 列对交换: Viana et al. (2009) 的高速算法，逐次优化两个变量
实务中ESE法或模拟退火是标配，pyDOE和SciPy都有现成实现。
Python实现示例
Python中怎样写？
SciPy的qmc.LatinHypercube是目前的标准。优化LHS（最大最小准则）也只需一行。
import numpy as np
from scipy.stats.qmc import LatinHypercube
from scipy.stats import norm, uniform
# --- 基础LHS（5变量, 50样本）---
d, N = 5, 50
sampler = LatinHypercube(d=d, seed=42, optimization="random-cd")
unit_samples = sampler.random(n=N) # [0,1]^d 上的LHS
# --- 逆CDF变换到目标分布 ---
# 变量1: 正态分布 N(210e3, 5e3) ← 弹性模量 [MPa]
# 变量2: 正态分布 N(0.3, 0.02) ← 泊松比
# 变量3: 均匀分布 U(1.5, 2.5) ← 厚度 [mm]
x1 = norm.ppf(unit_samples[:, 0], loc=210e3, scale=5e3)
x2 = norm.ppf(unit_samples[:, 1], loc=0.3, scale=0.02)
x3 = uniform.ppf(unit_samples[:, 2], loc=1.5, scale=1.0)
# ...以此类推optimization="random-cd"是什么意思？
是中心差异度（Centered Discrepancy）最小化的优化选项。还有"Lloyd"等其他选项。默认（None）是无优化的基础LHS。实务一定要加优化选项，因为素LHS的空间填充性有时不理想。
收敛性与样本数确定
LHS的收敛速度比蒙特卡罗快吗？
阶数都是$O(1/\sqrt{N})$，但常数不同。LHS的常数小，所以达到同样误差需要的样本数较少。给出实用指南：
| 用途 | 维数$d$ | 推荐样本数$N$ | 备注 |
|---|---|---|---|
| 均值·标准差估计 | 任意 | $N \geq 10d$ | 最低标准 |
| 概率分布估计（直方图） | $\leq 10$ | $N \geq 50d$ | 尾部精度确保 |
| 代理模型构建（Kriging） | $\leq 20$ | $N = 10d \sim 50d$ | 空间填充准则优化 |
| Sobol灵敏度指标估计 | 任意 | $N \geq 100(d+1)$ | 配合Saltelli方法 |
拉丁超立方体抽样（LHS）的实务应用
CAE分析中的LHS应用流程
在实际CAE项目中怎样应用LHS？告诉我整个流程。
标准流程5步：
- 选定输入变量: 决定哪些参数要变化（物性、尺寸、载荷等）。可用灵敏度筛选（如Morris法）事先排除影响小的变量，提高效率
- 设定概率分布: 确定各变量的分布类型（正态、对数正态、均匀、Weibull等）和参数。有实测数据就用分布拟合，没有就根据专家意见
- 生成LHS样本: 用优化LHS（最大最小距离推荐）生成$N$个点。有相关变量时应用Iman-Conover法
- 运行CAE分析: 对每个样本点执行求解器。用参数化脚本自动化工作流
- 后处理统计: 计算输出的统计量（均值、标准差、百分位）。必要时构建代理模型用于进一步分析
汽车零部件设计中具体是什么变量呢？
常见的有：
- 板厚（工差±0.1mm → 正态分布$N(t_0, 0.033)$，用3σ范围）
- 屈服应力（批量变异 → 对数正态分布）
- 焊点位置（±2mm → 均匀分布）
- 摩擦系数（0.1～0.3 → 均匀分布）
4～8个变量时，50个样本通常就有很好的精度。如果单个分析需要4小时，50个样本共200小时；并行运行10个，20小时完成。蒙特卡罗要达到同样精度可能需要200～300样本，计算时间就成了瓶颈。
样本数的确定方法
"$10d$以上"这个通用法则之外，有没有更实际的判断方法？
有几个办法：
- Bootstrap方法: 用现有$N$个点重新抽样，评估统计量的波动。若波动在可接受范围内就停止增加样本
- 逐次增加: 从$N_0 = 10d$开始，增加样本直到代理模型的交叉验证误差达到目标
- 计算预算方法: 根据"（单个分析时间 × 样本数 / 并行数）≤ 允许天数"反推$N$。这是最现实的方法
常见失误与应对
用LHS出错的情况都有哪些？
现场很常见的错误：
| 失误形式 | 原因 | 解决方案 |
|---|---|---|
| "用了LHS结果还是波动大" | 用的是无优化的素LHS，空间填充性不好 | 改用最大最小距离或CD准则优化 |
| "变量之间出现意外相关" | 独立生成置换导致伪相关出现 | 使用Iman-Conover法把相关控制到零 |
| "某些分析出现发散" | 尾部样本落在非物理的极端值范围 | 给输入分布设上下限截断（如±3σ） |
| "代理模型精度不理想" | 空间填充性差 | 用最大最小距离优化的LHS |
| "结果无法复现" | 没有固定随机数种子 | 明确设置种子并记录 |
| "生成样本特别慢" | 优化迭代次数过多 | N ≤ 200用默认参数，N > 1000改ESE法 |
| "追加样本后LHS性质崩坏" | 直接往现有LHS加点，破坏了层化 | 用PLHS（Progressive LHS）从设计初期就考虑扩展 |
分布尾部的样本导致分析发散，这个我之前没想到…
正是LHS的一个陷阱。LHS通过分层确保尾部有样本，这在统计上是好事。但对CAE来说，极端参数组合可能让求解器崩溃。标准做法是给正态分布设$[\mu - 3\sigma, \mu + 3\sigma]$的截断，对数正态分布用0.1%～99.9%分位数截断。
拉丁超立方体抽样（LHS）的软件对比
商用和开源工具对比
能用LHS的工具有哪些？
商用和开源都很丰富。
| 工具名 | 类型 | LHS功能 | 优化LHS | 相关控制 |
|---|---|---|---|---|
| Ansys optiSLang | 商用 | 标配 | 最大最小距离, CD | Iman-Conover |
| DAKOTA (Sandia) | 开源 | 标配 | 最大最小距离 | 秩相关控制 |
| modeFRONTIER | 商用 | 标配 | 最大最小距离, audze-eglais | 支持 |
| UQLab (ETH Zurich) | 学术开源 | 标配 | 最大最小距离, $\phi_p$ | Nataf变换 |
| SciPy (Python) | 开源 | qmc.LatinHypercube | random-cd, Lloyd | 需手动实现 |
| pyDOE2 (Python) | 开源 | lhs() | 最大最小距离, 相关 | 支持 |
| MATLAB | 商用 | lhsdesign() | 最大最小距离 | lhsnorm() |
| JMP | 商用 | Space Filling Design | 最大最小距离 | 支持 |
| LS-OPT | 商用 | 标配 | D-最优LHS | 支持 |
零成本开始的话选哪个？
Python的SciPy的qmc.LatinHypercube现在是最佳选择。只需NumPy/SciPy，有优化选项。要做大规模参数研究的话，DAKOTA（Sandia国家实验室开发）是业界标准，虽然没有GUI但命令行极其强大，完全免费。
拉丁超立方体抽样（LHS）的前沿研究
自适应LHS（Adaptive LHS）
最近的研究中LHS有什么发展？
大趋势是自适应抽样。传统LHS一开始就决定全部样本（"一次性设计"），自适应LHS则根据分析结果动态添加新样本。具体地，在代理模型预测不确定性大的地方重点抽样。
- EGO（全局高效优化）与集成: Kriging预测方差大的点添加LHS样本
- EIVAR（方差改进预期）: 在使方差降低幅度最大的位置添加
- 多保真度LHS: 低精度模型（粗网格）和高精度模型（细网格）使用不同的样本数
序列LHS（SLHS: Sequential Latin Hypercube Sampling）
之前提到的"追加样本破坏LHS性质"的问题，SLHS怎样解决？
SLHS的妙招是，初期设计时就为未来扩展预留好层的空间。Sheikholeslami & Razavi (2017) 的PLHS（渐进式LHS）是代表作。初期用$N_1$个样本的LHS，层划分时就按$N_1 + N_2$的方式设计，这样以后加$N_2$个新样本时，整体仍然保持LHS的分层性质。
高维问题的拓展
变量数很多（几十到几百）的时候，LHS还能用吗？
高维（$d > 20$）时空间填充的优化变困难，而按$10d$法则样本数会很大。应对方法：
- 降维与组合: 用活跃子空间（Active Subspace）或PCA降低维数，再应用LHS
- 两阶段设计: 先用Morris法或Sobol法筛选关键变量，再对这些变量做LHS
- 稀疏网格混合: 低维方向用稀疏网格，高维用LHS配合
- 准蒙特卡罗对比: Sobol或Halton数列在高维下差异度也很低。$d > 50$时有时QMC比LHS更优
高维的时候LHS和QMC，选谁更好？
没有绝对答案。$d \leq 20$且分布信息充分时LHS几乎总是最优。$d$很大或分布信息缺乏时，QMC（Sobol/Halton序列）通常更高效。现代框架（DAKOTA、UQLab）都支持两者，最好的做法是都试试，看结果。
拉丁超立方体抽样（LHS）故障排查
常见问题与解决方案
用LHS时遇到问题该怎么查？能给我个速查表吗？
给你总结一份。
| 现象 | 可能原因 | 处理办法 |
|---|---|---|
| 估计均值与真值偏差大 | 样本数不足或输入分布设置错误 | 增加N为2倍，检查分布参数 |
| 样本散布图呈"斜线" | 变量间产生非预期相关 | 用Iman-Conover法把相关控制到零 |
| 部分CAE分析发散 | 尾部非物理样本值 | 给分布设截断界（如±3σ） |
| 代理模型精度达不到 | 空间填充性差 | 改用最大最小距离或CD优化的LHS |
| 结果不可重复 | 随机数种子未固定 | 明确指定种子并记录 |
| LHS生成很慢 | 优化迭代过多 | N ≤ 200用默认，N > 1000改ESE方法 |
| 追加样本后LHS性质消失 | 无计划地往现有LHS加点 | 用PLHS从开始就规划扩展 |
总结一下，LHS的核心是什么？
"最实用的输入空间均匀覆盖策略"。相比蒙特卡罗总是更高效，实现也简单。在不确定性量化、DOE、代理模型构建等各个场景都是首选。理论卡壳了就回过头来看这篇文章。