text:
\begin{document}

\setlength{\parindent}{1.5em}
\setlength{\parskip}{1.2mm}
\renewcommand{\baselinestretch}{1.2}

Now we consider a generalization of the Z-shape graph matrix discussed in Section \ref{section z-graph matrix}.

\begin{defn}\label{defn:m-zshape}
\begin{figure}[h]
    \centering
    \includegraphics[scale=0.35]{m-layer-z-shape.png}
    \caption{The $m-$layer Z-shape $\mzshape$}
    \label{fig:m-zshape}
    \end{figure}

    Let $\mzshape$ be the bipartite shape with vertices $V(\mzshape)=\{u_1,\dots,u_m,v_1,\dots, v_m\}$ and edges $E(\mzshape)=\left\{\{u_i,v_i\}: i\in[m]\right\}\cup \left\{\{u_{i+1},v_{i}\}: i\in[m-1]\right\}$ with distinguished tuples of vertices $U_{\mzshape}=(u_1,\dots,u_m)$ and $V_{\mzshape}=(v_1,\dots, v_m)$. See Figure \ref{fig:m-zshape} for an illustration.

    We refer to $\mzshape$ as the {\bfseries\boldmath $m-$layer Z-shape}. Note that $\a_{Z(2)}$ is the Z-shape $\a_Z$ as in definition \ref{defn:z-shape}.
\end{defn}

\begin{rmk}\label{rmk num constr edges for m-Z}
    For the $m-$layer Z-shape $\mzshape$, the size of the minimum separator is $m$. By Lemma \ref{lem:min constr edges for bipartite shape}, for any nonzero-valued constraint graph $C\in\C_{\(\mzshape,2q\)}$, $\abs{E(C)}\geq m\cdot(q-1)$. By Corollary \ref{cor:dominantnumberofedges}, dominant constraint graphs $C\in\C_{\(\mzshape,2q\)}$ have $m\cdot(q-1)$ edges.
\end{rmk}

\begin{defn}\label{defn D(m,n)}
For $m,n$ positive integers,
    \begin{equation}
        D(m,n)=\frac{1}{m\cdot n+1}{(m+1)\cdot n\choose n}\,.
    \end{equation}
\end{defn}

\begin{rmk}
    The number $\displaystyle{D_n=\dfrac{1}{2n+1}{3n\choose n}}$ in Section \ref{section z-graph matrix} is $D(2,n)$ here. Also $D(m-1,n)=A_n(m,1)$ where the generalized Catalan number $\displaystyle{A_n(k,r)=\dfrac{r}{nk+r}{nk+r\choose n}}$ is defined in Remark \ref{rmk generalized Catalan num}.
\end{rmk}

Below is the main result of this section.
\begin{thm}\label{thm num of m-Z constr graph}
    Let $\mzshape$ be the $m-$layer Z-shape as in definition \ref{defn:m-zshape}. Then the number of dominant constraint graphs $C\in\C_{\(\mzshape,2q\)}$ is $D(m,q)$.
\end{thm}

\begin{rmk}
    When $m=2$, $D(m,n)=D(2,n)=D_n$, $\mzshape=\a_Z$, and this theorem is exactly Theorem \ref{thm num of Z constr graph}.
\end{rmk}

A direct corollary we get from the above theorem is the following.
\begin{cor}\label{cor M expected trace}
    Let $\mzshape$ be the m-layer Z-shape as in definition \ref{defn:m-zshape}. Let $M_{n,m}=\dfrac{1}{n^{m/2}}M_{\mzshape}(G)$ be the graph matrix where $G\sim G(n,1/2)$ and $\displaystyle r(n,m)=\dfrac{n!}{(n-m)!}$ be the dimension of $M_{n,m}$. Recall that $\displaystyle{D(m,q)=\dfrac{1}{mq+1}{(m+1)q\choose q}}$.
    Then
    \begin{equation}
        \lim_{n\to\infty}\dfrac{1}{r(n,m)}\Eb\[\trace\(\(M_{n,m}M_{n,m}^T\)^q\)\] = D(m,q).
    \end{equation}
\end{cor}

\begin{proof}
    Recall that Corollary \ref{cor:expected trace for bipartite shapes} says that for any bipartite shape $\alpha$,
    \begin{equation*}
        \lim_{n \to \infty}{\frac{1}{r_{approx}(n)}\Eb\[\trace\(\(\dfrac{M_{\alpha}M_{\alpha}^T}{n^{|V(\alpha)| - s_{\alpha}}}\)^q\)\]}
        = \abs{\left\{C \in \mathcal{C}_{\(\alpha,2q\)}: C \text{ is dominant}\right\}}.
    \end{equation*}

    Since $s_{\mzshape}=m$ and $\abs{V(\mzshape)}=2m$, $\displaystyle r_{approx}(n) = \dfrac{n!}{\(n - s_{\mzshape}\)!}=\dfrac{n!}{(n-m)!}=r(n,m)$ and $\displaystyle \frac{M_{\mzshape}M_{\mzshape}^T}{ n^{\abs{V(\mzshape)}-s_{\mzshape}}} = \dfrac{M_{\mzshape}M_{\mzshape}^T}{n^m}= M_{n,m}M_{n,m}^T$.

    By Theorem \ref{thm num of m-Z constr graph}, $\abs{\left\{C \in \mathcal{C}_{\(\mzshape,2q\)}: C \text{ is dominant}\right\}} = D(m,q)$ and the result follows.

\end{proof}

\subsection{Recurrence Relation for $D(m,n)$}

To prove the main result for this section, We need the following crucial recurrence relation for $D(m,n)$.
\begin{thm}\label{thm recurrence D(m,n)}
    \begin{equation}
        D(m,n+1)=\sum_{\substack{ i_0,\dots, i_{m}\geq 0: \\ i_0+\dots+i_{m}=n}} D(m,i_0)\dots D(m,i_{m})\,.
    \end{equation}
\end{thm}

\begin{proof}
    The proof is similar to the proof of Theorem \ref{thm recurrence Dk}. Let \emph{$W_{m,n}:=$ the set of all grid walks from $(0,0)$ to $(n,mn)$ that are weakly below the diagonal} and $d_{m,n}=\abs{W_{m,n}}$. We will prove that $d_{m,n}=D(m,n)$ and that $d_{m,n}$ satisfies the recurrence relation in the theorem.

    \begin{enumerate}[1.]

\setlength{\parindent}{1.5em}
\setlength{\parskip}{1.2mm}
\renewcommand{\baselinestretch}{1.2}

    \item $d_{m,n}=D(m,n)$:

    For $r\in\{0,1,\dots,mn\}$, let \emph{$W_{m,n}(r):=$ the set of grid walks from $(0,0)$ to $(n,mn)$ that are $r$ steps above the diagonal}. Then $\displaystyle{\bigcup_{r=0}^{mn} W_{m,n}(r)}$ is the set of all grid walks from $(0,0)$ to $(n,mn)$, which has cardinality $\displaystyle{(m+1)\cdot n\choose n}$. Also $|W_{m,n}(0)|=|W_{m,n}|=d_{m,n}$. By the same proof as in the Theorem \ref{thm recurrence Dk}, there is a bijection between $W_{m,n}(r)$ and $W_{m,n}(r-1)$ for each $r\in[mn]$. Thus $\displaystyle{d_{m,n}=\dfrac{1}{mn+1}{(m+1)\cdot n\choose n}}$.

        \begin{figure}[h]
            \centering
            \includegraphics[scale=0.3]{gridwalk-multi.png}
            \caption{Illustration of the proof part 2 for Theorem \ref{thm recurrence D(m,n)}.}
            \label{fig: D(m,n) recurrence}
        \end{figure}

    \item $\displaystyle{d_{m,n}=\sum_{\substack{ i_0,\dots, i_{m+1}\geq 0: \\ i_0+\dots+i_{m+1}=n-1}} d_{m,i_1}\dots d_{m,i_{m+1}}}$:

    Similar to the proof of Theorem \ref{thm recurrence Dk}, we now will find a bijection between $W_{m,n}$ and \\ $\displaystyle{\bigcup_{\substack{ i_0,\dots, i_{m}\geq 0: \\ i_0+\dots+i_{m}=n-1}} W_{m,i_0}\times\dots\times W_{m,i_{m}} }$.

    Let $d_k$ be the line that is shifted $k$ vertical grids down from the diagonal. i.e. $d_k$ is the line $y=m\cdot x - k$. Let $w=\(z_1,\dots,z_{n\cdot(m+1)}\)\in W_{m,n}$.

    \begin{enumerate}[-]
        \item Let $z_{i_0}=(a_0, m\cdot a_0)$ be the first point where $w$ touches the diagonal. Then $w_0:=\(z_{i_0},\dots,z_{n\cdot(m+1)}\)$ can be viewed as an element in $W_{m,n-a_0}$. Moreover, $w':=\(z_2,\dots,z_{i_0-1}\)$ is strictly below the diagonal $d_0$, thus weakly below $d_1$.

        \item Let $z_{i_1}=(a_1, m\cdot a_1-1)$ be the first point where $w'$ touches $d_1$. Then $w_1:=(z_{i_1},\dots,z_{i_0-1})$ can be viewed as an element in $W_{m,a_0-a_1}$. Moreover, $w':=(z_2,\dots,z_{i_1-1})$ is strictly below $d_1$, thus weakly below $d_2$.

        \item Continue this way, we get a sequence of points $z_{i_0},\dots,z_{i_{m-1}}$ and walks $w_0,\dots,w_{m-1}$ where each $z_{i_j}=(a_j,m\cdot a_j-j)$ is a point on $d_j$ and each $w_i$ can be viewed as an element in $W_{m,a_{i-1}-a_i}$.

        \item Since $z_{i_{m-1}}$ is the first point touching $d_{m-1}$, $w_{m}=(z_2,\dots, z_{i_{m-1}-1})$ is strictly below $d_{m-1}$ and thus weakly below $d_m$. Since $d_m$ crosses $(1,0)=z_{2}$, $w_m$ can be viewed as an element in $W_{m,a_{m-1}-1}$.
    \end{enumerate}

    Since $(n-a_0)+(a_0-a_1)+\dots+(a_{m-1}-1)=n-1$, we conclude that $$\displaystyle{(w_0,\dots,w_m)\in\bigcup_{\substack{ i_0,\dots, i_{m}\geq 0: \\ i_0+\dots+i_{m}=n-1}} W_{m,i_0}\times\dots\times W_{m,i_{m}}}\,.$$

    The other direction of the bijection can be constructed in a backward manner. It is not hard to prove this construction gives a bijection.

    \end{enumerate}

    Combining 1 and 2, we conclude that $D(m,n)$ satisfies the recurrence relation.

\end{proof}

\subsection{Properties of Dominant Constraint Graphs on $H\(\mzshape,2q\)$}

\begin{defn}\label{defn:m-zshape-H}
    Let $\mzshape$ be the multi-layer Z-shape as in Definition \ref{defn:m-zshape} and let $H\(\mzshape,2q\)$ be the multi-graph as in definition \ref{def:copies}. We label the vertices of $V_{(\a_Z)_i}$ as $\{a_{ij}, b_{ij}:j\in[m]\}$ and the vertices of $V_{(\a_Z)_i}^T$ as $\{b_{ij}, a_{(i+1)j}\}$. Let $V_i=\{a_{ij},b_{ij}:i\in[q]\}$. For $j\in[m]$, we call the induced subgraph of $H\(\mzshape,2q\)$ on vertices $V_i$ the \textit{$j^{th}$ wheel $W_j$}.

    We label the ``middle edges" of $H(\a,2q)$ in the following way: let $e_{2i-1,j}=\{a_{i(j+1)},b_{ij}\}$ and $e_{2i,j}=\{b_{ij},a_{(i+1)(j+1)}\}$ for $i=1,\dots,q$. For a fixed $j\in[m]$, we call the edges $e_{i,j}$'s the \textit{spokes between wheels $W_j$ and $W_{j+1}$} of $H\(\mzshape,2q\)$. See Figure \ref{fig:m-z-shape-H labeling} for an illustration.

    \begin{figure}[h]
        \centering
        \includegraphics[scale=0.3]{m-z-shape-H.png}
        \caption{$H\(\mzshape,2q\)$, here $m=3$}
        \label{fig:m-z-shape-H labeling}
    \end{figure}
\end{defn}

\begin{defn}\label{defn:m-z-shape induced constraint graph}
   Let $\mzshape$ be the multi-layer Z-shape. Let $C$ be a constraint graph on $H(\a_Z,2q)$. For $i=1,2$, We denote $C_i$ the induced subgraph of $C$ on vertices $V_i$. We call $C_i$ the \textit{induced constraint graph} of $C$ on $V_i$.
\end{defn}

Recall that a constraint graph $C\in\C_{(\a,2q)}$ is \emph{well-behaved} if whenever $u\=v$ in $C$, $u$ and $v$ are copies of the same vertex in $\a$ or $\a^T$.

\begin{thm}\label{thm:mzshape dominant are well behaved}
    All dominant constraint graphs in $\C_{\(\mzshape,2q\)}$ are well-behaved.
\end{thm}
\begin{proof}
    See appendix.
\end{proof}

We extend our definitions of wheel-respecting, parity-preserving, and non-crossing to $\C_{\(\mzshape,2q\)}$.

\begin{defn}\label{defn:well-defined m-z-shape constraint graph}
    Let $C$ be a constraint graph in $\C_{\(\mzshape,2q\)}$.
    \begin{enumerate}
        \item We say $C$ \textit{wheel-respecting} if for all $u\=v$, $u,v\in V_i$ for some $i\in[m]$. (i.e. no two vertices on different wheels are constrained together by $C$).
        \item We say $C$ is \textit{parity-preserving} if for each $i\in[m]$, the induced constraint graphs $C_i$ of $C$ on $V_i$ is parity-preserving.
        \item We say $C$ is \emph{non-crossing} if the induced constraint graphs $C_i$'s are non-crossing.
    \end{enumerate}

\end{defn}

\begin{rmk}
    If $C\in\C_{\(\mzshape,2q\)}$ is wheel-respecting and $G_i$ are representatives of $C_i$, then the graph $G$ with $V(G)=V(\mzshape,2q)$ and $E(G)=E(G_1)\sqcup\dots\sqcup E(G_m)$ is a representative of $C$.
\end{rmk}

\begin{prop}
    Let $C$ be a constraint graph in $\C_{\(\mzshape,2q\)}$. $C$ is well-behaved if and only if $C$ is wheel-respecting and parity-preserving.
\end{prop}

\begin{cor}\label{cor:mzshape dominant properties}
    If $C$ is a dominant constraint graph in $\C_{\(\mzshape,2q\)}$, then
    \begin{enumerate}
        \item $C$ is wheel-respecting, parity-preserving and non-crossing.
        \item The induced constraint graphs $C_i$ on wheels $W_i$ are dominant constraint graphs in $\C_{\(\a_0,2q\)}$.
    \end{enumerate}
\end{cor}

The same proofs for Lemma \ref{lem:z-shape pre-contract} yields the following statement.
\begin{lemma}\label{lem:m-zshape pre-contract}
    Let $\mzshape$ be as in definition \ref{defn:m-zshape} and let $C$ be a dominant constraint graph in $\C_{\(\mzshape,2q\)}$. If $a_{sj}\=a_{tj}$ for some $1\leq s < t\leq q$ and $j\in[m-1]$, then $a_{s(j+1)}\=a_{t(j+1)}$. Moreover, the spokes $\left\{e_{x,j}:x \in [2s-1,2t-2]\right\}$ can only be made equal to each other.

    Similarly, if $b_{sj}\=b_{tj}$ for some $1\leq s< t\leq q$ and $j\in \{2,3,\dots,m\}$, then $b_{s(j-1)}\=b_{t(j-1)}$. Moreover, the spokes $\left\{e_{x,j-1}:x \in [2s,2t-1]\right\}$ can only be made equal to each other.

\end{lemma}

\begin{cor}\label{cor:m-zshape split}
\begin{figure}[h]
    \centering

    \includegraphics[scale=0.3]{m-zshape-split-shade.png}

    \caption{Illustration of Corollary \ref{cor:m-zshape split}: $C\in\C_{\(\mzshape,2q\)}$ is dominant. $a_{ik}\longleftrightarrow a_{jk}$ and $b_{i(k-1)} \longleftrightarrow b_{j(k-1)}$. Here $q=7$, $m=4$, $k=3$ and $j-i=3$.}
    \label{fig:m-zshape split}
\end{figure}
    If $C$ is a dominant constraint graph in $\C_{\(\mzshape,2q\)}$, then the following statements are true:
    \begin{enumerate}
        \item If $a_{i1}\=a_{j1}$ for some $1\leq i<j\leq q$, then $a_{ik}\=a_{jk}$ for all $k\in[m]$.

        \item Similarly, if $b_{im}\=a_{jm}$ for some $1\leq i<j\leq q$, then $b_{ik}\=b_{jk}$ for all $k\in[m]$.

        \item More generally, if $a_{ik}\=a_{jk}$ and $b_{i(k-1)}\=b_{j(k-1)}$ for some $k\in[m+1]$, then $a_{is}\=a_{js}$ for all $k\leq s\leq m$ and $b_{it}\=b_{jt}$ for all $1\leq t\leq k-1$. Note that $k=1$ corresponds to case 1 above and $k=m+1$ corresponds to case 2 above. See Figure \ref{fig:m-zshape split} for an illustration.
    \end{enumerate}

    In all three cases, contracting the constrained vertices splits $H\(\mzshape,2q\)$ into $H\(\mzshape,2(j-i)\)$ and $H\(\mzshape,2(q-j+i)\)$, and the induced constraint graphs $C'\in\C_{\(\mzshape,2(j-i)\)}$ and $C''\in\C_{\(\mzshape,2(q-j+i)\)}$ are dominant.
\end{cor}
\begin{proof}[Proof sketch]
We focus on the third case as this is the trickiest case. To split $H\(\mzshape,2q\)$ into $H\(\mzshape,2(j-i)\)$ and $H\(\mzshape,2(q-j+i)\)$, imagine doing the following:
\begin{enumerate}
    \item Cut towards the center of the wheels through the vertices $b_{i1},b_{i2},\ldots,b_{i(k-1)}$, then cut along the spoke $\{a_{ik},b_{i(k-1)}\}$, and then cut towards the center of the wheels through the vertices $a_{ik},...a_{im}$.
    \item Similarly, cut towards the center of the wheels through the vertices $b_{j1},b_{j2},\ldots,b_{j(k-1)}$, then cut along the spoke $\{a_{jk},b_{j(k-1)}\}$, and then cut towards the center of the wheels through the vertices $a_{jk},...a_{jm}$.
\end{enumerate}
These cuts split $H\(\mzshape,2q\)$ into two halves. Taking each half and gluing it to itself along the cuts gives us $H\(\mzshape,2(j-i)\)$ and $H\(\mzshape,2(q-j+i)\)$.

To see that these the induced constraint graphs $C'$ and $C''$ are dominant, observe that except for the vertices along the cuts, each vertex appears in one half or the other but not both. Except for the spokes $\{a_{ik},b_{i(k-1)}\}$ and $\{a_{jk},b_{j(k-1)}\}$, each edge is incident to a vertex which is not part of the cut, so this implies that all of the edges in $H\(\mzshape,2q\)$ except for the spokes $\{a_{ik},b_{i(k-1)}\}$ and $\{a_{jk},b_{j(k-1)}\}$ (which are made equal to each other) appear in one half or the other but not both.

By Corollary \ref{cor:line-shape dividing constraint edge}, edges on wheels on one side of the split can only be made equal to the edges on the same side. By Lemma \ref{lem:z-shape pre-contract}, $b_{ir}\=b_{jr}$ implies that spokes between wheels $W_r$ and $W_{r-1}$ on one side of the split can only be made equal to the spokes on the same side. By the same lemma, $a_{ir}\=a_{jr}$ implies that spokes between wheels $W_{r}$ and $W_{r+1}$ on one side of the split can only be made equal to the spokes on the same side. By the assumptions on the constrained vertices, we conclude that other than the spokes between wheels $W_{k}$ and $W_{k-1}$, all other spokes are made equal to the spokes on its own side after splitting.

Thus so far, in both $H\(\mzshape,2(j-i)\)/{C'}$ and $H\(\mzshape,2(q-j+i)\)/{C''}$, each edge except for the spokes between $W_k$ and $W_{k-1}$ appears an even number of times.

For the spokes between $W_k$ and $W_{k-1}$, assume there is a spoke $e_s$ on one side of the split such that $e_s\not\=\{a_{ik},b_{i(k-1)}\}$ and $e_s\=e_t$ for some $e_t$ on the other side of the split. Then by Lemma \ref{lem:non-crossing spokes}, $e_s\=e_t\=e_{2i-1,k-1}\=e_{2j-1,k-1}$, a contradiction. Thus spokes on one side that are not made equal to $\{a_{ik},b_{i(k-1)}\}$ or $\{a_{jk},b_{j(k-1)}\}$ are only made equal to other spokes on the same side. Since the total number of spokes on each side is even, $\{a_{ik},b_{i(k-1)}\}$ needs to appear even number of times in $H\(\mzshape,2(j-i)\)/{C'}$ and $\{a_{jk},b_{j(k-1)}\}$ needs to appear even number of times in $H\(\mzshape,2(q-j+i)\)/{C''}$. Thus in both $H\(\mzshape,2(j-i)\)/{C'}$ and $H\(\mzshape,2(q-j+i)\)/{C''}$, each spoke between $W_k$ and $W_{k-1}$ appears an even number of times .

This implies that $C'$ and $C''$ are nonzero-valued constraint graphs on $H\(\mzshape,2(j-i)\)$ and $H\(\mzshape,2(q-j+i)\)$ and are thus dominant, as needed.
\end{proof}

\subsection{Proof of Theorem \ref{thm num of m-Z constr graph}}

Now we are ready to prove the main result of this section, Theorem \ref{thm num of m-Z constr graph}. The proof will be similar to the one for the Z-shape case. For the Z-shape, we split $H\(\a_Z,2q\)$ two times to get the recurrence relation, where first step is based on the constrained vertices on the inner wheel, and the second step is based on the constrained vertices on the outer wheel. For the multi-Z-shape with $m$ layers, we will split $H\(\mzshape,2q\)$ $m$ times for the recurrence relation, where the $r^{th}$ step will be based on the constrained vertices on the $(m-r+1)^{th}$ wheel, starting from the inner-most wheel.

\begin{defn}\label{defn:D(m,q,r)}
    For $r\in[m]$, let $\D_{m,q,r}$ denote the set of all dominant constraint graphs in $\C_{\(\mzshape,2q\)}$ such that $a_{1m}\=a_{2m}, a_{1(m-1)}\=a_{2(m-1)},\dots, a_{1(m-r+1)}\=a_{2(m-r+1)}$ in $C$. When $r=0$, define $\D_{m,q,0}=$ the set of all dominant constraint graphs in $\C_{\(\mzshape,2q\)}$. Define $D(m,q,r)=\abs{\D_{m,q,r}}$ for all $m,q,r$.
\end{defn}

We first consider the first step of the splitting.
\begin{lemma}\label{lem:m-zshape splitting first step}
    Let $\D_{m,q,r}$ be as defined in \ref{defn:D(m,q,r)}. There is a bijection between $\displaystyle\D_{m,q,0}$ and $\displaystyle\bigcup_{i=0}^{q-1} \D_{m,i,0}\times \D_{m,q-i,1}$. Thus
    \begin{equation}
        D(m,q,0)=\sum_{i=0}^{q-1} D(m,i,0)\cdot D(m,q-i,1).
    \end{equation}
\end{lemma}

\begin{proof}

    \begin{figure}[h]
        \centering
        \includegraphics[scale=0.3]{m-zshape-splitting.png}
        \caption{Illustration of Lemma \ref{lem:m-zshape splitting first step}: $a_{(i+1)m}$ is the first vertex that $a_{1m}$ is constrained to. Constrained vertices are indicated by dash lines. Here $m=3$.}
        \label{fig: m-zshape splitting step 1}
    \end{figure}

    Let $C\in\D_{m,q,0}$. i.e., $C$ is a dominant constraint graph on $H\(\mzshape,2q\)$. Let $i>0$ be the smallest index such that $a_{1m}$ is constrained to $a_{(i+1)m}$. If $a_{1m}$ is isolated then we take $i=q$.

    Since $C$ is dominant, by Corollary \ref{cor:mzshape dominant properties}, the induced constraint graph $C_m$ on $W_m$ is dominant. By Lemma \ref{lem:line shape first constrained vertex} and Lemma \ref{lem:line shape isolated vertex} (in the case when $a_{1m}$ is isolated), $b_{1m}\=b_{im}$. By Corollary \ref{cor:m-zshape split}, $H\(\mzshape,2q\)$ splits into two parts, $H\(\mzshape,2(i-1)\)$ and $H\(\mzshape,2(q-i+1)\)$ where $a_{1m}\=a_{(i+1)m}$ (see Figure \ref{fig: m-zshape splitting step 1} for an illustration). Moreover, the induced constraint graphs $C'\in\C_{\(\mzshape,2(i-1)\)}$ and $C''\in\C_{\mzshape,2(q-i+1)}$ are dominant.

    Since $C'$ is dominant on $H(\mzshape,2(i-1))$, we have $C'\in\D_{m,i-1,0}\,$. Since $C''$ is dominant on $H\(\mzshape,2(q-i+1)\)$ and $a_{1m}\=a_{(i+1)m}$ where $a_{1m}$ and $a_{(i+1)m}$ are adjacent in $H\(\mzshape,2(q-i+1)\)$ (see Figure \ref{fig: m-zshape splitting step 1} for an illustration), $C''\in\D_{m,q-i+1,1}$.

    So far we constructed a mapping $C\in\D_{m,q,0}\mapsto (C',C'')\in\D_{m,i-1,0}\times\D_{m,q-i+1,1}$ given that $i\in[q]$ is the smallest index such that $a_{1m}\=a_{(i+1)m}$ in $C$. Thus by considering the first vertex that $a_{1m}$ is constrained to in $C$, we can construct a map from $\D_{m,q,0}$ to $\displaystyle\bigcup_{i=0}^{q-1} \D_{m,i,0}\times\D_{m,q-i,1}$.

    The other direction of the map goes in reverse order of the contracting process. It can be verified that this is a bijection. Thus we proved that $\displaystyle D(m,q,0)=\sum_{i=0}^{q-1} D(m,i,0)\cdot D(m,q-i,1)$.
\end{proof}

Now we will consider the later steps of the splitting.
\begin{lemma}\label{lem:m-zshape intermediate splitting}
    For $r\geq 1$, there exists an bijection between $\D_{m,q,r}$ and $\displaystyle{\bigcup_{i=0}^{q-1}\D_{m,i,0}\times \D_{m,q-i,r+1}}$. Thus
    \begin{equation}
        D(m,q,r)=\sum_{i=0}^{q-1} D(m,i,0)\cdot D(m,q-i,r+1).
    \end{equation}
\end{lemma}

\begin{proof}
    \begin{figure}[h]
        \centering
        \includegraphics[scale=0.3]{m-zshape-intermediate-splitting-shaded.png}
        \caption{Illustration of proof for Lemma \ref{lem:m-zshape intermediate splitting}: On the left is a $C\in\D_{m,q,r}$, with constrained vertices indicated by dash lines. Here $q=7$, $m=3$ and $r=1$. $b_{i(m-r)}$ is the first vertex that $b_{1(m-r)}$ is constrained to in $C$. On the right is $H\(\mzshape,2(q-i+1)\)$ and $H\(\mzshape,2(i-1)\)$ after splitting.}
        \label{fig: m-zshape splitting step 2}
    \end{figure}

    Let $C\in\D_{m,q,r}$ for some $r\geq 1$. i.e. $C$ is dominant and $a_{1k}\=a_{2k}$ for all $m-r+1\leq k\leq m$. Let $i>1$ be the first index such that $b_{1(m-r)}\=b_{i(m-r)}$. If $b_{1(m-r)}$ is isolated then we take $i=q+1$. By Corollary \ref{cor:mzshape dominant properties}, since $C$ is dominant on $H\(\mzshape,2q\)$, the induced constraint graph $C_{m-r}$ is dominant, thus by Lemma \ref{lem:line shape first constrained vertex} and Lemma \ref{lem:line shape isolated vertex}, $a_{2(m-r)}\=a_{i(m-r)}$. By Lemma \ref{lem:m-zshape pre-contract}, $a_{2k)}\=a_{ik}$ for all $m-r\leq k\leq m$. Since $a_{1k}\=a_{2k}$ for all $m-r+1\leq k\leq m$ by assumption, $a_{1k}\= a_{ik}$ for all $m-r+1\leq k\leq m$. By Lemma \ref{lem:m-zshape pre-contract}, since $b_{1(m-r)}\=b_{i(m-r)}$, $b_{1k}\=b_{ik}$ for all $1\leq k\leq m-r$.

    So far we have that
    \begin{enumerate}
        \item $a_{1k}\= a_{2k}\= a_{ik}$ for all $m-r+1\leq k\leq m$.
        \item $b_{1k}\=b_{ik}$ for all $1\leq k\leq m-r$.
        \item $a_{2(m-r)}\=a_{i(m-r)}$.
    \end{enumerate}

    Since $b_{1(m-r)}\=b_{i(m-r)}$ and $a_{1(m-r+1)}\= a_{i(m-r+1)}$, by Corollary \ref{cor:m-zshape split}, $H\(\mzshape,2q\)$ splits into $H\(\mzshape,2(i-1)\)$ and $H\(\mzshape,2(q-i+1)\)$. Let $C'\in\C_{\(\mzshape,2(i-1)\)}$ and $C''\in\C_{\(\mzshape,2(q-i+1)\)}$ be the induced constraint graphs. By Corollary \ref{cor:line-shape dividing constraint edge}, the $C'$ and $C''$ are dominant. Notice that there is no extra required constrained vertices in $C''$, so $C''\in\D_{m,q-i+1,0}$. For $C'$, since $a_{2k}\=a_{ik}$ for all $m-r+1\leq k\leq m$ and $a_{2(m-r)}\=a_{i(m-r)}$, $C'\in\D_{m,i-1,r+1}$ (see Figure \ref{fig: m-zshape splitting step 2} for an illustration).

    So far we have constructed a mapping $C\in\D_{m,q,r}\mapsto (C',C'')\in\D_{m,q-i+1,0}\times\D_{m,i-1,r+1}$ given that $1<i<q+1$ is the first index such that $b_{1(m-r)}\=b_{i(m-r)}$ in $C$. Thus by considering the first vertex that $b_{1(m-r)}$ is constrained to in $C$, we can construct a mapping from $\D_{m,q,r}$ to $\displaystyle\bigcup_{i=0}^{q-1}\D_{m,i,0}\times\D_{m,q-i,r+1}$.

    The other direction of the map simply reverses the the splitting process. It can be verified that this is a bijection. Thus we proved that $\displaystyle D(m,q,r)=\sum_{i=0}^{q-1} D(m,i,0)\cdot D(m,q-i,r+1)$.
\end{proof}

Finally we arrive at the last step, where we identify the last item from the splitting process $\D_{m,q,m}$ with $\D_{m,q-1,0}$.
\begin{lemma}\label{lem:m-zshape contract all wheels}
    Let $\D_{m,q,r}$ be defined as in \ref{defn:D(m,q,r)}. For any $q\geq 1$ and $m\geq1$, there is a bijection between $\D_{m,q,m}$ and $\D_{m,q-1,0}$. Thus $D(m,q,m)=D(m,q-1,0)$.
\end{lemma}
\begin{proof}
    Let $C\in\D_{m,q,m}$. Then $a_{1j}\=a_{2j}$ for all $j\in[m]$. Contracting the constrained vertices we get the induced constraint graph $C'$ on $H(\mzshape,2(q-1))$, and $C'$ is dominant. i.e. $C'\in\D_{m,q-1,0}$.

    Given $C'\in\D_{m,q-1,0}$, we first relabel the indices by increasing all of them by $1$. We then expand $a_{2j}$ to $a_{1j}$ and $a_{2j}$ and make them constrained, for each $j\in[m]$. Adding these constraints $a_{1j}\=a_{2j}$ for all $j\in[m]$ to $C'$, we get a new constraint graph $C\in\D_{m,q,m}$.

    We have constructed maps between $\D_{m,q,m}$ and $\D_{m,q-1,0}$ in both directions. It is easy to see that this is a bijection. Thus we proved that $D(m,q,m)=D(m,q-1,0)$.
\end{proof}

With the above lemmas, we are ready to prove Theorem \ref{thm num of m-Z constr graph}.
\begin{proof}[Proof of Theorem \ref{thm num of m-Z constr graph}]

Recall that $\displaystyle{D(m,q)=\dfrac{1}{mq+1}{(m+1)q\choose q}}$ and we want to prove that $D(m,q,0)=$ the number of dominant constraint graphs in $\C_{\(\mzshape,2q\)}$ $=D(m,q)$. For $q=1$, we check that $D(m,1,0)=1=D(m,1)$. We will then show that
\begin{equation}\label{eqn M recurrence}
    D(m,q,0)=\sum_{\substack{i_1,\dots,i_{m+1}\geq 0\\ i_1+\dots+i_{m+1}=q-1}} D(m,i_1,0)\dots D(m,i_{m+1},0).
\end{equation}
which is the same recurrence relation for $D(m,q)$ shown in Theorem \ref{thm recurrence D(m,n)}. Thus by induction on $q$, we can prove $D(m,q,0)=D(m,q)$ for all $q\geq 1$.

Therefore it suffices to prove Equation \eqref{eqn M recurrence}.

Recall that from Lemma \ref{lem:m-zshape splitting first step}, Lemma \ref{lem:m-zshape intermediate splitting} and Lemma \ref{lem:m-zshape contract all wheels}, we have
\begin{enumerate}
    \item $\displaystyle D(m,q,0)=\sum_{i=0}^{q-1} D(m,i,0)\cdot D(m,q-i,1)$.
    \item $\displaystyle D(m,q,r) = \sum_{i=0}^{q-1} D(m,i,0)\cdot D(m,q-i,r+1)$.
    \item $\displaystyle D(m,q,m)=D(m,q-1,0)$.
\end{enumerate}

Thus
\begin{align*}
    D(m,q,0)
    &=\sum_{i_1=0}^{q-1} D(m,i_1,0)\cdot D(m,q-i_1,1) \\
    &= \sum_{i_1=0}^{q-1} D(m,i_1,0)\cdot\(\sum_{i_2=0}^{q-i_1-1} D(m,i_2,0)\cdot D(q-i_1-i_2,2)\) \\
    & \vdots \\
    & = \sum_{\substack{i_1,\dots,i_m\geq 0, i_{m+1}'\geq 1: \\ i_1+\dots+i_m+i_{m+1}'=q}} D(m,i_1,0)\cdot D(m,i_2,0)\dots D(m,i_{m},0)\cdot D(m,i_{m+1}',m)\\
    & = \sum_{\substack{i_1,\dots,i_m\geq 0, i_{m+1}'\geq 1: \\ i_1+\dots+i_m+i_{m+1}'=q}} D(m,i_1,0)\dots D(m,i_{m},0)\cdot D(m,i_{m+1}'-1,0) \\
    & = \sum_{i_j\geq 0: i_1+\dots+i_{m+1}=q-1} D(m,i_1,0)\dots D(m,i_{m},0)\cdot D(m,i_{m+1},0)
\end{align*}

This proves Equation \eqref{eqn M recurrence}, as needed.
\end{proof}

\end{document}


TITLE: How do we know the temperature on the planets?
QUESTION [12 upvotes]: I was watching a show and they were saying that the temperature of Pluto (I know it is not a planet) is about $-300 {}^\circ F$.  I know that depends where in the orbit Pluto is, but how do we determine the temperature?

REPLY [20 votes]: There's a few ways the temperature can be measured remotely.
The easiest way is to measure the amount, and for bonus, spectrum, of the radiated heat. All objects greater than Absolute Zero radiate a certain amount of energy. The wavelength spectrum can be determined by Planck's Law, and the amount of energy by
the Stefan-Boltzmann law, both of which are listed below.
$I'(\lambda,T) =\frac{2 hc^2}{\lambda^5}\frac{1}{ e^{\frac{hc}{\lambda kT}}-1}$
$j^{\star} = \sigma T^{4}$
So, to determine the temperature, all one needs to do is measure the energy over the waveband of interest, and one can get the temperature.
In the case of Pluto, an easier method can be done. We know the total amount of energy going toward Pluto, based off of its distance from the Sun. We can use this to measure the average temperature of an object receiving that amount of energy.
The last method involves knowing the physical properties of the object. For instance, Pluto has a atmosphere of Nitrogen, Carbon Dioxide, and Methane. The atmosphere will freeze as Pluto moves further out in its orbit. The temperature can be inferred by what materials haven't frozen out of the atmosphere yet, which can be determined with spectroscopy.
Typically, the first one is done, and the second when a more accurate method cannot be used (IE, no direct spectrum of the planet observed).


TITLE: Division with decimal in the divisor
QUESTION [2 upvotes]: I understand HOW to do division with a decimal in the divisor, but my question is MUST we remove the decimal in the divisor and if so, why?
Thanks.

REPLY [1 votes]: Remember that the word decimal is a shortened form of decimal fraction, meaning a fraction whose denominator is a power of $10$.  For example,
$$
3.14 = \frac{314}{100},
\qquad
0.0625 = \frac{625}{1000},
\qquad
5280.1 = \frac{52801}{10}.
$$

As an example, consider the division problem
$$
\frac{9.6}{3.84} = \frac{\frac{96}{10}}{\frac{384}{100}}
= \frac{96}{10} \cdot \frac{100}{384}
= \frac{96}{384} \cdot \frac{100}{10}
= \frac{96}{384} \cdot 10.
$$
This last form shows you that the original division calculation is equivalent the quotient of the whole numbers$96$ and $384$, as long as you adjust by a factor of $10$ afterwards.  (This is usually described as moving the decimal point.)
In order to actually do this calculation, you can multiply the numerator by a sufficiently large power of $10$ to make it a multiple of the denominator, making sure to divide by that same power of $10$ afterwards to compensate.  Note:  it won't always be the case that you get precisely a multiple, no matter how large a power of $10$ you multiply by.  In this case,
$$
\frac{96}{384} \cdot 10 = \frac{96}{384} \cdot \frac{100}{100} \cdot 10
= \frac{9600}{384} \cdot \frac{10}{100}
= 25 \cdot \frac{1}{10}
= \frac{25}{10}
= 2.5.
$$
The upshot:  the integer fact that $384 \cdot 25 = 9600$ is responsible for all of the following:
$$
\cdots
= \frac{0.0096}{0.00384}
= \frac{0.096}{0.0384}
= \frac{0.96}{0.384}
= \frac{9.6}{3.84}
= \frac{96}{38.4}
= \frac{960}{384}
= \frac{9600}{3840}
= \cdots
= 2.5
$$


TITLE: How to prove $e^x\left|\int_x^{x+1}\sin(e^t) \,\mathrm d t\right|\le 1.4$?
QUESTION [20 upvotes]: Related question asked by me on MathOverflow: How to prove $e^x\left|\int_x^{x+1}\sin(e^t) \,\mathrm d t\right|\le 1.4$?
This is a follow-up question to the question How to prove $ \mathrm{e}^x\left|\int_x^{x+1}\sin\mathrm e^t \mathrm d t\right|\leqslant 2$?, in which a weaker bound is proven using a nice trick.
Now my question is how to maximize and minimize $$f(x)=e^x\int_x^{x+1}\sin(e^t) \,\mathrm d t$$
or at least to prove $-1.4\le f(x)\le 1.4$.
Some observations, using the substitution $y=e^t$:
$$f(x)=e^x \int_{e^x}^{e^{x+1}} \frac{\sin(y)}y\,\mathrm dy=g(e^x),$$
where I have defined $$g(z)=z \int_z^{e z} \frac{\sin(y)}y\,\mathrm dy = z (\operatorname{Si}(e z)-\operatorname{Si}(z)).$$
($\operatorname{Si}$ is the Sine integral.)
So the question reduces to: What are the maxima/minima of $g(z)$ for $z\geq 0$ ?
Using the series of $\mathrm{Si}(z)$, we get
$$g(z)=\sum_{k=1}^\infty (-1)^{k-1} \frac{z^{2k}(e^{2k-1}-1)}{(2k-1)!\cdot(2k-1)}$$
and here is a plot of $g(z)$, which seems to be periodic:

REPLY [2 votes]: By the result of [1], we have
$$\left|\mathrm{Si}(y) - \frac{\pi}{2}\right| \le \frac{1}{y}, \ y > 0.\tag{1}$$
Thus, we have, for $z > 0$,
$$\frac{\pi}{2} - \frac{1}{z} \le \mathrm{Si}(z) \le \frac{\pi}{2} + \frac{1}{z}$$
and
$$\frac{\pi}{2} - \frac{1}{\mathrm{e}z} \le \mathrm{Si}(\mathrm{e}z) \le \frac{\pi}{2} + \frac{1}{\mathrm{e}z}.$$
Thus, we have, for $z > 0$,
$$-1 - \frac{1}{\mathrm{e}} \le z(\mathrm{Si}(\mathrm{e}z) - \mathrm{Si}(z)) \le 1 + \frac{1}{\mathrm{e}}.$$
Also, $1 + \frac{1}{\mathrm{e}} \approx 1.367879441$. We are done.
Reference
[1] Upper bound for the sine integral
Remark: The following stronger inequality holds, which implies (1):
$$\left|\mathrm{Si}(y) - \frac{\pi}{2}\right| \le \arctan \frac{1}{y}, \ y > 0.$$
See: Graham Jameson, Nick Lord and James McKee, An inequality for Si(x), Math. Gazette 99 (2015).
https://www.maths.lancs.ac.uk/jameson/siineqnotes.pdf


text:
\begin{document}

\maketitle

\begin{abstract}
Bayesian quadrature (BQ) is a method for solving numerical integration problems in a Bayesian manner, which allows users to quantify their uncertainty about the solution.
The standard approach to BQ is based on a Gaussian process (GP) approximation of the integrand. As a result, BQ is inherently limited to cases where GP approximations can be done in an efficient manner, thus often prohibiting very high-dimensional or non-smooth target functions.
This paper proposes to tackle this issue with a new Bayesian numerical integration algorithm based on Bayesian Additive Regression Trees (BART) priors, which we call BART-Int.
BART priors are easy to tune and well-suited for discontinuous functions.
We demonstrate that they also lend themselves naturally to a sequential design setting and that explicit convergence rates can be obtained in a variety of settings.
The advantages and disadvantages of this new methodology are highlighted on a set of benchmark tests including the Genz functions, and on a Bayesian survey design problem.
\end{abstract}

\section{Introduction}
\label{Introduction}

Numerical integration algorithms are key tools for modern statistics and machine learning. They allow us to tackle many of the integrals commonly arising in those fields and which cannot be computed in closed form. Examples include the computation of posterior expectations in Bayesian statistics, the marginalisation or conditioning of random variables, the computation of normalisation constants, calculations for the EM algorithm, or even the approximation of solutions to differential equations.

This paper will focus on the task of approximating the integral of some function $f:\mathcal{X} \rightarrow \mathbb{R}$ which is integrable with respect to some distribution $\Pi$ (whose density with respect to the Lebesgue measure is denoted $\pi$) over some domain $\mathcal{X}\subseteq \mathbb{R}^d$ ($d\in \mathbb{N}_+ = \mathbb{N}\backslash \{ 0\}$):
\begin{talign}
	\Pi[f] :=\int_{\mathcal{X}} f(x) d\Pi(x) = \int_{\mathcal{X}} f(x) \pi(x) dx.
\label{eq:integral}
\end{talign}
A large number of methods have been developed to tackle this problem.
Classical quadrature rules \cite{Davis2007} tend to be limited to low-dimensional problems and integrals against a small class of probability measures. Alternatively, Monte Carlo integration (MI) \cite{Robert2004} only requires sampling independently from $\Pi$ and enjoys a convergence rate of $\mathcal{O}(n^{-1/2})$. When $\pi$ is unnormalised, sequential Monte Carlo (SMC) \cite{Moral2006} or Markov chain Monte Carlo (MCMC) \cite{Robert2004} samplers can be used and enjoy a similar rate.  Quasi-Monte Carlo (QMC) \cite{Dick2013} can improve on this convergence rate, but is limited to integration against a uniform measure in a hyper-cube (or simple transformations of this problem). We note that all of these methods are \emph{quadrature rules}: they take the form $\hat{\Pi}[f]:=\sum_{i=1}^n w_i f(x_i)$ where $\{x_i \}_{i=1}^n\subset \mathcal{X}$ are design points and $\{ w_i\}_{i=1}^n\subset \mathbb{R}$ are weights. Finally, the Laplace approximation \cite{Rue2016} and variational inference \cite{Blei2017} can be used, but these cannot be guaranteed to provide accurate estimates as $n$ grows.

Although all of the methods above are commonly used in practice, they all lack a straightforward non-asymptotic approach to quantifying our uncertainty about the value of $\Pi[f]$ after a finite number of function evaluations. This is not necessarily a problem for applications where $n$ can be taken to be large relative to the difficulty of the problem. However, those methods will be sub-optimal when $f$ is both difficult to approximate and expensive to evaluate since the asymptotic results will not hold.

An alternative approach comes from the field of probabilistic numerics \cite{Suldin1959,Larkin1972,Diaconis1988,OHagan1992,Hennig2015}, and in particular Bayesian probabilistic numerical methods \cite{Cockayne2017BPNM}, which frame problems in numerical analysis as statistical estimation tasks. This allows for the quantification of uncertainty surrounding the value of a quantity of interest using probabilistic statements valid for finite $n$. For integration, the main approach is called Bayesian quadrature (BQ) \cite{OHagan1991,Diaconis1988,Rasmussen:2002:BMC:2968618.2968681,Briol2019PI}. Its name comes from the fact that the estimator takes the form of a quadrature rule, but also carries a Bayesian interpretation. However, we will see that not all Bayesian estimators have this property. For this reason, we prefer the general name of \textit{Bayesian probabilistic numerical integration (BPNI)}.

BPNI starts by positing a prior distribution for the integrand $f$. It then computes the posterior distribution given values of $f$ at some design points, and finally considers the implied distribution on $\Pi[f]$. The main advantage of BPNI is that all our knowledge about the problem can be straightforwardly implemented in a prior, and the posterior distribution on $\Pi[f]$ can allow us to quantify our uncertainty about the exact value of this quantity. BPNI critically relies upon a flexible statistical model for the integrand $f$, which is almost always chosen to be a Gaussian process (GP). In this case, we have a quadrature rule which we will call GP-BQ. GPs have the convenient property that the posterior distribution can be obtained in closed form for interpolation and regression with Gaussian noise. They also have well-studied concentration rates which carry over to the corresponding GP-BQ methods \cite{Briol2019PI,Kanagawa2017,Kanagawa2019,Wynne2020}.

GP-BQ has received a lot of attention in recent years. Efficient deterministic \cite{Karvonen2017symmetric,Karvonen2019},  adaptive \cite{Gunter2014,briol2015frank,Fisher2020,Jiang2019} and randomised \cite{belhadji2019kernel} point-selection schemes have been designed. The method was extended to cases where $\pi$ is not available \cite{Oates2017heart}, and when multiple integrals are computed simultaneously \cite{Xi2018MultiOutput,Gessner2019}. It has been applied to fields ranging from econometrics \cite{Oettershagen2017} to computer graphics \cite{Brouillat2009} and robotics \cite{Quann2018}. While GPs are virtually the only choice that has been explored in the BPNI literature, in practice they suffer from a number of challenges. These include:

\begin{enumerate}

    \item \textbf{Discontinuities:} GPs are not suited for integration problems where the integrand is discontinuous -- a common challenge in applications \cite{Cornford1999,Mohammadi2019}. As a result, they may lead to poor predictions and unreliable uncertainty estimates for these problems.

    \item \textbf{Computational cost:} A major issue with GPs is their cubic cost in the number of data points. This can be mitigated using approximate GPs, but these do not necessarily lead to closed-form estimators. Exceptions include the work of \cite{Karvonen2017symmetric} and \cite{Jagadeeswaran2018} which allow for near-linear computational cost, but require $\Pi$ and $\{x_i\}_{i=1}^n$ to exhibit symmetry properties which may not hold in practice.

    \item \textbf{High dimensions:} Few applications of GP-BQ exists where $d>10$ due to the curse of dimensionality. This issue is closely linked to the computational cost issue since the number of points needed to approximate functions well will grow exponentially with dimension. However, it is sometimes possible to use sparsity of $f$ to lower this cost \cite{Briol2019PI}.

\end{enumerate}

To tackle some of these issues, we propose to use an alternative model based on trees. The specific model is called Bayesian additive regression trees (BART) \cite{BART,Hill2020}, and as a result we call this BPNI method \textit{BART integration (BART-Int)}. BART is a sum-of-trees model similar in spirit and effectiveness to random forests. It has been successfully deployed in a number of settings including for causal inference \cite{Hill2011, dorie2019automated}, genomics \cite{Ding2012}, behavioural sciences \cite{yeager2019national} and Bayesian optimisation \cite{Chipman2012}. Extensive experimental results in the literature have shown that BART does not usually overfit, thanks to its ensemble structure, and that it only requires a limited amount of hyperparameter tuning.

\textit{\textbf{Contributions of the paper:}}
This paper derives a novel BPNI algorithm based on Bayesian priors with a tree structure. We show that BART is a natural choice for BPNI which can be preferable to GPs in any of the three settings highlighted. This is done through a mix of theory and experiments. On the theoretical side, we prove asymptotic convergence rates for BPNI with MCMC approximations. This result may be of independent interest since it is applicable for any nonparametric method for which a concentration rate is known. On the experimental side, we compare our method to GP-BQ on a range of problems including the Genz test functions, a rare-event simulation problem and a Bayesian survey design problem.

\section{Background: Bayesian Regression with GPs and BART}
\label{background}

\label{bart_functino}

Before presenting BART-Int, this section introduces background material on Bayesian regression with GPs and BART. Given data $(X,y)$ consisting of design points $X = (x_1,\ldots,x_n)^\top \subset \mathcal{X}^n$ and responses $y = (y_1,\ldots,y_n)^\top \subset \mathbb{R}^n$, the regression problem is to recover an unknown function $f:\mathcal{X}\rightarrow\mathbb{R}$ for which we have access to noisy measurements $y$. These are usually assumed to be independently and identically distributed (i.i.d.\@) Gaussian: $y_i = f(x_i) + \epsilon_i$, $\epsilon_i \sim \mathcal{N}(0,\sigma^{2})$ for $i=1,\ldots,n$, where $\mathcal{N}(a,b)$ denotes a Gaussian distribution with mean $a$ and variance $b$, and $N(\cdot |a, b)$ will denote its density.

Bayesian (nonparametric) regression \cite{Ghosal2017} consists of specifying a (nonparametric) prior distribution over the function $f$ and conditioning on the observations  to obtain a posterior distribution on $f$. This prior is often chosen to be some (real-valued) stochastic process $g:\mathcal{X} \times \Omega \rightarrow \mathbb{R}$ where $(\Omega,\Sigma,\mathbb{P})$ is some probability space. These can be thought of as a random function since $\forall \omega \in \Omega$, $g(\cdot;\omega)$ will be a function, and $\forall x \in \mathcal{X}$, $g(x;\cdot)$ will be a random variable.

A common choice of prior is a GP \cite{Rasmussen2006}. A GP is fully determined by its mean function $\mu:\mathcal{X}\rightarrow \mathbb{R}$ and covariance function (or kernel) $k:\mathcal{X}\times \mathcal{X}\rightarrow \mathbb{R}$, and hence often denoted by $\mathcal{GP}(\mu, k)$. Due to conjugacy properties of the Gaussian distribution, the posterior $f$ after conditioning on $(X,y)$ is once again a GP with posterior mean $\tilde{\mu}(x) = \mu(x) + k_{x,X} (k_{X,X}+\sigma^2 I)^{-1}(y - \mu_X)$ and covariance $\tilde{k}(x, x') = k(x,x') - k_{x, X} (k_{X,X}+\sigma^2 I)^{-1} k_{X, x'}$, where $k_{X,X}\coloneqq (k(x_i, x_j))_{i,j=1}^n$, $k_{x, X}\coloneqq (k(x,x_1), \ldots, k(x, x_n))$, $k_{X,x}=k_{x,X}^\top$, $\mu_X \coloneqq (\mu(x_1),\ldots,\mu(x_n))^\top$ and $I$ is the $n \times n$ identity matrix. The properties of GP posteriors are inherited from $\mu$ and $k$. Whenever these are continuous, the posterior mean will also be continuous. This is the case for the most common choices of kernels, such as the Gaussian or Mat\'ern kernels with sufficient smoothness. In those cases, discontinuities would have to be inputted manually or inferred from data \cite{Cornford1999}, making the modelling of functions with a large number of discontinuities a challenge for most GP models.

In this paper, we consider instead models based on tree structures. Popular examples include Bayesian CART \cite{Denison1998}, dynamic regression trees \cite{Taddy2011} and Mondrian forests \cite{Lakshminarayanan2016}. We will focus on BART \cite{BART,Hill2020} due to its strong empirical performance and theoretical results. BART is a model which consists of a combination of regression trees. A \textit{regression tree} is any step function $g_{\mathcal{T}, \beta}:\mathcal{X}\rightarrow\mathbb{R}$ with $K$ leaves/partitions: $g_{\mathcal{T}, \beta}(x)\coloneqq  \sum_{k=1}^K \beta_{k}\mathbbm{1}_{\chi_k}(x)$. Here, we denote by $\mathcal{T}\coloneqq\{ \chi_k\}_{k=1}^K$, where $\chi_k \subset \mathcal{X}$, the partition of the domain, and by ${\beta}\coloneqq(\beta_1,\ldots, \beta_K)^\top \in \mathbb{R}^K$ the leaf values. The function $\mathbbm{1}_{\chi_k}:\mathcal{X}\rightarrow \mathbb{R}$ is an indicator function taking value $1$ whenever $x \in \chi_k$ and $0$ otherwise.

Similarly to decision trees and random forest, a single tree comprising of a large number of leafs might overfit to the training data. As a result, it is common to use an ensemble of shallow trees. We will call $T$-\textit{additive regression tree} any function which takes the form of a sum of regression trees: $g_{\mathcal{E}, \mathcal{B}}(x)\coloneqq \sum_{t=1}^T g_{\mathcal{T}_t, {\beta}_t}(x),$ where $\mathcal{E}\coloneqq\{\mathcal{T}_t \}_{t=1}^T$ and $\mathcal{B}\coloneqq \{\beta_t\}_{t=1}^T$. Finally, we call \textit{Bayesian additive regression tree} (BART) any distribution on the family of $T$-additive regression trees. Such a distribution can be constructed by specifying a (prior) distribution on the
partition $\mathcal{E}$ and leaf values $\mathcal{B}$. BART is hence a stochastic process whose sample space $\Omega$ consists of the product space of $K_t$-partitions of $\mathcal{X}$ and $\mathbb{R}^{K_t}$. For simplicity, BART is usually restricted to an approximation domain $\mathcal{X} = [0,1]^d$, but this is not a requirement in full generality.

Denote by $p$ the density of this prior distribution. We will follow the majority of the BART literature \cite{CART, BART, rockova2019theory} and use prior models which factorise in the following way: $p(\mathcal{E}, \mathcal{B}) \coloneqq \prod_{t=1}^T p(\mathcal{T}_t)p(\beta_t|\mathcal{T}_t)$ and $ p(\beta_t|\mathcal{T}_t) \coloneqq \prod_{k=1}^{K_t} N(\beta_{t,k}|0, 1/(16T))$, where $T$ is the number of trees and $K_t$ is the number of leaves in the $t^{\text{th}}$ tree. The construction of the distribution on partitions is usually itself done via a tree generating stochastic process; see further details in \cite{BART} (and a full description in Appendix \ref{appendix:BART}).

Given this prior $\mathbb{P}$ on $T$-additive trees, we can condition on $(X,y)$ to obtain a posterior $\mathbb{P}_n$ (with density $p_n$). We focus on regression with i.i.d.\@ Gaussian noise, but several generalisations exist \cite{pratola2017}. The corresponding posterior distribution on parameters will then imply a posterior distribution on $T$-additive trees. BART is hence another Bayesian model which we can use to approximate the integrand, using for example the posterior mean: $g^n(x) = \mathbb{E}_{\mathbb{P}_n}[f(x)] = \mathbb{E}[f(x)|X,y]  =  \int_{\Omega} g_{\mathcal{E},\mathcal{B}}(x) p_n(\mathcal{E},\mathcal{B}) d\mathcal{E} d \mathcal{B}$.

Unfortunately, this is not available in closed form and needs to be approximated. This is usually done using Markov chain Monte Carlo (MCMC) methods \cite{BART,Pratola2016} (although optimisation methods can also be used \cite{he2019xbart}). At each iteration $j$, the MCMC algorithm produces $T$ regression trees. More precisely, for the $t^{\text{th}}$ regression tree, the algorithm returns a set of leaf values $\textstyle \{\beta_{t,k}^j\}_{k=1}^{K_{t,j}} $ and a partition of the domain $\textstyle \{\chi_{t,k}^j\}_{k=1}^{K_{t,j}}$ where $K_{t,j}$ is the number of leaves for tree $t$ at iteration $j$. This gives a T-additive regression tree which we will denote $g^n_j$. After $m$ iterations, the MCMC mean is given by:
\begin{talign}
    \hat{g}^n(x) & = \frac{1}{m}  \sum_{j=1}^m g^n_j(x)  = \frac{1}{m}  \sum_{j=1}^m \sum_{t=1}^T \sum_{k=1}^{K_{t,j}} \beta_{t,k}^j \mathbbm{1}_{\chi^j_{t, k}}(x). \label{eq:MCMCapprox_f}
\end{talign}
A similar expression can be obtained for the variance of our posterior on $f$ at $x \in \mathcal{X}$. However, when compared to GPs, the entire BART posterior can in fact provide us with a more complex, possibly multimodal, distribution on $f$.

\section{Numerical Integration with Bayesian Additive Regression Trees}
\label{BART}

The Bayesian approach to numerical integration can be succinctly summarised in the following steps:
\begin{enumerate}
    \item \emph{Posit a Bayesian prior on the integrand $f$}. This consists of a stochastic process $g:\mathcal{X}\times \Omega \rightarrow \mathbb{R}$ together with a prior measure $\mathbb{P}$ on $\Omega$.
    \item \emph{Condition this prior on the data $(X,y)$ to get a posterior on $f$}. This is given in terms of the posterior measure $\mathbb{P}_n$, and the posterior mean is $g^n := \mathbb{E}_{\mathbb{P}_n}[ g] = \mathbb{E}[g|X,y]$.
    \item \emph{Obtain the posterior distribution on $\Pi[f]$}. This is given by the marginal distribution obtained by integrating out $\Pi$, and a point estimate for the posterior mean $\Pi[f]$, given by $\hat{\Pi}_{\text{BPNI}}[f] = \Pi[g^n]$.
\end{enumerate}
We remark that the posterior distribution on $\Pi[f]$, or even just the estimate $\hat{\Pi}_{\text{BPNI}}[f]$, may not be available in closed-form. In that case, MCMC estimates can be used (see the following section). In the case of GP-BQ with a prior $\mathcal{GP}(\mu,k)$, the posterior on $\Pi[f]$ is available in closed-form \cite{Briol2019PI} and is given by a Gaussian with mean $\mathbb{E}[\Pi [f] |X,y] = \Pi[k_{\cdot,X}] (k_{X,X}+\sigma^2I)^{-1}(y-\mu_X)$ and variance $\mathbb{V}[\Pi [f] | X,y]  = \Pi\bar{\Pi}[k]- \Pi[k_{\cdot,X}](k_{X,X}+\sigma^2 I)^{-1}\Pi[k_{X,\cdot}],$ with $\Pi[k_{\cdot,X}] = (\Pi[k(\cdot,x_1)],\ldots,\Pi[k(\cdot,x_n)])$, $\Pi[k_{X,\cdot}] = \Pi[k_{\cdot,X}]^\top$ and $\bar{\Pi}$ denotes the integration functional with respect to the second input. The posterior mean can be expressed as a quadrature rule whose weights depend on $X$, $k$ and $\sigma$. This expression will be available in closed form whenever $\Pi[k(\cdot,x)]$, known as the \textit{kernel mean}, is available in closed form. This is potentially a restrictive condition, but these integrals can themselves be approximated using quadrature, or the kernels could be constructed with this condition in mind  \cite{Oates2017,Oates2016CF2}.

A significant challenge for GP-BQ is the computational cost, which is $\mathcal{O}(n^3)$ due to the need to invert an $n \times n$ matrix. This can possibly be alleviated using fast GP algorithms, but often at the cost of introducing approximations to the posterior. The GP-BQ method is hence more expensive than most standard MC methods, but this should be balanced with significantly faster convergence rates; see \cite{Briol2019PI,Kanagawa2019,Wynne2020}. Overall, this method should hence be preferred to MC methods whenever the number of quadrature points $n$ is small or moderate, or when the integrand is itself expensive.

We now derive a new BPNI method via BART. Suppose that $g$, $\mathbb{P}_0$ are a BART prior (as described in Sec.~\ref{bart_functino}) and we have data $(X,y)$; the posterior on $f$ implies a posterior on $\Pi[f]$. An important distinction with GP-BQ is that the posterior on $\Pi[f]$ is not available in closed-form, and must be approximated. It is also not fully specified by the mean and variance, and may, in fact, be multi-modal. However, the mean of this posterior can be useful if a point estimate for $\Pi[f]$ is required, and the variance can be a useful summary of uncertainty. Expressions for both of these quantities are included in the proposition below for completeness.
\vspace{2mm}
\begin{proposition}[\textbf{MCMC for BART-Int}]\label{eq:BARTBQ_mean&var}
The MCMC approximation of $\Pi[f]$ consists of samples:
\begin{talign*}
s_j = \Pi\left[g^n_j\right] = \sum_{t=1}^T\sum_{k=1}^{K_{t,j}} \beta^j_{t, k}\Pi\left[\mathbbm{1}_{\chi^j_{t, k}}\right], \qquad j\in\{1,\ldots,m\}.
\end{talign*}
These lead to estimates of the mean $ \textstyle \widehat{\mathbb{E}}_{\text{m}}[\Pi[f]|X,y] = \Pi[\hat{g}^n] =\frac{1}{m}\sum_{j=1}^m s_j$ and variance $\textstyle \widehat{\mathbb{V}}_{\text{m}}[\Pi[f]|X,y]  = \frac{1}{(m-1)} \sum_{j=1}^m (s_j - \Pi[\hat{g}^n])^2$.
\label{prop:mcmc}
\end{proposition}
The posterior samples $\{s_j\}_{j=1}^m$ can only be obtained in a closed form whenever probabilities of the form $\Pi(\chi)$ can be computed for any $\chi$ in the partition. This issue corresponds to the well-known issue of intractable kernel means for GP-BQ. The simplest case for which $\Pi(\chi)$ can be computed is when $\Pi$ is the uniform distribution. If $\pi$ is available in closed form, one option is to model $f \pi$ as the integrand, in which case the integral is once again against a uniform. However, this makes the specification of a BART prior challenging since very little will usually be known about $f \pi$, and as a result the uncertainty quantification may be unreliable.

Another solution, used in this paper and commonly used for intractable kernel means, is to approximate these probabilities using samples $\{\tilde{x}_i\}_{i=1}^l$ representative of $\Pi$. This leads to approximate MCMC samples and an estimate of $\Pi[f]$ of the form:
\begin{talign*}
&\hat{s}^l_j=\frac{1}{l}\sum_{i=1}^{l} g^n_j(\tilde{x}_i)= \frac{1}{l}\sum_{i=1}^{l} \sum_{t=1}^T\sum_{k=1}^{K_{t,j}} \beta^j_{t, k}\mathbbm{1}_{\chi^j_{t, k}}(\tilde{x}_i), \\
&\widehat{\mathbb{E}}^l_{\text{m}}[\Pi[f]|X,y] = \frac{1}{m}\sum_{j=1}^m \hat{s}_j^l.  \qquad  \widehat{\mathbb{V}}^l_{\text{m}}[\Pi[f]|X,y]  = \frac{1}{(m-1)} \sum_{j=1}^m \left(\hat{s}^l_j - \left( \frac{1}{m}\sum_{k=1}^m \hat{s}_k^l\right)\right)^2.
\end{talign*}
Using this approach might be counter-intuitive since it means that some of the uncertainty is not quantified in a Bayesian manner. However, when the integrand $f$ is expensive but a very large number of data points are available (i.e. $n \ll l$), the additional error can be made to be negligible relative to the overall integration error. These settings are common in practice, as will be shown in Sec. \ref{Performance}.

A significant advantage of BART-Int over GP-BQ is the computational cost. For BART, the cost is $\mathcal{O}(T m n)$ (or $\mathcal{O}(T m (n+l))$), although the constant depends on the properties of the trees (see \cite{Pratola2014} for a detailed analysis). When $n$ is large, this can be much cheaper than the $\mathcal{O}(n^3)$ for GP-BQ.
We note that $\{x_i\}_{i=1}^n$ can be selected through any quadrature rule. However, when $f$ is computationally expensive, it may be preferable to adaptively select $\{x_i\}_{i=1}^n$ using tools from experimental design and active learning \cite{Santner2018}. Recent advances in this direction in numerical integration focus on minimising some distance between $\Pi$ and $\{x_i\}_{i=1}^n$ (see e.g.\@ kernel herding \cite{Chen2010}, Stein Variational Gradient Descent \cite{Liu2016SVGD} and Stein points \cite{Chen2018}). All of these methods could be combined with our BART-Int estimators, but such an approach may be sub-optimal in the sense that our objective is not to approximate $\Pi$, but only $\Pi[f]$. We propose instead a method which focuses on improving the fit of the BART posterior mean. Our approach is hence closer to previous active learning strategies for GP-BQ, including the work of \cite{Gunter2014,briol2015frank}. It is summarised in Algorithm~\ref{alg:SQ}.

\begin{algorithm}[h!]
  \caption{Sequential Design for BART-Int}
  \label{alg:SQ}
\begin{algorithmic}[1]
  \State {\bfseries Inputs:} Initial point(s) $X$ (set of size $n_{\text{ini}}$), response(s) $y$, number of candidate design points $S$, total number of points $n$, acquisition criterion $J$ and number of MCMC samples $m$.
  \vspace{1mm}
  \For{$i \in \{ n_{\text{ini}},\ldots,n-1\}$ iterations}
  \State Obtain $m$ posterior samples (given $X^{i},y^{i}$) \& pick a candidate set $\mathcal{C} = \{c_1, \ldots, c_S\} \subset \mathcal{X}$.
  \State Find $\,c^* = \mbox{argmax}_{c \in \mathcal{C} } J(c)$ \& set $X^{i+1} \leftarrow X^i \cup \{c^*\}$, $\; y^{i+1}\leftarrow y^i\cup \{ y_c^{*}\}, \; y_{c}^*
  = f(c^*)+\epsilon^*$.
  \EndFor
\end{algorithmic}
\end{algorithm}
The second step consists of comparing the suitability of each point in $\mathcal{C}$ according to some criterion $J:\mathcal{X}\rightarrow \mathbb{R}$, then adding the best point to our design set. At iteration $i$, we propose to use
\begin{talign}
    J(c) =  \widehat{\mathbb{V}}_{m}[f(c)\pi(c)|X^i,y^i],
\end{talign}
which leads to selecting points where the uncertainty in $f$ is highest. Alternatively, one could choose $c$ which minimises the posterior variance: $J(c)=- \widehat{\mathbb{V}}_{m}[\Pi[f]|X^i \cup \{c\},y^i \cup \{y_c\}]$. However, this would require training a different MCMC sampler for each point in $\mathcal{C}$, which would significantly increase the computational cost.

\section{Theoretical Results}

We now introduce novel concentration results. All proofs are given in Appendix \ref{appendix:proofs}. These results hold for BART-Int, but are also of independent interest since they can lead to rates for other BPNI methods. The results will hold for an integrand in a normed subspace $\mathcal{H}$ of $L^2(\Pi):=\{h:\mathcal{X} \rightarrow \mathbb{R} \text{ s.t. } \|h\|_{L^2(\mathcal{X})}:= \int_{\mathcal{X}} h^2(x) \pi(x) dx < \infty\}$. They will depend on the contraction rate of the posterior onto the integrand, as measured in some empirical norm: $\textstyle \|f\|_n = (\frac{1}{n} \sum_{i=1}^n f(x_i)^2)^{1/2}$. Below, we will use $(X^n,y^n)$ to denote a dataset $(X,y)$ of size $n$.

\vspace{2mm}
\begin{theorem}[\textbf{Concentration Bound for BPNI}]\label{thm:convergence_rate}
Suppose we have a normed space $\mathcal{H} \subseteq L^2(\Pi)$ such that the integrand $f \in \mathcal{H}$. Furthermore, suppose that our BPNI posterior is based on a stochastic process $g$ (jointly measurable over $\mathcal{X} \times \Omega$) and a sequence of data $\{(X^n,y^n)\}_{n \in \mathbb{N}}$ such that $\forall \omega \in \Omega, g(\cdot,\omega) \in \mathcal{H}$ and $\exists N \in \mathbb{N}_+$ such that:
\begin{itemize}
\item [A1.]  $\exists \{\varepsilon_n\}_{n \geq N}$ such that $\lim_{n \rightarrow \infty} \mathbb{P}_n[\|f - g\|_n > A_n \varepsilon_n]=0$ for any $A_n \rightarrow \infty$ as $n \rightarrow \infty$.

\item [A2.] $\exists \{\gamma_n\}_{n \geq N}$ with $\gamma_n \rightarrow 0$ as $n \rightarrow \infty$ such that $\sup_{\|h\|_{\mathcal{H}} \leq 1} |\frac{1}{n}\sum_{i=1}^n h(x_i) - \Pi[h]| = O(\gamma_n)$.
\end{itemize}
Then, we have $\lim_{n\rightarrow \infty} \mathbb{P}_n[| \Pi[f] -  \Pi[g]| > C_n \max(\varepsilon_n,\gamma_n)] = 0$ for any $C_n \rightarrow \infty$ as $n \rightarrow \infty$.
\end{theorem}
We now make a number of remarks on the assumptions and statement of the theorem:
\begin{enumerate}
    \item To guarantee concentration of a BPNI posterior, we require two ingredients: (i) the Bayesian posterior should concentrate on $f$ as $n$ grows, and (ii) the sequence of point sets $X^n$ should give a sequence of quadrature rules which can integrate both $f$ and $g^n$ for any $n$ large enough. The rates at which (i) and (ii) occur then control the overall concentration rate of the BPNI posterior.
    \item The assumptions in this theorem are stated in a way such that the order in $n$ at which concentration of the BPNI posterior occurs is explicit. In particular, A1 is a standard assumption in Bayesian nonparametric and guarantees concentration of the posterior on $f$ at rate controlled by $O(\varepsilon_n)$. Similarly, A2 is a common assumption for quadrature rules and guarantees a quadrature rate of $O(\gamma_n)$. The resulting concentration rate of the BPNI posterior is then controlled by $O(\max(\varepsilon_n,\gamma_n))$, the slowest of these two rates.
    \item A number of quadrature rules can integrate functions in $L^2(\Pi)$ at a rate $O(n^{-1/2})$ (i.e. $\gamma_n =n^{-1/2}$), but have faster worst-case integration rates for subspaces of $L^2(\Pi)$ which are sufficiently regular. For this reason, we introduced the subspace $\mathcal{H}$ and the conditions that $f \in \mathcal{H}$ and $ \forall \omega \in \Omega$, $ g(\cdot,\omega) \in \mathcal{H}$, in order to allow for faster rates for $\gamma_n$.
    \item Results of the form in A1 are available for a variety of Bayesian nonparametric models \cite{Ghosal2017}, including GPs \cite{Vaart2011} and Bayesian neural networks \cite{Polson2018}. This is a rapidly evolving field and most results were only derived in recent years. The advantage of Thm. \ref{thm:convergence_rate} is that any contraction result derived in the future can be plugged in to understand the corresponding BPNI method.

    \item Assumption A2 will hold for most reasonable choices of design points. For example, it would hold for $V$-uniformly ergodic MCMC when $\mathcal{H}=\{h \in L^2(\Pi): \sup_{x \in \mathcal{X}} |h(x)|/V(x) < \infty \}$ (under regularity conditions given in \cite{Roberts2004}). Note that Assumption A1 may also depend on the design points, and a discussion of various cases can be found in \cite{jeong_art_2020}.

    \item As for BART-Int, we might sometimes need to approximate the BPNI using samples $\{\tilde{x}_i\}_{i=1}^l$. In that case, it is sufficient for $\{\tilde{x}_i\}_{i=1}^l$ to satisfy an assumption of the form of A2 for Thm. \ref{thm:convergence_rate} to hold. In general, it will be possible to take $l$ much larger than $n$, so we could potentially even allow for slower quadrature rates than $O(\gamma_n)$.
\end{enumerate}

Now that we have discussed our result on concentration of BPNI posteriors, we consider implications for BPNI with tree-based models. Most results are presented for the space of $\alpha$-H\"older continuous functions, denoted $\mathcal{H}^{\alpha}$, where $ \|f\|_{\mathcal{H}^\alpha} \coloneqq \sup_{x,y \in \mathcal{X}} |f(x)-f(y)|/\|x-y\|_2^{\alpha} < \infty$ (when $\alpha=1$, we recover Lipschitz functions). Several results are of interest:

\begin{itemize}
    \item \cite{rockova2017posterior,rockova2019theory} showed that Bayesian CART and BART concentrate on targets which are sums of functions in $\mathcal{H}^{\alpha}$ ($0 < \alpha \leq 1)$ that are constant in $d-d_0$ coordinates at a rate $\varepsilon_n = n^{-\frac{\alpha}{2\alpha+d_0}}\log^{\frac{1}{2}}n$. This does not require a-priori knowledge of the coordinates in which the function is constant.

    \item \cite{linero2018bayesian} showed that the soft-BART posterior raised to a fractional power concentrates at the rate $\varepsilon_n = n^{-\frac{\alpha}{2\alpha+d_0}}\log^\beta n+\sqrt{n^{-1}d_0 \log d}$ for any $\beta \geq \alpha(d_0+1)/(2\alpha+d_0)$, even when $\alpha>1$.
\end{itemize}
These results are not directly comparable with the most of the literature for GP-BQ \cite{Briol2019PI,Kanagawa2017,Kanagawa2019} since these consider the interpolation setting (where $y_i = f(x_i)$). The closest result is in \cite{Wynne2020}, which presents a result for regression with i.i.d. Gaussian noise and convergence in expectation. The rate in \cite{Wynne2020} is identical to those above except for for sparse high-dimensional functions where it is slower.

Recent Bernstein-von Mises results for BART presented in \cite{Rockova2019BvM} (Thm. 5.2) also provide guarantees on the uncertainty quantification provided by the BART-Int posterior distribution on $\Pi[f]$, and show that in the asymptotic limit of $n\rightarrow \infty$, the Bayesian credible intervals will coincide with frequentist confidence intervals. Similar results were obtained for sparse deep networks in \cite{Wang2020}.

Before concluding, we note that any MCMC approximation required to approximate the posterior will also impact the convergence. This is made more precise through the following result.

\vspace{2mm}
\begin{proposition}[\textbf{MCMC Approximation for BPNI}] \label{prop:MCMC_convergence}
Let $\hat{\mathbb{E}}_m[\Pi[g]]$ denote the MCMC approximation of $\mathbb{E}[\Pi[g]]$ where the expectation is with respect to some measure $\mathbb{P}$. Assume $\Pi[g] \in L^2(\mathbb{P})$, where the Markov chain targets $\mathbb{P}$ and is geometrically ergodic and reversible. Then, $\exists \sigma_{\text{MCMC}} >0$ such that
\begin{talign*}
\sqrt{m} | \mathbb{E}[\Pi[g]]- \hat{\mathbb{E}}_m[\Pi[g]]| \rightarrow \mathcal{N}(0,\sigma^2_{\text{MCMC}})
\end{talign*}
in distribution as $m \rightarrow \infty$.
\end{proposition}
Geometric ergodicity is a well-studied concept in MCMC theory; see Appendix \ref{sec:MCMC_proof}. Combining Thm. \ref{thm:convergence_rate} and Prop. \ref{prop:MCMC_convergence} gives us some intuition as to how to balance the computational cost of increasing $n$ and $m$ in order to obtain the smallest possible approximation error. Unfortunately, to the best of our knowledge, the ergodicity of MCMC samplers commonly used for BART has not yet been studied.

\section{Numerical Experiments}
\label{Performance}

We now illustrate BART-Int and GP-BQ on a range of synthetic problems and a Bayesian survey design problem. We emphasise that our goal is to compare methods which can provide a Bayesian quantification of uncertainty; as a result, it is very much possible that non-Bayesian methods could have better point estimate performance on some of these problems. For BART-Int, we used the default prior settings in \verb!dbarts! \cite{dorie2018dbarts}, whereas for GP-BQ we used a Mat\'ern kernel whose lengthscale was chosen through maximum likelihood. Further details can be found in Appendix \ref{appendix:experiments}.

\begin{table}[h!]
\centering
 \label{table:genz}
 \begin{minipage}[t]{\linewidth}
 \begin{adjustbox}{max width=\textwidth}
 \begin{tabular}[t]{| c | c | | c | c | c | c | c | c | c |}
    \hline
     Setup & Method  & Cont   & Copeak &Disc & Gaussian & Oscil & Prpeak & Step \\
     \hline
     \hline
    $d=1$  & BART-Int  & 1.21 &  \textbf{7.32e-01} & \textbf{6.06e-02} & \textbf{6.80e-01} & \textbf{228} & 2.58& \textbf{7.85e-03} \\
    $n_\text{ini}=20$ & MI  &
\textbf{ 7.11e-01} &
8.99e-01 &
1.55e-01 &
 8.41e-01 &
 251 &
 1.28 &
5.90e-02 \\
     $n_\text{seq}=20$& GP-BQ  &9.38e-01 &  3.41 & 1.61e-01 & 7.42e-01 & 229& \textbf{9.88e-01}& 2.88e-02\\
     \hline
      $d=10$ & BART-Int  & \textbf{9.10e-04} & 43.8  &\textbf{1.31e-02}
      &3.37e-03 &1.19e-02 &1.77e-02 &\textbf{4.33e-03} \\
     $n_{\text{ini}}=200$ & MI  &
     2.78e-03  &
     \textbf{2.12e-01} &
     1.33e-01&
     1.18e-02 &
     8.43e-02 &
     3.70e-03 &
     2.46e-02 \\
     $n_\text{seq}=200$&  GP-BQ   &
     1.52e-03  &
     2.79e-01 &
     1.84e-01&
     \textbf{3.02e-03} &
     \textbf{2.39e-03} &
     \textbf{6.59e-04} &
     8.36e-03 \\
     \hline
    \end{tabular}
    \end{adjustbox}
\caption{Integration of the six Genz test functions and a step function. The values correspond to the mean absolute percentage error (MAPE) over 20 separate runs.}
    \end{minipage}
    \vspace{-6mm}
\end{table}

\begin{wrapfigure}{r}{0.3\linewidth}
    \vspace{-7mm}
    \includegraphics[width=\linewidth,trim=0.5cm 0 0.5cm 0.2cm, clip]{combined_7.pdf}
    \vspace{-5mm}
    \caption{Integration of $f_\text{step}$ against a uniform, $d=1$, $n_{\text{ini}}=20$ and $n_{\text{seq}}=20$. }
    \label{figure:step}
    \vspace{-6mm}
\end{wrapfigure}
\paragraph{Genz Functions} We start with a standard benchmark set of multi-dimensional integrands proposed by Genz \cite{Genz} (see Appendix \ref{appendix:experiments}). This consists of six functions in $\mathcal{X}=[0,1]^d$, where $d$ can be varied to adjust the level of difficulty of the problem. We also add a step function:  $f_\text{step}(x) = \mathbbm{1}_{ \{x_1 \in (0.5, 1] \} }(x)$, where $x_1$ is the first component of $x$. We use $n_{\text{ini}}=20d$ design points, with $n_\text{seq}=20d$ additional points selected by the sequential scheme in Algorithm~\ref{alg:SQ}.

Table~1 shows the mean absolute percentage error (MAPE) of BART-Int ($m=1500$, $T=200$ $m=1000$, $T=50$, with a burn-in of 1000 and keeping every 5 samples afterwards), GP-BQ and Monte Carlo Integration (MI) for all six Genz functions in $d=1$ and $10$. The MAPE is given by given by $\frac{1}{r}\sum_{t=1}^r|\Pi[f]-\Hat{\Pi}_t[ f]|/ |\Pi[f]|$, where $\Hat{\Pi}_t[f]$ for $t=1,\ldots,r$, are estimates of $\Pi[f]$ for $r$ different initial i.i.d. uniform point sets.

BART-Int outperforms GP-BQ when $f$ is not continuous, e.g.\@ the Step functions and the Discontinuous functions. In those cases, the posterior distribution of BART (red dots) is also more concentrated around the truth than GP-BQ (whose 95\% credible intervals are in shaded blue); see Fig \ref{figure:step}. However, GP-BQ is strongest when estimating integrals of smooth continuous in $d=1$ (see all other integrands). This is to be expected since these functions are smooth enough to be well modelled by the GP (see  convergence rates in \cite{Wynne2020}), while the step-function nature of BART makes it more appropriate for non-smooth functions. Finally, in $d=10$, BART-Int outperforms GP-BQ for all integrands as it is adaptive to the important features of the integrand, with the largest gains being once again for discontinuous functions.

\vspace{-3mm}

\begin{wrapfigure}{r}{0.35\linewidth}
 \begin{center}
 \vspace{-8mm}
    \includegraphics[width=0.8\linewidth,trim=0.6cm 0.7cm 0.7cm 0, clip]{computational_complexity.pdf}
    \vspace{-1mm}
    \caption{Run-times of GP-BQ \& BART-Int for $\Pi[f_\text{step}]$ (w.o. seq. design). }
    \label{fig:computational_complexity}
    \end{center}
     \vspace{-8mm}
\end{wrapfigure}
\paragraph{Computational Complexity}
The computational complexity of BART-Int is $\mathcal{O}(Tmn)$, which is much cheaper than the $\mathcal{O}(n^3)$ for GP-BQ with even moderately large $n$. This behaviour can be seen empirically in Figure~\ref{fig:computational_complexity}, noting $T$ and $m$ have been fixed. The computational time of BART-Int is based on the assumption that the tree-based operations are constant, which is reasonable so long as the tree sizes are moderate \cite{Pratola2014}. Furthermore, it has been shown empirically that the run time of BART is almost independent of the dimension, $d$, of the data, see e.g. \cite{Pratola2014} and section~6 of \cite{BART}. However, with larger $d$, it would require a longer burn-in period for the parameter space to be fully explored, and thus the run-time for BART (thus also BART-Int) would increase for the purpose of approximation accuracy.

\paragraph{Rare Events Simulation}
Rare-event simulation is an important tool in risk management. It is used in a range of applications, such as for safety testing autonomous driving systems \cite{o2018scalable}, queuing systems \cite{rubino2009rare} and finance \cite{chan2012improved}. Traditionally, rare-event simulation is tackled using Monte Carlo methods \cite{rubino2009rare}, but we propose to use BART-Int instead.
\begin{wrapfigure}{r}{0.4\linewidth}
\begin{minipage}[t]{\linewidth}
\centering
\begin{adjustbox}{max width=0.99\textwidth}
\begin{tabular}[t]{| c | c | c | c |}
\hline
  & Method  & MAPE   & Std. Err. \\
 \hline
  & BART-Int  & 1.71e-01 & 2.56e-02 \\
 $d=5$  & MI  &1.95e-01 & 2.29e-02 \\
  $n=2500$&  GP-BQ  & \textbf{1.68e-01}& 2.09e-02 \\
 \hline
  & BART-Int  &\textbf{1.56e-02} &2.35e-03  \\
 $d=10$  & MI  & 9.98e-01  & 4.47e-04 \\
 $n=5000$ &  GP-BQ   & 2.72e-02& 5.20e-03 \\
 \hline
  & BART-Int & \textbf{8.40e-03}  &1.60e-03 \\
 $d=20$    & MI  &9.94e-01 & 6.34e-04\\
 $n=10000$ &  GP-BQ    &2.92e-02 & 4.90e-03 \\
 \hline
\end{tabular}
\end{adjustbox}
\captionof{table}{Performance on the portfolio loss event over 20 runs.}
\label{tab:high-dimensionality}
\vspace{-5mm}
\end{minipage}
\end{wrapfigure}
Let $\ell:\mathbb{R}^d\rightarrow\mathbb{R}$ be a measurable function representing our loss. The central quantity of interest here is the probability of obtaining a loss larger than some constant $\gamma \in \mathbb{R}$, which can be written as $p_\gamma = \int_{\mathcal{X}} \mathbbm{1}_{\{\ell(x)>\gamma\}} (x)\Pi(dx) $.

We consider a problem of high-dimensional portfolio management adapted from \cite{chan2012improved}. Suppose we have $d$ loans to obligors, each with value $c_i$ for $i=1,\ldots,d$. Let $x_i$ denote the financial strain on loan $i$, and suppose $d_i$ is a thresholds after which default occurs. We assume that the distribution of financial strains is $\text{Exp}(1)$. We can define the portfolio loss as $\ell(x) = \sum_{i=1}^d c_i \mathbbm{1}_{\{x_i > d_i\}} (x)$. Given a threshold $\gamma$, we can compute the probability of making a loss greater than $\gamma$ as $p_\gamma$.

We set $\gamma=2$, $d_i = 0.5i$ and $c_i=0.2i$. We use the same experimental settings as for the Genz functions, but do not use a sequential design and instead set $n=500d$. The number of post-burn-in samples is chosen to be $10^4$. We can see that for $d=5$ most of the algorithms perform similarly, but as $d$ increases BART-Int gradually performs better, especially when compared to GP-BQ.

\begin{wraptable}{r}{0.3\linewidth}
    \centering
    \vspace{-10mm}
    \includegraphics[width=0.99\linewidth,trim=0.5cm 0.5cm 1cm 0, clip]{resultsBin2.pdf}
    \captionof{figure}{Proportion of individuals with log-salary >$10$.}
    \label{fig:survey}
    \vspace{-8mm}
\end{wraptable}
\paragraph{Bayesian Survey Design} In surveys, response variables are collected from a subset of a population based on answers to a set of questions. Survey design concerns sampling strategies to obtain population-representative estimates from these samples. The standard approach is simple or stratified random (Monte Carlo) sampling. Bayesian hierarchical models are also often used in this setting to analyse survey data in order to stabilise estimates and to borrow strength when making sub-population estimates for underrepresented subgroups or locations \cite{gelman2006data}.

We propose a new approach we term ``Bayesian Survey Design'', using BART-Int to adaptively choose the next individual to survey.
Assume we have access to a small set of $n_{\text{ini}}$ survey responses with a mixture of continuous and categorical covariates (educational attainment, age, etc).
In addition, we assume there is a much larger set of $S$ individuals for whom demographic variables are known but the response variable is unknown but available for surveying. Such active surveying has previously also been explored in Bayesian active learning \cite{10.5555/3042573.3042683}. To demonstrate this approach, we use individual-level anonymised census data from the United States \cite{us_census_bureau_american_2018} and estimate the proportion of the population whose income is higher than around $\$$22,000, or log-income greater than $10$. This is equivalent to computing the integral of the following indicator function: $f(x) = \mathbbm{1}_{\{\text{log-income} >10 \}}(x)$ where the input consists of all other covariates in the survey and $\Pi$ is the distribution of these covariates in the population.

\begin{wraptable}{r}{0.4\linewidth}
\vspace{-5mm}
    \begin{adjustbox}{max width=0.9\textwidth}
    \begin{tabular}[t]{  c | c | c }
        \hline
          Method  & MAPE   & Std. Err. \\
         \hline
         BART-Int  & \textbf{4.66e-02} & 7.38e-03 \\
         MI  &5.68e-02 & 7.01e-03\\
          GP-BQ  & 9.73e-02& 1.44e-02 \\
         \hline
        \end{tabular}
    \end{adjustbox}
    \captionof{table}{Performance on the Survey Design problem over 20 runs.}
    \label{tab:survey_design}
    \vspace{-4mm}
\end{wraptable}

We consider a universe of $454,816$ possible respondents with $8$ demographic covariates. After a one-hot encoding of categorical variables, the dimensionality is $d=24$. We randomly select our initial set (of size $n_\text{ini}=20$) and candidate set (of size $S=10,000$). We compute ground truth using all $454,816$ observations and use Algorithm \ref{alg:SQ} to select $n_\text{seq}=200$ new individuals to survey. As seen in Table~\ref{tab:survey_design} and Figure~\ref{fig:survey}, BART-Int outperforms both MI and GP-BQ. As expected, the posterior distributions contract to the population mean as $n$ increases. Furthermore, the BART-Int posterior (red dots) is centered around the truth, whereas the 95\% credible intervals (blue shaded region) for GP-BQ seem overconfident at the wrong value around $n=100$.

\vspace{-1mm}

\section{Related Work}

Before concluding, we briefly comment on the connections between BART-Int and other tree-based algorithms for integration. Our work is most closely related to GP-BQ, but with tree-based priors. Although trees had not been considered in this literature, the nearest-neighbors method in \cite{llorente2020adaptive} uses GP kernels which depends on indicator functions at Voronoi partitions, and yield approximations to $\Pi[f]$ which closely resembles that in our paper.

We note that tree-based methods have received significant interest in the literature in numerical integration, but previous work does not come with a Bayesian interpretation. For example, the \texttt{VEGAS} algorithm \cite{PeterLepage1978} from the physics literature can be thought of as a tree-based method to reduce the variance of MI. More recent work also includes MCMC and SMC samplers with proposals estimated using trees \cite{felip2019tree, hanson2011polya}. Finally, \cite{foster2020model} also uses regression trees as surrogates, but again in a non-Bayesian framework.

\vspace{-2mm}

\section{Discussion}
\label{Discussion}

\vspace{-2mm}
We proposed a novel BPNI algorithm which uses BART instead of a GP to model the integrand. BART has several advantages in settings where GPs do not perform well. It is easy to tune, robust to overfitting and enjoys a low computational complexity. It is also a natural model for discontinuous functions and can automatically adapt to sparse high-dimensional functions. We demonstrated these advantages through theoretical results, including contraction of the posterior distribution onto the value of the integral, and through a set of benchmark tests and a survey design experiment. However, it is also important to note that BART-Int did not perform as well as GP-BQ in certain low dimensional and smooth settings. We therefore see BART-Int as a useful addition to the toolset of BPNI which complements, rather than replaces, GP-BQ.

There are a number of potential directions for future research. From an applications-viewpoint, there are interesting parallels to be made with algorithms used to estimate treatment effects in causal inference, and we foresee some possible use for BART-Int in this field. This was discussed in more details in an opinion piece published in \cite{hahn2020}. In terms of theory, little is known about the ergodicity of the MCMC sampler used to sample posterior trees. Further work in this area could help practitioners understand the impact of the MCMC approximation on the estimate of the integral for BART-Int. In terms of methodology, while the sequential design approach proposed in this paper uses the (estimated) posterior variance to select new query points, it would be interesting to explore how the use of other acquisition functions may improve the estimates. Similarly, extensions of BART, such as soft-BART \cite{linero2018bayesian} and heteroscedastic-BART \cite{pratola2017}, should also be of practical interest, as they may provide additional advantages over GP-based methods. Finally, some tree-based models can naturally handle categorical variables without the need for one-hot encoding; this could significantly improve the performance for problems such as the Bayesian survey design experiment in this paper.

\section*{Broader Impact}
Our paper is about numerical integration, a common computational problem in machine learning. It provides an approach to improve the accuracy of such approximations, as well as obtaining credible intervals representing our uncertainty about the value of the integral. The increased accuracy of the method may allow the users to reduce their computational requirements, which may, further down the line, have some impact on mitigating the impact of large computer clusters on climate change.

However, the broader impact of the method will mostly depend on the applications of the algorithms that it will be used to enhance: applications to ethical algorithms will have a positive impact, but application to unethical algorithms will also have a negative impact. For example, in the Bayesian survey design problem studied in the numerical experiments section, we have shown that there is potential for application in predicting population proportions resulting from binary outcome variables. The categorical variables for instance, however, may be unethically used as explanatory variables for certain problems, and this would be an issue depending on how the conclusions are drawn.

\begin{ack}
The authors would like to thank Matthew Fisher, Takuo Matsubara and Chris Oates for useful feedback on a previous draft of the paper. HZ, XL were supported by the EPSRC Centre for Doctoral Training in Modern Statistics and Statistical Machine Learning (EP/S023151/1) and the Department of Mathematics of Imperial College London. HZ was supported by Cervest Limited. XL was supported by the President’s PhD Scholarships of Imperial College London. RK and ZS were supported by the Undergraduate Research Opportunities Programme of Imperial College London. FXB was supported by the Lloyd's Register Foundation Programme on Data-Centric Engineering and the Alan Turing Institute under the EPSRC grant [EP/N510129/1], and through an Amazon Research Award on ``Transfer Learning for Numerical Integration in Expensive Machine Learning Systems''.
\end{ack}
\bibliographystyle{abbrv}

config: arxiv


TITLE: Write $0.2154154\overline{154}$ as a fraction
QUESTION [3 upvotes]: Let $x = 0.2154154\overline{154}$ , I have to prove that it is a rational number just by writing it as a fraction with the proper steps.
I note that the repeating part, $154$, is composed by 3 digits. Thus, using a trick that I have learnt, I can write an equation of this type:
$$1000x = ?$$
I am not sure about what goes on the right side of the equation because of this $0.2$
I mean, if, for instance, the number was $x = 0.154$ periodic, I would write the equation as:
$1000x = 154 + x$    and solve for $x$
I have tried various attempts, re-writing $0$, $2$ as $\frac{1}{5}$, but I don't get the number which equals our $x$.

REPLY [11 votes]: First Method
What you basically do in this method is find a multiple of $x$ that will get rid of the repeating part when you subtract it from x itself. (HINT: $10x$, $100x$, it all depends on the number of digits in the repeating part.)
$$x=0.2154154154...$$
$$10x=2.154154154...$$
$$10000x = 2154.154154...$$
$$10000x - 10x = 9990x = 2154 - 2 = 2152$$
$$x = \frac{2152}{9990}$$
Note we picked $10000x$ here because in the expression of $10x$, the repeating part has three digits.

Second Method
$$2.154154154 = 2 + 0.154 + 0.000154 + 0.000000154 + ...$$
$$= 2 + 0.154 + 0.154(0.001) + 0.154(0.001)^2 +  0.154(0.001)^3 + ...$$
$$2 + \frac{0.154}{(1-0.001)} = 2 + \frac{0.154}{0.999} = \frac{2.152}{0.999}$$
$$= \frac{2152}{999}$$
Thus $$2.\overline{154} = \frac{2152}{999}$$
Since the original number is $0.2\overline{154}$, the fractional form is:
$$ \frac{2152}{9990}$$


TITLE: Local mass function in spherically symmetric spacetime
QUESTION [0 upvotes]: I am studying the paper `Inflation and de Sitter Thermodynamics' at https://arxiv.org/abs/hep-th/0212327 . I have problems with the way they define a local mass function in a general spherically symmetric spacetime. I am quoting from Appendix B of the paper: "A general spherical spacetime is described by the metric
\begin{align}
ds^2 = g_{ab}dx^a dx^b + \rho^2 d\Omega^2 \qquad
\end{align}
where $g_{ab}$ is the metric on a two manifold with coordinates $(t,r)$ and $\rho$ is the physical radius of spherical slices."
Firstly, I do not understand the difference between $r$ and $\rho$.
Then, before eq. (35) they say, "In spherical symmetry, it is possible to define a local mass function $M(x^a)$ by
\begin{align}
1 - \frac{2GM}{\rho} = (\nabla \rho)^2~."
\end{align}
It is not at all obvious to me how the above relation comes about. Any insight into this will be very helpful.

REPLY [0 votes]: 1) You can choose whichever coordinates you like in GR, in order to aid physical intuition and/or to ease calculations, etc. In this case, the manifold is sliced into 2-spheres which is done by choosing $r=$ constant and $t=$ constant, to get the metric
$$ds^2 = \rho^2 d \Omega^2$$
which is exactly the metric for a 2-sphere of radius $\rho$. What they mean is that, in general, $r$ may be different from $\rho$, but it is $\rho$ that defines the sizes of the spherical slices.
2) That's a definition, as stated in the paper, so they're relating $\rho$ to a new function $M$
$$\frac{2GM}{\rho} = 2 \kappa \rho$$
$$\Rightarrow M = \frac{\kappa \rho^2}{G}$$


TITLE: Dimension in Whitney's theorem
QUESTION [9 upvotes]: There are some interesting examples of (classes of) manifolds for which I suspect that the Whitney theorem can be strenghtened. For example it is known that every (smooth) $n$-dimensional manifold can be embedded into $\mathbb{R}^{2n-1}$ for example if $M$ is orientable or is not of dimension of the form $2^k$. In particular I'm interested what is the minimal number $emb(n)$ with the property that each $n$-dimensional manifold can be embedded into $R^{emb(n)}$ if we assume additionally that $M$ is:
a) simultaneously orientable and not of dimension of the form $2^k$,
b) a Lie group, and
c) the boundary of some other manifold?

REPLY [8 votes]: The following paper of Elmer Rees addresses precisely your question for Lie groups:
Rees, Elmer
Some embeddings of Lie groups in Euclidean space.
Mathematika 18 (1971), 152–156.
The main result is that if a compact connected Lie group $G$ of rank $\ell$ admits a faithful representation on $\mathbb{R}^k$, then $G$ admits an embedding into Euclidean space whose codimension is $k-\ell$. This improves on Whitney's Theorem when $k-\ell<\operatorname{dim}G$, which is the case for many of the classical families of Lie groups.