\section{Introduction}
Consider the problem of maximizing an expensive-to-compute black-box objective function $f$ that depends on an \emph{optimization variable} $\mbf{x}$ and an \emph{environmental random variable} $\mbf{Z}$.
Due to the randomness in $\mbf{Z}$, the function evaluation $f(\mbf{x},\mbf{Z})$ of $f$ at $\mbf{x}$ is a random variable.
Though for such an objective function $f$, \emph{Bayesian optimization} (BO) can be naturally applied to maximize its expectation $\mbb{E}_{\mbf{Z}}[f(\mbf{x},\mbf{Z})]$ over $\mbf{Z}$ \cite{toscano2018bayesian}, this maximization objective overlooks the \emph{risks} of potentially undesirable function evaluations.
These risks can arise from either (a) the realization of an unknown distribution of $\mbf{Z}$ or (b) the realization of the random $\mbf{Z}$ given that the distribution of $f(\mbf{x},\mbf{Z})$ can be estimated well or that of $\mbf{Z}$ is known.
The issue (a) has been tackled by distributionally robust BO \cite{kirschner2020distributionally,nguyen2020distributionally} which maximizes $\mbb{E}_{\mbf{Z}}[f(\mbf{x},\mbf{Z})]$ under the worst-case realization of the distribution of $\mbf{Z}$.
To resolve the issue (b),
the risk from the uncertainty in $\mbf{Z}$ can be controlled via the mean-variance optimization framework \cite{iwazaki2020mean}, \emph{value-at-risk} ({\textsc{VaR}}), or \emph{conditional value-at-risk} ({\textsc{CVaR}}) \cite{borisk20,torossian2020bayesian}.
The work of \citet{bogunovic2018adversarially} has considered \emph{adversarially robust BO},
where $\mbf{z}$ is controlled by an adversary deterministically.\footnote{We use upper-case letter $\mbf{Z}$ to denote the environmental random variable and lower-case letter $\mbf{z}$ to denote its realization or a (non-random) variable.}
In this case, the objective is to find $\mbf{x}$ that maximizes the function under the worst-case realization of $\mbf{z}$, i.e., $\text{arg}\!\max_{\mbf{x}} \min_{\mbf{z}} f(\mbf{x},\mbf{z})$.

In this paper, we focus on case (b) where the distribution of
$\mbf{Z}$ is known (or well-estimated).
For example, in agriculture, although farmers cannot control the temperature of an outdoor farm,
its distribution
can be estimated from historical data and controlled in an indoor environment for optimizing the plant yield.
Given the distribution of $\mbf{Z}$, the objective is to control the risk
that the function evaluation $f(\mbf{x},\mbf{z})$,
for a $\mbf{z}$ sampled from $\mbf{Z}$,
is small.
One popular framework is to control the trade-off between the mean (viewed as reward) and the variance (viewed as risk) of the function evaluation with respect to $\mbf{Z}$ \cite{iwazaki2020mean}.
However, quantifying the risk
using variance implies indifference between positive and negative deviations from the mean, while people often have asymmetric risk attitudes \cite{goh2012portfolio}. In our problem of maximizing the objective function, it is reasonable to assume that people  are risk-averse towards only the negative deviations
from the mean, i.e.,
the risk of getting lower function evaluations.
Thus, it is more appropriate to adopt
risk measures
with this asymmetric property, such as \emph{value-at-risk} ({\textsc{VaR}})
which is a widely adopted risk measure in real-world applications (e.g., banking \cite{basel06}).
Intuitively, the risk
that the random $f(\mbf{x},\mbf{Z})$ is less than {\textsc{VaR}} at level $\alpha \in (0,1)$ does not exceed $\alpha$,
e.g., by specifying a small value of $\alpha$ as $0.1$, this risk is controlled to be at most $10\%$.
Therefore, to maximize the function $f$ while controlling the risk of undesirable (i.e., small) function evaluations,
we aim to maximize {\textsc{VaR}} of the random function $f(\mbf{x},\mbf{Z})$ over $\mbf{x}$.

The recent work of \citet{borisk20} has used BO to maximize {\textsc{VaR}} and has achieved state-of-the-art empirical performances.
They have assumed that we are able to select both $\mbf{x}$ and $\mbf{z}$ to query during BO,
which is motivated by fact that physical experiments can usually be studied by
simulation \cite{williams2000sequential}.
In the example on agriculture given above, we can control the temperature, light and water ($\mbf{z}$) in a small indoor environment to optimize the amount of fertilizer ($\mbf{x}$), which can then be used in an outdoor environment with random weather factors.
\citet{borisk20} have exploited the ability to select $\mbf{z}$ to model the function $f(\mbf{x},\mbf{z})$ as a GP, which allows them to retain the appealing closed-form posterior belief of the objective function.
To select the queries $\mbf{x}$ and $\mbf{z}$, they have designed a one-step lookahead approach based on the well-known \emph{knowledge gradient} (KG) acquisition function \cite{scott2011correlated}.
However, the one-step lookahead incurs an expensive nested optimization procedure, which is computationally expensive and hence requires approximations.
Besides, the acquisition function can only be approximated using samples of the objective function $f$ from the GP posterior and the environmental random variable $\mbf{Z}$.
While they have analysed the asymptotically unbiased and consistent estimator of the gradients, it is challenging to obtain a guarantee for the convergence of their algorithm.
Another recent work \cite{torossian2020bayesian} has also applied BO to maximize {\textsc{VaR}} using an asymmetric Laplace likelihood function and variational approximation of the posterior belief. However, in contrast to \citet{borisk20} and our work, they have focused on a different setting where the realizations of $\mbf{Z}$ are not observed.

In this paper, we adopt the setting of \citet{borisk20} which allows us to choose both $\mbf{x}$ and $\mbf{z}$ to query, and assume that the distribution of $\mbf{Z}$ is known or well-estimated.
Our contributions include:

\textbf{Firstly}, we propose a novel BO algorithm named \emph{Value-at-risk Upper Confidence Bound} ({V-UCB}) in Section~\ref{sec:vy}.
Unlike the work of \citet{borisk20}, {V-UCB} is equipped with a no-regret convergence guarantee and is more computationally efficient.
To guide its query selection and facilitate its proof of the no-regret guarantee, the classical GP-UCB algorithm~\cite{srinivas10ucb} constructs a \emph{confidence bound} of the objective function. Similarly, to maximize the {\textsc{VaR}} of a random function, we, for the first time to the best of our knowledge, construct a confidence bound of {\textsc{VaR}} (Lemma~\ref{lemma:confbound}).
The resulting confidence bound of {\textsc{VaR}} naturally gives rise to a strategy to select $\mbf{x}$.
However, it remains a major challenge to select $\mbf{z}$ to preserve the no-regret convergence of GP-UCB.
To this end, we firstly prove that our algorithm is no-regret as long as we ensure that at the selected $\mbf{z}$, the confidence bound of {\textsc{VaR}} \emph{lies within} the confidence bound of the objective function.
Next, we also prove that this query selection strategy is \emph{feasible}, i.e., such values of $\mbf{z}$, referred to as \emph{lacing values} (LV), exist

\textbf{Secondly}, although our theoretical no-regret property allows the selection of \emph{any} LV, we design a heuristic to select an LV such that it improves our empirical performance over random selection of LV (Section~\ref{subsec:improvedselectz}).
We also discuss the implications when $\mbf{z}$ cannot be selected by BO and is instead randomly sampled by the environment during BO (Remark~\ref{rmk:unknownZ}).
\textbf{Thirdly}, we show that adversarially robust BO \cite{bogunovic2018adversarially} can be cast as a special case of our V-UCB
when the risk level $\alpha$ of {\textsc{VaR}} approaches $0$ from the right and the domain of $\mbf{z}$ is the support of $\mbf{Z}$.
In this case,
adversarially robust BO \cite{bogunovic2018adversarially} selects the same input queries as those selected by V-UCB since the set of LV collapse into the set of minimizers of the lower bound of the objective function (Section~\ref{subsec:vy-stableopt}).
\textbf{Lastly}, we provide practical techniques for
implementing {V-UCB} with continuous random variable $\mbf{Z}$ (Section~\ref{sec:continuousz}): we (a) introduce \emph{local neural surrogate optimization} with the \emph{pinball loss} to optimize {\textsc{VaR}}, and (b) construct an objective function to search for an LV in the continuous support of $\mbf{Z}$.

The performance of our proposed algorithm is empirically demonstrated in optimizing several synthetic benchmark functions, a portfolio optimization problem, and a simulated robot task in Section~\ref{sec:experiments}.

\section{Problem Statement and Background}

Let the objective function be defined as $f: \mcl{D}_{\mbf{x}} \times \mcl{D}_{\mbf{z}} \rightarrow \mbb{R}$ where $\mcl{D}_{\mbf{x}} \subset \mbb{R}^{d_x}$ and $\mcl{D}_{\mbf{z}}\subset \mbb{R}^{d_z}$ are the bounded domain of the optimization variable $\mbf{x}$ and the support of the environmental random variable $\mbf{Z}$, respectively; $d_x$ and $d_z$ are the dimensions of $\mbf{x}$ and $\mbf{z}$, respectively. The support of $\mbf{Z}$ is defined as the smallest closed subset $\mcl{D}_{\mbf{z}}$ of $\mbb{R}^{d_z}$ such that $P(\mbf{Z} \in \mcl{D}_{\mbf{z}}) = 1$.
Let $\mbf{z} \in \mcl{D}_{\mbf{z}}$ denote a realization of the random variable $\mbf{Z}$.
Let $f(\mbf{x},\mbf{Z})$ denote a random variable whose randomness comes from $\mbf{Z}$.
The {\textsc{VaR}} of $f(\mbf{x},\mbf{Z})$ at \emph{risk level} $\alpha \in (0,1)$ is defined as:
\begin{equation}
V_{\alpha}(f(\mbf{x},\mbf{Z})) \triangleq \inf \{\omega: P(f(\mbf{x},\mbf{Z}) \le \omega) \ge \alpha\}
\label{eq:var}
\end{equation}
which implies the risk that $f(\mbf{x},\mbf{Z})$ is less than its {\textsc{VaR}} at level $\alpha$ does not exceed $\alpha$.

Our objective is to search for $\mbf{x} \in \mcl{D}_{\mbf{x}}$ that maximizes $V_\alpha(f(\mbf{x},\mbf{Z}))$ at a user-specified risk level $\alpha \in (0,1)$.
Intuitively, the goal is find $\mbf{x}$ where the evaluations of the objective function are as large as possible under most realizations of the environmental random variable $\mbf{Z}$ which is characterized by the probability of $1-\alpha$.

The unknown objective function $f(\mbf{x},\mbf{z})$ is modeled with a GP. That is, every finite subset of $\{f(\mbf{x},\mbf{z})\}_{(\mbf{x}, \mbf{z}) \in \mcl{D}_{\mbf{x}} \times \mcl{D}_{\mbf{z}}}$ follows a multivariate Gaussian distribution \cite{rasmussen06}.
The GP is fully specified by its \emph{prior} mean and covariance function $k_{(\mbf{x},\mbf{z}), (\mbf{x}',\mbf{z}')} \triangleq \text{cov}[f(\mbf{x},\mbf{z}), f(\mbf{x}',\mbf{z}')]$ for all $\mbf{x}, \mbf{x}'$ in $\mcl{D}_{\mbf{x}}$ and $\mbf{z}, \mbf{z}'$ in $\mcl{D}_{\mbf{z}}$. For notational simplicity (and w.l.o.g.), the former is assumed to be zero, while we use the \emph{squared exponential} (SE) kernel as its bounded maximum information gain can be used for later analysis \cite{srinivas10ucb}.

To identify the optimal $\mbf{x}_* \triangleq \text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_{\alpha}(f(\mbf{x},\mbf{Z}))$, BO algorithm selects an input query $(\mbf{x}_t,\mbf{z}_t)$ in the $t$-th iteration to obtain a noisy function evaluation $y_{(\mbf{x}_t,\mbf{z}_t)} \triangleq f(\mbf{x}_t,\mbf{z}_t) + \epsilon_t$ where $\epsilon_t \sim \mcl{N}(0,\sigma_n^2)$ are i.i.d. Gaussian noise with variance $\sigma_n^2$. Given noisy observations $\mbf{y}_{\mcl{D}_t} \triangleq (y_{(\mbf{x},\mbf{z})})_{(\mbf{x},\mbf{z}) \in \mcl{D}_t}^\top$ at observed inputs $\mcl{D}_t \triangleq \mcl{D}_{t-1} \cup \{(\mbf{x}_t, \mbf{z}_t)\}$ (and $\mcl{D}_0$ is the initial observed inputs), the GP posterior belief of function evaluation at any input $(\mbf{x},\mbf{z})$ is a Gaussian $p(f(\mbf{x},\mbf{z})|\mbf{y}_{\mcl{D}_t}) \triangleq \mcl{N}( f(\mbf{x},\mbf{z})|\mu_t(\mbf{x},\mbf{z}), \sigma_t^2(\mbf{x},\mbf{z}))$:
\begin{equation}
\begin{array}{r@{}l}
\mu_t(\mbf{x},\mbf{z}) &\triangleq \mbf{K}_{(\mbf{x},\mbf{z}),\mcl{D}_t} \bm{\Lambda}_{\mcl{D}_t\mcl{D}_t} \mbf{y}_{\mcl{D}_t}\ ,\\
\sigma_t^2(\mbf{x},\mbf{z}) &\triangleq k_{(\mbf{x},\mbf{z}),(\mbf{x},\mbf{z})}  - \mbf{K}_{(\mbf{x},\mbf{z}),\mcl{D}_t} \bm{\Lambda}_{\mcl{D}_t\mcl{D}_t} \mbf{K}_{\mcl{D}_t,(\mbf{x},\mbf{z})}
\label{eq:gppost}
\end{array}
\end{equation}
where $\bm{\Lambda}_{\mcl{D}_t\mcl{D}_t} \triangleq \left(
	\mbf{K}_{\mcl{D}_t\mcl{D}_t} + \sigma_n^2 \mbf{I}
\right)^{-1}$,
$\mbf{K}_{(\mbf{x},\mbf{z}),\mcl{D}_t} \triangleq (k_{(\mbf{x},\mbf{z}),(\mbf{x}',\mbf{z}')})_{(\mbf{x}',\mbf{z}') \in \mcl{D}_t}$,
$\mbf{K}_{\mcl{D}_t,(\mbf{x},\mbf{z})} \triangleq \mbf{K}_{(\mbf{x},\mbf{z}),\mcl{D}_t}^\top$,
$\mbf{K}_{\mcl{D}_t\mcl{D}_t} \triangleq (k_{(\mbf{x}',\mbf{z}'), (\mbf{x}'',\mbf{z}'')})_{(\mbf{x}',\mbf{z}'), (\mbf{x}'',\mbf{z}'') \in \mcl{D}_t}$,  $\mbf{I}$ is the identity matrix.

\section{BO of {\textsc{VaR}}}
\label{sec:vy}

Following the seminal work \cite{srinivas10ucb}, we use the \emph{cumulative regret} as the performance metric to quantify the performance of our BO algorithm. It is defined as $R_T \triangleq \sum_{t=1}^T r(\mbf{x}_t)$ where $r(\mbf{x}_t) \triangleq V_{\alpha}(f(\mbf{x}_*,\mbf{Z})) - V_{\alpha}(f(\mbf{x}_t,\mbf{Z}))$ is the \emph{instantaneous regret} and $\mbf{x}_* \triangleq \text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_{\alpha}(f(\mbf{x},\mbf{Z}))$. We would like to design a query selection strategy that incurs \emph{no regret}, i.e., $\lim_{T\rightarrow \infty} R_T / T = 0$.
Furthermore, we have that
$\min_{t \le T} r(\mbf{x}_t) \le R_T / T$, equivalently, $\max_{t \le T} V_{\alpha}(f(\mbf{x}_t,\mbf{Z})) \ge V_{\alpha}(f(\mbf{x}_*,\mbf{Z})) - R_T / T$.
Thus, $\text{lim}_{T \rightarrow \infty} \max_{t \le T} V_{\alpha}(f(\mbf{x}_t,\mbf{Z})) = V_{\alpha}(f(\mbf{x}_*,\mbf{Z}))$ for a no-regret algorithm.

The proof of the upper bound on the cumulative regret of GP-UCB is based on confidence bounds of the objective function \cite{srinivas10ucb}.
Similarly, in the next section, we start by constructing a confidence bound of $V_{\alpha}(f(\mbf{x},\mbf{Z}))$, which naturally leads to a query selection strategy for $\mbf{x}_t$.
\subsection{A Confidence Bound of $V_\alpha(f(\mbf{x},\mbf{Z}))$ and the Query Selection Strategy for $\mbf{x}_t$}
\label{subsec:confboundf}
Firstly, we adopt a confidence bound of the function $f(\mbf{x},\mbf{z})$ from \citet{chowdhury2017kernelized}, which assumes that $f$ belongs to a \emph{reproducing kernel Hilbert space} $\mcl{F}_k(B)$ such that its RKHS norm is bounded $\Vert f\Vert_k \le B$.

\begin{lemma}[\citet{chowdhury2017kernelized}]
Pick $\delta \in (0,1)$ and set
$\beta_t = (B + \sigma_n \sqrt{2(\gamma_{t-1} + 1 + \log1/ \delta)})^2$.
Then,
$f(\mbf{x},\mbf{z}) \in I_{t-1}[f(\mbf{x},\mbf{z})] \triangleq [l_{t-1}(\mbf{x},\mbf{z}), u_{t-1}(\mbf{x},\mbf{z})]$
$\forall \mbf{x} \in \mcl{D}_{\mbf{x}}, \mbf{z} \in \mcl{D}_{\mbf{z}}, t \ge 1 $ holds with probability $\ge 1 - \delta$ where
\begin{equation}
\begin{array}{r@{}l}
l_{t-1}(\mbf{x},\mbf{z})\ &\triangleq \mu_{t-1}(\mbf{x},\mbf{z}) - \beta_t^{1/2} \sigma_{t-1}(\mbf{x},\mbf{z})\\
u_{t-1}(\mbf{x},\mbf{z})\ &\triangleq \mu_{t-1}(\mbf{x},\mbf{z}) + \beta_t^{1/2} \sigma_{t-1}(\mbf{x},\mbf{z})\ .
\end{array}\vspace{-1mm}
\label{eq:fbound}
\end{equation}
\label{lemma:ucb51}
\end{lemma}
As the above lemma holds for both finite and continuous $\mcl{D}_x$ and $\mcl{D}_z$, it is
used
to analyse the regret in both cases. On the other hand, the confidence bound can be adopted to the Bayesian setting by changing only $\beta_t$ following the work of \citet{srinivas10ucb} as noted by \cite{bogunovic2018adversarially}.

Then, we exploit this confidence bound on the function evaluations (Lemma \ref{lemma:ucb51}) to formulate a confidence bound of $V_{\alpha}(f(\mbf{x},\mbf{Z}))$ as follows.
\begin{lemma}
\label{lemma:confbound}
Similar to the definition of $f(\mbf{x},\mbf{Z})$, let $l_{t-1}(\mbf{x}, \mbf{Z})$ and $u_{t-1}(\mbf{x},\mbf{Z})$ denote the random function over $\mbf{x}$ where the randomness comes from the random variable $\mbf{Z}$; $l_{t-1}$ and $u_{t-1}$ are defined in \eqref{eq:fbound}.
Then, $\forall \mbf{x} \in \mcl{D}_{\mbf{x}}$, $t \ge 1$,
\[
\begin{array}{r@{}l}
\displaystyle V_\alpha(f(\mbf{x},\mbf{Z})) \displaystyle &\displaystyle \in I_{t-1}[V_{\alpha}(f(\mbf{x},\mbf{Z}))]\\
	 &\displaystyle \triangleq [V_\alpha(l_{t-1}(\mbf{x},\mbf{Z})), V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))]
\end{array}
\]
holds with probability $\ge 1 - \delta$ for $\delta$ in Lemma~\ref{lemma:ucb51}, where $V_\alpha(l_{t-1}(\mbf{x},\mbf{Z}))$ and $V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$ are defined as \eqref{eq:var}.
\end{lemma}
The proof is in Appendix~\ref{app:proofvconfbound}.
Given the confidence bound $I_{t-1}[V_{\alpha}(f(\mbf{x},\mbf{Z}))] \triangleq [V_\alpha(l_{t-1}(\mbf{x},\mbf{Z})), V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$ in Lemma~\ref{lemma:confbound}, we follow the the well-known
``optimism in the face of uncertainty''
principle to select $\mbf{x}_t = \text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$. This query selection strategy for $\mbf{x}_t$ leads to an upper bound of $r(\mbf{x}_t)$:
\begin{align}
r(\mbf{x}_t) \le V_\alpha(u_{t-1}(\mbf{x}_t,\mbf{Z})) - V_\alpha(l_{t-1}(\mbf{x}_t,\mbf{Z}))\ \forall t \ge 1
	\label{eq:iregretbound1}
\end{align}
which holds with probability $\ge 1 - \delta$ for $\delta$ in Lemma~\ref{lemma:ucb51}, and is proved in Appendix~\ref{app:iregretbound1}.

As our goal is $\lim_{T\rightarrow \infty} R_T/T = 0$, given the selected query $\mbf{x}_t$, a reasonable query selection strategy of $\mbf{z}_t$ should gather informative observations at $(\mbf{x}_t, \mbf{z}_t)$ that improves the confidence bound $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$ (i.e., $I_{t}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$ is a proper subset of $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$ if $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))] \neq \emptyset$) which can be viewed as the uncertainty of $V_{\alpha}(f(\mbf{x}_t,\mbf{Z}))$.

Assume that there exists
$\mbf{z}_l \in \mcl{D}_{\mbf{z}}$ such that  $l_{t-1}(\mbf{x}_t,\mbf{z}_l) = V_\alpha(l_{t-1}(\mbf{x}_t,\mbf{Z}))$
and
$\mbf{z}_u \in \mcl{D}_{\mbf{z}}$ such that  $u_{t-1}(\mbf{x}_t,\mbf{z}_u) = V_\alpha(u_{t-1}(\mbf{x}_t,\mbf{Z}))$.
Lemma~\ref{lemma:confbound} implies that $V_\alpha(f(\mbf{x}_t,\mbf{Z})) \in I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))] = [l_{t-1}(\mbf{x}_t,\mbf{z}_l), u_{t-1}(\mbf{x}_t,\mbf{z}_u)]$ with high probability.
Hence, we may
na\"ively
want to query for observations at $(\mbf{x}_t, \mbf{z}_l)$ and $(\mbf{x}_t,\mbf{z}_u)$ to reduce $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$.
However, these observations may not always reduce $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$.
The insight is that $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$ changes (i.e., shrinks) when either of its boundary values (i.e., $l_{t-1}(\mbf{x}_t,\mbf{z}_l)$ or $u_{t-1}(\mbf{x}_t,\mbf{z}_u)$) changes.
Consider $u_{t-1}(\mbf{x}_t,\mbf{z}_u)$ and finite $\mcl{D}_{\mbf{z}}$ as an example, since $u_{t-1}(\mbf{x}_t,\mbf{z}_u) = V_\alpha(u_{t-1}(\mbf{x}_t,\mbf{Z}))$,
a natural cause for the change in $u_{t-1}(\mbf{x}_t,\mbf{z}_u)$ is when $\mbf{z}_u$ changes.
This happens if there exists $\mbf{z}' \neq \mbf{z}_u$ such that the \emph{ordering} of $u_{t-1}(\mbf{x}_t,\mbf{z}')$ relative to $u_{t-1}(\mbf{x}_t,\mbf{z}_u)$ changes given more observations. Thus, observations that are capable of reducing $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$ should be able to \emph{change the relative ordering} in this case.
We construct the following counterexample where observations at $\mbf{z}_u$ (and $\mbf{z}_l$) are not able to change the relative ordering, so they do not reduce $I_{t-1}[V_\alpha (f(\mbf{x}_t,\mbf{Z}))]$.

\begin{example}
\label{example:counter}
This example is described by Fig.~\ref{fig:orderinguncertainty}.
We reduce notational clutter by removing $\mbf{x}_t$ and $t$ since they are fixed in this example, i.e., we use $f(\mbf{z})$, $f(\mbf{Z})$, and $l(\mbf{z})$ to denote $f(\mbf{x}_t, \mbf{z})$, $f(\mbf{x}_t, \mbf{Z})$, and $l_{t-1}(\mbf{x}_t, \mbf{z})$ respectively.
We condition on the event $f(\mbf{z}) \in I[f(\mbf{z})] \triangleq [l(\mbf{z}), u(\mbf{z})]$ for all $\mbf{z} \in \mcl{D}_{\mbf{z}}$ which occurs with probability $\ge 1 - \delta$ in Lemma~\ref{lemma:ucb51}. In this example, $\mbf{z}_l = \mbf{z}_1$ and $l(\mbf{z}_1) = u(\mbf{z}_1)$, so there is no uncertainty in $f(\mbf{z}_l) = f(\mbf{z}_1)$. Similarly, there is no uncertainty in $f(\mbf{z}_u) = f(\mbf{z}_2)$. Thus, new observations at $\mbf{z}_l$ and $\mbf{z}_u$ change neither $l(\mbf{z}_l)$ nor $u(\mbf{z}_u)$, so these observations do not reduce the confidence bound $I[V_{\alpha=0.4}(f(\mbf{Z}))] = [l(\mbf{z}_l), u(\mbf{z}_u)]$ (plotted as the double-headed arrow in Fig.~\ref{fig:orderinguncertainty}b).
In fact, to reduce $I[V_{\alpha=0.4}(f(\mbf{Z}))]$, we should gather new observations at $\mbf{z}_0$ which potentially change the ordering of $u(\mbf{z}_0)$ relative to $u(\mbf{z}_2)$ (which is $u(\mbf{z}_u)$ without new observations). For example, after
getting
new observations at $\mbf{z}_0$, if $u(\mbf{z}_0)$ is improved to be in the white region between A and B ($u(\mbf{z}_0) > u(\mbf{z}_2)$ in Fig.~\ref{fig:orderinguncertainty}b changes to $u(\mbf{z}_0) < u(\mbf{z}_2)$ in Fig.~\ref{fig:orderinguncertainty}c), then $I[V_{\alpha=0.4}(f(\mbf{Z}))]$ is reduced to $[l(\mbf{z}_1), u(\mbf{z}_0)]$
because now $\mbf{z}_u=\mbf{z}_0$.
Thus, as the confidence bound $I[f(\mbf{z}_0)]$ is shortened with more and more observations at $\mbf{z}_0$, the confidence bound $I[V_{\alpha=0.4}(f(\mbf{Z}))]$ reduces  (the white region in Fig.~\ref{fig:orderinguncertainty} is `laced up').
\end{example}
\begin{figure*}[h!]
\centering
\begin{tabular}{@{}ccc@{}}
\includegraphics[width=0.25\textwidth]{img/uncertainty/ordering-uncertainty-top}
&
\includegraphics[width=0.25\textwidth]{img/uncertainty/ordering-uncertainty-bottom}
&
\includegraphics[width=0.25\textwidth]{img/uncertainty/ordering-uncertainty-extra}
\\
(a) & (b) & (c)
\end{tabular}
\caption{
A counterexample against selecting $\mbf{z}_u$ and $\mbf{z}_l$ as input queries.
Here $\mbf{z}$ follows a discrete uniform distribution over $\mcl{D}_{\mbf{z}} \triangleq \{\mbf{z}_0, \mbf{z}_1, \mbf{z}_2\}$.
(a) shows the mappings of $\mbf{z}$ to the upper bound $u(\mbf{z})$ and lower bound $l(\mbf{z})$.
The {\textsc{VaR}} at $\alpha=0.4$ of $u(\mbf{Z})$ and $l(\mbf{Z})$ are $u(\mbf{z}_2)$ and $l(\mbf{z}_1)$, respectively, i.e., $\mbf{z}_u = \mbf{z}_2$ and $\mbf{z}_l = \mbf{z}_1$.
(b) shows the values of $l(\mbf{z})$ and $u(\mbf{z})$ for all $\mbf{z}$ on the same axis, as well as the confidence bounds of $f(\mbf{z})$ and $V_{\alpha}(f(\mbf{Z}))$. The gray areas A and B indicate the intervals of values $\omega \in \mbb{R}$ where $\omega \le l(\mbf{z}_l) = l(\mbf{z}_1)$ and $\omega \ge u(\mbf{z}_u) = u(\mbf{z}_2)$, respectively.
(c) shows a hypothetical scenario when $I[f(\mbf{z}_0)]$ is shortened with more observations at $\mbf{z}_0$.
}
\label{fig:orderinguncertainty}
\end{figure*}

In the next section, we define a property of $\mbf{z}_0$ in Example~\ref{example:counter}
and prove the existence of $\mbf{z}$'s with this property.
Then, we prove that along with the optimistic selection of $\mbf{x}_t$, the selection of $\mbf{z}_t$ such that it satisfies this property leads to a no-regret algorithm.

\subsection{Lacing Value (LV) and the Query Selection Strategy for $\mbf{z}_t$}
\label{subsec:lv}

We note that in Example~\ref{example:counter}, as long as the confidence bound of the function evaluation at $\mbf{z}_0$ contains
the confidence bound
of {\textsc{VaR}}, observations at $\mbf{z}_0$ can
reduce the confidence bound of {\textsc{VaR}}. We name the values of $\mbf{z}$ satisfying this property as \emph{lacing values} (LV):

\begin{definition}[Lacing values]
\label{definition:lv}
\emph{Lacing values} (LV) with respect to $\mbf{x} \in \mcl{D}_{\mbf{x}}$ and $t \ge 1$ are $\mbf{z}_{\text{LV}} \in \mcl{D}_{\mbf{z}}$ that satisfies $l_{t-1}(\mbf{x},\mbf{z}_{\text{LV}}) \le V_{\alpha}(l_{t-1}(\mbf{x}, \mbf{Z})) \le V_{\alpha}(u_{t-1}(\mbf{x}, \mbf{Z})) \le u_{t-1}(\mbf{x},\mbf{z}_{\text{LV}})$, equivalently,
$I_{t-1}[V_{\alpha}(f(\mbf{x},\mbf{Z}))] \subset [l_{t-1}(\mbf{x},\mbf{z}_{\text{LV}}), u_{t-1}(\mbf{x},\mbf{z}_{\text{LV}})]\ .$
\end{definition}

Recall that the support $\mcl{D}_{\mbf{z}}$ of $\mbf{Z}$ is defined as the smallest closed subset $\mcl{D}_{\mbf{z}}$ of $\mbb{R}^{d_z}$ such that $P(\mbf{Z} \in \mcl{D}_{\mbf{z}}) = 1$ (e.g., $\mcl{D}_{\mbf{z}}$ is a finite subset of $\mbb{R}^{d_z}$ and $\mcl{D}_{\mbf{z}} = \mbb{R}^{d_z}$).
The following theorem guarantees the existence of lacing values and is proved in Appendix~\ref{app:prooflv}.

\begin{theorem}[Existence of lacing values]
\label{theorem:lv}
$\forall \alpha \in (0,1)$, $\forall \mbf{x} \in \mcl{D}_{\mbf{x}}$, $\forall t \ge 1$, there exists a lacing value in $\mcl{D}_{\mbf{z}}$ with respect to $\mbf{x}$ and $t$.
\end{theorem}

\begin{corollary}
\label{corollary:loclv}
Lacing values with respect to $\mbf{x} \in \mcl{D}_{\mbf{x}}$ and $t \ge 1$ are in $\mcl{Z}_l^{\le} \cap \mcl{Z}_u^{\ge}$ where $\mcl{Z}_l^{\le} \triangleq \{\mbf{z} \in \mcl{D}_{\mbf{z}}: l_{t-1}(\mbf{x},\mbf{z}) \le V_{\alpha}(l_{t-1}(\mbf{x}, \mbf{Z}))\}$ and $\mcl{Z}_u^{\ge} \triangleq \{\mbf{z} \in \mcl{D}_{\mbf{z}}: u_{t-1}(\mbf{x},\mbf{z}) \ge V_{\alpha}(u_{t-1}(\mbf{x}, \mbf{Z}))\}$.
\end{corollary}

As a special case, when $\mbf{z}_l = \mbf{z}_u$,  $I_{t-1}[V_\alpha(f(\mbf{x},\mbf{Z}))] = I_{t-1}[f(\mbf{x},\mbf{z}_l)]$ which means $\mbf{z}_l = \mbf{z}_u$ is an LV.
Based on Theorem~\ref{theorem:lv}, we can always select $\mbf{z}_t$ as an LV w.r.t $\mbf{x}_t$ defined in Definition~\ref{definition:lv}.
This strategy for the selection of $\mbf{z}_t$, together with the selection of $\mbf{x}_t = \text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$ (Section \ref{subsec:confboundf}), completes our algorithm: \emph{{\textsc{VaR}} Upper Confidence Bound} ({V-UCB}) (Algorithm~\ref{alg:v-ucb}).

\textbf{Upper Bound on Regret.}
As a result of the selection strategies for $\mbf{x}_t$ and $\mbf{z}_t$, our V-UCB algorithm enjoys the following upper bound on its instantaneous regret (proven in Appendix~\ref{app:iregretbound2}):
\begin{lemma}
\label{lemma:iregretbound2}
By selecting $\mbf{x}_t$ as a maximizer of $V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$ and selecting $\mbf{z}_t$ as an LV w.r.t $\mbf{x}_t$, the instantaneous regret is upper-bounded by:
\begin{equation*}
r(\mbf{x}_t) \le 2 \beta_t^{1/2} \sigma_{t-1}(\mbf{x}_t, \mbf{z}_t)\ \forall t \ge 1
\end{equation*}
with probability $\ge 1 - \delta$ for $\delta$ in Lemma~\ref{lemma:ucb51}.
\end{lemma}
Lemma~\ref{lemma:iregretbound2}, together with Lemma 5.4 from \citet{srinivas10ucb}, implies that the cumulative regret of our algorithm is bounded (Appendix~\ref{app:rtbound}): $R_T \le \sqrt{C_1 T \beta_T \gamma_T}$
where $C_1 \triangleq 8/\log(1 + \sigma_n^{-2})$, and $\gamma_T$ is the maximum information gain about $f$ that can be obtained from any set of $T$ observations.
\citet{srinivas10ucb} have analyzed $\gamma_T$ for several commonly used kernels such as SE and Mat\'ern kernels, and have shown that for these kernels, the upper bound on $R_T$ grows sub-linearly.
This implies that our algorithm is \emph{no-regret} because $\lim_{T\rightarrow \infty} R_T / T = 0$.

\begin{algorithm}[tb]
   \caption{The {V-UCB} Algorithm}
\begin{algorithmic}[1]
   \STATE {\bfseries Input:} $\mcl{D}_{\mbf{x}}$, $\mcl{D}_{\mbf{z}}$, prior $\mu_0=0, \sigma_0, k$
   \FOR{$i=1,2,\dots$}
   \STATE Select $\mbf{x}_t = \text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$
   \STATE Select $\mbf{z}_t$ as a \emph{lacing value} w.r.t. $\mbf{x}_t$ (Definition~\ref{definition:lv})
   \STATE Obtain observation $y_t \triangleq f(\mbf{x}_t,\mbf{z}_t) + \epsilon_t$
   \STATE Update the GP posterior belief to obtain $\mu_t$ and $\sigma_t$
   \ENDFOR
\end{algorithmic}
\label{alg:v-ucb}
\end{algorithm}

Inspired by \citet{bogunovic2018adversarially}, at the {$T$-th} iteration of {V-UCB}, we can recommend $\mbf{x}_{t_*(T)}$ as an estimate of the maximizer $\mbf{x}_*$ of $V_\alpha(f(\mbf{x},\mbf{Z}))$, where $t_*(T) \triangleq  \text{arg}\!\max_{t \in \{1,\dots,T\}} V_{\alpha}(l_{t-1}(\mbf{x}_t,\mbf{Z}))$. Then, the instantaneous regret $r(\mbf{x}_{t_*(T)})$ is upper-bounded by $\sqrt{C_1 \beta_T \gamma_T / T}$ with high probability as we show in Appendix~\ref{app:recommendxbound}.
In our experiments in Section~\ref{sec:experiments},
we recommend $\text{arg}\!\max_{\mbf{x} \in \mcl{D}_T} V_{\alpha}(\mu_{t-1}(\mbf{x}, \mbf{Z}))$ (where $\mu_{t-1}(\mbf{x}, \mbf{Z})$ is a random function defined in the same manner as $f(\mbf{x},\mbf{Z})$) as an estimate of $\mbf{x}_*$ due to its empirical convergence.

\textbf{Computational Complexity.}
To compare our computational complexity
with that of the {$\rho \text{KG}$} algorithm from \citet{borisk20}, we exclude the common part of training the GP model (line 6) and assume that $\mcl{D}_{\mbf{z}}$ is finite.
Then, the time complexity of {V-UCB} is dominated by that of the selection of $\mbf{x}_t$ (line 3) which includes the time complexity $\mcl{O}(|\mcl{D}_{\mbf{z}}| |\mcl{D}_{t-1}|^2)$ for the GP prediction at $\{\mbf{x}\} \times \mcl{D}_{\mbf{z}}$, and $\mcl{O}(|\mcl{D}_{\mbf{z}}|\log|\mcl{D}_{\mbf{z}}|)$ for the sorting of $u_{t-1}(\mbf{x},\mcl{D}_{\mbf{z}})$ and searching of {\textsc{VaR}}.
Hence, our overall complexity is $\mcl{O}(n |\mcl{D}_{\mbf{z}}|\ (|\mcl{D}_{t-1}|^2 + \log |\mcl{D}_{\mbf{z}}|))$, where $n$ is the number of iterations to maximize $V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$ (line 3).
Therefore, our V-UCB is more computationally efficient than {$\rho \text{KG}$} and its variant with approximation {$\rho \text{KG}^{apx}$}, whose complexities are $\mcl{O}(n_{\text{out}} n_{\text{in}} K |\mcl{D}_{\mbf{z}}|\ (|\mcl{D}_{t-1}|^2 + |\mcl{D}_{\mbf{z}}| |\mcl{D}_{t-1}| + |\mcl{D}_{\mbf{z}}|^2 + M |\mcl{D}_{\mbf{z}}|))$ of {$\rho \text{KG}$} and $\mcl{O}(n_{\text{out}} |\mcl{D}_{t-1}| K |\mcl{D}_{\mbf{z}}|\ (|\mcl{D}_{t-1}|^2 + |\mcl{D}_{\mbf{z}}| |\mcl{D}_{t-1}| + |\mcl{D}_{\mbf{z}}|^2 + M |\mcl{D}_{\mbf{z}}|))$, respectively.\footnote{$n_{\text{out}}$ and $n_{\text{in}}$ are the numbers of iterations for the outer and inner optimization respectively, $K$ is the number of fantasy GP models required for their one-step lookahead, and $M$ is the number of functions sampled from the GP posterior \cite{borisk20}.}

\subsection{On the Selection of $\mbf{z}_t$}
\label{subsec:improvedselectz}

Although Algorithm~\ref{alg:v-ucb} is guaranteed to be no-regret with any choice of LV as $\mbf{z}_t$, we would like to select the LV that can reduce a large amount of the uncertainty of $V_{\alpha}(f(\mbf{x}_t,\mbf{Z}))$.
However, relying on the information gain measure or the knowledge gradient method often incurs the expensive one-step lookahead.
Therefore, we use a simple heuristic by choosing the LV $\mbf{z}_{\text{LV}}$ with the maximum probability mass (or probability density if $\mbf{Z}$ is continuous) of $\mbf{z}_{\text{LV}}$.
We motivate this heuristic using an example with $\alpha = 0.2$, i.e., $V_{\alpha=0.2}(f(\mbf{x}_t,\mbf{Z})) = \inf \{\omega: P(f(\mbf{x}_t,\mbf{Z}) \le \omega) \ge 0.2\}$.
Suppose $\mcl{D}_{\mbf{z}}$ is finite and there are $2$ LV's $\mbf{z}_0$ and $\mbf{z}_1$ with $P(\mbf{z}_0) \ge 0.2$ and $P(\mbf{z}_1) = 0.01$.
Then, because $P(f(\mbf{x}_t,\mbf{Z}) \le f(\mbf{x}_t,\mbf{z}_0)) \ge P(\mbf{z}_0) \ge 0.2$, it follows that $V_{\alpha=0.2}(f(\mbf{x}_t,\mbf{Z})) \le f(\mbf{x}_t,\mbf{z}_0)$, i.e., querying $\mbf{z}_0$ at $\mbf{x}_t$ gives us information about an explicit upper bound on $V_{\alpha=0.2}(f(\mbf{x}_t,\mbf{Z}))$ to reduce its uncertainty.
In contrast, this cannot be achieved by querying $\mbf{z}_1$ due to its low probability mass.
This simple heuristic can also be implemented when $\mbf{Z}$ is a continuous random variable which we will introduce in Section~\ref{sec:continuousz}.
\begin{remark}
\label{rmk:unknownZ}
Although we assume that we can select both $\mbf{x}_t$ and $\mbf{z}_t$ during our
algorithm, Corollary~\ref{corollary:loclv} also gives us some insights about the scenario where we cannot select $\mbf{z}_t$.
In this case, in each iteration $t$, we select $\mbf{x}_t$ while $\mbf{z}_t$ is randomly sampled by the environment following the distribution of $\mbf{Z}$.
Next, we observe both $\mbf{z}_t$ and $y_{(\mbf{x}_t,\mbf{z}_t)}$ and then update the GP posterior belief of $f$.
Of note, Corollary~\ref{corollary:loclv} has shown that all LV lie in the set $\mcl{Z}_l^{\le} \cap \mcl{Z}_u^{\ge}$.
However, the probability of this set is usually small, because $P(\mbf{Z} \in \mcl{Z}_l^{\le} \cap \mcl{Z}_u^{\ge}) \le P(\mbf{Z} \in \mcl{Z}_l^{\le}) \le \alpha$ and small values of $\alpha$ are often used by real-world applications to manage risks.
Thus, the probability that the sampled $\mbf{z}_t$ is an LV is small.
As a result, we suggest sampling a large number of $\mbf{z}_t$'s from the environment to increase the chance that an LV is sampled.
On the other hand, the small probability of sampling an LV
motivates the need for us to select $\mbf{z}_t$.
\end{remark}

\subsection{{V-UCB} Approaches {\textsc{StableOpt}} as $\alpha \rightarrow 0^+$}
\label{subsec:vy-stableopt}

Recall that the objective of adversarially robust BO is to find $\mbf{x} \in \mcl{D}_{\mbf{x}}$ that maximizes
$\min_{\mbf{z} \in \mcl{D}_{\mbf{z}}} f(\mbf{x}, \mbf{z})$ \cite{bogunovic2018adversarially} by iteratively specifying input query $(\mbf{x}_t, \mbf{z}_t)$ to collect noisy observations $y_{\mbf{x}_t,\mbf{z}_t}$.
It is different from BO of {\textsc{VaR}} since its $\mbf{z}$ is not random but selected by an adversary who aims to minimize the function evaluation.
The work of \citet{bogunovic2018adversarially} has proposed a no-regret algorithm for this setting named {\textsc{StableOpt}}, which selects
\begin{equation}
\begin{array}{r@{}l}
\mbf{x}_t &= \text{arg}\!\max_{\mbf{x} \in \mcl{D}_x} \min_{\mbf{z} \in \mcl{D}_z} u_{t-1}(\mbf{x},\mbf{z})\ ,\\
\mbf{z}_t &= \text{arg}\!\min_{\mbf{z} \in \mcl{D}_z} l_{t-1}(\mbf{x}_t,\mbf{z})
\end{array}
\label{eq:stableopt}
\end{equation}
where $u_{t-1}$ and $l_{t-1}$ are defined in \eqref{eq:fbound}.

At first glance, BO of {\textsc{VaR}} and adversarially robust BO are seemingly different problems because $\mbf{Z}$ is a random variable in the former but not in the latter.
However, based on our key observation on the connection between the minimum value of a continuous function $h(\mbf{w})$ and the {\textsc{VaR}} of the random variable $h(\mbf{W})$ in the following theorem, these two problems and their solutions are connected as shown in Corollary~\ref{corollary:alpha0lv}, and~\ref{corollary:stableopt}.

\begin{theorem}
\label{theorem:0plus}
Let $\mbf{W}$ be a random variable with the support $\mcl{D}_w \subset \mbb{R}^{d_w}$ and dimension $d_w$.
Let $h$ be a continuous function mapping from $\mbf{w} \in \mcl{D}_w$ to $\mbb{R}$. Then, $h(\mbf{W})$ denotes the random variable whose realization is the function $h$ evaluation at a realization $\mbf{w}$ of $\mbf{W}$. Suppose $h(\mbf{w})$ has a minimizer $\mbf{w}_{\min} \in \mcl{D}_w$, then $\lim_{\alpha \rightarrow 0^+} V_{\alpha}(h(\mbf{W})) = h(\mbf{w}_{\min})\ .$
\end{theorem}

\begin{corollary}
\label{corollary:alpha0lv}
Adversarially robust BO which finds $\text{arg}\!\max_{\mbf{x}} \min_{\mbf{z}} f(\mbf{x},\mbf{z})$ can be cast as BO of {\textsc{VaR}} by letting (a) $\alpha$ approach $0$ from the right and (b) $\mcl{D}_{\mbf{z}}$ be the support of the environmental random variable $\mbf{Z}$, i.e., $\text{arg}\!\max_{\mbf{x}} \lim_{\alpha \rightarrow 0^+} V_\alpha(f(\mbf{x},\mbf{Z}))$.
\end{corollary}

Interestingly, from Theorem~\ref{theorem:0plus}, we observe that $\mcl{Z}_l^\le$ in Corollary~\ref{corollary:loclv}
approaches
the set of minimizers  $\text{arg}\!\min_{\mbf{z} \in \mcl{D}_{\mbf{z}}} l_{t-1}(\mbf{x}_t,\mbf{z})$ as $\alpha \rightarrow 0^+$.
Corollary \ref{corollary:stableopt} below shows that LV w.r.t $\mbf{x}_t$ becomes a minimizer of $l_{t-1}(\mbf{x}_t,\mbf{z})$ which is same as the selected $\mbf{z}_t$ of {\textsc{StableOpt}} in \eqref{eq:stableopt}.

\begin{corollary}
\label{corollary:stableopt}
The {\textsc{StableOpt}} solution to adversarially robust BO selects the same input query as that selected by the {V-UCB} solution to the corresponding BO of {\textsc{VaR}} in Corollary~\ref{corollary:alpha0lv}.
\end{corollary}

The proof of Theorem~\ref{theorem:0plus} and its two corollaries are shown in Appendix~\ref{app:0plus}.
We note that {V-UCB} is also applicable to the optimization of $V_{\alpha}(f(\mbf{x},\mbf{Z}))$ where the distribution of $\mbf{Z}$ is a conditional distribution given $\mbf{x}$. For example,
in robotics, if there exists noise/perturbation in the control, an optimization problem of interest is
$V_{\alpha}(f(\mbf{x} + \bm{\xi}(\mbf{x})))$ where $\bm{\xi}(\mbf{x})$ is the random perturbation that depends on $\mbf{x}$.

\subsection{Implementation of V-UCB with Continuous Random Variable $\mbf{Z}$}
\label{sec:continuousz}

The V-UCB algorithm involves two steps: selecting $\mbf{x}_t = \text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ and selecting $\mbf{z}_t$ as the LV $\mbf{z}_{\text{LV}}$ with the largest probability mass (or probability density). When $|\mcl{D}_{\mbf{z}}|$ is finite, given $\mbf{x}$,  $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ can be computed exactly.
The gradient of $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ with respect to $\mbf{x}$ can be obtained
easily (e.g., using
automatic differentiation provided in the Tensorflow library \cite{tensorflow2015-whitepaper})
to aid the selection of $\mbf{x}_t$.
In this case, the latter step can also be performed by constructing the set of all LV (checking the condition in the Definition~\ref{definition:lv} for all $\mbf{z} \in \mcl{D}_{\mbf{z}}$) and choosing the LV $\mbf{z}_{\text{LV}}$ with the largest probability mass.

{\bf{Estimation of {\textsc{VaR}}.}} When $\mbf{Z}$ is a continuous random variable, estimating {\textsc{VaR}} by an ordered set of samples (e.g., in \citet{borisk20}) may require a prohibitively large number of samples, especially for small values of $\alpha$. Thus, we employ the following popular pinball (or tilted absolute value) function in quantile regression \cite{koenker1978regression} to estimate {\textsc{VaR}} as a lower $\alpha$-quantile:
\begin{equation*}
\rho_{\alpha}(w) \triangleq \begin{cases}
    \alpha w &\text{if } w \ge 0\ ,\\
    (\alpha - 1) w &\text{if } w < 0
\end{cases}
\end{equation*}
where $w \in \mbb{R}$. In particular, to estimate $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ as $\nu \in \mbb{R}$, we find $\nu$ that minimizes:
\begin{equation}
    \mbb{E}_{\mbf{z} \sim p(\mbf{Z})}[\rho_{\alpha}(u_{t-1}(\mbf{x},\mbf{z}) - \nu)]\ .
    \label{eq:lossfindvar}
\end{equation}
A well-known example is when $\alpha = 0.5$ and $\rho_{\alpha}$ is the absolute value function, then the optimal $\nu$ is the median.
The loss in \eqref{eq:lossfindvar} can be optimized using stochastic gradient descent with a random batch of samples of $\mbf{Z}$ at each optimization iteration.

{\bf{Maximization of $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$.}} Unfortunately, to maximize $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ over $\mbf{x} \in \mcl{D}_{\mbf{x}}$, there is no gradient of $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ with respect to $\mbf{x}$ under the above approach. This situation resembles BO where there is no gradient information, but only noisy observations at input queries. Unlike BO, the observation (samples of $u_{t-1}(\mbf{x},\mbf{Z})$ at $\mbf{x}$) is not costly. Therefore, we propose the \emph{local neural surrogate optimization} (LNSO) algorithm to find $\text{arg}\!\max_{\mbf{x} \in \mcl{D}_{\mbf{x}}} V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ which is visualized in Fig.~\ref{fig:lnso}. Suppose the optimization is initialized at $\mbf{x} = \mbf{x}^{(0)}$, instead of maximizing $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ (whose gradient w.r.t. $\mbf{x}$ is unknown), we maximize a surrogate function $g(\mbf{x},\bm{\theta}^{(0)})$ (modeled by a neural network) that approximates $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ well in a local region of $\mbf{x}^{(0)}$, e.g., a ball $\mcl{B}(\mbf{x}^{(0)}, r)$ of radius $r$ in Fig.~\ref{fig:lnso}.
We obtain such parameters $\bm{\theta}^{(0)}$ by minimizing the following loss function:
\begin{equation}
\begin{array}{l}
\mcl{L}_g(\bm{\theta},\mbf{x}^{(0)})\\
\ \displaystyle\triangleq \mbb{E}_{\mbf{x} \in \mcl{B}(\mbf{x}^{(0)}, r)} \mbb{E}_{\mbf{Z} \sim p(\mbf{Z})} \left[ \rho_\alpha(u_{t-1}(\mbf{x},\mbf{z}) - g(\mbf{x};\bm{\theta}))\right]
\end{array}
\label{eq:losslnso}
\end{equation}
where the expectation $\mbb{E}_{\mbf{x} \in \mcl{B}(\mbf{x}^{(0)}, r)}$ is taken over a uniformly distributed $\mbf{X}$ in $\mcl{B}(\mbf{x}^{(0)}, r)$. Minimizing \eqref{eq:losslnso} can be performed with stochastic gradient descent.
If maximizing $g(\mbf{x},\bm{\theta}^{(0)})$ leads to a value $\mbf{x}^{(i)} \notin \mcl{B}(\mbf{x}^{(0)}, r)$ (Fig.~\ref{fig:lnso}), we update the local region to be centered at $\mbf{x}^{(i)}$ ($\mcl{B}(\mbf{x}^{(i)}, r)$) and find  $\bm{\theta}^{(i)} = \text{arg}\!\min_{\bm{\theta}} \mcl{L}_g(\bm{\theta}, \mbf{x}^{(i)})$ such that $g(\mbf{x},\bm{\theta}^{(i)})$ approximates $V_{\alpha}(u_{t-1}(\mbf{x},\mbf{Z}))$ well $\forall \mbf{x} \in \mcl{B}(\mbf{x}^{(i)},r)$. Then, $\mbf{x}^{(i)}$ is updated by maximizing $g(\mbf{x},\bm{\theta}^{(i)})$ for $\mbf{x} \in \mcl{B}(\mbf{x}^{(i)},r)$. The complete algorithm is described in Appendix~\ref{app:lnso}.

We prefer a small value of $r$ so that the ball $\mcl{B}$ is small. In such case, $V_\alpha(u_{t-1}(\mbf{x},\mbf{Z}))$ for $\mbf{x} \in \mcl{B}$ can be estimated well with a small neural network $g(\mbf{x},\bm{\theta})$ whose training requires a small number of iterations.

\begin{figure}
    \centering
    \includegraphics[width=0.27\textwidth]{img/lnso.pdf}
    \caption{Plot of a hypothetical optimization path (as arrows) of LNSO initialized at $\mbf{x}^{(0)}$. Input $\mbf{x}$ is $2$-dimensional. The boundary of a ball $\mcl{B}$ of radius $r$ is plotted as a dotted circle. When the updated $\mbf{x}$ moves out of $\mcl{B}$, the center of $\mcl{B}$ and $\bm{\theta}$ are updated.}
    \label{fig:lnso}
\end{figure}

{\bf{Search of Lacing Values.}}
Given a continuous random variable $\mbf{Z}$, to find an LV w.r.t $\mbf{x}_t$ in line 4 of Algorithm~\ref{alg:v-ucb}, i.e., to find a $\mbf{z}$ satisfying $d_u(\mbf{z}) \triangleq u_{t-1}(\mbf{x}_t, \mbf{z}) - V_\alpha(u_{t-1}(\mbf{x}_t,\mbf{Z})) \ge 0$ and $d_l(\mbf{z}) \triangleq V_\alpha(l_t(\mbf{x}_t,\mbf{Z})) - l_{t-1}(\mbf{x}_t, \mbf{z}) \ge 0$, we choose a $\mbf{z}$ that minimizes
\begin{equation}
\mcl{L}_{\text{LV}}(\mbf{z}) \triangleq \text{ReLU}(-d_u(\mbf{z})) + \text{ReLU}(-d_l(\mbf{z}))
\label{eq:llv}
\end{equation}
where $\text{ReLU}(\omega) = \max(\omega,0)$ is the rectified linear unit function ($\omega \in \mbb{R}$). To include the heuristic in Section~\ref{subsec:improvedselectz} which selects the LV with the highest probability density, we find $\mbf{z}$ that minimizes
\begin{equation*}
\mcl{L}_{\text{LV-P}}(\mbf{z}) \triangleq \mcl{L}_{\text{LV}}(\mbf{z}) - \mbb{I}_{d_u(\mbf{z}) \ge 0} \mbb{I}_{ d_l(\mbf{z}) \ge 0}\ p(\mbf{z})
\end{equation*}
where $\mcl{L}_{\text{LV}}(\mbf{z})$ is defined in \eqref{eq:llv}; $p(\mbf{z})$ is the probability density; $\mbb{I}_{d_u(\mbf{z}) \ge 0}$ and $\mbb{I}_{ d_l(\mbf{z}) \ge 0}$ are indicator functions.

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

In this section, we empirically evaluate the performance of {V-UCB}. The work of \citet{borisk20} has motivated the use of the approximated variant of their algorithm {$\rho \text{KG}^{apx}$} over its original version {$\rho \text{KG}$} by showing that {$\rho \text{KG}^{apx}$} achieves comparable empirical performances to {$\rho \text{KG}$} while incurring much less computational cost. Furthermore, {$\rho \text{KG}^{apx}$} has been shown to significantly outperform
other competing algorithms \cite{borisk20}.
Therefore, we choose {$\rho \text{KG}^{apx}$} as the main baseline to empirically compare with {V-UCB}. The experiments using {$\rho \text{KG}^{apx}$} is performed by
adding
new objective functions to the existing implementation of \citet{borisk20} at \texttt{https://github.com/saitcakmak/BoRisk}.

Regarding {V-UCB}, when $\mcl{D}_{\mbf{z}}$ is finite and the distribution of $\mbf{Z}$ is not uniform, we perform {V-UCB} by selecting $\mbf{z}_t$ as an LV at random, labeled as \emph{{V-UCB} Unif}, and by selecting $\mbf{z}_t$ as the LV with the maximum probability mass, labeled as \emph{{V-UCB} Prob}.

The performance metric
is defined as $V_\alpha(f(\mbf{x}_*,\mbf{Z})) - V_\alpha(f(\tilde{\mbf{x}},\mbf{Z}))$ where $\tilde{\mbf{x}}$ is the recommended input.
The evaluation of {\textsc{VaR}} is described in Section~\ref{sec:continuousz}.
The recommended input is $\text{arg}\!\max_{\mbf{x} \in \mcl{D}_T} V_{\alpha}(\mu_{t-1}(\mbf{x}, \mbf{Z}))$ for {V-UCB}, and $\text{arg}\!\min_{\mbf{x} \in \mcl{D}_{\mbf{x}}} \mbb{E}_{t-1}[V_{\alpha}(f(\mbf{x}, \mbf{Z}))]$ for {$\rho \text{KG}^{apx}$} \cite{borisk20}, where $\mbb{E}_{t-1}$ is the conditional expectation over the unknown $f$ given the observations $\mcl{D}_{t-1}$ (approximated by a finite set of functions sampled from the GP posterior belief).\footnote{While the work of \citet{borisk20} considers a minimization problem of {\textsc{VaR}}, our work considers a maximization problem of {\textsc{VaR}}. Therefore, the objective functions for {$\rho \text{KG}^{apx}$} are the negation of those for {V-UCB}. For {V-UCB} at risk level $\alpha$,
the risk level for {$\rho \text{KG}^{apx}$} is $1-\alpha$.}
We repeat each experiment $10$ times with different random initial observations $\mbf{y}_{\mcl{D}_0}$ and plot both the mean (as lines) and the $70\%$ confidence interval (as shaded areas) of the $\log10$ of the performance metric.
The detailed descriptions of experiments are deferred to Appendix~\ref{app:experiment}.

\subsection{Synthetic Benchmark Functions}
\label{subsec:syn}

We use $3$ synthetic benchmark functions: Branin-Hoo, Goldstein-Price, and Hartmann-3D functions to construct $4$ optimization problems: Branin-Hoo-$(1,1)$, Goldstein-Price-$(1,1)$, Hartmann-$(1,2)$, and Hartmann-$(2,1)$.
The tuples represent $(d_x, d_z)$ corresponding to the
dimensions of $\mbf{x}$ and $\mbf{z}$.
The noise variance $\sigma_n^2$ is set to $0.01$. The risk level $\alpha$ is $0.1$.
There are $2$ different settings: finite $\mcl{D}_{\mbf{z}}$ ($|\mcl{D}_{\mbf{z}}| = 64$ for Hartmann-$(1,2)$ and $|\mcl{D}_{\mbf{z}}|=100$ for the others) and continuous $\mcl{D}_{\mbf{z}}$. In the latter setting, $r=0.1$ and the surrogate function is a neural network of $2$ hidden layers with $30$ hidden neurons each, and sigmoid activation functions.

The results are shown in Fig.~\ref{fig:synfinite} and Fig.~\ref{fig:syncont} for the settings of discrete $\mcl{D}_{\mbf{z}}$ and continuous $\mcl{D}_{\mbf{z}}$, respectively. When $\mcl{D}_{\mbf{z}}$ is discrete (Fig.~\ref{fig:synfinite}), {V-UCB} Unif is on par with {$\rho \text{KG}^{apx}$} in optimizing Branin-Hoo-$(1,1)$ and Goldstein-Price-$(1,1)$, and outperforms {$\rho \text{KG}^{apx}$} in optimizing Hartmann-$(1,2)$ and Hartmann-$(2,1)$. {V-UCB} Prob is also able to exploit the probability distribution of $\mbf{Z}$ to outperform {V-UCB} Unif. When $\mcl{D}_{\mbf{z}}$ is continuous (Fig.~\ref{fig:syncont}), {V-UCB} Prob outperforms {$\rho \text{KG}^{apx}$}.
The unsatisfactory performance of {$\rho \text{KG}^{apx}$} in some experiments may be attributed to its approximation of the inner optimization problem in the acquisition function \cite{borisk20}, and the approximation of {\textsc{VaR}} using samples of $\mbf{Z}$ and the GP posterior belief.
\newcommand\figheight{0.201}
\begin{figure}[ht]
\centering
\begin{tabular}{@{}cc@{}}
\includegraphics[height=\figheight\textwidth]{img/negative_branin_discrete}
&
\includegraphics[height=\figheight\textwidth]{img/negative_goldstein_discrete}
\\
(a) Branin-Hoo-$(1,1)$
&
(b) Goldstein-Price-$(1,1)$
\\
\\
\includegraphics[height=\figheight\textwidth]{img/negative_hartmann3d_12_discrete}
&
\includegraphics[height=\figheight\textwidth]{img/negative_hartmann3d_21_discrete}
\\
(c) Hartmann-$(1,2)$
&
(d) Hartmann-$(2,1)$
\end{tabular}
\caption{Synthetic benchmark functions with finite $\mcl{D}_{\mbf{z}}$.}
\label{fig:synfinite}
\end{figure}
\begin{figure}[ht]
\centering
\begin{tabular}{@{}cc@{}}
\includegraphics[height=\figheight\textwidth]{img/negative_branin_continuous.pdf}
&
\includegraphics[height=\figheight\textwidth]{img/negative_goldstein_continuous.pdf}
\\
(a) Branin-Hoo-$(1,1)$
&
(b) Goldstein-Price-$(1,1)$
\\
\\
\includegraphics[height=\figheight\textwidth]{img/negative_hartmann3d_12_continuous.pdf}
&
\includegraphics[height=\figheight\textwidth]{img/negative_hartmann3d_21_continuous.pdf}
\\
(c) Hartmann-$(1,2)$
&
(d) Hartmann-$(2,1)$
\end{tabular}
\caption{Synthetic benchmark functions with continuous $\mcl{D}_{\mbf{z}}$.}
\label{fig:syncont}
\end{figure}
\subsection{Simulated Optimization Problems}
\label{subsec:simulated}

The first problem is portfolio optimization adopted by \cite{borisk20}.
There are $d_x = 3$ optimization variables (risk and trade aversion parameters, and the holding cost multiplier) and $d_z = 2$ environmental random variables (bid-ask spread and borrow cost). The variable $\mbf{Z}$ follows a discrete uniform distribution with $|\mcl{D}_{\mbf{z}}| = 100$. Hence, there is no difference between {V-UCB} Unif and {V-UCB} Prob. Thus, we only report the results of the latter.
The objective function is the posterior mean of a trained GP on the dataset in \citet{borisk20} of size $3000$ generated from CVXPortfolio. The noise variance $\sigma_n^2$ is set to $0.01$. The risk level $\alpha$ is set to $0.2$.

The second problem is a simulated robot pushing task for which we use the implementation from the work of \citet{wang17mes}.
The simulation is viewed as a $3$-dimensional function $\mbf{h}(r_x, r_y, t_p)$ returning the 2-D location of the pushed object, where $r_x,r_y \in [-5,5]$ are the robot location and $t_p \in [1,30]$ is the pushing duration.
The objective is to minimize the distance to a fixed goal location $\mbf{g}=(g_x,g_y)^\top$, i.e., the objective function of the maximization problem is $f_0(r_x, r_y, t_p) = -\Vert \mbf{h}(r_x, r_y, t_p) - \mbf{g}\Vert$. We assume that there are perturbations in specifying the robot location $W_x, W_y$ whose support $\mcl{D}_{\mbf{z}}$ includes $64$ equi-distant points in $[-1,1]^2$ and whose probability mass is proportional to $\exp(-(W_x^2 + W_y^2) / 0.4^2)$. It leads to a random objective function $f(r_x, r_y, t_p, W_x, W_y) \triangleq f_0(r_x + W_x, r_y + W_y, t_p)$.
We aim to maximize the {\textsc{VaR}} of $f$ which is more difficult than maximizing that of $f_0$.
Moreover, a random noise following $\mcl{N}(0,0.01)$ is added to the evaluation of $f$.
The risk level $\alpha$ is set to $0.1$.

The results are shown in Fig.~\ref{fig:portrobot}.
We observe that {V-UCB} outperforms {$\rho \text{KG}^{apx}$} in both problems. Furthermore, in comparison to our synthetic experiments, the difference between {V-UCB} Unif and {V-UCB} Prob is not significant in the robot pushing experiment. This is because the chance that a uniform sample of LV has a large probability mass is higher in the robot pushing experiment due to a larger region of $\mcl{D}_{\mbf{z}}$ having high probabilities.

\begin{figure}[ht]
\centering
\begin{tabular}{@{}cc@{}}
\includegraphics[height=\figheight\textwidth]{img/negative_portfolio_optimization_discrete}
&
\includegraphics[height=\figheight\textwidth]{img/robot_pushing_optimization_discrete}
\\
(a) Portfolio optimization $(3,2)$
&
(b) Robot pushing $(3,2)$
\end{tabular}
\caption{Simulated experiments.}
\label{fig:portrobot}
\end{figure}

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

To tackle the BO of {\textsc{VaR}} problem, we construct a no-regret algorithm, namely \emph{{\textsc{VaR}} upper confidence bound} (V-UCB), through the design of a confidence bound of {\textsc{VaR}} and a set of \emph{lacing values} (LV) that is guaranteed to exist. Besides, we introduce a heuristic to select an LV that improves the emprical performance of {V-UCB} over random selection of LV. We also draw
an elegant
connection between BO of {\textsc{VaR}} and adversarially robust BO in terms of both problem formulation and solutions. Lastly, we provide practical techniques for implementing {\textsc{VaR}} with continuous $\mbf{Z}$.
While {V-UCB} is more computationally efficient than the the state-of-the-art {$\rho \text{KG}^{apx}$} algorithm for BO of {\textsc{VaR}}, it also demonstrates competitive empirical performances in our experiments.


\section{Introduction}
Nowadays, the most powerful laser systems can reach intensity on the order $ 10^{22}~\mathrm{W/cm^{2}} $ and the dynamics of laser-plasma interaction becomes strongly relativistic \cite{Bahk2004,Yanovsky2008,Pirozhkov2017}.
One example is the onset of relativistic transparency in the interaction with the over-dense target \cite{Fernndez2017}.
The laser light cannot propagate in the over-dense target, as its frequency is lower or equal to the frequency of the electron plasma oscillations.
However, when the intense laser field is applied, the electron mass increases due to the relativistic motion reducing the plasma frequency.
It results in a decrease of the effective target density and thus the intense part of the laser field can propagate through the plasma.

Such an effect is especially important in the interaction with a thin over-dense foil, as it affects the temporal envelope of the laser pulse \cite{Vshivkov1998,Reed2009,Palaniyappan2012,Wei2017}.
Therefore, it has an impact on ion acceleration and radiation generation \cite{Henig2009,Kiefer2013}.
It was shown that in the case of an ultra-intense laser pulse, a relativistic plasma aperture is created in an ultra-thin foil and the laser pulse is diffracted \cite{GonzalezIzquierdo2016Optically}.
Since the plasma electrons are driven by the diffraction pattern of the laser field, the spatial structure of the accelerated electron beam can be controlled by varying the laser pulse parameters \cite{GonzalezIzquierdo2016Optically}.
As the electron dynamics consequently affects the spatial-intensity distribution of the accelerated ions, such an approach could be used for optically controlled ion acceleration \cite{GonzalezIzquierdo2016Towards,Higginson2018,Williamson2020}.

Here we study the increase of the laser field intensity caused by the interference of the light diffracted at the relativistic plasma aperture for the case of linearly polarized laser pulse.
It is investigated both analytically and using 2D and 3D Particle-In-Cell (PIC) simulations.
It is shown that due to the relativistic motion of the electrons, the laser intensity can be locally increased by almost an order of magnitude.
However, this effect is only accessible in the regime of relativistically induced transparency and thus it requires a careful choice of the laser and target parameters.
If the target is either under- or over-dense, the relativistic plasma aperture cannot be created and the diffraction structure responsible for the field intensity enhancement is not present.
The theoretical estimates for the maximum field strength and its spatial location are provided as a function of laser pulse and target parameters and compared with simulation results.

\section{Interaction with a plasma shutter}
In the following text, we study the process of laser intensity enhancement in the interaction of an ultra-intense laser pulse $( \sim10^{22}~\mathrm{W/cm^{2}} )$ with an ultra-thin $(\sim 10\mathrm{s}~\mathrm{nm})$ plasma layer having the electron plasma density $ n_{\mathrm{p}} $ greater than the value of the relativistic critical density $ n_{\mathrm{c}\gamma}=\gamma n_{\mathrm{c}} $ where $ \gamma $ is the relativistic Lorentz factor (of an electron having the charge $ e $ and mass $ m_{e} $) and $ n_\mathrm{c}= \omega_{0}^{2}m_{e}/\left(4\pi e^{2} \right) $ represents the non-relativistic critical plasma density.
We assume a laser pulse linearly polarized along the $ y $-axis propagating in the negative $ x $-direction focused to a spot of radius $ w_{0} $.
It has a Gaussian temporal envelope of full width at half maximum (FWHM) duration $ \tau $ in the laser intensity.
Its amplitude is characterized by the Lorentz invariant parameter $ a_{0}=eE_{0}/\left(m_{e}\omega_{0}c \right)  $ where $ E_{0} $ is the amplitude
of the electric field, $ \omega_{0} $ is the laser angular frequency and $ c $ is the speed of light \cite{Gibbon2005}.

When a relativistically over-dense shutter ($ n_{\mathrm{p}}>n_{\mathrm{c\gamma}} $) is assumed, the front part of the laser pulse is reflected from the foil surface.
However, as the laser pulse intensity grows to its maximum value, the laser pulse starts to penetrate the target.
The dominant mechanism allowing target penetration differs depending on the target thickness and density.
In the case of the ultra-thin foil, the electrons are pushed away due to the radiation pressure of the incoming laser pulse.
If the thicker target is assumed, it is rather the increasing relativistic mass of electrons that causes the foil to become relativistically transparent so that the rest of the laser pulse can pass through the foil.
Therefore, the laser field strength at which the plasma layer becomes transparent depends on the density and thickness of the foil as well as on the laser wavelength $ \lambda $.
In both cases, the laser pulse acquires a steep front edge as a result of the interaction \cite{Vshivkov1998,Reed2009,Palaniyappan2012}, that can be utilized, e.g., for enhanced photon emission due to the improved temporal contrast or for enhancement of ion acceleration from structured targets due to reduced development of transverse short-wavelength instabilities \cite{Jirka2020,Matys2020}.
Here we study the interaction with the ultra-thin foil.
For $ a_{0}\gg1 $  and an unperturbed plasma slab of thickness $ l \ll \lambda $, the laser pulse strength required for penetration of the foil (meaning that the laser pulse is able to push all electrons away from the initial position of the foil by the radiation pressure) can be estimated as \cite{Vshivkov1998,Bulanov2012,Bulanov2016}
\begin{equation}\label{a0reqMax}
a_{0}^{\mathrm{max}}=\dfrac{n_{\mathrm{p}}\pi l}{n_{\mathrm{c}}\lambda}.
\end{equation}

Naturally, the laser pulse penetrates the target near the laser axis at first.
With the increase of the laser pulse strength at later times, the target becomes transparent at all points where the condition \eqref{a0reqMax} for laser penetration is satisfied.

However, even for $  a_{0} < a_{0}^{\mathrm{max}}$, the laser pulse may penetrate through the foil at least partially due to the skin effect \cite{Gamaliy1990}.
The relativistically corrected skin depth is $ l_{\mathrm{skin}}=\sqrt{\gamma}c/\omega_{\mathrm{p}} $, where $ \omega_{\mathrm{p}}=\sqrt{4\pi n_{\mathrm{p}}e^{2}/m_{e} } $ is the plasma frequency \cite{Gibbon2005}.
The radiation pressure of the laser pulse will first squeeze the plasma layer of initial thickness $ l $ to a thin slab
of thickness $ d $ while the density will grow by the factor of $ l/d $.
The corresponding skin depth is $ c/\omega_{0}\sqrt{\gamma
n_{\mathrm{c}}d/(n_{\mathrm{p}}l)} $.
As is shown later in the text, the amplitude of the resulting field is enhanced due to the interference of the transmitted and scattered beams.
Since the strength of the scattered beam should be non-negligible in order to noticeably contribute to the interference, we consider the target to be transparent if the intensity of the transmitted (scattered) beam is at least of the order of $ \exp\left( -3/2\right)  $ of the incident intensity.
Thus, assuming the intensity of the incident laser pulse is reduced by a factor of $ \exp\left( -3/2\right)  $ after passing through the squeezed plasma layer of a thickness $ d $, we obtain
\begin{equation}\label{intensity_drop}
2k_{0}\sqrt{\dfrac{ldn_{\mathrm{p}}}{\gamma
n_{\mathrm{c}}}}=\dfrac{3}{2},
\end{equation}
where $ k_{0}=\omega_{0}/c $.
After finding the value of $ d $ from Eq.~\eqref{intensity_drop} and inserting into $ a_{0}^{\mathrm{t}}=n_{\mathrm{p}}\pi \left( l-d\right) /n_{\mathrm{c}}\lambda $ one may obtain the characteristic amplitude of the field required for penetration of the target
\begin{equation}\label{a0reqt}
a_{0}^{\mathrm{t}}=\dfrac{4/3\sqrt{2}k_{0}^{2}l^{2}n_{\mathrm{p}}}{\left( 3/4+8/3\sqrt{2}k_{0}l\right) n_{\mathrm{c}}}.
\end{equation}

For the ultra-thin foils with the thickness of the order of 10 nm assumed in this paper, $ a_{0}^{\mathrm{t}} \approx 0.5a_{0}^{\mathrm{max}} $.
While the most intense part of the laser pulse propagates through the thin plasma slab and a relativistic aperture is formed the electrons at the borders of the aperture oscillate in the direction of the resulting electric field, which is given by the interference of the incident laser wave with the wavefronts reflected from relativistically overdense plasma slab \cite{GonzalezIzquierdo2016Optically}.
In the case of s-polarization, target electrons are simultaneously driven towards the laser beam axis by the ponderomotive force of the two nodes created at the target front side twice per laser period \cite{GonzalezIzquierdo2016Optically}.
It results in an increase of the target density on the axis and therefore two separate relativistic plasma apertures are created in the vicinity of the laser beam axis at first, see Fig.~\ref{fig:schematic}(a) \cite{GonzalezIzquierdo2016Optically}.
It agrees with a double lobe profile of the transmitted laser field observed in three-dimensional simulations and the corresponding asymmetry measured in the distribution of accelerated electrons and ions \cite{Yin2011,Gray2014,GonzalezIzquierdo2016Towards,GonzalezIzquierdo2016Optically}.
We assume that at the radial position $ r $, i.e., at the outer edge of the relativistic plasma aperture, the laser pulse strength is equal to $ a_{0}^{\mathrm{t}}=a_{0}\exp\left( -r^{2}/w_{0}^{2} \right)  $.
As the high density along the beam axis prevents the transmission of the beam close to the axis, the strength of the field passing through one aperture can be estimated as $ a(r/2)=a_{0}\exp\left[ -r^{2}/\left( 2w_{0}\right)^{2} \right]  $.
The resulting field at the rear side of the target is created by the interference of two transmitted beams, each having the above mentioned strength.
Therefore, the amplitude of the diffracted field is approximately given by
\begin{equation}\label{a0resS}
a_{0}^{\mathrm{s}}=2a(r/2)=2a_{0}\exp\left[ -r^{2}/\left( 2w_{0}\right)^{2} \right].
\end{equation}
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{schematic}
\caption{Schematic of target geometry during the interaction with (a) s-polarized and (b) p-polarized laser pulse. (c) Minimal and (d) characteristic aperture radius for penetration of the foil.}
\label{fig:schematic}
\end{figure}

When p-polarization is assumed, the aperture is created on the axis at first, see Fig.~\ref{fig:schematic}(b).
Contrary to the previous case, dense electron bunches are driven alternatively from two opposite edges of the aperture towards the laser axis at the frequency $ 2\omega_{0} $  \cite{Duff2020}.
The photons of the transmitted laser pulse are thus scattered by these relativistic electron bunches generating a number of high harmonics.
The intensity of the $ n $th harmonic can be approximated by a power-law $ I_{n}\propto n^{-8/3} $ \cite{Baeva2006,Edwards2020}.
Assuming that laser field of strength $ a_{0}^{\mathrm{t}} $ is scattered by the electron bunches at two opposite sides of the aperture, the interference of the generated harmonics will create a field of strength $ 2g a_{0}^{\mathrm{t}} $, where $ g=\sqrt{\sum_{n=2}^{N}n^{-8/3} } $ is the fraction of the incoming field scattered into $ N $ harmonics.
Due to their further interference with the field propagating along the laser axis, the resulting field amplitude can be approximated as
\begin{equation}\label{a0resP}
a_{0}^{\mathrm{p}}= 2g a_{0}^{\mathrm{t}}+a_{0}.
\end{equation}
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{s}
\caption{The intensity of the s-polarized laser pulse and the density of target electrons for $ I_{0}=10^{22}~\mathrm{W/cm^{2}} $ and $ \lambda=805~\mathrm{nm} $. Target thickness is 30 nm. Black contours in electron density represent the relativistic critical density $ n_{\mathrm{c}\gamma}^{\mathrm{avg}} $.}
\label{fig:s}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{p}
\caption{The intensity of the p-polarized laser pulse and the density of target electrons for $ I_{0}=10^{22}~\mathrm{W/cm^{2}} $ and $ \lambda=805~\mathrm{nm} $. Target thickness is 30 nm. Black contours in electron density represent the relativistic critical density $ n_{\mathrm{c}\gamma}^{\mathrm{avg}} $.}
\label{fig:p}
\end{figure}

However, in both cases, the field amplification is expected to be limited by the target thickness that has to be smaller than the relativistically corrected skin depth $ l_{\mathrm{skin}}$.
Since Eq.~\eqref{a0reqt} is maximized at $ l= l_{\mathrm{skin}} $ for given target density and laser wavelength, such a choice presents through Eqs.~\eqref{a0resS} and ~\eqref{a0resP} the upper limit for the resulting field strength.

The longitudinal position of the maximum field strength at the rear side of the target can be estimated from the knowledge of the aperture size and propagation direction of the transmitted or scattered light, see Fig.~\ref{fig:schematic}.
The minimal aperture radius of the foil can be estimated as
\begin{equation}\label{w0Pinit}
w^\mathrm{min}=w_{0}\sqrt{\left|\ln\dfrac{a_{0}^{\mathrm{max}}}{a_{0}}\right| },
\end{equation}
as illustrated in Fig.~\ref{fig:schematic}(c).

In the case of s-polarization, the target surface is perturbed by the action of created electromagnetic nodes.
The propagation direction $ \theta $ (with respect to the laser axis) of the transmitted beam at the rear side of the target is given by the deformation of the target surface caused by radiation pressure and, therefore, it strongly depends on laser and target parameters \cite{GonzalezIzquierdo2018Radiation,McIlvenny2020}.
However,  $ \tan\theta $ can be estimated as the ratio of the laser pulse strengths at the axis and at the radial distance which are sufficiently strong to penetrate the target.
Below the strength limit required for target penetration, the foil deformation is given by the distribution of intensity iso-surfaces of the incoming laser beam.
In such a case, the angle can be estimated as $ \tan\theta= \tau c/\left( \sqrt{2\ln2}w_{0}\right)$.
When the plasma aperture is created, the ratio of longitudinal and radial distances given by iso-surfaces has to be weighted by the maximum $ a_{0}^{\mathrm{max}} $ and the characteristic value $ a_{0}^{\mathrm{t}} $ of field strengths required for penetration which results in
\begin{equation}\label{tanTheta}
\tan\theta=\dfrac{\tau c}{\sqrt{2\ln2}w_{0}}\dfrac{w^\mathrm{min}}{w^{\mathrm{avg}}},
\end{equation}
where
\begin{equation}\label{w0Pavg}
w^\mathrm{avg}=w_{0}\sqrt{\left|\ln\dfrac{a_{0}^{\mathrm{t}}}{a_{0}}\right| }.
\end{equation}
is considered, see Fig.~\ref{fig:schematic}(d).
Assuming, that the laser pulse starts to penetrate once $ a_{0}^\mathrm{max} $ is achieved, one can expect the interference maxima to occur at the distance
\begin{equation}\label{fs}
f_{\mathrm{s}}=w^\mathrm{min}\tan\theta
\end{equation}
from the initial target position.

In the case of p-polarization, the maximum field strength is created by the constructive interference of scattered harmonics with the field of the propagating laser beam.
Thus its location depends on the propagation direction of scattered photons which is a function of electron beam momentum.
However, the angle of the electron bunch propagation direction evolves with time.
At the beginning, the electrons are driven along the $ y $-axis.
At later times, the angle of propagation $ \psi $ can be estimated again from the geometrical properties of the Gaussian laser pulse as
\begin{equation}\label{tanPsi}
\tan\psi=\dfrac{\tau c}{\sqrt{2\ln2}w_{0}}.
\end{equation}
The laser photons are scattered into a number of harmonics, each having a different angle of propagation measured with respect to the electron beam direction.
Among them the second harmonic has a dominant contribution for the interference pattern and the corresponding angle can be estimated as \cite{Landau1980Classical}
\begin{equation}\label{theta2nd}
\theta_{\mathrm{2nd}}=\arccos\left[  \dfrac{1+\beta\cos\left( \psi/2\right) }{2\beta}\right],
\end{equation}
where $ \beta=\sqrt{1-1/\left( \gamma^{\mathrm{avg}}\right) ^{2}} $, $ \gamma^{\mathrm{avg}}=\sqrt{1+\left( a_{0}^{\mathrm{t}}\right) ^{2}/2}$ and $ \psi/2 $ is the average value of the angle at which the electron bunch propagates.
Thus, the maximum field strength is expected to be achieved at the distance of
\begin{equation}\label{fp}
f_{\mathrm{p}}=w^{\mathrm{avg}}\tan\left( \theta_{\mathrm{2nd}}+\psi/2\right) .
\end{equation}

\section{2D and 3D simulation results}
We compare our theoretical predictions with 2D/3D PIC simulations performed in the code {\sf EPOCH} \cite{Ridgers2014,Arber2015}.
The 2D simulation box resolved with 10 125 $ \times $ 4 500 cells is spanning from -30 $ \mathrm{\mu m} $ to 15 $\mathrm{\mu m} $ in the $ x $-direction and from -10 $\mathrm{\mu m} $ to  10 $\mathrm{\mu m} $ in the $ y $-direction.
The laser pulse of the peak intensity $ I_{0}=10^{22}~\mathrm{W/cm^{2}} $ ($ a_{0}=69 $) enters the simulation box at a boundary $ x=15~\mathrm{\mu m} $ and is focused to a spot of radius $ w_{0}=1.5~\mathrm{\mu m} $ located at $ x=0~\mathrm{\mu m} $.
It has a wavelength $ \lambda=805~\mathrm{nm} $ and a Gaussian temporal envelope of FWHM duration $ \tau=30~\mathrm{fs}$ in laser intensity.
We performed simulations for both s- and p-polarization of the laser pulse.
The $ \mathrm{Al}^{13+} $  solid foil with the corresponding free electron density 450$ n_{\mathrm{c}} $ and the thickness $ 10-60~\mathrm{nm}  $ is located at the focal spot.

The distributions of the transmitted laser intensity and the density of target electrons are shown in Figs.~\ref{fig:s} and \ref{fig:p} for s- and p-polarization, respectively.
The black contours correspond to the relativistic critical plasma density $ n_{\mathrm{c}\gamma}^{\mathrm{avg}}=\gamma^{\mathrm{avg}}n_{\mathrm{c}}  $.
In the first figure it is shown, that during the penetration of the target having a thickness of $ 30~\mathrm{nm} $, the laser field has a two-lobe structure (at $ x=-1.6~\mathrm{\mu m} $) creating later a symmetric field distribution at $ x=-2.4~\mathrm{\mu m} $.
This is a significant difference from the latter case, in which an asymmetric field diffraction pattern is created due to the interference of the penetrating laser field with the field scattered by the dense electron bunches which are driven alternately from two opposite sides of the aperture (see contours of the critical density near $ x= -0.6~\mathrm{\mu m}$, $ y=2~\mathrm{\mu m} $ and $ x= -1~\mathrm{\mu m}$, $ y=-2~\mathrm{\mu m} $ in Fig.~\ref{fig:p}).
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{order_and_focus}
\caption{(a) Longitudinal position of the maximum field strength at the rear side of the target given by Eqs.~\eqref{fs} and \eqref{fp} and obtained from PIC simulation. (b) The ratio of the first ten harmonics of the scattered light to the incident field of strength $ E_{0} $ obtained from PIC simulations and compared to the theory described in Refs.~\cite{Baeva2006,Edwards2020} for p-polarized laser beam.}
\label{fig:orderandfocus}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{test}
\caption{The electric field intensity distribution in $ k_{x}-k_{y} $ phase-space for (a) s- and (b) p-polarization. The laser of intensity $ I_{0}=10^{22}~\mathrm{W/cm^{2}} $ and wavelength $ \lambda=805~\mathrm{nm} $ irradiates the 30 nm thick target. The color scale is saturated to enhance the visibility of the higher frequency components.}
\label{fig:test}
\end{figure}

The longitudinal position of the maximum field strength generated by the interference pattern at the rear side of the target is shown in Fig.~\ref{fig:orderandfocus}(a).
The expected values calculated from Eqs.~\eqref{fs} and \eqref{fp} ($ x $-axis) are compared to the results obtained from PIC simulations ($ y $-axis) for different target thicknesses.
The solid line is added to guide the eye for the perfect agreement between the theory and simulation results.
As expected according to the theory, the focal distance becomes shorter as the target thickness grows.

In Fig.~\ref{fig:orderandfocus}(b) we compare the ratio of the first ten harmonics of the scattered light to the incident field of amplitude $ E_{0} $ in case of p-polarized laser beam.
The respective harmonics are obtained by filtering out all other frequency components in the Fourier transform of the laser electric field.
The results obtained from the simulation with $ 30~\mathrm{nm} $ target are in good agreement with the theory described in Ref.~\cite{Baeva2006,Edwards2020} which was used for calculating the resulting field amplitude $ a_{0}^{\mathrm{p}} $ in Eq.~\eqref{a0resP}.
It can be seen that the field strength for higher harmonics is small and it can be neglected.
This justifies our choice of $ N=10 $ in Eq.~\eqref{a0resP}.
The electric field distributions in $ k_{x}-k_{y} $ phase-space ($ \mathbf{k} $ is the wave-vector) for (a) s- and (b) p-polarization are shown in Fig.~\ref{fig:test} for the same laser and target parameters.
This spectrum is calculated only for the field which penetrates through the target.
There is a clear evidence of high harmonics propagating at an angle in case of p-polarization in contrast to s-polarization.
The $ k_{x} $ spectrum is broad as the maximum intensity due to interference is very well localized corresponding to a short spatial profile of the pulse in the longitudinal direction.
The $ k_{x} $ and $ k_{y} $ are normalized to the vacuum wavenumber of the laser $ k_{0}=2\pi/\lambda $.

\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{Fig_01}
\caption{Maximum intensity $ I $ achieved in the interaction of the s- and p-polarized Gaussian laser beam characterized by $ I_{0}=10^{22}~\mathrm{W/cm^{2}} $ with the target of different thickness. Lines represent the expected values given by Eqs.~\eqref{a0resS} and ~\eqref{a0resP}, markers show the results obtained from PIC simulations. The relativistically corrected skin depth is $ 42~\mathrm{nm} $.}
\label{fig:fig01}
\end{figure}

Figure \ref{fig:fig01} presents the comparison of resulting field amplitude as a function of target thickness for both types of a linearly polarized laser beam.
The simulation results show that the peak intensity of the diffracted field gradually grows with the target thickness as predicted by Eqs.~\eqref{a0resS} and Eqs.~\eqref{a0resP}.
The thicker the target is, the stronger field is required for penetration, and thus the resulting diffracted field achieves higher amplitude.
However, increasing the target thickness beyond the relativistically corrected skin depth $ l_{\mathrm{skin}}=42~\mathrm{nm} $ results in more pronounced laser energy absorption in the target and thus the intensity of the diffracted field becomes reduced.
The obtained results are therefore in good agreement with the theory as it  has a natural limit for target thickness given by the relativistically corrected skin depth.
This absorption is not taken into account in the theory and thus the theoretical curve stays constant.
In case of p-polarization, the first ten harmonics were considered in Eq.~\eqref{a0resP}.

In order to check the validity of our results and get more realistic information on the field amplification due to the interference of the laser field, we carried out an additional 3D simulation for different laser and target parameters ($ \lambda=1~\mathrm{\mu m} $, $ \sin^{2} $ temporal envelope, foil thickness 20~nm, $ n_{\mathrm{p}}=835n_{\mathrm{c}} $).
The simulation box of dimensions $ 44~\mathrm{\mu}m\times17~\mathrm{\mu}m\times17~\mathrm{\mu}m $ was resolved with $ 8~800\times680\times680 $ cells.
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{test1}
\caption{The distribution of the laser intensity in horizontal and vertical slices of the laser pulse after the interaction with a foil as obtained from 3D simulation.}
\label{fig:puls}
\end{figure}
As shown in Fig.~\ref{fig:puls}, the characteristic structures of the diffracted laser field in a plane perpendicular or parallel to the polarization axis (see Figs.~\ref{fig:s} and ~\ref{fig:p}) are present also in 3D case.
The horizontal and vertical slices of the laser beam correspond to p- and s-polarization from 2D, respectively.
The horizontal plane in 3D contains the electric field vector of the linearly polarized laser beam.
The laser field amplification due to the diffraction on an evolving plasma aperture in 3D case can be roughly estimated as the product of contributions from two planes which are perpendicular ('s') and parallel ('p') to the polarization direction and which were studied above in 2D case.
This approach predicts the maximum amplification of the laser intensity by a factor of 7.0 at a distance $ 3.7~\mathrm{\mu}m $ behind the shutter, which is comparable to the maximum value 7.1 obtained from 3D simulation at a distance $ 4.3~\mathrm{\mu}m $, see Fig.~\ref{fig:test2}, having in mind that in the simulation it is not guaranteed that the field interference maxima in both planes appear simultaneously at the same position.

\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{test2}
\caption{The distribution of the laser intensity in the horizontal slice of the laser pulse as obtained from 3D simulation at the moment of reaching the maximum value.}
\label{fig:test2}
\end{figure}

So far we have assumed an idealized laser pulse.
However, the laser systems are subject to a limited contrast \cite{Weber2017}.
Due to the interaction of the laser pulse pedestal with the target, the pre-plasma is formed before the main laser pulse arrives and thus it can significantly affect the interaction dynamics.
Considering the above mentioned laser intensity, the pedestal of the main laser pulse can completely destroy the nm-scale foil.
The foil expands to several microns along the propagation axis and therefore its density is below the relativistically corrected critical value for the main laser pulse.
As a result, the main laser pulse propagates in the expanded plasma without any significant distortion or focusing.
To avoid such a scenario in a real experiment, one might consider the interaction with two foils.
The pedestal of the main laser pulse interacts with the first foil only.
The main laser pulse then propagates trough the created under-dense plasma and interacts with the second foil.
Such a setup allows creating a relativistic plasma aperture in the second foil which is essential for laser pulse diffraction resulting in the increase of the laser intensity \cite{Nikl2021} which can be used, e.g., for enhancement of ion acceleration from an additional target \cite{Matys2021}.

\section{Conclusion}
In conclusion, we have studied the process of laser pulse diffraction on the relativistic plasma aperture created in the interaction with an over-dense ultra-thin foil.
Such a setup allows amplifying the laser intensity due to the interference of the transmitted/diffracted beam.
We provide the theoretical estimates for obtaining the maximum field strength and its spatial location as a function of the laser and target parameters.
Using 3D numerical simulations we have verified that the laser intensity can be increased at least seven-times when proper laser and target parameters are used.
This setup, based solely on relativistic plasma optics, thus presents a viable approach for obtaining localized laser intensity enhancement which might be interesting for applications like ion acceleration, generation of gamma rays and electron positron pair creation.

\section*{Acknowledgements}
Portions of this research were carried out at ELI Beamlines, a European user
facility operated by the Institute of Physics of the Academy of Sciences of the
Czech Republic.
This work is supported by the project High Field Initiative (HIFI) CZ.02.1.01/0.0/0.0/15\_003/0000449 from European
Regional Development Fund (ERDF). The support of Czech Science Foundation project No. 18-09560S is acknowledged. The support of Grant Agency of the Czech Technical University in Prague is appreciated, grant no. SGS19/192/OHK4/3T/14. This work was supported by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project "IT4Innovations National Supercomputing Center – LM2018140". Computational resources were supplied by the project "e-Infrastruktura CZ" (e-INFRA LM2018140) provided within the program Projects of Large Research, Development and Innovations Infrastructures.
The simulations were performed at the cluster ECLIPSE at ELI Beamlines.