Multivariate Statistical Divergences

• Multivariate Statistical Divergences are rigorous frameworks that measure dissimilarity between probability distributions using techniques like Bregman and KL divergences.
• They underpin applications in clustering, model selection, and robust inference across machine learning, information theory, and data science.
• Efficient estimators and geometric algorithms enable scalable computation and statistical modeling in high-dimensional and structured data settings.

Multivariate statistical divergences provide rigorous frameworks to quantify, analyze, and interpret the dissimilarity between multivariate probability distributions. They underpin a wide array of theoretical and applied methodologies in machine learning, information theory, statistics, and data-driven sciences. Modern research has established a rich landscape of such divergences, including parametric, nonparametric, geometric, information-theoretic, and optimization-driven forms, each tailored to specific properties of high-dimensional data, structured models, or inference schemes.

1. Bregman Divergences and Generalized Centroids

Bregman divergences constitute a unifying class of distortion measures generated by strictly convex, differentiable functions $F : \mathbb{R}^d \to \mathbb{R}$. For points $p, q$ in the domain, the divergence is defined as

$$D_F(p \| q) = F(p) - F(q) - \langle p-q, \nabla F(q) \rangle.$$

This generalizes Euclidean distance (quadratic generator) and information-theoretic divergences (e.g., Kullback–Leibler with $F(x) = x \log x - x$). Most Bregman divergences are not symmetric; consequently, centroid computation must distinguish between right-type (minimizing $\sum D_F(p_i \| c)$), left-type ($\sum D_F(c \| p_i)$), and symmetrized forms via symmetrized divergences $S_F(p, q) = \frac{1}{2}[D_F(p \| q) + D_F(q \| p)]$.

The paper (0711.3242) presents the following results:

• Right-type centroid:

$c_R^F = \frac{1}{n} \sum_{i=1}^n p_i$

(arithmetic mean; independent of $F$),

• Left-type centroid:

$$c_L^F = (\nabla F)^{-1}\left( \frac{1}{n} \sum_{i=1}^n \nabla F(p_i) \right)$$

(generalized mean under the mapping induced by $\nabla F$),

• Symmetrized centroid:

Absence of closed form; an efficient geodesic-walk algorithm (dichotomic search along the segment between $c_R^F$ and $c_L^F$) locates the unique intersection with the Bregman bisector.

This construction is exploited in domains where the objects of comparison are distributions over simplices (e.g., histograms for image analysis) or parameter spaces (e.g., natural parameters for multivariate normals), allowing for computation of cluster centroids under entropy-type measures, which is critical for scalable, geometry-aware clustering of high-dimensional and matrix-variate data. The method extends existing ad hoc schemes (e.g., Veldhuis' convex programs for KL) and avoids solving nonlinear matrix equations as in Riccati-based centroid computations.

2. Divergences in Exponential Families and Information Geometry

Within exponential families, divergence concepts acquire a strong geometric underpinning. The canonical divergence in a dually flat space (i.e., a statistical manifold endowed with affine connection and dual coordinates $\theta$ and $\eta$) is given by

$$D(P \| Q) = \psi(\theta(P)) + \varphi(\eta(Q)) - \langle \theta(P), \eta(Q) \rangle,$$

where $\psi, \varphi$ are convex potentials linked by Legendre transform. This specializes to Kullback–Leibler divergence in exponential or mixture families. The affine divergence—sum of canonical divergences in both directions—realizes the Jeffreys divergence (a symmetric measure): $D_A(P, Q) = D(P\|Q) + D(Q\|P).$

The paper (Nishiyama, 2018) introduces new divergence functions—affine, ψ-, and φ-divergences—encapsulating symmetrizations and skew variations such as the Bhattacharyya and Jensen–Shannon divergences. These divergences obey Euclidean-like relations: triangular laws, parallelogram identities, and a generalization of Lin's inequality for comparing potential- and coordinate-based divergences. Such geometric properties enable new algorithms for projection, clustering, and statistical modelling in multivariate spaces.

3. Alpha-Beta Log-Determinant and Spectral Divergences for SPD Matrices

The Alpha-Beta (AB) log-det divergence is a two-parameter family for comparing symmetric positive-definite (SPD) matrices

$$D_{AB}^{(\alpha, \beta)}(\mathbf{P}\|\mathbf{Q}) = \frac{1}{\alpha\beta} \log \det \left( \frac{\alpha\, (\mathbf{Q}^{-1}\mathbf{P})^\beta + \beta\, (\mathbf{Q}^{-1}\mathbf{P})^{-\alpha}}{\alpha+\beta} \right).$$

Relevant specializations recover Stein’s loss, the Jensen–Bregman LogDet (JBLD) divergence, Bhattacharyya divergence, the affine-invariant Riemannian metric (AIRM), KL divergence for Gaussian densities, and the S-divergence (square root leading to a metric structure on SPD cone) (Cichocki et al., 2014).

For multivariate Gaussian densities, gamma log-det divergences decompose into a covariance part (an AB divergence on covariances) and a mean part (Mahalanobis distance). These concepts naturally extend to multiway (Kronecker-structured) settings, crucial in multilinear tensor data and diffusion tensor imaging, with divergence measures decomposing additively over tensor modes.

Symmetrized forms (averaging AB divergences in both directions or via Jensen–Shannon-like expressions) further regularize these measures for applications such as covariance learning, SPD kernel design, or robust clustering.

4. Efficient Estimation and Computation in High Dimensions

Consistent and scalable estimation of f-divergences in high-dimensional settings is a critical challenge. Nonparametric approaches, such as optimally-weighted ensemble k–NN plug-in estimators, achieve the mean squared error (MSE) convergence rate $O(1/T)$ and asymptotic normality even in high dimensions (Moon et al., 2014). The approach constructs an ensemble of plug-in estimators via varying neighborhood size and combines them with bias-optimized weights solved via convex programming.

Information-theoretic divergences, such as the Kullback–Leibler, Rényi-α, and total correlation, now admit fast, minimax rate-optimal neural network estimators through variational representations constrained via finite neural architectures. The error decomposes into function approximation (controllable by the width of the network) and statistical estimation (controlled by sample size), with the total error matching parametric rates under compact support and sufficient smoothness (Sreekumar et al., 2021).

Moreover, the Radon–Nikodym-based k–NN estimator for the general graph divergence measure recasts mutual, conditional, and total information measures as estimable KL divergences between empirical joint and model-restricted distributions (Rahimzamani et al., 2018). This estimator is robust to mixtures, manifold-supported distributions, and mixed discrete-continuous variables.

5. Structure and Decomposition: Hierarchical and Spectral Analysis

For joint distributions compared to product reference models, as in the hierarchical KL decomposition,

$$KL(P_k \| Q^{\otimes k}) = \sum_{i=1}^k KL(P_i \| Q) + C(P_k),$$

where $C(P_k)$ is the total correlation (multi-information). Through Möbius inversion on the subset lattice, $C(P_k)$ further splits into hierarchical interaction information—distinctly identifying pairwise, triplet, and higher-order dependencies (Cook, 12 Apr 2025). This algebraically exact decomposition allows diagnostics into whether divergence arises from marginal mis-specification or statistical dependence among variables, clarifying sources of model mismatch.

Additionally, for location families with fixed scale, all f-divergences reduce to strictly increasing functions of the Mahalanobis distance between means. For scale families (fixed mean), the divergence depends only on the eigenvalues of scale matrix ratios—matrix spectral divergences (Nielsen et al., 2022). This correspondence holds for a broad class of generator functions and standard densities, facilitating computational tractability in high-dimensional problems, e.g., in model selection for Gaussian mixtures or robustified Cauchy models.

6. Applications: Clustering, Robust Inference, and Experimental Design

The discussed divergences are operationalized in multiple domains:

• Clustering: Entropic (symmetrized Bregman) centroids support center-based clustering of distributions (e.g., images, multivariate normals) (0711.3242), and symmetrized Hölder divergences empirically improve clustering accuracy over Cauchy–Schwarz divergences for Gaussians (Nielsen et al., 2017).
• Model Selection and Comparison: Shannon entropy and KL divergence facilitate maximum entropy model comparison and clustering in skew-normal/canonical fundamental skew-normal families (Muniz et al., 2014), with Monte Carlo estimation schemes available for non-closed-form divergences.
• Optimal Design and Testing: Bregman divergences constructed from optimal experimental design functionals (power means or simplicial measures) yield robust, moment-sensitive measures for distinguishing multivariate normals and testing identity of mean/covariance via empirically optimized ROC/AUC curves (Pronzato et al., 2018).
• Random Projection and Robust EM: For mixture models, projection-based estimation and divergence quantification via Kolmogorov–Smirnov on univariate projections provide robust alternatives to EM in the presence of outliers or singularity, and support statistical comparison of random partitions in clustering (Fraiman et al., 15 Mar 2025).

7. Geometry, Algebraic Decomposition, and Quantum Extensions

The information-geometric underpinnings—dually flatness, canonical divergences, and affine coordinates—facilitate generalized projection theorems (Pythagorean, triangular) and motivate the vector space analogy for divergence analysis (Nishiyama, 2018). Barycentric (Choquet-like) decompositions of extensive monotone divergences reveal that all classical (and quantum) divergences can be expressed as barycenters over a test spectrum of extremal divergences, fully characterized in the classical case and shown to cover the range of quantum relative entropies and matrix mean divergences in quantum information (Haapasalo, 23 Sep 2025). This identifies the core “building blocks” of divergence theory and formalizes the independence and variety observed among proposed quantum divergences.

This synthesis of multivariate statistical divergences—encompassing theoretical generality, practical computability, geometric structure, and robustness—forms the basis for modern, scalable statistical inference, model selection, and data analysis across disciplines involving high-dimensional, structured, or heterogeneous data.