\section*{Acknoledements}
This work was supported by MEXT KAKENHI 20H00601, 21H03498, 22H00300, JSPS KAKENHI Grant Number JP21J14673.

\section{Introduction}
\label{sec:intro}

\emph{Bayesian optimization} (BO) \citep{Mockus1978-Application} has become a widely-used framework for expensive black-box optimization problems.
To reduce the number of function evaluations, BO sequentially observes the noisy function value using an acquisition function computed from a Bayesian model.
BO has been applied to many different fields, including automatic machine learning \citep{Snoek2012-Practical,kandasamy2018-neural}, materials informatics \citep{ueno2016combo}, and drug design \citep{korovina2020-chembo,griffiths2020-constrained}.

Theoretical guarantees for BO have also been studied extensively \citep[e.g., ][]{Srinivas2010-Gaussian, Russo2014-learning}.
The Gaussian process upper confidence bound (GP-UCB) \citep{Srinivas2010-Gaussian}, which achieves a strong theoretical guarantee for cumulative regret bounds, is a seminal work in this field.
Based on this versatile framework for the analysis, while maintaining the theoretical guarantee, GP-UCB has been extended to various problems, which include multi-objective BO \citep{paria2020-flexible,Zuluga2016-epsilon}, multi-fidelity BO \citep{Kandasamy2016-Gaussian,Kandasamy2017-Multi}, parallel BO \citep{Contal2013-Parallel,Desautels2014-Parallelizing}, high-dimensional BO \citep{kandasamy2015-high}, and cascade BO \citep{Kusakawa2022-bayesian}.

However, the choice of a \emph{confidence parameter} $\beta_t$, which controls the exploitation-exploration trade-off, is a practical challenge for GP-UCB-based methods, where $t$ is the iteration of BO.
Since the theoretical choice of $\beta_t$ is considerably large and increases with $t$, particularly in the later iterations, GP-UCB with the theoretical $\beta_t$ causes over-exploration, which reduces the optimization performance.
Therefore, in practice, $\beta_t$ is often specified manually using some heuristics, which strongly affects the performance of GP-UCB.
Since this is a common problem with GP-UCB-based methods, mitigating the effect of the manual specification of $\beta_t$ is an important task.

To resolve this issue, \citet{berk2021-randomized} have proposed a \emph{randomized} GP-UCB (RGP-UCB), which uses a randomized confidence parameter $\zeta_t$, which follows the Gamma distribution.
Their experimental results demonstrated that RGP-UCB performed better than the standard GP-UCB.
Although \citet{berk2021-randomized} provided the regret analysis, it contains some technical issues, including an unknown convergence rate (See Appendix~\ref{app:issues_RGPUCB} for details).
Therefore, the theoretical justification of RGP-UCB is still an open problem.

\begin{table*}[t]
    \centering
    \begin{tabular}{c|c|c|c}
         & GP-UCB \citep{Srinivas2010-Gaussian} & RGP-UCB & IRGP-UCB (proposed) \\ \hline
         BCR (discrete) & $O(\sqrt{T \gamma_T \log T})$ & $O(\sqrt{T \gamma_T \log T})$ & * $O(\sqrt{T \gamma_T})$ \\
         BCR (continuous) & $O(\sqrt{T \gamma_T \log T})$ & $O(\sqrt{T \gamma_T \log T})$ & $O(\sqrt{T \gamma_T \log T})$ \\
         Sufficient condition  & \multirow{2}{*}{$\sqrt{ \gamma_T \log T / T} \lesssim \eta$} &  \multirow{2}{*}{$\sqrt{\gamma_T \log T / T} \lesssim \eta$} &  \multirow{2}{*}{* $\sqrt{\gamma_T  / T} \lesssim \eta / \sqrt{- \log \eta}$}\\
         of BSR $< \eta$ &&&
    \end{tabular}
    \caption{
        Summary of Bayesian regret of GP-UCB-based algorithms.
        %
        The first and second rows show the BCR bounds for discrete and continuous input domains, respectively, where $\gamma_T$ is the maximum information gain defined in Section~\ref{sec:background}.
        %
        The third row shows the sufficient conditions to achieve that BSR is lower than the predefined accuracy $\eta \in (0, 1)$ for both continuous and discrete input domains, where $\lesssim$ represents an inequality ignoring factors except for $T$ and $\eta$.
        %
        Star means better bounds than known results.
        %
        Note that the results in the right two columns are derived in this paper.
        }
    \label{tab:regret_summary}
\end{table*}

This study demonstrates a generalized regret bound of RGP-UCB, which holds for a wider class of distributions for $\zeta_t$ in the Bayesian setting, where an objective function $f$ is assumed to be a sample path from the GP.
However, this analysis requires an increase in $\zeta_t$, as $\mathbb{E} [\zeta_t] \approx \beta_t$, which implies the over-exploration cannot be avoided.
Furthermore, this generalization does not provide us the choice of the distribution of $\zeta_t$.

The aforementioned analyses motivated us to propose an improved RGP-UCB (IRGP-UCB), in which $\zeta_t$ follows a \emph{two-parameter exponential distribution}.
First, we show the sub-linear Bayesian cumulative regret (BCR) bounds for IRGP-UCB in the Bayesian setting.
In particular, when the input domain is finite, the BCR for IRGP-UCB achieves the best convergence rate, which shows a $O(\sqrt{\log T})$ multiplicative factor improvement from the known bounds, where $T$ is the entire time horizon.
Unfortunately, for the infinite domain, the convergence rate of IRGP-UCB is similar to that of the original GP-UCB and RGP-UCB.
On the other hand, IRGP-UCB achieves a tighter convergence rate with respect to $T$ for the Bayesian simple regret (BSR) than existing analyses for both finite and infinite domains.
These analyses also reveal that the increase in $\zeta_t$ for IRGP-UCB is unnecessary.
Therefore, the over-exploration of the original GP-UCB is theoretically alleviated via randomization using the two-parameter exponential distribution.
A key finding in our proof is a direct upper bound of $\mathbb{E}[\mathop{\rm argmax}\limits f(\*x)]$.
We show that $\mathbb{E}[\mathop{\rm argmax}\limits f(\*x)]$ can be bounded by the expectation of the UCB with our random non-increasing confidence parameter, which enables us to derive the tighter Bayesian regret bounds than known results.

Our main contributions are summarized as follows:
\begin{itemize}
    \item We provide the Bayesian regret bounds for RGP-UCB, which holds for a wider class of distributions for $\zeta_t$ and achieves the same convergence rate as GP-UCB.
    \item
        We propose yet another randomized variant of GP-UCB called IRGP-UCB, which sets the confidence parameter $\zeta_t$ using the two-parameter exponential distribution.
        %
        A notable advantage of IRGP-UCB is that increases in the confidence parameter in proportion to $t$ are unnecessary.
    \item
        We show Bayesian regret bounds for IRGP-UCB, which achieves the best convergence rate, in which $O(\sqrt{\log T})$ multiplicative factor is improved from known results.
        %
        This result suggests that the over-exploration of GP-UCB is theoretically alleviated.
    \item We provide the upper bound for $\mathbb{E}[\mathop{\rm argmax}\limits f(\*x)]$ along the way of the proof, which enables us to improve the Bayesian regret bounds.
\end{itemize}
The theoretical results are summarized in Table~\ref{tab:regret_summary}.
Finally, we demonstrate the effectiveness of IRGP-UCB through a wide range of experiments.

\subsection{Related Work}
\label{sec:related_work}

This study considers BO with the Bayesian setting, where the objective function $f$ is assumed to be a sample path from the GP.
Various BO methods have been developed in the literature, for example, expected improvement (EI) \citep{Mockus1978-Application}, entropy search (ES) \citep{Henning2012-Entropy}, and predictive entropy search (PES) \citep{Hernandez2014-Predictive}.
Although the regret analysis of EI for the noiseless and frequentist setting, in which $f$ is an element of reproducing kernel Hilbert space, is provided in \citep{bull2011convergence}, the Bayesian setting has not been considered.
Further, although the practical performance of ES and PES has been shown repeatedly, their regret analysis is an open problem.

GP-UCB is one of the prominent studies for theoretically guaranteed BO methods.
\citet{Srinivas2010-Gaussian} showed the high probability bound for cumulative regret.
Although Bayesian regret analysis for GP-UCB has not been shown explicitly, \citet{paria2020-flexible} have shown the BCR bound for a multi-objective extension of GP-UCB, which contains that of the original single-objective GP-UCB as the special case.
For completeness, we show the BCR bound of GP-UCB in Section~\ref{sec:GPUCB}.
Although many studies \citep[e.g., ][]{Srinivas2010-Gaussian,Chowdhury2017-on,janz2020-bandit} considered the frequentist setting, this study concentrates on the Bayesian setting.

\citet{berk2021-randomized} attempted to alleviate the GP-UCB hyperparameter tuning issue through randomization using the Gamma distribution.
Note that our result shown in Theorem~\ref{theo:BCR_RGPUCB_informal} differs from \citep{berk2021-randomized}.
We believe that these differences come from several technical issues in the proof in \citep{berk2021-randomized} (See Appendix~\ref{app:issues_RGPUCB} for details).
In addition, the final convergence rate in \citep{berk2021-randomized} includes the variables whose convergence rate is not proved.
In contrast, our Theorem~\ref{theo:BCR_RGPUCB_informal} fully clarifies the convergence rate without those unproved variables.
Moreover, Theorem~\ref{theo:BCR_RGPUCB_informal} is generalized to a wider class of distributions for confidence parameters that contain the Gamma distribution.
In addition, the two-parameter exponential distribution used in IRGP-UCB is not considered in \citep{berk2021-randomized}.

Another standard BO method with a theoretical guarantee is Thompson sampling (TS) \citep{Russo2014-learning,Kandasamy2018-Parallelised}.
TS achieves the sub-linear BCR bounds by sequentially observing the optimal point of the GP's posterior sample path.
Although TS does not require any hyperparameter, TS is often deteriorated by over-exploration, as discussed in \citep{Shahriari2016-Taking}.

\citet{Wang2016-Optimization} have shown the regret analysis of the GP estimation (GP-EST) algorithm, which can be interpreted as GP-UCB with the confidence parameter defined using $\hat{m}$, an estimator of $\mathbb{E}[\mathop{\rm argmax}\limits f(\*x)]$.
Their analysis requires an assumption $\hat{m} \geq \mathbb{E}[\mathop{\rm argmax}\limits f(\*x)]$, whose sufficient condition is provided in \citep[][Corollary~3.5]{Wang2016-Optimization}.
However, this sufficient condition does not typically hold, as discussed immediately after the corollary in \citep{Wang2016-Optimization}.
Furthermore, the final convergence rate contains $\hat{m}$ itself, whose convergence rate is not clarified.
In contrast, our Lemma~\ref{lem:bound_RGPUCB} shows the bound for $\mathbb{E}[\mathop{\rm argmax}\limits f(\*x)]$ under common regularity conditions.
\citet{Wang2017-Max} have shown the regret analysis of max-value entropy search (MES).
However, it is pointed out that their proof contains several technical issues \citep{takeno2022-sequential}.

\section{Background}
\label{sec:background}

\subsection{Bayesian Optimization}

We consider an optimization problem $\*x^* = \mathop{\text{argmax}}_{\*x \in {\cal X}} f(\*x)$, where $f$ is an unknown expensive-to-evaluate objective function, ${\cal X} \subset \mathbb{R}^d$ is an input domain, and $d$ is an input dimension.
BO sequentially observes the noisy function value aiming to minimize the number of function evaluations.
Thus, at each iteration $t$, we can query $\*x_t$ and obtain $y_t = f(\*x_t) + \epsilon_t$, where $\epsilon_t \sim {\cal N}(0, \sigma^2)$ is i.i.d. Gaussian noise with a positive variance $\sigma^2 > 0$.

We assume that $f$ is a sample path from a GP \citep{Rasmussen2005-Gaussian} with a zero mean and a stationary kernel function $k: {\cal X} \times {\cal X} \mapsto \mathbb{R}$ denoted as $f \sim {\mathscr G} {\cal P}(0, k)$.
At each iteration $t$, a dataset ${\cal D}_{t-1} \coloneqq \{ (\*x_i, y_i) \}_{i=1}^{t-1}$ is already obtained from the nature of BO.
Then, the posterior distribution $p(f \mid {\cal D}_{t-1})$ is a GP again.
The posterior mean and variance of $f(\*x)$ are derived as follows:
\begin{align*}
    \mu_{t-1}(\*x) &= \*k_{t-1}(\*x)^\top \bigl(\*K + \sigma^2 \*I_{t-1} \bigr)^{-1} \*y_{t-1}, \\
    \sigma_{t-1}^2 (\*x) &= k(\*x, \*x) - \*k_{t-1}(\*x) ^\top \bigl(\*K + \sigma^2 \*I_{t-1} \bigr)^{-1} \*k_{t-1}(\*x),
\end{align*}
where $\*k_{t-1}(\*x) \coloneqq \bigl( k(\*x, \*x_1), \dots, k(\*x, \*x_{t-1}) \bigr)^\top \in \mathbb{R}^{t-1}$, $\*K \in \mathbb{R}^{(t-1)\times (t-1)}$ is the kernel matrix whose $(i, j)$-element is $k(\*x_i, \*x_j)$, $\*I_{t-1} \in \mathbb{R}^{(t-1)\times (t-1)}$ is the identity matrix, and $\*y_{t-1} \coloneqq (y_1, \dots, y_{t-1})^\top \in \mathbb{R}^{t-1}$.
Hereafter, we denote that the probability density function (PDF) $p(\cdot \mid {\cal D}_{t-1}) = p_t (\cdot)$, the probability $\Pr (\cdot \mid {\cal D}_{t-1}) = \mathrm{Pr}_t (\cdot) $, and the expectation $\mathbb{E}[\cdot \mid {\cal D}_{t-1}] = \mathbb{E}_t[\cdot]$ for brevity.

\subsection{Preliminaries for Regret Analysis}

When ${\cal X}$ is infinite (continuous), the following regularity assumption is used in most analyses \citep[e.g., ][]{Srinivas2010-Gaussian,Kandasamy2018-Parallelised,paria2020-flexible}:
\begin{assumption}
    Let ${\cal X} \subset [0, r]^d$ be a compact and convex set, where $r > 0$.
    %
    Assume that the kernel $k$ satisfies the following condition on the derivatives of a sample path $f$.
    %
    There exists the constants $a, b > 0$ such that,
    \begin{align*}
        \Pr \left( \sup_{\*x \in {\cal X}} \left| \frac{\partial f}{\partial \*x_j} \right| > L \right) \leq a \exp \left( - \left(\frac{L}{b}\right)^2 \right),\text{ for } j \in [d],
    \end{align*}
    where $[d] = \{1, \dots, d\}$.
    \label{assump:continuous_X}
\end{assumption}
Our analysis also requires this assumption.

The convergence rates of regret bounds are characterized by \emph{maximum information gain} (MIG) \citep{Srinivas2010-Gaussian,vakili2021-information}.
MIG $\gamma_T$ is defined as follows:
\begin{definition}[Maximum information gain]
    Let $f \sim {\mathscr G} {\cal P} (0, k)$.
    %
    Let $A = \{ \*a_i \}_{i=1}^T \subset {\cal X}$.
    %
    Let $\*f_A = \bigl(f(\*a_i) \bigr)_{i=1}^T$, $\*\epsilon_A = \bigl(\epsilon_i \bigr)_{i=1}^T$, where $\forall i, \epsilon_i \sim {\cal N}(0, \sigma^2)$, and $\*y_A = \*f_A + \*\epsilon_A \in \mathbb{R}^T$.
    %
    Then, MIG $\gamma_T$ is defined as follows:
    \begin{align*}
        \gamma_T \coloneqq \mathop{\text{argmax}}_{A \subset {\cal X}; |A| = T} I(\*y_A ; \*f_A),
    \end{align*}
    where $I$ is the Shanon mutual information.
\end{definition}
MIG is known to be sub-linear for commonly used kernel functions, e.g., $\gamma_T = O\bigl( (\log T)^{d+1} \bigr)$ for RBF kernels and $\gamma_T = O\bigl( T^{\frac{d}{2\nu + d}} (\log T)^{\frac{2\nu}{2\nu + d}} \bigr)$ for Mat\`{e}rn-$\nu$ kernels \citep{Srinivas2010-Gaussian,vakili2021-information}.

\subsection{Bayesian Regret Analysis}

In this paper, we evaluate the performance of BO methods by \emph{Bayesian regret} \citep{Russo2014-learning,Kandasamy2018-Parallelised,paria2020-flexible}.
The BSR and BCR are defined as follows:
\begin{align}
    {\rm BSR}_T \coloneqq \mathbb{E}\bigl[ f(\*x^*) - \mathop{\rm argmax}\limits_{t \leq T} f(\*x_t) \bigr], \label{eq:BSR} \\
    {\rm BCR}_T \coloneqq \mathbb{E}\left[\sum_{t=1}^T f(\*x^*) - f(\*x_t) \right], \label{eq:BCR}
\end{align}
where the expectation is taken with all randomness, including $f, \epsilon_t$, and the randomness of BO algorithms.
We discuss the convergence of BSR by analyzing the BCR.
That is, if BCR is sub-linear, BSR converges to zero since ${\rm BSR}_T \leq \frac{1}{T} {\rm BCR}_T$.
Furthermore, we evaluate a required number of function evaluations to achieve ${\rm BSR}_T \leq \eta$, where $\eta > 0$ represents the desired accuracy of a solution.

According to \citep{Russo2014-learning}, Bayesian regret bounds directly provide high probability bounds.
Assume ${\rm BCR}_T = O(h(T))$ with some non-negative function $h$.
Then, the direct consequence of Markov's inequality implies that $\sum_{t=1}^T f(\*x^*) - f(\*x_t) = O \left( h(T) / \delta \right)$ holds with probability $1 - \delta$, where $\delta \in (0, 1)$.
Thus, the improvement of the convergence rate of BCR for $T$ shows that of high-probability bound although the rate for $\delta$ is worse compared to $\log (1/\delta)$ \citep{Srinivas2010-Gaussian}.

The existing Bayesian regret analyses typically use the technique called \emph{regret decomposition} \citep{Russo2014-learning}.
This approach decomposes BCR as follows:
\begin{align}
    {\rm BCR}_T = \sum_{t=1}^T \mathbb{E}\left[ f(\*x^*) - u_t(\*x_t) \right] + \mathbb{E}\left[ u_t(\*x_t) - f(\*x_t) \right],
    \label{eq:regret_depomsition}
\end{align}
where $u_t(\*x) \coloneqq \mu_{t-1}(\*x) + \beta^{1/2}_t \sigma_{t-1}(\*x)$.
If we use fixed constant for $\{\beta_t\}_{t \geq 1}$ (e.g., $\forall t\geq 1, \beta_t = 2$), bounding the first term $\sum_{t=1}^T \mathbb{E}\left[ f(\*x_*) - u_t(\*x_t) \right]$ is difficult.
Thus, in existing analyses, $\{\beta_t\}_{t \geq 1}$ is scheduled so that $\beta_t \propto \log t$.
Since $\beta_t$ remains in the upper bound, the regret decomposition using $u_t(\*x)$ deteriorates the final convergence rate, as we will show an example of GP-UCB in Theorem~\ref{theo:BCR_GPUCB_informal}.
IRGP-UCB avoids this problem by randomizing the confidence parameters, as shown in Section~\ref{sec:theorem}.

\subsection{GP-UCB}
\label{sec:GPUCB}

The GP-UCB \citep{Srinivas2010-Gaussian} selects the next evaluation by maximizing UCB as follows:
\begin{align*}
    \*x_t = \mathop{\text{argmax}}_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_t \sigma_{t-1}(\*x).
\end{align*}
Although BCR bound for the standard GP-UCB has not been shown explicitly, sub-linear regret bounds of GP-UCB can be shown in nearly the same way as \citep{paria2020-flexible}:
\begin{theorem}[Informal: BCR of GP-UCB]
    Suppose that $f \sim {\mathscr G} {\cal P}(0, k)$, where $k(\*x, \*x) = 1$.
    %
    Assume that ${\cal X}$ is finite or Assumption~\ref{assump:continuous_X} holds.
    %
    Then, by running GP-UCB with $\beta_t \propto \log(t)$, BCR can be bounded as follows:
    \begin{align}
        {\rm BCR}_T = O(\sqrt{T \beta_T \gamma_T}),
    \end{align}
    which implies sub-linearity.
    %
    \label{theo:BCR_GPUCB_informal}
\end{theorem}
See Appendix~\ref{app:GP_UCB_proof} for the proof\footnote{It is worth noting that the convergence rate for $a$ in Assumption~\ref{assump:continuous_X} is tightened in our proof compared with the proofs in \citep{paria2020-flexible,Kandasamy2018-Parallelised}. This is applied to all subsequent regret analyses for a continuous domain. See details for Appendixes~\ref{app:GP_UCB_proof} and \ref{app:lemmas}.}.
In the theorem, the confidence parameter $\beta_t$ must be scheduled as $\beta_t \propto \log(t)$.

\section{Algorithm}
\label{sec:algorithm}

Algorithm \ref{alg:RGPUCB} shows a simple procedure of IRGP-UCB.
The difference from the original GP-UCB is that $\zeta_t$ is a random variable, not a constant.
Using generated $\zeta_t$, the next evaluation is chosen as follows:
\begin{align*}
    \*x_t = \mathop{\text{argmax}}_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2}_t \sigma_{t-1}(\*x).
\end{align*}
Our choice of the distribution of $\zeta_t$ is a \emph{two-parameter exponential distribution} \citep[e.g., ][]{beg1980estimation,lam1994estimation}, whose PDF is written as follows:
\begin{align*}
    p(\zeta; s, \lambda) =
    \begin{cases}
    \lambda \exp\bigl( \lambda (\zeta - s) \bigr) & \text{if } \zeta \geq s, \\
    0 & \text{otherwise},
    \end{cases}
\end{align*}
where $\lambda$ is a rate parameter.
This can be seen as a distribution of the sum of $s$ and $Z \sim {\rm Exp} (\lambda)$.
Thus, the sampling from the two-parameter exponential distribution can be easily performed using an exponential distribution as follows:
\begin{align}
    \zeta_t \gets s + Z, \text{ where } Z \sim {\rm Exp} (\lambda).
    \label{eq:sampling_beta}
\end{align}
The theoretical choice of $s$ and $\lambda$ will be shown in Section~\ref{sec:theorem}.

\begin{algorithm}[!t]
    \caption{IRGP-UCB}\label{alg:RGPUCB}
    \begin{algorithmic}[1]
        \Require Input space ${\cal X}$, Parameters $s$ and $\lambda$ for $\zeta_t$, GP prior $\mu=0$ and $k$
        \State ${\cal D}_{0} \gets \emptyset $
        \For{$t = 1, \dots$}
            \State Fit GP to ${\cal D}_{t-1}$
            \State Generate $\zeta_t$ by Eq.~\eqref{eq:sampling_beta}
            \State $\*x_t \gets \mathop{\text{argmax}}_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta_t^{1/2} \sigma_{t-1}(\*x)$
            \State Observe $y_t = f(\*x_t) + \epsilon_t$ and ${\cal D}_{t} \gets {\cal D}_{t-1} \cup (\*x_t, y_t)$
        \EndFor
    \end{algorithmic}
\end{algorithm}

\section{Regret Analysis}
\label{sec:theorem}

First, we show the Bayesian regret analysis for a general RGP-UCB, which is not restricted to the two-parameter exponential distribution.
Next, we show the tight Bayesian regret bounds for IRGP-UCB.

\subsection{Regret Bound for General RGP-UCB}

Here, we provide a generalized theoretical analysis for RGP-UCB, which holds for a wider class of distributions.
\begin{theorem}[Informal: BCR of RGP-UCB]
    Suppose that $f \sim {\mathscr G} {\cal P}(0, k)$, where $k(\*x, \*x) = 1$.
    %
    Assume that ${\cal X}$ is finite or Assumption~\ref{assump:continuous_X} holds.
    %
    Let $\{\zeta_t\}_{t \geq 1}$ be a sequence of non-negative random variables, which satisfies $\mathbb{E}[\zeta_t] = O(\log t)$ and $\mathbb{E}[\exp ( - \zeta_t / 2)] = O (1 / t^2)$.
    %
    By running RGP-UCB with $\zeta_t$, BCR can be bounded as follows:
    \begin{align*}
        {\rm BCR}_T = O(\sqrt{T \gamma_T \mathbb{E}[\zeta_T]}),
    \end{align*}
    which implies sub-linearity.
    %
    \label{theo:BCR_RGPUCB_informal}
\end{theorem}
See Appendix~\ref{app:proof_general_RGPUCB} for the proof, in which we also rectify the technical issues in \citep{berk2021-randomized} listed in Appendix~\ref{app:issues_RGPUCB}.
A notable difference between Theorem~\ref{theo:BCR_RGPUCB_informal} and \citep[][Theorem~3]{berk2021-randomized} is that Theorem~\ref{theo:BCR_RGPUCB_informal} is applicable to a wider class of distributions of $\zeta_t$, not only Gamma distribution proposed in \citep{berk2021-randomized}.
Roughly speaking, Theorem~\ref{theo:BCR_RGPUCB_informal} only requires that the distribution has the parameters that can control the distribution so that $\mathbb{E}[\zeta_t] = O(\log t)$ and $\mathbb{E}[\exp ( - \zeta_t / 2)] = O (1 / t^2)$.
This condition can be satisfied by many distributions, such as Gamma, two-parameter exponential, and truncated normal distributions, as shown in Appendix~\ref{app:examples_RGPUCB_dist}.
See Appendix~\ref{app:examples_RGPUCB_dist} for a detailed discussion.

Although Theorem~\ref{theo:BCR_RGPUCB_informal} achieves sub-linear BCR bounds, the implication of Theorem~\ref{theo:BCR_RGPUCB_informal} is restricted.
In the proof, $\zeta_t$ must behave similar to $\beta_t$, that is, $\mathbb{E}[\zeta_t] \propto \beta_t$ and $\mathbb{E}[\exp ( - \zeta_t / 2)] \propto \exp ( - \beta_t / 2)$.
Therefore, based on this theorem, over-exploration of the GP-UCB cannot be avoided.
Furthermore, this theorem does not specify a distribution, that is suitable for the randomization of GP-UCB.

\subsection{Regret Bounds for IRGP-UCB}

First, for completeness, we provide a modified version of the basic lemma derived in \citep{Srinivas2010-Gaussian}:
\begin{lemma}
    Suppose that $f$ is a sample path from a GP with zero mean and a stationary kernel $k$ and ${\cal X}$ is finite.
    %
    Pick $\delta \in (0, 1)$ and $t \geq 1$.
    %
    Then, for any given ${\cal D}_{t-1}$,
    \begin{align*}
        \mathrm{Pr}_t \left( f(\*x) \leq \mu_{t-1}(\*x) + \beta^{1/2}_{\delta} \sigma_{t-1}(\*x), \forall \*x \in {\cal X} \right)
        \geq 1 - \delta,
    \end{align*}
    where $\beta_{\delta} = 2 \log (|{\cal X}| /( 2 \delta))$.
    \label{lem:bound_srinivas}
\end{lemma}
This lemma differs slightly from \citep[][Lemma~5.1]{Srinivas2010-Gaussian}, since, in Lemma~\ref{lem:bound_srinivas}, the iteration $t$ is fixed, and $\beta_{\delta}$ does not depend on $t$.

Based on Lemma~\ref{lem:bound_srinivas}, we show the following key lemma to show the improved convergence rate:
\begin{lemma}
    Let $f \sim {\mathscr G} {\cal P} (0, k)$, where $k$ is a stationary kernel and $k(\*x, \*x) = 1$, and ${\cal X}$ be finite.
    %
    Assume that $\zeta$ follows a two-parameter exponential distribution with $s = 2 \log (|{\cal X}| / 2)$ and $\lambda = 1/2$.
    %
    Then, the following inequality holds:
    \begin{align*}
        \mathbb{E}[f(\*x^*)] \leq \mathbb{E} \left[\mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2}_t \sigma_{t-1}(\*x) \right],
    \end{align*}
    for all $t \geq 1$.
    \label{lem:bound_RGPUCB}
\end{lemma}
\begin{proof}
    We here show the short proof of Lemma~\ref{lem:bound_RGPUCB} although detailed proof is shown in Appendix~\ref{app:proof_lemma_IRGPUCB}.
    %
    From the tower property of the expectation, it suffices to show the following inequality:
    \begin{align*}
        \mathbb{E}_t[f(\*x^*)] \leq \mathbb{E}_t \left[\mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2}_t \sigma_{t-1}(\*x) \right].
    \end{align*}
    %
    Since this inequality only considers the expectation given ${\cal D}_{t-1}$, we can fix ${\cal D}_{t-1}$ in the proof.
    %
    Furthermore, from Lemma~\ref{lem:bound_srinivas}, we can obtain the following inequality:
    \begin{align*}

        F_t^{-1}(1 - \delta) \leq \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{\delta} \sigma_{t-1}(\*x),
    \end{align*}
    where $F_t (\cdot) \coloneqq \mathrm{Pr}_t\left( f(\*x_*) < \cdot \right)$ is a cumulative distribution function of $f(\*x_*)$, and $F_t^{-1}$ is its inverse function.
    %
    Then, substituting $U \sim {\rm Uni}(0, 1)$ into $\delta$ and taking the expectation, we obtain the following inequality:
    \begin{align*}
        \mathbb{E}_t \left[ f(\*x_*) \right] \leq \mathbb{E}_U \left[ \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{U} \sigma_{t-1}(\*x) \right],
    \end{align*}
    which can be derived in a similar way to inverse transform sampling.
    %
    Hence, $\beta^{1/2}_{U} = 2 \log (|{\cal X}| / 2) - 2 \log(U)$ results in a random variable, which follows the two-parameter exponential distribution.
    %
    Consequently, we can obtain the desired bound.
\end{proof}

Lemma~\ref{lem:bound_RGPUCB} shows a direct upper bound of the expectation of $f(\*x^*)$.
A remarkable consequence of Lemma~\ref{lem:bound_RGPUCB} is an upper bound of BCR as follows:
\begin{align*}
    {\rm BCR}_T
    &= \sum_{t=1}^T \mathbb{E}\left[ f(\*x^*) - v_t(\*x_t) \right] + \mathbb{E}\left[ v_t(\*x_t) - f(\*x_t) \right] \\
    &\leq \sum_{t=1}^T \mathbb{E}\left[ v_t(\*x_t) - f(\*x_t) \right],
\end{align*}
where $v_t(\*x) \coloneqq \mu_{t-1}(\*x) + \zeta^{1/2}_t \sigma_{t-1}(\*x)$.
This upper bound eliminates the first term, which cannot be bounded without increasing the confidence parameter in the conventional regret decomposition Eq.~\eqref{eq:regret_depomsition}.
Importantly, $\mathbb{E}[\zeta_t] = 2 + s$ is a constant since $s$ and $\lambda$ for $\zeta$ does not depend on $t$ in Lemma~\ref{lem:bound_RGPUCB}.
Thus, using Lemma~\ref{lem:bound_RGPUCB}, we can obtain the tight BCR bound for the finite input domain:
\begin{theorem}[BCR bound for finite domain]
    Let $f \sim {\mathscr G} {\cal P} (0, k)$, where $k$ is a stationary kernel and $k(\*x, \*x) = 1$, and ${\cal X}$ be finite.
    %
    Assume that $\zeta_t$ follows a two-parameter exponential distribution with $s = 2 \log (|{\cal X}| / 2)$ and $\lambda = 1/2$ for any $t \geq 1$.
    %
    Then, by running IRGP-UCB with $\zeta_t$, BCR can be bounded as follows:
    \begin{align*}
        {\rm BCR}_T \leq \sqrt{C_1 C_2 T \gamma_T},
    \end{align*}
    where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$ and $C_2 \coloneqq 2 + s$.
    \label{theo:BCR_IRGPUCB_discrete}
\end{theorem}
See Appendix~\ref{app:proof_discrete} for detailed proof.

Theorem~\ref{theo:BCR_IRGPUCB_discrete} has two important implications.
First, Theorem~\ref{theo:BCR_IRGPUCB_discrete} shows the convergence rate $O(\sqrt{T \gamma_T})$, in which multiplicative $\sqrt{\log T}$ factor is improved compared with $O(\sqrt{T \gamma_T \log T})$ achieved by GP-UCB \citep{Srinivas2010-Gaussian}, RGP-UCB, and TS \citep{Russo2014-learning}.
Thus, IRGP-UCB achieves the best convergence rate for BCR.
Second, IRGP-UCB does not need to schedule the parameters of $\zeta_t$, i.e., $s$ and $\lambda$, in contrast to GP-UCB and RGP-UCB.
Therefore, through randomization, IRGP-UCB essentially alleviates the problem that the well-known $\beta_t \propto \log(t)$ strategy results in a practically too large confidence parameter.
Further, over-exploration in the later iterations can be avoided.
The IRGP-UCB is the first GP-UCB-based method that enjoys the above technical and practical benefits.

It is worth noting that our key lemma (Lemma~\ref{lem:bound_RGPUCB}) mainly requires Lemma~\ref{lem:bound_srinivas} only, which is expected to be satisfied in a wide range of exploration problems in which a GP is used as a surrogate model.
Therefore, we conjecture that the same proof technique can be applicable to more advanced problem settings, such as multi-objective BO \citep{paria2020-flexible}, for which further analyses are important future directions of our results.

Next, we show the BCR bound for the infinite (continuous) domain:
\begin{theorem}[BCR bound for infinite domain]
    Let $f \sim {\mathscr G} {\cal P} (0, k)$, where $k$ is a stationary kernel, $k(\*x, \*x) = 1$, and Assumption~\ref{assump:continuous_X} holds.
    %
    Assume that $\zeta_t$ follows a two-parameter exponential distribution with $s_t = 2d \log(bdr t^2 \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr)) - 2 \log 2$ and $\lambda = 1/2$ for any $t \geq 1$.
    %
    Then, by running IRGP-UCB, BCR can be bounded as follows:
    \begin{align*}
        {\rm BCR}_T \leq \frac{\pi^2}{6} + \sqrt{C_1 T \gamma_T (2 + s_T)},
    \end{align*}
    where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
    \label{theo:BCR_IRGPUCB_continuous}
\end{theorem}
See Appendix~\ref{app:proof_continuous} for the proof.
Unfortunately, in this case, $\mathbb{E}[\zeta_t] = O(\log t)$ and the resulting BCR bound is $O(\sqrt{T \gamma_T \log T})$.
This is because the discretization of input domain ${\cal X}_t$, which is refined as $|{\cal X}_t| = O(t^{2d})$, is required to obtain the BCR bound for the infinite domain.
Since we cannot avoid the dependence on $|{\cal X}|$ also in Theorem~\ref{theo:BCR_IRGPUCB_discrete}, the resulting BCR bound in Theorem~\ref{theo:BCR_IRGPUCB_continuous} requires the term $\log(|{\cal X}_t|) = O\bigl(2d \log(t) \bigr)$.

On the other hand, to bound the BSR, refining the discretization is not necessarily needed.
That is, we can obtain the following corollaries for finite and infinite input domains, respectively:
\begin{corollary}[BSR bound for finite domain]
    Fix a required accuracy $\eta > 0$.
    %
    Assume the same condition as in Theorem~\ref{theo:BCR_IRGPUCB_discrete}.
    %
    Then, by running IRGP-UCB, ${\rm BSR}_T \leq \eta$ by at most $T$ function evaluations, where $T$ is the smallest positive integer satisfying the following inequality:
    \begin{align*}
        \sqrt{\frac{C_1 C_2 \gamma_T}{T}} \leq \eta,
    \end{align*}
    where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$ and $C_2 \coloneqq 2 + s$.
    \label{coro:BSR_IRGPUCB_discrete}
\end{corollary}
\begin{corollary}[BSR bound for infinite domain]
    Fix a required accuracy $\eta > 0$.
    %
    Assume the same condition as in Theorem~\ref{theo:BCR_IRGPUCB_continuous} except for $s_\eta = 2d \log(2bdr \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr) / \eta) - 2 \log 2$.
    %
    Then, by running IRGP-UCB, ${\rm BSR}_T \leq \eta$ by at most $T$ function evaluations, where $T$ is the smallest positive integer satisfying the following inequality:
    \begin{align*}
        \sqrt{\frac{C_1 (2 + s_\eta) \gamma_T}{T}} \leq \frac{\eta}{2},
    \end{align*}
    where  $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
    \label{coro:BSR_IRGPUCB_continuous}
\end{corollary}
See Appendixes~\ref{app:proof_discrete} and \ref{app:proof_continuous} for the proof, respectively.
These corollaries suggest that, for both finite and infinite domains, an arbitrary accuracy solution can be obtained using finite function evaluations, whose rate is shown by the inequalities.
In Corollary~\ref{coro:BSR_IRGPUCB_continuous}, $s_\eta$ depends on $\eta$ instead of $t$.
Thus, in both corollaries, the convergence rate of the left-hand side with respect to $T$ is $O(\sqrt{\gamma_T / T})$, which is improved compared with the convergence rate $O(\sqrt{\gamma_T \log T / T})$ achieved by GP-UCB, RGP-UCB, and TS.
This faster convergence rate is an essential improvement derived based on Lemma~\ref{lem:bound_RGPUCB}.
Therefore, IRGP-UCB does not need to schedule $\zeta_t$ if we require optimization rather than cumulative regret minimization.

\section{Experiments}
\label{sec:experiment}

We demonstrate the experimental results on synthetic and benchmark functions and the materials dataset provided in \citep{liang2021benchmarking}.
As a baseline, we performed EI \citep{Mockus1978-Application}, TS \citep{Russo2014-learning}, MES \citep{Wang2017-Max}, joint entropy search (JES) \citep{hvarfner2022-joint}, and GP-UCB \citep{Srinivas2010-Gaussian}.
For the posterior sampling in TS, MES, and JES, we used random Fourier features \citep{Rahimi2008-Random}.
For Monte Carlo estimation in MES and JES, we used ten samples.
We evaluate the practical performance of BO by simple regret $f(\*x^*) - \mathop{\rm argmax}\limits_{t \leq T} f(\*x_t)$.

\subsection{Synthetic Function Experiment}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\linewidth]{fig/confidence_params_SynFun_pool.pdf}
    \caption{
        Confidence parameter of GP-UCB-based methods.
        %
        For GP-UCB, the solid line represents $\beta_t$.
        %
        For RGP-UCB and IRGP-UCB, the solid lines and shaded area represent the expectations $\mathbb{E}[\zeta_t]$ and $95 \%$ credible intervals, respectively.
        }
    \label{fig:synthetic-confidence-param}
\end{figure}
\hfill
\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\linewidth]{fig/Results_SynFun_pool_ell=0.1-d=3_log.pdf}
    \caption{Average and standard errors of simple regret.}
    \label{fig:synthetic-regret}
\end{figure}

We use synthetic functions generated from ${\mathscr G} {\cal P} (0, k)$, where $k$ is Gaussian kernel with a length scale parameter $\ell = 0.1$ and the input dimension $d = 3$.
We set the noise variance $\sigma^2 = 10^{-4}$.
The input domain consists of equally divided points in $[0, 0.9]$, i.e., ${\cal X} = \{0, 0.1, \dots, 0.9\}^d$ and $|{\cal X}| = 1000$.
Figure~\ref{fig:synthetic-confidence-param} shows the theoretical confidence parameters of GP-UCB, RGP-UCB with ${\rm Gamma}(\kappa_t=\log(|{\cal X}|t^2) / \log(1.5), \theta = 1)$, and IRGP-UCB with a two-parameter exponential distribution.
We can see that the confidence parameters for GP-UCB and RGP-UCB are considerably large due to logarithmic increase, particularly in the later iterations.
In contrast, the confidence parameter of IRGP-UCB is drastically smaller than that of GP-UCB and not changed with respect to the iteration.
Therefore, we can observe that IRGP-UCB alleviates the over-exploration.

Next, we report the simple regret in Figure~\ref{fig:synthetic-regret}.
Ten synthetic functions and ten initial training datasets are randomly generated.
Thus, the average and standard error for $10 \times 10$ trials are reported.
The hyperparameters for GP are fixed to those used to generate the synthetic functions.
The confidence parameters for GP-UCB, RGP-UCB, and IRGP-UCB are set as in Figure~\ref{fig:synthetic-confidence-param}.
We can see that the regrets of GP-UCB, RGP-UCB, TS, and JES converge slowly.
We conjecture that this slow convergence came from over-exploration.
On the other hand, EI, MES, and IRGP-UCB show fast convergence.
In particular, IRGP-UCB achieves the best average in most iterations.
These results suggest that IRGP-UCB bridges the gap between theory and practice in contrast to GP-UCB and RGP-UCB.

\subsection{Benchmark Function Experiments}

\begin{figure*}[!t]
    \centering
    \subfigure[Benchmark functions.]{
        \includegraphics[width=0.325\linewidth]{fig/Results_SimMax_Holder_table_log.pdf}
        \includegraphics[width=0.325\linewidth]{fig/Results_SimMax_Cross_in_tray_log.pdf}
        \includegraphics[width=0.325\linewidth]{fig/Results_SimMax_Ackley_log.pdf}
        }
    \label{fig:benchmark-regret}
    \hfill
    \subfigure[Materials datasets.]{
        \includegraphics[width=0.325\linewidth]{fig/Results_SimMax_Perovskite_log.pdf}
        \includegraphics[width=0.325\linewidth]{fig/Results_SimMax_P3HT_log.pdf}
        \includegraphics[width=0.325\linewidth]{fig/Results_SimMax_AgNP_log.pdf}
        }
    \label{fig:real-regret}
    \vspace{-8pt}
    \caption{Average and standard errors of simple regret.}
    \label{fig:benchmark-real-regret}
\end{figure*}

We employ three benchmark functions called Holder table ($d=2$), Cross in tray($d=2$), and Ackley ($d=4$) functions, whose analytical forms are shown at \url{https://www.sfu.ca/~ssurjano/optimization.html}.
For each function, we report the average and standard error for $10$ trails using ten random initial datasets ${\cal D}_0$, where $|{\cal D}_0| = 2^d$.
We set the noise variance $\sigma^2 = 10^{-4}$.
We used the Gaussian kernel with automatic relevance determination, whose hyperparameter was selected by marginal likelihood maximization per $5$ iterations \citep{Rasmussen2005-Gaussian}.

In this experiment, since the input domain is continuous, the theoretical choice of the confidence parameter contains an unknown variable.
Thus, we use the heuristic choice for confidence parameters.
For GP-UCB, we set the confidence parameter as $\beta_t = 0.2 d \log(2 t)$, which is the heuristics used in \citep{kandasamy2015-high,Kandasamy2017-Multi}.
For RGP-UCB, we set $\zeta_t \sim {\rm Gamma} (\kappa_t, \theta=1)$ with $\kappa_t = 0.2 d \log(2 t)$ since $\mathbb{E}[\zeta_t]$ must have the same order as $\beta_t$ (note that $\mathbb{E}[\zeta_t] = \theta \kappa_t$).
For IRGP-UCB, we set $s = d / 2$ and $\lambda = 1 / 2$.
Note that $\lambda$ is equal to the theoretical setting, and $s$ captures the dependence on $d$, as shown in Corollary~\ref{coro:BSR_IRGPUCB_continuous}.

Figure~\ref{fig:benchmark-regret} shows the results.
In the Holder table function, JES shows faster convergence until 40 iterations.
However, the regrets of all baseline methods stagnate in $[10^{-1}, 10^0]$.
In contrast, the regret of IRGP-UCB converges to $10^{-3}$ until 60 iterations.
In the Cross in tray function, IRGP-UCB showed rapid convergence at approximately 50 iterations.
In the Ackley function, IRGP-UCB constantly belonged to the top group and showed minimum regret after 125 iterations.
Since the input dimension is relatively large, TS and JES, which depends on the sampled maxima $\*x_*$, deteriorate by over-exploration.
Throughout the experiments, IRGP-UCB outperforms TS, GP-UCB, and RGP-UCB, which achieves sub-linear BCR bounds, and EI and MES, which are practically well-used.
This result supports the effectiveness of randomization using the two-parameter exponential distribution.

\subsection{Real-World Dataset Experiments}

This section provides the experimental results on the materials datasets provided in \citep{liang2021benchmarking}.
In the perovskite dataset \citep{sun2021data}, we optimize environmental stability with respect to composition parameters for halide perovskite ($d = 3$ and $|{\cal X}| = 94$).
In the P3HT/CNT dataset \citep{bash2021multi}, we optimize electrical conductivity with respect to composition parameters for carbon nanotube polymer blend ($d = 5$ and $|{\cal X}| = 178$).
In the AgNP dataset \citep{mekki2021two}, we optimize the absorbance spectrum of synthesized silver nanoparticles with respect to processing parameters for synthesizing triangular nanoprisms ($d = 5$ and $|{\cal X}| = 164$).
See \citep{liang2021benchmarking} for more details about each dataset.

We set the initial dataset size $|{\cal D}_0| = 2$ as with \citep{liang2021benchmarking}.
Since the dataset size is small at earlier iterations and the dataset contains fluctuations from real-world experiments, we observed that hyperparameter tuning could be unstable.
Thus, we optimized the hyperparameters of the RBF kernel in each iteration to avoid repeatedly obtaining samples using an inappropriate hyperparameter.
The other settings matched those used in the benchmark function experiments.

In the perovskite dataset, IRGP-UCB constantly belonged to the top group and showed the best performance after 20 iterations.
In the P3HT/CNT dataset, EI converged to the smallest value after 60 iterations.
On the other hand, IRGP-UCB shows faster convergence during the first 20 iterations.
In the AgNP dataset, IRGP-UCB found the optimal point until 42 iterations in all the trials earliest.
In this experiment, heuristic methods, EI, MES, and JES, showed worse performance and required at least 60 function evaluations to find the optimal point.
Consequently, we can observe the effectiveness of IRGP-UCB against real-world datasets.

\section{Conclusion}
\label{sec:conculusion}

First, by showing the generalized BCR bounds for RGP-UCB, this study showed that randomization without ingenuity does not improve the regret bounds.
Then, this study proposed an improved randomized GP-UCB (IRGP-UCB) using the two-parameter exponential distribution, whose Bayesian regret bound is tighter than known results.
Lemma~\ref{lem:bound_RGPUCB}, which directly bounds the expectation of the maximum value, plays a key role in the proof.
Additionally, we demonstrated the practical effectiveness of IRGP-UCB through extensive experiments, including the application to the materials dataset.

Several directions for future works can be considered.
First, whether we can show the tight BCR bounds for the continuous domain is of interest.
Second, since the Bayesian regret analysis of TS depends on the usual UCB, we may improve the BCR bound of TS using randomized UCB.
Last, we may be able to extend IRGP-UCB to the other various practical settings, where the usual UCB-based methods have been extended, as listed in Section~\ref{sec:intro}.

\section{Proof for BCR Bound of GP-UCB}
\label{app:GP_UCB_proof}

\citet{Srinivas2010-Gaussian} provided the high-probability bounds for cumulative regret, not BCR bounds.
On the other hand, although \citet{paria2020-flexible} have shown BCR bounds for the GP-UCB-based multi-objective BO method, they did not provide explicit proof for standard single-objective BO.
Thus, although the following result is an almost direct consequence of Theorem~2 in \citep{paria2020-flexible}, we here provide the explicit proof for standard GP-UCB for completeness.
On the other hand, for the case of (ii), the bound with respect to $a$ is tightened as $O\bigl( \log ( \log (a)) \bigr)$ compared with $O( \log(a))$ \citep{Kandasamy2018-Parallelised, paria2020-flexible}, which is mainly based on Lemma~\ref{lem:L_max_bound}.

\begin{theorem}
    Suppose that $f \sim {\mathscr G} {\cal P}(0, k)$, where $k$ is a stationary kernel and $k(\*x, \*x) = 1$.

    \vspace{-0.8\baselineskip}
    \begin{itemize}
        \setlength{\parskip}{0.05cm}
        \setlength{\itemsep}{0.05cm}
        \item[(i)] When ${\cal X}$ is a finite set, GP-UCB with $\beta_t = 2 \log(|{\cal X}| t^2 / \sqrt{2\pi} )$ achieves the following BCR bound:
        \begin{align}
            {\rm BCR}_T \leq \frac{\pi^2}{6} + \sqrt{ C_1 T \beta_T \gamma_T},
        \end{align}
        where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
        \item[(ii)] When ${\cal X}$ is an infinite set, let Assumption~\ref{assump:continuous_X} holds.
                    %
                    Then, GP-UCB with $\beta_t = 2d \log \bigl(bdr t^2 \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr) \bigr) + 2\log (t^2 / \sqrt{2\pi})$ achieves the following BCR bound:
        \begin{align}
            {\rm BCR}_T \leq \frac{\pi^2}{3} + \sqrt{ C_1 T \beta_T \gamma_T},
        \end{align}
        where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
    \end{itemize}
    \vspace{-0.8\baselineskip}
    \label{theo:BCR_GPUCB}
\end{theorem}
\begin{proof}
    First, we show the proof of (i).
    %
    For regret analysis of BCR, regret decomposition  \citep{Russo2014-learning} with $U_t(\*x) \coloneqq \mu_{t-1}(\*x) + \beta_t^{1/2} \sigma_{t-1}(\*x)$ is used as follows:
    \begin{align}
        {\rm BCR}_T
        &= \sum_{t=1}^T \mathbb{E} [f(\*x^*) - f(\*x_t)] \\
        &= \sum_{t=1}^T \mathbb{E} [f(\*x^*) - U_t(\*x^*) + U_t(\*x^*) - U_t(\*x_t) + U_t(\*x_t) - f(\*x_t)] \\
        &\leq \underbrace{\sum_{t=1}^T \mathbb{E} [f(\*x^*) - U_t(\*x^*)]}_{R_1} + \underbrace{\sum_{t=1}^T \mathbb{E} [ U_t(\*x_t) - f(\*x_t)]}_{R_2}. && \bigl( \because U_t(\*x^*) - U_t(\*x_t) \leq 0 \bigr)
    \end{align}

    First, we consider $R_1$.
    %
    We use the following fact for a Gaussian random variable $Z \sim {\cal N}(m, s^2)$ with $m \leq 0$:
    \begin{align}
        \mathbb{E}[Z_+] \leq \frac{s}{\sqrt{2 \pi}} \exp \left\{ - \frac{m^2}{2 s^2} \right\}, \label{eq:Gauss_plus_expectation_1}
    \end{align}
    where $Z_+ \coloneqq \mathop{\rm argmax}\limits \{0, Z\}$.
    %
    This fact has been frequently used in this literature \citep{Russo2014-learning,Kandasamy2018-Parallelised,paria2020-flexible}.
    %
    Then, we can extend $R_1$ as follows:
    %
    \begin{align}
        R_1
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}} [ \mathbb{E}_t[ f(\*x^*) - U_t(\*x^*)]] \\
        &\leq \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}} [ \mathbb{E}_t[ \bigl( f(\*x^*) - U_t(\*x^*) \bigr)_+ ]] && \bigl( \because f(\*x^*) - U_t(\*x^*) \leq \bigl( f(\*x^*) - U_t(\*x^*) \bigr)_+ \bigr) \\

        &\leq \sum_{t=1}^T \sum_{\*x \in {\cal X}} \mathbb{E}_{{\cal D}_{t-1}} \left[ \mathbb{E}_t \left[ \bigl(f(\*x) - U_t(\*x)\bigr)_+  \right] \right]. && \left( \because \bigl( f(\*x^*) - U_t(\*x^*) \bigr)_+ \leq \sum_{\*x \in {\cal X}} \bigl( f(\*x) - U_t(\*x) \bigr)_+ \right)
    \end{align}
    %
    Then, since $f(\*x) - U_t(\*x) \mid {\cal D}_{t-1} \sim {\cal N}( -\beta_t^{1/2}(\*x)\sigma_{t-1}(\*x), \sigma^2_{t-1}(\*x))$, we can apply Eq.~\eqref{eq:Gauss_plus_expectation_1}:
    \begin{align}
        R_1
        &\leq \sum_{t=1}^T \sum_{\*x \in {\cal X}} \mathbb{E}_{{\cal D}_{t-1}} \left[ \frac{\sigma_{t-1}(\*x)}{\sqrt{2 \pi}} \exp \left\{ - \frac{\beta_t}{2} \right\} \right] && \bigl( \because \text{Eq.~\eqref{eq:Gauss_plus_expectation_1}} \bigr) \\
        &\leq \sum_{t=1}^T |{\cal X}| \frac{1}{\sqrt{2 \pi}} \exp \left\{ - \frac{\beta_t}{2} \right\} && \bigl( \because \sigma_{t-1}(\*x) \leq \sigma_0(\*x) = 1 \bigr) \\
        &= \sum_{t=1}^T \frac{1}{t^2} \\
        &\leq \frac{\pi^2}{6}. && \left(\because \sum_{t=1}^\infty \frac{1}{t^2} = \frac{\pi^2}{6} \right) \label{eq:GP_UCB_R1_bound}
    \end{align}

    Then, we bound $R_2$ as follows:
    \begin{align}
        R_2
        &= \sum_{t=1}^T \mathbb{E} [ U_t(\*x_t) - f(\*x_t)] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}} [ \mathbb{E}_t [ U_t(\*x_t) - f(\*x_t)] ] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}} [ U_t(\*x_t) - \mu_{t-1}(\*x_t) ] \\
        &= \sum_{t=1}^T \mathbb{E} \left[ \beta^{1/2}_t \sigma_{t-1}(\*x_t) \right] \\
        &= \mathbb{E} \left[ \sum_{t=1}^T \beta^{1/2}_t \sigma_{t-1}(\*x_t) \right] \\
        &\leq \mathbb{E} \left[ \sqrt{ \sum_{t=1}^T \beta_t  \sum_{t=1}^T \sigma^2_{t-1}(\*x_t) } \right] && \bigl( \because \text{Cauchy--Schwarz  inequality} \bigr) \\
        &= \sqrt{ \sum_{t=1}^T \beta_t} \mathbb{E} \left[ \sqrt{ \sum_{t=1}^T \sigma^2_{t-1}(\*x_t) } \right] \\
        &\leq \sqrt{ T \beta_T} \mathbb{E} \left[ \sqrt{ \sum_{t=1}^T \sigma^2_{t-1}(\*x_t) } \right] && \bigl( \because \beta_1 < \beta_2 < \dots < \beta_T \bigr).
    \end{align}
    %
    Then, from Lemma~5.4 in \citep{Srinivas2010-Gaussian}, $\sqrt{\sum_{t=1}^T \sigma^2_{t-1}(\*x_t)} \leq \sqrt{C_1 \gamma_T}$ with probability $1$.
    %
    Thus, we obtain
    \begin{align}
        R_2 \leq \sqrt{ C_1 T \beta_T \gamma_T}. \label{eq:GP_UCB_R2_bound}
    \end{align}
    Consequently, by combining Eqs.~\eqref{eq:GP_UCB_R1_bound} and \eqref{eq:GP_UCB_R2_bound}, we obtain the theorem for (i).

    Next, we consider (ii).
    %
    Purely for the sake of analysis, we use a set of discretization ${\cal X}_t \subset {\cal X}$ for $t \geq 1$.
    %
    For any $t \geq 1$, let ${\cal X}_t \subset {\cal X}$ be a finite set with each dimension equally divided into $\tau_t = bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)$.
    %
    Thus, $|{\cal X}_t| = \tau_t^d$.
    %
    In addition, we define $[\*x]_t$ as the nearest points in ${\cal X}_t$ of $\*x \in {\cal X}$.

    Then, BCR can be decomposed as follows:
    \begin{align*}
        {\rm BCR}_T
        &= \sum_{t=1}^T \mathbb{E} \left[f(\*x^*) - f([\*x^*]_t) + f([\*x^*]_t) - U_t([\*x^*]_t) + U_t([\*x^*]_t) - U_t(\*x_t) + U_t(\*x_t) - f(\*x_t) \right] \\
        &\leq \sum_{t=1}^T \mathbb{E} \left[f(\*x^*) - f([\*x^*]_t) + f([\*x^*]_t) - U_t([\*x^*]_t) + U_t(\*x_t) - f(\*x_t) \right],
    \end{align*}
    where $\mathbb{E} \bigl[ U_t([\*x^*]_t) - U_t(\*x_t) \bigr] \leq 0$ since $\*x_t = \mathop{\text{argmax}}_{\*x \in {\cal X}} U_t (\*x)$.
    %
    Furthermore, we obtain the following:
    \begin{align*}
        \sum_{t=1}^T \mathbb{E} \left[f(\*x^*) - f([\*x^*]_t) \right]
        &\leq \sum_{t=1}^T \mathbb{E} \left[ \sup_{\*x \in {\cal X}} f(\*x) - f([\*x]_t) \right] \\
        &\leq \sum_{t=1}^T \frac{1}{t^2} && \left( \because \text{Lemma}~\ref{lem:discretized_error} \right) \\
        &\leq \frac{\pi^2}{6} && \left(\because \sum_{t=1}^\infty \frac{1}{t^2} = \frac{\pi^2}{6} \right)
    \end{align*}
    %
    The remaining terms $\sum_{t=1}^T \mathbb{E}\bigl[ f([\*x^*]_t) - U_t([\*x^*]_t) + U_t(\*x_t) - f(\*x_t) \bigr]$ can be bounded as with the case (i) by setting $\beta_t = 2 \log \bigl( |{\cal X}_t| t^2 / \sqrt{2\pi} \bigr)$.
    %
    That is,
    \begin{align}
        \sum_{t=1}^T \mathbb{E}\bigl[ f([\*x^*]_t) - U_t([\*x^*]_t) + U_t(\*x_t) - f(\*x_t) \bigr]
        &\leq \frac{\pi^2}{6} + \sqrt{C_1 T \beta_T \gamma_T}.
    \end{align}
    %
    By substituting $|{\cal X}_t| = \bigl(  bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)\bigr)^d$, we conclude the proof.
\end{proof}

\section{Technical Issues of Proof in \citep{berk2021-randomized}}
\label{app:issues_RGPUCB}

The proofs of Theorems~1, 2, and 3 in \citep{berk2021-randomized} contain many technical issues.
We enumerated them using the notation in \citep{berk2021-randomized} as follows:
\begin{itemize}
    \item
        Lemma~1 does not necessarily hold since $\alpha([\*x_t]_\tau)$ can be smaller than $\alpha([\*x^*]_\tau)$ although $\alpha(\*x_t) \geq \alpha(\*x^*)$.
    %
    \item
        To treat $f(\*x_t)$ as the random variable, which follows ${\cal N}(\mu_t(\*x), \sigma^2_t(\*x))$, the conditioning by $D_t$ is required.
        %
        However, they do not perform the conditioning in most equations, such as Eqs.~(11) and (12).
        %
        Hence, most of their treatment for conditional expectations is flawed.
    %
    \item
        From Eq.~(13) to the next equations, the term $|{\cal X}_{\rm dis}|$ is eliminated.
        %
        This term is imperatively required.
    %
    \item
        In Eqs.~(11), (12), (16), and (17), variables such as $\sigma_{t-1}([\*x_t]_r)$ are random variables that depend on $D_t$ and $[\*x_t]_r$.
        %
        However, these random variables are taken out from the expectation.
    %
    \item
        From Eq.(17) to Eq.~(19), a discretization for $\*x_t$ is eliminated.
    %
    \item
        In Eq.~(18), an approximately equal sign $\approx$ is used, which is inappropriate for theoretical proof.
    %
    \item
        The convergence rate of the term $F^{-1}(1 - 1/T)$ in the resulting bound is not clarified.
\end{itemize}
Hence, Theorems~1, 2, and 3 in \citep{berk2021-randomized} do not hold.
Furthermore, even if Theorem~3 holds, it would be insufficient to show the sub-linear BCR bound.

\section{Proof of Theorem for RGP-UCB}
\label{app:proof_general_RGPUCB}

Here, we rectify and generalize the theorems in \citep{berk2021-randomized}.
First, we show the generalized theorem, which shows the upper bound using the expectation and MGFs of $\{ \zeta_t \}_{t \geq 1}$.
Next, we provide some examples of the distributions, which can achieve the sub-linear regret bound by appropriately controlling the parameters.

\subsection{Generalized BCR bounds for RGP-UCB}

The following theorem shows the upper bound of BCR for setting $\{\zeta_t\}_{t \geq 1}$ as arbitrary non-negative random variables:
\begin{theorem}
    Let $f \sim {\mathscr G} {\cal P}(0, k)$, where $k$ is a stationary kernel $k$ and $k(\*x, \*x) = 1$.
    %
    Let $\{\zeta_t\}_{t \geq 1}$ be a sequence of non-negative random variables and their MGFs be $M_t (\cdot)$, where $M_t (-1/2) < \infty$ at least.
    \vspace{-0.8\baselineskip}
    \begin{itemize}
        \setlength{\parskip}{0.05cm}
        \setlength{\itemsep}{0.05cm}
        \item[(i)] When ${\cal X}$ is a finite set, by running RGP-UCB with $\zeta_t$, BCR can be bounded as follows:
            \begin{align*}
                {\rm BCR}_T \leq \frac{ |{\cal X}|}{\sqrt{2 \pi}} \sum_{t=1}^T M_t(-1/2) + \sqrt{ \sum_{t=1}^T \mathbb{E}[\zeta_t]  C_1 \gamma_T},
            \end{align*}
            where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
        \item[(ii)] When ${\cal X}$ is an infinite set, let Assumption~\ref{assump:continuous_X} holds.
                    %
                    By running RGP-UCB with $\zeta_t$, BCR can be bounded as follows:
            \begin{align}
                {\rm BCR}_T \leq \frac{\pi^2}{6} +  \frac{1}{\sqrt{2 \pi}} \sum_{t=1}^T (\tau_t)^d M_t(-1/2)  + \sqrt{ C_1 \sum_{t=1}^T \mathbb{E}[\zeta_t] \gamma_T},
            \end{align}
            where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$ and $\tau_t = bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)$.
    \end{itemize}
    \vspace{-0.8\baselineskip}
    \label{theo:BCR_RGPUCB}
\end{theorem}
\begin{proof}
    The procedure of the proof is similar to the proof for the standard GP-UCB shown in Appendix~\ref{app:GP_UCB_proof}.
    %
    BCR can be decomposed by $V_t(\*x) = \mu_{t-1}(\*x) + \zeta_t^{1/2} \sigma_{t-1}(\*x)$ as follows:
    \begin{align*}
        {\rm BCR}_T
        &= \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - V_t(\*x^*) + V_t(\*x^*) - V_t(\*x_t) + V_t(\*x_t) - f(\*x_t) \right] \\
        &\leq \underbrace{\sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - V_t(\*x^*) \right]}_{\eqqcolon R_1} + \underbrace{\sum_{t=1}^T \mathbb{E}\left[ V_t(\*x_t) - f(\*x_t) \right]}_{ \eqqcolon R_2}. && \bigl(\because V_t(\*x^*) - V_t(\*x_t) \leq 0 \bigr)
    \end{align*}
    %
    Then, we derive the upper bounds of $R_1$ and $R_2$, respectively.

    First, we consider $R_1$.
    %
    We use the following fact for Gaussian random variables $Z \sim {\cal N}(m, s^2)$ with $m \leq 0$:
    \begin{align}
        \mathbb{E}[Z_+] \leq \frac{s}{\sqrt{2 \pi}} \exp \left\{ - \frac{m^2}{2 s^2} \right\}, \label{eq:Gauss_plus_expectation}
    \end{align}
    where $Z_+ \coloneqq \mathop{\rm argmax}\limits \{0, Z\}$.
    %
    This fact is frequently used in the literature \citep{Russo2014-learning,Kandasamy2018-Parallelised,paria2020-flexible}.
    %
    Then, we can extend $R_1$ as follows:
    \begin{align}
        R_1
        &= \sum_{t=1}^T \mathbb{E}\left[ f(\*x^*) - V_t(\*x^*) \right] \\
        &\leq \sum_{t=1}^T \mathbb{E}\left[ \bigl( f(\*x^*) - V_t(\*x^*) \bigr)_+ \right] && \bigl(\because f(\*x^*) - V_t(\*x^*) \leq \bigl( f(\*x^*) - V_t(\*x^*) \bigr)_+ \bigr) \\
        &\leq \sum_{t=1}^T \mathbb{E}\left[ \sum_{\*x \in {\cal X}} \bigl( f(\*x) - V_t(\*x) \bigr)_+ \right] && \left(\because \bigl( f(\*x^*) - V_t(\*x^*) \bigr)_+ \leq \sum_{\*x \in {\cal X}} \bigl( f(\*x) - V_t(\*x) \bigr)_+ \right)\\

        &= \sum_{t=1}^T \sum_{\*x \in {\cal X}} \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} \left[ \mathbb{E} \left[ \bigl( f(\*x) - V_t(\*x) \bigr)_+ \mid {\cal D}_{t-1}, \zeta_t \right] \right].
    \end{align}
    %
    Here, $\bigl( f(\*x) - V_t(\*x) \bigr) \mid {\cal D}_{t-1}, \zeta_t$ follows ${\cal N}(-\zeta_t^{1/2} \sigma_{t-1}(\*x), \sigma_{t-1}^2(\*x))$.
    %
    Thus, by using Eq.~\eqref{eq:Gauss_plus_expectation}, we can obtain the following:
    \begin{align}
        R_1
        &\leq \sum_{t=1}^T \sum_{\*x \in {\cal X}} \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} \left[ \frac{\sigma_{t-1}(\*x)}{\sqrt{2 \pi}} \exp \left\{ - \frac{\zeta_t}{2} \right\} \right] && \bigl( \because \text{Eq.~\eqref{eq:Gauss_plus_expectation}} \bigr)\\

        &\leq \frac{|{\cal X}|}{\sqrt{2 \pi}} \sum_{t=1}^T \mathbb{E}_{\zeta_t} \left[ \exp \left\{ - \frac{\zeta_t}{2} \right\} \right] && \bigl( \because \sigma_{t-1}(\*x) \leq \sigma_0 (\*x) = 1 \bigr) \\
        &= \frac{|{\cal X}|}{\sqrt{2 \pi}} \sum_{t=1}^T M_t(- 1/2). \label{eq:R_1_bound} && \bigl(\because \text{Definition of MGF} \bigr)
    \end{align}

    Second, we expand $R_2$ as follows:
    %
    \begin{align}
        R_2
        &= \sum_{t=1}^T \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \right] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} \left[ \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \mid {\cal D}_{t-1}, \zeta_t \right] \right] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} \left[ V_t(\*x_t) - \mu_{t-1}(\*x_t) \right] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} \left[ \zeta_t^{1/2} \sigma_{t-1}(\*x_t) \right] \\
        &= \mathbb{E} \left[ \sum_{t=1}^T \zeta_t^{1/2} \sigma_{t-1}(\*x_t) \right] \\
        &\leq \mathbb{E} \left[ \sqrt{ \sum_{t=1}^T \zeta_t \sum_{t=1}^T \sigma_{t-1}^2(\*x_t) } \right]. && \bigl( \because \text{Cauchy--Schwarz  inequality} \bigr)
    \end{align}
    %
    Then, as with Lemma~5.4 in \citep{Srinivas2010-Gaussian}, $\sqrt{\sum_{t=1}^T \sigma_{t-1}^2(\*x_t) } \leq \sqrt{C_1 \gamma_T}$ with probability 1, where $C_1 = 2 / \log(1 + \sigma^{-2})$.
    %
    Hence,
    \begin{align}
        R_2
        &\leq \mathbb{E} \left[ \sqrt{ \sum_{t=1}^T \zeta_t } \right] \sqrt{C_1 \gamma_T} \\
        &\leq \sqrt{ \mathbb{E} \left[ \sum_{t=1}^T \zeta_t \right] } \sqrt{C_1 \gamma_T} && \bigl( \because \text{Jensen's  inequality} \bigr) \\
        &= \sqrt{ \sum_{t=1}^T \mathbb{E} \left[ \zeta_t \right] } \sqrt{C_1 \gamma_T}. \label{eq:R_2_bound}
    \end{align}
    %
    Consequently, combining Eqs.~\eqref{eq:R_1_bound} and \eqref{eq:R_2_bound}, we can obtain the theorem for (i).

    Then, we consider the case (ii).
    %
    As with the proof of Theorem~\ref{theo:BCR_GPUCB}, purely for the sake of analysis, we used a set of discretization ${\cal X}_t \subset {\cal X}$ for $t \geq 1$.
    %
    For any $t \geq 1$, let ${\cal X}_t \subset {\cal X}$ be a finite set with each dimension equally divided into $\tau_t = bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)$.
    %
    Thus, $|{\cal X}_t| = \tau_t^d$.
    %
    In addition, we define $[\*x]_t$ as the nearest points in ${\cal X}_t$ of $\*x \in {\cal X}$.

    Then, we decompose BCR as follows:
    \begin{align}
        {\rm BCR}_T
        &= \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f([\*x^*]_t) + f([\*x^*]_t) - V_t([\*x^*]_t) + V_t([\*x^*]_t) - V_t(\*x_t) + V_t(\*x_t) - f(\*x_t) \right] \\
        &\leq \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f([\*x^*]_t) + f([\*x^*]_t) - V_t([\*x^*]_t) + V_t(\*x_t) - f(\*x_t) \right],
    \end{align}
    where $\mathbb{E} \bigl[ V_t([\*x^*]_t) - V_t(\*x_t) \bigr] \leq 0$ since $\*x_t = \mathop{\text{argmax}}_{\*x \in {\cal X}} V_t(\*x)$.
    %
    Furthermore, we obtain the following:
    \begin{align*}
        \sum_{t=1}^T \mathbb{E} \left[f(\*x^*) - f([\*x^*]_t) \right]
        &\leq \sum_{t=1}^T \mathbb{E} \left[ \sup_{\*x \in {\cal X}} f(\*x) - f([\*x]_t) \right] \\
        &\leq \sum_{t=1}^T \frac{1}{t^2} && \left( \because \text{Lemma}~\ref{lem:discretized_error} \right) \\
        &\leq \frac{\pi^2}{6} && \left(\because \sum_{t=1}^\infty \frac{1}{t^2} = \frac{\pi^2}{6} \right)
    \end{align*}
    %
    The remaining terms $\sum_{t=1}^T \mathbb{E}\bigl[ f([\*x^*]_t) - V_t([\*x^*]_t) + V_t(\*x_t) - f(\*x_t) \bigr]$ can be bounded as with the case (i).
    %
    That is,
    \begin{align}
        \sum_{t=1}^T \mathbb{E}\bigl[ f([\*x^*]_t) - V_t([\*x^*]_t) + V_t(\*x_t) - f(\*x_t) \bigr]
        &\leq \frac{1}{\sqrt{2\pi}} \sum_{t=1}^T |{\cal X}_t| M_t(-1/2) + \sqrt{C_1 \sum_{t=1}^T \mathbb{E}[\zeta_t] \gamma_T}.
    \end{align}
    %
    By substituting $|{\cal X}_t| = (\tau_t)^d$, we conclude the proof.
\end{proof}

\subsection{Example of Distributions for Sub-linear BCR Bounds}
\label{app:examples_RGPUCB_dist}
Next, we show examples of the distributions that can achieve sub-linear regret bound using Theorem~\ref{theo:BCR_RGPUCB}.
Here, we discuss the case (ii) for infinite ${\cal X}$.
The derivation for case (i) can be obtained easily by replacing $|{\cal X}_t| = (\tau_t)^d$ with $|{\cal X}|$.

Theorem~\ref{theo:BCR_RGPUCB} shows that BCR can be bounded by the expectation and MGF.
Roughly speaking, if we can control MGF as $M_t (- 1/2) = O(1 / (t^2 |{\cal X}_t|) )$ and the expectation $\mathbb{E}[\zeta_t] = O(\log t)$ by scheduling the parameter of the distribution of $\zeta_t$, we can obtain sub-linear regret bounds, which is similar to Theorem~\ref{theo:BCR_GPUCB}.
This is because we can bound the first term by $ \frac{1}{\sqrt{2 \pi}} \sum_{t=1}^T 1 / t^2 = \pi^{2} / (6 \sqrt{2 \pi})$ if the following inequality holds:
\begin{align}
 |{\cal X}_t| M_t(- 1 / 2) \leq \frac{1}{t^2}.
\end{align}
This inequality can be rewritten as,
\begin{align}
    M_t(- 1 / 2) \leq \frac{1}{|{\cal X}_t| t^2}. \label{eq:condition_MGF}
\end{align}
Furthermore, the expectation  $\mathbb{E}[\zeta_t] = O(\log t)$ is required to bound the second term.

The above requirements can be satisfied not only by the Gamma distribution proposed in \citep{berk2021-randomized} but also by many other distributions, such as two-parameter exponential and truncated normal distributions.
Their PDFs, MGFs, and expectations are listed in Table~\ref{tab:dist_properties}.
Note that although we describe only those three distributions, we can use other distributions, such as chi-squared distribution and truncated normal with other truncation, to obtain sub-linear regret bound from Theorem~\ref{theo:BCR_RGPUCB}.
Each scheduling of the parameter and resulting regret bound are shown as follows:

\begin{table}[!t]
    \centering
    \caption{PDFs, MGFs, and expectations for distributions of $\zeta$, where $\Gamma$ is the Gamma function, and $\phi$ and $\Phi$ are PDF and cumulative distribution function of standard normal distribution, respectively.}
    \begin{tabular}{c|c|c|c}
         & PDF & MGF: $M(-1/2)$ & Expectation  \\ \hline
         Gamma: $ (k, \theta)$ & $\frac{\zeta^{k-1} e^{- \zeta / \theta}}{\Gamma(k) \theta^k}$ & $(1 + \theta / 2)^{- k}$ & $k \theta$ \\ \hline

         Two-parameter exponential: $(s, \lambda)$ & $\lambda e^{- \lambda (\zeta - s)}$ for $\zeta \geq s \geq 0$ & $\frac{\lambda}{\lambda + 1 / 2} e^{- s / 2}$& $s + 1 / \lambda$ \\ \hline
         Truncated normal: $(m, s^2 = 1)$ &  $\frac{\phi( (\zeta - m) / s)}{\Phi(1) - \Phi(-1)}$ & \multirow{2}{*}{$e^{- m / 2 + 1 / 8} \frac{\Phi(3/2) - \Phi(-1/2)}{\Phi(1) - \Phi(-1)} $} & \multirow{2}{*}{$m$} \\
         with truncation $[m - s, m + s]$ & for $m - s \leq \zeta \leq m + s$ &&
    \end{tabular}
    \label{tab:dist_properties}
\end{table}

 \paragraph{Gamma distribution:} If we use $\zeta_t$ following Gamma distribution with $k_t$ and $\theta$, we require the following:
 \begin{align*}
     (1 + \theta / 2)^{-k_t} \leq \frac{1}{|{\cal X}_t| t^2}.
 \end{align*}
 %
 Therefore, we obtain the following:
 \begin{align*}
     k_t \geq \frac{\log ( |{\cal X}_t| t^2  )}{\log(1 + \theta / 2)}.
 \end{align*}
 %
 The above derivation is mostly the same as the proof of Theorem~1 in \citep{berk2021-randomized}.
 %
 However, \citet{berk2021-randomized} bounded as $ \sqrt{\sum_{t=1}^T \mathbb{E}[\zeta_t]} \leq \sqrt{T \mathbb{E}[\mathop{\rm argmax}\limits_{t \leq T} \zeta_t]}$.
 %
 This operation is correct, but the evaluation of $\mathbb{E}[\mathop{\rm argmax}\limits_{t \leq T} \zeta_t]$ is slightly complicated.
 %
 To bound this term, \citet{berk2021-randomized} used nearly equal, which is forbidden.
 %
 On the other hand, we know the expectation $\mathbb{E}[\zeta_t] = k_t \theta$ analytically.
 %
 Let
 \begin{align}
     k_t
     &= \frac{\log ( |{\cal X}_t| t^2  )}{\log(1 + \theta / 2)} \\
     &= \frac{ d\log \bigl( bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)  \bigr) + 2\log (t)}{\log(1 + \theta / 2)}.
 \end{align}
 %
 Consequently, since $\mathbb{E}[\zeta_t] = k_t \theta$ is monotone increasing, we obtain the following bound:
 \begin{align}
     {\rm BCR}_T
     \leq \frac{\pi^2}{6} + \frac{\pi^2}{6\sqrt{2\pi}} + \sqrt{C_1 \gamma_t T \mathbb{E}[\zeta_T]}
     = \frac{\pi^2}{6} + \frac{\pi^2}{6\sqrt{2\pi}} + \sqrt{C_1 \gamma_t T \theta k_T}.
 \end{align}

\paragraph{Two-parameter exponential distribution:} If we use $\zeta_t$ following two-parameter exponential distribution with $s_t$ and $\lambda$, we requires the following:
 \begin{align*}
     \frac{\lambda}{\lambda + 1 / 2} e^{- s_t / 2} \leq e^{- s_t / 2} \leq \frac{1}{|{\cal X}_t| t^2},
 \end{align*}
 where we bound $\frac{\lambda}{\lambda + 1 / 2} < 1$ so that the condition for $\lambda$ to maintain $s_t \geq 0$ will be eliminated.
 %
 Therefore, we obtain the following:
 \begin{align*}
     s_t \geq 2 \log \left( |{\cal X}_t| t^2 \right).
 \end{align*}
 Hence, by setting $s_t = 2 \log \left( |{\cal X}_t| t^2 \right) = 2d\log \bigl( bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)  \bigr) + 4\log (t)$, since $\mathbb{E}[\zeta_t] = s_t + 1/\lambda$ is monotone increasing, we obtain the following bound:
 \begin{align}
     {\rm BCR}_T
     \leq \frac{\pi^2}{6} + \frac{\pi^2}{6\sqrt{2\pi}} + \sqrt{C_1 \gamma_t T \mathbb{E}[\zeta_T]}
     = \frac{\pi^2}{6} + \frac{\pi^2}{6\sqrt{2\pi}} + \sqrt{C_1 \gamma_t T (s_T + 1/\lambda)}.
 \end{align}

\paragraph{Truncated normal distribution:} For brevity, we use a truncated normal distribution with $s^2 = 1$ and truncation $[m - s, m + s]$.
We schedule $m_t$ so that the following inequality holds:
\begin{align*}
    e^{- m / 2 + 1 / 8} \frac{\Phi(3/2) - \Phi(-1/2)}{\Phi(1) - \Phi(-1)} < e^{- m_t / 2 + 1 / 8} < \frac{1}{|{\cal X}_t| t^2},
\end{align*}
where we bound $\frac{\Phi(1) - \Phi(0)}{\Phi(1/2) - \Phi(-1/2)} < 1$ for simplicity.
Then, we obtain the following condition:
\begin{align*}
    m_t \geq 2 \log \left( |{\cal X}_t| t^2 \right) + 1/4.
\end{align*}
Furthermore, the lower bound $m_t - s > 0$ should be satisfied.
Hence, by setting $m_t = 2 \log \left( |{\cal X}_t| t^2 \right) + 1 = 2d\log \bigl( bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)  \bigr) + 4\log (t) + 1$, since $\mathbb{E}[\zeta_t] = m_t$ is monotone increasing, we obtain the following bound:
\begin{align}
     {\rm BCR}_T
     \leq \frac{\pi^2}{6} + \frac{\pi^2}{6\sqrt{2\pi}} + \sqrt{C_1 \gamma_t T \mathbb{E}[\zeta_T]}
     = \frac{\pi^2}{6} + \frac{\pi^2}{6\sqrt{2\pi}} + \sqrt{C_1 \gamma_t T m_T}.
 \end{align}

\section{Detailed Proof of Lemma~\ref{lem:bound_RGPUCB}}
\label{app:proof_lemma_IRGPUCB}

This section provides detailed proof of the following Lemma~\ref{lem:bound_RGPUCB}:
\begin{replemma}{lem:bound_RGPUCB}
    Let $f \sim {\mathscr G} {\cal P} (0, k)$, where $k$ is a stationary kernel and $k(\*x, \*x) = 1$, and ${\cal X}$ be finite.
    %
    Assume that $\zeta$ follows a two-parameter exponential distribution with $s = 2 \log (|{\cal X}| / 2)$ and $\lambda = 1/2$.
    %
    Then, the following inequality holds:
    \begin{align*}
        \mathbb{E}[f(\*x^*)] \leq \mathbb{E} \left[\mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2}_t \sigma_{t-1}(\*x) \right],
    \end{align*}
    for all $t \geq 1$.
\end{replemma}
\begin{proof}
    It suffices to show that, for an arbitrary ${\cal D}_{t-1}$,
    \begin{align*}
        \mathbb{E}_t[f(\*x^*) | {\cal D}_{t-1}] \leq \mathbb{E}_t \left[\mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2}_t \sigma_{t-1}(\*x) \right],
    \end{align*}
    since the claim of Lemma~\ref{lem:bound_RGPUCB} can be rewritten as follows:
    \begin{align*}
        &\mathbb{E}_{{\cal D}_{t-1}}[ \mathbb{E}_t[f(\*x^*) ]] \\
        &\leq \mathbb{E}_{{\cal D}_{t-1}} \left[ \mathbb{E}_t \left[ \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2} \sigma_{t-1}(\*x) \right] \right].
    \end{align*}
    %
    Then, fix the dataset ${\cal D}_{t-1}$ and $\delta \in (0, 1)$.
    %
    From Lemma~\ref{lem:bound_srinivas}, we can see that
    \begin{align*}
        \mathrm{Pr}_t \left( f(\*x^*) \leq \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{\delta} \sigma_{t-1}(\*x) \right) \geq 1 - \delta,
    \end{align*}
    where $\beta_{\delta} =  2 \log (|{\cal X}| / (2\delta))$.
    %
    Note that the probability is taken only with $f(\*x^*)$ since ${\cal D}_{t-1}$ is fixed.
    %
    Using the cumulative distribution function $F_t(\cdot) \coloneqq \mathrm{Pr}_t(f(\*x^*) \leq \cdot)$ and its inverse function $F^{-1}_t$, we can rewrite
    \begin{align*}
        F_t \left( \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{\delta} \sigma_{t-1}(\*x) \right) &\geq 1 - \delta \\
        \Longleftrightarrow \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{\delta} \sigma_{t-1}(\*x) &\geq F^{-1}_t ( 1 - \delta ),
    \end{align*}
    since $F^{-1}_t$ is a monotone increasing function.
    %
    %
    Since $\delta \in (0, 1)$ is arbitrary, we can substitute $U \sim {\rm Uni}(0, 1)$ as,
    \begin{align*}
         \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{U} \sigma_{t-1}(\*x) \geq F^{-1}_t ( 1 - U ).
    \end{align*}
    %
    By taking the expectation with respect to $U$, we can obtain,
    \begin{align*}
         \mathbb{E}_{U} \left[ \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{U} \sigma_{t-1}(\*x)   \right]
         &\geq \mathbb{E}_{U} \bigl[ F^{-1}_t ( 1 - U )   \bigr] \\
         &= \mathbb{E}_{U} \bigl[ F^{-1}_t ( U )   \bigr],
    \end{align*}
    where we use the fact that $1 - U$ also follows ${\rm Uni}(0, 1)$.
    %
    Since $F^{-1}_t ( U )$ and $f(\*x^*)$ are identically distributed as with the inverse transform sampling, we obtain
    \begin{align*}
         \mathbb{E}_{U} \left[ \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \beta^{1/2}_{U} \sigma_{t-1}(\*x)   \right]
         \geq \mathbb{E}_t \bigl[ f(\*x^*) \bigr].
    \end{align*}
    %
    Then, $\beta_U$ can be decomposed as follows:
    \begin{align*}
        \beta_U = 2 \log (|{\cal X}| / 2) - 2 \log(U).
    \end{align*}
    %
    From the inverse transform sampling, $ - 2 \log(U) \sim {\rm Exp} (\lambda = 1/2)$.
    %
    Hence, $\beta_U$ follows a two-parameter exponential distribution with $s = 2 \log (|{\cal X}| / 2)$ and $\lambda = 1/2$.
    %
    Therefore, we obtain the following:
    \begin{align*}
         \mathbb{E}_t \left[ \mathop{\rm argmax}\limits_{\*x \in {\cal X}} \mu_{t-1}(\*x) + \zeta^{1/2} \sigma_{t-1}(\*x) \right]
         \geq \mathbb{E}_t \bigl[ f(\*x^*)  \bigr].
    \end{align*}
\end{proof}

\section{Proof of Theorem~\ref{theo:BCR_IRGPUCB_discrete}}
\label{app:proof_discrete}

Using Lemma~\ref{lem:bound_RGPUCB}, we can obtain the proof of the following theorem for the finite input domain:
\begin{reptheorem}{theo:BCR_IRGPUCB_discrete}
    Let $f \sim {\mathscr G} {\cal P} (0, k)$, where $k$ is a stationary kernel and $k(\*x, \*x) = 1$, and ${\cal X}$ be finite.
    %
    Assume that $\zeta_t$ follows a two-parameter exponential distribution with $s = 2 \log (|{\cal X}| / 2)$ and $\lambda = 1/2$ for any $t \geq 1$.
    %
    Then, by running IRGP-UCB with $\zeta_t$, BCR can be bounded as follows:
    \begin{align*}
        {\rm BCR}_T \leq \sqrt{C_1 C_2 T \gamma_T},
    \end{align*}
    where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$ and $C_2 \coloneqq 2 + s$.
    %
\end{reptheorem}
\begin{proof}
    From Lemma~\ref{lem:bound_RGPUCB}, we obtain
    \begin{align}
        {\rm BCR}_T
        &= \sum_{t=1}^T \mathbb{E}[f(\*x^*) - f(\*x_t)] \\
        &= \sum_{t=1}^T \mathbb{E}[f(\*x^*) - V_t(\*x_t) + V_t(\*x_t) - f(\*x_t)] \\
        &\leq \sum_{t=1}^T \mathbb{E}[V_t(\*x_t) - f(\*x_t)].
    \end{align}
    %
    Then, we can easily derive the bound as follows:
    \begin{align}
        {\rm BCR}_T
        &\leq \sum_{t=1}^T \mathbb{E}[V_t(\*x_t) - f(\*x_t)] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} [ \mathbb{E} [V_t(\*x_t) - f(\*x_t) | {\cal D}_{t-1}, \zeta_t]] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}} [ \mathbb{E}_{\zeta_t} [V_t(\*x_t) ] - \mu_{t-1}(\*x_t)] \\
        &= \sum_{t=1}^T \mathbb{E}_{{\cal D}_{t-1}, \zeta_t} [ \zeta_t^{1/2} \sigma_{t-1}(\*x_t) ] \\
        &= \mathbb{E} \left[ \sum_{t=1}^T \zeta_t^{1/2} \sigma_{t-1}(\*x_t) \right] \\
        &\leq \mathbb{E} \left[\sqrt{ \sum_{t=1}^T \zeta_t \sum_{t=1}^T \sigma_{t-1}^2 (\*x_t) } \right] && \bigl( \because \text{Cauchy--Schwarz  inequality} \bigr) \\
        &\leq \mathbb{E} \left[\sqrt{ \sum_{t=1}^T \zeta_t } \right] \sqrt{C_1 \gamma_T} && \bigl( \because \text{Lemma~5.4 in \citep{Srinivas2010-Gaussian}} \bigr) \\
        &\leq \sqrt{ \sum_{t=1}^T \mathbb{E} [ \zeta_t ]} \sqrt{C_1 \gamma_T} && \bigl( \because \text{Jensen's inequality} \bigr) \\
        &= \sqrt{C_1 C_2 T \gamma_T}, && \bigl( \because \text{Definition of $\zeta_t$} \bigr)
    \end{align}
    where $C_2 = \mathbb{E} [ \zeta_t ] = 2 + 2\log(|{\cal X}|/2)$.
    %
    This concludes the proof.
\end{proof}

\begin{corollary}
    Assume the same condition as in Theorem~\ref{theo:BCR_IRGPUCB_discrete}.
    %
    Fix a required accuracy $\eta > 0$.
    %
    Then, ${\rm BSR}_T \leq \eta$ by at most $T$ function evaluations, where $T$ is the smallest positive integer satisfying the following inequality:
    \begin{align*}
        \sqrt{\frac{C_1 C_2 \gamma_T}{T}} \leq \eta,
    \end{align*}
    where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$ and $C_2 \coloneqq 2 + 2 \log (|{\cal X}| / 2)$.
\end{corollary}
\begin{proof}
    This BCR bound in Theorem~\ref{theo:BCR_IRGPUCB_discrete} immediately suggests that the following:
    \begin{align*}
        {\rm BSR}_T \leq \frac{{\rm BCR}_T }{T} \leq \sqrt{ \frac{C_1 C_2 \gamma_T}{T}}.
    \end{align*}
    %
    Therefore, if $\sqrt{ C_1 C_2 \gamma_T / T} < \eta$, then ${\rm BSR}_T < \eta$.
\end{proof}

\section{Proof of Theorem~\ref{theo:BCR_IRGPUCB_continuous}}
\label{app:proof_continuous}

\begin{reptheorem}{theo:BCR_IRGPUCB_continuous}
    Let $f \sim {\mathscr G} {\cal P} (0, k)$, where $k$ is a stationary kernel, $k(\*x, \*x) = 1$, and Assumption~\ref{assump:continuous_X} holds.
    %
    Assume that $\zeta_t$ follows a two-parameter exponential distribution with $s_t = 2d \log(bdr t^2 \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr)) - 2 \log 2$ and $\lambda = 1/2$ for any $t \geq 1$.
    %
    Then, by running IRGP-UCB, BCR can be bounded as follows:
    \begin{align*}
        {\rm BCR}_T \leq \frac{\pi^2}{6} + \sqrt{C_1 T \gamma_T (2 + s_T)},
    \end{align*}
    where $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
\end{reptheorem}
\begin{proof}
    Purely for the sake of analysis, we used a set of discretization ${\cal X}_t \subset {\cal X}$ for $t \geq 1$.
    %
    For any $t \geq 1$, let ${\cal X}_t \subset {\cal X}$ be a finite set with each dimension equally divided into $\tau_t = bdr t^2 \bigl(\sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)$.
    %
    Thus, $|{\cal X}_t| = \tau_t^d$.
    %
    In addition, we define $[\*x]_t$ as the nearest points in ${\cal X}_t$ of $\*x \in {\cal X}$.

    Then, we decompose BCR as follows:
    \begin{align}
        {\rm BCR}_T
        &= \sum_{t=1}^T \mathbb{E} \Biggl[ f(\*x^*) - f( [\*x^*]_t ) + f( [\*x^*]_t ) - \mathop{\rm argmax}\limits_{\*x \in {\cal X}_t} V_t(\*x) + \mathop{\rm argmax}\limits_{\*x \in {\cal X}_t} V_t(\*x) - V_t(\*x_t) + V_t(\*x_t) - f(\*x_t) \Biggl].
    \end{align}
    %
    Obviously, $\mathbb{E}\bigl[\mathop{\rm argmax}\limits_{\*x \in {\cal X}_t} V_t(\*x) - V_t(\*x_t) \bigr] \leq 0$ from the definition of $\*x_t$.
    %
    Furthermore, we observe the following:
    \begin{align*}
        \mathbb{E} \left[ f( [\*x^*]_t ) - \mathop{\rm argmax}\limits_{\*x \in {\cal X}_t} V_t(\*x) \right] \leq \mathbb{E} \left[ \mathop{\rm argmax}\limits_{\*x \in {\cal X}_t} f( \*x ) - \mathop{\rm argmax}\limits_{\*x \in {\cal X}_t} V_t(\*x) \right],
    \end{align*}
    which can be bounded above by $0$ by setting $s_t = 2\log(|{\cal X}_t| / 2)$ and $\lambda = 1 / 2$ from Lemma~\ref{lem:bound_RGPUCB}.
    %
    Therefore, we obtain the following:
    \begin{align}
        {\rm BCR}_T
        &\leq \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f( [\*x^*]_t ) + V_t(\*x_t) - f(\*x_t) \right] \\
        &= \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f( [\*x^*]_t ) \right] + \sum_{t=1}^T \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \right]. \label{eq:BCR_continuous_two_term_bound}
    \end{align}

    First, we consider the first term $ \sum_{t=1}^T \mathbb{E} \bigl[ f(\*x^*) - f( [\*x^*]_t ) \bigr]$.
    %
    From Lemma~\ref{lem:discretized_error}, we can obtain the following:
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f( [\*x^*]_t ) \right]
        &\leq \sum_{t=1}^T \mathbb{E} \left[ \sup_{\*x \in {\cal X}} f(\*x) - f( [\*x]_t ) \right] \\
        &\leq \sum_{t=1}^T \frac{1}{t^2} && \bigl( \because \text{ Lemma~\ref{lem:discretized_error} } \bigr) \\
        &\leq \frac{\pi^2}{6}, && \left(\because \sum_{t=1}^\infty \frac{1}{t^2} = \frac{\pi^2}{6} \right) \label{eq:discrete_points_distance_bound}
    \end{align}
    where $u_t$ in Lemma~\ref{lem:discretized_error} corresponds to $t^2$.

    The second term is bounded as follows:
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \right]
        &\leq \sqrt{ \sum_{t=1}^T \mathbb{E}[\zeta_t] C_1 \gamma_T }, && \bigl(\because \text{See the proof of Theorem~\ref{theo:BCR_IRGPUCB_discrete}} \bigr) \label{eq:continouos_bound_width_bound}
    \end{align}
    where
    \begin{align}
        \mathbb{E}[\zeta_t]
        &= 2 + 2\log(|{\cal X}_t| / 2) \\
        &= 2 + 2d \log(bdr t^2 \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr)) - 2 \log 2,
    \end{align}
    which is monotone increasing.
    %
    Therefore, we obtain the following:
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \right]
        &\leq \sqrt{ C_1 T \gamma_T \mathbb{E}[\zeta_T] },
    \end{align}
    where $\mathbb{E}[\zeta_T] = 2 + s_T = 2 + 2d \log(bdr T^2 \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr)) - 2 \log 2$.
    %
    Consequently, combining Eqs.~\eqref{eq:BCR_continuous_two_term_bound}, \eqref{eq:discrete_points_distance_bound}, and \eqref{eq:continouos_bound_width_bound} concludes the proof.
\end{proof}

\begin{corollary}
    Assume the same condition as in Theorem~\ref{theo:BCR_IRGPUCB_continuous}.
    %
    Fix a required accuracy $\eta > 0$.
    %
    Then, by running IRGP-UCB with $s_\eta = 2d \log(2bdr \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr) / \eta) - 2 \log 2$ and $\lambda = 1 / 2$, ${\rm BSR}_T \leq \eta$ by at most $T$ function evaluations, where $T$ is the smallest positive integer satisfying the following inequality:
    \begin{align*}
        \sqrt{\frac{C_1 (2 + s_\eta) \gamma_T}{T}} \leq \frac{\eta}{2},
    \end{align*}
    where  $C_1 \coloneqq 2 / \log(1 + \sigma^{-2})$.
\end{corollary}
\begin{proof}
    We obtain the upper bound of BSR via the upper bound of BCR as follows:
    \begin{align}
        {\rm BSR}_T \leq \frac{{\rm BCR}_T}{T}.
    \end{align}
    %
    Let $\tau_t = 2bdr \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr) / \eta$.
    %
    Until the derivation of Eq.~\eqref{eq:BCR_continuous_two_term_bound}, the same derivation as in Theorem~\ref{theo:BCR_IRGPUCB_continuous} can be applied.
    %
    That is,
    \begin{align}
        {\rm BCR}_T
        &\leq \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f( [\*x^*]_t ) \right] + \sum_{t=1}^T \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \right].
    \end{align}
    %
    Then, for the first term, we can obtain the following
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ f(\*x^*) - f( [\*x^*]_t ) \right]

        &\leq \sum_{t=1}^T \frac{\eta}{2} && \bigl( \because \text{ Lemma~\ref{lem:discretized_error} } \bigr) \\
        &\leq \frac{\eta}{2} T.
    \end{align}
    %
    For the second term,
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ V_t(\*x_t) - f(\*x_t) \right]
        &\leq \sqrt{ C_1 \sum_{t=1}^T \mathbb{E}[\zeta_t] \gamma_T },
    \end{align}
    where, from the assumption, $\mathbb{E} [\zeta_t] = 2 + s_\eta$ for any $t \geq 1$.
    %
    Thus, we can derive the bound of BSR:
    \begin{align}
        {\rm BSR}_T \leq \frac{\eta}{2} + \sqrt{\frac{C_1 (2 + s_\eta) \gamma_T}{T}}.
    \end{align}
    %
    Hence, we obtain the smallest integer $T$ that satisfies ${\rm BSR}_T \leq \eta$ by arranging the following inequality:
    \begin{align*}
        \frac{\eta}{2} + \sqrt{\frac{C_1 (2 + s_\eta) \gamma_T}{T}} \leq \eta.
    \end{align*}
\end{proof}

\section{Auxiliary Lemmas}
\label{app:lemmas}

For convenience, we here show the assumption again:
\begin{repassumption}{assump:continuous_X}
    Let ${\cal X} \subset [0, r]^d$ be a compact and convex set, where $r > 0$.
    %
    Assume that the kernel $k$ satisfies the following condition on the derivatives of sample path $f$.
    %
    There exists the constants $a, b > 0$ such that,
    \begin{align*}
        \Pr \left( \sup_{\*x \in {\cal X}} \left| \frac{\partial f}{\partial \*x_j} \right| > L \right) \leq a \exp \left( - \left(\frac{L}{b}\right)^2 \right),\text{ for } j \in [d],
    \end{align*}
    where $[d] = \{1, \dots, d\}$.
\end{repassumption}

Then, from Assumption~\ref{assump:continuous_X}, we obtain several lemmas.
First, we show the upper bound of the supremum of the partial derivatives.
This result is tighter than Lemma~12 in \citep{Kandasamy2018-Parallelised}.
Note that, in the proof of Lemma~12 in \citep{Kandasamy2018-Parallelised}, there is a typo that $\Pr (L \geq t)$ is bounded above by $a \exp\{ t^2 / b^2 \}$, which should be $a d \exp\{ - t^2 / b^2 \}$.
\begin{lemma}
    Let $f \sim {\mathscr G} {\cal P} (0, k)$ and Assumption~\ref{assump:continuous_X} holds.
    %
    Let the supremum of the partial derivatives $L_{\rm max} \coloneqq \sup_{j \in [d]} \sup_{\*x \in {\cal X}} \left| \frac{\partial f}{\partial x_j} \right|$.
    %
    Then, $\mathbb{E}[L_{\rm max}]$ can be bounded above as follows:
    \begin{align}
        \mathbb{E}[L_{\rm max}] \leq b \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr).
    \end{align}
    \label{lem:L_max_bound}
\end{lemma}
\begin{proof}
    From the assumption, using union bound, we obtain the bound of the probability,
    \begin{align}
        \Pr\left( L_{\rm max} > L \right) \leq \sum_{j=1}^d a \exp \left( - \left(\frac{L}{b}\right)^2 \right) = a d \exp \left( - \left(\frac{L}{b}\right)^2 \right). \label{eq:L_max_prob}
    \end{align}
    %
    Then, using expectation integral identity, we can bound $\mathbb{E}[L_{\rm max}]$ as follows:
    \begin{align}
        \mathbb{E}[L_{\rm max}]
        &= \int_0^{\infty} \Pr\left( L_{\rm max} > L \right) {\rm d} L \\
        &\leq \int_0^{\infty} \mathop{\rm argmin}\limits \{ 1, ad e^{- (L / b)^2} \} {\rm d} L && \bigl(\because \text{Eq.~\eqref{eq:L_max_prob}} \bigr) \\
        &= b \sqrt{\log (ad)} + \int_{b\sqrt{\log (ad)}}^{\infty} ad e^{- (L / b)^2} {\rm d} L \\
        &= b \sqrt{\log (ad)} + a b d \sqrt{\pi} \int_{b\sqrt{\log (ad)}}^{\infty} \frac{1}{\sqrt{2\pi (b^2/2)}} e^{- (L / b)^2} {\rm d} L \\
        &= b \sqrt{\log (ad)} + a b d \sqrt{\pi} \left( 1 - \Phi \left( \frac{b\sqrt{\log (ad)}}{ b / \sqrt{2}} \right) \right) \\
        &\leq b \sqrt{\log (ad)} + \frac{ b \sqrt{\pi}}{2}, && \bigl(\because \text{Lemma~\ref{lem:Gauss_tail_bound}} \bigr)
    \end{align}
    where $\Phi$ is a cumulative distribution function of the standard normal distribution.
\end{proof}

Next, we show the bound for the difference of discretized outputs:
\begin{lemma}
    Let $f \sim {\mathscr G} {\cal P} (0, k)$ and Assumption~\ref{assump:continuous_X} holds.
    %
    Let ${\cal X}_t \subset {\cal X}$ be a finite set with each dimension equally divided into $\tau_t = bdr u_t \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2 \bigr)$ for any $t \geq 1$.
    %
    Then, we can bound the expectation of differences,
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ \sup_{\*x \in {\cal X}} f(\*x) - f( [\*x]_t ) \right] \leq \sum_{t=1}^T \frac{1}{u_t},
    \end{align}
    where $[\*x]_t$ is the nearest point in ${\cal X}_t$ of $\*x \in {\cal X}$.
    \label{lem:discretized_error}
\end{lemma}
\begin{proof}
    From the construction of ${\cal X}_t$, we can obtain the upper bound of L1 distance between $\*x$ and $[\*x]_t$ as follows:
    \begin{align}
        \sup_{\*x \in {\cal X}}\| \*x - [\*x]_t \|_1
        &\leq \frac{dr }{bdr u_t \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr) } \\
        &= \frac{1}{b u_t \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr)}.
        \label{eq:sup_L1dist_discretization}
    \end{align}
    %
    Note that this discretization does not depend on any randomness and is fixed beforehand.

    Then, we obtain the following:
    \begin{align}
        \sum_{t=1}^T \mathbb{E} \left[ \sup_{\*x \in {\cal X}} f(\*x) - f( [\*x]_t ) \right]
        &\leq \sum_{t=1}^T \mathbb{E} \left[ L_{\rm max} \sup_{\*x \in {\cal X}} \| \*x - [\*x]_t \|_1 \right] \\
        &\leq \sum_{t=1}^T \mathbb{E} \left[ L_{\rm max} \right] \frac{1}{b u_t \bigl( \sqrt{\log (ad)} + \sqrt{\pi} / 2\bigr)} && \bigl(\because \text{Eq.~\eqref{eq:sup_L1dist_discretization}} \bigr) \\
        &\leq \sum_{t=1}^T \frac{1}{u_t}. && \bigl(\because \text{Lemma~\ref{lem:L_max_bound}} \bigr)
    \end{align}
\end{proof}

We used the following useful lemma:
\begin{lemma}[in Lemma 5.2 of \citep{Srinivas2010-Gaussian}]
    For $c > 0$, the survival function of the standard normal distribution can be bounded above as follows:
    \begin{align}
        1 - \Phi(c) \leq \frac{1}{2} \exp (- c^2 / 2).
    \end{align}
    \label{lem:Gauss_tail_bound}
\end{lemma}
\begin{proof}
    Let $R \sim {\cal N}(0, 1)$.
    %
    Then, we obtain the following:
    \begin{align}
        \Pr(R > c) &= \int_{c}^\infty \frac{1}{\sqrt{2 \pi}} \exp ( - r^2 / 2) {\rm d}r \\
        &= \exp (- c^2 / 2) \int_{c}^\infty \frac{1}{\sqrt{2 \pi}} \exp ( - r^2 / 2 + c^2 / 2) {\rm d}r \\
        &= \exp (- c^2 / 2) \int_{c}^\infty \frac{1}{\sqrt{2 \pi}} \exp ( - (r - c)^2 / 2 - rc + c^2) {\rm d}r \\
        &= \exp (- c^2 / 2) \int_{c}^\infty \frac{1}{\sqrt{2 \pi}} \exp ( - (r - c)^2 / 2 - c(r - c)) {\rm d}r \\
        &\leq \exp (- c^2 / 2) \int_{c}^\infty \frac{1}{\sqrt{2 \pi}} \exp ( - (r - c)^2 / 2) {\rm d}r && \left( \because \text{$e^{-c(r-c)} \leq 1$ since $r \geq c > 0$} \right) \\
        &= \frac{1}{2} \exp (- c^2 / 2).
    \end{align}
    %

\end{proof}