BIDistances: From distance computation to theory-guided selection

library(BIDistances)

Why BIDistances?

Choosing a distance is part of the scientific model, not merely a technical preprocessing decision. A distance determines which observations are regarded as similar, which neighbourhoods are formed, and which structures a distance-based clustering algorithm can recover.

BIDistances connects four stages of this workflow:

  1. Compute distances through one dispatcher.
  2. Classify their mathematical properties and required conditions.
  3. Compare and select candidates through their distance distributions when distance-based cluster structures are the objective.
  4. Scale weighted Minkowski calculations across parallelDist, an internal multicore implementation, and optional OpenCL kernels.

The package also contains the Gene Ontology-derived TF-IDF distance and specialist functions for probability distributions, nearest-neighbour structures, time series, and toroidal coordinates.

Access to distance methods

The following release-specific count uses canonical method names accepted by DistanceMatrix() for numerical data. Aliases, the parallelDist user-defined custom route, direct-only functions, and the configurable mixed-data constructions in manydist are not counted.

MethodCount = data.frame(
  AccessLayer = c(
    "Package-defined DistanceMatrix routes",
    "Predefined parallelDist methods",
    "Additional philentropy identifiers"
  ),
  CanonicalRoutes = c(19, 41, 34),
  Availability = c(
    "Core",
    "Core import",
    "Suggested package"
  ),
  stringsAsFactors = FALSE
)

knitr::kable(MethodCount, row.names = FALSE)
AccessLayer CanonicalRoutes Availability
Package-defined DistanceMatrix routes 19 Core
Predefined parallelDist methods 41 Core import
Additional philentropy identifiers 34 Suggested package

The core installation therefore exposes 60 canonical named numerical-data routes. With philentropy installed, the total is 94. The 34 additional philentropy identifiers are obtained after removing exact overlaps and the chebyshev and squared_euclidean aliases already resolved by the dispatcher.

The package-wide choice is broader than this count because manydist supplies mixed-type presets and custom component combinations, while functions such as TWED_Distance() and ToroidalEuclidean_Distance() are intentionally called directly.

Data and Distance conventions

Most methods compare rows. The main exceptions are DTW_Distance() and MSMD_Distance(), which compare time series stored in columns, and EndresSchindelin_Distance(), which compares columns as samples of variables.

One dispatcher for different distance families

data("Hepta")
Data = Hepta$Data

D_euclidean = DistanceMatrix(
  Data,
  method = "euclidean"
)

D_minkowski_p3 = DistanceMatrix(
  Data,
  method = "minkowski",
  p = 3,
  threads = 2
)

D_cosine = DistanceMatrix(
  Data,
  method = "Cosine_Distance"
)

c(
  EuclideanRows = nrow(D_euclidean),
  MinkowskiRows = nrow(D_minkowski_p3),
  CosineRows = nrow(D_cosine)
)
#> EuclideanRows MinkowskiRows    CosineRows 
#>           212           212           212

The ordinary lower-case method names are first offered to parallelDist. Canonical names ending in _Distance select native package routes. When a name is unavailable in parallelDist, DistanceMatrix() can use philentropy if that package is installed. Mixed numerical and categorical data are handled through method = "manydist".

Theory-guided comparison of distance distributions

The theoretical motivation is the goal of finding distance-based cluster structures: distances inside a partition should be smaller than distances between partitions. Under this goal, a multimodal full distance distribution can indicate different regimes of small and large distances. DistanceDistributions() operationalizes this exploratory comparison before clustering.

DistanceDistributionAnalysis() reports the observed Hartigan dip statistic in SelectionStatistics$DipStatistic. Candidate vectors with at most 72,000 pairwise values retain diptest::dip.test(). Above that threshold, all analyzable distance columns are processed in one integrated C++ batch with an explicit asymptotic p-value calibration and a workspace reused between columns.

The full comparison and density visualization are intentionally shown without being executed during vignette construction because they are computationally intensive and exercise optional plotting internals. Run these chunks interactively to reproduce the result.

set.seed(42)

Selection = DistanceDistributions(
  Data = Data,
  DistanceMethods = c(
    "euclidean",
    "manhattan",
    "minkowski",
    "chord"
  ),
  CosineNonParallel = TRUE,
  CorrelationDist = TRUE,
  PlotIt = FALSE,
  PlotSampleSize = 5000
)

Selection$DistanceChoice
Selection$SelectionStatistics[, c(
  "Distance", "DipStatistic", "DipPValue", "BimodalityAmplitude", "Status"
)]

When minkowski is included, the current implementation examines several exponents in addition to the requested base method. Some candidates with 0 < p < 1 are fractional dissimilarities rather than metrics. Their inclusion can be useful for exploration, but their mathematical status must be kept separate from that of Minkowski metrics with p >= 1.

Selection$ggobject

The returned object contains:

This procedure addresses a specific modeling goal. Multimodality is evidence for distance-based structure under that goal; it is not a universal proof of cluster validity, biological meaning, or superiority for every downstream algorithm.

Mathematical properties are explicit

The word distance is used broadly in software. Some methods are metrics, others are pseudometrics, squared metric forms, non-metric dissimilarities, or divergences. DistanceProperties() makes these distinctions inspectable.

Properties = DistanceProperties()
Relevant = grepl(
  "Minkowski|Fractional|Tfidf|DTW|Mahalanobis|Cosine",
  Properties$Function
)

knitr::kable(
  Properties[Relevant, c("Function", "Classification", "Conditions")],
  row.names = FALSE
)
Function Classification Conditions
Minkowski_Distance metric or pseudometric p >= 1; all finite-p weights > 0 give a metric; zero weights can collapse distinct vectors; p = Inf is available only with unit weights
Fractional_Distance metric or non-metric distance/dissimilarity p >= 1 gives the lp metric; 0 < p < 1 gives Deza’s fractional lp distance, which generally violates the triangle inequality; the package does not return its p-th power
Mahalanobis_Distance quadrance (squared Mahalanobis metric) the covariance or precision matrix must be symmetric positive definite; the function returns (x-y)’ A^{-1} (x-y) without a square root; sqrt(output) is a Mahalanobis metric, up to Deza’s optional positive determinant normalization
Cosine_Distance non-metric distance/dissimilarity returns 1 - cosine similarity; positive scalar multiples collapse, and the triangle inequality can fail; two zero vectors have distance 0 and a zero/non-zero pair has distance 1
Tfidf_Distance / Tfidf_dist pseudometric on input rows the distance is the absolute difference between scalar weights; different rows may share a weight.
DTW_Distance non-metric distance/dissimilarity the implementation rejects non-finite or negative final pair costs, mirrors each calculated value, and sets the diagonal to zero; ordinary DTW can give zero for distinct series and generally violates the triangle inequality.

This is important when a downstream algorithm assumes the triangle inequality or identity of indiscernibles. For example:

Weighted Minkowski distances on CPU and OpenCL

For finite p >= 1, the implemented weighted Minkowski distance is

\[ d_{p,w}(x,y) = \left(\sum_{k=1}^{d} w_k |x_k-y_k|^p\right)^{1/p}. \]

The three backends have different practical roles:

Backend Role
parallelDist Established multithreaded baseline; weights are represented by coordinate scaling
multicore Internal RcppParallel implementation and CPU fallback
opencl Optional OpenCL implementation with full, batched-output, and blockwise execution
auto Tries OpenCL and falls back to the internal multicore implementation with a warning

The Euclidean entry points remain explicit rather than being hidden by the new generalization. DistanceMatrix(..., method = "euclidean") uses parallelDist, EuclideanMulticore_Distance() exposes the internal unweighted CPU implementation, EuclideanGPU_Distance() provides the established weighted CPU/OpenCL route, and ToroidalEuclidean_Distance() represents periodic two-dimensional geometry. Minkowski_Distance(..., p = 2) reuses the optimized Euclidean OpenCL implementation instead of replacing it with a generic power kernel.

Cross-backend agreement on the CPU

MinkowskiData = as.matrix(iris[1:40, 1:4])
Weights = c(1, 2, 0.5, 1)

D_parallel = Minkowski_Distance(
  Data = MinkowskiData,
  p = 3,
  Weights = Weights,
  backend = "parallelDist",
  threads = 2
)

D_multicore = Minkowski_Distance(
  Data = MinkowskiData,
  p = 3,
  Weights = Weights,
  backend = "multicore",
  threads = 2
)

c(
  MaximumAbsoluteDifference = max(abs(D_parallel - D_multicore)),
  EqualWithinTolerance = isTRUE(
    all.equal(D_parallel, D_multicore, tolerance = 1e-10)
  )
)
#> MaximumAbsoluteDifference      EqualWithinTolerance 
#>              7.771561e-16              1.000000e+00

The same route is available through the package dispatcher:

D_dispatch = DistanceMatrix(
  MinkowskiData,
  method = "Minkowski_Distance",
  p = 3,
  Weights = Weights,
  backend = "multicore",
  threads = 2
)

isTRUE(all.equal(D_dispatch, D_multicore, tolerance = 1e-10))
#> [1] TRUE

OpenCL route

D_opencl = Minkowski_Distance(
  Data = MinkowskiData,
  p = 3,
  Weights = Weights,
  backend = "opencl",
  Mem = 2
)

max(abs(D_opencl - D_multicore))

For p = 2, the OpenCL route delegates to the established EuclideanGPU_Distance() implementation. General exponents use separate Minkowski kernels, so the optimized Euclidean code remains isolated. For p = Inf, unit weights are required and the maximum distance is returned.

Before a large OpenCL calculation, the memory plan can be inspected without allocating the full matrix:

MemoryPlan = calculateMemoryDemandGPU(
  n = 70000,
  d = 784,
  mem = 6
)
MemoryPlan
#> $maxBatchSize
#> [1] 10720
#> 
#> $minNrBatches
#> [1] 7
#> 
#> $batchSizes
#> [1] 10720 10720 10720 10720 10720 10720  5680
#> 
#> $Version
#> [1] 1

The GPU backend still returns a dense pairwise matrix, whose storage grows quadratically with the number of cases. The memory planner chooses among complete, output-batched, and input-blocked execution, but it does not remove the quadratic size of the final result.

Gene Ontology-derived TF-IDF distances

Tfidf_Distance() converts a gene-by-GO-term feature matrix into one scalar TF-IDF weight per gene and then computes absolute pairwise differences between those weights.

data("Hearingloss_N109")
Gene2Term = Hearingloss_N109$FeatureMatrix_Gene2Term

GO_Result = Tfidf_Distance(
  Gene2Term,
  tf_fun = mean
)

c(
  Genes = nrow(Gene2Term),
  GOTerms = ncol(Gene2Term),
  DistanceRows = nrow(GO_Result$Distance)
)
#>        Genes      GOTerms DistanceRows 
#>          109          829          109

GO_Result$TfidfWeights[1:5]
#> [1] 1.269761 1.363305 1.661398 1.919593 1.363305
GO_Result$Distance[1:5, 1:5]
#>           [,1]      [,2]      [,3]      [,4]      [,5]
#> [1,] 0.0000000 0.0935443 0.3916371 0.6498323 0.0935443
#> [2,] 0.0935443 0.0000000 0.2980928 0.5562880 0.0000000
#> [3,] 0.3916371 0.2980928 0.0000000 0.2581952 0.2980928
#> [4,] 0.6498323 0.5562880 0.2581952 0.0000000 0.5562880
#> [5,] 0.0935443 0.0000000 0.2980928 0.5562880 0.0000000

The distance-only result can be obtained through DistanceMatrix():

GO_Distance = DistanceMatrix(
  Gene2Term,
  method = "Tfidf_Distance",
  tf_fun = mean
)

isTRUE(all.equal(GO_Distance, GO_Result$Distance))
#> [1] TRUE

This route illustrates a central advantage of the package: a domain-derived distance, its mathematical classification, and the general distance-distribution workflow are available in one software environment.

Time-series and mixed-data routes

MSMD_Distance() and DTW_Distance() interpret each column as one time series and use row order as temporal order. No time variable should be supplied.

TimeSeries = cbind(
  SeriesA = c(0, 1, 2, 1, 0),
  SeriesB = c(0, 1, 1.5, 1, 0),
  SeriesC = c(2, 1, 0, 1, 2)
)

D_msmd = DistanceMatrix(
  TimeSeries,
  method = "MSMD_Distance",
  ParameterC = 1
)
D_msmd
#>         SeriesA SeriesB SeriesC
#> SeriesA     0.0     0.5     6.0
#> SeriesB     0.5     0.0     5.5
#> SeriesC     6.0     5.5     0.0

TWED_Distance() requires explicit timestamps and is intentionally called directly. This avoids guessing a time variable from the standard matrix interface.

When manydist is installed, mixed numerical and categorical data can be passed through the same dispatcher:

MixedData = data.frame(
  expression = c(2.1, 1.7, 4.2, 3.8),
  age = c(42, 38, 61, 58),
  subtype = factor(c("A", "A", "B", "B"))
)

D_mixed = DistanceMatrix(
  MixedData,
  method = "manydist",
  preset = "gower"
)

Which entry point should be used?

Task Recommended entry point
Compute one of many pairwise matrices DistanceMatrix()
Compare candidate distance distributions for clustering DistanceDistributions()
Inspect metric or divergence status DistanceProperties()
Compute weighted Minkowski distances with selectable backends Minkowski_Distance()
Compute the optimized weighted Euclidean CPU/OpenCL route EuclideanGPU_Distance()
Obtain both GO-derived distances and TF-IDF weights Tfidf_Distance()
Use explicit timestamps for TWED TWED_Distance()
Compute toroidal coordinates ToroidalEuclidean_Distance()

References

[Thrun, 2021] Thrun, M. C.: The Exploitation of Distance Distributions for Clustering, International Journal of Computational Intelligence and Applications, Vol. 20(3), 2150016, DOI: 10.1142/S1469026821500164, 2021.

[Stier & Thrun, 2023] Stier, Q., & Thrun, M. C.: Deriving Homogeneous Subsets from Gene Sets by Exploiting the Gene Ontology, Informatica, Vol. 34(2), pp. 357-386, DOI: 10.15388/23-INFOR517, 2023.

[Aggarwal et al., 2001] Aggarwal, C. C., Hinneburg, A., & Keim, D. A.: On the Surprising Behavior of Distance Metrics in High Dimensional Space, Database Theory — ICDT 2001, pp. 420-434, 2001.

mirror server hosted at Truenetwork, Russian Federation.