\section*{Classification of short baseline neutrino oscillation experiments}

The short baseline anomalies come from 4 different classes of experiments~\cite{sterile}.
Based on this, we can classify short baseline neutrino oscillation experiments into
following 5 groups;
\begin{enumerate}
\item test of LSND signal,
\item test of MiniBooNE signal,
\item test of reactor antineutrino anomaly,
\item test of Gallium anomaly, and
\item others.
\end{enumerate}
We discuss each of the above group in following sections.
The last group is all other experiments,
they are mainly experiments motivated by 1~eV sterile neutrinos
and not short baseline anomalies themselves.
Therefore many experiments in (v) are not short baseline oscillation experiments.

The short baseline anomalies are unsolved mysteries in this community,
and they attract many theorists and experimentalists.
The search of 1~eV sterile neutrino is one of the big branches of
the neutrino experiment community~\cite{neutrino}.
Therefore, it is rather impossible to cover all experiments in this note,
however, we try to cover most of experiments planned in the near future.

\section{LSND signal and experiments designed to test to it}

\subsection*{LSND experiment}

The origin of $\Delta m^2_{sterile}\sim 1~\mbox{eV}^2$ is the LSND experiment~\cite{LSND_osc},
where muon antineutrino (0 to 53~MeV) are produced by pion decay-at-rest (DAR),
and detected by a liquid scintillator detector at 31~m from the target.
The LSND experiment measured $\bar\nu_e$ candidate events by
utilizing the coincidence of the prompt Cherenkov radiation
from the positron and the delayed neutron capture by a hydrogen.
\begin{eqnarray}
\bar\nu_\mu\stackrel{oscillation}{\longrightarrow}\bar\nu_e+p\rightarrow\e^+(Cherenkov)+n(capture)~.\nonumber
\end{eqnarray}

The LSND experiment observed an excess of $\bar\nu_e$ candidate events.
This small ($<1\%$ oscillation probability)
but statistically significant signal is
consistent with the presence of sterile neutrinos ($\bar\nu_\mu\rightarrow\nu_{sterile}\rightarrow\bar\nu_e$).

Meantime, the KARMEN experiment~\cite{KARMEN} excluded high $\Delta m^2$ region,
and the Bugey experiment~\cite{Bugey} excluded all low $\Delta m^2$ region
of the LSND signal region in $\Delta m^2 - sin^2 2\theta$ plane.
The combined result suggests the LSND signal is most likely due to sterile neutrinos around
1~eV region.

\subsection*{Experiments to test LSND signal}

The LSND experiment has  limited statistics,
also, the duty cycle (nominal run, $\sim$25~Hz) was high with a wide pulse.
This allows LSND to accept large amount of cosmic backgrounds.
Also, the detector was located to the direction of the primary beamline,
and neutrinos from pion decay-in-flight (DIF) made additional backgrounds.
Therefore, to test LSND signal, experiments are desired to have;
\begin{itemize}
\item LSND beam energy, baseline, and the detector which can tag $\bar\nu_e$ candidate events,
\item higher statistics,
\item known and narrow beam structure, and
\item detector located on off-axis.
\end{itemize}

The promising experiments are off-axis liquid scintillator experiments
at various neutron spallation sources in the world.
The high pion-DAR $\bar\nu_\mu$ flux is available (and free!)
from Oak Ridge National Laboratory (ORNL)~\cite{OscSNS},
J-PARC Materials and Life Science Experimental Facility (MLF)~\cite{OscJSNS},
and European spallation source (ESS)~\cite{ESSnuSB}.
On top of the high neutrino flux, proton pulses hitting the target
to produce neutrons have well known beam structure.
Therefore these experiments cover the desired features to test the LSND signal.
Presently, OscSNS~\cite{OscSNS} is about to write a proposal,
and J-PARC group~\cite{OscJSNS} submitted proposal to J-PARC.

\section{MiniBooNE signal and experiments designed to test it}

\subsection*{MiniBooNE experiment}

The MiniBooNE experiment is designed to test the LSND signal within
the two massive neutrino oscillation hypothesis.
However, muon (anti)neutrinos are now made by
pion DIF at the Fermilab Booster neutrino beamline~\cite{MB_beam},
and the baseline is 541~m from the target (pion decay length is $\sim$18~m).
The MiniBooNE detector is a spherical mineral oil based Cherenkov detector~\cite{MB_detec},
and the $\nu_e$($\bar\nu_e$) candidate signals are
measured as single isolated electron-like Cherenkov ring.
In this way, the systematics of MiniBooNE
is completely different from the LSND experiment,
but MiniBooNE can test 1~eV sterile neutrino hypothesis,
because of similar L/E with LSND.
\begin{eqnarray}
\nu_\mu&\stackrel{oscillation}{\longrightarrow}&\nu_e+n\rightarrow p+e^+(Cherenkov)~,\nonumber\\
\bar\nu_\mu&\stackrel{oscillation}{\longrightarrow}&\bar\nu_e+p\rightarrow n+\e^-(Cherenkov)~.\nonumber
\end{eqnarray}

The Booster neutrino beamline can run either in neutrino mode or in antineutrino mode,
by focussing either positive or negative mesons.
Since LSND signal was interpreted $\bar\nu_\mu\rightarrow\nu_{sterile}\rightarrow\bar\nu_e$ oscillations,
running in antineutrino mode is more interesting.
However, the antineutrino mode run suffers from lower statistics
and higher backgrounds~\cite{MB_ANTICCQE}
(especially from the muon-neutrino contamination in the antineutrino mode beam~\cite{MB_WS}),
therefore the experiment started in neutrino mode prior
to the antineutrino mode running,
and in the meantime systematics (the neutrino flux, neutrino interactions,
and the detector response) were studied.

Unlike the LSND experiment, expected signal to noise is much lower.
There are 2 dominant backgrounds of $\nu_e$($\bar\nu_e$) candidate events.
The first one is the intrinsic $\nu_e$($\bar\nu_e$) contaminated in the beam.
The majority of them are made by muon decay, therefore,
MiniBooNE constrains them by simultaneously measuring
$\nu_\mu$ charged current quasi-elastic (CCQE) events~\cite{MB_CCQE,MB_CCQEPRD},
where measured $\nu_\mu$ is related to intrinsic $\nu_e$ through
the pion decay chain ($\pi^+\rightarrow\nu_\mu+\mu^+,~\mu^+->\bar\nu_\mu+e^++\bar\nu_e$) in their simulation.

The second largest background is the misID of neutral current (NC) events,
mainly NC $\pi^{\circ}$ production.
Although a $\pi^{\circ}$ decays to two gamma rays which should be distinguishable from
an electron (positron) Cherenkov ring,
sometimes decay kinematics make two gamma rays look like one gamma ray
(asymmetric decay, gamma rays are too close).
Then, the Cherenkov ring from one gamma ray is indistinguishable from an electron (positron).
For this, MiniBooNE internally measured NC$\pi^{\circ}$ production,
and the measured information was used to correct $\pi^{\circ}$ production rates
in the simulation~\cite{MB_pi0}.

After the 10 years running in both neutrino and antineutrino mode,
the MiniBooNE experiment observed excesses in both neutrino
and antineutrino mode runs~\cite{MB_osc}.
The final result is shown in Figure~\ref{fig:MB_osc}.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=0.5]{Fig_osc_01.pdf}
\includegraphics[scale=0.5]{Fig_osc_02.pdf}
\end{center}
\caption{\label{fig:MB_osc}
(color online)
The final MiniBooNE oscillation results~\cite{MB_osc}.
The top plots are antineutrino mode results,
and bottom plots are neutrino mode results.
The left plots show the reconstructed neutrino energy distribution
of oscillation candidate events,
and the right plots show the allowed region in $\Delta m^2$-$sin^22\theta$,
where the best fit points are shown in black stars.
Both modes show excesses in the low energy region,
while the neutrino mode has higher statistical significance.
On the other hand,
the compatibility with the LSND signal is better in antineutrino mode.
 }
\end{figure}

\subsection*{Experiments to test MiniBooNE signal}

The measured signal, especially in neutrino mode,
does not quite agree with the expected sterile neutrino signal.
The MiniBooNE detector cannot distinguish an electron (positron) and a gamma ray,
therefore $\nu_\mu$NC interaction
with single gamma ray in the final state is a potential misID background.
Therefore, to test the MiniBooNE signal, experiments are desired to have;
\begin{itemize}
\item MiniBooNE beam energy and baseline,
\item ability to distinguish NC or CC interaction, or
\item ability to distinguish an electron (positron) and a gamma ray.
\end{itemize}

The MiniBooNE+ was proposed to fulfil these criteria~\cite{MB+}.
By doping scintillator (PPO) in the MiniBooNE detector mineral oil,
MiniBooNE+ can measure scintillation light from the neutron capture.
This allows statistical separation between $\nu_e$CCQE interaction
(higher proton multiplicity in the final state),
and $\nu_\mu$NC interactions
(protons and neutrons are half-and-half in the final state).
However, the proposal of MiniBooNE+ was not accepted by Fermilab recently.

The MicroBooNE experiment~\cite{uB} is a new experiment on
the Fermilab Booster neutrino beamline to test the MiniBooNE signal.
It is also an important project for the future large liquid argon (LAr) TPC experiment,
such as LBNE~\cite{LBNE}.
The MicroBooNE LArTPC detector
has an ability to separate an electron from single gamma ray,
by utilizing vertex-shower separation and dE/dx before developing the shower.
This clearly tells if MiniBooNE excess is by an electron (positron) or a
gamma ray~\cite{Fleming_uB}.
The MicroBooNE experiment is under commissioning stage,
and they expect beam data at the end of 2014.

Although the T2K experiment~\cite{T2K_osc} is designed to measure
the neutrino Standard Model ($\nu$SM) parameters and it does not use the Booster neutrino beam,
J-PARC neutrino beam~\cite{T2K_beam} has a similar beam peak ($\sim$600~MeV)
as the Booster neutrino beam but is narrower,
and the baseline to the near detector complex~\cite{T2K_detec}
is close (280~m) to what MiniBooNE has.
Therefore, T2K is sensitive to the MiniBooNE signals~\footnote{
Current sterile neutrino search analysis in T2K is looking for
$\nu_e$ disappearance~\cite{T2K_sterile}, instead of $\nu_e$ appearance.}.
The magnetic field in the near detector is a great advantage.
It can allow the sign selection of the signal.
The NC background (mostly ambient gamma rays) can be understood
from the internal measurement~\cite{T2K_nue}.

\section{Reactor antineutrino anomaly and experiments designed to test it}

\subsection*{Reactor antineutrino anomaly}

The re-evaluation of reactor electron antineutrino flux calculation provides
consistently lower rate than world reactor data (about 6\%)~\cite{ReactorAnomaly1}.
This, so called {\it reactor antineutrino anomaly}
can be interpreted as neutrino oscillations with
1~eV sterile neutrinos~\cite{ReactorAnomaly2}.

\subsection*{Experiments to test Reactor anomaly}

Reactor anomaly can be tested by small scale experiments,
by measuring neutrino flux with small detectors with very short baseline ($\sim$15~m).
To detect low energy reactor antineutrino ($\sim$4~MeV),
detector needs to be sensitive to low energy events.
The common choice is the liquid scintillator detector,
where large mass can be prepared at a relatively low cost.

There are number of such experiments designed for R\&D of
neutrino reactor monitoring for
nuclear non-proliferation (SCRAAM~\cite{SCRAAM}, Nucifer~\cite{Nucifer}, etc).
These experiments are naturally served to test reactor antineutrino anomaly.
Due to affordable cost of the experiments,
several new experiments are also proposed to test the reactor antineutrino anomaly
(DANSS~\cite{DANSS}, PROSPECT~\cite{PROSPECT}, STEREO~\cite{STEREO}, etc).

\section{Gallium anomaly and experiments designed to test}

\subsection*{Gallium anomaly}

The two of pp-solar neutrino experiments, SAGE~\cite{SAGE} and GALLEX~\cite{GALLEX},
used highly radioactive sources to calibrate gallium detectors.
\begin{eqnarray}
\nu_e + ^{71}Ga \rightarrow ^{71}Ge + e^-~.
\end{eqnarray}
But some of these measurements using mega-curie $^{51}$Cr or $^{37}$Ar sources
showed lower event rates than expected,
and this so-called {\it Gallium anomaly}
can be understood by neutrino oscillations with 1~eV sterile neutrinos~\cite{Giunti_Gallium}.

\subsection*{Experiments to test Gallium anomaly}

To test “Gallium anomaly”, a highly radio-active neutrino source
and a very sensitive detector are required.
Existing high sensitivity solar or reactor neutrino detectors
(Borexino as ``SOX''~\cite{SOX}, KamLAND as ``Ce-LAND''~\cite{Ce-LAND},
SNO+~\cite{SNO+}, DayaBay~\cite{DayaBay+}, etc)
are good candidates for this purpose.
Similarly, SAGE group proposed to build a new detector but
reuse the liquid gallium from old detectors~\cite{SAGE+}.
Proposed solar neutrino experiment (LENS~\cite{Lens-sterile}, etc) and coherent scattering
experiment (RICOCHET~\cite{RICOCHET}, etc) can look for 1~eV sterile neutrinos, too.

\section{Others, 1~eV sterile neutrino searches}

\subsection*{Tests by existing facilities}

Once we assume 1~eV sterile neutrinos, some existing facilities are also sensitive to signals,
even though the experiments are not originally designed to test 1~eV sterile neutrinos.
The IceCube experiment~\cite{IceCube_PeV} is designed
to measure astrophysical ultra high energy neutrinos,
however, 1~eV sterile neutrinos cause disappearance signals
for $>$100~GeV high energy atmospheric neutrinos~\cite{IceCube_sterile}.
MINOS+ experiment~\cite{MINOS+} is an extension run of
MINOS experiment~\cite{MINOS_osc,MINOS_sterile} during NOvA beam configuration era
(medium energy NuMI, $\sim$7~GeV peak at Sudan mine)~\cite{NOvA}.
One of physics goal of MINOS+ is to look for
$\nu_\mu$ disappearance signal due to 1 eV sterile neutrinos.

\subsection*{Tests by R\&D facilities}

Many R\&D experiments for other purposes often look for 1~eV sterile neutrinos.
For example IsoDAR experiment~\cite{IsoDAR} look for sterile neutrino oscillation using
$^8$Li isotope made by the high power cyclotron.
This cyclotron is a part of the R\&D for the DAE$\delta$ALUS experiment~\cite{DAEdALUS}.
The $\nu$STORM~\cite{nuSTORM} experiment uses the muon storage ring,
which is an important step for the future neutrino factory~\cite{VLENF}.
KDAR experiment~\cite{KDAR} uses mono-energetic kaon DAR muon neutrinos (236~MeV),
and the detector requires high resolution, such as LArTPC.
This can be considered a part of LArTPC technology R\&D,
and in fact, all LArTPC sterile neutrino searches, such as MicroBooNE~\cite{uB},
LAr1-ND~\cite{LAr1-ND}, and NESSiE~\cite{NESSiE},
have detector R\&D aspects for future large LArTPC experiments.

\subsection*{Ultimate 1~eV sterile neutrino search experiments}

Experiments including precise measurement of oscillation probability with
function of L/E can be considered in this group.
LSND reloaded~\cite{LSND-reloaded} was proposed to test short baseline neutrino oscillations by
measuring oscillations with function of L/E in a large detector,
such as gadolinium doped Super-Kamiokande.
Similar concept may be applied to reactor antineutrino anomaly experiments
and gallium anomaly experiments,
where neutrino sources are small and low energy.
In those experiments, precise L/E dependence measurement may be possible by
either using a large detector or multiple small detectors,
or moving sources and/or detectors.
The precise oscillation measured in this way
is a strong evidence of sterile neutrinos
and it is a missing part from past experiments.

\section*{Acknowledgement}

The author thanks Ranjan Dharmapalan for the careful reading of this manuscript.
The author thanks to the organisers of ``NuPhys2013, prospects in neutrino physics,
(Institute of Physics, London, UK)'' for the invitation to the conference.

\section*{References}


\section{Introduction}

Let $\ell$ be a prime and let $k$ be field of characteristic zero.
Let $a_1,\dots ,a_n\in k^\ast$ with symbol $\alpha:=\{a_1,\dots ,a_n\}\in K^M_n(k)$.
Associated to this there is a standard norm variety $X_{\alpha,\ell}$ (see \cite{rost02,SJ}) which played a crucial role in Voevodsky's proof of the Bloch--Kato conjecture \cite{Voe}.

 Asok proved that $X_{\alpha,\ell}$  is geometrically rationally connected, see \cite[Proposition 2.6 and Remark 2.9]{asok}.
In \cite[Theorem 1]{BHS}, the authors improve Asok's result by showing that  $X_{\alpha,\ell} $ is geometrically $\mathbb A^1$-connected.
Equivalently, $X_{\alpha,\ell} $  becomes universally $R$-trivial over the algebraic closure $\bar k$, which means that after any field extension $L/\bar k$, any two $L$-rational points of $X_{\alpha,\ell}\times_kL$ can be joined by a chain of rational curves.
(The result of Asok handled the case where $L$ is algebraically closed.)

Recall that a variety $X$ is retract rational if the identity on $X$ factors as rational map through some projective space.
Any retract rational variety is unirational, while the converse is not true in general.
Moreover,
Asok and Morel showed that retract rational varieties are $\mathbb A^1$-connected (at least in characteristic zero), see \cite[Theorem 2.3.6]{asok-morel}.
The converse to this is an open question, see \cite[Remark 2.3.11]{asok-morel}.
It is thus natural to ask whether standard norm varieties are not only geometrically $\mathbb A^1$-connected but in fact geometrically retract rational.
This question goes back to Asok \cite[Remark 2.9]{asok} and  Balwe--Hogadi--Sawant \cite[Introduction]{BHS}.
In this short note we answer this question positively.

\begin{theorem} \label{thm:main}
Let $k$ be a field of characteristic zero and let $\alpha:=\{a_1,\dots ,a_n\}\in K^M_n(k)$ be a symbol for some $n\geq 2$.
Then for any prime $\ell$, the standard norm variety $X_{\alpha,\ell}$ is geometrically retract rational.
\end{theorem}

The above theorem is formulated for fields of characteristic zero, because the standard norm varieties under consideration are mostly studied under this assumption, see e.g.\ \cite{SJ}.
However,  our proof does not make use of resolution of singularities and works in arbitrary characteristic.
In particular, standard norm varieties in positive characteristic, defined analogously to the case of characteristic zero,  are geometrically retract rational, despite the fact that it is not clear that these varieties have a smooth projective model.

We note that  $X_{\alpha,\ell}$ is in general not retract rational over the ground field $k$.
In fact, whenever $\alpha$ is nonzero in $K^M_n(k)/\ell$, the norm variety $X_{\alpha,\ell}$ is never retract rational (it is not even unirational, because it does not contain a rational point).

\subsection{Notation} \label{subsec:notation}
An algebraic scheme is a separated scheme of finite type over a field.
A $k$-variety is an integral algebraic scheme over a field $k$.
For a $k$-variety $X$ and a positive integer $m$, we denote by $\operatorname{Sym} ^mX=X^m/S_m$ the $m$-th symmetric power of $X$.
In general, $\operatorname{Sym}^mX$ is an algebraic scheme over $k$; it is integral, hence a variety in the above sense, if $X$ is geometrically integral.
We will use that $\operatorname{Sym}^m$ is a covariant functor on the category of geometrically integral $k$-varieties.

A $k$-variety $X$ is retract rational if there is a rational map $f:X\dashrightarrow \mathbb P^N$ for some $N\geq 1$ and a rational map $g:\mathbb P^N\dashrightarrow X$ such that $g\circ f$ is defined and agrees with the identity as a rational map.
In other words, there is a rational map $f:X\dashrightarrow \mathbb P^N$ which admits a rational section.
Clearly, if $X$ is retract rational, then so is any birational model (e.g.\ any non-empty open subset).
Moreover, $X$ is retract rational if and only if $X\times \mathbb P^n$ is retract rational for some $n\geq 0$.
In particular, any stably rational variety is retract rational.

By definition,  $X$ is retract rational if and only if  there are open dense subsets $U\subset X$ and $V\subset \mathbb P^N$ for some $N$ such that the identity on $U$ factors through $V$: $U\to V\to U$.
It  follows from this description that  retract rational varieties are geometrically integral (it is in the above notation enough to check that $U$ is geometrically integral, which follows from the fact that $V$ is geometrically integral).
In particular, if $X$ is a retract rational $k$-variety then any symmetric power $\operatorname{Sym} ^m X$ is a k-variety in the above sense.

\section{Three lemmas}

We say that a rational map $f:X\dashrightarrow Y$ between two $k$-varieties is generically injective if there are dense open subsets $U\subset X$ and $V\subset Y$ such that $f$ restricts to a closed embedding $U\hookrightarrow V$.
Equivalently, if $\eta_X\in X$ denotes the generic point, then $f$ induces an isomorphism between residue fields $\kappa(f(\eta_X))\stackrel{\cong}\to \kappa(\eta_X)$.

\begin{lemma} \label{lem:0}
Let $k$ be a field and let $Y$ be a retract rational $k$-variety.
Then for any rational map $f:Y\dashrightarrow \mathbb P^N$ which is generically injective, there is a rational section $g:\mathbb P^N\dashrightarrow Y$, i.e.\ a rational map $g$ such that $g\circ f$ is defined and coincides with the identity as a rational map.
\end{lemma}
\begin{proof}
By Saltman's criterion \cite[Theorem 3.4]{saltman} (see also \cite[Proposition 1.2]{CTS}), there is a dense open subset $V\subset Y$ such that for any local $k$-algebra $A$ with residue field $\kappa$, the map $V(A)\to V(\kappa)$ is surjective.
(In fact, this is the easy direction of Saltman's criterion: choose $V$ such that there is a dense open subset $U\subset \mathbb A^N$ such that the identity on $V$ factors via $V\to U\to V$.
Any map $\epsilon:\operatorname{Spec} \kappa\to \mathbb A^N$ lifts to a map $\epsilon_A:\operatorname{Spec} A\to \mathbb A^N$.
If $\operatorname{im}(\epsilon)\subset U$, then $\operatorname{im} (\epsilon_A)\subset U$, because $A$ is local.
Hence,   $U(A)\to U(\kappa)$ is surjective.
By the factorization $V\to U\to V$ of $\operatorname{id}_V$, $V(A)\to V(\kappa)$ is surjective as well.)

We apply this to the local ring $A:=\mathcal O_{\mathbb P^N,f(\eta_Y)}$ of $\mathbb P^N$ at the image of the generic point of $Y$.
Since $f$ is generically injective, the residue field of $A$ is isomorphic to $k(Y)$.
It follows that the diagonal point of $Y_{k(Y)}$ lifts to an $A$-valued point of $Y$.
In other words, there is a morphism $\operatorname{Spec} A\to Y$ which restricts to the inclusion of the generic point $\eta_Y$ on the closed point of $\operatorname{Spec} A$.
Since $A=\mathcal O_{\mathbb P^N,f(\eta_Y)}$, we conclude that there is a rational map $g:\mathbb P^N\dashrightarrow Y$ which is defined at $f(\eta_Y)$ and such that $g\circ f$ coincides as rational map with the identity on $Y$.
This concludes the proof of the lemma.
\end{proof}

\begin{remark}
Since any variety is birational to a hypersurface, it follows from Lemma \ref{lem:0} that any retract rational k-variety $Y$ admits a rational map $Y\dashrightarrow \mathbb P^{\dim Y+1}$ with a rational section.
\end{remark}

\begin{remark}
 The argument in Lemma \ref{lem:0}
 works verbatim if we replace $\mathbb P^N$ by any $k$-variety $Z$.
Hence, if $Y$ is retract rational, then any generically injective map $f:Y\dashrightarrow Z$ to a $k$-variety $Z$ admits a rational section.
 (Thanks to Remy van Dobben de Bruyn for pointing this out.)
In fact, this a priori stronger statement is   a formal consequence of Lemma \ref{lem:0} because any $k$-variety $Z$ admits a generically injective rational map $\iota:Z\dashrightarrow \mathbb P^N$ and so we may apply Lemma \ref{lem:0} to $\iota\circ f$ to produce a rational section of $\iota\circ f$ which restricts to a rational section of $f$, as we want.
 \end{remark}

\begin{lemma} \label{lem:1}
Let $Y$ be a retract rational $k$-variety over an infinite field $k$.
Then $\operatorname{Sym}^mY$ is retract rational over $k$ for any $m\geq 1$.
\end{lemma}
\begin{proof}
The statement is trivial if $\dim Y=0$, because retract rational varieties are geometrically integral and so $Y=\operatorname{Spec} k$ in this case.
We may therefore assume that $Y$ has positive dimension.
Up to replacing $Y$ by a projective model, we may assume that there is a closed embedding $f:Y\hookrightarrow  \mathbb P^{N}$ for some $N$.
Since $\dim Y\geq 1$,  composing the embedding $Y\hookrightarrow \mathbb P^N$ with some Veronese embedding of $\mathbb P^N$, we can ensure that
\begin{align}\label{eq:lem:1}
\dim   \left\langle f(Y)\right\rangle \geq m,
\end{align}
where $\left\langle f(Y)\right\rangle \subset \mathbb P^N$ denotes the smallest linear subspace that contains $f(Y)$.

Since $Y$ is retract rational,  Lemma \ref{lem:0} shows that $f$ admits a rational section, i.e.\ a rational map $g:\mathbb P^N\dashrightarrow Y$ such that $g\circ f$ is defined and coincides with the identity (as a rational map).
Taking the symmetric product, we get maps
$$
\operatorname{Sym}^m f:\operatorname{Sym}^m Y\longrightarrow \operatorname{Sym}^m\mathbb P^N\ \ \text{and}\ \ \operatorname{Sym}^mg:\operatorname{Sym}^m \mathbb P^N\dashrightarrow \operatorname{Sym}^m Y
$$
such that $\operatorname{Sym}^m g\circ \operatorname{Sym}^m f$ is defined and coincides with the identity.

By \cite{mattuck}, there is a birational map
$$
\varphi: \operatorname{Sym}^m\mathbb P^N\stackrel{\cong}\dashrightarrow \mathbb P^{m N} .
$$
Hence there are open dense subsets $U\subset \operatorname{Sym}^m\mathbb P^N$ and $V\subset \mathbb P^{mN}$ such that $\varphi$ induces an isomorphism $U\cong V$.

From now on we use that $k$ is infinite, which implies that the $k$-rational points of $U$ are Zariski dense.
We may therefore pick a general $k$-rational point $\xi\in U$ and note that $\xi$ corresponds to a set $\{x_1,\dots ,x_m\}$ of general $k$-rational points $x_i\in \mathbb P^N$.
By (\ref{eq:lem:1}), $\dim \left\langle f(Y)\right\rangle \geq m$ and so $N\geq m$ (because $\left\langle f(Y)\right\rangle $ is a subspace of $\mathbb P^N$).
It follows that the general points $x_1,\dots ,x_m$ of $\mathbb P^N$ are linearly independent in $\mathbb P^N$.
Similarly, a general $k$-rational point $\zeta\in \operatorname{Sym}^mY$ corresponds to a set $\{y_1,\dots ,y_m\}$ of general $k$-rational points $y_i\in Y$.
(Here we use that the set of $k$-rational points of $Y$ is Zariski dense because $k$ is infinite and $Y$ is retract rational, hence unirational over $k$.)
Since $\zeta$ is general, $f$ is defined at each of the points $y_1,\dots ,y_m$ and since $\dim \left\langle f(Y)\right\rangle \geq m$, the points $f(y_1),\dots ,f(y_m)$ of $\mathbb P^N$ are linearly independent.
It follows that there is a linear automorphism $h:\mathbb P^N\stackrel{\cong}\to \mathbb P^N$ defined over $k$ such that $h(f(y_i))=x_i$ for all $i$.
The induced automorphism
$$
\operatorname{Sym}^m h:\operatorname{Sym}^m\mathbb P^N\longrightarrow \operatorname{Sym}^m\mathbb P^N
$$
satisfies
$$
\operatorname{Sym}^m h(\operatorname{Sym}^m f(\zeta))=\xi\in U.
$$
It follows that the composition
$$
\operatorname{Sym}^m Y\stackrel{\operatorname{Sym}^m f}\dashrightarrow \operatorname{Sym}^m\mathbb P^N \stackrel{\operatorname{Sym}^m h} \longrightarrow \operatorname{Sym}^m\mathbb P^N \stackrel{\varphi}\dashrightarrow \mathbb P^{m N}
$$
is defined and admits a rational section, given by
$$
\mathbb P^{m N}\stackrel{\varphi^{-1}}\dashrightarrow \operatorname{Sym}^m\mathbb P^N  \stackrel{\operatorname{Sym}^m h^{-1}} \longrightarrow \operatorname{Sym}^m\mathbb P^N\stackrel{\operatorname{Sym}^m g}\dashrightarrow \operatorname{Sym}^m Y .
$$
(Here we use that $\operatorname{Sym}^m g$ is defined at $\operatorname{Sym}^m f(\zeta)$, because $\zeta\in \operatorname{Sym}^mY$ was chosen to be general and $g\circ f$ is defined as a rational map by assumption.)
This concludes the proof of the lemma.
\end{proof}

\begin{lemma} \label{lem:2}
Let $h:X\to Y$ be a dominant morphism of $k$-varieties whose generic fibre $X_\eta$ is irreducible.
If $X_\eta$ is retract rational over $k(Y)$ and $Y$ is retract rational over $k$, then $X$ is retract rational over $k$.
\end{lemma}
\begin{proof}
Since the generic fibre $X_\eta$ of $h$ is retract rational, there is a rational map
$$
f:X_\eta\dashrightarrow \mathbb P^r_{k(Y)}
$$
which admits a rational section $g:\mathbb P^r_{k(Y)}\dashrightarrow X_\eta$.
The maps $f$ and $g$ spread out to rational maps of $k$-varieties
$$
F:X\dashrightarrow \mathbb P^r\times Y\ \ \text{and}\ \ G:\mathbb P^r\times Y \dashrightarrow X ,
$$
such that $G\circ F$ is defined and agrees with the identity as a rational map.
Moreover, $F$ and $G$ are rational maps over $Y$, i.e.\ they are compatible with the corresponding maps to $Y$.
Since $Y$ is retract rational as a $k$-variety, there is a rational map $f':Y\dashrightarrow \mathbb P^N$ with a rational section $g':\mathbb P^N\dashrightarrow Y$.
We identify $ \mathbb P^r\times \mathbb P^N$ birationally with $\mathbb P^{r+N}$, by noting that both varieties are compactifications of $\mathbb A^{r+N}$.
The compositions
$$
(\operatorname{id}\times f')\circ F:X\dashrightarrow \mathbb P^{r+N}
$$
and
$$
G\circ (\operatorname{id}\times g'):\mathbb P^{r+N}\dashrightarrow X
$$
are then rational maps of $k$-varieties such that the composition
$$
G\circ (\operatorname{id}\times g')\circ (\operatorname{id}\times f')\circ F:X\dashrightarrow X
$$
is defined (this uses that $F$ and $G$ are rational maps over $Y$, hence defined at some point that dominates the generic point of $Y$) and coincides with the identity.
It follows that $X$ is retract rational, as claimed.
\end{proof}

\section{Proof of Theorem \ref{thm:main}}

We may replace $k$ by $\bar k$ and assume that $k$ is algebraically closed.
For $n=2$, $X_{\alpha,\ell}$ is a Severi--Brauer variety and so it is rational over the algebraically closed field $k$.
We may thus assume $n\geq 3$ and argue by induction.
By construction (see e.g.\ \cite[\S 2]{SJ} or \cite[\S 4]{BHS}), there is a birational model $X$ of $X_{\alpha,\ell}$ with  a dominant morphism $f:X\to \operatorname{Sym}^\ell Y$ such that
\begin{itemize}
\item $Y$ is birational to $X_{\alpha',\ell}$ with  $\alpha':=\{a_1,\dots ,a_{n-1}\}\in K^M_n(k)$;
\item the generic fibre of $f$ is birational to the norm hypersurface $\{N-a_n=0\} $ over $k(\operatorname{Sym}^\ell Y)$, where
$$
N:R_{L/K}\mathbb G_{m,L}\longrightarrow \mathbb G_{m, K}
$$
is the map of tori associated to the norm map of fields $N:L\to K$, where $L:=k(Y\times \operatorname{Sym}^{\ell-1}Y)$ and $K:=k(\operatorname{Sym}^\ell Y)$, and  $ R_{L/K}\mathbb G_{m,L}$ denotes the Weil restriction of $\mathbb G_{m,L}$ to $K$.
\end{itemize}
Here the models $X$ and $Y$ are not necessarily proper nor smooth.
(Of course we could always assume  that one of the two conditions is satisfied, but the argument does not need that.)
We could assume that the generic fibre of $f$ is actually isomorphic to the norm hypersurface  $\{N-a_n=0\} $ over $k(\operatorname{Sym}^\ell Y)$, but again this will not play any role in the argument.

Since $k=\bar k$, there is an element $\xi\in k$ with $\xi^\ell=a_n$.
We may regard $\xi$ as an element of $L=k(Y\times \operatorname{Sym}^{\ell-1}Y)$ and we note that $N(\xi)=a_n$.
Multiplication with $\xi$ induces an automorphism of $R_{L/K}\mathbb G_{m,L}$ which identifies the norm 1 torus $\{N-1=0\}$ with   $\{N-a_n=0\}$.
Hence, since $\{N-1=0\}$  is retract rational by \cite[Theorem 4.1]{endo},  $\{N-a_n=0\}$ is retract rational as well.
Theorem \ref{thm:main} therefore follows from Lemmas \ref{lem:1} and \ref{lem:2} above, where we use that $k$ is algebraically closed, hence an infinite field.

\begin{remark}
Stably rational varieties are retract rational and the converse is known to fail  in general over non-closed fields.
However, it is not known whether the two notions coincide over algebraically closed fields, cf.\  \cite{PS}.
While Lemma \ref{lem:2} remains true if we replace retract rationality by stable rationality, this does not seem to be
clear for
 Lemma \ref{lem:1}.
More importantly,  the norm 1 tori over $k(\operatorname{Sym}^\ell Y)$ that appear in the above argument are retract rational for all primes,  but they are not stably rational for primes $\ell\geq 5$,  see \cite[Theorem 4.1]{endo}.
In particular,   the norm varieties $X_{\alpha,\ell}$ for $\ell\geq 5$ yield by Theorem \ref{thm:main} natural candidates for varieties over algebraically closed fields that are retract rational but which may potentially be not stably rational in general.
\end{remark}

\section*{Acknowledgements}
Thanks to Remy van Dobben de Bruyn for comments and to Anand Sawant for explaining the results in \cite{BHS} to me.
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 948066 (ERC-StG RationAlgic).


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

\section{Introduction}
Preferential Bayesian optimization (BO) has been an attractive approach for solving problems where reliable feedback is limited to \emph{pairwise comparison}, the so-called \emph{duels}.
This preference setting often appears in \emph{human-in-the-loop optimization} problems such as visual design optimization~\citep{koyama2020sequential} and generative melody composition~\citep{zhou2021interactive} because it is easier for humans to judge which one is better than to give an absolute rating~\citep{kahneman2013prospect}.
The system (i.e., the optimization method) in the human-in-the-loop optimization presents choices and receives preferential feedback \emph{interactively}.
To reduce the waiting time for users, the system is required to \emph{quickly} present the new options to users by learning from the observed feedback information.

An important challenge in preferential BO, which uses the preferential Gaussian process (GP) model to represent preference structure, is that the posterior distribution is computationally intractable skew GP.
The existing approaches to this difficulty are twofold:
The first approach is the Gaussian approximation, such as Laplace approximation (LA) and expectation propagation (EP)~\citep{chu2005-preference,chu2005-extensions}.
Actually, most of the preferential BO algorithms proposed so far are based on this Gaussian approximation~\citep{brochu2010-tutorial,nielsen2015-perception,gonzalez2017preferential,siivola2021-preferential,fauvel2021-efficient}, leading to good computational efficiency.
However, the predictive accuracy can be poor, as discussed in Section~\ref{sec:comparison_PGP}, because the Gaussian approximation \emph{ignores} the skewness of the exact posterior represented as the skew GP~\citep{benavoli2020skew,benavoli2021-preferential,benavoli2021-unified}.
The second approach is to directly employ the skew GP model using the Markov chain Monte Carlo (MCMC) method~\citep{benavoli2020skew,benavoli2021-preferential,benavoli2021-unified}.
\citet{benavoli2020skew,benavoli2021-preferential,benavoli2021-unified} argued that the Gaussian approximation has poor estimation accuracy, and their MCMC-based estimation is effective in preferential BO.
However, the experimental evidence was insufficient to show why the Gaussian approximation degrade the preferential BO.

The aim of this study is to build truly \emph{practical} preferential BO with experimental evidence to support it.
To this end, we first evaluate the accuracy of the Gaussian approximation-based methods, which have been most used in the preferential BO literature.
Our experiments reveal that those methods fail to correctly predict the probability of duel outcomes, which is critical in the preferential BO setting.
This result implies that it is important to consider the skewness of the true posterior, and then we also improve MCMC for skew GP estimation.
Specifically, we present MC estimators that reduce the estimation variance for skew GP.
Furthermore, we empirically clarify the high computational efficiency of the Gibbs sampling compared with an MCMC method employed in \citep{benavoli2020skew,benavoli2021-unified}.
These improvements make the preferential BO algorithm with skew GP more attractive.
On the other hand, the preferential BO algorithms using MCMC can still be unrealistic in preferential BO, where the computational time may become a bottleneck in many cases.

To build practical preferential BO, we further develop a simple and computationally efficient preferential BO method that incorporates the skewness of the true posterior in a randomized manner.
The basic idea is to use a posterior additionally conditioned by a random sample from the original posterior itself, called \emph{hallucination}, by which we show that computing acquisition functions (AFs) in the proposed method is computationally efficient while reflecting the skewness of the original posterior.
Noteworthy, any powerful AFs developed so far in the standard BO literature (e.g., upper confidence bound (UCB)~\citep{Srinivas2010-Gaussian} and expected improvement (EI)~\citep{Mockus1978-Application}) can be integrated, which improves the flexibility of the proposed method.

Our contributions are summarized as follows:
\begin{enumerate}
    \item \textbfit{Rethinking existing preferential BO algorithms}: We reveal the poor accuracy of Gaussian approximation-based methods in the true posterior approximation, which implies the importance of considering \emph{skewness} in the preferential BO setting. Motivated by this observation, we improve MCMC for skew GP estimation.

    \item \textbfit{Building practical preferential BO algorithm}: We develop a simple and computationally efficient method, called \emph{hallucination believer}, that can reflect the skewness of the original posterior and integrate any powerful AFs developed so far in the standard BO literature.

    \item \textbfit{Extensive empirical validations}: Numerical experiments on 12 benchmark functions show that the proposed method achieves better or at least competitive performance in both terms of \emph{computational efficiency} and \emph{sample complexity} compared with Gaussian approximation-based preferential BO~\citep{gonzalez2017preferential,siivola2021-preferential,fauvel2021-efficient} and MCMC-based preferential BO~\citep{benavoli2021-preferential,benavoli2021-unified}, respectively.
\end{enumerate}

\section{Background}
\subsection{Preferential Bayesian Optimization Problem}

We consider that the preferential relation is modeled by a latent function  $f: {\cal X} \rightarrow \mathbb{R}$, where ${\cal X} \in \mathbb{R}^{d}$ is the input domain.
Our goal is to maximize the latent function as
\begin{align*}
    \*x_* = \mathop{\text{argmax}}_{\*x \in {\cal X}} f(\*x),
\end{align*}
through the dueling feedback, $\*x \succ \*x^\prime$, which implies $\*x$ is preferable to $\*x^\prime$.

\subsection{Gaussian Processes for Preferential Learning}
\label{sec:skewGP}
We assume that $f$ is a sample path of GP with zero mean and some stationary kernel $k: {\cal X} \times {\cal X} \rightarrow \mathbb{R}$.
Following \cite{chu2005-preference,chu2005-extensions}, we consider that the duel is determined as follows:
\begin{align*}
    \*x_{i, w} \succ \*x_{i, l}

    \Leftrightarrow
    f(\*x_{i, w}) + \epsilon_w > f(\*x_{i, l}) + \epsilon_l,
\end{align*}
where i.i.d. additive noise $\epsilon_w$ and $\epsilon_l$ follow the normal distribution ${\cal N}(0, \sigma^2_{\rm noise})$.
This is equivalent to assuming that the preferences are obtained by which a direct observation $y = f(\*x) + \epsilon$ is bigger or smaller, where $\epsilon \sim {\cal N}(0, \sigma^2_{\rm noise})$.
Therefore, the training data ${\cal D}_t$ can be written as,
\begin{align*}
    {\cal D}_t = \{ \*x_{i, w} \succ \*x_{i, l} \}_{i=1}^t \equiv \{v_i < 0\}_{i=1}^t,
\end{align*}
where $\*x_{i, w}$ and $\*x_{i, l}$ are a winner and a loser of $i$-th duel, respectively, and $v_i \coloneqq f(\*x_{i, l}) + \epsilon_l - f(\*x_{i, w}) - \epsilon_w$.
For brevity, we denote $\{v_i < 0\}_{i=1}^t$ as $\*v_t < \*0$, where $\*v_t \coloneqq (v_1, \dots, v_t)^\top$.

The exact posterior distribution $p(f \mid \*v_t < \*0)$ is skew GP, as shown in \citep{benavoli2021-preferential,benavoli2021-unified}.
That is, for any output vector $\*f_{\rm tes} \coloneqq \bigl( f(\*x_{1, {\rm tes}}), \dots, f(\*x_{m, {\rm tes}}) \bigr)^\top$, where $m \in \mathbb{N}$, we obtain the posterior distribution from Bayes' theorem as follows:
\begin{align}
    p\left(\*f_{\rm tes} \mid \*v_t < \*0 \right)
    = \frac{\Pr\left(\*v_t < \*0 \mid \*f_{\rm tes}\right)  p\left(\*f_{\rm tes}\right)}{\Pr(\*v_t < \*0)},
    \label{eq:posterior}
\end{align}
which is the distribution called multivariate unified skew normal \citep{azzalini2013-skewnormal,benavoli2021-preferential}.

Then, we briefly derive the components of Eq.~\eqref{eq:posterior}.
Hereafter, we denote $(i,j)$-th element of matrix as $[\cdot]_{i,j}$ and both $i$-th row of matrix and $i$-th element of vector as $[\cdot]_i$.
Suppoce that $\*X \in \mathbb{R}^{ (m + 2t) \times d}$ is
\begin{align*}
    \*X \coloneqq \! \bigl( \*x_{1, {\rm tes}}, \dots, \*x_{m, {\rm tes}}, \*x_{1, w}, \dots, \*x_{t, w}, \*x_{1, l}, \dots, \*x_{t, l} \bigr)^\top \!.
\end{align*}
Then, the prior for $\*f \coloneqq \bigl(f([\*X]_0), \dots, f([\*X]_{m + 2t}) \bigr)^\top \in \mathbb{R}^{m + 2t}$ is,
\begin{align*}
    \*f \sim {\cal N}(\*0, \*K),
\end{align*}
where $[\*K]_{i, j}$ is $k([\*X]_i, [\*X]_j)$.
Since $\*f_{\rm tes}$ and $\*v_t$ are linear combinations of $\*f$ and the noises $\epsilon_{w}$ and $\epsilon_{l}$, we see that,
\begin{align}
    \begin{bmatrix}
        \*f_{\rm tes} \\
        \*v_t
    \end{bmatrix}
    \sim {\cal N}(\*0, \*\Sigma),
    \label{eq:prior_tes_v}
\end{align}
where
\begin{align}
    \*\Sigma &\coloneqq \*A \bigl(\*K + \*B \bigr) \*A^\top \in \mathbb{R}^{ (m + t) \times (m + t) },
    \label{eq:prior_covariance} \\
    \*A &\coloneqq
    \begin{bmatrix}
        \*I_{m} & \*0 & \*0 \\
        \*0 & - \*I_{t} & \*I_{t}
    \end{bmatrix} \in \mathbb{R}^{ (m + t) \times (m + 2t) }, \\
    \*B &\coloneqq
    \begin{bmatrix}
        \*0 & \*0 \\
        \*0 & \sigma^2_{\rm noise} \*I_{2t}
    \end{bmatrix} \in \mathbb{R}^{(m + 2t) \times (m + 2t)},
\end{align}
and $\*I_{i} \in \mathbb{R}^{i \times i}$ is the identity matrix.
Thus, since the conditional distribution of multivariate normal (MVN) is MVN again \citep{Rasmussen2005-Gaussian}, $\*v_t \mid \*f_{\rm tes}$ follows MVN.
Consequently, $p\left(\*f_{\rm tes}\right)$ is a probability density function (PDF) of MVN and $\Pr(\*v_t < \*0)$ and $\Pr\left(\*v_t < \*0 \mid \*f_{\rm tes}\right)$ are cumulative distribution functions (CDF) of MVN.

Statistics for skew GP, such as expectation, variance, probability, and quantile, are hard to compute.
Specifically, they all need many times of computation of CDF of MVN \citep{azzalini2013-skewnormal}.
CDF of MVN \citep{Genz1992-numerical, Genton2017-Hierarchical} is computationally intractable and expensive.
Therefore, the posterior inference using CDF of MVN at each iteration for preferential BO can be unrealistic.

To approximate the statistics, \citet{benavoli2021-preferential,benavoli2021-unified} employed the MCMC, for which they showed that a sample path from skew GP $p(f | \*v_t < \*0)$ could be generated through the sampling from truncated MVN (TMVN) and MVN.
For example, they approximate the posterior mean $\mathbb{E}_{f | \*v_t < \*0} \bigl[ [\*f_{\rm tes}]_i \bigr]$ as,
\begin{align}
    \mathbb{E}_{f | \*v_t < \*0} \bigl[ [\*f_{\rm tes}]_i \bigr]
    \approx \frac{1}{M} \sum_{j=1}^M f_j (\*x_{i, {\rm tes}}),
    \label{eq:samplepath_estimation}
\end{align}
where $M$ is the number of MC samples and $f_j$ is a sample path from the posterior $f | \*v_t < \*0$.
Another popular approximation is Gaussian approximation using LA and EP \citep{chu2005-preference,chu2005-extensions}.
LA \citep{Mackay1996-bayesian} approximates the posterior by Gaussian distribution so that the mean is equal to the mode of the true posterior.
EP \citep{Minka2001-Minka} is an iterative moment-matching method, which aims to adjust the first and second moments.

\citet{benavoli2020skew} revealed that the true posterior of GP classification (GPC) with probit likelihood is also a skew GP.
Therefore, in principle, the approximate inferences for GPC can apply to preferential GP.
Actually, the survey of GPC \citep{Nickisch2008-approximations,Kuss2005-assessing} argues that EP exhibits excellent prediction performance \emph{in the GPC setting}.
However, the validity of using EP in the preferential BO has not been verified in detail so far.
Then, we will demonstrate in this study the predictive probability of the duel using EP can be inaccurate in Section~\ref{sec:comparison_PGP}.

\section{Rethinking Estimation of Preferential GP}
\label{sec:preferential GP_learning}
First, we investigate the limitation of the Gaussian approximation, which is the inaccurate predictive probability.
Second, we tackle the improvement of MCMC-based estimation for skew GP.
We show the practical effectiveness of Gibbs sampling for skew GP.
Furthermore, we newly derive MC estimators, showing that the estimation variance is lower than the estimators used in \citep{benavoli2021-preferential,benavoli2021-unified}.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\linewidth]{manuscripts/fig/preferentialGP_legend.pdf}
    \hfill
    \includegraphics[width=0.6\linewidth]{manuscripts/fig/preferentialGP_MCMC_acc.pdf}
    \hfill
    \includegraphics[width=0.6\linewidth]{manuscripts/fig/preferentialGP_LP.pdf}
    \hfill
    \includegraphics[width=0.6\linewidth]{manuscripts/fig/preferentialGP_EP.pdf}
    \caption{
        An illustrative example of preferential GP models fitted to 50 uniformly random duels, including $\*x_w \succ \*x_l$, from a sample path from GP.
        %
        We use the RBF kernel, which is used to generate the sample path, and set $\sigma^2_{\rm noise} = 10^{-4}$.
        %
        The top, middle, and bottom plots show the predictions of MCMC, LA, and EP, respectively.
        %
        The left plot shows the mean and 95\% credible interval of $f | \*v_t < \*0$.
        %
        The right plots show the estimator of $\Pr_t\bigl( f(\*x_1) \leq f(\*x_2) \bigr)$ and PDF $p(f(\*x_1), f(\*x_2) | \*v_t < \*0)$, for which an approximation by Eq.~\eqref{eq:posterior} using CDF of MVN is used in the top right plot.
        }
        \label{fig:examples_preferential GP}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_LA_mean.pdf}
    \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_LA_mode.pdf}
    \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_LA_prob.pdf}
    \subfigure[Truth vs. prediction plots for LA and EP.]{
        \centering
        \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_EP_mean.pdf}
        \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_EP_mode.pdf}
        \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_EP_prob.pdf}
    }

    \subfigure[RMSE of proposed MC estimator and full MC estimator \citep{benavoli2021-preferential} against ground truth.]{
        \includegraphics[width=0.2\linewidth]{manuscripts/fig/Ackley_mean_error.pdf}
        \includegraphics[width=0.19\linewidth]{manuscripts/fig/Ackley_mode_error.pdf}
        \includegraphics[width=0.19\linewidth]{manuscripts/fig/Ackley_prob_error.pdf}
    }

    \caption{
        Experimental results in the Ackley function ($d=4$).
        %
        We evaluate the estimators of the mean, mode, and duel probabilities $\mathrm{Pr}_t\bigl( f(\*x) \succ f(\*x^\prime) \bigr)$.
        %
        Other details are shown in Appendix~\ref{app:exp_settings}.
        %
        A similar tendency was confirmed in the Holder table ($d=2$) and Hartmann6 ($d=6$) functions, which are shown in Appendix~\ref{app:additional_exp}.
        }
    \label{fig:pgp_ackley}
\end{figure}

\subsection{Limitation of Gaussian Approximation-Based Preferential GP}
\label{sec:comparison_PGP}

To clarify the inaccuracy of Gaussian approximation, Figure~\ref{fig:examples_preferential GP} shows an illustrative example of the estimation of skew GP via MCMC, LA, and EP.
The MCMC-based estimation uses many (10000) MC samples with 1000 burn-in and 10 thinning generated by Gibbs sampling, whose details are discussed in Section~\ref{sec:Gibbs_sampling} and Appendix~\ref{app:detail_gibbs_sampling}.
Thus, we can expect that the MCMC is sufficiently precise and treat it as the ground truth.
Then, we will confirm the differences between the ground truth (MCMC), and LA and EP, respectively.

First, we can see that the LA is inaccurate since the mode can be very far away from the mean, particularly when $\sigma^2_{\rm noise}$ is small.
A similar problem for GPC is discussed in \citep[][Section~12.2.1]{Nickisch2008-approximations}.
Therefore, although LA is very fast, LA-based preferential BO will fail.

Second, EP is very accurate for estimating the mean and the credible interval of $f | \*v_t < \*0$, which matches the experimental observation in the literature of GPC \citep{Nickisch2008-approximations,Kuss2005-assessing}.
However, in the right plot, $\Pr_t (f(\*x_w) \leq f(\*x_l))$ is overestimated although $\*x_w \succ \*x_l$ is already obtained.
This implies the problem of EP that predictions for the statistics involving the joint distribution of several outputs (e.g., $f(\*x_w)$ and $f(\*x_l)$) can be inaccurate.
Particularly, since some AFs for preferential BO \citep{gonzalez2017preferential,nielsen2015-perception,fauvel2021-efficient} need to evaluate the joint distribution of $f(\*x^{(1)}_t)$ and $f(\*x^{(2)}_t)$, EP can degrade the performance of such AFs.
This is a crucial problem but has not been discussed in the literature on preferential GP \citep{chu2005-extensions,benavoli2021-preferential,benavoli2021-unified} and GPC \citep{Kuss2005-assessing,Nickisch2008-approximations}.

It should be noted that \citet{benavoli2021-unified} already provided the experimental comparison between EP and skew GP.
However, we argue that their setting is not common in the following aspects:
For example, \citep[][Figure~6 and 7]{benavoli2021-unified} considers that all the inputs in a certain interval lose, and all other inputs in another interval win.
Although the Gaussian approximation becomes inaccurate due to the heavy skewness in this setting, these biased training duels are unrealistic.
In contrast, Figure~\ref{fig:examples_preferential GP}, in which the usual uniformly random duels are used, considers a more common setting.

Figure~\ref{fig:pgp_ackley}~(a) further shows the truth (MCMC) and predictions of LA and EP plots using the Ackley function.
We can again confirm that LA underestimates all the statistics, and the mean and mode estimators of EP are very accurate.
However, for the probability of duels, EP underestimates the probabilities around 0 or 1.
Importantly, we can see that the estimator does not retain the relative relationship of the ground truth.
Thus, EP-based AF can select the different inputs from the inputs selected by the ground truth.
Therefore, Figure~\ref{fig:pgp_ackley}~(a) also suggests that EP-based preferential BO can deteriorate.

\subsection{Improving MCMC-based Preferential GP}
The above results motivated us to improve the MCMC-based estimation of the skew GP.
In this section, we show the practical effectiveness of Gibbs sampling and derive the low variance MC estimators.

\subsubsection{Gibbs Sampling for TMVN}

\label{sec:Gibbs_sampling}
In this study, for MCMC of skew GP that requires sampling from TMVN, we argue that the Gibbs sampling should be employed rather than linear elliptical slice sampling (LinESS) \citep{gessner2020-integrals}, which is used in the existing studies~\citep{benavoli2021-preferential,benavoli2021-unified}.
Gibbs sampling is a standard approach for the sampling from TMVN \citep{adler2008-efficient,Breslaw1994-random,geweke1991efficient,kotecha1999gibbs,li2015efficient,robert1995-simulation}.
In Gibbs sampling, we need to sample univariate truncated normal, for which we used the efficient rejection sampling \citep[][Section 2.1]{li2015efficient}, which uses several proposal distributions depending on the truncation.
Detailed procedure and pseudo-code of Gibbs sampling are shown in Appendix~\ref{app:detail_gibbs_sampling}.

Figure~\ref{fig:trace_ackley} shows the trace and autocorrelation plots of LinESS and Gibbs sampling, where we used the implementation by \citep{benavoli2021-preferential,benavoli2021-unified} for LinESS (\url{https://github.com/benavoli/SkewGP}).
The mixing time of Gibbs sampling is very fast compared with LinESS.
Furthermore, as shown in Table~\ref{tab:computation_time_MCMC_ackley}, Gibbs sampling is roughly 3x faster than LinESS.
Consequently, our experimental results suggest that Gibbs sampling is more suitable than LinESS for preferential GP.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.7\linewidth]{manuscripts/fig/Ackley_trace_00.pdf}
    \caption{Trace and autocorrelation plots of LinESS and Gibbs sampling for first 1000 samples used in the experiments in Figure~\ref{fig:pgp_ackley}.}
    \label{fig:trace_ackley}
\end{figure}

\begin{table}[!t]
    \caption{The mean and standard deviation of computational time of MCMC for Figure~\ref{fig:trace_ackley} with 10 random trials.}
    \label{tab:computation_time_MCMC_ackley}
    \centering
    \begin{tabular} {c|cccc}
        & LinESS & Gibbs sampling \\ \hline
        Time (sec) & $1.61 \pm 0.01$ & $0.54 \pm 0.01$
    \end{tabular}
\end{table}

\subsubsection{MC Estimation for Posterior Statistics}
\label{sec:procedure_inference}

To derive the low variance MC estimator, we use the following fact:
\begin{proposition}
    The exact posterior distribution additionally conditioned by $\*v_t$ is a GP.
    %
    That is, for an arbitrarily output vector $\*f_{\rm tes}$, we obtain
    \begin{align}
        \*f_{\rm tes} \mid \*v_t, \*v_t < \*0
        &\stackrel{\rm{d}}{=} \*f_{\rm tes} \mid \*v_t
        \sim {\cal N}(\*\mu_{{\rm tes} | \*v}, \*\Sigma_{{\rm tes} | \*v}), \\
        \*\mu_{{\rm tes} | \*v} &\coloneqq \*\Sigma_{{\rm tes}, \*v} \*\Sigma_{\*v, \*v}^{-1} \*v_t, \label{eq:mu_given_v} \\
        \*\Sigma_{{\rm tes} | \*v} &\coloneqq \*\Sigma_{{\rm tes}, {\rm tes}} - \*\Sigma_{{\rm tes}, \*v} \*\Sigma_{\*v, \*v}^{-1} \*\Sigma_{{\rm tes}, \*v}^\top,
    \end{align}
    where $\stackrel{\rm{d}}{=}$ implies the equality in distribution, and we denote the sub-matrices of Eq.~\eqref{eq:prior_covariance} as follows:
    \begin{align*}
        \*\Sigma
        \eqqcolon
        \begin{bmatrix}
            \*\Sigma_{{\rm tes}, {\rm tes}} & \*\Sigma_{{\rm tes}, \*v} \\
            \*\Sigma_{{\rm tes}, \*v}^\top & \*\Sigma_{\*v, \*v}
        \end{bmatrix},
    \end{align*}
    where $\*\Sigma_{{\rm tes}, {\rm tes}} \in \mathbb{R}^{m \times m}, \*\Sigma_{{\rm tes}, \*v} \in \mathbb{R}^{m \times t}, \*\Sigma_{\*v, \*v} \in \mathbb{R}^{t \times t}$.
    \label{prop:cond_SGP}
\end{proposition}
See Appendix~\ref{app:proof} for the proof.

One of the important consequences of Proposition~\ref{prop:cond_SGP} is that an expectation for the posterior inference can be partially calculated by an analytical form.
For example, the posterior mean $\mathbb{E}_{f | \*v_t < \*0 } [\*f_{\rm tes} ]$ can be derived as,
\begin{align}
    \mathbb{E}_{f | \*v_t < \*0 } [\*f_{\rm tes}]
    &= \mathbb{E}_t \bigl[ \mathbb{E}_{f | \*v_t, \*v_t < \*0} [\*f_{\rm tes} ] \bigr] \\

    &= \mathbb{E}_t [ \*\mu_{{\rm tes} | \*v} ] && (\because \text{Prop.~\ref{prop:cond_SGP}}) \\
    &= \*\Sigma_{{\rm tes}, \*v} \*\Sigma_{\*v, \*v}^{-1} \mathbb{E}_t [\*v_t],
    \label{eq:mean_estimator}
\end{align}
where we denote $\mathbb{E}_t [\cdot] \coloneqq \mathbb{E}_{\*v_t | \*v_t < \*0} [\cdot]$.
Note that, in the first and last equalities, we use the tower property of the expectation and Eq.~\eqref{eq:mu_given_v}, the definition of $\*\mu_{{\rm tes} | \*v}$, respectively.
Other statistics, such as variance and probability, can be estimated via the following derivation:
\begin{align}
    &\mathbb{V}_{f | \*v_t < \*0} [\*f_{\rm tes} ] = \mathbb{V}_t [\*\mu_{{\rm tes} | \*v}] + \*\Sigma_{{\rm tes} | \*v}, \\
    &\mathrm{Pr} \bigl( [\*f_{\rm tes}]_i  \leq c | \*v_t < \*0  \bigr)
    = \mathbb{E}_t \left[ \Phi \left(
    \frac{c - [\*\mu_{{\rm tes} | \*v}]_i}{ \sqrt{[\*\Sigma_{{\rm tes} | \*v}]_{i,i} }}
    \right) \right], \label{eq:CDF_estimator} \\

    &\mathrm{Pr}_t \bigl( [\*f_{\rm tes}]_i  \leq [\*f_{\rm tes}]_j | \*v_t < \*0  \bigr) \\
    &= \mathbb{E}_t \left[ \Phi \left(
    \frac{[\*\mu_{{\rm tes} | \*v}]_j - [\*\mu_{{\rm tes} | \*v}]_i}{ \sqrt{[\*\Sigma_{{\rm tes} | \*v}]_{i,i} + [\*\Sigma_{{\rm tes} | \*v}]_{j,j} - 2 [\*\Sigma_{{\rm tes} | \*v}]_{i,j} }}
    \right) \right],
\end{align}
where $\mathbb{V}$ means the (co)variance, and we denote $\mathbb{V}_t [\cdot] \coloneqq \mathbb{V}_{\*v_t | \*v_t < \*0} [\cdot]$, and $\Phi$ is the CDF of the standard normal distribution.
See Appendix~\ref{app:detailed_estimators} for detailed derivation.
As a result, we only need to apply MC estimation using the samples of $\*v_t$ to $\mathbb{V}_t$ or $\mathbb{E}_t$.
Note that after we obtain the samples of $\*v_t | \*v_t < \*0$ once, we can reuse it for all the estimations.

The remaining important statistic is the quantile of $[\*f_{\rm tes}]_i | \*v_t < \*0$, which is useful to compute the credible interval and UCB-based AFs.
Although this can be derived as an inverse function of Eq.~\eqref{eq:CDF_estimator}, the analytical form of this inverse function is unknown.
Thus, we compute the quantile by applying a binary search to the MC estimator of Eq.~\eqref{eq:CDF_estimator}.
See Appendix~\ref{app:detailed_estimators} for details.

The error of the MC estimator is in proportion to the variance of the MC estimator.
Then, our MC estimators are justified in terms of the error of the MC estimator as follows:
\begin{proposition}
    The variances of our MC estimators are must smaller than those of \citep{benavoli2021-preferential}, such as Eq.~\eqref{eq:samplepath_estimation}.
    %
    For example, the following inequality holds:
    \begin{align}
        \mathbb{V}_t \bigl[  [\*\mu_{{\rm tes} | \*v}]_{i} \bigr] \leq \mathbb{V}_{f | \*v_t < \*0} \bigl[ [\*f_{\rm tes}]_{i} \bigr],
    \end{align}
    for all $i \in \{1, \dots, m\}$.
    \label{prop:low_variance}
\end{proposition}
This proposition is a direct consequence of the low of total variance.
See Appendix~\ref{app:proof_variance} for details.
Thus, in principle, our MC estimators are always superior to those of \citep{benavoli2021-preferential,benavoli2021-unified}.

Figure~\ref{fig:pgp_ackley}~(b) shows the root mean squared error (RMSE) of our MC estimators and those of \citep{benavoli2021-preferential,benavoli2021-unified}, which we refer to as as a full MC estimator.
Note that the RMSE is computed using the ground truth.
Thus, if the ground truth is sufficiently accurate (equal to the true expectation), the RMSE can be interpreted as the standard deviation estimator of the MC estimator.
We can observe that the RMSE of our estimator is smaller than that of the full MC estimator in all the statistics and all the number of MC samples.
Therefore, the efficiency of our MC estimator was verified also experimentally.

\section{Preferential BO with Hallucination Believer}
\label{sec:proposed}

Although we tackled the speed-up of MCMC-based estimation, its computational time can still be unrealistic in preferential BO.
Furthermore, Gaussian approximation often suffers poor prediction performance ignoring the skewness of the posterior.
In this work, we propose a computationally efficient preferential BO method while considering the skewness in a randomized manner.

\subsection{Hallucination believer}

We proposed a randomized preferential BO method, called hallucination believer (HB), which uses the AFs for the standard BO calculated from GP $f | \*v_t = \tilde{\*v}_t$, where $\tilde{\*v}_t \sim p(\*v_t | \*v_t < \*0)$ is an i.i.d. sample and called \emph{hallucination}.
There are two key points: the conditioning on $\*v_t$ and using the hallucination $\tilde{\*v}_t$ for conditioning.
First, as shown in Proposition~\ref{prop:cond_SGP}, the conditioning on $\*v_t$ reduces the posterior to the GP $f | \*v_t = \tilde{\*v}_t$, which enables us to use any computationally efficient and powerful AF for the standard BO.
Second, by using hallucination $\tilde{\*v}_t$ from the true posterior, we aim to incorporate the skewness of the posterior in a randomized manner.
This is because, if we condition as $ \*v_t = \*c$, where $\*c$ is some constant vector such as the posterior mean, the proposed method ignores skewness, and the optimization performance will deteriorate.
Consequently, we need only one sample from TMVN $\tilde{\*v}_t$ and the usual GP-based AF calculations.
Therefore, the computational complexity is significantly cheaper than MCMC-based AFs \citep{benavoli2021-preferential}.
We will show the effectiveness of this hallucination-based randomization via extensive experiments in Section~\ref{sec:exp}.

\begin{algorithm}[t]
    \caption{Hallucination believer for preferential BO}\label{alg:HB}
    \begin{algorithmic}[1]
        \REQUIRE ${\cal D}_0 = \{\*x_{0,w}, \*x_{0,l} \}$, ${\cal X}$
        \FOR{$t = 1, \dots$}
            \STATE $\*x_{t}^{(1)} \gets \*x_{t-1,w}$
            \STATE Generate  $\tilde{\*v}_{t-1}$ from $p(\*v_{t-1}\mid \*v_{t-1} < \*0)$
            \STATE $\*x_{t}^{(2)} \gets \mathop{\text{argmax}}_{\*x \in {\cal X}} \alpha(\*x)$ based on GP $f \mid \tilde{\*v}_{t-1}$
            \STATE Set $\*x_{t,w}$ and $\*x_{t,l}$ as the winner and loser of the duel between $\*x_{t}^{(1)}$ and $\*x_{t}^{(2)}$, respectively
        \ENDFOR
    \end{algorithmic}
\end{algorithm}

Algorithm~\ref{alg:HB} shows the procedure of HB.
HB iteratively selects a pair of inputs by the following two steps:
(i) Select the winner of the past duels as the first input $\*x_{t}^{(1)}$ (line 2) as with \citep{benavoli2021-preferential}.
(ii) Using the posterior distribution conditioned by the hallucination, select the input that maximizes the AF $\alpha: {\cal X} \rightarrow \mathbb{R}$ as the second point $\*x_{t}^{(2)}$ (lines 3-4).

\paragraph{Choice of AF $\alpha$ for HB:}
For HB, we can use an arbitrary AF for the standard BO, e.g., EI \citep{Mockus1978-Application} and UCB \citep{Srinivas2010-Gaussian}.
However, HB combining Thompson sampling (TS) is meaningless since the sampling sequentially from $p(\*v_t | \*v_t < \*0)$ and $p(f | \*v_t)$ is just the TS.
On the other hand, AFs proposed for preferential BO, such as bivariate EI \citep{nielsen2015-perception} and maximally uncertain challenge (MUC) \citep{fauvel2021-efficient}, which aims for appropriate diversity by considering the correlation between $\*x^{(1)}_t$ and $\*x^{(2)}_t$, should not be integrated into HB.
This is because, in HB, diversity is achieved by hallucination-based randomization.
Thus, HB combining bivariate EI or MUC results in over-exploration.
Consequently, we use EI and UCB in our experiments.

\subsection{Relation to Existing Methods}
\label{sec:relation_TS}

HB and TS have a relationship in the perspective of using a random sample from the posterior.
HB used the hallucination from $p(\*v_t | \*v_t < \*0)$.
We can interpret that TS uses the most uncertain hallucination, i.e., sample path from $p(f | \*v_t < \*0)$.
Thus, the uncertainty of TS is very large, and TS often results in over-exploration.
On the other hand, the hallucination $\tilde{\*v}_t$ is the least uncertain random variable in the random variables, which reduces the posterior to the GP.
By this difference, HB has lower uncertainty than TS and alleviates over-exploration, particularly when the input dimension is large.

HB also has a relationship to the kriging believer \citep{Shahriari2016-Taking}, which is a well-known heuristic in the literature of parallel BO.
This relation is discussed in Appendix~\ref{app:diff_hallucination_kriging}.

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

Although dueling bandit \citep{yue2012-karmed} considers online learning from the duels, it does not consider the correlation between arms in general, as reviewed in \citep{sui2018-advancements}.
In contrast, preferential BO~\citep{brochu2010-tutorial} aims for more efficient optimization using the preferential GP, which can capture a flexible correlational structure in an input domain.
Since the true posterior is computationally intractable, most prior works employed Gaussian approximation-based preferential GP models \citep{chu2005-preference,chu2005-extensions}.
Typically, preferential BO sequentially duels $\*x_t^{(1)}$, which is the most exploitative point (e.g., maxima of the posterior mean), and $\*x_t^{(2)}$, which is the maxima of some AFs.

The design of AFs for $\*x_t^{(2)}$ is roughly twofold.
First approaches \citep{brochu2010-tutorial,siivola2021-preferential} use the AFs for standard BO, which assumes that the direct observation $y$ can be obtained (e.g., EI \citep{Mockus1978-Application}).
However, this approach selects similar duels repeatedly since the preferential GP model's variance is hardly decreased due to the weakness of the information obtained from the duel.
Second approach considers the relationship between $\*x_t^{(1)}$ and $\*x_t^{(2)}$.
For example, bivariate EI \citep{nielsen2015-perception} uses the improvement from $f(\*x_t^{(1)})$, and MUC \citep{gonzalez2017preferential,fauvel2021-efficient} employed the uncertainty of the probability that $\*x_t^{(2)}$ wins $\*x_t^{(1)}$.
Considering the correlation to $\*x_t^{(1)}$, these AFs are expected to explore appropriately.
However, as discussed in Section~\ref{sec:comparison_PGP}, Gaussian approximation (LA and EP) can be inaccurate for evaluation of the joint distribution of $f(\*x_t^{(1)})$ and $f(\*x_t^{(2)})$.

Instead of using Gaussian approximation, \citet{benavoli2021-preferential} directly estimates the true posterior using MCMC.
However, MCMC-based AFs require heavy computational time, which can be a bottleneck in practical applications.

Another well-used AF for preferential BO is TS \citep{Russo2014-learning}.
\citet{gonzalez2017preferential} proposed to select $\*x_t^{(1)}$ by TS, and \citet{siivola2021-preferential,sui2017multi} proposed to select both $\*x_t^{(1)}$ and $\*x_t^{(2)}$ by TS, both of which are based on Gaussian approximation.
\citet{benavoli2021-preferential} also proposed to select $\*x_t^{(2)}$ by TS, which is generated from skew GP.
Although the merit of Gaussian approximation is low computational time, the posterior sampling from skew GP only once is sufficiently fast (comparable to EP).
Therefore, the advantage of Gaussian approximation for TS-based AF is limited.
In Section~\ref{sec:exp}, we will show that TS can cause over-exploration when the input dimension is large in preferential BO, as with the standard BO reviewed in \citep{Shahriari2016-Taking}.

\section{Comparison of Preferential BO Methods}
\label{sec:exp}

We investigate the effectiveness of the proposed preferential BO method through comprehensive numerical experiments.
We employed the 12 benchmark functions.
In this section, we show the results for the 8 functions, called the Branin, Holder table, Bukin, Eggholder, Ackley, Harmann3, Hartmann4, and Hartmann6; others are shown in Appendix~\ref{app:additional_exp}.

We performed HB combining EI \cite{Mockus1978-Application} and UCB \cite{Srinivas2010-Gaussian}, denoted as HB-EI and HB-UCB, respectively.
We employed the baseline methods, LA-EI\citep{brochu2010-tutorial}, EP-EI \cite{siivola2021-preferential}, EP-MUC~\cite{fauvel2021-efficient}, DuelTS~\cite{benavoli2021-preferential}, DuelUCB~\cite{benavoli2021-preferential}, and EIIG\footnote{
\citet{benavoli2021-preferential,benavoli2021-unified} proposed EIIG that is EI minus information gain (IG), which is probably a typo since both EI and IG should be large. Thus, we used a modified one that EI plus IG.
}
\cite{benavoli2021-preferential}, where DuelTS, DuelUCB, and EIIG were based on skew GP and other methods with Suffix ``EP'' and ``LA'' employed EP and LA, respectively.
The first input $\*x^{(1)}_t$ for EP-EI, EP-MUC, and LA-EI is the maxima of the posterior mean over the training duels, and other methods, including HB, use the winner so far as $\*x^{(1)}_t$, following \citep{benavoli2020skew}.
For preferential GP models, we use RBF kernel with automatic relevance determination \citep{Rasmussen2005-Gaussian}, whose lengthscales are selected by marginal likelihood maximization per 10 iterations, and set fixed noise variance $\sigma^2_{\rm noise}=10^{-4}$.
We use marginal likelihood estimation by LA, which is sufficiently precise if $\sigma^2_{\rm noise}$ is small, as discussed in \citep[][Section~7.2]{Kuss2005-assessing}.
We show the accuracy of the marginal likelihood estimations of LA and EP in Appendix~\ref{app:additional_exp}.
Other detailed settings are shown in Appendix~\ref{app:exp_settings}.

Figure~\ref{fig:regret} show the results for (a) computational time and (b) iterations.
As a performance measure, we used the regret defined as $f(\*x_*) - f(\tilde{\*x}_t)$, where $\tilde{\*x}_t$ is a recommendation point at $t$-th iteration, for which we used $\*x^{(1)}_t$.
We report the mean and standard error of the regret over 10 random initialization, where the initial duel is $3d$ uniformly random input pair.

\begin{figure*}[!th]
    \centering
    \includegraphics[width=0.75\linewidth]{manuscripts/fig/Results_legend.pdf}

    \includegraphics[width=0.02\linewidth]{manuscripts/fig/regret_axis_label.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Branin.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Holder_table.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Bukin.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Eggholder.pdf}
    \hfill
    \subfigure[Computational time (sec)]{
        \includegraphics[width=0.02\linewidth]{manuscripts/fig/regret_axis_label.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Ackley.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Hartmann3.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Hartmann4.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Hartmann6.pdf}
    }
    \hfill
    \centering
    \includegraphics[width=0.02\linewidth]{manuscripts/fig/regret_axis_label.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Branin.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Holder_table.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Bukin.pdf}
    \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Eggholder.pdf}
    \hfill
    \subfigure[Iteration]{
        \includegraphics[width=0.02\linewidth]{manuscripts/fig/regret_axis_label.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Ackley.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Hartmann3.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Hartmann4.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Hartmann6.pdf}
    }

    \caption{
    Average and standard error of the regret of the proposed methods (HB-EI and HB-UCB) and the state-of-the-art preferential BO methods.
    The horizontal axis represents (a) the computational time (sec) and (b) the number of iterations.
    The vertical axis represents the regret, which is the smaller, the better it is.}

    \label{fig:regret}
\end{figure*}

First, we focus on the comparison against Gaussian approximation-based methods.
LA-EI is the worst method due to the poor approximation accuracy.
Note that, since we set the plotted interval for the small regret to focus the differences between HB and baseline methods except for LA-EI, LA-EI is out of the range in Hartmann3, 4, and 6 functions.
The regret of EP-EI, which selects the duel between the maximums of the posterior mean and EI, is often stagnated due to over-exploitation.
The regret of EP-MUC also stagnates in the Branin, Bukin, and Ackley functions since EP is often inaccurate in predicting the joint distribution of several outputs, as shown in Section~\ref{sec:comparison_PGP}.
Overall, incorporating the skewness of the true posterior, HB-EI and HB-UCB outperformed Gaussian approximation-based methods in Figure~\ref{fig:regret} (b).
Furthermore, the same tendency can be confirmed also in Figure~\ref{fig:regret} (a) since the computational time of HB-EI and HB-UCB is comparable with Gaussian approximation.

Second, DuelUCB and EIIG are MCMC-based AFs based on skew GP.
Therefore, although they often show rapid convergence in Figure~\ref{fig:regret}(b) because of the use of skew GP, they show slow convergence in Figure~\ref{fig:regret}(a) due to the heavy computational time of MCMC.
Hence, HB-EI and HB-UCB are comparable to and outperformed MCMC-based AFs in Figure~\ref{fig:regret}(b) and (a), respectively, which suggest the practical effectiveness of the proposed methods.

Last, since DuelTS requires only one sample path, its performance is not largely changed in Figure~\ref{fig:regret} (a) and (b).
Although the regret of DuelTS converges rapidly in low-dimensional functions, it converges considerably slowly in relatively high-dimensional functions, such as Hartmann4 and Hartmann6, in which the proposed methods clearly outperformed DuelTS.
In our experiments, we confirmed that TS results in over-exploration in relatively high-dimensional functions, which is the well-known problem of TS for the standard BO, as reviewed in \citep{Shahriari2016-Taking}.

\section{Conclusion}
Towards building a more practical preferential BO, this work has made the following important progresses.
We first revealed the poor performance of the Gaussian approximation for skew GP, such as Laplace approximation and expectation propagation, despite being the method most widely used in preferential BO.
This result motivated us to improve the MCMC-based estimation of the skew GP, for which we showed the practical effectiveness of Gibbs sampling and derived the low variance MC estimators.
Finally, we developed a new preferential BO method that achieves both high computational efficiency and low sample complexity, and then verified its effectiveness through extensive numerical experiments.
The regret analysis of the proposed method is one of the interesting future works.

\section{Proof of Proposition~\ref{prop:cond_SGP}}
\label{app:proof}

We here prove Proposition~\ref{prop:cond_SGP}.
The joint posterior distribution of $\*f_{\rm tes}$ and $\*v_t$ is truncated MVN, in which $\*v_t$ is truncated above at $\*0$, as follows:
\begin{align*}
    p(\*f_{\rm tes}, \*v_t | \*v_t < \*0)
    &= \frac{\Pr(\*v_t < \*0 | \*f_{\rm tes}, \*v_t) p(\*f_{\rm tes}, \*v_t) }{\Pr(\*v_t < \*0)}.
\end{align*}
Then, we can see that $\Pr(\*v_t < \*0 | \*f_{\rm tes}, \*v_t) = \mathbbm{1} \{ \*v_t < \*0 \}$, where $\mathbbm{1} \{ \*v_t < \*0 \} = 1$ if $\*v_t < \*0$, and 0 otherwise.
Hence, we can obtain
\begin{align*}
    p(\*f_{\rm tes}, \*v_t | \*v_t < \*0)
    &= \frac{\mathbbm{1} \{ \*v_t < \*0 \} p(\*f_{\rm tes}, \*v_t) }{\Pr(\*v_t < \*0)} \\
    &= p(\*f_{\rm tes} | \*v_t) \frac{\mathbbm{1} \{ \*v_t < \*0 \} p(\*v_t) }{\Pr(\*v_t < \*0)}.
\end{align*}
Furthermore, we know that $p(\*v_t | \*v_t < \*0) = \frac{\mathbbm{1} \{ \*v_t < \*0 \} p(\*v_t) }{\Pr(\*v_t < \*0)}$.
Therefore, we derive
\begin{align*}
    p(\*f_{\rm tes}, \*v_t | \*v_t < \*0)
    &= p(\*f_{\rm tes} | \*v_t) p(\*v_t | \*v_t < \*0),
\end{align*}
which shows $p(\*f_{\rm tes} | \*v_t) = p(\*f_{\rm tes} | \*v_t, \*v_t < \*0)$.
Since $\*f_{\rm tes}$ is arbitrary, $f \mid \*v_t$ is a GP.

We can obtain other proof.
From the property of truncated MVN \citep[][Conclusion~5]{Horrace2005-Some}, the conditional distribution of truncated MVN is truncated MVN keeping the original truncation, in which the parameters can be computed as with usual MVN.
Therefore, since remained variable $\*f_{\rm tes}$ is not truncated, the conditional distribution $p(\*f_{\rm tes} \mid \*v_t, \*v_t < \*0)$ is MVN:
\begin{align*}
    \*f_{\rm tes} \mid \*v_t, \*v_t < \*0
    \sim {\cal N}(\*\mu_{{\rm tes} | \*v }, \*\Sigma_{{\rm tes} | \*v}),
\end{align*}
which is equivalent to the distribution $p(\*f_{\rm tes} \mid \*v_t)$.
Hence, we can see that $f \mid \*v_t$ is a GP.

\section{Detailed Derivation of Statistics for Skew GP}
\label{app:detailed_estimators}

\paragraph{Variance:}
From the low of total variance, we can decompose the variance as follows:
\begin{align*}
    \mathbb{V}_{f | \*v_t < \*0}[\*f_{\rm tes}]
    &= \mathbb{V}_t[ \mathbb{E}_{f | \*v_t} [\*f_{\rm tes}]] + \mathbb{E}_t[\mathbb{V}_{f | \*v_t}[\*f_{\rm tes}]] \\
    &= \mathbb{V}_t[ \*\mu_{{\rm tes} | \*v }] + \mathbb{E}_t[\*\Sigma_{{\rm tes} | \*v}] && (\because \text{Proposition~\ref{prop:cond_SGP}}) \\
    &= \mathbb{V}_t[ \*\mu_{{\rm tes} | \*v }] + \*\Sigma_{{\rm tes} | \*v}. && (\because \text{$\*\Sigma_{{\rm tes} | \*v}$ does not depend on $\*v_t$.})
\end{align*}

\paragraph{CDF of $f(\*x)$:}
We can derive as follows:
\begin{align*}
    \mathrm{Pr} \bigl( [\*f_{\rm tes}]_i \leq c  | \*v_t < \*0 \bigr)
    &= \mathbb{E}_t \left[ \Pr \bigl( [\*f_{\rm tes}]_i \leq c | \*v_t, \*v_t < \*0 \bigr) \right] \\
    &= \mathbb{E}_t \left[ \Phi \left( \frac{c - [\*\mu_{{\rm tes} | \*v}]_i}{ \sqrt{[\*\Sigma_{{\rm tes} | \*v}]_{i,i} }} \right) \right].&& (\because \text{Proposition~\ref{prop:cond_SGP}})
\end{align*}

\paragraph{Probability of duel:}
Given $\*v_t$, $[\*f_{\rm tes}]_i$ and $[\*f_{\rm tes}]_j$ follow MVN jointly.
Therefore, $[\*f_{\rm tes}]_i - [\*f_{\rm tes}]_j | \*v_t$ also follows the Gaussian distribution ${\cal N}([\*\mu_{{\rm tes} | \*v}]_i - [\*\mu_{{\rm tes} | \*v}]_j, [\*\Sigma_{{\rm tes} | \*v}]_{i,i} + [\*\Sigma_{{\rm tes} | \*v}]_{j,j} - 2 [\*\Sigma_{{\rm tes} | \*v}]_{i,j})$.
Then, we can obtain the estimator as follows:
\begin{align*}
    \mathrm{Pr} \bigl( [\*f_{\rm tes}]_i  \leq [\*f_{\rm tes}]_j | \*v_t < \*0 \bigr)
    &= \mathrm{Pr} \bigl( [\*f_{\rm tes}]_i - [\*f_{\rm tes}]_j \leq 0 | \*v_t < \*0 \bigr) \\
    &= \mathbb{E}_t \left[ \mathrm{Pr} \bigl( [\*f_{\rm tes}]_i - [\*f_{\rm tes}]_j \leq 0 | \*v_t, \*v_t \leq \*0 \bigr) \right] \\
    &= \mathbb{E}_t \left[ \Phi \left(
    \frac{[\*\mu_{{\rm tes} | \*v}]_j - [\*\mu_{{\rm tes} | \*v}]_i}{ \sqrt{[\*\Sigma_{{\rm tes} | \*v}]_{i,i} + [\*\Sigma_{{\rm tes} | \*v}]_{j,j} - 2 [\*\Sigma_{{\rm tes} | \*v}]_{i,j} }}
    \right) \right]. && (\because \text{Proposition~\ref{prop:cond_SGP}})
\end{align*}

\paragraph{Quantile:}
As mentioned in the main paper, the quantile estimators can be defined as an inverse function of the MC estimator of Eq.~\eqref{eq:CDF_estimator}.
Let the estimator of Eq.~\eqref{eq:CDF_estimator} be $\hat{F}_i: \mathbb{R} \rightarrow [0, 1]$.
Unfortunately, since $\hat{F}_i$ is defined as the mean of the CDF of Gaussian distribution, the analytical form of $\hat{F}_i^{(-1)}$ is unknown.
Therefore, we apply the binary search to $\hat{F}_i$.
We define the initial search interval from the $M$ MC samples $\tilde{\*v}_{t, k}$, where $k \in \{1, \dots, M\}$.
From the construction of $\hat{F}_i$, we see that
\begin{align*}
    \hat{F}_i^{(-1)} (\alpha)
    \in [\ell, u]
    \coloneqq
        \left[
            \mathop{\rm argmin}\limits_k [\*\mu_{{\rm tes} | \tilde{\*v}_{t,k}}]_i + c \sqrt{[\*\Sigma_{{\rm tes} | \tilde{\*v}_{t,k}}]_{i,i} },
            \mathop{\rm argmax}\limits_k [\*\mu_{{\rm tes} | \tilde{\*v}_{t,k}}]_i + c \sqrt{[\*\Sigma_{{\rm tes} | \tilde{\*v}_{t,k}}]_{i,i} }
        \right],
\end{align*}
where $\alpha \in (0, 1)$, $c = \Phi^{-1}(\alpha)$, and $\ell$ and $u$ are the minimum and maximum values of the quantiles of $[\*f_{\rm tes}]_i | \tilde{\*v}_{t, k}$ over $k \in \{1, \dots, M\}$, respectively.
Therefore, we use $[\ell, u]$ as the initial interval.
Algorithm~\ref{alg:quantile} shows the pseudo-code.

\begin{algorithm}[t]
    \caption{Binary search for quantiles of skew GP} \label{alg:quantile}
    \begin{algorithmic}[1]
        \REQUIRE $\hat{F}_i: \mathbb{R} \rightarrow [0, 1]$, the desired percentage $\alpha \in (0, 1)$, tolerance $\epsilon$
        \STATE $c \gets \Phi^{-1}(\alpha)$
        \STATE $[\ell, u] \gets \left[ \mathop{\rm argmin}\limits_k [\*\mu_{{\rm tes} | \tilde{\*v}_{t,k}}]_i + c \sqrt{[\*\Sigma_{{\rm tes} | \tilde{\*v}_{t,k}}]_{i,i} }, \mathop{\rm argmax}\limits_k [\*\mu_{{\rm tes} | \tilde{\*v}_{t,k}}]_i + c \sqrt{[\*\Sigma_{{\rm tes} | \tilde{\*v}_{t,k}}]_{i,i} } \right]$
        \WHILE{$|\hat{F}_i(\gamma) - \alpha| > \epsilon$}
            \STATE $\gamma \gets \hat{F}_i\bigl( (\ell + u) / 2 \bigr)$
            \IF{$\gamma > \alpha$}
                \STATE $u \gets (\ell + u) / 2$
            \ELSE
                \STATE $\ell \gets (\ell + u) / 2$
            \ENDIF
        \ENDWHILE
        \ENSURE $\gamma$
    \end{algorithmic}
\end{algorithm}

\section{Proof of Proposition~\ref{prop:low_variance}}
\label{app:proof_variance}

We show the general fact for the MC estimator.
Let $X$ and $Y$ be random variables, where $Y$ depends on $X$.
Furthermore, we assume the conditional expectation $\mu_{Y|X} \coloneqq \mathbb{E}_{Y | X} [Y]$ can be analytically calculated from a realization of $X$.
Then, $Y$ itself and $\mu_{Y|X}$ are both unbiased estimators of $\mathbb{E}[Y]$, i.e., $\mathbb{E}_Y[Y] = \mathbb{E}_X[\mu_{Y|X}]$.
Intuitively, since $\mu_{Y|X}$ is computed analytically with respect to the expectation of $Y | X$, it should be a better estimator than $Y$.
This intuition is justified by the low of total variance shown below:
\begin{align*}
    \mathbb{V}[Y] = \mathbb{V}_X[\mu_{Y|X}] + \mathbb{E}_X[\mathbb{V}_{Y|X}[Y]].
\end{align*}
Since $\mathbb{E}_X[\mathbb{V}_{Y|X}[Y]] \geq 0$, we see $\mathbb{V}_X[\mu_{Y|X}] \leq \mathbb{V}[Y]$.
In the case of our mean estimator, $X$, $Y$, and $\mu_{Y|X}$ correspond to $\*v_t$, $f(\*x_{i, {\rm tes}})$, and $\mu_{{\rm tes} | \*v_t}$.
This concludes the proof.

\section{Detailed Procedure of Gibbs Sampling for TMVN}
\label{app:detail_gibbs_sampling}

Let us consider the sampling from TMVN $p(\*v | \*v < \*0)$, where original $\*v \in \mathbb{R}^{n}$ follows MVN below:
\begin{align*}
    \*v \sim {\cal N}(0, \*\Sigma),
\end{align*}
where $\*\Sigma \in \mathbb{R}^{n\times n}$ is an arbitrary covariance matrix.
Let $v_j$ and $\*v_{-j}$ be $j$-th element of $\*v$ and the vector consisting of the elements except for $v_j$, respectively.
In Gibbs sampling, we repeat the sampling $v_j \mid \*v_{-j}, v_j < 0$, which follows a univariate truncated normal distribution.
The conditional distribution $v_j \mid \*v_{-j} \sim {\cal N}(\mu_j, \sigma^{2}_j)$, where $\mu_j$ and $\sigma^2_j$ are computed efficiently by computing $\*\Sigma^{-1}$ once \citep[][Section~5.4.2]{Rasmussen2005-Gaussian}.
Algorithm~\ref{alg:GS} shows the procedure of Gibbs sampling.
For the sampling from the univariate truncated normal distribution, we employed the efficient rejection sampling \citep[][Section~2.1]{li2015efficient}, which uses several proposal distributions depending on the truncation.
Note that we did not employ the transformation of $\*v$, which is proposed in \citep[][Section 2.2]{li2015efficient} to improve the mixing time.
This is because we cannot confirm the effectiveness of this transformation due to additional computational time for transformation.

\begin{algorithm}[th]
    \caption{Gibbs sampling for TMVN}\label{alg:GS}
    \begin{algorithmic}
        \REQUIRE $\*v^{(0)}$, $\*\Sigma$
        \STATE Compute $\*\Sigma^{-1}$
        \FOR{$i = 1, \dots$}
            \STATE $\*v^{(i)} \gets \*v^{(i-1)}$
            \FOR{$j = 1, \dots, n$}
                \STATE $\mu^{(i)}_j \gets [\Sigma^{-1} \*v^{(i)}]_j / [\Sigma^{-1}]_{jj}$

                \STATE Set $v^{(i)}_j$ by the sampling from ${\cal N}(\mu^{(i)}_j, \sigma^{2}_j = 1 / [\Sigma^{-1}]_{jj})$ with truncation above at $\*0$
            \ENDFOR
        \ENDFOR
    \end{algorithmic}
\end{algorithm}

\section{Difference between hallucination and kriging believer}
\label{app:diff_hallucination_kriging}

Kriging believer (KB) \citep{Shahriari2016-Taking} is a well-known heuristic for parallel BO.
KB conditions on the ongoing function evaluation by the posterior mean, and then the next batch point is selected using this GP conditioned by the posterior mean.
Some variants, called constant liar, use some predefined constant instead of the posterior mean.
Furthermore, some studies averaged the resulting AF value by the sample of the posterior normal distribution \citep[e.g., ][]{Snoek2012-Practical}.
These studies aim to guarantee the diversity of the batch points via penalization by conditioning.

On the other hand, one of the important aims of HB is to reduce the skew GP to the standard GP.
For this purpose, we conditioned the latent truncated variable $\*v_t$, which is not related to parallel BO.
Furthermore, in preferential BO, if we conditioned the constant including the posterior mean, preferential BO methods cannot consider the skewness.
Thus, the conditioning by the constant results in poor performance as with Gaussian approximation.
On the other hand, although averaging by the samples from the posterior is promising, it requires huge computational time for MCMC with respect to skew GP.
Hence, we employed the conditioning by the hallucination of $\*v_t$.

\section{Detailed Experimental Settings}
\label{app:exp_settings}

All the details of benchmark functions are shown in \url{https://www.sfu.ca/~ssurjano/optimization.html}.

\paragraph{Comparison with respect to preferential GPs:}
We fitted LA-, EP-, and MCMC-based preferential GP with RBF kernel to 50 uniformly random duels in 10 random trials.
We set $\sigma^2_{\rm noise} = 10^{-4}$, and the hyperparameters of RBF kernel are determined by marginal likelihood maximization \citep{Rasmussen2005-Gaussian}, for which we use LA-based marginal likelihood estimation.
We evaluate the estimators of the mean, mode, and the duel probabilities $\mathrm{Pr}_t\bigl( f(\*x) \succ f(\*x^\prime) \bigr)$ at 100 uniformly random inputs and inputs included in the training duels.
The MCMC uses 10000 MC samples, 1000 burn-in, and 10 thinning.

\paragraph{Comparison of preferential BO methods:}
Although \cite{benavoli2021-preferential} originally employed LinESS for the sampling from truncated MVN, we employed Gibbs sampling in DuelTS, DuelUCB, and EIIG, as with our proposed methods, for a fair comparison.
For the parameters for Gibbs sampling, burn-in is 1000, and the MC sample size for DuelUCB and EIIG is 1000 (thinning is not performed).
Other settings for existing methods, such as the percentage for UCB, are set as with the suggestion from the original paper.
For HB-UCB, we use $\beta^{1/2} = 2$.

\section{Additional Experiments}
\label{app:additional_exp}
In this section, we provide additional experimental results.

\subsection{Marginal Likelihood Approximation}

The marginal likelihood of the preferential GP model is $\Pr(\*v_t < \*0)$, the CDF of MVN.
LA and EP can provide computationally efficient approximations of this.
In the literature of GPC, \citet{Kuss2005-assessing} showed the broad evaluations for the accuracy of these approximations, which suggest that EP is accurate.
Furthermore, \citet[][Section~7.2]{Kuss2005-assessing} described that LA is also accurate if the noise variance $\sigma^2_{\rm noise} < 1$.
We here show the examples using the Branin function in Figure~\ref{fig:MLE}, in which the marginal likelihood approximations by Scipy CDF of MVN, LA, and EP are shown.
We can see that LA and EP show reasonable approximations for the marginal likelihood.
In addition, the computational time shown in the title implies the efficiency of LA.
Therefore, we employed the LA-based marginal likelihood approximation throughout the paper.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{manuscripts/fig/log_cdf_test-50.pdf}
    \caption{Log marginal likelihood estimations for Branin function by Scipy CDF of MVN, LA, and EP. Vertical and horizontal axes represent lengthscales of RBF kernel. In titles, averages of computational times are shown.}
    \label{fig:MLE}
\end{figure}

\subsection{Approximation of Preferential GP for Other Benchmark Functions}
\label{app:comparison_sampling}

We here show the additional experiments regarding preferential GP, which include the comparisons between Gibbs sampling and LinESS,  MCMC and LA and EP, and our MC estimator and the full MC estimator, using the Holder table and Hartmann6 functions.
The experimental settings are the same as the experiments for the Ackley function.

Figures~\ref{fig:pgp_holder_table} and \ref{fig:pgp_hartmann6} show the results of the Holder table and Hartmann6 functions, respectively.
We can observe the same tendency as the Ackley function.
Only in Figure~\ref{fig:pgp_holder_table}~(b), the difference between the RMSE of MC estimators is lower than those of other functions.
This is because, since the dimension of the Holder table functions is relatively low, the remaining uncertainty of $f | \*v_t$ is smaller than that in other functions.
Furthermore, we can confirm the computational efficiency of Gibbs sampling also for the Holder table and Hartmann6 functions, as shown in Table~\ref{tab:computation_time_MCMC}.

\begin{table}[!t]
    \caption{The mean and standard deviation of computational time of MCMC for Figure~\ref{fig:pgp_holder_table} (a) and Figure~\ref{fig:pgp_hartmann6} (a) with 10 random trials.}
    \label{tab:computation_time_MCMC}
    \centering
    \begin{tabular} {c|cccc}
        Time (sec) & LinESS & Gibbs sampling \\ \hline
        Holder table & $1.71 \pm 0.00$ & $0.55 \pm 0.02$ \\
        Hartmann6 & $1.64 \pm 0.00$ & $0.56 \pm 0.00$
    \end{tabular}
\end{table}

\begin{figure}[t]
    \begin{minipage}{0.49\linewidth}
        \centering
        \subfigure[Trace and autocorrelation plots of proposed MC estimator and full MC estimator \citep{benavoli2021-preferential}.]{\includegraphics[width=\linewidth]{manuscripts/fig/Holder_table_trace_00.pdf}}

        \subfigure[RMSE of proposed MC estimator and full MC estimator \citep{benavoli2021-preferential} against ground truth.]{
            \centering
            \includegraphics[width=0.345\linewidth]{manuscripts/fig/Holder_table_mean_error.pdf}
            \includegraphics[width=0.317\linewidth]{manuscripts/fig/Holder_table_mode_error.pdf}
            \includegraphics[width=0.318\linewidth]{manuscripts/fig/Holder_table_prob_error.pdf}
        }

        \centering
        \includegraphics[width=0.32\linewidth]{manuscripts/fig/Holder_table_LA_mean.pdf}
        \includegraphics[width=0.32\linewidth]{manuscripts/fig/Holder_table_LA_mode.pdf}
        \includegraphics[width=0.32\linewidth]{manuscripts/fig/Holder_table_LA_prob.pdf}
        \subfigure[Truth vs. prediction plots for LA and EP.]{
            \centering
            \includegraphics[width=0.32\linewidth]{manuscripts/fig/Holder_table_EP_mean.pdf}
            \includegraphics[width=0.32\linewidth]{manuscripts/fig/Holder_table_EP_mode.pdf}
            \includegraphics[width=0.32\linewidth]{manuscripts/fig/Holder_table_EP_prob.pdf}
        }
        \caption{Reuslts of the Holder table function ($d=2$).}
        \label{fig:pgp_holder_table}
    \end{minipage}
    \hspace{0.02\linewidth}
    \begin{minipage}{0.49\linewidth}
        \centering
        \subfigure[Trace and autocorrelation plots of proposed MC estimator and full MC estimator \citep{benavoli2021-preferential}.]{\includegraphics[width=\linewidth]{manuscripts/fig/Hartmann6_trace_00.pdf}}

        \subfigure[RMSE of proposed MC estimator and full MC estimator \citep{benavoli2021-preferential} against ground truth.]{
            \centering
            \includegraphics[width=0.345\linewidth]{manuscripts/fig/Hartmann6_mean_error.pdf}
            \includegraphics[width=0.317\linewidth]{manuscripts/fig/Hartmann6_mode_error.pdf}
            \includegraphics[width=0.318\linewidth]{manuscripts/fig/Hartmann6_prob_error.pdf}
        }

        \centering
        \includegraphics[width=0.32\linewidth]{manuscripts/fig/Hartmann6_LA_mean.pdf}
        \includegraphics[width=0.32\linewidth]{manuscripts/fig/Hartmann6_LA_mode.pdf}
        \includegraphics[width=0.32\linewidth]{manuscripts/fig/Hartmann6_LA_prob.pdf}
        \subfigure[Truth vs. prediction plots for LA and EP.]{
            \centering
            \includegraphics[width=0.32\linewidth]{manuscripts/fig/Hartmann6_EP_mean.pdf}
            \includegraphics[width=0.32\linewidth]{manuscripts/fig/Hartmann6_EP_mode.pdf}
            \includegraphics[width=0.32\linewidth]{manuscripts/fig/Hartmann6_EP_prob.pdf}
        }
        \caption{Reuslts of the Hartmann6 function ($d=2$).}
        \label{fig:pgp_hartmann6}
    \end{minipage}
\end{figure}

\subsection{Regret Evaluation for Other Benchmark Functions}

Figure~\ref{fig:app-regret} shows the results for the Cross in tray, Langerman, Levy13, and Levy functions.
In these plots, we can still confirm that the proposed methods, HB-EI and HB-UCB, show superior performance in terms of both computational time and iteration.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.75\linewidth]{manuscripts/fig/Results_legend.pdf}

    \subfigure[Computational time (sec)]{
        \includegraphics[width=0.02\linewidth]{manuscripts/fig/regret_axis_label.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Cross_in_tray.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Langerman.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Levy13.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_wctime_Levy.pdf}
    }
    \hfill
    \centering
    \subfigure[Iteration]{
        \includegraphics[width=0.02\linewidth]{manuscripts/fig/regret_axis_label.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Cross_in_tray.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Langerman.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Levy13.pdf}
        \includegraphics[width=0.24\linewidth]{manuscripts/fig/Results_Levy.pdf}
    }

    \caption{Average and standard error of the proposed methods (HB-EI and HB-UCB) and the state-of-the-art preferential BO methods. The horizontal axis represents (a) the computational time (sec) and (b) the number of iterations.
    The vertical axis represents the regret, which is the smaller, the better it is.}
    \label{fig:app-regret}
\end{figure}