Approximating integrals with respect to stationary probability measures of iterated function systems

TL;DR

This paper introduces an efficient algorithm for approximating integrals with respect to stationary measures of iterated function systems, enabling better estimation of various mathematical properties like Hausdorff moments and Lyapunov exponents.

Contribution

The paper presents a novel algorithm for fast approximation of integrals related to stationary measures of iterated function systems under specific conditions.

Findings

Effective approximation of Hausdorff moments

Accurate estimation of Wasserstein distances

Reliable computation of Lyapunov exponents

Abstract

We study fast approximation of integrals with respect to stationary probability measures associated to iterated functions systems on the unit interval. We provide an algorithm for approximating the integrals under certain conditions on the iterated function system and on the function that is being integrated. We apply this technique to estimate Hausdorff moments, Wasserstein distances and Lyapunov exponents of stationary probability measures.

Tables

Table 1. Table 1 . Approximate values for the moments.

Order	Approx. moment value ($k=14$)	Actual moment value	Error
$n$		${\gamma }_n$
0	1	1	0
1	$0.\overline{6}$	2/3	0
2	$0.555555555555555555555555555555555555555555555554956$	$5/9$
3	$0.495726495726495726495726495726495726495726495726296$	$58/117$
4	$0.455270655270655270655270655270655270655270655272966$	$799/1755$
5	$0.424681939833454984970136485288000439515591030742856$	$54110/127413$
6	$0.400239922384444528966673488818010962533107055240313$	$662945/1656369$
7	$0.380103465368104155158312752269065347061865131565893$	$2064430846/5431233951$
8	$0.363174113035119080714419595189042981452860757191787$	$1213077397297/3340208879865$
9	$0.348713510641572814140979346643070945142831570520821$	$764170684622650/2191399705783431$
10	$0.336192206948823127238914429058461576988475565523927$	$16313445679660723325/48524163685162512633$

Table 2. Table 2 . Approximate values for the moments.

Order	Approximate moment value ($k=13$)
$n$
0	$1$
1	$0.330469717526485534080138479518406828981534429410127592033533774023242108\dots$
2	$0.1192803776960544798961200249581823359145663180489550186633549589883397537\dots$
3	$0.0461208401857915310881276274089274103313021776737494744720364563940134891\dots$
4	$0.0186956679319404288549585313723378148471406387052756144239642197703097910\dots$
5	$0.0078078479770635609245553780159130351213591283122475360490644426811585893\dots$
6	$0.0033201105732319037686859365646767147485917239274989027585387616786231639\dots$
7	$0.00142718211837850365241828109336276663400252833572475441921650238220944521\dots$
8	$0.000617598307785531412407227175919292282332872096283836068810742529448255746\dots$
9	$0.000268421862695922651075727295017088906810341162763633320631290756501684091\dots$
10	$0.000117017198360695628842316148053569471108768123919492595293008649132277667\dots$

Table 3. Table 3 . Approximate values for the Wasserstein distances.

Iteration	Approx. Wasserstein dsitance	Error
$k$	$w_k$
8	$0.3999999999583169972032649579562794534915$
9	$0.400000000000063388018611127635519504188$
10	$0.3999999999999999530973195805847764869926$
11	$0.4000000000000000000170204375320820851671$
12	$0.3999999999999999999999969551351875258863$
13	$0.4000000000000000000000000002694598038564$
14	$0.3999999999999999999999999999999881752292$
15	$0.4000000000000000000000000000000000001067$
16	$0.3999999999999999999999999999999999999999972$

Table 4. Table 4 . Approximate values for the Wasserstein distances.

Iter.	Approx. Wasserstein distance
$k$	$w_k$
8	$0.2210457542228009986686646214111284400538314616500475049250670839710830921157028$
9	$0.2210457542228009986686646648222322435005592621138097135184119027209843678362625$
10	$0.2210457542228009986686646648222083279218347168575015785464317336794134748041328$
11	$0.2210457542228009986686646648222083279244322328476185036762204108148528502575738$
12	$0.2210457542228009986686646648222083279244322327918382319334332352348209033215652$
13	$0.2210457542228009986686646648222083279244322327918382321725464968058605768813675$
14	$0.2210457542228009986686646648222083279244322327918382321725464966008628827356344$
15	$0.2210457542228009986686646648222083279244322327918382321725464966008628827708159$

Table 5. Table 5 . Approximate values for Lyapunov exponents.

Iter.	Approx. Lyapunov exponent
$k$	$w_k$
8	$1.73672081473731987719335668233174352405940161192454820335830004878102455942170283350510993376437455813744965215$
9	$1.73672081473731987719335669096051704652553164845149915406190548269505610944642757530134017955544773131075328248$
10	$1.73672081473731987719335669096051377335988973317409087754512488812667028438133069194736931707026244323079276022$
11	$1.73672081473731987719335669096051377336020590601023577825778088227884147013275372835954447415750580402943072648$
12	$1.73672081473731987719335669096051377336020590600637607989863098312523352089178183012079770679720814305305511767$
13	$1.73672081473731987719335669096051377336020590600637607991887362479762663771464637789116844202618748250354445394$
14	$1.73672081473731987719335669096051377336020590600637607991887362479193249845331374280452940254564000080656340808$
15	$1.73672081473731987719335669096051377336020590600637607991887362479193249845555716884110153252308930493603697019$
16	$1.73672081473731987719335669096051377336020590600637607991887362479193249845555716884109104549736936505689675419$
17	$1.73672081473731987719335669096051377336020590600637607991887362479193249845555716884109104549736965148549578668$
18	$1.73672081473731987719335669096051377336020590600637607991887362479193249845555716884109104549736965148549578625$

Equations

$$\Lambda =\bigcup _{i\in \mathcal{I}}{\varphi }_i(\Lambda ).$$

$$\int \phi \text{d}\mu =\sum _{i\in \mathcal{I}}{p}_i\int \phi \circ {\varphi }_i\text{d}\mu$$

$${\int }_0^{1}g\text{d}\mu$$

$${\mathcal{C}}_{\Phi }:=\bigcup _D\left\{g:[0,1]\to \mathbb{R}:g\text{mbox}hasanextensiontosomeboundedholomorphicfunctiononD\right\},$$

$${\sigma }^k\mathbf{i}:=$$

$${\tilde{\sigma }}^k\mathbf{i}:=$$

$$p_{\mathbf{i}}:=$$

$$\Pi (i_1{i}_2\dots ):=\underset{n\to \infty }{\lim }{\varphi }_{i_1}\circ \cdots \circ {\varphi }_{i_n}\left([0,1]\right).$$

$$t_m:=\sum _{\mathbf{i}\in {\mathcal{I}}^m}{p}_{\mathbf{i}}\frac{1}{1-{\varphi }_{\mathbf{i}}^{\prime }(z_{\mathbf{i}})}$$

$${\tau }_m:=\sum _{\mathbf{i}\in {\mathcal{I}}^m}{p}_{\mathbf{i}}\frac{g(z_{\mathbf{i}})+g(z_{\sigma \mathbf{i}})+\dots +g(z_{{\sigma }^m\mathbf{i}})}{1-{\varphi }_{\mathbf{i}}^{\prime }(z_{\mathbf{i}})}.$$

$${\alpha }_n:=\sum _{l=1}^n\frac{(-1{)}^l}{l!}\sum _{n_1+\dots +n_l=n}\sum _{j=1}^l\frac{{\tau }_{n_j}}{n_j}\prod _{1⩽m⩽l,m\ne j}\frac{t_{n_m}}{n_m}$$

$$a_n:=\sum _{l=1}^n\frac{(-1{)}^l}{l!}\sum _{n_1+\dots +n_l=n}\prod _{i=1}^l\frac{t_{n_i}}{n_i}.$$

$${\mu }_k(g):=\frac{{\sum }_{n=0}^k{\alpha }_n}{{\sum }_{n=0}^{k}n{a}_n}$$

$$\left|{\int }_0^{1}g\text{d}\mu -{\mu }_k(g)\right|<C\exp (-\lambda k^2)$$

$${\int }_0^{1}g(x)\text{d}\mu (x)$$

$$s_n(L)=\inf \{\Vert L-K\Vert :\text{rank}(K)⩽n-1\}.$$

$$\text{tr}(L)=\sum _{n=1}^{\infty }⟨L{e}_n,e_n{⟩}_H$$

$$\det (\text{Id}-zL)=\prod _{n=0}^{\infty }(1-z{\lambda }_n(L))$$

$$\det (\text{Id}-zL)=\sum _{n=0}^{\infty }{a}_n{z}^n.$$

$$a_n=\sum _{m=1}^n\frac{(-1{)}^m}{m!}\sum _{\begin{array}{c}{n}_1,\dots ,n_m\in {\mathbb{N}}^m\\ n_1+\dots +n_m=m\end{array}}\prod _{i=1}^m\frac{\text{tr}{L}^{n_i}}{n_i},$$

$$a_n=(-1{)}^n\sum _{i_1<\dots <i_n}{\lambda }_{i_1}(L)\dots {\lambda }_{i_n}(L)$$

$$|a_n|⩽\sum _{i_1<\dots <i_n}|{\lambda }_{i_1}(L)\dots {\lambda }_{i_n}(L)|.$$

$$A^2(D):=\left\{f:D\to \mathbb{C}\text{holomorphic}:{\int }_D|f(z){|}^2\text{d}V(z)<\infty \right\}$$

$$⟨f,g{⟩}_{A^2}={\int }_{D}f(z)\overline{g(z)}\text{d}V(z).$$

$$Lf(z)=\sum _{i=1}^k{w}_i(z)f({\psi }_i(z))$$

$${\lambda }_n(L)⩽C\left(\sum _{i=1}^k\Vert w_i{\Vert }_{\infty }\right)\exp (-\lambda n)$$

$$\text{tr}(L^n)=\sum _{\mathbf{i}\in {\mathcal{I}}^n}\frac{w_{\mathbf{i}}(z_{\mathbf{i}})}{1-{\psi }_{\mathbf{i}}^{\prime }(z_{\mathbf{i}})}$$

$$P_t=\frac{1}{2\pi i}{\int }_{\gamma }(s\text{Id}-L_t{)}^{-1}\text{d}s$$

$${\mathcal{L}}_{s}f(z)=\sum _{i=1}^N{p}_i\exp (sg({\varphi }_i(z)))f({\varphi }_i(z)).$$

$${\mathcal{L}}_s^{n}f(z)=\sum _{\mathbf{i}\in {\mathcal{I}}^n}{p}_{\mathbf{i}}\exp (s(g({\varphi }_{\mathbf{i}}z)+g({\varphi }_{\tilde{\sigma }\mathbf{i}}z)+\dots +g({\varphi }_{{\tilde{\sigma }}^{n-1}\mathbf{i}}z)))f({\varphi }_{\mathbf{i}}z).$$

$${\lambda }_n({\mathcal{L}}_s)⩽C\exp (-\lambda n^2)$$

Full text

1. Introduction

Let $\mathcal{I}:=\{1,\ldots,N\}$ for some $N\in\{2,3,\ldots\}$ and let $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$ be an iterated function system consisting of Lipschitz contractions $\phi_{i}:[0,1]\to[0,1]$. It is well-known that there exists a unique, non-empty and compact set $\Lambda\subset[0,1]$ such that

[TABLE]

This set is called the attractor of $\Phi.$ Given a probability vector $\mathbf{p}=(p_{1},p_{2},\ldots,p_{N})$ (meaning that each $0<p_{i}<1$ and $\sum_{i\in\mathcal{I}}p_{i}=1$), there exists a unique probability measure $\mu=\mu^{(\Phi,\mathbf{p})}$ such that

[TABLE]

for every continuous function $\varphi:[0,1]\to\mathbb{R}.$ This probability measure is called the stationary probability measure associated to $\Phi$ and $\mathbf{p}.$ Its support is the attractor of $\Phi.$

Iterated function systems were first studied by Hutchinson [17] who proved the existence and uniqueness of the attractor and stationary probability measure. Since the late eighties, the problem of algorithmically estimating stationary probability measures has been a topic of active research, see [24, 2, 12, 20, 13].

In these notes we focus on the related problem of algorithmically estimating the integrals of functions with respect to stationary measures. In many areas of analysis, dynamical systems and probability theory, important quantities can be expressed in terms of such integrals, see section 6. For example, in fractal geometry, the Hausdorff dimension of the attractor $\Lambda$ can be bounded or sometimes explicitly calculated in terms of an expression involving the Lyapunov exponent of particular stationary measures, see section 6.3. On the other hand, in dynamical systems, equilibrium measures often take the form of a stationary probability measure, and therefore it is a natural problem to study the approximation of integrals with respect to them. Apart from some special cases, for instance whenever the stationary measure is absolutely continuous with an explicit known density, it is impossible to analytically calculate these integrals, and so one must turn to approximating them. We use an operator theoretic approach similar to [21, Theorem 1], where Jenkinson and Pollicott constructed an algorithm yielding accelerated approximations of integrals of functions with respect to Lebesgue measure. Indeed, our main Theorem 2.3 can be seen as an analogue of their result. In particular we provide

• •

a sufficient condition on the iterated function system $\Phi,$ and

• •

a class $\mathcal{C}_{\Phi}$ of functions $g:[0,1]\to\mathbb{R},$

such that the integral

[TABLE]

can be approximated efficiently, where $\mu$ denotes any stationary probability measure associated to $\Phi$. In particular, in Theorem 2.3 we provide an algorithm which produces approximations $\mu_{k}(g)$ to $\int g\textup{d}\mu$ such that the approximations approach the integral at a super-exponential rate. We demonstrate the performance of this algorithm by using it to estimate Hausdorff moments, Wasserstein distances and Lyapunov exponents.

2. Main results

We introduce the following contraction condition motivated by the definition introduced in [3, Definition 2.2].

Definition 2.1**.**

A family of maps $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$, $\phi_{i}:[0,1]\to[0,1]$ is called complex contracting if there exists a non-empty, bounded, connected and open set $\mathcal{D}\subset\mathbb{C}$ such that $[0,1]\subset\mathcal{D}$ and

1. (1)

each $\phi_{i}$ extends holomorphically to $\mathcal{D}$ 2. (2)

$\phi_{i}^{\prime}$* is continuous on the boundary of $\mathcal{D}$* 3. (3)

$\sup_{i\in\mathcal{I}}\sup_{z\in\mathcal{D}}|\phi_{i}^{\prime}(z)|<1.$**

The significance of the above definition is that if $\Phi$ is a complex contracting family of maps then there exists some complex domain $D\subset\mathcal{D}$ such that $[0,1]\subset D$, each $\phi_{i}$ extends holomorphically to $D$ and $\overline{\phi_{i}(D)}\subset D$, see [3, Lemma 2.4]. In this case we say that $D$ is an admissable domain for $\Phi$. The existence of such a domain is what is actually required in order to construct our algorithm, however in many cases it is easier to check that a family is complex contracting than to check that such a domain $D$ exists.

Definition 2.2**.**

Given a complex contracting iterated function system $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$ on the unit interval, we define the class of functions

[TABLE]

where the union is taken over all admissable domains for $\Phi.$

In particular if $\Phi$ is complex contracting, then all sufficiently small Euclidean $\varepsilon$-neighbourhoods of $[0,1]$ are admissable, see [3, Remark 2.5]. Therefore any function which is bounded and holomorphic on some Euclidean $\varepsilon$-neighbourhood of $[0,1]$ belongs to $\mathcal{C}_{\Phi}$.

We introduce some notation. Let $\mathcal{I}^{n}=\{i_{1}\cdots i_{n}:i_{j}\in\mathcal{I}\}$ and $\phi_{i_{1}\ldots i_{n}}:=\phi_{i_{1}}\circ\cdots\circ\phi_{i_{n}}$ for $n\in\mathbb{N}.$ Define $\mathcal{I}^{\ast}=\bigcup_{n\in\mathbb{N}}\mathcal{I}^{n}$. For $\mathbf{i}=i_{1}\ldots i_{n}\in\mathcal{I}^{\ast}$ and $1\leqslant k\leqslant n-1$ we define

[TABLE]

Define $\Pi:\mathcal{I}^{\mathbb{N}}\to\Lambda$ by

[TABLE]

For each $\mathbf{i}\in\mathcal{I}^{\ast}$ let $(\mathbf{i})^{\infty}$ denote the periodic point $(\mathbf{i})^{\infty}=\mathbf{i}\mathbf{i}\mathbf{i}...\in\mathcal{I}^{\mathbb{N}}$ and let $z_{\mathbf{i}}$ denote the fixed point $z_{\mathbf{i}}=\Pi((\mathbf{i})^{\infty})$ of $\phi_{\mathbf{i}}$. Finally, we say that an iterated function system is non-overlapping if $\{\phi_{i}((0,1))\}_{i\in\mathcal{I}}$ are pairwise disjoint.

The following is our main result.

Theorem 2.3**.**

Let $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$ be a complex contracting iterated function system on the unit interval and $\mathbf{p}=(p_{1},\ldots,p_{N})$ a probability vector. Given $g\in\mathcal{C}_{\Phi},$ define

[TABLE]

and

[TABLE]

Also define

[TABLE]

and

[TABLE]

Then for

[TABLE]

and $\mu$ the stationary probability measure associated to $\Phi$ and $\mathbf{p},$

[TABLE]

for some constants $C,\lambda>0$ which are independent of $k.$

Note that while our approximations are non-effective, in some special cases the constants $C,\lambda$ appearing in (18) can be made explicit. For example, when each $\phi_{i}$ takes the form of a linear fractional transformation, that is, $\phi_{i}(x)=\frac{a_{i}x+b_{i}}{c_{i}x+d_{i}}$ for some constants $a_{i},b_{i},c_{i},d_{i}$ (which in particular includes the important class of self-similar iterated function systems), the constants $C$ and $\lambda$ can be bounded similarly to [19, 22]. Also note that since the time taken to process $n$ steps of the algorithm is exponential in $n$, the error decreases super-polynomially fast in time.

Finally, we observe that by using stationarity of the measure $\mu$, Theorem 2.3 can also be used to obtain approximations of integrals of some piecewise analytic functions. In particular, suppose that for some $K\in\mathbb{N}$, $\{\phi_{\mathbf{i}}([0,1])\}_{\mathbf{i}\in\mathcal{I}^{K}}$ are pairwise disjoint and let $g:[0,1]\to\mathbb{R}$ be a function (not necessarily continuous) such that each composition function $g\circ\phi_{\mathbf{i}}:[0,1]\to\mathbb{R}$ belongs to $\mathcal{C}_{\Phi}$ for every $\mathbf{i}\in\mathcal{I}^{K}.$ By continuous continuation there is a continuous function $\tilde{g}:[0,1]\to\mathbb{R}$ with $\tilde{g}|_{\phi_{\mathbf{i}}([0,1])}=g|_{\phi_{\mathbf{i}}([0,1])}$ for every $\mathbf{i}\in\mathcal{I}^{K}.$ Thus by stationarity,

[TABLE]

so since $\phi_{\mathbf{i}}(D)\subset D$, Theorem 2.3 can be applied to approximate each integral $\int g\circ\phi_{\mathbf{i}}\textup{d}\mu$.

3. Preliminaries

3.1. Trace class operators, determinants and approximation numbers

Given a compact operator $L:H\to H$ on a Hilbert space $H$, its $n$th approximation number is defined as

[TABLE]

A bounded linear operator on a complex separable Hilbert space $H$ is called trace-class if $\sum_{n=1}^{\infty}s_{n}(L)<\infty.$ Given a trace-class operator $L$, the trace is defined as

[TABLE]

where $\{e_{n}\}$ is any orthonormal basis and $\langle,\rangle_{H}$ is the inner product for the Hilbert space $H$.

Given a compact operator $L$, we denote by $\{\lambda_{n}(L)\}_{n\in\mathbb{N}}$ the monotone decreasing sequence of non-zero eigenvalues of $L$, listed with algebraic multiplicity. If $L$ is trace-class then it is compact and its sequence of eigenvalues $\lambda_{n}(L)$ is absolutely summable.

For a trace-class operator $L$, the Fredholm determinant of $L$ can be defined as

[TABLE]

which is an entire function of $z$ [5, Theorem 3.3], so in particular there exist $a_{n}\in\mathbb{C}$ such that

[TABLE]

Note that by (3) the roots of $\det(\textnormal{Id}-zL)$ are precisely the reciprocals of the eigenvalues of $L$, and the degree of each zero is given by the multiplicity of the corresponding eigenvalue. Moreover, each coefficient $a_{n}$ can be expressed in terms of the traces of $L^{m}$ for $1\leqslant m\leqslant n$:

[TABLE]

for a proof see for instance [19, Proposition 3.2].

On the other hand, by finding the coefficient of $z^{n}$ in (3) we see that

[TABLE]

therefore

[TABLE]

3.2. Bergman space

We will use results of Bandtlow and Jenkinson [3, 4] on operators acting on the Bergman space. Let $D$ be any connected, non-empty open subset of $\mathbb{C}$. We define the Bergman space $A^{2}(D)$ by

[TABLE]

where $V$ denotes 2-dimensional Lebesgue measure, normalised so that the unit ball has unit mass. $A^{2}(D)$ is a Hilbert space with inner product

[TABLE]

Consider the operator

[TABLE]

where $\psi_{i}$ are holomorphic functions on $D$ that satisfy $\bigcup_{i=1}^{k}\overline{\psi_{i}(D)}\subset D$ and $w_{i}$ are holomorphic bounded functions on $D$. Using the work of Bandtlow and Jenkinson [3, 4] we have the following result.

Proposition 3.1**.**

The operator $L$ preserves the Bergman space $A^{2}(D)$. Moreover

1. (1)

$L$* is a trace-class operator and there exist constants $C,\lambda>0$ which are independent of $N$ and $w_{i}$ such that its eigenvalues satisfy*

[TABLE]

where $\lVert w_{i}\rVert_{\infty}:=\sup_{z\in D}|w_{i}(z)|$, 2. (2)

the trace of $L^{n}$ is given by

[TABLE]

where $z_{\mathbf{i}}$ is the unique fixed point of $\psi_{\mathbf{i}}$ and $w_{\mathbf{i}}(z):=w_{i_{n}}(z)w_{i_{n-1}}(\phi_{i_{n}}z)\cdots w_{i_{1}}(\phi_{i_{2}\ldots i_{n}}z).$

Proof.

The first part is a special case of [4, Theorems 5.9 and 5.13]. The second part follows by [3, Theorem 4.2]. ∎

3.3. Analytic perturbation theory

We say that a bounded linear operator $L$ on a Banach space has spectral gap if $L=\lambda P+N$ where $P$ is a rank one projection (so $P^{2}=P$ and $\dim(\textnormal{Im}(P))=1$), $N$ is a bounded operator with spectral radius $\rho(N)<|\lambda|$ and $PN=NP=0$. $L$ does not need to be compact in order to have a spectral gap, however if the operator $L$ is compact and has a simple leading eigenvalue111Throughout the paper we say that an eigenvalue is simple if it is algebraically simple, that is, the eigenvalue has a one-dimensional generalised eigenspace. and no other eigenvalues with the same absolute value, it has a spectral gap.

We can use the standard techniques of perturbation theory [18] to relate $\int g\textup{d}\mu$ to the spectral properties of an appropriate operator. The following perturbation theorem is presented in a more general form in [16, Theorem 3.8].

Theorem 3.2** (Analytic perturbation theorem).**

Let $\{L_{t}\}_{t\in\mathbb{C}}$ be a family of bounded linear operators on a Banach space such that $t\mapsto L_{t}$ is holomorphic and $L_{0}$ has spectral gap. Then there exists an open neighbourhood $U\subset\mathbb{C}$ of 0 for which $L_{t}$ has spectral gap for all $t\in U$. Moreover there exist $\lambda(t),P_{t},N_{t}$ which are holomorphic families on $U$ such that:

1. (a)

$L_{t}=\lambda(t)P_{t}+N_{t}$, 2. (b)

$N_{t}P_{t}=P_{t}N_{t}=0$** 3. (c)

$P_{t}$* is a bounded rank one projection and has the form*

[TABLE]

for some small circle $\gamma$ around $\lambda$ which separates it from the rest of the spectrum of $L_{0}$, 4. (d)

$\rho(N_{t})<|\lambda(t)|-\varepsilon$* for some $\varepsilon>0$ which is independent of $t$.*

4. Transfer operator

Let $\Phi$ and $g$ satisfy the hypothesis of theorem 2.3. By assumption, there exists an admissable domain $D$ for $\Phi$ such that $g$ has a bounded and holomorphic extension to $D$. Abusing notation slightly, we also denote this extension by $g$. Fix a probability vector $\mathbf{p}$, $s\in\mathbb{C}$ and $f\in A^{2}(D)$. Define the operator $\mathcal{L}_{s}$ as

[TABLE]

Then by Proposition 3.1, $\mathcal{L}_{s}f\in A^{2}(D)$. Observe that the iterates of $\mathcal{L}_{s}$ are given by

[TABLE]

The following proposition summarises some of the spectral properties of $\mathcal{L}_{s}$.

Proposition 4.1**.**

$\mathcal{L}_{s}:A^{2}(D)\to A^{2}(D)$* is a trace-class operator with decreasing sequence of eigenvalues $\{\lambda_{n}(s)\}_{n\in\mathbb{N}}$. Moreover*

1. (1)

There exist constants $C,\lambda>0$ which are independent of $s$ such that the eigenvalues admit the bound

[TABLE]

for any $s\in\mathbb{C}$ with $|s|\leqslant 1$, 2. (2)

$\lambda_{1}(0)=1$* is a simple eigenvalue of $\mathcal{L}_{0}$ of maximum modulus,* 3. (3)

$\lambda_{1}(s)$* is analytic in $s$ in a neighbourhood of 0 and* 4. (4)

\int g\textup{d}\mu=\frac{\textup{d}}{\textup{d}s}\lambda_{1}(s)\biggr{|}_{s=0}.

Proof.

The fact that $\mathcal{L}_{s}$ is trace-class and the uniform bounds on $\lambda_{n}(\mathcal{L}_{s})$ follow from Proposition 3.1.

Next we prove (2). Clearly 1 is an eigenvalue of $\mathcal{L}_{0}$ with eigenfunction $\mathbf{1}$. To see that it is geometrically simple, suppose that $f\in A^{2}(D)$ is a fixed point of $\mathcal{L}_{0}$ and $f\neq 0$. We will show that $f$ must be a constant function. First, observe that

[TABLE]

where the right hand side is finite because $\bigcup_{i=1}^{k}\overline{\phi_{i}(D)}$ is a compact subset of $D$. Therefore,

[TABLE]

for some $z_{0}\in\bigcup_{i=1}^{k}\overline{\phi_{i}(D)}$. By the maximum-modulus principle, $f$ is constant on $D$. Using the same argument, we can establish that 1 is the only eigenvalue of modulus 1.

Therefore it remains to show that $1$ is an algebraically simple eigenvalue. We need to show that $\ker(\mathcal{L}_{0}-\textnormal{Id})^{2}$ is one dimensional (so only consists of the constant functions). For a contradiction suppose that there exists $f\in A^{2}(D)$ for which $(\mathcal{L}_{0}-\textnormal{Id})f\neq 0$ but $(\mathcal{L}_{0}-\textnormal{Id})f\in\ker(\mathcal{L}_{0}-\textnormal{Id})$. So in particular $(\mathcal{L}_{0}-\textnormal{Id})f=c\mathbf{1}$ for some constant $c$. In particular, $c\neq 0$ since $(\mathcal{L}_{0}-\textnormal{Id})f\neq 0$ and therefore by replacing $f$ by $c^{-1}f$ we obtain that $(\mathcal{L}_{0}-\textnormal{Id})f=\mathbf{1}$, that is, $\mathcal{L}_{0}f=\mathbf{1}+f$. By induction we see that

[TABLE]

On the other hand, define

[TABLE]

and define $\Gamma=\bigcap_{n=1}^{\infty}\Gamma^{n}$. Since $\Gamma^{n}$ is a nested sequence of closed subsets of $D$, $\Gamma$ is a compact subset of $D$. For any $z\in\Gamma$,

[TABLE]

By (6), $|\mathcal{L}_{0}^{n}f(z)|=|n+f(z)|\geqslant n-|f(z)|$ implying that

[TABLE]

which is clearly a contradiction since $f$ is bounded on $\Gamma$.

To see (3), since $\mathcal{L}_{0}$ is compact, $\lambda_{1}(0)$ is a simple eigenvalue of maximum modulus and $\{\mathcal{L}_{s}\}$ is clearly an analytic family of operators in $s$, we can apply the analytic perturbation theorem 3.2 to deduce that $\lambda_{1}(s)$ is analytic in a neighbourhood $U$ of 0.

Finally to prove (4) put

[TABLE]

as in (c) of Theorem 3.2. (a)-(c) of Theorem 3.2 imply that the image of $\mathcal{P}_{s}$ is an eigenspace for the eigenvalue $\lambda_{1}(s)$ and that $h_{s}=\mathcal{P}_{s}\mathbf{1}$ is an eigenfunction for the eigenvalue $\lambda_{1}(s)$. Note that $h_{0}=\mathbf{1}$. Since $s\mapsto\mathcal{P}_{s}$ is holomorphic it immediately follows that $s\mapsto h_{s}$ is also holomorphic. We write f_{0}=\frac{\textup{d}}{\textup{d}s}h_{s}\bigr{|}_{s=0}\in A^{2}(D).

Fix some $z_{0}\in\Gamma\cap(0,1)$. Observe that for each $n>1$ and $s\in U$,

[TABLE]

Denote $S_{n}g(\phi_{\mathbf{i}}z_{0})=g(\phi_{\mathbf{i}}z_{0})+\ldots g(\phi_{\tilde{\sigma}^{n}\mathbf{i}}z_{0})$. Differentiating at $s=0$ we obtain

[TABLE]

where we used that $\lambda_{1}(0)=1$ and $h_{0}=\mathbf{1}$. Therefore since $|\mathcal{L}_{0}^{n}f_{0}(z_{0})-f_{0}(z_{0})|\leqslant\sup_{z\in\Gamma}|f_{0}(z)|<\infty$

[TABLE]

The fact that $\lim_{n\to\infty}\frac{1}{n}\sum_{\mathbf{i}\in\mathcal{I}^{n}}p_{\mathbf{i}}S_{n}g(\phi_{\mathbf{i}}z_{0})=\lim_{n\to\infty}\sum_{\mathbf{i}\in\mathcal{I}^{n}}p_{\mathbf{i}}g(\phi_{\mathbf{i}}z_{0})$ follows because $\phi_{i}$ are uniformly contracting and $g$ has bounded derivative on $[0,1]$, therefore for any $\varepsilon,\delta>0$ one can choose $N$ sufficiently large so that for $n\geqslant N$, and all $k\leqslant(1-\varepsilon)n$,

[TABLE]

∎

5. Determinants and the algorithm

Since $\mathcal{L}_{s}$ is a trace-class operator for each $s\in\mathbb{C}$, we can define its determinant function $\det(\textnormal{Id}-z\mathcal{L}_{s})$. The following proposition summarises its properties.

Proposition 5.1**.**

For all $s\in\mathbb{C}$, $\det(\textnormal{Id}-z\mathcal{L}_{s})$ is an entire function of $z$ and can be written in the form

[TABLE]

for $b_{n}(s)\in\mathbb{C}$ where $b_{0}(s)=1$ for all $s$ and for $n\geqslant 1$ is defined as

[TABLE]

Moreover, the trace of $\mathcal{L}_{s}^{n}$ is given by

[TABLE]

Proof.

(9) follows from (4). To see (10) notice that by Proposition 3.1(2),

[TABLE]

Since for any $1\leqslant k\leqslant n$, $\phi_{\tilde{\sigma}^{k}\mathbf{i}}(z_{\mathbf{i}})$ is a fixed point of $\phi_{\sigma^{k}\mathbf{i}}$ it follows that $\phi_{\tilde{\sigma}^{k}\mathbf{i}}(z_{\mathbf{i}})=z_{\sigma^{k}\mathbf{i}}$ and therefore (10) follows. ∎

Proposition 5.2**.**

Let $\det(\textnormal{Id}-z\mathcal{L}_{s})=\sum_{n=0}^{\infty}b_{n}(s)z^{n}$ as before. Then there exist constants $C,\lambda>0$ such that $|b_{n}(s)|,|b_{n}^{\prime}(0)|\leqslant C\exp(-\lambda n^{2})$ for all $s\in\mathbb{C}$, $|s|\leqslant 1$. Moreover,

[TABLE]

Proof.

For the first part, note that by (5) and Proposition 4.1(1),

[TABLE]

where the last line follows by repeated geometric summation. The upper bound on $|b_{n}(s)|$ follows (for new constants $C$ and $\lambda$). By the Cauchy integral formula,

[TABLE]

therefore

[TABLE]

from which we obtain the bound on $|b_{n}^{\prime}(0)|$.

Let $U$ be the neighbourhood of 0 on which $\lambda_{1}(s)$ is analytic. Observe that since the zeroes of the determinant $\det(\textnormal{Id}-z\mathcal{L}_{s})$ are the reciprocals of the eigenvalues of $\mathcal{L}_{s}$,

[TABLE]

Since $|b_{n}(s)|=O(\exp(-\lambda n^{2}))$ uniformly on $U$ we can apply the Cauchy integral formula to deduce that the partial sums $\sum_{n=1}^{N}b_{n}^{\prime}(s)-nb_{n}(s)\lambda_{1}^{\prime}(s)$ converge uniformly on compact subsets of $U$ as $N\to\infty$. Therefore we can differentiate (12) and take derivatives inside the summation to obtain

[TABLE]

Since $(\lambda_{1}(0))^{-1}$ is a simple zero of $\det(\textnormal{Id}-z\mathcal{L}_{0})$, it follows that $\sum_{n=0}^{\infty}nb_{n}(0)\neq 0$ and so by rearranging (13) we obtain

[TABLE]

which completes the proof. ∎

We define the $k$th approximation

[TABLE]

Define

[TABLE]

and

[TABLE]

We are now ready to prove the main theorem.

Proof of Theorem 2.3. By (9), (14) and the definitions of $t_{n},\tau_{n}$ it follows that $\mu_{k}(g)$ is given by theorem 2.3. By the bounds on $b_{n}(0)$ and $b_{n}^{\prime}(0)$ given in Proposition 5.2 and (11) we obtain the desired bound on the error $|\int g\textup{d}\mu-\mu_{k}(g)|$. ∎

6. Applications

In this section we demonstrate Theorem 2.3 by estimating moments and Lyapunov exponents of stationary probability measures and the first Wasserstein distance between stationary probability measures.

6.1. Moments of stationary probability measures

The $n$-th moment of a positive Borel measure $\mu$ on $\mathbb{R}$ is defined by

[TABLE]

Moments are particularly important in analysis, probability and statistics. In analysis, the classical moment problem [1] is connected with characterising the image of the map

[TABLE]

which is connected to the problem of the extension of positive functionals. In probability and statistics, the method of moments [6] is useful for proving limit theorems and estimating distributions of samples.

Since moments of stationary probability measures can be computed analytically in some special cases (see Lemma 6.2), we can demonstrate the efficiency of the algorithm by using it to approximate moments.

Corollary 6.1**.**

Let $\Phi$ and $\mathbf{p}$ satisfy the assumptions of Theorem 2.3. Then, for each fixed $m\in\mathbb{N}$, Theorem 2.3 provides an algorithm with a super-exponential rate of convergence (18) that gives approximations $w_{k}\in\mathbb{R}^{m+1}$ to $\gamma_{[0,m]}:=(\gamma_{0},\gamma_{1},\ldots,\gamma_{m})$ for $k=1,2,\ldots.$

In the setting where $\Phi$ is a non-overlapping iterated function system of similarities it is possible to obtain the following analytic formula for the moments.

Lemma 6.2**.**

Let $\Phi=(\phi_{1},\ldots,\phi_{N}),$ $\phi_{i}(x)=\rho_{i}x+t_{i}$ with $0<\rho_{i}<1$ such that $\phi_{i}(0,1)\cap\phi_{j}(0,1)=\emptyset$ for every $i\neq j$ and $\mathbf{p}=(p_{1},\ldots,p_{N})$ a probability vector.

Then $\gamma_{0}=1$ and for every $n>0$

[TABLE]

Proof.

Directly from (17) and definition of $\Phi.$ ∎

We can use Theorem 2.3 to compute approximate values for the moments and compare these with the exact values given by Lemma 6.2.

Example 6.3**.**

Let $\Phi=\{\frac{1}{3}x,\frac{1}{3}x+\frac{2}{3}\}$ and $p=\left(\frac{1}{3},\frac{2}{3}\right).$ In table 1 we compare the approximate values for the moments with the exact values given by the formula.

On the other hand, when $\Phi$ is made up of non-affine contractions, there is no analytic formula available. In the following example we approximate the moments of a stationary probability measure in this setting.

Example 6.4**.**

Let $\Phi=\{\frac{1}{x+2},\frac{1}{x+4}\}$ and $\mathbf{p}=\left(\frac{1}{2},\frac{1}{2}\right).$ We can compute the approximate values for the moments. The results are in Table 2, where all of the digits are stable.222We say that a digit is stable if, from empirical observations, it appears to have converged to a stable value

6.2. First Wasserstein distance between stationary probability measures

Fraser [11] initiated the study of the interaction between fractal geometry and optimal transportation by proposing the problem of computing the Wasserstein distances between stationary probability measures. The Wasserstein distance is a metrization of the weak-$*$ topology of the space of probability Borel measures on a Polish space. For $m\in[1,\infty),$ the Wasserstein distance of order $m$ between two probability measures $\mu$ and $\nu$ is defined by

[TABLE]

In general, under no assumption on the iterated functions system and associated stationary probability measure, there is no hope of finding an analytic formula. Under some very restrictive assumptions it has been proved in [11, 8] that $W_{1}(\mu,\nu)$ can be obtained explicitly. Using results announced in [7] and Theorem 2.3 we can provide an estimate of $W_{1}(\mu,\nu)$ when $d=1$ and $\mu,\nu$ are stationary probability measures. The following result was obtained in [7].

Theorem 6.5**.**

Let $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$ be a non-overlapping iterated function system of differentiable Lipschitz contractions on the unit interval. Additionally, suppose that $\frac{d\phi_{i}}{dx}>0$ for each $i\in\mathcal{I}$ and that $(\mathbf{p},\mathbf{q})$ is a pair of probability vectors with the property that

[TABLE]

does not change sign for $i=1,\ldots,N$. Then

[TABLE]

Using the above result and Theorem 2.3, we obtain the following.

Corollary 6.6**.**

Let $\Phi$ and $\mathbf{p},\mathbf{q}$ satisfy the hypothesis of Theorem 6.5 and Theorem 2.3. Then Theorem 2.3 provides an algorithm that gives approximations $w_{k}$ $(k=1,2,\ldots)$ to $W_{1}(\mu^{(\Phi,\mathbf{p})},\mu^{(\Phi,\mathbf{q})})$ and which converges at a super-exponential rate (18).

Let $\Phi,\mathbf{p},\mathbf{q}$ satisfy the assumptions of Corollary 6.6 and additionally suppose that $\phi_{i}(x)=\rho_{i}x+t_{i}$ with $0<\rho_{i}<1$ for each $i=1,\ldots N$. In [7] it was proved that

[TABLE]

We can compare the approximate values for the Wasserstein distances given by Theorem 6.6 with the exact value given by the equation (15).

Example 6.7**.**

Let $\Phi=\left\{\frac{1}{3}x,\frac{1}{2}x+\frac{1}{2}\right\}$ and $p=\left(\frac{1}{3},\frac{2}{3}\right),$ $q=\left(\frac{3}{4},\frac{1}{4}\right).$ We compare the approximate value of the Wasserstein distance $w_{k}$ for $k=8,9,\ldots,16,$ with the exact value $W_{1}\left(\mu^{(\Phi,\mathbf{p})},\mu^{(\Phi,\mathbf{q})}\right)=\frac{2}{5}.$ The results are in table 3.

On the other hand, when the contractions in $\Phi$ are non-affine, an analytic formula is not available.

Example 6.8**.**

Let $\Phi=\left\{\frac{\sin(\pi x/4)}{6}+\frac{1}{4},\frac{\sin(\pi x/4)}{3}+\frac{2}{3}\right\}$ and $\mathbf{p}=\left(\frac{1}{7},\frac{6}{7}\right),$ $\mathbf{q}=\left(\frac{1}{2},\frac{1}{2}\right).$ We can compute the approximate values for the Wasserstein distance between the stationary probability measures associated to $(\Phi,\mathbf{p})$ and $(\Phi,\mathbf{q}).$ The results are in Table 4.

6.3. Lyapunov exponents

Let $\Phi$ be an iterated function system and $\mu$ be the stationary probability measure associated to $(\Phi,\mathbf{p}).$ The Lyapunov exponent of $\Phi$ with respect to $\mu$ can be defined by

[TABLE]

The Lyapunov exponent of $\mu$ describes the typical rate at which distances are contracted under the action of the maps in the iterated function system, from the point of view of the measure $\mu$. Therefore, Lyapunov exponents play an important role in fractal geometry; for example the Hausdorff dimension of $\mu$ is given in terms of its Lyapunov exponent, which in turn can provide bounds or even precise values for the Hausdorff dimension of the attractor $\Lambda$, see [9]. When $\Phi$ is non-overlapping, the Lyapunov exponent of $\mu$ has an equivalent interpretation as the $\mu$-typical rate of expansion of the dynamical system whose branches are given by the inverse images of the maps in $\Phi$, and therefore Lyapunov exponents also play an important role in ergodic theory, particularly in thermodynamic formalism. We can use Theorem 2.3 to approximate the Lyapunov exponent of $\Phi$ with respect to $\mu$.

Corollary 6.9**.**

Let $\mathbf{p}=(p_{1},\ldots,p_{N})$ a probability vector and let $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$ be a complex contracting iterated function system on the unit interval such that:

[TABLE]

Then Theorem 2.3 provides an algorithm that gives approximations $w_{k}$ $(k=1,2,\ldots)$ to the Lyapunov exponent of $\Phi$ with respect to $\mu$ which converges at the super-exponential rate (18).

Proof.

Follows by applying Theorem 2.3 to the function $g(x)=-\sum_{i=1}^{N}p_{i}\log|\phi_{i}^{\prime}(x)|$ which belongs to $\mathcal{C}_{\Phi}$ by (16) and the fact that $\Phi$ is complex contracting. ∎

Example 6.10**.**

Let $\Phi=\{\frac{\sin(\pi x/4)}{6}+\frac{1}{4},\frac{\sin(\pi x/4)}{3}+\frac{2}{3}\}$ and $\mathbf{p}=\left(\frac{1}{3},\frac{2}{3}\right).$ We compute the approximate values to the Lyapunov exponent of $\Phi$ with respect to $\mu$. The results are in Table 5.

7. Further remarks

7.1. Implementation

The algorithm was implemented in Python using the free library mpmath 333Python library for real and complex floating-point arithmetic with arbitrary precision.. To estimate the fixed points $z_{\mathbf{i}}$ with arbitrary precision the command

mpmath.findroot

was used and to estimate $\phi^{\prime}_{\mathbf{i}}(z_{\mathbf{i}})$ with arbitrary precision the command

mpmath.diff

was used.

7.2. Other applications

Our methods can be used to estimate other functions in probability theory, such as pointwise estimates of the Fourier-Stieltjes transform of stationary probability measures $\hat{\mu}(x):=\int_{0}^{1}e^{-2\pi ixy}\textup{d}\mu(y)$ and coefficients of the Fourier series of singular measures arising as stationary probability measures.

7.3. Extension to non constant weight functions

Let $\Phi$ be a complex contracting iterated function as before. We can extend our results to certain probability measures associated to non-constant weight functions.

In particular, suppose there exists an admissable domain $D$ for $\Phi$ and bounded holomorphic functions $\mathbf{p}=\{p_{i}\}_{i\in\mathcal{I}}$, $p_{i}:D\to\mathbb{C}$ with the property that

1. (1)

$\sum_{i\in\mathcal{I}}p_{i}(z)=1$ for all $z\in D$, 2. (2)

$p_{i}([0,1])\subset(0,1)$ for all $i\in\mathcal{I}$.

Then we can obtain an analogue of Theorem 2.3 for the unique probability measure $\mu=\mu^{(\Phi,\mathbf{p})}$ supported on the attractor of $\Phi$ with the property that

[TABLE]

for every continuous function $\varphi:[0,1]\to\mathbb{R}$ . This probability measure is called the stationary probability measure associated to $\Phi$ and $\mathbf{p}$, see [10]. It generalises the stationary probability measure considered for constant weight functions.

By following the proof of Theorem 2.3 with the operator

[TABLE]

one can obtain the following analogue of Theorem 2.3.

Theorem 7.1**.**

Let $\Phi=\{\phi_{i}\}_{i\in\mathcal{I}}$, $\mathbf{p}=\{p_{i}\}_{i\in\mathcal{I}}$ and $\mu=\mu^{\Phi,\mathbf{p}}$ be as above. Given $g\in\mathcal{C}_{\Phi},$ define

[TABLE]

and

[TABLE]

where

[TABLE]

Define $\alpha_{n}$, $a_{n}$ and $\mu_{k}(g)$ as in Theorem 2.3. We have

[TABLE]

for some constants $C,\lambda>0$ which are independent of $k.$

Acknowledgements.

IC would like to express his gratitude to Felipe Serrano and Franziska Schlösser for suggestions around the effective implementation of the algorithms. IC was partially supported by CONICYT PIA ACT172001. NJ was financially supported by the Leverhulme Trust (Research Project Grant number RPG-2016-194). Both authors would like to thank Ian Morris for helpful discussions and suggestions.