Implementations of measures

[sethaxen]

#3 proposes that we have some type to specify the base measure of a distribution. The goal of this issue is to reach a consensus on the following questions:

Should measures be explicitly or implicitly defined?

Explicit definitions increase flexibility at the cost of greater code complexity, while implicit definitions are simpler but less customizable.

Which measures will we define, and how will we define them?

Here are some popular measures:

• Hausdorff
• Invariant measure/Haar
• Lebesgue
• Volume measure

Any one of these can be normalized. In some cases these are equivalent. They also have different definitions. We should agree on a set of definitions for these and which should be included.

How will measure definitions relate to point representation?

e.g. the usual von Mises-Fisher density on the N-sphere is written wrt the unnormalized Hausdorff measure, and the usual representation is a vector with N+1 Cartesian coordinates. But what if a user wants to provide a vector of N spherical coordinates and compute the density wrt the Lebesgue measure on those coordinates? Other coordinates could in principle be used; how should the user specify which Lebesgue measure is being used?

Should we provide any functionality for changing measures?

This is related to the previous question.

For example, suppose a user is doing Monte Carlo sampling for a parameter on the positive N-hemisphere and wants to specify a prior. They then have 2 samplers. Sampler 1 proposes a move in the N-ball from an isotropic Gaussian in the first N coordinates, fills the N+1 entry by completing the norm, and then evaluates the density to accept or reject. Sampler 2 proposes a move on the geodesic of the N-sphere and then evaluates the density. Sampler 1 implicitly wants densities to be defined wrt the Lebesgue measure on the N-ball, while Sampler 2 implicitly wants densities to be defined wrt the Hausdorff measure on the N-sphere/positive N-hemisphere. Without a mechanism for changing measures, we would require the user manually compute the density correction if their sampler expects a different base measure than the one we implement.

[sethaxen]

Hausdorff measure

The Hausdorff measure $\mathcal{H}^k$ for a $k$-dimensional manifold $\mathcal{M}$ with Riemannian metric $g$ is generally defined as the limit of coverings of subsets $U_i \subset \mathcal{M}$. For $U \subset \mathcal{M}$, let $\mathrm{diam}(U) = \sup_{x,y \in U} d_g(x, y)$.
Then, for any $A \subset \mathcal{M}$, the $k$-dimensional Hausdorff measure is

$$\mathcal{H}^k(A) = \lim_{\delta \to 0} \inf \bigg\{ \sum_{i=1}^\infty \left(\frac{\mathrm{diam}(U_i)}{2}\right)^k : A \subseteq \bigcup_{i=1}^\infty U_i, U_i \subset \mathcal{M}, \mathrm{diam}(U_i) < \delta\bigg\},$$

where the infimum is over all countable covers of $A$ such that the constraints on $U_i$ for all $i$ are satisfied. Note that this construction does not require an embedding, only metric.

Defined this way, it's clear that if $\mathcal{M}=\mathbb{R}^k$, then $\mathcal{H}^k(\mathbb{B}^k)=1$, where $\mathbb{B}^k$ is the unit $k$-ball. So the relation between this and the usual Lebesgue measure $\lambda^k$ on $\mathbb{R}^k$ is

$$\lambda^k(A) = \mathrm{Vol}(\mathbb{B}^k) \mathcal{H}^k(A),$$

where $\mathrm{Vol}(\mathbb{B}^k) = \frac{\pi^{k/2}}{\Gamma(k/2+1)}$.

Different works will use other normalization conventions. Let

$$\tilde{\mathcal{H}}^k = \omega_k \mathcal{H}^k = \alpha_k \lambda^k$$

for constants $\omega_k$ and $\alpha_k = \frac{\omega_k}{\mathrm{Vol}(\mathbb{B}^k)}$.
Here are several different choices for $\omega_k$:

$\omega_k$	$\alpha_k$	$\tilde{\mathcal{H}}^k(\mathcal{M})$
$2^k$	$\frac{2^k}{\mathrm{Vol}(\mathbb{B}^k)}$	$2^k \mathcal{H}^k(\mathcal{M})$
$1$	$\frac{1}{\mathrm{Vol}(\mathbb{B}^k)}$	$\mathcal{H}^k(\mathcal{M})$
$\mathrm{Vol}(\mathbb{B}^k)$	$1$	$\mathrm{Vol}(\mathbb{B}^k) \mathcal{H}^k(\mathcal{M})$
$\mathcal{H}(\mathcal{M})^{-1}$	$\mathrm{Vol}(\mathcal{M})^{-1}$	$1$

The final choice only works for compact manifolds. It's convenient because a uniform distribution wrt this measure as a constant unit density. The 3rd is convenient because it coincides with the Lebesgue measure, and the uniform distribution on a compact manifold wrt this measure is the reciprocal of the "volume" of that manifold (e.g. this is the usual measure in directional statistics). The 1st and 2nd I think are mostly used in theoretical work.

So I'd be fine defining Hausdorff here as the 3rd definition and NormalizedHausdorff as the 4th.

[mateuszbaran]

Note that the Hausdorff measure (with the third normalization) has a nice equivalent definition for Riemannian manifolds, see https://link.springer.com/book/10.1007/978-3-7643-7480-8 (Chapter XII, Integration on manifolds, Eq. (1.1), which actually lets us compute densities in charts. volume_density in Manifolds.jl handles the density correction between Euclidean Lebesgue measure in normal charts and the Haudorff measure.

For example, suppose a user is doing Monte Carlo sampling for a parameter on the positive N-hemisphere and wants to specify a prior. They then have 2 samplers.

I guess users would prefer to not have to deal with adjusting prior distributions to what different samplers expect. Does Turing.jl currently have a mechanism for handling differences in base measures between samplers? Can we maybe wrap samplers in a decorator that adjusts densities and pretend they all have the same base measure?

[mateuszbaran]

But what if a user wants to provide a vector of N spherical coordinates and compute the density wrt the Lebesgue measure on those coordinates? Other coordinates could in principle be used; how should the user specify which Lebesgue measure is being used?

If a user works in a chart then they should specify which chart they work in. The chart then determines the in-chart Lebesgue measure. That should be sufficient I think?

[sethaxen]

Note that the Hausdorff measure (with the third normalization) has a nice equivalent definition for Riemannian manifolds, see https://link.springer.com/book/10.1007/978-3-7643-7480-8 (Chapter XII, Integration on manifolds, Eq. (1.1), which actually lets us compute densities in charts.

I don't see the explicit connection to the Hausdorff measure in Eq (1.1); What they're describing there is the volume measure, which should indeed be equivalent to the Hausdorff measure with the third normalization.

I guess users would prefer to not have to deal with adjusting prior distributions to what different samplers expect.

Agreed, I would say it's up to implementers of Monte Carlo proposal methods to specify the expected reference measure.

Does Turing.jl currently have a mechanism for handling differences in base measures between samplers?

For HMC sampling, Turing only supports parameters on manifolds for which one can have a bijective almost-everywhere map from R^n to the manifold. The choice of bijector and its Jacobian determinant implicitly define the base measure, and these have to match the implicit base measure of the Distributions implementation. Since the base measure is almost always Lebesgue in chart-free coordinates, it's mostly not been an issue. For full Turing compatibility, we would probably want to specify the "bijectors" ourselves in an extension/separate package. But in most cases we would want bijections with coordinate expansions, i.e. we have the map $f: \mathbb{R}^{m \times n} \to \mathcal{M} \times \mathcal{N}$, where $\mathcal{M}$ is the manifold we have a distribution on and $r \in \mathcal{N}$ are extra parameters we don't care about. If $\mathcal{N}$ is non-compact, then the user must also specify a prior on $r$ (but a default could be applied). A trivial example on the n-sphere is that it's usually more efficient to sample an unnormalized point $y$ and then map it to $(x=y/\lVert y \rVert, r=\lVert y \rVert)$ than to sample in spherical coordinates. Here $r$ needs a prior for the target distribution to have finite mass. This generalizes also to the Stiefel/Grassmann manifolds.

Can we maybe wrap samplers in a decorator that adjusts densities and pretend they all have the same base measure?

I think let's cross this bridge when we come to it. HMC with maps like above wouldn't require it but it would require Bijectors have some notion of parameter expansion or non-bijective distributions. Manifold HMC would be very advanced and we shouldn't plan for it yet. Manifold Gibbs proposals seem like they would be the first thing to work out, but I think we should still wait to work that out.

If a user works in a chart then they should specify which chart they work in. The chart then determines the in-chart Lebesgue measure. That should be sufficient I think?

I hadn't looked at the Atlas implementations in Manifolds for a while, but yes, they do seem sufficient now. Effectively we would want a way to decorate a distribution with an atlas; then the base measure is always Lebesgue, and the logpdf implementations could include the log-det-Jacobian of the chart (this would be the logdet of local_metric, right? is there a function for that?)

[mateuszbaran]

I don't see the explicit connection to the Hausdorff measure in Eq (1.1); What they're describing there is the volume measure, which should indeed be equivalent to the Hausdorff measure with the third normalization.

It looks like the connection is "elementary" for mathematicians studying geometric measure theory: https://mathoverflow.net/questions/380247/equality-of-hausdorff-measure-and-lebesgue-measure-on-manifolds-reference + linked questions.

I think the important part here is that on Riemannian manifolds there is a canonical choice of measure that can be constructed in different ways, so we probably should not use so many different names.

For HMC sampling, Turing only supports parameters on manifolds for which one can have a bijective almost-everywhere map from R^n to the manifold. The choice of bijector and its Jacobian determinant implicitly define the base measure, and these have to match the implicit base measure of the Distributions implementation. Since the base measure is almost always Lebesgue in chart-free coordinates, it's mostly not been an issue. For full Turing compatibility, we would probably want to specify the "bijectors" ourselves in an extension/separate package. But in most cases we would want bijections with coordinate expansions, i.e. we have the map

Manifold HMC would be very advanced and we shouldn't plan for it yet.

Somehow to me manifold HMC seems less confusing than using Euclidean HMC without charts on manifolds 😅 . Anyway, we can have both eventually.

I hadn't looked at the Atlas implementations in Manifolds for a while, but yes, they do seem sufficient now. Effectively we would want a way to decorate a distribution with an atlas; then the base measure is always Lebesgue, and the logpdf implementations could include the log-det-Jacobian of the chart (this would be the logdet of local_metric, right? is there a function for that?)

No, we don't have log-det-Jacobian function yet for charts. volume_density computes the det-Jacobian for normal (exponential) charts though.