--- title: "Introduction to Coreset" author: "Martin R. Smith" date: today format: html bibliography: ../inst/REFERENCES.bib csl: ../inst/apa.csl vignette: > %\VignetteIndexEntry{Introduction to Coreset} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ```{r} #| label: setup #| include: false knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library("Coreset") ``` The Coreset package selects a subset that represents a fixed pool of *N* items, based on approximate and exact solutions to a suite of objectives: The **Max-Min Diversity Problem** (MMDP, the discrete *p*-dispersion objective) selects $k$ elements such that the minimum distance between any pair of selected elements is as large as possible; the chosen elements are maximally separated. This can reward selections that leave the interior of the set unrepresented. This objective is suited to defining a representative sample from a fixed pool: picking biological specimens for sequencing that span available diversity, or choosing a representative subset of protein structures from a database. The **Max-Sum Diversity Problem** (the "maximum diversity problem") selects $k$ elements that maximize the *total* pairwise distance they contain. As every pair in the selection contributes to the max-sum score, the optimum tends to favour elements spread across the whole extent of the pool rather than elements that merely avoid a single close neighbour. The **Max-Mean Dispersion Problem** maximizes the *average* distance between selected elements; this differs from the Max Sum objective in that the number of elements to be selected is not specified in advance. The **discrete *k*-centre problem** selects $k$ elements such that the maximum distance from any element in the original set to a selected element is as small as possible. In ensuring that each point has a nearby representative, this objective can select points that reflect a central compromise, rather than selections that are closer to more local points. Whereas dispersion spreads to the extremes, covering gravitates into the interior. This objective is useful when selecting centres that represent each point in a dataset: for example, siting fire stations to guarantee that all buildings can be reached within a given response time. The **maximum-entropy (maxdet) sampling** objective selects the $k$ least redundant elements: that is, those that maximize the log-determinant of a similarity kernel built from the distances between them, and thus that together occupy the largest volume. ## Installation Install from CRAN with: ```{r} #| label: install #| eval: false install.packages("Coreset") ``` ## Quick start Our examples employ a built-in R `dist` object that contains road distances (in km) between 21 European cities. ```{r} #| label: quickstart # Load the `eurodist` dist object data("eurodist") # Set a seed for a reproducible selection set.seed(1) # Select 4 maximally dispersed cities ffPick <- FarFirst(4L, eurodist) # View distances between selected cities as.matrix(eurodist)[ffPick, ffPick] # Quickly extract the minimum distance for a given selection MinDist(eurodist, ffPick) ``` `FarFirst()` returned the indices of the four cities whose nearest-neighbour distance within the selection is largest. `MinDist()` reports that value explicitly. ## Methods at a glance Coreset provides solvers for each objective: | Function | Objective | Quality | Speed | Stochastic? | |---|---|---|---|---| | `FarFirst()` | max-min / k-centre | Good (2-approximation) | Very fast | No | | `DropAdd()` | MMDP (max-min) | High (≈ 99 % optimal) | Fast | No | | `Grasp()` | MMDP (max-min) | Highest | Moderate | Yes (`set.seed()`) | | `ExactMaxMin()` | MMDP (max-min) | Optimal (NP-hard) | Slow | No | | `ExactMaxSum()` | Max-Sum Diversity | Optimal (NP-hard) | Slow | No | | `MaxMean()` | max-mean | High | Time-budgeted | Yes (`set.seed()`) | | `KCentre()` | k-centre | Near-optimal | Fast | No | | `ExactKCentre()` | k-centre | Optimal (NP-hard) | Slow | No | | `MaxEntropy()` | maxdet | High / Optimal for small `k` | Fast | No | To compute the score for an arbitrary selection of points under each objective, use `MinDist()` (max-min), `MeanDist()` (max-mean) and `KCentreRadius()` (k-centre). `ExactMaxSum()` and `MaxEntropy()` report their own objective value (total pairwise distance and log-determinant respectively) as the `score` attribute of the selection they return. ## Fast greedy selection The @Gonzalez1985 algorithm builds the selection greedily: start from a seed point, then repeatedly add whichever unselected point is farthest from the current selection. This greedy rule guarantees a 2-approximation to the optimal T~k~ and runs in O(*N* · *m*) time. The choice of seed point influences the quality of the selection. Peripheral seeds are more likely than central seeds to represent extremes of the data, and hence to feature in the optimal selection. `FarFirst()` supports several methods for identifying starting points for the greedy search. The default is a two-step strategy that starts from a random selects a point at random, then starts greedy search from the point furthest from By default `FarFirst()` runs three `"random_furthest"` starts, each of which takes a randomly selected point, moves to the point furthest from it, and begins the farthest-first sweep from there. The results of the pass with the highest T~k~ are returned. Three starts captures most of the quality gain from restarting — on benchmarks across a wide range of datasets the improvement curve bends early (knee at n ≈ 3–4); for higher-quality results, prefer `DropAdd()`, whose tabu search escapes the plateau that restarts cannot. ```{r} #| label: gonzalez-random set.seed(1) FarFirst(6L, eurodist) # default: best of three random starts ``` The number of random starts can be configured via the `nSeeds` argument. ```{r} #| label: gonzalez-ensemble set.seed(1) # Fewer starts provides a faster run, but the solution found may be inferior FarFirst(6L, eurodist, nSeeds = 2L) ``` Other strategies to select peripheral seeds are available via the `strategy` argument; these are described at `PickPoint()`. The best solution from the requested strategies is used. ```{r} #| label: gonzalez-seeds ffPick <- FarFirst(6L, eurodist, strategy = c("diameter", "anti_medoid")) MinDist(eurodist, ffPick) attr(ffPick, "winning_strategy") ``` With large datasets of points that are not associated with natural coordinates, it may not be feasible to compute and store all N×N distances. In such cases, a distance function may take the place of the distance matrix. This function is passed one index `i` at a time, and must return the distance from `i` to each other point. `N`, the number of objects, must also be specified. ```{r} #| label: gonzalez-oracle data("USArrests") arrestTypes <- USArrests[, c("Murder", "Assault", "Rape")] StateDist <- function(i) { diffs <- sweep(arrestTypes, 2, unlist(arrestTypes[i, ]), "-") sqrt(rowSums(diffs ^ 2)) } idx <- FarFirst(4L, StateDist, N = nrow(arrestTypes), strategy = 1L) arrestTypes[idx, ] ``` ## DropAdd tabu search The **DropAdd** heuristic [@Porumbel2011] refines an initial selection by alternately dropping and adding points from the selection; it typically reaches ≈ 99 % of the optimal T~k~. ```{r} #| label: dropadd daPick <- DropAdd(6L, eurodist, plateau = 500L) daPick labels(eurodist)[daPick] ``` The algorithm terminates after `plateau` iterations do not improve T~k~. Where *N* is too large for a distance matrix to fit in memory (roughly N > 46 000), pass a coordinate matrix via `DropAdd(points = ...)`. ## GRASP with path relinking **GRASP + path relinking** [@Resende2010] combines a randomised greedy construction phase with extended local search, and then refines an "elite" set of good solutions by interpolating between elite-pair trajectories (path relinking). It achieves the highest T~k~ of the three heuristics, at a proportionally higher cost. ```{r} #| label: grasp set.seed(0) grPick <- Grasp(6L, eurodist, plateau = 50L) grPick labels(eurodist)[grPick] attr(grPick, "pr_calls") # path-relinking calls performed ``` `plateau` controls how many consecutive non-improving GRASP iterations trigger termination; `maxSeconds` is available as an absolute time cap. ## Comparing methods on a simulated example To see how the methods relate visually, we generate 50 points in two dimensions and select *k* = 8 from each. ```{r} #| label: sim-data set.seed(1) # Seed selected such that FarFirst < DropAdd < Grasp pts <- matrix(rnorm(100), ncol = 2) # 50 points, 2 dimensions d50 <- dist(pts) k <- 8L ``` ```{r} #| label: compare-run set.seed(1) # Seed selected such that FarFirst < DropAdd < Grasp ffPick <- FarFirst(k, d50) da50Pick <- DropAdd(k, d50, plateau = 500L) gr50Pick <- Grasp(k, d50, plateau = 50L) ``` Even a small difference in T~k~ can correspond to a meaningfully more dispersed selection. ```{r} #| label: compare-scores scores <- c( FarFirst = attr(ffPick, "score"), DropAdd = attr(da50Pick, "score"), Grasp = attr(gr50Pick, "score") ) round(scores, 3) ``` Plotting the selections against the point cloud makes the differences concrete. When two methods select the same index, their symbols overlap; the T~k~ table above captures the quality distinction even when the visual overlap is high. ```{r} #| label: compare-plot #| fig-width: 6.5 #| fig-height: 6 #| fig-cap: "Selections returned by each method on 50 random 2-D points #| (k = 8). Coloured symbols mark selected points; grey circles are the #| full candidate set." methods <- list( FarFirst = ffPick, DropAdd = da50Pick, Grasp = gr50Pick ) cols <- c(FarFirst = "#E41A1C", DropAdd = "#377EB8", Grasp = "#4DAF4A") pchs <- c(FarFirst = 1L, DropAdd = 3L, Grasp = 4L) plot(pts, pch = 1L, col = "grey75", asp = 1L, xlab = "x", ylab = "y", frame.plot = FALSE, main = "Coreset method comparison") for (nm in names(methods)) { sel <- methods[[nm]] points(pts[sel, 1L], pts[sel, 2L], pch = pchs[nm], col = cols[nm], cex = 1.6) } legend_labels <- lapply(seq_along(methods), function(i) { bquote(.(names(methods)[i]) ~ (T[k] == .(sprintf("%.3f", scores[i])))) }) legend("topleft", legend = as.expression(legend_labels), pch = pchs, col = cols, pt.bg = cols, pt.cex = 1.4, bty = "n") ``` ## Exact solution For small instances (roughly N ≤ 25–30), `ExactMaxMin()` solves the problem to proven optimality via a node-packing integer programme [@Sayyady2016], using the **highs** solver (which we must first install). ```{r} #| label: install-highs #| eval: false install.packages("highs") ``` ```{r} #| label: exact-data set.seed(1L) pts30 <- matrix(rnorm(60L), ncol = 2L) d30 <- dist(pts30) ``` ```{r} #| label: exact #| eval: !expr requireNamespace("highs", quietly = TRUE) exPick <- ExactMaxMin(6L, d30, maxSeconds = 30L) attr(exPick, "proven") # TRUE ⟹ objective is the global optimum attr(exPick, "score") # Compare to the greedy heuristic on the same instance ff30Pick <- FarFirst(6L, d30) c(exact = attr(exPick, "score"), farFirst = MinDist(d30, ff30Pick)) ``` `$proven = TRUE` certifies that no selection can achieve a higher T~k~. `ExactMaxMin()` is NP-hard; it is wise to set a time budget once the number of points exceeds ~30 candidates, or to switch to a heuristic method. ## Scoring `MinDist()` computes the T~k~ objective for any index set. It accepts a `dist` object, a square distance matrix, or a coordinate matrix via the `points` argument: ```{r} #| label: MinDist MinDist(d50, ffPick) # from dist MinDist(as.matrix(d50), ffPick) # from square matrix MinDist(points = pts, idx = ffPick) # from coordinates ``` ## Max-sum diversity and Max-mean dispersion The Max-Sum Diversity Problem selects the $k$-subset with the highest total pairwise distance. `ExactMaxSum()` finds an optimal solution via per-node integer-program linearisation [@Kuo1993]. ```{r} #| label: exact-maxsum #| eval: !expr requireNamespace("highs", quietly = TRUE) smPick <- ExactMaxSum(6L, d30, maxSeconds = 30L) attr(smPick, "proven") # TRUE ⟹ objective is the global optimum attr(smPick, "score") # total pairwise distance within the selection ``` Where the number of elements is not specified _a priori_, Max-Sum diversity generalizes to the Max-Mean Dispersion Problem, which maximizes the sum of pairwise distances divided by the number of selected elements. `MaxMean()` implements reinforcement-learning tabu search [@Dieudonne2020]: each restart constructs a candidate selection, randomly at first, then guided by a *Q*-learning memory of which elements proved valuable. The selection is refined with a one-flip tabu search that adds or removes a single element at a time. The search continues until its time budget expires. The objective is only interesting when negative distances occur: otherwise the optimal selection tends to include all points. ```{r} #| label: maxmean-data set.seed(1) affinity <- matrix(runif(30L * 30L, min = -10, max = 10), nrow = 30L) affinity <- (affinity + t(affinity)) / 2 # symmetric diag(affinity) <- 0 ``` ```{r} #| label: maxmean set.seed(1) mmPick <- MaxMean(affinity, maxSeconds = 2) mmPick attr(mmPick, "size") # the algorithm chose this subset size attr(mmPick, "score") # achieved mean-dispersion objective ``` `MaxMean()` chooses the subset size to retaining only the elements that raise the average separation. `MeanDist()` scores any index set under the same objective, so a hand-picked selection can be compared directly: ```{r} #| label: meandist MeanDist(affinity, mmPick) # matches attr(mmPick, "score") MeanDist(affinity, 1:30) # the full set scores lower ``` ## Covering: the k-centre problem The above methods seek to spread the selection such that its members are mutually far apart. The k-centre problem instead minimizes the covering radius $R$, the largest distance from any point to its nearest chosen centre, such that no point of the pool is left far from a representative [@Gonzalez1985; @Hochbaum1985], typically resulting in selections that reach further into the interior of a sample. ### Heuristic approach `FarFirst()` is the quickest approximation to the K-centres problem. The CDSh algorithm implemented in `KCentres()` gives a more sophisticated heuristic [@GarciaDiaz2019; @GarciaDiaz2017]; its solutions are typically within 1–3.5 % of the optimum, compared to ~10% for `FarFirst()`. ```{r} #| label: kcentre centres <- KCentre(4L, eurodist) labels(eurodist)[centres] centres ``` `KCentreRadius()` scores any centre set by its covering radius (lower is better). CDSh covers at least as tightly as the Gonzalez 2-approximation baseline: ```{r} #| label: kcentre-radius ff <- FarFirst(4L, eurodist, strategy = "peripheral") c(KCentre = KCentreRadius(eurodist, centres), FarFirst = KCentreRadius(eurodist, ff)) ``` Like `MinDist()`, `KCentreRadius()` accepts a `points` coordinate matrix, allowing it to score a selection on a set of points whose distance matrix would be too large to fit in memory. ### Exact solver For small instances, `ExactKCentre()` finds the optimal solution, using a minimum-set-cover integer program approach akin to `ExactMaxMin()`, and using again the **highs** solver. ```{r} #| label: exact-kcentre #| eval: !expr requireNamespace("highs", quietly = TRUE) kc <- ExactKCentre(4L, eurodist) kc attr(kc, "proven") # TRUE ⟹ radius is the global covering optimum ``` The covering optimum is sometimes attained by fewer centres (once every point is covered, extra centres cannot lower the radius); `indices` then has length below the requested number, and the reported `radius` is still the proven optimum. `ExactKCentre()` is NP-hard, so — like `ExactMaxMin()` — it is a ground-truth reference for small instances, not a scalable method. The dispersion and covering optima differ even on this small example: dispersion selects cities at the rim of the map, while covering pulls inward to keep every city near a centre. ```{r} #| label: dispersion-vs-covering #| eval: !expr requireNamespace("highs", quietly = TRUE) disp <- ExactMaxMin(4L, eurodist) labels(eurodist)[disp] # dispersion: pushed to the extremes labels(eurodist)[kc] # covering: pulled toward the interior ``` ## Maximum-entropy (maxdet) selection The maximum entropy selection seeks to select $k$ elements that contain as much of the information of the original set as possible. The route to doing so is dropping elements that are redundant to selected elements. If redundancy between two points is equated to the overlap of volumes centred on each point, then maximum entropy can be cast as finding a selection that maximizes the log-determinant of its kernel block, $\log\det K_S$ [@Shewry1987], equivalently the maximum-a-posteriori mode of a determinantal point process [@Kulesza2012]. `MaxEntropy()` builds a radial-basis kernel from the distance matrix. This is repaired to be positive semi-definite where needed, since an arbitrary distance is not guaranteed to be of negative type. ```{r} #| label: maxentropy mePick <- MaxEntropy(4L, eurodist) labels(eurodist)[mePick] attr(mePick, "score") # the achieved log-determinant attr(mePick, "exact") # TRUE if certified by exact enumeration ``` The selection is built greedily by pivoted Cholesky, substituting exact enumeration when $\binom{n}{k}$ does not exceed `maxCombos`. `negMass` reports the fraction of spectral mass removed by the positive-semidefinite repair; as this fraction increases, so the result must be considered more approximate. ```{r} #| label: maxentropy-negmass attr(mePick, "negMass") ``` ## When to use which method For **dispersion** (spread the selection; maximize T~k~): | Scenario | Recommended | |---|---| | Speed matters most | `FarFirst()` (ensemble default) | | Deterministic, good quality valued | `DropAdd()` | | Best quality, `set.seed()` for reproducibility | `Grasp()` | | N > 46 000 (distance matrix infeasible) | `DropAdd(points = ...)` or `FarFirst(points = ...)` | | Arbitrary metric with no coordinate embedding | `FarFirst(, N = ...)` | | Proven optimum, N ≤ ~ 25–30, **highs** installed | `ExactMaxMin()` | | Score a selection's T~k~ | `MinDist()` | For **total dispersion** (fixed-size subset maximizing the *sum* of pairwise distances): | Scenario | Recommended | |---|---| | Proven optimum, N ≤ ~ 25–30, **highs** installed | `ExactMaxSum()` | For **average dispersion** (select a subset of the size that maximizes mean separation): | Scenario | Recommended | |---|---| | Maximize mean pairwise distance, size unfixed | `MaxMean()` | | Distances may be negative | `MaxMean()` | | Score a selection's max-mean objective | `MeanDist()` | For **covering** (minimise the radius; no point far from a centre): | Scenario | Recommended | |---|---| | Near-optimal covering, fast and deterministic | `KCentre()` (CDSh) | | Quick 2-approximation baseline | `FarFirst(strategy = "peripheral")` | | Proven optimum, small N | `ExactKCentre()` | | Score a centre set's covering radius (matrix-free at large N) | `KCentreRadius(points = ...)` | For **maximum-entropy (maxdet) selection** (density-blind volume maximization): | Scenario | Recommended | |---|---| | Maximize spanned volume / avoid near-duplicate selections | `MaxEntropy()` | | Proven optimum, `choose(n, k)` small | `MaxEntropy(exact = TRUE)` | ## Related problems `Coreset` selects a subset from a given set of elements. Several established packages solve neighbouring objectives: - **k-medoids / k-median** selects elements that minimize the mean distance from each element to its nearest centre. Implementations include: * [`cluster::pam()`](https://cran.r-project.org/package=cluster): generates the full *O(N²)* dissimilarity matrix, and hence caps at *n* ≤ 65 536; * [`banditpam`](https://cran.r-project.org/package=banditpam), a matrix-free $O(N \log N)$ implementation restricted to coordinate data; * `cluster::clara()` (PAM / FastPAM / FasterPAM); * [`ClusterR::Cluster_Medoids()`](https://cran.r-project.org/package=ClusterR). - **k-means** ([`stats::kmeans()`](https://rdrr.io/r/stats/kmeans.html)) selects elements so as to minimize the within-cluster sum of squares around centres that are coordinate means, not data points; as such, it applies only to Euclidean coordinates. k-means++ ([`TreeDist::KMeansPP()`]( https://ms609.github.io/TreeDist/reference/KMeansPP.html)) initializes its selection using D²-weighted seeding, a randomized relative of `FarFirst()`'s farthest-first traversal. - [`maximin`](https://cran.r-project.org/package=maximin) solves the related design problem of adding *new* points at positions that maximize the minimum inter-point distance. ## References