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

This file is documentation for the SIAM \LaTeX\ style, including how
to typeset the main document, the {\scshape Bib}\TeX\xspace\ file, and any supplementary
material. More information
about SIAM's editorial style can be found in the style manual, available
at \url{https://www.siam.org/journals/pdf/stylemanual.pdf}.
The major changes in the SIAM online-only class are summarized in \cref{sec:changes}.
The SIAM \LaTeX\@ files can be found at
\url{https://www.siam.org/journals/auth-info.php}. The files that
are distributed are given below. 
\begin{itemize}
\item \texttt{siamonline190516.cls} (required): Main SIAM online-only \LaTeX\ class file.
\item \texttt{siamplain.bst} (required): Bibliographic style file for
  {\scshape Bib}\TeX\xspace.
\item \texttt{docsiamonline.tex}: Produces this documentation.
\item \texttt{references.bib}: {\scshape Bib}\TeX\xspace\ database for this
  documentation and examples.
\item \texttt{ex\_article.tex}: Template for article.
\item \texttt{ex\_supplement.tex}: Template for supplement.
\item \texttt{ex\_shared.tex}: Template for shared information for
  article and supplement.
\end{itemize}
To use these files, put \texttt{siamonline190516.cls} and
\texttt{siamplain.bst} in the directory with your
paper or, alternatively, into your \LaTeX\@ and {\scshape Bib}\TeX\xspace\@ paths,
respectively. 
The outline of a SIAM \LaTeX\ article is shown in
\cref{ex:outline}. Templates are provided and discussed in more detail
in \cref{sec:template}.

\section{Class options}
\label{sec:class-options}

Class options can be included in the bracketed argument of the
command, separated by commas. The possible class options are:
\begin{itemize}
\item \code{review} --- Recommended for submitting your manuscript to
  a SIAM journal.  Adds line numbers as well as the statement ``This
  manuscript is for review purposes only'' to the bottom of each page.
\item \code{final} --- Turns off the black boxes that help authors
  identify lines that are too long. The final published version will
  have this option on.
\item \code{supplement} --- Specifies that the file is a supplement
  and not the main document, causing changes in the appearance of the
  title and numbering; see \cref{sec:supplement} for details.
\end{itemize}

\begin{example}[label={ex:outline},listing only,%
  listing options={style=siamlatex,{morekeywords=[1]{maketitle},
      morekeywords=[2]{siamonline190516}},}]%
  {Document outline}
\documentclass{siamonline190516}
\begin{document}
\maketitle
\end{document}
\end{example}

\section{Front matter}
\label{sec:front}
The title and author parts are formatted using the standard
\code{\title}, \code{\author}, and \code{\maketitle} commands as
described in Lamport \cite{La86}. The title and author should be
declared in the preamble. The title and author names are automatically
converted to uppercase in the document.
If there is more than one author, each additional author should be preceded by
the \code{\and} command.  
The addresses and support acknowledgments are added via
\code{\thanks}. Each author's thanks should specify their address.
The support acknowledgment should be put in the title thanks,
unless specific support needs to be specified for individual authors,
in which case it should follow the author address.  
The header for this file was produced by the code in \cref{ex:header}, including an
example of a shared footnote. Each thanks produces a footnote, so the
footnote of the second author is \#3.
The command
\code{\headers{title}{authors}} command, with
the title (possibly shortened to fit) and the authors' names, creates
the page headers, automatically converted to uppercase.

~

\examplefile[label={ex:header},listing only,%
listing options={style=siamlatex,%
deletetexcs={and,thanks,title,author},%
{moretexcs=[2]{and,thanks,title,author,maketitle,headers,email}}}
]{Title and authors in preamble}{tmp_\jobname_header.tex}

Following the author and title is the abstract, key words listing, and AMS subject
classifications, designated using the \code{abstract}, \code{keywords}, and \code{AMS}
environments. 
Authors are responsible for providing AMS numbers which can be found
on the AMS web site \cite{AMSMSC2010}. The abstract, keywords, and AMS subject classifications
for this document are specified in \cref{ex:abstract}.

\examplefile[label={ex:abstract},%
before upper={\preamble{\boldsymbol newcommand\{\boldsymbol BibTeX\}\{\{\boldsymbol scshape Bib\}\boldsymbol TeX\boldsymbol xspace\}}},
listing only,%
listing options={style=siamlatex,%
{morekeywords=[2]{abstract,keywords,AMS}}}
]{Abstract, keywords, and AMS classifications}{tmp_\jobname_abstract.tex}

A more complete example, including a PDF supplement, that uses the
included files \texttt{ex\_article.tex}, \texttt{ex\_supplement.tex},
and \texttt{ex\_shared.tex} is discussed in
\cref{sec:template}. The example files can be used as a starting point
for producing a document.

\section{Cross references and hyperlinks}
\label{sec:cr+hyp}

SIAM now supports cross references and hyperlinks via the
\texttt{cleveref} and \texttt{hyperef} packages, which are loaded by
the class file.

\subsection{Cleveref}
\label{sec:cleveref}

SIAM strongly recommends using the commands provided by 
the \texttt{cleveref} package for cross referencing.
The package is automatically loaded and already customized
to adhere to SIAM's style guidelines. 
To create a cross reference, use the command \code{\cref} (inside
sentence) or \code{\Cref} (beginning of a sentence) in place of the
object name and \code{\ref}. 
The \texttt{cleveref}
package enhances \LaTeX's cross-referencing features, allowing
the format of cross references to be determined automatically
according to the ``type" of cross reference (equation, section, etc.)
and the context in which the cross reference is used.
So, the package
\emph{automatically} inserts the object name as well as the
appropriate hyperlink; see \cref{ex:cref}.
It may require two \LaTeX\@ compilations for the references to show up
correctly.
Additional examples are shown in the sections below
for equations, tables, figures, sections, etc.

\begin{example}[label=ex:cref,bicolor,listing options={style=siamlatex,%
{morekeywords=[2]{cref,ref}}}]{Advantage of using cleveref}
The normal way to get a cross reference with a hyperlink requires a
lot of typing: \hyperref[thm:mvt]{Theorem~\ref*{thm:mvt}}.
The \texttt{cleveref} package gets both the name and hyperlink
automatically using a single macro: \cref{thm:mvt}.
It also handles multiple references with the same macro, such as
\cref{thm:mvt,fig:pgfplots,fig:testfig}.
\end{example}

\subsection{Hyperef}
\label{sec:hyperef}

Hyperlinks are created with the \code{\href} and \code{\url} commands,
as shown in \cref{ex:href}.
SIAM has also defined the \code{\email} command, as shown in
\cref{ex:header}.

\begin{example}[label={ex:href},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{href,url}}}]{Creating hyperlinks}
The \href{https://www.siam.org}{SIAM homepage} has general information.
There are times when the author may want to specify
the location explicitly instead by using \url{https://www.siam.org}.
\end{example}

Note that homepage links via \code{\url} in the \code{\thanks} environment require
special formatting for the tilde (\string~) character. The formatting
is used in the template and shown in \cref{ex:shared}.


\section{Math and equations}
\label{sec:math}
Here we show some example equations, with numbering, and examples of
referencing the equations. SIAM now includes the package 
\texttt{amsmath} by default, and we include some of its features as
well, although the reader should consult the package user manual for
further guidance \cite{amsmath,shortmath}.
 Several of the example are adapted 
from Mittlebach and Goossen's guide to \LaTeX~\cite{MiGo04}.

\Cref{ex:textmath} is a straightforward example of inline mathematics
equations that does not use any special packages or features.

\begin{example}[label={ex:textmath},bicolor]{Inline math}
The following shows an example of math in text: 
Let $S=[s_{ij}]$ ($1\leq i,j\leq n$) be a $(0,1,-1)$-matrix of order $n$.
\end{example}

In \cref{ex:bbm}, we show the recommended method for getting
blackboard fonts using the \texttt{amsfonts} package. This is not
loaded by default and must be included in the preamble.

\begin{example}[label={ex:bbm},bicolor,before upper={\preamble{\boldsymbol
      usepackage\{amsfonts\}}},%
listing options={style=siamlatex,%
{morekeywords=[2]{mathbb}}}]{Blackboard math}
Blackboard bold characters, such as $\mathbb{C}$ and $\mathbb{R}$,
should be created with the \texttt{amsfonts} package, although this
is not included by default.
\end{example}

\newpage
\Cref{ex:smallmatrix} shows the \code{smallmatrix} environment for an
inline matrix from the \texttt{amsmath} package, which is included by
default.

\begin{example}[label={ex:smallmatrix},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{smallmatrix}}}]{Inline matrix}
Matrices of no more than two rows appearing in text can be created
as shown in the next example:
$B = \bigl[ \begin{smallmatrix} B_{11} & B_{12} \\ 
  B_{21} & B_{22} \end{smallmatrix} \bigr]$.
\end{example}

Bigger matrices can be rendered with environments from
the \texttt{amsmath} package, such as \code{bmatrix} and
\code{pmatrix} used in
\cref{ex:matrices}.

\begin{example}[label={ex:matrices},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{bmatrix,pmatrix}}}]{Creating matrices}
Display matrices can be rendered using environments from \texttt{amsmath}:
\begin{equation}\label{eq:matrices}
  S=\begin{bmatrix}1&0\\0&0\end{bmatrix}
  \quad\text{and}\quad
  C=\begin{pmatrix}1&1&0\\1&1&0\\0&0&0\end{pmatrix}.
\end{equation}
\Cref{eq:matrices} shows some example matrices.
\end{example}

\Cref{ex:dmo} shows how to use the \code{\DeclareMathOperator} command
from the \texttt{amsopn} package to declare the \code{\Range} macro. 
(This example also uses the \texttt{braket} package for the
\code{\set} macro, but this is not necessarily recommended by SIAM.)
\begin{example}[label={ex:dmo},%
before upper={\preamble{\boldsymbol usepackage\{braket,amsfonts,amsopn\}}\\
\noindent\preamble{\boldsymbol DeclareMathOperator\{\boldsymbol Range\}\{Range\}}},%
bicolor,%
listing options={style=siamlatex,%
{moretexcs=[2]{Range}}}
]{Declaring math operators}
An example of a math operator:
\begin{equation}\label{eq:range}
  \Range(A) = \set{ y \in \mathbb{R}^n | y = Ax }.
\end{equation}
\end{example}

\newpage
\Cref{ex:foo} shows how to use the \code{align} environment from
\texttt{amsmath} to easily align multiple equations.

\begin{example}[label={ex:foo},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{align}}}]{Aligned equations}
\Cref{eq:a,eq:b,eq:c} show three aligned equations.
\begin{align} 
  f &= g, \label{eq:a} \\
  f' &= g', \quad\text{and} \label{eq:b} \\
  \mathcal{L}f &= \mathcal{L}g \label{eq:c}.
\end{align}
\end{example}

Another way to number a set of equations is the
\code{subequations} environment from \texttt{amsmath}, as shown in \cref{ex:aligned}.

\begin{example}[label={ex:aligned},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{subequations}}}]{Subequations}
We calculate the Fr\'{e}chet derivative of $F$ as follows:
\begin{subequations}
\begin{align}
  F'(U,V)(H,K) 
  &= \langle R(U,V),H\Sigma V^{T} + U\Sigma K^{T} -
  P(H\Sigma V^{T} + U\Sigma K^{T})\rangle \label{eq:aa} \\
  &= \langle R(U,V),H\Sigma V^{T} + U\Sigma K^{T}\rangle 
  \nonumber \\
  &= \langle R(U,V)V\Sigma^{T},H\rangle + 
  \langle \Sigma^{T}U^{T}R(U,V),K^{T}\rangle. \label{eq:bb}
\end{align}
\end{subequations}
\Cref{eq:aa} is the first line, and \cref{eq:bb} is the last line.
\end{example}

For an equation split over multiple lines, \cref{ex:ml} shows the
usage of the \code{multline} environment provided by \texttt{amsmath}.
\begin{example}[label={ex:ml},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{multline}}}]{Equation split across lines}
We claim that the projection $g(U,V)$ is given by the pair of matrices:
\begin{multline} \label{eq:ml}
  g(U,V) = \biggl( \frac{R(U,V)V\Sigma^{T}U^{T} 
    - U\Sigma V^{T}R(U,V)^{T}}{2}U,\\
  \frac{R(U,V)^{T}U\Sigma V^{T}-V \Sigma^{T}U^{T}R(U,V)}{2}V \biggr).
\end{multline}
\end{example}

\section{Theorem-like environments}
\label{sec:thm}
SIAM loads \texttt{ntheorem} package and uses it to define the following
theorem-like environments:
\code{theorem}, 
\code{lemma}, 
\code{corollary},
\code{definition}, and
\code{proposition}.
SIAM also defines a \code{proof} environment that automatically
inserts the symbol ``$\,\proofbox\,$'' at the end of any proof, even if it ends in an
equation environment. \emph{Note that the document may need to be
  compiled twice for the mark to appear.}
Some of the calculus examples were adapted from \cite{CalcI}. 
\Cref{ex:theorem} shows usage of the \code{theorem} environment. 
An optional argument can be used to name the theorem.
\Cref{ex:cor} illustrates a corollary, without a name, and the
proof environment.

\begin{example}[label=ex:theorem,bicolor,parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{theorem}}}]{Theorem}
\begin{theorem}[Mean Value Theorem]\label{thm:mvt}
  Suppose $f$ is a function that is continuous on the closed interval
  $[a,b]$.  and differentiable on the open interval $(a,b)$.
  Then there exists a number $c$ such that $a < c < b$ and
  \begin{displaymath}
    f'(c) = \frac{f(b)-f(a)}{b-a}.
  \end{displaymath}
  In other words, $f(b)-f(a) = f'(c)(b-a)$.
\end{theorem}
\end{example}

\begin{example}[label=ex:cor,bicolor,parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{corollary,proof}}}]%
{Corollary and proof}
\begin{corollary}
  Let $f(x)$ be continuous and differentiable everywhere. If $f(x)$
  has at least two roots, then $f'(x)$ must have at least one root.
\end{corollary}
\begin{proof}
  Let $a$ and $b$ be two distinct roots of $f$.
  By \cref{thm:mvt}, there exists a number $c$ such that
  \begin{displaymath}
    f'(c) = \frac{f(b)-f(a)}{b-a} = \frac{0-0}{b-a} = 0.
  \end{displaymath}
\end{proof}
\end{example}

\newpage
SIAM also defines commands to create your own theorem- and
remark-like environments:
\begin{itemize}
\item \code{newsiamthm} --- Small caps header, italized body. 
\item \code{newsiamremark} --- Italics header, roman body.
\end{itemize}
Each command takes two arguments. The first is the environment name,
and the second is the name to show in the document. These commands should
be used instead of \code{\newtheorem}.
\Cref{ex:claim,ex:ref} shows how to use the commands above, including how to
specify the plural version for \texttt{cleveref} if it is unusual.

\begin{example}[label=ex:claim,bicolor,%
before upper={\preamble{\boldsymbol newsiamthm\{claim\}\{Claim\}}\\
\noindent\preamble{\boldsymbol newsiamremark\{hypothesis\}\{Hypothesis\}}\\
\noindent\preamble{\boldsymbol crefname\{hypothesis\}\{Hypothesis\}\{Hypotheses\}}},%
parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{claim,proof,hypothesis}}}]{New theorem-like environment}
\begin{claim}\label{cl:constant}
  If $f'(x) = 0$ for all $x \in (a,b)$ then $f(x)$ is constant on $(a,b)$.
\end{claim}
\begin{hypothesis}\label{hyp1}
The function $f$ is continuously differentiable.
\end{hypothesis}
\begin{hypothesis}\label{hyp2}
The random variable is normally distributed.
\end{hypothesis}
\end{example}


\begin{example}[label=ex:ref,bicolor,listing options={style=siamlatex,%
{morekeywords=[2]{cref}}}]{References}
We can reference multiple types of objects with a single reference:
\cref{cl:constant,thm:mvt,hyp1}.
\end{example}


\section{Tables}
\label{sec:tab}
Table captions should go above the tables.
\Cref{ex:simpletable} shows the code to generate a
\cref{tab:simpletable}.
A more complicated example is shown in 
\cref{ex:table}, which generates \cref{tab:KoMa14}. This
example uses subfloats via the \texttt{subfig} package, as well as
special column options from the \texttt{array} package.

\begin{tcbverbatimwrite}{tmp_\jobname_simpletable.tex}
\begin{table}[tbhp]
{\footnotesize
  \caption{Example table}\label{tab:simpletable}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
}
\end{table}
\end{tcbverbatimwrite}

\examplefile[label={ex:simpletable},%
listing only, listing options={style=siamlatex}]%
{Example table.}{tmp_\jobname_simpletable.tex}

\input{tmp_\jobname_simpletable.tex}

\begin{tcbverbatimwrite}{tmp_\jobname_table.tex}
\newcolumntype{R}{>{$}r<{$}}  %
\newcolumntype{V}[1]{>{[\;}*{#1}{R@{\;\;}}R<{\;]}} %
\begin{table}[tbhp]
{\footnotesize
 \captionsetup{position=top}
  \caption{Example table adapted from Kolda and Mayo \rm{\cite{KoMa14}}.}\label{tab:KoMa14}
\begin{center}
  \subfloat[$\beta=1$]{
    \begin{tabular}{|r|R|V{3}|c|r@{\,$\pm$\,}l|} \hline
occ. & \multicolumn{1}{c|}{$\lambda$} & \multicolumn{4}{c|}{$\mathbf{x}$} & 
fevals &  \multicolumn{2}{c|}{time (sec.)}\\ \hline
 718 & 11.3476 &  0.5544 &  0.3155 &  1.2018 &  0.0977 &   45 & 0.17 & 0.06 \\ \hline 
 134 &  3.7394 &  0.2642 & -1.1056 &  0.2657 & -0.3160 &   31 & 0.12 & 0.05 \\ \hline 
   4 & \multicolumn{6}{c|}{\emph{--- Failed to converge ---}} & 0.21 & 0.10 \\ \hline 
 \end{tabular}}

  \subfloat[$\beta=-1$]{
    \begin{tabular}{|r|R|V{3}|c|r@{\,$\pm$\,}l|} \hline
occ. & \multicolumn{1}{c|}{$\lambda$} & \multicolumn{4}{c|}{$\mathbf{x}$} & 
fevals &  \multicolumn{2}{c|}{time (sec.)}\\ \hline
  72 & -1.1507 &  0.2291 &  0.6444 &  0.3540 & -0.8990 &   34 & 0.14 & 0.06 \\ \hline 
 624 & -6.3985 &  0.1003 &  0.1840 &  0.5305 &  1.2438 &   48 & 0.19 & 0.08 \\ \hline 
   2 & \multicolumn{6}{c|}{\emph{--- Failed to converge ---}} & 0.23 & 0.02 \\ \hline 
 \end{tabular}}
\end{center}
}
\end{table}
\end{tcbverbatimwrite}

\examplefile[label={ex:table},%
before upper={\preamble[\scriptsize]{\boldsymbol usepackage\{array\}}\\[-0.4em]
\noindent\preamble[\scriptsize]{\boldsymbol usepackage[caption=false]\{subfig\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example table with subtables.}{tmp_\jobname_table.tex}

\input{tmp_\jobname_table.tex}

\newpage
\section{Figures}
\label{sec:fig}

It is recommended that all figures be generated in high resolution.
In the past, SIAM has required encapsulated
postscript (EPS) format for final production. This is still an
acceptable format, but SIAM also now allows high-resolution 
PDF, JPEG, and PNG figures. 
If working with EPS images and using \texttt{pdflatex}, we
recommend the package \texttt{epstopdf} to automatically convert EPS
images to PDF for inclusion in PDF documents created by
\texttt{pdflatex}.
\Cref{ex:fig} shows the code to generate \cref{fig:testfig}.  This
example uses the \texttt{graphicx} package for the
\code{\includegraphics} command. 


\begin{tcbverbatimwrite}{tmp_\jobname_fig.tex}
\begin{figure}[tbhp]
  \centering
  \subfloat[$\epsilon_{\max}=5$]{\label{fig:a}\includegraphics{lexample_fig1}}
  \subfloat[$\epsilon_{\max}=0.5$]{\label{fig:b}\includegraphics{lexample_fig2}}
  \caption{Example figure using external image files.}
  \label{fig:testfig}
\end{figure}
\end{tcbverbatimwrite}

\examplefile[label={ex:fig},%
before upper={\preamble[\scriptsize]{\boldsymbol usepackage\{graphicx,epstopdf\}}\\[-0.4em]
\noindent\preamble[\scriptsize]{\boldsymbol usepackage[caption=false]\{subfig\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example figure with subfigures and external files}{tmp_\jobname_fig.tex}

\input{tmp_\jobname_fig.tex}

Another option for figures is a graphics-generator that is platform- and
format-independent.
PGF is a TeX macro package for generating  such graphics and works together with the most important TeX backend drivers, including pdftex and dvips. 
The user-friedly syntax layer called TikZ.
Here we show an example using \texttt{PGFPLOTS}, useful for drawing high-quality 
plots directly in \LaTeX.
\Cref{ex:data} and \cref{ex:pgfplots} shows the data and code, respectively, to generate \cref{fig:pgfplots}, adapted from \cite{pgfplots}.

\examplefile[label={ex:data},listing only,
listing options={style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example data file (data.dat)}{data.dat}

\begin{tcbverbatimwrite}{tmp_\jobname_tikz.tex}
\begin{figure}[tbhp]
  \centering
  \begin{tikzpicture}
    \begin{loglogaxis}[height=2.75in, grid=major,
      xlabel={Degrees of Freedom}, ylabel={$L_2$ Error},
      legend entries={$d=2$,$d=3$}]
      \addplot table [x=d2_dof,y=d2_l2_err] {data.dat};
      \addplot table [x=d3_dof,y=d3_l2_err] {data.dat};
    \end{loglogaxis}
  \end{tikzpicture}
  \caption{Example \texttt{PGFPLOTS} figure.}
  \label{fig:pgfplots}
\end{figure}
\end{tcbverbatimwrite}

\examplefile[label={ex:pgfplots},%
before upper={\preamble[\scriptsize]{\boldsymbol usepackage\{pgfplots\}}},%
listing only, listing options={%
  style=siamlatex}]%
{Example TikZ/PGF for platform-independent graphics.}{tmp_\jobname_tikz.tex}

\input{tmp_\jobname_tikz.tex}

  
\newpage
\section{Algorithms}
\label{sec:algs}

SIAM automatically includes the \texttt{algorithm} package in the
class definition. This provides the float environment.
Users have the choice of \texttt{algpseudocode},
\texttt{algorithmic}, and other packages for actually formatting the algorithm.
For example, \cref{alg:buildtree} is produced by the code in \cref{ex:alg}.
In order to reference lines within the algorithm, we need to tell the
\texttt{cleveref} package how to do the referencing, which is the
second line of \cref{ex:alg}. Then we can use the code
\code{\cref{line3}} to produce \cref{line3}.

\begin{tcbverbatimwrite}{tmp_\jobname_alg.tex}
\begin{algorithm}
\caption{Build tree}
\label{alg:buildtree}
\begin{algorithmic}[1]
\STATE{Define $P:=T:=\{ \{1\},\ldots,\{d\}$\}}
\WHILE{$\#P > 1$}
\STATE\label{line3}{Choose $C^\prime\in\mathcal{C}_p(P)$ with $C^\prime := \operatorname{argmin}_{C\in\mathcal{C}_p(P)} \varrho(C)$}
\STATE{Find an optimal partition tree $T_{C^\prime}$ }
\STATE{Update $P := (P{\setminus} C^\prime) \cup \{ \bigcup_{t\in C^\prime} t \}$}
\STATE{Update $T := T \cup \{ \bigcup_{t\in\tau} t : \tau\in T_{C^\prime}{\setminus} \mathcal{L}(T_{C^\prime})\}$}
\ENDWHILE
\RETURN $T$
\end{algorithmic}
\end{algorithm}
\end{tcbverbatimwrite}

\examplefile[float=htpb,label={ex:alg},%
before upper={\preamble[\scriptsize]{\boldsymbol usepackage\{algorithmic\}}\\[-0.4em]
\preamble[\scriptsize]{\boldsymbol Crefname\{ALC@unique\}\{Line\}\{Lines\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example algorithm}{tmp_\jobname_alg.tex}

\input{tmp_\jobname_alg.tex}

\section{Sections}
\label{sec:sec}
Sections are denoted using standard \LaTeX\ section commands, i.e.,
\code{\section}, \code{\subsection}, etc. 
If you wish to end the section title with something other that a period (the default), you have to add the command \code{\nopunct} at the end of the title.

Appendices are created with the normal sectioning commands, following
the command \code{


\section{A detailed example}

Here we include some equations and theorem-like environments to show
how these are labeled in a supplement and can be referenced from the
main text.
Consider the following equation:
\begin{equation}
  \label{eq:suppa}
  a^2 + b^2 = c^2.
\end{equation}
You can also reference equations such as \cref{eq:matrices,eq:bb} 
from the main article in this supplement.

\lipsum[100-101]

\begin{theorem}
  An example theorem.
\end{theorem}

\lipsum[102]
 
\begin{lemma}
  An example lemma.
\end{lemma}

\lipsum[103-105]

Here is an example citation: \cite{KoMa14}.

\section[Proof of Thm]{Proof of \cref{thm:bigthm}}
\label{sec:proof}

\lipsum[106-112]

\section{Additional experimental results}
\Cref{tab:foo} shows additional
supporting evidence. 

\begin{table}[htbp]
{\footnotesize
  \caption{Example table.}  \label{tab:foo}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
}
\end{table}


\bibliographystyle{siamplain}

\section{Introduction}
\label{s:Introduction}
We consider the maximization of an expensive-to-evaluate oracle $f(\mathbf{x}, t)$,
where $\mathbf{x}\in \bm{\mathcal{X}} \subset \mathbb{R}^d$ are the inputs or \emph{action} parameters
in a compact domain and $t$ is the \emph{context} from a (possibly infinite)
set of contexts $\mathcal{T}$. An observation $\hat{y} \in \mathbb{R}$ made by observing
$f$ at the action-context pair $(\mathbf{x}, t)$ is assumed to be corrupted by
additive stochastic noise $\epsilon$; that is, $\hat{y} = f(\mathbf{x},t) + \epsilon$.
Furthermore, we assume a context space $\mathcal{T}$ whose members have a unique
ordering (hereafter ``time''). We seek an action $\mathbf{x}^*_T$ at a given future time $T \in \mathcal{T}$ that solves
\begin{equation}
\operatornamewithlimits{arg\,max}_{\mathbf{x} \in \bm{\mathcal{X}}} f(\mathbf{x}, T)
\label{eq:prob}    
\end{equation}
in as few stochastic observations of $f$ as possible.
Such a problem fundamentally differs from general
context-dependent bandit problems \cite{krause2011contextual,srinivas2009gaussian} in that we are interested solely in
the action that maximizes the oracle value $f$ 
at time $t=T$, as opposed to maximizing the
\emph{cumulative oracle value} until $T$. Additionally, we
assume that only one observation can be made per time 
$t$.
The time-dependent maximization that we address arises, for example, in quantum
computing and control~\cite{magann2020digital,zhu2018training}, where 
values for the control parameter $\mathbf{x}$ are sought by using
noisy observations $f(\mathbf{x},t)$ 
to achieve a
desired task at future time $T$; see \cref{ss:qoc}. Another application  relevant to our work
is multistage stochastic programming for portfolio optimization \cite[Ch. 1]{shapiro2014lectures}, where the goal is to find a 
policy $\{(\mathbf{x}_i, t_i) :~i=0,\ldots,n-1\}$ for redistributing assets 
to maximize the expected profit at $T=t_n$ and the profit at each stage is
realized by querying a noisy oracle. 
In many settings, $\bm{\mathcal{X}}$ is a simple domain; here 
we assume that $\bm{\mathcal{X}}:= [0,1]^d$ so that the inputs $\mathbf{x}$ are normalized.


The main challenges of the problem \eqref{eq:prob} that we seek to address are (i) the absence of structure (e.g., convexity) in $f$ that
could potentially be exploited, 
(ii) the high cost
associated with each
observation of $f$, and 
(iii) the time-varying setting where observations of $f$ cannot be made in the past or future.
Bayesian optimization (BO)~\cite{brochu2010tutorial, shahriari2015taking} with
Gaussian process (GP) priors~\cite{gramacy2020surrogates,rasmussen:williams:2006}  suits the
specific
challenges posed by our problem (mainly (i) and (ii) mentioned
above),
where observations at each round $i$, 
$\hat{y}_i =
f(\mathbf{x}_i, t_i) + \epsilon_i$, 
are made judiciously by leveraging information
gained from previous observations, as a means of coping with the high costs of
each observation. The key idea  is to specify GP prior distributions on the
oracle value and the noise. That is, $f(\mathbf{x},t) \sim \mathcal{GP}\left(0, k( (\mathbf{x}, t),
(\mathbf{x}',t'))\right)$ and $\epsilon \sim \mathcal{GP}(0,\sigma^2_\epsilon)$, where $\sigma^2_\epsilon$ is a constant noise variance. 

The covariance function $k$ captures the correlation between the observations
in the joint $(\mathbf{x},t)$ space; here we use the product composite form
given by $k((\mathbf{x}, t), (\mathbf{x}', t')) = k_{\mathbf{x}}(\mathbf{x},\mathbf{x}'; \theta_\mathbf{x}) \times k_t(t,t';
\theta_t) $, where $\theta_\mathbf{x}$ and $\theta_t$ 
parameterize the covariance functions for $\mathbf{x}$ and $t$, respectively. We estimate the
GP hyperparameters 
$\bm{\Omega}= \lbrace \theta_x, \theta_t, \sigma^2_\epsilon \rbrace$
from data by maximizing the marginal
likelihood. We use the squared-exponential  kernels defined as $k_\mathbf{x}(\mathbf{x}, \mathbf{x}') = \textnormal{exp} \left( -\frac{\|\mathbf{x} -
\mathbf{x}'\|^2}{2 \theta^2_\mathbf{x}}\right)$ and $k_t(t, t') = \textnormal{exp} \left( -\frac{\|t -
t'\|^2}{2 \theta^2_t}\right)$.
The
posterior predictive distribution of the output $Y$, conditioned on available observations from the oracle, is given by 
\begin{equation}
    \begin{split}
        Y(\mathbf{x}, t |
      \bm{\mathcal{D}}_n, \bm{\Omega}) &\sim \mathcal{GP}(\mu_n(\mathbf{x}, t), \sigma^2_n (\mathbf{x}, t)),~ \quad \bm{\mathcal{D}}_n :=  \lbrace
(\mathbf{x}_i, t_i), \hat{y}_i \rbrace _{i=1}^n \\
\mu_n(\mathbf{x},t) &= \mathbf{k}_n^\top [\mathbf{K}_n + \sigma_\epsilon^2 \mathbf{I}]^{-1} \mathbf{y}_n\\
\sigma^2_n(\mathbf{x}, t) &=  k((\mathbf{x},t), (\mathbf{x},t)) - \mathbf{k}_n^\top [\mathbf{K}_n + \sigma_\epsilon^2 \mathbf{I}]^{-1} \mathbf{k}_n,
    \end{split}
    \label{e:GP}
\end{equation}
where $\mathbf{k}_n$ is a vector of covariances between $(\mathbf{x},t)$ and all observed
points in $\bm{\mathcal{D}}_n$, $\mathbf{K}_n$ is a sample covariance matrix of observed points in
$\bm{\mathcal{D}}_n$, $\mathbf{I}$ is the identity matrix, and $\mathbf{y}_n$ is the vector of all
observations in $\bm{\mathcal{D}}_n$; \eqref{e:GP} is then used as a surrogate model for
$f$ in Bayesian optimization. Note that $\mu_n$ and $\sigma_n^2$
are the posterior mean and variance of the GP, respectively, where the
subscript $n$ implies the conditioning based on $n$ past observations. BO then
proceeds by defining an \emph{acquisition function}, $\alpha(\mathbf{x},t)$, in terms of the GP posterior that is
optimized in lieu of the expensive $f$ to select the next
point $\mathbf{x}_{n+1}$ at $t_{n+1}$, and the process continues until a budget of $q$ observations is reached; see Algorithm~\ref{a:BO}. We denote by $\mathbf{x}_i$ and $\mathbf{x}_t$ a decision made by maximizing the acquisition function at $t_i$ and $t$, respectively. We use $\mathbf{x}_t^*$ to denote a maximizer of $f(\cdot, t)$.


\begin{algorithm2e}[t]
\textbf{Given:} $\bm{\mathcal{D}}_n = \lbrace
(\mathbf{x}_i, t_i), \hat{y}_i \rbrace _{i=1}^n$, 
 total budget $q$, 
 schedule $\lbrace t_{n+1},\ldots, t_{q}= T \rbrace$,
 and GP hyperparameters $\bm{\Omega}$ \\
\KwResult{$(\mathbf{x}_T, \hat{y}_T)$}
  \For{$i=n+1, \ldots, q$, }{
  Find $\mathbf{x}_i \in \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\operatornamewithlimits{arg\,max}}~ \alpha(\mathbf{x}, t_i)$ \qquad (acquisition function maximization based on $\bm{\Omega}$)\\
    Observe $\hat{y}_i$ = $f(\mathbf{x}_i, t_i) + \epsilon_i$\\ 
    Append $\bm{\mathcal{D}}_i = \bm{\mathcal{D}}_{i-1} \cup \lbrace (\mathbf{x}_i, t_i), \hat{y}_i \rbrace$\\ 
    Update GP hyperparameters $\bm{\Omega}$ \\
 }
 \caption{Generic Bayesian optimization}
 \label{a:BO}
\end{algorithm2e}

Typical acquisition functions in BO take a \emph{greedy} (``myopic'')
approach, where each $\mathbf{x}_i$ 
is selected as though it were the
last, without accounting for the potential impact on future selections.
Such myopic acquisition functions include the
probability of improvement (\texttt{PI})~\cite{kushner1964new}, the expected improvement
(\texttt{EI}) \cite{jones1998efficient, mockus1978application}, and the GP upper
confidence bound (\texttt{UCB})~\cite{srinivas2009gaussian}, which are, respectively,
\begin{equation}
    \begin{split}
        \alpha_{\textrm{PI}}(\mathbf{x}, t) &= \Phi \left( z \right) \\
        \alpha_{\textrm{EI}}(\mathbf{x}, t) &= z\sigma_n(\mathbf{x}, t) \Phi\left(z \right) + \sigma_n(\mathbf{x}, t) \phi \left( z \right) \\
        \alpha_{\textrm{UCB}}(\mathbf{x},t) &= \mu_n(\mathbf{x}, t) + \beta^{1/2} \sigma_n(\mathbf{x}, t),
    \end{split}
    \label{e:myopic_acqs}
\end{equation}
where $z = \frac{\mu_n(\mathbf{x}, t) - \xi}{\sigma_n(\mathbf{x}, t)}$; $\Phi$ and $\phi$ denote
the standard normal cumulative and probability density functions, respectively;
$\xi$ is a user-specified target; and $\beta$ is a confidence parameter.
Whereas the \texttt{PI} and \texttt{EI} acquisition functions are probabilistic measures of
improvement over a user-specified target, \texttt{UCB} is an optimistic
estimate. 
Under suitable regularity conditions, the \texttt{PI}, \texttt{EI}, and
\texttt{UCB} acquisition functions are guaranteed to achieve asymptotic
consistency~\cite{bull2011convergence, mockus1978application,
srinivas2009gaussian, vazquez2010convergence} in BO. For a limited budget $q$, however, 
such acquisition functions can be suboptimal~\cite{gonzalez2016glasses}.
Furthermore, the myopic nature of these basic acquisition functions is such that they seek to maximize the 
oracle value at $t=T$ only when they arrive at $T$, meaning that the past ($t<T$) decisions 
are not necessarily made to facilitate making the best decision for 
$t=T$. On the other hand, we are interested in a strategy that seeks to
maximize $f(\mathbf{x},T)$. 
Since the time horizon $t=T<\infty$, this is a \emph{finite-horizon} problem.

For time-dependent objectives, an optimal decision may be horizon-dependent,
which poses a challenge for improvement-based acquisition functions such as \texttt{PI} and \texttt{EI}. Because the target $\xi$ is typically chosen as the best observed value in $\{\hat{y}_i\}_{i=1}^n$, specifying an appropriate target when the maximum oracle value $\max f(\mathbf{x},t)$
changes with time $t$ is challenging. For example, what if the maximum oracle value at $T$ is
much smaller than anything observed until current time $t$? A further challenge in our context is that observations at the target time $T$ are not realizable until we conclude the optimization. Therefore, in such situations, myopic 
acquisition functions that seek
the best oracle value at the current $t$ may not necessarily be  useful. Instead,
we are interested in acquisition functions that  optimize the
\emph{long-sighted} decision (at $t=T$) even if it means incurring a suboptimal
oracle value at intermediate $t<T$. 
In the BO literature, strategies that make decisions by accounting for the
decisions remaining in the future are generally referred to as \emph{lookahead}
approaches, which we adopt to solve our  finite-horizon problem. 



We now illustrate the benefit of a lookahead acquisition function in the
context of time-dependent BO. \Cref{fig:EI_v_2LEY} shows the difference
between the decisions made via EI and the lookahead acquisition function
\texttt{2LEY} introduced in \cref{s:mLBO}. The
top left of the figure shows contours of the true (noise-free) oracle with a target horizon of 
$T=3.0$. We assume we have made $n$ decisions before $t_{n+1}=2.625$, and we
are interested in making a decision at current $t_{n+1}$ (denoted by the horizontal
dashed line). The top right shows the $\texttt{2LEY}$ and $\texttt{EI}$
acquisition functions and the corresponding decisions made by maximizing each
of them (denoted by vertical dashed lines). In the bottom row, we show the GP
predictions at $T$, with each of the decisions made via the myopic and
lookahead acquisition functions incorporated, along with the noise-free $f(\mathbf{x},
T)$ for reference, and the final decisions at $T$ (denoted by vertical dashed
lines). Notice that the lookahead decision results in a more accurate
prediction of the global maximum at $T$ compared with the myopic decision. This illustration shows that 
lookahead acquisition functions have the potential to make better decisions
than do myopic acquisition functions when the oracle value at a target (future) horizon is of interest. One could ask whether a randomly made decision that offers a \emph{space filling} effect (e.g., $0\leq \mathbf{x} \leq 0.4$ at $t=2.625$ ) could better learn $f(\mathbf{x}, T)$ and hence lead to a better final decision at $T$. We show that this is not the case, with a randomly chosen decision $(0.2474, 2.625)$ that does not outperform the lookahead method; see bottom right of \Cref{fig:EI_v_2LEY}.


\begin{figure}[htbp!]
    \centering
    \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{EI_v_2LEY.pdf} 
        \caption{Left: contour of $f(\mathbf{x}, t)$ in $\bm{\mathcal{X}}\times \mathcal{T}$ with $T=3$. Red circles denote the locations of $n$ previous observations. Horizontal line denotes the time $t_{n+1}=2.625$ at which the next observation is to be made. Right: myopic (\texttt{EI}) and lookahead (\texttt{2LEY}) acquisition functions. }
    \end{subfigure}\\
    \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{EI_v_2LEY_v_Random_T.pdf}   
        \caption{GP prediction at $T$ with decisions incorporated from myopic (\texttt{EI}) and lookahead (\texttt{2LEY}) decisions as well as a decision chosen uniformly at random from $\boldsymbol{\mathcal{X}}$ at $t$. Shaded region represents 95\% confidence region.}
    \end{subfigure}
    \caption{Comparison of the decisions made via a myopic (\texttt{EI}) 
    and a lookahead (\texttt{2LEY}) acquisition function. The top right shows the acquisition functions and the respective decisions (maximizers). The GP predictions, incorporating the decision made via either of the acquisition functions, are shown in the bottom. One can see that if the final decision at $T$ is made as   $\operatornamewithlimits{arg\,max}_{\mathbf{x}\in\bm{\mathcal{X}}}~\mu_{n+1}(\mathbf{x}, T)$, then the lookahead approach predicts the maximizer at $T$ more accurately than does the myopic approach. A random decision $(0.2474, 2.625)$, even though it is \emph{space filling}, does not outperform the lookahead decision.}
    \label{fig:EI_v_2LEY}
\end{figure}

In foundational work on lookahead BO, Osborne \cite{osborne2010bayesian}
showed that proper Bayesian reasoning can be used to define an acquisition
function where the ($n+1$)th observation is made by marginalizing a loss function
at the final ($q$th) observation  over all the remaining observations $(\mathbf{x}_i, \hat{y}_i)$,
$i=n+2,\ldots, q$. Others have proposed
\emph{lookahead EI} in the context of finite-budget BO. For example, Ginsbourger and Le
Riche \cite{ginsbourger2010towards} showed that the lookahead EI is a dynamic
program that chooses the EI maximizer in expectation, considering all possible
strategies of the same budget. 
A practical challenge with lookahead strategies is the computational
tractability (as demonstrated in, for example, \cite{ginsbourger2010towards,osborne2010bayesian}), which can be partially mitigated by various
approximation methods \cite{gonzalez2016glasses,lam2016bayesian,wu2019practical}. Yue and Kontar \cite{yue2020non} explored the
theoretical properties of the approximate dynamic programming approach for
multistep lookahead acquisition functions and discussed when these strategies are guaranteed
to perform better than their myopic counterparts. In this work
we draw inspiration from \cite{osborne2010bayesian} to address the
limited budget aspect of our problem. However, we propose a general framework
for lookahead BO that readily applies to the finite-horizon,
time-dependent BO; we are unaware of any previous work on this topic.

In the context of
time-dependent BO, the only works we are aware of are
\cite{nyikosabayesian} and \cite{bogunovic2016time}. 
We fundamentally differ
from their approach in the following ways: (i) their goal is to track a
time-varying optimum whereas ours is to predict it at the target horizon $T$), and (ii) they do not address
the limited budget constraint as we do. Furthermore, setting $k_t(t_i,
t_j)=(1-\epsilon)^{|t_i - t_j|/2}$ (where $\epsilon$ is the \emph{forgetting}
factor as defined in \cite{bogunovic2016time}), our method is equivalent to
\cite{bogunovic2016time} if we myopically make decisions at each $t$ via \texttt{UCB}.  



Another plausible way to solve the finite-horizon (target $t=T$) BO problem
is to
sequentially choose points that are most
\emph{informative}~\cite{cover2012elements} about the maximizer
at the target horizon; see, for example, 
\cite{hennig2012entropy, hernandez2014predictive, wang2017max}. Although such
information-theoretic approaches are known to involve an intractable form of
the information entropies, they are also highly reliant on an accurate underlying GP model~\cite{mackay1992information}. 


In this work, we extend myopic acquisition functions in order to make
nonmyopic decisions for
finite-horizon BO. Specifically, we propose a generalized
framework that specifies a \emph{value function} 
that quantifies an objective at the target horizon $T$ and that can then be used in a lookahead acquisition
function. Our main contributions are as follows. 
\begin{enumerate}
    \item We present the first strategy to solve finite-horizon BO for time-dependent 
   
    oracles.
   
    \item We propose lookahead versions of myopic acquisition functions~(\cref{s:mLBO}) that
      are applicable to both finite-budget BO and finite-horizon, time-dependent settings. 
    \item We present a practical algorithm that recursively applies a two-step
      lookahead acquisition function~(\cref{s:r2L}) 
     
   
   
   
      to make new observations time-dependent oracles.
    \item We establish the utility of our method by demonstrating it on carefully constructed synthetic test functions (\cref{ss:synthetic_1d,ss:synthetic_hd}) and a real-world quantum optimal control problem~(\cref{ss:qoc}). 
\end{enumerate}
Software that implements our lookahead framework for limited horizon, finite-budget time-dependent BO problems is available upon request and will be released publicly with the final version of this manuscript.

The remainder of the paper is organized as follows. In  \cref{s:mLBO}
we present our generalized lookahead acquisition framework for generic value
functions. In \cref{s:r2L} we present a recursive two-step lookahead
approach that is more practical than a more-than-two-step lookahead approach.
We discuss theoretical properties of our approach in \cref{s:theory}.
We then
present the results of our numerical experiments on synthetic test functions
(\cref{ss:synthetic_1d,ss:synthetic_hd}) as well as a
real-world quantum optimal control problem in \cref{ss:qoc}. We conclude the paper with an outlook on future work.

\section{Generalized m-step lookahead acquisition function}
\label{s:mLBO}
Given the data from $n$ observations $\bm{\mathcal{D}}_n$, let us define a value function $\upsilon$ by 
\begin{equation}
    \upsilon_n(\mathbf{x}, t; \boldsymbol{\omega}) := \mathbb{E}_{y\sim Y|\bm{\mathcal{D}}_n}\left[
   
    h(y(\mathbf{x}, t); \boldsymbol{\omega})     
    \right],
    \label{e:value_function}
\end{equation}
where $h$ is a scalar-valued function, $\boldsymbol{\omega} \in \mathbb{R}^p$ are
parameters that parameterize $h$, 
$Y|\bm{\mathcal{D}}_n$ is the GP posterior 
given $n$ observations, and $y(\mathbf{x}, t)$ is a sample path realized from $Y(\mathbf{x},
t)$. For example, the \texttt{PI}, \texttt{EI}, and \texttt{UCB} acquisition functions can be obtained from defining $h$ as 

\begin{equation}
    \begin{split}
   
    h(y(\mathbf{x}, t); \xi) 
    &= y(\mathbf{x}, t) - \xi \implies \upsilon(\mathbf{x},t;\xi) = \Phi \left( z \right) \qquad \text{(PI)}\\
   
    h(y(\mathbf{x}, t); \xi) 
    &= [y(\mathbf{x}, t) - \xi]^+ \implies \upsilon(\mathbf{x},t;\xi) = z\sigma(\mathbf{x}, t) \Phi\left(z \right) + \sigma_n(\mathbf{x}, t) \phi \left( z \right) \qquad \text{(EI)}\\
       
        h(y(\mathbf{x}, t); \beta)
        &= y(\mathbf{x}, t) +  \beta^{1/2}  \sigma(\mathbf{x},t) \implies \upsilon(\mathbf{x},t;\beta) = \mu(\mathbf{x}, t) + \beta^{1/2}\sigma(\mathbf{x}, t), \qquad \text{(UCB)} 
        \end{split}
    \label{e:value_functions}
\end{equation}
where recall that $\xi$ and $\beta$ are user-specified target and confidence parameters, respectively.
Similarly, when $h$ is the identity function $h(y(\mathbf{x}, t); \boldsymbol{\omega})  = y(\mathbf{x}, t)$, then $\upsilon(\mathbf{x},t)= \mu(\mathbf{x}, t)$. 

We use the following notation for conciseness:
\begin{align*}
&\mathbf{x}^t_{n+1} := (\mathbf{x}, t_{n+1}) && \bm{\mathcal{D}}_{n,j} := \bm{\mathcal{D}}_n \bigcup_{i=1,\ldots,j} \{ \mathbf{x}^t_{n+i}, y_{n+i} \}\\ 
&y_{n+1} := y(\mathbf{x}^t_{n+1}) &&\mathbb{E}_{n,j} := \mathbb{E}_{y\sim Y|\bm{\mathcal{D}}_{n,j}} \\
&\mu_{n,j} := \mathbb{E}_{n,j}Y, ~~\sigma^2_{n,j} := \mathbb{E}_{n,j} \left( Y - \mu_{n,j}\right)^2 && \upsilon^T_{n,j} := \upsilon_{n,j}(\mathbf{x}, T; \boldsymbol{\omega}) = \mathbb{E}_{n,j}~\left[h(y(\mathbf{x}, T); \boldsymbol{\omega}) \right],
\end{align*}
where we note
that whenever we explicitly use $y(\mathbf{x},t)$, we refer to a sample path drawn from
$Y(\mathbf{x},t)$ and therefore $y_{n+1} = y(\mathbf{x}_{n+1},t_{n+1})$ is a draw from $y(\mathbf{x},t)$ at $(\mathbf{x}_{n+1}, t_{n+1})$. 
On the other hand, when we use $\hat{y}_i$, we refer to an observation of the oracle at
$(\mathbf{x}_i,t_i)$. Similarly, note that $\bm{\mathcal{D}}_{n,k}$ contains $n$ oracle evaluations 
and $k$ candidate points $\mathbf{x}$ and the $k$ corresponding draws from the GP at times 
$t_{n+1}$, \ldots, $t_{n+k}$. Furthermore, we note here that $\mathbf{x}^t_{n+1}$ is $d$-dimensional.

We define our $m$-step lookahead acquisition function as
\begin{equation}
    \begin{split}
      \alpha_m(\mathbf{x}^t_{n+1}) = \int \ldots \int  & \left[\max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{m-1}(\tilde{\mathbf{x}}, T; \boldsymbol{\omega})|\bm{\mathcal{D}}_{n,m-1}\right] ~ p(y_{n+m-2}|\bm{\mathcal{D}}_{n,m-2})\\ &p(\mathbf{x}^t_{n+m-1}|\bm{\mathcal{D}}_{n,m-2})\cdots 
    p(y_{n}|\bm{\mathcal{D}}_n)~dy_{n+m-2}\cdots dy_{n}~d\mathbf{x}^t_{n,m-1} \cdots d\mathbf{x}^t_{n+2}.
    \end{split}
    \label{e:mstep}
\end{equation}
Note that in \eqref{e:mstep}, the \emph{inner maximization} is always with respect to $\mathbf{x}$ at $t=T$ and that after evaluating the integral, the resulting expression is a function of $\mathbf{x}^t_{n+1}$. Furthermore, since we are interested in choosing the maximizer of $\alpha_m(\mathbf{x}^t_{n+1})$ for $t_{n+1}$, we do not marginalize $\mathbf{x}^t_{n+1}$. Computing \eqref{e:mstep} involves the computation of $m-1$ nested expectations. Additionally, when the expectations are approximated via Monte Carlo sampling, there is a total of $N\times (m-1)$ $d$-dimensional  optimization solves due to the inner maximization in \eqref{e:mstep}. The computation or even the Monte Carlo approximation of \eqref{e:mstep} becomes intractable as $m$ increases, warranting further approximation, as is done, for instance, in \cite{gonzalez2016glasses, lam2016bayesian, yue2020non}. When $m=2$, however, the computation of \eqref{e:mstep} is much more straightforward, as we discuss in the following section.









\section{Generalized two-step lookahead acquisition function}
\label{s:r2L}
In this work, as a practical alternative to the computationally intractable multistep ($m>2$) 
lookahead acquisition function, we recursively optimize the two-step lookahead acquisition function~\eqref{e:E_twostep} at each time step of the time-dependent BO problem; this process is graphically described in \cref{f:method}.
Given $\bm{\mathcal{D}}_n$ and with $m=2$, we define our \emph{two-step} lookahead acquisition function by
\begin{equation}
\begin{array}{rl}
 
 


  \alpha_2(\mathbf{x},t_{n+1})
  
&=  \displaystyle \int \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~ \left[ \upsilon_{n,1}(\tilde{\mathbf{x}}, T; \boldsymbol{\omega}) | \, 
\bm{\mathcal{D}}_n \bigcup \{ (\mathbf{x},t_{n+1}), y_{n+1} \}
\right]~ p(y_{n}|\bm{\mathcal{D}}_{n}) ~dy_{n}, 
\end{array}
\label{e:twostep}
\end{equation}
where $p(y_{n}|\bm{\mathcal{D}}_{n})$ is the probability density of the GP posterior distribution
conditioned on $\bm{\mathcal{D}}_n$.
We write \eqref{e:twostep} concisely as 
\begin{equation}
  \alpha_2(\mathbf{x}^t_{n+1}) = \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}(\tilde{\mathbf{x}},T;\boldsymbol{\omega})|\bm{\mathcal{D}}_{n,1} \right].
    \label{e:E_twostep}
\end{equation}
In \eqref{e:E_twostep}, recall that $\bm{\mathcal{D}}_{n,1}$ is the union of $\bm{\mathcal{D}}_n$ and the $n+1$th observation being a draw from the GP $Y_n$, and we compute $\upsilon_{n,1}(\tilde{\mathbf{x}}, T; \boldsymbol{\omega})$ by taking the expectation with respect to $Y|\bm{\mathcal{D}}_{n,1}$. In general, it is not guaranteed that
\eqref{e:E_twostep} admits a closed-form expression. Hence we make the
Monte Carlo approximation
\begin{equation}
    \alpha_2(\mathbf{x}^t_{n+1}) \approx \hat{\alpha}_2(\mathbf{x}^t_{n+1}) := \frac{1}{N} \sum_{j=1}^N \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}^{T,j},
    \label{e:twostep_mc}
\end{equation}
where $\upsilona{n}{1}^{T,j} = \upsilona{n}{1}(\tilde{\mathbf{x}}, T)|y^j_{n+1}$ with the $N$ iid samples $\{y^j_{n+1}: \,j=1,\ldots,N\}$ drawn from $Y_{n}$. In practice, the computation of \eqref{e:E_twostep} involves the following steps (given $\mathbf{x}_{n+1}$):
\begin{enumerate}
    \item Generate $N$ iid samples $\{ y^j_{n+1}\}_{i=1}^N$ drawing from $\mathcal{N}(\mu_n(\mathbf{x}_{n+1}, t_{n+1}), \sigma_n^2(\mathbf{x}_{n+1}, t_{n+1}))$.
    \item For each $y^j_{n+1},~j=1,\ldots,N$
    \begin{enumerate}
        \item[2.1] update $\bm{\mathcal{D}}_n$ to $\bm{\mathcal{D}}_{n,1}$ to get GP posterior $Y|\bm{\mathcal{D}}_{n,1}$ and
        \item[2.2] maximize a draw of the projected value function $\nu^j = \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}(\tilde{\mathbf{x}}, T)|\bm{\mathcal{D}}_{n,1}$.
    \end{enumerate}
    \item Compute $\hat{\alpha}_2(\mathbf{x}_{n+1}) \approx \frac{1}{N} \sum_{j=1}^N \nu^j$.
\end{enumerate}

In order to make a decision at $t_{n+1}$, \eqref{e:E_twostep} needs to be maximized. We show how the maximization of the acquisition function can be done efficiently by constructing unbiased estimators of the acquisition function's gradient. Under suitable regularity conditions, we can interchange the expectation and gradient
operators (see \cref{thm:gradexp}
) to obtain
\begin{equation}
    \nabla_{\mathbf{x}}\alpha_{2}(\mathbf{x}^t_{n+1}) = \mathbb{E}_{n} \left[ \nabla_{\mathbf{x}}~ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}(\tilde{\mathbf{x}}, T) | \bm{\mathcal{D}}_{n,1} \right].
    \label{e:twostep_T_grad}
\end{equation}
Similar to the Monte Carlo approximation for the acquisition function presented
in \eqref{e:twostep_mc}, \eqref{e:twostep_T_grad} can be approximated as
\begin{equation}
    \nabla_{\mathbf{x}}\alpha_{2}(\mathbf{x}^t_{n+1}) \approx \tilde{\nabla}_{\mathbf{x}}\alpha_{2}(\mathbf{x}^t_{n+1}) = \frac{1}{N} \sum_{i=1}^{N} d\nu^i :=  \frac{1}{N} \sum_{i=1}^{N}  \nabla_{\mathbf{x}}~ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}^{T,i}.
    \label{e:twostep_T_grad_mc}
\end{equation}
In this way, sample-averaged gradients can be used in a gradient-based optimizer
to efficiently optimize \eqref{e:E_twostep}. We clarify the issues associated
with computing the gradients of functions involving the $\max()$ operator as in
\eqref{e:twostep_T_grad_mc} in \Cref{s:theory}. 

We now illustrate in detail the two-step
acquisition function for a specific value function.
\subsection*{Two-step lookahead expected oracle value 
(\texttt{2LEY}) acquisition function}
\label{ss:r2LEY}
\begin{figure}
    \centering
    \includegraphics[width=\textwidth, trim=0.5cm 4cm 0.5cm 6cm,clip]{Sequential 2-step Lookahead Strategy.pdf}
    \caption{Recursive two-step lookahead method for time-dependent problems.
    Here, we start with $n$ observations up to $t_n$ and have $q-n$ 
   
    remaining
    observations to make, where $q$ is the total budget. We sequentially select
    one point $\mathbf{x}_{n+i}$ at each step that maximizes the expected value
    function evaluated at $T$; this is denoted by the black dashed lines
    at $t_{n+i}$ and $T$, whereas the gray dashed lines denote the
    remaining decisions to be made. We repeat the same process until the
    decision is made at the penultimate step $i=q-1$. The final decision at $T$
    is made by maximizing the value function at $T$.}
    \label{f:method}
\end{figure}
As an illustration, we set $h$ to be the identity function, and hence $\upsilona{n}{1}(\mathbf{x}, T) = \mu_{n,1}(\mathbf{x},
T)|\bm{\mathcal{D}}_{n,1}$.
We define our two-step lookahead expected oracle value (\texttt{2LEY}) acquisition function as
\begin{equation}
    \begin{split}
    \alpha_{2LEY}(\mathbf{x}^t_{n+1}) := & \mathbb{E}_n \left[  \underset{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}{\max}~\left.\mu_{n,1}(\tilde{\mathbf{x}}, T) \right| \bm{\mathcal{D}}_{n,1} \right] \\
        =& \mathbb{E}_n \left[ \mu^*_{n,1}( \mathbf{x}^t_{n+1}, T) \right],
    \end{split}
    \label{e:2LEY_1}
\end{equation}
where $\mu^*_{n,1} = \mu_{n,1}(\mathbf{x}^*, T)$ with $\mathbf{x}^*$ being a maximizer of $\mu_{n,1}(\mathbf{x},
T)|\bm{\mathcal{D}}_{n,1}$ and 
$y_{n+1}$ is a draw from $\mathcal{N} \left(\mu_n(\mathbf{x}^t_{n+1}), \sigma_n^2(\mathbf{x}^t_{n+1}) \right)$. 
The dependence of the second line of \eqref{e:2LEY_1} on
$\mathbf{x}^t_{n+1}$ can
be seen by realizing that
\begin{equation}
    \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T)  = \mathbf{k}_{n+1}^\top \mathbf{K}^{-1}_{n+1} [\mathbf{y}_n^\top, y_{n+1}]^\top,
\end{equation}
where $ \mathbf{k}_{n+1} = [ k\left((\mathbf{x}^*,T), (\mathbf{x}_1, t_1)\right), \ldots,
k\left((\mathbf{x}^*,T), (\mathbf{x}_n, t_n)\right), k\left((\mathbf{x}^*,T),
(\mathbf{x}_{n+1},t_{n+1})\right)]^\top$, and $\mathbf{K}_{n+1, ij} = k\left((\mathbf{x}_i, t_i), (\mathbf{x}_j, t_j) \right);~\forall i=1,\ldots,n+1, ~ j=1,\ldots,n+1$. The expectation in \eqref{e:2LEY_1} is approximated as
\begin{equation}
    \alpha_{2LEY}(\mathbf{x}^t_{n+1}) = \mathbb{E}_n \left[ \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right] \approx \frac{1}{N} \sum_{j=1}^N  \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) | y^j_{n+1}.
    \label{e:2LEY_MC}
\end{equation}
The gradient of $\mu_{n,1}^*(\mathbf{x}^t_{n+1}, T)$ with respect to $\mathbf{x}^t_{n+1}$ is given by
\begin{equation}
    \begin{split}
        \nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) = &\frac{\partial \mathbf{k}_{n+1}^\top}{\partial \mathbf{x}_{n+1}} \mathbf{K}^{-1}_{n+1} + \mathbf{k}_{n+1}^\top \frac{\partial \mathbf{K}^{-1}_{n+1}}{\partial \mathbf{x}_{n+1}} \\
        = &\frac{\partial \mathbf{k}_{n+1}^\top}{\partial \mathbf{x}_{n+1}} \mathbf{K}^{-1}_{n+1} +
        \mathbf{k}_{n+1}^\top~\mathbf{K}^{-1}_{n+1} \frac{\partial \mathbf{K}_{n+1}}{\partial \mathbf{x}_{n+1}}\mathbf{K}^{-1}_{n+1},
    \end{split}
    \label{e:2LEY-grad_a}
\end{equation}
where $ \frac{\partial \mathbf{K}_{n+1}}{\partial \mathbf{x}_{n+1}}$ is a matrix of
elementwise derivatives and the second line follows from a well-known lemma on the derivative of a matrix inverse~\cite[pp.~201--202]{rasmussen:williams:2006}.

With a continuously differentiable kernel $k(\cdot, \cdot)$, we have that
\begin{equation}
    \nabla_{\mathbf{x}} \alpha_{2LEY}(\mathbf{x}^t_{n+1}) = \nabla_{\mathbf{x}} \mathbb{E}_n \left[ \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right] = \mathbb{E}_n \left[\nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right],
    \label{e:2LEY_gr_ex_interchange}
\end{equation}
where the interchange of the gradient and expectation operators is via \cref{thm:gradexp}. Then the gradient is approximated as 
\begin{equation}
    \nabla_{\mathbf{x}} \alpha_{2LEY}(\mathbf{x}^t_{n+1}) = \mathbb{E}_n \left[\nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right] \approx \frac{1}{N} \sum_{j=1}^N \nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) | y^j_{n+1}.
    \label{e:2LEY_grad_MC}
\end{equation}
Furthermore, \eqref{e:2LEY_MC} and
\eqref{e:2LEY_grad_MC} are unbiased estimators for $\alpha_{2LEY}$ and $\nabla_{\mathbf{x}}
\alpha_{2LEY}$, respectively; see \cref{thm:gradexp}. This gradient can then be used in a gradient-based optimizer to
maximize our two-step lookahead acquisition function in an efficient manner.

Similarly, the two-step lookahead  extensions for
\texttt{EI}, \texttt{PI}, and \texttt{UCB} are obtained by picking the
corresponding value functions as $\alpha_{\textrm{EI}}(\mathbf{x},T)$, $\alpha_{\textrm{PI}}(\mathbf{x},T)$ and $\alpha_{\textrm{UCB}}(\mathbf{x},T)$ and are
named respectively, $\alpha_{2LEI}$, $\alpha_{2LPI}$, and $\alpha_{2LUCB}$. Our
overall algorithm recursively maximizes the two-step acquisition function until
the penultimate step and maximizes the actual value function at the final
step. The overall algorithm is graphically depicted in \cref{f:method} and is presented in
Algorithm~\ref{a:r2LBO}. We refer to Algorithm~\ref{a:r2LBO} with the value functions chosen as $\mu_{n,1}(\mathbf{x}, T)$, $\alpha_{EI}$, $\alpha_{PI}$, and $\alpha_{UCB}$ as \texttt{r2LEY}, \texttt{r2LEI}, \texttt{r2LPI}, and \texttt{r2LUCB}, respectively.


\begin{algorithm2e}[H]
\SetAlgoLined
\textbf{Given:} $\bm{\mathcal{D}}_n$,  
total budget $q$, 
schedule $\lbrace t_{n+1},\ldots, t_{q}:=  T \rbrace$, and GP hyperparameters $\bm{\Omega}$ \\
\KwResult{$(\mathbf{x}_T, \hat{y}_T)$}
  \For{$i=n+1, \ldots, q-1$, }{
  Define value function $\upsilon^T_{i-1, 1}|\bm{\mathcal{D}}_{i-1, 1}$ \\
  Find $\mathbf{x}_i = \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\operatornamewithlimits{arg\,max}}~ \hat{\alpha}_2(\mathbf{x}, t_i)$, \qquad (acquisition function maximization)\\
  \quad with $\hat{\alpha}_2$ and $\nabla_{\mathbf{x}}~\hat{\alpha}_2$ evaluated via \eqref{e:twostep_mc} and \eqref{e:twostep_T_grad_mc} resp. \\
    Observe $\hat{y}_i$ = $f(\mathbf{x}_i, t_i) + \epsilon_i$\\ 
    Append $\bm{\mathcal{D}}_i = \bm{\mathcal{D}}_{i-1} \cup \lbrace (\mathbf{x}_i, t_i), \hat{y}_i \rbrace$\\ 
    Update GP hyperparameters $\bm{\Omega}$ \\
 }
 Find $\mathbf{x}_T \in \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\operatornamewithlimits{arg\,max}}~\upsilon_{q}(\mathbf{x}, T)|\bm{\mathcal{D}}_{q-1}$ \qquad (value function maximization)\\
 Observe $\hat{y}_T$ = $f(\mathbf{x}_T, T) + \epsilon_T$\\ 
 \caption{Recursive two-step lookahead Bayesian optimization} 
 \label{a:r2LBO}
\end{algorithm2e}











\section{Theoretical properties}
\label{s:theory}

We now present the theoretical foundation of our method and a few properties.  
  Let $\mathcal{Y}$
  be the support of $y_{n+1}$, whose density is $p(y_{n+1})$. As before, let the
  domain of $\mathbf{x}^t_{n+1}$ be $\bm{\mathcal{X}}$. 





We make the following remark about the differentiability of our acquisition function; this justifies the interchange of the gradient and expectation operators.

\begin{remark}
We write the inner term of our two-step lookahead acquisition function in \eqref{e:E_twostep} as $\max_{\mathbf{x}\in\bm{\mathcal{X}}} \mathbb{E}_{n,1}h(y)$, where we omit the arguments of $y$ and $\boldsymbol{\omega}$ for the sake of conciseness. The GP posterior mean $\mu(\mathbf{x}, t)$ and standard deviation $\sigma(\mathbf{x}, t)$ are continuously differentiable with the sufficient condition that we choose a continuously differentiable prior covariance kernel $k$~\cite{smith1995differentiation, wang2016parallel}. We then use the reparameterization trick~\cite{wilson2017reparameterization,
wilson2018maximizing} to write the GP sample path as $y(\mathbf{x}, T) := \mu(\mathbf{x}, T) + \sigma(\mathbf{x},
T)\times \gamma$, where $\gamma \sim \mathcal{N}(0, 1)$, which makes the differentiability of our acquisition function transparent. 


From a purely implementation standpoint, we use the gradient $\nabla
\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_{n, 1}~h(y)}$ if it exists and a
\emph{subgradient}~\cite{boyd2004convex} otherwise. In our implementation, the
gradients are computed by backpropagating the output of the
$\max()$ operation; in \Cref{app:grad_nd} we provide examples of how (sub)gradients are calculated for nonsmooth functions.
\end{remark}








Let $\mathbf{x}^*$ be a maximizer of $\mathbb{E}_{n,1} h(y(\mathbf{x}, T); \boldsymbol{\omega}) = \upsilon_{n,1}(\mathbf{x}, T)|\bm{\mathcal{D}}_{n,1}$. In what follows, we write
  $g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) := \upsilon_{n,1}(\mathbf{x}^*, T)|\bm{\mathcal{D}}_{n,1}$ to denote the maximum value of the value function at $T$ after having updated $\bm{\mathcal{D}}_n$ with $y_{n+1}$, and $\nabla_{\mathbf{x}}g^* := \partial g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})/\partial \mathbf{x}^t_{n+1} $. Note again that the dependence of $g^*$ on $\mathbf{x}^t_{n+1}$ is due to $\bm{\mathcal{D}}_{n,1}$.
  
\begin{theorem}[Interchange of gradient and expectation operators] 
\label{thm:gradexp}
 Let $\nabla_{\mathbf{x}}g^*$ exist, and let $g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})\times
  p(y_{n+1})$ and $\nabla_{\mathbf{x}}g^* \times
  p(y_{n+1})$ be continuous on $\bm{\mathcal{X}} \times \mathcal{Y}$. 
  Further assume
  that there exist functions $q_0$ and $q_1$ such that
\begin{equation*}
\begin{array}{rl}
  |g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) \times p(y_{n+1})| &\leq q_0(y_{n+1}) \\
  \|
  \nabla_{\mathbf{x}} g^* \times  p(y_{n+1})\| &\leq q_1(y_{n+1})
  \end{array}
  \qquad \forall (\mathbf{x}^t_{n+1}, y_{n+1}) \in \bm{\mathcal{X}} \times \mathcal{Y}, 
  \end{equation*}
where
  $\int_{\mathcal{Y}}q_1(y_{n+1}) dy_{n+1} < \infty$ and
  $\int_{\mathcal{Y}}q_2(y_{n+1}) dy_{n+1} < \infty$. Then,
\[ \nabla_{\mathbf{x}} \mathbb{E}_j \left[ g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) \right] = \mathbb{E}_n \left[\nabla_{\mathbf{x}} g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) \right].\]

Furthermore, a realization of $\nabla_{\mathbf{x}} g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})$ yields
an unbiased estimate of the gradient of $\alpha_2(\mathbf{x}_{n+1}^t)$ in \eqref{e:E_twostep}.
\end{theorem}

\begin{proof}
Using $g^*(\mathbf{x}^t_{n+1})$ to denote  $g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})$ for the sake of brevity, we have
\begin{equation}
\begin{split}
    \nabla_\mathbf{x} \alpha_2(\mathbf{x}^t_{n+1}) := \frac{\partial}{\partial \mathbf{x}^t_{n+1}} \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}) \right] &= \underset{h\rightarrow 0}{\textnormal{lim}}~\frac{1}{h} \left \lbrace \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}+h) \right] -  \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}) \right] \right \rbrace \\
    &= \underset{h\rightarrow 0}{\textnormal{lim}}~ \left \lbrace \mathbb{E}_n \frac{1}{h} \left[ g^*(\mathbf{x}^t_{n+1}+h) - g^*(\mathbf{x}^t_{n+1}) \right] \right \rbrace \\
    &= \underset{h\rightarrow 0}{\textnormal{lim}}~ \left \lbrace \mathbb{E}_n \frac{\partial g^*(\bar{\mathbf{x}}^h_{n+1})}{\partial \mathbf{x}_{n+1}}  \right \rbrace ,\\
\end{split} 
\label{e:grad_limit}
\end{equation}
where $\bar{\mathbf{x}}^h_{n+1} = \lambda \mathbf{x}^t_{n+1} + (1-\lambda)(\mathbf{x}^t_{n+1} + h)$ for some
$\lambda \in [0,1]$ and the last line follows from the first-order mean value
theorem.
We bring the limit inside the integral by Lebesgue's dominated convergence theorem to obtain
\begin{equation}
    \frac{\partial}{\partial \mathbf{x}^t_{n+1}} \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}) \right] =  \mathbb{E}_n~\underset{h\rightarrow 0}{\textnormal{lim}} \frac{\partial g^*(\bar{\mathbf{x}}^h_{n+1})}{\partial \mathbf{x}^t_{n+1}}   = \mathbb{E}_n \frac{\partial g^*(\mathbf{x}^t_{n+1})}{\partial \mathbf{x}^t_{n+1}}.
    \label{e:LDC}
\end{equation}
Since \eqref{e:grad_limit} holds, a realization $\partial g^*(\mathbf{x}^t_{n+1}) / \partial \mathbf{x}^t_{n+1}$ yields an unbiased estimate of the gradient $\nabla_\mathbf{x} \alpha_2(\mathbf{x}^t_{n+1})$.
\end{proof}



\begin{proposition}[Special case when $m=2$ and $\upsilon^T = \mu(\mathbf{x}, T)$] 
When $m=2$, our \texttt{2LEY} acquisition function is equivalent to the knowledge gradient
(KG)~\cite{frazier2008knowledge, wu2016parallel} acquisition function offset by a constant factor.

\begin{proof}
The KG acquisition function, in the context of our problem, is given as
\begin{equation}
    \alpha_{KG}(\mathbf{x}^t_{n+1}) = \mathbb{E}_{n}\left[ \underset{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}{\max}~ \mu_{n,1}(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] - \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\max}~ \mu_{n}(\mathbf{x},T).
    \label{e:KG_time}
\end{equation}

With $m=2$, the \texttt{2LEY} acquisition function at $\mathbf{x}_{n+1}$ is defined as follows:
\begin{equation}\label{eq:nm1}
    \begin{split}
    \alpha_{2LEY}(\mathbf{x}^t_{n+1}) & = \int_{y_{n}} \left[ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}} ~\mathbb{E}_{n,1} ~y(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] p(y_{n}|\bm{\mathcal{D}}_{n})~dy_{n} \\
    & = \int_{y_{n}}\left[  \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~ \mu_{n,1}(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] p(y_{n}|\bm{\mathcal{D}}_{n})~dy_{n} \\
    & = \mathbb{E}_{n}\left[ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~ \mu_{n,1}(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] \\
    & = \alpha_{KG}(\mathbf{x}^t_{n+1}) + \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\mu_{n}(\tilde{\mathbf{x}},T), 
    \end{split}
\end{equation}
where the last line follows from \cref{e:KG_time}. Note that $\alpha_{2LEY}$ and $\alpha_{KG}$ share the same maximizer. 
\end{proof}
\label{prop:special_cases}
\end{proposition}

\begin{proposition}[Special case when $|T-t| \gg\theta_t$] 
With the prior assumption $f(\mathbf{x},t)\sim \mathcal{GP}(0, k(\cdot,\cdot))$, where $k(\cdot,\cdot) := k_{\mathbf{x}}(\cdot,\cdot; ~\theta_{\mathbf{x}})\times k_t(\cdot,\cdot; ~\theta_{t})$, with $k_{\mathbf{x}}$ and $k_t$ being squared-exponential kernels,
our two-step lookahead acquisition function approaches $\alpha_2(\mathbf{x}^t_{n+1}) = \max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\upsilon_n(\tilde{\mathbf{x}}, T)}$ as $|T - t_{n+1}| \rightarrow \infty$. In other words, far away from the target horizon $T$, $\alpha_2(\mathbf{x}^t_{n+1})$ approaches a constant function.


\begin{proof}
When $|T-t_{n+1}| \gg \theta_t$, then $k_t(t_{n+1}, T) \approx 0$; this yields $k\left( (\mathbf{x},t_{n+1}), (\mathbf{x},T) \right) \approx 0$. Let $\mu_n$ and $\sigma_n^2$ be the posterior mean and variance given $\bm{\mathcal{D}}_n$, respectively. Given another observation 
$y_{n+1}$ drawn from the GP $Y_n$ at $(\mathbf{x}, t_{n+1})$, where $|T-t_{n+1}| \gg \theta_t$, the updated posterior mean and variance are given by
\begin{equation}
    \begin{split}
        \mu_{n,1} =& \mathbf{k}^\top_{n+1}\mathbf{K}^{-1}_{n+1}[\mathbf{y}^\top_n, y_{n+1}]^\top  \approx [\mathbf{k}^\top_n, 0] \begin{bmatrix}
        \mathbf{K}_n & \mathbf{0} \\
        \mathbf{0} & 1
        \end{bmatrix}^{-1}[\mathbf{y}^\top_n, y_{n+1}]^\top = \mu_n \\
        \sigma^2_{n,1} =& 1-\mathbf{k}^\top_{n+1}\mathbf{K}^{-1}_{n+1}\mathbf{k}_{n+1} \approx 1 - [\mathbf{k}^\top_n, 0] \begin{bmatrix}
        \mathbf{K}_n & \mathbf{0} \\
        \mathbf{0} & 1
        \end{bmatrix}^{-1}[\mathbf{k}^\top_n, 0]^\top = \sigma_n^2,
    \end{split}
\end{equation}
where we have assumed, without loss of generality, that $k((\mathbf{x},t), (\mathbf{x},t))=1$ and have used the matrix inversion lemma.
Therefore, writing our acquisition function using the reparameterization trick, we see that
\begin{equation}
    \begin{split}
    \alpha_2(\mathbf{x}^t_{n+1}) =& \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_{n, 1}~h(y(\tilde{\mathbf{x}}, T))} \right]\\ =& \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\mu_{n,1}(\tilde{\mathbf{x}}, T) + \sigma_{n,1}(\tilde{\mathbf{x}}, T)\gamma)} \right]\\ \approx& 
    \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\mu_{n}(\tilde{\mathbf{x}}, T) + \sigma_{n}(\tilde{\mathbf{x}}, T)\gamma)} \right]\\
    =& \max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\upsilon_n(\tilde{\mathbf{x}}, T)},
    \end{split}
\end{equation}
and thus, in the limit $|T - t_{n+1}| \rightarrow \infty$,
$\alpha_2(\mathbf{x}^t_{n+1}) = \max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\upsilon_n(\tilde{\mathbf{x}}, T)}$.
\end{proof}
\label{prop:t_llt_T}


\begin{remark}
    Following \Cref{prop:t_llt_T}, one can  see that if $k_t(t_{n+1}, T) \approx 0$, then $k_t(t_i, T) \approx 0,~\forall i=1,\ldots, n$ since $t_{i} > t_{i-1}$. This further leads to $\mathbf{k}_{n+1} \approx \mathbf{0}$ and hence $\mu_{n,1} \approx 0$ and $\sigma_{n,1}^2 \approx 1$. This results in (for $|T-t_{n+1}| \gg \theta_t$),
    \begin{equation}
    \begin{split}
    \alpha_{2}(\mathbf{x}^t_{n+1}) =& \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\mu_{n}(\tilde{\mathbf{x}}, T) + \sigma_{n}(\tilde{\mathbf{x}}, T)\gamma)} \right] \\
    \approx & \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\gamma)} \right] \\
    =& c,
    \end{split}
    \label{e:2_inner}
\end{equation}
where $c$ is some constant. Suppose the $(n+1)$th decision is made by maximizing $\alpha_{2}(\mathbf{x}^t_{n+1})$ with a multistart local optimization procedure started from a set of starting points $\{\mathbf{x}\} \in \bm{\mathcal{X}}$, where $\{\mathbf{x}\}$ are sampled from distribution $p_\mathbf{x}$. Then, in the limit $|T - t_{n+1}| \rightarrow \infty$, decisions made via our lookahead acquisition function are equivalent to samples drawn randomly from $p_\mathbf{x}$.\footnote{Note that this assumes that the local optimizer seeks a critical point and returns the starting point as the maximizer.} 

We take advantage of this property and set $p_\mathbf{x} = \mathcal{U}(\bm{\mathcal{X}})$, where $\mathcal{U}$ is the uniform distribution. In practice, this could result in a decisions that \emph{explore} $\bm{\mathcal{X}}$ when the $T$ is far away and hence can lead to improved learning of $f(\mathbf{x},t)$, thus facilitating a more accurate prediction of the maximum at $T$.
\end{remark}

\label{prop:special_cases2}
\end{proposition}





\section{Numerical experiments}
We demonstrate our proposed approach on synthetic test cases and a quantum optimal control problem. 

We compare our method with
widely used myopic approaches in BO, namely,  \texttt{EI}, \texttt{PI}, and \texttt{UCB}.
We refer to ``\texttt{EImumax}'' and ``\texttt{PImumax}'' as the \texttt{EI} and \texttt{PI} strategies,
respectively, with the target $\xi$ set as the maximum of the GP posterior mean at the
current step, that is, $\underset{\mathbf{x} \in \bm{\mathcal{X}}}{\max}~ \mu(\mathbf{x}, t)$. We set 
the confidence parameter for \texttt{UCB}  to $\beta=2$. Additionally, we
 compare our method with ``\texttt{mumax},'' which selects
points by maximizing the GP posterior mean at the current time. We also make direct comparisons
between the myopic acquisition functions and their lookahead counterparts.
Moreover, we compare our method with a strategy that selects points uniformly at random
from $\mathbf{x} \sim \mathcal{U}(\bm{\mathcal{X}})$; we call this acquisition function ``\texttt{Random}.'' 

Our
metric for comparison is the oracle value $f(\mathbf{x}_T, T)$ at the target time $T$. Each query to $f(\mathbf{x}, t)$ includes additive noise with mean zero and variance $\sigma_\epsilon^2 = 10^{-3}$. Each repetition of each experiment is provided with \emph{starting} samples that include oracle evaluations at a set of points distributed with uniform spacing in $[0, t_{\textnormal{start}}]$ with each point sampled uniformly at random in $\bm{\mathcal{X}}$, where $t_\textnormal{start}<T$ is arbitrarily set for each experiment. Within each repetition the same starting samples are given to all acquisition functions. 


\subsection{Synthetic one-dimensional test cases}
\label{ss:synthetic_1d}
Our synthetic test functions take the
form $f(\mathbf{x}, t) = f_{\mathbf{x}}(\mathbf{x}) + f_{\mathbf{x} t}(\mathbf{x}, t)$ or  $f(\mathbf{x}, t) = f_{\textnormal{ref}}(\mathbf{x} (t)) $ 
in order to  induce time dependence to the
maximizer of $f(\mathbf{x},t)$.

We first consider quadratic one-dimensional functions $f_{\mathbf{x}}$
over the domain $\bm{\mathcal{X}} = [0,1]$ with
\begin{equation} 
f_{\mathbf{x}}(\mathbf{x}) = -4 (\mathbf{x} - 0.5)^2.
\label{eq:1dq}
\end{equation}
The time-dependent component takes the following four forms:
\begin{equation*}
    \begin{split}
    \textnormal{Quadratic-a}: \quad f_{\mathbf{x} t} = &\sin \left(\pi(\mathbf{x}+t)\right) + \cos \left(\pi(\mathbf{x}+t)\right) \\
    \textnormal{Quadratic-b}: \quad f_{\mathbf{x} t} = &\sin \left(\pi(\mathbf{x} t)\right) + \textnormal{cos}\left(\pi(\mathbf{x} t)\right) \\
    \textnormal{Quadratic-c}: \quad f_{\mathbf{x} t} = &\sin \left(\pi(\mathbf{x}[t-3]^+)\right) +
    \textnormal{cos}\left(\pi(\mathbf{x}[t-3]^+)\right)\\
    \textnormal{Quadratic-d}: \quad f_{\mathbf{x} t} = &2\mathbf{x}\textnormal{sin}(t) - \textnormal{sin}^2(t),  
    \end{split}
\end{equation*}
where $[t-3]^+ \equiv \max(0, t-3)$ is specified to induce movement of
$\mathbf{x}_t^*$ for $t\geq 3$. 
\cref{f:1d_synthetic} shows the trajectories of $\mathbf{x}_t^*$, the maximizer at time $t$, for each of these test cases.
\begin{figure}
    \centering
    \includegraphics[trim={100 0 100 0}, clip, width=1\linewidth]{quadratics.pdf}
    \caption{One-dimensional synthetic test cases. The dashed line represents the trajectory of $\mathbf{x}_t^*$.}
   
    \label{f:1d_synthetic}
\end{figure}
These test cases are deliberately designed to reflect a variety of situations, such as discontinuous change of $\mathbf{x}^*_t$ (Quadratic-a and Quadratic-b) and sudden dynamics (Quadratic-c). 

\cref{f:quadratic_exp} compares the points selected by our lookahead approach with various value functions and the corresponding myopic approaches. In each plot
20 repetitions 
are overlaid, with samples in $\bm{\mathcal{X}}\times[0, 1]$ used to start all of the methods;  white
circles denote the points chosen by each method at times $t>1$. 
The lookahead BO does not try to track a moving maximizer, which results in points being more spread out in a space-filling fashion than with myopic acquisition functions. 
This approach is particularly  beneficial in handling oracles whose maximizers can go through a discontinuous change (e.g., Quadratic-a and Quadratic-b).
In the case of Quadratic-c, we demonstrate (see further details in \cref{sf:c}) that our algorithm is able
to handle oracles where $\mathbf{x}_t^*$ may undergo sudden dynamics ($\approx t=3.0$ in
this case). 
For Quadratic-d, where the challenges of the other
one-dimensional cases do not exist, the algorithm performs much better in
predicting $\mathbf{x}_T^*$.
\begin{figure}
    \centering
    \begin{subfigure}{.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip, width=1\linewidth]{qa_mumax_2LEY.pdf}
        \caption{Quadratic-a}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip,width=1\linewidth]{qb_EI_2LEI.pdf}
        \caption{Quadratic-b}
    \end{subfigure} \\
    \begin{subfigure}{.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip,width=1\linewidth]{qc_PI_2LPI.pdf}
        \caption{Quadratic-c}
    \end{subfigure}%
    \begin{subfigure}{.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip,width=1\linewidth]{qd_UCB_2LUCB.pdf}
        \caption{Quadratic-d}
    \end{subfigure}    
    \caption{Comparison of our two-step lookahead (right) and myopic (left) acquisition functions on the one-dimensional synthetic test cases. The plot shows
    an overlay of 20 repetitions of the algorithm with (shared) random starting points. The
    white circles denote the decision made at each step,
    and the contours in the background
    represents the true noise-free oracle $f(\mathbf{x}, t)$}
    \label{f:quadratic_exp}
\end{figure}

    


\cref{f:quad_myopic_v_lookahead} compares the time history of the average oracle value (over 20 replications
with randomized starting points) resulting from decisions made via the myopic and lookahead acquisition functions. We present results for only
Quadratic-a and Quadratic-c for the sake of brevity. One can see that
the lookahead acquisition functions (dashed line plots) make decisions with
lower oracle value during the early rounds. While the acquisition function per se does not
\emph{seek} to make decisions with low current oracle value, the consequence of
looking ahead at $T$ is that decisions made at current $t$ may incur a low oracle value. This is a fundamental distinction between our proposed lookahead and the myopic acquisition functions. Also, the average oracle value at target
time $T$ (denoted by the empty circle) is consistently higher than the myopic
counterparts (denoted by the square). Furthermore, the predictions by the
lookahead acquisition functions are with higher confidence, as visualized by
the error bars (one standard deviation in the predictions) shown at $T$; note that we highlight the error bars for the
lookahead acquisition functions in black.

The average oracle values at $T$
predicted by all the acquisition functions are somewhat consistent, with the
\texttt{r2LEY} being one of the most competitive, including the rest of the test cases.
For this reason, moving forward we predominately
use the \texttt{r2LEY} as the representative lookahead acquisition function.
\begin{figure}
    \centering
    \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=.8\linewidth]{qa_myopic_v_lookahead.pdf}
        \caption{Quadratic-a}
        \label{sf:a}        
    \end{subfigure}\\
        \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=.8\linewidth]{qc_myopic_v_lookahead.pdf}
        \caption{Quadratic-c}
        \label{sf:c}        
    \end{subfigure}
\caption{Average of 20 replications with random
starting points of the myopic (lookahead) acquisition function on 
one-dimensional test cases. The final point at $T=4$ is highlighted with square (circle) symbols and  variability at $T$ is visualized via
colored (black) error bars.
The vertical dashed line represents the time until which $n=15$ starting points
are collected and the total budget is set to $q=45$.}
\label{f:quad_myopic_v_lookahead}
\end{figure}


\subsection{Synthetic higher dimensional test cases}
\label{ss:synthetic_hd}
We now demonstrate our method on higher-dimensional test functions, up to $d=6$.
We begin by presenting a modified Griewank function
($d=2$)~\cite{simulationlib}, where we multiply the original
Griewank function $f_{\textnormal{ref}}$ by a Gaussian weight function in order to create a
unique global maximum. Then, we modify the input $\mathbf{x}
\in \bm{\mathcal{X}}=[-5,5]^2$ with a time-dependent rotation to yield the function
\begin{equation}
    f(\mathbf{x},t) = f_{\textnormal{ref}}(\mathbf{x}(t)) \times \textnormal{exp}~\left( - \| \mathbf{x} - \mathbf{x}_0 \|^2/\ell\right), 
\end{equation}
where $\mathbf{x}_0 = [3, 0]$, $\ell \in \mathbb{R}_+$  is a scaling
constant\footnote{We set
$\ell = 160$ to ensure the weighting is somewhat mild and does not
significantly change the nature of the true Griewank function.} and $\mathbf{x}(t) = R(\zeta_t) \times \mathbf{x},$
where
\[ R(\zeta_t) = 
\begin{bmatrix}
\cos(\zeta_t) & -\sin(\zeta_t) \\
\sin(\zeta_t) & \cos(\zeta_t)
\end{bmatrix} \]
is a $2\times 2$ rotation matrix and $\zeta_t = \pi t /4$. Time snapshots of
this test function are shown in \cref{f:mod_griewank} and illustrate a 
counterclockwise rotation of the unique global maximizer with time. As before, all oracle observations include additive mean zero noise with variance $\sigma_\epsilon^2 = 10^{-3}$.

The performance of all the myopic acquisition functions is compared against \texttt{r2LEY} in \cref{f:mod_griewank_result}, where the
contour plot at $T$ is shown with the predictions from 20 repetitions of each
acquisition function (black $\times$ symbols) overlaid. Note that the total
budget is fixed at $q=90$, where $n=60$ is used in $\bm{\mathcal{X}} \times [2,3]$ to
start the algorithm, and the target time is $T=4$.  The
lookahead acquisition function clearly predicts the global maximum at $T$ with greater
consistency than the rest do; all but two repetitions predicted $\mathbf{x}_T$ such that $\| \mathbf{x}_T - \mathbf{x}^*_T\| < 0.26$.

\begin{figure}[htb!]
    \centering
    \includegraphics[trim={50 50 50 50}, clip, width=.8\linewidth]{modified_griewank.pdf}
    \caption{Time snapshots of the modified Griewank function. Notice the unique global maximum at each time and the counterclockwise rotation of the global maximizer with time.}
    \label{f:mod_griewank}
\end{figure}
\begin{figure}[htb!]
    \centering
    \includegraphics[width=.8\linewidth]{mod_griewank_result.pdf}
    \caption{Predictions at $T=4.0$ on the modified Griewank test case. Symbols represent predictions at $T$ after 20 repetitions of each algorithm, and the contours represent the noise-free $f(\mathbf{x}, T)$. The total budget is set to $q=90$, where $n=60$ samples in $2.0\leq t \leq 3.0$ are used to start each algorithm.}
    \label{f:mod_griewank_result}
\end{figure}


We also demonstrate our method on the Hartmann-3d ($d=3$) and Hartmann-6d
($d=6$) test functions. In both cases, the Hartmann
function~\cite{simulationlib} is summed with $f_{\mathbf{x} t} (\mathbf{x}, t) = \mathbf{1}^\top [2\textnormal{sin}(t)\mathbf{x}
- \textnormal{sin}^2(t)\mathbf{1}]$ (where $\mathbf{1}$ is a vector of ones of length $d$) to induce movement of the maximum with time.
For Hartmann-3d, the total budget is set to $q=120$ with $n=90$
samples queried in $[0, 1]^3 \times [2,3]$ to start the algorithm. For
Hartmann-6d, we set $q=150$ with $n=120$ samples queried in $[0, 1]^3 \times
[2,3]$ to start the algorithm, and our target time is $T=4$. 
\Cref{f:hartmann} shows the time histories of $f(\mathbf{x}_i,t)$
averaged over 20 replications. The predictions at $T$ are highlighted by adding
a marker in the plot. In both  test functions, the
lookahead acquisition function results in a higher oracle value $f(\mathbf{x}_T, T)$ at the
target time $T$ compared with other myopic acquisition functions. 

\begin{figure}[htb!]
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{hart3d.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{hart6d.pdf}
    \end{subfigure}
    \caption{Hartmann-3d ($d=3$) and Hartmann-6d ($d=6$) test functions. Average of 20 repetitions. Left of vertical line corresponds to starting samples. The error bars for the oracle value at $T$ for the best two acquisition functions (\texttt{r2LEY} and \texttt{mumax}) shown are slightly shifted in time for better visibility. The \texttt{r2LEY} shows comparatively less variability.}
   
    \label{f:hartmann}
\end{figure}

\subsection{Quantum optimal control}
\label{ss:qoc}
We next compare methods for optimizing time-dependent oracles on a challenging real-world problem.
Optimally shaping electromagnetic fields to control quantum systems is a widespread
application~\cite{brif2010control}; similar problems appear when controlling molecular transformations relevant to chemical,
biological, and materials science applications. Quantum optimal control (QOC) is
a method to shape electromagnetic fields to steer a quantum system toward a
desired control target; see
\cref{f:quantum_chemistry} where the target is a specific outcome from a chemical reaction. 

\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.5\linewidth]{Quantum_chemistry.PNG}
    \caption{Illustration of QOC to control chemical reactions. A femtosecond
    laser controlled by QOC is used to provide control over the selective
    dissociation of molecules. Figure source:~\cite{magann2020digital}}
    \label{f:quantum_chemistry}
\end{figure}

The problem we consider is based on a diatomic molecule \ce{HF}, whose vibrational energy mode is modeled as a nonrotating Morse oscillator on a 1D grid that varies with time and is illuminated with a laser
field to induce dissociation. The electromagnetic field is parameterized by a set of coefficients $\mathbf{x} \in \bm{\mathcal{X}} \subset \mathbb{R}^5$ that are the control parameters; more specifically, $\mathbf{x}$ contains the frequency components of the electromagnetic field. The digital quantum simulation of the Morse 
oscillator $f(\mathbf{x}, t)$ takes as input the control parameter $\mathbf{x}$ and a \emph{terminal} time $t
\in [0, T]$---which is not to be confused with our target time $T$---whose
 output is the observable control objective that needs to be minimized with respect to $\mathbf{x}$. 
The cost of each query of the simulation increases monotonically with the terminal time $t$, and
we are interested in finding the optimal control at a specific target time $T$. Since searching for the optimal control at $T$ is expensive,  we seek to predict 
 the optimal control at $T$ with cheaper evaluations of $f(\mathbf{x}, t)$ for $t<T$.

In this regard, we use BO to find the control parameters $\mathbf{x}_T^*$ that maximize the negative value of $f(\mathbf{x},T),
\forall \mathbf{x} \in \bm{\mathcal{X}}$. We begin by collecting observations of their system
$\bm{\mathcal{D}}_n= \{(\mathbf{x}_i, t_i), \hat{y}_i \},~i=1,\ldots,n$, which we use to fit
a GP that emulates $f$ in $\bm{\mathcal{X}} \times {\mathcal{T}}$. 
Then, we sequentially
make observations $\{(\mathbf{x}_i, t_i), \hat{y}_i \},~i=n+1,\ldots,q$, where 
each $\mathbf{x}_i$ is selected by maximizing our proposed two-step lookahead acquisition function. 
\Cref{f:quantum_problem} provides a schematic of
this problem. To conform to our problem setting, we allow evaluations only at increasing values of $t$. 
See~\cite{magann2020digital} for more information about this problem. 



We set the total budget to $q=180$, where $n=150$ samples are queried over $4000
\leq t \leq 5000$ to start the algorithm and $T=6000$. The time histories of
each acquisition function for the QOC problem are shown in
\cref{f:quantum_results}, where we include results from \texttt{r2LEY} and \texttt{r2LPI} among the lookahead approaches. 
The best-known  maximum $f(\mathbf{x}^*,
T)=0.0$ at $T=6000$, which is shown in the figure as the black square. In terms of the average performance, the
lookahead \texttt{r2LEY} outperforms all other methods, with the \texttt{mumax}
being a  close competitor, as shown in the top-left of \Cref{f:quantum_results}. We visualize the variability in the predictions at $T$ via boxplots in the top right of \Cref{f:quantum_results}, where we observe that the \texttt{r2LEY} has the highest median oracle value to interquartile distance ratio. This shows that although \texttt{r2LEY} does not result in the best median oracle value, it provides the most confident prediction, while still being competitive with \texttt{UCB} and \texttt{mumax}. 
Even though the lookahead method \texttt{r2LPI} does not necessarily outperform all the other myopic methods, it shows better performance---in terms of both average and median oracle value---compared with its myopic counterpart \texttt{PImumax}. This behavior is consistent with  what was observed with the one-dimensional test cases presented in \cref{ss:synthetic_1d}.



\begin{figure}[tb!]
    \centering
    \includegraphics[width=1\linewidth]{Quantum_problem3.pdf}
    \caption{Overview of the quantum optimal control problem recreated from \cite{magann2020digital}.
    The quantum system dynamics, driven by a field
    parameterized by $\mathbf{x}$, are simulated for the time range $[0, t]$ (where $t$
    is the terminal time), after which the control objective function
   
    $f(\mathbf{x},t)$ is calculated. A recent framework \cite{magann2020digital} is proposed
    for efficiently evaluating the control
    objective function 
   
    by using a future quantum computer.
    During an initial sample generation phase, the control objective function
   
    is evaluated $n$
    times to generate initial samples $\bm{\mathcal{D}}_n$ that are used by the GP to
    learn $f(\mathbf{x}, t),~ \mathbf{x}\in \bm{\mathcal{X}}, t\in [0, t]$. Then, in the optimization
    phase, the trained GP 
   
    is used in the BO framework to find the optimal
    control at target time $T$. The inputs at $(\mathbf{x}_i,t_i),~i=,n+1,\ldots,q$ are
    chosen at \emph{optimal} locations as guided by the BO acquisition function
    optimization.}
    \label{f:quantum_problem}
\end{figure}
\begin{figure}
    \centering
    \begin{subfigure}[t]{.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{quantum.pdf}
        \caption{Average $f(\mathbf{x},t)$ over 20 repetitions with randomized starting points. Left of vertical line denotes starting points.}
    \end{subfigure}%
        \begin{subfigure}[t]{.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{quantum_boxplot.pdf}
        \caption{Variability in the oracle value at $T$, $f(\mathbf{x}, T)$. }
    \end{subfigure}\\
   
   
   
   
   
   
   

    \caption{Demonstration on the quantum optimal control problem.}
    \label{f:quantum_results}
\end{figure}






















\subsection{Implementation details}
We implemented our recursive lookahead Bayesian optimization framework
in PyTorch~\cite{pytorch_bib} using the GP and BO libraries GPyTorch~\cite{gardner2018gpytorch} and
BoTorch, respectively~\cite{balandat2019botorch}. For
finite-horizon time-dependent BO with limited budget, the lookahead
acquisition functions outperform their myopic counterparts in securing the best
oracle value at the horizon. This improved accuracy comes at the cost of additional
computational overhead in terms of optimizing the acquisition functions. Our
Monte Carlo approximation to the acquisition function and its gradients, as presented in Line 5 of 
Algorithm~\ref{a:r2LBO}, lends
itself well to be used with a gradient-based optimizer.
For further computational improvements, however, we
also adopt the \emph{one-shot} optimization of the lookahead acquisition
function as implemented in BoTorch. The approach formulates line 4 of Algorithm~\ref{a:r2LBO} as a deterministic optimization problem over the joint $\mathbf{x}$ and $\mathbf{x}'$, where $\mathbf{x}' := \{ \mathbf{x}^j\}_{j=1}^N$ are the candidate maximizers or \emph{fantasy} points. While this increases the dimensionality of the optimization problem from $d$ to $N+d$, we observe consistent results (with our Monte Carlo approach) with $N$ set to as small as 32, resulting in significant computational speedup. Moreover, the one-shot approach solves a deterministic optimization problem as opposed to a stochastic optimization problem with the Monte Carlo approach.
We include
both the Monte Carlo and one-shot approaches in our software and recommend the
latter as the dimensionality $d$ of the problem increases; all the results presented in this work are based on the one-shot approach. 
For the $d=5$
QOC problem, the acquisition function optimization takes approximately
30 seconds of wall-clock time when run in serial. Because many of the time-dependent oracles we are targeting require significantly more computational cost and time, we consider the per-iterate cost of our
approach to be trivial. Furthermore, the PyTorch framework allows
for leveraging advanced computation hardware such as graphics processing units, which may allow for further computational efficiency.

\section{Conclusions and future work}
Bayesian optimization with a lookahead acquisition
function outperforms well-known myopic acquisition functions in solving the finite-budget, finite-horizon, time-dependent maximization of an
expensive stochastic oracle.
Specifically, for such problems, the lookahead approach
makes probabilistically sound decisions that give better average-case
returns at the target horizon, while considering 
strategies that
lead up to the target horizon. In this regard, the lookahead acquisition
functions take advantage of the epistemic uncertainty---explained by the GP---in
the predictions about the oracle. In this work we introduce a generalized
framework that offers lookahead extensions to widely known myopic acquisition
functions. We demonstrate our method on 
illustrative 
synthetic test functions (that involve discontinuities in the
oracle maximizer) as well as a quantum optimal control application problem.
Because of the lack of analytically closed-form expressions for our lookahead
acquisition functions, their construction (and optimization) incurs more
computational overhead than do the myopic acquisition functions. Such an
overhead is insignificant, however, when optimizing expensive stochastic oracles that take, for example, hours of wall-clock time per evaluation. 
Overall, we demonstrate that our
proposed acquisition function results in more confident and accurate
predictions of the global maxima at the target horizon compared with Bayesian optimization using myopic acquisition functions.

The primary challenge with handling time-dependent optimization problems with
BO is that the effect of the dynamics of the oracle has to be learned by a GP model. When the
distance to the horizon is large, it is particularly challenging to train a
\emph{global} GP that emulates the oracle in $\bm{\mathcal{X}} \times \mathcal{T}$. A natural
alternative is to work with \emph{local} GPs within
local
regions that
are dynamically updated. This is one avenue for future work that we
are considering.

A related problem is multifidelity
BO (see, e.g., \cite{song2019general, kandasamy2017multi, takeno2019multi}) in which one seeks to solve
$    \max_{\mathbf{x} \in \bm{\mathcal{X}} }~f(\mathbf{x})$ using observations of $\hat{f}(\mathbf{x}, s)$. The \emph{fidelity} is given by the
parameter $s\in[0,1]$, with $f(\mathbf{x}):= \hat{f}(\mathbf{x}, 1)$ corresponding to the
highest fidelity. Typically, there is a cost function $c(s)\in \mathbb{R}$, which is monotonically increasing in the fidelity
level $s$. 
By placing a GP prior on the function $\hat{f}(\mathbf{x}, s)$ directly, 
one can
sequentially optimize acquisition functions such as
\[ \max_{(\mathbf{x},s) \in \bm{\mathcal{X}}\times [0,1]}~\frac{\alpha(\mathbf{x}, s)}{c(s)}.\]

Our problem is equivalent to the multifidelity BO problem with $s$
replaced by $t\in[0,T]$ with the following key differences:
\begin{enumerate}
    \item In our case, $\hat{f}(\mathbf{x},T)$ can never be observed when $t< T$. More generally, any
      set $\{t_i \}, i=1, \ldots, n$ has a unique ordering $t_1 < \ldots < t_n$;
      and while we are at any $t_i$, evaluation of $f(\mathbf{x}, t \neq t_i)$ is not possible.
    \item We make only one observation $f(\mathbf{x}, t)$ for each $t\in[0,T]$.
    \item The cost of each evaluation $c(s) = c$ is a constant. However, similar to the general multifidelity framework, we
    are still interested in maximizing $f(\mathbf{x}, T)$.
\end{enumerate}
Another area of future work is defining a time-dependent cost
function that can be used within the multifidelity framework to also select the
time schedule of observations automatically; this approach  is of specific interest in
handling continuous-time systems.


\section{Introduction}
\label{s:Introduction}
We consider the maximization of an expensive-to-evaluate oracle $f(\mathbf{x}, t)$,
where $\mathbf{x}\in \bm{\mathcal{X}} \subset \mathbb{R}^d$ are the inputs or \emph{action} parameters
in a compact domain and $t$ is the \emph{context} from a (possibly infinite)
set of contexts $\mathcal{T}$. An observation $\hat{y} \in \mathbb{R}$ made by observing
$f$ at the action-context pair $(\mathbf{x}, t)$ is assumed to be corrupted by
additive stochastic noise $\epsilon$; that is, $\hat{y} = f(\mathbf{x},t) + \epsilon$.
Furthermore, we assume a context space $\mathcal{T}$ whose members have a unique
ordering (hereafter ``time''). We seek an action $\mathbf{x}^*_T$ at a given future time $T \in \mathcal{T}$ that solves
\begin{equation}
\operatornamewithlimits{arg\,max}_{\mathbf{x} \in \bm{\mathcal{X}}} f(\mathbf{x}, T)
\label{eq:prob}    
\end{equation}
in as few stochastic observations of $f$ as possible.
Such a problem fundamentally differs from general
context-dependent bandit problems \cite{krause2011contextual,srinivas2009gaussian} in that we are interested solely in
the action that maximizes the oracle value $f$ 
at time $t=T$, as opposed to maximizing the
\emph{cumulative oracle value} until $T$. Additionally, we
assume that only one observation can be made per time 
$t$.
The time-dependent maximization that we address arises, for example, in quantum
computing and control~\cite{magann2020digital,zhu2018training}, where 
values for the control parameter $\mathbf{x}$ are sought by using
noisy observations $f(\mathbf{x},t)$ 
to achieve a
desired task at future time $T$; see \cref{ss:qoc}. Another application  relevant to our work
is multistage stochastic programming for portfolio optimization \cite[Ch. 1]{shapiro2014lectures}, where the goal is to find a 
policy $\{(\mathbf{x}_i, t_i) :~i=0,\ldots,n-1\}$ for redistributing assets 
to maximize the expected profit at $T=t_n$ and the profit at each stage is
realized by querying a noisy oracle. 
In many settings, $\bm{\mathcal{X}}$ is a simple domain; here 
we assume that $\bm{\mathcal{X}}:= [0,1]^d$ so that the inputs $\mathbf{x}$ are normalized.


The main challenges of the problem \eqref{eq:prob} that we seek to address are (i) the absence of structure (e.g., convexity) in $f$ that
could potentially be exploited, 
(ii) the high cost
associated with each
observation of $f$, and 
(iii) the time-varying setting where observations of $f$ cannot be made in the past or future.
Bayesian optimization (BO)~\cite{brochu2010tutorial, shahriari2015taking} with
Gaussian process (GP) priors~\cite{gramacy2020surrogates,rasmussen:williams:2006}  suits the
specific
challenges posed by our problem (mainly (i) and (ii) mentioned
above),
where observations at each round $i$, 
$\hat{y}_i =
f(\mathbf{x}_i, t_i) + \epsilon_i$, 
are made judiciously by leveraging information
gained from previous observations, as a means of coping with the high costs of
each observation. The key idea  is to specify GP prior distributions on the
oracle value and the noise. That is, $f(\mathbf{x},t) \sim \mathcal{GP}\left(0, k( (\mathbf{x}, t),
(\mathbf{x}',t'))\right)$ and $\epsilon \sim \mathcal{GP}(0,\sigma^2_\epsilon)$, where $\sigma^2_\epsilon$ is a constant noise variance. 

The covariance function $k$ captures the correlation between the observations
in the joint $(\mathbf{x},t)$ space; here we use the product composite form
given by $k((\mathbf{x}, t), (\mathbf{x}', t')) = k_{\mathbf{x}}(\mathbf{x},\mathbf{x}'; \theta_\mathbf{x}) \times k_t(t,t';
\theta_t) $, where $\theta_\mathbf{x}$ and $\theta_t$ 
parameterize the covariance functions for $\mathbf{x}$ and $t$, respectively. We estimate the
GP hyperparameters 
$\bm{\Omega}= \lbrace \theta_x, \theta_t, \sigma^2_\epsilon \rbrace$
from data by maximizing the marginal
likelihood. We use the squared-exponential  kernels defined as $k_\mathbf{x}(\mathbf{x}, \mathbf{x}') = \textnormal{exp} \left( -\frac{\|\mathbf{x} -
\mathbf{x}'\|^2}{2 \theta^2_\mathbf{x}}\right)$ and $k_t(t, t') = \textnormal{exp} \left( -\frac{\|t -
t'\|^2}{2 \theta^2_t}\right)$.
The
posterior predictive distribution of the output $Y$, conditioned on available observations from the oracle, is given by 
\begin{equation}
    \begin{split}
        Y(\mathbf{x}, t |
      \bm{\mathcal{D}}_n, \bm{\Omega}) &\sim \mathcal{GP}(\mu_n(\mathbf{x}, t), \sigma^2_n (\mathbf{x}, t)),~ \quad \bm{\mathcal{D}}_n :=  \lbrace
(\mathbf{x}_i, t_i), \hat{y}_i \rbrace _{i=1}^n \\
\mu_n(\mathbf{x},t) &= \mathbf{k}_n^\top [\mathbf{K}_n + \sigma_\epsilon^2 \mathbf{I}]^{-1} \mathbf{y}_n\\
\sigma^2_n(\mathbf{x}, t) &=  k((\mathbf{x},t), (\mathbf{x},t)) - \mathbf{k}_n^\top [\mathbf{K}_n + \sigma_\epsilon^2 \mathbf{I}]^{-1} \mathbf{k}_n,
    \end{split}
    \label{e:GP}
\end{equation}
where $\mathbf{k}_n$ is a vector of covariances between $(\mathbf{x},t)$ and all observed
points in $\bm{\mathcal{D}}_n$, $\mathbf{K}_n$ is a sample covariance matrix of observed points in
$\bm{\mathcal{D}}_n$, $\mathbf{I}$ is the identity matrix, and $\mathbf{y}_n$ is the vector of all
observations in $\bm{\mathcal{D}}_n$; \eqref{e:GP} is then used as a surrogate model for
$f$ in Bayesian optimization. Note that $\mu_n$ and $\sigma_n^2$
are the posterior mean and variance of the GP, respectively, where the
subscript $n$ implies the conditioning based on $n$ past observations. BO then
proceeds by defining an \emph{acquisition function}, $\alpha(\mathbf{x},t)$, in terms of the GP posterior that is
optimized in lieu of the expensive $f$ to select the next
point $\mathbf{x}_{n+1}$ at $t_{n+1}$, and the process continues until a budget of $q$ observations is reached; see Algorithm~\ref{a:BO}. We denote by $\mathbf{x}_i$ and $\mathbf{x}_t$ a decision made by maximizing the acquisition function at $t_i$ and $t$, respectively. We use $\mathbf{x}_t^*$ to denote a maximizer of $f(\cdot, t)$.


\begin{algorithm2e}[t]
\textbf{Given:} $\bm{\mathcal{D}}_n = \lbrace
(\mathbf{x}_i, t_i), \hat{y}_i \rbrace _{i=1}^n$, 
 total budget $q$, 
 schedule $\lbrace t_{n+1},\ldots, t_{q}= T \rbrace$,
 and GP hyperparameters $\bm{\Omega}$ \\
\KwResult{$(\mathbf{x}_T, \hat{y}_T)$}
  \For{$i=n+1, \ldots, q$, }{
  Find $\mathbf{x}_i \in \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\operatornamewithlimits{arg\,max}}~ \alpha(\mathbf{x}, t_i)$ \qquad (acquisition function maximization based on $\bm{\Omega}$)\\
    Observe $\hat{y}_i$ = $f(\mathbf{x}_i, t_i) + \epsilon_i$\\ 
    Append $\bm{\mathcal{D}}_i = \bm{\mathcal{D}}_{i-1} \cup \lbrace (\mathbf{x}_i, t_i), \hat{y}_i \rbrace$\\ 
    Update GP hyperparameters $\bm{\Omega}$ \\
 }
 \caption{Generic Bayesian optimization}
 \label{a:BO}
\end{algorithm2e}

Typical acquisition functions in BO take a \emph{greedy} (``myopic'')
approach, where each $\mathbf{x}_i$ 
is selected as though it were the
last, without accounting for the potential impact on future selections.
Such myopic acquisition functions include the
probability of improvement (\texttt{PI})~\cite{kushner1964new}, the expected improvement
(\texttt{EI}) \cite{jones1998efficient, mockus1978application}, and the GP upper
confidence bound (\texttt{UCB})~\cite{srinivas2009gaussian}, which are, respectively,
\begin{equation}
    \begin{split}
        \alpha_{\textrm{PI}}(\mathbf{x}, t) &= \Phi \left( z \right) \\
        \alpha_{\textrm{EI}}(\mathbf{x}, t) &= z\sigma_n(\mathbf{x}, t) \Phi\left(z \right) + \sigma_n(\mathbf{x}, t) \phi \left( z \right) \\
        \alpha_{\textrm{UCB}}(\mathbf{x},t) &= \mu_n(\mathbf{x}, t) + \beta^{1/2} \sigma_n(\mathbf{x}, t),
    \end{split}
    \label{e:myopic_acqs}
\end{equation}
where $z = \frac{\mu_n(\mathbf{x}, t) - \xi}{\sigma_n(\mathbf{x}, t)}$; $\Phi$ and $\phi$ denote
the standard normal cumulative and probability density functions, respectively;
$\xi$ is a user-specified target; and $\beta$ is a confidence parameter.
Whereas the \texttt{PI} and \texttt{EI} acquisition functions are probabilistic measures of
improvement over a user-specified target, \texttt{UCB} is an optimistic
estimate. 
Under suitable regularity conditions, the \texttt{PI}, \texttt{EI}, and
\texttt{UCB} acquisition functions are guaranteed to achieve asymptotic
consistency~\cite{bull2011convergence, mockus1978application,
srinivas2009gaussian, vazquez2010convergence} in BO. For a limited budget $q$, however, 
such acquisition functions can be suboptimal~\cite{gonzalez2016glasses}.
Furthermore, the myopic nature of these basic acquisition functions is such that they seek to maximize the 
oracle value at $t=T$ only when they arrive at $T$, meaning that the past ($t<T$) decisions 
are not necessarily made to facilitate making the best decision for 
$t=T$. On the other hand, we are interested in a strategy that seeks to
maximize $f(\mathbf{x},T)$. 
Since the time horizon $t=T<\infty$, this is a \emph{finite-horizon} problem.

For time-dependent objectives, an optimal decision may be horizon-dependent,
which poses a challenge for improvement-based acquisition functions such as \texttt{PI} and \texttt{EI}. Because the target $\xi$ is typically chosen as the best observed value in $\{\hat{y}_i\}_{i=1}^n$, specifying an appropriate target when the maximum oracle value $\max f(\mathbf{x},t)$
changes with time $t$ is challenging. For example, what if the maximum oracle value at $T$ is
much smaller than anything observed until current time $t$? A further challenge in our context is that observations at the target time $T$ are not realizable until we conclude the optimization. Therefore, in such situations, myopic 
acquisition functions that seek
the best oracle value at the current $t$ may not necessarily be  useful. Instead,
we are interested in acquisition functions that  optimize the
\emph{long-sighted} decision (at $t=T$) even if it means incurring a suboptimal
oracle value at intermediate $t<T$. 
In the BO literature, strategies that make decisions by accounting for the
decisions remaining in the future are generally referred to as \emph{lookahead}
approaches, which we adopt to solve our  finite-horizon problem. 


We now illustrate the benefit of a lookahead acquisition function in the
context of time-dependent BO. \Cref{fig:EI_v_2LEY} shows the difference
between the decisions made via EI and the lookahead acquisition function
\texttt{2LEY} introduced in \cref{s:mLBO}. The
top left of the figure shows contours of the true (noise-free) oracle with a target horizon of 
$T=3.0$. We assume we have made $n$ decisions before $t_{n+1}=2.625$, and we
are interested in making a decision at current $t_{n+1}$ (denoted by the horizontal
dashed line). The top right shows the $\texttt{2LEY}$ and $\texttt{EI}$
acquisition functions and the corresponding decisions made by maximizing each
of them (denoted by vertical dashed lines). In the bottom row, we show the GP
predictions at $T$, with each of the decisions made via the myopic and
lookahead acquisition functions incorporated, along with the noise-free $f(\mathbf{x},
T)$ for reference, and the final decisions at $T$ (denoted by vertical dashed
lines). Notice that the lookahead decision results in a more accurate
prediction of the global maximum at $T$ compared with the myopic decision. This illustration shows that 
lookahead acquisition functions have the potential to make better decisions
than do myopic acquisition functions when the oracle value at a target (future) horizon is of interest. One could ask whether a randomly made decision that offers a \emph{space filling} effect (e.g., $0\leq \mathbf{x} \leq 0.4$ at $t=2.625$ ) could better learn $f(\mathbf{x}, T)$ and hence lead to a better final decision at $T$. We show that this is not the case, with a randomly chosen decision $(0.2474, 2.625)$ that does not outperform the lookahead method; see bottom right of \Cref{fig:EI_v_2LEY}.

\begin{figure}[htbp!]
    \centering
    \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{EI_v_2LEY.pdf} 
        \caption{Left: contour of $f(\mathbf{x}, t)$ in $\bm{\mathcal{X}}\times \mathcal{T}$ with $T=3$. Red circles denote the locations of $n$ previous observations. Horizontal line denotes the time $t_{n+1}=2.625$ at which the next observation is to be made. Right: myopic (\texttt{EI}) and lookahead (\texttt{2LEY}) acquisition functions. }
    \end{subfigure}\\
    \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{EI_v_2LEY_v_Random_T.pdf}   
        \caption{GP prediction at $T$ with decisions incorporated from myopic (\texttt{EI}) and lookahead (\texttt{2LEY}) decisions as well as a decision chosen uniformly at random from $\boldsymbol{\mathcal{X}}$ at $t$. Shaded region represents 95\% confidence region.}
    \end{subfigure}
    \caption{Comparison of the decisions made via a myopic (\texttt{EI}) 
    and a lookahead (\texttt{2LEY}) acquisition function. The top right shows the acquisition functions and the respective decisions (maximizers). The GP predictions, incorporating the decision made via either of the acquisition functions, are shown in the bottom. One can see that if the final decision at $T$ is made as   $\operatornamewithlimits{arg\,max}_{\mathbf{x}\in\bm{\mathcal{X}}}~\mu_{n+1}(\mathbf{x}, T)$, then the lookahead approach predicts the maximizer at $T$ more accurately than does the myopic approach. A random decision $(0.2474, 2.625)$, even though it is \emph{space filling}, does not outperform the lookahead decision.}
    \label{fig:EI_v_2LEY}
\end{figure}

In foundational work on lookahead BO, Osborne \cite{osborne2010bayesian}
showed that proper Bayesian reasoning can be used to define an acquisition
function where the ($n+1$)th observation is made by marginalizing a loss function
at the final ($q$th) observation  over all the remaining observations $(\mathbf{x}_i, \hat{y}_i)$,
$i=n+2,\ldots, q$. Others have proposed
\emph{lookahead EI} in the context of finite-budget BO. For example, Ginsbourger and Le
Riche \cite{ginsbourger2010towards} showed that the lookahead EI is a dynamic
program that chooses the EI maximizer in expectation, considering all possible
strategies of the same budget. 
A practical challenge with lookahead strategies is the computational
tractability (as demonstrated in, for example, \cite{ginsbourger2010towards,osborne2010bayesian}), which can be partially mitigated by various
approximation methods \cite{gonzalez2016glasses,lam2016bayesian,wu2019practical}. Yue and Kontar \cite{yue2020non} explored the
theoretical properties of the approximate dynamic programming approach for
multistep lookahead acquisition functions and discussed when these strategies are guaranteed
to perform better than their myopic counterparts. In this work
we draw inspiration from \cite{osborne2010bayesian} to address the
limited budget aspect of our problem. However, we propose a general framework
for lookahead BO that readily applies to the finite-horizon,
time-dependent BO; we are unaware of any previous work on this topic.

In the context of
time-dependent BO, the only works we are aware of are
\cite{nyikosabayesian} and \cite{bogunovic2016time}. 
We fundamentally differ
from their approach in the following ways: (i) their goal is to track a
time-varying optimum whereas ours is to predict it at the target horizon $T$), and (ii) they do not address
the limited budget constraint as we do. Furthermore, setting $k_t(t_i,
t_j)=(1-\epsilon)^{|t_i - t_j|/2}$ (where $\epsilon$ is the \emph{forgetting}
factor as defined in \cite{bogunovic2016time}), our method is equivalent to
\cite{bogunovic2016time} if we myopically make decisions at each $t$ via \texttt{UCB}.  



Another plausible way to solve the finite-horizon (target $t=T$) BO problem
is to
sequentially choose points that are most
\emph{informative}~\cite{cover2012elements} about the maximizer
at the target horizon; see, for example, 
\cite{hennig2012entropy, hernandez2014predictive, wang2017max}. Although such
information-theoretic approaches are known to involve an intractable form of
the information entropies, they are also highly reliant on an accurate underlying GP model~\cite{mackay1992information}. 


In this work, we extend myopic acquisition functions in order to make
nonmyopic decisions for
finite-horizon BO. Specifically, we propose a generalized
framework that specifies a \emph{value function} 
that quantifies an objective at the target horizon $T$ and that can then be used in a lookahead acquisition
function. Our main contributions are as follows. 
\begin{enumerate}
    \item We present the first strategy to solve finite-horizon BO for time-dependent 
   
    oracles.
    \item We propose lookahead versions of myopic acquisition functions~(\cref{s:mLBO}) that
      are applicable to both finite-budget BO and finite-horizon, time-dependent settings. 
    \item We present a practical algorithm that recursively applies a two-step
      lookahead acquisition function~(\cref{s:r2L}) 
     
   
   
   
      to make new observations time-dependent oracles.
    \item We establish the utility of our method by demonstrating it on carefully constructed synthetic test functions (\cref{ss:synthetic_1d,ss:synthetic_hd}) and a real-world quantum optimal control problem~(\cref{ss:qoc}). 
\end{enumerate}
Software that implements our lookahead framework for limited horizon, finite-budget time-dependent BO problems is available upon request and will be released publicly with the final version of this manuscript.

The remainder of the paper is organized as follows. In  \cref{s:mLBO}
we present our generalized lookahead acquisition framework for generic value
functions. In \cref{s:r2L} we present a recursive two-step lookahead
approach that is more practical than a more-than-two-step lookahead approach.
We discuss theoretical properties of our approach in \cref{s:theory}.
We then
present the results of our numerical experiments on synthetic test functions
(\cref{ss:synthetic_1d,ss:synthetic_hd}) as well as a
real-world quantum optimal control problem in \cref{ss:qoc}. We conclude the paper with an outlook on future work.

\section{Generalized m-step lookahead acquisition function}
\label{s:mLBO}
Given the data from $n$ observations $\bm{\mathcal{D}}_n$, let us define a value function $\upsilon$ by 
\begin{equation}
    \upsilon_n(\mathbf{x}, t; \boldsymbol{\omega}) := \mathbb{E}_{y\sim Y|\bm{\mathcal{D}}_n}\left[
   
    h(y(\mathbf{x}, t); \boldsymbol{\omega})     
    \right],
    \label{e:value_function}
\end{equation}
where $h$ is a scalar-valued function, $\boldsymbol{\omega} \in \mathbb{R}^p$ are
parameters that parameterize $h$, 
$Y|\bm{\mathcal{D}}_n$ is the GP posterior 
given $n$ observations, and $y(\mathbf{x}, t)$ is a sample path realized from $Y(\mathbf{x},
t)$. For example, the \texttt{PI}, \texttt{EI}, and \texttt{UCB} acquisition functions can be obtained from defining $h$ as 

\begin{equation}
    \begin{split}
   
    h(y(\mathbf{x}, t); \xi) 
    &= y(\mathbf{x}, t) - \xi \implies \upsilon(\mathbf{x},t;\xi) = \Phi \left( z \right) \qquad \text{(PI)}\\
   
    h(y(\mathbf{x}, t); \xi) 
    &= [y(\mathbf{x}, t) - \xi]^+ \implies \upsilon(\mathbf{x},t;\xi) = z\sigma(\mathbf{x}, t) \Phi\left(z \right) + \sigma_n(\mathbf{x}, t) \phi \left( z \right) \qquad \text{(EI)}\\
       
        h(y(\mathbf{x}, t); \beta)
        &= y(\mathbf{x}, t) +  \beta^{1/2}  \sigma(\mathbf{x},t) \implies \upsilon(\mathbf{x},t;\beta) = \mu(\mathbf{x}, t) + \beta^{1/2}\sigma(\mathbf{x}, t), \qquad \text{(UCB)} 
        \end{split}
    \label{e:value_functions}
\end{equation}
where recall that $\xi$ and $\beta$ are user-specified target and confidence parameters, respectively.
Similarly, when $h$ is the identity function $h(y(\mathbf{x}, t); \boldsymbol{\omega})  = y(\mathbf{x}, t)$, then $\upsilon(\mathbf{x},t)= \mu(\mathbf{x}, t)$. 

We use the following notation for conciseness:
\begin{align*}
&\mathbf{x}^t_{n+1} := (\mathbf{x}, t_{n+1}) && \bm{\mathcal{D}}_{n,j} := \bm{\mathcal{D}}_n \bigcup_{i=1,\ldots,j} \{ \mathbf{x}^t_{n+i}, y_{n+i} \}\\ 
&y_{n+1} := y(\mathbf{x}^t_{n+1}) &&\mathbb{E}_{n,j} := \mathbb{E}_{y\sim Y|\bm{\mathcal{D}}_{n,j}} \\
&\mu_{n,j} := \mathbb{E}_{n,j}Y, ~~\sigma^2_{n,j} := \mathbb{E}_{n,j} \left( Y - \mu_{n,j}\right)^2 && \upsilon^T_{n,j} := \upsilon_{n,j}(\mathbf{x}, T; \boldsymbol{\omega}) = \mathbb{E}_{n,j}~\left[h(y(\mathbf{x}, T); \boldsymbol{\omega}) \right],
\end{align*}
where we note
that whenever we explicitly use $y(\mathbf{x},t)$, we refer to a sample path drawn from
$Y(\mathbf{x},t)$ and therefore $y_{n+1} = y(\mathbf{x}_{n+1},t_{n+1})$ is a draw from $y(\mathbf{x},t)$ at $(\mathbf{x}_{n+1}, t_{n+1})$. 
On the other hand, when we use $\hat{y}_i$, we refer to an observation of the oracle at
$(\mathbf{x}_i,t_i)$. Similarly, note that $\bm{\mathcal{D}}_{n,k}$ contains $n$ oracle evaluations 
and $k$ candidate points $\mathbf{x}$ and the $k$ corresponding draws from the GP at times 
$t_{n+1}$, \ldots, $t_{n+k}$. Furthermore, we note here that $\mathbf{x}^t_{n+1}$ is $d$-dimensional.

We define our $m$-step lookahead acquisition function as
\begin{equation}
    \begin{split}
      \alpha_m(\mathbf{x}^t_{n+1}) = \int \ldots \int  & \left[\max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{m-1}(\tilde{\mathbf{x}}, T; \boldsymbol{\omega})|\bm{\mathcal{D}}_{n,m-1}\right] ~ p(y_{n+m-2}|\bm{\mathcal{D}}_{n,m-2})\\ &p(\mathbf{x}^t_{n+m-1}|\bm{\mathcal{D}}_{n,m-2})\cdots 
    p(y_{n}|\bm{\mathcal{D}}_n)~dy_{n+m-2}\cdots dy_{n}~d\mathbf{x}^t_{n,m-1} \cdots d\mathbf{x}^t_{n+2}.
    \end{split}
    \label{e:mstep}
\end{equation}
Note that in \eqref{e:mstep}, the \emph{inner maximization} is always with respect to $\mathbf{x}$ at $t=T$ and that after evaluating the integral, the resulting expression is a function of $\mathbf{x}^t_{n+1}$. Furthermore, since we are interested in choosing the maximizer of $\alpha_m(\mathbf{x}^t_{n+1})$ for $t_{n+1}$, we do not marginalize $\mathbf{x}^t_{n+1}$. Computing \eqref{e:mstep} involves the computation of $m-1$ nested expectations. Additionally, when the expectations are approximated via Monte Carlo sampling, there is a total of $N\times (m-1)$ $d$-dimensional  optimization solves due to the inner maximization in \eqref{e:mstep}. The computation or even the Monte Carlo approximation of \eqref{e:mstep} becomes intractable as $m$ increases, warranting further approximation, as is done, for instance, in \cite{gonzalez2016glasses, lam2016bayesian, yue2020non}. When $m=2$, however, the computation of \eqref{e:mstep} is much more straightforward, as we discuss in the following section.









\section{Generalized two-step lookahead acquisition function}
\label{s:r2L}
In this work, as a practical alternative to the computationally intractable multistep ($m>2$) 
lookahead acquisition function, we recursively optimize the two-step lookahead acquisition function~\eqref{e:E_twostep} at each time step of the time-dependent BO problem; this process is graphically described in \cref{f:method}.
Given $\bm{\mathcal{D}}_n$ and with $m=2$, we define our \emph{two-step} lookahead acquisition function by
\begin{equation}
\begin{array}{rl}
 
 


  \alpha_2(\mathbf{x},t_{n+1})
  
&=  \displaystyle \int \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~ \left[ \upsilon_{n,1}(\tilde{\mathbf{x}}, T; \boldsymbol{\omega}) | \, 
\bm{\mathcal{D}}_n \bigcup \{ (\mathbf{x},t_{n+1}), y_{n+1} \}
\right]~ p(y_{n}|\bm{\mathcal{D}}_{n}) ~dy_{n}, 
\end{array}
\label{e:twostep}
\end{equation}
where $p(y_{n}|\bm{\mathcal{D}}_{n})$ is the probability density of the GP posterior distribution
conditioned on $\bm{\mathcal{D}}_n$.
We write \eqref{e:twostep} concisely as 
\begin{equation}
  \alpha_2(\mathbf{x}^t_{n+1}) = \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}(\tilde{\mathbf{x}},T;\boldsymbol{\omega})|\bm{\mathcal{D}}_{n,1} \right].
    \label{e:E_twostep}
\end{equation}
In \eqref{e:E_twostep}, recall that $\bm{\mathcal{D}}_{n,1}$ is the union of $\bm{\mathcal{D}}_n$ and the $n+1$th observation being a draw from the GP $Y_n$, and we compute $\upsilon_{n,1}(\tilde{\mathbf{x}}, T; \boldsymbol{\omega})$ by taking the expectation with respect to $Y|\bm{\mathcal{D}}_{n,1}$. In general, it is not guaranteed that
\eqref{e:E_twostep} admits a closed-form expression. Hence we make the
Monte Carlo approximation
\begin{equation}
    \alpha_2(\mathbf{x}^t_{n+1}) \approx \hat{\alpha}_2(\mathbf{x}^t_{n+1}) := \frac{1}{N} \sum_{j=1}^N \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}^{T,j},
    \label{e:twostep_mc}
\end{equation}
where $\upsilona{n}{1}^{T,j} = \upsilona{n}{1}(\tilde{\mathbf{x}}, T)|y^j_{n+1}$ with the $N$ iid samples $\{y^j_{n+1}: \,j=1,\ldots,N\}$ drawn from $Y_{n}$. In practice, the computation of \eqref{e:E_twostep} involves the following steps (given $\mathbf{x}_{n+1}$):
\begin{enumerate}
    \item Generate $N$ iid samples $\{ y^j_{n+1}\}_{i=1}^N$ drawing from $\mathcal{N}(\mu_n(\mathbf{x}_{n+1}, t_{n+1}), \sigma_n^2(\mathbf{x}_{n+1}, t_{n+1}))$.
    \item For each $y^j_{n+1},~j=1,\ldots,N$
    \begin{enumerate}
        \item[2.1] update $\bm{\mathcal{D}}_n$ to $\bm{\mathcal{D}}_{n,1}$ to get GP posterior $Y|\bm{\mathcal{D}}_{n,1}$ and
        \item[2.2] maximize a draw of the projected value function $\nu^j = \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}(\tilde{\mathbf{x}}, T)|\bm{\mathcal{D}}_{n,1}$.
    \end{enumerate}
    \item Compute $\hat{\alpha}_2(\mathbf{x}_{n+1}) \approx \frac{1}{N} \sum_{j=1}^N \nu^j$.
\end{enumerate}

In order to make a decision at $t_{n+1}$, \eqref{e:E_twostep} needs to be maximized. We show how the maximization of the acquisition function can be done efficiently by constructing unbiased estimators of the acquisition function's gradient. Under suitable regularity conditions, we can interchange the expectation and gradient
operators (see \cref{thm:gradexp}
) to obtain
\begin{equation}
    \nabla_{\mathbf{x}}\alpha_{2}(\mathbf{x}^t_{n+1}) = \mathbb{E}_{n} \left[ \nabla_{\mathbf{x}}~ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}(\tilde{\mathbf{x}}, T) | \bm{\mathcal{D}}_{n,1} \right].
    \label{e:twostep_T_grad}
\end{equation}
Similar to the Monte Carlo approximation for the acquisition function presented
in \eqref{e:twostep_mc}, \eqref{e:twostep_T_grad} can be approximated as
\begin{equation}
    \nabla_{\mathbf{x}}\alpha_{2}(\mathbf{x}^t_{n+1}) \approx \tilde{\nabla}_{\mathbf{x}}\alpha_{2}(\mathbf{x}^t_{n+1}) = \frac{1}{N} \sum_{i=1}^{N} d\nu^i :=  \frac{1}{N} \sum_{i=1}^{N}  \nabla_{\mathbf{x}}~ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\upsilona{n}{1}^{T,i}.
    \label{e:twostep_T_grad_mc}
\end{equation}
In this way, sample-averaged gradients can be used in a gradient-based optimizer
to efficiently optimize \eqref{e:E_twostep}. We clarify the issues associated
with computing the gradients of functions involving the $\max()$ operator as in
\eqref{e:twostep_T_grad_mc} in \Cref{s:theory}. 

We now illustrate in detail the two-step
acquisition function for a specific value function.
\subsection*{Two-step lookahead expected oracle value 
(\texttt{2LEY}) acquisition function}
\label{ss:r2LEY}
\begin{figure}
    \centering
    \includegraphics[width=\textwidth, trim=0.5cm 4cm 0.5cm 6cm,clip]{Sequential_2-step_Lookahead_Strategy.pdf}
    \caption{Recursive two-step lookahead method for time-dependent problems.
    Here, we start with $n$ observations up to $t_n$ and have $q-n$ 
   
    remaining
    observations to make, where $q$ is the total budget. We sequentially select
    one point $\mathbf{x}_{n+i}$ at each step that maximizes the expected value
    function evaluated at $T$; this is denoted by the black dashed lines
    at $t_{n+i}$ and $T$, whereas the gray dashed lines denote the
    remaining decisions to be made. We repeat the same process until the
    decision is made at the penultimate step $i=q-1$. The final decision at $T$
    is made by maximizing the value function at $T$.}
    \label{f:method}
\end{figure}
As an illustration, we set $h$ to be the identity function, and hence $\upsilona{n}{1}(\mathbf{x}, T) = \mu_{n,1}(\mathbf{x},
T)|\bm{\mathcal{D}}_{n,1}$.
We define our two-step lookahead expected oracle value (\texttt{2LEY}) acquisition function as
\begin{equation}
    \begin{split}
    \alpha_{2LEY}(\mathbf{x}^t_{n+1}) := & \mathbb{E}_n \left[  \underset{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}{\max}~\left.\mu_{n,1}(\tilde{\mathbf{x}}, T) \right| \bm{\mathcal{D}}_{n,1} \right] \\
        =& \mathbb{E}_n \left[ \mu^*_{n,1}( \mathbf{x}^t_{n+1}, T) \right],
    \end{split}
    \label{e:2LEY_1}
\end{equation}
where $\mu^*_{n,1} = \mu_{n,1}(\mathbf{x}^*, T)$ with $\mathbf{x}^*$ being a maximizer of $\mu_{n,1}(\mathbf{x},
T)|\bm{\mathcal{D}}_{n,1}$ and 
$y_{n+1}$ is a draw from $\mathcal{N} \left(\mu_n(\mathbf{x}^t_{n+1}), \sigma_n^2(\mathbf{x}^t_{n+1}) \right)$. 
The dependence of the second line of \eqref{e:2LEY_1} on
$\mathbf{x}^t_{n+1}$ can
be seen by realizing that
\begin{equation}
    \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T)  = \mathbf{k}_{n+1}^\top \mathbf{K}^{-1}_{n+1} [\mathbf{y}_n^\top, y_{n+1}]^\top,
\end{equation}
where $ \mathbf{k}_{n+1} = [ k\left((\mathbf{x}^*,T), (\mathbf{x}_1, t_1)\right), \ldots,
k\left((\mathbf{x}^*,T), (\mathbf{x}_n, t_n)\right), k\left((\mathbf{x}^*,T),
(\mathbf{x}_{n+1},t_{n+1})\right)]^\top$, and $\mathbf{K}_{n+1, ij} = k\left((\mathbf{x}_i, t_i), (\mathbf{x}_j, t_j) \right);~\forall i=1,\ldots,n+1, ~ j=1,\ldots,n+1$. The expectation in \eqref{e:2LEY_1} is approximated as
\begin{equation}
    \alpha_{2LEY}(\mathbf{x}^t_{n+1}) = \mathbb{E}_n \left[ \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right] \approx \frac{1}{N} \sum_{j=1}^N  \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) | y^j_{n+1}.
    \label{e:2LEY_MC}
\end{equation}
The gradient of $\mu_{n,1}^*(\mathbf{x}^t_{n+1}, T)$ with respect to $\mathbf{x}^t_{n+1}$ is given by
\begin{equation}
    \begin{split}
        \nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) = &\frac{\partial \mathbf{k}_{n+1}^\top}{\partial \mathbf{x}_{n+1}} \mathbf{K}^{-1}_{n+1} + \mathbf{k}_{n+1}^\top \frac{\partial \mathbf{K}^{-1}_{n+1}}{\partial \mathbf{x}_{n+1}} \\
        = &\frac{\partial \mathbf{k}_{n+1}^\top}{\partial \mathbf{x}_{n+1}} \mathbf{K}^{-1}_{n+1} +
        \mathbf{k}_{n+1}^\top~\mathbf{K}^{-1}_{n+1} \frac{\partial \mathbf{K}_{n+1}}{\partial \mathbf{x}_{n+1}}\mathbf{K}^{-1}_{n+1},
    \end{split}
    \label{e:2LEY-grad_a}
\end{equation}
where $ \frac{\partial \mathbf{K}_{n+1}}{\partial \mathbf{x}_{n+1}}$ is a matrix of
elementwise derivatives and the second line follows from a well-known lemma on the derivative of a matrix inverse~\cite[pp.~201--202]{rasmussen:williams:2006}.

With a continuously differentiable kernel $k(\cdot, \cdot)$, we have that
\begin{equation}
    \nabla_{\mathbf{x}} \alpha_{2LEY}(\mathbf{x}^t_{n+1}) = \nabla_{\mathbf{x}} \mathbb{E}_n \left[ \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right] = \mathbb{E}_n \left[\nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right],
    \label{e:2LEY_gr_ex_interchange}
\end{equation}
where the interchange of the gradient and expectation operators is via \cref{thm:gradexp}. Then the gradient is approximated as 
\begin{equation}
    \nabla_{\mathbf{x}} \alpha_{2LEY}(\mathbf{x}^t_{n+1}) = \mathbb{E}_n \left[\nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) \right] \approx \frac{1}{N} \sum_{j=1}^N \nabla_{\mathbf{x}} \mu_{n,1}^*(\mathbf{x}^t_{n+1}, T) | y^j_{n+1}.
    \label{e:2LEY_grad_MC}
\end{equation}
Furthermore, \eqref{e:2LEY_MC} and
\eqref{e:2LEY_grad_MC} are unbiased estimators for $\alpha_{2LEY}$ and $\nabla_{\mathbf{x}}
\alpha_{2LEY}$, respectively; see \cref{thm:gradexp}. This gradient can then be used in a gradient-based optimizer to
maximize our two-step lookahead acquisition function in an efficient manner.

Similarly, the two-step lookahead  extensions for
\texttt{EI}, \texttt{PI}, and \texttt{UCB} are obtained by picking the
corresponding value functions as $\alpha_{\textrm{EI}}(\mathbf{x},T)$, $\alpha_{\textrm{PI}}(\mathbf{x},T)$ and $\alpha_{\textrm{UCB}}(\mathbf{x},T)$ and are
named respectively, $\alpha_{2LEI}$, $\alpha_{2LPI}$, and $\alpha_{2LUCB}$. Our
overall algorithm recursively maximizes the two-step acquisition function until
the penultimate step and maximizes the actual value function at the final
step. The overall algorithm is graphically depicted in \cref{f:method} and is presented in
Algorithm~\ref{a:r2LBO}. We refer to Algorithm~\ref{a:r2LBO} with the value functions chosen as $\mu_{n,1}(\mathbf{x}, T)$, $\alpha_{EI}$, $\alpha_{PI}$, and $\alpha_{UCB}$ as \texttt{r2LEY}, \texttt{r2LEI}, \texttt{r2LPI}, and \texttt{r2LUCB}, respectively.


\begin{algorithm2e}[H]
\SetAlgoLined
\textbf{Given:} $\bm{\mathcal{D}}_n$,  
total budget $q$, 
schedule $\lbrace t_{n+1},\ldots, t_{q}:=  T \rbrace$, and GP hyperparameters $\bm{\Omega}$ \\
\KwResult{$(\mathbf{x}_T, \hat{y}_T)$}
  \For{$i=n+1, \ldots, q-1$, }{
  Define value function $\upsilon^T_{i-1, 1}|\bm{\mathcal{D}}_{i-1, 1}$ \\
  Find $\mathbf{x}_i = \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\operatornamewithlimits{arg\,max}}~ \hat{\alpha}_2(\mathbf{x}, t_i)$, \qquad (acquisition function maximization)\\
  \quad with $\hat{\alpha}_2$ and $\nabla_{\mathbf{x}}~\hat{\alpha}_2$ evaluated via \eqref{e:twostep_mc} and \eqref{e:twostep_T_grad_mc} resp. \\
    Observe $\hat{y}_i$ = $f(\mathbf{x}_i, t_i) + \epsilon_i$\\ 
    Append $\bm{\mathcal{D}}_i = \bm{\mathcal{D}}_{i-1} \cup \lbrace (\mathbf{x}_i, t_i), \hat{y}_i \rbrace$\\ 
    Update GP hyperparameters $\bm{\Omega}$ \\
 }
 Find $\mathbf{x}_T \in \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\operatornamewithlimits{arg\,max}}~\upsilon_{q}(\mathbf{x}, T)|\bm{\mathcal{D}}_{q-1}$ \qquad (value function maximization)\\
 Observe $\hat{y}_T$ = $f(\mathbf{x}_T, T) + \epsilon_T$\\ 
 \caption{Recursive two-step lookahead Bayesian optimization} 
 \label{a:r2LBO}
\end{algorithm2e}











\section{Theoretical properties}
\label{s:theory}

We now present the theoretical foundation of our method and a few properties.  
  Let $\mathcal{Y}$
  be the support of $y_{n+1}$, whose density is $p(y_{n+1})$. As before, let the
  domain of $\mathbf{x}^t_{n+1}$ be $\bm{\mathcal{X}}$. 




We make the following remark about the differentiability of our acquisition function; this justifies the interchange of the gradient and expectation operators.

\begin{remark}
We write the inner term of our two-step lookahead acquisition function in \eqref{e:E_twostep} as $\max_{\mathbf{x}\in\bm{\mathcal{X}}} \mathbb{E}_{n,1}h(y)$, where we omit the arguments of $y$ and $\boldsymbol{\omega}$ for the sake of conciseness. The GP posterior mean $\mu(\mathbf{x}, t)$ and standard deviation $\sigma(\mathbf{x}, t)$ are continuously differentiable with the sufficient condition that we choose a continuously differentiable prior covariance kernel $k$~\cite{smith1995differentiation, wang2016parallel}. We then use the reparameterization trick~\cite{wilson2017reparameterization,
wilson2018maximizing} to write the GP sample path as $y(\mathbf{x}, T) := \mu(\mathbf{x}, T) + \sigma(\mathbf{x},
T)\times \gamma$, where $\gamma \sim \mathcal{N}(0, 1)$, which makes the differentiability of our acquisition function transparent. 


From a purely implementation standpoint, we use the gradient $\nabla
\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_{n, 1}~h(y)}$ if it exists and a
\emph{subgradient}~\cite{boyd2004convex} otherwise. In our implementation, the
gradients are computed by backpropagating the output of the
$\max()$ operation; in \Cref{app:grad_nd} we provide examples of how (sub)gradients are calculated for nonsmooth functions.
\end{remark}








Let $\mathbf{x}^*$ be a maximizer of $\mathbb{E}_{n,1} h(y(\mathbf{x}, T); \boldsymbol{\omega}) = \upsilon_{n,1}(\mathbf{x}, T)|\bm{\mathcal{D}}_{n,1}$. In what follows, we write
  $g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) := \upsilon_{n,1}(\mathbf{x}^*, T)|\bm{\mathcal{D}}_{n,1}$ to denote the maximum value of the value function at $T$ after having updated $\bm{\mathcal{D}}_n$ with $y_{n+1}$, and $\nabla_{\mathbf{x}}g^* := \partial g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})/\partial \mathbf{x}^t_{n+1} $. Note again that the dependence of $g^*$ on $\mathbf{x}^t_{n+1}$ is due to $\bm{\mathcal{D}}_{n,1}$.
  
\begin{theorem}[Interchange of gradient and expectation operators] 
\label{thm:gradexp}
 Let $\nabla_{\mathbf{x}}g^*$ exist, and let $g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})\times
  p(y_{n+1})$ and $\nabla_{\mathbf{x}}g^* \times
  p(y_{n+1})$ be continuous on $\bm{\mathcal{X}} \times \mathcal{Y}$. 
  Further assume
  that there exist functions $q_0$ and $q_1$ such that
\begin{equation*}
\begin{array}{rl}
  |g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) \times p(y_{n+1})| &\leq q_0(y_{n+1}) \\
  \|
  \nabla_{\mathbf{x}} g^* \times  p(y_{n+1})\| &\leq q_1(y_{n+1})
  \end{array}
  \qquad \forall (\mathbf{x}^t_{n+1}, y_{n+1}) \in \bm{\mathcal{X}} \times \mathcal{Y}, 
  \end{equation*}
where
  $\int_{\mathcal{Y}}q_1(y_{n+1}) dy_{n+1} < \infty$ and
  $\int_{\mathcal{Y}}q_2(y_{n+1}) dy_{n+1} < \infty$. Then,
\[ \nabla_{\mathbf{x}} \mathbb{E}_j \left[ g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) \right] = \mathbb{E}_n \left[\nabla_{\mathbf{x}} g^*(\mathbf{x}^t_{n+1}, T, y_{n+1}) \right].\]

Furthermore, a realization of $\nabla_{\mathbf{x}} g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})$ yields
an unbiased estimate of the gradient of $\alpha_2(\mathbf{x}_{n+1}^t)$ in \eqref{e:E_twostep}.
\end{theorem}

\begin{proof}
Using $g^*(\mathbf{x}^t_{n+1})$ to denote  $g^*(\mathbf{x}^t_{n+1}, T, y_{n+1})$ for the sake of brevity, we have
\begin{equation}
\begin{split}
    \nabla_\mathbf{x} \alpha_2(\mathbf{x}^t_{n+1}) := \frac{\partial}{\partial \mathbf{x}^t_{n+1}} \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}) \right] &= \underset{h\rightarrow 0}{\textnormal{lim}}~\frac{1}{h} \left \lbrace \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}+h) \right] -  \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}) \right] \right \rbrace \\
    &= \underset{h\rightarrow 0}{\textnormal{lim}}~ \left \lbrace \mathbb{E}_n \frac{1}{h} \left[ g^*(\mathbf{x}^t_{n+1}+h) - g^*(\mathbf{x}^t_{n+1}) \right] \right \rbrace \\
    &= \underset{h\rightarrow 0}{\textnormal{lim}}~ \left \lbrace \mathbb{E}_n \frac{\partial g^*(\bar{\mathbf{x}}^h_{n+1})}{\partial \mathbf{x}_{n+1}}  \right \rbrace ,\\
\end{split} 
\label{e:grad_limit}
\end{equation}
where $\bar{\mathbf{x}}^h_{n+1} = \lambda \mathbf{x}^t_{n+1} + (1-\lambda)(\mathbf{x}^t_{n+1} + h)$ for some
$\lambda \in [0,1]$ and the last line follows from the first-order mean value
theorem.
We bring the limit inside the integral by Lebesgue's dominated convergence theorem to obtain
\begin{equation}
    \frac{\partial}{\partial \mathbf{x}^t_{n+1}} \mathbb{E}_n \left[ g^*(\mathbf{x}^t_{n+1}) \right] =  \mathbb{E}_n~\underset{h\rightarrow 0}{\textnormal{lim}} \frac{\partial g^*(\bar{\mathbf{x}}^h_{n+1})}{\partial \mathbf{x}^t_{n+1}}   = \mathbb{E}_n \frac{\partial g^*(\mathbf{x}^t_{n+1})}{\partial \mathbf{x}^t_{n+1}}.
    \label{e:LDC}
\end{equation}
Since \eqref{e:grad_limit} holds, a realization $\partial g^*(\mathbf{x}^t_{n+1}) / \partial \mathbf{x}^t_{n+1}$ yields an unbiased estimate of the gradient $\nabla_\mathbf{x} \alpha_2(\mathbf{x}^t_{n+1})$.
\end{proof}



\begin{proposition}[Special case when $m=2$ and $\upsilon^T = \mu(\mathbf{x}, T)$] 
When $m=2$, our \texttt{2LEY} acquisition function is equivalent to the knowledge gradient
(KG)~\cite{frazier2008knowledge, wu2016parallel} acquisition function offset by a constant factor.

\begin{proof}
The KG acquisition function, in the context of our problem, is given as
\begin{equation}
    \alpha_{KG}(\mathbf{x}^t_{n+1}) = \mathbb{E}_{n}\left[ \underset{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}{\max}~ \mu_{n,1}(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] - \underset{\mathbf{x} \in \bm{\mathcal{X}}}{\max}~ \mu_{n}(\mathbf{x},T).
    \label{e:KG_time}
\end{equation}

With $m=2$, the \texttt{2LEY} acquisition function at $\mathbf{x}_{n+1}$ is defined as follows:
\begin{equation}\label{eq:nm1}
    \begin{split}
    \alpha_{2LEY}(\mathbf{x}^t_{n+1}) & = \int_{y_{n}} \left[ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}} ~\mathbb{E}_{n,1} ~y(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] p(y_{n}|\bm{\mathcal{D}}_{n})~dy_{n} \\
    & = \int_{y_{n}}\left[  \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~ \mu_{n,1}(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] p(y_{n}|\bm{\mathcal{D}}_{n})~dy_{n} \\
    & = \mathbb{E}_{n}\left[ \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~ \mu_{n,1}(\tilde{\mathbf{x}},T)|\bm{\mathcal{D}}_{n,1} \right] \\
    & = \alpha_{KG}(\mathbf{x}^t_{n+1}) + \max_{\tilde{\mathbf{x}} \in \bm{\mathcal{X}}}~\mu_{n}(\tilde{\mathbf{x}},T), 
    \end{split}
\end{equation}
where the last line follows from \cref{e:KG_time}. Note that $\alpha_{2LEY}$ and $\alpha_{KG}$ share the same maximizer. 
\end{proof}
\label{prop:special_cases}
\end{proposition}

\begin{proposition}[Special case when $|T-t| \gg\theta_t$] 
With the prior assumption $f(\mathbf{x},t)\sim \mathcal{GP}(0, k(\cdot,\cdot))$, where $k(\cdot,\cdot) := k_{\mathbf{x}}(\cdot,\cdot; ~\theta_{\mathbf{x}})\times k_t(\cdot,\cdot; ~\theta_{t})$, with $k_{\mathbf{x}}$ and $k_t$ being squared-exponential kernels,
our two-step lookahead acquisition function approaches $\alpha_2(\mathbf{x}^t_{n+1}) = \max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\upsilon_n(\tilde{\mathbf{x}}, T)}$ as $|T - t_{n+1}| \rightarrow \infty$. In other words, far away from the target horizon $T$, $\alpha_2(\mathbf{x}^t_{n+1})$ approaches a constant function.


\begin{proof}
When $|T-t_{n+1}| \gg \theta_t$, then $k_t(t_{n+1}, T) \approx 0$; this yields $k\left( (\mathbf{x},t_{n+1}), (\mathbf{x},T) \right) \approx 0$. Let $\mu_n$ and $\sigma_n^2$ be the posterior mean and variance given $\bm{\mathcal{D}}_n$, respectively. Given another observation 
$y_{n+1}$ drawn from the GP $Y_n$ at $(\mathbf{x}, t_{n+1})$, where $|T-t_{n+1}| \gg \theta_t$, the updated posterior mean and variance are given by
\begin{equation}
    \begin{split}
        \mu_{n,1} =& \mathbf{k}^\top_{n+1}\mathbf{K}^{-1}_{n+1}[\mathbf{y}^\top_n, y_{n+1}]^\top  \approx [\mathbf{k}^\top_n, 0] \begin{bmatrix}
        \mathbf{K}_n & \mathbf{0} \\
        \mathbf{0} & 1
        \end{bmatrix}^{-1}[\mathbf{y}^\top_n, y_{n+1}]^\top = \mu_n \\
        \sigma^2_{n,1} =& 1-\mathbf{k}^\top_{n+1}\mathbf{K}^{-1}_{n+1}\mathbf{k}_{n+1} \approx 1 - [\mathbf{k}^\top_n, 0] \begin{bmatrix}
        \mathbf{K}_n & \mathbf{0} \\
        \mathbf{0} & 1
        \end{bmatrix}^{-1}[\mathbf{k}^\top_n, 0]^\top = \sigma_n^2,
    \end{split}
\end{equation}
where we have assumed, without loss of generality, that $k((\mathbf{x},t), (\mathbf{x},t))=1$ and have used the matrix inversion lemma.
Therefore, writing our acquisition function using the reparameterization trick, we see that
\begin{equation}
    \begin{split}
    \alpha_2(\mathbf{x}^t_{n+1}) =& \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_{n, 1}~h(y(\tilde{\mathbf{x}}, T))} \right]\\ =& \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\mu_{n,1}(\tilde{\mathbf{x}}, T) + \sigma_{n,1}(\tilde{\mathbf{x}}, T)\gamma)} \right]\\ \approx& 
    \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\mu_{n}(\tilde{\mathbf{x}}, T) + \sigma_{n}(\tilde{\mathbf{x}}, T)\gamma)} \right]\\
    =& \max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\upsilon_n(\tilde{\mathbf{x}}, T)},
    \end{split}
\end{equation}
and thus, in the limit $|T - t_{n+1}| \rightarrow \infty$,
$\alpha_2(\mathbf{x}^t_{n+1}) = \max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\upsilon_n(\tilde{\mathbf{x}}, T)}$.
\end{proof}
\label{prop:t_llt_T}


\begin{remark}
    Following \Cref{prop:t_llt_T}, one can  see that if $k_t(t_{n+1}, T) \approx 0$, then $k_t(t_i, T) \approx 0,~\forall i=1,\ldots, n$ since $t_{i} > t_{i-1}$. This further leads to $\mathbf{k}_{n+1} \approx \mathbf{0}$ and hence $\mu_{n,1} \approx 0$ and $\sigma_{n,1}^2 \approx 1$. This results in (for $|T-t_{n+1}| \gg \theta_t$),
    \begin{equation}
    \begin{split}
    \alpha_{2}(\mathbf{x}^t_{n+1}) =& \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\mu_{n}(\tilde{\mathbf{x}}, T) + \sigma_{n}(\tilde{\mathbf{x}}, T)\gamma)} \right] \\
    \approx & \mathbb{E}_n \left[\max_{\tilde{\mathbf{x}}\in\bm{\mathcal{X}}}{\mathbb{E}_\gamma~h(\gamma)} \right] \\
    =& c,
    \end{split}
    \label{e:2_inner}
\end{equation}
where $c$ is some constant. Suppose the $(n+1)$th decision is made by maximizing $\alpha_{2}(\mathbf{x}^t_{n+1})$ with a multistart local optimization procedure started from a set of starting points $\{\mathbf{x}\} \in \bm{\mathcal{X}}$, where $\{\mathbf{x}\}$ are sampled from distribution $p_\mathbf{x}$. Then, in the limit $|T - t_{n+1}| \rightarrow \infty$, decisions made via our lookahead acquisition function are equivalent to samples drawn randomly from $p_\mathbf{x}$.\footnote{Note that this assumes that the local optimizer seeks a critical point and returns the starting point as the maximizer.} 

We take advantage of this property and set $p_\mathbf{x} = \mathcal{U}(\bm{\mathcal{X}})$, where $\mathcal{U}$ is the uniform distribution. In practice, this could result in a decisions that \emph{explore} $\bm{\mathcal{X}}$ when the $T$ is far away and hence can lead to improved learning of $f(\mathbf{x},t)$, thus facilitating a more accurate prediction of the maximum at $T$.
\end{remark}

\label{prop:special_cases2}
\end{proposition}





\section{Numerical experiments}
We demonstrate our proposed approach on synthetic test cases and a quantum optimal control problem. 

We compare our method with
widely used myopic approaches in BO, namely,  \texttt{EI}, \texttt{PI}, and \texttt{UCB}.
We refer to ``\texttt{EImumax}'' and ``\texttt{PImumax}'' as the \texttt{EI} and \texttt{PI} strategies,
respectively, with the target $\xi$ set as the maximum of the GP posterior mean at the
current step, that is, $\underset{\mathbf{x} \in \bm{\mathcal{X}}}{\max}~ \mu(\mathbf{x}, t)$. We set 
the confidence parameter for \texttt{UCB}  to $\beta=2$. Additionally, we
 compare our method with ``\texttt{mumax},'' which selects
points by maximizing the GP posterior mean at the current time. We also make direct comparisons
between the myopic acquisition functions and their lookahead counterparts.
Moreover, we compare our method with a strategy that selects points uniformly at random
from $\mathbf{x} \sim \mathcal{U}(\bm{\mathcal{X}})$; we call this acquisition function ``\texttt{Random}.'' 

Our
metric for comparison is the oracle value $f(\mathbf{x}_T, T)$ at the target time $T$. Each query to $f(\mathbf{x}, t)$ includes additive noise with mean zero and variance $\sigma_\epsilon^2 = 10^{-3}$. Each repetition of each experiment is provided with \emph{starting} samples that include oracle evaluations at a set of points distributed with uniform spacing in $[0, t_{\textnormal{start}}]$ with each point sampled uniformly at random in $\bm{\mathcal{X}}$, where $t_\textnormal{start}<T$ is arbitrarily set for each experiment. Within each repetition the same starting samples are given to all acquisition functions. 


\subsection{Synthetic one-dimensional test cases}
\label{ss:synthetic_1d}
Our synthetic test functions take the
form $f(\mathbf{x}, t) = f_{\mathbf{x}}(\mathbf{x}) + f_{\mathbf{x} t}(\mathbf{x}, t)$ or  $f(\mathbf{x}, t) = f_{\textnormal{ref}}(\mathbf{x} (t)) $ 
in order to  induce time dependence to the
maximizer of $f(\mathbf{x},t)$.

We first consider quadratic one-dimensional functions $f_{\mathbf{x}}$
over the domain $\bm{\mathcal{X}} = [0,1]$ with
\begin{equation} 
f_{\mathbf{x}}(\mathbf{x}) = -4 (\mathbf{x} - 0.5)^2.
\label{eq:1dq}
\end{equation}
The time-dependent component takes the following four forms:
\begin{equation*}
    \begin{split}
    \textnormal{Quadratic-a}: \quad f_{\mathbf{x} t} = &\sin \left(\pi(\mathbf{x}+t)\right) + \cos \left(\pi(\mathbf{x}+t)\right) \\
    \textnormal{Quadratic-b}: \quad f_{\mathbf{x} t} = &\sin \left(\pi(\mathbf{x} t)\right) + \textnormal{cos}\left(\pi(\mathbf{x} t)\right) \\
    \textnormal{Quadratic-c}: \quad f_{\mathbf{x} t} = &\sin \left(\pi(\mathbf{x}[t-3]^+)\right) +
    \textnormal{cos}\left(\pi(\mathbf{x}[t-3]^+)\right)\\
    \textnormal{Quadratic-d}: \quad f_{\mathbf{x} t} = &2\mathbf{x}\textnormal{sin}(t) - \textnormal{sin}^2(t),  
    \end{split}
\end{equation*}
where $[t-3]^+ \equiv \max(0, t-3)$ is specified to induce movement of
$\mathbf{x}_t^*$ for $t\geq 3$. 
\cref{f:1d_synthetic} shows the trajectories of $\mathbf{x}_t^*$, the maximizer at time $t$, for each of these test cases.
\begin{figure}
    \centering
    \includegraphics[trim={100 0 100 0}, clip, width=1\linewidth]{quadratics.pdf}
    \caption{One-dimensional synthetic test cases. The dashed line represents the trajectory of $\mathbf{x}_t^*$.}
    \label{f:1d_synthetic}
\end{figure}
These test cases are deliberately designed to reflect a variety of situations, such as discontinuous change of $\mathbf{x}^*_t$ (Quadratic-a and Quadratic-b) and sudden dynamics (Quadratic-c). 

\cref{f:quadratic_exp} compares the points selected by our lookahead approach with various value functions and the corresponding myopic approaches. In each plot
20 repetitions 
are overlaid, with samples in $\bm{\mathcal{X}}\times[0, 1]$ used to start all of the methods;  white
circles denote the points chosen by each method at times $t>1$. 
The lookahead BO does not try to track a moving maximizer, which results in points being more spread out in a space-filling fashion than with myopic acquisition functions. 
This approach is particularly  beneficial in handling oracles whose maximizers can go through a discontinuous change (e.g., Quadratic-a and Quadratic-b).
In the case of Quadratic-c, we demonstrate (see further details in \cref{sf:c}) that our algorithm is able
to handle oracles where $\mathbf{x}_t^*$ may undergo sudden dynamics ($\approx t=3.0$ in
this case). 
For Quadratic-d, where the challenges of the other
one-dimensional cases do not exist, the algorithm performs much better in
predicting $\mathbf{x}_T^*$.
\begin{figure}
    \centering
    \begin{subfigure}{.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip, width=1\linewidth]{qa_mumax_2LEY.pdf}
        \caption{Quadratic-a}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip,width=1\linewidth]{qb_EI_2LEI.pdf}
        \caption{Quadratic-b}
    \end{subfigure} \\
    \begin{subfigure}{.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip,width=1\linewidth]{qc_PI_2LPI.pdf}
        \caption{Quadratic-c}
    \end{subfigure}%
    \begin{subfigure}{.5\textwidth}
        \centering
        \includegraphics[trim={30 0 50 0}, clip,width=1\linewidth]{qd_UCB_2LUCB.pdf}
        \caption{Quadratic-d}
    \end{subfigure}    
    \caption{Comparison of our two-step lookahead (right) and myopic (left) acquisition functions on the one-dimensional synthetic test cases. The plot shows
    an overlay of 20 repetitions of the algorithm with (shared) random starting points. The
    white circles denote the decision made at each step,
    and the contours in the background
    represents the true noise-free oracle $f(\mathbf{x}, t)$}
    \label{f:quadratic_exp}
\end{figure}

    


\cref{f:quad_myopic_v_lookahead} compares the time history of the average oracle value (over 20 replications
with randomized starting points) resulting from decisions made via the myopic and lookahead acquisition functions. We present results for only
Quadratic-a and Quadratic-c for the sake of brevity. One can see that
the lookahead acquisition functions (dashed line plots) make decisions with
lower oracle value during the early rounds. While the acquisition function per se does not
\emph{seek} to make decisions with low current oracle value, the consequence of
looking ahead at $T$ is that decisions made at current $t$ may incur a low oracle value. This is a fundamental distinction between our proposed lookahead and the myopic acquisition functions. Also, the average oracle value at target
time $T$ (denoted by the empty circle) is consistently higher than the myopic
counterparts (denoted by the square). Furthermore, the predictions by the
lookahead acquisition functions are with higher confidence, as visualized by
the error bars (one standard deviation in the predictions) shown at $T$; note that we highlight the error bars for the
lookahead acquisition functions in black.

The average oracle values at $T$
predicted by all the acquisition functions are somewhat consistent, with the
\texttt{r2LEY} being one of the most competitive, including the rest of the test cases.
For this reason, moving forward we predominately
use the \texttt{r2LEY} as the representative lookahead acquisition function.
\begin{figure}
    \centering
    \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=.8\linewidth]{qa_myopic_v_lookahead.pdf}
        \caption{Quadratic-a}
        \label{sf:a}        
    \end{subfigure}\\
        \begin{subfigure}{1\textwidth}
        \centering
        \includegraphics[width=.8\linewidth]{qc_myopic_v_lookahead.pdf}
        \caption{Quadratic-c}
        \label{sf:c}        
    \end{subfigure}
\caption{Average of 20 replications with random
starting points of the myopic (lookahead) acquisition function on 
one-dimensional test cases. The final point at $T=4$ is highlighted with square (circle) symbols and  variability at $T$ is visualized via
colored (black) error bars.
The vertical dashed line represents the time until which $n=15$ starting points
are collected and the total budget is set to $q=45$.}
\label{f:quad_myopic_v_lookahead}
\end{figure}


\subsection{Synthetic higher dimensional test cases}
\label{ss:synthetic_hd}
We now demonstrate our method on higher-dimensional test functions, up to $d=6$.
We begin by presenting a modified Griewank function
($d=2$)~\cite{simulationlib}, where we multiply the original
Griewank function $f_{\textnormal{ref}}$ by a Gaussian weight function in order to create a
unique global maximum. Then, we modify the input $\mathbf{x}
\in \bm{\mathcal{X}}=[-5,5]^2$ with a time-dependent rotation to yield the function
\begin{equation}
    f(\mathbf{x},t) = f_{\textnormal{ref}}(\mathbf{x}(t)) \times \textnormal{exp}~\left( - \| \mathbf{x} - \mathbf{x}_0 \|^2/\ell\right), 
\end{equation}
where $\mathbf{x}_0 = [3, 0]$, $\ell \in \mathbb{R}_+$  is a scaling
constant\footnote{We set
$\ell = 160$ to ensure the weighting is somewhat mild and does not
significantly change the nature of the true Griewank function.} and $\mathbf{x}(t) = R(\zeta_t) \times \mathbf{x},$
where
\[ R(\zeta_t) = 
\begin{bmatrix}
\cos(\zeta_t) & -\sin(\zeta_t) \\
\sin(\zeta_t) & \cos(\zeta_t)
\end{bmatrix} \]
is a $2\times 2$ rotation matrix and $\zeta_t = \pi t /4$. Time snapshots of
this test function are shown in \cref{f:mod_griewank} and illustrate a 
counterclockwise rotation of the unique global maximizer with time. As before, all oracle observations include additive mean zero noise with variance $\sigma_\epsilon^2 = 10^{-3}$.

The performance of all the myopic acquisition functions is compared against \texttt{r2LEY} in \cref{f:mod_griewank_result}, where the
contour plot at $T$ is shown with the predictions from 20 repetitions of each
acquisition function (black $\times$ symbols) overlaid. Note that the total
budget is fixed at $q=90$, where $n=60$ is used in $\bm{\mathcal{X}} \times [2,3]$ to
start the algorithm, and the target time is $T=4$.  The
lookahead acquisition function clearly predicts the global maximum at $T$ with greater
consistency than the rest do; all but two repetitions predicted $\mathbf{x}_T$ such that $\| \mathbf{x}_T - \mathbf{x}^*_T\| < 0.26$.

\begin{figure}[htb!]
    \centering
    \includegraphics[trim={50 50 50 50}, clip, width=.8\linewidth]{modified_griewank.pdf}
    \caption{Time snapshots of the modified Griewank function. Notice the unique global maximum at each time and the counterclockwise rotation of the global maximizer with time.}
    \label{f:mod_griewank}
\end{figure}
\begin{figure}[htb!]
    \centering
    \includegraphics[width=.8\linewidth]{mod_griewank_result.pdf}
    \caption{Predictions at $T=4.0$ on the modified Griewank test case. Symbols represent predictions at $T$ after 20 repetitions of each algorithm, and the contours represent the noise-free $f(\mathbf{x}, T)$. The total budget is set to $q=90$, where $n=60$ samples in $2.0\leq t \leq 3.0$ are used to start each algorithm.}
    \label{f:mod_griewank_result}
\end{figure}


We also demonstrate our method on the Hartmann-3d ($d=3$) and Hartmann-6d
($d=6$) test functions. In both cases, the Hartmann
function~\cite{simulationlib} is summed with $f_{\mathbf{x} t} (\mathbf{x}, t) = \mathbf{1}^\top [2\textnormal{sin}(t)\mathbf{x}
- \textnormal{sin}^2(t)\mathbf{1}]$ (where $\mathbf{1}$ is a vector of ones of length $d$) to induce movement of the maximum with time.
For Hartmann-3d, the total budget is set to $q=120$ with $n=90$
samples queried in $[0, 1]^3 \times [2,3]$ to start the algorithm. For
Hartmann-6d, we set $q=150$ with $n=120$ samples queried in $[0, 1]^3 \times
[2,3]$ to start the algorithm, and our target time is $T=4$. 
\Cref{f:hartmann} shows the time histories of $f(\mathbf{x}_i,t)$
averaged over 20 replications. The predictions at $T$ are highlighted by adding
a marker in the plot. In both  test functions, the
lookahead acquisition function results in a higher oracle value $f(\mathbf{x}_T, T)$ at the
target time $T$ compared with other myopic acquisition functions. 

\begin{figure}[htb!]
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{hart3d.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{hart6d.pdf}
    \end{subfigure}
    \caption{Hartmann-3d ($d=3$) and Hartmann-6d ($d=6$) test functions. Average of 20 repetitions. Left of vertical line corresponds to starting samples. The error bars for the oracle value at $T$ for the best two acquisition functions (\texttt{r2LEY} and \texttt{mumax}) shown are slightly shifted in time for better visibility. The \texttt{r2LEY} shows comparatively less variability.}
    \label{f:hartmann}
\end{figure}

\subsection{Quantum optimal control}
\label{ss:qoc}
We next compare methods for optimizing time-dependent oracles on a challenging real-world problem.
Optimally shaping electromagnetic fields to control quantum systems is a widespread
application~\cite{brif2010control}; similar problems appear when controlling molecular transformations relevant to chemical,
biological, and materials science applications. Quantum optimal control (QOC) is
a method to shape electromagnetic fields to steer a quantum system toward a
desired control target; see
\cref{f:quantum_chemistry} where the target is a specific outcome from a chemical reaction. 

\begin{figure}[tb!]
    \centering
    \includegraphics[width=0.5\linewidth]{Quantum_chemistry.PNG}
    \caption{Illustration of QOC to control chemical reactions. A femtosecond
    laser controlled by QOC is used to provide control over the selective
    dissociation of molecules. Figure source:~\cite{magann2020digital}}
    \label{f:quantum_chemistry}
\end{figure}

The problem we consider is based on a diatomic molecule \ce{HF}, whose vibrational energy mode is modeled as a nonrotating Morse oscillator on a 1D grid that varies with time and is illuminated with a laser
field to induce dissociation. The electromagnetic field is parameterized by a set of coefficients $\mathbf{x} \in \bm{\mathcal{X}} \subset \mathbb{R}^5$ that are the control parameters; more specifically, $\mathbf{x}$ contains the frequency components of the electromagnetic field. The digital quantum simulation of the Morse 
oscillator $f(\mathbf{x}, t)$ takes as input the control parameter $\mathbf{x}$ and a \emph{terminal} time $t
\in [0, T]$---which is not to be confused with our target time $T$---whose
 output is the observable control objective that needs to be minimized with respect to $\mathbf{x}$. 
The cost of each query of the simulation increases monotonically with the terminal time $t$, and
we are interested in finding the optimal control at a specific target time $T$. Since searching for the optimal control at $T$ is expensive,  we seek to predict 
 the optimal control at $T$ with cheaper evaluations of $f(\mathbf{x}, t)$ for $t<T$.

In this regard, we use BO to find the control parameters $\mathbf{x}_T^*$ that maximize the negative value of $f(\mathbf{x},T),
\forall \mathbf{x} \in \bm{\mathcal{X}}$. We begin by collecting observations of their system
$\bm{\mathcal{D}}_n= \{(\mathbf{x}_i, t_i), \hat{y}_i \},~i=1,\ldots,n$, which we use to fit
a GP that emulates $f$ in $\bm{\mathcal{X}} \times {\mathcal{T}}$. 
Then, we sequentially
make observations $\{(\mathbf{x}_i, t_i), \hat{y}_i \},~i=n+1,\ldots,q$, where 
each $\mathbf{x}_i$ is selected by maximizing our proposed two-step lookahead acquisition function. 
\Cref{f:quantum_problem} provides a schematic of
this problem. To conform to our problem setting, we allow evaluations only at increasing values of $t$. 
See~\cite{magann2020digital} for more information about this problem. 


We set the total budget to $q=180$, where $n=150$ samples are queried over $4000
\leq t \leq 5000$ to start the algorithm and $T=6000$. The time histories of
each acquisition function for the QOC problem are shown in
\cref{f:quantum_results}, where we include results from \texttt{r2LEY} and \texttt{r2LPI} among the lookahead approaches. 
The best-known  maximum $f(\mathbf{x}^*,
T)=0.0$ at $T=6000$, which is shown in the figure as the black square. In terms of the average performance, the
lookahead \texttt{r2LEY} outperforms all other methods, with the \texttt{mumax}
being a  close competitor, as shown in the top-left of \Cref{f:quantum_results}. We visualize the variability in the predictions at $T$ via boxplots in the top right of \Cref{f:quantum_results}, where we observe that the \texttt{r2LEY} has the highest median oracle value to interquartile distance ratio. This shows that although \texttt{r2LEY} does not result in the best median oracle value, it provides the most confident prediction, while still being competitive with \texttt{UCB} and \texttt{mumax}. 
Even though the lookahead method \texttt{r2LPI} does not necessarily outperform all the other myopic methods, it shows better performance---in terms of both average and median oracle value---compared with its myopic counterpart \texttt{PImumax}. This behavior is consistent with  what was observed with the one-dimensional test cases presented in \cref{ss:synthetic_1d}.



\begin{figure}[tb!]
    \centering
    \includegraphics[width=1\linewidth]{Quantum_problem3.pdf}
    \caption{Overview of the quantum optimal control problem recreated from \cite{magann2020digital}.
    The quantum system dynamics, driven by a field
    parameterized by $\mathbf{x}$, are simulated for the time range $[0, t]$ (where $t$
    is the terminal time), after which the control objective function
    $f(\mathbf{x},t)$ is calculated. A recent framework \cite{magann2020digital} is proposed
    for efficiently evaluating the control
    objective function 
    by using a future quantum computer.
    During an initial sample generation phase, the control objective function
    is evaluated $n$
    times to generate initial samples $\bm{\mathcal{D}}_n$ that are used by the GP to
    learn $f(\mathbf{x}, t),~ \mathbf{x}\in \bm{\mathcal{X}}, t\in [0, t]$. Then, in the optimization
    phase, the trained GP 
    is used in the BO framework to find the optimal
    control at target time $T$. The inputs at $(\mathbf{x}_i,t_i),~i=,n+1,\ldots,q$ are
    chosen at \emph{optimal} locations as guided by the BO acquisition function
    optimization.}
    \label{f:quantum_problem}
\end{figure}
\begin{figure}
    \centering
    \begin{subfigure}[t]{.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{quantum.pdf}
        \caption{Average $f(\mathbf{x},t)$ over 20 repetitions with randomized starting points. Left of vertical line denotes starting points.}
    \end{subfigure}%
        \begin{subfigure}[t]{.5\textwidth}
        \centering
        \includegraphics[width=1\linewidth]{quantum_boxplot.pdf}
        \caption{Variability in the oracle value at $T$, $f(\mathbf{x}, T)$. }
    \end{subfigure}\\
   
   
   
   
   
   
   

    \caption{Demonstration on the quantum optimal control problem.}
    \label{f:quantum_results}
\end{figure}






















\subsection{Implementation details}
We implemented our recursive lookahead Bayesian optimization framework
in PyTorch~\cite{pytorch_bib} using the GP and BO libraries GPyTorch~\cite{gardner2018gpytorch} and
BoTorch, respectively~\cite{balandat2019botorch}. For
finite-horizon time-dependent BO with limited budget, the lookahead
acquisition functions outperform their myopic counterparts in securing the best
oracle value at the horizon. This improved accuracy comes at the cost of additional
computational overhead in terms of optimizing the acquisition functions. Our
Monte Carlo approximation to the acquisition function and its gradients, as presented in Line 5 of 
Algorithm~\ref{a:r2LBO}, lends
itself well to be used with a gradient-based optimizer.
For further computational improvements, however, we
also adopt the \emph{one-shot} optimization of the lookahead acquisition
function as implemented in BoTorch. The approach formulates line 4 of Algorithm~\ref{a:r2LBO} as a deterministic optimization problem over the joint $\mathbf{x}$ and $\mathbf{x}'$, where $\mathbf{x}' := \{ \mathbf{x}^j\}_{j=1}^N$ are the candidate maximizers or \emph{fantasy} points. While this increases the dimensionality of the optimization problem from $d$ to $N+d$, we observe consistent results (with our Monte Carlo approach) with $N$ set to as small as 32, resulting in significant computational speedup. Moreover, the one-shot approach solves a deterministic optimization problem as opposed to a stochastic optimization problem with the Monte Carlo approach.
We include
both the Monte Carlo and one-shot approaches in our software and recommend the
latter as the dimensionality $d$ of the problem increases; all the results presented in this work are based on the one-shot approach. 
For the $d=5$
QOC problem, the acquisition function optimization takes approximately
30 seconds of wall-clock time when run in serial. Because many of the time-dependent oracles we are targeting require significantly more computational cost and time, we consider the per-iterate cost of our
approach to be trivial. Furthermore, the PyTorch framework allows
for leveraging advanced computation hardware such as graphics processing units, which may allow for further computational efficiency.

\section{Conclusions and future work}
Bayesian optimization with a lookahead acquisition
function outperforms well-known myopic acquisition functions in solving the finite-budget, finite-horizon, time-dependent maximization of an
expensive stochastic oracle.
Specifically, for such problems, the lookahead approach
makes probabilistically sound decisions that give better average-case
returns at the target horizon, while considering 
strategies that
lead up to the target horizon. In this regard, the lookahead acquisition
functions take advantage of the epistemic uncertainty---explained by the GP---in
the predictions about the oracle. In this work we introduce a generalized
framework that offers lookahead extensions to widely known myopic acquisition
functions. We demonstrate our method on 
illustrative 
synthetic test functions (that involve discontinuities in the
oracle maximizer) as well as a quantum optimal control application problem.
Because of the lack of analytically closed-form expressions for our lookahead
acquisition functions, their construction (and optimization) incurs more
computational overhead than do the myopic acquisition functions. Such an
overhead is insignificant, however, when optimizing expensive stochastic oracles that take, for example, hours of wall-clock time per evaluation. 
Overall, we demonstrate that our
proposed acquisition function results in more confident and accurate
predictions of the global maxima at the target horizon compared with Bayesian optimization using myopic acquisition functions.

The primary challenge with handling time-dependent optimization problems with
BO is that the effect of the dynamics of the oracle has to be learned by a GP model. When the
distance to the horizon is large, it is particularly challenging to train a
\emph{global} GP that emulates the oracle in $\bm{\mathcal{X}} \times \mathcal{T}$. A natural
alternative is to work with \emph{local} GPs within 
local
regions that
are dynamically updated. This is one avenue for future work that we
are considering.

A related problem is multifidelity
BO (see, e.g., \cite{song2019general, kandasamy2017multi, takeno2019multi}) in which one seeks to solve
$    \max_{\mathbf{x} \in \bm{\mathcal{X}} }~f(\mathbf{x})$ using observations of $\hat{f}(\mathbf{x}, s)$. The \emph{fidelity} is given by the
parameter $s\in[0,1]$, with $f(\mathbf{x}):= \hat{f}(\mathbf{x}, 1)$ corresponding to the
highest fidelity. Typically, there is a cost function $c(s)\in \mathbb{R}$, which is monotonically increasing in the fidelity
level $s$. 
By placing a GP prior on the function $\hat{f}(\mathbf{x}, s)$ directly, 
one can
sequentially optimize acquisition functions such as
\[ \max_{(\mathbf{x},s) \in \bm{\mathcal{X}}\times [0,1]}~\frac{\alpha(\mathbf{x}, s)}{c(s)}.\]

Our problem is equivalent to the multifidelity BO problem with $s$
replaced by $t\in[0,T]$ with the following key differences:
\begin{enumerate}
    \item In our case, $\hat{f}(\mathbf{x},T)$ can never be observed when $t< T$. More generally, any
      set $\{t_i \}, i=1, \ldots, n$ has a unique ordering $t_1 < \ldots < t_n$;
      and while we are at any $t_i$, evaluation of $f(\mathbf{x}, t \neq t_i)$ is not possible.
    \item We make only one observation $f(\mathbf{x}, t)$ for each $t\in[0,T]$.
    \item The cost of each evaluation $c(s) = c$ is a constant. However, similar to the general multifidelity framework, we
    are still interested in maximizing $f(\mathbf{x}, T)$.
\end{enumerate}
Another area of future work is defining a time-dependent cost
function that can be used within the multifidelity framework to also select the
time schedule of observations automatically; this approach  is of specific interest in
handling continuous-time systems.

