# Old Phase Remnants in First Order Phase Transitions

Philip Lu<EMAIL_ADDRESS>Center for Theoretical Physics, Department of
Physics and Astronomy, Seoul National University, Seoul 08826, Korea  Kiyoharu
Kawana<EMAIL_ADDRESS>Center for Theoretical Physics, Department of Physics
and Astronomy, Seoul National University, Seoul 08826, Korea  Ke-Pan Xie
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Nebraska,
Lincoln, NE 68588, USA

###### Abstract

First order phase transitions (FOPTs) are usually described by the nucleation
and expansion of new phase bubbles in the old phase background. While the
dynamics of new phase bubbles have been extensively studied, a comprehensive
treatment of the shrinking old phase remnants remained undeveloped. We present
a novel formalism for remnant statistics in FOPTs and perform the first
analytical calculations of their distribution. By shifting to the reverse time
description, we identify the shrinking remnants with expanding old phase
bubbles, allowing a quantitative evolution and determination of the population
statistics. Our results not only provide essential input for cosmological
FOPT-induced soliton/primordial black hole formation scenarios, but can also
be readily applied to generic FOPTs.

Introduction– First order phase transitions (FOPTs) are found across
disciplines as diverse as biology Lu et al. (2021); Narayanan et al. (2017),
condensed matter physics Blume (1966); Sachdev (2011), and cosmology Sato
(1981). In the cosmological context, FOPTs are a natural consequence of many
Beyond Standard Model theories, and could play a crucial role in generating
the matter-antimatter asymmetry Kuzmin et al. (1985); Joyce et al. (1996,
1995); Cohen et al. (1993); Morrissey and Ramsey-Musolf (2012), forming dark
matter Baker and Kopp (2017); Baker et al. (2020); Chway et al. (2020); Witten
(1984); Krylov et al. (2013); Huang and Li (2017); Bai and Long (2018); Bai et
al. (2019); Hong et al. (2020); Asadi et al. (2021); Marfatia and Tseng
(2021a) and primordial black holes (PBHs) Crawford and Schramm (1982); Hawking
et al. (1982); La and Steinhardt (1989); Moss (1994); Konoplich et al. (1998,
1999); Kodama et al. (1982); Lewicki and Vaskonen (2020); Kusenko et al.
(2020); Gross et al. (2021); Baker et al. (2021a, b); Kawana and Xie (2022);
Marfatia and Tseng (2021b); Huang and Xie (2022); Liu et al. (2021);
Davoudiasl et al. (2021); Jung and Okui (2021); Hashino et al. (2021); Maeso
et al. (2021), and leave detectable signals in current or near-future
gravitational wave detectors Xue et al. (2021); Arzoumanian et al. (2020);
Romero et al. (2021); Caprini et al. (2016, 2020); Liang et al. (2022).

Cosmological FOPTs happen through the nucleation and growth of new true vacuum
(TV) phase bubbles in the old false vacuum (FV) phase background. More
attention has been focused on the calculation and estimation of the properties
of TV bubbles, whose statistics are relevant for electroweak baryogenesis
Kuzmin et al. (1985); Joyce et al. (1996, 1995); Cohen et al. (1993);
Morrissey and Ramsey-Musolf (2012), inflationary models, and for the
production of gravitational waves Hogan (1986); Maggiore (2000); Kawana
(2022). As a result, analytic methods were developed to estimate the TV
nucleation rate, wall velocity, and bubble distribution Callan and Coleman
(1977); Guth and Tye (1980); Guth and Weinberg (1981); Affleck (1981); Dine et
al. (1992); Espinosa et al. (2010); Megevand and Ramirez (2017); Ellis et al.
(2019); Wang et al. (2020); Bodeker and Moore (2009, 2017); Höche et al.
(2021); Gouttenoire et al. (2021). On the other hand, FV remnants are more
relevant for the mechanisms involving trapping particles in the FV to realize
baryogenesis Arakawa et al. (2021), dark matter Witten (1984); Krylov et al.
(2013); Huang and Li (2017); Bai and Long (2018); Bai et al. (2019); Hong et
al. (2020); Asadi et al. (2021); Marfatia and Tseng (2021a) and PBHs Baker et
al. (2021a); Kawana and Xie (2022); Baker et al. (2021b); Marfatia and Tseng
(2021b); Huang and Xie (2022). Lacking an equivalent detailed description of
the FV remnants, previous studies had to either use naive estimations of the
average remnant size and density Krylov et al. (2013); Hong et al. (2020), or
take those observables as free parameters Baker et al. (2021a, b).

Figure 1: Progression of the FOPT: with the flow of time from left to right,
the initial FV (white) space is nucleated with TV (blue) bubbles, which
populate and surround the shrinking FV bubbles. In the reversed time
description, FV bubbles nucleate and grow in a TV dominated space, eventually
resulting in shrinking TV bubbles.

In this Letter, we develop a method for calculating remnant statistics from
first principles. In the existing framework, the FV phase acquires a decay
probability to the energetically favorable TV phase below the critical
temperature. The vacuum pressure causes the TV bubbles to expand, filling up
the space and leaving shrinking pockets of disconnected FV remnants (as
illustrated in Fig. 1). We extend the TV bubble nucleation formalism to
include FV bubble nucleation by considering the phase transition in reverse.
From this reverse time description, the centers of the collapsing remnants
with time flowing forward can be viewed as the nucleation sites of FV bubbles
with time flowing backwards. Thus, the methods used for TV bubbles can be
adapted for the dynamics of FV bubbles. We perform the first calculation of
the properties and evolution of the FV bubbles based on the reverse time
description of the FOPT. Although we present our results within a cosmological
context, they can be easily adapted to general FOPTs in fields such as
condensed matter physics and biology.

True vacuum bubbles– We review the existing TV nucleation formalism for
cosmological phase transitions. We assume thin walls and constant velocity
$v_{w}$ throughout, which is a good approximation for a range of moderately
strong FOPTs Megevand and Ramirez (2017). Consider a Universe initially in an
FV phase at high temperatures. Below the critical temperature $T_{c}$, the TV
phase becomes energetically favorable, giving a non-zero probability for the
FV space to tunnel to the TV. The TV nucleation rate per unit volume and unit
time is Affleck (1981)

$\Gamma(t)=A(t)e^{-S(t)},$ (1)

where $S(t)$ is the smaller of the two instanton bounce actions $S_{3}/T$
Linde (1983) and $S_{4}$ Coleman (1977) ($S_{3}/T$ is typically smaller
Megevand and Ramirez (2017)).

Assuming the bubbles grow spherically outwards, the radius of a bubble
nucleated at time $t^{\prime}$ is

$R(t,t^{\prime})=v_{w}\int_{t^{\prime}}^{t}dt^{\prime\prime}\frac{a(t)}{a(t^{\prime\prime})}~{},$
(2)

where $a(t)$ is the scale factor of the FLRW metric. The fraction of space in
the FV is Guth and Tye (1980); Guth and Weinberg (1981)

$f_{\rm fv}(t)=e^{-I(t)}~{},$ (3)

where

$I(t)=\int_{t_{c}}^{t}dt^{\prime}\Gamma(t^{\prime})\frac{a^{3}(t^{\prime})}{a^{3}(t)}\frac{4\pi}{3}R^{3}(t,t^{\prime})~{},$
(4)

with $t_{c}$ being the cosmic time corresponding to $T_{c}$. It will be
convenient to define the average TV nucleation rate, $\bar{\Gamma}(t)=f_{\rm
fv}(t)\Gamma(t)$, which contains a factor of $f_{\rm fv}$ since TV bubbles can
only nucleate in patches of FV.

The distribution of TV bubbles at time $t$ with size $R$ must equal the
average nucleation rate at a time $t^{\prime}$ at which $R(t,t^{\prime})=R$
(when the scale factor is negligible, $t^{\prime}=t-R/v_{w}$),

$\frac{dn_{\rm
tv}}{dR}(t)=\frac{1}{v_{w}}\bar{\Gamma}(t^{\prime})\frac{a^{4}(t^{\prime})}{a^{4}(t)}~{}.$
(5)

Integrating this equation, the total bubble number density is given by

$n_{\rm
tv}(t)=\int_{t_{c}}^{t}dt^{\prime}\bar{\Gamma}(t^{\prime})\frac{a^{3}(t^{\prime})}{a^{3}(t)}~{},$
(6)

where we have used Eq. (2).

False vacuum bubbles– In the forward evolution of the phase transition,
spherical TV bubbles percolate when the FV volume fraction drops below $f_{\rm
fv}=0.71$ Rintoul and Torquato (1997a), forming an infinite connected cluster.
As $f_{\rm fv}$ decreases, FV regions are separated into shrinking remnants.
When the process is considered in reverse, FV bubbles nucleate in a TV
background, forming an infinite connected cluster above $f_{\rm fv}\approx
1-0.71=0.29$. This is the reverse time description of the phase transition, in
which we identify shrinking remnants forward in time with nucleating bubbles
backward in time.

We develop an analogous formalism for FV bubbles in this reverse time
description and calculate a FV bubble nucleation rate. To do so, the
nucleation point of a FV bubble can be identified as the projected center of a
collapsing remnant, with the nucleation rate equal to the collapse rate. In
the final stages of collapse, the size of the remnant is much smaller than the
radii of the enveloping TV bubbles, so that the collapsing walls can be
treated as flat or planar.

The projection interpretation, in which we extrapolate the wall trajectories
to the point of collapse, suggests a counterpart to the reverse time
description (the projection interpretation is extensively discussed in the
last part of the Letter). In the reverse time description, the FV bubbles
nucleated at point $C$ and time $t^{\prime}>t$ grow to a size
$R_{r}(t,t^{\prime})$ at time $t$. In the alternative forward time
description, the FV remnants of size $R_{r}(t,t^{\prime})$ at time $t$ are
those which are projected to collapse at point $C$ and time $t^{\prime}>t$.
These two viewpoints are complementary and fundamentally equivalent.

In one spatial dimension, i.e. on a line, it takes two TV walls to close a FV
remnant at a point. For a final collapse at time $t$, there needs to be
nucleation at some time $t^{\prime}<t$ a distance $v_{w}(t-t^{\prime})$ away
(neglecting scale factors), and another nucleation at a time
$t^{\prime\prime}<t^{\prime}$ at a distance $v_{w}(t-t^{\prime\prime})$ away
and opposite of the first nucleation. Moving to higher dimensions, the only
form of collapse with non-vanishing probability in two (three) dimensions is
by three (four) walls forming a triangle (tetrahedron).

The probability of FV remnant collapse per unit volume per unit time,
$d^{2}P_{r}/dVdt$, at a point $C$ can be found by integrating along the past
wall cone (similar to the light cone but with wall velocity $v_{w}$) and
finding the four points of TV nucleation, so that the four resulting walls
would meet at point $C$ and time $t$. We order the TV walls by their proximity
in time to $t$, with wall 1 being the most recently nucleated wall and thereby
avoid unnecessary combinatoric factors. The nucleation point of wall 1 is then
only constrained to lie on the past wall cone of point $C$, but subsequent
walls have to obey angular restrictions to form a closed FV remnant. The
angular factors become more complicated with higher spatial dimensions, so we
first introduce the formalism in two dimensions before developing the three
dimensional case.

Two spatial dimensions– For simplicity, we omit scale factors in this
exposition and restore them in the full three-dimensional expression Eq. (14).
In its final stage, the infinitesimal shrinking FV remnant is a triangle
collapsing towards its incenter surrounded by three TV walls. The collapse
probability per unit area per unit time, $d^{2}P_{r}/dAdt$ can be found by
integrating along the past wall cone of the collapse point $C$. The radial and
angular integrations can be separated, and the radial part is

$v_{w}f_{\rm
fv}(t)\int_{t_{c}}^{t}dt_{1}v_{w}(t-t_{1})\Gamma(t_{1})\int_{t_{c}}^{t_{1}}dt_{2}v_{w}(t-t_{2})\Gamma(t_{2})\\\
\times\int_{t_{c}}^{t_{2}}dt_{3}v_{w}(t-t_{3})\Gamma(t_{3})~{},$ (7)

with the integration times ordered as $t_{c}<t_{3}<t_{2}<t_{1}<t$,
corresponding to the nucleation time of each successive wall. The factors of
$\Gamma$ gives the TV nucleation rate of each wall, and the factor of the FV
fraction, $f_{\rm fv}$, is required since only FV points are eligible to
collapse to TV. Since we are integrating along the past wall cone, the
integration space is in the FV and no additional factors of $f_{\rm fv}$ show
up inside the integrals. Otherwise, the TV would have already spread to point
$C$ before time $t$.

Figure 2: Geometry of a collapsing FV remnant. Left: In 2-dimensional space,
to form a closed triangle, $-\hat{\textbf{r}}_{3}$ must lie in the angular
range bounded by $\hat{\textbf{r}}_{1}$ and $\hat{\textbf{r}}_{2}$. Right: In
3-dimensional space, to form a closed tetrahedron, $-\hat{\textbf{r}}_{4}$
must lie in the solid angle delimited by $\hat{\textbf{r}}_{1}$,
$\hat{\textbf{r}}_{2}$ and $\hat{\textbf{r}}_{3}$ (the four walls are omitted
for clarity).

For the angular factor, we denote the normal vector of wall $i$ by
$\hat{\textbf{r}}_{i}$, i.e. the unit vector pointing from the TV nucleation
point towards $C$. For the three walls to form a closed triangle,
$-\hat{\textbf{r}}_{3}$ must lie in the angular range bounded by
$\hat{\textbf{r}}_{1}$ and $\hat{\textbf{r}}_{2}$, as illustrated in the left
panel of Fig. 2. To integrate over all the allowed TV bubble configurations,
we first choose $\hat{\textbf{r}}_{1}$ along the $\hat{x}$ axis and
parameterize the other two normal vectors as

$\hat{\textbf{r}}_{2}=(\cos\phi_{2},\sin\phi_{2}),\quad\hat{\textbf{r}}_{3}=(-\cos\phi_{3},-\sin\phi_{3})~{}.$
(8)

Then for $0<\phi_{2}<\pi$, the closure of the FV triangle requires
$0<\phi_{3}<\phi_{2}$; while the case $\pi<\phi_{2}<2\pi$ is just a reflection
of the case $0<\phi_{2}<\pi$ across the $\hat{x}$ axis. Therefore, the angular
integral reads

$2\pi\int_{0}^{\pi}2d\phi_{2}\int_{0}^{\phi_{2}}d\phi_{3}=2\pi^{3}~{},$ (9)

where the first “$2\pi$” factor represents the integral over the arbitrary
$\hat{\textbf{r}}_{1}$ angle, while the “2” factor in $\phi_{2}$ integral
accounts for the $\pi<\phi_{2}<2\pi$ region.

In the reverse time description, the FV bubble nucleation rate is defined as
the remnant collapse probability per unit area per unit time, $\Gamma_{r}^{\rm
2d}(t)=d^{2}P_{r}/dAdt$. Combining Eqs. (7) and (9), the FV bubble nucleation
rate in two dimensions is

$\Gamma_{r}^{\rm 2d}(t)=2\pi^{3}v_{w}^{4}f_{\rm
fv}(t)\int_{t_{c}}^{t}dt_{1}(t-t_{1})\Gamma(t_{1})\times\\\
\int_{t_{c}}^{t_{1}}dt_{2}(t-t_{2})\Gamma(t_{2})\int_{t_{c}}^{t_{2}}dt_{3}(t-t_{3})\Gamma(t_{3})~{}.$
(10)

Three spatial dimensions– Next we build on our two dimensional results and
apply the method to three dimensions. Integrating over the past wall cone, the
radial factor is similar, but the angular factor is more complicated, as the
three-dimensional remnant is now tetrahedral with four collapsing walls.
Ordering the walls temporally and labeling the normal vectors as before, the
condition that the four walls form a tetrahedron is that
$-\hat{\textbf{r}}_{4}$ should lie in the solid angle delimited by
$\hat{\textbf{r}}_{1}$, $\hat{\textbf{r}}_{2}$ and $\hat{\textbf{r}}_{3}$, as
illustrated in the right panel of Fig. 2.

With the collapse point $C$ at the origin, we set the $\hat{z}$ axis in the
direction of $\hat{\textbf{r}}_{1}$ and $\hat{y}$ axis along the
$\hat{\textbf{r}}_{1}\times\hat{\textbf{r}}_{2}$ direction, defining our
spherical coordinate system. The normal vectors can then be parameterized as

$\begin{split}&\hat{\textbf{r}}_{2}=(\sin\theta_{2},0,\cos\theta_{2}),\\\
&\hat{\textbf{r}}_{3}=(\sin\theta_{3}\cos\phi_{3},\sin\theta_{3}\sin\phi_{3},\cos\theta_{3}),\\\
&\hat{\textbf{r}}_{4}=(-\sin\theta_{4}\cos\phi_{4},-\sin\theta_{4}\sin\phi_{4},-\cos\theta_{4})~{}.\end{split}$
(11)

When $0<\phi_{3}<\pi$, the closure condition requires $0<\phi_{4}<\phi_{3}$
and $0<\theta_{4}<\bar{\theta}(\phi_{4})$, where

$\bar{\theta}(\phi_{4})={\rm arccot}\Big{(}\cot\theta_{2}\cos\phi_{4}\\\
+\frac{\cot\theta_{3}-\cot\theta_{2}\cos\phi_{3}}{\sin\phi_{3}}\sin\phi_{4}\Big{)}~{}$
(12)

is determined by the intersection of plane $C\hat{\textbf{r}}_{4}\hat{z}$
[$(\sin\phi_{4})x-(\cos\phi_{4})y=0$] and plane
$C\hat{\textbf{r}}_{2}\hat{\textbf{r}}_{3}$
[$(\cot\theta_{2})x+(\sin\phi_{3})^{-1}(\cot\theta_{3}-\cot\theta_{2}\cos\phi_{3})y-z=0$].
Here the range of the arccotangent function is limited to $(0,\pi)$, as the
case $\pi<\phi_{3}<2\pi$ is just a reflection of the case $0<\phi_{3}<\pi$
over the $C\hat{x}\hat{z}$ plane. Therefore, the angular integral reads

$4\pi\int_{0}^{\pi}d\theta_{2}\int_{0}^{\pi}d\theta_{3}\int_{0}^{\bar{\theta}(\phi_{4})}d\theta_{4}\int_{0}^{2\pi}d\phi_{2}\times\\\
\int_{0}^{\pi}2d\phi_{3}\int_{0}^{\phi_{3}}d\phi_{4}\sin\theta_{2}\sin\theta_{3}\sin\theta_{4}=32\pi^{4}~{},$
(13)

where the first “$4\pi$” factor represents the integral over the arbitrary
$\hat{\textbf{r}}_{1}$, while the factor of “2” in the $\phi_{3}$ integral
accounts for the $\pi<\phi_{3}<2\pi$ region. A more intuitive derivation of
Eq. (13) can be found in the Appendix.

With Eq. (13) in hand, the FV bubble nucleation rate at time $t$ is

$\Gamma_{r}(t)\equiv\frac{d^{2}P_{r}}{dVdt}=32\pi^{4}v_{w}^{9}f_{\rm
fv}(t)\int_{t_{c}}^{t}dt_{1}\left(\int_{t_{1}}^{t}dt_{1}^{\prime}\frac{a(t_{1})}{a(t_{1}^{\prime})}\right)^{2}\Gamma(t_{1})\frac{a(t_{1})}{a(t)}\int_{t_{c}}^{t_{1}}dt_{2}\left(\int_{t_{2}}^{t}dt_{2}^{\prime}\frac{a(t_{2})}{a(t_{2}^{\prime})}\right)^{2}\Gamma(t_{2})\frac{a(t_{2})}{a(t)}\\\
\times\int_{t_{c}}^{t_{2}}dt_{3}\left(\int_{t_{3}}^{t}dt_{3}^{\prime}\frac{a(t_{3})}{a(t_{3}^{\prime})}\right)^{2}\Gamma(t_{3})\frac{a(t_{3})}{a(t)}\int_{t_{c}}^{t_{3}}dt_{4}\left(\int_{t_{4}}^{t}dt^{\prime}_{4}\frac{a(t_{4})}{a(t^{\prime}_{4})}\right)^{2}\Gamma(t_{4})\frac{a(t_{4})}{a(t)}~{},$
(14)

where the integrals over the scale factors come from the integration element
of the radial direction in spherical coordinates. The integral simplifies to

$\left(\int_{t_{i}}^{t}dt_{i}^{\prime}\frac{a(t_{i})}{a(t_{i}^{\prime})}\right)^{2}\to(t-t_{i})^{2}~{},$
(15)

when the scale factors can be taken as constant.

False vacuum bubble distribution– In the reverse time description, we can use
the FV bubble nucleation rate, Eq. (14), to find the FV bubble distribution.
Analogously to the TV bubble case, a FV bubble which nucleates at
$t^{\prime}>t$ has a radius

$R_{r}(t,t^{\prime})=v_{w}\int_{t}^{t^{\prime}}dt^{\prime\prime}\frac{a(t)}{a(t^{\prime\prime})}~{},$
(16)

at time $t$ (note that FV bubbles are nucleated in the future).

Since FV bubbles can only nucleate in patches of TV, the average FV bubble
nucleation rate is $\bar{\Gamma}_{r}(t)=(1-f_{\rm fv}(t))\Gamma_{r}(t)$. The
FV bubble size distribution is then

$\frac{n_{\rm
fv}}{dR_{r}}(t)=\frac{1}{v_{w}}\bar{\Gamma}_{r}(t^{\prime})\frac{a^{4}(t^{\prime})}{a^{4}(t)}~{},$
(17)

where $t^{\prime}$ is resolved using Eq. (16), reducing to
$t^{\prime}=t+R_{r}/v_{w}$ in the limit of constant scale factor. Integrating
Eq. (17) over $R_{r}$ yields the overall FV bubble number density

$n_{\rm
fv}(t)=\int_{t}^{t_{e}}dt^{\prime}\bar{\Gamma}_{r}(t^{\prime})\frac{a^{3}(t^{\prime})}{a^{3}(t)}~{},$
(18)

where $t_{e}$ is the ending time of the FOPT, which can be effectively taken
as $+\infty$ and $f_{\rm fv}(t_{e})=0$.

Exponential nucleation– So far our results are rather general and apply to any
FOPT scenario as long as the TV nucleation rate $\Gamma(t)$ is available. We
now evaluate our results in the exponential nucleation rate approximation,

$\Gamma(t)=\Gamma_{\ast}e^{\beta(t-t_{\ast})}~{},$ (19)

expanded around an arbitrary time $t_{*}\in(t_{c},t_{e})$, where
$\beta=-dS(t)/dt|_{t_{*}}$ can be treated as the inverse of the FOPT duration.
This exponential approximation is accurate if the FOPT proceeds rapidly
compared to the Hubble time scale, i.e. $\beta/H(t_{*})\gg 1$. Hence, the
scale factor $a(t)$ is approximately constant over the transition and can be
neglected.

There are many motivated choices for the time $t_{\ast}$ in the exponential
case such as the nucleation time ($n_{\rm tv}H^{-3}=1$) or TV bubble
percolation time ($f_{\rm fv}=0.71$) Holcomb and Rehr (1969); Rintoul and
Torquato (1997b); Li et al. (2020), etc. The exponential approximation fits
numerical simulations the best when $t_{\ast}$ is chosen to be the time at
which the bubble statistics are computed. For the remnant distribution, we
choose $t_{*}$ to be the FV bubble percolation time, defined such that $f_{\rm
fv}(t_{\ast})=0.29$ and $I_{\ast}=-\ln(0.29)=1.238$.

We can explicitly solve for the FV filling fraction, Eq. (3) with

$I(t)=\frac{8\pi v_{w}^{3}\Gamma_{\ast}}{\beta^{4}}e^{\beta(t-t_{\ast})}\equiv
I_{\ast}e^{\beta(t-t_{\ast})}~{},$ (20)

where $\beta(t_{*}-t_{c})\gg 1$ is used. Integrating Eq. (14), the average FV
bubble nucleation rate is

$\bar{\Gamma}_{r}(t)\approx\frac{I_{*}^{4}\beta^{4}}{192v_{w}^{3}}e^{4\beta(t-t_{*})}e^{-I(t)}\left(1-e^{-I(t)}\right)~{}.$
(21)

This FV nucleation rate has the necessary properties of a solution. The
factors of $f_{\rm fv}=e^{-I}$ and $(1-f_{\rm fv})=(1-e^{-I})$ constrain the
majority of FV bubbles to nucleate during the phase transition,
$t_{c}<t<t_{e}$, and the factor of $e^{4\beta(t-t_{*})}$ biases the FV
nucleation towards later times, $f_{\rm fv}\ll 1$, which is to be expected of
collapsing FV remnants. By using Eq. (17), the FV bubble distribution is found
to be

$\displaystyle\frac{dn_{\rm fv}}{dR_{r}}(t_{*})$
$\displaystyle\sim\frac{1}{v_{w}}\Gamma_{r}\left(t_{*}+R/v_{w}\right)$ (22)
$\displaystyle=\frac{I_{*}^{4}\beta^{4}}{192v_{w}^{3}}e^{4\beta
R_{r}/v_{w}}e^{-I_{*}e^{\beta R_{r}/v_{w}}}\left(1-e^{-I_{*}e^{\beta
R_{r}/v_{w}}}\right)~{}.$

False vacuum bubble normalization– We determine the approximate shape of the
FV bubbles by normalizing the total volume contained in the FV bubbles to the
filling fraction $f_{\rm fv}$ at the remnant percolation time $t_{*}$. While
the projected collapse is tetrahedral, which would suggest a volume of
$8\sqrt{3}R^{3}$ for a regular tetrahedron, with $R$ the inradius, the FV
bubble is expected to become more rounded as additional TV walls partially
cover the bubble. This mimics the effects of tension, which is neglected in
this treatment, but which tends to smooth the bubble edges. We expect FV
bubbles at the percolation threshold to be roughly spherical but abnormally
shaped, as depicted in Fig. 1. We therefore use the volume formula
$A(4\pi/3)R^{3}$, where $(A-1)$ measures the departure from sphericality.
Summing up the FV bubble volumes and normalizing to the FV fraction,

$\int dR_{r}A\frac{4\pi R_{r}^{3}}{3}\frac{dn_{\rm
fv}}{dR_{r}}(t_{*})=0.29~{},$ (23)

yields $A=1.15$. This suggests that the FV bubbles/remnants are approximately
spherical but with a non-negligible deviation. From another point of view, the
factor $A=1.15$ is the required normalization for the formalism to be self-
consistent. Although overlap between adjacent FV bubbles is ignored here, the
volume $V\approx 1.15(4\pi/3)R^{3}$ is the average volume that effectively
“belongs” to each collapsing remnant of size $R$ at the percolation time.

Note that the value of the FV percolation chosen here, $f_{\rm fv}=0.29$, is
strictly valid only for spherical bubbles of equal size. Since the FV bubble
distribution is not monochromatic, nor are they spherical, the exact
percolation threshold would have to be determined by a numerical simulation.
We do not expect the threshold to deviate significantly from $f_{\rm
fv}=0.29$.

Application to Primordial Black Holes– As a concrete example, we apply our
method to the Fermi-ball/PBH formation scenario proposed in Refs. Hong et al.
(2020); Kawana and Xie (2022) and subsequently studied in Refs. Marfatia and
Tseng (2021a, b); Huang and Xie (2022). During the FOPT, an asymmetric
population of dark fermions $\chi$-$\overline{\chi}$ is trapped between the
expanding TV bubble walls into the collapsing FV remnants due to a large mass
differential for $\chi$ in the two phases. As the remnants shrink, the
fermions and anti-fermions annihilate, leaving only the asymmetrical portion
supported by degeneracy pressure. The total number of $\chi$ fermions trapped
in a remnant with size $R_{*}$ at the FV percolation time $t_{*}$ is Hong et
al. (2020); Kawana and Xie (2022)

$Q_{\rm FB}=\frac{\eta_{\chi}s(t_{*})}{f_{\rm fv}(t_{*})}A\frac{4\pi
R_{*}^{3}}{3}~{},$ (24)

where $\eta_{\chi}=(n_{\chi}-n_{\bar{\chi}})/s$ is the $\chi$-asymmetry with
$s$ the entropy density.

Figure 3: A benchmark $f(M_{\rm PBH})$ distribution (red curve) for the Fermi-
ball/PBH formation scenario of Ref. Kawana and Xie (2022), where $f$ is the
PBH mass distribution probability function. The FOPT parameters are shown in
the figure, with $U_{0}$ the energy density difference between the TV and FV.
The vertical black thin line and blue dashed line represent $\left\langle
M_{\rm PBH}\right\rangle$ from our work and the estimate of $M_{\rm PBH}$ from
Ref. Kawana and Xie (2022), respectively.

Since the Fermi-ball/PBH mass $M_{\rm PBH}\propto Q_{\rm FB}$ Hong et al.
(2020), the distribution of $R_{*}$ is key to deriving the Fermi-ball/PBH mass
profile. Lacking methods to compute the $R_{*}$ distribution, Refs. Hong et
al. (2020); Kawana and Xie (2022) estimated the average size to be

$R_{*}\sim\left(\frac{3v_{w}}{4\pi
A\Gamma_{*}}\right)^{1/4}=1.43\frac{v_{w}}{\beta}~{},$ (25)

resulting in a monochromatic Fermi-ball/PBH mass distribution. With the
technique developed in this Letter, we are able to derive the $R_{*}$
distribution using Eq. (22), and obtain

$\left\langle R\right\rangle_{*}=\int dR_{r}R_{r}\frac{dn_{\rm
fv}}{dR_{r}}(t_{*})=1.12\frac{v_{w}}{\beta}~{}.$ (26)

Although this seems quite close to the estimated value, a continuous $R_{*}$
distribution results in an extended Fermi-ball/PBH mass profile, which greatly
impacts experimental constraints Green and Kavanagh (2021). In Fig. 3, we
compute an example PBH distribution at present time which comprises all of
dark matter and lie within the PBH mass window, $10^{17}~{}{\rm
g}<M<10^{23}~{}{\rm g}$, where there are no reliable bounds on PBH dark matter
Carr et al. (2021).

Discussion– In our derivation of the FV bubble nucleation rate Eq. (14), we
made a few simplifying assumptions. First, the walls were assumed to be
infinitesimally thin and the wall velocity constant. TV bubble mergers can
alter the effective location of the TV wall nucleation and expansion
direction. Furthermore, surface tension tends to shape the collapsing remnants
to be more spherical, whereas the collapsing remnant is treated as tetrahedral
with four planar walls. All of these effects are exacerbated near the end of
the phase transition when TV bubbles inevitably merge, surface tension becomes
more prominent as the surface area to volume ratio increases, and particle
trapping may stop or slow the collapse.

We resolve these issues by interpreting the method and formalism of computing
$\Gamma_{r}(t)$ developed in this Letter as a projected description of the
FOPT rather than a physical description. In other words, beyond the remnant
percolation time $t_{*}$ at which we evaluate Eq. (22), the future evolution
of the FOPT is irrelevant to the remnant statistics evaluated at $t_{*}$.
Causally, the remnant size distribution and number density at $t_{*}$ cannot
be affected by events at later times $t>t_{*}$. Thus, the remnant statistics
at the remnant percolation time depend only on the history of the FOPT up to
the that time, before which these four assumptions are only mildly violated.
The calculated observables will be the same whether, in the later stages of
the transition, our idealized collapse scenario is applied or a more physical
scenario with surface tension effects is used instead. Hence, the projection
interpretation is that our formalism essentially traces the walls of the
collapsing remnants forward past time $t_{*}$ to find the collapse point, in
order to then trace the collapse backwards in the reverse time description and
infer the size of the remnant at time $t_{*}$.

We offer an analogy: the shadow of a falling apple can be used to infer its
instantaneous position. Whether or not the apple eventually lands directly on
its shadow, is perturbed by a gust of wind, or drops on an unwitting head, is
immaterial to the determination of its instantaneous position. Likewise,
whether the remnant eventually collapses spherically or is stopped by
degeneracy pressure is irrelevant to its size distribution at remnant
percolation. Therefore, in the derivation of the remnant distribution, our
tetrahedral collapse model is more appropriate than a physically realistic
spherical collapse model because the shadow (or projection) of the collapsing
walls is tetrahedral and not spherical.

In summary, we have performed the first analytical calculation of the FV
remnant distribution and evolution in FOPTs. By identifying the center of a
collapsing remnant as a FV bubble nucleation in reverse, the established TV
bubble nucleation formalism can be adapted to derive remnant statistics. Our
results provide a more sophisticated treatment of FOPT remnants than
previously available, and are directly applicable to many new physics
mechanisms involving trapped particles Arakawa et al. (2021); Witten (1984);
Krylov et al. (2013); Huang and Li (2017); Bai and Long (2018); Bai et al.
(2019); Hong et al. (2020); Asadi et al. (2021); Marfatia and Tseng (2021a);
Baker et al. (2021a); Kawana and Xie (2022); Baker et al. (2021b); Marfatia
and Tseng (2021b); Huang and Xie (2022). The novel formalism developed in this
Letter can be readily generalized and applied outside the cosmological
context.

Acknowledgments.– We thank Sunghoon Jung, Hyung Do Kim, Taehun Kim, and Howard
Chen for useful discussions. The work of P.L. was supported by the Grant Korea
NRF-2019R1C1C1010050. K.P.X. is supported by the University of Nebraska-
Lincoln. The work of K.K. is supported by Grant Korea NRF-2019R1C1CC1010050,
2019R1A6A1A10073437.

## References

  * Lu et al. (2021) Y. Lu, Y. Lu, J. Ma, J. Li, X. Huang, Q. Jia, D. Ma, M. Liu, H. Zhang, X. Yu, et al., bioRxiv (2021), URL https://www.biorxiv.org/content/early/2021/04/20/2021.04.15.439963.
  * Narayanan et al. (2017) A. Narayanan, A. B. Meriin, M. Y. Sherman, and I. I. Cissé, bioRxiv (2017), URL https://www.biorxiv.org/content/early/2017/06/19/148395.
  * Blume (1966) M. Blume, Physical Review 141, 517 (1966).
  * Sachdev (2011) S. Sachdev, _Quantum Phase Transitions_ (2011).
  * Sato (1981) K. Sato, Mon. Not. Roy. Astron. Soc. 195, 467 (1981).
  * Kuzmin et al. (1985) V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov, Phys. Lett. B 155, 36 (1985).
  * Joyce et al. (1996) M. Joyce, T. Prokopec, and N. Turok, Phys. Rev. D 53, 2958 (1996), eprint hep-ph/9410282.
  * Joyce et al. (1995) M. Joyce, T. Prokopec, and N. Turok, Phys. Rev. Lett. 75, 1695 (1995), [Erratum: Phys.Rev.Lett. 75, 3375 (1995)], eprint hep-ph/9408339.
  * Cohen et al. (1993) A. G. Cohen, D. B. Kaplan, and A. E. Nelson, Ann. Rev. Nucl. Part. Sci. 43, 27 (1993), eprint hep-ph/9302210.
  * Morrissey and Ramsey-Musolf (2012) D. E. Morrissey and M. J. Ramsey-Musolf, New J. Phys. 14, 125003 (2012), eprint 1206.2942.
  * Baker and Kopp (2017) M. J. Baker and J. Kopp, Phys. Rev. Lett. 119, 061801 (2017), eprint 1608.07578.
  * Baker et al. (2020) M. J. Baker, J. Kopp, and A. J. Long, Phys. Rev. Lett. 125, 151102 (2020), eprint 1912.02830.
  * Chway et al. (2020) D. Chway, T. H. Jung, and C. S. Shin, Phys. Rev. D 101, 095019 (2020), eprint 1912.04238.
  * Witten (1984) E. Witten, Phys. Rev. D 30, 272 (1984).
  * Krylov et al. (2013) E. Krylov, A. Levin, and V. Rubakov, Phys. Rev. D 87, 083528 (2013), eprint 1301.0354.
  * Huang and Li (2017) F. P. Huang and C. S. Li, Phys. Rev. D 96, 095028 (2017), eprint 1709.09691.
  * Bai and Long (2018) Y. Bai and A. J. Long, JHEP 06, 072 (2018), eprint 1804.10249.
  * Bai et al. (2019) Y. Bai, A. J. Long, and S. Lu, Phys. Rev. D 99, 055047 (2019), eprint 1810.04360.
  * Hong et al. (2020) J.-P. Hong, S. Jung, and K.-P. Xie, Phys. Rev. D 102, 075028 (2020), eprint 2008.04430.
  * Asadi et al. (2021) P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer, and J. Smirnov, Phys. Rev. Lett. 127, 211101 (2021), eprint 2103.09822.
  * Marfatia and Tseng (2021a) D. Marfatia and P.-Y. Tseng, JHEP 11, 068 (2021a), eprint 2107.00859.
  * Crawford and Schramm (1982) M. Crawford and D. N. Schramm, Nature 298, 538 (1982).
  * Hawking et al. (1982) S. W. Hawking, I. G. Moss, and J. M. Stewart, Phys. Rev. D 26, 2681 (1982).
  * La and Steinhardt (1989) D. La and P. J. Steinhardt, Phys. Lett. B 220, 375 (1989).
  * Moss (1994) I. G. Moss, Phys. Rev. D 50, 676 (1994).
  * Konoplich et al. (1998) R. Konoplich, S. Rubin, A. Sakharov, and M. Y. Khlopov, Astronomy Letters 24, 413 (1998).
  * Konoplich et al. (1999) R. V. Konoplich, S. G. Rubin, A. S. Sakharov, and M. Y. Khlopov, Phys. Atom. Nucl. 62, 1593 (1999).
  * Kodama et al. (1982) H. Kodama, M. Sasaki, and K. Sato, Prog. Theor. Phys. 68, 1979 (1982).
  * Lewicki and Vaskonen (2020) M. Lewicki and V. Vaskonen, Phys. Dark Univ. 30, 100672 (2020), eprint 1912.00997.
  * Kusenko et al. (2020) A. Kusenko, M. Sasaki, S. Sugiyama, M. Takada, V. Takhistov, and E. Vitagliano, Phys. Rev. Lett. 125, 181304 (2020), eprint 2001.09160.
  * Gross et al. (2021) C. Gross, G. Landini, A. Strumia, and D. Teresi, JHEP 09, 033 (2021), eprint 2105.02840.
  * Baker et al. (2021a) M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht (2021a), eprint 2105.07481.
  * Baker et al. (2021b) M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht (2021b), eprint 2110.00005.
  * Kawana and Xie (2022) K. Kawana and K.-P. Xie, Phys. Lett. B 824, 136791 (2022), eprint 2106.00111.
  * Marfatia and Tseng (2021b) D. Marfatia and P.-Y. Tseng (2021b), eprint 2112.14588.
  * Huang and Xie (2022) P. Huang and K.-P. Xie (2022), eprint 2201.07243.
  * Liu et al. (2021) J. Liu, L. Bian, R.-G. Cai, Z.-K. Guo, and S.-J. Wang (2021), eprint 2106.05637.
  * Davoudiasl et al. (2021) H. Davoudiasl, P. B. Denton, and J. Gehrlein (2021), eprint 2109.01678.
  * Jung and Okui (2021) T. H. Jung and T. Okui (2021), eprint 2110.04271.
  * Hashino et al. (2021) K. Hashino, S. Kanemura, and T. Takahashi (2021), eprint 2111.13099.
  * Maeso et al. (2021) D. N. Maeso, L. Marzola, M. Raidal, V. Vaskonen, and H. Veermäe (2021), eprint 2112.01505.
  * Xue et al. (2021) X. Xue et al., Phys. Rev. Lett. 127, 251303 (2021), eprint 2110.03096.
  * Arzoumanian et al. (2020) Z. Arzoumanian et al. (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), eprint 2009.04496.
  * Romero et al. (2021) A. Romero, K. Martinovic, T. A. Callister, H.-K. Guo, M. Martínez, M. Sakellariadou, F.-W. Yang, and Y. Zhao, Phys. Rev. Lett. 126, 151301 (2021), eprint 2102.01714.
  * Caprini et al. (2016) C. Caprini et al., JCAP 04, 001 (2016), eprint 1512.06239.
  * Caprini et al. (2020) C. Caprini et al., JCAP 03, 024 (2020), eprint 1910.13125.
  * Liang et al. (2022) Z.-C. Liang, Y.-M. Hu, Y. Jiang, J. Cheng, J.-d. Zhang, and J. Mei, Phys. Rev. D 105, 022001 (2022), eprint 2107.08643.
  * Hogan (1986) C. J. Hogan, Mon. Not. Roy. Astron. Soc. 218, 629 (1986).
  * Maggiore (2000) M. Maggiore, Phys. Rept. 331, 283 (2000), eprint gr-qc/9909001.
  * Kawana (2022) K. Kawana (2022), eprint 2201.00560.
  * Callan and Coleman (1977) C. G. Callan, Jr. and S. R. Coleman, Phys. Rev. D 16, 1762 (1977).
  * Guth and Tye (1980) A. H. Guth and S. H. H. Tye, Phys. Rev. Lett. 44, 631 (1980), [Erratum: Phys.Rev.Lett. 44, 963 (1980)].
  * Guth and Weinberg (1981) A. H. Guth and E. J. Weinberg, Phys. Rev. D 23, 876 (1981).
  * Affleck (1981) I. Affleck, Phys. Rev. Lett. 46, 388 (1981).
  * Dine et al. (1992) M. Dine, R. G. Leigh, P. Y. Huet, A. D. Linde, and D. A. Linde, Phys. Rev. D 46, 550 (1992), eprint hep-ph/9203203.
  * Espinosa et al. (2010) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 06, 028 (2010), eprint 1004.4187.
  * Megevand and Ramirez (2017) A. Megevand and S. Ramirez, Nucl. Phys. B 919, 74 (2017), eprint 1611.05853.
  * Ellis et al. (2019) J. Ellis, M. Lewicki, and J. M. No, JCAP 04, 003 (2019), eprint 1809.08242.
  * Wang et al. (2020) X. Wang, F. P. Huang, and X. Zhang, JCAP 05, 045 (2020), eprint 2003.08892.
  * Bodeker and Moore (2009) D. Bodeker and G. D. Moore, JCAP 05, 009 (2009), eprint 0903.4099.
  * Bodeker and Moore (2017) D. Bodeker and G. D. Moore, JCAP 05, 025 (2017), eprint 1703.08215.
  * Höche et al. (2021) S. Höche, J. Kozaczuk, A. J. Long, J. Turner, and Y. Wang, JCAP 03, 009 (2021), eprint 2007.10343.
  * Gouttenoire et al. (2021) Y. Gouttenoire, R. Jinno, and F. Sala (2021), eprint 2112.07686.
  * Arakawa et al. (2021) J. Arakawa, A. Rajaraman, and T. M. P. Tait (2021), eprint 2109.13941.
  * Linde (1983) A. D. Linde, Nucl. Phys. B 216, 421 (1983), [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * Coleman (1977) S. R. Coleman, Phys. Rev. D 15, 2929 (1977), [Erratum: Phys.Rev.D 16, 1248 (1977)].
  * Rintoul and Torquato (1997a) M. D. Rintoul and S. Torquato, Journal of physics a: mathematical and general 30, L585 (1997a).
  * Holcomb and Rehr (1969) D. F. Holcomb and J. J. Rehr, Phys. Rev. 183, 773 (1969), URL https://link.aps.org/doi/10.1103/PhysRev.183.773.
  * Rintoul and Torquato (1997b) M. D. Rintoul and S. Torquato, Journal of Physics A: Mathematical and General 30, L585 (1997b), URL https://doi.org/10.1088/0305-4470/30/16/005.
  * Li et al. (2020) M. Li, H. Chen, and J. Lin, Computer Methods in Applied Mechanics and Engineering 361, 112815 (2020), ISSN 0045-7825, URL https://www.sciencedirect.com/science/article/pii/S0045782519307078.
  * Green and Kavanagh (2021) A. M. Green and B. J. Kavanagh, J. Phys. G 48, 043001 (2021), eprint 2007.10722.
  * Carr et al. (2021) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rept. Prog. Phys. 84, 116902 (2021), eprint 2002.12778.

Appendix - Angular Contribution

We present an alternative derivation of the angular factor in three dimensions
(Eq. (13)). Here we order the walls and set up the spherical coordinate system
in the same manner as in the main body of the text. Requiring
$-\hat{\textbf{r}}_{4}$ to be in the solid angle established by
$\hat{\textbf{r}}_{1}$,$\hat{\textbf{r}}_{2}$,and $\hat{\textbf{r}}_{3}$
results in the two conditions

$\pi\leq\phi_{4}\leq\pi+\phi_{3}~{},$ (27)

and

$\max(\theta_{2},\theta_{3},\theta_{4})\geq\pi-\min(\theta_{2},\theta_{3},\theta_{4})~{}.$
(28)

The $\phi$ condition can be understand as the requirement that walls 2, 3 and
4 form a triangle in the plane of wall 1, and so is analogous to the two-
dimensional case (Eq. (9)). The $\theta$ condition can be interpreted as the
need to form a closed tetrahedron in the $\hat{z}$ direction. To understand
this condition, first sort the angles
$\theta_{2},\theta_{3},\theta_{4}\xrightarrow{}\theta_{l}\leq\theta_{m}\leq\theta_{h}$.
Take the limiting case where wall $h$ is arranged opposite of wall $l$, so
that the corresponding azimuthal angles satisfy $\phi_{h}=\pi+\phi_{l}$. We
see that the closure requirement is satisfied if
$\pi-\theta_{l}\leq\theta_{h}\leq\pi$, where at the lower bound wall $h$ lies
directly opposite of wall $l$ and at the upper bound, wall $h$ is directly
opposite wall $1$.

The, the $\phi$ contribution to the three-dimension probability integral is

$\int_{0}^{2\pi}d\phi_{1}\int_{0}^{2\pi}d\phi_{2}\int_{0}^{\pi}2d\phi_{3}\int_{\pi}^{\pi+\phi_{3}}d\phi_{4}=4\pi^{4}~{},$
(29)

and the $\theta$ contribution is

$\displaystyle\begin{split}&\int_{0}^{\pi}d\theta_{1}\sin\theta_{1}~{}3!\int_{0}^{\pi}d\theta_{l}\sin\theta_{l}\int_{\max(\pi-\theta_{l},\theta_{l})}^{\pi}d\theta_{h}\sin\theta_{h}\\\
&\times\int_{\theta_{l}}^{\theta_{h}}d\theta_{m}\sin\theta_{m}=8~{}.\end{split}$
(30)

The combinatoric factor of $3!$ comes from the sorting of
$\theta_{2}$,$\theta_{3}$, and $\theta_{4}$ into the low, medium, and high
angles. As in Eq. (13), the combined angular contribution is $32\pi^{4}$.


# A Comparison of Machine Learning Methods for Data with High-Cardinality
Categorical Variables

Fabio Sigrist
Lucerne University of Applied Sciences and Arts Email<EMAIL_ADDRESS>Address: Lucerne University of Applied Sciences and Arts, Suurstoffi 1, 6343
Rotkreuz, Switzerland.

###### Abstract

High-cardinality categorical variables are variables for which the number of
different levels is large relative to the sample size of a data set, or in
other words, there are few data points per level. Machine learning methods can
have difficulties with high-cardinality variables. In this article, we
empirically compare several versions of two of the most successful machine
learning methods, tree-boosting and deep neural networks, and linear mixed
effects models using multiple tabular data sets with high-cardinality
categorical variables. We find that, first, machine learning models with
random effects have higher prediction accuracy than their classical
counterparts without random effects, and, second, tree-boosting with random
effects outperforms deep neural networks with random effects.

Keywords: Mixed effects machine learning, random effects, high-cardinality
categorical variables, tree-boosting, deep neural networks

##  1 Introduction

High-cardinality categorical variables are variables for which the number of
different levels is large relative to the sample size of a data set or,
equivalently, there is little data per level. High-cardinality categorical
variables can pose difficulties for machine learning methods such as deep
neural networks and tree-based models.

A simple strategy for dealing with categorical variables is to use one-hot
encoding or dummy variables. But this approach often does not work well for
high-cardinality categorical variables due to the reasons described below. For
neural networks, a frequently adopted solution is to use entity embeddings
(Guo and Berkhahn, 2016) that map every level of a categorical variable into a
low-dimensional Euclidean space. For tree-boosting, an alternative to one-hot
encoding is to assign a number to every level of a categorical variable, and
then consider this as a one-dimensional numeric variable. Another solution
implemented in the LightGBM boosting library (Ke et al., 2017) works by
partitioning all levels into two subsets using an approximate approach
(Fisher, 1958) when finding splits in the tree-building algorithm. Further,
the CatBoost boosting library (Prokhorenkova et al., 2018) implements an
approach based on ordered target statistics calculated using random partitions
of the training data for handling categorical predictor variables.

Random effects (Laird et al., 1982; Pinheiro and Bates, 2006) can also be used
as a tool for handling high-cardinality categorical variables. In a random
effects model, it is assumed that a (potentially transformed) parameter
$\mu\in\mathbb{R}^{n}$ of the response variable distribution equals the sum of
fixed $F(X)$ and random effects $Zb$:

$\mu=F(X)+Zb,~{}~{}~{}~{}b\sim\mathcal{N}(0,\Sigma),$

where $F(X)$ is the row-wise evaluation of a function
$F(\cdot):\mathbb{R}^{p}\rightarrow\mathbb{R}$,
$F(X)=(F(X_{1}),\dots,F(X_{n}))^{T}$,
$X_{i}=(X_{i1}\dots,X_{ip})^{T}\in\mathbb{R}^{p}$ is the $i$-th row of the
fixed effects variables matrix $X\in\mathbb{R}^{n\times p}$, $i=1,\dots,n$,
$b\in\mathbb{R}^{m}$, and $Z\in\mathbb{R}^{n\times m}$. These models are
called mixed effects models since they contain both fixed effects $F(X)$ and
random effects $Zb$. If the conditional response variable distribution is
Gaussian and there is a single high-cardinality categorical variable, such a
mixed effects model can also be written as

$y_{ij}=F(x_{ij})+b_{i}+\epsilon_{ij},~{}~{}b_{i}\overset{iid}{\sim}\mathcal{N}(0,\sigma^{2}_{1}),~{}~{}\epsilon_{ij}\overset{iid}{\sim}\mathcal{N}(0,\sigma^{2}),$
(1)

where $j=1,\dots,n_{i}$ is the sample index within level $i$ with $n_{i}$
being the number of samples for which the categorical variable attains level
$i$, $i=1,\dots,q$ is the level index with $q$ being the number of levels of
the categorical variable, and $x_{ij}$ are the fixed effects predictor
variables for observation $ij$. The total number of samples is
$n=\sum_{i=1}^{q}n_{i}$. Further, the random effects $b_{i}$ and
$\epsilon_{ij}$ are assumed to be independent. For this model, the matrix $Z$
is simply a binary incidence matrix that maps every random effect $b_{i}$ to
its corresponding observations and $\Sigma=\sigma^{2}I_{m}$.

In (generalized) linear mixed effects models it is assumed that $F(\cdot)$ is
a linear function: $F(X)=X\beta$. In the last years, linear mixed effects
models have been extended to non-linear ones using single trees (Hajjem et
al., 2011; Sela and Simonoff, 2012; Fu and Simonoff, 2015), random forest
(Hajjem et al., 2014), tree-boosting (Sigrist, 2022, 2023), and most recently
(in terms of first public preprint) deep neural networks (Simchoni and Rosset,
2021, 2023; Avanzi et al., 2023). In contrast to classical independent machine
learning models, the random effects introduce dependence among samples.

###  1.1 Why are random effects useful for high-cardinality categorical
variables?

For high-cardinality categorical variables, there is little data for every
level. Intuitively, if the response variable has a different (conditional)
mean for many levels, traditional machine learning models (with, e.g., one-hot
encoding, embeddings, or simply one-dimensional numeric variables) may have
problems with over- or underfitting for such data. From the point of view of a
classical bias-variance trade-off, independent machine learning models may
have difficulties balancing this trade-off and finding an appropriate amount
of regularization. For instance, overfitting may occur which means that a
model has a low bias but high variance.

Broadly speaking, random effects act as a prior, or regularizer, which models
the difficult part of a function, i.e., the part whose “dimension” is similar
to the total sample size, and, in doing so, provide an effective way for
finding a balance between over- and underfitting or bias and variance. For
instance, for a single categorical variable, random effects models will shrink
estimates of group intercept effects towards the global mean. This process is
sometimes also called “information pooling”. It represents a trade-off between
completely ignoring the categorical variable (= underfitting / high bias and
low variance) and giving every level in the categorical variable “complete
freedom” in estimation (= overfitting / low bias and high variance).
Importantly, the amount of regularization, which is determined by the variance
parameters of the model, is learned from the data. Specifically, in the above
single-level random effects model in (1), a (point) prediction $\hat{y}^{p}$
for the response variable for a sample with predictor variables $x^{p}$ and
categorical variable having level $i$ is given by

$\hat{y}^{p}=\hat{F}(x^{p})+\frac{\hat{\sigma}^{2}_{1}}{\hat{\sigma}^{2}/n_{i}+\hat{\sigma}^{2}_{1}}(\bar{y}_{i}-\bar{F}_{i}),$

where $\hat{F}(x^{p})$ is the trained function evaluated at $x^{p}$,
$\hat{\sigma}^{2}_{1}$ and $\hat{\sigma}^{2}$ are variance estimates, and
$\bar{y}_{i}$ and $\bar{F}_{i}$ are sample means of $y_{ij}$ and $F(x_{ij})$,
respectively, for level $i$. Ignoring the categorical variable would give the
prediction $\hat{y}^{p}=\hat{F}(x^{p})$, and a fully flexible model without
regularization gives $\hat{y}^{p}=\hat{F}(x^{p})+(\bar{y}_{i}-\bar{F}_{i})$.
I.e., the difference between these two extreme cases and the random effects
model is the shrinkage factor
$\frac{\hat{\sigma}^{2}_{1}}{\hat{\sigma}^{2}/n_{i}+\hat{\sigma}^{2}_{1}}$
(which goes to zero if the number of samples $n_{i}$ for level $i$ is large).
Related to this, random effects models allow for more efficient (i.e., lower
variance) estimation of the fixed effects function $F(\cdot)$ (Sigrist, 2022).
See also Sigrist (2023, Section 1.1) for a discussion on why random effects
are useful for modeling high-cardinality categorical variables.

In line with the above argumentation, Sigrist (2023, Section 4.1) find in
empirical experiments that tree-boosting combined with random effects
outperforms traditional independent tree-boosting the more, the lower the
number of samples per level of a categorical variable, i.e., the higher the
cardinality of a categorical variable.

##  2 Methods, data sets, and experimental settings

In the following, we compare several methods using multiple real-world data
sets with high-cardinality categorical variables. We use all the publicly
available tabular data sets from Simchoni and Rosset (2021, 2023) and also the
same experimental setting as in Simchoni and Rosset (2021, 2023). In addition,
we include the Wages data set analyzed in Sigrist (2022).

###  2.1 Methods

We consider the following methods:

  * •

‘Linear’: linear mixed effects models

  * •

‘NN_Embed’: deep neural networks with embeddings

  * •

‘LMMNN’: combining deep neural networks and random effects (Simchoni and
Rosset, 2021, 2023)

  * •

‘LGBM_Num’: tree-boosting by assigning a number to every level of categorical
variables and considering these as one-dimensional numeric variables

  * •

‘LGBM_Cat’: tree-boosting with the approach of LightGBM (Ke et al., 2017) for
categorical variables

  * •

‘CatBoost’: tree-boosting with the approach of CatBoost (Prokhorenkova et al.,
2018) for categorical variables

  * •

‘GPBoost’: combining tree-boosting and random effects (Sigrist, 2022, 2023)

The MERF algorithm of Hajjem et al. (2014) is another promising mixed effects
machine learning method, but its current implementation in the form of a
Python package111https://github.com/manifoldai/merf is prohibitively slow for
the sample sizes of the data sets considered here. Also note that, recently
(starting with version 1.6), the XGBoost library (Chen and Guestrin, 2016) has
also implemented the same approach as LightGBM for handling categorical
variables.222See
https://xgboost.readthedocs.io/en/stable/tutorials/categorical.html#optimal-
partitioning (retrieved on June 30, 2023) We do not consider this as a
separate approach here.

###  2.2 Data sets

Table 1 gives an overview of the data sets. For more details on the data, we
refer to Simchoni and Rosset (2021, 2023) and Sigrist (2022).

Dataset | $n$ | $p$ | $K$ | Cat. var. | Nb. levels | Response var.
---|---|---|---|---|---|---
Airbnb |  $50$K | $196$ | $1$ | host |  $39$K | price (log)
IMDb |  $88$K | $159$ | $2$ | director |  $38$K | avg. movie score
|  |  |  | movie type |  $1.7$K |
Spotify |  $28$K | $14$ | $4$ | artist |  $10$K | song danceability
|  |  |  | album |  $22$K |
|  |  |  | playlist |  $2.3$K |
|  |  |  | subgenre | $553$ |
News |  $81$K | $176$ | $2$ | source |  $5.4$K | nb. shares (log)
|  |  |  | title |  $72$K |
InstEval |  $73$K | $3$ | $3$ | student |  $2.9$K | teacher rating
|  |  |  | teacher |  $1.1$K |
|  |  |  | department | $14$ |
Rossmann |  $33$K | $23$ | $1$ | store |  $1.1$K | total sales (in 100K)
AUimport |  $125$K | $8$ | $1$ | commodity |  $5$K | total import (log)
Wages |  $28$K | $52$ | $1$ | person |  $4.7$K | wage (log)
Table 1: Summary of data sets. $n$ is the number of samples, $p$ is the number
of predictor variables (excl. high-cardinality categorical variables), $K$ is
the number of high-cardinality categorical variables, and ’Cat. var.’
describes the categorical variable(s).

For all methods with random effects, we include random effects for every
categorical variable mentioned in Table 1 with no prior correlation among
random effects. The Rossmann, AUImport, and Wages data sets are longitudinal
data sets. This means that the samples for every level or the categorical
variables are repeated measurements over time $t$. As in Simchoni and Rosset
(2023), we additionally include random coefficients (= random slopes), for the
time variables $t$ and $t^{2}$ besides intercept random effects with no prior
correlation among the random effects. I.e., these random effects models are
given by

$y_{ij}=F(x_{ij})+b_{0,i}+b_{1,i}\cdot t_{ij}+b_{2,i}\cdot
t_{ij}^{2}+\epsilon_{ij},~{}~{}b_{k,i}\overset{iid}{\sim}\mathcal{N}(0,\sigma^{2}_{k}),k\in\\{1,2,3\\},~{}~{}\epsilon_{ij}\overset{iid}{\sim}\mathcal{N}(0,\sigma^{2}),$

where $b_{1,i}$ and $b_{2,i}$ are the random slopes.

###  2.3 Experimental setting

We use the same experimental setting as in Simchoni and Rosset (2023) to
compare the different methods. This means that we perform 5-fold cross-
validation with the test mean squared error (MSE) to measure prediction
accuracy. For the longitudinal data sets, we use the “random mode” 5-fold CV
setting of Simchoni and Rosset (2023) where in every fold, 80% of the data is
used for training and 20% of the data for validation irrespective of the time
variable $t$. The results for neural networks with embeddings and also neural
networks with random effects (LMMNN) are taken from Simchoni and Rosset (2021,
2023); see Simchoni and Rosset (2021, 2023) for more details on the
specifications of these two modeling approaches. For linear mixed effects
models, LGBM_Num, LGBM_Cat, GPBoost, and GPBoost_I, we use the GPBoost library
version 1.2.1 which is built upon the LightGBM library. Further, we use
version 1.1.1 of the CatBoost library.

For all tree-boosting methods, we choose tuning parameters on every of the
five training sets in the $5$-fold CV by randomly splitting the training data
into inner training data containing 80% of the outer training data and
validation data consisting of the remaining 20%. We use a deterministic grid
search with the mean squared error as a selection criterion and consider the
following parameter combinations: number of boosting iterations
$M\in\\{1,\dots,1000\\}$, learning rate $\in\\{1,0.1,0.01\\}$, maximal tree-
depth $\in\\{1,2,3,5,10\\}$, minimal number of samples per leaf
$\in\\{10,100,1000\\}$, and L2 penalty on leaf values $\in\\{0,1,10\\}$. For
the machine learning models with random effects, one can either exclude or
include the high-cardinality categorical variables in the fixed effects
function. When additionally including them, one allows for potential
interaction between the categorical variables and other predictor variables.
For GPBoost, we consider this as a tuning parameter option. However, the
difference between these two modeling options is minor which is an indication
that there is no interaction present; see Table 4 in the appendix, where we
report the results when separately either excluding or including the high-
cardinality categorical variables in the fixed effects tree ensemble function.

Code for pre-processing the data with instructions on how to download the data
and code for running the experiments can be found at
https://github.com/fabsig/Compare_ML_HighCardinality_Categorical_Variables.
Pre-processed data for modeling can also be found on the above webpage for
data sets for which the license of the original source permits it.

##  3 Results

The results are reported in Figure 1 and Table 2. Table 2 reports test mean
squared errors (MSE) and corresponding standard errors. To summarize the
prediction accuracy over the different data sets, we report average relative
differences to the best result and average ranks. The former is obtained by
calculating the relative difference of a test MSE of a method compared to the
lowest MSE for every data set and then taking the average over all data sets.
In Figure 1, we also report the average relative difference to the best
result. In addition, Table 3 in the appendix reports average wall clock times.

Figure 1: Average relative difference (in %) to the lowest test MSE. The Wages data set is not included for calculating this since not all methods were run on it. Dataset | Linear | LGBM_Cat | LGBM_Num | CatBoost | GPBoost | NN_Embed | LMMNN
---|---|---|---|---|---|---|---
Airbnb | 0.268 | 0.131 | 0.132 | 0.131 | 0.125 | 0.158 | 0.142
| (0.0748) | (0.00283) | (0.00285) | (0.00267) | (0.00264) | (0.00158) | (0.00167)
IMDb | 1.02 | 0.954 | 1.05 | 0.922 | 0.850 | 1.26 | 0.974
| (0.00730) | (0.00721) | (0.00830) | (0.0117) | (0.00825) | (0.124) | (0.00933)
Spotify | 0.0111 | 0.00858 | 0.00880 | 0.00851 | 0.00910 | 0.0164 | 0.00926
| (9.20e-05) | (6.23e-05) | (5.79e-05) | (4.69e-05) | (0.000164) | (0.000540) | (8.34e-05)
News | 1.88 | 1.89 | 2.37 | 1.85 | 1.72 | 1.89 | 1.81
| (0.0146) | (0.0147) | (0.0157) | (0.0160) | (0.0131) | (0.0215) | (0.0190)
InstEval | 1.44 | 1.44 | 1.53 | 1.46 | 1.47 | 1.50 | 1.45
| (0.00672) | (0.00746) | (0.00549) | (0.00711) | (0.0230) | (0.00669) | (0.00492)
Rossmann | 0.0153 | 0.00627 | 0.0124 | 0.01000 | 0.00877 | 0.0516 | 0.0105
| (0.000552) | (0.000188) | (0.000848) | (0.000267) | (0.000284) | (0.00600) | (0.000162)
AUimport | 0.743 | 1.26 | 4.56 | 2.15 | 0.651 | 3.35 | 0.713
| (0.0138) | (0.0269) | (0.0177) | (0.0477) | (0.00396) | (0.455) | (0.00700)
Wages | 0.321 | 0.102 | 0.0992 | 0.0920 | 0.0830 |  |
| (0.0222) | (0.00259) | (0.00234) | (0.00251) | (0.00200) |  |
Avg. rank | 4.71 | 2.71 | 5.57 | 2.86 | 2.14 | 6.43 | 3.57
Avg. rel. dif. | 47.6 | 17.4 | 111. | 44.5 | 7.02 | 189. | 17.3
Table 2: Average test mean squared error (MSE) and standard errors in
parentheses. Best results are bold. ’Avg. rank’ denotes the average rank of a
method. ’Avg. rel. dif.’ denotes the average relative difference in % of a
method compared to the best method. The Wages data set is not included for
calculating ranks and relative differences since not all methods were run on
it.

The results in Table 2 show that combined tree-boosting and random effects
(GPBoost) has the highest prediction accuracy. GPBoost has an average relative
difference to the best result of $7.02\%$ and an average rank of $2.14$.
Combined neural networks and random effects (LMMNN) have an average relative
difference to the best result of $17.3\%$ and an average rank of $3.57$.
LightGBM (LGBM_Cat) has a similar average relative difference to the best
result of $17.4\%$ and an average rank of $2.71$. Next, CatBoost has an
average relative difference to the best method of $44.5\%$ and an average rank
of $2.86$. Linear mixed effects models perform worse having an average
relative difference of approximately $47.6\%$ and an average rank of $4.71$.
Overall worst perform neural networks with embeddings having an average
relative difference to the best result of $189\%$. Tree-boosting with the
categorical variables transformed to one-dimensional numeric variables
(LGBM_Num) performs slightly better with an average relative difference to the
best result of $111\%$.333In their online documentation, LightGBM recommends
“For a categorical feature with high cardinality, it often works best to treat
the feature as numeric …”; see
https://lightgbm.readthedocs.io/en/latest/Advanced-Topics.html#categorical-
feature-support (retrieved on June 30, 2023). We clearly come to a different
conclusion.

##  4 Conclusion

We empirically compare various versions of tree-boosting and deep neural
networks as well as linear mixed effects models on multiple tabular data sets
with high-cardinality variables. We find that, first, machine learning models
with random effects have higher prediction accuracy than their classical
independent counterparts without random effects and, second, tree-boosting
with random effects performs better than deep neural networks with random
effects. While there may be several reasons for the latter finding, this is in
line with the recent work of Grinsztajn et al. (2022) who find that tree-
boosting outperforms deep neural networks (and also random forest) on tabular
data without high-cardinality categorical variables. Similarly, Shwartz-Ziv
and Armon (2022) conclude that tree-boosting “outperforms deep models on
tabular data.”

##  Appendix A Additional results

Dataset | Linear | LGBM_Cat | LGBM_Num | CatBoost | GPBoost | NN_Embed | LMMNN
---|---|---|---|---|---|---|---
Airbnb | 1.93 | 6.41 | 0.952 | 16.8 | 14.4 | 825. | 95.4
IMDb | 183. | 1.82 | 1.63 | 19.8 | 236. | 55.8 | 144.
Spotify | 99.1 | 5.06 | 1.10 | 7.07 | 90.6 | 12.7 | 57.1
News | 138. | 1.02 | 1.54 | 93.1 | 63.9 | 31.1 | 83.9
InstEval | 95.3 | 4.64 | 0.773 | 25.5 | 682. | 51.1 | 92.5
Rossmann | 8.22 | 3.61 | 1.32 | 54.6 | 1.45 | 53.6 | 76.3
AUimport | 37.7 | 8.87 | 0.604 | 17.7 | 45.9 | 107. | 175.
Wages | 19.6 | 3.61 | 2.41 | 43.6 | 34.2 |  |
Table 3: Average wall-clock time in seconds. Note that except for the linear model, the wall-clock time depends on the chosen tuning parameters. The experiments for ’NN_Embed’ and ’LMMNN’ were run on Google Colab with NVIDIA Tesla V100 GPU machines (see Simchoni and Rosset, 2023), and all for all other methods, a laptop with an Intel i7-12800H processor was used. Dataset | Linear | LGBM_Cat | LGBM_Num | CatBoost | GPBoost_E | GPBoost_I | NN_Embed | LMMNN
---|---|---|---|---|---|---|---|---
Airbnb | 0.268 | 0.131 | 0.132 | 0.131 | 0.124 | 0.125 | 0.158 | 0.142
| (0.0748) | (0.00283) | (0.00285) | (0.00267) | (0.00272) | (0.00265) | (0.00158) | (0.00167)
IMDb | 1.02 | 0.954 | 1.05 | 0.922 | 0.883 | 0.850 | 1.26 | 0.974
| (0.00730) | (0.00721) | (0.00830) | (0.0117) | (0.00837) | (0.00825) | (0.124) | (0.00933)
Spotify | 0.0111 | 0.00858 | 0.00880 | 0.00851 | 0.00894 | 0.00986 | 0.0164 | 0.00926
| (9.20e-05) | (6.23e-05) | (5.79e-05) | (4.69e-05) | (6.99e-05) | (0.000227) | (0.000540) | (8.34e-05)
News | 1.88 | 1.89 | 2.37 | 1.85 | 1.78 | 1.72 | 1.89 | 1.81
| (0.0146) | (0.0147) | (0.0157) | (0.0160) | (0.0143) | (0.0131) | (0.0215) | (0.0190)
InstEval | 1.44 | 1.44 | 1.53 | 1.46 | 1.44 | 1.47 | 1.50 | 1.45
| (0.00672) | (0.00746) | (0.00549) | (0.00711) | (0.00685) | (0.0230) | (0.00669) | (0.00492)
Rossmann | 0.0153 | 0.00627 | 0.0124 | 0.01000 | 0.00910 | 0.00877 | 0.0516 | 0.0105
| (0.000552) | (0.000188) | (0.000848) | (0.000267) | (0.000229) | (0.000284) | (0.00600) | (0.000162)
AUimport | 0.743 | 1.26 | 4.56 | 2.15 | 0.714 | 0.651 | 3.35 | 0.713
| (0.0138) | (0.0269) | (0.0177) | (0.0477) | (0.00550) | (0.00396) | (0.455) | (0.00700)
Wages | 0.321 | 0.102 | 0.0992 | 0.0920 | 0.0829 | 0.0830 |  |
| (0.0222) | (0.00259) | (0.00234) | (0.00251) | (0.00191) | (0.00205) |  |
Avg. rank | 5.71 | 3.29 | 6.43 | 3.71 | 2.43 | 2.71 | 7.43 | 4.29
Avg. rel. dif. | 47.7 | 17.4 | 111. | 44.5 | 9.63 | 8.36 | 189. | 17.3
Table 4: Results when either excluding (’GPBoost_E’) or including
(’GPBoost_I’) the categorical variables in the fixed effects tree ensemble
function for GPBoost.

## References

  * Avanzi et al. [2023] B. Avanzi, G. Taylor, M. Wang, and B. Wong.  Machine Learning with High-Cardinality Categorical Features in Actuarial Applications.  _arXiv preprint arXiv:2301.12710_ , 2023.
  * Chen and Guestrin [2016] T. Chen and C. Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining_ , pages 785–794. ACM, 2016.
  * Fisher [1958] W. D. Fisher.  On grouping for maximum homogeneity.  _Journal of the American statistical Association_ , 53(284):789–798, 1958.
  * Fu and Simonoff [2015] W. Fu and J. S. Simonoff.  Unbiased regression trees for longitudinal and clustered data.  _Computational Statistics & Data Analysis_, 88:53–74, 2015.
  * Grinsztajn et al. [2022] L. Grinsztajn, E. Oyallon, and G. Varoquaux.  Why do tree-based models still outperform deep learning on typical tabular data?  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 507–520. Curran Associates, Inc., 2022.
  * Guo and Berkhahn [2016] C. Guo and F. Berkhahn.  Entity embeddings of categorical variables.  _arXiv preprint arXiv:1604.06737_ , 2016.
  * Hajjem et al. [2011] A. Hajjem, F. Bellavance, and D. Larocque.  Mixed effects regression trees for clustered data.  _Statistics & probability letters_, 81(4):451–459, 2011.
  * Hajjem et al. [2014] A. Hajjem, F. Bellavance, and D. Larocque.  Mixed-effects random forest for clustered data.  _Journal of Statistical Computation and Simulation_ , 84(6):1313–1328, 2014.
  * Ke et al. [2017] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu.  LightGBM: A highly efficient gradient boosting decision tree.  In _Advances in Neural Information Processing Systems_ , pages 3149–3157, 2017.
  * Laird et al. [1982] N. M. Laird, J. H. Ware, et al.  Random-effects models for longitudinal data.  _Biometrics_ , 38(4):963–974, 1982.
  * Pinheiro and Bates [2006] J. Pinheiro and D. Bates.  _Mixed-effects models in S and S-PLUS_.  Springer Science & Business Media, 2006.
  * Prokhorenkova et al. [2018] L. Prokhorenkova, G. Gusev, A. Vorobev, A. V. Dorogush, and A. Gulin.  CatBoost: unbiased boosting with categorical features.  In _Advances in Neural Information Processing Systems_ , pages 6638–6648, 2018.
  * Sela and Simonoff [2012] R. J. Sela and J. S. Simonoff.  RE-EM trees: a data mining approach for longitudinal and clustered data.  _Machine learning_ , 86(2):169–207, 2012.
  * Shwartz-Ziv and Armon [2022] R. Shwartz-Ziv and A. Armon.  Tabular data: Deep learning is not all you need.  _Information Fusion_ , 81:84–90, 2022.
  * Sigrist [2022] F. Sigrist.  Gaussian Process Boosting.  _The Journal of Machine Learning Research_ , 23(1):10565–10610, 2022.
  * Sigrist [2023] F. Sigrist.  Latent Gaussian Model Boosting.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 45(2):1894–1905, 2023.
  * Simchoni and Rosset [2021] G. Simchoni and S. Rosset.  Using random effects to account for high-cardinality categorical features and repeated measures in deep neural networks.  _Advances in Neural Information Processing Systems_ , 34:25111–25122, 2021.
  * Simchoni and Rosset [2023] G. Simchoni and S. Rosset.  Integrating Random Effects in Deep Neural Networks.  _Journal of Machine Learning Research_ , 24(156):1–57, 2023.