\section{Introduction}

Black-box optimization of an expensive-to-evaluate function based on point queries is a ubiquitous problem in machine learning. Bayesian optimization (or Gaussian process optimization) refers to a class of methods using Gaussian processes (GPs), whose main idea is to place a prior over the unknown function and update the posterior according to point query results. Bayesian optimization has a wide range of applications including parameter tuning \cite{Sno12}, experimental design \cite{griffiths2020constrained}, and robotics \cite{Liz07}.

In the literature on Bayesian optimization, the problem of optimizing a real-valued black-box function is usually studied under two settings: (1) Bayesian setting: the target function is sampled from a known GP prior, and (2) non-Bayesian setting: the target function has a low norm in reproducing kernel Hilbert space (RKHS).

In this work, we consider a setting falling ``in between'' the white-box setting (where the full definition of the target function is known) and the black-box setting, namely, a grey-box setting in which the algorithms can leverage partial internal information of the target function beyond merely the final outputs or even slightly modify the target function  \cite{astudillo2021thinking}.  Existing grey-box optimization methods exploit internal information such as observations of intermediate outputs for composite functions \cite{astudillo2019bayesian} and lower fidelity but faster approximation of the final output for modifiable functions \cite{huang2006sequential}. Numerical experiments show that the grey-box methods significantly outperform standard black-box methods \cite{astudillo2021thinking}.

Though any real-valued network can be treated as a single black-box function and solved with classical Bayesian optimization methods, we explore the benefits offered by utilizing the network structure information and the exact intermediate results under the noise-free grey-box setting.

\subsection{Related Work}

Numerous works have proposed Bayesian optimization algorithms for optimizing a single real-valued black-box function under the RKHS setting. For the noisy setting, \cite{Sri09,Cho17,gupta2022regret} provide a typical cumulative regret ${\widetilde{O}}(\sqrt{T}\gamma_T)$, where $T$ is the time horizon, $\gamma_T$ is the maximum information gain associated to the underlying kernel, and the ${\widetilde{O}}(\cdot)$ notation hides the poly-logarithmic factors. Recently, \cite{Jamieson21,Salgia2020,li2022gaussian} achieved cumulative regret ${\widetilde{O}}(\sqrt{T\gamma_T})$, which nearly matches algorithm-independent lower bounds for the squared exponential and Mat\'ern kernels \cite{Sca17a,Cai20}. For the noise-free setting, \cite{bull2011convergence} achieves a nearly optimal simple regret $O(T^{-\nu/d})$ for the Mat\'ern kernel with smoothness $\nu$. This result implies a two-batch algorithm that uniformly selects $T^\frac{d}{\nu+d}$ points in the first batch and repeatedly picks the returned point in the second batch, has cumulative regret $O(T^\frac{d}{\nu+d})$. In addition, \cite{lyu2019efficient} provides deterministic cumulative regret $O(\sqrt{T\gamma_T})$, with the rough idea being to substitute zero noise into the analysis of GP-UCB \cite{Sri09}.

Recently, several Bayesian optimization algorithms for optimizing a composition of multiple functions under the noise-free setting have been proposed. \cite{nguyen2016cascade} provides a method for cascade Bayesian optimization; \cite{astudillo2019bayesian} studies optimizing a composition of a black-box function and a known cheap-to-evaluate function; and \cite{astudillo2021bayesian} studies optimizing a network of functions sampled from a GP prior under the grey-box setting. Both \cite{astudillo2019bayesian} and \cite{astudillo2021bayesian} prove the asymptotic consistency of their expected improvement sampling based methods.

The most related work to ours is \cite{kusakawa2021bayesian}, which introduces two confidence bound based algorithms along with their regret guarantees. A detailed comparison regarding problem setup and theoretical performance is provided in \cref{sec:comparison}.  In short, we significantly improve certain dependencies in their regret for a UCB-type approach, and we study two new directions -- non-adaptive sampling and algorithm-independent lower bounds -- that were not considered therein (see Section \ref{sec:contributions} for a summary).

For function networks with $m$ layers, our cumulative regret bounds are expressed in terms of
\begin{align*}
    \Sigma_T=\max_{i\in[m]}\max_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}}\sum_{t=1}^T \sig{i}_{t-1}({\mathbf{z}}_t),
\end{align*}
where $\mathcal{X}^{(i)}$ and $\sig{i}_{t-1}(\cdot)$ denote the domain and posterior standard deviation  of the $i$-th layer respectively. This is a term for general layer-composed networks associated to domains $\mathcal{X}^{(1)},\dots,\mathcal{X}^{(m)}$ and kernel $k$. When $m=1$, the term $\max_{{\mathbf{x}}_1,\dots,{\mathbf{x}}_T\in\mathcal{X}}\sum_{t=1}^T \sigma_{t-1}({\mathbf{x}}_t)$ often appears in the cumulative regret analysis of classic black-box optimization \cite{Sri09,lyu2019efficient,vakili2022open}. Explicit upper bounds on $\Sigma_T$ will be discussed in \cref{sec:ucb_disc}

\subsection{Contributions} \label{sec:contributions}

We study the problem of optimizing an $m$-layer function network in the noise-free grey-box setting, where the exact intermediate results are observable. We focus on three types of network structures: (1) chain: a cascade of scalar-valued functions, (2) multi-output chain: a cascade of vector-valued functions, and (3) feed-forward network: a fully connected feed-forward network of scalar-valued functions.  Then:

\begin{itemize} \itemsep0ex
	\item We propose a fully sequential upper confidence bound based algorithm GPN-UCB along with its upper bound on cumulative regret for each network structure.
	\item We introduce a non-adaptive sampling based method, and provide its theoretical upper bound on the simple regret for the Mat\'ern kernel.
	\item We provide algorithm-independent lower bounds on the simple regret and cumulative regret for an arbitrary algorithm optimizing any chain, multi-output chain, or feed-forward nework associated to the Mat\'ern kernel.  We show that the upper and lower bounds nearly match in broad cases of interest.
\end{itemize}

Let $d$ denote the dimension of the domain, $d_{\max}$ denote the maximum dimension among all the $m$ layers, and $D_{2,m}$ denote the product of dimensions from the second layer to the last layer. With $B>0$ restricting the magnitude and smoothness of each layer and $L>1$ restricting the slope of each layer, a partial summary of the proposed cumulative regret bounds for the Mat\'ern kernel with smoothness $\nu\ge 1$ when $d\ge\nu$ and $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$ for $c=\Theta(1)$ is displayed in \cref{tab:summary_cml}, showing that the upper bounds and the lower bound share the same $BL^{m-1}T^{1-\nu/d}$ dependence when $d=d_{\max}$.

We highlight that to our knowledge, we are the first to attain provably near-optimal scaling (in broad cases of interest), and doing so requires both improving the existing upper bounds (e.g., \cite{kusakawa2021bayesian}) and attaining novel lower bounds.  A full summary of our theoretical results is provided in \cref{sec:summary}.

\begin{table*}[t!]
	\centering
	\begin{tabular}{ll}
		\toprule
		Lower bound & $\Omega(B(cL)^{m-1}T^{1-\nu/d})$ \\
		\midrule
		Upper bound (chains) & $2^{m+1} BL^{m-1}T^{1-\nu/d}$ \\
		Upper bound (multi-output chains) & $5^m BL^{m-1}T^{1-\nu/d_{\max}}$ \\
		Upper bound (feed-forward networks) & $2^{m+1} \sqrt{D_{2,m}} BL^{m-1}T^{1-\nu/d_{\max}}$ \\
		\bottomrule
	\end{tabular}
	\caption{Summary of proposed cumulative regret bounds for the Mat\'ern kernel when $d\ge\nu\ge1$ and $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$.}
	\label{tab:summary_cml}
\end{table*}

\section{Problem Setup}\label{sec:setup}

\begin{figure*}[t!]
\centering
\includegraphics[width=0.8\textwidth]{figs/setup.pdf}
\caption{A function network $g$ of $m$ layers with input ${\mathbf{x}}$ and output $y$.}
\label{fig:setup}
\end{figure*}

We consider optimizing a real-valued grey-box function $g$ on $\mathcal{X}=[0,1]^d$ based on noise-free point queries.  (We briefly comment on noisy settings in Section \ref{sec:conclusion}.) As shown in \cref{fig:setup}, the target function $g$ is known to be a network of $m$ unknown layers $\ff{i}$ with $i\in[m]$. In general, for any input ${\mathbf{x}}\in\mathcal{X}$, the network $g$ has a structure such that
\begin{align*}
    \xx{1}&={\mathbf{x}}, \\
    \xx{i+1}&=\ff{i}(\xx{i})\text{ for }i\in[m-1],\\
    y&=g({\mathbf{x}})=\ff{m}(\xx{m})\in{\mathbb{R}},
\end{align*}
where $\xx{i}$ has dimension $d_i$ for each $i\in[m]$. The domain of $\ff{i}$ is $\mathcal{X}^{(i)}$, and the range of $\ff{i}$ is $\mathcal{X}^{(i+1)}$.\footnote{$\mathcal{X}^{(1)}$ is an alternative notation for $\mathcal{X}$, and $d_1$ is an alternative notation for $d$.} For any ${\mathbf{z}}\in\mathcal{X}^{(i)}$, there exists a corresponding ${\mathbf{x}}\in\mathcal{X}$ such that $\xx{i}={\mathbf{z}}$.

We aim to find ${\mathbf{x}}^\ast=\argmax_{{\mathbf{x}}\in\mathcal{X}}g({\mathbf{x}})$ based on a sequence of point queries up to time horizon $T$. When we query $g$ with input ${\mathbf{x}}_t$ at time step $t$, the intermediate noise-free results $\xx{2}_t,\dots,\xx{m}_t$ and the final noise-free output $y_t$ are accessible. We measure the performance using the following two notions:
\begin{itemize} \itemsep0ex
    \item \textbf{Simple regret}: With ${\mathbf{x}}_T^\ast$ being  additional point returned after $T$ rounds, the simple regret is defined as $r_T^\ast=g({\mathbf{x}}^\ast)-g({\mathbf{x}}_T^\ast)$;
    \item \textbf{Cumulative regret}: The cumulative regret incurred over $T$ rounds is defined as $R_T=\sum_{i=1}^T r_t$ with $r_t=g({\mathbf{x}}^\ast)-g({\mathbf{x}}_t)$.
\end{itemize}

\subsection{Kernelized Bandits}

We assume $g$ is a composition of multiple constituent functions, for which we consider both scalar-valued functions and vector-valued functions based on a given kernel.

\textbf{Scalar-valued functions.}
For a given scalar-valued kernel $k:\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}$, consider a function space $\mathcal{S}_k:=\{f(\cdot)=\sum_{i=1}^{n_0} a_i k(\cdot,{\mathbf{x}}_i): n_0\in{\mathbb{N}}, a_i\in{\mathbb{R}}, {\mathbf{x}}_i\in \mathcal{X}\}$. Then, the reproducing kernel Hilbert space (RKHS) corresponding to kernel $k$, denoted by $\mathcal{H}_k$, can be obtained by forming the completion of $\mathcal{S}_k$, and the elements in $\mathcal{H}_k$ are called scalar-valued kernelized bandits or GP bandits.
$\mathcal{H}_k$ is equipped with the inner product
\begin{align}
    \langle f, f' \rangle_k = \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} a_i b_j k({\mathbf{x}}_i,{\mathbf{x}}'_j)\label{eq:norm}
\end{align}
for $f=\sum_{i=1}^{n_1} a_i k(\cdot,{\mathbf{x}}_i)$ and $f'=\sum_{j=1}^{n_2} b_j k(\cdot,{\mathbf{x}}'_j)$. This inner product satisfies the reproducing property, such that $\langle f(\cdot), k(\cdot,{\mathbf{x}})\rangle_k=f({\mathbf{x}}),\forall {\mathbf{x}}\in\mathcal{X},\forall f\in\mathcal{H}_k$. The RKHS norm of $f$ is $\|f\|_k=\sqrt{\langle f, f\rangle_k}$, and we use $\mathcal{H}_k(B):=\{f\in\mathcal{H}_k:\|f\|_k\le B\}$ to denote the set of functions whose RKHS norm is upper bounded by some known constant $B>0$. In this work, we mainly focus on the Mat\'ern kernel:
\begin{align*}
	k_\text{Mat\'ern}({\mathbf{x}},{\mathbf{x}}')&=\frac{2^{1-\nu}}{\Gamma(\nu)}\Bigg(\frac{\sqrt{2\nu}d_{{\mathbf{x}},{\mathbf{x}}'}}{l}\Bigg)^\nu B_\nu \Bigg(\frac{\sqrt{2\nu}d_{{\mathbf{x}},{\mathbf{x}}'}}{l}\Bigg),
\end{align*}
where $d_{{\mathbf{x}},{\mathbf{x}}'}=\| {\mathbf{x}}-{\mathbf{x}}'\|_2$, $l>0$ denotes the length-scale, $\nu>0$ is a smoothness parameter, $\Gamma$ is the Gamma function, and $B_\nu$ is the modified Bessel function.

\textbf{Vector-valued functions.}
An operator-valued kernel $\Gamma:\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}^{n\times n}$ is called a multi-task kernel on $\mathcal{X}$ if $\Gamma(\cdot,\cdot)$ is symmetric positive definite. Moreover, a single-task kernel with $n=1$ recovers a scalar-valued kernel. For a given multi-task kernel $\Gamma$, similarly to the scalar-valued kernels, there exists an RKHS of vector-valued functions $\mathcal{H}_\Gamma$, which is the completion of $\mathcal{S}_\Gamma:=\{f(\cdot)=\sum_{i=1}^{n_0} \Gamma(\cdot,{\mathbf{x}}_i){\mathbf{a}}_i: n_0\in{\mathbb{N}}, {\mathbf{a}}_i\in{\mathbb{R}}^n, {\mathbf{x}}_i\in \mathcal{X}\}$. The elements in $\mathcal{H}_\Gamma$ are called vector-valued kernelized bandits, and $\mathcal{H}_\Gamma$ is equipped with the inner product \cite{carmeli2006vector}
\begin{align}
    \langle f,f'\rangle_\Gamma=\sum_{i=1}^{n_1}\sum_{j=1}^{n_2}\langle \Gamma({\mathbf{x}}_i,{\mathbf{x}}_j') {\mathbf{a}}_i,{\mathbf{b}}_j\rangle\label{eq:mul_norm}
\end{align}
for $f=\sum_{i=1}^{n_1}\Gamma(\cdot,{\mathbf{x}}_i){\mathbf{a}}_i$ and $f'=\sum_{j=1}^{n_2}\Gamma(\cdot,{\mathbf{x}}'_j){\mathbf{b}}_i$. This inner product satisfies the reproducing property, such that $\langle f(\cdot),\Gamma(\cdot,{\mathbf{x}}){\mathbf{v}}\rangle_\Gamma=\langle f({\mathbf{x}}),{\mathbf{v}}\rangle_2,\forall{\mathbf{x}}\in\mathcal{X},\forall {\mathbf{v}}\in{\mathbb{R}}^n,\forall f\in\mathcal{H}_\Gamma.$ The RKHS norm of $f\in\mathcal{H}_\Gamma$ is $\|f\|_\Gamma=\sqrt{\langle f,f\rangle_\Gamma}$, and we focus on $\mathcal{H}_\Gamma(B):=\{f\in\mathcal{H}_\Gamma:\|f\|_\Gamma\le B\}$ containing functions with norm at most $B>0$.  We will often pay particular attention to the Mat\'ern kernel $\Gamma_\text{Mat\'ern}(\cdot,\cdot)=k_\text{Mat\'ern}(\cdot,\cdot){\mathbf{I}}$, where $k_\text{Mat\'ern}$ is the scalar-valued Mat\'ern kernel and ${\mathbf{I}}$ is the identity matrix of size $n$.

\textbf{Surrogate GP model.}
As is common in kernelized bandit problems, our algorithms employ a surrogate Bayesian GP model for $f\in\mathcal{H}_k(B)$. For $f$ sampled from the prior with zero mean and kernel $k$, given a sequence of points $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ and their noise-free observations $\mathbf{y}_t=(y_1,\dots,y_t)$ up to time $t$, the posterior distribution of $f$ is a GP with mean and variance given by \cite{rasmussen06}
\begin{align}
    \mu_t({\mathbf{x}}) &= \mathbf{k}_t({\mathbf{x}})^T \mathbf{K}_t^{-1}\mathbf{y}_t \label{eq:mean}\\
    \sigma_t({\mathbf{x}})^2  &= k({\mathbf{x}},{\mathbf{x}}) - \mathbf{k}_t({\mathbf{x}})^T \mathbf{K}_t^{-1}\mathbf{k}_t({\mathbf{x}}),\label{eq:var}
\end{align}
where $\mathbf{k}_t({\mathbf{x}}) = [k({\mathbf{x}}, {\mathbf{x}}_i)]_{i=1}^{t}\in{\mathbb{R}}^{t\times 1}$ and $\mathbf{K}_t = [k({\mathbf{x}}_i, {\mathbf{x}}_j)]_{i,j=1}^t\in{\mathbb{R}}^{t\times t}$. The following lemma shows that the posterior confidence region defined with parameter $B$ is always deterministically valid.
\begin{lem}\textup{\citep[Corollary 3.11]{kanagawa2018gaussian}}
\label{lem:conf}
For $f\in\mathcal{H}_k(B)$, let $\mu_t({\mathbf{x}})$ and $\sigma_t({\mathbf{x}})^2$ denote the posterior mean and variance based on $t$ points $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ and their noise-free observations $(y_1,\dots,y_t)$ using \eqref{eq:mean} and \eqref{eq:var}. Then, it holds for all ${\mathbf{x}}\in\mathcal{X}$ that
\begin{align*}
    |f({\mathbf{x}})-\mu_t({\mathbf{x}})|\le B\sigma_t({\mathbf{x}}).
\end{align*}
\end{lem}
We also impose a surrogate GP model for functions in $\mathcal{H}_\Gamma(B)$. The posterior mean and variance matrix based on $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ and their noise-free observations $Y_t=({\mathbf{y}}_1,\dots, {\mathbf{y}}_t)$ are given by \cite{chowdhury2021no}
\begin{align}
    \mu_t({\mathbf{x}})&=G_t({\mathbf{x}})^T G_t^{-1}Y_t,\label{eq:multi_mean}\\
    \Gamma_t({\mathbf{x}},{\mathbf{x}})&=\Gamma({\mathbf{x}},{\mathbf{x}})-G_t({\mathbf{x}})^T G_t^{-1}G_t({\mathbf{x}})\label{eq:multi_var},
\end{align}
where $G_t({\mathbf{x}})=[\Gamma({\mathbf{x}},{\mathbf{x}}_i)]_{i=1}^t \in {\mathbb{R}}^{nt\times n}$, $G_t=[\Gamma({\mathbf{x}}_i,{\mathbf{x}}_j)]_{i,j=1}^t \in {\mathbb{R}}^{nt\times nt}$, and $Y_t=[{\mathbf{y}}_i]_{i=1}^t\in{\mathbb{R}}^{nt\times 1}$. With $\|\cdot\|_2$ denoting the spectral norm of a matrix, the following lemma provides a deterministic confidence region.

\begin{restatable}{lem}{mulconf}
\label{lem:mul}
For $f\in\mathcal{H}_\Gamma(B)$, let $\mu_t({\mathbf{x}})$ and $\Gamma_t({\mathbf{x}},{\mathbf{x}})$ denote the posterior mean and variance based on $t$ points $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ and their noise-free observations $({\mathbf{y}}_1,\dots,{\mathbf{y}}_t)$ using \eqref{eq:multi_mean} and \eqref{eq:multi_var}. Then, it holds for all ${\mathbf{x}}\in\mathcal{X}$ that
\begin{align*}
    \|f({\mathbf{x}})-\mu_t({\mathbf{x}})\|_2\le B \|\Gamma_t({\mathbf{x}},{\mathbf{x}})\|_2^{1/2}.
\end{align*}
\end{restatable}
The proof is given in \cref{sec:cr_mul}.

\subsection{Lipschitz Continuity}

We also assume that each constituent function in the network $g$ is Lipschitz continuous. For a constant $L> 1$, we denote by $\mathcal{F}(L)$ the set of functions such that
\begin{align*}
    \mathcal{F}(L)=\{f:\|f({\mathbf{x}})-f({\mathbf{x}}')\|_2\le L\|{\mathbf{x}}-{\mathbf{x}}'\|_2,\forall {\mathbf{x}},{\mathbf{x}}'\},
\end{align*}
where $L$ is called the Lipschitz constant. This is a mild assumption, as \cite{lee2022multi} has shown that Lipschitz continuity is a guarantee for functions in $\mathcal{H}_k(B)$ for the commonly-used squared exponential kernel and Mat\'ern kernel with smoothness $\nu>1$.

\subsection{Network structures.}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\linewidth]{figs/chain.pdf}
    \caption{A chain with $m=3$, where $\{\ff{i}\}_{i=1}^3$ are scalar-valued functions.}
    \label{fig:chain}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\linewidth]{figs/mul.pdf}
    \caption{A multi-output chain with $m=3$, where $\{\ff{i}\}_{i=1}^3$ are vector-valued functions.}
    \label{fig:mul}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\linewidth]{figs/net.pdf}
    \caption{A feed-forward network with $m=3$, where $\ff{i,j}$ is a scalar-valued function for each $i\in[3],j\in[d_{i+1}]$.}
    \label{fig:net}
\end{figure}

In this work, we consider three types of network structure for $g$:
\begin{itemize}
    \item \textbf{Chain}: For a scalar-valued kernel $k$, a chain is a cascade of scalar-valued functions. Specifically, $d_1\ge 1$, $d_2=d_3=\cdots=d_m=1$, and $\ff{i}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m]$. An example is given in \cref{fig:chain}.
    \item \textbf{Multi-output chain}: For an operator-valued kernel $\Gamma$, a multi-output chain is a cascade of vector-valued functions. Specifically, $d_i\ge 1$ and $\ff{i}\in\mathcal{H}_\Gamma(B)\cap\mathcal{F}(L)$ for each $i\in[m]$. An example is given in \cref{fig:mul}.
    \item \textbf{Feed-forward network}: For a scalar-valued kernel $k$, a feed-forward network is a fully-connected feed-forward network of scalar-valued functions. Specifically, $d_i\ge 1$ and $\ff{i}({\mathbf{z}})=[\ff{i,j}({\mathbf{z}})]_{j=1}^{d_{i+1}}$ with $\ff{i,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m],j\in[d_{i+1}]$. An example is given in \cref{fig:net}.
\end{itemize}
For all the cases above, the overall network $g$ is scalar-valued with the dimension of the final output $y$ being $d_{m+1}=1$.

\section{GPN-UCB Algorithm and Regret Bounds}

\begin{algorithm}[!t]
\caption{GPN-UCB (Gaussian Process Network - Upper Confidence Bound)}
\label{algo:ucb}
\begin{algorithmic}[1]

    \FORALL{$t \gets 1, 2, \dots, T$}
        \STATE Select ${\mathbf{x}}_t\gets\argmax_{{\mathbf{x}}\in\mathcal{X}}\mathrm{UCB}_{t-1}({\mathbf{x}})$
        \STATE Obtain observations $\xx{2}_t,\dots,\xx{m}_t,$ and $y_t$.
        \STATE Compute $\mathrm{UCB}_t$ using \eqref{eq:ucb_chain}, \eqref{eq:ucb_mul}, or \eqref{eq:ucb_net} based on $\{\xx{1}_s,\dots,\xx{m}_s,y_s\}_{s=1}^t$.
    \ENDFOR
\end{algorithmic}
\end{algorithm}

In this section, we propose a fully sequential algorithm GPN-UCB (see \cref{algo:ucb}) for chains, multi-output chains, and feed-forward networks. The algorithm works with structure-specific upper confidence bounds. Similar to GP-UCB for scalar-valued functions \cite{Sri09}, the proposed algorithm repeatedly queries the point with the highest posterior upper confidence bound, while the posterior upper confidence bound $\mathrm{UCB}_{t-1}$ used here is computed based on not only the historical final outputs $\{y_s\}_{s=1}^{t-1}$ but also the intermediate results $\{\xx{2}_s,\dots,\xx{m}_s\}_{s=1}^{t-1}$.

We note that our goal in this paper is establishing what regret bounds are fundamentally attainable, and some components of our algorithm may be impractical as stated.  We comment more on computation in Remark \ref{rem:comp}.

\subsection{GPN-UCB for Chains}

A chain is a cascade of scalar-valued functions. For each $i\in[m]$, we denote by $\muu{i}_t$ and $\sig{i}_t$ the posterior mean and standard deviation of $\ff{i}$ computed using \eqref{eq:mean} and \eqref{eq:var} based on $\{\xx{i}_s,\xx{i+1}_s\}_{s=1}^{t}$.\footnote{$\xx{1}$ is an alternative notation for ${\mathbf{x}}$, and $\xx{m+1}$ is an alternative notation for $y$.} Then, based on \cref{lem:conf}, the upper confidence bound and lower confidence bound of $\ff{i}({\mathbf{z}})$ based on $t$ exact observations are defined as follows:
\begin{align}
    \UCBB{i}_t({\mathbf{z}}) &= \muu{i}_{t}({\mathbf{z}})+B\sig{i}_{t}({\mathbf{z}}), \label{eq:UCB_orig}\\
    \LCBB{i}_t({\mathbf{z}}) &= \muu{i}_{t}({\mathbf{z}})-B\sig{i}_{t}({\mathbf{z}}). \label{eq:LCB_orig}
\end{align}
Since $\ff{i}\in\mathcal{F}(L)$, we have for any ${\mathbf{z}},{\mathbf{z}}'$ that
\begin{align}
    \UCBB{i}_t({\mathbf{z}}')+L\|{\mathbf{z}}-{\mathbf{z}}'\|_2
    &\ge \ff{i}({\mathbf{z}}')+L\|{\mathbf{z}}-{\mathbf{z}}'\|_2 \nonumber\\
    &\ge \ff{i}({\mathbf{z}}),\label{eq:lip_ucb}\\
    \LCBB{i}_t({\mathbf{z}}')-L\|{\mathbf{z}}-{\mathbf{z}}'\|_2
    &\le \ff{i}({\mathbf{z}}')-L\|{\mathbf{z}}-{\mathbf{z}}'\|_2 \nonumber\\
    &\le \ff{i}({\mathbf{z}}).
\end{align}
It follows that
\begin{align}
    \UCBBB{i}_t({\mathbf{z}})&:=\min_{{\mathbf{z}}'}\big(\UCBB{i}_t({\mathbf{z}}')+L\|{\mathbf{z}}-{\mathbf{z}}'\|_2\big),\label{eq:ucb_i}\\
    \LCBBB{i}_t({\mathbf{z}})&:=\max_{{\mathbf{z}}'}\big(\LCBB{i}_t({\mathbf{z}}')-L\|{\mathbf{z}}-{\mathbf{z}}'\|_2\big)\label{eq:lcb_i}
\end{align}
are also valid confidence bounds for $\ff{i}({\mathbf{z}})$. $\UCBBB{i}_t({\mathbf{z}})$ is the lower envelope of a collection of upper bounds for $\ff{i}({\mathbf{z}})$, which can be obtained by considering multiple values of ${\mathbf{z}}'$ in \eqref{eq:lip_ucb}. Then, since $g$ is a cascade of $\ff{i}$'s, for any input ${\mathbf{x}}$, we can recursively construct a confidence region of $\xx{i+1}$ based on the confidence region of $\xx{i}$, and the following upper confidence bound on $g({\mathbf{x}})$ is valid:
\begin{align}
    \mathrm{UCB}_t({\mathbf{x}})&=\max_{{\mathbf{z}}\in\del{m}_t({\mathbf{x}})}\UCBBB{m}_t({\mathbf{z}}),\label{eq:ucb_chain}

\end{align}
where $\del{i}_t({\mathbf{x}})$ denotes the confidence region of $\xx{i}$:
\begin{align*}
    \del{1}_t({\mathbf{x}})&=\{{\mathbf{x}}\},\\
    \del{i+1}_t({\mathbf{x}})&=\Bigg[\min_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\LCBBB{i}_t({\mathbf{z}}),\max_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\UCBBB{i}_t({\mathbf{z}})\Bigg]
\end{align*}
\hfill for $i\in[m-1]$.

The theoretical performance of \cref{algo:ucb} for chains using the upper confidence bound in \eqref{eq:ucb_chain} is provided in the following theorem.

\begin{restatable}[GPN-UCB for chains]{thm}{ucbchain}
\label{thm:ucb_chain}
Under the setup of \cref{sec:setup}, for $B>0$ and $L> 1$, a scalar-valued kernel $k$, and a chain $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m]$, \cref{algo:ucb} achieves
\begin{align*}
    R_T\le 2^{m+1} BL^{m-1}\Sigma_T,
\end{align*}
where $\Sigma_T=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}} \sum_{t=1}^T \sig{i}_{t-1}({\mathbf{z}}_t)$.
\end{restatable}

The proof is given in \cref{sec:ucb_chain}, and upper bounds on $\Sigma_T$ will be discussed in \cref{sec:ucb_disc}.  Regardless of such upper bounds, we note that $B\Sigma_T$ serves as a noise-free regret bound for standard GP optimization \cite{vakili2022open}, and thus, the key distinction here is the multiplication by $L^{m-1}$.  See Section \ref{sec:lower} for a study of the extent to which this dependence is unavoidable.

\begin{rem} \label{rem:comp}
    As stated above, our algorithm may be difficult to implement exactly in practice; in particular:
    \begin{itemize} \itemsep0ex
        \item Since $\mathcal{X}^{(2)},\dots,\mathcal{X}^{(m)}$ are not known, \eqref{eq:ucb_i} and \eqref{eq:lcb_i} are computed based on all ${\mathbf{z}}'\in{\mathbb{R}}^{d_i}$, which may be difficult.  On the other hand, we view this as a ``theoretical trick'', and we expect that in practical scenarios the effect of skipping this step and using \eqref{eq:UCB_orig}--\eqref{eq:LCB_orig} in place of \eqref{eq:ucb_i}--\eqref{eq:lcb_i} would be minimal.
        \item Recursively computing \eqref{eq:ucb_chain} may also be challenging, but some practical methods for approximating this kind of acquisition function are discussed in \citep[Sec.~4.3]{astudillo2021bayesian}.
    \end{itemize}
The same goes for the variations given below.
\end{rem}

\subsection{GPN-UCB for Multi-Output Chains}

A multi-output chain is a cascade of vector-valued functions. For any input ${\mathbf{z}}$ of the multi-output function $\ff{i}$, we define the confidence region of $\ff{i}({\mathbf{z}})$ as
\begin{align}
    \overbar{\mathcal{C}}^{(i)}_t({\mathbf{z}})&=\bigcap_{{\mathbf{z}}'}\mathcal{C}^{(i)}_t({\mathbf{z}},{\mathbf{z}}'),\label{eq:cbarz}
\end{align}
with
\begin{align}
    &\mathcal{C}^{(i)}_t({\mathbf{z}}')=\{\muu{i}_t({\mathbf{z}}')+ {\mathbf{u}}:\nonumber\\
    &\hspace{4ex}{\mathbf{u}}\in{\mathbb{R}}^{d_{i+1}}, \|{\mathbf{u}}\|_2\le B\|\Gamma_t^{(i)}({\mathbf{z}},{\mathbf{z}})\|_2^{1/2} \},\label{eq:cz}\\
    &\mathcal{C}^{(i)}_t({\mathbf{z}},{\mathbf{z}}')=\{{\mathbf{v}}+{\mathbf{w}}:\nonumber\\
    &\hspace{4ex}{\mathbf{v}}\in\mathcal{C}^{(i)}_t({\mathbf{z}}'),\|{\mathbf{w}}\|_2\le L\|{\mathbf{z}}-{\mathbf{z}}'\|_2\}.\label{eq:czz}
\end{align}
\cref{lem:mul} shows that $\mathcal{C}^{(i)}_t({\mathbf{z}}')$ is a valid deterministic confidence region for $\ff{i}({\mathbf{z}}')$. Assuming $\ff{i}\in\mathcal{F}(L)$, $\{\ff{i}({\mathbf{z}}')+{\mathbf{w}}: {\mathbf{w}}\in{\mathbb{R}}^{d_{i+1}}, \|{\mathbf{w}}\|_2\le L\|{\mathbf{z}}-{\mathbf{z}}'\|_2\}$ containing all the points satisfying the Lipchitz property is a valid confidence region for $\ff{i}({\mathbf{z}})$, and therefore its superset $\mathcal{C}^{(i)}_t({\mathbf{z}},{\mathbf{z}}')$ is also a valid confidence region for $\ff{i}({\mathbf{z}})$. Since $\ff{i}({\mathbf{z}})$ must belong to the intersection of all its confidence regions, $\overbar{\mathcal{C}}^{(i)}_t({\mathbf{z}})$ is again a valid deterministic confidence region for $\ff{i}({\mathbf{z}})$.
Hence, noting that $\overbar{\mathcal{C}}^{(m)}_{t}$ is a subset of ${\mathbb{R}}$ (unlike the vector-valued layers), the upper confidence bound for $g({\mathbf{x}})$ for any input ${\mathbf{x}}$ based on $t$ observations is
\begin{align}
    \mathrm{UCB}_{t}({\mathbf{x}})&=\max_{{\mathbf{z}}\in\del{m}_{t}({\mathbf{x}})} \overbar{\mathcal{C}}^{(m)}_{t}({\mathbf{z}}),\label{eq:ucb_mul}

\end{align}
where $\del{i}_t({\mathbf{x}})$ denotes the confidence region of $\xx{i}$ based on $t$ observations such that
\begin{align}
    \del{1}_t({\mathbf{x}})&=\{{\mathbf{x}}\},\\
    \del{i+1}_t({\mathbf{x}})&=\bigcup_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})} \overbar{\mathcal{C}}^{(i)}_t({\mathbf{z}})\text{ for }i\in[m-1].\label{eq:del}
\end{align}

The cumulative regret achieved by \cref{algo:ucb} for multi-output chains using the upper confidence bound in \eqref{eq:ucb_mul} is provided in the following theorem.

\begin{restatable}[GPN-UCB for multi-output chains]{thm}{ucbmul}
\label{thm:ucb_mul}
Under the setup of \cref{sec:setup}, given $B>0$ and $L> 1$, an operator-valued kernel $\Gamma$, and a multi-output chain $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}\in\mathcal{H}_\Gamma(B)\cap\mathcal{F}(L)$ for each $i\in[m]$, \cref{algo:ucb} achieves
\begin{align*}
    R_T\le 5^m BL^{m-1}\Sigma^\Gamma_T,
\end{align*}
where $\Sigma^\Gamma_T=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}} \sum_{t=1}^T\|\Gamma^{(i)}_{t-1}({\mathbf{z}}_t,{\mathbf{z}}_t)\|_2^{1/2}$.
\end{restatable}
The proof is given in \cref{sec:ucb_mul}.

\begin{rem}\label{rem:mul}
Fix $\mathcal{X}^{(1)},\mathcal{X}^{(2)},\dots,\mathcal{X}^{(m)}$ with dimension $d_1,d_2,\dots,d_m\ge 1$ respectively. For $\Gamma(\cdot,\cdot)=k(\cdot,\cdot){\mathbf{I}}$ with $k$ being a scalar-valued kernel and ${\mathbf{I}}$ being the identity matrix of size $d_{i+1}$, recalling that
\begin{align*}
	\Sigma_T&=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}} \sum_{t=1}^T \sig{i}_{t-1}({\mathbf{z}}_t),\\
	\Sigma^\Gamma_T&=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}} \sum_{t=1}^T\|\Gamma^{(i)}_{t-1}({\mathbf{z}}_t,{\mathbf{z}}_t)\|_2^{1/2},
\end{align*}
it follows from \cref{lem:multi_var} (see \cref{sec:multi_var}) that $\Sigma^\Gamma_T=\Sigma_T$.
\end{rem}

The upper bound on $\Sigma^\Gamma_T$ for general operator-valued kernels will be discussed in \cref{sec:ucb_disc}.

\subsection{GPN-UCB for Feed-Forward Networks}

In the feed-forward network structure, $\ff{i}({\mathbf{z}})=[\ff{i,j}({\mathbf{z}})]_{j=1}^{d_{i+1}}$ and each $\ff{i,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ is a scalar-valued function. Similar to \eqref{eq:ucb_i} and \eqref{eq:lcb_i}, with $\muu{i,j}_t({\mathbf{z}})$ and $\sig{i,j}_t({\mathbf{z}})^2$ denoting the posterior mean and variance of $\ff{i,j}({\mathbf{z}})$ using \eqref{eq:mean} and \eqref{eq:var}, the following confidence bounds on $\ff{i,j}({\mathbf{z}})$ based on $\{(\xx{i}_s,\xx{i+1,j}_s)\}_{s=1}^t$ are valid:
\begin{align}
    \UCBBB{i,j}_t({\mathbf{z}})&=\min_{{\mathbf{z}}'}(\UCBB{i,j}_t({\mathbf{z}}')+L\|{\mathbf{z}}-{\mathbf{z}}'\|_2), \nonumber \\
    \LCBBB{i,j}_t({\mathbf{z}})&=\max_{{\mathbf{z}}'}(\LCBB{i,j}_t({\mathbf{z}}')-L\|{\mathbf{z}}-{\mathbf{z}}'\|_2), \label{eq:ucb_ij_bar}
\end{align}
where
\begin{align}
    \UCBB{i,j}_t({\mathbf{z}}) &= \muu{i,j}_{t}({\mathbf{z}})+B\sig{i,j}_{t}({\mathbf{z}}), \nonumber \\
    \LCBB{i,j}_t({\mathbf{z}}) &= \muu{i,j}_{t}({\mathbf{z}})-B\sig{i,j}_{t}({\mathbf{z}}). \label{eq:ucb_ij_final}
\end{align}
Then, the upper confidence bound of $g({\mathbf{x}})$ based on $t$ observations is
\begin{align}
    \mathrm{UCB}_t({\mathbf{x}})&=\max_{{\mathbf{z}}\in\del{m}_t({\mathbf{x}})}\UCBBB{m,1}_t({\mathbf{z}}),\label{eq:ucb_net}

\end{align}
where
\begin{align}
    \del{1}_t({\mathbf{x}})=&\{{\mathbf{x}}\},\nonumber\\
    \del{i+1,j}_t({\mathbf{x}})=&\Bigg[\min_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\LCBBB{i,j}_t({\mathbf{z}}) , \max_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\UCBBB{i,j}_t({\mathbf{z}}) \Bigg]\nonumber\\
    &\hspace{12ex}\text{ for }i\in[m-1],j\in[d_{i+1}],\nonumber\\
    \del{i}_t({\mathbf{x}})=& \del{i,1}_t({\mathbf{x}})\times\cdots\times\del{i,d_i}_t({\mathbf{x}})\text{ for }i\in[m].\label{eq:del_d}
\end{align}

The following theorem provides the theoretical performance of \cref{algo:ucb} for feed-forward networks using the upper confidence bound in \eqref{eq:ucb_net}.

\begin{restatable}[GPN-UCB for feed-forward networks]{thm}{ucbnet}
\label{thm:ucb_net}
Under the setup of \cref{sec:setup}, given $B>0$ and $L> 1$, a scalar-valued kernel $k$, and a feed-forward network $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}({\mathbf{z}})=[\ff{i,j}({\mathbf{z}})]_{j=1}^{d_{i+1}}$ and $\ff{i,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m],j\in[d_{i+1}]$, \cref{algo:ucb} achieves
\begin{align*}
    R_T\le 2^{m+1} \sqrt{D_{2,m}} B L^{m-1} \Sigma_T,
\end{align*}
where $D_{2,m}=\prod_{i=2}^{m} d_i$ and $\Sigma_T=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}}\sum_{t=1}^T \sig{i,1}_{t-1}({\mathbf{z}}_t)$.
\end{restatable}
The proof is given in \cref{sec:ucb_net}.

Since $\sig{i,1}$ in the feed-forward network setting is computed using \eqref{eq:var}, which is exactly the same as how $\sig{i}$ is computed in the chain setting, despite the slightly different superscripts, the $\Sigma_T$ term in \cref{thm:ucb_net} represents the same quantity as in \cref{thm:ucb_chain} (depending on $\mathcal{X}^{(1)},\mathcal{X}^{(2)},\dots,\mathcal{X}^{(m)}$). In particular, when $d_2=\cdots=d_m=1$, \cref{thm:ucb_net} recovers \cref{thm:ucb_chain}.

\subsection{Upper Bounds on $\Sigma_T$ and $\Sigma^\Gamma_T$}
\label{sec:ucb_disc}

A simple way to establish an upper bound on the $\Sigma_T$ term in \cref{thm:ucb_chain}, \cref{rem:mul}, and \cref{thm:ucb_net} is to essentially set the noise term to be zero in a known result for the noisy setting. For the scalar-valued function (GP bandit) optimization problem under the noisy setting, most existing upper bounds on cumulative regret are expressed in terms of the maximum information gain corresponding to the kernel defined as \cite{Sri09}
\begin{align*}
    \gamma_t=\max_{{\mathbf{x}}_1,\dots,{\mathbf{x}}_t}\frac{1}{2}\log\det(\mathbf{I}_t+\lambda^{-1}\mathbf{K}_t)
\end{align*}
for a free parameter $\lambda>0$, and \cite{Sri09} has shown that the sum of posterior variances in the noisy setting satisfies $\sum_{t=1}^T\sigma'_{t-1}({\mathbf{x}}_t)^2=O(\gamma_T)$, where
\begin{align*}
    \sigma'_t({\mathbf{x}})^2  &= k({\mathbf{x}},{\mathbf{x}}) - \mathbf{k}_t({\mathbf{x}})^T (\mathbf{K}_t+\lambda\mathbf{I}_t)^{-1}\mathbf{k}_t({\mathbf{x}}).
\end{align*}
Using $\sigma_t({\mathbf{x}})\le\sigma'_t({\mathbf{x}})$ and the Cauchy-Schwartz inequality, we have $\Sigma_T=O\big(\sqrt{T\gamma_T}\big)$. An existing upper bound on $\gamma_T$ for the Mat\'ern kernel with dimension $d$ and smoothness $\nu$ on a fixed compact domain is \cite{Vak20a}
\begin{align*}
    \gamma_T^\text{Mat\'ern}={\widetilde{O}}\big(T^\frac{d}{2\nu+d}\big).
\end{align*}
In our setting, a simple sufficient condition for this bound to apply is that  $d_{\max}=\max_{i\in[m]}d_i$, $L$, and $m$ are constant, since then the Lipschitz assumption implies that each domain $\mathcal{X}^{(1)},\mathcal{X}^{(2)},\dots,\mathcal{X}^{(m)}$ is also compact/bounded.  More generally, we believe that uniformly bounded domains is a mild assumption, and when it holds, the bound  $\Sigma_T=O\big(\sqrt{T\gamma_T}\big)$ simplifies to $\Sigma_T=O\big(T^\frac{\nu+d_{\max}}{2\nu+d_{\max}}\big)$.

In addition, \cite{vakili2022open} provides the following conjecture on the upper bound on $\Sigma_T$ for the Mat\'ern kernel
\begin{align*}
    \Sigma^\text{Mat\'ern}_T=\begin{cases}
    O(T^{1-\nu/d_{\max}}) &\text{ when }d_{\max}>\nu,\\
    O(\log T)&\text{ when }d_{\max}=\nu,\\
    O(1)&\text{ when }d_{\max}<\nu.
    \end{cases}
\end{align*}
We will discuss in Section \cref{sec:discussion} how if this conjecture is true, we can deduce the near-optimality of GPN-UCB for the Mat\'ern kernel.

For an arbitrary operator-valued kernel $\Gamma:\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}^{n\times n}$ and a free parameter $\lambda$, the maximum information gain is defined as \cite{chowdhury2021no}
\begin{align*}
    \gamma^\Gamma_t=\max_{{\mathbf{x}}_1,\dots,{\mathbf{x}}_t}\frac{1}{2}\log\det(\mathbf{I}_{nt}+\lambda^{-1}G_t).
\end{align*}
\cite{chowdhury2021no} has shown that $\sum_{t=1}^T\|\Gamma_{t-1}({\mathbf{x}}_t,{\mathbf{x}}_t)\|_2=O(\gamma^\Gamma_T)$, and therefore $\Sigma^\Gamma_T=O\big(\sqrt{T\gamma^\Gamma_T}\big)$ by similar reasoning to that above.

\section{Non-Adaptive Sampling Based Method}

In this section, we propose a simple non-adaptive sampling based method (see \cref{algo:fd}) for each structure, and provide the corresponding theoretical simple regret for the Mat\'ern kernel. For a set of $T$ sampled points $\{{\mathbf{x}}_s\}_{s=1}^T$, its fill distance is defined as the largest distance from a point in the domain to the closest sampled point \cite{wendland2004scattered}:
\begin{align*}
    \delta_T=\max_{{\mathbf{x}}\in\mathcal{X}}\min_{s\in[T]}\|{\mathbf{x}}-{\mathbf{x}}_s\|_2.
\end{align*}
\cref{algo:fd} samples $T$ points with $\delta_T=O(T^{-1/d})$. For $\mathcal{X}=[0,1]^d$, a simple way to construct such a sample is to use a uniform $d$-dimensional grid with step size $T^{-1/d}$. The algorithm observes the selected points in parallel, computes a structure-specific ``composite mean'' $\mu_T^g$ (to be defined shortly) for the overall network $g$, and returns the point that maximizes $\mu_T^g$.

\begin{algorithm}[!t]
\caption{Non-Adaptive Sampling Based Method}
\label{algo:fd}
\begin{algorithmic}[1]
    \STATE Choosing $\{{\mathbf{x}}_s\}_{s=1}^T$ such that $\delta_T=O(T^{-\frac{1}{d}})$.
    \STATE Obtain observations $\{\xx{2}_s,\dots,\xx{m}_s,y_s\}_{s=1}^T$.
    \STATE Compute $\mu_T^g$
    based on $\{\xx{2}_s,\dots,\xx{m}_s,y_s\}_{s=1}^T$.
    \ENSURE ${\mathbf{x}}^\ast_T=\argmax_{{\mathbf{x}}\in\mathcal{X}}\mu_T^g({\mathbf{x}})$.
\end{algorithmic}
\end{algorithm}

The composite posterior mean of $g({\mathbf{x}})$ with chain structure is defined as
\begin{align}
    \mu_T^g({\mathbf{x}})=(\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}),\label{eq:mean_chain}
\end{align}
where $\muu{i}_T$ denotes the posterior mean of $\ff{i}$ computed using \eqref{eq:mean} based on $\{(\xx{i}_s,\xx{i+1}_s)\}_{s=1}^T$ for each $i\in[m]$.

Then, the following theorem provides the theoretical upper bound on the simple regret of \cref{algo:fd} using \eqref{eq:mean_chain}.  Note that the notation ${\widetilde{O}}(\cdot)$ hides poly-logarithmic factors \emph{with respect to the argument}, e.g., ${\widetilde{O}}(\sqrt{T}) = O( \sqrt{T} \cdot (\log T)^{O(1)})$ and ${\widetilde{O}}(2^{n}) = O(2^n \cdot n^{O(1)})$.

\begin{restatable}[Non-adaptive sampling method for chains]{thm}{fdchain}
\label{thm:fd_chain}
Under the setup of \cref{sec:setup}, given $B=\Theta(L)$,  $k=k_\text{Mat\'ern}$ with smoothness $\nu$, and a chain $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m]$, we have
\begin{itemize}
	\item When $\nu\le 1$, \cref{algo:fd} achieves
	\begin{align*}
		r^\ast_T={\widetilde{O}}(BL^{m-1}T^{-\nu^m/d});
	\end{align*}
	\item When $\nu>1$, \cref{algo:fd} achieves
	\begin{align*}
		r^\ast_T={\widetilde{O}}\big(\max\big\{&BL^{(m-1)\nu}T^{-\nu/d},\\ &B^{1+\nu+\nu^2+\cdots+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu^2/d} \big\}\big).
	\end{align*}
\end{itemize}
\end{restatable}
The proof is given in \cref{sec:fd_chain}, and the optimality will be discussed in \cref{sec:discussion}.

When $\nu>1$, the simple regret upper bound takes the maximum of two terms. The first term has a smaller constant factor, while the second term has a smaller $T$-dependent factor. By taking the highest-order constant factor and the highest-order $T$-dependent factor, we can deduce the weaker but simpler bound  $r^\ast_T=O(B^{1+\nu+\nu^2+\cdots+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu/d})$.

We also consider two more restrictive cases, where we remove the assumption of $B=\Theta(L)$, but have additional assumptions on $g$ as follows:
\begin{itemize}
	\item \textbf{Case 1}: We additionally assume that $\muu{i}_T\circ\cdots\circ\muu{1}_T({\mathbf{x}}^\ast)\in\mathcal{X}^{(i+1)}$ and $\muu{i}_T\circ\cdots\circ\muu{1}_T({\mathbf{x}}^\ast_T)\in \mathcal{X}^{(i+1)}$ for all $i\in[m-1]$.
	\item \textbf{Case 2}: We additionally assume that all the domains $\mathcal{X}^{(i)}$ are known. Defining
	\begin{align*}
		\muutil{i}_T({\mathbf{z}}) &= \argmin_{{\mathbf{z}}'\in\mathcal{X}^{(i+1)}}|\muu{i}_T({\mathbf{z}})-{\mathbf{z}}'|,
	\end{align*}
	we slightly modify the algorithm to return
	\begin{align*}
		{\mathbf{x}}^\ast_T=\argmax_{{\mathbf{x}}\in\mathcal{X}}\muutil{m}_T\circ\cdots\circ\muutil{1}_T({\mathbf{x}}).
	\end{align*}
\end{itemize}
\begin{rem}
	Under the assumptions of either Case 1 or Case 2, \cref{algo:fd} achieves for chains that
	\begin{align}
		r^\ast_T=
		\begin{cases}
			O(BL^{m-1}T^{-\nu/d})&\text{ when $\nu\le 1$,}\\
			O(BL^{(m-1)\nu}T^{-\nu/d})&\text{ when $\nu> 1$.}
		\end{cases}
        \label{eq:more_restrictive}
	\end{align}
\end{rem}
The proof is given in \cref{sec:rest_chain}.

The composite posterior means and simple regret upper bounds for multi-output chains and feed-forward networks are provided in \cref{sec:fd_mul} and \cref{sec:fd_net} repectively, where the simple regret upper bounds are stated only for the case that the domain of each layer is a hyperrectangle. Removing this restrictive assumption is left for future work.

\section{Algorithm-Independent Lower Bounds}
\label{sec:lower}
In this section, we provide algorithm-independent lower bounds on the simple regret and cumulative regret for any algorithm optimizing chains, multi-output chains, or feed-forward networks for the scalar-valued kernel $k_\text{Mat\'ern}$ or the operator-valued $\Gamma_\text{Mat\'ern}(\cdot,\cdot)=k_\text{Mat\'ern}(\cdot,\cdot)\mathbf{I}$ with smoothness $\nu$.

\begin{restatable}[Lower bound on simple regret]{thm}{lowersim}
\label{thm:lower_sim}
Fix $\epsilon\in(0,1)$, sufficiently large $B>0$, $k=k_\text{Mat\'ern}$, and $\Gamma=\Gamma_\text{Mat\'ern}$ with smoothness $\nu\ge 1$. Suppose that there exists an algorithm (possibly randomized) that achieves average simple regret ${\mathbb{E}}[r_T^\ast]\le\epsilon$ after $T$ rounds for any $m$-layer chain, multi-output chain, or feed-forward network on $[0,1]^d$ with some $L=\Theta(B)$. Then, provided that $\frac{\epsilon}{B}$ is sufficiently small, it is necessary that
\begin{align*}
    T=\Omega\Bigg(\Big(\frac{B(cL)^{m-1}}{\epsilon}\Big)^{d/\nu}\Bigg),
\end{align*}
for some $c=\Theta(1)$.
\end{restatable}
The proof is given in \cref{sec:lower_sim}, and the high-level steps are similar to \cite{bull2011convergence}, but the main differences are significant. For each structure, we consider a collection of $M$ hard functions $\mathcal{G}=\{g_1,\dots,g_M\}$, where each $g_j$ is obtained by shifting a base function ${\overbar{g}}$ of the specified structure and cropping the shifted function into $[0,1]^d$. Then, we show that there exists a worst-case function in $\mathcal{G}$ with the provided lower bound. Different from \cite{bull2011convergence}, the hard functions we construct here are function networks. We define the first layer as a ``needle'' function with much smaller height and width than \cite{bull2011convergence}. For subsequent layers, we construct a function with corresponding RKHS norm such that the output is always larger than the input. As a consequence, the ``needle'' function gets higher and higher when being fed into subsequent layers, and the composite function is a ``needle'' function with some specified height but a much smaller width.

The lower bound on simple regret readily implies the following lower bound on cumulative regret.
\begin{restatable}[Lower bound on cumulative regret]{thm}{lowercml}
\label{thm:lower_cml}
Fix sufficiently large $B>0$, $k=k_\text{Mat\'ern}$, and $\Gamma=\Gamma_\text{Mat\'ern}$ with smoothness $\nu\ge 1$. Suppose that there exists an algorithm (possibly randomized) that achieves average cumulative regret ${\mathbb{E}}[R_T]$ after $T$ rounds for any $m$-layer chain, multi-output chain, or feed-forward network on $[0,1]^d$ with some $L=\Theta(B)$. Then, it is necessary that
\begin{align*}
    {\mathbb{E}}[R_T]=\begin{cases}
    	\Omega\big(\min\{T, B(cL)^{m-1}T^{1-\nu/d}\}\big)&\text{ when }d>\nu,\\
    	\Omega\big(\min\{T, \big(B(cL)^{m-1}\big)^{d/\nu}\}\big)&\text{ when }d\le\nu,
    \end{cases}
\end{align*}
for some $c=\Theta(1)$.
\end{restatable}
The proof is given in \cref{sec:lower_cml}.

\section{Comparison of Upper and Lower Bounds}
\label{sec:discussion}
In this section, we compare the algorithmic upper bounds of GPN-UCB (\cref{algo:ucb}) and non-adaptive sampling (\cref{algo:fd}) to the algorithmic-independent lower bounds in \cref{sec:lower}. A full summary of the proposed regret bounds for the Mat\'ern kernel is provided in \cref{sec:summary}.

For GPN-UCB, the cumulative regret upper bound for chains (\cref{thm:ucb_chain}) matches the lower bound (\cref{thm:lower_cml}) up to a $2^m$ factor when $d\ge\nu\ge1$ and $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$. The upper bound for multi-output chains (\cref{thm:ucb_mul}) is similarly optimal (up to a $5^m$ term) when $d_{\max}=d\ge\nu\ge1$, while there is always an $O(\sqrt{D_{2,m}})$ gap for the $T$-independent factor of feed-forward networks (\cref{thm:ucb_net}). When $d<\nu$ and $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$, the cumulative regret lower bound for all the three structures is $\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$, while the upper bound always contains an $O(BL^{m-1})$ factor; hence, the terms behave similarly but there remains more room to close them further.

We expect that the discrepancies for multi-output chains and feed-forward networks are due to the looseness of the proposed lower bound. Since the hard functions $\mathcal{G}$ used in analysis always produce a single-entry vector output for intermediate layers, for a fixed value of $B$, there might exist a worse hard function network with more nonzero entries for intermediate outputs and a probably higher final regret.

For non-adaptive sampling, when $\nu=1$, the upper bound for chains (\cref{thm:fd_chain}) matches the lower bound (\cref{thm:lower_sim}) up to a $c^{m-1}$ factor. When $\nu>1$, \cref{thm:fd_chain} shows that the simple regret upper bound takes the maximum of two terms, where the first term has a matched $T$-dependent factor. However, both terms have a larger $T$-independent factor than the lower bound when $\nu>1$; in our analysis, this arises from magnifying the uncertainty from each layer to the next.

\section{Conclusion} \label{sec:conclusion}
We have proposed an upper confidence bound based method GPN-UCB and a non-adaptive sampling based method for optimizing chains, multi-output chains, and feed-forward networks in the noise-free grey-box setting. We provided cumulative regret upper bounds for the former and simple regret upper bounds for the latter. We also provided algorithm-independent lower bounds on simple regret and cumulative regret, showing that our algorithms are near-optimal in broad cases of interest.
\section{Confidence Region for Vector-Valued Functions (Proof of \cref{lem:mul})}
\label{sec:cr_mul}
In this section, we prove \cref{lem:mul}, which is restated as follows.
\mulconf*
\begin{proof}
We first review the GP posterior and confidence region for operator-valued kernel $\Gamma:\mathcal{X}\times\mathcal{X}\to{\mathbb{R}}^{n\times n}$ in the noisy setting.

\begin{restatable}{lem}{confmul}
\textup{\citep[Theorem 1]{chowdhury2021no}}
For $f\in\mathcal{H}_\Gamma(B)$, given a sequence of points $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ and their noisy observations $Y'_t=({\mathbf{y}}'_1,\dots,{\mathbf{y}}'_t)$, where ${\mathbf{y}}'_i=f({\mathbf{x}}_i)+\epsilon_i$ with $\epsilon_i$ being i.i.d. $\sigma$-sub-Gaussian for each $i\in[t]$ for some $\sigma>0$, let $\mu'_t({\mathbf{x}})$ and $\Gamma'_t({\mathbf{x}},{\mathbf{x}})$ denote the posterior mean and variance computed using
\begin{align}
    \mu'_t({\mathbf{x}})&=G_t({\mathbf{x}})^T (G_t+\lambda {\mathbf{I}}_{nt})^{-1}Y'_t,\label{eq:mu_noisy}\\
    \Gamma'_t({\mathbf{x}},{\mathbf{x}}')&=\Gamma({\mathbf{x}},{\mathbf{x}}')-G_t({\mathbf{x}})^T (G_t+\lambda {\mathbf{I}}_{nt})^{-1}G_t({\mathbf{x}}'),\label{eq:gamma_noisy}
\end{align}
where $G_t({\mathbf{x}})=[\Gamma({\mathbf{x}},{\mathbf{x}}_i)]_{i=1}^t \in {\mathbb{R}}^{nt\times n}$, $G_t=[\Gamma({\mathbf{x}}_i,{\mathbf{x}}_j)]_{i,j=1}^t \in {\mathbb{R}}^{nt\times nt}$, $Y'_t=[{\mathbf{y}}'_i]_{i=1}^t\in{\mathbb{R}}^{nt\times 1}$, and $\lambda>0$ is a regularization parameter.
Then, for any $\lambda>0$ and $\delta\in(0,1]$, with probability at least $1-\delta$, it holds for all ${\mathbf{x}}\in\mathcal{X}$ that
\begin{align}
    \|f({\mathbf{x}})-\mu'_t({\mathbf{x}})\|_2\le \alpha_t\|\Gamma'_t({\mathbf{x}},{\mathbf{x}})\|_2^{1/2}
\end{align}
with $\alpha_t=B+\frac{\sigma}{\sqrt{n}}\sqrt{2\log(1/\delta)+\log\det(I_{nt}+\lambda^{-1}G_t)}$.
\end{restatable}
Since zero noise is $\sigma$-sub-Gaussian for any $\sigma>0$, by setting $\sigma=\lambda=\frac{1}{a}$ and then taking $a\to\infty$, \eqref{eq:mu_noisy} and \eqref{eq:gamma_noisy} converge to
\begin{align}
    \mu_t({\mathbf{x}})&=G_t({\mathbf{x}})^T G_t^{-1}Y_t,\\
    \Gamma_t({\mathbf{x}},{\mathbf{x}}')&=\Gamma({\mathbf{x}},{\mathbf{x}}')-G_t({\mathbf{x}})^T G_t^{-1}G_t({\mathbf{x}}'),
\end{align}
with $Y_t=({\mathbf{y}}_1,\dots, {\mathbf{y}}_t)$, thus yielding the posterior mean and variance based on $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ and their noise-free observations $Y_t$.

Then, by setting $\sigma=\delta=\lambda=\frac{1}{a}$ for $a\to\infty$, with $\{\lambda_i\}_{i=1}^{nt}$ being the eigenvalues of $G_t$, we obtain
\begin{align}
	\lim\limits_{a\to\infty}\alpha_t&=B+\frac{1}{\sqrt{n}} \cdot \lim\limits_{a\to\infty} \sqrt{\frac{2\log a}{a^2}+\frac{\sum_{i=1}^{nt} \log(a\lambda_i+1)}{a^2}} = B.
\end{align}
Hence, we obtain \cref{lem:mul} for the deterministic confidence region based on noise-free observations.
\end{proof}

\section{Posterior Variance for Vector-Valued Functions}
\label{sec:multi_var}

In this section, we state and prove \cref{lem:multi_var}.

\begin{restatable}{lem}{multivar}
\label{lem:multi_var}
For a scalar-valued kernel $k$, define $\Gamma({\mathbf{x}},{\mathbf{x}}')=k({\mathbf{x}},{\mathbf{x}}'){\mathbf{I}}_n$ with ${\mathbf{I}}_n$ being the $n\times n$ identity matrix. For $f_1\in\mathcal{H}_k(B)$ and $f_2\in\mathcal{H}_\Gamma(B)$ with domain $\mathcal{X}$ and constant $B>0$, let $\sigma_t({\mathbf{x}})^2$ and $\Gamma_t({\mathbf{x}},{\mathbf{x}})$ denote the posterior variance (matrix) for $f_1$ and $f_2$ based on $t$ points $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ in the noise-free setting computed using \eqref{eq:var} and \eqref{eq:multi_var}. Then, it holds for all ${\mathbf{x}}\in\mathcal{X}$ that
\begin{align}
    \Gamma_t({\mathbf{x}},{\mathbf{x}})=\sigma_t({\mathbf{x}})^2{\mathbf{I}}_n.
\end{align}
\end{restatable}
\begin{proof}
Let $\sigma'_t({\mathbf{x}})^2$ and $\Gamma'_t({\mathbf{x}},{\mathbf{x}})$ denote the posterior variance (matrix) for $f_1$ and $f_2$ based on $t$ points $({\mathbf{x}}_1,\dots,{\mathbf{x}}_t)$ in the noisy setting. For any $\lambda>0$, we have
\begin{align}
    \sigma'_t({\mathbf{x}})^2  &= k({\mathbf{x}},{\mathbf{x}}) - \mathbf{k}_t({\mathbf{x}})^T (\mathbf{K}_t+\lambda\mathbf{I}_t)^{-1}\mathbf{k}_t({\mathbf{x}}),\\
    \Gamma'_t({\mathbf{x}},{\mathbf{x}})&=\Gamma({\mathbf{x}},{\mathbf{x}})-G_t({\mathbf{x}})^T (G_t+\lambda {\mathbf{I}}_{nt})^{-1}G_t({\mathbf{x}}'),\label{eq:multi_var_noisy}
\end{align}
where $\mathbf{k}_t({\mathbf{x}}) = [k({\mathbf{x}}, {\mathbf{x}}_i)]_{i=1}^{t}\in{\mathbb{R}}^{t\times 1}$, $\mathbf{K}_t = [k({\mathbf{x}}_i, {\mathbf{x}}_j)]_{i,j=1}^t\in{\mathbb{R}}^{t\times t}$, $G_t({\mathbf{x}})=[\Gamma({\mathbf{x}},{\mathbf{x}}_i)]_{i=1}^t \in {\mathbb{R}}^{nt\times n}$, and $G_t=[\Gamma({\mathbf{x}}_i,{\mathbf{x}}_j)]_{i,j=1}^t \in {\mathbb{R}}^{nt\times nt}$.
With $\otimes$ denoting the Kronecker product, we have
$G_t({\mathbf{x}})= {\mathbf{k}}_t({\mathbf{x}})\otimes{\mathbf{I}}_n$ and $G_t= {\mathbf{K}}_t\otimes{\mathbf{I}}_n$. Then, it follows from \eqref{eq:multi_var_noisy} that
\begin{align}
    \Gamma'_t({\mathbf{x}},{\mathbf{x}})&=\Gamma({\mathbf{x}},{\mathbf{x}})-G_t({\mathbf{x}})^T(G_t+\lambda {\mathbf{I}}_{nt})^{-1}G_t({\mathbf{x}})\\
    &=k({\mathbf{x}},{\mathbf{x}}){\mathbf{I}}_n-\big({\mathbf{k}}_t({\mathbf{x}})^T\otimes{\mathbf{I}}_n\big)(\mathbf{K}_t\otimes {\mathbf{I}}_{n}+\lambda{\mathbf{I}}_{nt})^{-1}\big({\mathbf{k}}_t({\mathbf{x}})\otimes{\mathbf{I}}_n\big)\\
    &=k({\mathbf{x}},{\mathbf{x}}){\mathbf{I}}_n-\big({\mathbf{k}}_t({\mathbf{x}})^T\otimes{\mathbf{I}}_n\big)\big((\mathbf{K}_t+\lambda {\mathbf{I}}_t)\otimes {\mathbf{I}}_{n}\big)^{-1}\big({\mathbf{k}}_t({\mathbf{x}})\otimes{\mathbf{I}}_n\big)\\
    &=k({\mathbf{x}},{\mathbf{x}}){\mathbf{I}}_n-\big({\mathbf{k}}_t({\mathbf{x}})^T\otimes{\mathbf{I}}_n\big)\big(({\mathbf{K}}_t+\lambda {\mathbf{I}}_t)^{-1}\otimes {\mathbf{I}}_{n}\big)\big({\mathbf{k}}_t({\mathbf{x}})\otimes{\mathbf{I}}_n\big)\\
    &=k({\mathbf{x}},{\mathbf{x}}){\mathbf{I}}_n-\big({\mathbf{k}}_t({\mathbf{x}})^T({\mathbf{K}}_t+\lambda{\mathbf{I}}_t)^{-1}{\mathbf{k}}_t({\mathbf{x}})\big)\otimes{\mathbf{I}}_n\\
    &=k({\mathbf{x}},{\mathbf{x}}){\mathbf{I}}_n-\big({\mathbf{k}}_t({\mathbf{x}})^T({\mathbf{K}}_t+\lambda{\mathbf{I}}_t)^{-1}{\mathbf{k}}_t({\mathbf{x}})\big){\mathbf{I}}_n\\
    &=\big(k({\mathbf{x}},{\mathbf{x}})-({\mathbf{k}}_t({\mathbf{x}})^T({\mathbf{K}}_t+\lambda{\mathbf{I}}_t)^{-1}{\mathbf{k}}_t({\mathbf{x}})\big){\mathbf{I}}_n\\
    &=\sigma'_t({\mathbf{x}})^2{\mathbf{I}}_n
\end{align}
Taking $\lambda\to0$, we obtain in the noise-free setting that
\begin{align}
    \Gamma_t({\mathbf{x}},{\mathbf{x}})=\sigma_t({\mathbf{x}})^2{\mathbf{I}}_n.
\end{align}
\end{proof}

\section{Analysis of GPN-UCB (\cref{algo:ucb})}
\subsection{Proof of \cref{thm:ucb_chain} (Chains)}
\label{sec:ucb_chain}
In this section, we prove \cref{thm:ucb_chain}, which is restated as follows.
\ucbchain*

\begin{proof}
With ${\overbar{\mathbf{x}}}_t^{(i)}=\argmax_{{\mathbf{z}}\in\del{i}_{t-1}({\mathbf{x}}_t)}\UCBBB{i}({\mathbf{z}})$ and ${\widetilde{\mathbf{x}}}_t^{(i)}=\argmin_{{\mathbf{z}}\in\del{i}_{t-1}({\mathbf{x}}_t)}\LCBBB{i}({\mathbf{z}})$, the simple regret is upper bounded as follows:
\begin{align}
    r_t&=g({\mathbf{x}}^\ast)-g({\mathbf{x}}_t) \label{eq:sim_0} \\
    &\le \mathrm{UCB}_{t-1}({\mathbf{x}}_t)-\LCBBB{m}_{t-1}({\widetilde{\mathbf{x}}}_t^{(m)})\label{eq:sim_1}\\
    &= \UCBBB{m}_{t-1}({\overbar{\mathbf{x}}}_t^{(m)})-\LCBBB{m}_{t-1}({\widetilde{\mathbf{x}}}_t^{(m)}) \label{eq:sim_2}\\
    &\le \UCBB{m}_{t-1}(\xx{m}_t)+L\|{\overbar{\mathbf{x}}}_t^{(m)}-\xx{m}_t\|_2 -\LCBB{m}_{t-1}(\xx{m}_t)+L\|\xx{m}_t-{\widetilde{\mathbf{x}}}_t^{(m)}\|_2\label{eq:sim_3}\\
    &\le 2B\sig{m}_{t-1}(\xx{m}_t)+2L\cdot\dia{\del{m}_{t-1}({\mathbf{x}}_t)}\label{eq:sim_4}\\
    &\le 2B\Big(
    \sig{m}_{t-1}(\xx{m}_t)
    +(2L)\sig{m-1}_{t-1}(\xx{m-1}_t) +\cdots+(2L)^{m-1}\sig{1}_{t-1}(\xx{1}_t)
    \Big)\label{eq:sim_5}\\
    &= 2B\sum_{i=1}^m(2L)^{m-i}\sig{i}_{t-1}(\xx{i}_t),\label{eq:sim_6}
\end{align}
where:
\begin{itemize}
    \item \eqref{eq:sim_1} follows since $g({\mathbf{x}}^\ast)\le \mathrm{UCB}_{t-1}({\mathbf{x}}^\ast)\le \mathrm{UCB}_{t-1}({\mathbf{x}}_t)$ (due to the algorithm maximizing the UCB score) and $g({\mathbf{x}}_t)\ge\LCBBB{m}_{t-1}({\mathbf{x}}_t^{(m)})\ge \LCBBB{m}_{t-1}({\widetilde{\mathbf{x}}}_t^{(m)})$;
    \item \eqref{eq:sim_2} follows from \eqref{eq:ucb_chain} and the definition of ${\overbar{\mathbf{x}}}_t^{(m)}$;
    \item \eqref{eq:sim_3} follows from \eqref{eq:ucb_i} and \eqref{eq:lcb_i};
    \item \eqref{eq:sim_4} follows by defining $\dia{\del{i}_{t-1}({\mathbf{x}}_t)}$ as the diameter of confidence region $\del{i}_{t-1}({\mathbf{x}}_t)$;
    \item \eqref{eq:sim_5} follows from the following recursion:
    \begin{align}
    \dia{\del{i+1}_{t-1}({\mathbf{x}}_t)}&= \UCBBB{i}_{t-1}({\overbar{\mathbf{x}}}^{(i)}_t)-\LCBBB{i}_{t-1}({\widetilde{\mathbf{x}}}^{(i)}_t)\\
    &\le 2B\sig{i}_{t-1}(\xx{i}_t)+2L\cdot\dia{\del{i}_{t-1}({\mathbf{x}}_t)},
    \end{align}
    where the inequality is obtained by following \eqref{eq:sim_2}-\eqref{eq:sim_4} with $m$ replaced by $i$.
\end{itemize}
Then, the cumulative regret is
\begin{align}
    R_T&=\sum_{t=1}^T r_t\\
    &\le 2B \sum_{t=1}^T \sum_{i=1}^m(2L)^{m-i}\sig{i}_{t-1}(\xx{i}_t)\\
    &= 2B \frac{(2L)^m-1}{2L-1} \Sigma_T\\
    &\le 2^{m+1} BL^{m-1}\Sigma_T,
\end{align}
where $\Sigma_T=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}}\sum_{t=1}^T\sig{i}_{t-1}({\mathbf{z}}_t)$.
\end{proof}

\subsection{Proof of \cref{thm:ucb_mul} (Multi-Output Chains)}
\label{sec:ucb_mul}
In this section, we prove \cref{thm:ucb_mul}, which is restated as follows.
\ucbmul*
\begin{proof}
For each $i\in [m-1]$, there must exist ${\overbar{\mathbf{x}}}^{(i)}_t,{\widetilde{\mathbf{x}}}^{(i)}_t\in\del{i}_{t-1}({\mathbf{x}}_t)$ such that the following upper bound on $\dia{\del{i+1}_{t-1}({\mathbf{x}}_t)}$ in terms of $\dia{\del{i}_{t-1}({\mathbf{x}}_t)}$ holds:
\begin{align}
    \dia{\del{i+1}_{t-1}({\mathbf{x}}_t)} &\le \|\ff{i}({\overbar{\mathbf{x}}}^{(i)}_t)-\ff{i}({\widetilde{\mathbf{x}}}^{(i)}_t)\|_2+\dia{\overbar{\mathcal{C}}^{(i)}_{t-1}({\overbar{\mathbf{x}}}^{(i)}_t)}+\dia{\overbar{\mathcal{C}}^{(i)}_{t-1}({\widetilde{\mathbf{x}}}^{(i)}_t)}\label{eq:multi_del_1}\\
    &\le L\cdot\dia{\del{i}_{t-1}({\mathbf{x}}_t)}+\dia{\mathcal{C}^{(i)}_{t-1}({\overbar{\mathbf{x}}}^{(i)}_t,\xx{i}_t)}+\dia{\mathcal{C}^{(i)}_{t-1}({\widetilde{\mathbf{x}}}^{(i)}_t,\xx{i}_t)}\label{eq:multi_del_2}\\
    &\le L\cdot\dia{\del{i}_{t-1}({\mathbf{x}}_t)}+2\cdot\dia{\mathcal{C}^{(i)}_{t-1}(\xx{i}_t)}+2L\|{\overbar{\mathbf{x}}}^{(i)}_t-\xx{i}_t\|_2+2L\|{\widetilde{\mathbf{x}}}^{(i)}_t-\xx{i}_t\|_2\label{eq:multi_del_3}\\
    &\le 4B\|\Gamma_{t-1}^{(i)}(\xx{i}_t,\xx{i}_t)\|_2^{1/2}+5L\cdot\dia{\del{i}_{t-1}({\mathbf{x}}_t)}, \label{eq:multi_del_4}
\end{align}
where:
\begin{itemize}
    \item \eqref{eq:multi_del_1} holds since there must exist ${\overbar{\mathbf{x}}}^{(i)}_t,{\widetilde{\mathbf{x}}}^{(i)}_t\in\del{i}_{t-1}({\mathbf{x}}_t)$ such that $\dia{\del{i+1}_{t-1}({\mathbf{x}}_t)}=\|{\overbar{\mathbf{z}}}-{\widetilde{\mathbf{z}}}\|_2$ for some ${\overbar{\mathbf{z}}}\in\overbar{\mathcal{C}}^{(i)}_{t-1}({\overbar{\mathbf{x}}}^{(i)}_t)$ and some ${\widetilde{\mathbf{z}}}\in\overbar{\mathcal{C}}^{(i)}_{t-1}({\widetilde{\mathbf{x}}}^{(i)}_t)$ (see \eqref{eq:del}), and the right hand side follows from the triangle inequality along with $\ff{i}({\overbar{\mathbf{x}}}^{(i)}_t)\in \overbar{\mathcal{C}}^{(i)}_{t-1}({\overbar{\mathbf{x}}}^{(i)}_t)$ and $\ff{i}({\widetilde{\mathbf{x}}}^{(i)}_t)\in \overbar{\mathcal{C}}^{(i)}_{t-1}({\widetilde{\mathbf{x}}}^{(i)}_t)$ (by \cref{lem:mul}). Here $\dia{\overbar{\mathcal{C}}^{(i)}_{t-1}(\cdot)}$ denotes the diameter of $\overbar{\mathcal{C}}^{(i)}_{t-1}(\cdot)$ (i.e., the Euclidean distance between the most distant pair of points in $\overbar{\mathcal{C}}^{(i)}_{t-1}(\cdot)$);
    \item \eqref{eq:multi_del_2} holds since $\ff{i}$ has Lipschitz constant $L$, and since ${\overbar{\mathbf{x}}}^{(i)}_t,{\widetilde{\mathbf{x}}}^{(i)}_t\in\del{i}_{t-1}({\mathbf{x}}_t)$ and $\overbar{\mathcal{C}}^{(i)}_{t-1}({\mathbf{z}})\subseteq\mathcal{C}^{(i)}_{t-1}({\mathbf{z}},\xx{i}_t) $ for all ${\mathbf{z}}\in\del{i}_{t-1}({\mathbf{x}}_t)$ by the definition in \eqref{eq:cbarz};
    \item \eqref{eq:multi_del_3} holds since $\dia{\mathcal{C}^{(i)}_{t-1}({\mathbf{z}},\xx{i}_t)}\le \dia{\mathcal{C}^{(i)}_{t-1}(\xx{i}_t)}+2L\|{\mathbf{z}}-\xx{i}_t\|_2$ by the definition in \eqref{eq:czz};
    \item \eqref{eq:multi_del_4} holds since $\dia{\mathcal{C}^{(i)}_{t-1}(\xx{i}_t)}=2B\|\Gamma_{t-1}^{(i)}(\xx{i}_t,\xx{i}_t)\|_2^{1/2}$ (see \eqref{eq:cz}) and  $\xx{i}_t,{\overbar{\mathbf{x}}}^{(i)}_t,{\widetilde{\mathbf{x}}}^{(i)}_t\in\del{i}_{t-1}({\mathbf{x}}_t)$.
\end{itemize}

Analogous to the UCB, we define $\mathrm{LCB}_t({\mathbf{x}})=\min_{{\mathbf{z}}\in\del{m}_t({\mathbf{x}})}\overbar{\mathcal{C}}_t^{(m)}({\mathbf{z}})$.  Moreover, we define
${\overbar{\mathbf{x}}}^{(m)}_t=\argmax_{{\mathbf{z}}\in\del{m}_{t-1}({\mathbf{x}}_t)}\big(\max \overbar{\mathcal{C}}^{(m)}_{t-1}({\mathbf{z}})\big)$, and ${\widetilde{\mathbf{x}}}^{(m)}_t=\argmin_{{\mathbf{z}}\in\del{m}_{t-1}({\mathbf{x}}_t)} \big(\min\overbar{\mathcal{C}}^{(m)}_{t-1}({\mathbf{z}})\big)$.  Then, we have
\begin{align}
    r_t={}&g({\mathbf{x}}^\ast)-g({\mathbf{x}}_t)\\
    \le{}& \mathrm{UCB}_{t-1}({\mathbf{x}}^\ast) - \mathrm{LCB}_{t-1}({\mathbf{x}}_t)\\
    \le{}& \mathrm{UCB}_{t-1}({\mathbf{x}}_t) - \mathrm{LCB}_{t-1}({\mathbf{x}}_t)\\
    ={}& \max\overbar{\mathcal{C}}^{(m)}_{t-1}({\overbar{\mathbf{x}}}^{(m)}_t) - \min \overbar{\mathcal{C}}^{(m)}_{t-1}({\widetilde{\mathbf{x}}}^{(m)}_t) \\
    \le{}& \Big(\muu{m}_{t-1}(\xx{m}_t)+B\|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}+L\|{\overbar{\mathbf{x}}}^{(m)}_t-\xx{m}_t\|_2\Big)\nonumber\\
    &-\Big(\muu{m}_{t-1}(\xx{m}_t)-B\|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}-L\|{\widetilde{\mathbf{x}}}^{(m)}_t-\xx{m}_t\|_2\Big)\label{eq:rt_ucb_1}\\
    \le{}& 2B \|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}+2L\cdot\dia{\del{m}_{t-1}({\mathbf{x}}_t)}\label{eq:rt_ucb_1a}\\
    \le{}& 2B \|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}+2L\Big( 4B \|\Gamma_{t-1}^{(m-1)}(\xx{m-1}_t,\xx{m-1}_t)\|_2^{1/2}+ 5L\cdot\dia{\del{m-1}_{t-1}({\mathbf{x}}_t)}\Big)\label{eq:rt_ucb_2}\\
    \le{}& 2B \|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}
    + \sum_{i=1}^{m-1} 2L \cdot (5L)^{m-1-i} \cdot 4B \|\Gamma_{t-1}^{(i)}(\xx{i}_t,\xx{i}_t)\|_2^{1/2}, \label{eq:rt_ucb_3}
\end{align}
where:
\begin{itemize}
    \item \eqref{eq:rt_ucb_1} holds since
    \begin{align}
        \max\overbar{\mathcal{C}}^{(m)}_{t-1}({\overbar{\mathbf{x}}}^{(m)}_t)
        &\le \max\mathcal{C}^{(m)}_{t-1}({\overbar{\mathbf{x}}}^{(m)}_t,\xx{m}_t)&&\text{(by \eqref{eq:cbarz})}\\
        &\le \max\mathcal{C}^{(m)}_{t-1}( \xx{m}_t) + L\|{\overbar{\mathbf{x}}}^{(m)}_t-\xx{m}_t\|_2 &&\text{(by \eqref{eq:czz})}\\
        &= \muu{m}_{t-1}(\xx{m}_t)+B\|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}+L\|{\overbar{\mathbf{x}}}^{(m)}_t-\xx{m}_t\|_2 &&\text{(by \eqref{eq:cz})}.
    \end{align}
    Similarly, we have
    \begin{align}
        \min\overbar{\mathcal{C}}^{(m)}_{t-1}({\widetilde{\mathbf{x}}}^{(m)}_t)
        &\ge \min\mathcal{C}^{(m)}_{t-1}({\widetilde{\mathbf{x}}}^{(m)}_t,\xx{m}_t)\\
        &\ge \min\mathcal{C}^{(m)}_{t-1}( \xx{m}_t) - L\|{\widetilde{\mathbf{x}}}^{(m)}_t-\xx{m}_t\|_2\\
        &= \muu{m}_{t-1}(\xx{m}_t)-B\|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}-L\|{\widetilde{\mathbf{x}}}^{(m)}_t-\xx{m}_t\|_2.
    \end{align}
    \item \eqref{eq:rt_ucb_1a} holds since all of the ${\mathbf{x}}$ vectors in \eqref{eq:rt_ucb_1a} lie in $\del{m}_{t-1}({\mathbf{x}}_t)$;
    \item \eqref{eq:rt_ucb_2} and \eqref{eq:rt_ucb_3} follow from the recursive relation in \eqref{eq:multi_del_4}.
\end{itemize}
Then, the cumulative regret is
\begin{align}
    R_T={}&\sum_{t=1}^T r_t\\
    \le{}& 2B \sum_{t=1}^T\|\Gamma_{t-1}^{(m)}(\xx{m}_t,\xx{m}_t)\|_2^{1/2}
    + 2L \cdot 4B \sum_{t=1}^T\|\Gamma_{t-1}^{(m-1)}(\xx{m-1}_t,\xx{m-1}_t)\|_2^{1/2}\nonumber\\
    {}&+ \sum_{i=1}^{m-2} 2L \cdot (5L)^{m-1-i} \cdot 4B \sum_{t=1}^T\|\Gamma_{t-1}^{(i)}(\xx{i}_t,\xx{i}_t)\|_2^{1/2} \\
    \le{}& 4B\frac{(5L)^m-1}{5L-1} \Sigma^\Gamma_T \\
    \le{}& 5^m BL^{m-1}\Sigma^\Gamma_T ,
\end{align}
where $\Sigma^\Gamma_T=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}} \sum_{t=1}^T\|\Gamma^{(i)}_{t-1}({\mathbf{z}}_t,{\mathbf{z}}_t)\|_2^{1/2}$.
\end{proof}

\subsection{Proof of \cref{thm:ucb_net} (Feed-Forward Networks)}
\label{sec:ucb_net}
In this section, we prove \cref{thm:ucb_net}, which is restated as follows.
\ucbnet*
\begin{proof}
Recall the UCB and LCB definitions in \eqref{eq:ucb_ij_bar}--\eqref{eq:ucb_ij_final}.
For a fixed input ${\mathbf{x}}\in\mathcal{X}$, defining ${\overbar{\mathbf{x}}}^{(i,j)}=\argmax_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\UCBBB{i,j}_t({\mathbf{z}})$ and ${\widetilde{\mathbf{x}}}^{(i,j)}=\argmin_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\LCBBB{i,j}_t({\mathbf{z}})$, the diameter of the confidence region of $\xx{i+1,j}$ is
\begin{align}
    \dia{\del{i+1,j}_t({\mathbf{x}})}&= \max_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\UCBBB{i,j}_t({\mathbf{z}}) - \min_{{\mathbf{z}}\in\del{i}_t({\mathbf{x}})}\LCBBB{i,j}_t({\mathbf{z}}) \\
    &= \UCBBB{i,j}_{t}({\overbar{\mathbf{x}}}^{(i,j)})-\LCBBB{i,j}_{t}({\widetilde{\mathbf{x}}}^{(i,j)}) \\
    &\le \UCBB{i,j}_{t}(\xx{i})+L\|{\overbar{\mathbf{x}}}^{(i,j)}-\xx{i}\|_2 -\LCBB{i,j}_{t}(\xx{i})+L\|\xx{i}-{\widetilde{\mathbf{x}}}^{(i,j)}\|_2\\
    &\le 2B\sig{i,j}_{t}(\xx{i})+2L\cdot\dia{\del{i}_{t}({\mathbf{x}})},
\end{align}
and therefore, the squared diameter of the confidence region of $\xx{i+1}=[\xx{i+1,1},\dots,\xx{i+1, d_{i+1}}]^T$ is
\begin{align}
    \dia{\del{i+1}_t({\mathbf{x}})}^2&=\sum_{j=1}^{d_{i+1}} \dia{\del{i+1,j}_t({\mathbf{x}})}^2\\
    &\le \sum_{j=1}^{d_{i+1}} \Big( 2B\sig{i,j}_{t}(\xx{i})+2L\cdot\dia{\del{i}_{t}({\mathbf{x}})} \Big)^2\\
    &= d_{i+1} \Big( 2B\sig{i,1}_{t}(\xx{i})+2L\cdot\dia{\del{i}_{t}({\mathbf{x}})} \Big)^2,
\end{align}
where the last step follows since $\sig{i,j}_t(\cdot)=\sig{i,1}_t(\cdot)$ for each $j\in[d_{i+1}]$ by the symmetry of our setup (each function is associated with the same kernel).

Then, by recursion, the diameter of the confidence region of $\xx{m}$ is
\begin{align}
    \dia{\del{m}_t({\mathbf{x}})}
    \le{}& 2\sqrt{d_m}B\sig{m-1,1}_t(\xx{m-1})+2\sqrt{d_m}L\cdot\dia{\del{m-1}_t({\mathbf{x}})}\\
    \le{}& 2\sqrt{d_m}B\sig{m-1,1}_t(\xx{m-1})+(2\sqrt{d_m}L)(2\sqrt{d_{m-1}}B)\sig{m-2,1}_t(\xx{m-2}) \nonumber \\
    & +(2\sqrt{d_m}L)(2\sqrt{d_{m-1}}L)\dia{\del{m-2}_t({\mathbf{x}})}\\
    \le{}& \sum_{i=1}^{m-1} \Big(\prod_{s=i+2}^m (2\sqrt{d_s}L) \Big) (2\sqrt{d_{i+1}}B) \sig{i,1}_t(\xx{i}),
\end{align}
where we set $\prod_{s=m+1}^m (2\sqrt{d_s}L)=1$.

Similarly to \eqref{eq:sim_0}--\eqref{eq:sim_6} for the case of chains, with ${\overbar{\mathbf{x}}}_t^{(m)}=\argmax_{{\mathbf{z}}\in\del{m}_{t-1}({\mathbf{x}}_t)}\UCBBB{m,1}_{t-1}({\mathbf{z}})$ and ${\widetilde{\mathbf{x}}}_t^{(m)}=\argmin_{{\mathbf{z}}\in\del{m}_{t-1}({\mathbf{x}}_t)}\LCBBB{m,1}_{t-1}({\mathbf{z}})$, the simple regret is upper bounded as follows:
\begin{align}
    r_t&=g({\mathbf{x}}^\ast)-g({\mathbf{x}}_t)\\
    &\le \mathrm{UCB}_{t-1}({\mathbf{x}}_t)-\LCBBB{m,1}_{t-1}({\widetilde{\mathbf{x}}}_t^{(m)})\\
    &= \UCBBB{m,1}_{t-1}({\overbar{\mathbf{x}}}_t^{(m)})-\LCBBB{m,1}_{t-1}({\widetilde{\mathbf{x}}}_t^{(m)})\\
    &\le \UCBB{m,1}_{t-1}(\xx{m}_t)+L\|{\overbar{\mathbf{x}}}_t^{(m)}-\xx{m}_t\|_2 -\LCBB{m,1}_{t-1}(\xx{m}_t)+L\|\xx{m}_t-{\widetilde{\mathbf{x}}}_t^{(m)}\|_2\\
    &\le 2B\sig{m,1}_{t-1}(\xx{m}_t)+2L\cdot\dia{\del{m}_{t-1}({\mathbf{x}}_t)}\\
    &\le 2B\sig{m,1}_{t-1}(\xx{m}_t)+ 2L\sum_{i=1}^{m-1} \Big(\prod_{s=i+2}^m (2\sqrt{d_s}L) \Big) (2\sqrt{d_{i+1}}B) \sig{i,1}_{t-1}(\xx{i}_t)\\
    &\le \sum_{i=1}^{m} 2^{m-i+1} B L^{m-i} \Big(\prod_{s=i+1}^{m} \sqrt{d_s} \Big)  \sig{i,1}_{t-1}(\xx{i}_t),
\end{align}
where we use the convention $\prod_{s=m+1}^m \sqrt{d_s}=1$.

Hence, the cumulative regret is
\begin{align}
    R_T&= \sum_{t=1}^T r_t\\
    &\le \sum_{i=1}^{m} 2^{m-i+1} B L^{m-i} \Big(\prod_{s=i+1}^{m} (\sqrt{d_s}) \Big) \sum_{t=1}^T  \sig{i,1}_{t-1}(\xx{i}_t) \\
    &\le 2B\frac{(2L)^m-1}{2L-1} \sqrt{D_{2,m}} \Sigma_T \\
    &\le 2^{m+1} B L^{m-1} \sqrt{D_{2,m}} \Sigma_T ,
\end{align}
where $D_{2,m}=\prod_{i=2}^{m} d_i$ and $\Sigma_T=\max\limits_{i\in[m]}\max\limits_{{\mathbf{z}}_1,\dots,{\mathbf{z}}_T\in\mathcal{X}^{(i)}}\sum_{t=1}^T \sig{i,1}_{t-1}({\mathbf{z}}_t)$.
\end{proof}

\section{Analysis of Non-Adaptive Sampling (\cref{algo:fd})}

As mentioned in the main body, the analysis in this appendix is restricted to the case that all domains in the network are \emph{hyperrectangular}.  This is a somewhat restrictive assumption, but we note that the primary reason for assuming this is to be able to apply Lemma \ref{lem:fill} below.  Hence, if Lemma \ref{lem:fill} can be generalized, then our results also generalize in the same way.

The lemma, stated as follows, provides an upper bound on the posterior standard deviation on a hyperrectangular domain for the Mat\'ern kernel in terms of the fill distance of the sampled points and the kernel parameter.
\begin{restatable}{lem}{fillmat}\textup{\citep[Theorem 5.4]{kanagawa2018gaussian}}
\label{lem:fill}
For the Mat\'ern kernel with smoothness $\nu$, a hyperrectangular domain $\mathcal{X}$, and any set of $T$ points $\mathcal{X}_T\subset \mathcal{X}$ with fill distance $\delta_T$, define
\begin{align}
	{\overbar{\sigma}}_T:=\max_{{\mathbf{x}}\in\mathcal{X}}\sigma_T({\mathbf{x}}),
\end{align}
where $\sigma_T(\cdot)$ is the posterior standard deviation computed based on $\mathcal{X}_T$ using $\eqref{eq:var}$. There exists a constant $\delta_0$ depending only on the kernel parameters such that if $\delta_T\le \delta_0$ then
\begin{align}
    {\overbar{\sigma}}_T = O\big( (\delta_T)^\nu \big).
\end{align}
\end{restatable}
Recalling that $\{\xx{i}_s\}_{s=1}^T$ denote the intermediate outputs of $\{{\mathbf{x}}_s\}_{s=1}^T$ right after $\ff{i-1}$, we define the fill distance of the $\mathcal{X}^{(i)}$ with respect to $\{\xx{i}_s\}_{s=1}^T$ as
\begin{align}
	\delta_T^{(i)}&=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\min_{s\in[T]}\|{\mathbf{z}}-\xx{i}_s\|_2.\label{eq:delta_i_T}
\end{align}
We also provide corollaries on the posterior standard deviation of a single layer for each network structure.
\begin{restatable}[Chains]{cor}{fillchain}
	\label{cor:fill_chain}
	With $k$ being the Mat\'ern kernel with smoothness $\nu$, consider a chain $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m]$.
	For any set of $T$ points $\mathcal{X}_T=\{{\mathbf{x}}_s\}_{s=1}^T\subset\mathcal{X}$ with fill distance $\delta_T$, let $\mathcal{X}^{(i)}_T=\{\xx{i}_s\}_{s=1}^T$ be the noise-free observations of $\ff{i-1}\circ\cdots\circ\ff{1}$.
	Define
	\begin{align}
		{\overbar{\sigma}}^{(i)}_T:=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\sig{i}_T({\mathbf{z}}),
	\end{align}
	where $\sig{i}_T(\cdot)$ is the posterior standard deviation computed based on $\mathcal{X}^{(i)}_T$ using $\eqref{eq:var}$. Then, for any $i\in[m]$, there exists a constant $\delta_0$ depending only on the kernel parameters such that if $\delta_T\le \delta_0$ then
	\begin{align}
		{\overbar{\sigma}}^{(i)}_T =  O\big((L^{i-1} \delta_T)^\nu\big).
	\end{align}
\end{restatable}
\begin{proof}
	First, it is straightforward that $\ff{i-1}\circ\cdots\circ\ff{1}$ has Lipschitz constant $L^{i-1}$. For any input ${\mathbf{x}}, {\mathbf{x}}_s\in\mathcal{X}$, $\ff{i-1}\circ\cdots\circ\ff{1}$ outputs $\xx{i},\xx{i}_s$ such that
	\begin{align}
		|\xx{i}-\xx{i}_s|=|\ff{i-1}\circ\cdots\circ\ff{1}({\mathbf{x}})-\ff{i-1}\circ\cdots\circ\ff{1}({\mathbf{x}}_s)|\le L^{i-1} \|{\mathbf{x}}-{\mathbf{x}}_s\|_2.
	\end{align}

	If $\{{\mathbf{x}}_s\}_{s=1}^T$ has fill distance $\delta_T$, then the fill distance of $\{\xx{i}_s\}_{s=1}^T$ is
	\begin{align}
		\delta_T^{(i)}&=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\min_{s\in[T]}|{\mathbf{z}}-\xx{i}_s|=\max_{{\mathbf{x}}\in\mathcal{X}}\min_{s\in[T]}|\xx{i}-\xx{i}_s|
		\le \max_{{\mathbf{x}}\in\mathcal{X}}\min_{s\in[T]} L^{i-1} \|{\mathbf{x}}-{\mathbf{x}}_s\|_2
		= L^{i-1} \delta_T.\label{eq:fill_i}
	\end{align}
	Then, by \cref{lem:fill} we have
	\begin{align}
		{\overbar{\sigma}}^{(i)}_T = \max_{{\mathbf{z}}\in\mathcal{X}^{(i)}} \sig{i}_T({\mathbf{z}}) = O\big((\delta_T^{(i)})^\nu\big) = O \big((L^{i-1} \delta_T)^\nu\big).
	\end{align}
\end{proof}

\begin{restatable}[Multi-output chains]{cor}{fillmul}
	\label{cor:fill_mul}
	For $\Gamma(\cdot,\cdot)=k(\cdot,\cdot){\mathbf{I}}_n$, with $k$ being the Mat\'ern kernel with smoothness $\nu$ and ${\mathbf{I}}_n$ being the identity matrix of size $n$, consider a chain $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}\in\mathcal{H}_\Gamma(B)\cap\mathcal{F}(L)$ and $\mathcal{X}^{(i)}$ being a hyperrectangle for each $i\in[m]$ and $n=d_{i+1}$.
	For any set of $T$ points $\mathcal{X}_T=\{{\mathbf{x}}_s\}_{s=1}^T\subset\mathcal{X}$ with fill distance $\delta_T$, let $\mathcal{X}^{(i)}_T=\{\xx{i}_s\}_{s=1}^T$ be the noise-free observations of $\ff{i-1}\circ\cdots\circ\ff{1}$.
	Define
	\begin{align}
		\overline{\Gamma}^{(i)}_T:=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\|\Gamma^{(i)}_T({\mathbf{z}},{\mathbf{z}})\|_2^{1/2},
	\end{align}
	where $\Gamma^{(i)}_T({\mathbf{z}},{\mathbf{z}})$ is the posterior variance matrix computed based on $\mathcal{X}^{(i)}_T$ using $\eqref{eq:multi_var}$, and $\|\cdot\|_2$ denotes the matrix spectral norm. Then, for any $i\in[m]$, there exists a constant $\delta_0$ depending only on the kernel parameters such that if $\delta_T\le \delta_0$ then
	\begin{align}
		\overline{\Gamma}^{(i)}_T =  O\big((L^{i-1} \delta_T)^\nu\big).
	\end{align}
\end{restatable}
\begin{proof} Recalling that $\|\cdot\|_2$ with a matrix argument denotes the spectral norm, we have
	\begin{align}
		\overbar{\Gamma}^{(i)}_T&=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\|\Gamma^{(i)}_T({\mathbf{z}},{\mathbf{z}})\|_2^{1/2}\\
		&=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}} \|\big(\sig{i}_T({\mathbf{z}})\big)^2 {\mathbf{I}}_n \|_2^{1/2} && \text{(by \cref{lem:multi_var})} \\
		&=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}} \sig{i}_T({\mathbf{z}}) && \text{(by $\|{\mathbf{I}}_n\|_2=1$)} \\
		&=O\big((\delta_T^{(i)})^\nu\big) && \text{(by \cref{lem:fill})}\label{eq:fill_mul}\\
		&= O \big((L^{i-1} \delta_T)^\nu\big). &&\text{(by \eqref{eq:fill_i})}
	\end{align}
\end{proof}

\begin{restatable}[Feed-forward networks]{cor}{fillnet}
	\label{cor:fill_net}
	With $k$ being the Mat\'ern kernel with smoothness $\nu$, consider a feed-forward network $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$
	with $\ff{s}(\cdot)=[\ff{s,j}(\cdot)]_{j=1}^{d_{s+1}}$ and $\ff{s,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ and $\mathcal{X}^{(s)}$ being a hyperrectangle for each $s\in[m],j\in[d_{s+1}]$.
	For any set of $T$ points $\mathcal{X}_T=\{{\mathbf{x}}_s\}_{s=1}^T\subset\mathcal{X}$ with fill distance $\delta_T$, let $\mathcal{X}^{(i)}_T=\{\xx{i}_s\}_{s=1}^t$ be the noise-free observations of $\ff{i-1}\circ\cdots\circ\ff{1}$.
	Define
	\begin{align}
		{\overbar{\sigma}}^{(i,j)}_T:=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\sig{i,j}_T({\mathbf{z}}),
	\end{align}
	where $\sig{i,j}_T(\cdot)$ is the posterior standard deviation computed based on $\mathcal{X}^{(i)}_T$ using $\eqref{eq:var}$. Then, for each $i\in[m]$ and $j\in[d_{i+1}]$, there exists a constant $\delta_0$ depending only on the kernel parameters such that if $\delta_T\le \delta_0$ then
	\begin{align}
		{\overbar{\sigma}}^{(i,j)}_T =  O\big((\sqrt{D_{2,i}}L^{i-1} \delta_T)^\nu\big), \label{eq:sigbar}
	\end{align}
	where $D_{2,i}=\prod_{s=2}^i d_s$.
\end{restatable}
\begin{proof}
	Since $\ff{s,j}\in\mathcal{F}(L)$ for each $s\in[m]$, we have for any ${\mathbf{z}},{\mathbf{z}}'\in\mathcal{X}^{(s)}$ that
	\begin{align}
		\|\ff{s}({\mathbf{z}})-\ff{s}({\mathbf{z}}')\|_2
		&=\sqrt{\sum_{j=1}^{d_{s+1}}\|\ff{s,j}({\mathbf{z}})-\ff{s,j}({\mathbf{z}}')\|_2^2}\\
		&\le \sqrt{\sum_{j=1}^{d_{s+1}}L^2\|{\mathbf{z}}-{\mathbf{z}}'\|_2^2}\\
		&= \sqrt{d_{s+1}}L\|{\mathbf{z}}-{\mathbf{z}}'\|_2,
	\end{align}
	and therefore $\ff{i-1}\circ\cdots\circ\ff{1}$ is Lipschitz continuous with constant $\sqrt{D_{2,i}}L^{i-1}$, where $D_{2,i}=\prod_{s=1}^i d_s$.
	If $\{{\mathbf{x}}_s\}_{s=1}^T$ has fill distance $\delta_T$, then the fill distance of $\{\xx{i}_s\}_{s=1}^T$ is
	\begin{align}
		\delta_T^{(i)}&=\max_{{\mathbf{z}}\in\mathcal{X}^{(i)}}\min_{s\in[T]}\|{\mathbf{z}}-\xx{i}_s\|_2\\
		&=\max_{{\mathbf{x}}\in\mathcal{X}}\min_{s\in[T]}\|\xx{i}-\xx{i}_s\|_2\\
		&\le \max_{{\mathbf{x}}\in\mathcal{X}}\min_{s\in[T]}\sqrt{D_{2,i}} L^{i-1} \|{\mathbf{x}}-{\mathbf{x}}_s\|_2\\
		&\le \sqrt{D_{2,i}} L^{i-1} \delta_T.\label{eq:fill_i_net}
	\end{align}
	Hence, by \cref{lem:fill} we have
	\begin{align}
		{\overbar{\sigma}}^{(i,j)}_T = \max_{{\mathbf{z}}\in\mathcal{X}^{(i)}} \sig{i,j}_T({\mathbf{z}}) = O\big((\delta_T^{(i)})^\nu\big) = O \big((\sqrt{D_{2,i}}L^{i-1} \delta_T)^\nu\big).
	\end{align}
\end{proof}

\begin{rem}
	\cref{cor:fill_mul} and \cref{cor:fill_net} require $\mathcal{X}^{(i)}$ being a hyperrectangle. This is because they are the consequences of \cref{lem:fill}, which only holds for a hyperrectangular domain. For chains, the requirement of hyperrectangular domain is trivial, since $\mathcal{X}=[0,1]^d$ and $\mathcal{X}^{(2)},\dots,\mathcal{X}^{(m)}$ have single dimension.
\end{rem}

\subsection{Proof of \cref{thm:fd_chain} (Chains)}
\label{sec:fd_chain}
In this section, we prove \cref{thm:fd_chain}, which is restated as follows.
\fdchain*

\begin{proof}
Since $\ff{i}\in\mathcal{F}(L)$ implies that $\ff{m}\circ\cdots\circ\ff{i}$ is Lipschitz continuous with constant $L^{m-i+1}$, defining
\begin{align}
	{\widetilde{\sigma}}^{(i)}_T({\mathbf{x}}):=(\sig{i}_T\circ\muu{i-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}),
\end{align}
we have for all ${\mathbf{x}}\in\mathcal{X}$ that
\begin{align}
	g({\mathbf{x}})={}&(\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{2})\big(\ff{1}({\mathbf{x}})\big)\label{eq:fd_start}\\
	\le{}& (\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{2})\big(\muu{1}_T({\mathbf{x}})\big) + L^{m-1} B\sig{1}_T({\mathbf{x}})\\
	={}& (\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{3})\big((\ff{2}\circ\muu{1}_T)({\mathbf{x}})\big) + L^{m-1} B\sig{1}_T({\mathbf{x}})\\
	\le{}& (\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{3})\big((\muu{2}_T\circ\muu{1}_T)({\mathbf{x}})\big) + L^{m-2}B\sig{2}_T\big(\muu{1}_T({\mathbf{x}})\big) + L^{m-1} B\sig{1}_T({\mathbf{x}})\\

	\le{}& (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})+B\sig{m}_T\big((\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})\big)\nonumber\\
	&+LB\sig{m-1}_T\big((\muu{m-2}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})\big)+\cdots+ L^{m-1} B\sig{1}_T({\mathbf{x}})\\
	={}& (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}) + B\sum_{i=1}^m L^{m-i} (\sig{i}_T\circ\muu{i-1}_T\circ\cdots\circ\muu{1})({\mathbf{x}})\\
	={}& (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}) + B\sum_{i=1}^m L^{m-i} {\widetilde{\sigma}}_T^{(i)}({\mathbf{x}}),\label{eq:fd_1}
\end{align}
where the first two inequalities use the confidence bounds and Lipschitz assumption, and the third inequality follows by continuing recursively.  Similarly, we also have
\begin{align}
	g({\mathbf{x}})\ge (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}) - B\sum_{i=1}^m L^{m-i} {\widetilde{\sigma}}_T^{(i)}({\mathbf{x}}).
\end{align}
Hence, defining ${\widetilde{\sigma}}_T^{(i)}:=\max_{{\mathbf{x}}\in\mathcal{X}}{\widetilde{\sigma}}_T^{(i)}({\mathbf{x}})$, it follows that
\begin{align}
	r^\ast_T={}&g({\mathbf{x}}^\ast)-g({\mathbf{x}}^\ast_T)\\
	\le{}& \Big( (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}^\ast) + B\sum_{i=1}^m L^{m-i} {\widetilde{\sigma}}_T^{(i)} ({\mathbf{x}}^\ast) \Big)\nonumber\\
	&-\Big( (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}^\ast_T) - B\sum_{i=1}^m L^{m-i} {\widetilde{\sigma}}_T^{(i)} ({\mathbf{x}}^\ast_T) \Big) \\

	\le{}& 2B \sum_{i=1}^m L^{m-i} \max\{{\widetilde{\sigma}}_T^{(i)} ({\mathbf{x}}^\ast), {\widetilde{\sigma}}_T^{(i)} ({\mathbf{x}}^\ast_T)\}\label{eq:fd_end} \\
	\le{}& 2B \sum_{i=1}^m L^{m-i} {\widetilde{\sigma}}_T^{(i)},
\end{align}
where \eqref{eq:fd_end} uses the fact that ${\mathbf{x}}^\ast_T$ maximizes the composite posterior mean function.  It remains to derive an upper bound on ${\widetilde{\sigma}}_T^{(i)}$ for $i\in[m]$. Firstly, by \cref{lem:fill}, we have
\begin{align}
	{\widetilde{\sigma}}^{(1)}_T({\mathbf{x}})=\sig{1}_T({\mathbf{x}})=O\big((\delta_T)^\nu\big)=O(T^{-\nu/d}).\label{eq:sigtil_1}
\end{align}
We now represent the upper bound on ${\widetilde{\sigma}}^{(i+1)}_T$ in terms of $\{{\widetilde{\sigma}}^{(j)}_T:j\in[i]\}$ for $i\in[m-1]$. For all ${\mathbf{x}}\in\mathcal{X}$, the distance between $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})$ and its closest point in $\mathcal{X}^{(i+1)}$ is
\begin{align}
	\delta_{i+1}({\mathbf{x}})&=\min_{{\mathbf{z}}\in\mathcal{X}^{(i+1)}}|(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})-{\mathbf{z}}|\\
	&\le|\muu{i}_T\circ\cdots\circ\muu{1}_T({\mathbf{x}})-\ff{i}\circ\cdots\circ\ff{1}({\mathbf{x}})|\\
	&= B\sum_{j=1}^{i}L^{i-j}{\widetilde{\sigma}}^{(j)}_T({\mathbf{x}}), \label{eq:delta_i_1}
\end{align}
where the last step follows from \eqref{eq:fd_1}.
Then, recalling the definition of $\delta^{(i)}_T$ in \eqref{eq:delta_i_T}, we have that the distance between $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})$ and its closest point in $\mathcal{X}^{(i+1)}_T$ is
\begin{align}
	\delta'_{i+1}({\mathbf{x}})&=\min_{s\in[T]}|(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})-\xx{i+1}_s|\\
	&\le\delta^{(i+1)}_T+\delta_{i+1}({\mathbf{x}}) \label{eq:delta_prime_2}\\
	&\le O(L^i T^{-1/d})+B\sum_{j=1}^{i}L^{i-j}{\widetilde{\sigma}}^{(j)}_T({\mathbf{x}}) \label{eq:delta_prime_3}\\
	&= O\Big(\max\Big\{L^iT^{-1/d}, i B \cdot \max_{j\in[i]}L^{i-j} {\widetilde{\sigma}}^{(j)}_T({\mathbf{x}}) \Big\}\Big),
\end{align}
where \eqref{eq:delta_prime_2} applies the triangle inequality, and \eqref{eq:delta_prime_3} uses \eqref{eq:fill_i} and \eqref{eq:delta_i_1}.

Now, we extend $\mathcal{X}^{(i+1)}$ to $\widetilde{\mathcal{X}}^{(i+1)}$, the shortest interval that covers both $\{\muu{i}_T\circ\cdots\circ\muu{1}_T({\mathbf{x}}):{\mathbf{x}}\in\mathcal{X}\}$ and $\mathcal{X}^{(i+1)}$. Then, the fill distance of the extended domain with respect to $\mathcal{X}^{(i+1)}_T$ is
\begin{align}
	\delta'_{i+1}:=\max_{{\mathbf{x}}\in\mathcal{X}}\delta'_{i+1}({\mathbf{x}})=O\Big(\max\Big\{L^iT^{-1/d},  i B \cdot\max_{j\in[i]}L^{i-j} {\widetilde{\sigma}}^{(j)}_T \Big\}\Big),
\end{align}
and by \cref{lem:fill}, we have
\begin{align}
	{\widetilde{\sigma}}^{(i+1)}_T=\max_{{\mathbf{x}}\in\mathcal{X}}{\widetilde{\sigma}}^{(i+1)}_T({\mathbf{x}})=O\big((\delta'_{i+1})^\nu\big)=O\Bigg(\Big(\max\Big\{L^iT^{-1/d},  i B\cdot\max_{j\in[i]}L^{i-j} {\widetilde{\sigma}}^{(j)}_T \Big\}\Big)^\nu\Bigg).\label{eq:sigtil}
\end{align}
We now split the analysis into two cases.

\underline{The case that $\nu \le 1$}:
Recall that we assume $B=\Theta(L)$, and using the fact that ${\widetilde{O}}(L^{c\nu})$ is equivalent to $O\big( (L^{c\nu} {\rm poly}(c\log L)) \big)$ for fixed $\nu$, it follows from \eqref{eq:sigtil_1} and \eqref{eq:sigtil} that when $\nu \le 1$, we have
\begin{align}
	{\widetilde{\sigma}}_T^{(2)}&=O\Big(\max\Big\{L^\nu T^{-\nu/d},B^\nu T^{-\nu^2/d}\Big\}\Big)=O(L^\nu T^{-\nu^2/d}), \\
	{\widetilde{\sigma}}^{(3)}_T&=O\Big(\max\Big\{L^{2\nu}T^{-\nu/d}, 2^\nu B^\nu L^\nu T^{-\nu^2/d}, 2^\nu B^\nu L^{\nu^2} T^{-\nu^3/d}\Big\}\Big)=O(2^\nu B^\nu L^\nu T^{-\nu^3/d})={\widetilde{O}}( B^\nu L^\nu T^{-\nu^3/d}),\\
	&\vdotswithin{=}\nonumber\\
	{\widetilde{\sigma}}^{(i)}_T&={\widetilde{O}}\big((i-1)^\nu B^\nu L^{(i-2)\nu}T^{-\nu^i/d}\big)={\widetilde{O}}(B^\nu L^{(i-2)\nu}T^{-\nu^i/d}),
\end{align}
and the simple regret is
\begin{align}
	r^\ast_T&\le 2B\sum_{i=1}^m L^{m-i}{\widetilde{\sigma}}^{(i)}_T={\widetilde{O}}(BL^{m-1}T^{-\nu^m/d}).
\end{align}

\underline{The case that $\nu \le 1$}:
With $B=\Theta(L)$, using the property of ${\widetilde{O}}(\cdot)$ similarly to the first case, it follows from \eqref{eq:sigtil_1} and \eqref{eq:sigtil} that when $\nu > 1$, we have
\begin{align}
	{\widetilde{\sigma}}_T^{(2)}&=O\Big(\max\Big\{L^\nu T^{-\nu/d},B^\nu T^{-\nu^2/d}\Big\}\Big)=O(L^\nu T^{-\nu/d}), \\
	{\widetilde{\sigma}}^{(3)}_T&=O\Big(\max\Big\{L^{2\nu}T^{-\nu/d}, 2^\nu B^\nu L^{\nu^2} T^{-\nu^2/d}\Big\}\Big)=O(2^\nu B^\nu L^{\nu^2} T^{-\nu/d})={\widetilde{O}}(B^\nu L^{\nu^2} T^{-\nu/d}),\\
	&\vdotswithin{=}\nonumber\\
	{\widetilde{\sigma}}^{(i)}_T&={\widetilde{O}}\Big(\max\Big\{L^{(i-1)\nu}T^{-\nu/d}, (i-1)^\nu B^{\nu+\cdots+\nu^{i-2}}L^{\nu^{i-1}}T^{-\nu^2/d}\Big\}\Big)\nonumber\\
	&= {\widetilde{O}}\Big(\max\Big\{L^{(i-1)\nu}T^{-\nu/d}, B^{\nu+\cdots+\nu^{i-2}}L^{\nu^{i-1}}T^{-\nu^2/d}\Big\}\Big),
\end{align}
and the simple regret is
\begin{align}
	r^\ast_T&\le 2B\sum_{i=1}^m L^{m-i}{\widetilde{\sigma}}^{(i)}_T
	={\widetilde{O}}\Big(\max\Big\{BL^{(m-1)\nu}T^{-\nu/d}, B^{1+\nu+\nu^2+\cdots+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu^2/d} \Big\}\Big).
\end{align}
\end{proof}

\subsection{Two More Restrictive Cases for Chains}
\label{sec:rest_chain}
Recall the following two more restrictive cases introduced in the main body:
\begin{itemize}
	\item \textbf{Case 1}: We additionally assume that $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}^\ast)\in\mathcal{X}^{(i+1)}$ and $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}^\ast_T)\in \mathcal{X}^{(i+1)}$ for $i\in[m-1]$.
	\item \textbf{Case 2}: We additionally assume that all the $\mathcal{X}^{(i)}$'s are known. Defining
	\begin{align}
		\muutil{i}_T({\mathbf{z}}) &= \argmin_{{\mathbf{z}}'\in\mathcal{X}^{(i+1)}}|\muu{i}_T({\mathbf{z}})-{\mathbf{z}}'|, \label{eq:mutil}
	\end{align}
	we let the algorithm return
	\begin{align}
		{\mathbf{x}}^\ast_T=\argmax_{{\mathbf{x}}\in\mathcal{X}}(\muutil{m}_T\circ\cdots\circ\muutil{1}_T)({\mathbf{x}}).
	\end{align}
\end{itemize}
In Case 1, it follows from \cref{cor:fill_chain} (along with \eqref{eq:sigbar} and $\delta_T= \Theta(T^{-\frac{1}{d}})$) that for each ${\mathbf{x}}'\in\{{\mathbf{x}}^\ast,{\mathbf{x}}^\ast_T\}$
\begin{align}
	{\widetilde{\sigma}}^{(i+1)}_T({\mathbf{x}}')&=(\sig{i+1}_T\circ\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}')\le{\overbar{\sigma}}^{(i+1)}_T=O\big((L^i \delta_T)^\nu\big)=O\big((L^i T^{-1/d})^\nu\big).
\end{align}
By substituting these upper bounds into \eqref{eq:fd_end}, it holds that
\begin{align}
	r^\ast_T&\le 2B \sum_{i=1}^m L^{m-i} \max\{{\widetilde{\sigma}}_T^{(i)} ({\mathbf{x}}^\ast), {\widetilde{\sigma}}_T^{(i)} ({\mathbf{x}}^\ast_T)\} \\
	&\le 2B \sum_{i=1}^m L^{m-i}  O\big((L^i T^{-1/d})^\nu\big)\\
	&=\begin{cases}
		O(BL^{m-1}T^{-\nu/d})&\text{ when $\nu\le 1$,}\\
		O(BL^{(m-1)\nu}T^{-\nu/d})&\text{ when $\nu> 1$.}
	\end{cases}
\end{align}
In Case 2, for all ${\mathbf{z}}\in\mathcal{X}^{(i)}$, we have $\muutil{i}_T({\mathbf{z}})\in\mathcal{X}^{(i+1)}$, and
\begin{align}
	|\ff{i}({\mathbf{z}})-\muutil{i}_T({\mathbf{z}})|&\le  |\ff{i}({\mathbf{z}})-\muu{i}_T({\mathbf{z}})| + |\muu{i}_T({\mathbf{z}})-\muutil{i}_T({\mathbf{z}})|\\
	&\le |\ff{i}({\mathbf{z}})-\muu{i}_T({\mathbf{z}})| + |\muu{i}_T({\mathbf{z}})-\ff{i}({\mathbf{z}})| \\
	&\le 2B\sig{i}_T({\mathbf{z}}),
\end{align}
where we used \eqref{eq:mutil} and the confidence bounds.
By recursion, we also have $(\muutil{i}_T\circ\cdots\circ\muutil{1}_T)({\mathbf{x}})\in\mathcal{X}^{(i+1)}$ for all ${\mathbf{x}}\in\mathcal{X}$, and therefore
\begin{align}
	(\sig{i+1}_T\circ\muutil{i}_T\circ\cdots\circ\muutil{1}_T)({\mathbf{x}})\le{\overbar{\sigma}}^{(i+1)}_T.
\end{align}

Hence, replacing $\mu$ with ${\widetilde{\mu}}$ in \eqref{eq:fd_start}--\eqref{eq:fd_end}, it follows from \cref{cor:fill_chain} that
\begin{align}
	r^\ast_T\le 4B\sum_{i=1}^m L^{m-i}{\overbar{\sigma}}_T^{(i)}=
	\begin{cases}
		O(BL^{m-1}T^{-\nu/d})&\text{ when $\nu\le 1$,}\\
		O(BL^{(m-1)\nu}T^{-\nu/d})&\text{ when $\nu> 1$.}
	\end{cases}
\end{align}

\subsection{Non-Adaptive Sampling Method for Multi-Output Chains}
\label{sec:fd_mul}
For multi-output chains, the composite posterior mean of $g({\mathbf{x}})$ is
\begin{align}
	\mu_T^g({\mathbf{x}})=(\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}),\label{eq:mean_mul}
\end{align}
where $\muu{i}_T$ denotes the posterior mean of $\ff{i}$ computed using \eqref{eq:multi_mean} based on $\{(\xx{i}_s,\xx{i+1}_s)\}_{s=1}^T$ for each $i\in[m]$.

We again assume that each $\mathcal{X}^{(i)}$ is a hyperrectangle of dimension $d_i$.  Then, the upper bound on simple regret of \cref{algo:fd} using \eqref{eq:mean_mul} is provided in the following theorem.

\begin{restatable}[Non-adaptive sampling method for multi-output chains]{thm}{fdmul}
	\label{thm:fd_mul}
	Under the setup of \cref{sec:setup}, given $B=\Theta(L)$, $k=k_\text{Mat\'ern}$, $\Gamma(\cdot,\cdot)=k(\cdot,\cdot){\mathbf{I}}$, and a multi-output chain  $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}\in\mathcal{H}_\Gamma(B)\cap\mathcal{F}(L)$ and $\mathcal{X}^{(i)}$ being a hyperrectangle for each $i\in[m]$,
	\begin{itemize}
		\item when $\nu\le 1$, \cref{algo:fd} achieves
		\begin{align*}
			r^\ast_T={\widetilde{O}}(BL^{m-1}T^{-\nu^m/d});
		\end{align*}
		\item when $\nu>1$, \cref{algo:fd} achieves
		\begin{align*}
			r^\ast_T={\widetilde{O}}\big(\max\big\{BL^{(m-1)\nu}T^{-\nu/d}, B^{1+\nu+\nu^2+\cdots+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu^2/d} \big\}\big).
		\end{align*}
	\end{itemize}
\end{restatable}

\begin{proof}
The analysis is similar to the case of chains, so we omit some details and focus on the main differences.
Defining $\widetilde{\Gamma}_T^{(i)}({\mathbf{x}})=\|\Gamma_T^{(i)}\circ\muu{i-1}_T\circ\cdots\circ\muu{1}_T({\mathbf{x}})\|_2^{1/2}$, similarly to the case of chains, we have
\begin{align}
	|g({\mathbf{x}}) - (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})|\le B\sum_{i=1}^m L^{m-i} \widetilde{\Gamma}_T^{(i)}({\mathbf{x}}),
\end{align}
and
\begin{align}
	r^\ast_T=g({\mathbf{x}}^\ast)-g({\mathbf{x}}^\ast_T)&\le 2 B\sum_{i=1}^m L^{m-i} \max\{ \widetilde{\Gamma}_T^{(i)}({\mathbf{x}}^\ast),\widetilde{\Gamma}_T^{(i)}({\mathbf{x}}^\ast_T) \}.
\end{align}
For each ${\mathbf{x}}'\in\{{\mathbf{x}}^\ast,{\mathbf{x}}^\ast_T\}$, we extend $\mathcal{X}^{(i+1)}$ to the smallest hyperrectangle $\widetilde{\mathcal{X}}^{(i+1)}$ that covers $(\muu{i}_T\circ\cdots\circ\muu{1})({\mathbf{x}}')$, and the original domain $\mathcal{X}^{(i+1)}$ (see \cref{fig:mul_fill}). Recalling the definition of $\delta^{(i+1)}_T$ in \eqref{eq:delta_i_T}, with $\delta_{i+1}({\mathbf{x}})=\min_{{\mathbf{z}}\in\mathcal{X}^{(i+1)}}\|(\muu{i}_T\circ\cdots\circ\muu{1})({\mathbf{x}})-{\mathbf{z}}\|_2$, the fill distance of $\widetilde{\mathcal{X}}^{(i+1)}$ with respect to $\{\xx{i+1}_s\}_{s=1}^T$ is
\begin{align}
	\delta'_{i+1}&=\max_{{\mathbf{z}}'\in\widetilde{\mathcal{X}}^{(i+1)}}\min_{s\in[T]}\|{\mathbf{z}}'-\xx{i+1}_s\|_2\\
	&\le \delta_T^{(i+1)}+\max_{{\mathbf{z}}'\in\widetilde{\mathcal{X}}^{(i+1)}}\min_{{\mathbf{z}}\in\mathcal{X}^{(i+1)}}\|{\mathbf{z}}'-{\mathbf{z}}\|_2\\
	&\le \delta_T^{(i+1)}+\delta_{i+1}({\mathbf{x}}')\\
	&\le \delta_T^{(i+1)}+B\sum_{j=1}^i L^{i-j} \widetilde{\Gamma}_T^{(i)}({\mathbf{x}}').
\end{align}
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.6\linewidth]{figs/mul_fill.pdf}
	\caption{Extended domain for multi-output chains.}
	\label{fig:mul_fill}
\end{figure}
Since $\widetilde{\mathcal{X}}^{(i+1)}$ is a hyperrectangle and $\muu{i}_T\circ\cdots\circ\muu{1}({\mathbf{x}}')\in\widetilde{\mathcal{X}}^{(i+1)}$, it follows from \eqref{eq:fill_mul} that
\begin{align}
	\widetilde{\Gamma}_T^{(i+1)}({\mathbf{x}}')&=O\big((\delta'_{i+1})^\nu\big)\\
	&=O\Bigg(\Big(\max\Big\{L^iT^{-1/d},  i B\cdot\max_{j\in[i]}L^{i-j} \widetilde{\Gamma}_T^{(j)} ({\mathbf{x}}') \Big\}\Big)^\nu\Bigg).
\end{align}
By \cref{lem:multi_var}, we also have
\begin{align}
	\widetilde{\Gamma}_T^{(1)}({\mathbf{x}}')=\|\Gamma_T^{(1)}({\mathbf{x}}')\|_2^{1/2}=\sigma_T^{(1)}({\mathbf{x}}')=O\big((\delta_T)^\nu)=O(T^{-\nu/d}).
\end{align}

This recursive relation is exactly the same as that of chains in \eqref{eq:sigtil}, showing \cref{thm:fd_chain} extends to multi-output chains.

For the two more restrictive cases, \cref{cor:fill_mul}, with the same posterior standard deviation upper bound as \cref{cor:fill_chain}, implies that the results in \cref{sec:rest_chain} also extend to multi-output chains.

\end{proof}

\subsection{Non-Adaptive Sampling Method for Feed-Forward Networks}
\label{sec:fd_net}
For feed-forward networks of scalar-valued functions, let $\muu{i,j}_T$ denote the posterior mean of $\ff{i,j}$ computed using \eqref{eq:mean} based on $\{(\xx{i}_s,\xx{i+1,j}_s)\}_{s=1}^T$. The composite posterior mean of $g({\mathbf{x}})$ with feed-forward network structure is
\begin{align}
	\mu_T^g({\mathbf{x}})=\muu{m,1}_T(\zz{m})\label{eq:mean_net}
\end{align}
with
\begin{align}
	\zz{1}&={\mathbf{x}},\\
	\zz{i+1,j}&=\muu{i,j}_T(\zz{i})&&\text{ for }i\in[m-1],j\in[d_{i+1}],\\
	\zz{i+1}&=\muu{i}_T(\zz{i})=[\zz{i+1,1},\dots,\zz{i+1,d_i}]&&\text{ for }i\in[m-1].
\end{align}

Aussiming each $\mathcal{X}^{(i)}$ is a hyperrectangle of dimension $d_i$, the following theorem provides the upper bound on simple regret of \cref{algo:fd} using \eqref{eq:mean_net} for feed-forward networks.

\begin{restatable}[Non-adaptive sampling method for feed-forward networks]{thm}{fdnet}
	\label{thm:fd_net}
	Under the setup of \cref{sec:setup}, given $B=\Theta(L)$,  $k=k_\text{Mat\'ern}$, and  a feed-forward network $g=\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{1}$ with $\ff{i}({\mathbf{z}})=[\ff{i,j}({\mathbf{z}})]_{j=1}^{d_{i+1}}$, $\ff{i,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$, and $\mathcal{X}^{(i)}$ being a hyperrectangle for each $i\in[m],j\in[d_{i+1}]$,
	\begin{itemize}
		\item when $\nu\le 1$, \cref{algo:fd} achieves
		\begin{align}
			r^\ast_T={\widetilde{O}}(\sqrt{D_{2,m}}BL^{m-1}T^{-\nu^m/d});
		\end{align}
		\item when $\nu>1$, \cref{algo:fd} achieves
		\begin{align}
			r^\ast_T={\widetilde{O}}\Big(\max\Big\{(D_{2,m})^{\nu/2}BL^{(m-1)\nu}T^{-\nu/d}, {\widetilde{D}}_{2,m}^\nu B^{1+\nu+\nu^2+\cdots+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu^2/d} \Big\}\Big),
		\end{align}
	\end{itemize}
where $D_{2,m}=\prod_{i=2}^m d_i$ and ${\widetilde{D}}_{2,m}^\nu=\prod_{i=2}^m (d_i)^{\nu^{m+1-i}/2}$.

\end{restatable}

\begin{proof}
With $\muu{i}_T(\cdot)=[\muu{i,j}_T(\cdot)]_{j=1}^{d_{i+1}}$, we have
\begin{align}
    \|\ff{i}({\mathbf{z}})-\muu{i}_T({\mathbf{z}})\|_2
    &=\sqrt{\sum_{j=1}^{d_{i+1}}\big(\ff{i,j}({\mathbf{z}})-\muu{i,j}_T({\mathbf{z}})\big)^2}\\
    &\le\sqrt{\sum_{j=1}^{d_{i+1}}\big(B\sig{i,j}_T({\mathbf{z}})\big)^2}\\
    &=\sqrt{d_{i+1}}B\sig{i,1}_T({\mathbf{z}}),
\end{align}
where we use the confidence bounds and the fact that $\sig{i,j}_T({\mathbf{z}}) = \sig{i,1}_T({\mathbf{z}})$ by the symmetry of our setup.

Since each $\ff{i,j}$ is Lipschitz continuous with parameter $L$, we have that $\ff{i}$ is Lipschitz continuous with parameter $\sqrt{d_{i+1}}L$, and $\ff{m}\circ\cdots\circ\ff{i}$ is Lipschitz continuous with parameter $\sqrt{D_{i+1,m+1}}L^{m-i+1}$, where $D_{i,j}=\prod_{s=i}^j d_s$ and $d_{m+1}=1$.
Then, defining ${\widetilde{\sigma}}_T^{(i,1)}({\mathbf{x}})=(\sig{i,1}_T\circ\muu{i-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})$, we can follow \eqref{eq:fd_start}--\eqref{eq:fd_1} to obtain
\begin{align}
    g({\mathbf{x}})={}&(\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{2})\big(\ff{1}({\mathbf{x}})\big)\\
    \le{}& (\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{2})\big(\muu{1}_T({\mathbf{x}})\big) + \sqrt{D_{3,m+1}} L^{m-1} \|\ff{1}({\mathbf{x}})-\muu{1}_T({\mathbf{x}})\|_2\\
    ={}& (\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{2})\big(\muu{1}_T({\mathbf{x}})\big) + \sqrt{D_{2,m+1}} L^{m-1} B\sig{1,1}_T({\mathbf{x}})\\
    \le{}& (\ff{m}\circ\ff{m-1}\circ\cdots\circ\ff{3})\big((\muu{2}_T\circ\muu{1}_T)({\mathbf{x}})\big) + \sqrt{D_{3,m+1}}L^{m-2}B\sig{2,1}_T\big(\muu{1}_T({\mathbf{x}})\big)\nonumber\\
    &+ \sqrt{D_{2,m+1}}L^{m-1} B\sig{1,1}_T({\mathbf{x}})\\

    \le{}& (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})+B\sig{m,1}\big((\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})\big)\nonumber\\
    &+LB\sig{m-1}\big((\muu{m-2}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})\big)+\cdots+ \sqrt{D_{2,m+1}}L^{m-1} B\sig{1,1}_T({\mathbf{x}})\\
    \le{}& (\muu{m}_T\circ\muu{m-1}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}) + B\sum_{i=1}^m \sqrt{D_{i+1,m+1}} L^{m-i} {\widetilde{\sigma}}_T^{(i,1)}({\mathbf{x}}).\label{eq:fd_2}
\end{align}
Hence, we have
\begin{align}
    r^\ast_T&=g({\mathbf{x}}^\ast)-g({\mathbf{x}}^\ast_T)\le 2 B\sum_{i=1}^m \sqrt{D_{i+1,m}} L^{m-i} \max\{{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast),{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast_T)\}.
\end{align}
It remains to upper bound $\max\{{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast),{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast_T)\}$ for $i\in[m]$. Firstly, by \cref{lem:fill}, we have
\begin{align}
	{\widetilde{\sigma}}^{(1,1)}_T({\mathbf{x}})=\sig{1,1}_T({\mathbf{x}})=O\big((\delta_T)^\nu\big)=O(T^{-\nu/d}).\label{eq:sigtil_net_1}
\end{align}
We now represent the upper bound on ${\widetilde{\sigma}}^{(i+1,1)}_T(\cdot)$ in terms of $\{{\widetilde{\sigma}}^{(j,1)}_T(\cdot):j\in[i]\}$ for $i\in[m-1]$. For all ${\mathbf{x}}\in\mathcal{X}$, the distance between $\muu{i}_T\circ\cdots\circ\muu{1}_T({\mathbf{x}})$ and its closest point in $\mathcal{X}^{(i+1)}$ is
\begin{align}
	\delta_{i+1}({\mathbf{x}})&=\min_{{\mathbf{z}}\in\mathcal{X}^{(i+1)}}|(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})-{\mathbf{z}}|\\
	&\le|(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})-\ff{i}\circ\cdots\circ\ff{1}({\mathbf{x}})|\\
	&\le B\sum_{j=1}^{i}\sqrt{D_{j+1,i+1}}L^{i-j}{\widetilde{\sigma}}^{(j,1)}_T({\mathbf{x}}),\label{eq:delta_i_net}
\end{align}
where the last step follows from \eqref{eq:fd_2}.
Then, the distance between $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})$ and its closest point in $\mathcal{X}^{(i+1)}_T$ is
\begin{align}
	\delta'_{i+1}({\mathbf{x}})&=\min_{s\in[T]}|(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}})-\xx{i+1}_s|\\
	&\le\delta^{(i+1)}_T+\delta_{i+1}({\mathbf{x}})\label{eq:delta_prime_net_1}\\
	&\le O(\sqrt{D_{2,i+1}}L^i T^{-1/d})+B\sum_{j=1}^{i}\sqrt{D_{j+1,i+1}}L^{i-j}{\widetilde{\sigma}}^{(j,1)}_T({\mathbf{x}})\label{eq:delta_prime_net_2}\\
	&= O\Big(\max\Big\{\sqrt{D_{2,i+1}}L^iT^{-1/d}, iB\cdot \max_{j\in[i]}\sqrt{D_{j+1,i+1}}L^{i-j} {\widetilde{\sigma}}^{(j,1)}_T({\mathbf{x}}) \Big\}\Big),
\end{align}
where \eqref{eq:delta_prime_net_1} applies the triangle inequality, and \eqref{eq:delta_prime_net_2} uses \eqref{eq:fill_i_net} and \eqref{eq:delta_i_net}.

Now, for each ${\mathbf{x}}'\in\{{\mathbf{x}}^\ast,{\mathbf{x}}^\ast_T\}$, we extend $\mathcal{X}^{(i+1)}$ to $\widetilde{\mathcal{X}}^{(i+1)}$, the smallest hyperrectangle that covers $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}')$ and $\mathcal{X}^{(i+1)}$. Then, the fill distance of the extended domain with regard to $\mathcal{X}^{(i+1)}_T$ is
\begin{align}
	\delta'_{i+1}=O\Big(\max\Big\{\sqrt{D_{2,i+1}}L^iT^{-1/d}, iB\cdot\max_{j\in[i]}\sqrt{D_{j+1,i+1}}L^{i-j} {\widetilde{\sigma}}^{(j,1)}_T({\mathbf{x}}') \Big\}\Big),
\end{align}
and by \cref{lem:fill}, we have
\begin{align}
	{\widetilde{\sigma}}^{(i+1,1)}_T({\mathbf{x}}')=O\big((\delta'_{i+1})^\nu\big)=O\Bigg(\Big(\max\Big\{\sqrt{D_{2,i+1}}L^iT^{-1/d}, iB\cdot\max_{j\in[i]}\sqrt{D_{j+1,i+1}}L^{i-j} {\widetilde{\sigma}}^{(j,1)}_T ({\mathbf{x}}') \Big\}\Big)^\nu\Bigg).\label{eq:sigtil_net}
\end{align}
We now consider two cases;

\underline{The case that $\nu \le 1$}:
Recalling the assumption $B=\Theta(L)$, and using the fact that ${\widetilde{O}}(L^{c\nu})$ is equivalent to $O\big( (L^{c\nu} {\rm poly}(c\log L)) \big)$ for fixed $\nu$, it follows from \eqref{eq:sigtil_net_1} and \eqref{eq:sigtil_net} that when $\nu\le 1$, we have
\begin{align}
	{\widetilde{\sigma}}_T^{(2,1)}({\mathbf{x}}')&=O\Big(\max\Big\{(d_2)^{\nu/2}L^\nu T^{-\nu/d}, (d_2)^{\nu/2} B^\nu T^{-\nu^2/d}\Big\}\Big)=O\big((d_2)^{\nu/2}L^\nu T^{-\nu^2/d}\big), \\
	{\widetilde{\sigma}}^{(3,1)}_T({\mathbf{x}}')&=O\Big(\max\Big\{(D_{2,3})^{\nu/2}L^{2\nu}T^{-\nu/d}, 2^\nu (D_{2,3})^{\nu/2} B^\nu  L^\nu T^{-\nu^2/d}, 2^\nu (d_3)^{\nu/2}  B^\nu (d_2)^{\nu^2/2} L^{\nu^2}  T^{-\nu^3/d}\Big\}\Big)\nonumber\\
	&=O\big(2^\nu(D_{2,3})^{\nu/2} B^\nu L^\nu T^{-\nu^3/d}\big)={\widetilde{O}}\big( (D_{2,3})^{\nu/2} B^\nu L^\nu T^{-\nu^3/d}\big),\\
	&\vdotswithin{=}\nonumber\\
	{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}')&={\widetilde{O}}\big((i-1)^\nu(D_{2,i})^{\nu/2} B^\nu L^{(i-2)\nu}T^{-\nu^i/d}\big)={\widetilde{O}}\big((D_{2,i})^{\nu/2} B^\nu L^{(i-2)\nu}T^{-\nu^i/d}\big),
\end{align}
and the simple regret is
\begin{align}
	r^\ast_T&\le 2B\sum_{i=1}^m \sqrt{D_{i+1,m}} L^{m-i}\max\{{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast),{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast_T)\}\\
	&={\widetilde{O}}(\sqrt{D_{2,m}}BL^{m-1}T^{-\nu^m/d}).
\end{align}

\underline{The case that $\nu > 1$}:
With $B=\Theta(L)$, and using the property of ${\widetilde{O}}(\cdot)$ similar to the first case, it follows from \eqref{eq:sigtil_net_1} and \eqref{eq:sigtil_net} that when $\nu > 1$, we have
\begin{align}
	{\widetilde{\sigma}}_T^{(2,1)}({\mathbf{x}}')&=O\Big(\max\Big\{(d_2)^{\nu/2}L^\nu T^{-\nu/d}, (d_2)^{\nu/2} B^\nu T^{-\nu^2/d}\Big\}\Big)=O((d_2)^{\nu/2}L^\nu T^{-\nu/d}), \\
	{\widetilde{\sigma}}^{(3,1)}_T({\mathbf{x}}')&=O\Big(\max\Big\{(D_{2,3})^{\nu/2}L^{2\nu}T^{-\nu/d}, (d_3)^{\nu/2}  B^\nu (d_2)^{\nu^2/2} L^{\nu^2} T^{-\nu^2/d}\Big\}\Big)\nonumber\\
	&=O\big(2^\nu(d_3)^{\nu/2} (d_2)^{\nu^2/2} B^\nu L^{\nu^2} T^{-\nu/d}\big)={\widetilde{O}}\big((d_3)^{\nu/2} (d_2)^{\nu^2/2} B^\nu L^{\nu^2} T^{-\nu/d}\big),\\
	&\vdotswithin{=}\nonumber\\
	{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}')&={\widetilde{O}}\Big(\max\Big\{(D_{2,i})^{\nu/2}L^{(i-1)\nu}T^{-\nu/d}, (i-1)^\nu (d_i)^{\nu/2}(d_{i-1})^{\nu^2/2}\cdots(d_2)^{\nu^{i-1}/2}B^{\nu+\cdots+\nu^{i-2}}L^{\nu^{i-1}}T^{-\nu^2/d}\Big\}\Big)\nonumber\\
	&={\widetilde{O}}\Big(\max\Big\{(D_{2,i})^{\nu/2}L^{(i-1)\nu}T^{-\nu/d}, (d_i)^{\nu/2}(d_{i-1})^{\nu^2/2}\cdots(d_2)^{\nu^{i-1}/2}B^{\nu+\cdots+\nu^{i-2}}L^{\nu^{i-1}}T^{-\nu^2/d}\Big\}\Big),
\end{align}
and the simple regret is
\begin{align}
	r^\ast_T&\le 2B\sum_{i=1}^m \sqrt{D_{i+1,m}} L^{m-i}\max\{{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast),{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast_T)\}\\
	&={\widetilde{O}}\Big(\max\Big\{(D_{2,m})^{\nu/2}BL^{(m-1)\nu}T^{-\nu/d}, {\widetilde{D}}_{2,m}^\nu B^{1+\nu+\nu^2+\cdots+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu^2/d} \Big\}\Big),
\end{align}
where ${\widetilde{D}}_{2,m}^\nu=\prod_{i=2}^m (d_i)^{\nu^{m+1-i}/2}$.
\end{proof}

Next, we recall the following two restrictive cases introduced in the main body:
\begin{itemize}
	\item \textbf{Case 1}: We additionally assume that $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}^\ast)\in\mathcal{X}^{(i+1)}$ and $(\muu{i}_T\circ\cdots\circ\muu{1}_T)({\mathbf{x}}^\ast_T)\in \mathcal{X}^{(i+1)}$ for $i\in[m-1]$.
	\item \textbf{Case 2}: We additionally assume that all the $\mathcal{X}^{(i)}$'s are known. Defining
	\begin{align}
		\muutil{i}_T({\mathbf{z}}) &= \argmin_{{\mathbf{z}}'\in\mathcal{X}^{(i+1)}}|\muu{i}_T({\mathbf{z}})-{\mathbf{z}}'|,
	\end{align}
	we let the algorithm return
	\begin{align}
		{\mathbf{x}}^\ast_T=\argmax_{{\mathbf{x}}\in\mathcal{X}}(\muutil{m}_T\circ\cdots\circ\muutil{1}_T)({\mathbf{x}}).
	\end{align}
\end{itemize}
Similarly to chains and multi-output chains, in either case, it follows that
\begin{align}
	r^\ast_T\le 2 B\sum_{i=1}^m \sqrt{D_{i+1,m}} L^{m-i} \max\{{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast),{\widetilde{\sigma}}^{(i,1)}_T({\mathbf{x}}^\ast_T)\} =O\Big( B\sum_{i=1}^m \sqrt{D_{i+1,m}} L^{m-i} {\overbar{\sigma}}_T^{(i,1)} \Big).
\end{align}
Moreover, by \cref{cor:fill_net}, we can further upper bound this by
\begin{align}
	r^\ast_T\le O\Big( B\sum_{i=1}^m \sqrt{D_{i+1,m}} L^{m-i} {\overbar{\sigma}}_T^{(i,1)} \Big)=
	\begin{cases}
		O(\sqrt{D_{2,m}}BL^{m-1}T^{-\nu/d})&\text{ when $\nu\le 1$,}\\
		O\big((D_{2,m})^{\nu/2}BL^{(m-1)\nu}T^{-\nu/d}\big)&\text{ when $\nu> 1$.}
	\end{cases}
\end{align}

\section{Lower Bound on Simple Regret (Proof of \cref{thm:lower_sim})}
\label{sec:lower_sim}

Recall the high-level intuition behind our lower bound outlined in Section \ref{sec:lower}: We use the idea from \cite{bull2011convergence} of having a small bump that is difficult to locate, but unlike \cite{bull2011convergence}, we exploit Lipschitz functions at the intermediate layers to ``amplify'' that bump (which means the original bump can have a much narrower width for a given RKHS norm).

We first show that, for fixed $\epsilon>0$, we can construct a base function ${\overbar{g}}$ with height $2\epsilon$ and support radius $w=\Theta\big((\frac{\epsilon}{ BL^{m-1}})^{1/\nu}\big)$ for each network structure.

\subsection{Hard Function for Chains}

To construct a base function with chain structure, we define the following two scalar-valued functions.

Considering the ``bump'' function $h({\mathbf{x}})=\exp\big(\frac{-1}{1-\|{\mathbf{x}}\|_2^2}\big)\mathds{1}\{\|{\mathbf{x}}\|_2< 1\}$, for some $\epsilon_1>0$ and width $w>0$, we define \cite{bull2011convergence}
\begin{align}
    {\widetilde{h}}({\mathbf{x}},\epsilon_1,w)&=\frac{2\epsilon_1}{h(\mathbf{0})}h\Big(\frac{{\mathbf{x}}}{w}\Big),\label{eq:htil}
\end{align}
which is a scaled bump function with height $2\epsilon_1$ and compact support $\{{\mathbf{x}}\in\mathcal{X}:\|{\mathbf{x}}\|_2<w\}$.

For a fixed $u>0$ and ${\widetilde{L}}=\sqrt{2}B/\sqrt{k(0)-k(2u)}$, where $k$ is the Mat\'ern kernel on ${\mathbb{R}}$ with smoothness $\nu$ and $k(|x-x'|)=k(x,x')$, we define
\begin{align}
    {\widetilde{g}}_k(\cdot)& = \frac{{\widetilde{L}}}{2} \big(k(\cdot,u)-k(\cdot,-u)\big).\label{eq:gtilk}
\end{align}

\begin{restatable}[Hard function for chains]{thm}{hardchain}
\label{thm:hard_chain}
For the Mat\'ern kernel $k$ with smoothness $\nu\ge1$, sufficiently small $\epsilon>0$, and sufficiently large $B$, there exists $\epsilon_1>0,L=\Theta(B), c=\Theta(1), w=\Theta\big((\frac{\epsilon}{ B(cL)^{m-1}})^{1/\nu}\big),$ and $u>0$ such that ${\overbar{g}}:={\overbar{f}}^{(m)}\circ{\overbar{f}}^{(m-1)}\circ\cdots\circ{\overbar{f}}^{(1)}$ with ${\overbar{f}}^{(1)}(\cdot)={\widetilde{h}}(\cdot,\epsilon_1,w)$ and ${\overbar{f}}^{(s)}(\cdot)={\widetilde{g}}_k(\cdot)$ for $s\in [2,m]$ has the following properties:
\begin{itemize}
    \item ${\overbar{f}}^{(i)}\in\mathcal{H}_k(B)\cap\mathcal{F}(L)$ for each $i\in[m]$,
    \item $\max_{{\mathbf{x}}} {\overbar{g}}({\mathbf{x}})= 2\epsilon$,
    \item ${\overbar{g}}({\mathbf{x}})>0$ when $\|{\mathbf{x}}\|_2<w$ , and ${\overbar{g}}({\mathbf{x}})=0$ otherwise.
\end{itemize}
\end{restatable}
\begin{proof}
For $k$ being the Mat\'ern kernel with smoothness $\nu$, recalling that
\begin{align}
    {\widetilde{h}}({\mathbf{x}},\epsilon_1,w)&=\frac{2\epsilon_1}{h(\mathbf{0})}h\Big(\frac{{\mathbf{x}}}{w}\Big)
\end{align}
with $h$ being the bump function, \citep[Section A.2]{bull2011convergence} has shown that for some constant $C_1$,
\begin{align}
    \|{\widetilde{h}}\|_k\le C_1 \frac{2\epsilon_1}{h(\mathbf{0})} \Big(\frac{1}{w}\Big)^\nu \|h\|_k.
\end{align}
Hence, we have $\|{\widetilde{h}}\|_k\le B$ when $w=\big(\frac{2C_1\|h\|_k\epsilon_1}{h(\mathbf{0})B}\big)^{1/\nu}$. When $\frac{\epsilon_1}{B}$ is sufficiently small, the diameter of the support satisfies $2w\ll 1$. Since the bump function is infinitely differentiable, ${\widetilde{h}}$ is Lipschitz continuous with some constant $L'=\Theta(\frac{\epsilon_1}{w})$, and our assumption of $\nu\ge 1$ further implies $L'=O(B)$.

Recall that for a fixed $u>0$ and the Mat\'ern kernel $k$ on a one-dimension domain, we define
\begin{align}
    {\widetilde{g}}_k(\cdot)& = \frac{{\widetilde{L}}}{2} \big(k(\cdot,u)-k(\cdot,-u)\big)
\end{align}
with ${\widetilde{L}}=\frac{\sqrt{2}B}{\sqrt{k(0)-k(2u)}}$. Applying \eqref{eq:norm} with ${\widetilde{g}}_k$ replacing both $f$ and $f'$, we have
\begin{align}
    \|{\widetilde{g}}_k\|_k=\frac{{\widetilde{L}}}{2}\sqrt{k(u,u)+2k(u,-u)+k(-u,-u)}=\frac{{\widetilde{L}}}{2}\sqrt{2k(0)-2k(2u)}=B,
\end{align}
where $k(d_{{\mathbf{x}},{\mathbf{x}}'})=k({\mathbf{x}},{\mathbf{x}}')$ with $d_{{\mathbf{x}},{\mathbf{x}}'}=\|{\mathbf{x}}-{\mathbf{x}}'\|_2$.

We choose the two constants $u>{\widetilde{u}}>0$ to satisfy the following:
\begin{enumerate}
    \item $k(d_{{\mathbf{x}},{\mathbf{x}}'})$ is non-increasing when $d_{{\mathbf{x}},{\mathbf{x}}'}\in[u-{\widetilde{u}},u+{\widetilde{u}}]$.
    \item $k(u-{\widetilde{u}})-k(u+{\widetilde{u}})\ge\frac{2}{{\widetilde{L}}}=\frac{\sqrt{2}\sqrt{k(0)-k(2u)}}{B}$.
    \item Defining $r_{\max}=\sup_{z\in(0,{\widetilde{u}}]}\frac{k(u-z)-k(u+z)}{2z}$ and $r_{\min}=\inf_{z\in(0,{\widetilde{u}}]}\frac{k(u-z)-k(u+z)}{2z}$, it holds that $r_{\max}=\Theta(1)$ and $r_{\min}=\Theta(1)$.
\end{enumerate}
Whenever $k(d_{{\mathbf{x}},{\mathbf{x}}'})$ is continuous and non-increasing in $d_{{\mathbf{x}},{\mathbf{x}}'}$ on ${\mathbb{R}}^+$ (e.g., the Mat\'ern kernel), condition 1 is satisfied. This condition guarantees that $k(\cdot, u)$ is non-decreasing on $[0,{\widetilde{u}}]$, and $k(\cdot, -u)$ is non-increasing on $[0,{\widetilde{u}}]$, which together implies that ${\widetilde{g}}_k$ is non-decreasing on $[0,{\widetilde{u}}]$. Under condition 1, condition 2 is satisfied as long as $B$ is sufficiently large, and it implies that ${\widetilde{g}}_k({\widetilde{u}})\ge 1$. Examples of ${\widetilde{g}}_k$ for $k=k_\text{Mat\'ern}$ are given in Figure \ref{fig:gtil}.

Let $L=\max\{L',\sup_{z\in(0,{\widetilde{u}}]}\frac{{\widetilde{g}}_k(z)}{z}\}=\max\{L',r_{\max}{\widetilde{L}}\}$ and $\alpha=\inf_{z\in(0,{\widetilde{u}}]}\frac{{\widetilde{g}}_k(z)}{zL}=\frac{r_{\min}{\widetilde{L}}}{L}$. Then, it holds for all $z\in(0,{\widetilde{u}}]$ that
\begin{align}
    1<\alpha L\le \frac{{\widetilde{g}}_k(z)}{z} \le L. \label{eq:slope}
\end{align}
Clearly $\alpha$ cannot exceed $1$, and moreover, condition 3 above along with $L'=O(B)$ implies that $\alpha=\min\big\{r_{\min}\sqrt{\frac{2}{k(0)-k(2u)}}\frac{B}{L'},\frac{r_{\min}}{r_{\max}}\big\}=\Theta(1)$.

\begin{figure*}[t!]
	\minipage[t]{0.33\textwidth}
	\includegraphics[width=\linewidth]{figs/g_100.pdf}
	\caption*{$\nu=1$}
	\endminipage
	\minipage[t]{0.33\textwidth}
	\includegraphics[width=\linewidth]{figs/g_150.pdf}
	\caption*{$\nu=3/2$}
	\endminipage
	\minipage[t]{0.33\textwidth}
	\includegraphics[width=\linewidth]{figs/g_200.pdf}
	\caption*{$\nu=2$}
	\endminipage
	\caption{$k=k_\text{Mat\'ern}$ with lengthscale $l=1$ and smoothness $\nu$, and ${\widetilde{g}}_k$ with $B=5$, $u=0.5$, and ${\widetilde{u}}=0.3$.}
	\label{fig:gtil}
\end{figure*}

Defining ${\overbar{g}}_i={\overbar{f}}^{(i)}\circ{\overbar{f}}^{(i-1)}\circ\dots\circ{\overbar{f}}^{(1)}$, we denote the height of ${\overbar{g}}_i$ by $h_i=\max_{\mathbf{x}} {\overbar{g}}_i({\mathbf{x}})$. If $h_{i-1}\le{\widetilde{u}}$, \eqref{eq:slope} implies
\begin{align}
    h_{i-1}<\alpha L h_{i-1}\le h_i \le L h_{i-1}.\label{eq:height}
\end{align}
Then, for a given $\epsilon\in(0, {\widetilde{g}}_k({\widetilde{u}})]$, where condition 2 implies ${\widetilde{g}}_k({\widetilde{u}})\ge 1$, we choose $\epsilon_1$ to satisfy
\begin{align}
    ({\overbar{f}}^{(m)} \circ \cdots \circ {\overbar{f}}^{(2)}) (2\epsilon_1)=({\widetilde{g}}_k \circ \cdots \circ {\widetilde{g}}_k) (2\epsilon_1)=2\epsilon.
\end{align}
By $h_1=2\epsilon_1$ and \eqref{eq:height}, this choice of $\epsilon_1$ must also satisfy
\begin{align}
    2 (\alpha L)^{m-1} \epsilon_1 = (\alpha L)^{m-1} h_1 \le h_m = 2\epsilon \le L^{m-1} h_1 = 2 L^{m-1} \epsilon_1,
\end{align}
implying (via \eqref{eq:slope}) that
\begin{align}
    \frac{\epsilon}{L^{m-1}} \le \epsilon_1 \le \frac{\epsilon}{(\alpha L)^{m-1}}<\epsilon.
\end{align}
Since ${\widetilde{h}}$ has a compact support with radius $w$ and ${\widetilde{g}}(0)=0$, ${\overbar{g}}$ also has a compact support with radius
\begin{align}
    w = \Theta\bigg(\Big(\frac{\epsilon_1}{B}\Big)^{1/\nu}\bigg) = \Theta \bigg(\Big(\frac{\epsilon}{B(cL)^{m-1}}\Big)^{1/\nu}\bigg),
\end{align}
for some constant $c=\Theta(1)$.

Lastly, \eqref{eq:slope} implies ${\widetilde{g}}_k$ on $[0,{\widetilde{u}}]$ is a member of $\mathcal{F}(L)$, and $L\ge L'$ guarantees ${\widetilde{h}}\in\mathcal{F}(L)$.
\end{proof}

\subsection{Hard Function for Multi-Output Chains}

In this section, we consider the operator-valued Mat\'ern kernel $\Gamma^{(i)}({\mathbf{x}},{\mathbf{x}}')=k^{(i)}({\mathbf{x}},{\mathbf{x}}'){\mathbf{I}}^{(i+1)}$, where $k^{(i)}$ is the scalar-valued Mat\'ern kernel on ${\mathbb{R}}^{d_i}$ and ${\mathbf{I}}^{(i+1)}$ is the identity matrix of size $d_{i+1}$. Then, for a fixed ${\mathbf{u}}^{(i)}\in{\mathbb{R}}^{d_i}$ and ${\widetilde{L}}^{(i)}=\sqrt{2}B/\sqrt{k^{(i)}(0)-k^{(i)}(2\|{\mathbf{u}}^{(i)}\|_2)}$ where $k^{(i)}(\|{\mathbf{x}}-{\mathbf{x}}'\|_2)=k^{(i)}({\mathbf{x}},{\mathbf{x}}')$, we define
\begin{align}
    {\widetilde{g}}^{(i)}(\cdot)& = \frac{{\widetilde{L}}^{(i)}}{2} \big(\Gamma^{(i)}(\cdot,{\mathbf{u}}^{(i)})-\Gamma^{(i)}(\cdot,-{\mathbf{u}}^{(i)})\big).
\end{align}

\begin{restatable}[Hard function for multi-output chains]{thm}{hardmul}
\label{thm:hard_mul}
Let ${\mathbf{I}}^{(i)}$ denote the identity matrix of size $d_i$, let ${\mathbf{e}}^{(i)}_1$ denote the first column of ${\mathbf{I}}^{(i)}$, and let $k^{(i)}$ denote the scalar-valued Mat\'ern kernel on ${\mathbb{R}}^{d_i}$ with smoothness $\nu\ge1$. For $\Gamma^{(i)}({\mathbf{x}},{\mathbf{x}}')=k^{(i)}({\mathbf{x}},{\mathbf{x}}'){\mathbf{I}}^{(i+1)}$, sufficiently small $\epsilon>0$, and sufficiently large $B$, there exists $\epsilon_1>0,L=\Theta(B),c=\Theta(1),w=\Theta\big((\frac{\epsilon}{ B(cL)^{m-1}})^{1/\nu}\big),$ and $u\in{\mathbb{R}}$ such that ${\overbar{g}}:={\overbar{f}}^{(m)}\circ{\overbar{f}}^{(m-1)}\circ\cdots\circ{\overbar{f}}^{(1)}$ with ${\overbar{f}}^{(1)}(\cdot)={\widetilde{h}}(\cdot,\epsilon_1,w){\mathbf{e}}_1^{(2)}$ and ${\overbar{f}}^{(s)}(\cdot)={\widetilde{g}}^{(s)}(\cdot){\mathbf{e}}_1^{(s+1)}$ with ${\mathbf{u}}^{(s)}=u{\mathbf{e}}_1^{(s)}$ for $s\in [2,m]$ has the following properties:
\begin{itemize}
    \item ${\overbar{f}}^{(i)}\in\mathcal{H}_{\Gamma^{(i)}}(B)\cap\mathcal{F}(L)$ for each $i\in[m]$,
    \item $\max_{{\mathbf{x}}} {\overbar{g}}({\mathbf{x}})= 2\epsilon$,
    \item ${\overbar{g}}({\mathbf{x}})>0$ when $\|{\mathbf{x}}\|_2<w$ , and ${\overbar{g}}({\mathbf{x}})=0$ otherwise.
\end{itemize}
\end{restatable}

\begin{proof}
The rough idea is to reduce to the case of regular chains by only making use of a single coordinate throughout the network.  We leave open the question as to whether the lower bound can be improved by utilizing all coordinates.

With ${\mathbf{I}}$ denoting the identity matrix of size $d_2$ and ${\mathbf{e}}_1$ denoting the first column of ${\mathbf{I}}$, we aim to show that ${\overbar{f}}^{(1)}\in\mathcal{H}_{\Gamma^{(1)}}(B)$ by showing that if $\Gamma(\cdot,\cdot)=k(\cdot,\cdot){\mathbf{I}}\in{\mathbb{R}}^{d_2\times d_2}$ and ${\widetilde{h}}\in\mathcal{H}_k(B)$ for some scalar-valued kernel $k$ and some constant $B$, then the function ${\widetilde{h}}'(\cdot)={\widetilde{h}}(\cdot){\mathbf{e}}_1$ satisfies ${\widetilde{h}}' \in\mathcal{H}_\Gamma(B)$. Since ${\widetilde{h}}\in\mathcal{H}_k(B)$, there exists a sequence $\{(a_i,{\mathbf{x}}_i)\}_{i=1}^\infty$ such that ${\widetilde{h}}(\cdot)=\sum_{i=1}^\infty a_i k(\cdot,{\mathbf{x}}_i)$. Then, using the definition of RKHS norm for vector-valued functions in \eqref{eq:mul_norm}, we have
\begin{align}
    \|{\widetilde{h}}'\|_\Gamma^2&=\lim_{n\to\infty}\Big\|\sum_{i=1}^n a_i k(\cdot,{\mathbf{x}}_i){\mathbf{e}}_1\Big\|_\Gamma^2\\
    &=\lim_{n\to\infty}\Big\|\sum_{i=1}^n a_i \Gamma(\cdot,{\mathbf{x}}_i){\mathbf{e}}_1\Big\|_\Gamma^2\\
    &=\sum_{i,j=1}^\infty \langle \Gamma({\mathbf{x}}_i, {\mathbf{x}}_j)(a_i{\mathbf{e}}_1),a_j{\mathbf{e}}_1\rangle \\
    &=\sum_{i,j=1}^\infty a_i a_j k({\mathbf{x}}_i,{\mathbf{x}}_j)\\
    &=\|{\widetilde{h}}\|_k^2 \\
    &\le B^2,
\end{align}
and therefore ${\widetilde{h}}'\in\mathcal{H}_\Gamma(B)$. Next, for $i\in[2,m]$, with ${\mathbf{u}}^{(i)}=u{\mathbf{e}}_1^{(i)}$ we have
\begin{align}
    \|{\overbar{f}}^{(i)}\|_\Gamma &= \frac{{\widetilde{L}}^{(i)}}{2}\sqrt{ \langle \Gamma^{(i)}({\mathbf{u}}^{(i)},{\mathbf{u}}^{(i)}){\mathbf{e}}_1,{\mathbf{e}}_1\rangle - 2\langle \Gamma^{(i)}({\mathbf{u}}^{(i)},-{\mathbf{u}}^{(i)}){\mathbf{e}}_1,{\mathbf{e}}_1\rangle + \langle \Gamma^{(i)}(-{\mathbf{u}}^{(i)},-{\mathbf{u}}^{(i)}){\mathbf{e}}_1,{\mathbf{e}}_1\rangle } \\
    &= \frac{{\widetilde{L}}^{(i)}}{2}\sqrt{k^{(i)}({\mathbf{u}}^{(i)},{\mathbf{u}}^{(i)})-2k^{(i)}({\mathbf{u}}^{(i)},-{\mathbf{u}}^{(i)})+k^{(i)}(-{\mathbf{u}}^{(i)},-{\mathbf{u}}^{(i)})}\\
    &= \frac{{\widetilde{L}}^{(i)}}{2}\sqrt{2k^{(i)}(0)-2k^{(i)}(2u)}\\
    &=B.
\end{align}
Since ${\widetilde{h}}\in\mathcal{F}(L')$ for some $L'>1$, we also have ${\widetilde{h}}'\in\mathcal{F}(L')$.

We reuse the choice of $u$ and ${\widetilde{u}}$ in the previous case. Then, for any $z\in{\mathbb{R}}$, with ${\mathbf{z}}^{(i)}=z{\mathbf{e}}_1^{(i)}$ and ${\mathbf{u}}^{(i)}=u{\mathbf{e}}_1^{(i)}$, we have
\begin{align}
    {\overbar{f}}^{(i)}({\mathbf{z}})&={\widetilde{g}}^{(i)}({\mathbf{z}}^{(i)},{\mathbf{u}}^{(i)}){\mathbf{e}}_1^{(i+1)}\\
    &=\frac{{\widetilde{L}}^{(i)}}{2}\big(\Gamma^{(i)}({\mathbf{z}}^{(i)},{\mathbf{u}}^{(i)})-\Gamma^{(i)}({\mathbf{z}}^{(i)},-{\mathbf{u}}^{(i)})\big){\mathbf{e}}_1^{(i+1)}\\
    &=\frac{{\widetilde{L}}^{(i)}}{2}\big(k^{(i)}({\mathbf{z}}^{(i)},{\mathbf{u}}^{(i)})-k^{(i)}({\mathbf{z}}^{(i)},-{\mathbf{u}}^{(i)})\big){\mathbf{e}}_1^{(i+1)}\\
    &=\frac{{\widetilde{L}}^{(i)}}{2}\big(k(\|{\mathbf{z}}^{(i)}-{\mathbf{u}}^{(i)}\|)-k(\|{\mathbf{z}}^{(i)}+{\mathbf{u}}^{(i)}\|)\big){\mathbf{e}}_1^{(i+1)}\\
    &=\frac{{\widetilde{L}}}{2}\big(k(|z-u|)-k(|z+u|)\big){\mathbf{e}}_1^{(i+1)}\\
    &={\widetilde{g}}_k(z){\mathbf{e}}_1^{(i+1)},
\end{align}
where ${\widetilde{g}}_k(\cdot)$ depending on $u$ is defined in \eqref{eq:gtilk}. Hence, as illustrated in \cref{fig:mul_lower}, for any input ${\mathbf{x}}\in[0,1]^d$ of ${\overbar{g}}$, we have
\begin{align}
    \xx{2,1}&={\widetilde{h}}({\mathbf{x}}),\\
    \xx{2,j}&=0 &&\hspace{-25ex}\text{ for }j\ge 2, \\
    \xx{i+1,1}&={\widetilde{g}}_k(\xx{i,1})&&\hspace{-25ex}\text{ for }i\ge 2,\\
    \xx{i+1,j}&=0&&\hspace{-25ex}\text{ for }i\ge 2\text{ and }j\ge 2.
\end{align}
\begin{figure}[t!]
    \centering
    \includegraphics[width=0.6\textwidth]{figs/mul_lower.pdf}
    \caption{Illustration of ${\overbar{g}}$ for multi-output chain.}
    \label{fig:mul_lower}
\end{figure}
By a similar argument to the case of single-output chains, there exists $L\ge L'$ and $\alpha=\Theta(1)$ such that for all $z\in[0,{\widetilde{u}}]$,
\begin{align}
    1<\alpha L\le \frac{{\widetilde{g}}_k(z)}{z}\le L.\label{eq:mul_slope}
\end{align}

For $\epsilon\le{\widetilde{g}}_k({\widetilde{u}})$, we choose $\epsilon_1$ to satisfy
\begin{align}
    ({\overbar{f}}^{(m)}\circ{\overbar{f}}^{(m-1)}\circ\cdots\circ{\overbar{f}}^{(2)})(2\epsilon_1{\mathbf{e}}_1^{(2)})=({\widetilde{g}}_k\circ\cdots\circ{\widetilde{g}}_k)(2\epsilon_1)=2\epsilon,
\end{align}
and ${\overbar{g}}$ has a compact support with radius
\begin{align}
    w = \Theta\Bigg(\Big(\frac{\epsilon_1}{B}\Big)^{1/\nu}\Bigg) = \Theta \Bigg(\Big(\frac{\epsilon}{B(cL)^{m-1}}\Big)^{1/\nu}\Bigg),
\end{align}
for some constant $c=\Theta(1)$.

Lastly, since $L\ge L'$, we immediately deduce that ${\overbar{f}}^{(i)}\in\mathcal{F}(L)$ on its domain for each $i\in[m]$.
\end{proof}

\subsection{Hard Function for Feed-Forward Networks}

For a fixed ${\mathbf{u}}^{(i)}\in{\mathbb{R}}^{d_i}$ and ${\widetilde{L}}^{(i)}=\sqrt{2}B/\sqrt{k^{(i)}(0)-k^{(i)}(2\|{\mathbf{u}}^{(i)}\|_2)}$, where $k^{(i)}(\|{\mathbf{x}}-{\mathbf{x}}'\|_2)=k^{(i)}({\mathbf{x}},{\mathbf{x}}')$, we define
\begin{align}
    {\widetilde{g}}_k^{(i)}(\cdot)& = \frac{{\widetilde{L}}^{(i)}}{2} \big(k^{(i)}(\cdot,{\mathbf{u}}^{(i)})-k^{(i)}(\cdot,-{\mathbf{u}}^{(i)})\big).
\end{align}

\begin{restatable}[Hard function for feed-forward networks]{thm}{hardnet}
\label{thm:hard_net}
Let ${\mathbf{e}}_1^{(i)}$ denote the first column of identity matrix of size $d_i$. For the Mat\'ern kernel $k$ with smoothness $\nu\ge 1$, sufficiently small $\epsilon>0$, and sufficiently large $B$, there exists $\epsilon_1>0,L=\Theta(B),c=\Theta(1), w=\Theta\big((\frac{\epsilon}{ B(cL)^{m-1}})^{1/\nu}\big),$ and $u>0$ such that ${\overbar{g}}:={\overbar{f}}^{(m)}\circ{\overbar{f}}^{(m-1)}\circ\cdots\circ{\overbar{f}}^{(1)}$ with ${\overbar{f}}^{(i)}(\cdot)=[{\overbar{f}}^{(i,j)}(\cdot)]_{j=1}^{d_{i+1}}$ for $i\in[m]$, ${\overbar{f}}^{(1,1)}(\cdot)={\widetilde{h}}(\cdot,\epsilon_1,w)$, ${\overbar{f}}^{(s,1)}(\cdot)={\widetilde{g}}_k^{(s)}(\cdot)$ with ${\mathbf{u}}^{(s)}=u{\mathbf{e}}_1^{(s)}$, and ${\overbar{f}}^{(s,r)}(\cdot)=0$ for $s\in [2,m]$ and $r\neq 1$ has the following properties:
\begin{itemize}
    \item ${\overbar{f}}^{(i,j)}\in\mathcal{H}_{k}(B)\cap\mathcal{F}(L)$ for each $i\in[m],j\in[d_{i+1}]$;
    \item $\max_{{\mathbf{x}}} {\overbar{g}}({\mathbf{x}})= 2\epsilon$;
    \item ${\overbar{g}}({\mathbf{x}})>0$ when $\|{\mathbf{x}}\|_2<w$ , and ${\overbar{g}}({\mathbf{x}})=0$ otherwise.
\end{itemize}
\end{restatable}

\begin{proof}
We adopt a similar general approach to the case of chains and multi-output chains, but with some different details.

As noted in our analysis of chains, we have ${\widetilde{h}}\in\mathcal{H}_{k}(B)$ and ${\widetilde{h}}\in\mathcal{F}(L')$ for some constant $L'>1$, and we also have
\begin{align}
    \|{\widetilde{g}}^{(i)}_k\|_k&=\frac{{\widetilde{L}}^{(i)}}{2}\sqrt{k^{(i)}({\mathbf{u}}^{(i)},{\mathbf{u}}^{(i)})-2k^{(i)}({\mathbf{u}}^{(i)},-{\mathbf{u}}^{(i)})+k^{(i)}(-{\mathbf{u}}^{(i)},-{\mathbf{u}}^{(i)})}\\
    &=\frac{{\widetilde{L}}^{(i)}}{2}\sqrt{2k^{(i)}(0)-2k^{(i)}(2\|{\mathbf{u}}^{(i)}\|_2)}\\
    &= B.
\end{align}
Reusing the choice of $u$ and ${\widetilde{u}}$ in the case of chains, with ${\mathbf{z}}^{(i)}=z{\mathbf{e}}_1^{(i)}$ and ${\mathbf{u}}^{(i)}=u{\mathbf{e}}_1^{(i)}$, we have
\begin{align}
    {\widetilde{g}}^{(i)}_k({\mathbf{z}}^{(i)})&= \frac{{\widetilde{L}}^{(i)}}{2} \big(k^{(i)}({\mathbf{z}}^{(i)},{\mathbf{u}}^{(i)})-k^{(i)}({\mathbf{z}}^{(i)},-{\mathbf{u}}^{(i)})\big)\\
    &= \frac{{\widetilde{L}}^{(i)}}{2} \big(k^{(i)}(\|{\mathbf{z}}^{(i)}-{\mathbf{u}}^{(i)}\|_2)-k^{(i)}(\|{\mathbf{z}}^{(i)}+{\mathbf{u}}^{(i)}\|_2)\big)\\
    &= \frac{{\widetilde{L}}^{(i)}}{2} \big(k^{(i)}(|z-u|)-k^{(i)}(|z+u|)\big)\\
    &= {\widetilde{g}}_k(z)
\end{align}
Hence, as illustrated in \cref{fig:net_lower}, for any input ${\mathbf{x}}\in[0,1]^d$ of ${\overbar{g}}$, we have
\begin{align}
    \xx{2,1}&={\widetilde{h}}({\mathbf{x}}),\\
    \xx{2,j}&=0 &&\hspace{-25ex}\text{ for }j\ge 2, \\
    \xx{i+1,1}&={\widetilde{g}}_k(\xx{i,1})&&\hspace{-25ex}\text{ for }i\ge 2,\\
    \xx{i+1,j}&=0&&\hspace{-25ex}\text{ for }i\ge 2\text{ and }j\ge 2.
\end{align}
\begin{figure}[t!]
    \centering
    \includegraphics[width=0.6\textwidth]{figs/net_lower.pdf}
    \caption{Illustration of ${\overbar{g}}$ for feed-forward network.}
    \label{fig:net_lower}
\end{figure}
By a similar argument to the previous cases, there exists $L\ge L'$ and $\alpha=\Theta(1)$ such that for all $z\in[0,{\widetilde{u}}]$,
\begin{align}
    1<\alpha L\le \frac{{\widetilde{g}}_k(z)}{ z}\le L.\label{eq:net_slope}
\end{align}

For $\epsilon\le{\widetilde{g}}_k({\widetilde{u}})$, we choose $\epsilon_1$ satisfying
\begin{align}
    ({\overbar{f}}^{(m)}\circ{\overbar{f}}^{(m-1)}\circ\cdots\circ{\overbar{f}}^{(2)})(2\epsilon_1{\mathbf{e}}_1^{(2)})=({\widetilde{g}}_k\circ\cdots\circ{\widetilde{g}}_k)(2\epsilon_1)=2\epsilon,
\end{align}
and ${\overbar{g}}$ has a compact support with radius
\begin{align}
    w = \Theta\Bigg(\Big(\frac{\epsilon_1}{B}\Big)^{1/\nu}\Bigg) = \Theta \Bigg(\Big(\frac{\epsilon}{B(cL)^{m-1}}\Big)^{1/\nu}\Bigg),
\end{align}
for some constant $c=\Theta(1)$.

Lastly, due to $L\ge L'$, ${\overbar{f}}^{(i,j)}\in\mathcal{F}(L)$ on its domain for each $i\in[m],j\in[d_{i+1}]$.
\end{proof}

\subsection{Lower Bound on Simple Regret}

With the preceding ``hard functions'' established, the final step is to essentially follow that of \cite{bull2011convergence}.  We provide the details for completeness.

By splitting the domain $\mathcal{X}=[0,1]^d$ into a grid of dimension $d$ with spacing $2w$, we construct $M=(\lfloor \frac{1}{2w} \rfloor)^d$ functions with disjoint supports by shifting the origin of ${\overbar{g}}$ to the center of each cell and cropping the shifted function into $[0,1]^d$, which is denoted by $\mathcal{G}=\{g_1,\dots,g_M\}$. For $g$ sampled uniformly from $\mathcal{G}$, we first show that the expected simple regret ${\mathbb{E}}[r_T^\ast]$ of an arbitrary algorithm is lower bounded, then there must exist a function in $\mathcal{G}$ that has the same lower bound.

Now, we prove \cref{thm:lower_sim}, which is restated as follows.
\lowersim*
\begin{proof}
    Since the theorem concerns worst-case RKHS functions, it suffices to establish the same result when the function $g$ is drawn uniformly from the hard subset $\mathcal{G}$ introduced above.  Given that $g$ is random, Yao's minimax principle implies that it suffices to consider deterministic algorithms.

    For any given deterministic algorithm, let $\mathcal{X}'_T=\{{\mathbf{x}}'_1,{\mathbf{x}}'_2,\cdots,{\mathbf{x}}'_T\}$ be the set of points that it would sample if the function were zero everywhere.  We observe that if $g = g_i$ for some $g_i$ satisfying $g({\mathbf{x}}'_1) = g({\mathbf{x}}'_2) = \dotsc = g({\mathbf{x}}'_T) = 0$ (i.e., the bump in $g_i$ does not cover any of the points in $\mathcal{X}'_T$), then the algorithm will precisely sample ${\mathbf{x}}'_1,{\mathbf{x}}'_2,\cdots,{\mathbf{x}}'_T$.  Moreover, if there are multiple such functions $g_i$, then the final returned point ${\mathbf{x}}_T^*$ can only be below $2\epsilon$ (i.e., the bump height) for at most one of those functions, since their supports are disjoint by construction.

    Now suppose that $T \le \frac{M}{2} - 1$, where $M$ is the number of functions in $\mathcal{G}$.  This means that regardless of the values ${\mathbf{x}}'_1,{\mathbf{x}}'_2,\cdots,{\mathbf{x}}'_T$, there are at least $\frac{M}{2}+1$ functions in $\mathcal{G}$ such that the sampled values are all zero.  Hence, there are at least $\frac{M}{2}$ functions where a simple regret of at least $2\epsilon$ is incurred, meaning that the average simple regret is at least $\epsilon$.

    With this result in place, \cref{thm:lower_sim} immediately follows by substituting $M=\Theta\big((1/w)^d\big)$ and $w=\Theta \big((\frac{\epsilon}{B(cL)^{m-1}})^{1/\nu}\big)$.
\end{proof}

\section{Lower Bound on Cumulative Regret (Proof of \cref{thm:lower_cml})}
\label{sec:lower_cml}
In this section, we prove \cref{thm:lower_cml}, which is restated as follows.
\lowercml*

\begin{proof}
By rearranging \cref{thm:lower_sim}, we have
$\epsilon=\Omega(B(cL)^{m-1}T^{-\nu/d})$, which implies that the lower bound on cumulative regret is
\begin{align}
	{\mathbb{E}}[R_T]=\Omega(\epsilon T)=\Omega(B(cL)^{m-1}T^{1-\nu/d}).
\end{align}
However, this lower bound is loose when $d<\nu$. \cref{thm:lower_sim} implies that to have simple regret at most $\epsilon=\Theta(1)$ requires $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$. When $T= \Theta\big((B(cL)^{m-1})^{d/\nu}\big)$, we have ${\mathbb{E}}[r^\ast_T]=\Omega(1)$ and ${\mathbb{E}}[R_T]=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$. Since cumulative regret is always non-decreasing in $T$, when $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$, we also have ${\mathbb{E}}[R_T]=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$.

For $T=o\big((B(cL)^{m-1})^{d/\nu}\big)$ (which is equivalent to $T=o(B(cL)^{m-1}T^{1-\nu/d})$), we show by contradiction that ${\mathbb{E}}[R_T]=\Omega(T)$. Suppose on the contrary that there exists an algorithm guaranteeing ${\mathbb{E}}[R_T]=o(T)$ when $T=T_0$ for some $T_0 = o\big((B(cL)^{m-1})^{d/\nu}\big)$. Then, by repeatedly selecting the best point among the first $T_0$ time steps, the algorithm attains ${\mathbb{E}}[R_T]=o\big((B(cL)^{m-1})^{d/\nu}\big)$ when $T=\Theta\big((B(cL)^{m-1})^{d/\nu}\big)$, which contradicts the lower bound for $T=\Omega\big((B(cL)^{m-1})^{d/\nu}\big)$.

Hence, \cref{thm:lower_cml} follows by combining the two cases.

\end{proof}

\section{Summary of Regret Bounds}
\label{sec:summary}
A detailed summary of our regret bounds for the Mat\'ern kernel is given in \cref{tab:summary}.
\begin{table}[t!]
	\centering
	\begin{tabular}{ll}
		\toprule
		\textbf{Algorithm-Independent Cumulative Regret Lower Bound} & \\
		Chains/Multi-Output Chains/Feed-Forward Networks & $\begin{cases}
			\Omega\big(\min\{T, B(cL)^{m-1}T^{1-\nu/d}\}\big)&\text{ when }d>\nu\ge1,\\
			\Omega\big(\min\{T, (B(cL)^{m-1})^{d/\nu}\}\big)&\text{ when } d\le\nu.
		\end{cases}$ \\
		\textbf{Algorithmic Cumulative Regret Upper Bound} \\ \textbf{(If Conjecture of \cite{vakili2022open} Holds)}\tablefootnote{For bounds not requiring this conjecture, see Section \ref{sec:ucb_disc}.} & \\
		Chains & $\begin{cases}
			O(2^mBL^{m-1}T^{1-\nu/d})&\text{ when }d>\nu,\\
			{\widetilde{O}}(2^mBL^{m-1})&\text{ when }d\le\nu.
		\end{cases}$
		\\
		Multi-Output Chains & $\begin{cases}
			O(5^mBL^{m-1}T^{1-\nu/d_{\max}})&\text{ when }d_{\max}>\nu,\\
			{\widetilde{O}}(5^mBL^{m-1})&\text{ when }d_{\max}\le\nu.
		\end{cases}$
		\\
		Feed-Forward Networks & $\begin{cases}
			O(2^m\sqrt{D_{2,m}}BL^{m-1}T^{1-\nu/d_{\max}})&\text{ when }d_{\max}>\nu,\\
			{\widetilde{O}}(2^m\sqrt{D_{2,m}}BL^{m-1})&\text{ when }d_{\max}\le\nu.
		\end{cases}$
		\\
		\midrule
		\textbf{Algorithm-Independent Simple Regret Lower Bound} & \\
		Chains/Multi-Output Chains/Feed-Forward Networks & $\Omega(B(cL)^{m-1}T^{-\nu/d})$\hspace{4ex}when $\nu\ge1.$ \\ [3mm]
		\textbf{Algorithmic Simple Regret Upper Bound} & \\
		Chains/Multi-Output Chains & $\begin{cases}
			{\widetilde{O}}(BL^{m-1}T^{-\nu^m/d})&\text{when $\nu\le 1$,}\\
			{\widetilde{O}}(B^{1+\nu+\nu^2+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu/d})&\text{when $\nu> 1$.}
		\end{cases}$ \\
	   Chains/Multi-Output Chains (Restrictive Cases) & $\begin{cases}
	   	O(BL^{m-1}T^{-\nu/d})&\text{ when $\nu\le 1$,}\\
	   	O(BL^{(m-1)\nu}T^{-\nu/d})&\text{ when $\nu> 1$.}
	   \end{cases}$ \\
   		Feed-Forward Networ
   		ks &   $\begin{cases}
   			{\widetilde{O}}(\sqrt{D_{2,m}}BL^{m-1}T^{-\nu^m/d})&\text{ when $\nu\le 1$,}\\
   			{\widetilde{O}}( {\widetilde{D}}_{2,m}^\nu B^{1+\nu+\nu^2+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu/d} )&\text{ when $\nu> 1$.}
   		\end{cases}$ \\
   		Feed-Forward Networks (Restrictive Cases) & $\begin{cases}
   		O(\sqrt{D_{2,m}}BL^{m-1}T^{-\nu/d})&\text{ when $\nu\le 1$,}\\
   		O\big((D_{2,m})^{\nu/2}BL^{(m-1)\nu}T^{-\nu/d}\big)&\text{ when $\nu> 1$.}
   		\end{cases}$ \\
		\bottomrule
	\end{tabular}
	\caption{Summary of regret bounds for the Mat\'ern kernel. (The algorithmic simple regret upper bounds are valid when the domain of each layer is a hyperrectangle.) $T$ denotes the time horizon; $m$ denotes the number of layers; $B$ denotes the RKHS norm upper bound of each layer; $L$ denotes the Lipschitz constant upper bound of each layer; $\nu$ denotes the smoothness of the Mat\'ern kernel; $d$ denotes the domain dimension of $g$; $d_{\max}=\max_{i\in[m]}d_i$ denotes the maximum dimension among all the $m$ layers; $D_{2,m}=\prod_{i=2}^m d_i$ denotes the product of the dimensions from the second to the last layer; ${\widetilde{D}}_{2,m}^\nu=\prod_{i=2}^m (d_i)^{\nu^{m+1-i}/2}$. The lower bounds hold for some $c=\Theta(1)$.}
	\label{tab:summary}
\end{table}

\section{Comparison to \cite{kusakawa2021bayesian}}
\label{sec:comparison}
In this section, we compare our work to \cite{kusakawa2021bayesian}, which proposed two confidence bound based algorithms cascade UCB (cUCB) and optimistic improvement (OI) under the noise-free setting. Both algorithms utilize a novel posterior standard deviation defined using the Lipschitz constant.

\textbf{Problem Setup.}
\begin{itemize}
	\item In each layer, \cite{kusakawa2021bayesian} assumes the entries in the same layer are mutually independent, which is equivalent to our feed-forward network structure.
	\item Different from our assumption of $\ff{i,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L, \|\cdot\|_2)$ for each $i\in[m]$ and $j\in[d_{i+1}]$, where the Lipschitz continuity is with respect to $\|\cdot\|_2$, \cite{kusakawa2021bayesian} assumes $\ff{i,j}\in\mathcal{H}_k(B)\cap\mathcal{F}(L_f, \|\cdot\|_1)$, where the Lipschitz continuity is with respect to $\|\cdot\|_1$. Therefore, for a fixed scalar-valued function $\ff{i,j}\in\mathcal{F}(L, \|\cdot\|_2)\cap\mathcal{F}(L_f, \|\cdot\|_1)$ with the smallest possible $L$ and $L_f$, it is satisfied that $L_f\le L\le\sqrt{d_i}L_f$.
	\item \cite{kusakawa2021bayesian} has an additional assumption that $\sig{i,j}_t\in\mathcal{F}(L_\sigma,\|\cdot\|_1)$ for all $t\ge 1$ for some constant $L_\sigma$. The constant $L_\sigma$ is not used in the algorithms, and only appears in the regret bounds.
\end{itemize}

\textbf{Regret Bounds. }
\begin{itemize}
	\item cUCB achieves cumulative regret $R_T=O\big(B(BL_\sigma+L_f)^{m-1}D_{2,m}(\sum_{i=2}^m d_i) \sqrt{T\gamma_T}\big)$, with a significantly larger $T$-independent factor than that in $R_T=O(2^mBL^{m-1}\sqrt{D_{2,m}T\gamma_T})$ of our GPN-UCB (\cref{algo:ucb}).
	\item OI achieves simple regret $r^\ast_T=O\big(m^{m^2+m+\frac{1}{2}}B^{m+1}L_f^{m^3}(BL_\sigma+L_f)^{3m^3}(D_{2,m})^{2m^2}(\sum_{i=2}^m d_i)^{m^2+1}T^{-\frac{\nu}{2\nu+d}}\big)$ for the Mat\'ern kernel with smoothness $\nu$. When $\nu\le 1$, our non-adaptive sampling (\cref{algo:fd}) achieves $r^\ast_T={\widetilde{O}}(\sqrt{D_{2,m}}BL^{m-1}T^{-\nu^m/d})$, which has a significantly smaller $T$-independent factor. When $\nu>1$, our \cref{algo:fd} achieves $r^\ast_T={\widetilde{O}}( {\widetilde{D}}_{2,m}^\nu B^{1+\nu+\nu^2+\nu^{m-2}}L^{\nu^{m-1}}T^{-\nu/d} )$, which has a smaller $T$-dependent factor.  The $T$-independent factors are more difficult to compare, though ours certainly become preferable under the more restrictive scenarios discussed leading up to \eqref{eq:more_restrictive}.
\end{itemize}

\textbf{Generality.} \cite{kusakawa2021bayesian} allows additional (multi-dimensional) input independent of previous layers for each layer. GPN-UCB and non-adaptive sampling can also be adapted to accept additional input. Let $\widehat{\mathcal{X}}^{(i)}\subset{\mathbb{R}}^{{\widehat{d}}_i}$ denote the domain of the additional input of $\ff{i}$, and let $\widehat{{\mathbf{x}}}_{[i:j]}=[\widehat{{\mathbf{x}}}^{(i)},\dots,\widehat{{\mathbf{x}}}^{(j)}]$ denote the concatenation of the addition inputs from $\ff{i}$ to $\ff{j}$.
\begin{itemize}
	\item For GPN-UCB (\cref{algo:ucb}), we simply modify the upper confidence bound to be
	\begin{align}
		\mathrm{UCB}_t({\mathbf{x}},\widehat{{\mathbf{x}}})&=\max_{{\mathbf{z}}\in\Delta^{(m)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:m-1]})}\UCBBB{m,1}_t([{\mathbf{z}},\widehat{{\mathbf{x}}}^{(m)}]),
	\end{align}
	where
	\begin{align}
		\Delta^{(1)}_t({\mathbf{x}})=&\{{\mathbf{x}}\},\\
		\Delta^{(i+1,j)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i]})=&\Bigg[\min_{{\mathbf{z}}\in\Delta^{(i)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i-1]})}\LCBBB{i,j}_t([{\mathbf{z}},\widehat{{\mathbf{x}}}^{(i)}]) , \max_{{\mathbf{z}}\in\Delta^{(i)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i-1]})}\UCBBB{i,j}_t([{\mathbf{z}},\widehat{{\mathbf{x}}}^{(i)}]) \Bigg]\nonumber\\
		&\hspace{42ex}\text{ for }i\in[m-1],j\in[d_{i+1}],\\
		\Delta^{(i)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i-1]})=& \Delta^{(i,1)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i-1]})\times\cdots\times\Delta^{(i,d_i)}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i-1]}) \nonumber\\
		&\hspace{56ex}\text{ for }i\in[m].
	\end{align}
	The cumulative regret upper bound will remain the same, since $\dia{\del{i}_t({\mathbf{x}},\widehat{{\mathbf{x}}}_{[2:i-1]})}$ remains the same.
	\item For non-adaptive sampling (\cref{algo:fd}), with ${\widehat{d}}=d+\sum_{i=2}^m {\widehat{d}}_i$, we choose $\{{\mathbf{x}}_s, \widehat{{\mathbf{x}}}^{(2)}_s,\dots,\widehat{{\mathbf{x}}}^{(m)}_s\}_{s=1}^T$ such that
	\begin{align}
		\max_{[{\mathbf{x}},\widehat{{\mathbf{x}}}^{(2)},\dots,\widehat{{\mathbf{x}}}^{(m)}]\in \mathcal{X}\times\widehat{\mathcal{X}}^{(2)}\times\cdots\times\widehat{\mathcal{X}}^{(m)}} \min_{s\in[T]} \| [{\mathbf{x}},\widehat{{\mathbf{x}}}^{(2)},\dots,\widehat{{\mathbf{x}}}^{(m)}] - [{\mathbf{x}}_s,\widehat{{\mathbf{x}}}^{(2)}_s,\dots,\widehat{{\mathbf{x}}}^{(m)}_s] \|_2=O(T^{-1/{\widehat{d}}}).
	\end{align}
	Then, modify the composite mean to be
	\begin{align}
		\mu_T^g({\mathbf{x}},\widehat{{\mathbf{x}}})=\muu{m,1}_T([\zz{m},\widehat{{\mathbf{x}}}^{(m)}])
	\end{align}
	with
	\begin{align}
		\zz{1}&={\mathbf{x}},\\
		\zz{i+1,j}&=\muu{i,j}_T([\zz{i},\widehat{{\mathbf{x}}}^{(i)}])&&\text{ for }i\in[m-1],j\in[d_{i+1}],\\
	\zz{i+1}&=\muu{i}_T([\zz{i},\widehat{{\mathbf{x}}}^{(i)}])=[\zz{i+1,1},\dots,\zz{i+1,d_i}]&&\text{ for }i\in[m-1].
	\end{align}
	Then, the only change in simple regret upper bound is that $d$ is replaced by ${\widehat{d}}$.
\end{itemize}