Reference

DataFrame(comparison::SplitComparison) -> DataFrame

One row per splitter: method, kernel, ratio, subset sizes, convergence report, and the discrepancy score (energy_distance or mmd, lower is better).

AbstractSplitter

Supertype of every splitting method usable with datasplit. Concrete splitters expose kernel and ratio fields, which compare/DataFrame(::SplitComparison) read.

DiscrepancyEstimator

Supertype of estimators selected by the estimator keyword: Exact, Subsample, RandomSlices, RandomFeatures.

EnergyKernel()

The kernel k(x, y) = −‖x − y‖, whose maximum mean discrepancy between two samples is the energy distance. Support points under this kernel are those of Mak & Joseph (2018), optimized by their closed-form majorization–minimization update.

Exact()

Evaluate every pairwise term (block-wise, threaded; O(nm) kernel evaluations).

FourierFeatureMap(kernel::GaussianKernel{Float64}, p, D, rng)

Random Fourier features (Rahimi & Recht 2007) for a p-dimensional Gaussian kernel: z(x) = √(2/D) cos(Wx + b) with W ~ N(0, σ⁻² I) and b ~ U[0, 2π], so that E[z(x)ᵀz(y)] = kernelvalue(kernel, x, y). Drawn once per call from rng; callable on a column vector.

GaussianKernel(bandwidth = :median)

The Gaussian (RBF) kernel k(u, v) = exp(−‖u − v‖² / (2σ²)). Support points under this kernel minimize the squared maximum mean discrepancy (MMD², Gretton et al. 2012) by projected gradient descent with Armijo backtracking.

bandwidth is either a positive number σ or :median, in which case σ is set at fit time to the median pairwise Euclidean distance of (a sample of) the standardized data. The resolved value is stored in result.method.kernel. :median fails with an ArgumentError when at least half of all row pairs coincide (e.g. a single binary categorical column): pass a numeric bandwidth in that case. Choose σ on the scale of the standardized data (:median does this), since a bandwidth far below the row spacing makes the objective flat and the optimizer stops at the initial sample.

HerdingSplitter(; kernel = GaussianKernel(), ratio = 0.2,
                  n_threads = Threads.nthreads(), rng = Random.default_rng())

Split by greedy kernel herding (Chen, Welling & Smola 2010): the smaller subset is chosen row by row to minimize the MMD² (energy distance for EnergyKernel) to the whole data. The data term is always computed exactly (threaded, O(N²)); see herd for why an approximate data term is not offered. Deterministic given the data and a numeric kernel; rng only feeds a :median bandwidth. See herd for the rule and cost.

Examples

data = randn(MersenneTwister(1), 200, 3)
result = datasplit(HerdingSplitter(rng = MersenneTwister(2)), data)

KernelThinningSplitter(; kernel = EnergyKernel(), ratio = 0.2, delta = 0.5,
                         compress = :auto, n_threads = Threads.nthreads(),
                         rng = Random.default_rng())

Split by generalized kernel thinning with the target kernel (Dwivedi & Mackey 2022; kernel halving from Dwivedi & Mackey 2024): the smaller side is chosen by kernel_thinning, so it minimizes the MMD² (energy distance for EnergyKernel) to the data without continuous optimization or a nearest-neighbor step, with the papers' high-probability MMD guarantee of order √(log n / n) for the KT-SPLIT candidates and a KT-SWAP result never worse than a uniform random subset. Cost is O(N²) kernel evaluations like HerdingSplitter, unless compress selects Compress++ for near-linear time.

• kernel: EnergyKernel() (default) or GaussianKernel(σ); a :median bandwidth is resolved at datasplit time and stored in result.method.
• ratio: fraction of rows assigned to the test set, in (0, 1); above one half the selection is the complement of a kernel-thinning selection of the other side (see kernel_thinning).
• delta: the failure probability δ of the kernel-thinning guarantees: the papers' δ_i = δ/L, applied as δ_i/m at every halving step (the experiments use δ = 0.5).
• compress: :auto (default) runs Compress++ when n ≪ N makes it cheaper than plain kernel thinning and the target is the data itself. It does not fire at the default 20% ratio, so datasplit with the default splitter is unchanged; it can fire below roughly a 10% ratio for N ≥ 10⁴, where n is small enough for Compress++ to be cheaper. Pass compress = :never to keep the plain path there; :always forces Compress++.
• rng: the input shuffle, the halving coin flips, and the baseline draw.

SplitResult.converged is always true; iterations is the number of KT-SWAP replacements. weights and reference act on the KT-SWAP objective only (see kernel_thinning for the differences from the paper).

Examples

data = randn(MersenneTwister(1), 1_000, 5)
result = datasplit(KernelThinningSplitter(rng = MersenneTwister(2)), data)

MatrixSearch(Xt)

Nearest/k-nearest column search over Xt (p × m, one point per column) by explicit @inbounds @simd distance loops rather than a NearestNeighbors tree. Xt is a plain Matrix{Float64}, so the search compiles once no matter how large p is, unlike KDTree/BruteTree, which specialize on SVector{p, Float64} and become impractical to compile above a few thousand columns (issue #72).

Preprocessor

Fitted preprocessing: per-column encoding rules, which encoded columns are kept (those not constant on both the fit set and, when given, the extra set), and the mean and scale of every kept column. Built by fit_preprocessor, applied by apply_preprocessor. Internal.

RandomFeatures(D)

Estimate a Gaussian-kernel quantity with D random Fourier features (Rahimi & Recht 2007) drawn with rng: unbiased, O((n + m) D p). Defined for GaussianKernel only. With sample weights the feature means are weighted means.

RandomSlices(k)

Estimate an energy distance from k random one-dimensional projections (rng), each evaluated exactly by sorting: unbiased, O(k (n + m) log(n + m)). Defined for the energy kernel only. See the Methods page for the identity E_θ |⟨θ, u⟩| = κ_p ‖u‖ behind it. With sample weights the per-direction one-dimensional energy distance is weighted (prefix sums of the sorted weights).

SplitComparison

Result of compare: the fitted splitters (with any :median bandwidth resolved from the data), their results, their splitquality values, and the kernel they were scored under, index-aligned.

SplitKernel

Abstract supertype of kernels usable with SupportPointSplitter and HerdingSplitter. Each kernel defines the discrepancy between the selected subset and the data distribution that the splitter minimizes.

SplitResult

Outcome of datasplit: index partition plus an honest report of the optimizer's convergence. selected is the side (:test or :train) that holds the rows the splitter chose; the other side is the complement.

Subsample(m, repeats = 8)

Average the exact statistic over repeats random size-m row subsets drawn with rng. Carries a positive bias of order 1/m (the zero diagonal of the V-statistic); meant for comparing splits, not as an absolute value.

SupportPointSplitter(; kernel = EnergyKernel(), ratio = 0.2, kappa = nothing,
                      max_iterations = 500, tolerance = 1e-10,
                      n_threads = Threads.nthreads(),
                      rng = Random.default_rng(), verbose = false)

Configuration for optimal data splitting via support points (Joseph & Vakayil 2022): the smaller side is computed as a set of support points and mapped to data rows by sequential nearest-neighbor selection.

• kernel: EnergyKernel() (default) or GaussianKernel(σ); a :median bandwidth is resolved from the data at datasplit time and the resolved kernel is stored in result.method.
• ratio: fraction of rows assigned to the test set, in (0, 1).
• kappa: absolute per-iteration subsample size for stochastic optimization; nothing uses all rows every iteration. Stochastic mode runs only when kappa is below the number of rows of the target (the data, or the reference when one is given).
• tolerance: convergence when the largest squared displacement of any support point in one iteration is below this value. In stochastic mode the running-average weight decays with the iteration count, so convergence there partly reflects that step-size decay rather than the objective flattening out. For GaussianKernel on full data, convergence never fires before the second iteration, and also triggers when the relative objective decrease falls below an internal rtol = 1e-8 (not exposed here); in stochastic mode the Gaussian kernel uses the MM sweep and the displacement rule only.
• rng: source of all randomness (initialization, jitter, stochastic subsampling); pass a seeded RNG for reproducible splits.
• verbose: print per-iteration progress.

In high dimension the optimized points may move less than the spacing between rows, in which case the selection is the initial random sample; prefer HerdingSplitter there (see the Benchmarks page).

TwinningSplitter(; ratio = 0.2, start = :farthest, rng = Random.default_rng())

Split by twinning (Vakayil & Joseph 2022): the standardized rows are covered by n disjoint groups, each a row u_i and its nearest ungrouped neighbors, formed as a chain where u_{i+1} is the ungrouped row nearest to the farthest member of group i; the smaller side is {u_1, …, u_n}. The objective is the energy distance (the paper's Proposition 1 ties it to the SPlit objective), so kernel is a fixed EnergyKernel(). No optimizer, no kappa, serial; average cost O(pN log N).

• ratio: fraction of rows assigned to the test set, in (0, 1).
• start: the starting row u_1: :farthest (default) is the row farthest from the centroid of the standardized data and consumes no randomness, :random draws it with rng, an integer names the row.
• rng: used only by start = :random.

SplitResult.converged is always true and iterations is the number of groups formed. weights and reference are not defined for twinning and raise an ArgumentError.

Differences from the paper

The paper assumes 1/ratio = r is an integer and takes n = ⌈N/r⌉, with the last group absorbing the remainder. Here n follows the generic rule of datasplit (or the caller's n in selectrows) and r = ⌊N/n⌋, the N − rn leftover rows forming groups of r + 1 spread evenly along the chain; when N = rn the two agree exactly. Grouped rows are masked in the search structure (a k-d tree below TWINNING_BRUTE_FORCE_DIMENSION columns, otherwise a matrix brute-force search), rebuilt on the remaining rows once more than half of its rows are masked (the paper masks without rebuilding).

Examples

data = randn(MersenneTwister(1), 1_000, 5)
result = datasplit(TwinningSplitter(), data)          # deterministic
folds = multiplet(TwinningSplitter(), data, 5)        # 5 distribution-balanced folds

_compress(kernel, X, seq, g, δ_halve, rng; n_threads) -> Vector{Int}

Compress (Shetty, Dwivedi & Mackey 2022, Alg. 1) of the row sequence seq with oversampling g: sequences of at most 4^g rows are returned as they are; longer ones are split into four consecutive parts, compressed recursively, concatenated, and halved by _symmetrized_halve. Returns about 2^g √(length(seq)) rows.

_compress_plus_plus(kernel, X, n; delta, rng, n_threads) -> (rows, swaps)

Compress++ (Shetty, Dwivedi & Mackey 2022, Alg. 2): Compress a random permutation of the rows with g = max(4, ⌈log₂(2n/√N)⌉), then THIN the compressed set to n rows with kernel_thinning (whose data term and swap candidates are the compressed rows, as in the paper). delta is split evenly between the halvings (each gets delta / 2K, K the number of HALVE calls) and THIN (delta / 2), a union bound in place of the paper's per-call schedule. If a compressed set has no more than n rows, g is increased by one and Compress is rerun.

_kernel_halving(kernel, Xt, seq, δ_step, rng; n_threads) -> (S₁, S₂)

Kernel halving (Dwivedi & Mackey 2024, Alg. 2) of the row sequence seq (indices into the columns of Xt = permutedims(X)): rows are taken two at a time, one goes to each half, and the assignment is swapped with probability min(1, ½(1 − α/a)₊) where α = Σ_{z∈S₂} f(z) − Σ_{z∈S₁} f(z), f = k(x,·) − k(x′,·), and a is the paper's adaptive threshold for the per-step failure probability δ_step. A trailing unpaired row is dropped.

_knn(s::MatrixSearch, q, k, skip) -> (idxs, dists)

The k columns of s.Xt nearest to q with !skip(j), by increasing distance. dists are Euclidean, matching NearestNeighbors' knn. k must not exceed the number of available (unskipped) columns; this method throws an ArgumentError where NearestNeighbors' knn would instead return fewer results.

For k <= 32, this selects by k passes of argmin over the scratch distances, allocating only the two length-k output vectors (idxs and dists); larger k falls back to partialsortperm, which allocates additional O(m) scratch.

_kt_swap(kernel, Xt, candidates, baseline, d, n_threads) -> (rows, swaps)

KT-SWAP (Dwivedi & Mackey 2022, Alg. 1b) with the target measure encoded in the data term d (d[z] = Σ_l v̄_l k(z, r_l)). Keeps the candidate (the baseline included) with the smallest (1/n²) Σ_{a,b∈S} k(a,b) − (2/n) Σ_{a∈S} d(a), then makes one pass over its positions, replacing each row by the row of the data outside the coreset that lowers the objective most (only if it lowers it; ties to the lowest row index). Returns the refined coreset in position order and the number of replacements.

_median_bandwidth(data, rows) -> GaussianKernel

Build the :median-heuristic GaussianKernel from the pairwise Euclidean distances among data[rows, :]: the bandwidth is the median of those distances. Throws ArgumentError when that median is zero (e.g. at least half of all row pairs coincide), since a zero bandwidth degenerates the kernel. Shared by the unweighted and weighted resolve methods, which differ only in how rows is drawn.

_nn(s::MatrixSearch, q, skip) -> (idx, dist)

Nearest column of s.Xt to q with !skip(j), by an explicit distance loop; allocation-free. dist is Euclidean, matching NearestNeighbors' nn.

_twin_groups(X, n, start, rng; search) -> Vector{Vector{Int}}

The n groups of Algorithm 1 on the standardized rows of X, in formation order; each group lists u_i first and then its neighbors by increasing distance to u_i. Grouped rows are masked in the search structure; when more than half of the structure's rows are masked it is rebuilt on the alive rows (an implementation detail, not from the paper; total rebuild work is O(N log N)). search is :kdtree, :brute_tree (a NearestNeighbors BruteTree), or :matrix (a MatrixSearch, which compiles once no matter how wide X is, unlike the NearestNeighbors trees); it defaults to :matrix at or above TWINNING_BRUTE_FORCE_DIMENSION columns and :kdtree below it.

apply_preprocessor(prep::Preprocessor, data) -> Matrix{Float64}

Encode data with the rules of prep and standardize every kept column with prep's mean and scale. A categorical level prep has not seen, a column-count or column-name mismatch, and a numeric/categorical kind mismatch are ArgumentErrors. Internal.

best(comparison::SplitComparison) -> (method, result)

The splitter/result pair with the lowest discrepancy.

compare(methods, data; kernel = EnergyKernel(), kwargs...) -> SplitComparison

Run datasplit with each splitter in methods on data and score every split with splitquality under kernel. kernel selects the scoring discrepancy; a :median bandwidth is resolved once on the preprocessed data with rng and the resolved kernel is stored in the comparison. Remaining keyword arguments are forwarded to splitquality.

weights is forwarded to datasplit, to the :median resolution of the scoring kernel, and to splitquality.

reference and reference_weights are forwarded to datasplit and splitquality, and a :median scoring kernel is then resolved on the encoded reference.

standardize = false uses a numeric matrix or vector as it is (no centering, scaling, or constant-column removal), for cosine-normalized embeddings; a DataFrame then raises an ArgumentError.

datasplit(splitter::AbstractSplitter, data) -> SplitResult

Split data (matrix, DataFrame, or vector; observations in rows) into train and test sets whose distributions are as similar as possible, using the method's own procedure; see the splitter types (SupportPointSplitter, HerdingSplitter).

weights (one non-negative entry per row; nothing for uniform) makes the split target the weighted empirical distribution Σ w̄ᵢ δ(xᵢ): the smaller subset is chosen to approximate it, preprocessing standardizes with the weighted mean and variance, and a :median bandwidth is resolved from rows drawn in proportion to the weights: this only changes the resolved bandwidth for datasets above 1000 rows; below that every row enters the median and the weights do not change it. The train/test labeling rule is unchanged. Weights proportional to duplication counts are equivalent to duplicating rows, up to the common column rescaling of the weighted standardization, which changes nothing under EnergyKernel or a :median bandwidth but does matter for a fixed numeric Gaussian bandwidth.

reference (same kind and columns as data; optionally weighted by reference_weights) makes the chosen side approximate the distribution of reference instead of the data: preprocessing is fit on reference and applied to both sets, a :median bandwidth is resolved on the encoded reference, and candidates remain the rows of data. weights cannot be combined with reference. The train/test labeling rule is unchanged; result.selected names the side that holds the chosen rows. See selectrows for the indices alone.

standardize = false uses a numeric matrix or vector as it is (no centering, scaling, or constant-column removal), for cosine-normalized embeddings; a DataFrame then raises an ArgumentError.

energydistance(X, Y; estimator = Exact(), rng = Random.default_rng(),
               n_threads = Threads.nthreads(), subsample = nothing, repeats = 8)

Energy distance between two samples whose rows are observations.

estimator selects how it is computed: Exact (default) evaluates every pairwise term block-wise, threaded over n_threads; Subsample averages the exact statistic over repeats random size-m row subsets drawn with rng: this estimate carries a positive bias of order 1/m, so use it to compare splits rather than as an absolute value; RandomSlices averages k random one-dimensional projections drawn with rng (unbiased). The subsample = m, repeats = r keywords are a compatibility path equivalent to estimator = Subsample(m, r).

weights_x and weights_y (one non-negative entry per row; nothing means uniform) give the energy distance between the weighted empirical distributions Σᵢ w̄ᵢ δ(xᵢ) and Σⱼ v̄ⱼ δ(yⱼ), with the weights scaled to sum one: 2 Σ w̄ᵢ v̄ⱼ ‖xᵢ − yⱼ‖ − Σ w̄ᵢ w̄ₖ ‖xᵢ − xₖ‖ − Σ v̄ⱼ v̄ₗ ‖yⱼ − yₗ‖. Weights proportional to duplication counts are equivalent to duplicating rows. A constant weight vector is treated as nothing, so uniform weights take the unweighted path and reproduce its value exactly.

Vectors are treated as single-column samples.

fit_preprocessor(data; weights = nothing, extra = nothing) -> Preprocessor

Learn the preprocessing on data: categorical columns are Helmert-encoded over the canonical-order union of their levels in data and in extra (the set the preprocessor will also be applied to). An encoded column constant on both data and extra (or constant on data when extra is not given) carries no information and is dropped. A column constant on data but varying on extra is kept instead of dropped: it is centered at data's constant value and scaled by extra's own (unweighted) spread, so it standardizes data's rows to exactly 0 and penalizes extra rows away from that value. Every other kept column gets the mean and scale of data (weighted forms when weights is given, as in preprocess); weights never applies to extra. With extra === nothing, this is exactly the previous behavior and preprocess stays bit-identical. Internal.

helmert_matrix(nlevels::Int) -> Matrix{Float64}

Helmert contrast matrix of size (nlevels, nlevels - 1). Column j assigns 1/j to the first j levels and -1 to level j + 1.

herd(kernel, X, n; weights = nothing, target = nothing, target_weights = nothing,
     n_threads = Threads.nthreads()) -> Vector{Int}

Select n rows of X by kernel herding: the first row maximizes the mean kernel value to the data, and each later row maximizes mean_l k(x, x_l) − (1/(T+1)) Σ_t k(x, s_t) over rows not yet selected (Chen, Welling & Smola 2010, Eq. 8). The data term is the exact mean over all N rows, including the self-term, computed threaded (O(N²)); approximating it with a DiscrepancyEstimator was evaluated and rejected: the Design experiments page records the measurement: the noise such estimators introduce is correlated across candidate rows and makes greedy selection unreliable. The procedure is deterministic given the data and a numeric kernel; ties go to the lowest row index. Cost O(N² + nN).

weights (one non-negative entry per row, nothing for uniform) replaces the data term by Σₗ w̄ₗ k(x, xₗ) with w̄ scaled to sum one, so the selection targets the weighted empirical distribution; the selected-set term is unchanged. A constant weight vector is treated as nothing, so uniform weights take the unweighted path and reproduce it exactly.

target (a matrix with the same columns as X) replaces the data term by the mean kernel value to the rows of target, weighted by target_weights when given; the candidates stay the rows of X. Cost O(NM + nN). weights and target are mutually exclusive.

Examples

X = randn(MersenneTwister(1), 200, 3)
herd(GaussianKernel(1.0), X, 40)

isresolved(kernel) -> Bool

Whether every kernel parameter is numeric (no data-dependent placeholders).

kernel_thinning(kernel, X, n; delta = 0.5, weights = nothing, target = nothing,
                target_weights = nothing, compress = :never,
                n_threads = Threads.nthreads(),
                rng = Random.default_rng()) -> (rows, swaps)

Select n rows of X by generalized kernel thinning with the target kernel (Dwivedi & Mackey 2022). For n ≤ N/2: with m = ⌊log₂(N/n)⌋, the first L = n·2^m rows of a random permutation are split by m rounds of kernel halving into 2^m candidate subsets of size n (KT-SPLIT), and KT-SWAP keeps the candidate (or a uniform random baseline) with the smallest MMD² to the target measure and refines it by one pass of best single-row swaps over all N rows. For n > N/2, the result is the complement of a kernel-thinning selection of the N - n rows not chosen (same rng order, delta, compress, and target measure; see "Differences from the paper" for why this is a reasonable rule). delta is the failure probability δ of the kernel-thinning guarantees: the papers' δ_i = δ/L, applied as δ_i/m at every halving step; weights, target, target_weights define the target measure as in herd and act on KT-SWAP only. Rough kernel-evaluation estimates, used only to choose the path: O(L²) for KT-SPLIT, O(N²) for the data term, O(nN) for KT-SWAP, all threaded. Deterministic given rng and independent of n_threads.

compress selects Compress++ (Shetty, Dwivedi & Mackey 2022) in place of plain kernel thinning: :never (default) never runs it; :always always runs it and requires weights === target === nothing, since Compress++ is defined for the data's own distribution, not a weighted or reference target; :auto runs it exactly when weights === target === nothing and _compress_pays_off judges it cheaper than plain kernel thinning at this N and n, which needs n far below N: not at a 20% split ratio, but below roughly 10% for N ≥ 10⁴.

Differences from the paper

Target-kernel thinning is used (no square-root kernel). The papers thin N to ⌊N/2^m⌋; here n is given and only L = n·2^m shuffled rows enter KT-SPLIT, the rest take part through KT-SWAP (equal to the paper when N/n = 2^m). Swap candidates exclude rows already in the coreset so the result is a set of distinct rows. weights/target change only the KT-SWAP objective. For n > N/2 the two halves of kernel halving are symmetric: for a subset S of size n_c = N - n with complement C, the mean-embedding identity μ_N = (n_c/N) μ_S + ((N - n_c)/N) μ_C gives MMD(C, P_N) = (n_c/(N - n_c)) · MMD(S, P_N) ≤ MMD(S, P_N) when n_c ≤ N/2, so the complement is at least as close to the data as the thinned set; for a weighted or reference target this identity does not hold exactly, and the rule is applied the same way regardless. Compress++'s HALVE is kernel thinning of the block's own rows rather than the paper's generic halving primitive, δ is split evenly across the halvings and the final THIN rather than following the paper's per-call schedule, the four-way splits of Compress are of possibly-uneven sizes rather than exactly ℓ/4, and the oversampling g is tied to the requested n (via _compress_g) rather than fixed across a run.

kernelgrad!(g, kernel, u, v) -> g

Write the gradient of kernelvalue(kernel, u, v) with respect to u into g.

kernelvalue(kernel, u, v) -> Float64

Evaluate the kernel at two points (vectors of equal length). For EnergyKernel the value is −‖u − v‖, the (conditionally negative definite) kernel whose MMD² is the energy distance.

mmd(X, Y, kernel; estimator = Exact(), rng = Random.default_rng(),
    n_threads = Threads.nthreads(), subsample = nothing, repeats = 8)

Squared maximum mean discrepancy (Gretton et al. 2012) between two samples whose rows are observations, under kernel.

estimator selects how it is computed: Exact (default) accumulates the V-statistic mean k(X,X) + mean k(Y,Y) − 2 mean k(X,Y) block-wise, threaded over n_threads; Subsample averages the exact statistic over repeats random size-m row subsets drawn with rng: this estimate carries a positive bias of order 1/m, so use it to compare splits rather than as an absolute value; RandomFeatures estimates it with D random Fourier features drawn with rng (unbiased), and is defined for GaussianKernel only. For EnergyKernel() this delegates to energydistance with the same estimator (so RandomSlices works and RandomFeatures raises ArgumentError). A :median bandwidth is resolved on the pooled rows of X and Y with rng. The subsample = m, repeats = r keywords are a compatibility path equivalent to estimator = Subsample(m, r).

weights_x and weights_y (one non-negative entry per row of X and Y; nothing means uniform) turn each sample into a weighted empirical distribution, Σᵢ w̄ᵢ δ(xᵢ) with w̄ scaled to sum one, and the statistic becomes Σ w̄ᵢ w̄ₖ k(xᵢ, xₖ) + Σ v̄ⱼ v̄ₗ k(yⱼ, yₗ) − 2 Σ w̄ᵢ v̄ⱼ k(xᵢ, yⱼ). Subsample draws rows uniformly and rescales the drawn weights to sum one. Weights proportional to duplication counts are equivalent to duplicating rows. A constant weight vector is treated as nothing, so uniform weights take the unweighted path and reproduce its value exactly. A :median bandwidth is still resolved on the pooled rows without weights; pass a numeric bandwidth or resolve it with datasplit/splitquality for a weighted one.

multiplet(splitter::AbstractSplitter, data, k; strategy = :sequential,
          weights = nothing, reference = nothing, reference_weights = nothing)
  -> Vector{Vector{Int}}

Partition the rows of data into k folds whose distributions each resemble the whole data (Vakayil & Joseph 2022, Section 5). Returns k sorted index vectors with sizes differing by at most one; 2 ≤ k ≤ N.

• :sequential (the paper's S1, default): select ⌊N/k⌋ rows with the splitter, then ⌊N'/(k−1)⌋ of the remaining N' rows, and so on; the last fold is what remains. Each run preprocesses the remaining rows afresh (or, with reference, fits the preprocessing on the reference as selectrows does). The splitter's own ratio is ignored. Any splitter.
• :halving (S2): split every part into its selected half (⌊N_part/2⌋ rows) and the complement, level by level, until k parts exist; k must be a power of two. Any splitter.
• :single (S3): one twinning run with ⌊N/k⌋ groups of at least k rows; fold j collects the j-th member (by neighbor rank) of every group, and the N mod k-style leftovers above rank k go one each to the first folds. TwinningSplitter only; one pass, and the paper measures it slightly behind S1 and S2.

weights, reference, and reference_weights are forwarded to selectrows (a TwinningSplitter rejects them). An integer start on a TwinningSplitter is interpreted within each run's remaining rows for :sequential and :halving (and must be a valid position there), so prefer :farthest or :random with those strategies.

standardize = false uses a numeric matrix or vector as it is (no centering, scaling, or constant-column removal), for cosine-normalized embeddings; a DataFrame then raises an ArgumentError.

Examples

data = randn(MersenneTwister(1), 500, 4)
folds = multiplet(TwinningSplitter(), data, 5)
holdout = data[folds[1], :]

optimal_split_ratio(x, y; method = :simple, degree = 2) -> Float64

Optimal fraction of data to assign to the test set, following the γ* = 1/(√p + 1) result of Joseph (2022, Eq. 11).

• method = :simple: p is taken as the number of encoded predictor columns of x (after preprocess: categorical columns Helmert-encoded, constant columns dropped) plus one, for the intercept.
• method = :regression: the paper's practical strategy (Section 3) for when the model is unknown: expand x into a larger feature set (e.g. polynomial terms up to degree), fit a linear regression on the full data, and select p via a model-selection criterion such as Cp or AIC. Not implemented in this release; calling with this method raises an error rather than silently falling back to :simple.

x may be an AbstractVector, AbstractMatrix, or DataFrame; y must have the same number of observations as x.

preprocess(data) -> Matrix{Float64}
preprocess(data, weights) -> Matrix{Float64}

Validate and transform data for splitting: reject missing values, encode categorical columns with Helmert contrasts, drop constant columns, and standardize every remaining column. Accepts AbstractMatrix, DataFrame, and AbstractVector inputs. Equivalent to fitting a Preprocessor on data and applying it to data.

With weights (one non-negative entry per row), standardization uses the weighted mean μⱼ = Σ w̄ᵢ xᵢⱼ and the unbiased weighted variance σⱼ² = Σ w̄ᵢ (xᵢⱼ − μⱼ)² / (1 − Σ w̄ᵢ²) with w̄ the weights scaled to sum one, which reduces to the n − 1 denominator of std for uniform weights; the encoding steps are the same. weights = nothing is the unweighted method. A constant weight vector is treated as nothing, so uniform weights take the unweighted path and reproduce it exactly.

resolve(kernel, data, rng) -> kernel with numeric parameters

Replace data-dependent placeholders. For GaussianKernel(:median) the bandwidth becomes the median pairwise Euclidean distance over min(size(data, 1), 1_000) rows drawn with rng (Gretton et al. 2012, "median heuristic"). Numeric kernels are returned unchanged.

resolve(kernel, data, rng, weights)

Weighted form of resolve: for GaussianKernel(:median) the rows behind the median heuristic are drawn with probability proportional to weights (without replacement), so the bandwidth reflects the weighted distribution. This only changes the resolved bandwidth for datasets above 1000 rows; below that every row enters the median and the weights do not change it. weights = nothing is the unweighted method; numeric kernels and EnergyKernel are returned unchanged. A constant weight vector is treated as nothing, so uniform weights take the unweighted path and reproduce it exactly.

select_nearest(data, points; search = _nearest_search(size(data, 2))) -> Vector{Int}

Indices of the data rows selected for the support points, in support-point order. data and points have observations in rows.

search is :kdtree (a NearestNeighbors k-d tree, doubling k on a query that returns only claimed rows) or :matrix (one MatrixSearch pass per point); the default follows NEAREST_BRUTE_FORCE_DIMENSION. Both implement the same claim-the-nearest-unclaimed-row rule and agree except on exact ties, which continuous data does not produce.

When the optimized points sit closer to their own starting row than to any other row, this returns the initial random sample unchanged. That is expected in high dimension, where the displacement is small relative to the row spacing (see the Benchmarks page).

selectrows(splitter::AbstractSplitter, data, n; weights = nothing,
       reference = nothing, reference_weights = nothing,
       standardize = true) -> Vector{Int}

Indices of the n rows of data the splitter chooses, in selection order (support-point order for SupportPointSplitter, greedy order for HerdingSplitter, formation order for TwinningSplitter, coreset position order for KernelThinningSplitter), without building a train/test partition. The chosen rows approximate the data's own distribution (weighted by weights) or, when reference is given, the distribution of reference (weighted by reference_weights): preprocessing is then fit on reference and applied to both, candidates stay the rows of data, and weights may not be given. Convergence diagnostics are reported by datasplit.

standardize = false uses a numeric matrix or vector as it is (no centering, scaling, or constant-column removal), for cosine-normalized embeddings; a DataFrame then raises an ArgumentError.

sphere_constant(p) -> Float64

κ_p = E_θ |⟨θ, e₁⟩| for θ uniform on the unit sphere of ℝ^p: κ_p = Γ(p/2) / (√π Γ((p+1)/2)), computed by the recursion κ_1 = 1, κ_2 = 2/π, κ_{p+2} = κ_p · p / (p + 1).

splitquality(data, result::SplitResult;
             kernel = EnergyKernel(), estimator = nothing,
             exact_threshold = 20_000, subsample = nothing, repeats = 8,
             rng = Random.default_rng(), n_threads = Threads.nthreads()) -> Float64

Discrepancy between the train and test rows of data under the same preprocessing datasplit applied: the energy distance by default, or mmd under kernel. Smaller is better.

estimator = nothing (the default) computes exactly (Exact) when the total row count is at most exact_threshold, and otherwise falls back to a fixed DiscrepancyEstimator chosen by the selection experiment on the Design experiments page (currently RandomSlices(64) for EnergyKernel and RandomFeatures(512) for GaussianKernel; see _fallback_estimator). Pass any DiscrepancyEstimator to override. The old subsample = m, repeats = r keywords are a compatibility path: when subsample is given explicitly, it always wins over estimator and maps to Subsample(m, r).

weights (one non-negative entry per row of data) applies the weighted preprocessing datasplit used and compares the weighted train rows with the weighted test rows, each side's weights rescaled to sum one.

reference (with optional reference_weights) scores the rows on the result.selected side against reference instead of comparing train with test: preprocessing is fit on reference and applied to both, and the discrepancy is between the selected rows (uniform) and the reference (weighted). weights cannot be combined with reference.

standardize = false uses a numeric matrix or vector as it is (no centering, scaling, or constant-column removal), for cosine-normalized embeddings; a DataFrame then raises an ArgumentError.

support_points(kernel::EnergyKernel, data, n; kwargs...) -> (points, converged, iterations)

Compute n support points for data (rows are observations) under the energy kernel by the closed-form majorization–minimization (MM) sweep of Mak & Joseph (2018). The sweep costs one pass over the data and the point set and never increases the objective on full data. Returns the points, whether the point-movement tolerance was reached, and the number of iterations actually used.

Convergence compares the largest squared displacement of any support point in one iteration to tolerance. In stochastic mode (kappa !== nothing), the running-average weight for iteration i is n0 / (i + n0) with n0 = 0.2n, which decays toward zero as iterations proceed, so convergence there partly reflects this step-size decay rather than the objective flattening out. n0 = 0.2n is an implementation constant, not from the papers, chosen by a small convergence experiment (see _n0_factor, an internal tuning knob not exposed on SupportPointSplitter).

weights (one non-negative entry per row, nothing for uniform) makes the points approximate the weighted empirical distribution Σ w̄ᵢ δ(xᵢ): the data sums in the MM update carry ŵᵢ = N w̄ᵢ, which is exactly 1.0 for uniform weights. In stochastic mode _subsampling (internal) selects how the kappa rows are drawn: :uniform draws them uniformly and rescales their weights to mean one within the subsample; :proportional draws them with probability proportional to the weights and treats the subsample as uniform (this needs at least kappa rows with positive weight). The default was chosen by the weighted-kappa experiment on the Design experiments page. A constant weight vector is treated as nothing, so uniform weights take the unweighted path and reproduce it exactly.

target (a matrix with the same columns as data) makes the points approximate the empirical distribution of target instead of data: the data term of the objective runs over the rows of target, weighted by target_weights (sum-one normalized, nothing for uniform; a constant vector is treated as nothing), while the initial points and the bounding box come from data, whose rows the points are later rounded to. In stochastic mode kappa subsamples the rows of target. weights is only for the case without a target; giving both is an ArgumentError. A target with duplicate rows is jittered by 1e-3 of its column range like the data, so weighting a reference is equivalent to duplicating its rows only up to that jitter.

support_points(kernel::GaussianKernel, data, n; kwargs...)

Support points under a Gaussian kernel. On full data (kappa === nothing, or kappa at or above the number of target rows), minimize the squared MMD between the point set and the data by projected gradient descent with Armijo backtracking on the projected step (the objective never increases across accepted steps). In stochastic mode (kappa below the number of target rows), support points instead come from the Gaussian MM sweep (mean-shift data term, majorized repulsion, see the Methods page), with the energy path's running-average blend and the displacement rule only, no line search. The kernel must be resolved (numeric bandwidth); datasplit resolves it.

The first trial step of the full-data path is scale-aware: t0 = 0.1 * scale / max(‖∇f‖, eps()), with scale the median per-dimension data range and ‖∇f‖ the largest gradient row norm, so the initial move is a tenth of the data scale regardless of the point count or gradient magnitude (Armijo backtracking and the 2t warm start on later iterations are unchanged). Convergence never fires before the second iteration, and then when either the largest squared displacement is below tolerance or the relative objective decrease |f_{t-1} - f_t| / max(|f_t|, 1e-12) is below rtol. f is _mmd_objective, which omits the constant data self-term and is bounded in [-1, 1] for a Gaussian kernel, so rtol acts as an absolute per-iteration tolerance on that bounded objective, not a relative tolerance on the (orders-of-magnitude smaller) true MMD². rtol applies to the full-data path only; in stochastic mode convergence uses only the displacement rule, as on the energy path.

weights (one non-negative entry per row, nothing for uniform) makes the points minimize the MMD² to the weighted empirical distribution Σ w̄ᵢ δ(xᵢ): the data term of the objective and of the gradient carries w̄, with ŵ = N w̄ (exactly 1.0 for uniform weights) inside the gradient loop, so unweighted results are unchanged.

target (a matrix with the same columns as data) makes the points approximate the empirical distribution of target instead of data: the data term of the objective runs over the rows of target, weighted by target_weights (sum-one normalized, nothing for uniform; a constant vector is treated as nothing), while the initial points and the bounding box come from data, whose rows the points are later rounded to. weights is only for the case without a target; giving both is an ArgumentError. A target with duplicate rows is jittered by 1e-3 of its column range like the data, so weighting a reference is equivalent to duplicating its rows only up to that jitter.

test_indices(result::SplitResult) -> Vector{Int}

train_indices(result::SplitResult) -> Vector{Int}