# Vectorial analogues
of Cauchy’s surface area formula

Daniel Hug and Rolf Schneider

###### Abstract

Cauchy’s surface area formula says that for a convex body $K$ in
$n$-dimensional Euclidean space the mean value of the $(n-1)$-dimensional
volumes of the orthogonal projections of $K$ to hyperplanes is a constant
multiple of the surface area of $K$. We prove an analogous formula, with the
volumes of the projections replaced by their moment vectors. This requires to
introduce a new vector-valued valuation on convex bodies.
Keywords: Cauchy’s surface area formula, moment vector, tensor-valued
valuation
2020 Mathematics Subject Classification: Primary 52A20, Secondary 53C65

##  1 Introduction

Cauchy’s surface area formula is probably the first result in a branch later
named integral geometry. The $n$-dimensional version of Cauchy’s formula says
that

$S(K)=\frac{1}{\kappa_{n-1}}\int_{S^{n-1}}V_{n-1}(K|u^{\perp})\,{\rm d}u$

for each convex body $K$ (a nonempty, compact, convex set) in Euclidean space
${\mathbb{R}}^{n}$. Here $K|u^{\perp}$ is the image of $K$ under orthogonal
projection to the linear subspace $u^{\perp}$ orthogonal to the unit vector
$u$. By $S(K)$ we denote the surface area of $K$ and by $V_{k}$ the
$k$-dimensional volume. Further, $\int_{S^{n-1}}f(u)\,{\rm d}u$ always denotes
integration over the unit sphere $S^{n-1}$ with respect to the spherical
Lebesgue measure, and $\kappa_{k}$ is the volume of the $k$-dimensional unit
ball ($\omega_{k}=k\kappa_{k}$ is its surface area). For proofs of more
general projection formulas we refer to [14, p. 301] and [16, p. 222] or [6,
Chap. 5].

Classical integral geometry, as it can be found, e.g., in the books by
Hadwiger [3, Chap. 6] and Santaló [11] (see also [14, Sect. 4.4]), exhibits,
for example, intersection theorems such as the kinematic formulas and the
Crofton formulas, mainly for real-valued valuations on convex bodies. Integral
geometric formulas for vector-valued functionals were investigated in [4, 10,
12]. An extensive theory of integral geometric intersection formulas for
tensor valued valuations, also with values in spaces of measures, was
developed in [2, 5, 7, 8, 9, 13, 15]. We also refer to the survey article [1].

While the classical integral geometry of real-valued functionals also
considered projections, the mentioned investigations of tensor valuations are
restricted to intersections. It seems that not even such a simple question as
the following vector-valued analogue of Cauchy’s surface area formula has been
answered. First we recall (e.g., from [14, Sect. 5.4.1]) that for a convex
body $K\subset{\mathbb{R}}^{n}$ and for $k\in\\{{\rm dim}\,K,\dots,n\\}$ the
moment vector $z_{k+1}(K)$ is defined by

$z_{k+1}(K):=\int_{K}x\,\mathcal{H}^{k}({\rm d}x),$

where $\mathcal{H}^{k}$ is the $k$-dimensional Hausdorff measure (the lower
index indicates the degree of homogeneity). The question is then whether the
integral

$\int_{S^{n-1}}z_{n}(K|u^{\perp})\,{\rm d}u$

can be expressed by familiar functionals of $K$. It turns out that, compared
to the cases dealing with intersections, additional tensor valuations are
required. After introducing these, we shall give an answer to the question in
Theorem 1.

##  2 Preliminaries, and a result

We work in $n$-dimensional Euclidean vector space ${\mathbb{R}}^{n}$ (with
origin $o$) and use its scalar product $\langle\cdot\,,\cdot\rangle$ to
identify the space with its dual space. Thus, we identify a vector
$x\in{\mathbb{R}}^{n}$ with the linear functional $\langle x,\cdot\rangle$. By
${\mathbb{T}}^{r}$ we denote the real vector space (with its usual topology)
of symmetric $r$-tensors (tensors of rank $r\in{\mathbb{N}}_{0}$) on
${\mathbb{R}}^{n}$. The elements of ${\mathbb{T}}^{r}$ are symmetric
$r$-linear functionals on ${\mathbb{R}}^{n}$. The symmetric tensor product of
$a\in{\mathbb{T}}^{r}$ and $b\in{\mathbb{T}}^{s}$ is denoted by $a\cdot b=ab$
and is an element of ${\mathbb{T}}^{r+s}$. In particular, for
$x\in{\mathbb{R}}^{n}$ we write $x^{r}:=x\cdots x$ ($r$ factors); we then have
$x^{r}(a_{1},\dots,a_{r})=\langle x,a_{1}\rangle\cdots\langle x,a_{r}\rangle$
for $a_{1},\dots,a_{r}\in{\mathbb{R}}^{n}$.

Let $K\in{\mathcal{K}}^{n}$, where ${\mathcal{K}}^{n}$ denotes the set of
convex bodies in ${\mathbb{R}}^{n}$. Its volume can be written as

$V_{n}(K)=\int_{K}{\rm d}x,$ (1)

where $\int f(x)\,{\rm d}x$ indicates integration with respect to the Lebesgue
measure on ${\mathbb{R}}^{n}$. This motivates one to define a symmetric
$r$-tensor by

$\Psi_{r}(K):=\frac{1}{r!}\int_{K}x^{r}{\rm d}x,\quad K\in{\mathcal{K}}^{n}.$

The factor before the integral is for convenience, and we have
$\Psi_{1}(K)=z_{n+1}(K)$. For $z_{n+1}$, we get the polynomial expansion

$z_{n+1}(K+\lambda B^{n})=\sum_{j=0}^{n}\binom{n}{j}\lambda^{j}q_{j}(K)$

(see [14, (5.95)]), where $B^{n}$ is the unit ball of ${\mathbb{R}}^{n}$ and
$\lambda\geq 0$. The vector-valued coefficients $q_{j}$ can be represented by

$q_{j}(K)=\frac{1}{n}\int_{\partial K}x\,C_{n-j}(K,{\rm d}x)$

for $j\geq 1$ (whereas $q_{0}=z_{n+1}$). See [14, (5.98)], and note that the
curvature measure $C_{m}(K,\cdot)$ is concentrated on the boundary of $K$. In
particular (by [14, (4.31)]), if $K$ has nonempty interior, then

$q_{1}(K)=\frac{1}{n}\int_{\partial K}x\,\mathcal{H}^{n-1}({\rm d}x).$ (2)

Let ${\mathcal{K}}^{n}_{o}$ denote the set of convex bodies
$K\in{\mathcal{K}}^{n}$ with $o\in{\rm int}\,K$. Alternatively to (1), for
$K\in{\mathcal{K}}^{n}_{o}$ the volume of $K$ can also be written as the
integral of the cone-volume measure of $K$ over the unit sphere $S^{n-1}$. The
cone-volume measure of $K\in{\mathcal{K}}^{n}_{o}$ is defined by

$V_{K}(\omega):=\mathcal{H}^{n}\left(\bigcup_{x\in\tau(K,\omega)}[o,x]\right),\quad\omega\in\mathcal{B}(S^{n-1}).$

Here $\tau(K,\omega)$ is the set of boundary points of $K$ at which there
exists an outer unit normal vector of $K$ belonging to $\omega$, $[o,x]$ is
the closed segment with endpoints $o$ and $x$, and $\mathcal{B}(S^{n-1})$
denotes the $\sigma$-algebra of Borel sets in $S^{n-1}$. It is known (see,
e.g., [14, p. 501]) that

$\displaystyle V_{K}(\omega)$
$\displaystyle=\frac{1}{n}\int_{\omega}h_{K}(u)\,S_{n-1}(K,{\rm d}u)$
$\displaystyle=\frac{1}{n}\int_{{\mathbb{R}}^{n}\times\omega}\langle
x,u\rangle\,\Theta_{n-1}(K,{\rm d}(x,u)),$

where $h_{K}$ is the support function of $K$, $S_{n-1}(K,\cdot)$ is its
surface area measure and $\Theta_{n-1}(K,\cdot)$ is a support measure of $K$
(see [14, Sect. 4.2]). For the last equality, we refer to [14, (4.11)] and the
simple observation that $h_{K}(u)=\langle x,u\rangle$ for $(x,u)$ in the
support of the measure $\Theta_{n-1}(K,\cdot)$. From these representations it
is clear that $V_{K}$ is a measure-valued valuation (see [14, Thm. 4.2.1]).
Moreover, the integral representations are well defined for all convex bodies.

The preceding discussion suggests to define a symmetric $r$-tensor functional
by

$\Upsilon_{r}(K):=\int_{S^{n-1}}u^{r}\,V_{K}({\rm d}u),\quad
K\in{\mathcal{K}}^{d}_{o}.$

For general convex bodies $K\in{\mathcal{K}}^{d}$, it is consistent to define
$\Upsilon_{r}:{\mathcal{K}}^{n}\to{\mathbb{T}}^{r}$ by

$\displaystyle\Upsilon_{r}(K):$
$\displaystyle=\frac{1}{n}\int_{S^{n-1}}h_{K}(u)u^{r}\,S_{n-1}(K,{\rm d}u)$
(3) $\displaystyle=\frac{1}{n}\int_{{\mathbb{R}}^{n}\times S^{n-1}}\langle
x,u\rangle u^{r}\,\Theta_{n-1}(K,{\rm d}(x,u)).$

Just as the $\Psi_{r}$, also the $\Upsilon_{r}$ are rotation covariant,
continuous tensor-valued valuations. But the translation behavior is
different. We know (see [14, (5.104)]) that, for $t\in{\mathbb{R}}^{n}$,

$\Psi_{r}(K+t)=\sum_{j=0}^{r}\frac{1}{j!}\Psi_{r-j}(K)t^{j},$

where $\Psi_{r-j}(K)t^{j}$ is a symmetric tensor product. On the other hand,
if we define a translation invariant $r$-tensor by

$\Xi_{r}(K):=\frac{1}{n}\int_{S^{n-1}}u^{r}\,S_{n-1}(K,{\rm d}u),$

for $K\in{\mathcal{K}}^{n}$, then

$\Upsilon_{r}(K+t)=\Upsilon_{r}(K)+\Xi_{r+1}(K)(t),$

where the symmetric $(r+1)$-tensor $\Xi_{r+1}(K)$ is applied to the vector
$t$, which results in a tensor of rank $r$.

Considering that $\Upsilon_{0}(K)$ is the volume of $K$ and that

$\frac{1}{n}\int_{S^{n-1}}h_{K_{1}}(u)\,S(K_{2},\dots,K_{n},{\rm
d}u)=V(K_{1},K_{2},\dots,K_{n})$

(where $S(K_{2},\dots,K_{n},\cdot)$ is the mixed area measure of
$K_{2},\dots,K_{n}$) is the mixed volume of $K_{1},\dots,K_{n}$ and hence is
symmetric in its arguments, one may ask whether also

$\int_{S^{n-1}}h_{K_{1}}(u)u^{r}S(K_{2},\dots,K_{n},{\rm d}u)$ (4)

is symmetric in $K_{1},\dots,K_{n}$. That this is in general not the case, can
be seen from the fact that, for $r=1$, the vector

$\int_{S^{n-1}}h_{K}(u)u\,S(B^{n},\dots,B^{n},{\rm d}u)$

is proportional to the Steiner point of $K$ (as follows from [14, (1.31)]),
whereas

$\int_{S^{n-1}}h_{B^{n}}(u)u\,S(K,B^{n},\dots,B^{n},{\rm d}u)$

is always equal to $o$, by [14, (5.30)].

Since (4) is not symmetric in $K_{1},\ldots,K_{n}$, we introduce a symmetric
tensor functional $\Upsilon^{(r)}:({\mathcal{K}}^{n})^{n}\to{\mathbb{T}}^{r}$
by

$\Upsilon^{(r)}(K_{1},\ldots,K_{n}):=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{n}\int_{S^{n-1}}h_{K_{i}}(u)u^{r}\,S(K_{1},\dots,\check{K_{i}},\ldots,K_{n},{\rm
d}u),$

where $\check{K_{i}}$ means that $K_{i}$ is omitted.

We have mentioned tensor valuations of rank higher than necessary for the
following, since we hope to come back to them later.

We can now state the vectorial counterpart to Cauchy’s surface area formula.

######  Theorem 1.

For $K\in{\mathcal{K}}^{n}$,

$\int_{S^{n-1}}z_{n}(K|u^{\perp})\,{\rm
d}u=\frac{n\kappa_{n-1}}{n+1}(nq_{1}(K)-\Upsilon_{1}(K)).$

The functional $q_{1}$ is a special mixed moment vector, that is, we have

$\frac{n}{n+1}q_{1}(K)=z(K[n],B^{n}).$

The mixed moment vectors $z:({\mathcal{K}}^{n})^{n+1}\to{\mathbb{R}}^{n}$ are
symmetric and multilinear in each of their $n+1$ arguments (see [14, Section
5.4.1]. The choice $K=\lambda_{1}K_{1}+\dots+\lambda_{n}K_{n}$ with
$\lambda_{1},\ldots,\lambda_{n}\geq 0$ in Theorem 1, polynomial expansion and
comparison of the coefficients of $\lambda_{1}\cdots\lambda_{n}$ yields the
following consequence.

######  Corollary 1.

Let $K_{1},\ldots,K_{n}\in{\mathcal{K}}^{n}$. Then

$\displaystyle\int_{S^{n-1}}z_{n}(K_{1}|u^{\perp},\ldots,K_{n}|u^{\perp})\,{\rm
d}u$
$\displaystyle=n\kappa_{n-1}\left(z(K_{1},\ldots,K_{n},B^{n})-\frac{1}{n+1}\Upsilon^{(1)}(K_{1},\ldots,K_{n})\right).$

##  3 Proof of Theorem 1

Both sides of the asserted equation depend continuously on $K$. Hence we can
assume that $K\in{\mathcal{K}}^{n}$ has nonempty interior. Let
$\partial^{\prime}K$ denote the set of all $x\in\partial K$ at which the outer
unit normal vector $\nu_{K}(x)$ of $K$ at $x$ is unique. Then
$\partial^{\prime}K$ is a Borel set, $\nu_{K}$ is continuous on
$\partial^{\prime}K$ and $\partial K\setminus\partial^{\prime}K$ has
$\mathcal{H}^{n-1}$-measure zero. Moreover, for $y\in\partial^{\prime}K$ and
$u\in S^{n-1}$ the projection map $p_{u}:\partial K\to u^{\perp}$, $y\mapsto
y-\langle y,u\rangle u$, has the Jacobian
$Jp_{u}(y)=|\langle\nu_{K}(y),u\rangle|$. For $x\in{\rm
relint}\,(K|u^{\perp})$ we have
$p_{u}^{-1}(\\{x\\})=\\{y_{u}^{-}(K,x),y_{u}^{+}(K,x)\\}$ and
$p_{u}(y^{\pm}_{u}(K,x))=x$. Therefore, we can apply the coarea formula and
then Fubini’s theorem to obtain

$\displaystyle 2\int_{S^{n-1}}z_{n}(K|u^{\perp})\,{\rm d}u$ $\displaystyle=$
$\displaystyle 2\int_{S^{n-1}}\int_{K|u^{\perp}}x\,\mathcal{H}^{n-1}({\rm
d}x)\,{\rm d}u$ (5) $\displaystyle=$
$\displaystyle\int_{S^{n-1}}\int_{K|u^{\perp}}\int_{p_{u}^{-1}(\\{x\\})}p_{u}(y)\,\mathcal{H}^{0}({\rm
d}y)\,\mathcal{H}^{n-1}({\rm d}x)\,{\rm d}u$ $\displaystyle=$
$\displaystyle\int_{S^{n-1}}\int_{\partial
K}|\langle\nu_{K}(y),u\rangle|(y-\langle y,u\rangle u)\,\mathcal{H}^{n-1}({\rm
d}y)\,{\rm d}u$ $\displaystyle=$ $\displaystyle\int_{\partial
K}y\int_{S^{n-1}}|\langle\nu_{K}(y),u\rangle|\,{\rm
d}u\,\mathcal{H}^{n-1}({\rm d}y)$ $\displaystyle-\int_{\partial
K}\int_{S^{n-1}}|\langle\nu_{K}(y),u\rangle|\langle y,u\rangle u\,{\rm
d}u\,\mathcal{H}^{n-1}({\rm d}y).$

For the first integral in (5) we can use that

$\int_{S^{n-1}}|\langle v,u\rangle|\,{\rm d}u=2\kappa_{n-1}$

for all $v\in S^{n-1}$ and obtain

$\int_{\partial K}y\int_{S^{n-1}}|\langle\nu_{K}(y),u\rangle|\,{\rm
d}u\,\mathcal{H}^{n-1}({\rm d}y)=2n\kappa_{n-1}q_{1}(K),$

by (2). For the second integral, we need a lemma. For this, we denote by $Q$
the metric $2$-tensor on ${\mathbb{R}}^{n}$, that is, $Q(x,y)=\langle
x,y\rangle$ for $x,y\in{\mathbb{R}}^{n}$.

######  Lemma 1.

Let $n\geq 2$ and $v\in S^{n-1}$. Then

$\int_{S^{n-1}}|\langle v,u\rangle|u^{2}\,{\rm
d}u=\frac{2\kappa_{n-1}}{n+1}(v^{2}+Q).$

###### Proof.

We use a decomposition of spherical Lebesgue measure to get

$\displaystyle\int_{S^{n-1}}|\langle v,u\rangle|u^{2}\,\mathcal{H}^{n-1}({\rm
d}u)$ $\displaystyle=\int_{-1}^{1}\int_{S^{n-1}\cap
v^{\perp}}(1-\tau^{2})^{\frac{n-3}{2}}\left|\left\langle v,\tau
v+\sqrt{1-\tau^{2}}\,w\right\rangle\right|\left(\tau
v+\sqrt{1-\tau^{2}}\,w\right)^{2}\mathcal{H}^{n-2}({\rm d}w)\,{\rm d}\tau$
$\displaystyle=\int_{-1}^{1}\int_{S^{n-1}\cap
v^{\perp}}(1-\tau^{2})^{\frac{n-3}{2}}|\tau|\left(\tau^{2}v^{2}+2\tau\sqrt{1-\tau^{2}}\,vw+(1-\tau^{2})w^{2}\right)\mathcal{H}^{n-2}({\rm
d}w)\,{\rm d}\tau$
$\displaystyle=\int_{-1}^{1}|\tau|\tau^{2}(1-\tau^{2})^{\frac{n-3}{2}}{\rm
d}\tau\,\omega_{n-1}v^{2}+\int_{-1}^{1}|\tau|(1-\tau^{2})^{\frac{n-1}{2}}{\rm
d}\tau\int_{S^{n-1}\cap v^{\perp}}w^{2}\,\mathcal{H}^{n-2}({\rm d}w),$ (6)

since

$\int_{S^{n-1}\cap v^{\perp}}w\,\mathcal{H}^{n-2}({\rm d}w)=o.$

Let $B(\cdot\,,\cdot)$ denote the Beta function. Then

$\displaystyle\int_{-1}^{1}|\tau|\tau^{2}(1-\tau^{2})^{\frac{n-3}{2}}{\rm
d}\tau$ $\displaystyle=$ $\displaystyle
2\int_{0}^{1}\tau^{3}(1-\tau^{2})^{\frac{n-3}{2}}{\rm
d}\tau=\int_{0}^{1}s(1-s)^{\frac{n-3}{2}}{\rm d}s$ (7) $\displaystyle=$
$\displaystyle
B\left(2,\frac{n-1}{2}\right)=\frac{\Gamma(2)\Gamma(\frac{n-1}{2})}{\Gamma(\frac{n+3}{2})}=\frac{4}{n^{2}-1}$

and

$\int_{-1}^{1}|\tau|(1-\tau^{2})^{\frac{n-1}{2}}{\rm
d}\tau=\int_{0}^{1}(1-s)^{\frac{n-1}{2}}{\rm d}s=\frac{2}{n+1}.$

It is known from [15, (24)] that

$\int_{S^{n-1}\cap v^{\perp}}w^{2}\,\mathcal{H}^{n-2}({\rm
d}w)=2\frac{\omega_{n+1}}{\omega_{3}}Q_{v^{\perp}}=\kappa_{n-1}(Q-v^{2}),$ (8)

where $Q_{v^{\perp}}=Q-v^{2}$ is the metric tensor of $v^{\perp}$. Combination
of (6)–(8) yields the assertion. ∎

Using Lemma 1, we can write the second integral in (5) as

$\displaystyle\int_{\partial K}\int_{S^{n-1}}|\langle
u,\nu_{K}(y)\rangle|\langle y,u\rangle u\,{\rm d}u\,\mathcal{H}^{n-1}({\rm
d}y)\allowdisplaybreaks$ $\displaystyle=\int_{\partial
K}\int_{S^{n-1}}|\langle u,\nu_{K}(y)\rangle|u^{2}(y)\,{\rm
d}u\,\mathcal{H}^{n-1}({\rm d}y)\allowdisplaybreaks$
$\displaystyle=\frac{2\kappa_{n-1}}{n+1}\int_{\partial
K}(\nu_{K}(y)^{2}+Q)(y)\,\mathcal{H}^{n-1}({\rm d}y)$
$\displaystyle=\frac{2\kappa_{n-1}}{n+1}\int_{\partial
K}\langle\nu_{K}(y),y\rangle\nu_{K}(y)\,\mathcal{H}^{n-1}({\rm
d}y)+\frac{2\kappa_{n-1}}{n+1}\int_{\partial K}y\,\mathcal{H}^{n-1}({\rm d}y)$
$\displaystyle=\frac{2\kappa_{n-1}}{n+1}\int_{S^{n-1}}h_{K}(u)u\,S_{n-1}(K,{\rm
d}u)+\frac{2n\kappa_{n-1}}{n+1}q_{1}(K)$
$\displaystyle=\frac{2n\kappa_{n-1}}{n+1}\left[\Upsilon_{1}(K)+q_{1}(K)\right].$

We have used that the surface area measure $S_{n-1}(K,\cdot)$ is the image
measure of $\mathcal{H}^{n-1}$ under the Gauss map, together with the
transformation theorem for integrals, and we have again used (2), as well as
(3). Taking both integrals of (5) together, we complete the proof of Theorem
1.

##  4 Another analogue of Cauchy’s formula

We mention briefly another way of obtaining an analogue of Cauchy’s surface
area formula. For this, we consider a convex body $K\in{\mathcal{K}}^{n}_{o}$
and a measurable, bounded function $f:\partial K\to{\mathbb{T}}^{r}$. For
$x\in K|u^{\perp}$ the intersection $K\cap(x+{\mathbb{R}}u)$ is a segment
$[y_{u}^{-}(K,x),y_{u}^{+}(K,x)]$ with boundary points
$y_{u}^{-}(K,x),y_{u}^{+}(K,x)\in\partial K$ satisfying, say, $\langle
y_{u}^{-}(K,x),u\rangle\leq\langle y_{u}^{+}(K,x),u\rangle$. Defining the
$r$-tensor

$F_{f}(K,u):=\int_{K|u^{\perp}}f(y_{u}^{+}(K,x))\,\mathcal{H}^{n-1}({\rm
d}x),$

we state that

$\int_{S^{n-1}}F_{f}(K,u)\,{\rm d}u=\kappa_{n-1}\int_{\partial
K}f(x)\,\mathcal{H}^{n-1}({\rm d}x).$ (9)

For $r=0$ and $f=1$ this is Cauchy’s surface area formula. While (9) involves
integration over projections, the relation is rather a statement about real-
valued than tensor-valued functions.

For the proof of (9), we write

$\displaystyle\partial_{u}^{\pm}K$ $\displaystyle:=$
$\displaystyle\\{y_{u}^{\pm}(K,x):x\in K|u^{\perp}\\}.$

By the coarea formula we then have

$\int_{\partial_{u}^{\pm}K}f(x)|\langle\nu_{K}(x),u\rangle|\,\mathcal{H}^{n-1}({\rm
d}x)=\int_{K|u^{\perp}}f(y_{u}^{\pm}(K,x))\,\mathcal{H}^{n-1}({\rm d}x).$ (10)

Since $\langle\nu_{K}(x),u\rangle=0$ for $\mathcal{H}^{n-1}$-almost all $x\in
K\cap({\rm relbd}\,(K|u^{\perp})+{\mathbb{R}}u)$ (here ${\rm relbd}\,$ refers
to the boundary relative to $u^{\perp}$), we have

$\int_{K\cap({\rm
relbd}\,(K|u^{\perp})+{\mathbb{R}}u)}f(x)|\langle\nu_{K}(x),u\rangle|\,\mathcal{H}^{n-1}({\rm
d}x)=0.$ (11)

It follows from $y_{-u}^{+}(K,x)=y_{u}^{-}(K,x)$, (10) and (11) that

$F_{f}(K,u)+F_{f}(K,-u)=\int_{\partial
K}f(x)|\langle\nu_{K}(x),u\rangle|\,\mathcal{H}^{n-1}({\rm d}x).$

Fubini’s theorem gives

$\displaystyle\int_{S^{n-1}}F_{f}(K,u)\,{\rm d}u$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{S^{n-1}}\int_{\partial
K}f(x)|\langle\nu_{K}(x),u\rangle|\,\mathcal{H}^{n-1}({\rm d}x)\,{\rm d}u$
$\displaystyle=$ $\displaystyle\frac{1}{2}\int_{\partial
K}f(x)\int_{S^{n-1}}|\langle\nu_{K}(x),u\rangle|\,{\rm
d}u\,\mathcal{H}^{n-1}({\rm d}x)$ $\displaystyle=$
$\displaystyle\kappa_{n-1}\int_{\partial K}f(x)\,\mathcal{H}^{n-1}({\rm d}x).$

This completes the proof of (9).

We add a final remark. Tsukerman and Veomett [17], in their Section 3, aim at
an extension of Cauchy’s formula to moment vectors, but in an irritating way.
Their notion $z_{n}(K_{u})$ has different meanings in the formulation of
Theorem 3 and in its proof. According to the definitions given on page 927,
$z_{n}(K_{u})=\int_{K_{u}}x\,\mathcal{H}^{n-1}({\rm d}x)$, where $K_{u}$ is a
part of the boundary of $K$. But this contradicts the third displayed formula
on page 928, which should read

$D_{[o,u]}(z_{n+1})(K)=\lim_{\epsilon\to
0^{+}}\frac{z_{n+1}((K+\epsilon[o,u])\setminus
K)}{\epsilon}=(n+1)z(K,\dots,K,[o,u]),$

and this is different from $z_{n}(K_{u})$. With the notations used above, we
have

$(n+1)z(K,\dots,K,[o,u])=\int_{K|u^{\perp}}y_{u}^{+}(K,x)\,\mathcal{H}^{n-1}({\rm
d}x).$

Therefore, with corrections, the argument on p. 928 of [17] can be interpreted
to yield the special case $f(x)=x$ of formula (9).

Acknowledgements. D. Hug was supported by DFG research grant HU 1874/5-1 (SPP
2265).

## References

  * [1] Bernig, A., Hug, D.: Integral geometry and algebraic structures for tensor valuations. In: Kiderlen, M., Vedel Jensen, E.B. (eds.) Tensor Valuations and their Applications in Stochastic Geometry and Imaging, Lecture Notes in Mathematics 2177, pp. 79–109, Springer, Cham (2017)
  * [2] Bernig, A., Hug, D.: Kinematic formulas for tensor valuations. J. Reine Angew. Math. 736, 141–191 (2018)
  * [3] Hadwiger, H.: Vorlesungen über Inhalt, Oberfläche und Isoperimetrie. Springer, Berlin (1957)
  * [4] Hadwiger, H., Schneider, R.: Vektorielle Integralgeometrie. Elem. Math. 26, 49–57 (1971)
  * [5] Hug, D., Schneider, R., Schuster, R.: Integral geometry of tensor valuations. Adv. in Appl. Math. 41, 482–509 (2008)
  * [6] Hug, D., Weil, W.: Lectures on Convex Geometry. Springer, Cham (2020)
  * [7] Hug, D., Weis, J.A.: Crofton formulae for tensor-valued curvature measures. In: Kiderlen, M., Vedel Jensen, E.B. (eds.) Tensor Valuations and their Applications in Stochastic Geometry and Imaging, Lecture Notes in Mathematics 2177, pp. 111–156, Springer, Cham (2017)
  * [8] Hug, D., Weis, J.A.: Crofton formulae for tensorial curvature measures: the general case. In: Bianchi, G., Colesanti, A., Gronchi, P. (eds.) Analytic aspects of convexity, Springer INdAM Series 25, pp. 39–60, Springer, Cham (2018)
  * [9] Hug, D., Weis, J.A.: Kinematic formulas for tensorial curvature measures. Annali Mat. Pura Appl. 197, 1349–1384 (2018)
  * [10] Müller, H.R.: Über Momente ersten und zweiten Grades in der Integralgeometrie. Rend. Circ. Mat. Palermo, II. Ser., 2, 119–140 (1953)
  * [11] Santaló, L.A.: Integral Geometry and Geometric Probability. Addison–Wesley, Reading, Ma (1976)
  * [12] Schneider, R.: Krümmungsschwerpunkte konvexer Körper II. Abh. Math. Sem. Univ. Hamburg 37, 204–217 (1972)
  * [13] Schneider, R.: Tensor valuations on convex bodies and integral geometry. Rend. Circ. Mat. Palermo, II. Ser., 65, 295–316 (2000)
  * [14] Schneider, R.: Convex Bodies: The Brunn–Minkowski Theory. 2nd edn., Encyclopedia of Mathematics and Its Applications, vol. 151, Cambridge University Press, Cambridge (2014)
  * [15] Schneider, R., Schuster, R.: Tensor valuations on convex bodies and integral geometry, II. Rend. Circ. Mat. Palermo, II. Ser., Suppl. 70, 295–314 (2002)
  * [16] Schneider, R., Weil, W.: Stochastic and Integral Geometry. Springer, Berlin (2008)
  * [17] Tsukerman, E., Veomett, E.: Brunn–Minkowski theory and Cauchy’s surface area formula. Amer. Math. Monthly 124, 922–929 (2017)

Authors’ addresses:
Daniel Hug
Karlsruhe Institute of Technology (KIT)
Department of Mathematics
D-76128 Karlsruhe, Germany
E-mail<EMAIL_ADDRESS>
Rolf Schneider
Mathematisches Institut
Albert-Ludwigs-Universität
D-79104 Freiburg i. Br., Germany
E-mail<EMAIL_ADDRESS>


# Solutions of the Multivariate Inverse Frobenius–Perron Problem

Colin Fox  Department of Physics, University of Otago, Dunedin, New Zealand.
<EMAIL_ADDRESS>Li-Jen Hsiao  System Manufacturing Center,
National Chung-Shan Institute of Science & Technology, Taiwan.
<EMAIL_ADDRESS>Jeong Eun (Kate) Lee  Department of
Statistics, The University of Auckland, Auckland, New Zealand.
<EMAIL_ADDRESS>
###### Abstract

We address the inverse Frobenius–Perron problem: given a prescribed target
distribution $\rho$, find a deterministic map $M$ such that iterations of $M$
tend to $\rho$ in distribution. We show that all solutions may be written in
terms of a factorization that combines the forward and inverse Rosenblatt
transformations with a uniform map, that is, a map under which the uniform
distribution on the $d$-dimensional hypercube as invariant. Indeed, every
solution is equivalent to the choice of a uniform map. We motivate this
factorization via $1$-dimensional examples, and then use the factorization to
present solutions in $1$ and $2$ dimensions induced by a range of uniform
maps.

##  I Introduction

A basic question in the theory of discrete dynamical systems, and in
statistical mechanics, is whether a chaotic iterated function $M:X\rightarrow
X$, that maps a space $X\subseteq\mathbb{R}^{d}$ back onto $X$, has an
equilibrium distribution, with probability density function (PDF) $\rho(x)$
111The density function is with respect to some underlying measure, typically
Lebesgue. Throughout this paper we will use the same symbol for the
distribution and associated PDF, with meaning taken from context.. A necessary
condition is that the distribution $\rho$ is _invariant_ under $M$, i.e. if
$x\sim\rho$ ($x$ is distributed as $\rho$) then so is $M(x)$, and further that
the orbit of almost all points $x\in X$ defined as
$\mathcal{O}^{+}(x)=\left\\{x,M(x),M^{2}(x),M^{3}(x),\ldots\right\\}$ tends in
distribution to $\rho$. Then, under mild conditions, the map is ergodic for
$\rho$, that is, expectations with respect to $\rho$ may be replaced by
averages over the orbit dorfman1999introduction ; LasotaMackey .

For example, it is well known UlamvonNeumann1947 ; may1976simple ;
dorfman1999introduction ; LasotaMackey  that the logistic map
$M_{\text{log}}(x)=4x(1-x)$, for $x\in[0,1]$, is chaotic with the equilibrium
distribution having PDF

$\rho_{\text{log}}(x)=\frac{1}{\pi\sqrt{x(1-x)}},$

implying that $M_{\text{log}}$ is ergodic for $\rho_{\text{log}}$.

Our motivating interest is in using this ergodic property to implement sample-
based inference for Bayesian analysis or machine learning. In those settings,
the target distribution $\rho$ is defined by the application. Generating a
sequence $\left\\{x_{0},x_{1},x_{2},x_{3},\ldots\right\\}$ that is ergodic for
$\rho$ is useful because then expectations of any quantity of interest can be
computed as averages over the sequence, i.e.,

$\lim_{N\to\infty}\frac{1}{N}\sum_{i=0}^{N-1}g\left(x_{i}\right)=\int_{X}g(x)\rho(x)\,\mathrm{d}x.$
(1)

Commonly in statistics, such ergodic sequences are generated using a
_stochastic_ iteration that simulates a Markov chain targeting $\rho$ mcmc_liu
; here we explore _deterministic_ iterations that generate an orbit that is
ergodic for $\rho$.

The equilibrium distribution for a given iterated function, if it exists, can
be approximated by computing the orbit of the map for some starting point then
performing kernel density estimation, or theoretically by seeking the
stationary distribution of the Frobenius–Perron (FP) operator that is the
transfer KlusKoltaiSchutte2016 , aka push-forward, operator induced by a
deterministic map dorfman1999introduction ; LasotaMackey ; we present the FP
equation in Section II.

The inverse problem that we consider here, of determining a map that gives a
prescribed equilibrium distribution, is the _inverse Frobenius–Perron problem_
(IFPP) and has been studied extensively grossmann1977invariant ;
ershov1988solution ; diakonos1996construction ; diakonos1999stochastic ;
pingel1999theory ; bollt2000controlling ; nie2013new ; nie2018matrix .
Summaries of previous approaches to the IFPP are presented in
bollt2000controlling ; rogers2008synthesis  that characterize approaches as
based on conjugate functions (see grossmann1977invariant  for details) or the
matrix method (see rogers2008synthesis  for details), and wei2015solutions
that also lists the differential equation approach. Existing work almost
solely considers the IFPP in one-variable, $d=1$, the exception being the
development of a two-dimensional solution in bollt2000controlling  that is
also presented in rogers2008synthesis .

The matrix method, first suggested by Ulam Ulam1960 , solves the IFPP for a
piecewise-constant approximation to the target density, using a transition
matrix representation of the approximated FP operator. Convergence of the
discrete approximation is related to Ulam’s conjecture, and has been proved
for the multidimensional problem; see wei2015solutions  and references
therein. While the matrix method allows construction of solutions, at least in
the limit, existing methods only offer a limited class of very non-smooth
solutions, so are not clearly useful for characterizing all solutions, as we
do here. We do not further consider the matrix methods.

The development in this paper starts with the differential equation approach
in which the IFPP for restricted forms of distributions and maps is written as
a differential equation that may be solved. We re-derive some existing
solutions to the IFPP this way in Sec. III. The main contribution of this
paper is to show that the form of these solutions may be generalized to give
the general solution of the IFPP for any probability distribution in any
dimension $d$, as presented in the factorization theorem of Sec. IV. This
novel factorization represents solutions of the IFPP in terms of the
Rosenblatt transformation rosenblatt-1952  and a _uniform map_ , that is a map
on $[0,1]^{d}$ that leaves the uniform distribution invariant. In particular,
we show that the conjugating functions in grossmann1977invariant  are exactly
the inverse Rosenblatt transformations. For a given Rosenblatt transformation,
there is a one-to-one correspondence between solution of the IFPP and choice
of a uniform map.

This reformulation of the IFPP in terms of two well-studied constructs leads
to practical analytic and numerical solutions by exploiting existing, well-
developed methods for Rosenblatt transformations, and for deterministic
iterations that target the uniform distribution. The factorization also allows
us to establish equivalence of solutions of the IFPP and other methods that
employ a deterministic map within the generation of ergodic sequences. This
standardizes and simplifies existing solution methods by showing that they are
special cases of constructing the Rosenblatt transformation (or its inverse),
and selection of a uniform map.

This paper starts with definitions of the IFPP and the Lyuapunov exponent in
Sec. II. Solutions of the IFPP in one-dimension, $d=1$, are developed in Sec.
III. These solutions for $d=1$ motivate the factorization theorem in Sec. IV
that presents a general solution to the IFPP for probability distributions
with domains in $\mathbb{R}^{d}$ for any $d$. Sec. V presents further examples
of univariate, $d=1$, solutions to the IFPP based on the factorization theorem
in Sec. IV. Two 2-dimensional numerical examples are presented in Section VI.3
to demonstrate that the theoretical constructs may be implemented in practice.
A summary and discussion of results is presented in Section VII, including a
discussion of some existing computational methods that can be viewed as
implicitly implementing the factorization solution of the IFPP presented here.

##  II Inverse Frobenius–Perron Problem and Lyapunov Exponent

In this section we define the forward and inverse Frobenius–Perron problems,
and also the Lyuapunov exponent that measures chaotic behaviour.

###  II.1 Frobenius–Perron Operator

A deterministic map $x_{n+1}=M(x_{n})$ defines a map on probability
distributions over state, called the transfer operator KlusKoltaiSchutte2016 .
Consider the case where the initial state $x_{0}\sim\rho_{0}(\cdot)$ ($x_{0}$
is distributed as $\rho_{0}$) for some distribution $\rho_{0}$ and let
$\rho_{n}$ denote the $n$-step distribution, i.e., the distribution over
$x_{n}=M^{n}(x_{0})$ at iteration $n$. The transfer operator that maps
$\rho_{n}\mapsto\rho_{n+1}$ induced by $M$ is given by the Frobenius–Perron
operator associated with $M:x\mapsto y$ dorfman1999introduction ; LasotaMackey
; Gaspard1998

$\rho_{n+1}(y)=\sum_{x\in M^{-1}(y)}\frac{\rho_{n}\left(x\right)}{|J(x)|}$ (2)

where $|J(x)|$ denotes the Jacobian determinant 222We have used the language
of differential maps, as all the maps that we display in this paper are
differentiable almost everywhere varberg1971change. More generally,
$|J(x)|^{-1}$ denotes the density of $\rho_{n}M^{-1}$ with respect to
$\rho_{n+1}$, see, e.g., (LasotaMackey, , Remark 3.2.4.). of $M$ at $x$, and
the sum is over inverse images of $y$.

The equilibrium distribution $\rho$ of $M$ satisfies

$\rho(y)=\sum_{x\in M^{-1}(y)}\frac{\rho(x)}{|J(x)|}$ (3)

and we say that $\rho$ is invariant under $M$.

###  II.2 Inverse Frobenius–Perron Problem

The inverse problem that we address is finding an iterative map $M$ that has a
given distribution $\rho$ as its equilibrium distribution. We do this by
performing the inverse Frobenius–Perron problem (IFPP) of finding $M$ that
satisfies (3) to ensure that $\rho$ is an invariant distribution of $M$.
Establishing chaotic hence ergodic behaviour is a separate calculation.

We will assume throughout that $\rho$ is absolutely continuous with respect to
the underlying measure so that a probability density function $\rho(x)$
exists, and, furthermore, that $\rho(x)>0$, $\forall x\in X$.

###  II.3 Lyapunov Exponent

The Lyapunov exponent $h$ of an iterative map gives the average exponential
rate of divergence of trajectories. We define the (maximal) Lyapunov exponent
$h$ as dorfman1999introduction ; LasotaMackey

$h=\lim_{N\to\infty}\frac{1}{N}\log\left|\frac{\mathrm{d}x_{N}}{\mathrm{d}x_{0}}\right|=\lim_{N\to\infty}\frac{1}{N}\sum^{N-1}_{n=0}\log|J(x_{n})|$
(4)

that features the starting value $x_{0}$. For ergodic maps the dependency on
$x_{0}$ is lost as $N\to\infty$ and the Lyapunov exponent may be written

$h=\int_{X}\log|J(x)|\rho(x)\,\mathrm{d}x$ (5)

where $\rho(x)$ is the invariant density. A positive Lyapunov exponent $h$
indicates that the map is chaotic.

The theoretical value for the Lyapunov exponent may be obtained using (5),
while (4) provides an empirical value obtained by iterating the map $M$. For
example, the Lyapunov exponent of the logistic map evaluated by (5) is
$h_{\text{log}}=\log 2\approx 0.693147$ while (4) evaluated over an orbit with
$10000$ iterations gave $h_{\text{log}}\approx 0.693140$.

For chaotic maps, any uncertainty in initial value means that an orbit cannot
be precisely predicted, since initial states with any separation become
arbitrarily far apart, within $X$, as iterations increase. Hence it is useful
to characterize the orbit statistically, in terms of the equilibrium
distribution over states in the orbit.

It is interesting to note that theoretical chaotic and ergodic behaviour does
not necessarily occur when iterations of a map are implemented on a finite-
precision computer. For example, when the logistic map, above, on $[0,1]$ is
iterated on a binary computer the multiplication by $4$ corresponds to a shift
left by two bits and all subsequent operations maintain lowest order bits that
are zero. Repeated iterations eventually produce the number zero, no matter
the starting value. While it is simple to correct this non-ideal behaviour, as
was done to give the numerical Lyapunov exponent, above, it is important to
note that computer implementation can have very different dynamics to the
mathematical model.

##  III Solution of the IFPP in $1$-dimension

In this Section we develop solutions of the IFPP in one-dimension, $d=1$.
Without loss of generality we consider distributions on the unit interval
$X=[0,1]$, as the domain of any univariate distribution may be transformed by
a change of variable to $X=[0,1]$, including when the domain is the whole real
line $(-\infty,\infty)$.

For distributions over a scalar random variable, the FP equation for the
invariant density (3) simplifies to

$\rho(y)=\sum_{x\in M^{-1}(y)}\frac{\rho\left(x\right)}{|M^{\prime}(x)|}.$ (6)

###  III.1 The Simplest Solution

We first note, almost trivially, that the identity map $M=I$, where $I(x)=x$,
has $\rho$ as an invariant distribution, so solves the IFPP for any $\rho$.
Somewhat less trivial is to derive this simplest solution by assuming that $M$
is monotonic increasing and $M(0)=0$, so that there is just one inverse image
in (6). Writing $|M^{\prime}(x)|=\mathrm{d}M/\mathrm{d}x$ gives the
differential equation with separated variables

$\rho(M)\mathrm{d}M=\rho(x)\mathrm{d}x$ (7)

that has solution

$F(M)=F(x)$

where $F(x)=\int_{0}^{x}\rho(x^{\prime})\,\mathrm{d}x^{\prime}$ is the
cumulative distribution function (CDF) for $\rho$. If $F$ is invertible denote
the inverse by $F^{-1}$, called the the inverse distribution function (IDF),
otherwise let $F^{-1}$ denote the generalized inverse distribution function,
$\displaystyle F^{-1}(p)=\inf\\{x\in X:F(x)\geq p\\}$. Then,
$M=F^{-1}(F(x))=x$, or $M(x)=x$, almost everywhere. Hence the identity map is
the _unique_ monotonic increasing map that has $\rho$ as its invariant
distribution. Clearly, the identity map is not ergodic for $\rho$.

We may generalize this solution by setting $M(0)=k$, for some $k\in[0,1)$, and
also only requiring $M$ to be piecewise continuous. Allowing one discontinuity
in $M$ we write the integral of the separated differential equation (7) as

$F(M)=F(x)+k\mod 1$

giving the solution to the IFPP

$M(x)=\left(F^{-1}\circ T^{c}\circ F\right)(x)$ (8)

where $c=F^{-1}(k)$. Here $T^{c}$ denotes the operator that translates by $c$
with wrap-around on $[0,1)$ 333Hence $T^{c}$ is the translation operator on
the unit circle $S^{1}$

$T^{c}(y)=y+c-\lfloor y+c\rfloor,$ (9)

where $\displaystyle\lfloor\,\ \rfloor$ denotes the floor function; see Figure
1 (left). The identity map is recovered when $k=c=0$.

Figure 1: Translation operator $T^{c}$ for $c=0.3$ (left) and triangle map
$t_{1}$ (right).

It is easily seen that the Lyapunov exponent of the map in (8) is $h=0$, so
the map is not chaotic. However, interestingly, this map can generate a
sequence of states that produce a numerical integration rule with respect to
$\rho$, since appropriate choices of $c$ and number of iterations $N$ can
produce a rectangle rule quadrature or a quasi Monte Carlo lattice rule; see
Sec. IV.3.

###  III.2 Exploiting Symmetry in $\rho(x)$

When the PDF $\rho$ has reflexive symmetry about $1/2$, i.e.
$\rho(x)=\rho(1-x)$, we can simplify the FP equation (6) by assuming that the
map $M$ has the same symmetry. Specifically, we write the triangle map (see
Figure 1 (right))

$t_{1}(x)=1-2\left|x-1/2\right|$ (10)

that has reflexive symmetry about $1/2$, and write

$M(x)=m(t_{1}(x))$ (11)

where $m(x):[0,1]\to[0,1]$ is a monotonic increasing map with $m(0)=0$ (and,
it will turn out, $m(1)=1$). Hence the FP equation simplifies to

$\rho(y)=2\frac{\rho\left(x\right)}{|M^{\prime}(x)|},\quad x\in M^{-1}(y),$
(12)

that we can write as the separated equations

$\displaystyle\rho(M)\,\mathrm{d}M=2\rho(x)\mathrm{d}x,\quad$ $\displaystyle
x<1/2,$ $\displaystyle\rho(M)\,\mathrm{d}M=-2\rho(x)\mathrm{d}x,\quad$
$\displaystyle x>1/2,$

that has the continuous solution

$F(M)=t_{1}(F(x))$

giving the continuous solution to the IFPP

$M(x)=\left(F^{-1}\circ t_{1}\circ F\right)(x).$ (13)

One can solve for $m(x)=F^{-1}(2F(x/2))$, though we do not further consider
the function $m$.

The approach we have used here simplifies the approach in
diakonos1996construction , while ‘doubly symmetric’ maps of the form (11) were
considered in gyorgyi1984fully , and again in pingel1999theory ;
diakonos1999stochastic .

###  III.3 Symmetric Triangular Distribution

To give a concrete example of the solution in (13), we consider the symmetric
triangular distribution on $[0,1]$ with PDF

$\rho_{\text{tri}}(x)=2-4\left|x-\frac{1}{2}\right|$ (14)

that has reflexive symmetry about $x=\frac{1}{2}$. The CDF is

$F_{\text{tri}}(x)=\begin{cases}2x^{2}&0\leq x\leq\frac{1}{2}\\\
\frac{1}{2}+4x-2x^{2}&\frac{1}{2}\leq x\leq 1\end{cases}$

giving the unimodal map, after substituting into (13),

$M_{\text{tri}}(x)=\left\\{\begin{array}[]{rl}\sqrt{2}x&0\leq
x\leq\frac{1}{\sqrt{8}}\\\ 1-\sqrt{\frac{1}{2}-2x^{2}}&\frac{1}{\sqrt{8}}\leq
x\leq\frac{1}{2}\\\ 1-\sqrt{\frac{1}{2}-2(1-x)^{2}}&\frac{1}{2}\leq x\leq
1-\frac{1}{\sqrt{8}}\\\ \sqrt{2}(1-x)&1-\frac{1}{\sqrt{8}}\leq x\leq
1\end{array}\right.$ (15)

shown in Figure 2 (left). The same map was derived in grossmann1977invariant .
Figure 2 (right) shows a normalized histogram of $10^{6}$ iterates of
$M_{\text{tri}}$ starting at $x=0.3$, confirming that the orbit of
$M_{\text{tri}}$ converges to the desired triangular distribution. The
numerical implementation avoids finite-precision effects, as discussed later.

Figure 2: Iterative map $M_{\text{tri}}$ in (15) (left) and a histogram of
$1\times 10^{6}$ iterates of the map $M_{\text{tri}}$ (right).

The theoretical Lyapunov exponent for $M_{\text{tri}}$ is $h_{\text{tri}}=\log
2\approx 0.693147$ while (4) evaluated over an orbit with $10^{6}$ iterations
gave $h_{\text{tri}}\approx 0.693148$.

##  IV Solutions of the IFPP for General Multi-Variate Target Distributions

The solutions to the $1$-dimensional IFPP with special structure in Eqs (8)
and (13) are actually examples of a general solution to the IFPP for multi-
variate probability distributions with no special structure. We state that
connection via a theorem that establishes a factorization of all possible
solutions to the IFPP, and that also provides a practical means of solving the
IFPP.

We first introduce the forward and inverse Rosenblatt transformations, that is
the multi-variate generalization the CDF and IDF for univariate distributions.

###  IV.1 Forward and Inverse Rosenblatt Transformations

A simple transformation of an absolutely continuous $d$-variate distribution
into the uniform distribution on the $d$-dimensional hypercube was introduced
by Rosenblatt rosenblatt-1952 , as follows. The joint PDF can be written as a
product of conditional densities,

$\rho(x_{1},\ldots,x_{d})=\rho_{1}(x_{1})\rho_{2}(x_{2}|x_{1})\cdots\rho_{d}(x_{d}|x_{1}\ldots,x_{d-1}),$

where $\rho_{k}(x_{k}|x_{1}\ldots,x_{k-1})$ is a conditional density given by

$\rho_{k}(x_{k}|x_{1}\ldots,x_{k-1})=\frac{p_{k}(x_{1},\ldots,x_{k})}{p_{k-1}(x_{1}\ldots,x_{k-1})},$
(16)

in terms of the marginal densities,

$p_{k}=\int\rho(x_{1},\ldots,x_{d})\,\mathrm{d}x_{k+1}\cdots\mathrm{d}x_{d},$
(17)

where $k=1,\ldots,d$.

Let $z=(z_{1},\ldots,z_{d})=R(x_{1},\ldots,x_{d})$ where $R$ is the Rosenblatt
transformation rosenblatt-1952  from the state-space
$X\subseteq\mathbb{R}^{d}$ of $\rho$ to the $d$-dimensional unit cube
$[0,1]^{d}$, defined in terms of the (cumulative) distribution function $F$ by

$\displaystyle z_{1}$
$\displaystyle=F_{1}(x_{1})=\int_{-\infty}^{x_{1}}\rho_{1}(x_{1}^{\prime})\,\mathrm{d}x_{1}^{\prime},$
$\displaystyle z_{2}$
$\displaystyle=F_{2}(x_{2}|x_{1})=\int_{-\infty}^{x_{2}}\rho_{2}(x_{2}^{\prime}|x_{1})\,\mathrm{d}x_{2}^{\prime},$
$\displaystyle\vdots$ $\displaystyle z_{d}$
$\displaystyle=F_{d}(x_{2}|x_{1},\ldots,x_{d-1})=\int_{-\infty}^{x_{d}}\rho_{d}(x_{d}^{\prime}|x_{1},\ldots,x_{d-1})\,\mathrm{d}x_{d}^{\prime}.$

As noted in rosenblatt-1952 , there are $d!$ transformations of this type,
corresponding to the $d!$ ways of ordering the coordinates. Further
multiplicity is introduced by considering coordinate transformations, such as
rotations.

Notice that in $1$-dimension, the Rosenblatt transformation $R(x)$ is simply
the CDF $F(x)$.

It follows that if $x\sim\rho$ then
$z=R(x)\sim\operatorname{Unif}([0,1]^{d})$, i.e., $z$ is uniformly distributed
on the $d$-dimensional unit cube rosenblatt-1952 . When $\rho(x)>0$, $\forall
x\in X$ the distribution functions are strictly monotonic increasing and the
inverse of the Rosenblatt transformation $R^{-1}$ is well defined, otherwise
let $R^{-1}$ denote the generalized inverse as in Sec. III.1. Then, if
$z\sim\operatorname{Unif}([0,1]^{d})$ it follows that $x=R^{-1}(z)\sim\rho$,
i.e., $x$ is distributed as the desired target distribution $\rho$
johnson1987multivariate ; this is the basis of the _conditional distribution
method_ for generating multi-variate random variables, that generalizes the
inverse cumulative transformation method for univariate distributions devroye-
rvgen-1986 ; johnson1987multivariate ; hormann-rvgen-2004 ;
dolgov2020approximation . These results may also be established by
substituting $R$ or $R^{-1}$ into the the FP equation (2), noting that there
is a single inverse image and that the Jacobian determinant of $R$ equals the
target PDF $\rho(x)$.

In the remainder of this paper, we refer to _any_ map $R$ satisfying
$x\sim\rho\Rightarrow R(x)\sim\operatorname{Unif}([0,1]^{d})$ as a Rosenblatt
transformation, with (generalized) inverse as defined above.

###  IV.2 Factorization Theorem

The following theorem characterizes solutions to the IFPP.

######  Theorem 1.

Given a probability distribution $\rho$ in $d$-dimensions, a map $M(x)$ is a
solution of the IFPP, that is, $M(x)$ satisfies the FP equation (3), if and
only if

$M(x)=(R^{-1}\circ U\circ R)(x),$ (18)

where $R$ is a Rosenblatt transformation and $U$ is a ‘uniform map’ on the
unit $d$-dimensional hypercube, i.e. a map that has
$\operatorname{Unif}([0,1]^{d})$ as invariant distribution 444When $U$
satisfies the stronger condition that $\operatorname{Unif}([0,1]^{d})$ is the
equilibrium distribution, $U$ is called an exact map LasotaMackey ..

Proof: We will show that $\rho$ is an invariant distribution of $M$ iff $M$
has the form (18). ($\Rightarrow$) Assume $M$ has the form (18). If
$x\sim\rho$ then $R(x)\sim\operatorname{Unif}([0,1]^{d})$ hence
$U(R(x))\sim\operatorname{Unif}([0,1]^{d})$, as
$\operatorname{Unif}([0,1]^{d})$ in invariant under $U$, hence
$M(x)=R^{-1}(U(R(x)))\sim\rho$, as desired. ($\Leftarrow$) If $\rho$ is
invariant under $M$ then $U=R\circ M\circ R^{-1}$ is a uniform map, since if
$z\sim\operatorname{Unif}([0,1]^{d})$ then $R^{-1}(z)\sim\rho$,
$M(R^{-1}(z))\sim\rho$, and
$R(M(R^{-1}(z)))\sim\operatorname{Unif}([0,1]^{d})$. Inserting this $U$ into
(18) gives the desired factorization. ∎

The first part of the proof shows that any uniform map $U$ induces a solution
to the IFPP, though the particular solution depends on the particular
Rosenblatt transformation. The second part of the proof shows that different
solutions to the IFPP effectively differ only by the choice of the uniform map
$U$, once the Rosenblatt transformation is determined, that is, a coordinate
system is chosen with an ordering of those coordinates.

Grossmann and Thomae grossmann1977invariant  called dynamical systems $M$ and
$U$ related by a formula of the form (18) as ‘related by conjugation’, or just
‘conjugate’, and the map $R^{-1}$ in (18) is a ‘conjugating function’. Thus,
in the language of grossmann1977invariant , Theorem 1 shows that the IFPP for
any distribution $\rho$ has a solution (actually, it shows that there are
infinitely many solutions), every solution map is conjugate to a uniform map,
and the conjugating function is precisely the inverse Rosenblatt
transformation.

Notice that both the translation operator $T^{c}$ in (9) and the triangle map
in (10) are uniform maps on the unit interval $[0,1]$. Thus, the solutions to
the IFPP given in Eqs (8) and (13) are examples of the general solution form
in (18). In particular, while the solution to the IFPP in (13) was derived
assuming symmetry of the target density $\rho(\cdot)$, (13) actually gives a
solution of the IFPP for _any_ density $\rho(\cdot)$. Unimodal maps of this
form were derived in ershov1988solution .

Computed examples of solutions to the IFPP given by the factorization (18) are
presented in Sec. V, in one dimension, and in Sec VI in two dimensions. High
dimensional calculations are discussed in Sec. VII.

###  IV.3 Properties of $M$ from $U$

Many properties of the map $M$ are inherited from the uniform map $U$.

When $R$ and $R^{-1}$ are continuous, $M$ is continuous iff $U$ is continuous.
In one dimension, monotonicity of the CDF and IDF imply that the number of
modes of $U$ equals the number of modes of $M$; in particular, $M$ is unimodal
iff $U$ is unimodal.

Constructing iterative maps with specific periodicity of the orbit is possible
through the use of translation operators $T^{c}$ as uniform maps, defined in
Eq. (9). Consider, first, maps in one dimension on $[0,1]$. If the shift
$c\neq 0$ and $c\notin\mathbb{Q}$, the map is aperiodic. However, in the case
that $c\neq 0$ and $c\in\mathbb{Q}$ such that

$c=\frac{N}{D}$ (19)

with $N,D\in\mathbb{N}$ and $\gcd(N,D)=1$, then the map is periodic with
periodicity $D$, and iterative maps constructed with $U=T^{c}$ exhibit the
same periodicity. These properties may be extended to multi-dimensional
settings when the translation constant $c$ is a vector of shifts in each
coordinate direction, as used in rank-one lattice rules for quasi-Monte Carlo
integration dick2013high .

The factorization in Theorem 1 also shows that performing an iteration
$x_{n+1}=M(x_{n})$ with an iterative map $M$ on the space $X$ is equivalent to
applying the corresponding uniform map $z_{n+1}=U(z_{n})$ on the space
$[0,1]^{d}$ through the transformations $R$ and $R^{-1}$, as indicated in the
following (commuting) diagram.

${x_{n}}$${z_{n}}$${x_{n+1}}$${z_{n+1}}$$\scriptstyle{\displaystyle
R}$$\scriptstyle{\displaystyle M}$$\scriptstyle{\displaystyle
U}$$\scriptstyle{\displaystyle R^{-1}}$$\scriptstyle{\displaystyle
R}$$\scriptstyle{\displaystyle R^{-1}}$

Using this commuting property, it is straightforward to prove the following
lemma.

######  Lemma 1.

For given distribution $\rho$, let $R$ be a Rosenblatt transformation for
$\rho$. Let $M=R^{-1}\circ U\circ R$ be a solution of the IFFP, as guaranteed
by Theorem 1, where $U$ is a uniform map. Then

$\mathcal{O}^{+}_{M}(x_{0})=R^{-1}\left(\mathcal{O}^{+}_{U}(R(x_{0}))\right).$
(20)

Hence, instead of iterating $M$ on the space $X$ to produce the sequence
$\\{x_{1},x_{2},x_{3},\ldots\\}$, Eq. (20) shows that one can iterate the map
$U$ on the space $[0,1]^{d}$ to produce the sequence
$\\{z_{1},z_{2},z_{3},\ldots\\}$ and then evaluate $x_{n}=R^{-1}z_{n}$,
$n=1,2,\ldots$, to produce exactly the same sequence on $X$. Since the map $M$
is mixing or ergodic iff the uniform map $U$ is mixing or ergodic,
respectively, in this sense mixing and ergodicity of $M$ is inherited from
$U$.

Using the expansion in Eq. (20), we see that $M$ is deterministic or
stochastic iff $U$ is deterministic or stochastic, respectively. Even though
we are not considering stochastic maps here, we note that, for stochastic
maps, iterations of $M$ are correlated or independent iff iterations of $U$
are correlated or independent, respectively.

Some other properties that are, and are not, preserved by the transformation
from $U$ to $M$ are discussed in grossmann1977invariant .

##  V Examples in One Dimension

###  V.1 Uniform Maps on $[0,1]$

We have already encountered three uniform maps on the interval $[0,1]$, namely
the identity map $I(x)=x$ (Figure 3 (top, left)) and the translation operator
(9) (Figure 1 (left)), that have Lyapunov exponent $h=0$, and the triangle map
(Figure 1 (right)) with Lyapunov exponent $h=\log 2$.

Some further elementary uniform maps on $[0,1]$, and associated Lyapunov
exponents, are:

  * •

$\ell$ periods of a sawtooth function on $[0,1]$ (Figure 3, top-right, for
$l=3$ periods) 555The two period sawtooth map $s_{2}$ is also called the
Bernoulli map, and its orbit $\mathcal{O}^{+}(x)$ is the dyadic
transformation.

$s_{\ell}(x)=\ell*x-\lfloor\ell*x\rfloor,$ (21)

with Lyapunov exponent $h=\log\ell$,

  * •

$\ell$ periods of a triangle function on $[0,1]$ (Figure 3 (bottom, left) for
$l=3$ periods) 666This is the ‘broken linear transformation’ in
grossmann1977invariant  of order $p=2\ell$.

$t_{\ell}(x)=1-2\left|s_{\ell}(x)-1/2\right|.$ (22)

with Lyapunov exponent $h=\log 2\ell$,

  * •

and the asymmetric triangle, for $c\in(0,1)$ (Figure 3 (bottom, right) for
$c=0.3$)

$t^{c}(x)=\begin{cases}\displaystyle\frac{x}{c}&0\leq x\leq c\\\\[10.00002pt]
\displaystyle\frac{1-x}{1-c}&c\leq x\leq 1\end{cases}$ (23)

with Lyapunov exponent $0\leq h=-c\log c-(1-c)\log(1-c)\leq\log 2$.

Figure 3: Four examples of uniform maps on $[0,1]$. Top-left: identity map.
Top-right: $\ell$ periods of a sawtooth wave for $\ell=3$. Bottom-left: $\ell$
periods of a triangle wave for $\ell=3$. Bottom-right: asymmetric triangle for
$c=0.3$.

Obviously, many more uniform maps are possible. Further examples can be formed
by partitioning the domain and range of any uniform map, and then permuting
the subintervals. Many existing ‘matrix-based’ methods for constructing
solutions to the IFPP, can be viewed as examples of such a partition-and-
permute of an elementary uniform map rogers2004synthesizing . Uniform maps of
other forms are developed in huang2005characterizing  from two-segmental
Lebesgue processes, producing uniform maps that are curiously non-linear.

######  Lemma 2.

The composition of uniform maps is also a uniform map, i.e., if $U_{1}$ and
$U_{2}$ are uniform maps then so is $U=U_{1}\circ U_{2}$.

An example is $t_{\ell}$, that can be constructed as the composition
$t_{\ell}=s_{\ell}\circ t_{1}$.

We mentioned the numerical artifacts that can occur with finite-precision
arithmetic, particularly when implementing maps on a binary computer and when
the endpoints of the interval $X$ and constants in the maps have exact binary
representations. We avoided these artifacts by composing the stated uniform
map with the translation $T^{c}$ for $c=1/3\times 10^{-9}$ that does not have
finite binary representation 777Computation was performed in MatLab that
implements IEEE Standard 754 for double-precision binary floating-point
format.. This small shift is indiscernible in the graphs of the maps.

###  V.2 Ramp Distribution

To give a concrete example of the solution in (13) for a distribution without
reflexive symmetry, we consider the ramp distribution with PDF

$\rho_{\text{ramp}}(x)=2x$ (24)

that has CDF

$F_{\text{ramp}}(x)=x^{2}.$

We produce a unimodal, continuous map by choosing the uniform map $t_{1}$, as
in (13), to give

$M_{\text{ramp}}=\begin{cases}\sqrt{2}x&0\leq x\leq\frac{1}{\sqrt{2}}\\\
\sqrt{2}\sqrt{1-x^{2}}&\frac{1}{\sqrt{2}}\leq x\leq 1\end{cases}$ (25)

shown in Figure 4 (left). Figure 2 (right) shows a normalized histogram of
$10^{6}$ iterates of $M_{\text{ramp}}$ starting at $x=0.3$, confirming that
the orbit of $M_{\text{ramp}}$ converges to the desired ramp distribution, as
guaranteed by Theorem 1. The numerical implementation avoids finite-precision
effects, as discussed earlier.

Figure 4: Iterative map $M_{\text{ramp}}$ in (25) (left) and a normalized
histogram of $1\times 10^{6}$ iterates that approximates the equilibrium PDF
(right).

The estimated Lyapunov exponent for this map is $h\approx 1.040035$, which is
greater than the Lyapunov exponent for the inducing triangular map $t_{1}$
which is $\log 2\approx 0.693147$.

Using a different uniform map gives a different solution to the IFPP. For
example, choosing $s_{3}$ gives the map $M=F^{-1}\circ s_{3}\circ F$ shown in
Figure 5 (left). A normalized histogram over an orbit of $10^{6}$ iterations
is shown in Figure 5 (right), confirming that this map is also ergodic for
$\rho_{\text{ramp}}$.

Figure 5: Iterative map $M_{\text{ramp}}$ in (25) (left) and a normalized
histogram of $1\times 10^{6}$ iterates that approximates the equilibrium PDF
(right).

The estimated Lyapunov exponent for this map is $h\approx 1.098612$, which is
the same numerical value as the Lyapunov exponent for the sawtooth map $s_{3}$
which is $\log 3\approx 1.098612$. Comparing with the result for the map
induced by $t_{1}$, this shows that the Lyapunov exponent of a map is not
generally equal to the Lyapunov exponent of the inducing uniform map.

###  V.3 The Logistic Map and Alternatives

The logistic map, mentioned in the introduction, is

$M_{\text{log}}(x)=4x(1-x).$ (26)

The equilibrium distribution of this map

$\rho_{\text{log}}(x)=\frac{1}{\pi\sqrt{x(1-x)}},$ (27)

can be easily verified by substituting into the FP equation (6). The CDF of
$\rho_{\text{log}}(x)$ is

$F_{\text{log}}(x)=\int_{0}^{x}\rho_{\text{log}}(x^{\prime})dx^{\prime}=\frac{2}{\pi}\arcsin(\sqrt{x}),$
(28)

and the IDF is

$F_{\text{log}}^{-1}(x)=\sin^{2}(\frac{\pi x}{2})=\frac{1}{2}(1-\cos(\pi x)).$
(29)

The logistic map (26) is induced by the factorization (18) by choosing the
triangle map $t_{1}$ as uniform map, i.e., substituting the CDF (28) and IDF
(29) into (13); see Figure 6 (left). Equivalently, one may note that the
logistic map (26) is transformed into the triangle map $t_{1}$ by the change
of variables $z=F^{-1}(x)$; in the language of grossmann1977invariant ,
$M_{\text{log}}$ and $t_{1}$ are _conjugate_ dynamical laws.

Other iterative maps which preserves the same equilibrium distribution (27)
can be constructed by choosing another uniform map, such as $\ell$ periods of
a triangle function (22). This gives the iterative maps

$M_{\ell}=F_{\text{log}}^{-1}\circ t_{\ell}\circ
F_{\text{log}}=\sin^{2}(2\ell\arcsin(\sqrt{x})),\ell\geq 1,$ (30)

that coincide with the $n^{\text{th}}$ power of the logistic map (26) for
$\ell=2^{n-1}$. Figure 6 (right) shows the map which preserves the same
equilibrium distribution as the logistic map but induced by the uniform map
$t_{3}$. Since $3$ is not of the form $2^{n-1}$, this map is not simply a
power of the logistic map.

Figure 6: Logistic map, that equals $M_{\text{log}}=F_{\text{log}}^{-1}\circ
t_{1}\circ F_{\text{log}}$ (left), and the map $M=F_{\text{log}}^{-1}\circ
t_{3}\circ F_{\text{log}}$ that has the same equilibrium distribution (right).

The theoretical value of the Lyapunov exponent of the map in (30) is $\log
2\ell$, using (5). Table 1 gives the theoretical values of the Lyapunov
exponent and experimentally calculated values using $10000$ iterations, as in
(4), for some values of $\ell$.

$\ell$ | $h_{\text{e}}$ | $h_{\text{t}}$
---|---|---
$1$ | $0.692819$ | $0.693147$
$2$ | $1.386284$ | $1.386294$
$4$ | $2.079430$ | $2.079442$
$2^{24}$ | $17.328594$ | $17.328680$
$2^{39}$ | $27.725713$ | $27.725887$
Table 1: Experimental $h_{\text{e}}$ and theoretical $h_{\text{t}}$ Lyapunov
exponents for the maps in (6) for a range of $\ell$, given to $6$ decimal
places.

##  VI Two Examples in Two Dimensions

###  VI.1 Uniform Maps on $[0,1]^{2}$

Two well-known examples of uniform maps in the two-dimensional unit square,
$X=[0,1]^{2}$, are the baker’s map

$U_{\text{b}}(x_{1},x_{2})=\left(2x_{1}\mod
1,\frac{x_{2}+u\left(x_{1}-\frac{1}{2}\right)}{2}\right),$ (31)

where $u$ is the unit step function, and the Arnold cat map

$U_{\text{A}}(x_{1},x_{2})=\left((2x_{1}+x_{2})\mod 1,(x_{1}+x_{2})\mod
1\right).$ (32)

Other uniform maps in $d>1$ dimensions may be formed by $1$-dimensional
uniform maps acting on each coordinate, giving the coordinate-wise uniform map

$U(x)=\left(U_{1}(x_{1}),U_{2}(x_{2}),\ldots,U_{n}(x_{d})\right)$ (33)

where $U_{i}(x)$, $i=1,2,\ldots,d$, are uniform maps in one dimension. We will
use the baker’s map (31) and a coordinate-wise uniform map in the
$2$-dimensional examples that follow.

###  VI.2 Checker-Board Distribution

This example demonstrates construction of a map in two-dimensions that targets
the checker board distribution, shown in Figure 8 (bottom-left) and Figure 9
(bottom-left), using the factorization (18).

The first step in constructing a solution to the IFPP for this distribution is
to construct the forward and inverse Rosenblatt transformations, which
requires the marginal density functions (17) that may be evaluated
analytically in this case. A plot of the two components of the functions $R$
and $R^{-1}$ is shown in Figure 7.

Figure 7: Plots of the Rosenblatt transformation $R$ and its inverse $R^{-1}$
for the checker-board distribution. Top row shows the components of $R$: first
component (left) and second component (right). Bottom row shows the components
of $R^{-1}$: first component (left) and second component (right).

We construct two solutions to the IFPP, each induced by choosing a particular
uniform map: the first is the baker’s map (31); the second is a component-wise
uniform map (33) with an asymmetric triangle map (23) acting on each
component,

$U(x_{1},x_{2})=\left(U_{1}(x_{1}),U_{2}(x_{2})\right)$ (34)

where $U_{1}=t^{c}$ for $c=0.3$ and $U_{2}=t^{c}$ with $c=0.9$.

Figure 8 (top row) shows the two components of the map induced by the baker’s
map, the checker-board distribution (bottom-left), and a histogram of $10^{6}$
iterates (bottom-right) showing that the map does indeed converge in
distribution to the desired distribution.

Figure 8: Iterated function constructed with $U$ being the baker’s
transformation, a histogram of iterates, and the checker board distribution.
Shown is: the $x_{1}$ part of the constructed map (top-left), the $x_{2}$ part
of the constructed map (top-right), the checker board distribution (bottom-
left), and a histogram of iterates of the constructed map (bottom-right).

Figure 9 (top row) shows the two components of the map constructed using the
two component-wise asymmetric triangular maps, the checker-board distribution
(bottom-left), and a histogram of $10^{6}$ iterates (bottom-right) showing
that the map also converges in distribution to the desired distribution, and
hence is also a solution to the IFPP.

Figure 9: Same as Figure 8 except $M$ is constructed with $U$ being component-
wise asymmetric triangular maps.

###  VI.3 A Numerical Construction

Numerical implementation of the factorized solution (18) is not difficult in
few dimensions. In this section we present an example of numerical
implementation using a normalized greyscale image of a coin 888A pre-2006 New
Zealand $50$ cent coin., piecewise constant over pixels, as the target
distribution; see Figure 10 (left). The marginal distributions (17) are
evaluated as linear interpolation of cumulative sums over pixel values, and
hence the CDF and then forward and inverse Rosenblatt transformations follow
as in section IV.1. The uniform map was produced as component-wise univariate
translation maps, specifically

$U(x_{1},x_{2})=\left(U_{1}(x_{1}),U_{2}(x_{2})\right)$

where $U_{1}=T^{c}$ for $c=0.6$ and $U_{2}=T^{c}$ with $c=0.2$. The resulting
map is given by Eq. (18).

Figure 10: A normalized greyscale image of a coin used as the target
distribution, and a normalized histogram of iterations of the map targeting
this distribution. Original image (left), and normalized histogram of $10^{6}$
iterations (right).

Figure 10 (right) shows a normalized histogram, binned to pixels, of $10^{6}$
iterations of this map. As can be seen, the estimated PDF from the orbit of
this map does reproduce the image of the coin. However, there are also obvious
artifacts near the edge of the image showing that mixing could be better. We
conjecture that a chaotic uniform map would produce better mixing and fewer
numerical artifacts.

##  VII Summary and Discussion

We have shown how the solution of the IFPP, of finding an iterative map with a
given invariant distribution, can be constructed from uniform maps through the
factorization established in Theorem 1,

$M=R^{-1}\circ U\circ R$

where $R$ denotes the Rosenblatt transformation that has Jacobian determinant
equal to density function of the invariant distribution. In one-dimension, $R$
is exactly the CDF of the given distribution, so the factorization generalizes
existing one-dimensional solutions to the setting of arbitrary multi-variate
distributions. The factorization also shows the relationship between arbitrary
iterative maps and uniform maps, i.e., given a Rosenblatt transformation the
solution of the IFPP is equivalent to the choice of a uniform map that has
$\operatorname{Unif}([0,1]^{d})$ as invariant distribution.

We find the factorization (18) appealing as it shows that solution of the IFPP
for arbitrary distributions, and in multi dimensions, is reduced to two
standard and well-studied problems, i.e., constructing the Rosenblatt
transformation (or CDF in one dimension) and designing a uniform map. It is
therefore surprising, to us, that the factorization (18), and more generally
the Rosenblatt transformation, appears to not be widely used in the study of
chaotic iterated functions and the IFPP. Grossmann and Thomae
grossmann1977invariant , in one of the earliest studies of the IFPP,
essentially derived the factorization (18) by introducing conjugate maps and
establishing the relation (in their notation) that
$\rho^{*}(x)=\mathrm{d}h^{-1}(x)/\mathrm{d}x$ where $\rho^{*}$ is the
invariant distribution and $h$ is the conjugating function; see
(grossmann1977invariant, , Figure 3). It is a small step to identify that $h$
is the IDF, generalized in multi-dimensions by the inverse Rosenblatt
transformation. However, the connection was not made in grossmann1977invariant
, despite the Rosenblatt transformation having been already known, in
statistics, for some decades rosenblatt-1952 .

We constructed solutions to the IFPP for distributions with special reflexive
symmetry structure, and then with no special structure, by constructing the
Rosenblatt transformation, and its inverse, for some examples in one and two
dimensions. For simple distributions with analytic form the Rosenblatt
transformation may be constructed analytically, while numerically-defined
distributions required calculation of the marginal distributions (17) using
numerical techniques.

Although this factorization and construction is applicable to high-dimensional
problems, the main difficulty is obtaining all necessary marginal densities,
which requires the high-dimensional integral over $x_{k+1}\ldots x_{d}$ in
(17). In general, this calculation can be extremely costly. Even a simple
discretization of the PDF $\rho$, or of the argument of the marginal densities
(17), leads to exponential cost with dimension.

To overcome this cost, Dolgov _et al._ dolgov2020approximation  precomputed an
approximation of $\rho(x_{1},\ldots,x_{d})$ in a compressed tensor train
representation that allows fast computation of integrals in (17), and
subsequent simulation of the inverse Rosenblatt transformation $R^{-1}$ from
the conditionals in (17), and showed that computational cost scales _linearly_
with dimension $d$. Practical examples presented in dolgov2020approximation ,
in dimension $d\leq 32$, demonstrate that operation by the forward and inverse
Rosenblatt transformations is computationally feasible for multivariate
problems with no special structure.

Finding a solution of the IFPP with desired properties is reduced to a
standard problem of designing a uniform map on $[0,1]^{d}$, for which there
are many existing efficient options. For example, standard computational
uniform random number generators, that produce pseudo-random sequences of
numbers, are one such existing uniform map, as are the quasi-Monte Carlo rules
that we mentioned earlier dick2013high . These induce pseudo-random and quasi-
Monte Carlo sequences, respectively, on the space $X$ via the inverse
Rosenblatt transformation $R^{-1}$ gentle2003random . Both these schemes were
demonstrated in practical high-dimensional settings in dolgov2020approximation
.

The RHS of Eq. (20) in $d=1$ dimension is exactly the standard computational
route for implementing inverse cumulative transformation sampling from $\rho$,
since computational uniform pseudo-random number generators perform a
deterministic iteration on $[0,1]$ to implement a uniform map gentle2003random
. For $d>1$, the RHS of Eq. (20) is the conditional distribution method that
generalizes the inverse cumulative transformation method, as mentioned above
devroye-rvgen-1986 ; johnson1987multivariate ; hormann-rvgen-2004 ;
dolgov2020approximation . Hence Lemma 1 shows that standard computational
implementation of both the inverse cumulative transformation method in $d=1$
and the conditional distribution method in $d>1$ is equivalent to implementing
a solution to the IFFP. In this sense, computational inverse cumulative
transformation sampling from $\rho$ can be viewed as the prototype for all
iterative maps that target the distribution $\rho$, with each ergodic sequence
corresponding to a particular choice of uniform map.

We mentioned that the Rosenblatt transformation associated with a given
distribution $\rho$ is not unique. Actually, any two Rosenblatt
transformations for $\rho$ are related by a uniform map, as shown in the
following Lemma.

######  Lemma 3.

If $R_{1}$ is a Rosenblatt transformation for $\rho$ then $R_{2}$ is a
Rosenblatt transformation for $\rho$ if and only if

$R_{2}=U\circ R_{1}$

for some uniform map $U$.

Proof: ($\Rightarrow$) Since $R_{1}$ and $R_{2}$ are Rosenblatt
transformations for $\rho$, then $U=R_{2}\circ R_{1}^{-1}$ is a uniform map
and $R_{2}=U\circ R_{1}$. ($\Leftarrow$) If $R_{2}=U\circ R_{1}$ then if
$x\sim\rho$, $R_{1}(x)\sim\operatorname{Unif}([0,1]^{d})$ and $U\circ
R_{1}(x)\sim\operatorname{Unif}([0,1]^{d})$ so $R_{2}$ is a Rosenblatt
transformation. ∎

Hence, any Rosenblatt transformation $R$ may be written as $R=U\circ R_{0}$
for some uniform map $U$ and a fixed Rosenblatt transformation $R_{0}$.

The Rosenblatt transformations that map any distribution to the uniform
distribution on the hypercube may also be used to understand mappings between
spaces that are designed to transform one distribution to another, such as the
‘transport maps’ developed in parno2018transport . Consider distributions
$\rho_{A}$ and $\rho_{B}$, with Rosenblatt transformations $R_{A}$ and
$R_{B}$, respectively, that may be related as in the following diagram.

${\rho_{A}}$${\operatorname{Unif}([0,1]^{d})}$${\rho_{B}}$$\scriptstyle{\displaystyle
R_{A}}$$\scriptstyle{\displaystyle R_{A}^{-1}}$$\scriptstyle{\displaystyle
R_{B}^{-1}}$$\scriptstyle{\displaystyle U}$$\scriptstyle{\displaystyle R_{B}}$

The diagram suggests a proof of the following lemma, that generalizes the
factorization Theorem 1.

######  Lemma 4.

A map $M$ satisfies $x\sim\rho_{A}\Rightarrow M(x)\sim\rho_{B}$ iff it can be
written as $M=R_{B}^{-1}\circ U\circ R_{A}$, where $R_{A}$ and $R_{B}$ are
Rosenblatt transformations for $\rho_{A}$ and $\rho_{B}$, respectively, and
$U$ is a uniform map.

Hence, for given Rosenblatt transformations, the choice of a map, that maps
samples from $\rho_{A}$ to samples from $\rho_{B}$, is equivalent to the
choice of a uniform map. Alternatively, if a fixed uniform map is selected,
such as the identity map, the choice of map $M$ is completely equivalent to
the choice of Rosenblatt transformations. This factorization also shows that
the equivalence class of conjugate maps, noted in grossmann1977invariant , for
each dimension $d$, is generated by the uniform maps, and each member of the
equivalence class contains maps that target each distribution, when the
associated Rosenblatt transformation satisfies the mild conditions to be a
conjugating function as defined in grossmann1977invariant .

## References

  * [1] The density function is with respect to some underlying measure, typically Lebesgue. Throughout this paper we will use the same symbol for the distribution and associated PDF, with meaning taken from context.
  * [2] J.R. Dorfman.  Cambridge Lecture Notes in Physics: An introduction to chaos in nonequilibrium statistical mechanics, volume 14.  Cambridge University Press, 1999.
  * [3] Andrzej Lasota and Michael C. Mackey.  Chaos, fractals, and noise.  Springer-Verlag, New York, second edition, 1994.
  * [4] S. Ulam and J. von Neumann.  On combination of stochastic and deterministic processes.  Bulletin of the American Mathematical Society, 53(11):1120, 1947\.
  * [5] R.M. May.  Simple mathematical models with very complicated dynamics.  Nature, 261(5560):459–467, 1976.
  * [6] J. S. Liu.  Monte Carlo Strategies in Scientific Computing.  Springer-Verlag, 2001.
  * [7] Stefan Klus, Péter Koltai, and Christof Schütte.  On the numerical approximation of the Perron–Frobenius and Koopman operator.  Journal of Computational Dynamics, 3(1):51–79, 2016.
  * [8] Siegfried Grossmann and Stefan Thomae.  Invariant distributions and stationary correlation functions of one-dimensional discrete processes.  Zeitschrift für Naturforschung a, 32(12):1353–1363, 1977.
  * [9] S. V. Ershov and Georgii Gennad’evich Malinetskii.  The solution of the inverse problem for the Perron–Frobenius equation.  USSR Computational Mathematics and Mathematical Physics, 28(5):136–141, 1988.
  * [10] F.K. Diakonos and P. Schmelcher.  On the construction of one-dimensional iterative maps from the invariant density: The dynamical route to the beta distribution.  Physics Letters A, 211(4):199–203, 1996.
  * [11] FK Diakonos, D. Pingel, and P. Schmelcher.  A stochastic approach to the construction of one-dimensional chaotic maps with prescribed statistical properties.  Physics Letters A, 264(2-3):162–170, 1999.
  * [12] D. Pingel, P. Schmelcher, and F.K. Diakonos.  Theory and examples of the inverse Frobenius–Perron problem for complete chaotic maps.  Chaos, 9(2):357–366, 1999.
  * [13] Erik M Bollt.  Controlling chaos and the inverse Frobenius–Perron problem: global stabilization of arbitrary invariant measures.  International Journal of Bifurcation and Chaos, 10(05):1033–1050, 2000.
  * [14] Xiaokai Nie and Daniel Coca.  A new approach to solving the inverse Frobenius–Perron problem.  In 2013 European Control Conference (ECC), pages 2916–2920. IEEE, 2013.
  * [15] Xiaokai Nie and Daniel Coca.  A matrix-based approach to solving the inverse Frobenius–Perron problem using sequences of density functions of stochastically perturbed dynamical systems.  Communications in Nonlinear Science and Numerical Simulation, 54:248–266, 2018.
  * [16] A. Rogers, R. Shorten, D.M. Heffernan, and D. Naughton.  Synthesis of piecewise-linear chaotic maps: Invariant densities, autocorrelations, and switching.  International Journal of Bifurcation and Chaos, 18(8):2169–2189, 2008.
  * [17] Nijun Wei.  Solutions of the inverse Frobenius–Perron problem.  Master’s thesis, Concordia University, 2015.
  * [18] Stanisław Marcin Ulam.  A Collection of Mathematical Problems.  Interscience Publishers, 1960.
  * [19] M. Rosenblatt.  Remarks on a multivariate transformation.  Ann. Math. Stat., 23(3):470–472, 1952.
  * [20] Pierre Gaspard.  Chaos, scattering and statistical mechanics.  Cambridge University Press, 1998.
  * [21] We have used the language of differential maps, as all the maps that we display in this paper are differentiable almost everywhere [varberg1971change]. More generally, $|J(x)|^{-1}$ denotes the density of $\rho_{n}M^{-1}$ with respect to $\rho_{n+1}$, see, e.g., [3, Remark 3.2.4.].
  * [22] Hence $T^{c}$ is the translation operator on the unit circle $S^{1}$.
  * [23] G Györgyi and P Szépfalusy.  Fully developed chaotic 1-d maps.  Zeitschrift für Physik B Condensed Matter, 55(2):179–186, 1984\.
  * [24] M.E. Johnson.  Multivariate Statistical Simulation.  John Wiley & Sons, 1987.
  * [25] L. Devroye.  Non-Uniform Random Variate Generation.  Springer-Verlag, 1986.
  * [26] W. Hörmann, J. Leydold, and G. Derflinger.  Automatic Nonuniform Random Variate Generation.  Springer-Verlag, 2004.
  * [27] Sergey Dolgov, Karim Anaya-Izquierdo, Colin Fox, and Robert Scheichl.  Approximation and sampling of multivariate probability distributions in the tensor train decomposition.  Statistics and Computing, 30(3):603–625, 2020.
  * [28] When $U$ satisfies the stronger condition that $\operatorname{Unif}([0,1]^{d})$ is the equilibrium distribution, $U$ is called an exact map [3].
  * [29] Josef Dick, Frances Y Kuo, and Ian H Sloan.  High-dimensional integration: the quasi-Monte Carlo way.  Acta Numerica, 22:133, 2013.
  * [30] The two period sawtooth map $s_{2}$ is also called the Bernoulli map, and its orbit $\mathcal{O}^{+}(x)$ is the dyadic transformation.
  * [31] This is the ‘broken linear transformation’ in [8] of order $p=2\ell$.
  * [32] Alan Rogers, Robert Shorten, and Daniel M Heffernan.  Synthesizing chaotic maps with prescribed invariant densities.  Physics Letters A, 330(6):435–441, 2004.
  * [33] W. Huang.  Characterizing chaotic processes that generate uniform invariant density.  Chaos, Solitons & Fractals, 25(2):449–460, 2005.
  * [34] Computation was performed in MatLab that implements IEEE Standard 754 for double-precision binary floating-point format.
  * [35] A pre-2006 New Zealand $50$ cent coin.
  * [36] James E Gentle.  Random number generation and Monte Carlo methods, volume 381\.  Springer, 2003.
  * [37] Matthew D Parno and Youssef M Marzouk.  Transport map accelerated Markov chain Monte Carlo.  SIAM/ASA Journal on Uncertainty Quantification, 6(2):645–682, 2018\.