G)_{*}\overline{X})$

(where $\overline{X}$ is the $\boldsymbol{X}$-controlled path with trace $X$
and zero Gubinelli derivatives) and the conclusion is implied by 2. below.

The proof of 2. is straightforward (with the second claim deduced from the Faà
di Bruno formula 59).

As for 3., we have

$\displaystyle((\overline{J}\cdot\overline{K})*\overline{H}))_{\boldsymbol{\gamma}}^{a\choose
b}$
$\displaystyle=\\!\\!\\!\\!\sum_{\begin{subarray}{c}m=1,\ldots,|{\boldsymbol{\gamma}}|\\\
({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\end{subarray}}\\!\\!\\!\\!(\overline{J}\cdot\overline{K})^{a\choose
b}_{\boldsymbol{k}}\overline{H}^{k_{1}}_{{\boldsymbol{\gamma}}^{1}}\cdots\overline{H}^{k_{m}}_{{\boldsymbol{\gamma}}^{m}}$
$\displaystyle=\\!\\!\\!\\!\sum_{\begin{subarray}{c}m=1,\ldots,|{\boldsymbol{\gamma}}|\\\
({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\\\
({\boldsymbol{i}},{\boldsymbol{j}})\in\mathrm{Sh}^{-1}({\boldsymbol{k}})\end{subarray}}\\!\\!\\!\\!\overline{J}^{a\choose
c}_{{\boldsymbol{i}}}\overline{K}^{c\choose
b}_{\boldsymbol{j}}\overline{H}^{k_{1}}_{{\boldsymbol{\gamma}}^{1}}\cdots\overline{H}^{k_{m}}_{{\boldsymbol{\gamma}}^{m}}$
$\displaystyle=\\!\\!\\!\\!\sum_{\begin{subarray}{c}l+q=1,\ldots,|{\boldsymbol{\gamma}}|\\\
({\boldsymbol{\alpha}},{\boldsymbol{\beta}})\in\mathrm{Sh}^{-1}({\boldsymbol{\gamma}})\\\
({\boldsymbol{\alpha}}^{1},\ldots,{\boldsymbol{\alpha}}^{l})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\alpha}})\\\
({\boldsymbol{\beta}}^{1},\ldots,{\boldsymbol{\beta}}^{q})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\beta}})\end{subarray}}\\!\\!\\!\\!\overline{J}^{a\choose
c}_{{\boldsymbol{i}}}\overline{K}^{c\choose
b}_{\boldsymbol{j}}\overline{H}^{i_{1}}_{{\boldsymbol{\alpha}}^{1}}\cdots\overline{H}^{i_{l}}_{{\boldsymbol{\alpha}}^{l}}\overline{H}^{j_{1}}_{{\boldsymbol{\beta}}^{1}}\cdots\overline{H}^{j_{q}}_{{\boldsymbol{\beta}}^{q}}$
$\displaystyle=(\overline{J}*\overline{H})\cdot(\overline{K}*\overline{H})$

where the second last identity follows from Section 1 (with $m=2$).

As for 4., it is easy to show, using associativity of composition and of
$\Delta_{\shuffle}$, that both sides coincide with
$\times^{3}\circ(\overline{J}\boxtimes\overline{K}\boxtimes\overline{H})\circ\Delta_{\shuffle}^{3}$
where $\times^{3}$ denotes composition of three linear maps.

Finally, 5. is shown as follows:

$\displaystyle F^{*}G^{*}\overline{H}$
$\displaystyle=(G^{*}\overline{H}*\overline{F}(X))\cdot\overline{DF}(X)$
$\displaystyle=(((\overline{H}*\overline{G}(F(X)))\cdot\overline{DG}(F(X)))*\overline{F}(X))\cdot\overline{DF}(X)$
$\displaystyle=(\overline{H}*\overline{G}(F(X))*\overline{F}(X))\cdot(\overline{DG}(F(X))*\overline{F}(X))\cdot\overline{DF}(X)$
$\displaystyle=((\overline{H}*\overline{G\circ F}(X))\cdot\overline{D(G\circ
F)}(X)$ $\displaystyle=(G\circ F)^{*}\overline{H}.$

Here, we have used the previous points 2., 3. and 4. in the proposition, as
well as the fact that

$\overline{G}(F(X))*\overline{F}(X)=G_{*}\overline{F}(X)=G_{*}F_{*}\overline{X}=(G\circ
F)_{*}\overline{X}=\overline{G\circ F}(X)$

and similarly that

$\displaystyle(\overline{DG}(F(X))*\overline{F}(X))\cdot\overline{DF}(X)$
$\displaystyle=\overline{DG\circ F}(X)\cdot\overline{DF}(X)$
$\displaystyle=\overline{DG\circ F(\cdot)DF(\cdot)}(X)$
$\displaystyle=\overline{D(G\circ F)}(X)$

where we have used Section 2. ∎

In the next theorem we prove the property, well-known in both ordinary and
stochastic calculus, which allows to substitute the differential. This will be
especially convenient when manipulating RDEs. We have only introduced the
theory necessary to handle weakly geometric rough paths, and the theorem is
therefore stated in this context, but one can expect this type of result to
also hold true in other settings, such as Itô calculus and branched rough
paths.

######  Theorem  (Associativity of the rough integral).

Let $\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],V)$,
$\overline{H}\in\mathscr{D}_{\boldsymbol{X}}(\mathcal{L}(V,W))$,
$I\coloneqq\int\overline{H}\emph{d}\boldsymbol{X}$, $\overline{I}$ and
$\boldsymbol{I}$ respectively the canonical controlled and rough paths above
$I$, $\overline{K}\in\mathscr{D}_{\boldsymbol{I}}(\mathcal{L}(W,U))$. Then

$\bigg{(}\overline{\int}\overline{K}\emph{d}\boldsymbol{I}\bigg{)}*\overline{I}=\overline{\int}(\overline{K}*\overline{I})\cdot\overline{H}\emph{d}\boldsymbol{X}$
(69)

and therefore

$\boldsymbol{\int}\overline{K}\emph{d}\boldsymbol{I}=\boldsymbol{\int}(\overline{K}*\overline{I})\cdot\overline{H}\emph{d}\boldsymbol{X}.$
(70)

###### Proof.

In this proof we will denote, for a tuple ${\boldsymbol{\gamma}}$,
$\gamma^{\cdot}$ its last entry and ${\boldsymbol{\gamma}}^{-}$ the tuple
obtained by removing $\gamma^{\cdot}$, so
${\boldsymbol{\gamma}}=({\boldsymbol{\gamma}}^{-},\gamma^{\cdot})$. Moreover,
we will interchangeably use the two indexing notations for controlled
integrands, e.g.
$\overline{H}^{k}_{\boldsymbol{\gamma}}=\overline{H}^{k\choose\gamma^{\cdot}}_{{\boldsymbol{\gamma}}^{-}}$.
For $|{\boldsymbol{\gamma}}|\geq 1$ we then have

$\displaystyle\bigg{(}\overline{\int}(\overline{K}*\overline{I})\cdot\overline{H}\mathrm{d}\boldsymbol{X}\bigg{)}^{c}_{\boldsymbol{\gamma}}$
$\displaystyle=((\overline{K}*\overline{I})\cdot\overline{H})^{c}_{\boldsymbol{\gamma}}$
$\displaystyle=\sum_{({\boldsymbol{\alpha}},{\boldsymbol{\beta}})\in\mathrm{Sh}^{-1}({\boldsymbol{\gamma}}^{-})}(\overline{K}*\overline{I})^{c\choose
h}_{\boldsymbol{\alpha}}\overline{H}^{h\choose\gamma^{\cdot}}_{\boldsymbol{\beta}}$
$\displaystyle=\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\alpha}},{\boldsymbol{\beta}})\in\mathrm{Sh}^{-1}({\boldsymbol{\gamma}}^{-})\\\
({\boldsymbol{\alpha}}^{1},\ldots,{\boldsymbol{\alpha}}^{n})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\alpha}})\end{subarray}}\\!\\!\overline{K}^{c\choose
h}_{\boldsymbol{k}}\overline{H}^{k_{1}}_{{\boldsymbol{\alpha}}^{1}}\cdots\overline{H}^{k_{n}}_{{\boldsymbol{\alpha}}^{n}}\overline{H}^{h\choose\gamma^{\cdot}}_{\boldsymbol{\beta}}$
$\displaystyle=\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\alpha}},{\boldsymbol{\beta}})\in\mathrm{Sh}^{-1}({\boldsymbol{\gamma}}^{-})\\\
({\boldsymbol{\alpha}}^{1},\ldots,{\boldsymbol{\alpha}}^{n})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\alpha}})\end{subarray}}\\!\\!\overline{K}^{c}_{({\boldsymbol{k}},h)}\overline{H}^{k_{1}}_{{\boldsymbol{\alpha}}^{1}}\cdots\overline{H}^{k_{n}}_{{\boldsymbol{\alpha}}^{n}}\overline{H}^{h}_{({\boldsymbol{\beta}},\gamma^{\cdot})}$
$\displaystyle=\\!\\!\sum_{\begin{subarray}{c}m=1,\ldots,|{\boldsymbol{\gamma}}|\\\
({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\end{subarray}}\\!\\!\overline{K}^{c}_{{\boldsymbol{k}}}\overline{H}^{k_{1}}_{{\boldsymbol{\gamma}}^{1}}\cdots\overline{H}^{k_{m}}_{{\boldsymbol{\gamma}}^{m}}$
$\displaystyle=\\!\\!\sum_{\begin{subarray}{c}m=1,\ldots,|{\boldsymbol{\gamma}}|\\\
({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\end{subarray}}\\!\\!\bigg{(}\overline{\int}\overline{K}\mathrm{d}\boldsymbol{I}\bigg{)}^{c}_{\boldsymbol{k}}\overline{I}^{k_{1}}_{{\boldsymbol{\gamma}}^{1}}\cdots\overline{I}^{k_{m}}_{{\boldsymbol{\gamma}}^{m}}$
$\displaystyle=\bigg{(}\bigg{(}\overline{\int}\overline{K}\mathrm{d}\boldsymbol{I}\bigg{)}*\overline{I}\bigg{)}_{\boldsymbol{\gamma}}^{c}.$

At the trace level, we have, through a similar argument

$\displaystyle\bigg{(}\overline{\int}(\overline{K}*\overline{I})\cdot\overline{H}\mathrm{d}\boldsymbol{X}\bigg{)}^{c}_{{\boldsymbol{\gamma}};st}$
$\displaystyle\approx\sum_{|{\boldsymbol{\gamma}}|=1}^{\lfloor
p\rfloor}((\overline{K}*\overline{I})\cdot\overline{H})^{c}_{{\boldsymbol{\gamma}};s}\boldsymbol{X}^{\boldsymbol{\gamma}}_{st}$
$\displaystyle=\\!\\!\sum_{\begin{subarray}{c}|{\boldsymbol{\gamma}}|=1,\ldots{\lfloor
p\rfloor}\\\ m=1,\ldots,|{\boldsymbol{\gamma}}|\\\
({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\end{subarray}}\\!\\!\overline{K}^{c}_{{\boldsymbol{k}};s}\overline{H}^{k_{1}}_{{\boldsymbol{\gamma}}^{1};s}\cdots\overline{H}^{k_{m}}_{{\boldsymbol{\gamma}}^{m};s}\boldsymbol{X}^{\boldsymbol{\gamma}}_{st}$
$\displaystyle\approx\int\overline{K}^{c}\mathrm{d}\boldsymbol{I}.$

As for the statement at the level of rough paths, we have

$\displaystyle\boldsymbol{\int}\overline{K}\mathrm{d}\boldsymbol{I}$
$\displaystyle=\Big{\uparrow}_{\\!\boldsymbol{I}}\overline{\int}\overline{K}\mathrm{d}\boldsymbol{I}$
$\displaystyle=\Big{\uparrow}_{\\!\boldsymbol{X}}\bigg{[}\bigg{(}\overline{\int}\overline{K}\mathrm{d}\boldsymbol{I}\bigg{)}*\overline{I}\bigg{]}$
$\displaystyle=\Big{\uparrow}_{\\!\boldsymbol{X}}\bigg{[}\overline{\int}(\overline{K}*\overline{I})\cdot\overline{H}\mathrm{d}\boldsymbol{X}\bigg{]}$
$\displaystyle=\boldsymbol{\int}(\overline{K}*\overline{I})\cdot\overline{H}\mathrm{d}\boldsymbol{X}$

where we have used 1. in Section 2 and the previous statement. ∎

The next result, for which geometricity is essential, tells us that $F_{*}$
and $F^{*}$ behave as adjoint operators under the rough integral pairing. Its
proof is an immediate consequence of Section 2 and Section 2.

######  Theorem  (Pushforward-pullback adjunction).

Let $\boldsymbol{X},\overline{H},F$ be as above, then

$\bigg{(}\overline{\int}\overline{H}\emph{d}F_{*}\boldsymbol{X}\bigg{)}*\overline{F}(X)=\overline{\int}F^{*}\overline{H}\emph{d}\boldsymbol{X}$
(71)

and therefore

$\boldsymbol{\int}\overline{H}\emph{d}F_{*}\boldsymbol{X}=\boldsymbol{\int}F^{*}\overline{H}\emph{d}\boldsymbol{X}.$
(72)

Next we move on to the topic of rough differential equations (RDEs). We will
introduce two equivalent notions of solution to an RDE. Given a field of
linear maps $F\in C^{\infty}(W,\mathcal{L}(V,W))$ and a smooth map $g\in
C^{\infty}(W,U)$ we define, for $y\in W$

$Fg(y)\coloneqq Dg(y)\circ F(y)\in\mathcal{L}(V,U)$ (73)

and inductively

$F^{n}g(y)\coloneqq F(\eta\mapsto
F^{n-1}g(\eta))|_{\eta=y}\in\mathcal{L}(V,\mathcal{L}(V^{\otimes
n-1},U))=\mathcal{L}(V^{\otimes n},U).$ (74)

In coordinates we denote

$F_{\gamma}g^{c}(y)\coloneqq
Fg(y)^{c}_{\gamma}=\partial_{k}g(y)F^{k}_{\gamma}(y)\ \Rightarrow\
F^{n}g(y)_{(\gamma_{1},\ldots,\gamma_{n})}^{c}=F_{\gamma_{1}}\cdots
F_{\gamma_{n}}g^{c}(y).$

We will also use the compact notation $F_{\boldsymbol{\gamma}}g^{c}(y)$ for
the latter.

######  Remark .

Note that $F_{\gamma_{1}}\cdots F_{\gamma_{n}}g^{c}(y)$ can be read right to
left as well as left to right, i.e. it is equal to $F_{\gamma_{1}}\cdots
F_{\gamma_{n-1}}(F_{\gamma_{n}}g^{c})(y)$ for $n\geq 2$. This can be seen by
induction on $n$ (with the quantifier $\forall g$ inside the inductive
hypothesis). For $n=2$ the statement is tautological. For the inductive step
we have

$\displaystyle F_{\gamma_{1}}\cdots F_{\gamma_{n+1}}g^{c}(y)$
$\displaystyle=F_{\gamma_{1}}(F_{\gamma_{2}}\cdots F_{\gamma_{n+1}}g^{c})(y)$
$\displaystyle=F_{\gamma_{1}}(F_{\gamma_{2}}\cdots
F_{\gamma_{n}}(F_{\gamma_{n+1}}g^{c}))(y)$ $\displaystyle=F_{\gamma_{1}}\cdots
F_{\gamma_{n}}(F_{\gamma_{n+1}}g^{c})(y)$

where in the second identity we have used the inductive hypothesis.

######  Definition  (Davie solution to an RDE).

Let $F\in C^{\infty}(W,\mathcal{L}(V,W))$. A _solution_ to the RDE

$\mathrm{d}Y=F(Y)\mathrm{d}\boldsymbol{X},\quad Y_{0}=y_{0}$ (75)

is a path $Y\in C([0,T],W)$ starting at $y_{0}$ with the property that for all
$g\in C^{\infty}(W)$.

$g(Y)_{st}\approx\sum_{n=1}^{\lfloor p\rfloor}\langle
F^{n}g(Y_{s}),\boldsymbol{X}^{n}_{st}\rangle.$ (76)

######  Proposition  (Gubinelli solution to an RDE).

Let $Y$ be a solution to 75. Then

$\overline{Y}\coloneqq(Y,F(Y)=F\mathbbm{1}(Y),\ldots,F^{{\lfloor
p\rfloor}-1}\mathbbm{1}(Y))\in\mathscr{D}_{\boldsymbol{X}}(W)$ (77)

and moreover

$\overline{Y}=y_{0}+\overline{\int}F_{*}\overline{Y}\mathrm{d}\boldsymbol{X}.$
(78)

Conversely, if an $\boldsymbol{X}$-controlled controlled path satisfying the
above identity, its trace satisfies Section 2.

In order to prove this proposition we will make use of the following

######  Lemma .

The $\boldsymbol{X}$-controlled path $\overline{Y}$ has the form 77 if and
only if $\overline{Y}_{n}=(F_{*}\overline{Y})_{n-1}$ for $n\geq 1$. Moreover,
in this case

$((Fg)_{*}\overline{Y})_{n-1}=F^{n}g(Y),\quad n\geq 1$

with $g$ as above.

###### Proof.

If we prove

$(Fg)_{*}(Y,F(Y),\ldots,F^{{\lfloor
p\rfloor}-1}\mathbbm{1}(Y))_{n-1}=F^{n}g(Y)$

we will have shown both the second statement and the only if part of the first
(just choose $g=\mathbbm{1}$). In order to show this, we first establish the
identity

$F_{\boldsymbol{\gamma}}g^{c}(y)=\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\\\
|{\boldsymbol{\gamma}}^{1}|,\ldots,|{\boldsymbol{\gamma}}^{m}|\geq
1\end{subarray}}\partial_{\boldsymbol{k}}g^{c}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}(y).$ (79)

Note that this is not a closed form formula for $F^{n}g$, since iterated
compositions of the vector fields $F_{\boldsymbol{\gamma}}$ also appear on the
right hand side, but it will be useful for us nonetheless. We proceed by
induction on $|{\boldsymbol{\gamma}}|$. For $|{\boldsymbol{\gamma}}|=1$ there
is nothing to show. For the inductive step, using Section 2 we have

$\displaystyle\mathrel{\phantom{=}}F_{(\gamma_{1},\ldots,\gamma_{n+1})}g^{c}(y)$
$\displaystyle=F_{\gamma_{1}}(F_{\gamma_{2},\ldots,\gamma_{n+1}}g^{c})(y)$
$\displaystyle=F_{\gamma_{1}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}(\gamma_{2},\ldots,\gamma_{n+1})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{\boldsymbol{k}}g^{c}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}(y)$
$\displaystyle=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}(\gamma_{2},\ldots,\gamma_{n+1})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\big{(}\partial_{h,{\boldsymbol{k}}}g^{c}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}$
$\displaystyle\mathrel{\phantom{=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}(\gamma_{2},\ldots,\gamma_{n+1})}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!}}+\sum_{l=1}^{m}\partial_{{\boldsymbol{k}}}g^{c}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots\partial_{h}(F_{{\boldsymbol{\gamma}}^{l}}\mathbbm{1}^{k_{l}})\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}\big{)}F^{h}_{\gamma_{1}}(y)$
$\displaystyle=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}(\gamma_{2},\ldots,\gamma_{n+1})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\big{(}\partial_{h,{\boldsymbol{k}}}g^{c}F^{h}_{\gamma_{1}}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}$
$\displaystyle\mathrel{\phantom{=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}(\gamma_{2},\ldots,\gamma_{n+1})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!}}+\sum_{l=1}^{m}\partial_{{\boldsymbol{k}}}g^{c}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{\gamma_{1},{\boldsymbol{\gamma}}^{l}}\mathbbm{1}^{k_{l}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}\big{)}(y)$
$\displaystyle=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}(\gamma_{1},\ldots,\gamma_{n+1})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{\boldsymbol{k}}g^{c}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}(y).$

Now, for $n=|{\boldsymbol{\gamma}}|\geq 1$ we have

$\displaystyle((Fg)_{*}\overline{Y})^{c}_{{\boldsymbol{\gamma}}}$
$\displaystyle=\\!\\!\\!\\!\sum_{({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in\overline{\mathrm{Sh}}^{-1}(\gamma_{1},\ldots,\gamma_{n-1})}\\!\\!\\!\\!\partial_{\boldsymbol{k}}(F_{\gamma_{n}}g^{c})F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{m}}\mathbbm{1}^{k_{m}}(Y)$
$\displaystyle=F_{(\gamma_{1},\ldots,\gamma_{n-1})}(F_{\gamma_{n}}g^{c})(y)$
$\displaystyle=F_{\boldsymbol{\gamma}}g^{c}(y)$

where we have used Section 2 and 79.

We now show the if implication of the first statement. Namely, we need to show
that if $\overline{Y}\in\mathscr{D}_{\boldsymbol{X}}(W)$ has the property that
$\overline{Y}_{n}=(F_{*}\overline{Y})_{n-1}$ for $n=1,\ldots,{\lfloor
p\rfloor}$ then $\overline{Y}_{n}=F^{n}\mathbbm{1}(Y)$. We show this by
induction on $n$. For $n=1$ the assertion is obvious. For the inductive step
we have

$\displaystyle\overline{Y}^{h}_{\boldsymbol{\gamma}}$
$\displaystyle=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{\boldsymbol{k}}F^{h}(Y)\overline{Y}^{k_{1}}_{{\boldsymbol{\gamma}}^{1}}\cdots\overline{Y}^{k_{1}}_{{\boldsymbol{\gamma}}^{1}}$
$\displaystyle=\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\begin{subarray}{c}({\boldsymbol{\gamma}}^{1},\ldots,{\boldsymbol{\gamma}}^{m})\in
A_{1}\cap\overline{\mathrm{Sh}}^{-1}({\boldsymbol{\gamma}})\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{\boldsymbol{k}}F^{h}F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}\cdots
F_{{\boldsymbol{\gamma}}^{1}}\mathbbm{1}^{k_{1}}(Y)$
$\displaystyle=F_{\boldsymbol{\gamma}}\mathbbm{1}^{h}(Y)$

where we have used 79 and the inductive hypothesis. ∎

###### Proof of Section 2.

Let $Y$ be a Davie solution to the RDE. Taking $g$ in 76 to be
$\mathbbm{1},F,\ldots,F^{{\lfloor p\rfloor}-2}\mathbbm{1}$ proves that
$\overline{Y}$ defined in 77 is indeed an element of
$\mathscr{D}_{\boldsymbol{X}}(W)$. By Section 2 we then have
$F_{*}\overline{Y}=(F(Y),\ldots,F^{{\lfloor p\rfloor}}\mathbbm{1}(Y))$ and by
Section 2

$\int_{s}^{t}F_{*}\overline{Y}\mathrm{d}\boldsymbol{X}\approx\langle(F_{*}\overline{Y})_{s},\boldsymbol{X}^{\geq
1}_{st}\rangle\approx Y_{st}$

again by the Davie definition. Since both the left and right hand sides are
increments of paths, we conclude by [Lyo98, Theorem 3.3.1] that identity must
hold. Therefore, since $Y_{0}=y_{0}$, 78 holds at the trace level, and for
$n\geq 1$

$\bigg{(}y_{0}+\overline{\int}F_{*}\overline{Y}\mathrm{d}\boldsymbol{X}\bigg{)}_{n}=(F_{*}\overline{Y})_{n-1}=\overline{Y}_{n}.$

Conversely, assume that there exists some
$\overline{Y}\in\mathscr{D}_{\boldsymbol{X}}(W)$ s.t. 78 holds: this implies
that for $m\geq 0$

$(F_{*}\overline{Y})_{m}=\bigg{(}y_{0}+\overline{\int}F_{*}\overline{Y}\mathrm{d}\boldsymbol{X}\bigg{)}_{m+1}=\overline{Y}_{m+1}$

and therefore by Section 2 $\overline{Y}$ must have the form 77. Finally,
$Y_{0}=y_{0}$ and for $g\in C^{\infty}(W)$ and
$\boldsymbol{Y}\coloneqq\uparrow_{\\!\boldsymbol{X}}\\!\\!\overline{Y}$

$\displaystyle g(Y)_{st}$
$\displaystyle=\int_{s}^{t}\overline{Dg}(Y)\mathrm{d}\boldsymbol{Y}$
$\displaystyle=\int_{s}^{t}(\overline{Dg}(Y)*\overline{Y})\cdot(\overline{F}(Y)*\overline{Y})\mathrm{d}\boldsymbol{X}$
$\displaystyle=\int_{s}^{t}(\overline{Dg}(Y)\cdot\overline{F}(Y))*\overline{Y}\mathrm{d}\boldsymbol{X}$
$\displaystyle=\int_{s}^{t}\overline{(Dg(\cdot)F(\cdot))}(Y)*\overline{Y}\mathrm{d}\boldsymbol{X}$
$\displaystyle=\int_{s}^{t}(Fg)_{*}\overline{Y}\mathrm{d}\boldsymbol{X}$
$\displaystyle\approx\sum_{n=1}^{\lfloor p\rfloor}\langle
F^{n}g(Y_{s}),\boldsymbol{X}^{n}_{st}\rangle$

where we have used Section 2, Section 2 and Section 2. This concludes the
proof. ∎

The above proposition tells us that once we have the solution in the sense of
Section 2 we can obtain an $\boldsymbol{X}$-controlled path, and thus by
Section 2 a rough path. If we want to emphasise the existence of these
superstructures we will write

$\mathrm{d}\overline{Y}=F(Y)\mathrm{d}\boldsymbol{X}\quad\text{and}\quad\mathrm{d}\boldsymbol{Y}=F(Y)\mathrm{d}\boldsymbol{X},\quad
Y_{0}=y_{0}$ (80)

i.e. $\boldsymbol{Y}\coloneqq\uparrow_{\\!\boldsymbol{X}}\\!\\!\overline{Y}$.
Notice that the initial condition only involves the trace.

The next result will be instrumental in defining RDEs on manifolds in a
coordinate-invariant manner.

######  Theorem  (Change of variable formula for RDE solutions).

Let $\boldsymbol{X},F,\boldsymbol{Y}$ be as above, $g\in C^{\infty}(W,U)$.
Then $(\boldsymbol{Y},g_{*}\boldsymbol{Y})$ jointly solve the RDE

$\mathrm{d}\begin{pmatrix}\boldsymbol{Y}\\\
\boldsymbol{Z}\end{pmatrix}=\begin{pmatrix}F(Y)\\\
Dg(Y)F(Y)\end{pmatrix}\mathrm{d}\boldsymbol{X}.$ (81)

In particular, if $g$ is invertible, Defining
$C^{\infty}(U,\mathcal{L}(V,U))\ni F_{g}(z)\coloneqq
Dg(g^{-1}(Z))F(g^{-1}(Z))$, $g_{*}\boldsymbol{Y}$ coincides with the rough
path solution to

$\mathrm{d}\boldsymbol{Z}=F_{g}(Z)\mathrm{d}\boldsymbol{X}.$ (82)

###### Proof.

Using Section 2, Section 2 and Section 2 we have

$\displaystyle\mathrm{d}(g_{*}\boldsymbol{Y})$
$\displaystyle=y_{0}+\boldsymbol{\int}\overline{Dg}(Y)\mathrm{d}\boldsymbol{Y}$
$\displaystyle=y_{0}+\boldsymbol{\int}Dg_{*}\overline{Y}\cdot
F_{*}\overline{Y}\mathrm{d}\boldsymbol{X}$
$\displaystyle=y_{0}+\boldsymbol{\int}(Dg(\cdot)F(\cdot))_{*}\overline{Y}\mathrm{d}\boldsymbol{X}.$

This proves the first claim; as for the second, we continue

$\displaystyle\mathrm{d}(g_{*}\boldsymbol{Y})$
$\displaystyle=\boldsymbol{\int}(Dg(\cdot)F(\cdot))_{*}\overline{Y}\mathrm{d}\boldsymbol{X}$
$\displaystyle=\boldsymbol{\int}(Dg(g^{-1}(\cdot))F(g^{-1}(\cdot)))_{*}g_{*}\overline{Y}\mathrm{d}\boldsymbol{X}$

where we have again used Section 2. This concludes the proof. ∎

The following theorem is proved in [CDL15, Corollary 2.17, Theorem 4.2] in the
case of $2\leq p<3$ and the proof carries over to the general case:

######  Theorem  (Local existence and uniqueness).

Precisely one of the following two possibility holds w.r.t. 75

  1. 1.

A solution on $[0,T]$ exists;

  2. 2.

There exists an $S\leq T$ and a solution on $[0,S)$, with $Y_{[0,S)}$ not
contained in any compact set of $\mathbb{R}^{e}$.

Moreover, in either case, the solution is unique on the interval on which it
is defined.

######  Example  (Non-autonomous RDEs).

We can define RDEs that also depend on the driving signal, by doubling the
variables i.e.

$\mathrm{d}\boldsymbol{Y}=F(Y,X)\mathrm{d}\boldsymbol{X}\
\stackrel{{\scriptstyle\text{def}}}{{\Longleftrightarrow}}\
\mathrm{d}\begin{pmatrix}\boldsymbol{X}\\\
\boldsymbol{Y}\end{pmatrix}=\begin{pmatrix}\mathbbm{1}\\\
F(Y,X)\end{pmatrix}\mathrm{d}\boldsymbol{X}.$ (83)

This will be important when defining RDEs on manifolds, driven by a manifold-
valued rough path.

##  3 Weakly geometric rough paths on manifolds

In this section we show how our algebraic framework for weakly geometric rough
paths can be deployed to transfer the theory to the manifold setting. In this
section $M$ will be a smooth $m$-dimensional manifold and $TM$ its tangent
bundle. The following definition is similar in spirit to the the one provided
in [BL15].

######  Definition  (Manifold-valued rough path).

Given a smooth atlas $\\{\varphi\colon
A_{\varphi}\to\mathbb{R}^{m}\\}_{\varphi}$ of $M$, an $M$-valued $p$-_weakly
geometric rough path_ controlled by $\omega$ on $[0,T]$,
$\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],M)$, consists of a collection
of rough paths
${{}^{\varphi}\\!}\boldsymbol{X}|_{[a_{\varphi},b_{\varphi}]}\in\mathscr{C}^{p}_{\omega}([a_{\varphi},b_{\varphi}],\mathbb{R}^{m})$
for all possible collections of intervals $[a_{\varphi},b_{\varphi}]$ (indexed
by the chart $\varphi$ in the atlas) s.t.
$\text{Im}({{}^{\varphi}\\!}X|_{[a_{\varphi},b_{\varphi}]})\subseteq\text{Range}(\varphi)$,
and that for each pair of intervals $[a_{\varphi},b_{\varphi}]$,
$[a_{\psi},b_{\psi}]$

$(\psi\circ\varphi^{-1})_{*}{{}^{\varphi}\\!}\boldsymbol{X}={{}^{\psi}\\!}\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([a_{\varphi},b_{\varphi}]\cap[a_{\psi},b_{\psi}],\mathbb{R}^{m}).$
(84)

The _trace_ of $\boldsymbol{X}$ is the path $t\mapsto
X_{t}\coloneqq\varphi^{-1}({{}^{\varphi}\\!}X_{t})\in M$ whenever
$t\in[a_{\varphi},b_{\varphi}]$ (independently of $\varphi$),
$X\in\mathcal{C}^{p}([0,T],M)$.

It makes sense to allow the mappings
$\varphi\mapsto{{}^{\varphi}\\!}\boldsymbol{X}$ and
$\varphi\mapsto[a_{\varphi},b_{\varphi}]$ to be multi-valued, so that the same
chart can be used multiple times (e.g. if the trace $X$ goes back and forth
between charts). The definition only depends on the smooth structure of $M$
(i.e. an equivalence class of atlases), since a rough path according to a
particular atlas is uniquely extended to the maximal atlas for the smooth
structure by reading 84 as a definition of the right hand side. Also note that
the definition is already fixed once we have
${{}^{\varphi}\\!}\boldsymbol{X}|_{[a_{\varphi},b_{\varphi}]}$ on a set of
intervals, one for each chart, s.t.
$\bigcup_{\varphi}(a_{\varphi},b_{\varphi})=(0,T)$: if the atlas if finite we
call this a _finite representation_ of $\boldsymbol{X}$. The following example
shows how our theory applies to the case of Stratonovich calculus of manifold-
valued semimartingales.

######  Example  (Stratonovich rough path).

Let $X$ be an $M$-valued continuous semimartingale, i.e. $f(X)$ is a real-
valued semimartingale for all $f\in C^{\infty}(M)$. Since semimartingales are
a.s. bounded $p$-variation for any $p>2$ we only need to define a rough path
above $X$ up to level $2$: we then define its _Stratonovich rough path_ in
coordinates by

$\boldsymbol{X}^{\alpha\beta}_{st}\coloneqq\int_{s}^{t}X^{\alpha}_{su}\circ\mathrm{d}X^{\beta}_{u}$
(85)

where the integral is intended in the Stratonovich sense. This is well-known
to a.s. define a rough path in the linear setting, but since the above
coordinate expression is taken accoding to a chart, we need the following
lemma to establish that it defines a (stochastic) rough path in the sense of
Section 3.

######  Lemma .

Let $X$ be an $\mathbb{R}^{m}$-valued continuous semimartingale, $f\in
C^{\infty}(\mathbb{R}^{m},\mathbb{R}^{n})$ and $\boldsymbol{X}$ be defined as
in 85. Then

$(f_{*}\boldsymbol{X})^{ij}_{st}=\int_{s}^{t}f^{i}(X)\circ\mathrm{d}f^{j}(X).$
(86)

###### Proof.

Using that, in the linear setting, Stratonovich integrals a.s. coincide with
rough integrals against the Stratonovich rough path [FH14, Theorem 9.1]

$\displaystyle\int_{s}^{t}f^{i}(X)\circ\mathrm{d}f^{j}(X)$
$\displaystyle=\int_{s}^{t}f^{i}(X)\partial_{\gamma}f^{j}(X)\circ\mathrm{d}X$
$\displaystyle=\int_{s}^{t}\overline{f}^{i}\overline{\partial}_{\gamma}\overline{f}^{j}(X)\mathrm{d}\boldsymbol{X}^{\gamma}$
$\displaystyle\approx
f^{i}\partial_{\gamma}f^{j}(X_{s})X_{st}^{\gamma}+(\partial_{\alpha}f^{i}\partial_{\beta}f^{j}+f^{i}\partial_{\alpha\beta}f^{j})(X_{s})\boldsymbol{X}_{st}^{\alpha\beta}$
$\displaystyle\approx
f^{i}(X_{s})f^{j}(X)_{st}+\partial_{\alpha}f^{i}\partial_{\beta}f^{j}(X_{s})\boldsymbol{X}^{\alpha\beta}_{st}.$

We then have

$\int_{s}^{t}f^{i}(X)_{su}\mathrm{d}f^{j}(X_{u})=\int_{s}^{t}f^{i}(X)\mathrm{d}f^{j}(X)-f^{i}(X_{s})f^{j}(X)_{st}\approx\partial_{\alpha}f^{i}\partial_{\beta}f^{j}(X_{s})\boldsymbol{X}^{\alpha\beta}_{st}$

concluding the proof. ∎

We can define a path valued in a fixed vector space $V$ and controlled by a
manifold-valued rough path as follows:

######  Definition .

Let $\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],M)$. We define an
$V$-valued _$\boldsymbol{X}$ -controlled path_
$\overline{H}\in\mathscr{D}_{\boldsymbol{X}}(V)$ to be a collection
${{}^{\varphi}\\!}\overline{H}\in\mathscr{D}_{{{}^{\varphi}\\!}\boldsymbol{X}}(V)$
with $\varphi,a_{\varphi},b_{\varphi}$ as in Section 3 and

${{}^{\varphi}\\!}\overline{H}*\overline{(\varphi\circ\psi^{-1})}({{}^{\psi}\\!}X)={{}^{\psi}\\!}\overline{H}.$

As for rough path, this only depends on the smooth structure of $M$. If we
wish to define a controlled integrand, the trace of $H$ must lie in the bundle
$\mathcal{L}(TM,W)$ for some vector space $W$, above $X$: for this we need a
separate definition:

######  Definition  (Controlled integrand).

Let $\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],M)$. We define an
$W$-valued _$\boldsymbol{X}$ -controlled integrand_
$\overline{H}\in\mathscr{D}_{\boldsymbol{X}}(\mathcal{L}(TM,W))$ to be a
collection
${{}^{\varphi}\\!}\overline{H}\in\mathscr{D}_{{{}^{\varphi}\\!}\boldsymbol{X}}(\mathcal{L}(\mathbb{R}^{m},W))$
with $\varphi,a_{\varphi},b_{\varphi}$ as in Section 3 and

$(\varphi\circ\psi^{-1})^{*}{{}^{\varphi}\\!}\overline{H}={{}^{\psi}\\!}\overline{H}.$
(87)

The _trace_ of $\overline{H}$ is the path $H\coloneqq{{}^{\varphi}\\!}H\circ
T_{X}\varphi$, which is valued in the fibre of $X$ of the bundle
$\mathcal{L}(TM,W)$ (if $W=\mathbb{R}^{e}$ then this bundle is $(T^{*}M)^{e}$.

We can now define rough integration on manifolds.

######  Definition  (Rough integral on manifolds).

Let $\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],M)$ and
$\overline{H}\in\mathscr{D}_{\boldsymbol{X}}(\mathcal{L}(TM,W))$. We define
the _rough integral_

$\int_{0}^{\cdot}\overline{H}\mathrm{d}\boldsymbol{X}\coloneqq\sum_{[s_{\varphi},t_{\varphi}]}\int_{s_{\varphi}}^{t_{\varphi}}{{}^{\varphi}\\!}\overline{H}\mathrm{d}{{}^{\varphi}\\!}\boldsymbol{X}\in\mathcal{C}^{p}_{\omega}([0,T],W)$
(88)

where we are summing over a finite partition of $[0,\cdot]$ whose intervals
$[s_{\varphi},t_{\varphi}]$ are indexed by charts $\varphi$ with the property
that each $X_{[s_{\varphi},t_{\varphi}]}$ is contained in the domain of
$\varphi$.

This definition does not depend on the subdivision (since the integral taken
w.r.t. two different subdivisions coincides with that take w.r.t. to their
common refinement) and does not depend on the charts used, since, for another
chart $\psi$ and $X_{[s,t]}$ contained in the domains of $\varphi$ and $\psi$

$\displaystyle\int_{s}^{t}{{}^{\psi}\\!}\overline{H}\mathrm{d}{{}^{\psi}\\!}\boldsymbol{X}$
$\displaystyle=\int_{s}^{t}(\varphi\circ\psi^{-1})^{*}{{}^{\varphi}\\!}\overline{H}\mathrm{d}(\psi\circ\varphi^{-1})_{*}{{}^{\varphi}\\!}\boldsymbol{X}$
$\displaystyle=\int_{s}^{t}(\psi\circ\varphi^{-1})^{*}(\varphi\circ\psi^{-1})^{*}{{}^{\varphi}\\!}\overline{H}\mathrm{d}{{}^{\varphi}\\!}\boldsymbol{X}$
$\displaystyle=\int_{s}^{t}{{}^{\varphi}\\!}\overline{H}\mathrm{d}{{}^{\varphi}\\!}\boldsymbol{X}$

by Section 2 and Section 2.

We proceed to the topic of RDEs on manifolds: following the approach of [É89],
and for maximum generality, we will consider both the driving rough path and
the solution to be manifold-valued; for this purpose we let $N$ be an
$n$-dimensional manifold.

######  Definition  (RDEs on manifolds).

Let $\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],M)$ and $F$ a section of
the vector bundle $\mathcal{L}(TM,TN)$ (i.e. the vector bundle over $N\times
M$ with fibres $\mathcal{L}(TM,TN)_{y,x}\coloneqq\mathcal{L}(T_{x}M,T_{y}N)$).
We will say that $\boldsymbol{Y}\in\mathscr{C}^{p}_{\omega}([0,T],M)$ is a
solution to the RDE

$\mathrm{d}\boldsymbol{Y}=F(Y,X)\mathrm{d}\boldsymbol{X},\quad Y_{0}=y_{0}$
(89)

if for all charts $\varphi$ and $\psi$ on $M$ and $N$ respectively we have
$Y_{0}=y_{0}$ and on all intervals $[s,t]$ s.t. $X_{[s,t]}$ is contained in
the domain of $\varphi$ and $Y_{[s,t]}$ is contained in the domain of $\psi$

$\mathrm{d}{{}^{\psi}\\!}\boldsymbol{Y}=F({{}^{\psi}\\!}Y,{{}^{\varphi}\\!}X)\mathrm{d}({{}^{\varphi}\\!}\boldsymbol{X})\quad\text{on
}[s,t]$

in the sense of Section 2, with initial condition ${{}^{\psi}\\!}Y_{s}$.

This does not depend on the chart: given other chats $\overline{\varphi}$,
$\overline{\psi}$ we have

$\displaystyle\mathrm{d}\begin{pmatrix}{{}^{\overline{\varphi}}\\!}\boldsymbol{X}\\\
{{}^{\overline{\psi}}\\!}\boldsymbol{Y}\end{pmatrix}=\mathrm{d}\begin{pmatrix}(\overline{\varphi}\circ\varphi^{-1})_{*}{{}^{\varphi}\\!}\boldsymbol{X}\\\
(\overline{\psi}\circ\psi^{-1})_{*}{{}^{\psi}\\!}\boldsymbol{Y}\end{pmatrix}=\begin{pmatrix}\mathbbm{1}\\\
F({{}^{\overline{\psi}}\\!}Y,{{}^{\overline{\varphi}}\\!}X)\end{pmatrix}\mathrm{d}({{}^{\overline{\varphi}}\\!}\boldsymbol{X})$

by Section 2 applied to the change of variable
$(\overline{\varphi}\circ\varphi^{-1},\overline{\psi}\circ\psi^{-1})$. Local
existence can be inferred from Section 2, which in particular implies global
existence if $M$ is compact. More general conditions that guarantee global
existence can be found in [Wei18, Dri18].

We end with a few brief remarks which link the topic of this paper to the
existing literature, without elaborating on the details.

######  Remark  (Stratonovich calculus).

We can further expand on Section 3 to include the following: if $H$ is a
$\mathcal{L}(TM,W)$-valued semimartingale above the $M$-valued semimartingale
$X$ (i.e. $H_{t}\in\mathcal{L}(T_{X_{t}}M,W)$) the Stratonovich integral $\int
H\mathrm{d}X$ coincides a.s. with the rough integral
$\int\overline{H}\mathrm{d}\boldsymbol{X}$ against the Stratonovich rough path
$\boldsymbol{X}$ for any choice of $\overline{H}$ with trace $H$, and the
solution to the Stratonovich differential equation
$\mathrm{d}Y=F(Y,X)\mathrm{d}X$ coincides a.s. with the RDE 89.

######  Remark  (The extrinsic viewpoint).

In [CDL15] the topic of manifold-valued theory of rough paths, rough
integration (specifically of 1-forms) and RDEs was treated from the extrinsic
point of view. Here
$\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],\mathbb{R}^{d})$ is defined
in [CDL15, Definition 3.17] to be _constrained_ to a smoothly embedded
manifold $M$ if its trace is $M$-valued and for all 1-forms
$F\in\Gamma\mathcal{L}(T\mathbb{R}^{d},W)$ ($\Gamma$ denoting the space of
sections)

$\forall x\in M\ F(x)|_{T_{x}M}=0\ \Rightarrow\
\boldsymbol{\int}\overline{F}(X)\mathrm{d}\boldsymbol{X}=0\in\mathscr{C}^{p}_{\omega}([0,T],W).$
(90)

In [CDL15, Corollary 3.32, Proposition 3.35] this is shown to be equivalent to
the trace $X$ being $M$-valued and $(I\otimes
Q(X_{s}))\boldsymbol{X}^{2}_{st}\approx 0$, or equivalently to
$(P(X_{s})\otimes
P(X_{s}))\boldsymbol{X}^{2}_{st}\approx\boldsymbol{X}^{2}_{st}$, where for
$x\in M$ $P(x)$ is the orthogonal projection
$T_{x}\mathbb{R}^{d}\twoheadrightarrow T_{x}M$ and $Q\coloneqq\mathbbm{1}-P$.
Moreover, one may replace $\boldsymbol{X}^{2}_{st}$ with its antisymmetric
part $(\wedge\boldsymbol{X}^{2})_{st}$ in these identities (because
$\odot\boldsymbol{X}^{2}$ is already fixed by the trace).

This approach carries over to the case of higher $p$ considered here. If
$\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],\mathbb{R}^{d})$ we may say
that it is _constrained_ to the smoothly embedded manifold $M$ if
$\pi_{*}\boldsymbol{X}=\boldsymbol{X}$, where $\pi$ is the Riemannian
projection of a tubular neighbourhood $U$ of $M$ onto $M$ (i.e. it maps a
point in $U$ to the unique point on $M$ closest to it - this is well-defined
and smooth on a thin enough tubular neighbourhood). This extends the
definition of [CDL15] since $D_{x}\pi=P(x)$. In order to generalise the
equivalent condition $(I\otimes Q(X_{s}))(\wedge\boldsymbol{X})_{st}$ we can
take the log of our original condition, i.e.
$\log\pi_{*}\boldsymbol{X}=\log\boldsymbol{X}$: this has the advantage of
eliminating all the redundancies of the former (as explained in [LS06,
p.767]), and its precise coordinate expression can be derived by using [FV10,
Definition 7.20], but at higher orders cannot be described in terms of
antisymmetric tensors. The Chen-Strichartz formula [Bau04, Theorem 1.1],
however, expresses $\log\boldsymbol{X}$ as a Lie polynomial; the task of
expressing it in a basis of the Lie algebra is more complex still [Rei17].

Care must be taken when defining the rough integral of a controlled integrand,
since it is no longer the case that for an $\boldsymbol{X}$-controlled
integrand $\overline{H}$, $\int\overline{H}\mathrm{d}\boldsymbol{X}$ (defined
in the ordinary sense, where $\boldsymbol{X}$ is considered an element of
$\mathscr{C}^{p}_{\omega}([0,T],\mathbb{R}^{d})$) does not always only depend
on the trace $H$ of $\boldsymbol{X}$ restricted to $TM$, although this is
indeed the case when $\overline{H}$ is given by a 1-form: this is because if
$F\in\Gamma\mathcal{L}(T\mathbb{R}^{d},W)$ ($\Gamma$ denoting the space of
sections) vanishes on $TM$ (i.e. $F(x)P(x)=0$ for $x\in M$) by Section 2 we
have

$\begin{split}&\int\overline{F}(X)\mathrm{d}\boldsymbol{X}=\int\overline{F}(X)\mathrm{d}\pi_{*}\boldsymbol{X}=\int\pi^{*}\overline{F}(X)\mathrm{d}\boldsymbol{X}=0\\\
\text{since}\quad&\pi^{*}\overline{F}(X)=(\overline{F}(X)*\overline{\pi}(X))\cdot\overline{D\pi}(X)=\overline{F}(X)\cdot\overline{P}(X)=\overline{F(\cdot)P(\cdot)}(X)=0\end{split}$
(91)

where we have used Section 2. This then implies that if
$F,G\in\Gamma\mathcal{L}(T\mathbb{R}^{d},W)$ restrict to the same element of
$\Gamma\mathcal{L}(TM,W)$ then
$\int\overline{F}(X)\mathrm{d}\boldsymbol{X}=\int\overline{G}(X)\mathrm{d}\boldsymbol{X}$.
Similarly, we have that if
$\overline{H},\overline{K}\in\mathcal{D}_{\boldsymbol{X}}(\mathcal{L}(TM,W))$
are such that $\pi^{*}\overline{H}=\pi^{*}\overline{K}$ then
$\int\overline{H}\mathrm{d}\boldsymbol{X}=\int\overline{K}\mathrm{d}\boldsymbol{X}$,
but this involves conditions on all levels of $\overline{H}$, not just the
trace (for a simple counterexample where this identity fails when only
assuming $HP(X)=KP(X)$ see [ABCRF20, Example 4.3]).

Finally, the original definition of constrained rough path given by
integration 90 also carries over to higher $p$. The fact that this is implied
by $\pi_{*}\boldsymbol{X}=\boldsymbol{X}$ was shown in 91. For the converse,
we rewrite the identity as $(\mathbbm{1}-\pi)_{*}\boldsymbol{X}=0$: the trace
level is implied by the fact that $X$ is $M$-valued, and at orders $\geq 1$
the identity
$\boldsymbol{\int}\overline{Q}(X)\mathrm{d}\boldsymbol{X}=(\mathbbm{1}-\pi)_{*}\boldsymbol{X}$
is straightforward to check.

######  Remark  (Parallel transport and Cartan development).

If $M$ has a connection $\nabla$ (with Christoffel symbols $\Gamma$), given
$\boldsymbol{X}\in\mathscr{C}^{p}_{\omega}([0,T],M)$ we can define parallel
transport of vectors above its trace $X$ as the solution to the RDE driven by
$\boldsymbol{X}$, valued in $TM$, defined by the horizontal lift based at $X$:
in local coordinates the parallel frame path $t\mapsto A_{t}\in T_{X_{t}}M$
satisfies

$\mathrm{d}A^{\gamma}=-\Gamma^{\gamma}_{\alpha\beta}(X)A^{\beta}\mathrm{d}\boldsymbol{X}^{\alpha}.$
(92)

The Cartan development $\boldsymbol{Y}$ of a
$\boldsymbol{Z}\in\mathscr{C}^{p}_{\omega}([0,T],T_{o}M)$ (with $o\in M$ a
fixed basepoint) can be viewed as the projection onto $M$ of an $FM$-valued
RDE driven by $\boldsymbol{Z}$, defined the fundamental horizontal vector
fields (the full solution additionally consists of a parallel frame above the
developed path): in local coordinates

$\begin{cases}\mathrm{d}\boldsymbol{Y}^{k}=A^{k}_{\gamma}\mathrm{d}\boldsymbol{Z}^{\gamma}\\\
\mathrm{d}A^{k}_{\gamma}=-\Gamma^{k}_{ij}(Y)A^{i}_{\alpha}A^{j}_{\gamma}\mathrm{d}\boldsymbol{Z}^{\alpha}\end{cases}.$
(93)

Here we have emphasised how development actually defines a (possibly
explosive) $M$-valued rough path $\boldsymbol{Y}$, and not just its trace.

## References

  * [ABCRF20] John Armstrong, Damiano Brigo, Thomas Cass, and Emilio Rossi Ferrucci.  Non-geometric rough paths on manifolds.  arXiv:2007.06970, 2020.  https://arxiv.org/abs/2007.06970.
  * [Bau04] Fabrice Baudoin.  An Introduction to the Geometry of Stochastic Flows.  Published by Imperial College Press and distributed by World Scientific Publishing Co., 2004.
  * [BL15] Youness Boutaib and Terry Lyons.  A new definition of rough paths on manifolds.  arXiv:1510.07833v2, 2015.  https://arxiv.org/abs/1510.07833v2.
  * [CDL15] Thomas Cass, Bruce K. Driver, and Christian Litterer.  Constrained rough paths.  Proc. Lond. Math. Soc. (3), 111(6):1471–1518, 2015.
  * [CDLL16] Thomas Cass, Bruce K. Driver, Nengli Lim, and Christian Litterer.  On the integration of weakly geometric rough paths.  J. Math. Soc. Japan, 68(4):1505–1524, 10 2016.
  * [Dri18] Bruce Driver.  Global existence of geometric rough flows.  arXiv:1810.03708, 2018.  https://arxiv.org/abs/1810.03708.
  * [É89] Michel Émery.  Stochastic calculus in manifolds.  Universitext. Springer-Verlag, Berlin, 1989.  With an appendix by P.-A. Meyer.
  * [EFP15] Kurusch Ebrahimi-Fard and Frédéric Patras.  Cumulants, free cumulants and half-shuffles.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2176):20140843, 2015.
  * [FH14] Peter K. Friz and Martin Hairer.  A course on rough paths.  Universitext. Springer, Cham, 2014.  With an introduction to regularity structures.
  * [Foi13] Löıc Foissy.  An introduction to hopf algebras of trees.  2013\.
  * [FV10] Peter K. Friz and Nicolas B. Victoir.  Multidimensional stochastic processes as rough paths, volume 120 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, 2010.  Theory and applications.
  * [Gub04] Massimiliano Gubinelli.  Controlling rough paths.  Journal of Functional Analysis, 216(1):86 – 140, 2004.
  * [Gub10] Massimiliano Gubinelli.  Ramification of rough paths.  J. Differential Equations, 248(4):693–721, 2010.
  * [HK15] Martin Hairer and David Kelly.  Geometric versus non-geometric rough paths.  Ann. Inst. H. Poincaré Probab. Statist., 51(1):207–251, 02 2015\.
  * [LCL07] Terry Lyons, Michael Caruana, and Thierry Lévy.  Differential equations driven by rough paths ecole d’eté de probabilités de saint-flour xxxiv-2004, 2007.
  * [LS06] Terry J. Lyons and Nadia Sidorova.  On the radius of convergence of the logarithmic signature.  Illinois J. Math., 50(1-4):763–790, 2006.
  * [Lyo98] Terry Lyons.  Differential equations driven by rough signals.  Revista Matemática Iberoamericana, 14(2):215–310, 1998.
  * [Lyo02] Terry Lyons.  System control and rough paths.  Oxford mathematical monographs. Clarendon Press, Oxford, 2002.
  * [Man06] Dominique Manchon.  Hopf algebras, from basics to applications to renormalization.  arXiv:math/0408405v2, 2006.  https://arxiv.org/abs/math/0408405.
  * [Rei17] Jeremy Reizenstein.  Calculation of iterated-integral signatures and log signatures.  arXiv:1712.02757, 2017.  https://arxiv.org/abs/1712.02757.
  * [Reu93] Christophe. Reutenauer.  Free lie algebras.  London Mathematical Society monographs. New series ; 7. Clarendon, Oxford, 1993.
  * [Wei18] Martin Weidner.  A geometric view on rough differential equations.  PhD thesis, Imperial College London, 8 2018.  https://spiral.imperial.ac.uk:8443/handle/10044/1/62658.


# Precision-based attacks and interval refining: how to break, then fix,
differential privacy on finite computers

Samuel Haney
Tumult Labs
<EMAIL_ADDRESS>Damien Desfontaines
Tumult Labs
<EMAIL_ADDRESS>Luke Hartman
Tumult Labs
<EMAIL_ADDRESS>Ruchit Shrestha
Tumult Labs
<EMAIL_ADDRESS>Michael Hay
Tumult Labs
Colgate University
<EMAIL_ADDRESS>
###### Abstract

Despite being raised as a problem over ten years ago, the imprecision of
floating point arithmetic continues to cause privacy failures in the
implementations of differentially private noise mechanisms. In this paper, we
highlight a new class of vulnerabilities, which we call _precision-based
attacks_ , and which affect several open source libraries. To address this
vulnerability and implement differentially private mechanisms on floating-
point space in a safe way, we propose a novel technique, called _interval
refining_. This technique has minimal error, provable privacy, and broad
applicability. We use interval refining to design and implement a variant of
the Laplace mechanism that is equivalent to sampling from the Laplace
distribution and rounding to a float. We report on the performance of this
approach, and discuss how interval refining can be used to implement other
mechanisms safely, including the Gaussian mechanism and the exponential
mechanism.

##  1 Introduction

There are many issues that can arise when translating abstract differentially
private algorithms to real-world implementations. One such issue, as observed
by Mironov [24] over ten years ago, is the use of floating point arithmetic.
Mironov demonstrated this issue using the Laplace distribution, commonly used
as a building block for DP mechanisms: he showed that sampling Laplace noise
in a naive way creates “holes” in the output distribution, and that these
holes can be leveraged to leak private data and nullify the privacy guarantee.
Later, similar issues were found in other common primitives, like the
exponential mechanism [19] and the Gaussian mechanism [20]. The researchers
exposing these vulnerabilities also proposed mitigations, but privacy
libraries tend not to use these mitigations: some are complex, not
generalizable, and in some cases come with a large utility cost. Instead, the
creators of libraries have invented other techniques to mitigate floating
point attacks.

  * •

Google’s differential privacy team developed alternative procedures to add
Laplace and Gaussian noise by sampling _discrete_ distributions, and rounding
to a fixed granularity [17]. These are implemented in Google’s differential
privacy libraries [6].

  * •

Holohan et al. [18] proposed sampling Laplace and Gaussian noise in a way that
makes it computationally harder for an attacker to reverse-engineer the
output, without any rounding. This approach is implemented in diffprivlib [3],
IBM’s differential privacy library.

  * •

The authors of SmartNoise Core [11] and OpenDP [9] implemented yet another
approach, which consists of using MPFR [5], an arbitrary precision arithmetic
library, to generate a _hole-free_ unit Laplace or Gaussian distribution [10].

Unfortunately, we show that the approaches taken by diffprivlib, SmartNoise
Core and OpenDP are susceptible to a new class of floating point
vulnerabilities, which we call _precision-based attacks_ , and describe in
Section 2.

Google’s approach has other drawbacks: (1) the technique cannot be easily
generalized to a safe implementation of other mechanisms, like the exponential
mechanism, (2) the rounding precision has to be set in advance, which can be a
usability issue and limits the mechanism safety to a strict subset of the
floating-point range, and (3) it involves a small $\delta$ in the privacy
guarantee, which makes it difficult to use with privacy accounting based on
pure or zero-concentrated DP.

In this paper, we propose a new technique for safely implementing DP
mechanisms in floating-point space. Our technique, which we call _interval
refining_ , has three key properties.

Minimal error.

It has the same error as if we were sampling the desired distribution exactly,
and rounding to the nearest float. The only error comes from this last
rounding step, which is unavoidable for any mechanism returning a float.

Provable privacy.

It provides the same privacy guarantees as the abstract algorithm it
implements. As such, it can be used in contexts where the privacy accounting
is based on pure DP, zero-concentrated DP, or other definitions forbidding
infinite privacy loss.

Broad applicability.

It can be used to implement any mechanism in floating-point space. In this
paper, we describe our implementation of the Laplace mechanism, but we also
outline how the technique can be applied to other mechanisms, such as the
Gaussian and exponential mechanisms.

Our technique is conceptually simple, and follows the same idea as inverse
transform sampling, where a sample from a uniform distribution is transformed
into a sample of an arbitrary probability distribution by passing it as input
to the inverse cumulative distribution function (CDF) of this distribution.
But instead of simply generating a sample (which may be an irrational number),
we sample an _interval_ , initially large and then iteratively refined, until
the entire interval would be rounded to the same floating point number, and we
can return this float. Using this technique, we implement the Laplace
mechanism (Algorithm 1) in a way that is identical to sampling from the real-
valued Laplace distribution and rounding the sample to the next highest float.
This process is illustrated in Fig. 1.

The rest of the paper is structured as follows: In Section 2, we describe
novel floating point attacks with applications to existing differential
privacy libraries. In Section 3, we outline the key ideas underlying our
technique, and explain how to simulate sampling from a real-valued
distribution and rounding the result to floating-point. In Section 4, we
instantiate our technique and give a precise algorithm implementing sampling
from the Laplace distribution, and rounding to the next highest float; we
describe our implementation in the Tumult platform [13] and report on
empirical results. In Section 5, we discuss extensions, limitations, and
future work.

##  2 Precision-based attacks

In this section, we outline the vulnerabilities we found in existing
libraries. These vulnerabilities, which we call _precision-based attacks_ ,
rely on two simple observations. First, if a double-precision floating-point
number $x$ is such that $2^{k}\leq|x|<2^{k+1}$ for some integer $k$, then $x$
is a multiple of $2^{k-52}$. This quantity $2^{k-52}$ is known as the _unit in
the last place_ [16] (or _ulp_) of $x$, we denote it by $\mathrm{ulp}_{x}$.
Second, when adding any double to a fixed double $x$, then the output is
always a multiple of $\mathrm{ulp}_{x}/2$. The proof of these facts can be
found in Appendix A.

This creates a simple vulnerability affecting implementations that add noise
to floating-point numbers without any rounding: when adding noise to e.g. $1$,
all possible outputs are multiples of $2^{-53}$. But when adding noise to 0,
the output is exactly the value of the sampled noise. If it is possible to
sample noise that is _not_ a multiple of $2^{-53}$, then this creates a
distinguishing event between inputs $0$ and $1$. This is exactly what happens
in additive noise mechanisms in diffprivlib, SmartNoise Core, and OpenDP: if
the noise value $r$ is such that $|r|<0.5=2^{-1}$, then it might not be a
multiple of $2^{-53}$. If we return it as is (after adding it to 0), the
attacker can deduce that the true value was not $1$. This can happen
arbitrarily often as the noise scale gets smaller; with Laplace noise of scale
1 (corresponding to a counting query with $\varepsilon=1$), approximately 25%
of outputs are distinguishing events in diffprivlib, SmartNoise Core, and
OpenDP.

Perhaps surprisingly, precision-based attacks can create vulnerabilities in
other algorithms, besides simple additive noise mechanisms. Consider the
mechanism to compute quantiles based on the exponential mechanism, introduced
in [25]. This algorithm works in three steps: it splits the output space in
intervals based on the input data, privately chooses one interval using the
exponential mechanism, and returns a uniform number from this interval. To
implement the last step and generate a uniformly random number in an interval
$[x,y)$, a naive approach consists of generating a uniformly random number $r$
in $[0,1)$, and returning $x+(y-x)\cdot r$. The possible precision of the
output depends on the value of $x$, which creates the opportunity for
precision-based attacks. Finding a pair of input databases demonstrating such
distinguishing events is a little more involved than for additive noise
mechanisms. In Appendix A, we provide more detail on this class of
vulnerabilities, and show pairs of datasets which lead to distinguishing
events for quantile mechanisms in both diffprivlib and SmartNoise Core.

This novel class of attacks shows that mitigating floating-point
vulnerabilities is more difficult than it seems. It suggests that approaches
implemented in production-grade software should be formally documented, and
provide a proof that they satisfy the desired privacy guarantees.

These vulnerabilities were communicated to, and acknowledged by, the authors
of diffprivlib, SmartNoise Core, and OpenDP via personal correspondence and a
public bug report [4] in November and December 2021.

##  3 Overview of interval refining

In this section, we give an overview of our interval refining technique. Our
goal is to simulate the following process: sample $X\sim\mathcal{D}$, where
$\mathcal{D}$ is a distribution over sample space $\Omega=\mathbb{R}$, and
then round the value to a 64 bit floating point number. Any rounding scheme is
acceptable; for the purposes of exposition, we round up. If we call
$\mathsf{float}:\Omega\rightarrow S$ a function that maps any real number to
the next largest 64-bit floating point number, our goal is to sample
$X\sim\mathcal{D}$ and output $\mathsf{float}(X)$.

The main complication of this approach is that prior to rounding, the sampled
value may be an irrational number and thus cannot be represented with finite
memory. To solve this issue, we need three key ideas.

The first key idea starts with the observation that the float function
partitions the sample space into _intervals_ , where all elements in each
interval are assigned to the same float value: $\mathsf{float}$ maps elements
of the (infinite) sample space to a finite space $S$ (the set of 64-bit
floats). Therefore, our goal is to sample an element $s\in S$ with probability
equal to the probability of sampling an element from $\mathcal{D}$ that maps
to $s$. That is, $\Pr[\text{sampling
}s]=\Pr_{X\sim\mathcal{D}}[X\in\mathsf{float}^{-1}(s)]$. This means that
instead of sampling $X$, we can sample a progressively finer interval around
$X$: we start from a large interval, and iteratively refine it until all
values within the interval map to the same float. Conceptually, we never
sample $X$ directly; instead, we sample enough information about $X$ to
determine to which float it should be mapped.

This gives us the following high-level process.

  1. 1.

Set the current interval $I$ to be entire sample space:
$I=\Omega=[-\infty,\infty]$

  2. 2.

Partition the current interval $I$ into disjoint intervals $I_{1},I_{2},\dots$
such that $\bigcup_{i}I_{i}=I$, and sample interval $I_{i}$ with probability
proportional to $Pr_{X\sim\mathcal{D}}[X\in I_{i}|X\in I]$.

  3. 3.

Set the current interval $I$ to the selected interval $I_{i}$.

  4. 4.

If all values in $I$ round to the same float $s$ (i.e., $\exists
s:I\subseteq\mathsf{float}^{-1}(s)$), then return $s$.

  5. 5.

Otherwise, repeat from Step 2.

For step 2, there are a number of reasonable ways to partition the space. For
reasons that will be explained, we choose to partition $I$ into two equi-
probable sub-intervals $I_{1}$, $I_{2}$ such that $Pr_{X\sim\mathcal{D}}[X\in
I_{1}]=Pr_{X\sim\mathcal{D}}[X\in I_{2}]$. One can show that the process
outlined above samples $s$ with the correct probability.

One challenge remains: the interval boundaries might be irrational numbers,
which cannot be exactly represented on a computer with finite memory.

This brings us to our second key idea: using _inverse transform sampling_ [8].
The intuition is simple: sampling $U\sim\mathrm{Uniform}(0,1)$ and computing
$X=F^{-1}_{\mathcal{D}}(U)$, where $F^{-1}_{\mathcal{D}}$ is the inverse
cumulative distribution function for distribution $\mathcal{D}$, is the same
as sampling $X\sim\mathcal{D}$. Applying this observation to our setting,
instead of sampling intervals in the sample space of the distribution of
interest, we can instead sample intervals uniformly in $[0,1]$, the sample
space of $\mathrm{Uniform}(0,1)$, and use the inverse CDF to map these
intervals to the sample space of the distribution we want to sample from. A
nice feature of this is that the probability of a uniform random variable
being in an interval is equal to the interval’s width, so if we want to sample
equiprobable intervals, we simply divide the current interval in half and pick
a half at random.

The process now looks like this:

  1. 1.

Set the current interval to $[0,1]$.

  2. 2.

Divide the current interval in half, and select one of the halves uniformly at
random.

  3. 3.

Set the current interval to be the selected half, and denote the endpoints of
the current interval as $(a,b)$.

  4. 4.

If all values in interval $[F_{\mathcal{D}}^{-1}(a),F_{\mathcal{D}}^{-1}(b)]$
round to the same float $s$, then return $s$.

  5. 5.

Otherwise, repeat from Step 2.

Note that in this updated process, the interval end points $a,b$ are dyadic
rationals and therefore can be represented exactly on a computer with finite
memory.

There is one more complication to tackle: irrational numbers might still occur
in the termination step (4). This brings us to our third and final key idea:
approximating the inverse CDF to find an interval $I$ that contains
$[F_{\mathcal{D}}^{-1}(a),F_{\mathcal{D}}^{-1}(b)]$, and whose endpoints are
rational numbers. Then, the termination step becomes:

  1. 4.

Compute an interval $I$ with rational endpoints such that
$[F_{\mathcal{D}}^{-1}(a),F_{\mathcal{D}}^{-1}(b)]\subseteq I$. If all values
in interval $I$ round to the same float $s$, then return $s$.

This approximate interval may be too wide and thus the algorithm may fail to
terminate at the appropriate iteration, or fail to terminate entirely (e.g. if
the approximation $I$ is always the entire real number line). Our method of
approximating $I$, discussed in Section 4, has a parameter that controls the
precision of the approximation: by increasing the precision in each iteration,
we can show that the algorithm eventually terminates with probability 1. The
fact that the approximation error may cause the algorithm to run for
additional iterations is not a problem: once the algorithm has arrived at
endpoints $(a,b)$ such that
$[F_{\mathcal{D}}^{-1}(a),F_{\mathcal{D}}^{-1}(b)]\subseteq\mathsf{float}^{-1}(s)$
for some $s$, further iterations will not change the outcome, since any
subinterval of $[a,b]$ will still map to the same outcome $s$.

##  4 Laplace mechanism: algorithm and implementation

In this section, we describe our algorithm for safely sampling from a Laplace
distribution, explain how it can be the basis for the Laplace mechanism, and
report on our experience implementing it.

Algorithm 1 is our proposed technique for sampling a Laplace random variable
rounded to a float. Let $\mathcal{D}_{\mu,\beta}$ denote the Laplace
distribution with a given location $\mu$ and scale $\beta$. This algorithm
uses a subroutine, IntervalInvCDF, for computing the inverse CDF of
${\mathcal{D}_{\mu,\beta}}$ on an interval $[a,b]$. The subroutine takes as
input a tuple $\langle a,b\rangle$ of interval endpoints and a parameter
$prec$ that controls the working precision of the function, which influences
how closely the returned interval approximates the true interval.

A little more formally, IntervalInvCDF satisfies the following properties:

$F_{{\mathcal{D}_{\mu,\beta}}}^{-1}(x)\in\textsf{{IntervalInvCDF}}_{\mu,\beta}(\langle
a,b\rangle,prec)\text{ for }x\in[a,b]$ (1)

and

$\lim_{prec\to\infty}\textsf{{IntervalInvCDF}}_{\mu,\beta}(\langle
a,b\rangle,prec)=[F_{{\mathcal{D}_{\mu,\beta}}}^{-1}(a),F_{{\mathcal{D}_{\mu,\beta}}}^{-1}(b)]$
(2)

Algorithm 1 starts with the $[0,1]$ interval. On each iteration, it picks one
half of this interval uniformly at random and computes the inverse CDF (10) on
it to obtain an interval $[s,t]$. The precision is increased in each iteration
to ensure that even samples near a boundary between two floats can eventually
be distinguished. The algorithm terminates when $s$ and $t$ both round up to
the same floating point number. The check for termination happens on 11 using
the NextFloat function, which takes as input an arbitrary-precision floating
point number and rounds it to the next highest 64-bit number. Fig. 1
visualizes a run of this algorithm.

Algorithm 1 SampleLaplace$(\mu,\beta)$

1:$a,b\leftarrow\langle 0,1\rangle$

2:$prec\leftarrow 0$

3:while True do

4: if $\textsf{{RandBit}}()=0$ then

5: $a\leftarrow\frac{a+b}{2}$

6: else

7: $b\leftarrow\frac{a+b}{2}$

8: end if

9: $prec\leftarrow prec+1$

10: $\langle s,t\rangle\leftarrow\textsf{{IntervalInvCDF}}_{\mu,\beta}(\langle
a,b\rangle,prec)$ $\triangleright$ For any $r\in[a,b]$,
$F_{{\mathcal{D}_{\mu,\beta}}}^{-1}(r)\in[s,t]$

11: if $\textsf{{NextFloat}}(s)=\textsf{{NextFloat}}(t)$ then

12: return $\textsf{{NextFloat}}(s)$

13: end if

14:end while

We now formally state the claim that Algorithm 1 is arbitrarily close to
sampling from a Laplace and rounding to the next float. Let
$\mathcal{D}_{\text{{{SampleLaplace}}}({\mu,\beta,k})}$ be the probability
distribution over the outputs of Algorithm Algorithm 1 on inputs $\mu$,
$\beta$ after running for at most $k$ iterations; we define the output to be
$\bot$ if the algorithm hasn’t terminated after $k$ rounds. Let
$\mathcal{D}^{\mathsf{float}}_{\mu,\beta}$ be the probability distribution of
$\mathsf{float}(X)$ where $X\sim\mathcal{D}_{\mu,\beta}$.

######  Theorem 1.

For any $\mu$, $\beta>0$ where $\mu$ and $\beta$ are arbitrary-precision
floats,

$\lim_{k\rightarrow\infty}\mathsf{TVD}(\mathcal{D}_{\text{{{SampleLaplace}}}({\mu,\beta,k})},\mathcal{D}^{\mathsf{float}}_{\mu,\beta})=0$

where $\mathsf{TVD}$ is the total variation distance.

Proof of this theorem appears in Appendix B.

This algorithm can be the basis for a variant of the Laplace mechanism [15],
which noisily evaluates some target function $f$ on the private data $x$. To
do so, it is essential to invoke SampleLaplace with $\mu=f(x)$ rather than
invoking it with $\mu=0$ and adding the result to $f(x)$. The precision-based
attacks of Section 2 illustrate one of the potential vulnerabilities with the
latter approach.

Figure 1: A visualization of Algorithm 1. The green, red, blue and orange
lines correspond to the four iterations before the algorithm terminates, and
the ticks along the x-axis represent floating point numbers (not drawn to
scale). The colored lines below the x-axis represent intervals output by
IntervalInvCDF on the corresponding input intervals $(a,b)$ shown on the
y-axis. Note that each interval contains
$\set{F^{-1}_{{\mathcal{D}_{\mu,\beta}}}(x)\mid a\leq x\leq b}$, shown by the
sub-interval within the dotted lines but may be slightly wider due to the
approximation error. It is this approximation that prevents the algorithm from
terminating in the third iteration (corresponding to the blue intervals). By
the fourth iteration (orange intervals), the output interval lies completely
between two consecutive floating point numbers which suffices for termination
on 11.

We have implemented Algorithm 1 in the Tumult platform [13]. Our
implementation is in Python. All steps in the algorithm can be computed
_exactly_ with the exception of IntervalInvCDF, which is approximate in that
it returns an interval that may be _wider_ than the “true” interval (i.e.,
$[F_{\mathcal{D}}^{-1}(a),F_{\mathcal{D}}^{-1}(b)]$). Nevertheless, for the
correctness of the algorithm, the returned interval must be guaranteed to
contain the true interval.

To ensure this guarantee, the implementation relies on Arb [22], a C library
for arbitrary-precision ball arithmetic. Ball arithmetic enables computing
with real numbers by explicitly and automatically tracking error bounds
throughout the computation. With ball arithmetic, a function
$f:\mathcal{X}\rightarrow\mathcal{Y}$ that ordinarily takes a number and
returns a number is implemented as a function $F$ that takes in an interval
(represented by a midpoint and radius, a “ball”) and returns an interval with
the property that the returned interval contains the correct answer for any
number in the input interval (i.e., $\forall x\in[x_{1},x_{2}]:f(x)\in
F([x_{1},x_{2}])$). This correctness guarantee composes nicely and allows us
to construct a function like IntervalInvCDF from the base functions in Arb and
be confident the interval it returns contains the true interval.

Arb allows one to increase the working precision which generally speaking
reduces the size of the resulting interval, though this is not a universal
guarantee and depends on the computation. For the inverse CDF of the Laplace
distribution, which only involves basic arithmetic and a log function, we
believe, based on the Arb documentation, that the error decreases
exponentially with working precision (cf. Sec. 2.3.6 of [21]).

A difference between the implementation in the Tumult platform and Algorithm 1
is that instead of sampling a single bit per iteration, we sample 63 bits,
thereby choosing among $2^{63}$ equiprobable intervals in each iteration
rather than just two. This improves the runtime of the algorithm by a factor
of $\approx$ 40x.

We simulated 20 million samples. The throughput is 30,000 samples per second.
Almost all samples terminate in a single iteration, but 2% ($\approx$400,000)
took two iterations, and we never witnessed a sample requiring 3 or more
iterations.

##  5 Extensions, limitations, and future work

While in this paper we describe a technique for safely sampling from a target
distribution and rounding to a float, it should be cautioned that this cannot
necessarily be used as a drop in replacement in more complex algorithms that
rely on sampling from specific distributions as a subroutine. For some of
these algorithms (e.g., the Sparse Vector Technique [15], PrivTree [26]), the
proof of correctness leverages properties of the target distribution that may
not hold for its rounded variant.

Nevertheless, we believe that our technique, interval refining, can be
usefully extended to some of these more complex algorithms. To illustrate this
point, we briefly describe how it can be used for sampling a noisy argmax, a
subroutine that can be used inside the exponential mechanism [23] (via the
Gumbel-max trick [14]) or as a mechanism itself (e.g., the Report Noisy Max
function [15]). The noisy argmax problem is as follows: given $n$ private
values $x_{1},\dots,x_{n}$, report $\operatorname*{arg\,max}_{i}x_{i}+Z_{i}$
where $Z_{i}$ is noise sampled independently from some target distribution
(e.g., Gumbel). Note that we cannot simply sample _rounded_ noise values as
this would not be equivalent. But we can still use the interval refining
technique: we maintain intervals around the noisy value of each element and we
terminate once we have found the largest noisy value (i.e., there is one
interval that is strictly larger than all others). The interval refining
techniques gives us flexibility to chooose different termination conditions to
suit the particular application.

We are in the process of implementing the interval refining technique for the
exponential mechanism [23] and an algorithm for noisy quantiles [25]. A formal
description of the algorithm and a proof of its correctness are left as future
work.

A limitation of the interval refinement technique is that the runtime can be
difficult to analyze as the number of iterations required depends on the shape
of the inverse CDF and the approximation technique used to calculate it. In
practice, with our implementation of the Laplace mechanism, we never witnessed
a sample fail to reach the termination condition, with 98% terminating after a
single iteration. A more careful characterization of the runtime is the
subject of future work.

## References

  * [1] Chorus.  https://github.com/uvm-plaid/chorus.  Accessed 2022-04-13.
  * [2] diffpriv: easy differential privacy in R.  https://github.com/brubinstein/diffpriv.  Accessed 2022-04-13.
  * [3] Diffprivlib: the IBM differential privacy library.  https://github.com/IBM/differential-privacy-library.  Accessed 2022-04-13.
  * [4] Floating-point issue in noise samplers.  https://github.com/opendp/opendp/issues/414.  Accessed 2022-04-13.
  * [5] The GNU MPFR library.  https://mpfr.org.  Accessed 2022-04-13.
  * [6] Google’s differential privacy libraries.  https://github.com/google/differential-privacy.  Accessed 2022-04-13.
  * [7] IEEE 754.  https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules.  Accessed 2022-04-13.
  * [8] Inverse transform sampling.  https://en.wikipedia.org/wiki/Inverse_transform_sampling.  Accessed 2022-05-02.
  * [9] OpenDP.  https://github.com/opendp/opendp.  Accessed 2022-04-13.
  * [10] samplers.rs – opendp.  https://github.com/opendp/opendp/blob/bf0456891805ea68c99ba30be226cf7067c3af22/rust/src/samplers.rs#L543.  Accessed 2022-04-13.
  * [11] SmartNoise core.  https://github.com/opendp/smartnoise-core.  Accessed 2022-04-13.
  * [12] SmartNoise SQL.  https://github.com/opendp/smartnoise-sdk/tree/main/sql.  Accessed 2022-04-13.
  * [13] The tumult platform.  https://tmlt.io/platform.  Accessed 2022-04-13.
  * [14] Ryan Adams.  The gumbel-max trick for discrete distributions.  https://lips.cs.princeton.edu/the-gumbel-max-trick-for-discrete-distributions/.  Accessed 2022-05-02.
  * [15] Cynthia Dwork and Aaron Roth.  The algorithmic foundations of differential privacy.  Found. Trends Theor. Comput. Sci., 9(3–4):211–407, aug 2014\.
  * [16] David Goldberg.  What every computer scientist should know about floating-point arithmetic.  ACM Comput. Surv., 23(1):5–48, mar 1991.
  * [17] Google Differential Privacy Team.  Secure noise generation.  https://github.com/google/differential-privacy/blob/main/common_docs/Secure_Noise_Generation.pdf.  Accessed 2022-04-13.
  * [18] Naoise Holohan and Stefano Braghin.  Secure Random Sampling in Differential Privacy.  In Elisa Bertino, Haya Shulman, and Michael Waidner, editors, Computer Security – ESORICS 2021, Lecture Notes in Computer Science, pages 523–542, Cham, 2021. Springer International Publishing.
  * [19] Christina Ilvento.  Implementing the exponential mechanism with base-2 differential privacy.  In Proceedings of the 2020 ACM SIGSAC Conference on Computer and Communications Security, pages 717–742, 2020.
  * [20] Jiankai Jin, Eleanor McMurtry, Benjamin IP Rubinstein, and Olga Ohrimenko.  Are we there yet? Timing and floating-point attacks on differential privacy systems.  arXiv preprint arXiv:2112.05307, 2021.
  * [21] F. Johansson.  Arb - a c library for arbitrary-precision ball arithmetic.  https://arblib.org/arb.pdf.  Accessed 2022-05-02.
  * [22] F. Johansson.  Arb: efficient arbitrary-precision midpoint-radius interval arithmetic.  IEEE Transactions on Computers, 66:1281–1292, 2017.
  * [23] Frank McSherry and Kunal Talwar.  Mechanism design via differential privacy.  In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS’07), pages 94–103. IEEE, 2007.
  * [24] Ilya Mironov.  On significance of the least significant bits for differential privacy.  In Proceedings of the 2012 ACM Conference on Computer and Communications Security, CCS ’12, pages 650–661, New York, NY, USA, October 2012. Association for Computing Machinery.
  * [25] Adam Smith.  Privacy-preserving statistical estimation with optimal convergence rates.  In Proceedings of the forty-third annual ACM symposium on Theory of computing, pages 813–822, 2011.
  * [26] Jun Zhang, Xiaokui Xiao, and Xing Xie.  Privtree: A differentially private algorithm for hierarchical decompositions.  In Proceedings of the 2016 International Conference on Management of Data, SIGMOD ’16, page 155–170, New York, NY, USA, 2016. Association for Computing Machinery.

##  Appendix A Further details on precision-based attacks

In this section, we provide additional details on the precision-based attacks
described in Section 2.

We consider double-precision floating-point numbers as specified by IEEE 754
[7], called _doubles_ for short. A double uses 64 bits: 1 sign bit $s$, 11
exponent bits $e$, and 52 mantissa bits $d_{1}\dots d_{52}$; the corresponding
floating-point number is:

${\left(-1\right)}^{s}\cdot{\left(1.d_{1}\dots d_{52}\right)}\cdot
2^{e-1023}.$

An exponent value of $e=0$ is used to represent 0 and subnormal numbers, while
an exponent value of $e=2047$ is used to represent infinities and NaN values;
we will ignore these edge cases. In IEEE 754, arithmetic operations between
doubles are _correctly rounded_ : they must be computed exactly, and rounded
to the closest double. (In case the result of an operation falls exactly
between two successive doubles, special rounding rules apply.) In the rest of
this paper, we denote by $\oplus$, $\ominus$, and $\otimes$ the floating-point
analogues of addition, substraction, and multiplication.

The _unit in the last place_ , or _ulp_ , of a double $x$ is the size of the
interval between $x$ and the next consecutive double. We denote it as
$\mathrm{ulp}_{x}$. For example, the ulp of $x=1$ (represented by $s=0$,
$e=1023$, and $d_{i}=0$ for all $i$) is $2^{-52}$, while the ulp of
$x=1-2^{-53}$ (represented by $s=0$, $e=1022$, and $d_{i}$=1 for all $i$) is
$2^{-53}$. The first fact in Section 2 follows directly from the definition of
doubles: if a double $x$ is such that $2^{k}\leq x<2^{k+1}$, then
$\mathrm{ulp}_{x}=2^{k-52}$.

Let us now formalize and prove the second observation.

######  Theorem 2.

Let $x$ and $y$ be two doubles, with $x\neq 0$. Then $x\oplus y$ is a multiple
of $\mathrm{ulp}_{x}/2$.

###### Proof.

Let us assume, without loss of generality, that $x>0$. Let $k$ be such that of
$2^{k}\leq x<2^{k+1}$; note that $\mathrm{ulp}_{x}=2^{k-52}$. There are two
cases, based on whether $y<-2^{k-1}$.

  * •

If $y<-2^{k-1}$, then $\mathrm{ulp}_{y}$ is at least
$2^{k-1-52}=\mathrm{ulp}_{x}/2$, so $y$ is a multiple of $\mathrm{ulp}_{x}/2$.
Since $x$ is a multiple of $\mathrm{ulp}_{x}$, it’s also a multiple of
$\mathrm{ulp}_{x}/2$, and $x\oplus y$ is a multiple of $\mathrm{ulp}_{x}/2$.

  * •

If $y\geq-2^{k-1}$, then $x\oplus y\geq-2^{k-1}+2^{k}=2^{k-1}$. Then $x\oplus
y$ has a ulp of at least $2^{k-1-52}=\mathrm{ulp}_{x}/2$, so regardless of the
value of $y$, the sum will be a multiple of $\mathrm{ulp}_{x}/2$.

∎

Precision-based attacks on additive noise mechanisms follow directly from this
fact.

### Additive noise mechanisms

Distinguishing events can be found simply by adding noise to $0$ or $1$, which
can be two outputs of a counting query evaluated on two neighboring datasets.
All noisy values obtained from $1$ are multiples of $2^{-53}$, while a large
number of noisy values obtained from $0$ are not. This affects systems which
do not attempt to mitigate floating-point vulnerabilities, like Chorus [1] or
diffpriv [2]. More interestingly, the vulnerability also affects diffprivlib
(Gaussian, Analytic Gaussian, and Staircase mechanisms, as well as all
variants of the Laplace mechanism except the snapping mechanism), SmartNoise
Core (Laplace, Gaussian, and Truncated Gaussian mechanisms), and OpenDP
(Laplace, Gaussian, and Analytic Gaussian mechanisms), even though these
libraries attempt to mitigate floating-point attacks.

Note that in SmartNoise Core and OpenDP, users need to explicitly opt in to
using floating-point primitives using a flag, and are warned that doing so is
potentially risky. Therefore, it could be argued that this problem is
technically not a vulnerability. However, at the time of writing of this
paper, some tools relying on these libraries set this flag by default, and do
not warn users who are using these primitives; this is the case for e.g.
SmartNoise SQL [12].

### Quantiles mechanisms

It is a little less straightforward to find distinguishing events for quantile
mechanisms based on the exponential mechanism. Recall that the naive mechanism
first splits the output space in intervals based on the data, then choosing an
interval $[x,y)$ using the exponential mechanism, sampling a uniform number
$r$ in $[0,1)$, and returning $x\oplus(y\ominus x)\otimes r$.

In diffprivlib, $r$ is generated using the SystemRandom.random() function from
Python’s standard library, which generates multiples of $2^{-53}$. This means
that with $D_{1}=[0,0,1]$, all outputs will be multiples of $2^{-53}$. With
$D_{2}=[0,0.25,1]$ however, whenever the interval $[0,0.25)$ is selected by
the exponential mechanism, and the output of SystemRandom.random() is not a
multiple of $2^{-51}$, then the returned number will not be a multiple of
$2^{-53}$, creating distinguishing events.

SmartNoise Core attempts to prevent vulnerabilities by using MPFR [5] to
generate $r$ in a hole-free way: all possible doubles in $[0,1)$ can be
generated. As a consequence, the previous choice of $D_{1}$ and $D_{2}$ does
not create distinguishing events. However, if we use $D_{1}=[-1,1,1]$, then
the only possible sampled interval is $[-1,1)$: because of the addition with
$1$, all outputs will be multiples of $2^{-53}$. Using $D_{2}=[-1,0,1]$, the
interval $[0,1)$ might be sampled, in which case the output might not be a
multiple of $2^{-53}$.

Interestingly, SmartNoise Core does _not_ block the use of this quantiles
mechanism on the user-specified floating-point option which is necessary to
access floating-point additive noise mechanisms. This underscores that
floating-point vulnerabilities can occur in places where they would not be
expected, and thereby evade scrutiny.

##  Appendix B Proof of Theorem 1

We begin by recalling the statement of the theorem:

See 1

To simplify notation we let
$P^{(k)}=\mathcal{D}_{\text{{{SampleLaplace}}}({\mu,\beta,k})}$ and
$Q=\mathcal{D}^{\mathsf{float}}_{\mu,\beta}$. Recall that $S$ is the set of
rational numbers representable as floats and that $P^{(k)}$ is a distribution
over $S\cup\\{\bot\\}$111While $\bot$ is not a possible outcome for
distribution $Q$, we can consider it as a possible outcome with probability 0,
and thus the TVD is well-defined., where $s\in S$ is the event that the
algorithm terminates in the first $k$ rounds and outputs $s$, and $\bot$ is
the event that the algorithm does not terminate after $k$ rounds. $Q$ is the
distribution over $S$ that results from picking a real number from
$\mathcal{D}_{\mu,\beta}$ and rounding it to the nearest float.

The proof the theorem has the following structure:

  1. 1.

We show that for all $k\in[1,\infty)$ and for all $s\in S$, $P^{(k)}(s)\leq
Q(s)$ (Claim 3).

  2. 2.

We show that $\lim_{k\rightarrow\infty}P^{(k)}(\bot)=0$ (Claim 4). That is,
the probability that the algorithm does not terminate within $k$ rounds goes
to zero as $k$ goes to infinity.

  3. 3.

We show that (1) and (2) together imply that the theorem holds.

Before we can prove the first of these claims, we need to show some properties
about the intervals produced in the intermediate stages of the algorithm. We
call the interval that the algorithm sets to $[a_{k},b_{k}]$ at the beginning
of each iteration the _initial interval_ for that iteration. We let
$P_{init}^{(k)}$ denote the probability distribution over choices of initial
interval that result in the algorithm terminating in that round. That is,
$P^{(k)}_{init}(I)$ is the probability that the algorithm picks $I$ as its
initial interval in some round, and then terminates at the end of the round.
Note that whether the algorithm terminates on a given initial interval is
deterministic, so the algorithm always terminates in the first round that it
chooses an initial interval in
$\operatorname{supp}\left(P^{(k)}_{init}\right)$. We let
$P^{(k)}_{init}(\bot)$ be the probability that the algorithm does not
terminate in the first $k$ rounds. See Figure 2 for an illustration of the
possible initial intervals chosen by the Algorithm.

Figure 2: An illustration of the possible initial intervals chosen by the
4-round version of the algorithm. In each round, the algorithm picks one of
the two sub intervals of the interval it chose in the previous round. The
solid green intervals represent intervals that, if chosen, cause the algorithm
to terminate in the given round. Dotted gray intervals represent intervals
that could be chosen as the initial interval in a given round, but will never
cause the algorithm to terminate. Red intervals with x’s represent intervals
that can never be chosen by the algorithm because they are sub intervals of
the green intervals. Finally, blue intervals marked with squares represent
intervals in the last round that may be chosen and cause the algorithm to
output $\bot$. The distribution $P^{(4)}_{init}$ is a distribution over
intervals $a$ through $e$ and $\bot$, with probabilities
$1/4,1/4,1/8,1/8,1/16$ and $3/16$ respectively. Note that the green intervals
are disjoint (see Claim 1 for a proof).

Next, we want to consider the distributions that result when this initial
interval is transformed by the inverse CDF function. We let $F^{-1}$ denote
the interval version of the inverse CDF function for distribution
$\mathcal{D}_{\mu,\beta}$. That is $F^{-1}(I)$ is inverse CDF applied to each
endpoint of interval $I$. Additionally, we let $\tilde{F}^{-1}_{k}$ denote the
imprecise interval inverse CDF function with precision $k$,
$\textsf{{IntervalInvCDF}}_{\mu,\beta}(\cdot,k)$. We consider two
distributions: the one resulting from applying $F^{-1}$ to $P^{(k)}_{init}$,
denoted $P^{(k)}_{F^{-1}}$, and the one resulting from applying
$\tilde{F}^{-1}_{k}$ to $P^{(k)}_{init}$, denoted $P^{(k)}_{\tilde{F}^{-1}}$
(we drop the subscript $k$ on $\tilde{F}^{-1}$ since it is already present in
the superscript). The probability of $\bot$ in each distribution is the same,
i.e.
$P^{(k)}_{init}(\bot)=P^{(k)}_{F^{-1}}(\bot)=P^{(k)}_{\tilde{F}^{-1}}(\bot)$.

The following two claims (Claim 1 and 2) help to prove Claim 3. We first show
that intervals in the support of $P^{(k)}_{init}$ are disjoint, and so are the
intervals in the support of $P^{(k)}_{F^{-1}}$ (however, the same is not
necessarily true of $P^{(k)}_{\tilde{F}^{-1}}$).

######  Claim 1.

All intervals in $\operatorname{supp}(P^{(k)}_{init})$ are pairwise disjoint.
Additionally, all intervals in $\operatorname{supp}(P^{(k)}_{F^{-1}})$ are
pairwise disjoint.

###### Proof.

Suppose a pair of intervals
$I,I^{\prime}\in\operatorname{supp}\left(P^{(k)}_{int}\right)$ is overlapping.
Note that each interval only has the possibility of being chosen by the
algorithm on some fixed round (since the length of the intervals the algorithm
considers halves every round). We let $round(I)$ denote this round for
interval $I$.

Next, we claim that if the algorithm chooses interval $I$ as its initial
interval in round $round(I)$, then it must be the case that the algorithm
terminates in that round. Since the interval is in the support of the
distribution, the algorithm must terminate with some positive probability
after choosing $I$. However, whether the algorithm terminates after choosing
an interval in a given round is deterministic, and therefore the algorithm
either always terminates for a given interval, or never terminates.

Finally, note that any two intervals chosen by the algorithm can only
intersect if one is strictly contained within the other. Assume WLOG that
$I^{\prime}\subseteq I$, and therefore $round(I)<round(I^{\prime})$. On any
run of the algorithm that terminates in round $round(I^{\prime})$ after
choosing $I^{\prime}$, it must be the case that on round $round(I)$, the
algorithm chose $I$ as its initial interval. Therefore, the algorithm should
have terminated in $round(I)$, giving a contradiction.

The same property holds for $P^{(k)}_{F^{-1}}$ since the inverse CDF function
is monotonically increasing, and therefore a pair of intervals in
$\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$ that the corresponding pair
of intervals in $\operatorname{supp}\left(P^{(k)}_{init}\right)$ would
overlap. ∎

######  Claim 2.

For all $I\in\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$, $I\neq\bot$,

$P^{(k)}_{F^{-1}}(I)=\mathcal{D}_{\mu,\beta}(I).$ (3)

###### Proof.

First, we claim that for all
$J\in\operatorname{supp}\left(P^{(k)}_{init}\right)$, $P^{(k)}_{init}(J)=|J|$.
Note that for the algorithm to choose interval $J$, it must have chosen the
only interval in each of the previous rounds that contains $J$, i.e. there is
only one sequence of choices that leads to $J$ being chosen. Let $round(J)$ be
defined as in the proof of Claim 1. Let $y=round(J)$, and let
$J_{1},\ldots,J_{y-1}$ be the sequence of initial intervals the algorithm
chose in the first $y-1$ rounds. Then, the algorithm will never terminate
after choosing any $J_{i}$ (by Claim 1), and the probability that the
algorithm chooses $J_{i}$ given that it chose $J_{i-1}$ is $1/2$. The claim
follows from induction over $J_{i}$.

Next, fix $I\neq\bot\in\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$. Let
$J=F(I)$ be the corresponding initial interval that the algorithm chooses.
Then,

$\displaystyle P^{(k)}_{F^{-1}}(I)$ $\displaystyle=P^{(k)}_{init}(J)$
$\displaystyle=|J|$ $\displaystyle=F(I[1])-F(I[0])$ (4)
$\displaystyle=\mathcal{D}_{\mu,\beta}(I).$

The notation $F$ in equation 4 denotes the non-interval CDF function. ∎

######  Claim 3.

For all $k\in[1,\infty)$ and for all $s\in S$, $P^{(k)}(s)\leq Q(s)$.

###### Proof.

Let $\mathcal{I}_{s}$ be the subset of
$\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$ that results in the
algorithm producing output $s$. That is, if we let $float^{-1}(s)$ be the
preimage of float $s$ under function $float$, then
$\mathcal{I}_{s}=\left\\{I\in\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)\text{
s.t. }\tilde{F}^{-1}(F(I))\subseteq float^{-1}(s)\right\\}$. Then,

$\displaystyle P^{(k)}(s)$ $\displaystyle=P^{(k)}_{F^{-1}}(\mathcal{I}_{s})$
(by definition of $\mathcal{I}_{s}$)
$\displaystyle=\sum_{I\in\mathcal{I}_{s}}P^{(k)}_{F^{-1}}(I)$
$\displaystyle=\sum_{I\in\mathcal{I}_{s}}\mathcal{D}_{\mu,\beta}(I)$ (by Claim
2)
$\displaystyle=\mathcal{D}_{\mu,\beta}\left(\bigcup_{I\in\mathcal{I}_{s}}I\right)$
(by Claim 1) $\displaystyle\leq Q(s).$ (by definition of $\mathcal{I}_{s}$)

∎

######  Claim 4.

$\lim_{k\rightarrow\infty}P^{(k)}(\bot)=0$. That is, the probability that the
algorithm does not terminate within $k$ rounds goes to zero as $k$ goes to
infinity.

###### Proof.

We bound the probability that the algorithm does not terminate in $k$ rounds,
given that the algorithm reached round $k$. This is an upper bound on the
probability that the algorithm does not terminate in $k$ rounds,
$P^{(k)}(\\{\bot\\})$.

Let $I$ denote the algorithm’s choice of an initial interval for round $k$.
The algorithm does not terminate if and only if $s\in\tilde{F}^{-1}_{k}(I)$
for some $s\in S$. We consider two cases:

Case 1: $s\in F^{-1}(I)$ for some $s\in S$. There are only $|S|$ intervals in
$\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$ for which this is the case
(by Claim 1), and for all such intervals $J$, $P^{(k)}_{F^{-1}}(J)=2^{-k}$.
Therefore, letting $\mathcal{I}_{1}$ denote this set of intervals,

$P^{(k)}_{F^{-1}}(\mathcal{I}_{1})\leq 2^{-k}|S|.$ (5)

Case 2: $s\not\in F^{-1}(I)$ for any $s\in S$ (but $\tilde{F}^{-1}_{k}$ does
contain some $s\in S$). Then, we use our assumption about the precision of
$\tilde{F}^{-1}_{k}$ (Equation 2). In particular, Equation 2 implies that
there is a function $\delta(k)$ such that

$\lim_{k\rightarrow\infty}\delta(k)=0,$

and for any interval $J$,

$d\left(J,\tilde{F}^{-1}_{k}(I)\right)\leq\delta(k),$ (6)

where $d(I,J)$ is defined as $|I[0]-J[0]|+|I[1]-J[1]|$. In particular, this
true of $J=F^{-1}(I)$. Then, the condition for case 2 is only possible if at
least one of the endpoints of $F^{-1}(I)$ is within distance $\delta(k)$ of
some $s\in S$. We split this into two subcases.

Case 2(a): Only one of the endpoints of $F^{-1}(I)$ is within distance
$\delta(k)$ of some $s\in S$ (the other endpoint may be far from all elements
of $S$, or it may be close to some other $t\in S,t\not=s$). By Claim 1, this
can be true of at most two intervals in
$\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$ per element of $S$. Denote
this set of intervals $\mathcal{I}_{2}$. There are at most $2|S|$ of these
intervals, and each has probability $2^{-k}$ of being chosen, so

$P^{(k)}_{F^{-1}}(\mathcal{I}_{2})\leq 2^{1-k}|S|.$ (7)

Case 2(b): There exists an $s\in S$ such that both endpoints of $F^{-1}(I)$
are within $\delta(k)$ of $s$. Let $\mathcal{I}_{3}$ denote the set of the
intervals of $\operatorname{supp}\left(P^{(k)}_{F^{-1}}\right)$ meeting this
condition. Then, since the intervals in $\mathcal{I}_{3}$ are disjoint, the
total length of the intervals in $\mathcal{I}_{3}$ is at most
$2\cdot\delta(k)\cdot|S|$. So,

$\displaystyle P^{(k)}_{F^{-1}}(\mathcal{I}_{3})$
$\displaystyle=\mathcal{D}_{\mu,\beta}\left(\bigcup_{J\in\mathcal{I}_{3}}J\right)$
$\displaystyle\leq\Omega_{\mathcal{D}_{\mu,\beta}}(2\cdot\delta(k)\cdot|S|),$

where $\Omega_{\mathcal{D}}(x)$ denotes the maximum probability that
distribution $\mathcal{D}$ assigns to a subset of reals of total length at
most $x$. For the Laplace distribution, $lim_{|x|\rightarrow
0}\Omega_{\mathcal{D}}(x)=0$.

Combining these cases we get

$\displaystyle P^{(k)}(\bot)$ $\displaystyle\leq
P^{(k)}_{F^{-1}}(\mathcal{I}_{1})+P^{(k)}_{F^{-1}}(\mathcal{I}_{2})+P^{(k)}_{F^{-1}}(\mathcal{I}_{3})$
$\displaystyle\leq
2^{-k}|S|+2^{1-k}|S|+\Omega_{\mathcal{D}_{\mu,\beta}}(2\cdot\delta(k)\cdot|S|)$
$\displaystyle\operatorname{\stackrel{{\scriptstyle\text{def}}}{{=}}}\alpha(k).$

Then the claim follows since $\lim_{k\rightarrow\infty}\alpha(k)=0$. ∎

Finally, we show how Claims 3 and 4 imply that
$\lim_{k\rightarrow\infty}TVD(P^{(k)},Q)=0$.

$\displaystyle\mathsf{TVD}(P^{(k)},Q)$ $\displaystyle=\frac{1}{2}\sum_{s\in
S\cup\\{\bot\\}}\left|P^{(k)}(s)-Q(s)\right|$
$\displaystyle=\frac{1}{2}\sum_{s\in
S}\left|P^{(k)}(s)-Q(s)\right|+\frac{1}{2}P^{(k)}(\bot)$ (since $\bot$ is not
in the support of Q) $\displaystyle=\frac{1}{2}\sum_{s\in
S}\left(Q(s)-P^{(k)}(s)\right)+\frac{1}{2}P^{(k)}(\bot)$ (by Claim 3)
$\displaystyle=\frac{1}{2}-\frac{1}{2}\sum_{s\in
S}P^{(k)}(s)+\frac{1}{2}P^{(k)}(\bot)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}P^{(k)}(S)+\frac{1}{2}P^{(k)}(\bot)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\left(1-P^{(k)}(\bot)\right)+\frac{1}{2}P^{(k)}(\bot)$
$\displaystyle=P^{(k)}(\bot)$

Therefore, by Claim 4,

$\lim_{k\rightarrow\infty}\mathsf{TVD}(P^{(k)},Q)=0$ (8)

which completes the proof.