Skip to content

Top-level API

The top-level gmm_divergence namespace is intentionally small. It exposes the classes and functions most users reach for first, plus namespace modules for configuration objects.

What Lives At The Root

Symbol Purpose
Gaussian, GaussianMixture Core Gaussian-family classes
combine_gaussians Build a combined Gaussian mixture
kl_divergence, symmetric_kl_divergence, jensen_shannon_divergence Main divergence helpers
component_kl_matrix Pairwise component KL diagnostics
fit_mixture_weights, prune_mixture Main fitting helpers
DivergenceResult, FitResult Result containers

What Lives In Namespaces

Configuration and specialized helpers are grouped by namespace:

Namespace Contains
gmm_divergence.divergence Estimator configuration such as MonteCarlo
gmm_divergence.fitting Objectives, optimizers, and candidate selectors
gmm_divergence.sampling Sampling specifications such as Draw and Samples
gmm_divergence.covariance Covariance regularizers and epsilon heuristics
gmm_divergence.distributions Gaussian-family details and combination metadata

Example

import gmm_divergence as gd

p = gd.Gaussian.univariate(mean=0.0, variance=1.0)
q = gd.Gaussian.univariate(mean=1.0, variance=2.0)
method = gd.divergence.MonteCarlo(sampling=gd.sampling.Draw(50_000, rng=0))
result = gd.kl_divergence(p, q, method=method)

The sections below document the root exports explicitly. Namespace-only configuration objects are documented on their namespace pages.

Distributions

Gaussian dataclass 

Gaussian(mean: FloatArray, covariance: Covariance)

covariance instance-attribute 

covariance: Covariance

Covariance array of shape (n_features, n_features).

dim property 

dim: int

Dimensionality of the Gaussian.

mean instance-attribute 

mean: FloatArray

Mean array of shape (n_features,).

chol 

chol() -> FloatArray

Compute or retrieve the Cholesky factor of the covariance.

Source code in src/gmm_divergence/distributions/_gaussian.py 
92
93
94
95
96
97
98
99
def chol(self) -> FloatArray:
    """Compute or retrieve the Cholesky factor of the covariance."""
    if self._chol is not None:
        return self._chol

    chol = np.linalg.cholesky(self.covariance).astype(np.float64)
    object.__setattr__(self, "_chol", chol)
    return chol

component_arrays 

component_arrays() -> tuple[Weights, FloatArray, Covariances]

Return the mean and covariance as component arrays.

Source code in src/gmm_divergence/distributions/_gaussian.py 
134
135
136
def component_arrays(self) -> tuple[Weights, FloatArray, Covariances]:
    """Return the mean and covariance as component arrays."""
    return np.array([1.0]), self.mean[None, :], self.covariance[None, :, :]

from_arrays classmethod 

from_arrays(mean: ArrayLike, covariance: ArrayLike) -> Gaussian

Create a Gaussian instance from array-like inputs.

Source code in src/gmm_divergence/distributions/_gaussian.py 
28
29
30
31
@classmethod
def from_arrays(cls, mean: npt.ArrayLike, covariance: npt.ArrayLike) -> Gaussian:
    """Create a Gaussian instance from array-like inputs."""
    return cls(mean=cast("FloatArray", mean), covariance=cast("Covariance", covariance))

from_regularized_arrays classmethod 

from_regularized_arrays(
    mean: ArrayLike,
    covariance: ArrayLike,
    *,
    regularization: CovarianceRegularizer = "diagonal_loading",
) -> Gaussian

Create a Gaussian after explicitly regularizing its covariance.

This constructor keeps from_arrays strict while providing a convenient path for estimated or nearly singular covariances.

Source code in src/gmm_divergence/distributions/_gaussian.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
@classmethod
def from_regularized_arrays(
    cls,
    mean: npt.ArrayLike,
    covariance: npt.ArrayLike,
    *,
    regularization: CovarianceRegularizer = "diagonal_loading",
) -> Gaussian:
    """Create a Gaussian after explicitly regularizing its covariance.

    This constructor keeps `from_arrays` strict while providing a convenient
    path for estimated or nearly singular covariances.
    """
    regularized = regularize_covariance(covariance, method=regularization, batched=False)
    return cls(mean=cast("FloatArray", mean), covariance=regularized)

isotropic classmethod 

isotropic(mean: ArrayLike, variance: float = 1.0) -> Gaussian

Create an isotropic Gaussian instance.

Source code in src/gmm_divergence/distributions/_gaussian.py 
57
58
59
60
61
62
63
@classmethod
def isotropic(cls, mean: npt.ArrayLike, variance: float = 1.0) -> Gaussian:
    """Create an isotropic Gaussian instance."""
    mean = np.asarray(mean, dtype=np.float64)
    n_features = mean.shape[0]
    covariance = np.eye(n_features) * variance
    return cls(mean=mean, covariance=covariance)

log_det 

log_det() -> float

Compute or retrieve the log-determinant of the covariance.

Source code in src/gmm_divergence/distributions/_gaussian.py 
101
102
103
104
105
106
107
108
109
def log_det(self) -> float:
    """Compute or retrieve the log-determinant of the covariance."""
    if self._log_det is not None:
        return self._log_det

    chol = self.chol()
    log_det = 2 * np.sum(np.log(np.diag(chol)))
    object.__setattr__(self, "_log_det", log_det)
    return log_det

logpdf 

logpdf(x: ArrayLike) -> FloatArray

Evaluate the log-density of the Gaussian at given points.

Source code in src/gmm_divergence/distributions/_gaussian.py 
117
118
119
120
121
122
123
124
125
126
127
128
def logpdf(self, x: npt.ArrayLike) -> FloatArray:
    """Evaluate the log-density of the Gaussian at given points."""
    x = as_points(x, n_features=self.dim, name="x")

    d = self.mean.shape[0]
    chol = self.chol()
    diff = x - self.mean
    rhs = np.transpose(diff)
    y = np.transpose(np.linalg.solve(chol, rhs))
    mahalanobis = np.sum(y**2, axis=-1)
    log_det = self.log_det()
    return -0.5 * (d * np.log(2 * np.pi) + log_det + mahalanobis)

moments 

moments() -> tuple[FloatArray, Covariance]

Return the mean and covariance of the Gaussian.

Source code in src/gmm_divergence/distributions/_gaussian.py 
138
139
140
def moments(self) -> tuple[FloatArray, Covariance]:
    """Return the mean and covariance of the Gaussian."""
    return self.mean, self.covariance

pdf 

pdf(x: ArrayLike) -> FloatArray

Evaluate the density of the Gaussian at given points.

Source code in src/gmm_divergence/distributions/_gaussian.py 
130
131
132
def pdf(self, x: npt.ArrayLike) -> FloatArray:
    """Evaluate the density of the Gaussian at given points."""
    return np.exp(self.logpdf(x))

sample 

sample(n_samples: int, rng: Generator | int | None = None) -> FloatArray

Draw samples from the Gaussian.

Source code in src/gmm_divergence/distributions/_gaussian.py 
111
112
113
114
115
def sample(self, n_samples: int, rng: np.random.Generator | int | None = None) -> FloatArray:
    """Draw samples from the Gaussian."""
    n_samples = as_positive_sample_count(n_samples)
    rng = np.random.default_rng(rng)
    return rng.multivariate_normal(mean=self.mean, cov=self.covariance, size=n_samples)

standard classmethod 

standard(dim: int) -> Gaussian

Create a standard Gaussian instance with zero mean and identity covariance.

Source code in src/gmm_divergence/distributions/_gaussian.py 
65
66
67
68
@classmethod
def standard(cls, dim: int) -> Gaussian:
    """Create a standard Gaussian instance with zero mean and identity covariance."""
    return cls(mean=np.zeros(dim, dtype=np.float64), covariance=np.eye(dim, dtype=np.float64))

univariate classmethod 

univariate(mean: float = 0.0, variance: float = 1.0) -> Gaussian

Create a univariate Gaussian instance.

Source code in src/gmm_divergence/distributions/_gaussian.py 
49
50
51
52
53
54
55
@classmethod
def univariate(cls, mean: float = 0.0, variance: float = 1.0) -> Gaussian:
    """Create a univariate Gaussian instance."""
    return cls(
        mean=np.array([mean], dtype=np.float64),
        covariance=np.array([[variance]], dtype=np.float64),
    )

GaussianMixture dataclass 

GaussianMixture(weights: Weights, means: FloatArray, covariances: Covariances)

covariances instance-attribute 

covariances: Covariances

Covariance array of shape (n_components, n_features, n_features).

dim property 

dim: int

Dimensionality of the Gaussian mixture.

means instance-attribute 

means: FloatArray

Mean array of shape (n_components, n_features).

n_components property 

n_components: int

Number of components in the Gaussian mixture.

weights instance-attribute 

weights: Weights

Weight array of shape (n_components,).

as_gaussian 

as_gaussian(*, require_single: Literal[True]) -> Gaussian | None

as_gaussian(*, require_single: Literal[False] = False) -> Gaussian

as_gaussian(*, require_single: bool = False) -> Gaussian | None

Return a Gaussian approximation of the mixture using moment matching.

Source code in src/gmm_divergence/distributions/_mixture.py 
210
211
212
213
214
215
216
217
def as_gaussian(self, *, require_single: bool = False) -> Gaussian | None:
    """Return a Gaussian approximation of the mixture using moment matching."""
    if require_single and self.n_components > 1:
        return None
    if self.n_components == 1:
        return self.get_component(0)
    mean, covariance = self.moments()
    return Gaussian(mean=mean.astype(np.float64, copy=False), covariance=covariance)

chol 

chol() -> FloatArray

Compute or retrieve the cached Cholesky factors.

Source code in src/gmm_divergence/distributions/_mixture.py 
143
144
145
146
147
148
149
150
def chol(self) -> FloatArray:
    """Compute or retrieve the cached Cholesky factors."""
    if self._chol is not None:
        return self._chol

    chol = np.linalg.cholesky(self.covariances).astype(np.float64)
    object.__setattr__(self, "_chol", chol)
    return chol

component_arrays 

component_arrays() -> tuple[Weights, FloatArray, Covariances]

Return the weights, means, and covariances as arrays.

Source code in src/gmm_divergence/distributions/_mixture.py 
219
220
221
def component_arrays(self) -> tuple[Weights, FloatArray, Covariances]:
    """Return the weights, means, and covariances as arrays."""
    return self.weights, self.means, self.covariances

diagnostics 

diagnostics() -> MixtureDiagnostics

Return diagnostics for the Gaussian mixture.

Source code in src/gmm_divergence/distributions/_mixture.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
def diagnostics(self) -> MixtureDiagnostics:
    """Return diagnostics for the Gaussian mixture."""
    weight_sum = float(np.sum(self.weights))
    max_weight = float(np.max(self.weights))
    min_weight = float(np.min(self.weights))
    weights_entropy = float(-np.sum(self.weights * np.log(self.weights + 1e-12)))
    covar_condition_numbers = [float(np.linalg.cond(cov)) for cov in self.covariances]
    return MixtureDiagnostics(
        n_components=self.n_components,
        dim=self.dim,
        weight_sum=weight_sum,
        max_weight=max_weight,
        min_weight=min_weight,
        weights_entropy=weights_entropy,
        covar_condition_numbers=covar_condition_numbers,
    )

from_arrays classmethod 

from_arrays(
    weights: ArrayLike, means: ArrayLike, covariances: ArrayLike
) -> GaussianMixture

Create a Gaussian mixture from array-like parameters.

Source code in src/gmm_divergence/distributions/_mixture.py 
58
59
60
61
62
63
64
65
66
67
@classmethod
def from_arrays(
    cls, weights: npt.ArrayLike, means: npt.ArrayLike, covariances: npt.ArrayLike
) -> GaussianMixture:
    """Create a Gaussian mixture from array-like parameters."""
    return cls(
        weights=cast("Weights", weights),
        means=cast("FloatArray", means),
        covariances=cast("Covariances", covariances),
    )

from_components classmethod 

from_components(
    components: Sequence[Gaussian], weights: ArrayLike | None = None
) -> GaussianMixture

Create a Gaussian mixture from a sequence of Gaussian components and optional weights.

Source code in src/gmm_divergence/distributions/_mixture.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
@classmethod
def from_components(
    cls, components: Sequence[Gaussian], weights: npt.ArrayLike | None = None
) -> GaussianMixture:
    """Create a Gaussian mixture from a sequence of Gaussian components and optional weights."""
    means = np.array([comp.mean for comp in components], dtype=np.float64)
    covariances = np.array([comp.covariance for comp in components], dtype=np.float64)
    if weights is None:
        weights = np.ones(len(components), dtype=np.float64) / len(components)
    else:
        weights = np.asarray(weights, dtype=np.float64)
    return cls(weights=weights, means=means, covariances=covariances)

from_regularized_arrays classmethod 

from_regularized_arrays(
    weights: ArrayLike,
    means: ArrayLike,
    covariances: ArrayLike,
    *,
    regularization: CovarianceRegularizer = "diagonal_loading",
) -> GaussianMixture

Create a Gaussian mixture after explicitly regularizing covariances.

This constructor keeps from_arrays strict while providing a convenient path for estimated or nearly singular component covariances.

Source code in src/gmm_divergence/distributions/_mixture.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
@classmethod
def from_regularized_arrays(
    cls,
    weights: npt.ArrayLike,
    means: npt.ArrayLike,
    covariances: npt.ArrayLike,
    *,
    regularization: CovarianceRegularizer = "diagonal_loading",
) -> GaussianMixture:
    """Create a Gaussian mixture after explicitly regularizing covariances.

    This constructor keeps `from_arrays` strict while providing a convenient
    path for estimated or nearly singular component covariances.
    """
    regularized = regularize_covariance(covariances, method=regularization, batched=True)
    return cls(
        weights=cast("Weights", weights),
        means=cast("FloatArray", means),
        covariances=regularized,
    )

get_component 

get_component(index: int) -> Gaussian

Return the Gaussian component at the specified index.

Source code in src/gmm_divergence/distributions/_mixture.py 
170
171
172
173
174
175
176
def get_component(self, index: int) -> Gaussian:
    """Return the Gaussian component at the specified index."""
    if index < 0 or index >= self.n_components:
        msg = f"Component index {index} is out of bounds for {self.n_components} components."
        raise IndexError(msg)

    return Gaussian(mean=self.means[index], covariance=self.covariances[index])

log_dets 

log_dets() -> FloatArray

Compute or retrieve the cached log-determinants of the covariances.

Source code in src/gmm_divergence/distributions/_mixture.py 
152
153
154
155
156
157
158
159
160
def log_dets(self) -> FloatArray:
    """Compute or retrieve the cached log-determinants of the covariances."""
    if self._log_dets is not None:
        return self._log_dets

    chol = self.chol()
    log_dets = 2.0 * np.sum(np.log(np.diagonal(chol, axis1=1, axis2=2)), axis=1)
    object.__setattr__(self, "_log_dets", log_dets)
    return log_dets

logpdf 

logpdf(x: ArrayLike) -> FloatArray

Evaluate the log-density of the Gaussian mixture at given points.

Source code in src/gmm_divergence/distributions/_mixture.py 
137
138
139
140
141
def logpdf(self, x: npt.ArrayLike) -> FloatArray:
    """Evaluate the log-density of the Gaussian mixture at given points."""
    x = as_points(x, n_features=self.dim, name="x")

    return gmm_logpdf(x=x, gmm=self)

moments 

moments() -> tuple[FloatArray, Covariance]

Return the mean and covariance of the Gaussian mixture.

Source code in src/gmm_divergence/distributions/_mixture.py 
250
251
252
253
254
255
256
257
258
259
260
261
def moments(self) -> tuple[FloatArray, Covariance]:
    """Return the mean and covariance of the Gaussian mixture."""
    weights, means, covariances = self.component_arrays()
    mean = np.sum(weights[:, None] * means, axis=0)
    mean_delta = means - mean
    covariance = np.sum(
        weights[:, None, None]
        * (covariances + mean_delta[:, :, None] * mean_delta[:, None, :]),
        axis=0,
    )
    covariance = 0.5 * (covariance + covariance.T)
    return mean.astype(np.float64, copy=False), covariance

pdf 

pdf(x: ArrayLike) -> FloatArray

Evaluate the density of the Gaussian mixture at given points.

Source code in src/gmm_divergence/distributions/_mixture.py 
166
167
168
def pdf(self, x: npt.ArrayLike) -> FloatArray:
    """Evaluate the density of the Gaussian mixture at given points."""
    return np.exp(self.logpdf(x))

sample 

sample(n_samples: int, rng: Generator | int | None = None) -> FloatArray

Draw samples from the Gaussian mixture.

Source code in src/gmm_divergence/distributions/_mixture.py 
162
163
164
def sample(self, n_samples: int, rng: np.random.Generator | int | None = None) -> FloatArray:
    """Draw samples from the Gaussian mixture."""
    return sample_gmm(self, n_samples=n_samples, rng=rng)

select_components 

select_components(indices: ArrayLike) -> GaussianMixture

Return a new Gaussian mixture containing only the specified components.

Source code in src/gmm_divergence/distributions/_mixture.py 
178
179
180
181
182
183
184
185
186
187
188
189
def select_components(self, indices: npt.ArrayLike) -> GaussianMixture:
    """Return a new Gaussian mixture containing only the specified components."""
    indices = np.asarray(indices, dtype=np.intp)
    if np.any(indices < 0) or np.any(indices >= self.n_components):
        msg = f"Component indices must be in the range [0, {self.n_components})."
        raise IndexError(msg)

    return GaussianMixture(
        weights=as_weights(self.weights[indices], expected_length=indices.shape[0]),
        means=self.means[indices],
        covariances=self.covariances[indices],
    )

combine_gaussians 

combine_gaussians(
    sources: Sequence[Gaussian | GaussianMixture],
    weights: ArrayLike | None = None,
    *,
    include_mapping: Literal[False] = False,
) -> GaussianMixture

combine_gaussians(
    sources: Sequence[Gaussian | GaussianMixture],
    weights: ArrayLike | None = None,
    *,
    include_mapping: Literal[True],
) -> CombinedGaussianMixture

combine_gaussians(
    sources: Sequence[Gaussian | GaussianMixture],
    weights: ArrayLike | None = None,
    *,
    include_mapping: bool = False,
) -> GaussianMixture | CombinedGaussianMixture

Combine Gaussian mixtures and/or Gaussians into a single GaussianMixture.

A single Gaussian is internally converted to a GaussianMixture with one component. The resulting GaussianMixture has a number of components equal to the sum of the number of components in the input.

Parameters:

  • sources (Sequence[Gaussian | GaussianMixture]) –

    The input distributions to combine.

  • weights (ArrayLike, default: None ) –

    Weights for each input distribution. If None, equal weights are used.

  • include_mapping (bool, default: False ) –

    Whether to return a CombinedMixture with mapping information. If False, only the combined GaussianMixture is returned. The mapping information includes the index of the original input distribution and the local component index within that distribution for each component in the combined mixture.

Returns:

  • GaussianMixture or CombinedMixture –

    The combined GaussianMixture. If return_mapping is True, a CombinedMixture containing the combined mixture and mapping information is returned.

Source code in src/gmm_divergence/distributions/_combine.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def combine_gaussians(
    sources: Sequence[Gaussian | GaussianMixture],
    weights: npt.ArrayLike | None = None,
    *,
    include_mapping: bool = False,
) -> GaussianMixture | CombinedGaussianMixture:
    """Combine Gaussian mixtures and/or Gaussians into a single GaussianMixture.

    A single Gaussian is internally converted to a GaussianMixture with one
    component. The resulting GaussianMixture has a number of components equal to
    the sum of the number of components in the input.

    Parameters
    ----------
    sources : Sequence[Gaussian | GaussianMixture]
        The input distributions to combine.
    weights : ArrayLike, optional
        Weights for each input distribution. If None, equal weights are used.
    include_mapping : bool, default=False
        Whether to return a CombinedMixture with mapping information. If False,
        only the combined GaussianMixture is returned. The mapping information
        includes the index of the original input distribution and the local
        component index within that distribution for each component in the
        combined mixture.

    Returns
    -------
    GaussianMixture or CombinedMixture
        The combined GaussianMixture. If `return_mapping` is True, a
        CombinedMixture containing the combined mixture and mapping information
        is returned.

    """
    if len(sources) == 0:
        msg = "components must contain at least one distribution"
        raise ValueError(msg)

    weights_: list[float]
    if weights is None:
        weights_ = [1 / len(sources)] * len(sources)
    else:
        weights_arr = as_weights(weights, expected_length=len(sources), name="weights")
        weights_ = weights_arr.tolist()

    mixtures = [as_mixture(component) for component in sources]
    mixture = GaussianMixture(
        weights=np.concatenate([w * m.weights for w, m in zip(weights_, mixtures, strict=True)]),
        means=np.concatenate([m.means for m in mixtures], axis=0),
        covariances=np.concatenate([m.covariances for m in mixtures], axis=0),
    )

    if not include_mapping:
        return mixture

    mapping = MixtureMapping(
        source_index=np.concatenate([
            np.full(len(m.weights), i, dtype=np.intp) for i, m in enumerate(mixtures)
        ]),
        local_component_index=np.concatenate([
            np.arange(len(m.weights), dtype=np.intp) for m in mixtures
        ]),
    )
    return CombinedGaussianMixture(mixture=mixture, mapping=mapping)

Divergence Helpers

kl_divergence 

kl_divergence(
    p: GaussianLike,
    q: GaussianLike,
    /,
    *,
    method: KLMethod = "monte_carlo",
    prefer_closed_form: bool = True,
) -> DivergenceResult

Compute the Kullback--Leibler divergence between two Gaussian-family distributions.

Computes

\[ D_{\mathrm{KL}}(p \| q) = \mathbb{E}_{x \sim p} \left[ \log p(x) - \log q(x) \right]. \]

Here, p is treated as the reference distribution and q as the approximating distribution.

Parameters:

  • p (Gaussian or GaussianMixture) –

    The two Gaussian-family distributions to compare. They must have the same dimensionality.

  • q (Gaussian or GaussianMixture) –

    The two Gaussian-family distributions to compare. They must have the same dimensionality.

  • method (str or KL method configuration, default: "monte_carlo" ) –

    Method used to compute or estimate the KL divergence. Passing a string runs that method with its defaults. Use a method configuration object, such as divergence.MonteCarlo(sampling=sampling.Draw(50_000, rng=0)), for method-specific options.

  • prefer_closed_form (bool, default: True ) –

    If True, the function will attempt use closed form if both inputs are Gaussian, even if the user specified a different method.

Returns:

  • DivergenceResult –

    Result object containing the estimated KL divergence and metadata about the computation, such as the method used and whether the result is exact or approximate.

Notes

The KL divergence is asymmetric:

\[ D_{\mathrm{KL}}(p \| q) \neq D_{\mathrm{KL}}(q \| p). \]

Therefore, swapping p and q generally gives a different result.

Examples:

Compute the KL divergence using the default Monte Carlo estimator:

result = kl_divergence(p, q)
print(result.value)

Require a closed-form expression:

result = kl_divergence(p, q, method="closed_form")

Configure Monte Carlo sampling:

result = kl_divergence(
    p,
    q,
    method=divergence.MonteCarlo(sampling=sampling.Draw(50_000, rng=0)),
)

Use precomputed samples:

samples = p.sample(50_000, rng=0)
result = kl_divergence(
    p, q, method=divergence.MonteCarlo(sampling=sampling.Samples(samples))
)
Source code in src/gmm_divergence/divergence/_api.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def kl_divergence(
    p: GaussianLike,
    q: GaussianLike,
    /,
    *,
    method: KLMethod = "monte_carlo",
    prefer_closed_form: bool = True,
) -> DivergenceResult:
    r"""Compute the Kullback--Leibler divergence between two Gaussian-family distributions.

    Computes

    $$
    D_{\mathrm{KL}}(p \| q)
    =
    \mathbb{E}_{x \sim p}
    \left[
        \log p(x) - \log q(x)
    \right].
    $$

    Here, `p` is treated as the reference distribution and `q` as the
    approximating distribution.

    Parameters
    ----------
    p, q : Gaussian or GaussianMixture
        The two Gaussian-family distributions to compare. They must have the
        same dimensionality.
    method : str or KL method configuration, default="monte_carlo"
        Method used to compute or estimate the KL divergence. Passing a string
        runs that method with its defaults. Use a method configuration object,
        such as `divergence.MonteCarlo(sampling=sampling.Draw(50_000, rng=0))`, for
        method-specific options.
    prefer_closed_form : bool, default=True
        If `True`, the function will attempt use closed form if both inputs are
        Gaussian, even if the user specified a different method.

    Returns
    -------
    DivergenceResult
        Result object containing the estimated KL divergence and metadata about
        the computation, such as the method used and whether the result is exact
        or approximate.

    Notes
    -----
    The KL divergence is asymmetric:

    $$
    D_{\mathrm{KL}}(p \| q) \neq D_{\mathrm{KL}}(q \| p).
    $$

    Therefore, swapping `p` and `q` generally gives a different result.

    Examples
    --------
    Compute the KL divergence using the default Monte Carlo estimator:

    ```python
    result = kl_divergence(p, q)
    print(result.value)
    ```

    Require a closed-form expression:

    ```python
    result = kl_divergence(p, q, method="closed_form")
    ```

    Configure Monte Carlo sampling:

    ```python
    result = kl_divergence(
        p,
        q,
        method=divergence.MonteCarlo(sampling=sampling.Draw(50_000, rng=0)),
    )
    ```

    Use precomputed samples:

    ```python
    samples = p.sample(50_000, rng=0)
    result = kl_divergence(
        p, q, method=divergence.MonteCarlo(sampling=sampling.Samples(samples))
    )
    ```
    """
    _validate_same_dimension(p, q)
    spec, options = KL_REGISTRY.resolve(method)
    if prefer_closed_form and isinstance(p, Gaussian) and isinstance(q, Gaussian):
        return kl_closed_form(p, q)

    match spec.name:
        case "monte_carlo":
            options = cast_options(options, MonteCarlo)
            return kl_monte_carlo(p, q, sampling=options.sampling)
        case "unscented":
            return kl_unscented(p, q)
        case "gaussian_approximation":
            options = cast_options(options, MomentMatchedGaussian)
            return kl_gaussian_approximation(p, q, approximation=options.approximation)
        case "closed_form":
            p, q = _require_gaussian_pair(p, q, spec.name)
            return kl_closed_form(p, q)
        case "variational":
            return kl_variational(p, q)
        case _:
            msg = "Unhandled KL method registry entry."
            raise AssertionError(msg)

symmetric_kl_divergence 

symmetric_kl_divergence(
    p: GaussianLike,
    q: GaussianLike,
    /,
    *,
    method: KLMethod = "monte_carlo",
    prefer_closed_form: bool = True,
) -> DivergenceResult

Compute the symmetric KL divergence between two Gaussian-family distributions.

Computes

\[ D_{\mathrm{SKL}}(p, q) = \frac{1}{2} \left[ D_{\mathrm{KL}}(p \| q) + D_{\mathrm{KL}}(q \| p) \right]. \]

The same KL estimation method is used in both directions.

Parameters:

  • p (Gaussian or GaussianMixture) –

    The two Gaussian-family distributions to compare. They must have the same dimensionality.

  • q (Gaussian or GaussianMixture) –

    The two Gaussian-family distributions to compare. They must have the same dimensionality.

  • method (str or KL method configuration, default: "monte_carlo" ) –

    Method used for each directed KL estimate.

  • prefer_closed_form (bool, default: True ) –

    If True, each directed estimate will attempt to use closed form when both inputs are Gaussian.

Returns:

  • DivergenceResult –

    Result containing the symmetric KL value. For sampled methods, num_samples is the total sample count across both directed estimates when both counts are available.

Source code in src/gmm_divergence/divergence/_api.py 
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def symmetric_kl_divergence(
    p: GaussianLike,
    q: GaussianLike,
    /,
    *,
    method: KLMethod = "monte_carlo",
    prefer_closed_form: bool = True,
) -> DivergenceResult:
    r"""Compute the symmetric KL divergence between two Gaussian-family distributions.

    Computes

    $$
    D_{\mathrm{SKL}}(p, q)
    =
    \frac{1}{2}
    \left[
        D_{\mathrm{KL}}(p \| q) + D_{\mathrm{KL}}(q \| p)
    \right].
    $$

    The same KL estimation method is used in both directions.

    Parameters
    ----------
    p, q : Gaussian or GaussianMixture
        The two Gaussian-family distributions to compare. They must have the
        same dimensionality.
    method : str or KL method configuration, default="monte_carlo"
        Method used for each directed KL estimate.
    prefer_closed_form : bool, default=True
        If `True`, each directed estimate will attempt to use closed form when
        both inputs are Gaussian.

    Returns
    -------
    DivergenceResult
        Result containing the symmetric KL value. For sampled methods,
        `num_samples` is the total sample count across both directed estimates
        when both counts are available.
    """
    forward = kl_divergence(p, q, method=method, prefer_closed_form=prefer_closed_form)
    reverse = kl_divergence(q, p, method=method, prefer_closed_form=prefer_closed_form)
    return DivergenceResult(
        value=0.5 * (forward.value + reverse.value),
        method="symmetric_kl",
        num_samples=_sum_num_samples(forward, reverse),
    )

jensen_shannon_divergence 

jensen_shannon_divergence(
    p: GaussianLike,
    q: GaussianLike,
    /,
    *,
    method: KLMethod = "monte_carlo",
    prefer_closed_form: bool = True,
) -> DivergenceResult

Compute the Jensen-Shannon divergence between two Gaussian-family distributions.

Computes

\[ D_{\mathrm{JS}}(p, q) = \frac{1}{2}D_{\mathrm{KL}}(p \| m) + \frac{1}{2}D_{\mathrm{KL}}(q \| m), \qquad m = \frac{1}{2}p + \frac{1}{2}q. \]

For Gaussian and Gaussian-mixture inputs, the midpoint distribution m is represented as a Gaussian mixture using existing component-combination logic.

Parameters:

  • p (Gaussian or GaussianMixture) –

    The two Gaussian-family distributions to compare. They must have the same dimensionality.

  • q (Gaussian or GaussianMixture) –

    The two Gaussian-family distributions to compare. They must have the same dimensionality.

  • method (str or KL method configuration, default: "monte_carlo" ) –

    Method used for the two KL estimates against the midpoint mixture.

  • prefer_closed_form (bool, default: True ) –

    Passed through to kl_divergence. Note that the midpoint is generally a Gaussian mixture, so closed form is only available for methods and inputs that support it.

Returns:

  • DivergenceResult –

    Result containing the Jensen-Shannon divergence value. For sampled methods, num_samples is the total sample count across both directed estimates when both counts are available.

Source code in src/gmm_divergence/divergence/_api.py 
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
def jensen_shannon_divergence(
    p: GaussianLike,
    q: GaussianLike,
    /,
    *,
    method: KLMethod = "monte_carlo",
    prefer_closed_form: bool = True,
) -> DivergenceResult:
    r"""Compute the Jensen-Shannon divergence between two Gaussian-family distributions.

    Computes

    $$
    D_{\mathrm{JS}}(p, q)
    =
    \frac{1}{2}D_{\mathrm{KL}}(p \| m)
    +
    \frac{1}{2}D_{\mathrm{KL}}(q \| m),
    \qquad
    m = \frac{1}{2}p + \frac{1}{2}q.
    $$

    For Gaussian and Gaussian-mixture inputs, the midpoint distribution `m` is
    represented as a Gaussian mixture using existing component-combination
    logic.

    Parameters
    ----------
    p, q : Gaussian or GaussianMixture
        The two Gaussian-family distributions to compare. They must have the
        same dimensionality.
    method : str or KL method configuration, default="monte_carlo"
        Method used for the two KL estimates against the midpoint mixture.
    prefer_closed_form : bool, default=True
        Passed through to `kl_divergence`. Note that the midpoint is generally
        a Gaussian mixture, so closed form is only available for methods and
        inputs that support it.

    Returns
    -------
    DivergenceResult
        Result containing the Jensen-Shannon divergence value. For sampled
        methods, `num_samples` is the total sample count across both directed
        estimates when both counts are available.
    """
    _validate_same_dimension(p, q)
    midpoint = combine_gaussians([p, q], weights=[0.5, 0.5])
    p_to_midpoint = kl_divergence(p, midpoint, method=method, prefer_closed_form=prefer_closed_form)
    q_to_midpoint = kl_divergence(q, midpoint, method=method, prefer_closed_form=prefer_closed_form)
    return DivergenceResult(
        value=0.5 * (p_to_midpoint.value + q_to_midpoint.value),
        method="jensen_shannon",
        num_samples=_sum_num_samples(p_to_midpoint, q_to_midpoint),
    )

component_kl_matrix 

component_kl_matrix(p: GaussianLike, q: GaussianLike) -> FloatArray

Return pairwise Gaussian-component KL divergences.

The returned matrix has shape (p_components, q_components), where entry (i, j) is

\[ D_{\mathrm{KL}}\!\left(p_i \| q_j\right) \]

for component i of p and component j of q. A single Gaussian is treated as a one-component Gaussian mixture.

Parameters:

Returns:

  • FloatArray –

    Pairwise component KL matrix.

Source code in src/gmm_divergence/divergence/_api.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def component_kl_matrix(p: GaussianLike, q: GaussianLike, /) -> FloatArray:
    r"""Return pairwise Gaussian-component KL divergences.

    The returned matrix has shape `(p_components, q_components)`, where entry
    `(i, j)` is

    $$
    D_{\mathrm{KL}}\!\left(p_i \| q_j\right)
    $$

    for component `i` of `p` and component `j` of `q`. A single `Gaussian` is
    treated as a one-component Gaussian mixture.

    Parameters
    ----------
    p, q : Gaussian or GaussianMixture
        Gaussian-family distributions whose components are compared.

    Returns
    -------
    FloatArray
        Pairwise component KL matrix.
    """
    _validate_same_dimension(p, q)
    _, p_means, p_covariances = p.component_arrays()
    _, q_means, q_covariances = q.component_arrays()
    return pairwise_gaussian_kl(p_means, p_covariances, q_means, q_covariances)

Fitting Helpers

fit_mixture_weights 

fit_mixture_weights(
    p: Gaussian | GaussianMixture,
    q_i: Sequence[Gaussian | GaussianMixture],
    /,
    *,
    method: FitMethod = "softmax_lbfgsb",
    objective: FitObjective = "forward",
    x0: ArrayLike | None = None,
    candidate_selector: CandidateSelector | None = None,
) -> FitResult

Fit weights for a mixture of fixed candidate distributions.

Fits nonnegative weights w_i such that the weighted mixture

\[ q_w(x) = \sum_i w_i q_i(x) \]

approximates the reference distribution p according to the selected objective.

Parameters:

  • p (Gaussian or GaussianMixture) –

    Reference distribution.

  • q_i (sequence of Gaussian or GaussianMixture) –

    Candidate distributions whose weights are fitted.

  • method (str or optimizer configuration, default: 'softmax_lbfgsb' ) –

    Optimizer used for the weights. Passing a string runs that optimizer with defaults. Use SoftmaxLBFGSB(...) or SimplexSLSQP(...) for optimizer-specific options.

  • objective (str or WeightFitObjective configuration, default: 'forward' ) –

    Objective used for fitting. Passing a string runs that objective with defaults. Use ForwardKL(...), ReverseKL(...), BidirectionalKL(...), JensenShannon(...), or MomentMatching(...) for objective-specific options.

  • x0 (array - like, default: None ) –

    Initial weights for the optimized variables. If None, the optimizer's default initialization is used.

Returns:

  • FitResult –

    Result containing the fitted weights, fitted mixture, fit objective, objective value, forward/reverse KL diagnostics, and optimizer metadata.

Source code in src/gmm_divergence/fitting/_api.py 
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def fit_mixture_weights(
    p: Gaussian | GaussianMixture,
    q_i: Sequence[Gaussian | GaussianMixture],
    /,
    *,
    method: FitMethod = "softmax_lbfgsb",
    objective: FitObjective = "forward",
    x0: npt.ArrayLike | None = None,
    candidate_selector: CandidateSelector | None = None,
) -> FitResult:
    r"""Fit weights for a mixture of fixed candidate distributions.

    Fits nonnegative weights `w_i` such that the weighted mixture

    $$
    q_w(x) = \sum_i w_i q_i(x)
    $$

    approximates the reference distribution `p` according to the selected
    objective.

    Parameters
    ----------
    p : Gaussian or GaussianMixture
        Reference distribution.
    q_i : sequence of Gaussian or GaussianMixture
        Candidate distributions whose weights are fitted.
    method : str or optimizer configuration, optional
        Optimizer used for the weights. Passing a string runs that optimizer
        with defaults. Use `SoftmaxLBFGSB(...)` or `SimplexSLSQP(...)` for
        optimizer-specific options.
    objective : str or WeightFitObjective configuration, optional
        Objective used for fitting. Passing a string runs that objective with
        defaults. Use `ForwardKL(...)`, `ReverseKL(...)`, `BidirectionalKL(...)`,
        `JensenShannon(...)`, or `MomentMatching(...)` for objective-specific
        options.
    x0 : array-like, optional
        Initial weights for the optimized variables. If `None`, the optimizer's
        default initialization is used.

    Returns
    -------
    FitResult
        Result containing the fitted weights, fitted mixture, fit objective,
        objective value, forward/reverse KL diagnostics, and optimizer metadata.

    """
    method_spec, optimizer = _OPTIMIZER_REGISTRY.resolve(method)
    _objective_spec, objective_config = _OBJECTIVE_REGISTRY.resolve(objective)

    match method_spec.name:
        case "softmax_lbfgsb":
            optimizer = cast_options(optimizer, SoftmaxLBFGSB)
            objective_config = _cast_fit_objective(objective_config)
            return wfit.fit_mixture_weights(
                p=p,
                q_i=q_i,
                objective=objective_config,
                optimizer=optimizer,
                x0=x0,
                candidate_selection=candidate_selector,
            )
        case "simplex_slsqp":
            optimizer = cast_options(optimizer, SimplexSLSQP)
            objective_config = _cast_fit_objective(objective_config)
            return wfit.fit_mixture_weights(
                p=p,
                q_i=q_i,
                objective=objective_config,
                optimizer=optimizer,
                x0=x0,
                candidate_selection=candidate_selector,
            )
        case _:
            msg = "Unhandled fit optimizer registry entry."
            raise AssertionError(msg)

prune_mixture 

prune_mixture(
    mixture: GaussianMixture, *, min_weight: float = 0.0001
) -> GaussianMixture

Prune components of a Gaussian mixture with small weights.

This is a common post-processing step after fitting to remove components that contribute negligibly to the mixture, improving efficiency and interpretability.

Parameters:

  • mixture (GaussianMixture) –

    The mixture to prune.

  • min_weight (float, default: 0.0001 ) –

    Minimum weight threshold for keeping components. Components with weights below this threshold will be removed. Default is 1e-4.

Returns:

  • GaussianMixture –

    The pruned mixture with weights normalized to sum to one.

Raises:

  • ValueError –

    If all components are pruned.

Source code in src/gmm_divergence/fitting/_api.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def prune_mixture(mixture: GaussianMixture, *, min_weight: float = 1e-4) -> GaussianMixture:
    """Prune components of a Gaussian mixture with small weights.

    This is a common post-processing step after fitting to remove components
    that contribute negligibly to the mixture, improving efficiency and
    interpretability.

    Parameters
    ----------
    mixture : GaussianMixture
        The mixture to prune.
    min_weight : float, optional
        Minimum weight threshold for keeping components. Components with weights
        below this threshold will be removed. Default is 1e-4.

    Returns
    -------
    GaussianMixture
        The pruned mixture with weights normalized to sum to one.

    Raises
    ------
    ValueError
        If all components are pruned.
    """
    validate_nonnegative_finite(min_weight, name="min_weight")

    weights = mixture.weights
    keep_mask = weights >= min_weight
    if not np.any(keep_mask):
        msg = "All components were pruned, increase min_weight threshold."
        raise ValueError(msg)

    return mixture.select_components(np.nonzero(keep_mask)[0])

Results

DivergenceResult dataclass 

DivergenceResult(
    value: float,
    method: str | None = None,
    num_samples: int | None = None,
    monte_carlo_stats: MonteCarloStatistics | None = None,
)

Result of a divergence estimation.

method class-attribute instance-attribute 

method: str | None = None

The method used for estimation, if applicable.

monte_carlo_stats class-attribute instance-attribute 

monte_carlo_stats: MonteCarloStatistics | None = None

Monte Carlo statistics related to the estimation, if applicable.

Contains the sample mean, sample variance, standard error, and effective sample size of the pointwise divergence estimates when using Monte Carlo estimation.

num_samples class-attribute instance-attribute 

num_samples: int | None = None

The number of samples used for estimation, if applicable.

value instance-attribute 

value: float

The estimated divergence.

FitResult dataclass 

FitResult(
    weights: Weights,
    objective_value: float,
    scipy_result: OptimizeResult | None,
    fitted_mixture: CombinedGaussianMixture,
    fit_objective: FitObjective,
    fit_method: FitMethod,
    alpha: float | None = None,
    iterations: int | None = None,
    converged: bool | None = None,
    used_candidate_indices: list[int] | None = None,
)

Result of fitting a Gaussian mixture.

Primary result when using fit_mixture_weights and related functions.

alpha class-attribute instance-attribute 

alpha: float | None = None

Forward-objective weight used by bidirectional fitting, if applicable.

converged class-attribute instance-attribute 

converged: bool | None = None

Whether the optimization routine reported convergence, if applicable.

fit_method instance-attribute 

fit_method: FitMethod

The optimization method used to fit the mixture weights.

fit_objective instance-attribute 

fit_objective: FitObjective

The objective used to fit the mixture weights.

fitted_mixture instance-attribute 

fitted_mixture: CombinedGaussianMixture

The full combined Gaussian mixture corresponding to the fitted weights.

iterations class-attribute instance-attribute 

iterations: int | None = None

The number of iterations taken by the optimization routine, if applicable.

objective_value instance-attribute 

objective_value: float

The final scalar objective value minimized by the optimizer.

scipy_result instance-attribute 

scipy_result: OptimizeResult | None

The full result object returned by the optimization routine, if used.

weights instance-attribute 

weights: Weights

The fitted mixture weights as a 1D array.

candidate_weight_dict 

candidate_weight_dict() -> dict[int, float]

Return fitted weights keyed by original candidate index.

Source code in src/gmm_divergence/results.py 
92
93
94
def candidate_weight_dict(self) -> dict[int, float]:
    """Return fitted weights keyed by original candidate index."""
    return dict(self.candidate_weights())

candidate_weights 

candidate_weights() -> list[tuple[int, float]]

Return fitted weights paired with original candidate indices.

If candidate selection was used, indices refer to the original q_i sequence passed to fit_mixture_weights. Otherwise they are simply 0, 1, ..., n_candidates - 1.

Source code in src/gmm_divergence/results.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def candidate_weights(self) -> list[tuple[int, float]]:
    """Return fitted weights paired with original candidate indices.

    If candidate selection was used, indices refer to the original `q_i`
    sequence passed to `fit_mixture_weights`. Otherwise they are simply
    `0, 1, ..., n_candidates - 1`.
    """
    indices = (
        range(self.weights.shape[0])
        if self.used_candidate_indices is None
        else self.used_candidate_indices
    )
    return sorted(
        [
            (int(index), float(weight))
            for index, weight in zip(indices, self.weights, strict=True)
        ],
        key=operator.itemgetter(1),
        reverse=True,
    )