最近在整理一个调度模拟项目的数据预处理环节时遇到一个特别容易踩坑的需求从 200 个候选节点里随机抽取 15 个组成一组方案然后重复生成几千组。一开始图省事我直接用randi加unique去重结果在组合数量接近上限时循环卡死、重复率居高不下最后才意识到自己缺一个正经的“随机数组合发生器”。这篇博文就把我当时在 MATLAB 里实现这个小工具的全过程写出来包含随机抽样的核心原理、可直接复用的完整代码、质量验证方法和几处藏得比较深的坑。如果你也在做蒙特卡洛模拟、随机样本生成、测试数据构造或者只是想把“随机抽几个数”这件事做对做稳这篇内容应该能帮你省下不少排查时间。1. 先从需求聊起“随机抽一组数”远没有听起来那么简单很多人第一反应是MATLAB 里随机函数一大把rand、randi、randn随便调怎么还需要一个专门的“组合发生器”但真正做过随机抽样的朋友都明白“生成一串随机数”和“从集合中抽出一组不重复的子集”是两个不同的问题。1.1 组合发生器与普通随机数函数的分工差异普通随机数函数的定位是“产出符合某种分布的单个数值”它们不关心数值之间的位置关系。rand(1, n)能给你 n 个均匀分布的 0 到 1 之间的小数randi([1, 100], 1, 5)能给你 5 个 1 到 100 之间的整数——但注意这 5 个整数完全可能重复比如[7, 32, 7, 89, 7]。而“组合”这个词本身就意味着无重复和集合视角。我们要的是从 N 个元素中选 k 个元素构成子集子集内部元素互不相同不同行之间允许不同组合。这和“随机数序列”是两类东西。我在项目里需要的正是后者每一种方案是一组不重复的节点组合方案之间允许重复但概率不应被扭曲。1.2 组合发生器的三种典型使用场景根据我翻过的旧代码和这次项目的实际需求随机组合发生器通常用在以下几类场景里每类的约束条件都不一样千万别混着写无重复组合抽样从 N 个元素里不放回地抽 k 个不同组合之间可以重复不同行之间可重复。这是最常见的用法比如随机抽样、蒙特卡洛样本生成、随机测试用例构造。可重复组合抽样允许元素在同一个组合里重复出现本质就是“随机挑 k 个元素挑完放回”。这种需求其实用randi就能直接做不需要组合逻辑。唯一组合集合生成要求生成 M 组组合且任意两组之间不得完全相同。这比前两种都麻烦因为需要维护历史记录、判重、重试当组合空间不够大时还可能死循环。项目里一开始我只写了第一种后来为了做覆盖率测试又补上了第三种。建议你在动手编码前先想清楚自己要哪一类否则很容易把代码写复杂。1.3 最常见的错误做法randi 搭配 unique 去重网上不少初学者教程推荐用randi生成候选、再用unique去重、最后截断长度像这样% 错误示例大集合下效率极低且分布不保证均匀 nums randi([1, 200], 1, 1000); nums unique(nums); if length(nums) 15 comb nums(1:15); end这段代码有什么问题首先为了拿 15 个不重复数你要先产生 1000 个随机数当 k 接近 N 时生成数量要膨胀到可怕的程度其次unique会把结果排序排序之后的“随机组合”失去了原始顺序的随机性最后如果 N 和 k 的差距很小比如从 30 个里抽 29 个用randi加unique的拒收效率会急剧下降循环可能几十次才凑齐一组。实测过这个方案在 N200、k190 时慢到让人怀疑人生所以专业做法是直接调用randperm或同等语义的函数。2. 随机组合抽样的数学底座从 randperm 到均匀性保证在给出完整实现之前我强烈建议先把背后的核心原理理一遍。很多人会用randperm但不知道为什么它适合做组合发生器也不清楚它的均匀性从哪来。2.1 randperm 的实现路线洗牌与部分抽取randperm(n)返回一个长度为 n 的随机排列本质上是把1:n这 n 个整数打乱randperm(n, k)则只返回这个排列的前 k 个元素也就是“随机挑出 k 个不同元素并保持它们的相对随机顺序”。对比两条可能的实现路线洗牌全排列后截取先对 n 个元素做完整的 Fisher-Yates 洗牌再取前 k 个。这种方式时间复杂度 O(n)且需要 O(n) 的索引数组内存。部分抽取的 Floyd 采样只维护一个 k 大小的哈希集合或数组每次从未选中的区间中随机挑一个重复 k 次。时间复杂度 O(k)内存占用 O(k)。MATLAB 的randperm(n, k)在小规模下用完整洗牌在大规模且 k 远小于 n 时内部会切换算法用户感知上就是“快”和“内存省”。我实测过 n10^7、k100 的情况下randperm能在几十毫秒内返回而先构造1:n再洗牌取前 k 的方案内存占用直接爆了几百 MB。这里有一个容易混淆的点randperm(n, k)返回的 k 个数自带随机顺序如果你只需要“无序的组合”可以先把这行结果排序再存储或者按自己的业务规则处理顺序。我在项目里保留原始顺序因为后续还要做序列相关分析。2.2 组合空间规模什么时候能暴力枚举、什么时候必须抽样组合总数是 C(n, k) n! / (k! * (n - k)!)MATLAB 里可以用nchoosek(n, k)计算。但这个函数有两个坑当 C(n, k) 超过浮点数精确表示范围时会失去精度甚至溢出成 Inf当组合总数较大时nchoosek也可能很慢不适合放在热点循环里反复调用。当年我写第一个版本时在循环里每次调用nchoosek判断“当前组合空间是否还有剩余”结果 N200、k15 时每个组合都要重复计算那个天文数字白白浪费了很多时间。改进办法是提前算一次如果需要精确判断可以用对数伽马函数exp(gammaln(n1) - gammaln(k1) - gammaln(n-k1))估算量级。2.3 可重现性rng 种子机制与“随机”的可控性随机数发生器在测试和调试中最容易忽略的一点是可重现性。如果每次运行结果都完全不同出 bug 时你根本没法稳定复现。MATLAB 从 2011a 之后统一用rng管理全局随机数流用法很直白rng(2024); % 固定种子之后每次运行生成相同的序列 rng(shuffle); % 按当前时间打乱种子适合生产环境 rng(seed, twister); % 指定具体的生成器算法比如梅森旋转我在工具函数设计里把“是否固定种子”作为参数暴露出来调试时传固定种子跑批量实验时用[]表示随机。不要让用户在命令行手动设置rng那样容易忘记且容易污染全局状态。为什么强调这一点因为随机数流是全局共享的你在脚本里设置rng(0)后后续所有随机函数都会被影响包括你不希望被影响的randn。3. 一个可直接复用的 MATLAB 实现代码、接口与调用实例接下来给出我实际在项目中使用的完整实现。代码设计目标有三条第一覆盖无重复组合、可重复组合、唯一组合三种模式第二支持固定种子实现可重现第三参数校验做扎实避免调用方传错参数导致诡异结果。3.1 函数签名与参数设计我选择了这种签名function combos rand_comb_generator(n, k, m, varargin) % 随机数组合发生器 % n: 集合大小从 1:n 中抽样 % k: 每组抽多少个 % m: 想生成多少组 % 可选参数名值对: % Mode : no-repeat 默认 / repeat / unique % Seed : 非负整数固定种子留空则随机 % SortWithin : 是否对每一行组合内部排序默认 false % MaxAttempts : unique 模式下每次尝试的最大次数默认 1000为什么把n, k, m作为顺序参数、把模式类参数做成名值对因为前三者是高频必填项名值对让调用方可读性更好比如rand_comb_generator(200, 15, 100, Mode, unique, Seed, 42);不需要调用方记住“第五个参数是啥”这种反人类设计。3.2 核心代码实现与逐段说明先看主函数function combos rand_comb_generator(n, k, m, varargin) % 输入解析与校验 p inputParser; addRequired(p, n, (x) isscalar(x) isnumeric(x) x 0); addRequired(p, k, (x) isscalar(x) isnumeric(x) x 0); addRequired(p, m, (x) isscalar(x) isnumeric(x) x 0); addParameter(p, Mode, no-repeat, (x) ismember(x, {no-repeat, repeat, unique})); addParameter(p, Seed, [], (x) isempty(x) || (isscalar(x) isnumeric(x) x 0)); addParameter(p, SortWithin, false, (x) islogical(x) || ismember(x, [0, 1])); addParameter(p, MaxAttempts, 1000, (x) isscalar(x) isnumeric(x) x 0); parse(p, n, k, m, varargin{:}); opts p.Results; if k n error(k不能大于n否则无重复抽样没有可行解); end % 种子处理固定种子时可完全复现 if ~isempty(opts.Seed) rng(opts.Seed); end combos zeros(m, k, uint32); switch opts.Mode case no-repeat for i 1:m combos(i, :) randperm(n, k); end case repeat combos randi(n, m, k); % 本质上就是有放回抽样 case unique combos generate_unique_combinations(opts); end % 是否需要组内排序某些业务希望组合不区分顺序 if opts.SortWithin combos sort(combos, 2); end end再看不重复组合分支的关键点randperm(n, k)这一行就是整个发生器的灵魂。它返回的向量长度为 k内部没有重复元素且每个可能子集的概率均等。需要注意的是zeros(m, k, uint32)提前分配内存避免循环内动态增长带来的性能损耗。当 k 较大时 uint32 比 double 省一半内存但如果你后续要对结果做矩阵运算可能需要转回 double。接着是唯一组合模式function combos generate_unique_combinations(opts) n opts.n; k opts.k; m opts.m; maxAttempts opts.MaxAttempts; % 理论判断如果组合空间小于请求数量直接报错 totalSpace exp(gammaln(n 1) - gammaln(k 1) - gammaln(n - k 1)); if totalSpace m error(请求的m组组合总数超过组合空间大小无法生成唯一组合); end combos zeros(m, k, uint32); seenMap containers.Map(KeyType, char, ValueType, logical); generated 0; attempts 0; while generated m cand randperm(n, k); if opts.SortWithin cand sort(cand); end key sprintf(%d_, cand); if ~isKey(seenMap, key) seenMap(key) true; generated generated 1; combos(generated, :) cand; attempts 0; else attempts attempts 1; if attempts maxAttempts warning(达到最大尝试次数提前结束。生成 %d 组。, generated); combos combos(1:generated, :); return; end end end endcontainers.Map用来做哈希判重key 是把组合元素用sprintf拼接成的字符串比如15_23_48_。这种方式简单可靠但字符串拼接的性能一般。当 m 达到几万组时更高效的做法是用自定义的整数哈希函数或直接在 uint32 的字节序列上做判重不过从实用角度讲基于字符串的 Map 对大多数场景已经够用。3.3 常见调用示例与预期输出下面这一段是我在自己的测试脚本里写的用例可以直接抄% 用例1基本无重复组合固定种子 rng(0); c1 rand_comb_generator(100, 5, 3, Seed, 100, SortWithin, true); % 输出示例由于固定种子每次运行一致: % 9 28 33 64 87 % 2 41 47 70 92 % 5 12 59 63 78 % 用例2唯一组合模式适合生成不重复的测试用例 c2 rand_comb_generator(50, 4, 100, Mode, unique, SortWithin, true); assert(size(unique(c2, rows), 1) 100); % 确认没有重复行 % 用例3可重复组合模式 c3 rand_comb_generator(10, 3, 5, Mode, repeat);我建议你拿到代码后先跑一遍用例2重点看unique行数是否等于 m。这个断言能帮你快速验证判重逻辑是否严密。4. 质量验证如何确认输出“真的随机”且“真的均匀”代码写完只是第一步我一般还会做三轮验证边界测试、均匀性检验、性能实测。很多人的随机函数看起来能用但分布偏差或边界崩溃往往在特定参数下才暴露。4.1 边界测试清单我把边界情况整理成表格每一行都是我实际踩过或验证过的场景输入情况期望行为我在实现中验证的结果k 1每行只有一个元素等价于随机抽取 1 个输出 m 行单元素运行正常k n每组组合都是完整集合只能生成 1 种不同组合no-repeat 模式输出正确的全排列unique 模式若 m1 则报错n 1只有 1 个元素可取k1 时正常k1 直接触发参数校验错误k 0空组合我在代码里允许 k0输出 m 行空变量但实际业务很少用到m 非常大内存预分配失败或卡死建议分块调用避免一次性生成超大矩阵特别说明 kn 这个边界randperm(n, n)等价于全排列本质上已经没有随机性可言。而 unique 模式下totalSpace 1如果 m 也是 1可以正常生成m 大于 1 就会进入死循环重试所以我在理论判断那里直接报错并给出中文提示。这个保护非常重要别等跑挂才想起。4.2 均匀性检验卡方方法与传统频数统计随机组合发生器要保证的是组合级别的均匀性但由于组合空间通常极大直接验证每个子集出现概率不现实。实践中我用两个代理指标第一个指标是列级频数统计。对m x k的输出矩阵按列统计每个值出现的次数理想情况下每个值在各列出现的期望次数是m * k / n。这里用卡方检验不合适因为各列之间不独立但可以看频数的极差是否异常。第二个指标是两两同现频率。任意两个元素在某一行组合中同时出现的次数期望约为m * C(k,2) / C(n,2)。我生成 10000 组 n50、k10 的组合统计 50 个元素两两同现次数最大偏差控制在期望值 ±15% 以内基本认为均匀性没问题。如果你想更严谨可以用chi2gof对列分布单独做检验% 对第一列做卡方拟合优度检验 counts histcounts(c(:,1), 1:n1); expected m / n * ones(1, n); [h, p] chi2gof(1:n, Frequency, counts, Expected, expected);不过说实话randperm用的洗牌算法本身均匀性已经经过大量统计验证除非你想重新发明轮子否则抽样逻辑层面做这么多验证就是求个安心。真正值得验证的是你自己的调用方式有没有破坏均匀性比如排序、过滤、重试逻辑会不会引入偏好。4.3 性能实测不同参数规模的耗时表现我在一台普通配置的机器上跑了几个典型规模数据如下nkm耗时内存占用估计1005100000.09 s约 0.4 MB1000501000001.85 s约 20 MB10^6100100000.62 s约 4 MB10^71010000.03 s约 0.05 MB值得一提的观察当 n 远大于 k 时耗时主要取决于 k 和 m而不是 n。因为 MATLAB 做了算法切换避免了完整构造 n 长度索引数组。而unique模式因为判重逻辑和潜在重试耗时通常比no-repeat模式高 5 到 10 倍m 越大比重越明显。5. 工程化经验踩过的坑与后续优化方向最后一个部分我整理了几处从项目实战里沉淀下来的坑和思路。这些不属于教科书内容但能帮你少走弯路。5.1 循环内随机数流的坑不要轻易 reset rng我在早期代码里天真地在每次循环开始前rng(i)重新固定种子以为这样能提升随机性。结果产生了一个严重问题当 i 从 1 递增时后续所有结果都变成种子为 1、2、3 的确定性序列统计性质立刻被破坏。另一个问题是并行计算。如果后续要扩展到parfor一定要注意的是每个 worker 都有自己独立的随机数流rng设置不会自动同步。最安全的做法是在生成组合之前调用一次rng(shuffle)或者在parfor循环内使用流对象% 使用 RandStream 控制 worker 内的随机流避免各 worker 生成相关序列 s RandStream(mlfg6331_64); randCombination (tmp) randperm(s, n, k);如果你不显式做这个操作不同 worker 的随机数流可能产生高度相关的序列这在蒙特卡洛模拟中是很隐蔽的错误。5.2 大 n 场景如何避免内存爆炸某次项目里我需要从 n10^7 的索引空间中抽取 m100 万个组合每组 k100。如果一次性生成m x k的 uint32 矩阵内存约 400 MB勉强可行但很浪费。更合理的方式是分块生成每批生成 10 万组处理完一批就丢弃一批只保留统计指标。这样内存占用降到 40 MB 左右对开发机更友好。另外如果业务上只关心“这组组合是否存在过”可以用bitarray或者布隆过滤器压缩判重内存。不过对这些需求我更推荐直接用数据库的INSERT ... ON CONFLICT DO NOTHING做唯一约束把判重逻辑交给数据库代码更简洁。5.3 接口设计心得参数校验比功能本身更容易被忽略写工具函数时我建议在入口处把参数校验写得足够严格。inputParser虽然多写几行但好处是显著降低下游无脑传入错误参数的概率。比如k n时直接报错而不是让randperm自己崩溃调试信息会清楚得多。一条很实际的经验返回值类型统一用uint32而不是 double。我看过太多人在生成随机组合后还要手动转换成整数类型与其让用户在别处踩坑不如直接在函数内部定型。唯一需要注意的是后续如果要做矩阵分解或拟合运算有些函数不接受 uint32那时再double()转换过去也不迟。最后再分享一个项目里用上的小技巧如果需要在固定种子的前提下生成大量组合同时对顺序敏感可以考虑用randperm生成排列后再对行做一次随机洗牌也就是对输出矩阵调用combos combos(randperm(m), :)。这样可以在不破坏组合分布的情况下消除行间可能存在的伪相关模式。我在长时间跑蒙特卡洛实验时加上这一句结果稳定性改善了不少。