简介压缩包内含三个MATLAB脚本mode.m、rank_sort_new.m、hypervolume.m聚焦超体积指标在多目标优化中的应用主要面向进化算法研究者、研究生以及需要量化帕累托前沿质量的开发人员。hypervolume.m 是核心代码负责根据参考点计算解集覆盖体积rank_sort_new.m 实现非支配排序为超体积计算提供有序解集mode.m 用于提取多目标解集的模式或中心参考点辅助确定计算基准三者配合可快速搭建“解集生成—排序—超体积评估”链路。压缩包共3个文件总大小仅6KB代码轻量易读便于在此基础上进行二次扩展。目前已有954人学习或下载说明该工具在MOEA实验、算法对比和论文复现中有一定参考价值。通过这份资源读者能获得超体积指标的标准实现思路、非支配排序与参考点设置的关键算法并可直接运行脚本评估不同解集的收敛性与分布性节省从头编码的时间。1. hypervolume_一个把“解集质量”从玄学变成数字的小工程跑多目标进化算法的人都有过这种经历一代算法跑完弹回来 50 个非支配解看着 Pareto 前沿的形状感觉“好像比上一代好”但要你说出好在哪又只能支支吾吾。我当年也这样直到把 hypervolume_ 这个指标做成一个能直接跑的工程模块才把这种“感觉”换成了稳定输出的小数。hypervolume_ 要做的事很简单给一组解和一个参考点算出一个数这个数越大说明这一代解集离真实前沿越近、分布越均匀。它能解决的是多目标优化里最头疼的评估问题——多组解摆在一起到底哪组更好做算法调参、论文对比实验、线上模型效果回归都绕不开它。2. 先搞懂它在算什么参考点、支配区和一坨超体积2.1 一个三点例子讲透超体积的几何含义假设目标空间是二维两个目标都是最小化三个解分别是 (0.2, 0.9)、(0.4, 0.5)、(0.7, 0.3)参考点取 (1, 1)。每个解都向参考点的方向撑出一个矩形第一个解撑出 x∈[0.2,1]、y∈[0.9,1] 这块区域。三个矩形取并集并集的面积就是这个解集的超体积。手算一遍比背公式更管用x 从 0.2 到 0.4只有第一个矩形覆盖宽度是 0.2y 方向覆盖 [0.9,1]面积贡献 0.02x 从 0.4 到 0.7前两个矩形都覆盖y 方向并集下界是 0.5贡献 0.15x 从 0.7 到 1.0三个矩形都覆盖y 方向下界是 0.3贡献 0.21。总面积 0.38。这段手工推演可以写成几行代码验证顺便帮自己确认对“支配”的理解没跑偏points [(0.2, 0.9), (0.4, 0.5), (0.7, 0.3)] ref (1.0, 1.0) xs sorted(p[0] for p in points) ys_by_x {p[0]: p[1] for p in points} ys [ys_by_x[x] for x in xs] area 0.0 min_y float(inf) prev_x 0.0 for x, y in zip(xs, ys): if min_y ! float(inf): area (x - prev_x) * (ref[1] - min_y) min_y min(min_y, y) prev_x x area (ref[0] - prev_x) * (ref[1] - min_y) print(area) # 0.38这段逻辑的核心是按第一个目标轴排序后扫描过程中维护“当前已遇到点的最小第二目标值”。因为多目标最小化里解 p 能支配的区域是 [p1, r1]×[p2, r2]在 x 轴的任意一个切片上y 方向被覆盖的下界是“所有 x 值小于等于该切片的解里最小的 y”。逐个推进就能把并集面积一点点累出来。2.2 为什么同行宁可算它也不用间距和世代距离做多目标评估时最常见的替代品有三个世代距离GD、反世代距离IGD、间距指标Spacing。GD 度量“求解结果离真实前沿的平均距离”IGD 度量“真实前沿上的点离求解结果的平均距离”间距度量“解在目标空间分布得均不均匀”。它们各自都有明显软肋IGD 必须知道真实 Pareto 前沿实际工程问题里这个前沿根本不存在只能拿一个近似前沿顶上近似前沿本身有偏指标跟着就偏。间距指标只看均匀性一组解挤在角落但点间距均匀它照样给高分。超体积指标不需要真实前沿它只需要一个参考点。另外它有个很漂亮的单调性把解集整体往理想点方向推一步超体积严格变大。这意味着“收敛更好”这件事会直接反应在数字上。加上它对分布均匀性也有惩罚——点都堆在局部的话覆盖面积长不上去。我遇到过实际案例同一个双目标问题我用 IGD 对比两个算法A 算法得分 0.031B 算法得分 0.047按 IGD 越小越好A 胜出但用超体积算A 是 0.52B 是 0.61B 反而赢。后来把两个解集画出来才发现A 只是贴着一小段真实前沿分布得特别近但覆盖范围窄B 虽然离前沿稍远覆盖了整个目标空间。对使用者来说B 的可用性显然更好。这就是我不再用 IGD 做单一结论的原因。2.3 参考点怎么选一个参数决定整个数字的可靠性参考点选不好超体积计算出来的数字再精确也没有意义。参考点的语义是“目标空间里的最差点”必须保证所有解的所有目标值都不比参考点差否则对应维度会出现负贡献。常见做法是取每个目标维度上所有解的最大值再乘一个 1.1 之类的松弛系数避免某个解刚好压在参考点上时出现零面积。三种常用选法如下选法具体操作适用场景注意事项动态最大值用当前解集各目标最大值单代内部比较不同代之间不可比静态松弛值用预估的工程边界值乘系数多代纵向比较必须固定换了就不可比归一化固定点先把目标值映射到 [0,1]参考点固定为 (1,1,...,1)写论文、跑 benchmark归一化时注意异常点多数论文实验里用第三种因为超体积对参考点位置敏感这件事是公认的坑参考点拉得越远靠近理想点的解贡献占比越小超体积对不同解集的区分度会下降。我自己做横向对比时参考点一旦定下来整个实验周期内不会动它并且会记录在实验配置里这是能让别人复现你数字的前提。3. 把 hypervolume_ 写出来从最小实现到通用接口3.1 二维精确计算先剔被支配点再排序扫描二维超体积是最容易写对、也最适合当回归测试基准的实现。工程上不能直接套前文的简单扫描因为真实算法跑出来的种群解往往混着被支配的点直接扫描会重复计算很多区域数字虚高。所以我习惯先做一个非支配筛选再进扫描流程。import numpy as np def is_non_dominated(points): n len(points) keep np.ones(n, dtypebool) for i in range(n): if not keep[i]: continue for j in range(n): if i j or not keep[j]: continue # 如果 j 在每个目标上都不劣于 i且至少一个目标严格优于 i则 i 被支配 if np.all(points[j] points[i]) and np.any(points[j] points[i]): keep[i] False break return points[keep] def hv_2d(points, ref): points is_non_dominated(np.asarray(points, dtypefloat)) order np.argsort(points[:, 0]) points points[order] r1, r2 ref hv 0.0 min_y float(inf) prev_x 0.0 for x, y in points: if min_y ! float(inf): hv (x - prev_x) * (r2 - min_y) min_y min(min_y, y) prev_x x hv (r1 - prev_x) * (r2 - min_y) return hv这里的is_non_dominated用的是两两比较O(n²)对一两千个点没问题几万个点就要换成排序或 KD 树思路后面会提。hv_2d里的扫描算法来自 HSOHypervolume by Slicing Objectives思想的二维特例按第一个目标轴排序从左往右扫维护前缀最小第二目标值。这个“前缀最小”是二维超体积最容易写反的地方我见过好几个开源实现把min写成max结果算出来虚高一大截。3.2 三维精确计算把第一个维度切掉降维打击三维精确计算不用发明新算法用一个很直接的降维思路先把第一目标轴切成若干段每段上用二维超体积函数去算剩余两个维度上的投影覆盖面积再乘上段的厚度。这里的“段”取决于第一目标值排序后相邻解之间的间隙。def hv_3d(points, ref): points is_non_dominated(np.asarray(points, dtypefloat)) order np.argsort(points[:, 0]) points points[order] r1, r2, r3 ref hv 0.0 collected [] prev_x 0.0 for p in points: if len(collected) 0: proj np.array([[q[1], q[2]] for q in collected]) hv (p[0] - prev_x) * hv_2d(proj, (r2, r3)) collected.append(p) prev_x p[0] proj np.array([[q[1], q[2]] for q in collected]) hv (r1 - prev_x) * hv_2d(proj, (r2, r3)) return hv判断第一段为什么不需要算第一段是从 0 到第一个点的第一目标值这段里没有任何解在 x 轴方向上覆盖它投影集合为空hv_2d对空集返回 0所以直接跳过。最后一段从最后一个点的第一目标值到参考点的第一目标值所有点都在覆盖所以要用完整集合。中间每段的投影集合是随着扫描不断累加的因为只有当某个解的第一目标值已经“进入”当前切片轴范围它才可能参与覆盖。3.3 高维近似蒙特卡洛采样当高维测谎仪超过三个目标以后精确计算的代价上升得非常快工程上更常用的是蒙特卡洛采样近似。思路很简单在参考点框定的超立方体里随机撒点统计有多少点被解集支配比例乘以立方体体积就是超体积的估计值。def hv_mc(points, ref, n_samples50000, seed42): rng np.random.default_rng(seed) points np.asarray(points, dtypefloat) ref np.asarray(ref, dtypefloat) samples rng.uniform(0.0, ref, size(n_samples, len(ref))) # points[:, None, :] samples 得到 (n_points, n_samples, n_dims) dominated np.all(points[:, None, :] samples, axis2).any(axis0) ratio dominated.mean() return ratio * np.prod(ref)np.all(points[:, None, :] samples, axis2)这一步的含义是对每个采样点判断是否存在一个解在每个目标上都比它小如果是说明这个采样点落在了解集的支配区域内。没有任何一个解支配它则它不在超体积区域内。np.prod(ref)是参考点盒子的总体积。维度升高以后这个近似的方差会变大因为支配区域在超立方体里的占比会越来越低大量采样点落在支配区外。所以要么加大采样量要么做分层采样。固定seed是必须的否则两次跑同一个数据得到不同数字实验记录里你根本解释不清是算法变了还是采样波动。3.4 模块接口怎么设计一个能长期维护的 hypervolume_功能写完后接口设计决定这个工具能不能陪你把整个实验周期走完。我习惯用函数式而不是类式接口因为多目标实验的每个环节拿到的都是 numpy 数组函数式接口改造成本最低。模块名就叫hypervolume_目录里放一个核心文件加两个示例脚本。def compute(points, ref, methodauto, seed42): points np.asarray(points, dtypefloat) ref np.asarray(ref, dtypefloat) if method auto: if points.shape[1] 2: return hv_2d(points, ref) elif points.shape[1] 3: return hv_3d(points, ref) else: return hv_mc(points, ref, seedseed) if method 2d: return hv_2d(points, ref) if method 3d: return hv_3d(points, ref) if method mc: return hv_mc(points, ref, seedseed) raise ValueError(funknown method: {method})methodauto根据目标维度自动派发日常用起来最省事。有一点值得强调传入的points一定是“每行一个解、每列一个目标”的二维数组这是整个模块最容易出错的地方。我把行列搞反过一次当时所有解的第一目标值被当成不同维度的目标计算结果比实际大了好几倍之后我在接口处加了维度校验才彻底杜绝这个问题。4. 算不动的时候怎么办性能瓶颈和工程提速三板斧4.1 复杂度从哪里来从几千个解到几十个维度的差距超体积精确计算的复杂度主要由两个因素决定解的数量 n 和目标维度 k。二维排序扫描是 O(n log n)三维降维切片的复杂度已经是 O(n²) 量级因为每个切段都要对投影集合做一次二维计算而每个投影集合的大小在递增。到了四维五维精确算法涉及的递归划分和支配关系检索会进一步膨胀工程实测里超过五个目标还做精确计算数据量稍大就能把一次实验拖到分钟级甚至更久。理论界对这个问题的结论也很激进精确计算超体积在高维场景下本质上是把高维空间切成无数小块再逐块判断维度涨上去后这个划分的规模是指数级增长的。所以工程上有个默认的分界线三维以内用精确计算四到六维要看数据量决定六维以上直接换近似方法不要硬撑。4.2 提速第一板斧先把被支配点清掉被支配点对超体积的数值没有任何贡献但会增加排序、比较、递归的计算量。尤其是高维场景支配关系更稀疏混入的无效点比例很高。我之前跑一个五目标实验种群大小 200非支配筛选后只剩 37 个点后续所有计算量直接砍掉八成。def is_non_dominated_fast(points): n, k points.shape order np.argsort(points[:, 0]) sorted_pts points[order] keep np.ones(n, dtypebool) min_rest np.full(k - 1, np.inf) for i in range(n): # 当前点其余维度是否比已知前缀都差 if np.any(sorted_pts[i, 1:] min_rest): keep[order[i]] False min_rest np.minimum(min_rest, sorted_pts[i, 1:]) return points[keep]这个版本的思路是按第一维排序后只扫一遍如果当前点的所有剩余维度都已经“比某个前缀点差或相等”那它一定是被支配的。它牺牲了一点准确性来换速度严格场景还是用两两比较那个版本。加速的本质是减少进入精确计算或采样流程的解数量这一刀砍下去最直接。4.3 提速第二板斧增量维护环境选择里的边际贡献在基于超体积的进化算法里比如 SMS-EMOA每一代要淘汰一个解标准做法是删掉“超体积贡献最小”的那个。朴素实现是每次移除一个点就重算一遍整体超体积时间复杂度立刻翻倍。工程上更聪明的做法是只算边际贡献每次删一个候选解计算它被移走后超体积下降了多少然后删掉下降最小的那个。def contribution(points, idx, ref, method3d): mask np.ones(len(points), dtypebool) mask[idx] False return compute(points, ref, methodmethod) - compute(points[mask], ref, methodmethod) def select_worst(points, ref, method3d): candidates [contribution(points, i, ref, method) for i in range(len(points))] return int(np.argmin(candidates))这种方式仍然是 O(n) 次重算但配合第一板斧的非支配筛选实际运行时间比完全不裁剪快很多。再往后还有用 R2 指标替代超体积做环境选择的做法那个改的是候选解对权重向量的聚合值计算更快但不是同一个指标别混用概念。4.4 参考点归一化一个解决精度问题的隐藏加速归一化到 [0,1] 不只是为了统一参考点它还顺带解决了浮点精度和采样效率问题。目标值如果分布在 1e5 这个量级参考点取 1.2e5蒙特卡洛采样在 [0, 1.2e5] 的超立方体里撒点支配区域占比会小到离谱采样估计的方差也随之变大。归一化以后所有目标值都压在单位立方体里数值稳定性好得多。def normalize(points, ref): pts np.asarray(points, dtypefloat) return pts / ref除法完成两点一是把参考点变成全 1 向量二是把解的目标值按比例映射到 [0,1] 区间。注意参考点必须严格大于所有目标值否则归一化后某些维度的值会超过 1这时对应解的支配区域超出了参考点盒子超体积计算就失去意义了。这个坑在下一章细讲。5. 超体积计算避坑实录五个我踩过的翻车现场5.1 参考点不固定两次结果根本不可比现象同一个算法跑两次第一次超体积 0.61第二次 0.58看起来第二次退步了。仔细查才发现两次实验的参考点分别来自各代解集的最大值第二次迭代早期出现了一个异常大目标值的解把参考点拉远了所有解的超体积整体变小。原因动态参考点让每次计算的“盒子”大小不同数字的绝对大小没有可比性。解决实验启动前就把参考点定死用上一章说的静态松弛值或归一化固定点。参考点一旦定下写进实验配置不许中途改。做参数敏感性分析时可以专门扫参考点取值但那是另一件事不能和正式实验混在一起。5.2 把最大化问题直接喂进去体积算出负值现象一个收益最大化的双目标问题解的目标值越大越好我按最小化流程直接算结果超体积是负数。原因超体积的全部几何推导都建立在“目标越小越好”的假设上。最大化问题里解往右上角移动才叫收敛支配方向完全反了。参考点取最大值时解比参考点还大矩形区域变成负宽度。解决最大化问题先对目标值取负号变成最小化问题再计算。转换后参考点取负值空间里的合适上界。更推荐的做法是在优化阶段就把问题建模成最小化避免评估时再做一层变换引入浮点误差。5.3 被支配点没过滤面积虚高得离谱现象一组包含 200 个解的数据直接算二维超体积是 0.47过滤掉 143 个被支配点后再算只有 0.22。被支配点对真实前沿没有任何贡献数字却大了一倍。原因被支配点也会撑出支配矩形它们的矩形落在已有点的支配区域内部按排序扫描的累加逻辑这些点的 y 值如果比前缀最小值小就会错误地拉大覆盖宽度。解决计算结果前先执行非支配筛选。我现在的流程是筛选、归一化、排序三步固定链条任何一步都不能跳。高维场景下筛选本身也有成本但那笔账值得付。5.4 参考点设得太大浮点精度把面积冲没了现象一个目标的范围是 [0, 100]另一个是 [0, 1e6]参考点直接取各自最大值。计算三维超体积时数值一直不稳定同一份数据两次运行差 4%。原因目标维度之间的量级差了几个数量级。小数值维度的覆盖宽度只有 0.1 量级乘上大数值维度后浮点数在小数的低位上出现舍入误差。蒙特卡洛采样时采样点在短维度上几乎总落在支配区域内长维度上的覆盖决策被浮点误差干扰。解决先按目标维度做量级归一化或者把所有目标值缩放到 [0,1] 再计算。归一化的缩放系数要记录下来报告超体积时写清是归一化后的数值。同一实验内的多个解集必须用同一套缩放系数。5.5 蒙特卡洛采样没固定种子A/B 结论翻车现象两个算法各跑十次算法 A 的超体积均值 0.523算法 B 是 0.518看起来 A 赢了。但每次单独对比时B 有三次比 A 高且这三次恰好发生在不同种子下。原因蒙特卡洛采样本身就是随机过程不固定种子每次采样的点集都不一样估计值的波动可比算法差异还大。样本量 5 万在高维下远远不够方差依然显著。解决固定随机种子并加大采样量到 1e5 以上。更严谨的做法是同一个数据重复采样 10 次报告均值和标准差。我还给自己加了一条纪律用近似方法做出的排序结论一定要用三维精确计算或更大的采样量复核一遍避免把采样噪声当成算法改进。6. 把 hypervolume_ 接进你的进化算法一个实用小技巧6.1 用 HV 曲线替代“肉眼判断收敛”算法调参时我除了看最终一代的非支配前沿还会把每一代的 hypervolume_ 值拉成一条曲线。曲线变平说明收敛进入平台期曲线还在爬说明优化仍在有效推进。这个曲线也可以直接作为早停条件连续 20 代超体积提升都小于 1e-4停止迭代。history [] for gen in range(max_gen): offspring algorithm.ask() fitness evaluate(offspring) algorithm.tell(fitness) points algorithm.get_non_dominated() hv compute(points, ref, methodauto) history.append(hv) if len(history) 20 and abs(history[-1] - history[-21]) 1e-4: break这个技巧的关键点是methodauto在低维是精确值、高维是采样近似前后要统一曲线的波动才会小到能看出趋势。注意参考点必须和正式实验的参考点一致否则早停条件会误触发。6.2 我现在的验证习惯写了这个模块之后我养成了三个习惯。第一任何新数据集接入前先用二维精确计算跑一组已知答案的样例确认算法实现没回归。第二每次实验记录参考点和归一化系数报告结果时附上这两个参数。第三高维结果永远附带不确定性说明采样方法给出的数字从来不是“精确值”。超体积这个指标的边界我也越来越清楚它能告诉你解集整体质量好不好但不会告诉你解集在哪个局部区域特别优秀。我吃过一次亏高维问题上只依赖近似超体积做早停结果算法停在了一个局部覆盖很烂但整体体积数字还行的状态。后来改成超体积曲线配合各目标箱线图一起看问题才解决。做多目标优化指标是导航仪不是目的地。希望这些踩坑记录能帮你少走几段弯路把 hypervolume_ 真正用起来。本文还有配套的精品资源点击获取
