* Richardson et al. (2019) Richardson A, Filos-Ratsikas A, Faltings B (2019) Rewarding high-quality data via influence functions. CoRR abs/1908.11598, 1908.11598
  * Richardson et al. (2020) Richardson A, Filos-Ratsikas A, Faltings B (2020) Budget-bounded incentives for federated learning. In: Yang Q, Fan L, Yu H (eds) Federated Learning - Privacy and Incentive, Lecture Notes in Computer Science, vol 12500, Springer, pp 176–188, DOI 10.1007/978-3-030-63076-8“˙13, URL https://doi.org/10.1007/978-3-030-63076-8_13
  * Rousseeuw (1987) Rousseeuw PJ (1987) Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of computational and applied mathematics 20:53–65
  * Schomm et al. (2013) Schomm F, Stahl F, Vossen G (2013) Marketplaces for data: an initial survey. ACM SIGMOD Record 42(1):15–26
  * Shah and Zhou (2015) Shah NB, Zhou D (2015) Double or nothing: Multiplicative incentive mechanisms for crowdsourcing. In: Cortes C, Lawrence ND, Lee DD, Sugiyama M, Garnett R (eds) Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7-12, 2015, Montreal, Quebec, Canada, pp 1–9, URL https://proceedings.neurips.cc/paper/2015/hash/c81e728d9d4c2f636f067f89cc14862c-Abstract.html
  * Shah and Zhou (2016) Shah NB, Zhou D (2016) No oops, you won’t do it again: Mechanisms for self-correction in crowdsourcing. In: Balcan M, Weinberger KQ (eds) Proceedings of the 33nd International Conference on Machine Learning, ICML 2016, New York City, NY, USA, June 19-24, 2016, JMLR.org, JMLR Workshop and Conference Proceedings, vol 48, pp 1–10, URL http://proceedings.mlr.press/v48/shaha16.html
  * Shah et al. (2015) Shah NB, Zhou D, Peres Y (2015) Approval voting and incentives in crowdsourcing. In: Bach FR, Blei DM (eds) Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015, JMLR.org, JMLR Workshop and Conference Proceedings, vol 37, pp 10–19, URL http://proceedings.mlr.press/v37/shaha15.html
  * Shapley (1953) Shapley LS (1953) A value for n-person games. Contributions to the Theory of Games 2:307–317
  * Shnayder et al. (2016a) Shnayder V, Agarwal A, Frongillo RM, Parkes DC (2016a) Informed truthfulness in multi-task peer prediction. In: Conitzer V, Bergemann D, Chen Y (eds) Proceedings of the 2016 ACM Conference on Economics and Computation, EC ’16, Maastricht, The Netherlands, July 24-28, 2016, ACM, pp 179–196, DOI 10.1145/2940716.2940790, URL https://doi.org/10.1145/2940716.2940790
  * Shnayder et al. (2016b) Shnayder V, Frongillo RM, Parkes DC (2016b) Measuring performance of peer prediction mechanisms using replicator dynamics. In: Kambhampati S (ed) Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, IJCAI 2016, New York, NY, USA, 9-15 July 2016, IJCAI/AAAI Press, pp 2611–2617, URL http://www.ijcai.org/Abstract/16/371
  * Sim et al. (2020) Sim RHL, Zhang Y, Chan MC, Low BKH (2020) Collaborative machine learning with incentive-aware model rewards. In: Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event, PMLR, Proceedings of Machine Learning Research, vol 119, pp 8927–8936, URL http://proceedings.mlr.press/v119/sim20a.html
  * Singer (2010) Singer Y (2010) Budget feasible mechanisms. In: 51th Annual IEEE Symposium on Foundations of Computer Science, FOCS 2010, October 23-26, 2010, Las Vegas, Nevada, USA, IEEE Computer Society, pp 765–774, DOI 10.1109/FOCS.2010.78, URL https://doi.org/10.1109/FOCS.2010.78
  * Snowflake (2021) Snowflake (2021) Snowflake data marketplace. https://www.snowflake.com/data-marketplace/, accessed: 2021-05-09
  * Stahl and Vossen (2016) Stahl F, Vossen G (2016) Data quality scores for pricing on data marketplaces. In: Nguyen NT, Trawinski B, Fujita H, Hong T (eds) Intelligent Information and Database Systems - 8th Asian Conference, ACIIDS 2016, Da Nang, Vietnam, March 14-16, 2016, Proceedings, Part I, Springer, Lecture Notes in Computer Science, vol 9621, pp 215–224, DOI 10.1007/978-3-662-49381-6“˙21, URL https://doi.org/10.1007/978-3-662-49381-6_21
  * Tang et al. (2013) Tang R, Wu H, Bao Z, Bressan S, Valduriez P (2013) The price is right - models and algorithms for pricing data. In: Decker H, Lhotská L, Link S, Basl J, Tjoa AM (eds) Database and Expert Systems Applications - 24th International Conference, DEXA 2013, Prague, Czech Republic, August 26-29, 2013\. Proceedings, Part II, Springer, Lecture Notes in Computer Science, vol 8056, pp 380–394, DOI 10.1007/978-3-642-40173-2“˙31, URL https://doi.org/10.1007/978-3-642-40173-2_31
  * Tang et al. (2014) Tang R, Amarilli A, Senellart P, Bressan S (2014) Get a sample for a discount - sampling-based XML data pricing. In: Decker H, Lhotská L, Link S, Spies M, Wagner RR (eds) Database and Expert Systems Applications - 25th International Conference, DEXA 2014, Munich, Germany, September 1-4, 2014. Proceedings, Part I, Springer, Lecture Notes in Computer Science, vol 8644, pp 20–34, DOI 10.1007/978-3-319-10073-9“˙3, URL https://doi.org/10.1007/978-3-319-10073-9_3
  * Tramèr et al. (2016) Tramèr F, Zhang F, Juels A, Reiter MK, Ristenpart T (2016) Stealing machine learning models via prediction apis. In: Holz T, Savage S (eds) 25th USENIX Security Symposium, USENIX Security 16, Austin, TX, USA, August 10-12, 2016, USENIX Association, pp 601–618, URL https://www.usenix.org/conference/usenixsecurity16/technical-sessions/presentation/tramer
  * Upadhyaya et al. (2016) Upadhyaya P, Balazinska M, Suciu D (2016) Price-optimal querying with data apis. Proc VLDB Endow 9(14):1695–1706, DOI 10.14778/3007328.3007335, URL http://www.vldb.org/pvldb/vol9/p1695-upadhyaya.pdf
  * Vaughan (2017) Vaughan JW (2017) Making better use of the crowd: How crowdsourcing can advance machine learning research. J Mach Learn Res 18:193:1–193:46, URL http://jmlr.org/papers/v18/17-234.html
  * Wang et al. (2020) Wang T, Rausch J, Zhang C, Jia R, Song D (2020) A principled approach to data valuation for federated learning. In: Yang Q, Fan L, Yu H (eds) Federated Learning - Privacy and Incentive, Lecture Notes in Computer Science, vol 12500, Springer, pp 153–167, DOI 10.1007/978-3-030-63076-8“˙11, URL https://doi.org/10.1007/978-3-030-63076-8_11
  * Yan and Procaccia (2021) Yan T, Procaccia AD (2021) If you like shapley then you’ll love the core. In: Thirty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2021, Thirty-Third Conference on Innovative Applications of Artificial Intelligence, IAAI 2021, The Eleventh Symposium on Educational Advances in Artificial Intelligence, EAAI 2021, Virtual Event, February 2-9, 2021, AAAI Press, pp 5751–5759, URL https://ojs.aaai.org/index.php/AAAI/article/view/16721
  * Yang et al. (2012) Yang D, Xue G, Fang X, Tang J (2012) Crowdsourcing to smartphones: incentive mechanism design for mobile phone sensing. In: Akan ÖB, Ekici E, Qiu L, Snoeren AC (eds) The 18th Annual International Conference on Mobile Computing and Networking, Mobicom’12, Istanbul, Turkey, August 22-26, 2012, ACM, pp 173–184, DOI 10.1145/2348543.2348567, URL https://doi.org/10.1145/2348543.2348567
  * Yoon et al. (2020) Yoon J, Arik SÖ, Pfister T (2020) Data valuation using reinforcement learning. In: Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event, PMLR, Proceedings of Machine Learning Research, vol 119, pp 10842–10851, URL http://proceedings.mlr.press/v119/yoon20a.html
  * Yu and Zhang (2017) Yu H, Zhang M (2017) Data pricing strategy based on data quality. Comput Ind Eng 112:1–10, DOI 10.1016/j.cie.2017.08.008, URL https://doi.org/10.1016/j.cie.2017.08.008
  * Yu et al. (2020a) Yu H, Liu Z, Liu Y, Chen T, Cong M, Weng X, Niyato D, Yang Q (2020a) A sustainable incentive scheme for federated learning. IEEE Intell Syst 35(4):58–69, DOI 10.1109/MIS.2020.2987774, URL https://doi.org/10.1109/MIS.2020.2987774
  * Yu et al. (2020b) Yu H, Yang K, Zhang T, Tsai Y, Ho T, Jin Y (2020b) Cloudleak: Large-scale deep learning models stealing through adversarial examples. In: 27th Annual Network and Distributed System Security Symposium, NDSS 2020, San Diego, California, USA, February 23-26, 2020, The Internet Society, URL https://www.ndss-symposium.org/ndss-paper/cloudleak-large-scale-deep-learning-models-stealing-through-adversarial-examples/
  * Zhang and Beltran (2020) Zhang M, Beltran F (2020) A survey of data pricing methods. SSRN
  * Zhang et al. (2020) Zhang M, Beltrán F, Liu J (2020) Selling data at an auction under privacy constraints. In: Adams RP, Gogate V (eds) Proceedings of the Thirty-Sixth Conference on Uncertainty in Artificial Intelligence, UAI 2020, virtual online, August 3-6, 2020, AUAI Press, Proceedings of Machine Learning Research, vol 124, pp 669–678, URL http://proceedings.mlr.press/v124/zhang20b.html
  * Zhang et al. (2016) Zhang X, Yang Z, Sun W, Liu Y, Tang S, Xing K, Mao X (2016) Incentives for mobile crowd sensing: A survey. IEEE Commun Surv Tutorials 18(1):54–67, DOI 10.1109/COMST.2015.2415528, URL https://doi.org/10.1109/COMST.2015.2415528
  * Zhou and Zheng (2009) Zhou X, Zheng H (2009) Trust: A general framework for truthful double spectrum auctions. In: IEEE INFOCOM 2009, IEEE, pp 999–1007
  * Zhou et al. (2014) Zhou Y, Porwal U, Zhang C, Ngo HQ, Nguyen L, Ré C, Govindaraju V (2014) Parallel feature selection inspired by group testing. In: Ghahramani Z, Welling M, Cortes C, Lawrence ND, Weinberger KQ (eds) Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada, pp 3554–3562, URL https://proceedings.neurips.cc/paper/2014/hash/fb8feff253bb6c834deb61ec76baa893-Abstract.html


The Base Measure Problem and its Solution

Alexey Radul Boris Alexeev

Google Reseach Google Research

###### Abstract

Probabilistic programming systems generally compute with probability density
functions, leaving the base measure of each such function implicit. This
mostly works, but creates problems when densities with respect to different
base measures are accidentally combined or compared. Mistakes also happen when
computing volume corrections for continuous changes of variables, which in
general depend on the support measure. We motivate and clarify the problem in
the context of a composable library of probability distributions and bijective
transformations. We solve the problem by standardizing on Hausdorff measure as
a base, and deriving formulas for comparing and combining mixed-dimension
densities, as well as updating densities with respect to Hausdorff measure
under diffeomorphic transformations. We also propose a software architecture
that implements these formulas efficiently in the common case. We hope that by
adopting our solution, probabilistic programming systems can become more
robust and general, and make a broader class of models accessible to
practitioners.

##  1 Introduction

Suppose we are designing a composable library of software to represent
probability distributions and transformations thereof which, following The TFP
Team, (2019), we will call “bijectors”. TensorFlow Probability (The TFP Team,,
2019) and PyTorch Distributions (Paszke et al.,, 2019) are such libraries in
their own right, targeting the probabilistic machine learning space. General-
purpose probabilistic programming languages such as Stan (Carpenter et al.,,
2017), BLOG (Li and Russell,, 2013), Anglican (Tolpin et al.,, 2016), or
Venture (Mansinghka et al.,, 2014) must of necessity also include such
libraries, to implement their primitive distributions and deterministic
functions.

Suppose furthermore that we are designing this library to operate on
explicitly vector-valued probability distributions.111No distinction need be
made for our purposes between vectors, matrices, tuples, or other structures,
as long as joint distributions over non-trivial powers of $\mathbb{R}$ are in
scope. This is the case for TensorFlow Probability and PyTorch Distributions,
for instance, to take advantage of vectorized hardware; and studying the
general case is instructive even for a scalar design.

In this setting, it’s conventional to represent a probability distribution $P$
as a function $p:\mathbb{R}^{n}\to\mathbb{R}$ computing the probability
density of $P$ at points in $\mathbb{R}^{n}$, together with a sampler $x\sim
s()$ drawing random variates $x\in\mathbb{R}^{n}$ distributed according to
$P$. Given a map $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, the pushforward $fP$ is
then sampled as $x^{\prime}\sim f(s())$, and its density is computed as

$fp(x^{\prime})=\frac{1}{|\det(Jf_{x})|}p(x),\ \ \mathrm{for}\
x=f^{-1}(x^{\prime}).$ (1)

The Jacobian-determinant correction $1/|\det(Jf_{x})|$ accounts for the
possibility that $f$ changes the volume of an infinitesimal volume element
near $x$. The Jacobian determinant can be computed by forming the Jacobian of
$f$, for example with automatic differentiation; but for many functions $f$,
it’s available more efficiently. Thus a conventional choice is to package such
$f$, together with their inverse $f^{-1}$, in a Bijector class with a method
for computing said Jacobian determinant.

This conventional architecture admits a serious bug. We name this bug the
_Base Measure Problem_ , because it consists of neglecting the base measure
with respect to which we are computing densities. Our contributions answer
these questions:

  * •

What’s the problem? We give a clear and intuitive example of the Base Measure
Problem in Section 2;

  * •

How common is it? We briefly survey several areas where the Base Measure
Problem recurs in different guises in Section 3;

  * •

What’s the right answer? We propose a more nuanced standard base measure in
Section 4, and derive complete density manipulation formulas for that choice
using standard results; and

  * •

How do we fix our software? We detail common special cases of our formulas in
Section 5, as well as how to arrange a software library to optimize them.

We stress that while we do propose to explicitly represent information about
measures, all the information we will need will be local to a single point
(the point itself and various directional derivatives thereat), and thus
require no symbolic algebra to compute with.

##  2 Motivating Example

Consider the uniform distribution $P$ on the unit circle in $\mathbb{R}^{2}$.
The natural probability density function to write down for this is

$p(x,y)=\begin{cases}\frac{1}{2\pi}&\mathrm{when}\ x^{2}+y^{2}=1,\\\
0&\mathrm{otherwise}.\end{cases}$

Of course, if we were sticklers we would note that the base measure implied by
the type of the samples is Lebesgue measure on $\mathbb{R}^{2}$; and with
respect to this base measure the density of $P$ is $+\infty$ on points on the
unit circle. But that’s clearly less helpful to our users than $1/2\pi$, and
$1/2\pi$ is _a_ density for $P$, just with respect to222Technically, the
density function is the Radon–Nikodym derivative of our probability measure
with respect to the base measure. Lebesgue measure along the circle. So let’s
wing it and go with that.

Figure 1: 3000 samples from the uniform distribution on the unit circle
(blue), and the same samples scaled non-isotropically by 2x along the $x$ axis
and 20x along the $y$ axis (orange). The transformed distribution is clearly
not uniform, even though the original distribution is, and the Jacobian of the
transformation is constant over $\mathbb{R}^{2}$.

Now, a user looking at a density computed for a sample has no way to know
which base measure we meant, so we have created an instance of the base
measure problem. To see how it bites us, let $f$ be a somewhat contrived non-
isotropic scaling of $\mathbb{R}^{2}$, given by $f(x,y)=(2x,20y)$. What
happens when we try to compose our uniform distribution $P$ with our non-
isotropic scaling $f$? The sampler is fine, but the conventional density rule
(1) gives

$\displaystyle fp_{\mathrm{bad}}(x^{\prime},y^{\prime})$
$\displaystyle=\frac{1}{40}p(f^{-1}(x^{\prime},y^{\prime}))$
$\displaystyle=\begin{cases}\frac{1}{80\pi}&\mathrm{when}\
(x^{\prime}/2)^{2}+(y^{\prime}/20)^{2}=1,\\\ 0&\mathrm{otherwise}.\end{cases}$

Now we have definitely made a mistake. First of all, $fp_{\mathrm{bad}}$
doesn’t integrate to 1, because the perimeter of the ellipse
$(x^{\prime}/2)^{2}+(y^{\prime}/20)^{2}=1$ is approximately $81.28$, which is
considerably less than $80\pi\approx 251.32$. Second, as we can see by drawing
a few samples and plotting them in Figure 1, the true distribution $fp$ isn’t
uniform! It’s clearly denser near $(0,20)$ than $(2,0)$.333With respect to arc
length as the base measure.

###  2.1 What went wrong?

The problem is that when computing the volume change induced by the change of
variables $f$, we forgot that the base measure for $P$ wasn’t Lebesgue on
$\mathbb{R}^{2}$. It’s actually Lebesgue along the circle, and $f$ changes arc
length differently at different points. Indeed, locally near the point
$(x,y)$, the circle is the line in the direction $(-y,x)$. The directional
derivative of $f$ is $(-2y,20x)$, and the change in arclength that $f$ induces
is therefore $\sqrt{4y^{2}+400x^{2}}$. Hence the correct density is

$\displaystyle fp_{\mathrm{good}}(x^{\prime},y^{\prime})=\ $ $\displaystyle
1/2\pi\sqrt{y^{\prime 2}/100+100x^{\prime 2}}$ $\displaystyle\mathrm{when}\
(x^{\prime}/2)^{2}+(y^{\prime}/20)^{2}=1,\ \mathrm{and}$ $\displaystyle 0\
\mathrm{otherwise}.$

##  3 How common is this problem?

While we chose to introduce the base measure problem on a continuous example,
it actually occurs most often when transforming discrete distributions
embedded in $\mathbb{R}^{n}$. The density function of, say, the Bernoulli
distribution is $1/2$ at $0$ or $1$ and $0$ elsewhere, but this is a density
with respect to counting measure, not Lebesgue measure on $\mathbb{R}$.
Therefore, when transforming this distribution with a bijector, we should not
apply any Jacobian correction, because all bijections are counting-volume-
preserving.

The same problem also shows up when dealing with distributions on symmetric
matrices, such as the Wishart or LKJ distributions. The $k\times k$ symmetric
matrices occupy a lower-dimensional sub-manifold of $\mathbb{R}^{k\times k}$,
and the density of the Wishart distribution is with respect to Lebesgue
measure on that submanifold rather than all of $\mathbb{R}^{k\times k}$.
Bijections will in general change the volume under that measure differently
than the volume under Lebesgue measure on $\mathbb{R}^{k\times k}$.

Sometimes changes of variables occur in the inference algorithm rather than
the model. For example, the Jacobian of the deterministic transformation $h$
in the reversible-jump MCMC framework (Green and Hastie,, 2009) is due to
interpreting $h$ as a change of variables. One must therefore compute this
density correction with respect to the proper base measure whenever jumping
between submanifolds of $\mathbb{R}^{n}$. For example, such jumps would occur
in model selection involving a model with a latent von Mises-Fisher-
distributed random variable (which is defined on the unit sphere).

The same problem can _also_ show up without any changes of variables at all.
It suffices to compare the density of different points under the same
distribution, or the density of the same point under different distributions.
For instance, the Indian GPA problem discussed in e.g. Wu et al., (2018) boils
down to treating a density under counting measure as comparable to a density
under Lebesgue measure on $\mathbb{R}$, even though the former represents an
infinitely larger mass of probability. This incarnation of the base measure
problem is more widely studied in the literature, albeit without a name. Our
solution smoothly covers this scenario, as we will see in Section 4.3.

##  4 What’s the right answer?

The root of the base measure problem is that we didn’t want to compute with
measures directly, but lost the base measure when representing probability
distributions with densities. It’s not actually possible to infer the correct
base measure from the data type representing the sample: a point on the unit
circle in $\mathbb{R}^{2}$ is represented with two floating-point numbers, but
using Lebesgue measure on $\mathbb{R}^{2}$ as the base is not helpful.

We propose standardizing on Hausdorff measure as a universal base measure for
probability densities. Specifically, for $P$ a distribution on a
$d$-dimensional manifold $\mathcal{M}$ embedded in $\mathbb{R}^{n}$, let’s
define the density to be with respect to the $d$-dimensional Hausdorff measure
$H^{d}$. We propose $d$-Hausdorff measure because it provides a coherent
definition for “$d$-dimensional volume” of a surface embedded in
$\mathbb{R}^{n}$, even when the surface is curved.444We propose standardizing
at all in order to minimize the amount of information carried by the base
measure, so it can be as implicit as possible. The only degree of freedom in
Hausdorff measure is the dimension, and that is unavoidable. We scale $H^{d}$
to agree with the standard volume of straight surfaces, and thus with Lebesgue
measure when $d=n$. We begin with fixed-dimension distributions for
simplicity, and briefly address mixed-dimension distributions in Section 4.2.

###  4.1 Manifolds of dimension $d$

What do we need to represent about a sample $x$ besides its density $p(x)$ in
order to compute in this scheme? Clearly we need to represent $d$; but, as the
example of the unit circle shows, that’s not enough to transform densities
correctly. Recall that our scaling transformation had a constant Jacobian when
viewed as a function from $\mathbb{R}^{2}$ to $\mathbb{R}^{2}$, but we have to
apply a position-dependent correction to account for its effect on the unit
circle density. The additional thing we need to represent is the tangent space
of $\mathcal{M}$ at the sample point $x$. Notably, we do not need to
computationally represent any non-local information about $\mathcal{M}$ or
$P$—just the tangent space at $x$, its dimensionality $d$, and the density
$p(x)$.

$\mathcal{M}$$f$$f\mathcal{M}$$x$$\delta v_{0}$$\delta v_{1}$$f(x)$$\delta
v_{0}^{\prime}$$\delta v_{1}^{\prime}$ Figure 2: Change of local 2-volume
under a transformation from $\mathbb{R}^{3}$ to $\mathbb{R}^{3}$. On the left,
the small parallelogram with sides $\delta v_{0}$ and $\delta v_{1}$ is
tangent to $\mathcal{M}$ at $x$ and forms a volume element. On the right, $f$
takes $\mathcal{M}$ to $f\mathcal{M}$, $x$ to $f(x)$, and the volume element
to the small parallelogram tangent to $f\mathcal{M}$ at $f(x)$ with sides
$\delta v_{0}^{\prime}$ and $\delta v_{1}^{\prime}$. When computing the
density at $f(x)$ of the pushforward distribution $fP$ supported on
$f\mathcal{M}$, we have to correct for the change in volume of this
parallelogram, but we have to disregard any stretching or compression $f$ may
do in the perpendicular direction. Theorem 1 formalizes this for arbitrary
dimensions.

Let us now formalize what we are actually trying to compute and how we can
compute it.

######  Definition 1.

Given a probability distribution $P$ over any space $\Omega$, and any function
$f:\Omega\to\Omega^{\prime}$, the _pushforward_ distribution $fP$ is the
measure

$fP(S)=P(f^{-1}(S))$

for all events $S\subset\Omega^{\prime}$ for which $f^{-1}(S)$ is measurable
in $\Omega$.

The sampler for $fP$ is just to apply $f$ to a sample from $P$. We seek a
tractable density function for $fP$, and fortunately the following theorem
gives it to us for nice $f$.

######  Theorem 1.

Consider

  * •

A probability distribution $P$ over a $d$-dimensional manifold $\mathcal{M}$
in $\mathbb{R}^{n}$,

  * •

with density $p$ with respect to $d$-dimensional Hausdorff measure $H^{d}$;

  * •

a diffeomorphism
$f:\mathbb{R}^{n}\to(f\mathbb{R}^{n}\subseteq\mathbb{R}^{m})$; and

  * •

a point $x\in\mathcal{M}$.

  * •

Let $v_{i}$ be an arbitrary basis for the tangent space to $\mathcal{M}$ at
$x$, and pack it into a $d$-by-$n$ matrix $V$.

Then

  1. 1.

The pushforward $fP$ is supported on the pushforward manifold $f\mathcal{M}$
in $\mathbb{R}^{m}$.

  2. 2.

The directional derivatives $v_{i}^{\prime}=df(x+\varepsilon
v_{i})/d\varepsilon$ are a basis for the tangent space to $f\mathcal{M}$ at
$f(x)$. Let us pack them into a $d$-by-$m$ matrix $V^{\prime}$.

  3. 3.

The density $fp$ of $fP$ with respect to $H^{d}$ at $x^{\prime}=f(x)$ is

$fp(x^{\prime})=p(x)\frac{\sqrt{\det(VV^{T})}}{\sqrt{\det(V^{\prime}V^{\prime
T})}}.$ (2)

Claims 1 and 2 are standard, and give us the rule for propagating a basis for
the tangent space at $x$ to $x^{\prime}$. Claim 3 amounts to a restatement in
probability terms of standard notions of change of volume. The intuition is
that the relevant volume is volume in the tangent space to $\mathcal{M}$ at
$x$ (see Figure 2). We choose as a volume element a small parallelepiped:
centered at $x$ and with sides $\delta v_{i}$ for $\delta\to 0$.555Since the
$v_{i}$ are vectors in the tangent space to $\mathcal{M}$ at $x$, it would be
more natural to let $x$ be the vertex of the parallelepiped. However, we want
to end up with a neighborhood of $x$, so we shift the parallelepiped to have
$x$ in the center. The volume of this element is
$\delta^{d}\sqrt{\det(VV^{T})}$. The function $f$ transforms this volume
element to the parallelepiped with center $x^{\prime}$ and sides $\delta
v^{\prime}_{i}$. Therefore, the volume correction we are looking for is
$\sqrt{\det(V^{\prime}V^{\prime T})}/\sqrt{\det(VV^{T})}$. We reproduce a
formal proof of Claim 3 in Appendix A.

As an aside, Theorem 1 applies without modification to discrete spaces: $d=0$,
0-Hausdorff measure is, up to scaling, identical with the usual counting
measure, and we are free to define all functions of a discrete domain to be
vacuously diffeomorphic. The matrices $V$ and $V^{\prime}$ have zero rows, and
we recover the standard result that $fp(x^{\prime})=p(x)$ for $p$ a
probability mass function on a discrete space.

Observe that the inputs required for Theorem 1 are all local: the point $x$, a
basis $v_{i}$ for the tangent space to $\mathcal{M}$ at $x$, the density
$p(x)$ of $P$ at $x$, and the ability to compute the value and directional
derivatives of $f$ at $x$. We do not need any non-local information about $P,\
\mathcal{M}$, or $f$, except that $f$ is invertible. Observe also that the
outputs given by Theorem 1 give enough information about $fP$ at
$x^{\prime}=f(x)$ to apply the same theorem again with a new diffeomorphism
$g:\mathbb{R}^{m}\to(g\mathbb{R}^{m}\subseteq\mathbb{R}^{k})$. We invite the
interested reader to verify that pushing distributions forward by Theorem 1
commutes with function composition, that is, the double pushforward $g(fP)$
computes the same densities as a single pushforward by the composition
$(g\circ f)P$. This equivalence lends itself to modularity in software,
permitting, for example, the distribution $g(fP)$ to be represented taking
advantage of special structure in $f$ or $g$, without requiring the user to
derive that special structure for $g\circ f$.

###  4.2 Mixed-dimension manifolds

Theorem 1 has given us a formula for transforming densities of probability
distributions on $d$-dimensional manifolds embedded in $\mathbb{R}^{n}$, and
justifies standardizing on $d$-dimensional Hausdorff measure as a universal
base measure for densities of such distributions.

The results generalize directly to finite unions of manifolds, of potentially
different dimensions. Namely, let $P$ be a sum of measures $P_{j}$ on
manifolds $\mathcal{M}_{j}$ of dimensions $d_{j}$, where without loss of
generality we require that the intersection of any $\mathcal{M}_{j}$ and
$\mathcal{M}_{k}$ have strictly lower dimension than at least one of $d_{j}$
and $d_{k}$. In this situation, we use the sum of the corresponding Hausdorff
measures $H^{d_{j}}$ on those manifolds as the base measure. Then the density
$p(x)$ representing $P$ is given by the density $p_{j}(x)$ of the $P_{j}$
corresponding to the lowest-dimensional manifold $\mathcal{M}_{j}$ that $x$ is
on. If in addition to $\mathcal{M}_{j}$, $x$ is on manifolds of higher
dimension, any densities $p_{j^{\prime}}$ corresponding to measures on those
higher-dimensional manifolds become zero with respect to $H^{d_{j}}$.

If $x$ occurs on multiple support manifolds of minimal dimension, then by
assumption $x$ is part of a lower-dimensional intersection, whose total mass
under $P$ is necessarily zero. We are therefore formally free to assign
whatever probability density we wish to $x$. We conjecture that reporting the
sum $\sum_{\\{j|d_{j}\mathrm{minimal}\\}}P_{j}(x)$ will follow the principle
of least surprise, but are prepared to be corrected by future experience in
probabilistic programming.

###  4.3 Multiple points

When doing probabilistic inference, we concern ourselves with more than just
one point $x$. We must compare the probability densities of different points
to assess which are the best (and by how much), and we must add the
probability densities of different points to form an empirical mixture. In
MCMC we do the former when computing the acceptance ratio; and in SMC and
likelihood weighting we do both when we resample or form the self-normalizing
estimate of the evidence. These comparisons and additions are a priori
vulnerable to the base measure problem; so to fix them, let us derive measure-
aware formulas.

When computing a probability density ratio $p(y)/p(x)$, we are actually
computing the limiting ratio of finite probability masses

$\displaystyle\operatorname{Rat}_{P}(y,x)=\lim_{\delta\to
0}P(B^{\delta}_{y})/P(B^{\delta}_{x}),$ (3)

where $B^{\delta}_{x}$ denotes the radius-$\delta$ ball about $x$.666For this
to be well-defined, $x$ and $y$ must live in a single metric space. For
purposes of computational probability, it suffices to assume we are operating
in a disjoint sum of as many copies of $\mathbb{R}^{n}$, for as many values of
$n$, as needed. We use the Euclidean distance as a metric within each
$\mathbb{R}^{n}$, and posit that the distance between points in different
$\mathbb{R}^{n}$ is infinite. As $\delta$ becomes small, $P(B^{\delta}_{x})$
becomes $Cp(x)\delta^{d_{x}}$, where $d_{x}$ is the (Hausdorff) dimension of
the support of $P$ at $x$, and $C$ is a constant that depends on the metric
but not on $\delta$ or $x$. Our limit is therefore

$\displaystyle\operatorname{Rat}_{P}(y,x)=\lim_{\delta\to
0}\delta^{d_{y}-d_{x}}p(y)/p(x).$ (4)

In the common case where $d_{x}=d_{y}$, this of course reduces to $p(y)/p(x)$;
but in the general case, (4) becomes

$\operatorname{Rat}_{P}(y,x)=\begin{cases}p(y)/p(x)&\mathrm{when}\
d_{y}=d_{x};\\\ 0&\mathrm{when}\ d_{y}>d_{x},p(y)<\infty,\\\ &\mathrm{and}\
p(x)>0;\\\ \infty&\mathrm{when}\ d_{y}<d_{x},p(x)<\infty,\\\ &\mathrm{and}\
p(y)>0;\\\ \mathrm{undefined}&\mathrm{otherwise}.\end{cases}$

Note that in the “otherwise” case, the limit (3) may well exist, but the
quantities $d_{x}$, $d_{y}$, $p(x)$, and $p(y)$ do not suffice to determine
it; so for computational purposes, we throw up our hands and say “undefined”.

Summation is more subtle, because it is not dimensionless. Let us define the
two-output function

$\displaystyle\operatorname{Sum}_{P}(x,y)=d,\lim_{\delta\to
0}\frac{P(B_{x}^{\delta})+P(B_{y}^{\delta})}{C\delta^{d}},$ (5)

where $d$ is chosen to make the limit finite and positive.777This $d$ is
unique when it exists. Then, using the same notation for $d_{x}$ and $d_{y}$,
we derive

$\operatorname{Sum}_{P}(y,x)=\begin{cases}d_{x},p(x)+p(y)&\mathrm{when}\
d_{x}=d_{y};\\\ d_{x},p(x)&\mathrm{when}\ d_{x}<d_{y},p(y)<\infty,\\\
&\mathrm{and}\ p(x)>0;\\\ d_{y},p(y)&\mathrm{when}\
d_{x}>d_{y},p(x)<\infty,\\\ &\mathrm{and}\ p(y)>0;\\\
\mathrm{undefined}&\mathrm{otherwise}.\end{cases}$

Again, we let the “otherwise” case be “undefined”, because while there may
still be a $d$ for which the limit in (5) exists, we cannot compute it from
just $d_{x}$, $d_{y}$, $p(x)$, and $p(y)$.

For completeness, we mention multiplying probabilities, even though it is
straightforward:

$\operatorname{Prod}_{P}(x,y)=d_{x}+d_{y},p(x)p(y).$

Note that to compute any of $\operatorname{Rat}_{P}$,
$\operatorname{Sum}_{P}$, or $\operatorname{Prod}_{P}$, the only information
we need about the measure $P$ is local: the density $p(\cdot)$ and support
dimension $d_{(\cdot)}$ at the points $x$ and $y$. The support dimension $d$
tells us which Hausdorff measure $H^{d}$ serves as the base measure for the
density $p(\cdot)$, and that’s enough to add, multiply, and divide, even when
$P$ spans manifolds of many dimensionalities in different places. This $d$ is
the number of basis vectors spaning the tangent space required to compute
volume corrections per Section 4.1; but for $\operatorname{Rat}_{P}$,
$\operatorname{Sum}_{P}$, and $\operatorname{Prod}_{P}$ we do not actually
need the basis itself.

##  5 How do we fix the software?

We discuss how to adjust current probabilistic programming software to fix the
base measure problem. The largest change is needed to enable correct changes
of variables (Section 5.1); given that, updating common inference algorithms
is a local adjustment, which we explicate for completeness in Section 5.2.

###  5.1 Correct changes of variables

Correctly updating densities under changes of variables requires two new
architectural features of probabilistic programming software: First, querying
a probability distribution $P$ at a point $x$ must produce a local measure,
which includes the density $p(x)$, the local dimension $d_{x}$, and a
representation of the tangent space at $x$ to the support manifold
$\mathcal{M}$. Second, changes of variable induced by deterministic
transformations must take the tangent space into account. The current standard
architecture can be seen as the specialization of the above to the case when
the dimension $d$ is fixed, and the support manifold $\mathcal{M}$ is the
entire host space $\mathbb{R}^{n}$.

To elaborate, we know from Section 4.1 that to compute the density of a
pushforward $fP$ at a point $x^{\prime}$ with respect to $d$-Hausdorff measure
it suffices to compute

  * •

The preimage point $x=f^{-1}(x^{\prime})$;

  * •

The density $p(x)$ of $P$ at $x$;

  * •

A basis $v_{i}$ for the tangent space to the support $\mathcal{M}$ of $P$ at
$x$;

  * •

The directional derivatives $v_{i}^{\prime}$ of $f$ at $x$ in directions
$v_{i}$, forming a basis for the tangent space at $x^{\prime}$ to the support
of $fP$; and

  * •

The correction term $\sqrt{\det(VV^{T})}/\sqrt{\det(V^{\prime}V^{\prime T})}$
from (2).

  * •

For composition, we also need to return the pushforward tangent basis
$v_{i}^{\prime}$

At first glance this may seem prohibitively expensive: even if we have an
automatic differentiation system ready to hand that lets us compute the
needful directional derivatives of $f$, do we really have to do that, and then
perform two full matrix multiplies and determinants, for every probability
density evaluation?

Fortunately, there are several special cases we can take advantage of to save
work, and to recover the performance of the conventional architecture when it
gives the correct answer:

  * •

If our distribution $P$ is actually discrete, the tangent space is
0-dimensional, $VV^{T}$ is the 0x0 matrix, and its determinant is vacuously 1.

  * •

If $P$ is supported on all of $\mathbb{R}^{n}$, the tangent space is also
$\mathbb{R}^{n}$, and we may take $v_{i}$ to be the standard basis, so the
$\sqrt{\det(VV^{T})}$ term is unity and we need not explicitly compute it.

  * •

If in addition $f$ maps $\mathbb{R}^{n}$ to $\mathbb{R}^{n}$ and not to
$\mathbb{R}^{m}$ for $m>n$, then the matrix $V^{\prime}$ is the Jacobian of
$f$ at $x$ and itself has a determinant. We can save a matrix multiply and a
square root because in this case $\sqrt{\det(V^{\prime}V^{\prime T})}=|\det
V^{\prime}|$. This is where we recover the conventional Jacobian-determinant
correction to probability densities.

  * •

Distributions on structured-sparse matrices, such as lower-triangular
matrices, are supported on manifolds with axis-aligned tangent spaces. In such
a case, we need only keep track of a mask representing the present dimensions,
taking the basis $v_{i}$ to be a subset of the standard basis. The
$\sqrt{\det(VV^{T})}$ term again disappears, though if $f$ does not preserve
the sparsity, the $\sqrt{\det(V^{\prime}V^{\prime T})}$ term may need to be
computed.

  * •

If in addition $f$ is a univariate transformation applied coordinate-wise, the
sparse structure will be preserved, and the correction will just be the
product of partial derivatives of $f$, $\sqrt{\det(V^{\prime}V^{\prime
T})}=\prod_{i}\frac{\partial}{\partial x_{i}}f(x)$, with $i$ here ranging over
the dimensions present in the manifold but not the others.

  * •

Finally, simplices and symmetric matrices may also merit special treatment,
because the basis $v_{i}$ can again be implicit.

This list of special cases suggests ad-hoc polymorphism as a representation
strategy. Specifically, a typical probabilistic programming system probably
already has a Bijector class (hierarchy) for representing transformations $f$.
To implement this strategy for volume corrections, we can add a TangentSpace
class hierarchy for the tangent space to $\mathcal{M}$ at $x$. This second
hierarchy can have dedicated classes for efficient cases, such as the zero
space tangent to a discrete support, the full space tangent to an
$n$-dimensional support in $\mathbb{R}^{n}$, and so on. To take advantage of
special cases based both on the tangent space and the bijector, we can define
the density correction function with two-argument dispatch, or a Visitor-like
pattern (Gamma et al.,, 1995) to emulate it in languages where only single-
argument dispatch is available.

When working with distributions on high-dimensional spaces, or with large
batches of distributions, we expect the cost of the dispatch to be small
compared with the cost of the linear algebra the dispatch avoids. Further, on
a tracing platform like JAX (Bradbury et al.,, 2019), TensorFlow (Abadi et
al.,, 2015), or TorchScript (PyTorch Contributors,, 2018), the dispatch only
occurs once during tracing, avoiding repeated linear algebra during execution.

###  5.2 Measure-aware inference algorithms

Once a probabilistic programming system correctly propagates local measures
per Section 5.1, defending standard inference algorithms against the base
measure problem is just a matter of using the formulas from Section 4.3 for
operations that compare or combine probability densities. We note again that
only the density and local dimension are needed here; the tangent basis is
only used for changing variables.

We spell out two standard inference algorithms explicitly. To resample in SMC,
just throw out particles of non-minimally-dimensioned weights, as in Algorithm
1. In the one-stage case this reduces to the dimensionality-aware likelihood
weighting algorithm of Wu et al., (2018).

Algorithm 1 Measure-aware resampling for SMC

Input: Particles $p_{j}$, dimensioned weights $(d_{j},w_{j})$

Input: Desired number of new particles $N$.

Compute $d=\min d_{j}$.

if All $w_{j}=0$ where $d_{j}=d$ then

Error: The resampling is undefined.

end if

if Any $w_{j}=+\infty$ where $d_{j}>d$ then

Error: The resampling is undefined.

end if

for $i=1$ to $N$ do

Sample index $k_{i}\in\\{j|d_{j}=d\\}$ proportional to $w_{j}$.

Output $p_{k_{i}}$.

end for

Return Dimensioned self-normalization estimate $(d,\sum_{j|d_{j}=d}w_{j})$.

To compute the Metropolis-Hastings acceptance or rejection for MCMC, just
compare the target density dimension-major, as in Algorithm 2. Note that the
proposed state $y$ must live on a same-dimension manifold as the current state
$x$, otherwise the proposal $q$ is not reversible. This can be arranged with
the usual reversible-jump MCMC technique from Green and Hastie, (2009).
However, for general probabilitic programs $\pi$, $x$ and $y$ may yield
differently-dimensioned likelihoods, and M-H must account for this.

Algorithm 2 Measure-aware Metropolis-Hastings

Input: Target $\pi(\cdot)$, proposal $q(\cdot|\cdot)$, current state $x$

Sample proposed state $y\sim q(\cdot|x)$.

Evaluate dimensioned probabilities

$(d_{x},p_{x})=\pi(x)$

$(d_{y},p_{y})=\pi(y)$

if $d_{y}>d_{x}$ then

if $p_{x}\neq 0$ and $p_{y}<+\infty$ then

Reject: $d_{x},p_{x}$ dominates $d_{y},p_{y}$.

else

Error: The comparison is undefined.

end if

end if

if $d_{y}<d_{x}$ then

if $p_{y}\neq 0$ and $p_{x}<+\infty$ then

Accept: $d_{y},p_{y}$ dominates $d_{x},p_{x}$.

else

Error: The comparison is undefined.

end if

end if

Sample $u$ uniformly between 0 and 1.

Accept if $u<\frac{p_{y}q(x|y)}{p_{x}q(y|x)}$, else Reject.

For example, in the much-discussed Indian GPA problem (see Wu et al., (2018)
for one restatement), we might do MCMC on a boolean latent variable modeling
whether a student is American or Indian. The dimensioned likelihood of the
observed GPA of 4.0 for an American student would be something like $(0,0.1)$,
indicating a point-mass of size 0.1 at the maximum GPA that American high
schools report. Similarly, the dimensioned likelihood of a 4.0 GPA for an
Indian student would be $(1,0.1)$, being the density of a scalar continuous
uniform distribution of GPAs between 0 and 10. Algorithm 2 would then
correctly accept every proposal to change the latent from “Indian” to
“Ameican”, and reject every proposal to change back.

##  6 Related Work

Broadly, all general-purpose probabilistic programming systems are related to
the present work, in that they must either produce incorrect results when the
base measure problem arises (Warfield,, 2020; The Venture Team,, 2017), avoid
the problem entirely, or somehow solve it. Three specific systems are worth
mentioning:

Stan (Carpenter et al.,, 2017) is an example of a system that avoids the base
measure problem. There are two places in a Stan model where the base measure
problem might occur: in the automatic unconstraining transformations and in
applying densities to transformed parameters. The former Stan can get right
because both the base measure and unconstraining transformation are implied by
the type of the parameter being unconstrained. The base measure problem
occurs, in a sense, but the remedy is both local and internal to the
implementation of Stan. The latter Stan pushes on the user: if a user writes a
Stan model with a parameter $x$ and a transformed parameter $f(x)$, and wishes
to code a density on $f(x)$ that corresponds to the pushforward of some known
density $p(x)$, then that user must think about base measures and partial
derivatives themselves. Stan does not help; but then again, it also does not
claim it would help, so Stan doesn’t compute anything visibly “incorrect”.

BLOG discusses addressing the base measure problem (not by that name) in Wu et
al., (2018). Relative to the present work, that treatment is simpler in two
ways: they only discuss scalar random variables, where every tangent space is
either $\mathbb{R}$ or the zero vector space; and they do not explicitly
discuss transformations of random variables.

Some cases of the base measure problem are addressed by Hakaru’s
disintegration transform (Narayanan and Shan,, 2020), which generalizes
density and also needs correct base measures. Hakaru defines a restricted
language of base measures with respect to which it can symbolically compute
(unnormalized) disintegrations of s-finite measures. As this restricted set
includes discrete-continuous mixtures, it is sufficient to correctly handle
the Indian GPA problem; but for the present example of the unit circle, it
would need a richer notion of constraints and their tangent or normal spaces
(Narayanan,, 2020).

For example, the following model of the ellipse example from Section 2 is
expressible in the system from (Narayanan and Shan,, 2020) but gives an
incorrect disintegration:

    x ~ lebesgue
    y ~ lebesgue
    observe x**2 + y**2 = 1
    return (2*x, 20*y)

The issue is that the constraint restricts the distribution to a lower-
dimensional manifold $\mathcal{M}$, and the tangent space to that manifold is
not taken into account when applying the Jacobian correction due to the
transformation $f(x,y)=(2x,20y)$. This tangent space could be recovered, since
it is perpendicular to the gradient of the constraint, so we hope that this
manifestation of the base measure problem will be relatively easy to fix in
Hakaru.

There is also a line of work on formal semantics of probabilistic programming
languages, for example Borgström et al., (2013) and Staton et al., (2016).
That work taken on its own terms is generally not vulnerable to the base
measure problem at all, because the semantics are invariably in given terms of
measures, not densities. The present work can be viewed as a bridge to
implementation, providing a complete and efficient local representation of a
measure.

##  7 Future Work

A common thread in the above related work is to address base measures
dimension-wise. In other words, reduce to the scalar-oriented setting and
track one “discrete or $\mathbb{R}$” bit per scalar random variable. This is
appealing because no explicit tangent spaces are needed, since they are all
either trivial or $\mathbb{R}$. A vector-valued distribution $P$ on a
$d$-dimensional manifold embedded in $\mathbb{R}^{n}$ can then be modeled as
$n$ deterministic scalar functions $f_{i}$ from $d$ scalar random variables.
Computing the density of $P$ involves constraining the outputs of the $f_{i}$
one by one, with a density correction given by the derivative of $f_{i}$ as
appropriate. After $d$ irredundant constraints, the remaining distribution is
discrete, and no further derivatives are needed.

It would be interesting to work out whether the preceding strategy actually
works. The difficulty we predict is that, even if a joint function $f$ is
invertible as a vector-valued function, its coordinates $f_{i}$ may not be
invertible individually. In this case, treating them separately seems to
require explicitly manipulating their preimages, which may be non-trivial
discrete sets. For example, if we represent a distribution on the unit circle
as a distribution on $[0,2\pi)$ followed by the functions
$(\sin\cdot,\cos\cdot)$ then the preimage of $\sin$ doesn’t actually give us a
unique point at which to evaluate its derivative; whereas the preimage of
$\cos$, while no longer contributing a derivative per se, is still needed to
identify the correct branch of $\sin$.

##  8 Conclusion

We have named the _Base Measure Problem_ and provided a solution to it.
Implementing the solution in probabilistic programming systems should cause
negligible loss of performance for cases that were already correctly handled,
and expand the set of models in which the system can compute correct
probability densities. Implementation does carry a code complexity cost, but
that cost is minimized by using two-argument dispatch, or emulating it with a
Visitor pattern. Despite correctly accounting for measures, no non-local
information is required.

##  9 Acknowledgments

The authors thank Srinivas Vasudevan, Wynn Vonnegut, and Praveen Narayanan for
enlightening discussion. The authors also thank the anonymous reviewers for
valuable constructive criticism.

## References

  * Abadi et al., (2015) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G. S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., and Zheng, X. (2015).  TensorFlow: Large-scale machine learning on heterogeneous systems.  Software available from tensorflow.org.
  * Borgström et al., (2013) Borgström, J., Gordon, A., Greenberg, M., Margetson, J., and Van Gael, J. (2013).  Measure transformer semantics for bayesian machine learning.  Logical Methods in Computer Science, 9(3).
  * Bradbury et al., (2019) Bradbury, J., Frostig, R., Hawkins, P., Johnson, M., Leary, C., Maclaurin, D., and Wanderman-Milne, S. (2017–2019).  JAX.  Specifically the vmap functionality.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of statistical software, 76(1).
  * Gamma et al., (1995) Gamma, E., Helm, R., Johnson, R., and Vlissides, J. (1995).  Design Patterns: Elements of Reusable Object-Oriented Software.  Addison Wesley.
  * Green and Hastie, (2009) Green, P. J. and Hastie, D. I. (2009).  Reversible jump mcmc.  Genetics, 155(3):1391–1403.
  * Li and Russell, (2013) Li, L. and Russell, S. J. (2013).  The blog language reference.  tech. rep., Technical Report UCB/EECS-2013–51.
  * Mansinghka et al., (2014) Mansinghka, V., Selsam, D., and Perov, Y. (2014).  Venture: a higher-order probabilistic programming platform with programmable inference.  arXiv preprint arXiv:1404.0099.
  * Narayanan, (2020) Narayanan, P. (2020).  Personal communication.
  * Narayanan and Shan, (2020) Narayanan, P. and Shan, C.-c. (2020).  Symbolic disintegration with a variety of base measures.  ACM Trans. Program. Lang. Syst., 42(2).
  * Paszke et al., (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019).  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems, pages 8026–8037.
  * PyTorch Contributors, (2018) PyTorch Contributors (2018).  Torch script.
  * Staton et al., (2016) Staton, S., Yang, H., Wood, F., Heunen, C., and Kammar, O. (2016).  Semantics for probabilistic programming.  Proceedings of the 31st Annual ACM/IEEE Symposium on Logic in Computer Science.
  * The TFP Team, (2019) The TFP Team (2018–2019).  TensorFlow Probability.
  * The Venture Team, (2017) The Venture Team (2017).  Venture suffers from the base measure problem.  Personal communication.
  * Tolpin et al., (2016) Tolpin, D., van de Meent, J.-W., Yang, H., and Wood, F. (2016).  Design and implementation of probabilistic programming language anglican.  In Proceedings of the 28th Symposium on the Implementation and Application of Functional Programming Languages, IFL 2016, New York, NY, USA. Association for Computing Machinery.
  * Warfield, (2020) Warfield, R. (2020).  Tensorflow probability issue #761: Incorrect probabilities after transforming a quantizeddistribution.
  * Wu et al., (2018) Wu, Y., Srivastava, S., Hay, N., Du, S. S., and Russell, S. (2018).  Discrete-continuous mixtures in probabilistic programming: Generalized semantics and inference algorithms.  In International Conference on Machine Learning.

##  Appendix A Proof of Theorem 1, Claim 3

To prove Theorem 1, Claim 3 formally, we go through the measures $P$ and $fP$,
starting with the following lemma. The only technical trick is to enlarge our
parallelepipeds to open sets in $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$, so that
measuring them with $P$ and $fP$ captures all the mass near $x$ and
$x^{\prime}$, respectively, despite any curvature of $\mathcal{M}$ or
$f\mathcal{M}$. Then we will note that in the limit the only thing we care
about is the projections of those open sets back to the respective tangent
spaces, which are the parallelepipeds we started with.

######  Lemma 1.

Let $P$, $\mathcal{M}$, $p$, $x$, $v_{i}$, and $V$ be as in the statement of
Theorem 1. Let $G_{V}$ be the $d$-parallelepiped with center $x$ and sides
$v_{i}$, and let $G_{V+}$ be any $n$-parallelepiped centered at $x$ with
section $G_{V}$. Let $G_{V}^{\delta}$ and $G^{\delta}_{V+}$ be $G_{V}$ and
$G_{V+}$, respectively, scaled by $\delta$ about $x$. Then the density $p$ of
$P$ with respect to $d$-dimensional Hausdorff measure is given by

$p(x)=\frac{1}{\sqrt{\det(VV^{T})}}\lim_{\delta\to
0}P(G^{\delta}_{V+})/\delta^{d}.$

###### Proof of lemma.

The definition of $p$ being a density for $P$ with respect to $H^{d}$ is that
for all measurable subsets $B\subset\mathbb{R}^{n}$,

$P(B)=\int_{B}p\,dH^{d}.$

$P$ is only supported on $\mathcal{M}$, so we may take the integral on the
right to be over the intersection $B\cap\mathcal{M}$. We choose for $B$ the
neighborhood $G^{\delta}_{V+}$.

As $G^{\delta}_{V+}$ becomes sufficiently small, we may assume $p$ is constant
on $G^{\delta}_{V+}\cap\mathcal{M}$, so

$P(G^{\delta}_{V+})\approx p(x)H^{d}(G^{\delta}_{V+}\cap\mathcal{M}).$

As $\delta\to 0$, $\mathcal{M}$ near $x$ approaches its tangent space, so
$G^{\delta}_{V+}\cap\mathcal{M}$ becomes $G^{\delta}_{V}$. The right-hand term
becomes the $d$-dimensional Hausdorff measure of $G^{\delta}_{V}$, which given
our choice of scaling for the Hausdorff measure is just the volume thereof,
namely $\delta^{d}\sqrt{\det(VV^{T})}$. ∎

###### Proof of Claim 3.

Let $G_{V^{\prime}}$ be the $d$-parallelepiped in $\mathbb{R}^{m}$ with center
$x^{\prime}$ and sides $v^{\prime}_{i}$. Let $G_{V^{\prime}+}$ be any
$m$-parallelepiped in $\mathbb{R}^{m}$ with section $G_{V^{\prime}}$ centered
at $x^{\prime}$. Let $G^{\delta}_{V^{\prime}+}$ be $G_{V^{\prime}+}$ scaled by
$\delta$ about $x^{\prime}$. By the lemma, we have

$fp(x^{\prime})=\frac{1}{\sqrt{\det(V^{\prime}V^{\prime T})}}\lim_{\delta\to
0}fP(G^{\delta}_{V^{\prime}+})/\delta^{d}.$

The preimage $f^{-1}(G^{\delta}_{V^{\prime}+})$ approaches an
$n$-parallelepiped centered at $x$ with section $G^{\delta}_{V}$. This is
because the basis vectors $v_{i}^{\prime}$ forming $V^{\prime}$ are the
directional derivatives of $f$ in the directions given by the basis vectors
$v_{i}$ forming $V$. In directions normal to $\mathcal{M}$, we just need the
derivatives of $f$ to be finite. We may thus apply the lemma again to write

$\displaystyle\lim_{\delta\to 0}fP(G^{\delta}_{V^{\prime}+})/\delta^{d}$
$\displaystyle=\lim_{\delta\to
0}P(f^{-1}(G^{\delta}_{V^{\prime}+}))/\delta^{d}$
$\displaystyle=\sqrt{\det(VV^{T})}p(x),$

giving the desired result. ∎