# Scale symmetry and composition of compact star matter

Long-Qi Shao  College of Physics, Jilin University, Changchun 130012, China
School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute
for Advanced Study, UCAS, Hangzhou, 310024, China  Yong-Liang Ma
<EMAIL_ADDRESS>School of Fundamental Physics and Mathematical Sciences,
Hangzhou Institute for Advanced Study, UCAS, Hangzhou, 310024, China  College
of Physics, Jilin University, Changchun 130012, China

###### Abstract

The dense compact star matter is studied by using the skyrmion crystal
approach. The chiral effective theory used includes the lightest scalar meson,
the lowest-lying vector mesons as well as pions. Consistency with the vector
manifestation and the dilaton limit fixed point at high density constrains the
anomalous dimension of the gluon field $1.0\lesssim|\gamma_{G^{2}}|\lesssim
2.0$ and leads to the significance of the scale symmetry breaking in the
intrinsic parity-odd part of the effective theory. The speed of sound
$v_{s}^{2}\simeq 1/3$ and the polytropic index $\gamma\simeq 1$—both satisfy
the conformal limits—after the dilaton limit fixed point at high density but
the matter is still in the hadronic phase. This means that neither the
conformal speed of sound nor the smallness of the polytropic index can be used
as a criterion of the onset of quark matter. These conclusions are significant
for constructing the equation of state of nuclear matter.

##  I introduction

Dense nuclear matter relevant to the core of compact star is an uncharted
domain although it has been studied for several decades Brown and Rho (2002);
Li  _et al._ (2008); Holt  _et al._ (2016); Drews and Weise (2017); Li  _et
al._ (2019); Ma and Rho (2020a); Kim  _et al._ (2020); Ma and Rho (2021).
Several fundamental but significant questions in dense nuclear matter, such as
the pattern of symmetry, the constituents and so on are not yet well
understood. The recent observations of massive neutron stars with about two
times solar mass and detections of the gravitational waves from binary neutron
star mergers and neutron star-black hole mergers open new laboratories for
dense nuclear matter (see, e.g., Li  _et al._ (2019); Ma and Rho (2020a); Kim
_et al._ (2020); Ma and Rho (2021) and references therein).

The access to dense nuclear matter is mainly based on the models or effective
theories anchored on the symmetries of QCD. The parameters are fixed by using
the nuclear matter properties around saturation density. The properties of
dense nuclear matter are predicted by extrapolating such obtained equation of
state (EoS) (e.g., Refs. Drews and Weise (2017); Furnstahl (2020)). Whether
this extrapolation is reasonable is questionable. An alternative access is to
implement the conjectured symmetry/constituent to the construction of EoS
(e.g. Alford and Sedrakian (2017)). This construction is highly model
dependent. Therefore, it is valuable to obtain some information of dense
nuclear matter using a unified model.

It is widely accepted that when considered at the large $N_{c}$ limit, baryons
can be described as skyrmions Skyrme (1961); Witten (1979) and consequently
skyrmion matter can be regarded as nuclear matter Klebanov (1985). The
skyrmion approach provides a model independent way to study nucleon, nuclei
and nuclear matter in a unified model in the sense of large $N_{c}$ limit
Zahed and Brown (1986); Naya and Sutcliffe (2018); Ma and Rho (2017). As a
result, it is interesting to use the skyrmion approach to extract some model
independent information of dense nuclear matter where the large $N_{c}$
argument is applicable McLerran and Pisarski (2007); McLerran and Reddy
(2019); Zhao and Lattimer (2020).

Since the isoscalar scalar meson sigma and isoscalar vector meson omega are
indispensable ingredients for the nuclear force Zel’dovich (1961); Serot and
Walecka (1986), we include both in the calculation. The scalar meson is
introduced as the dilaton $\chi$, the Nambu-Goldstone boson of scale symmetry
breaking Crewther and Tunstall (2015); Crewther (2020). The omega meson and
its flavor partners rho mesons are included through the hidden local flavor
symmetry Bando  _et al._ (1985, 1988); Harada and Yamawaki (2003). Both
symmetries are hidden in the matter free space and are expected to emerge in
the (super-)highly condensed matter. In summary, the effective theory (denoted
as $s$HLS) used in this work includes the lightest scalar meson, omega, rho
and the pseudo-Nambu Goldstone bosons pions.

Based on the Wilsonian renormalization group (RG) approach, it is found that
at high energy scale $f_{\pi}^{\ast}\to 0,m_{\rho}^{\ast}\to m_{\pi}\to 0$,
i.e., there is a vector manifestation (VM) fixed point in the hidden local
symmetry (HLS) Harada and Yamawaki (2001, 2003). It is not strange to expect
that the VM appears at (super-)high density. In addition, in the approach to
the baryonic matter using the dilaton compensated chiral effective theory,
people found that there is a dilaton limit (DL) fixed point which states that
the medium modified decay constant $f_{\sigma}^{\ast}\to 0$ in the theory at
high density Beane and van Kolck (1994); Paeng  _et al._ (2012). We expect
both the VM fixed point and the DL fixed point to be realized in the skyrmion
crystal approach.

The purpose of this work is to study how the scale symmetry is realized in
dense nuclear matter and what is the possible constituents in the cores of the
massive neutron stars using skyrmion crystal approach based on $s$HLS. The
$s$HLS to be specified later includes the leading order of the chiral-scale
counting, the intrinsic parity-odd terms as well as a dilaton potential which
breaks the scale symmetry both explicitly and spontaneously. Although the
homogeneous Wess-Zumino (hWZ) term, the intrinsic parity-odd term is scale
invariant, we couple dilaton to it to obtain the finite energy of the skyrmion
matter Park  _et al._ (2008). Therefore, the dilaton compensated hWZ term does
not conserve the scale symmetry. Different from Park  _et al._ (2008) where
the dilaton couples to the hWZ term in an arbitrary way, the coupling between
dilaton and hWZ term has a solid physical argument in the genuine dilaton
approach Crewther and Tunstall (2015); Crewther (2020).

We found that the requirement $2\lesssim n_{\rm DLFP}/n_{0}\lesssim 10$ with
$n_{0}$ being the saturation density of nuclear matter and $n_{\rm DLFP}$
being the onset of the DL fixed point (the same as the VM fixed point in the
present approach) yields the magnitude of the anomalous dimension of gluon
field $1.0\lesssim|\gamma_{G^{2}}|\lesssim 2.0$ which is consistent with Park
_et al._ (2008); Ma and Rho (2018). More interestingly, we found that after
$n_{\rm DLFP}$, the speed of sound $v_{s}^{2}\to 1/3$ and the polytropic index
$\gamma\to 1$ at density relevant to the cores of massive neutron stars, both
satisfy the constraints from conformal invariance. This indicates that the
speed of sound of the compact star matter could be $v_{s}^{2}\approx 1/3$ Ma
and Rho (2019a); Zhao and Lattimer (2020); Kapusta and Welle (2021); Margueron
_et al._ (2021), in contrast to the usual conclusion in the literature Bedaque
and Steiner (2015); Tews  _et al._ (2018); Moustakidis  _et al._ (2017). Note
that the conformal speed of sound of the dense matter after $n_{\rm DLFP}$
does not mean that the matter is scale invariant since the trace of the
energy-momentum tensor (TEMT) is not zero but a density independent constant.
Therefore, the matter after $n_{\rm DLFP}$ is pseudoconformal Paeng  _et al._
(2017); Ma and Rho (2019b); Ma  _et al._ (2019); Ma and Rho (2020a). In
contrast to the usual picture of dense nuclear matter in which the conformal
speed of sound $v_{s}^{2}=1$ and polytropic index $\gamma\to 1$ are regarded
as the criterions of the emergence of the deconfined quark Annala  _et al._
(2020), there is no deconfined quark in the core of massive stars in this
pseudoconformal nuclear matter Ma and Rho (2020b).

We organize this paper as follows: In Sec. II we present the key points of the
$s$HLS used in this work. In Sec. III we discuss the nuclear matter properties
calculated from the skyrmion crystal in detail. Our conclusions and
discussions are presented in Sec. IV.

##  II Hidden scale and Hidden local symmetric Lagrangian

The Skyrme type model we are going to consider in the following is based on
the hidden local symmetry (HLS) Bando  _et al._ (1985, 1988); Harada and
Yamawaki (2003). Among a variety of approaches to the model of the lightest
scalar meson Crewther and Tunstall (2015); Golterman and Shamir (2016);
Crewther (2020); Del Debbio and Zwicky (2021), we here use the genuine dilaton
approach Crewther and Tunstall (2015); Crewther (2020) in a minimal way such
that both the chiral and scale symmetries can (partially) restore in matter.

The dilaton compensated HLS (denoted as $s$HLS) can be written as Crewther and
Tunstall (2015); Li  _et al._ (2017)

$\displaystyle\mathcal{L}_{s\mathrm{HLS}}$ $\displaystyle=$
$\displaystyle\mathcal{L}_{(2)}^{s\mathrm{HLS}}+\tilde{\mathcal{L}}_{\mathrm{hWZ}}+\mathcal{V}(\chi).$
(1)

The leading chiral-scale counting order Lagrangian
$\mathcal{L}_{(2)}^{s\mathrm{HLS}}$ is

$\displaystyle\mathcal{L}_{(2)}^{s\mathrm{HLS}}$ $\displaystyle=$
$\displaystyle f_{\pi}^{2}\left(\frac{\chi}{f_{\sigma}}\right)^{2}{\rm
Tr}\left[\hat{a}_{\perp\mu}\hat{a}_{\perp}^{\mu}\right]+af_{\pi}^{2}\left(\frac{\chi}{f_{\sigma}}\right)^{2}{\rm
Tr}\left[\hat{a}_{\|\mu}\hat{a}_{\|}^{\mu}\right]$ (2)
$\displaystyle{}-\frac{1}{2g^{2}}\operatorname{Tr}\left[V_{\mu\nu}V^{\mu\nu}\right]+\frac{1}{2}\partial_{\mu}\chi\partial^{\mu}\chi$

where the pion fields are $U(x)=\exp(2i\pi/f_{\pi})=\xi_{L}^{\dagger}\xi_{R}$
and

$\displaystyle\hat{a}_{\perp\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2i}\left(D_{\mu}\xi_{R}\cdot\xi_{R}^{\dagger}-D_{\mu}\xi_{L}\cdot\xi_{L}^{\dagger}\right),$
$\displaystyle\hat{a}_{\parallel\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2i}\left(D_{\mu}\xi_{R}\cdot\xi_{R}^{\dagger}+D_{\mu}\xi_{L}\cdot\xi_{L}^{\dagger}\right).$
(3)

The covariant derivative
$D_{\mu}\xi_{L,R}=\partial_{\mu}\xi_{L,R}-iV_{\mu}\xi_{L,R}$ with $V_{\mu}$
being the gauge field of the HLS. In this work, we take hidden local symmetry
as $H_{\rm local}=U(2)$ and identify

$\displaystyle V_{\mu}=g\left(\rho_{\mu}+\omega_{\mu}\right),$ (4)

where $g$ is the gauge coupling constant of HLS. In Eq. (2), the conformal
compensator field $\chi(x)=f_{\sigma}e^{\sigma/f_{\sigma}}$ with $\sigma$
being the dilaton field and $f_{\sigma}$ being the decay constant of dilaton.
Through out this work, we take the chiral limit.

The potential $\mathcal{V}(x)$ in our approach is responsible for both the
spontaneous and explicit breaking of scale symmetry. After taking the saddle
point equation we can write the dilaton potential in terms of dilaton mass as

$\displaystyle
V_{\chi}={}-\frac{f_{\sigma}^{2}m_{\sigma}^{2}}{4\beta^{\prime}}\left(\frac{\chi}{f_{\sigma}}\right)^{4}+\frac{f_{\sigma}^{2}m_{\sigma}^{2}}{\beta^{\prime}\left(4+\beta^{\prime}\right)}\left(\frac{\chi}{f_{\sigma}}\right)^{4+\beta^{\prime}},$
(5)

where $\beta^{\prime}=|\gamma_{G^{2}}|$ is the nonperturbative anomalous
dimension of gluon field which cannot be calculated from QCD. Note that the
dilaton potential (5) is different from the logarithm form used in Ma and Rho
(2018). The logarithm form in Ma and Rho (2018) can be reduced to from (5) by
taking $\beta^{\prime}\ll 1$ so than Eq. (5) is a more general form of the
potential. We will find that the general expression of the dilaton potential
(5) leads to a more stringent constraint on the nonperturbative anomalous
dimension of gluon field.

How the scale symmetry breaking effect is involved in the intrinsic parity-odd
part $\mathcal{L}_{\rm hWZ}$ is a subtle question. Here, for simplicity and
without lose of generality, we take a universal scale symmetry breaking form

$\displaystyle\tilde{\mathcal{L}}_{\mathrm{hWZ}}=\left(c_{h}+\left(1-c_{h}\right)\left(\frac{\chi}{f_{\sigma}}\right)^{\beta^{\prime}}\right)\mathcal{L}_{\mathrm{hWZ}},$
(6)

where $c_{h}$ is a free parameter. When $c_{h}=1$ or
$\langle\chi\rangle=f_{\sigma}$, the dilaton compensator disappears and
$\tilde{\mathcal{L}}_{\rm hWZ}$ reduces to the scale invariant form
$\mathcal{L}_{\rm hWZ}$. The hWZ term $\mathcal{L}_{\rm hWZ}$ takes the form

$\displaystyle\mathcal{L}_{\mathrm{hWZ}}$ $\displaystyle=$
$\displaystyle\frac{N_{c}}{16\pi^{2}}\sum_{i=1}^{3}c_{i}\mathcal{L}_{i}$ (7)

where in terms of the 1-form and 2-form notations

$\displaystyle\mathcal{L}_{1}$ $\displaystyle=$ $\displaystyle
i\operatorname{Tr}\left[\hat{\alpha}_{\mathrm{L}}^{3}\hat{\alpha}_{\mathrm{R}}-\hat{\alpha}_{\mathrm{R}}^{3}\hat{\alpha}_{\mathrm{L}}\right],$
(8a) $\displaystyle\mathcal{L}_{2}$ $\displaystyle=$ $\displaystyle
i\operatorname{Tr}\left[\hat{\alpha}_{\mathrm{L}}\hat{\alpha}_{\mathrm{R}}\hat{\alpha}_{\mathrm{L}}\hat{\alpha}_{\mathrm{R}}\right],$
(8b) $\displaystyle\mathcal{L}_{3}$ $\displaystyle=$
$\displaystyle\operatorname{Tr}\left[F_{\mathrm{V}}\left(\hat{\alpha}_{\mathrm{L}}\hat{\alpha}_{\mathrm{R}}-\hat{\alpha}_{\mathrm{R}}\hat{\alpha}_{\mathrm{L}}\right)\right],$
(8c)

with

$\displaystyle\hat{\alpha}_{L}$
$\displaystyle=\hat{\alpha}_{\|}-\hat{\alpha}_{\perp},$
$\displaystyle\hat{\alpha}_{R}$
$\displaystyle=\hat{\alpha}_{\|}+\hat{\alpha}_{\perp},$ $\displaystyle F_{V}$
$\displaystyle=dV-iV^{2}.$ (9)

Since in the medium $\left\langle\chi\right\rangle^{\ast}/f_{\sigma}<1$, when
$c_{h}\neq 1$ and $\beta^{\prime}\neq 0$, the scale symmetry breaking is
included in the intrinsic parity-odd part. In this note, we take $c_{h}$ as
well as the anomalous dimension $\beta^{\prime}$ as unknown parameters and fix
them by requiring the consistency to the VM and DL fixed points. It should be
noted that when $c_{h}=0$ and $\beta=3$, a specific choice of the low energy
constants in HLS, one reduces to the model of Ref. Park  _et al._ (2008).

##  III  Structure of compact star matter

Now we are ready to calculate the nuclear matter by using the skyrmion crystal
approach. This is accessed by putting the skyrmions calculated from Lagrangian
(1) to a specific crystal lattice Klebanov (1985). Here we take the face-
centered cubic (FCC) lattice which is known yields the lowest ground state
energy so far Kugler and Shtrikman (1988).

In the numerical calculation, we take the typical values $m_{\sigma}=640$ MeV
and $f_{\sigma}=240$ MeV and run $\beta^{\prime}\in(0,4)$ and $c_{h}\in(0,1)$.
The low energy constants $c_{i}$ cannot be determined phenomenologically so
far. Here we choose $c_{1}=-c_{2}=2/3$ and $c_{3}=0$ such that

$\displaystyle\tilde{\mathcal{L}}_{\mathrm{hWZ}}=g\omega_{\mu}B^{\mu}\left(c_{h}+\left(1-c_{h}\right)\left(\frac{\chi}{f_{\sigma}}\right)^{\beta^{\prime}}\right),$
(10)

with $B_{\mu}$ being the topology (baryon) current—the minimal model Meissner
_et al._ (1987). To illustrate the parameter independence of the existence of
DL fixed point, we compare the results with another set of values
$c_{1}=2/3,c_{2}=-1/3$ and $c_{3}=1$ which are fixed by the possible Seiberg
duality that identifies the $\omega$ field with the Chern-Simons field
together with the consistence of vector meson dominance Karasik (2021).

###  III.1 Dilaton limit fixed point and vector manifestation

We first study the constraints on $c_{h}$ and $\beta^{\prime}$ by requiring
the consistency to the VM/DL fixed point, typically, at density $2\lesssim
n_{\rm DLFP}/n_{0}\lesssim 10$. From Lagrangian (1), one can obtain

$\displaystyle f_{\pi}^{\ast}\propto m_{\rho}^{\ast}\propto
f_{\sigma}^{\ast}.$ (11)

Therefore, the VM fixed point is locked to the DL fixed point. Note that since
the general dilaton potential (5) is different from that used in Ma and Rho
(2018) which is valid only $\beta^{\prime}\ll 1$, the constraints are found to
be $0\lesssim c_{h}\lesssim 0.2$ and $1.0\lesssim\beta^{\prime}\lesssim 2.0$,
consistent to but slightly different from that obtained in Ma and Rho (2018).

Figure 1: Density dependence of the medium modified dilaton decay constant
$f_{\sigma}^{\ast}$.

The density dependence of $f_{\sigma}^{\ast}/f_{\sigma}$ is shown in Fig. 1
with some typical values of $c_{h}$ and $\beta^{\prime}$. One can easily see
that, for a fixed $\beta^{\prime}$, the bigger $c_{h}$, the larger $n_{\rm
DLFP}$. This is because, the bigger $c_{h}$, the less omega repulsive force is
affected by density and the less omega force is screened therefore the more
difficult $n_{\rm DLFP}$ can be arrived. When $c_{h}$ is bigger than a certain
typical value—$\approx 0.2$—the omega force is not screened enough so that the
scale restoration cannot happen. In addition, for a fixed $c_{h}$, the smaller
$\beta^{\prime}$ the larger $n_{\rm DLFP}$ because of the larger repulsive
force from omega meson. When $\beta^{\prime}$ is smaller than a typical
value—$\approx 0.5$—the DL/VM fixed point cannot appear and when
$\beta^{\prime}$ exceeds a typical value—$\approx 2.0$—the DL/VM fixed point
appears at a low density that inconsistent with nature.

The result from both $c_{h}\neq 1$ and $\beta^{\prime}\neq 0$ implies the
significance of the scale symmetry breaking in the hWZ part Park  _et al._
(2008); Ma and Rho (2018). Moreover, it tells us that in between the low
density—which cannot be accessed in the present approach—and the DL fixed
point, the symmetry is not visible, hence hidden Ma and Rho (2020c).

###  III.2 Pseudocomformal dense nuclear matter

In the skyrmion crystal approach, the per-skyrmion energy $E/A$ in the medium
can be calculated at a fixed crystal size, or equivalently, density Klebanov
(1985). Interestingly we found that when $n\gtrsim n_{\rm DLFP}$, $E/A$ can be
nicely fitted by the following function

$\displaystyle
E/A=a\left(\frac{n}{n_{0}}\right)^{1/3}+b\left(\frac{n}{n_{0}}\right)^{-1}.$
(12)

The fitted curves compared to the crystal data points are shown in Fig. 2.

Figure 2: Comparison of the crystal data and fitted curves.

From function (12) one can easily conclude that the speed of sound satisfies
the conformal limit

$\displaystyle v_{s}^{2}$ $\displaystyle=$ $\displaystyle\frac{\partial
P(n)}{\partial n}/\frac{\partial\epsilon(n)}{\partial n}\to 1/3,$ (13)

where $P(n)$ and $\epsilon(n)$ are, respectively, the pressure density and
energy density. In addition, one finds that the polytropic index $\gamma$ in
the nuclear matter approaches to

$\displaystyle\gamma=\frac{d\ln P}{d\ln\epsilon}\to 1,$ (14)

when $n\gtrsim n_{\rm DLFP}$. It also satisfies the constraint from conformal
invariance.

Figure 3: Trace of energy-momentum tensor as a function of density.

How to understand $v_{s}^{2}\to 1/3$ in the dense matter? Does it mean that
the matter becomes a conformal invariant theory at high density? To have a
deeper understanding of the structure of the dense nuclear matter, let us
check the density dependence of the trace of the energy-momentum tensor (TEMT)
using the fitting function (12). Our result is plotted in Fig. 3. From this
figure, one can easily see that the TEMT becomes a density independent nonzero
constant after $n_{\rm DLFP}$. This means that, the conformal symmetry is NOT
restored after $n_{\rm DLFP}$. Using the relation between TEMT
$\theta_{\mu}^{\mu}$ and the speed of sound

$\displaystyle\frac{\partial\theta_{\mu}^{\mu}}{\partial
n}=\frac{\partial\epsilon(n)}{\partial n}(1-3v_{s}^{2}),$ (15)

one can see that, since $\theta_{\mu}^{\mu}$ is independent of density, the
left hand side of this equation is zero. Because $\epsilon(n)$ varies with
density, $\frac{\partial\epsilon(n)}{\partial n}\neq 0$, one has
$1-3v_{s}^{2}=0$, i.e., $v_{s}^{2}=1/3$. Since the TEMT is not zero after
$n_{\rm DLFP}$ but the speed of sound satisfies the conformal limit, the
matter is pseudoconformal. Therefore, one can conclude neither the emergent of
the conformal symmetry nor an onset of the deconfined quark matter from the
conformal speed of sound and unit of polytropic index.

###  III.3 Minimal model vs Seiberg duality

What discussed above are based on the minimal model with parameter choice
$c_{1}=-c_{2}=2/3$ and $c_{3}=0$ Meissner  _et al._ (1987). So far, no
confidential phenomena can uniquely fix these parameters. Here, to show that
the existence of DL fixed point is independent of the choice of $c_{i}$, we
choose another set of parameter $c_{1}=2/3,c_{2}=-1/3$ and $c_{3}=1$ Karasik
(2021)(Karasik model) which is consistent with the vector meson dominance and
the vector meson-gluon duality that realizes the Seiberg duality in the strong
interaction sector Seiberg (1995).

Figure 4: Medium modified dilaton decay constant $f_{\sigma}^{\ast}$ from the
minimal model (solid) and Karasik model (dashed).

We plot $f_{\sigma}^{\ast}$ from the minimal model and the Karasik model in
Fig. 4 for $c_{h}=0.2$ and $\beta^{\prime}=1.0$. One can see that, the value
of $n_{\rm DLFP}$ is smaller in the Karasik model than that in the minimal
model. This is because, compared to the minimal model, there is an extra
attractive interaction from rho meson in the Karasik model. Effectively, the
role of this extra rho meson interaction is identical to increase
$\beta^{\prime}$ in the minimal model.

More importantly, this comparison indicates that, in a dilaton compensated
Skyrme type model with a proper choice of the parameters, the existence of DL
and VM fixed points is a genuine conclusion.

###  III.4 The cores of massive stars

We finally confront the EoS calculated in the present work with that
constrained from astrophysics including the gravitational wave detection. The
purpose is to simply give a qualitative idea since only in the large $N_{c}$
limit a skyrmion can be regarded as a baryon. The constraint (green area) is
from Annala  _et al._ (2018) with upper bound of neutron star $\gtrsim
2.0M_{\odot}$. Since the EoS from the crystal approach has negative pressure
at $n\lesssim 5.7n_{0}$, we only consider the EoS in the region
$5.7n_{0}\lesssim n\lesssim 10n_{0}$.

Figure 5: Comparison of the EoS calculated in this work upto $10n_{0}$ and
that from the constraint from astrophysics.

The comparison of the presently calculated EoS with the constraint is shown in
Fig. 5. We find that at density $n\gtrsim 7.5n_{0}>n_{\rm DLFP}$ the EoS falls
into the region of the constraint. This is not difficult to understand if one
respects the fact that a baryon can be regarded as a skyrmion in the sense of
large $N_{c}$ limit and the nuclear matter at low density is not a crystal
solid therefore the present numerical results make sense at large density.

##  IV Conclusion and Discussion

We study in this paper the properties of dense nuclear matter using the
skyrmion crystal approach based on $s$HLS including the lightest scalar meson
$\sigma$, isovector vector meson $\rho$, isoscalar vector meson $\omega$, in
addition to the pseudo-Nambu-Goldstone boson $\pi$.

We found that the requirement of the scale symmetry restoration in dense
nuclear matter, a natural expectation from physics, constrains the magnitudes
of the scale symmetry breaking as $0\lesssim c_{h}\lesssim 0.2$ and
$1.0\lesssim\beta^{\prime}\lesssim 2.0$. We also found that after $n_{\rm
DLFP}$, both the speed of sound and the polytropic index approach to the
conformal limit $v_{s}^{2}\to 1/3$ and $\gamma\to 1$ while the TEMT does not
vanish but becomes a density independent constant. Therefore, the dense
nuclear matter is a pseudoconformal one. To have a qualitative idea, we
compare the EoS calculated in the present approach and the constraint from
heavy ion physic and astrophysics.

Since the TEMT is a nonzero constant, our matter is still in the hadronic
phase after DL fixed point, not in the quark phase. This means that, in
contrast to Annala  _et al._ (2020), in nuclear matter, neither the conformal
speed of sound nor the smallness of the polytropic index can be used as a
sufficient criterion of the onset of the quark matter. The constant TEMT after
DL fixed point may closely related to chiral invariant mass of nucleon Detar
and Kunihiro (1989) which has been found significant in constructing the EoS
of nuclear matter Motohiro  _et al._ (2015). The chiral invariant mass of
nucleon may attribute to the multiquark condensate and gluon condensate in the
present approach.

###### Acknowledgements.

Y. L. M would like to thank Mannque Rho for his valuable discussions and
comments. The work of Y. L. M. was supported in part by National Science
Foundation of China (NSFC) under Grant No. 11875147 and No. 12147103.

## References

  * Brown and Rho (2002) G. E. Brown and M. Rho, Phys. Rept. 363, 85 (2002), arXiv:hep-ph/0103102 .
  * Li  _et al._ (2008) B.-A. Li, L.-W. Chen, and C. M. Ko, Phys. Rept. 464, 113 (2008), arXiv:0804.3580 [nucl-th] .
  * Holt  _et al._ (2016) J. W. Holt, M. Rho, and W. Weise, Phys. Rept. 621, 2 (2016), arXiv:1411.6681 [nucl-th] .
  * Drews and Weise (2017) M. Drews and W. Weise, Prog. Part. Nucl. Phys. 93, 69 (2017), arXiv:1610.07568 [nucl-th] .
  * Li  _et al._ (2019) B.-A. Li, P. G. Krastev, D.-H. Wen, and N.-B. Zhang, Eur. Phys. J. A 55, 117 (2019), arXiv:1905.13175 [nucl-th] .
  * Ma and Rho (2020a) Y.-L. Ma and M. Rho, Prog. Part. Nucl. Phys. 113, 103791 (2020a), arXiv:1909.05889 [nucl-th] .
  * Kim  _et al._ (2020) M. Kim, C.-H. Lee, Y.-M. Kim, K. Kwak, Y. Lim, and C. H. Hyun, Int. J. Mod. Phys. E 29, 2030007 (2020).
  * Ma and Rho (2021) Y.-L. Ma and M. Rho, Symmetry 13, 1888 (2021), arXiv:2009.09219 [nucl-th] .
  * Furnstahl (2020) R. J. Furnstahl, Eur. Phys. J. A 56, 85 (2020), arXiv:1906.00833 [nucl-th] .
  * Alford and Sedrakian (2017) M. G. Alford and A. Sedrakian, Phys. Rev. Lett. 119, 161104 (2017), arXiv:1706.01592 [astro-ph.HE] .
  * Skyrme (1961) T. H. R. Skyrme, Proc. Roy. Soc. Lond. A 260, 127 (1961).
  * Witten (1979) E. Witten, Nucl. Phys. B 160, 57 (1979).
  * Klebanov (1985) I. R. Klebanov, Nucl. Phys. B 262, 133 (1985).
  * Zahed and Brown (1986) I. Zahed and G. E. Brown, Phys. Rept. 142, 1 (1986).
  * Naya and Sutcliffe (2018) C. Naya and P. Sutcliffe, Phys. Rev. Lett. 121, 232002 (2018), arXiv:1811.02064 [hep-th] .
  * Ma and Rho (2017) Y.-L. Ma and M. Rho, Sci. China Phys. Mech. Astron. 60, 032001 (2017), arXiv:1612.06600 [nucl-th] .
  * McLerran and Pisarski (2007) L. McLerran and R. D. Pisarski, Nucl. Phys. A 796, 83 (2007), arXiv:0706.2191 [hep-ph] .
  * McLerran and Reddy (2019) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019), arXiv:1811.12503 [nucl-th] .
  * Zhao and Lattimer (2020) T. Zhao and J. M. Lattimer, Phys. Rev. D 102, 023021 (2020), arXiv:2004.08293 [astro-ph.HE] .
  * Zel’dovich (1961) Y. B. Zel’dovich, Zh. Eksp. Teor. Fiz. 41, 1609 (1961).
  * Serot and Walecka (1986) B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16, 1 (1986).
  * Crewther and Tunstall (2015) R. J. Crewther and L. C. Tunstall, Phys. Rev. D 91, 034016 (2015), arXiv:1312.3319 [hep-ph] .
  * Crewther (2020) R. J. Crewther, Universe 6, 96 (2020), arXiv:2003.11259 [hep-ph] .
  * Bando  _et al._ (1985) M. Bando, T. Kugo, S. Uehara, K. Yamawaki, and T. Yanagida, Phys. Rev. Lett. 54, 1215 (1985).
  * Bando  _et al._ (1988) M. Bando, T. Kugo, and K. Yamawaki, Phys. Rept. 164, 217 (1988).
  * Harada and Yamawaki (2003) M. Harada and K. Yamawaki, Phys. Rept. 381, 1 (2003), arXiv:hep-ph/0302103 .
  * Harada and Yamawaki (2001) M. Harada and K. Yamawaki, Phys. Rev. Lett. 86, 757 (2001), arXiv:hep-ph/0010207 .
  * Beane and van Kolck (1994) S. R. Beane and U. van Kolck, Phys. Lett. B 328, 137 (1994), arXiv:hep-ph/9401218 .
  * Paeng  _et al._ (2012) W.-G. Paeng, H. K. Lee, M. Rho, and C. Sasaki, Phys. Rev. D 85, 054022 (2012), arXiv:1109.5431 [hep-ph] .
  * Park  _et al._ (2008) B.-Y. Park, M. Rho, and V. Vento, Nucl. Phys. A 807, 28 (2008), arXiv:0801.1374 [hep-ph] .
  * Ma and Rho (2018) Y.-L. Ma and M. Rho, Phys. Rev. D 97, 094017 (2018), arXiv:1612.04079 [nucl-th] .
  * Ma and Rho (2019a) Y.-L. Ma and M. Rho, Phys. Rev. D 100, 114003 (2019a), arXiv:1811.07071 [nucl-th] .
  * Kapusta and Welle (2021) J. I. Kapusta and T. Welle, Phys. Rev. C 104, L012801 (2021), arXiv:2103.16633 [nucl-th] .
  * Margueron  _et al._ (2021) J. Margueron, H. Hansen, P. Proust, and G. Chanfray, Phys. Rev. C 104, 055803 (2021), arXiv:2103.10209 [nucl-th] .
  * Bedaque and Steiner (2015) P. Bedaque and A. W. Steiner, Phys. Rev. Lett. 114, 031103 (2015), arXiv:1408.5116 [nucl-th] .
  * Tews  _et al._ (2018) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018), arXiv:1801.01923 [nucl-th] .
  * Moustakidis  _et al._ (2017) C. C. Moustakidis, T. Gaitanos, C. Margaritis, and G. A. Lalazissis, Phys. Rev. C 95, 045801 (2017), [Erratum: Phys.Rev.C 95, 059904 (2017)], arXiv:1608.00344 [nucl-th] .
  * Paeng  _et al._ (2017) W.-G. Paeng, T. T. S. Kuo, H. K. Lee, Y.-L. Ma, and M. Rho, Phys. Rev. D 96, 014031 (2017), arXiv:1704.02775 [nucl-th] .
  * Ma and Rho (2019b) Y.-L. Ma and M. Rho, Phys. Rev. D 99, 014034 (2019b), arXiv:1810.06062 [nucl-th] .
  * Ma  _et al._ (2019) Y.-L. Ma, H. K. Lee, W.-G. Paeng, and M. Rho, Sci. China Phys. Mech. Astron. 62, 112011 (2019), arXiv:1804.00305 [nucl-th] .
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Ma and Rho (2020b) Y.-L. Ma and M. Rho, (2020b), arXiv:2006.14173 [nucl-th] .
  * Golterman and Shamir (2016) M. Golterman and Y. Shamir, Phys. Rev. D 94, 054502 (2016), arXiv:1603.04575 [hep-ph] .
  * Del Debbio and Zwicky (2021) L. Del Debbio and R. Zwicky, (2021), arXiv:2112.11363 [hep-ph] .
  * Li  _et al._ (2017) Y.-L. Li, Y.-L. Ma, and M. Rho, Phys. Rev. D 95, 114011 (2017), arXiv:1609.07014 [hep-ph] .
  * Kugler and Shtrikman (1988) M. Kugler and S. Shtrikman, Physics Letters B 208, 491 (1988).
  * Meissner  _et al._ (1987) U. G. Meissner, N. Kaiser, and W. Weise, Nucl. Phys. A 466, 685 (1987).
  * Karasik (2021) A. Karasik, SciPost Phys. 10, 138 (2021), arXiv:2010.10544 [hep-th] .
  * Ma and Rho (2020c) Y.-L. Ma and M. Rho, Phys. Rev. Lett. 125, 142501 (2020c), arXiv:2002.03310 [nucl-th] .
  * Seiberg (1995) N. Seiberg, Nucl. Phys. B 435, 129 (1995), arXiv:hep-th/9411149 .
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Detar and Kunihiro (1989) C. E. Detar and T. Kunihiro, Phys. Rev. D 39, 2805 (1989).
  * Motohiro  _et al._ (2015) Y. Motohiro, Y. Kim, and M. Harada, Phys. Rev. C 92, 025201 (2015), [Erratum: Phys.Rev.C 95, 059903 (2017)], arXiv:1505.00988 [nucl-th] .


# Optimal designs for the development of personalized treatment rules

David Azriela, Yosef Rinottb and Martin Poschc Corresponding author: Martin
Posch; email<EMAIL_ADDRESS>
###### Abstract

In the present paper, personalized treatment means choosing the best treatment
for a patient while taking into account certain relevant personal covariate
values. We study the design of trials whose goal is to find the best treatment
for a given patient with given covariates. We assume that the subjects in the
trial represent a random sample from the population, and consider the
allocation, possibly with randomization, of these subjects to the different
treatment groups in a way that depends on their covariates.

We derive approximately optimal allocations, aiming to minimize expected
regret, assuming that future patients will arrive from the same population as
the trial subjects. We find that, for the case of two treatments, an
approximately optimal allocation design does not depend on the value of the
covariates but only on the variances of the responses. In contrast, for the
case of three treatments the optimal allocation design does depend on the
covariates as we show for specific scenarios. Another finding is that the
optimal allocation can vary a lot as a function of the sample size, and that
randomized allocations are relevant for relatively small samples, and may not
be needed for very large studies.

a Faculty of Industrial Engineering and Management, The Technion
b Department of Statistics and Federmann Center for the Study of Rationality,
The Hebrew University
c Center for Medical Statistics, Informatics and Intelligent Systems, Medical
University of Vienna

Key words: personalized medicine; minimal regret; optimal allocation;
experimental design.

##  1 Introduction

This work concerns the study of optimal designs of clinical trials with the
goal of estimating personalized treatment rules. We consider multi-armed
clinical trials for treatments, where the primary outcome may depend on
patients’ covariates. Modeling the outcome a function of the covariates does
not only reduce the error variance, but allows one also to estimate
“personalized” treatment rules. These rules assign the treatment with the best
estimated outcome for a patient having given values of the covariates.

In the related context of bandit problems with covariates, Goldenshluger and
Zeevi (2013) studied a sequential allocation scheme where at each stage one
out of two treatments needs to be assigned under a minimax framework. A high-
dimensional version of this problem was studied in Bastani and Bayati (2015).
The goal in the Bandit formulation is to determine the best treatment while
minimizing some loss or regret function for subjects in the trial sample. This
is a different setting than the one considered here, since rather than be
concerned with optimal treatment for subjects in the trial sample, we look for
designs that will allocate patients to treatments efficiently for the goal of
finding the best treatment for future patients who may require one of the
treatments under study. Personalized design as a function of covariates
appears in many papers, such as Qian and Murphy (2011), Tian et al. (2014),
Ballarini et al. (2018), or Pan and Zhao (2020). However, these papers are
concerned with inference about the optimal treatment rather than the design of
the trial, which is the goal of this paper.

Problem statement: Consider $K$ possible treatments $T_{1},\ldots,T_{K}$. Let
$Y$ be a one-dimensional continuous response variable and let ${\bf
X}\in\mathbb{R}^{p}$ be a vector of a subject’s covariates. We assume that the
joint distribution of the covariates is continuous with density denoted by
$f({\bf x})$. The expected response is $E(Y|{\bf X},T_{k})=g_{k}({\bf X})$,
where $g_{k}$ is an unknown function. The optimal treatment for a subject with
covariate ${\bf x}$ is $\delta^{*}({\bf
x})=\arg\max_{k\in\\{1,\ldots,K\\}}g_{k}({\bf x})$ (assuming that a higher
response is better). If the above arg max contains more than one $k$, one
arbitrary treatment is selected.

Suppose that a clinical trial with $n$ subjects is performed in order to
estimate $\delta^{*}$. Let ${\bf X}_{1},\ldots,{\bf X}_{n}$ denote a sample
according to $f$ of the covariates of the $n$ subjects. The design we study
here consists of allocating subjects to treatments, taking account of their
covariates. Each subject $i$ is allocated independently to treatment $k$ with
probability $\pi_{k}({\bf X}_{i})$, where $\pi_{1}({\bf
x}),\ldots,\pi_{K}({\bf x})$ are non-negative functions satisfying for each
${\bf x}$, $\sum_{k=1}^{K}\pi_{k}({\bf x})=1$. The allocation functions define
densities for the covariates in the $K$ treatment groups. Specifically, in
treatment $k$, ${\bf X}$ is sampled from the density $f_{k}({\bf x}):=f({\bf
x})\pi_{k}({\bf x})/\nu_{k}$ where $\nu_{k}:=\int f({\bf x})\pi_{k}({\bf
x})d{\bf x}$. Our purpose is to find the optimal allocation functions
$\pi_{1}({\bf x}),\ldots,\pi_{K}({\bf x})$ that minimize the regret, which is
defined next.

Let $\widehat{\delta}$ denote the final estimate of the optimal rule
$\delta^{*}$ based on the covariates, responses, and allocations of the
subjects in the trial. Then the regret for a randomly chosen future patient is
defined as

$R(\pi_{1},\ldots,\pi_{K}):=E\left[g_{\delta^{*}(\widetilde{\bf
X})}(\widetilde{\bf X})-g_{\widehat{\delta}(\widetilde{\bf X})}(\widetilde{\bf
X})\right],$ (1)

where $\widetilde{\bf X}$ denotes the covariate vector of the chosen patient
and the expectation is over $\widetilde{\bf X}$, $\widehat{\delta}$, and the
randomizations in the allocations. In words, the regret is the difference
between the expected response to the optimal treatment $\delta^{*}$ and its
estimate $\widehat{\delta}$ for any independent future subject arising from
the same population as those in the trial. A similar criterion appears in Qian
and Murphy (2011). Our goal is to minimize the regret over all feasible
allocations, that is, non-negative functions $\pi_{k}(\cdot)$ that satisfy
$\sum_{k=1}^{K}\pi_{k}({\bf x})=1$ for every ${\bf x}$.

In this study we assume that we have some preliminary knowledge about the
parameters of the problem, such as the functions $g_{k}$, and the variances of
$Y$ given ${\bf X},T_{k}$ and the density $f$ of the covariates. This approach
concurs with existing literature on locally optimal designs; see, e.g.,
Chernoff (1953) and Silvey (2013, Chapter 6), where optimality is achieved for
a given set of parameters. Our goal is to find optimal allocation functions
relative to the assumed parameter values, and obtain a treatment rule
$\widehat{\delta}$ based on an experiment according to our design. Knowledge
about the parameters can arise from previous experience or theory, from a
pilot study, or in some situations from earlier phases of the study. Other
approaches, like sequential or Bayesian designs will not be discussed here.
However, our approach suggests a natural adaptive formulation where after
observing an initial part of the sample one estimates the parameters and aims
at an optimal design relative to these estimates, and the process can be
repeated.

A trial of the kind we study ends with recommendations on treatments,
amounting to a claim of causality. Here we allow the allocation to treatments
to be a function of the covariate vector $\bf X$ only, and thus conditioned on
$\bf X$, the allocation and the response are obviously independent. In this
paper we shall assume that a linear model holds true (see Section 2); however,
assuming any tractable parametric model on $E(Y|{\bf X})$ can be handled in
the same way, using non-linear least squares. Assuming a model, linear or more
complex, rather than taking the approach that “all models are wrong”, is
necessary if one wants to avoid a separate trial for each level of $\bf X$,
which seems inefficient. A linear model assumption may be reasonable if the
range of $\bf X$ is limited, which is often very natural. Another possibility
is to consider a polynomial linear model of some degree rather than a simple
linear model (see Section 3.3). The parameters of any model for $E(Y|{\bf X})$
obviously depend only on the conditional distribution of $Y$ given $\bf X$,
allowing us to use them for causal inference given $X$; see, e.g., Hernán and
Robins (2010). This remains true even if the allocation functions are
“deterministic”, taking only the values 0 and 1. On the other hand, suppose we
use the best linear predictors without assuming a linear model. In this case,
the parameters of the best linear predictor (projection coefficients) depend
on the distribution of $X$. In this case, causality cannot be claimed.
Starting in Section 2 we shall assume a linear model.

Main results: We study the regret via an approximation that assumes “ideal”
conditions, such as normality instead of asymptotic normality of linear
regression coefficient estimators. We show that this yields a useful
approximation for quite small sample sizes. We also study an asymptotic
regret. It turns out that when there are three or more treatments, optimal
designs for moderate sample sizes may be totally different from the asymptotic
design, and in this sense one should beware of asymptotics.

For the case of two treatments, we show that the optimal trial design
allocates patients in proportion to the standard deviations of the response
under each treatment, and the optimal allocation probabilities do not depend
on the covariates. This is shown for the case where the outcome in each
treatment group depends on the covariates according to a $p$-dimensional
linear regression model, as well as for the case of a single covariate and a
polynomial regression model.

For the case of three or more treatments and a linear regression model in a
single covariate, we combine theoretical arguments and numerical calculations
to study some examples of optimal designs. This includes a study comparing
three treatments (see Ebbeling et al. (2018)) based on data that are
reconstructed from a published paper; we compare the design in this paper to
an optimal design obtained by our approach. We show that the optimal
allocation in general depends on the covariate. For specific cases, we show
that the asymptotically optimal allocation is deterministic, i.e., the range
of the covariate is partitioned into intervals and in each interval the
allocation probability to a particular treatment is 1. Such allocation rules
arise typically in very large studies.

##  2 The regret under “ideal” conditions

Let $(Y_{1},{\bf X}_{1}),\ldots,(Y_{n},{\bf X}_{n})$ be a sample obtained by
given allocation functions $\pi_{1},\ldots,\pi_{n}$, along with the treatment
allocated for each ${\bf X}_{i}$. Assume $K\geq 2$ and ${\bf
X}\in\mathbb{R}^{p}$, and under treatment $T_{k}$ we have $Y=g_{k}({\bf
X})+\varepsilon$, where $g_{k}({\bf X})=E(Y|{\bf
X},T_{k})=\alpha_{k}+{\boldsymbol{\beta}}_{k}^{t}{\bf X}$, and
$\sigma_{k}^{2}:=Var(\varepsilon|{\bf X},T_{k})$. Let
$\widehat{\alpha}_{k},{\widehat{\boldsymbol{\beta}}_{k}}$ denote the OLS
estimators based on data from treatment $T_{k}$, $\widehat{g}_{k}({\bf X})$
the corresponding estimated regression functions, and $\widehat{\delta}({\bf
x})=\arg\max_{k}\widehat{g}_{k}({\bf X})$ the decision rule. Then,

$\displaystyle R(\pi_{1},\ldots,\pi_{K})$
$\displaystyle=\sum_{k=1}^{K}\int_{\mathbb{R}^{p}}P\left({\widehat{\delta}({\bf
x})}=k\right)\left[g_{\delta^{*}({\bf x})}({\bf x})-g_{k}({\bf
x})\right]f({\bf x})d{\bf x}.$ (2)

For $k=1,\ldots,K$, let

${\boldsymbol{Q}}_{k}:=\int_{\mathbb{R}^{p}}\binom{1}{\bf x}({1},{\bf
x}^{t})f_{k}({\bf x})d{\bf
x},\quad{\boldsymbol{\Sigma}}_{k}:=\frac{1}{\nu_{k}}\sigma^{2}_{k}{\boldsymbol{Q}}_{k}^{-1},\quad\xi_{k}^{2}({\bf
x}):=({1},{\bf x}^{t}){\boldsymbol{\Sigma}}_{k}\binom{1}{\bf x}.$ (3)

It is well known that
$\sqrt{n}[(\widehat{\alpha}_{k},{\widehat{\boldsymbol{\beta}}_{k}})-({\alpha}_{k},{\boldsymbol{\beta}}_{k})]\to
N(0,{\boldsymbol{\Sigma}}_{k})$, and in particular ${\boldsymbol{\Sigma}}_{k}$
is the asymptotic variance of the OLS estimators.

In order to approximate the probability $P\left({\widehat{\delta}({\bf
x})}=k\right)$, we first impose “ideal” conditions, to be dropped later: 1.
The asymptotic variance ${\boldsymbol{\Sigma}}_{k}$ is the exact variance
matrix of the OLS estimators and thus $\xi_{k}^{2}({\bf
x})/n=Var(\widehat{g}_{k}({\bf x}))$; and 2. The error $\varepsilon|{\bf
X},T_{k}\sim N(0,\sigma_{k}^{2})$. Under the ideal conditions the estimators
$\widehat{g}_{k}({\bf x})$ are jointly normal, and independent conditionally
on ${\bf X}_{1},\ldots{\bf X}_{n},T_{1},\ldots T_{K}$, with constat
expectation $\alpha_{k}+{\boldsymbol{\beta}}_{k}^{t}{\bf x}$, and by normality
they are independent. Therefore, with $P_{I}$ and $E_{I}$ denoting the
probability and expectation under the ideal conditions, we have

$\displaystyle P_{I}\left({\widehat{\delta}({\bf
x})}=k\right)=P_{I}\left(\max_{l=1,\ldots,K,l\not=k}\widehat{g}_{l}({\bf
x})<\widehat{g}_{k}({\bf x})\right)$ $\displaystyle=E_{I}\prod_{\ell\neq
k}P_{I}\left(\frac{\sqrt{n}[\widehat{g}_{l}({\bf x})-{g}_{l}({\bf
x})]}{\xi_{\ell}({\bf x})}<\frac{\sqrt{n}[\widehat{g}_{k}({\bf
x})-{g}_{\ell}({\bf x})]}{\xi_{\ell}({\bf x})}\,\Big{|}\,\widehat{g}_{k}({\bf
x})\right)$ $\displaystyle=E_{I}\prod_{\ell\neq
k}P\left(Z\leq\frac{\sqrt{n}[\widehat{g}_{k}({\bf x})-{g}_{\ell}({\bf
x})]}{\xi_{\ell}({\bf x})}\Big{|}\,\widehat{g}_{k}({\bf
x}\right)=E\prod_{\ell\neq k}P\left(\widetilde{Z}\leq\frac{{Z}\xi_{l}({\bf
x})+\sqrt{n}[g_{k}({\bf x})-g_{l}({\bf x})]}{\xi_{k}({\bf x})}\Big{|}Z\right)$
$\displaystyle=\int\prod_{l=1,\ldots,K,l\not=k}\Phi\left(\frac{{z}\xi_{l}({\bf
x})+\sqrt{n}[g_{k}({\bf x})-g_{l}({\bf x})]}{\xi_{k}({\bf
x})}\right)\varphi(z)dz,$ (4)

where the $E$ and $P$ in the penultimate expression are with respect to $Z$
and $\widetilde{Z}$, which are independent $N(0,1)$, and $\widetilde{Z}$
represents $\frac{\sqrt{n}[\widehat{g}_{k}({\bf x})-{g}_{k}({\bf
x})]}{\xi_{k}({\bf x})}$; $\varphi$ and $\Phi$ denote the standard normal
density and cumulative distribution function. The product in the second line
was justified above by independence under the ideal conditions, which are used
also in both equalities of the third line. We define the ideal regret, denoted
by $R_{I}$, by

$R_{I}(\pi_{1},\ldots,\pi_{K}):=\sum_{k=1}^{K}\int_{\mathbb{R}^{p}}P_{I}\left({\widehat{\delta}({\bf
x})}=k\right)\left[g_{\delta^{*}({\bf x})}({\bf x})-g_{k}({\bf
x})\right]f({\bf x})d{\bf x},$ (5)

where $P_{I}\left({\widehat{\delta}({\bf x})}=k\right)$ is given in (2).

##  3 The case of two treatments

###  3.1 Regret approximations and optimal design

We assume that $K=2$ and $Y=g_{k}({\bf X})+\varepsilon$ under treatment
$T_{k}$ for $k=1,2$, where $g_{k}({\bf X})$ denotes the conditional mean of
$Y$ given ${\bf X},T_{k}$, ${\bf X}\in\mathbb{R}^{p}$ and

$g_{1}({\bf X})=\alpha_{1}+{\boldsymbol{\beta}}_{1}^{t}{\bf X},~{}g_{2}({\bf
X})=\alpha_{2}+{\boldsymbol{\beta}}_{2}^{t}{\bf X},~{}Var(\varepsilon|{\bf
X},T_{1})=\sigma_{1}^{2}~{}\text{and}~{}Var(\varepsilon|{\bf
X},T_{2})=\sigma_{2}^{2}.$ (6)

We do not assume normality of the errors. However, we assume the existence of
the moment generating function of $\varepsilon$ when conditioned on ${\bf
X},T_{k}$. This assumption is needed for some large deviation exponential
bounds used below, and could be relaxed to assuming finiteness of some high
order moments instead. We also assume that ${\bf X}$ is continuous with a
bounded density $f({\bf x})$ supported on $[0,1]^{p}$.

Under the ideal conditions of Section 2 with $K=2$, Equation (2) becomes

$P_{I}\left({\widehat{\delta}({\bf
x})}=k\right)=\Phi\left(\frac{\sqrt{n}[g_{k}({\bf x})-g_{l}({\bf
x})]}{\sqrt{V({\bf x})}}\right)\\!\\!,\text{where }V({\bf
x}):={\xi^{2}_{1}({\bf x})+\xi^{2}_{2}({\bf x})}=(1,{\bf
x}^{t})\left({\boldsymbol{\Sigma}}_{1}+{\boldsymbol{\Sigma}}_{2}\right)\binom{1}{{\bf
x}}$ (7)

for $1\leq k\neq l\leq 2$. The ideal regret (5) is

$\displaystyle R_{I}(\pi_{1},\pi_{2})=$ $\displaystyle\int_{\\{{\bf
x}:g_{1}({\bf x})>g_{2}({\bf x})\\}}\Phi\left(\frac{\sqrt{n}[g_{2}({\bf
x})-g_{1}({\bf x})]}{\sqrt{V({\bf x})}}\right)[g_{1}({\bf x})-g_{2}({\bf
x})]f({\bf x})d{\bf x}$ $\displaystyle+\int_{\\{{\bf x}:g_{1}({\bf
x})<g_{2}({\bf x})\\}}\Phi\left(\frac{\sqrt{n}[g_{1}({\bf x})-g_{2}({\bf
x})]}{\sqrt{V({\bf x})}}\right)[g_{2}({\bf x})-g_{1}({\bf x})]f({\bf x})d{\bf
x}$ $\displaystyle=\int_{[0,1]^{p}}\Phi\left(-\frac{\sqrt{n}|g_{1}({\bf
x})-g_{2}({\bf x})|}{\sqrt{V({\bf x})}}\right)|g_{2}({\bf x})-g_{1}({\bf
x})|f({\bf x})d{\bf x}.$ (8)

By (3.1), minimization of $R_{I}(\pi_{1},\pi_{2})$ amounts to finding the
allocation functions $\pi_{1},\pi_{2}$ that minimize $V({\bf x})$. The
following theorem states that this can be done by minimizing $V({\bf x})$
uniformly over ${\bf x}$.

######  Theorem 3.1.

The allocation functions $\pi^{0}_{1}({\bf
x}):=\frac{\sigma_{1}}{\sigma_{1}+\sigma_{2}}$ and $\pi^{0}_{2}({\bf
x}):=\frac{\sigma_{2}}{\sigma_{1}+\sigma_{2}}$ minimize both $V({\bf x})$
uniformly over ${\bf x}$, and $R_{I}(\pi_{1},\pi_{2})$.

Recalling that $V({\bf x})$ can be written as a quadratic form, see (7),
Theorem 3.1 follows immediately from Lemma 3.2 below whose proof, as all other
proofs, is given in the Appendix. Given matrices ${A}$ and ${B}$, we write
${A}\succeq{B}$ if ${A}-{B}$ is positive semi-definite. Let
${\boldsymbol{Q}}:=\int_{[0,1]^{p}}\binom{1}{\bf x}({1},{\bf x})f({\bf
x})d{\bf x}$, and recall the notation of (3). We have

######  Lemma 3.2.

With the above definitions,

${\boldsymbol{\Sigma}}_{1}+{\boldsymbol{\Sigma}}_{2}\succeq(\sigma_{1}+\sigma_{2})^{2}\mathbb{\boldsymbol{Q}}^{-1}.$
(9)

Moreover, the lower bound is attained when $\pi_{1}({\bf x})=\pi^{0}_{1}({\bf
x})$ and $\pi_{2}({\bf x})=\pi^{0}_{2}({\bf x})$.

To see the latter statement note that when $\pi_{k}({\bf x})$ do not depend on
$\bf x$, both matrices ${\boldsymbol{Q}}_{k}$ are proportional to
${\boldsymbol{Q}}$ and therefore both matrices ${\boldsymbol{\Sigma}}_{k}$ are
proportional to ${\boldsymbol{Q}}^{-1}$.

We now discuss the regret (2). In Theorem 3.3 we approximate the regret by the
ideal regret, and then we use the above results for the ideal regret to find
an asymptotically optimal design for the regret itself. We later demonstrate
numerically that the latter design provides a good approximation for finite
sample sizes. See also Theorem 3.3. The regret is given by

$R(\pi_{1},\pi_{2})=\int_{\\{{\bf x}:g_{1}({\bf x})>g_{2}({\bf
x})\\}}P(\widehat{g}_{2}({\bf x})>\widehat{g}_{1}({\bf x}))(g_{1}({\bf
x})-g_{2}({\bf x}))f({\bf x})d{\bf x}\\\ +\int_{\\{{\bf x}:g_{2}({\bf
x})>g_{1}({\bf x})\\}}P(\widehat{g}_{1}({\bf x})>\widehat{g}_{2}({\bf
x}))(g_{2}({\bf x})-g_{1}({\bf x}))f({\bf x})d{\bf x}.$ (10)

######  Theorem 3.3.

Under model (6) with $\nu_{1},\nu_{2}>0$, we have for any $\varepsilon>0$

$n^{3/2-\varepsilon}\left|R(\pi_{1},\pi_{2})-R_{I}(\pi_{1},\pi_{2})\right|\to
0\text{ as }n\to\infty.$

To facilitate the representation, we state the next result in the one-
dimensional case first. In this case $g_{1}$ and $g_{2}$ are functions of a
single variable, and we assume that their intersection point
$\theta:=\frac{\alpha_{1}-\alpha_{2}}{\beta_{2}-\beta_{1}}$ is in $[0,1]$.
Assume without loss of generality, that $\beta_{2}>\beta_{1}$. It follows that
for $x<\theta$, $g_{1}(x)>g_{2}(x)$, that is, $T_{1}$ is the better treatment,
and when $x>\theta$, $T_{2}$ is better. In the case $p=1$ the limit of the
regret is simple, given by the following theorem.

######  Theorem 3.4.

Under model (6) with $p=1$, if $\nu_{1},\nu_{2}>0$, then
$\displaystyle\lim_{n\to\infty}nR(\pi_{1},\pi_{2})=\frac{V(\theta)f(\theta)}{2(\beta_{2}-\beta_{1})}$.

It follows that in order to minimize $\lim_{n\to\infty}nR(\pi_{1},\pi_{2})$ we
have to minimize $V(\theta)$, which by Theorem 3.1 is achieved by taking
$\pi_{1}({\bf x})=\pi^{0}_{1}({\bf x})$ and $\pi_{2}({\bf x})=\pi^{0}_{2}({\bf
x})$.

We remark that by a delta method calculation, the asymptotic variance of
$\widehat{\theta}=\frac{\widehat{\alpha}_{1}-\widehat{\alpha}_{2}}{\widehat{\beta}_{2}-\widehat{\beta}_{1}}$
is $\frac{V(\theta)}{(\beta_{2}-\beta_{1})^{2}}$, which is proportional to the
limit of $nR(\pi_{1},\pi_{2})$. Thus, minimizing this limit is equivalent to
minimizing the asymptotic variance of the estimator of the intersection point
$\theta$.

We return to general dimension $p$. Set
${\boldsymbol{\beta}}_{k}^{t}:=(\beta_{k,1},\ldots,\beta_{k,p})$, $k=1,2$,
${\boldsymbol{\beta}}_{k,-1}^{t}:=(\beta_{k,2},\ldots,\beta_{k,p})$ and
without loss of generality let ${\beta}_{2,1}>{\beta}_{1,1}$. For given ${\bf
x}_{-1}=(x_{2},\ldots,x_{p})$ define $\theta_{1}:=\theta_{1}({\bf
x}_{-1}):=\frac{\alpha_{1}-\alpha_{2}+({\boldsymbol{\beta}}_{1,-1}-{\boldsymbol{\beta}}_{2,-1})^{t}{\bf
x}_{-1}}{\beta_{2,1}-\beta_{1,1}}$. If $\theta_{1}\in[0,1]$, treatment $T_{1}$
is better for a covariate vector $\bf x$ satisfying
$x_{1}<\theta_{1}=\theta_{1}({\bf x}_{-1})$, and otherwise $T_{2}$ is better.
The limit of the regret is given in the following theorem.

######  Theorem 3.5.

Assume model (6) and $\nu_{1},\nu_{2}>0$, then

$\lim_{n\to\infty}nR(\pi_{1},\pi_{2})=\frac{1}{2(\beta_{2,1}-\beta_{1,1})}\int_{[0,1]^{p-1}}V(\theta_{1},{\bf
x}_{-1})f(\theta_{1},{\bf x}_{-1})I\\{\theta_{1}\in[0,1]\\}d{\bf x}_{-1}.$
(11)

Thus, asymptotically, $nR(\pi_{1},\pi_{2})$ is proportional to the integral of
the variance of the estimate of the intersection curve $\theta_{1}({\bf
x}_{-1})$, with weights proportional to the density at the intersection points
and by (11) $\lim_{n\to\infty}nR(\pi_{1},\pi_{2})$ is minimized by the
allocation $\pi^{0}_{1}({\bf x}),\pi^{0}_{2}({\bf x})$ since this alloction
minimizes $V$ uniformly.

The following corollary of Theorem 3.3 shows that $\pi_{1}^{0},\pi_{2}^{0}$,
which minimize the ideal, regret are also approximately optimal for the regret
$R(\pi_{1},\pi_{2})$ itself.

######  Corollary 3.1.

For any $\varepsilon>0$ there exists $C>0$ such that for any allocation
functions $\pi_{1},\pi_{2}$ with $\nu_{1},\nu_{2}>0$

$R(\pi_{1},\pi_{2})\geq
R(\pi_{1}^{0},\pi_{2}^{0})-\frac{C}{n^{3/2-\varepsilon}}.$

In particular, the optimal design that minimizes $R(\pi_{1},\pi_{2})$ can
improve the allocation $\pi_{1}^{0},\pi_{2}^{0}$ by at most
${C}/{n^{3/2-\varepsilon}}$.

The error term is meaningful as it is of smaller order than the regret, which
is $1/n$ according to Theorem 3.5.

###  3.2 A numerical example

To illustrate the asymptotic approximations of the regret and the optimal
allocations derived in Section 3.1 consider an example of model (6) with
$p=1$,
$\alpha_{1}=0.2,\,\alpha_{2}=0,\,\beta_{1}=0.5,\,\beta_{2}=1,\,\sigma_{1}^{2}=0.1,\,\sigma_{2}^{2}=0.2$,
and $X\sim U[0,1]$. The optimal allocation is $\nu_{1}=0.414$. The $\sigma$’s
are chosen such that $R^{2}$ in both regressions is about $0.2$. Figure 1(a)
shows the regression lines. For $x$ smaller than $\theta=0.4$, treatment 1 is
better and otherwise treatment 2 is preferred. The regret (10) is evaluated
using simulations (with $10^{5}$ replications) and is compared to the ideal
regret (3.1). Two scenarios are considered for the residuals in the regression
model: normal and centered exponential. Figures (b) and (c) show the regret
and ideal regret for $n=100$ under the allocation function
$\pi_{1}(x)=\nu_{1}$ where $\nu_{1}$ vary from 0.2 to 0.6. The allocation
optimizing the ideal regret is marked with a vertical line. While there is a
slight deviation between the ideal and actual regret, it seems that the
optimal allocations, with respect to both, are very close. Thus, the
approximation of Theorem 3.3 works well also for small $n$ in this example.
Figure 1 (d) shows $nR_{I}(\pi_{1},\pi_{2})$ of the balanced design
($\pi_{1}=\pi_{2}=1/2$) for large $n$ and its limit, see Theorem 3.4. It is
shown that the regret slowly converges to its limit, and it is close to the
limit only when $n$ is quite large.

(a) The regression lines

(b) Actual and ideal regret - normal residual

(c) Actual and ideal regret - exponential residual

(d) Ideal regret and asymptotic regret

Figure 1: Figure (a): the regression lines. Figures (b) and (c) show the
regret (computed by a simulation with 95% confidence intervals), and the ideal
regret for $n=100$ and different allocation ratios, where in (b) the residual
is normal and in (c) it is exponential (centered). The regret is calculated
for the allocation function $\pi_{1}(x)=\nu_{1}$ where $\nu_{1}$ varies from
0.2 to 0.6. The optimal allocation is marked by the vertical line. Figure (d)
shows $n$ times the ideal regret (blue line) and its limit (gray line).

###  3.3 Polynomial Regression

We now consider the model of (6) where the functions $g_{k}(x)$ are assumed to
be polynomials of a single continuous covariate $X$ with density $f$, i.e.,
$g_{k}(X)=\alpha_{k}+\sum_{j=1}^{J}\beta_{jk}X^{j}$, where $J$ denotes the
degree of the polynomial. Let $\theta_{1}<\cdots<\theta_{L}$, where $L\leq J$,
be the crossing points of $g_{1}(x)$ and $g_{2}(x)$, and define $\theta_{0}=0$
and $\theta_{L+1}=1$. Assume without loss of generality that
$g_{1}(x)>g_{2}(x)$ for $x\in(\theta_{\ell},\theta_{\ell+1})$ if $\ell$ is
odd, and the reverse inequality holds when $\ell$ is even. Therefore, the
regret is

$R(\pi_{1},\pi_{2})=\sum_{\ell\text{ is
odd}}\int_{\theta_{\ell}}^{\theta_{\ell+1}}P(\widehat{g}_{2}(x)>\widehat{g}_{1}(x))(g_{1}(x)-g_{2}(x))f(x)dx\\\
+\sum_{\ell\text{ is
even}}\int_{\theta_{\ell}}^{\theta_{\ell+1}}P(\widehat{g}_{1}(x)>\widehat{g}_{2}(x))(g_{2}(x)-g_{1}(x))f(x)dx.$

Notice that Theorems 3.3 and 3.5, which provide results for multivariate ${\bf
X}$, do not apply in this case because the random vector ${\bf
X}=(X,X^{2},\ldots,X^{J})$ does not have a joint density. However, using the
arguments of Theorem 3.3 (under similar notation and regularity conditions) it
can be shown that
$n^{3/2-\varepsilon}|R(\pi_{1},\pi_{2})-R_{I}(\pi_{1},\pi_{2})|\to 0$, for any
$\varepsilon>0$, where

$R_{I}(\pi_{1},\pi_{2})=\int_{0}^{1}\Phi\left(\frac{-\sqrt{n}|g_{2}(x)-g_{1}(x)|}{\sqrt{V(x)}}\right)|g_{2}(x)-g_{1}(x)|f(x)dx,$
(12)

and now ${\bf x}^{t}=(1,x,\ldots,x^{J})$,
${\boldsymbol{Q}}_{k}=\int_{\mathbb{R}^{p}}{\bf x}{\bf x}^{t}f_{k}({x})d{x}$
and
${\boldsymbol{\Sigma}}_{k}=\frac{1}{\nu_{k}}\sigma^{2}_{k}{\boldsymbol{Q}}_{k}^{-1},~{}~{}k=1,2$,
and $V(x)={\bf
x}^{t}\left({\boldsymbol{\Sigma}}_{1}+{\boldsymbol{\Sigma}}_{2}\right){\bf
x}$, which is the asymptotic variance of
$\widehat{g}_{1}(x)-\widehat{g}_{2}(x)$.

Lemma 3.2 implies that the design where
$\pi_{1}(x)=\frac{\sigma_{1}}{\sigma_{1}+\sigma_{2}}$ minimizes $V(x)$
uniformly over all $x$’s. It follows that this design is asymptotically
optimal also for this problem.

By Proposition .1 (in the appendix) we have

$\lim_{n\to\infty}nR(\pi_{1},\pi_{2})=\sum_{\ell=1}^{L}\frac{f(\theta_{\ell})V(\theta_{\ell})}{2|({\boldsymbol{\beta}}_{2}-{\boldsymbol{\beta}}_{1})^{t}{\boldsymbol{\zeta}}_{\ell}|},$

where
${\boldsymbol{\zeta}}_{\ell}:=(1,2\theta_{\ell},...,J\theta_{\ell}^{J-1})^{t}$.

A careful inspection of the proofs of the above results shows that the
optimality of the design with
$\pi_{1}(x)=\frac{\sigma_{1}}{\sigma_{1}+\sigma_{2}}$ generalizes to
regression functions of the form $g_{k}(X)=\sum_{j}\beta_{jk}h_{j}(X)$ for any
functions $h_{j}$ having bounded derivatives.

##  4 $K$ Treatments and one Covariate

###  4.1 Regret and ideal regret

We consider the case $p=1$ and
$g_{k}({X})=\alpha_{k}+{\beta}_{k}{X},~{}Var(Y|{X},T_{k})=\sigma_{k}^{2},~{}k=1,\ldots,K$.
As in Section 3, we assume the existence of the moment generating function of
$\varepsilon$ when conditioned on $({\bf X},T_{k})$. We also assume that ${X}$
is continuous with density $f({x})$ and is supported on $[0,1]$. Let
$\theta_{k}$, $k=0,\ldots,K$ be in increasing order, and such that treatment
$k$ is best in the interval $(\theta_{k-1},\theta_{k})$; we assume that each
treatment is best in some open interval, or equivalently that the intervals
are nonempty. Equations (2) and (5) for the regret and ideal regret imply
immediately

$R(\pi_{1},\ldots,\pi_{K})=\sum_{k=1}^{K}\sum_{m=1}^{K}\int_{\theta_{m-1}}^{\theta_{m}}P\left(\widehat{g}_{k}(x)>\max_{\ell\neq
k}\widehat{g}_{\ell}(x)\right)\left[g_{m}({x})-g_{k}({x})\right]f({x})d{x}$

and the ideal regret is

$\displaystyle R_{I}(\pi_{1},\ldots,\pi_{K})$
$\displaystyle=\sum_{k=1}^{K}\sum_{m=1}^{K}\int_{\theta_{m-1}}^{\theta_{m}}\left(\int\prod_{l=1,\ldots,K,l\not=k}\Phi\left(\frac{z\xi_{l}(x)+\sqrt{n}[g_{k}({x})-g_{l}({x})]}{\xi_{k}(x)}\right)\varphi(z)dz\right)\left[g_{m}({x})-g_{k}({x})\right]f({x})d{x}.$
(13)

Similar to the results in Section 3, we give asymptotic results on the rate of
convergence of the regret to the ideal regret $R_{I}$ as well as the limit of
the regret. Given allocation functions $\pi_{k}(x)$, the distribution of the
estimated regression coefficients $(\widehat{\alpha}_{k},\widehat{\beta}_{k})$
is approximately bivariate normal with means $(\alpha_{k},\beta_{k})$ and
covariance

${\boldsymbol{\Sigma}}_{k}=\frac{\sigma^{2}_{k}}{\nu_{k}}\left(\begin{array}[]{ll}\frac{\tau^{2}_{k}+\mu_{k}^{2}}{\tau^{2}_{k}}&\frac{-\mu_{k}}{\tau^{2}_{k}}\\\
\frac{-\mu_{k}}{\tau^{2}_{k}}&\frac{1}{\tau^{2}_{k}}\end{array}\right)\,,$
(14)

where we assume positivity of $\nu_{k}:=\int f({x})\pi_{k}({x})d{x}$, and
denote the mean and variance of the covariate in group $k$ by
$\mu_{k}:=\int_{0}^{1}xf_{k}(x)\,dx$ and
$\tau_{k}^{2}:=\int_{0}^{1}(x-\mu_{k})^{2}f_{k}(x)\,dx$, respectively,
$k=1,\ldots,K$. Parallel to Theorems 3.3 and 3.4 we have

######  Theorem 4.1.

Under the assumptions in the beginning of Section 4.1 we have for any
$\varepsilon>0$

$\lim_{n\to\infty}n^{3/2-\varepsilon}\left|R(\pi_{1},\ldots,\pi_{K})-R_{I}(\pi_{1},\ldots,\pi_{K})\right|=0.$

######  Theorem 4.2.

Under the assumptions in the beginning of Section 4.1, we have

$n\lim_{n\to\infty}R(\pi_{1},\ldots,\pi_{K})=\sum_{m=1}^{K-1}\frac{V_{m}(\theta_{m})f(\theta_{m})}{2|\beta_{m+1}-\beta_{m}|},$
(15)

where
$V_{m}(x)=(1,x)\left({\boldsymbol{\Sigma}}_{m}+{\boldsymbol{\Sigma}}_{m+1}\right)\binom{1}{x}$.

Using (14) it is possible to write $V_{m}(\theta_{m})$ explicitly as follows

$V_{m}(\theta_{m})=\frac{\sigma_{m}^{2}}{\nu_{m}}\left[1+\frac{(\theta_{m}-\mu_{m})^{2}}{\tau_{m}^{2}}\right]+\frac{\sigma_{m+1}^{2}}{\nu_{m+1}}\left[1+\frac{(\theta_{m}-\mu_{m+1})^{2}}{\tau_{m+1}^{2}}\right].$

Theorem 4.2 implies that asymptotically the optimal allocation problem reduces
to minimize a weighted average of the variances of
$\widehat{\theta}_{1},\ldots,\widehat{\theta}_{K-1}$ (see the discussion after
Theorem 3.4), which are the estimates of the intersection points. The weights
depend on the $\beta$’s and on the density $f$ at the intersection points.
Note also that $V_{m}(\theta_{m})$ is the sum of the asymptotic variances of
$\widehat{\alpha}_{m}+\widehat{\beta}_{m}\theta_{m}$ and
$\widehat{\alpha}_{m+1}+\widehat{\beta}_{m+1}\theta_{m}$. Corollary 3.1
continues to hold for $K$ treatments and thus, optimizing the ideal regret
approximately optimizes the regret itself.

###  4.2 A lower bound for the asymptotic regret (15)

By (4.1), (14), and (15) (see also (3)) both the ideal and the asymptotic
regret depend on the functions $\pi_{k}(x)$ only via the quantities $\nu_{k}$,
$\mu_{k},\tau_{k}^{2},k=1,\ldots,K$. Therefore, we can minimize (4.1) and (15)
in these parameters to obtain a lower bound for the regret. Let
$\mu:=\int_{0}^{1}xf(x)\,dx$, $\tau^{2}:=\int_{0}^{1}(x-\mu)^{2}f(x)\,dx$ and
recall that $f(x)=\sum_{k=1}^{K}\nu_{k}f_{k}(x)$. Then, by the relations for
central moments of mixture distributions we have

$\sum_{k=1}^{K}\nu_{k}=1,\quad\sum_{k=1}^{K}\nu_{k}\mu_{k}=\mu,\quad\sum_{k=1}^{K}\nu_{k}(\tau_{k}^{2}+\mu_{k}^{2})=\tau^{2}+\mu^{2}.$
(16)

Minimizing (4.1) or (15) in $\nu_{k}$, $\mu_{k},\tau_{k},k=1,\ldots,K$ subject
to the constraints (16) and $\nu_{k},\tau_{k}\geq 0,k=1,\ldots,K$ yields a
lower bound for the achievable regret. To obtain the minimum regret (instead
of a lower bound) one needs to add additional constraints on
$\mu_{k},\tau_{k},k=1,\ldots,K$ to restrict optimization to values for which
there exist mixture components $f_{k}$ with weighted sum $f$ that assume these
moments.

We now consider an example with uniform $f$ and $K=3$, and minimize the
asymptotic regret (15). Since now $f(x)\leq 1$ the constraint
$f(x)=\sum_{k=1}^{3}\nu_{k}f_{k}(x)$ implies that $\nu_{k}f_{k}(x)\leq
1,\,k=1,2,3,\,x\in[0,1]$. By Lemma 1 below this implies that the variance of
$f_{k}(x),k=1,2,3$ is bounded from below by $\nu_{k}^{2}/12$.

As an example consider the setting where
$(\alpha_{1},\alpha_{2},\alpha_{3})=(0.0,-0.1,-1.2),\,(\beta_{1},\beta_{2},\beta_{3})=(0.2,0.5,2.0)$,
and the residual variances $\sigma_{k}^{2},\,k=1,2,3$ are equal across
treatment groups. Figure 2 shows the resulting scenario and optimized
allocation probabilities (whose computation is explained below). Then,
numerically (see details below) minimizing (15) in $\nu_{k}$,
$\mu_{k},\tau_{k},k=1,\ldots,K$ subject to the constraints (16) as well as
$\tau_{k}^{2}\geq\nu_{k}^{2}/12,\,k=1,2,3$ gives
$\nu_{k}=0.346,\,0.444,\,0.210,\,\mu_{k}=0.342,\,0.512,\,0.735,\,\tau^{2}_{k}=0.00997,\,0.132,\,0.00369$,
respectively, and the minimized value of (15) is $12.128$. We note that
$\tau_{k}^{2}=\nu_{k}^{2}/12$, for $k=1,3$. Applying Lemma 1 for $c=1/\nu_{k}$
for $k=1,3$, the only density $f_{k}(x)$ with mean $\mu_{k}$ and variance
$\tau_{k}^{2}$ and satisfying the constraint $f_{k}(x)\leq
1/\nu_{k},\,x\in[0,1]$ is the uniform distributions on
$[\mu_{k}-\nu_{k}/2,\,\mu_{k}+\nu_{k}/2]$ with densities
$f_{k}(x)=1_{[\mu_{k}-\nu_{k}/2,\,\,\mu_{k}+\nu_{k}/2]}(x)/\nu_{k}$.
Furthermore, we set

$f_{2}(x)=\frac{1-\nu_{1}f_{1}(x)-\nu_{3}f_{3}(x)}{\nu_{2}}\,.$

Because the supports of $f_{1},f_{3}$ are disjoint and $f_{k}(x)\leq
1/\nu_{k},\,k=1,3$, $f_{2}(x)$ is a valid density. By construction, the mean
and variance of $f_{2}(x)$ are $\mu_{2},\tau^{2}_{2}$ and the three densities
satisfy $f(x)=\sum_{k=1}^{3}\nu_{k}f_{k}(x)$. Thus, the resulting allocation
probabilities $\pi_{k}(x)=f_{k}(x)/\nu_{k}$ achieve the lower bound computed
above, and hence, assuming the minimization algorithm yielded the correct
minimum (which we verified in several ways, including the fact that it is
unique), it is an optimal allocation (see Figure 2 for a plot). It follows
that in this case the above constraints suffice. Furthermore, as $f_{1},f_{3}$
are unique, they are a unique optimal allocation (under the assumption that
the minimization in $\nu_{k}$, $\mu_{k},\tau_{k},k=1,\ldots,K$ gives a unique
solution). It also follows that the optimal solution leads to a deterministic
allocation.

|
---|---
Figure 2: The functions $g_{k}$ for the three treatments (left) and the
optimal allocation probabilities minimizing (15) (as a stacked area chart,
right) for the example in Section 4.2. For instance, for $x$’s under the black
area, treatment 1 is assigned with probability $\pi_{1}(x)=1$.

Some comments: (1) To implement the numerical optimization of (15) subject to
the constraints, we use the R command optim with the following
reparametrization, which allows one to apply optimization algorithms for
unrestricted optimization. Let $A_{i}=(a_{i,1},a_{i,2})\in R^{2},\,i=1,2,3$.
First we transform the vectors $A_{i}$ into proportions setting
$a^{\prime}_{ik}=\exp(a_{i,k})/\sum_{m=1}^{3}\exp(a_{i,m}),\,i,k=1,2,3$ where
$a_{i,3}\equiv 0$, and set
$\nu_{k}=a^{\prime}_{1k},\,\mu_{k}=a^{\prime}_{2k}\mu/\nu_{k}$, and

$\tau^{2}_{k}=\frac{\nu_{k}^{2}}{12}+\frac{a^{\prime}_{3,k}\left(\mu^{2}+\tau^{2}-\sum_{m=1}^{3}(\nu_{m}\mu_{m}^{2}+\nu_{m}^{3}/12)\right)}{\nu_{k}},\,\,k=1,2,3.$

As long as the second summand remains positive, this parametrization ensures
that $\tau^{2}_{k}\geq\nu_{k}^{2}/12$ and that the constraints of (16) hold.
To confirm the robustness of the solution, the optimization is performed
repeatedly with randomly chosen initial conditions. (2) If the intersection
points $\theta_{k,k+1}$ of the lines $g_{k}$ move closer together so will the
intervals of $x$’s where all patients are allocated to Treatment 1 or
Treatment 3 (the blue and black regions in Figure 2). For scenarios where they
overlap, we can no longer use the above approach to obtain an optimal solution
and the solutions seem to become more complex. This is the case, for example,
for $\alpha_{k}=0.0,\,-0.1,\,-.8,\,\beta_{k}=0.2,\,0.5,\,2.0,\,\,k=1,2,3$.

##### Lemma 1

Let $c>0$. Among all continuous distributions $f$ with (i) mean $\mu$ such
that (ii) $f(x)\leq c$ for all $x\in\mathbb{R}$, the uniform distribution on
$[\mu-1/(2c),\mu+1/(2c)]$ has the smallest variance. This minimum variance is
$1/(c^{2}12)$.

##  5 A clinical trial to determine personalized diets

To demonstrate our approach, we consider a study by Ebbeling et al. (2018) in
which three diets were compared in a parallel group design. We fitted a model
to the data of that trial and derived an alternative optimal design based on
the estimated model parameters.

Patients were randomly assigned to three diet groups. The diets differ in
their carbohydrate content, high (60%), moderate (40%) or low (20%) and were
given for 20 weeks. The total number of subjects allocated to these three
treatment groups was 54, 53 and 57, respectively, so the design is (almost)
balanced. The main covariate is insulin secretion (insulin concentration 30
minutes after oral glucose). The primary outcome in the original trial was
averaged total energy expenditure over two measurements, in the middle and the
end of the trial. Total energy expenditure is measured by doubly labeled
water; for details on this measure see Hills et al. (2014).

Of the 1685 subjects initially screened, 234 participated in a run-in phase
which preceded the trial itself. Of these, 164 achieved a target of 12% ($\pm$
2%) weight loss and qualified for the trial, and were randomly assigned to one
of the three diets.

To account for the cost of the diets, either in terms of money or in terms of
inconvenience to the dieter, we consider the utility of the diets as an
outcome variable which is defined as the energy expenditure minus the diet
cost. We assume that the lower the carbon content of the diet, the higher the
cost (see Hagberg et al. (2019) for a similar definition) but otherwise our
choice of the cost for this example is arbitrary. On the basis of Figure 4 in
Ebbeling et al. (2018), we assume that under treatment $T_{k}$, $k=1,2,3$,

$\text{utility}=\alpha_{k}+\beta_{k}X+\varepsilon_{k}-\text{cost}_{k},$ (17)

where costk= 0, 150, 300 for $k=1,2,3$, respectively, $X$ is insulin secretion
and its distribution is Gamma(3.12, 0.02), and the standard deviations of
$\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}$ are 190, 150, 130,
respectively. The design we propose is based on assumed values of
$\alpha_{k}$, which are 40, -80, -240, and of $\beta_{k}$, which are -0.8,
0.1, and 0.8 for $k=1,2,3$, respectively. These values and the parameters of
the gamma distribution of $X$ and the standard deviation of $\varepsilon$’s
are all based on Figure 4 of Ebbeling et al. (2018), which plots the energy
expenditure as a function of the insulin secretion for the three diets in the
data. In particular, the parameters of the Gamma distribution are chosen to
match empirical quantiles of the $X$’s.

In order to numerically find the optimal allocation probabilities $\pi_{k}(x)$
that minimize the regret, we use a certain parameterization. Let ${\bf
a}=(a_{0},\ldots,a_{m})$, and set $h_{\bf
a}(x)=\exp(\sum_{j=0}^{m}a_{j}x^{j})$. Furthermore, let ${\bf A}$ denote a
$2\times(m+1)$ matrix of real numbers. Let
$\pi_{k}(x)=\frac{h_{A_{k}}(x)}{\sum_{j=0}^{3}h_{A_{j}}(x)}$, $k=1,2,3$, where
$A_{k}$ denotes the $k$-th row of ${\bf A}$ and we set $h_{A_{3}}(x)\equiv 1$.
Now we can approximately minimize the (ideal) regret by plugging $\pi_{k}(x)$
into (4.1) and minimizing with respect to ${\bf A}$. In this example we chose
$m=4$ and use the R function optim for the numerical optimization.

Figure 3(a) plots the regression lines (17) for the three diets and the
assumed distribution of the covariate $x$. The high, moderate, and low
carbohydrate diets are optimal when $X$ belongs to the intervals
$(0,133),(133,229),(229,\infty)$, respectively. The difference between the
treatments is more pronounced for extreme values of $X$. Hence, people with
very small or very large $X$ would benefit more from a personalized treatment
choice compared to those with medium values of $X$.

Figures 3(b)-(d) present the optimal allocation probabilities, as a stacked
area chart, when $n=164,\,1000$, and $\infty$. When the sample size is larger,
the optimal allocation rule is closer to being deterministic. For $n=164$,
which was the actual trial size, the optimal trial allocation rule tends to
allocate subjects to $T_{3}$ (respectively, $T_{2}$) for large (respectively,
small) values of $X$. Compared to equal allocation, which was the actual
design used, there is a reduction of about $11\%$ of the regret. A saving of
11% may be quite significant if the treatments are applied to a large number
of future patients. Furthermore, the regret of the balanced design with
$n=164$ can be achieved under the optimal design with $n=143$, which
represents a reduction of 13% of the sample size. Taking into account that in
this experiment only about 10% of recruited patients passed the screening
process, a saving of about 20 subjects amounts to a reduction of about 200
subjects that would have to be recruited and pretested.

To investigate if the reduction in the regret is due to unbalanced allocation
or the covariate dependent allocation, we also computed the regret for the
optimal unbalanced allocation that is restricted to allocation probabilities
that do not depend on the covariate. We refer to this design as unbalanced.
This can be achieved as above setting $m=0$. We found numerically the optimal
unbalanced design for $n=164$ leads to reduction of the regret by about 8%,
whereas optimizing with allocation functions that depend on the covariate,
leads to a reduction by an additional 3%. For larger sample sizes, it appears
that the gain by optimized covariate dependent allocation probabilities
increases (see Table 1). As the algorithm only gives an approximate optimal
solution, we also used the approach in Section 4.2 to compute a lower bound
for the limit $nR_{I}$ for $n\to\infty$. The obtained lower bound of 735.1 is
close to 738.2, the corresponding regret obtained by the derived allocation
function depicted in Figure 3d.

It should be noticed that all these designs are based on the true values of
the parameters and therefore require a preliminary study, or an adaptive
design.

|  | % reduction | optimized allocation
---|---|---|---
$n$ | $nR_{I}$ | optimal | unbalanced | high | moderate | low
$164$ | 842.8 | 10.9 | 7.9 | 45.3 | 36.7 | 18.3
$1000$ | 829.6 | 13.1 | 4.4 | 36.8 | 39.8 | 23.7
$\infty$ | 738.2 | 17.1 | 3.2 | 34.6 | 40.2 | 25.5
Table 1: Diets example: $nR_{I}$, the reductions of the regret of the optimal
design compared to a balanced design, the reductions of the regret of a design
with optimized allocation ratios only (where allocation probabilities do not
depend on the $X$) compared to the balanced design, and the optimized
allocation probabilities in percent (for the design where allocation
probabilities do not depend on $X$).

(a) the regression models

(b) optimal allocation for $n=164$

(c) optimal allocation for $n=1000$

(d) optimal allocation for $n\to\infty$

Figure 3: Figure (a) plots the regression models (17) and the assumed density
of the covariate. Figures (b)-(d) present the optimal allocation probabilities
(as a stacked area chart) when $n=164,1000,\infty$.

## Appendix: Proofs

#### Proof of Lemma 3.2

We use the convexity relation $\lambda A^{-1}+(1-\lambda)B^{-1}\succeq[\lambda
A+(1-\lambda)B]^{-1}$ (Moore, 1973) where $A$ and $B$ are positive definite
matrices. We have

${\boldsymbol{\Sigma}}_{1}+{\boldsymbol{\Sigma}}_{2}={\sigma_{1}}\left(\frac{{\nu_{1}}{\boldsymbol{Q}}_{1}}{\sigma_{1}}\right)^{-1}+{\sigma_{2}}\left(\frac{{\nu_{2}}{\boldsymbol{Q}}_{2}}{\sigma_{2}}\right)^{-1}\\\
\succeq(\sigma_{1}+\sigma_{2})\Big{[}\frac{\sigma_{1}}{\sigma_{1}+\sigma_{2}}\left(\frac{{\nu_{1}}{\boldsymbol{Q}}_{1}}{\sigma_{1}}\right)+\frac{\sigma_{2}}{\sigma_{1}+\sigma_{2}}\left(\frac{{\nu_{2}}{\boldsymbol{Q}}_{2}}{\sigma_{2}}\right)\Big{]}^{-1}=(\sigma_{1}+\sigma_{2})^{2}\mathbb{\boldsymbol{Q}}^{-1},$

where in the last equality we used the fact that
$\nu_{1}{\boldsymbol{Q}}_{1}+{\nu_{2}}{\boldsymbol{Q}}_{2}=\mathbb{\boldsymbol{Q}}$.
Equality holds when the two matrices to the left of the $\succeq$ sign are
equal, which happens when $\pi_{1}({\bf
x})=\frac{\sigma_{1}}{\sigma_{1}+\sigma_{2}}$ and $\pi_{2}({\bf
x})=\frac{\sigma_{2}}{\sigma_{1}+\sigma_{2}}$. ∎

#### Proof of Theorem 3.3

We start with the case of $p=1$. It is enough to show (as the other part of
the integral, from $\theta$ to 1, is symmetric) that

$n^{3/2-\varepsilon}\int_{0}^{\theta}\left|P(\widehat{g}_{2}(x)>\widehat{g}_{1}(x))-\Phi\left(\frac{-\sqrt{n}[g_{1}(x)-g_{2}(x)]}{\sqrt{V(x)}}\right)\right|[g_{1}(x)-g_{2}(x)]f(x)dx\to
0.$ (18)

The idea of the proof is quite simple. For $x\in\theta\pm 1/\sqrt{n}$ the
quantity in absolute value converges to zero, and the range of the integral
and $g_{1}(x)-g_{2}(x)$ are of order $1/\sqrt{n}$ each. For other $x$’s both
the probability and $\Phi$ in the integrand are exponentially small in $n$.

We will prove (18) by showing that

$n^{3/2-\varepsilon}\int_{0}^{\theta}\left|P(\widehat{g}_{2}(x)>\widehat{g}_{1}(x))-E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)\right|[g_{1}(x)-g_{2}(x)]f(x)dx\to
0\text{ and}$ (19)
$n^{3/2-\varepsilon}\int_{0}^{\theta}\left|E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)-\Phi\left(\frac{-\sqrt{n}[g_{1}(x)-g_{2}(x)]}{\sqrt{V(x)}}\right)\right|[g_{1}(x)-g_{2}(x)]f(x)dx\to
0,$ (20)

where

$S^{2}(x)=n(1,x)\left\\{\sigma_{1}^{2}\left(\sum_{i\in
T_{1}}\binom{1}{X_{i}}(1,X_{i})\right)^{-1}+\sigma_{2}^{2}\left(\sum_{i\in
T_{2}}\binom{1}{X_{i}}(1,X_{i})\right)^{-1}\right\\}\binom{1}{x},$ (21)

and the expectation in (19) and (20) is with respect to $X_{1},\ldots,X_{n}$
which appear in $S^{2}(x)$ and the random sets $T_{k}$. We start with (20).
Since $\theta$ is the intersection point,
$\alpha_{1}+\beta_{1}\theta=\alpha_{2}+\beta_{2}\theta$; therefore,
$g_{1}(x)-g_{2}(x)=\alpha_{1}+\beta_{1}x-(\alpha_{2}+\beta_{2}x)=(\beta_{2}-\beta_{1})(\theta-x)$.
Hence, (20) can be written as

$n^{3/2-\varepsilon}\int_{0}^{\theta}\left|E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)-\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)\right|(\beta_{2}-\beta_{1})(\theta-x)f(x)dx\to
0.$

Fix $c_{n}=n^{\varepsilon/5}$. We will show that

$n^{3/2-\varepsilon}\int_{0}^{\theta-
c_{n}/\sqrt{n}}\left|E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)-\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)\right|(\beta_{2}-\beta_{1})(\theta-x)f(x)dx\to
0,$ (22)

and that

$n^{3/2-\varepsilon}\int_{\theta-
c_{n}/\sqrt{n}}^{\theta}\left|E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)-\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)\right|(\beta_{2}-\beta_{1})(\theta-x)f(x)dx\to
0.$ (23)

We show (22) by arguing that the probabilities in (22) are exponentially
small. With notation defined in Section 2 we have that

$V(x)\leq(1+x^{2})\left(\frac{\sigma_{1}^{2}}{\nu_{1}\lambda_{min}({\boldsymbol{Q}_{1}})}+\frac{\sigma_{2}^{2}}{\nu_{2}\lambda_{min}({\boldsymbol{Q}_{2}})}\right)\leq
C,$

for a constant $C>0$ (since $x$ is bounded). With the bound
$\Phi(-t)\leq\exp(-t^{2}/2)$ for $t>1$ and $(\theta-x)\sqrt{n}\leq c_{n}$ we
have

$\Phi\left(\frac{-(\beta_{2}-\beta_{1})(\theta-x)\sqrt{n}}{\sqrt{V(x)}}\right)\leq\exp\left(-\frac{(\beta_{2}-\beta_{1})^{2}c_{n}^{2}}{2C}\right),$

which shows that the second normal probability $\Phi(\cdot)$ in (22) is
exponentially small. We now argue that the expectation in (22) is also
exponentially small. Consider the event

$A_{n}:=\left\\{\lambda_{min}\left(\frac{1}{n}\sum_{i\in
T_{1}}\binom{1}{X_{i}}(1,X_{i})\right)\leq\nu_{1}\frac{\lambda_{min}({\boldsymbol{Q}_{1}})}{2}\right\\}\\\
\bigcup\left\\{\lambda_{min}\left(\frac{1}{n}\sum_{i\in
T_{2}}\binom{1}{X_{i}}(1,X_{i})\right)\leq\nu_{2}\frac{\lambda_{min}({\boldsymbol{Q}_{2}})}{2}\right\\}.$
(24)

Writing $\frac{1}{n}\sum_{i\in
T_{1}}\binom{1}{X_{i}}(1,X_{i})=\frac{1}{n}\sum_{i=1}^{n}\binom{1}{X_{i}}(1,X_{i})I(i\in
T_{1})$, a sum of iid terms, we can apply Theorem 5.1 of Tropp (2012) and the
union bound to conclude that $P(A_{n})\leq C\exp(-Cn)$ for some constant
$C>0$. Going back to (22), write

$E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)\\\
=E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)I(A_{n})+E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)I(A_{n}^{c})\\\
\leq
P(A_{n})+E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)I(A_{n}^{c}).$

When $A_{n}^{c}$ occurs, the denominator $S^{2}(x)$ is bounded and the same
argument as above shows that this probability is exponentially small. We
conclude that both probabilities in (22) are exponentially small, which
implies that (22) is correct.

We now prove (23). By the mean value theorem, for $x\in(\theta-
c_{n}/\sqrt{n},\,\theta)$ and for some $V^{*}(x)$ between $S^{2}(x)$ and
$V(x)$, we have for some $C>0$

$\left|\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{S^{2}(x)}}\right)-\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)\right|\\\
=|S^{2}(x)-V(x)|\varphi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V^{*}(x)}}\right)\frac{\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{2}(V^{*}(x))^{-3/2}\\\
\leq{C}|S^{2}(x)-V(x)|c_{n}(\min\\{S^{2}(x),V(x)\\})^{-3/2},\text{ and}$
$\\{S^{2}(x)-V(x)\\}=(1,x)\left\\{\frac{\sigma_{1}^{2}}{\nu_{1}}\left({\bf
M}_{1}^{-1}-{\boldsymbol{Q}}_{1}^{-1}\right)+\frac{\sigma_{2}^{2}}{\nu_{2}}\left({\bf
M}_{2}^{-1}-{\boldsymbol{Q}}_{2}^{-1}\right)\right\\}\binom{1}{x},$ (25)

where ${\bf M}_{k}=\frac{1}{n\nu_{k}}\sum_{i\in
T_{k}}\binom{1}{X_{i}}(1,X_{i})$, $k=1,2$. We have that ${\bf
M}_{1}^{-1}-{\boldsymbol{Q}}_{1}^{-1}={\bf
M}_{1}^{-1}\left({\boldsymbol{Q}}_{1}-{\bf
M}_{1}\right){\boldsymbol{Q}}_{1}^{-1}$, implying that the difference of the
inverses is small when the difference of the matrices is small and the inverse
is bounded. Define the event

$\widetilde{A}_{n}:=\\{\max\\{\|{\boldsymbol{Q}}_{1}-{\bf
M}_{1}\|_{\infty},\|{\boldsymbol{Q}}_{2}-{\bf M}_{2}\|_{\infty}\\}\geq
c_{n}/\sqrt{n}\\},$

where $\|\cdot\|_{\infty}$ denotes the maximum (entry-wise) norm. By
Hoeffding’s inequality (since the random variables are bounded),
$P(\widetilde{A}_{n})\leq C\exp(-Cc_{n}^{2})$ (which is exponentially small).
Therefore, as before, this event can be ignored. On the complement of
$\widetilde{A}_{n}$, ${\bf M}_{1}^{-1}$, ${\bf M}_{2}^{-1}$, $S^{2}(x)$ are
all bounded, and by similar arguments $S^{2}(x)$ and $V(x)$ are bounded below.
Hence, when $\widetilde{A}_{n}^{c}$ occurs,

$\left|\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)-\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)\right|\leq
C\frac{c_{n}^{2}}{\sqrt{n}}.$

Therefore the integral in (23) is bounded by

$C{n}^{1-\varepsilon}c_{n}^{2}\int_{\theta-
c_{n}/\sqrt{n}}^{\theta}(\theta-x)f(x)dx\leq
Cn^{1/2-\varepsilon}c_{n}^{3}\int_{\theta-c_{n}/\sqrt{n}}^{\theta}f(x)dx\leq
Cc_{n}^{4}/{n^{\varepsilon}}=n^{-\varepsilon/5}\to 0,$

which completes the proof of (23) and hence of (20).

We now show (19) by dividing the integral into two ranges: (a) $[\theta-
c_{n}/\sqrt{n},\theta]$, and (b) $[0,\theta-c_{n}/\sqrt{n}]$.

We have that

$P\left(\widehat{g}_{2}(x)>\widehat{g}_{1}(x)\right)=P\left(\left[\binom{\widehat{\alpha}_{1}-\alpha_{1}}{\widehat{\beta}_{1}-\beta_{1}}-\binom{\widehat{\alpha}_{2}-\alpha_{2}}{\widehat{\beta}_{2}-\beta_{2}}\right]^{t}\binom{1}{x}>(\beta_{2}-\beta_{1})(x-\theta)\right)\\\
=P\left(\left[\frac{{\bf M}_{1}^{-1}}{n\nu_{1}}\sum_{i\in
T_{1}}\binom{\varepsilon_{i}}{X_{i}\varepsilon_{i}}-\frac{{\bf
M}_{2}^{-1}}{n\nu_{2}}\sum_{i\in
T_{2}}\binom{\varepsilon_{i}}{X_{i}\varepsilon_{i}}\right]^{t}\binom{1}{x}>(\beta_{2}-\beta_{1})(x-\theta)\right)\\\
=P\left(\sum_{i=1}^{n}a_{i}\varepsilon_{i}>(\beta_{2}-\beta_{1})(x-\theta)\right),$

where $a_{i}=\left\\{\begin{array}[]{cc}(1,X_{i})\frac{{\bf
M}_{1}^{-1}}{n\nu_{1}}\binom{1}{x}&i\in T_{1}\\\ -(1,X_{i})\frac{{\bf
M}_{2}^{-1}}{n\nu_{2}}\binom{1}{x}&i\in T_{2}\end{array}\right.$, ${\bf
M}_{1},{\bf M}_{2}$ are defined in (25) and $\varepsilon_{i}$ are the
regression error terms. Notice that $Var(\varepsilon_{i})$ is equal to
$\sigma_{1}^{2}$ if $i\in T_{1}$ and to $\sigma_{2}^{2}$ if $i\in T_{2}$. We
are going to condition on $X_{1},\ldots,X_{n}$ obtained in the experiment and
their allocation, denoted together by $D$. Then
$Var(\sum_{i=1}^{n}a_{i}\varepsilon_{i}|D)=S^{2}(x)/n$, where $S^{2}(x)$ is
defined in (21).

We aim to apply the Berry Esseen Theorem, for which we need some bounds. The
values $na_{i}$ are bounded below since the matrices ${\bf M}_{k}$ are bounded
above. In order to obtain upper bounds, recall the event $A_{n}$ defined in
(24). On the set $A_{n}^{c}$ we have that $na_{i}$ are bounded above. Since
$P(A_{n})$ is exponentially small, we can ignore it, and assume that $na_{i}$
are bounded with probability one. We can apply the Berry-Esseen Theorem for
the non-identically distributed case (for a convenient reference see Chen et
al. (2010) (3.27)) to
$\sum_{i=1}^{n}\frac{\sqrt{n}a_{i}}{S(x)}\varepsilon_{i}\equiv\sum_{i=1}^{n}V_{i}$
where $\sum_{i=1}^{n}Var(V_{i}|D)=1$ and $E(|V_{i}|^{3}\mid D)\leq
max\\{na_{i}/S(x)E(|\varepsilon_{i}|^{3}|D):1\leq i\leq n\\}n^{-3/2}$ to
obtain by unconditioning on $D$

$\left|P\left(\sum_{i=1}^{n}a_{i}\varepsilon_{i}>(\beta_{2}-\beta_{1})(x-\theta)\right)-E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)\right|\leq
C/\sqrt{n}.$

Therefore,

$n^{3/2-\varepsilon}\int_{\theta-
c_{n}/\sqrt{n}}^{\theta}\left|P(\widehat{g}_{2}(x)>\widehat{g}_{1}(x))-E\Phi\left(\frac{-\sqrt{n}(\beta_{2}-\beta_{1})(\theta-x)}{S(x)}\right)\right|(\beta_{2}-\beta_{1})(\theta-x)f(x)dx\\\
\leq n^{3/2-\varepsilon}\frac{C}{\sqrt{n}}\int_{\theta-
c_{n}/\sqrt{n}}^{\theta}(\beta_{2}-\beta_{1})(\theta-x)f(x)dx\leq
n^{3/2-\varepsilon}\frac{C}{\sqrt{n}}\cdot\frac{c_{n}}{\sqrt{n}}\int_{\theta-
c_{n}/\sqrt{n}}^{\theta}f(x)dx\\\ \leq
C\frac{c_{n}^{2}}{{n}^{\varepsilon}}=C\frac{n^{2\varepsilon/5}}{{n}^{\varepsilon}}\to
0.$ (26)

By a standard large deviation bound, the integral of (19) from 0 to $\theta-
c_{n}/\sqrt{n}$ multiplied by $n^{3/2-\varepsilon}$ is easily shown to be of
order $n^{3/2-\varepsilon}O(exp(-c_{n}\sqrt{n}))\to 0$. This proves the
theorem for covariates in $\mathbb{R}$ ($p$=1).

We now consider the case of covariates in $\mathbb{R}^{p}$ for $p\geq 2$.
Suppose, without loss of generality, that ${\beta}_{2,1}>{\beta}_{1,1}$. For
given ${\bf x}_{-1}=(x_{2},\ldots,x_{p})$ define $\theta_{1}({\bf
x}_{-1}):=\frac{\alpha_{1}-\alpha_{2}+({\boldsymbol{\beta}}_{1,-1}-{\boldsymbol{\beta}}_{2,-1})^{t}{\bf
x}_{-1}}{\beta_{2,1}-\beta_{1,1}}$. For ${\bf x}$ such that $\theta_{1}({\bf
x}_{-1})\notin[0,1]$ the probability of making a mistake (i.e.,
$\widehat{g}_{2}({\bf x})-\widehat{g}_{1}({\bf x})$ has the wrong sign) is
exponentially small. Therefore,

$R(\pi_{1},\pi_{2})=\int_{[0,1]^{p-1}}\int_{0}^{\theta_{1}}P(\widehat{g}_{2}({\bf
x})>\widehat{g}_{1}({\bf x}))(g_{1}({\bf x})-g_{2}({\bf x}))f({\bf
x})dx_{1}I(\theta_{1}\in[0,1])d{{\bf x}}_{-1}\\\
+\int_{[0,1]^{p-1}}\int_{\theta_{1}}^{1}P(\widehat{g}_{2}({\bf
x})>\widehat{g}_{1}({\bf x}))(g_{1}({\bf x})-g_{2}({\bf x}))f({\bf
x})dx_{1}I(\theta_{1}\in[0,1])d{{\bf x}}_{-1}+a_{n},$

where $a_{n}$ is exponentially small. By a slight variation on the above one-
dimensional case, the inner integral equals

$\int_{0}^{1}\Phi\left(\frac{-\sqrt{n}|g_{2}({\bf x})-g_{1}({\bf
x})|}{\sqrt{V({\bf x})}}\right)|g_{2}({\bf x})-g_{1}({\bf x})|f({\bf
x})dx_{1}I(\theta_{1}\in[0,1])+o(1/n^{3/2-\varepsilon}),$ (27)

where the error term $o$ is uniform in ${\bf x}_{-1}$. Taking the outer
integral, the theorem follows. ∎

#### Proof of Theorem 3.4

As computed two lines below (21),
$g_{1}(x)-g_{2}(x)=(\beta_{2}-\beta_{1})(\theta-x)$. We will show that for
continuous $t$

$\lim_{t\to\infty}t\int_{0}^{\theta}\Phi\left(\frac{-\sqrt{t}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)(\beta_{2}-\beta_{1})(\theta-x)f(x)dx=\frac{V(\theta)f(\theta)}{4(\beta_{2}-\beta_{1})}.$
(28)

The integral from $\theta$ to 1 is similar, yielding the desired limit. By
L’Hôpital’s rule the limit in (28) equals

$\lim_{t\to\infty}\frac{t^{2}}{2}\int_{0}^{\theta}\varphi\left(\frac{\sqrt{t}(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\right)\frac{(\beta_{2}-\beta_{1})(\theta-x)}{\sqrt{V(x)}}\frac{1}{\sqrt{t}}(\beta_{2}-\beta_{1})(\theta-x)f(x)dx.$

Substitution of $x$ by $y=(\theta-x)\sqrt{t}$ in the integral yields

$\lim_{t\to\infty}\frac{1}{2}\int_{0}^{\theta\sqrt{t}}\varphi\left(\frac{(\beta_{2}-\beta_{1})y}{\sqrt{V(\theta-y/\sqrt{t})}}\right)\frac{(\beta_{2}-\beta_{1})^{2}}{\sqrt{V(\theta-y/\sqrt{t})}}y^{2}f(\theta-y/\sqrt{t})dy.$

The limit is

$\frac{1}{2}\int_{0}^{\infty}\varphi\left(\frac{(\beta_{2}-\beta_{1})y}{\sqrt{V(\theta)}}\right)\frac{(\beta_{2}-\beta_{1})^{2}}{\sqrt{V(\theta)}}y^{2}f(\theta)dy.$

Substitution of $y$ with $z=\frac{(\beta_{2}-\beta_{1})y}{\sqrt{V(\theta)}}$
in the integral yields the claimed limit
$\frac{f(\theta)V(\theta)}{2(\beta_{2}-\beta_{1})}\int_{0}^{\infty}\varphi\left(z\right)z^{2}dz=\frac{f(\theta)V(\theta)}{4(\beta_{2}-\beta_{1})}.$
∎

#### Proof of Theorem 3.5

Consider the inner integral in (27). By Theorem 3.4, the limit of the inner
integral times $n$ is

$\frac{V(\theta_{1},{\bf x}_{-1})f(\theta_{1},{\bf
x}_{-1})}{2(\beta_{2,1}-\beta_{1,1})}.$

By taking the outer integral the result follows. ∎

### The limit of the regret for polynomial regression

Suppose that
$g_{1}(x)=\alpha_{1}+(x,\ldots,x^{J})^{t}{\boldsymbol{\beta}}_{1}$ and
$g_{2}(x)=\alpha_{2}+(x,\ldots,x^{J})^{t}{\boldsymbol{\beta}}_{2}$. Let
$\theta$ be an intersection point, i.e., $g_{1}(\theta)=g_{2}(\theta)$. We
have

$g_{2}(x)-g_{1}(x)=g_{2}(x)-g_{2}(\theta)-[g_{1}(x)-g_{1}(\theta)]={\boldsymbol{\eta}}^{t}({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}}),$

where
${\boldsymbol{\eta}}:=(x-\theta,x^{2}-\theta^{2},\ldots,x^{J}-\theta^{J})^{t}$.
The following proposition is parallel to Theorem 3.4 (in the one-dimensional
case) and provides the limit of the regret; see (12). The proof is similar,
although some modifications are required.

######  Proposition .1.

For continuous $t$

$\lim_{t\to\infty}t\int_{0}^{\theta}\Phi\left(\frac{-\sqrt{t}({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\eta}}}{\sqrt{V(x)}}\right)({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\eta}}f(x)dx=\frac{f(\theta)V(\theta)}{4|({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\zeta}}|},$
(29)

where ${\boldsymbol{\zeta}}:=(1,2\theta,...,J\theta^{J-1})^{t}$.

Proof. By L’Hôpital’s rule the limit in (29) equals

$\lim_{t\to\infty}\frac{t^{2}}{2}\int_{0}^{\theta}\varphi\left(\frac{\sqrt{t}({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\eta}}}{\sqrt{V(x)}}\right)\frac{({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\eta}}}{\sqrt{V(x)}}\frac{1}{\sqrt{t}}({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\eta}}f(x)dx.$

Substituting $y=(\theta-x)\sqrt{t}$ we get

$\sqrt{t}{\boldsymbol{\eta}}=\sqrt{t}\begin{pmatrix}x-\theta\\\
x^{2}-\theta^{2}\\\ \vdots\\\
x^{J}-\theta^{J}\end{pmatrix}=\sqrt{t}\begin{pmatrix}x-\theta\\\
(x-\theta)(x+\theta)\\\ \vdots\\\
(x-\theta)(x^{J-1}+x^{J-2}\theta+\ldots+\theta^{J-1})\end{pmatrix}=y\widetilde{\boldsymbol{\eta}},$

where

$\widetilde{\boldsymbol{\eta}}:=(1,x+\theta,\ldots,x^{J-1}+x^{J-2}\theta+\ldots+\theta^{J-1})^{t}\\\
=(1,2\theta-y/\sqrt{t},\ldots,(\theta-y/\sqrt{t})^{J-1}+(\theta-y/\sqrt{t})^{J-2}\theta+\ldots+\theta^{J-1})^{t}.$

The integral reads

$\lim_{t\to\infty}\frac{1}{2}\int_{0}^{\theta\sqrt{t}}\varphi\left(y\frac{({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}\widetilde{\boldsymbol{\eta}}}{\sqrt{V(\theta-y/\sqrt{t})}}\right)\frac{y^{2}[({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}\widetilde{\boldsymbol{\eta}}]^{2}}{\sqrt{V(\theta-y/\sqrt{t})}}f(\theta-y/\sqrt{t})dy.$

The limit is

$\frac{1}{2}\int_{0}^{\infty}\varphi\left(\frac{y({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\zeta}}}{\sqrt{V(\theta)}}\right)\frac{[y({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\zeta}}]^{2}}{\sqrt{V(\theta)}}f(\theta)dy.$

Substitution of $y$ with
$z=\frac{y({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\zeta}}}{\sqrt{V(\theta)}}$
yields

$\frac{f(\theta)V(\theta)}{2|({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\zeta}}|}\int_{0}^{\infty}\varphi\left(z\right)z^{2}dz=\frac{f(\theta)V(\theta)}{4|({\boldsymbol{\beta}_{2}}-{\boldsymbol{\beta}_{1}})^{t}{\boldsymbol{\zeta}}|},$

which completes the proof. ∎

For the proofs of Theorems 4.1 and 4.2 we need the lemma below and some
notation. For every $m$ and $k$, let
$I_{m,k}:=\int_{\theta_{m-1}}^{\theta_{m}}P\left(\widehat{g}_{k}(x)>\max_{\ell\neq
k}\widehat{g}_{\ell}(x)\right)\left[g_{m}({x})-g_{k}({x})\right]f({x})d{x}$.
We then have $R(\pi_{1},\ldots,\pi_{K})=\sum_{k=1}^{K}\sum_{m=1}^{K}I_{m,k}$.

######  Lemma .2.

Under the assumptions in the beginning of Section 4.1, the integral $I_{m,k}$
is exponentially small for $k\neq m-1,m+1$, and for every $\varepsilon>0$

$\displaystyle\lim_{n\to\infty}n^{3/2-\varepsilon}\left|I_{m,m-1}-\int_{\theta_{m-1}}^{\theta_{m-1}+\frac{c_{n}}{\sqrt{n}}}\Phi\left(\frac{-\sqrt{n}\\{g_{m}(x)-g_{m-1}(x)\\}}{\sqrt{V_{m-1}(x)}}\right)\left[g_{m}({x})-g_{m-1}({x})\right]f({x})d{x}\right|=0,$
(30)
$\displaystyle\lim_{n\to\infty}n^{3/2-\varepsilon}\left|I_{m,m+1}-\int_{\theta_{m}-\frac{c_{n}}{\sqrt{n}}}^{\theta_{m}}\Phi\left(\frac{-\sqrt{n}\\{g_{m}(x)-g_{m+1}(x)\\}}{\sqrt{V_{m}(x)}}\right)\left[g_{m}({x})-g_{m+1}({x})\right]f({x})d{x}\right|=0,$

where $c_{n}=n^{\varepsilon/5}$ and
$V_{m}(x)=(1,x)\left({\boldsymbol{\Sigma}}_{m}+{\boldsymbol{\Sigma}}_{m+1}\right)\binom{1}{x}$.
Moreover, for every $m$ and $k$

$n^{3/2-\varepsilon}\left|I_{m,k}-\int_{\theta_{m-1}}^{\theta_{m}}\left(\int\prod_{l=1,\ldots,K,l\not=k}\Phi\left(\frac{z\xi_{k}(x)+\sqrt{n}[g_{k}({x})-g_{l}({x})]}{\xi_{l}(x)}\right)\varphi(z)dz\right)\left[g_{m}({x})-g_{k}({x})\right]f({x})d{x}\right|\to
0,$

as $n\to\infty$.

#### Proof of Lemma .2

The proof of the first claim in Lemma .2 is similar to that of Theorem 3.3.
Here is a sketch. For $x\in I_{m,k}$ and bounded away from $\theta_{m-1}$ and
$\theta_{m}$, and $k\neq m$ both the quantities

$P\left(\widehat{g}_{k}(x)>\max_{\ell\neq k}\widehat{g}_{\ell}(x)\right)\text{
and
}\int\prod_{l=1,\ldots,K,l\not=k}\Phi\left(\frac{z\xi_{k}(x)+\sqrt{n}[g_{k}({x})-g_{l}({x})]}{\xi_{l}(x)}\right)\varphi(z)dz$

are exponentially small (for $k=m$ the regret is zero). The same holds true
for every $x$ and when $k\neq m-1,m+1$.

We now prove (30) concerning $I_{m,m-1}$ (where $m>1$). The other relation for
$I_{m,m+1}$ is similar. Fix $c_{n}=n^{\varepsilon/5}$ and consider $x$ such
that $x\in(\theta_{m-1},\theta_{m-1}+c_{n}/\sqrt{n})$ . For such $x$ we have
by a standard large deviations argument that

$P\left(\widehat{g}_{m-1}(x)>\max_{\ell\neq{m-1}}\widehat{g}_{\ell}(x)\right)=P(\widehat{g}_{m-1}(x)>\widehat{g}_{m}(x))+a_{n},\text{
and }\\\
\int\prod_{l=1,\ldots,K,l\not=m-1}\Phi\left(\frac{z\xi_{m-1}(x)+\sqrt{n}[g_{m-1}({x})-g_{l}({x})]}{\xi_{l}(x)}\right)\varphi(z)dz=\Phi\left(\frac{-\sqrt{n}[g_{m}(x)-g_{m-1}(x)]}{\sqrt{V_{m-1}(x)}}\right)+b_{n},$
(31)

where $a_{n}$ and $b_{n}$ are exponentially small (uniformly in $x$). By a
Berry-Esseen type bound applied to the first part of the integrand, and the
smallness of the other part for $x$ is the given range of the integral, and
the smallness of the range itself , as in (26), we obtain

$\int_{\theta_{m-1}}^{\theta_{m-1}+\frac{c_{n}}{\sqrt{n}}}\left[P(\widehat{g}_{m-1}(x)>\widehat{g}_{m}(x))-\Phi\left(\frac{-\sqrt{n}[g_{m}(x)-g_{m-1}(x)]}{\sqrt{V_{m}(x)}}\right)\right]\left[g_{m}({x})-g_{m-1}({x})\right]f({x})d{x}=o(1/n^{3/2-\varepsilon}),$

and (30) follows. This, together with (31), implies the last part of the
theorem when $k=m-1$. ∎

#### Proof of Theorem 4.1

This follows directly from the last part of Lemma .2 and the fact that
$R(\pi_{1},\ldots,\pi_{K})=\sum_{k,m=1}^{K}I_{m,k}$. Note that we can replace
that latter sum by $\sum_{m=1}^{K}I_{m,m-1}$ since $I_{m,k}$ is exponentially
small for $k\neq m-1,m+1$.

### Proof of Theorem 4.2

The first part of Lemma .2 implies that

$\sum_{k=1}^{K}\sum_{m=1}^{K}I_{m,k}=\sum_{m=1}^{K-1}\int_{\theta_{m}-\frac{c_{n}}{\sqrt{n}}}^{\theta_{m}+\frac{c_{n}}{\sqrt{n}}}\Phi\left(\frac{-\sqrt{n}|g_{m}(x)-g_{m-1}(x)|}{\sqrt{V_{m-1}(x)}}\right)\left|g_{m}({x})-g_{m-1}({x})\right|f({x})d{x}+o(1/n).$

The limit of the latter integrals can be calculated using Theorem 3.4 implying
the result. ∎

## References

  * Ballarini et al. (2018) Ballarini, N. M., G. K. Rosenkranz, T. Jaki, F. König, and M. Posch (2018).  Subgroup identification in clinical trials via the predicted individual treatment effect.  PloS one 13(10), e0205971.
  * Bastani and Bayati (2015) Bastani, H. and M. Bayati (2015).  Online decision-making with high-dimensional covariates.  Available at SSRN 2661896.
  * Chen et al. (2010) Chen, L. H., L. Goldstein, and Q.-M. Shao (2010).  Normal approximation by Stein’s method.  Springer Science & Business Media.
  * Chernoff (1953) Chernoff, H. (1953).  Locally optimal designs for estimating parameters.  Ann. Math. Statist. 24(4), 586–602.
  * Ebbeling et al. (2018) Ebbeling, C. B., H. A. Feldman, G. L. Klein, J. M. W. Wong, L. Bielak, S. K. Steltz, P. K. Luoto, R. R. Wolfe, W. W. Wong, and D. S. Ludwig (2018).  Effects of a low carbohydrate diet on energy expenditure during weight loss maintenance: randomized trial.  BMJ 363.
  * Goldenshluger and Zeevi (2013) Goldenshluger, A. and A. Zeevi (2013).  A linear response bandit problem.  Stochastic Systems 3(1), 230–261.
  * Hagberg et al. (2019) Hagberg, L., A. Winkvist, H. K. Brekke, F. Bertz, E. H. Johansson, and E. Huseinovic (2019).  Cost-effectiveness and quality of life of a diet intervention postpartum: 2-year results from a randomized controlled trial.  BMC Public Health 19(1), 38.
  * Hernán and Robins (2010) Hernán, M. A. and J. M. Robins (2010).  Causal inference.
  * Hills et al. (2014) Hills, A. P., N. Mokhtar, and N. M. Byrne (2014).  Assessment of physical activity and energy expenditure: an overview of objective measures.  Frontiers in Nutrition 1, 5.
  * Moore (1973) Moore, M. (1973).  A convex matrix function.  The American Mathematical Monthly 80(4), 408–409.
  * Pan and Zhao (2020) Pan, Y. and Y.-Q. Zhao (2020).  Improved doubly robust estimation in learning optimal individualized treatment rules.  Journal of the American Statistical Association 0(0), 1–12.
  * Qian and Murphy (2011) Qian, M. and S. A. Murphy (2011).  Performance guarantees for individualized treatment rules.  Annals of statistics 39(2), 1180.
  * Silvey (2013) Silvey, S. (2013).  Optimal design: an introduction to the theory for parameter estimation, Volume 1.  Springer Science & Business Media.
  * Tian et al. (2014) Tian, L., A. A. Alizadeh, A. J. Gentles, and R. Tibshirani (2014).  A simple method for estimating interactions between a treatment and a large number of covariates.  Journal of the American Statistical Association 109(508), 1517–1532.
  * Tropp (2012) Tropp, J. A. (2012).  User-friendly tail bounds for sums of random matrices.  Foundations of Computational Mathematics 12(4), 389–434.