Electroweak phase transition in the presence of hypermagnetic field and the
generation of gravitational waves
H. Abedia<EMAIL_ADDRESS>M. Ahmadvandb
<EMAIL_ADDRESS>and S. S. Gousheha<EMAIL_ADDRESS>
aDepartment of Physics, Shahid Beheshti University, G. C., Evin, Tehran 19839,
Iran
bSchool of Physics, Institute for Research in Fundamental Sciences (IPM), P.
O. Box 19395-5531, Tehran, Iran

We investigate the effect of a large-scale background hypermagnetic field on
the electroweak phase transition. We propose an effective weak angle which is
varying during the electroweak phase transition and upon its use, show that
for a strong enough hypermagnetic field the phase transition occurs in two
steps and becomes first-order. We obtain all of the important quantities
characterizing the details of the phase transition, including the phase
transition latent heat, temperature and duration. We then explore one of the
consequences of this model which is the generation of gravitational waves. We
calculate the gravitational wave spectrum generated during the first-order
electroweak phase transition and find that, for strong enough background
hypermagnetic fields, these signals can be detected by the Ultimate-DECIGO and
BBO correlated interferometers.

##  1 Introduction

Magnetic fields are ubiquitous in the Universe. The presence of magnetic
fields in stars, galaxies and intergalactic spaces raises the question about
the origin of these fields [1]. One possibility is that these fields have been
produced by some astrophysical and dynamo processes [2]. However, using these
mechanisms it is hard to explain large-scale magnetic fields observed in the
intergalactic spaces. The other scenario is that these long-range magnetic
fields can be relics of the early Universe. In this case, several mechanisms
for magnetogenesis have been proposed to be taking place during the early
times such as inflation and cosmological phase transitions [3].
Considering such primordial magnetic fields, one can study their effects on
significant events taking place after the big bang, e.g. the Electroweak Phase
Transition (EWPT). Assuming these magnetic fields are produced before the
EWPT, long-range component of these fields, i.e. the hypermagnetic fields, can
survive in the plasma with high conductivity. Due to the chiral Abelian
anomaly, these hypermagnetic fields can play an important role in the EW
baryogenesis scenarios [4, 5]. Furthermore, a background hypermagnetic field,
$B^{\mathrm{bg}}_{Y}$, can influence the nature of EWPT. Without such a field,
the Standard Model (SM) predicts that the EWPT is a crossover. However, the
presence of sufficiently strong hypermagnetic fields could change the
situation and make it first-order [6, 7]. On the other hand, due to the
coupling between these fields and magnetic dipole moment of sphalerons, the
energy barrier of sphalerons decreases so that sphaleron processes can
threaten the EW baryogenesis scenarios in this context [8]. In our previous
paper [9], based on gravitational anomaly and chiral gravitational waves (GWs)
sourced by helical magnetic fields, we found the possibility to violate $B-L$
symmetry, where $B$ and $L$ are baryon and lepton numbers, respectively. Then,
relying on sphaleron processes, we presented a possible mechanism for
generation of the matter-antimatter asymmetry of the Universe.
The influence of hypermagnetic fields on the EWPT has initially been addressed
in [6] and [7]. Considering only the direct effect of hypermagnetic fields on
the Gibbs free energy, they showed analytically and also numerically by
lattice simulation [10] that the EWPT becomes first-order for strong enough
hypermagnetic fields. In this work, starting from an appropriate Lagrangian
containing the effect of a constant $B^{\mathrm{bg}}_{Y}$ and calculating the
effective potential of symmetric and broken phases, we show that by
introducing an effective weak angle which is varying during the PT, one can
describe the PT dynamically. By analyzing this effective potential we find
that as the temperature falls below a temperature $T_{0}$, the vacuum of the
symmetric phase becomes unstable and the usual EW crossover transition to a
stable vacuum occurs. However, as the temperature decreases further, the
presence of a large $B^{\mathrm{bg}}_{Y}$ produces a second vacuum whose
stability eventually surpasses the first. Consequently, a first-order phase
transition occurs between these two well seperated vacuua. Then, by
calculating the bounce solution from the bounce action we find the transition
temperature at which the bubbles nucleate. At this temperature, the latent
heat of the PT and its relation to the strength of the hypermagnetic fields
are obtained. We also calculate the duration of the PT.
A first-order EWPT can have many profound consequenses, in particular for
matter-antimatter asymmetry generation as well as for the generation of
gravitational waves (GWs). Here we concentrate on the latter. The GWs are
useful probes, providing valuable information about the early Universe, partly
because they have the least attenuation during the propagation. In general,
any first-order PT in the early Universe is regarded as a source for the
generation of GWs111The relation of hypermagnetic fields to the GW radiation
at the EWPT is initiated in [11] [12, 13]. In this case during the evolution
of bubbles three sources can contribute to the production of the GW spectrum:
bubble collisions [14], Magnetohydrodynamic (MHD) turbulence [15], and sound
waves [16]. In this paper, by calculating the required quantities associated
with the PTs, we find the GW energy spectrum. We show that bubbles can run
away and the generated GWs can be in the sensitivity range of some future
space-based GW experiments, i.e. The Deci-Hertz Interferometer Gravitational
Wave Observatory (DECIGO) and the Big Bang Observer (BBO), whose the primary
objective is to track primordial GWs [17, 18].
In Section 2, we introduce the one-loop effective action in the presence of an
$B^{\mathrm{bg}}_{Y}$ and obtain the details of the PT including the bounce
solution, transition temperature, duration of the PT and latent heat. In
Section 3, we calculate the GWs spectra using the parameters obtained in the
previous section. We conclude in Section 4.

##  2 The electroweak phase transition

###  2.1 Model

The presence of a large-scale background hypermagnetic field,
$B^{\mathrm{bg}}_{Y}$, at the EWPT, regardless of its origin, can have
important consequences on the dynamics on the PT. To study these effects, one
can decompose the total hypercharge field strength as
$f^{t}_{\mu\nu}=f_{\mu\nu}+f^{bg}_{\mu\nu}$ and consequently the Lagrangian
can be written as

$\mathcal{L}=\mathcal{L}_{SM}-\frac{1}{2}(f^{bg})_{\mu\nu}f^{\mu\nu}-\frac{1}{4}(f^{bg})_{\mu\nu}(f^{bg})^{\mu\nu},$
(1)

where $\mathcal{L}_{SM}$ is the Lagrangian of the standard model. As shown in
[7], one can solve the equations of motion for both symmetric and broken
phases and then obtain the corresponding Gibbs free energies in terms of
$B^{\mathrm{bg}}_{Y}$. In the high conductive plasma of the symmetric phase,
only the large-scale hypermagnetic field can survive while other fields are
screened. Moreover, in the case of a first-order PT and existing of bubbles,
in order for the Z-component not to penetrate the bubbles, the wall width
should be larger than the correlation length of the field, i.e.
$l_{wall}>l_{Z}$. Only the Maxwellian component of $B^{\mathrm{bg}}_{Y}$
remains unscreened. Also, there is a criterion to avoid the W-condensation,
i.e. $eB^{\mathrm{bg}}_{Y}<(m_{W}\phi/v)^{2}$ [7]. We will investigate these
constraints in the next section.
The Gibbs free energies corresponding to the effective potentials in the
symmetric and broken phases are given by $V(0,T)$ and
$V(\phi,T)+1/2(B^{\mathrm{bg}}_{Y})^{2}\sin^{2}\theta_{w}$, respectively. Here
$\theta_{w}$ is the Weinberg or weak mixing angle. These free energies
correspond to those of [7] which have been shifted by
$1/2(B^{\mathrm{bg}}_{Y})^{2}$ due to the last term of Eq. (1), although the
physics of the system will be unchanged. Similar to the Meissner effect, one
can study effects of the extra energy difference between two phases,
$1/2(B^{\mathrm{bg}}_{Y})^{2}\sin^{2}\theta_{w}$, by a single form effective
potential. To obtain this potential, one can consider an effective weak angle
which varies during the EWPT from zero in the symmetric phase to its final
nonzero value, i.e., $\theta_{w}$, in the broken phase. An analogous scheme
for the weak angle varying with time has been studied in [19]. Since the Higgs
VEV varies with temperature, we propose a phenomenological Ansatz for the weak
effective angle as below,

$\Theta_{w}=\theta_{w}\Big{[}\frac{1}{2}+\frac{1}{2}\tanh\Big{(}\frac{\phi-m}{s}\Big{)}\Big{]},$
(2)

where $m$ and $s$ are two constants which determine the mid-value of $\phi$
and the smoothness of $\Theta_{w}$ profile, respectively. Consequently, we use
the following potential as the finite temperature effective potential which
describes the system during the phase transition and is valid in the symmetric
and broken phases

$V_{eff}(\phi,T,B^{\mathrm{bg}}_{Y})=V(\phi,T)+\frac{(B^{\mathrm{bg}}_{Y})^{2}}{2}\sin^{2}\Theta_{w}.$
(3)

The phenomenological Ansatz we have chosen for $\Theta_{w}(\phi)$ is a natural
smooth choice consistent with mentioned requirements. As a result, our
proposed effective potential coincides with the corresponding Gibbs free
energies in the symmetric and broken phases, and can describe the dynamics of
the PT continuously.
In Fig. (1), $\Theta_{w}$ is shown for $m$ = 95 and $s$ = 10. These parameters
are chosen in such a way as to not only avoid the Z-penetration of the bubbles
and W-condensation, but also make baryogenesis possible222In order for the
sphalerons to be active outside the bubbles before the transition temperature,
$v(T)/T$ must be sufficiently less than one. This condition can be satisfied
for $m<100$. On the other hand, as for avoiding W-condensation, $m>90$ is
acceptable. Hence, considering these two constraints, we can take $90<m<100$..

Figure 1:  Our proposed model for the weak angle as a function of $\phi$
during the EWPT for $m$ = 95 and $s$ = 10.

Conductivity of primordial plasma is proportional to its temperature [10], and
hence is very high. In such a plasma, the only long-range or background field
that can survive is $B^{\mathrm{bg}}_{Y}$. We assume such long-range
$B^{\mathrm{bg}}_{Y}$ exist and that their scale is much larger than the
typical size of bubbles during the PT. Therefore we take this field to be
constant on the scale of the bubbles and further assume
$B^{\mathrm{bg}}_{Y}=bT^{2}$ where $b$ is approximately constant during the
EWPT.
In Eq. (3), $V(\phi,T)$ consists of the following terms

$V(\phi,T)=V_{0}+V_{1}+V_{th},$ (4)

where the tree level potential of the Higgs field, $V_{0}$, is given by

$V_{0}(\phi)=-\frac{1}{2}\mu^{2}\phi^{2}+\frac{\lambda}{4}\phi^{4},$ (5)

where $\lambda=\mu^{2}/\nu^{2}$ is fixed by the Higgs mass,
$2\lambda\nu^{2}=125\,\mathrm{GeV}$, and Higgs VEV at zero temperature,
$\nu=246\,\mathrm{GeV}$. One-loop quantum correction can be written as [20]

$\displaystyle V_{1}(\phi)$ $\displaystyle=$
$\displaystyle\frac{6}{64\pi^{2}}\Big{[}m_{W}^{4}(\phi)\Big{(}\log\frac{m_{W}^{2}(\phi)}{m_{W}^{2}(\nu)}-\frac{5}{6}\Big{)}+2m_{W}^{2}(\phi)m_{W}^{2}(\nu)\Big{]}$
(6) $\displaystyle+$
$\displaystyle\frac{3}{64\pi^{2}}\Big{[}m_{Z}^{4}(\phi)\Big{(}\log\frac{m_{Z}^{2}(\phi)}{m_{Z}^{2}(\nu)}-\frac{5}{6}\Big{)}+2m_{Z}^{2}(\phi)m_{Z}^{2}(\nu)\Big{]}$
$\displaystyle-$
$\displaystyle\frac{12}{64\pi^{2}}\Big{[}m_{t}^{4}(\phi)\Big{(}\log\frac{m_{t}^{2}(\phi)}{m_{t}^{2}(\nu)}-\frac{3}{2}\Big{)}+2m_{t}^{2}(\phi)m_{t}^{2}(\nu)\Big{]},$

where $m_{W}(\phi)=g_{2}\phi/2$, $m_{Z}(\phi)=\sqrt{(g_{2}^{2}+g^{\prime
2})}\phi/2$ and $m_{t}(\phi)=y_{t}\phi/\sqrt{2}$ are masses of the massive
gauge fields and the top quark, respectively. The thermal correction term is
as follows [20]

$V_{th}(\phi,T)\equiv\sum_{i=W,Z,t}\pm\frac{n_{i}T^{4}}{2\pi^{2}}J_{B,F}\Big{(}\frac{m^{2}_{i}(\phi)}{T^{2}}\Big{)},$
(7)

where $n_{i}$ is the number of degrees of freedom of the mentioned particles
and

$J_{B,F}(x)=\int_{0}^{\infty}dy~{}y^{2}\log\Big{[}1\mp\exp\Big{(}-\sqrt{y^{2}+x}\Big{)}\Big{]}.$
(8)

High temperature expansion of the thermal corrections is as follows

$J_{B}(x)=\frac{\pi^{2}}{12}x-\frac{\pi}{6}x^{\frac{3}{2}}-\frac{x^{2}}{32}\log\frac{x}{a_{b}}+\mathcal{O}(x^{3}),$
(9)
$J_{F}(x)=-\frac{\pi^{2}}{24}x-\frac{x^{2}}{32}\log\frac{x}{a_{f}}+\mathcal{O}(x^{3}),$
(10)

where $\log a_{b}=5.4076$ and $\log a_{f}=2.6351$. In the next subsection,
within the above model, we calculate all of the major physical quantities
characterizing the phase transion.

###  2.2 Dynamics of the phase transition

In a typical PT, the existence of a barrier between two minima of the
effective potential produces a first-order PT. In [6, 7], using only the
difference of the Gibbs free energies of two phases and without considering
the PT dynamics, it is argued that the presence of a strong enough
$B^{\mathrm{bg}}_{Y}$ in the effective potential increases the area under the
barrier and delays the PT222It should be noted that in [6, 7] the value of
$B^{\mathrm{bg}}_{Y}$ required for the strongly first-order EWPT is calculated
before the observation of the Higgs particle.. In this paper, we use $V_{eff}$
as given by Eqs. (3, 2) to calculate the dynamics of the EWPT. In Fig. (2) we
plot $V_{eff}$ as a function of $\phi$ for a few important temperatures, all
at a fixed value of $B^{\mathrm{bg}}_{Y}$. In the symmetric phase and well
above the EWPT the system has a stable vacuum at $\phi=0$. As the temperature
decreases to $T_{0}\approx 163\,\mathrm{GeV}$, the second derivative at
$\phi=0$ becomes zero and the vacuum becomes metastable. As the temperature
decreases further, a stable minimum starts to form and the VEV continuously
increases from zero through the usual EW crossover. The new feature in the
presence of a strong $B^{\mathrm{bg}}_{Y}$ is the formation of a second
minimum, when the temperature decreases further. At the critical temperature,
$T_{c}$, the two vacua become degenerate. Shortly after $T_{c}$, at nucleation
or transition temperature, $T_{*}$, the newly formed vacuum becomes the
absolute minimum and tunneling to the true vacuum can be fulfilled. Finally,
at temperature $T_{f}\approx 118\,\mathrm{GeV}$ the barrier between the minima
disappears, i.e. only one vacuum remains. This temperature can be considered
as the end of the EWPT.

Figure 2:  The effective potential at different temperatures for
$B^{\mathrm{bg}}_{Y}=0.5\,T^{2}_{*}$ is displayed. The green dotted curve is
at $T=200\,\mathrm{GeV}$ far above the PT with one minimum at $\phi=0$, the
purple dashed curve at $T_{0}=163\,\mathrm{GeV}$ shows the beginning of
crossover, the blue dotdashed curve $T=145\,\mathrm{GeV}$ depicts a typical
case where the second local minimum is forming, and the red solid curve
corresponds to the nucleation temperature $T_{*}=135.7\,\mathrm{GeV}$. Since
$T_{c}=137.2\,\mathrm{GeV}$ is very close to $T_{*}$, we do not show its
curve. Also not shown is the curve for $T_{f}\approx 118\,\mathrm{GeV}$, where
only one minimum remains.

The probability of true vacuum bubble nucleation per unit Hubble spacetime
volume at finite temperature is given by [21]

$P\simeq\frac{M_{\mathrm{pl}}^{4}}{T^{4}}\exp\Big{(}-\frac{S_{3}(T)}{T}\Big{)},$
(11)

where $M_{\mathrm{pl}}$ is the Planck mass and $S_{3}(T)$ is the three-
dimensional Euclidean bounce action333It should also be noted that the
hypermagnetic field may deform the bubbles and degrade the spherical symmetry
to an axial symmetry. Therefore the calculations presented are an
approximation. An attempt to improve the calculations should include the back-
reaction of the plasma which might produce a restoring force.

$S_{3}(T)=\int_{0}^{\infty}4\pi
r^{2}~{}dr\Big{[}\frac{1}{2}\Big{(}\frac{d\phi}{dr}\Big{)}^{2}+V_{eff}(\phi,T,B^{\mathrm{bg}}_{Y})\Big{]}.$
(12)

As shown below, the presence of a strong hypermagnetic field can significantly
affect the bounce action, and in particular the EWPT is delayed. In fact, the
transition and bubble nucleation occurs when the bubble formation probability
is of the order of one, $P\sim 1$. From this condition, we can find the
transition temperature, $T_{*}$, [21]

$\frac{S_{3}(T_{*})}{T_{*}}=4\ln\Big{(}\frac{T_{*}}{H_{*}}\Big{)},$ (13)

where $H\simeq T^{2}/M_{\mathrm{pl}}$ is the Hubble expansion parameter. Thus,
we should first obtain the function $S_{3}(T)/T$. Extremizing the bounce
action leads to the bounce equation which is shown below along with the
appropriate boundary conditions,

$\frac{d^{2}\phi}{dr^{2}}+\frac{2}{r}\frac{d\phi}{dr}=\frac{\partial
V_{eff}}{\partial\phi},~{}~{}~{}~{}~{}~{}~{}\frac{d\phi}{dr}(r=0)=0,~{}~{}~{}~{}\phi(r=\infty)=\phi_{\mathrm{false}}.$
(14)

From these equations, we can find the radial profile of the Higgs, i.e. the
bounce solution, which, as we shall see, connects the two phases through a
step-like but smooth function. To obtain the solution, we use the “any bubble”
code [22] which takes advantage of a multiple-shooting method. The result is
shown in Fig. (3). The expectation value of the Higgs field inside the bubble
is completely consistent with the absolute minimum of the red solid curve in
Fig. (2). Moreover, as shown in Fig. (3), the wall width of the true vacuum
bubbles, $l_{w}\sim(0.04-0.06)\,\mathrm{GeV}^{-1}$, is larger than the
correlation length of Z-field, $l_{Z}\sim 1/m_{Z}=2/(\sqrt{g_{2}^{2}+g^{\prime
2}}v(T_{*}))\sim 0.01-0.02\,\mathrm{GeV}^{-1}$, which is compatible with our
assumptions.444To prevent the W-condensation at the transition temperature,
$eB^{\mathrm{bg}}_{Y}<g_{2}^{2}v(T_{*})^{2}/4$ which sets an upper bound for
the values of $B^{\mathrm{bg}}_{Y}$ that we study in this paper.

Figure 3:  The Higgs bubble profile or the bounce solution that separates the
two phases is shown for three different values of $B^{\mathrm{bg}}_{Y}$.
Figure 4:  The effective potential is shown at transition temperature for
various values of the background hypermagnetic field.

Putting the bounce solution into the Eq. (12) and computing the integral, we
obtain $S_{3}(T)/T$ as a function of temperature. Then, using Eq. (13), we
calculate the transition temperature, $T_{*}$, for different hypermagnetic
field values and show the results in Fig. (4). As seen from Fig. (4), by
increasing $B^{\mathrm{bg}}_{Y}$, the barrier between the vacua is increased
and PT occurs at a lower temperature.
Other important PT characteristics such as latent heat and duration of the PT
are computed at $T_{*}$. The latent heat is given by

$L=-T_{c}\frac{d\Delta V_{eff}(T)}{dT}\Bigg{|}_{T_{c}},$ (15)

from which one can define the vacuum energy density as

$\epsilon_{*}=\Big{(}\Delta V_{eff}(T)-T\frac{d\Delta
V_{eff}(T)}{dT}\Big{)}\Bigg{|}_{T=T_{*}},$ (16)

where $\Delta
V_{eff}(T_{*})=V_{eff}(\phi_{\mathrm{false}}(T_{*}),T_{*})-V_{eff}(\nu(T_{*}),T_{*})$.
Moreover, the duration of the PT, $\tau^{-1}$, can be obtained by the
following procedure. We assume bubbles nucleate with a rate per spacetime
volume which is given by555Here we adhere to the usual definition, with the
understanding that $t$ terminates at $t_{*}$ corresponding to $T_{*}$.
$P=P_{0}\exp(\tau t)$. As a result, we have $\tau=\dot{P}/P$. Then, according
to $dT/dt\simeq-HT$ and Eq. (11), we can find $\tau$ from $S_{3}(T)/T$
function at $T_{*}$ obtained in the previous section [23]

$\frac{\tau}{H_{*}}=T_{*}\frac{d}{dT}\Big{(}\frac{S_{3}(T)}{T}\Big{)}\Big{|}_{T_{*}}.$
(17)

In Table 1, we show values of various quantities characterizing the EWPT,
including $\epsilon_{*}$ and $\tau$, for three values of
$B^{\mathrm{bg}}_{Y}$. It is interesting to note that $\tau^{-1}\sim
10^{-3}(t_{f}-t_{*})$, where $t_{f}$ and $t_{*}$ are the times associated with
$T_{f}$ and $T_{*}$, respectively. In the next section, we will use these
quantities for computing the spectrum of the GWs produced during the EWPT.

$B^{\mathrm{bg}}_{Y}/T^{2}$ | $T_{c}[\mathrm{GeV}]$ | $T_{*}[\mathrm{GeV}]$ | $S_{3}(T_{*})/T_{*}$ | $L[\mathrm{GeV}^{4}]$ | $\epsilon_{*}[\mathrm{GeV}^{4}]$ | $\tau/H_{*}$
---|---|---|---|---|---|---
0.3 | 144.8 | 144.2 | 144.5 | $1.02\times 10^{8}$ | $1.02\times 10^{8}$ | $6.6\times 10^{4}$
0.5 | 137.2 | 138.9 | 144.8 | $1.51\times 10^{8}$ | $1.50\times 10^{8}$ | $2.5\times 10^{4}$
0.7 | 130.4 | 128.1 | 145.0 | $1.87\times 10^{8}$ | $1.84\times 10^{8}$ | $1.5\times 10^{4}$
Table 1: The values of the critical temperature, the transition temperature,
the on-shell bounce action, the latent heat, the vacuum energy density and the
inverse of PT duration are shown for three different values of
$B^{\mathrm{bg}}_{Y}$.

##  3 Gravitational wave generation

In this section, we calculate the spectrum of the GW generated during the
first-order EWPT. Considering the ever increasing detection capabilities of GW
detectors, GW can be used as an effective probe of the early Universe. In
particular, we investigate whether the GW produced in our model falls within
the detection range of Ultimate-DECIGO and BBO correlated interferometers.
During cosmological first-order PTs and evolution of the bubbles of true
vacuum, three processes can give rise to GW radiation. Indeed, when bubbles
nucleate and grow, because of their collisions, part of the latent heat
released during the transition is converted to GWs. Moreover, a fraction of
the energy is transferred to the plasma and causes the plasma motion which in
turn puts forward two other GW sources: MHD turbulence and sound waves.
The contribution of the first mentioned source to the GW frequency spectrum is
calculated by numerical simulations using the envelope approximation and
expressed in terms of the PT parameters [23]:

$h^{2}\Omega_{en}(f)=1.67\times
10^{-5}\Big{(}\frac{0.11v_{b}^{3}}{0.42+v_{b}^{2}}\Big{)}\Big{(}\frac{H_{*}}{\tau}\Big{)}^{2}\Big{(}\frac{\kappa\alpha}{1+\alpha}\Big{)}^{2}\Big{(}\frac{100}{g_{*}}\Big{)}^{\frac{1}{3}}S_{en}(f),$
(18)

where $h$ is the present Hubble parameter $H_{0}$ in units of 100 km
$\mathrm{sec}^{-1}\mathrm{Mpc}^{-1}$, $v_{b}$ is the bubble wall velocity, the
factor $\kappa$ stands for the fraction of the vacuum energy which is
converted into the kinetic energy of the bubbles, $g_{*}\simeq 106$ is the
number of effective relativistic degrees of freedom at the EWPT and $\alpha$
denotes the ratio of the vacuum energy density to the thermal energy density

$\alpha=\frac{\epsilon_{*}}{\frac{\pi^{2}}{30}g_{*}T_{*}^{4}}.$ (19)

The Expression for $\epsilon_{*}$ is given by Eq. (16). The spectral shape of
the GW is given by the following analytic fit [24]

$S_{en}(f)=\frac{3.8(\frac{f}{f_{en}})^{2.8}}{1+2.8(\frac{f}{f_{en}})^{3.8}},$
(20)

where the present-day red-shifted peak frequency is given by the following
relation

$f_{en}=16.5\times
10^{-6}\Big{(}\frac{0.62}{1.8-0.1v_{b}+v_{b}^{2}}\Big{)}\Big{(}\frac{\tau}{H_{*}}\Big{)}\Big{(}\frac{T_{*}}{100~{}\mathrm{GeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}}.$
(21)

Furthermore, we should take into account the other two sources contributing to
the GW energy density. Gravitational wave contributions form sound waves,
which is numerically calculated in [25], and MHD turbulence as a Kolmogorov-
type turbulence modeled by [26] are given by

$h^{2}\Omega_{sw}(f)=2.65\times
10^{-6}\Big{(}\frac{H_{*}}{\tau}\Big{)}\Big{(}\frac{\kappa_{sw}\alpha}{1+\alpha}\Big{)}^{2}\Big{(}\frac{100}{g_{*}}\Big{)}^{\frac{1}{3}}v_{b}~{}S_{sw}(f),$
(22)

and

$h^{2}\Omega_{tu}(f)=3.35\times
10^{-4}\Big{(}\frac{H_{*}}{\tau}\Big{)}\Big{(}\frac{\kappa_{tu}\alpha}{1+\alpha}\Big{)}^{\frac{3}{2}}\Big{(}\frac{100}{g_{*}}\Big{)}^{\frac{1}{3}}v_{b}~{}S_{tu}(f),$
(23)

where their spectral shapes are as follows [27]

$\displaystyle S_{sw}(f)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{f}{f_{sw}}\Big{)}^{3}\Big{(}\frac{7}{4+3(\frac{f}{f_{sw}})^{2}}\Big{)}^{\frac{7}{2}},$
(24) $\displaystyle S_{tu}(f)$ $\displaystyle=$
$\displaystyle\frac{(\frac{f}{f_{tu}})^{3}}{(1+\frac{f}{f_{tu}})^{\frac{11}{3}}(1+\frac{8\pi
f}{h_{*}})},$ (25)

with

$h_{*}=16.5\times
10^{-6}[\mathrm{Hz}]\Big{(}\frac{T_{*}}{100~{}\mathrm{GeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}},$
(26)

as the red-shifted Hubble parameter. The red-shifted peak frequencies in the
spectral shapes of these GW spectra are given by

$\displaystyle f_{sw}=1.9\times
10^{-5}\Big{(}\frac{1}{v_{b}}\Big{)}\Big{(}\frac{\tau}{H_{*}}\Big{)}\Big{(}\frac{T_{*}}{100~{}\mathrm{GeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}},$
$\displaystyle f_{tu}=2.7\times
10^{-5}\Big{(}\frac{1}{v_{b}}\Big{)}\Big{(}\frac{\tau}{H_{*}}\Big{)}\Big{(}\frac{T_{*}}{100~{}\mathrm{GeV}}\Big{)}\Big{(}\frac{g_{*}}{100}\Big{)}^{\frac{1}{6}}.$
(27)

As can be easily seen from Eqs. (18-3), the bubble growth velocity, $v_{b}$,
has an important role in the frequency distribution of GWs originating from
each of the three sources. An important parameter which affects $v_{b}$ is
$\alpha$ as defined in Eq. (19). In particular, a critical value of $\alpha$
is given by [27]

$\alpha_{\infty}=\frac{30}{24\pi^{2}}\frac{\sum_{i}c_{i}\Delta
m_{i}^{2}(\phi_{*})}{g_{*}T_{*}^{2}},$ (28)

where $c_{i}=n_{i}\,(c_{i}=n_{i}/2)$ and $n_{i}$ is the number of degrees of
freedom for boson (fermion) species and $\Delta m_{i}^{2}$ is the squared mass
difference of particles between two phases. Thus, the main contribution to
$\alpha_{\infty}$ comes from the particles becoming heavy during the PT, i.e.
$W$, $Z$ gauge bosons and $t$ quark. For $\alpha<\alpha_{\infty}$, $v_{b}$
remains subluminal and the available energy is transformed into fluid motion.
hence the dominant contributions to GW come from sound waves and MHD
turbulance, i.e. $h^{2}\Omega(f)\simeq h^{2}\Omega_{sw}+h^{2}\Omega_{tu}$.

For $\alpha>\alpha_{\infty}$ the bubbles are “runaway”, that is the excess
vacuum energy density leads to bubble acceleration and $v_{b}$ is bounded only
by the speed of light, $v_{b}=1$. In this case, all three sources contribute
to the GW spectrum, i.e. $h^{2}\Omega(f)\simeq
h^{2}\Omega_{en}+h^{2}\Omega_{sw}+h^{2}\Omega_{tu}$.

In Table 2, we show values of $\alpha$ and $\alpha_{\infty}$ for three
different values of $B^{\mathrm{bg}}_{Y}$. In all of these cases
$\alpha>\alpha_{\infty}$ so that bubbles are runaway, i.e. $v_{b}=1$. In this
case, the efficiency factor for the bubble collision source, which enters
directly in Eq. (18), can be expressed as [27, 28]

$\displaystyle\kappa=1-\frac{\alpha_{\infty}}{\alpha}.$ (29)

In fact the fraction $\alpha_{\infty}/\alpha$ is transformed into the fluid
motion and thermal energy. A fraction $\kappa_{\infty}$ of
$\alpha_{\infty}/\alpha$ is converted to the plasma motion with the efficiency
factor $\kappa_{v}$, and the remaining is spent for reheating the plasma by
the efficiency factor $\kappa_{th}$ [27, 28],

$\displaystyle\kappa_{v}=\frac{\alpha_{\infty}}{\alpha}\kappa_{\infty},~{}~{}~{}~{}~{}~{}~{}~{}\kappa_{th}=\frac{\alpha_{\infty}}{\alpha}(1-\kappa_{\infty}),$
(30)

where

$\displaystyle\kappa_{\infty}=\frac{\alpha_{\infty}}{0.73+0.083\sqrt{\alpha_{\infty}}+\alpha_{\infty}}.$
(31)

In addition, the fraction of plasma motion which is turbulence,
$\varepsilon=\kappa_{tu}/\kappa_{v}$, can be of the order of
$\varepsilon=0.05$ [25]. Therefore, the dominant source coming from the plasma
motion is attributed to the sound waves,
$\kappa_{sw}=(1-\varepsilon)\kappa_{v}$.

$B^{\mathrm{bg}}_{Y}/T_{*}^{2}$ | $\alpha$ | $\alpha_{\infty}$
---|---|---
0.3 | 0.007 | 0.003
0.5 | 0.013 | 0.005
0.7 | 0.020 | 0.008
Table 2: For three different values of $B^{\mathrm{bg}}_{Y}$, $\alpha$ and
$\alpha_{\infty}$ which are necessary to specify GW signals are displayed.

Now, by obtaining the key parameters and determining the contribution of each
GW source, Eqs. (18), (22) and (23), we compute the GW spectrum generated from
the EWPT and show the results for $B^{\mathrm{bg}}_{Y}/T_{*}^{2}=0.7$ in Fig.
(5). As can be seen in Fig. (5), the main contribution comes from the bubble
collisions, followed by that of the sound waves. In Fig. (6), we show the GW
spectrum for the three values of $B^{\mathrm{bg}}_{Y}$ shown in Table 2. As
listed in Table 2, larger $B^{\mathrm{bg}}_{Y}$ gives rise to greater $\alpha$
and so stronger first-order EWPT (see Fig. (4)), leading to GWs with higher
energy density but lower peak frequency (see Fig. (6)).
As can be seen in Fig. (6), these GWs can be detected by future space-based
Ultimate-DECIGO and BBO correlated detectors [17]. Ultimate-DECIGO and BBO
correlated are two future-planned interferometers which will be able to detect
GWs around $0.1-10\,\mathrm{Hz}$ at which noises raised from irresolvable
gravitational wave signals are negligible [29]. Therefore, these detectors
have higher sensitivity and can cover the gap frequency band between LISA and
ground-based detectors. Moreover, the interferometers which consist of several
detectors can enhance their sensitivity by a few orders of magnitude by making
correlation analysis between independent detectors. In fact, Ultimate-DECIGO
and BBO correlated interferometers can reach a sensitivity level of the order
of $10^{-20}$ and $10^{-17}$, respectively, around $0.1\,\mathrm{Hz}$ [18].
Finally, we expect that these GW spectra predicted by our model with peak
frequency between $0.1-1\,\mathrm{Hz}$ can be captured by the most sensitive
frequency range of these two detectors, provided $B^{\mathrm{bg}}_{Y}$ is
large enough.

Figure 5:  We plot the contribution of each source producing GWs during the
EWPT in the presence of the hypermagnetic field
$B^{\mathrm{bg}}_{Y}/T_{*}^{2}=0.7$. The Black solid curve shows the total GW
spectrum. The significant contribution belongs to the bubble collision source,
the red dashed curve. The blue dotdashed and green doted curves show the sound
wave and MHD turbulence contributions, respectively. Figure 6:  We display the
energy density of GWs generated from the first-order EWPT within our model,
for the three values of $B^{\mathrm{bg}}_{Y}$ shown in Tables 1 and 2. Note
that as the value of $B^{\mathrm{bg}}_{Y}$ increases, the energy density of
GWs increases, while the peak frequency decreases. As can be seen, for large
enough $B^{\mathrm{bg}}_{Y}$, these GWs can be in the sensitivity range of the
Ultimate-DECIGO (U-DECIGO) and BBO correlated (BBO Corr.) interferometers.

##  4 Conclusion

We have shown that the presence of a large-scale background hypermagnetic
field, $B^{\mathrm{bg}}_{Y}$, can have profound effects on the EWPT. In
particular, we have shown that as the temperature drops, $B^{\mathrm{bg}}_{Y}$
produces a second minimum at a yet larger value of $\phi$ which eventually
becomes the true vacuum. This makes a two-step transition: a crossover to a
temporary vacuum followed by a first-order to the final vacuum. The barrier
between the two minima also depends on the strength of $B^{\mathrm{bg}}_{Y}$.
We have computed all of the important characteristic quantities for the first-
order part, including the details of the formation and evolution of bubbles of
true vacuum associated with the broken phase, the Higss profile, the duration
of the PT, the latent heat and all of the characteristic temperatures. We find
that by enhancing $B^{\mathrm{bg}}_{Y}$ the transition temperature is lowered
and VEV is increased. We have then explored one of the consequences of the
first-order part of EWPT, which is the generation of GWs. We have shown that
in our model the bubbles are “runaway”, i.e. the bubble wall speed $v_{b}=1$
and thereby all three sources of GWs, i.e. bubble collisions, sound waves and
MHD turbulence contribute. Particularly, we have shown that the bubble
collision is the dominated source of the GWs. We have obtained the energy
density spectrum of GWs as a function of frequency and shown that their
detection is within the range of Ultimate-DECIGO interferometer for
$B^{\mathrm{bg}}_{Y}\gtrsim 0.5T_{*}^{2}$ and BBO correlated interferometer
for $B^{\mathrm{bg}}_{Y}\gtrsim 0.7T_{*}^{2}$, all for peak frequencies of
about $0.1-1\,\mathrm{Hz}$.

Acknowledgments

HA and SSG would like to thank the research council of the Shahid Beheshti
University for financial support. MA would like to thank Yasaman Farzan for
useful comments. MA is supported by Saramadan grant (c/o M. M. Sheikh-Jabbari)
no. ISEF/M/97219.

## References

  * [1] A. Neronov and I. Vovk, Science 328, 73 (2010) [arXiv:1006.3504 [astro-ph.HE]].
  * [2] A. Brandenburg and K. Subramanian, Phys. Rept. 417, 1 (2005) [astro-ph/0405052].
  * [3] M. S. Turner and L. M. Widrow, Phys. Rev. D 37, 2743 (1988); T. Vachaspati, Phys. Lett. B 265, 258 (1991).
  * [4] M. Giovannini and M. E. Shaposhnikov, Phys. Rev. Lett. 80, 22 (1998) [hep-ph/9708303]; M. Giovannini, Phys. Rev. D 61, 063502 (2000) [hep-ph/9906241].
  * [5] M. Giovannini, Phys. Rev. D 92, no. 12, 121301 (2015) [arXiv:1511.00138 [astro-ph.CO]]; S. Rostam Zadeh and S. S. Gousheh, Phys. Rev. D 95, no. 5, 056001 (2017) [arXiv:1607.00650 [hep-ph]].
  * [6] M. Giovannini and M. E. Shaposhnikov, Phys. Rev. D 57, 2186 (1998) [hep-ph/9710234].
  * [7] P. Elmfors, K. Enqvist and K. Kainulainen, Phys. Lett. B 440, 269 (1998) [hep-ph/9806403].
  * [8] D. Comelli, D. Grasso, M. Pietroni and A. Riotto, Phys. Lett. B 458, 304 (1999) [hep-ph/9903227].
  * [9] H. Abedi, M. Ahmadvand and S. S. Gousheh, Phys. Lett. B 786, 35 (2018) [arXiv:1805.10645 [hep-ph]].
  * [10] K. Kajantie, M. Laine, J. Peisa, K. Rummukainen and M. E. Shaposhnikov, Nucl. Phys. B 544, 357 (1999) [hep-lat/9809004].
  * [11] M. Giovannini, Phys. Rev. D 61, 063004 (2000) [hep-ph/9905358]; M. Giovannini, Class. Quant. Grav. 34, no. 13, 135010 (2017) [arXiv:1702.04293 [gr-qc]].
  * [12] A. Ashoorioon and T. Konstandin, JHEP 0907, 086 (2009) [arXiv:0904.0353 [hep-ph]]; R. G. Cai, M. Sasaki and S. J. Wang, JCAP 1708, no. 08, 004 (2017) [arXiv:1707.03001 [astro-ph.CO]]; A. Mazumdar and G. White, arXiv:1811.01948 [hep-ph].
  * [13] M. Ahmadvand and K. Bitaghsir Fadafan, Phys. Lett. B 772, 747 (2017) [arXiv:1703.02801 [hep-th]]; M. Ahmadvand and K. Bitaghsir Fadafan, Phys. Lett. B 779, 1 (2018) [arXiv:1707.05068 [hep-th]].
  * [14] A. Kosowsky, M. S. Turner and R. Watkins, Phys. Rev. Lett. 69, 2026 (1992).
  * [15] A. Kosowsky, A. Mack and T. Kahniashvili, Phys. Rev. D 66, 024030 (2002) [astro-ph/0111483].
  * [16] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014) [arXiv:1304.2433 [hep-ph]].
  * [17] K. Yagi and N. Seto, Phys. Rev. D 83, 044011 (2011) Erratum: [Phys. Rev. D 95, no. 10, 109901 (2017)] [arXiv:1101.3940 [astro-ph.CO]].
  * [18] A. Cooray, Mod. Phys. Lett. A 20, 2503 (2005) [astro-ph/0503118].
  * [19] K. Kamada and A. J. Long, Phys. Rev. D 94, no. 12, 123509 (2016) [arXiv:1610.03074 [hep-ph]].
  * [20] M. Quiros, hep-ph/9901312.
  * [21] A. D. Linde, Phys. Lett. 100B, 37 (1981); A. D. Linde, Nucl. Phys. B 216, 421 (1983) Erratum: [Nucl. Phys. B 223, 544 (1983)].
  * [22] A. Masoumi, K. D. Olum and B. Shlaer, JCAP 1701, no. 01, 051 (2017) [arXiv:1610.06594 [gr-qc]].
  * [23] M. Kamionkowski, A. Kosowsky and M. S. Turner, Phys. Rev. D 49, 2837 (1994) [astro-ph/9310044].
  * [24] S. J. Huber and T. Konstandin, JCAP 0809, 022 (2008) [arXiv:0806.1828 [hep-ph]].
  * [25] M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. D 92, no. 12, 123009 (2015) [arXiv:1504.03291 [astro-ph.CO]].
  * [26] C. Caprini, R. Durrer and G. Servant, JCAP 0912, 024 (2009) [arXiv:0909.0622 [astro-ph.CO]].
  * [27] C. Caprini et al., JCAP 1604, no. 04, 001 (2016) [arXiv:1512.06239 [astro-ph.CO]].
  * [28] J. R. Espinosa, T. Konstandin, J. M. No and G. Servant, JCAP 1006, 028 (2010) [arXiv:1004.4187 [hep-ph]].
  * [29] S. Sato et al., J. Phys. Conf. Ser. 840, no. 1, 012010 (2017).


# Grafting: Making Random Forests Consistent

Nicholas Waltz

###### Abstract

Despite their performance and widespread use, little is known about the theory
of Random Forests. A major unanswered question is whether, or when, the Random
Forest algorithm is consistent. The literature explores various variants of
the classic Random Forest algorithm to address this question and known short-
comings of the method. This paper is a contribution to this literature.
Specifically, the suitability of grafting consistent estimators onto a shallow
CART is explored. It is shown that this approach has a consistency guarantee
and performs well in empirical settings.

##  1 Introduction

This paper studies a variant of the canonical Random Forest [9] that addresses
some of the short-comings commonly associated with the classical method.
Specifically, its inconsistency in some settings and its inability to deal
with high-frequency periodic patterns. It is shown that this variant has a
consistency guarantee and performs well in empirical studies.

Random Forests are a non-parametric ensemble learning method that work by
averaging a large number of tree learners. The standard approach due to
Breiman [9] uses CARTs (Classification and Regression Trees) [11]. They are
known for their combination of fast computation time, interpretability and
good performance across a wide set of applications [16]. In particular, due to
the underlying feature selection in the CART algorithm, they are particularly
suited to high-dimensional settings [27, 12, 18], theoretically even when the
ambient dimension exceeds the sample size. This last aspect sets them apart
from other non-parametric methods and explains their popularity in disciplines
that deal with high-dimensional data, such as public health [21], labour
economics [15, 30], and policy evaluation [1, 22]. It is also one of the
reasons why Hal Varian, (ex) Google chief Economist, singles out Random
Forests in his “Tricks for Econometrics” paper [31].

Yet, despite its popularity, “little is known about the mathematical
properties of the method” [6], and a large body of literature has examined the
algorithm or variants of it under a variety of assumptions [19, 13, 23]. This
paper serves as a contribution to this literature. In particular, an open
question is when and if the Random Forest is consistent. There exists one
major result [27] on the consistency of CART-based Random Forests. This
result, and all variations on it, exists exclusively for the case of additive
models, and when the resampling mechanism of the Random Forest is replaced
with subsampling. It works by showing consistency of the CART algorithm under
strong regularity assumptions. Indeed, further exploration of this shows that
Random Forests perform catastrophically for such models, as the CART has
difficulty learning the additive model well [15] and approaches that
explicitly exploit the additive nature of the distribution would be directly
preferred (see [29]). As of yet, the question, whether the Breiman Random
Forest is universally consistent is an open one. We know, however, that under
the configurations studied in [27] the CART-based random forest is not
universally consistent, which can be shown by means of a clever example [5]
(see Figure 1).

By consistency, here, we refer to $L^{2}$ consistency for non-parametric
estimators, namely, that for a given distribution, the estimator of the
regression function $\hat{m}_{n}(x)$, which is constructed from the sample (in
this case a Random Forest), converges to the conditional expectation function
$m(x)=\mathbb{E}\left[Y|X=x\right]$ in the sense that

$\mathbb{E}\left[m(X)-\hat{m}_{n}(X)\right]^{2}\to 0\text{.}$

Using fairly standard notation, the square here applies to the term inside the
expectation operator. Universal consistency then refers to the case, where
there is such a consistency guarantee regardless of the distribution.

Consistency of the method is important, as it is vital for inference,
especially in settings, where Random Forests are applied to extract causal
relationships (e.g., for estimating treatment effects). For this reason,
consistency is often subsumed in applications. Random Forests, even though the
method is two decades old, are seeing an increased uptake in recent
applications specifically due to their suitability for the relatively novel
Double Machine Learning approach [12], where in fact Random Forests are
provided as the example methodology, because of their fast computation time
and suitability for high dimensions.

Simultaneously, there is a large literature on the short-comings of the
method. The most significant short-coming, known since the inception of the
CART tree learner [11], is the so-called “end cut preference” of the algorithm
that leads to large bias of the algorithm along the edges of the feature space
[17]. In particular, CARTs are unable to pick up on high-frequency periodic
patterns, which is usually not a problem with other non-parametric methods.
Also, inherent to the method is a difficulty with picking up smooth,
continuous patterns, due to the fact that tree learners are piece-wise
constant functions and thus comparable to, e.g., histogram estimators.

Figure 1: Contours of the Biau CEF, the Random Forest estimate, and the
estimate of the variant studied in this paper. The bottom-left and top-right
squares contain infinitely many stripes of decreasing width, so the CART
algorithm will never split the checkerboard pattern in the middle square.

For these reasons, a large literature studies variants of the method that
address these problems. Generally, these fall into one of two categories, (1)
toy models used as a metaphor for the CART-based algorithm and (2) variants
with specific use cases. The most studied variants are the Centered Forest
[10], where splits are performed at the median or midpoint of a node along a
feature chosen at random, the Purely Random Forest [4], where both the feature
and split points are fully random, Mondrian Forests [23], where the tree
learners are constructed using Mondrian partitions, and Kernel Random Forests
[24], where the splitting is replaced with an asymmetric kernel function. This
list goes on (e.g., [4, 7, 20, 32]).

Two other variants of the algorithm that are noteworthy are the Iterative
Random Forest [3], where the method is iterated on iteratively reweighted
features, which allows for the discovery of high-dimensional patterns. And,
while not a Random Forest directly, although a direct connection has been
conjectured [9], Boosting [8], where tree stumps are trained on residuals
until convergence.

This paper explores the suitability of “grafting” consistent estimators onto
the leaf nodes of a shallow CART step, where a particular emphasis will be
placed on grafting on Centered Trees (or “scions”, as it were). As will be
seen,

  1. 1.

This comes with a consistency guarantee of the entire algorithm,

  2. 2.

It outperforms Random Forests in some settings, and otherwise delivers on-par
performance,

  3. 3.

It allows traditionally unsuitable estimators like Kernel Regression to adapt
to high-dimensional settings.

Composition of the paper: After some setup, the main theoretical results are
presented in Section 2, including some discussion over variations of the
algorithm. In Section 3, the algorithm is applied in an empirical setting.
Section 4 presents various numerical experiments, including discussion of the
parameters of the method. Proofs and derivations are presented at the end in
an Appendix for readability.

##  2 Theoretical Results

###  2.1 Tools, Assumptions and Notation

Throughout this paper, we will assume a setting, in which we have some
(unknown) probability measure space
$([0,1]^{p}\times\mathbb{R},\mathcal{F},\mathbb{P})$, where $\mathbb{P}$ is
defined by a distribution, such that

$Y=m(X)+\varepsilon$

with $m(x)=\mathbb{E}\left[Y|X=x\right]$,
$\varepsilon\sim\text{i.i.d.}(0,\sigma^{2})$, and $X\sim\mathcal{U}[0,1]^{p}$.
From this, we draw an i.i.d. sample
$\mathcal{D}_{n}:=\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ of $n$ independent copies of
$(X,Y)$.

It should be noted that these assumptions are fairly standard and that we
impose them for theoretical elegance. Also, it should be pointed out that
under these assumptions, consistency of the Breiman Random Forest has not been
proved and thus cannot be assumed.
Further, we assume,
Assumption 1: $m_{j}(x)=\mathbb{E}\left[Y|X^{(j)}=x\right]$ is
$L_{j}$-Lipschitz.
Assumption 2: We allow $p_{n}\to\infty$, however, we impose sparsity, namely
that $Y$ depends only in $S$ relevant features, such that
$\mathbb{E}\left[Y|X\right]=\mathbb{E}\left[Y|X^{s}\right]$ where $X^{s}$ is
$S$-dimensional and only contains the $S$ relevant features.

Note that Assumption 2 implies that the projection of $m$ onto any irrelevant
feature is
$m_{j}(x)=\mathbb{E}\left[Y|X^{(j)}=x\right]=\mathbb{E}\left[Y\right]$.
Further we assume a functional form $Y(X)$, where $Y(X_{i})=Y_{i}$. Also, let
$\hat{\mathbb{E}}_{n}\left[F(X)\right]=\frac{1}{n}\sum_{i}F(X_{i})$. Lastly,
we refer to the $j$th component of a vector $x$ using the notation $x^{(j)}$.

###  2.2 Random Forests

Given $\mathcal{D}_{n}$, a random forest is an ensemble of recursively grown
partition estimators, so-called trees or tree learners,
$t_{n}(x;\theta,\mathcal{D}_{n})$, where $\theta$ is the (random)
configuration used to grow the tree, which usually contains information on
resampling and the order in which splits are performed. Let
$\mathbf{\Theta}_{n}^{M}:=(\theta^{1}_{n},\theta^{2}_{n},\ldots,\theta^{M}_{n})$
be a vector of $M$ random draws of the configuration vector. The random forest
regressor is then

$\hat{m}_{n,M}(x;\mathbf{\Theta}^{M}_{n},\mathcal{D}_{n})=\frac{1}{M}\sum_{l=1}^{M}{t_{n}(x;\theta^{l}_{n},\mathcal{D}_{n})}\text{.}$

$t_{n}(x;\theta_{n},\mathcal{D}_{n})$ partitions the feature space into so-
called leaf nodes, which are typically hyperrectangles. To each leaf node,
$t_{n}$ assigns a constant value, which serves as the prediction at that leaf
node: specifically, the prediction is the average of the samples that fall
into the leaf node. Let $A(x;\theta_{n},\mathcal{D}_{n})$ be the leaf node
that $x$ falls into. Further, included in $\theta_{n}$ is usually a resampling
step, where a new sample of size $a_{n}\leq n$ is drawn from
$\mathcal{D}_{n}$. We account for this explicitly. Let
$s_{i}(\theta_{n},\mathcal{D}_{n})$ be the sampling function, which is an
integer factor between $0$ and $a_{n}$ measuring how often $(X_{i},Y_{i})$ was
sampled for the construction of a specific tree
$t_{n}(x;\theta_{n},\mathcal{D}_{n})$. It satisfies
$\sum_{i}s_{i}(\theta_{n},\mathcal{D}_{n})=a_{n}$. For sampling without
replacement, $s_{i}(\theta_{n},\mathcal{D}_{n})\in\\{0,1\\}$, for sampling
with replacement $s_{i}(\theta_{n},\mathcal{D}_{n})\in\\{0,1,\ldots,a_{n}\\}$.
Then, it can be shown that the random forest has a local averaging form

$\hat{m}_{n,M}(x;\mathbf{\Theta}^{M}_{n},\mathcal{D}_{n})=\sum_{i}^{n}{W_{i}(x)Y_{i}}\text{,}$

where $W_{i}(x)=\frac{s_{i}\mathbf{1}\\{X_{i}\in
A(x)\\}}{\sum_{j}s_{j}\mathbf{1}\\{X_{j}\in A(x)\\}}$.

For purposes of notational decluttering, where it is obvious, we will often
omit indices and the explicit dependence on the data $\mathcal{D}_{n}$ and
randomisation vector $\theta_{n}$ throughout the paper.

In classification settings, where $Y_{i}\in\\{0,1\\}$, $\hat{m}_{n,M}(x)$ is
then used to construct a classifier

$\hat{g}_{n,M}(x)=\begin{cases}1&\text{if
}\hat{m}_{n,M}(x)>\frac{1}{2}\text{,}\\\ 0&\text{if
}\hat{m}_{n,M}(x)\leq\frac{1}{2}\text{.}\end{cases}$

To show consistency of the classifier, it is sufficient to show consistency of
$\hat{m}_{n,M}(x)$.

As pointed out by Breiman originally [9] and in [25], random forest algorithms
interpolate well. In particular, adding more trees improves performance and we
can show by a strong law of large numbers argument that

$\hat{m}_{n,\infty}(x)=\lim_{M\to\infty}{\hat{m}_{n,M}(x)}=\mathbb{E}_{\theta_{n}}\left[t_{n}(x;\theta_{n},\mathcal{D}_{n})\right]\eqqcolon\hat{m}_{n}(x)\text{.}$

$\mathbb{E}_{\theta_{n}}\left[\cdot\right]$ here refers to the expectation
over the distribution of the configuration vector $\theta_{n}$.
$\hat{m}_{n}(x)$ is the so-called infinite forest. It admits the local
averaging representation

$\hat{m}_{n}(x)=\sum_{i}w_{i}(x)Y_{i}\text{,}$

where we write $w_{i}(x)=\mathbb{E}_{\theta_{n}}\left[W_{i}(x)\right]$. This
is the starting point of our analysis.

####  2.2.1 Tree algorithms

We consider three tree algorithms in this paper.
Algorithm (A): CART [9]:

  1. 1.

Draw $a_{n}\leq n$ samples from $\mathcal{D}_{n}$ with replacement.

  2. 2.

Set the first node $A=[0,1]^{p}$ (root node).

  3. 3.

Split each node $A$ into two daughter nodes $A_{L}$ and $A_{R}$. For this,
choose a set $T$ of $\nu$ test-coordinates uniformly at random from
$\\{1,2,...,p\\}$ and split the node at

$(j^{*},z^{*})=\arg\max_{\\{(j,z):j\in T,z\in
Z^{j}(A,\mathcal{D}_{n})\\}}{L_{n}(j,z;A)}$

Where $Z^{j}(A,\mathcal{D}_{n})$ is the set of midpoints between the samples
in $A$ along feature $j$ and

$\displaystyle L_{n}(j,z;A)$
$\displaystyle=\hat{\mathbb{E}}_{n}\left[(Y(X)-\bar{Y}_{A})^{2}|X\in A\right]$
$\displaystyle-\frac{N(A_{L})}{N(A)}\hat{\mathbb{E}}_{n}\left[(Y(X)-\bar{Y}_{A_{L}})^{2}|X\in
A_{L}\right]$
$\displaystyle-\frac{N(A_{R})}{N(A)}\hat{\mathbb{E}}_{n}\left[(Y(X)-\bar{Y}_{A_{R}})^{2}|X\in
A_{R}\right]$

is the (im)purity gain. Here, $A_{L}=\\{x\in A:x^{(j)}\leq z\\}$ and
$A_{R}=A\backslash A_{L}$,
$N(A)=\sum_{i}\mathbf{1}\\{X_{i}\in A\\}$ and
$\bar{Y}_{A}=\hat{\mathbb{E}}_{n}\left[Y(X)|X\in A\right]$.

  4. 4.

Continue splitting the resulting daughter nodes until (1) all samples in a
node have the same value, (2) there is only one sample in a node, or (3) as
long as the resultant nodes contain more than $q_{n}$ samples, whichever
happens first.

  5. 5.

The final prediction is the mean at each leaf node.

This is the original algorithm due to Breiman. $q_{n}$ here corresponds to the
minpoints variable and $\nu$ to the mtry variable in the original
implementation of the algorithm. Note that in the original paper, $q_{n}=1$,
i.e., we have fully grown trees. Some discussion of $L_{n}(j,z;A)$, the CART-
criterion, is in order at this juncture. Note that inspection of the the CART
criterion reveals an alternate form (see Appendix):

$L_{n}(j,z;A)=\frac{N(A_{L})N(A_{R})}{N(A)^{2}}\left[\bar{Y}_{A_{L}}-\bar{Y}_{A_{R}}\right]^{2}\text{,}$

which is often the preferred implementation, given it requires less
computation. Note that the CART prediction is simply the average over the
samples in each leaf node, i.e., $\bar{Y}_{A}$, where $A$ using $L_{n}$. This
shows that the CART criterion explicitly optimises over the difference between
the resulting predictions of the tree, subject to a reweighting that punishes
splits close to the edge.

We know by a law of large numbers argument that this converges to its
population version almost surely:

$L_{n}(j,z;A)\to\lambda(A_{L})\lambda(A_{R})\left[\mathbb{E}\left[Y|X\in
A_{L}\right]-\mathbb{E}\left[Y|X\in A_{R}\right]\right]^{2}\equiv
L(j,z;A)\text{,}$

where $\lambda$ is the Lebesgue measure. If $L(j,z;A)$ has a unique maximum
$(j^{*},z^{*})$, by reference to empirical process theory we expect
$(j_{n}^{*},z_{n}^{*})\to^{p}(j^{*},z^{*})$, where $(j_{n}^{*},z_{n}^{*})$
refers to the maximiser of $L_{n}(j,z;A)$. This is investigated further in
[2]. This form also reveals that the CART criterion performs splits in the
direction that results in the least amount of variation in the leaf nodes.
Algorithm (B): Centered (median) trees [10]:

  1. 1.

Draw $a_{n}\leq n$ samples from $\mathcal{D}_{n}$ without replacement.

  2. 2.

Set the first node $A=[0,1]^{p}$ (root node).

  3. 3.

Split each node $A$ into two daughter nodes. For this, take the samples that
fall into $A$ and compute the sample median in each feature $\\{1,2,...,p\\}$.
Then, choose a feature uniformly at random (i.e., with probability $1/p$), and
split the node at the sample median in this feature.

  4. 4.

Continue splitting the resultant nodes until (1) all samples in a node have
the same value, (2) there is only one sample in a node, or (3) as long as
there are $q_{n}$ or more samples in the resultant nodes, whichever happens
first.

  5. 5.

The prediction at each node is the mean of the node sample.

This is the algorithm that has been introduced by Breiman in 2004 for its
theoretical elegance, to study a setting, in which Random Forests are
consistent. It is perhaps the most studied variant of the Random Forest [26].
Note how the bootstrap regime is replaced by subsampling in this case. This is
a common feature in most theoretical work on random forests, even in the study
of the CART-based Breiman Random Forests [27], as, unlike bootstrapping, where
the resampling changes the distribution, the subsampled distribution is the
same as that of the original data.

In [13] it is derived that when $m$ is Lipschitz, $X$ is uniform on the
$p$-hypercube, and $Y\in[-D,D]$ almost surely for some $D$, the algorithm is
consistent with rate $n^{-\frac{\log{(1-3/(4p))}}{\log{2}-\log{(1-3/(4p)})}}$,
whenever trees are grown to depth
$\frac{1}{\log{2}-\log{(1-3/(4p))}}\left(\log{n}+C_{1}\right)$ and
$a_{n}>C_{2}n^{\frac{\log{2}}{\log{2}-\log{(1-3/(4p))}}}$, for some constants
$C_{1}$ and $C_{2}$. Here, $\log(\cdot)$ refers to the natural logarithm. Also
note that here, the Centered Forest is parametrised via its depth $k_{n}$
rather than leaf sample size $q_{n}$. That is, we control how often we may
split to get to a leaf node. For Centered Forests, these two notions are
equivalent, as each split effectively halves the sample, so $q_{n}\approx
2^{-k_{n}}a_{n}$.

We now move on to the algorithm studied in this paper:
Algorithm (C): Grafted Trees:

  1. 1.

Draw $a_{n}\leq n$ samples from without replacement.

  2. 2.

Grow a tree using the CART criterion, such that the minimum number of leaf
samples is $\alpha_{n}q_{n}$, where $\alpha_{n}\geq 1$.

  3. 3.

Grow a Centered Forest on the leaf nodes of he CART step as long as there are
at least $q_{n}$ in the resultant nodes.

Algorithm (C*): Generalized Grafting:

  1. 1.

Draw $a_{n}\leq n$ samples from $\mathcal{D}_{n}$ without replacement.

  2. 2.

Grow a tree using the CART criterion, such that the minimum number of leaf
samples is $\alpha_{n}q_{n}$, , where $\alpha_{n}\geq 1$.

  3. 3.

Train a non-parametric regressor $\hat{\mu}_{n}(x)$ on the leaf nodes of the
CART step.

Note that due to its recursive nature, Algorithm (C) has the same computation
time benefits as the algorithm using full CARTs, namely by using
multithreading. In fact, we even expect faster computation, as the random
splits can be generated once beforehand and then just accessed from memory.

(C) and (C*) are the focus and contribution of this paper and we will now turn
to them. Given our tree learner, we are interested in showing
$L^{2}(\mathbb{P})$ consistency:

$\mathbb{E}\left[\hat{m}_{n}(X)-m(X)\right]^{2}\to 0\text{.}$

###  2.3 Grafting, theoretical results

We first show a general bias-variance decomposition theorem for (C) and our
assumptions.

Theorem 1: Under Assumption 1 the $L^{2}(\mathbb{P})$ error of the random
forest regressor defined by Algorithm (C) has the bound

$\mathbb{E}\left[\hat{m}_{n}(X)-m(X)\right]^{2}\leq
p\sum_{j}{L_{j}^{2}\mathbb{E}\left[\ell_{j}(A(X))^{2}\right]}+\tilde{C}\sigma^{2}\frac{a_{n}/q_{n}}{n}\text{,}$
(1)

where $\ell_{j}(A(x))=\sup_{y,z\in A}|y^{(j)}-z^{(j)}|$ is the $j$th side-
length of the leaf node $A(x)$ containing $x$.

Remark 1: In fact, under Assumption 2, we get

$\mathbb{E}\left[\hat{m}_{n}(X)-m(X)\right]^{2}\leq
S\sum_{j\in\mathcal{S}}{L_{j}^{2}\mathbb{E}\left[\ell_{j}(A(X))^{2}\right]}+\tilde{C}\sigma^{2}\frac{a_{n}/q_{n}}{n}\text{,}$

as $L_{j}=0$ for non-relevant features.

We note two things about this. Firstly, the second term, which is the training
error term, i.e., how quickly we “denoise” the sample decays with $n^{-1}$
rather than $a_{n}^{-1}$, which is the rate for a single tree. This is the
consequence of subsampling multiple (infinite) trees. Also, note how equation
(1) gives us a potential starting point for understanding how the Grafted
Forest with Centered “scions” is able to pick up on high-frequency patterns in
the data better than Breiman Random Forests. Through the process of repeatedly
halving the nodes, we have a guarantee that the final side-lengths of the
hyperrectangular partition sets decrease in all features. This is a guarantee
that we don’t generally have with CART [5].

####  2.3.1 Algorithm (C), Grafted Trees

We now move on to deriving consistency for the Grafted Tree algorithm defined
in (C).

Theorem 2: Consistency of Algorithm (C). Under Assumption 1, for $p$ fixed,
and for any choice of $\alpha_{n}$, and $a_{n}/q_{n}$, we have a worst bound
on the $L^{2}(\mathbb{P})$ error of

$\mathbb{E}\left[\hat{m}_{n}(X)-m(X)\right]^{2}\leq
C_{1}p^{2}\left(1-\frac{5}{8p}\right)^{\log_{2}{\alpha_{n}}}+C_{2}\frac{a_{n}/q_{n}}{n}\text{,}$

where $C_{1}$ and $C_{2}$ are independent of $p$ and $n$.

A suitable choice of the parameters $a_{n}$, $q_{n}$ and $\alpha_{n}$
consistency then requires $\frac{a_{n}/q_{n}}{n}\to 0$ and
$\alpha_{n}\to\infty$. E.g. define $\alpha_{n}=2^{\kappa_{n}}$ and
$q_{n}=a_{n}/2^{2\kappa_{n}}$, where $\kappa_{n}=c/2\log_{2}a_{n}$ for $c<1$.
Then provided that $a_{n}^{c}/n\to 0$, this yields consistency, which will
hold for any legal choice of increasing $a_{n}$ with $a_{n}\leq n$.

Note that the rate achieved by this algorithm, albeit a loose bound, beats the
best rate we know for CART-based Random Forests [19], which is of order
$\log{n}^{-1}$.

####  2.3.2 Algorithm (C*)

Theorem 3: Let $\hat{\mu}_{n}$ be consistent with some rate $r(n_{K})\to 0$
over all compact subsets $K$ of $[0,1]^{p}$, where $n_{K}$ is the number of
samples in $K$, such that

$\mathbb{E}\left[(\hat{\mu}_{n}(X)-m(X))^{2}|X\in K\right]\leq
r(n_{K})\text{.}$

Then, Algorithm C* is consistent under our assumptions if
$\alpha_{n}q_{n}\to\infty$.

$r(n)$ here refers to the convergence rate of our estimator of choice, e.g.
for a Nadaraya-Watson estimator with bandwidth $h_{n}\asymp
n^{-\frac{1}{4+\beta}}$ for $\beta>0$, under our assumptions, we would have
$r(n)=n^{-\frac{\beta}{4+\beta}}$ [28]. A note here is that we should expect
the rate of convergence $r(n)$ to decay faster on the leaf nodes in general,
given that the CART step systematically minimises the node variance, which we
in turn expect to lead to a quicker convergence in the training error on the
leaf node estimators.

####  2.3.3 Sparse setting

Another interesting question regards the sparsity setting implied by
Assumption 2. The initial CART step serves both as a detector of discontinuity
and a variance (or variation) reduction tool. We know that it can also serve
as a feature selection tool. Indeed, under the sparsity setting in Assumption
2, it can be shown under some additional smoothness assumptions [27], that for
arbitrary $\xi$, $k$, and some corresponding $N_{\xi}$, whenever the sample
size $n>N_{\xi}$, we know that with probability $1-\xi$, for the splits
$j^{*}_{l}$ with $l<k$ we have

$j^{*}_{l}\in\mathcal{S}\text{.}$

One direct way of using this, would be to train the grafted estimators only on
the features selected in the CART step. If the CART selects the correct
features, we know by Remark 1, that this would imply stable convergence of the
algorithm independent from $p_{n}$ which may increase arbitrarily quickly.
This is explored in Section 4.

##  3 Empirical Application

In this section, the Grafted Tree algorithm is applied to the Boston Housing
dataset, which is a widely used dataset in machine learning and statistics.
It’s a relatively small dataset containing information about homes in the
Boston suburbs. Each entry in the dataset corresponds to aggregated data about
the $n=506$ neighborhoods or housing zones in the Boston area. For each
suburb, 13 features, both numeric and categorical, and the outcome, the median
value of owner-occupied homes, are documented. In more detail, the features
are:

  1. 1.

CRIM: per capita crime rate by town.

  2. 2.

ZN: proportion of residential land zoned for lots over 25,000 sq. ft.

  3. 3.

INDUS: proportion of non-retail business acres per town.

  4. 4.

CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).

  5. 5.

NOX: nitrogen oxide concentration (parts per 10 million).

  6. 6.

RM: average number of rooms per dwelling.

  7. 7.

AGE: proportion of owner-occupied units built prior to 1940.

  8. 8.

DIS: weighted distances to five Boston employment centers.

  9. 9.

RAD: index of accessibility to radial highways.

  10. 10.

TAX: full-value property tax rate per $10,000.

  11. 11.

PTRATIO: pupil-teacher ratio by town.

  12. 12.

B: $1000(Bk-0.63)^{2}$ where $Bk$ is the proportion of people of African
American descent by town.

  13. 13.

LSTAT: % lower status of the population.

These features are then used to predict the outcome MEDV, which is the median
value of owner-occupied homes in $1000s.

For training, the dataset was split into a training set of 404 samples
(roughly 80% of the sample) and test set of 102 samples. We will refer to the
indices of the test set as $I_{test}$ and correspondingly, the indices of the
training set as $I_{train}$.

For all three algorithms a bootstrap sample size of $a_{n}=400$ was chosen and
$M=100$ trees were grown. Three models were considered, the Breiman Random
Forest (Algorithm (A)), and a bagged version, i.e. sampling with replacement
of the Centered Forest (Algorithm (B)) and the Grafted Trees (Algorithm (C)).
The parameters $q_{n}$ and $\alpha_{n}$ (for the Grafted Trees) were
determined via 50-fold cross-validation, using scikit’s RandomSearchCV with
the maximum amount of iterations. For the Breiman Random Forest $\nu$ was set
to 13, i.e. the full feature set, based on the recommendation in [14]. The
result of this was an optimal choice of $q_{n}=1$ for the Breiman Random
Forest (which corresponds to the original method in [9]), $q_{n}=3$ for the
Centered Forest and $q_{n}=5$, $\alpha_{n}=4$ for the Grafted Trees.

Lastly, using the cross-validated parameters, the models were trained and the
test error

$Err_{test}(\hat{m}_{n,M},m)=\frac{1}{|I_{test}|}\sum_{i\in
I_{test}}\left(Y_{i}-\hat{m}_{n,M}(X_{i})\right)^{2}$

was computed. The computed test error was $10.46$ for the Breiman Random
Forest, $11.23$ for the Grafted Tree algorithm and $26.45$ for the Centered
Forest algorithm. From this, we can clearly see that the Grafted Trees
estimator vastly outperforms the Centered Forest, which is a widely-studied
version of the Random Forest which due to its consistency guarantee. Also, the
performance appears to be on-par with the Breiman Random Forest, albeit
slightly worse in the specific test set.

##  4 Experiments

To gain a better understanding of the performance of the algorithm, we now
turn to simulated data. For a conditional expectation function $m(X)$ of
choice, we will sample $n$ data points $(X_{i},Y_{i})$ from a model

$Y=m(X)+\varepsilon\text{,}$

where $\varepsilon\sim N(0,1)$ and $X\sim\mathcal{U}[0,1]^{p}$.

For analysis, the $L^{2}$ error is then computed over a mesh
$\mathcal{X}_{mesh}$ of $[0,1]^{p}$ using

$Err(\hat{m}_{n,M},m)=\frac{1}{|\mathcal{X}_{mesh}|}\sum_{x\in\mathcal{X}_{mesh}}(\hat{m}_{n,M}(x)-m(x))^{2}\text{,}$

where $\hat{m}_{n,M}$ is the estimator of interest.

###  4.1 $L^{2}$ error as a function of $n$

First, we wish to gain an understanding of the performance of the estimator
defined by Algorithm (C) as the sample size increases. To this end, samples of
increasing size were generated for different choices of $m(x)$, and $p=3$ was
fixed. The algorithm was trained using $M=100$ trees and $a_{n}=\lceil
n/1.3\rceil$. $q_{n}$ and $\alpha_{n}$ were chosen via 50-fold cross-
validation. For comparison, the Breiman Random Forest was also trained and its
corresponding $q_{n}$ was also chosen via cross-validation. Lastly, the
$L^{2}$ errors were plotted, alongside the optimal choice of $\alpha_{n}$ and
the respective choices for $q_{n}$ for the Grafted Trees and the Breiman
Random Forest. Some of these experiments are presented now.

Figure 2: Plots for $m(x)=100\sin(200x^{(1)}x^{(2)})$ Figure 3: Plots for
$m(x)=100(x^{(1)})^{4}$ Figure 4: Plots for $m(x)=\cos(30(x^{(3)})^{3})$
Figure 5: Plots for $m(x)=\cos(200x^{(1)}+x^{(2)})+x^{(3)}$ Figure 6: Plots
for $m(x)=x^{(1)}x^{(2)}x^{(3)}$ Figure 7: Plots for
$m(x)=(x^{(1)})^{2}+(x^{(2)})^{3}$

We find that the $L^{2}$ error decreases with sample size and that Grafted
Trees demonstrate similar performance to Breiman Random Forests across the set
of functions considered.

###  4.2 Role of $\alpha_{n}$

If we inspect the optimal value for $\alpha_{n}$ chosen in the cross-
validation step of the experiments above, a common feature is that the optimal
choice of $\alpha_{n}$ increases with sample size. This is promising, given
that this is a necessary condition for consistency of the method. Further, we
regularly encounter values of $\alpha_{n}$ above 10 or more, which seems to
suggest that grafting is preferred over Breiman Random Forests in the
scenarios studied.

To inspect this further, the $L^{2}$ error was examined for different values
of $\alpha_{n}$ for seven functions of choice. For this, $q_{n}$ was set to
$10$; $10,000$ samples were drawn and $100$ trees were grown with $a_{n}=8000$
and only $\alpha_{n}$ was varied. The results of this are presented in Figure
8.

Figure 8: $L^{2}$ error for different values of $\alpha_{n}$, here,
$q_{n}=10$. The list of functions considered is $\sin(200x^{(1)}x^{(2)})$,
$x^{(1)}x^{(2)}$, $x^{(1)}x^{(2)}+x^{(3)}$,
$\sin(x^{(1)}x^{(2)})+\sin(x^{(3)})$,
$-\cos((x^{(1)})^{4}(x^{(2)})^{5}x^{(3)})+0.2(x^{(2)})^{3}$, and
$\sin(x^{(1)}x^{(2)})+x^{(1)}x^{(2)}x^{(3)}$.

Looking at Figure 8, it appears that in the settings studied, we prefer a
shallower CART step as performance gets worse for deeper CART trees.

###  4.3 Feature selection

In this section, the viability of using the CART step to select the relevant
features is examined. To this end, 1000 training samples were drawn from the
distribution defined by

$Y=X^{(1)}X^{(2)}+\varepsilon$

with $\varepsilon\sim N(0,0.1)$ and $X\sim\mathcal{U}[0,1]^{p}$, with $p\geq
2$. Note that this reproduces the sparse setting in Assumption 2, namely with
$S=2$. All values for $p$ ranging from $2$ to $102$ were considered, or
equivalently, between $0$ and $100$ noise features were added. Then, the
$L^{2}$ error was computed for Breiman’s original Random Forests (i.e.
$q_{n}=1$), Grafted Trees with $q_{n}=10$ and $\alpha_{n}=10$, and Centered
Forests with $q_{n}=10$, all of which were trained on the same sample. For
each, $M=100$ trees were grown and $a_{n}=\lceil n/1.3\rceil$ was used. The
results are plotted in Figure 9.

Figure 9: $L^{2}$ error as $p$ increases for the Grafted Trees method vs
Centered Forests and Breiman Random Forests.

We can clearly see that the $L^{2}$ error for both the Breiman Random Forest
and the Grafted Trees is independent of the number of noise features. Indeed,
the Grafted Trees outperform the Breiman Random Forest in this particular
setup. We can also see that the Centered Forest $L^{2}$ error increases
significantly with the number of redundant features. This confirms the
usefulness of using the CART step for feature selection in Algorithm (C).

####  4.3.1 Feature selection in Algorithm (C*)

Similarly, we may interested in the viability of using feature selection for
making consistent algorithms that are unsuitable for high-dimensional settings
adapt to high dimensions by integrating them in an approach that uses
Algorithm (C*). To explore this, the same setup as before was used but with
Algorithm (C*) using a Kernel Ridge estimator on the leaf nodes,
$\alpha_{n}q_{n}=100$, $M=100$ and $a_{n}=\lceil n/1.3\rceil$. This was then
compared to the $L^{2}$ performance of the plain Kernel Ridge estimator
applied to the full dataset and Breiman’s original Random Forest. The results
are presented in Figure 10.

Figure 10: $L^{2}$ error as $p$ increases for the Grafted Kernel Ridge method
vs plain Kernel Ridge Regression and Breiman Random Forests.

We can clearly see that the grafted version of the Kernel Ridge estimator
outperforms the plain version in the high dimensional setting considered and
that its $L^{2}$ error is independent of the number of noise features added.
Indeed, it also outperforms Breiman’s original Random Forest in this setup.

##  5 Conclusion and future work

This paper explored the suitability of using Grafted Trees in prediction
settings. Consistency was proved and the algorithm was empirically examined. A
feature of this method that was not explored in this paper is the fact that in
the limit, we expect tree construction to be independent from the data. Given
that consistency has been shown, this makes the Grafted Trees an attractive
candidate for extracting causal relatitionships, such as treatment effects, as
this independence would obviate the need for, e.g., sample splitting.

## References

  * [1] Susan Athey and Guido Imbens.  Recursive partitioning for heterogeneous causal effects.  Proceedings of the National Academy of Sciences of the United States of America, 113(27):7353–7360, July 2016.
  * [2] Moulinath Banerjee and Ian W. McKeague.  Confidence sets for split points in decision trees.  The Annals of Statistics, 35(2), April 2007.  arXiv:0708.1820 [math, stat].
  * [3] Sumanta Basu, Karl Kumbier, James B. Brown, and Bin Yu.  Iterative random forests to discover predictive and stable high-order interactions.  Proceedings of the National Academy of Sciences, 115(8):1943–1948, February 2018.  Publisher: Proceedings of the National Academy of Sciences.
  * [4] Gerard Biau.  Analysis of a Random Forests Model.  April 2012.
  * [5] Gérard Biau, Luc Devroye, and Gábor Lugosi.  Consistency of Random Forests and Other Averaging Classifiers.  Journal of Machine Learning Research, 9(66):2015–2033, 2008.
  * [6] Gérard Biau and Erwan Scornet.  A random forest guided tour.  TEST, 25(2):197–227, June 2016.
  * [7] Gérard Biau, Erwan Scornet, and Johannes Welbl.  Neural Random Forests, April 2018.  arXiv:1604.07143 [cs, math, stat].
  * [8] Leo Breiman.  Bias, Variance, and Arcing Classifiers.  1996\.
  * [9] Leo Breiman.  Random Forests.  Machine Learning, 45(1):5–32, October 2001.
  * [10] Leo Breiman.  Consistency for a simple model of random forests.  2004\.  Publisher: Citeseer.
  * [11] Leo Breiman, Jerome Friedman, Charles J. Stone, and R. A. Olshen.  Classification and Regression Trees.  CRC Press LLC, Boca Raton, UNITED STATES, 1984.
  * [12] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68, February 2018.
  * [13] Roxane Duroux and Erwan Scornet.  Impact of subsampling and tree depth on random forests.  ESAIM: Probability and Statistics, 22:96–128, 2018.  Publisher: EDP Sciences.
  * [14] Robin Genuer, Jean-Michel Poggi, and Christine Tuleau-Malot.  Variable selection using Random Forests.  Pattern Recognition Letters, 31(14):2225–2236, October 2010.  Publisher: Elsevier.
  * [15] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  Random Forests.  In Trevor Hastie, Robert Tibshirani, and Jerome Friedman, editors, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Springer Series in Statistics, pages 587–604. Springer, New York, NY, 2009.
  * [16] Jeremy Howard and Mike Bowles.  The Two Most Important Algorithms in Predictive Modeling Today. Strata Conference Presentation, New York, February 28 2012. - References - Scientific Research Publishing, 2012.
  * [17] Hemant Ishwaran.  The effect of splitting on random forests.  Machine Learning, 99(1):75–118, April 2015.
  * [18] Jason Klusowski.  Sparse Learning with CART.  In Advances in Neural Information Processing Systems, volume 33, pages 11612–11622. Curran Associates, Inc., 2020.
  * [19] Jason M. Klusowski.  Universal Consistency of Decision Trees in High Dimensions, June 2021.  arXiv:2104.13881 [cs, math, stat].
  * [20] Yi Lin and Yongho Jeon.  Random Forests and Adaptive Nearest Neighbors.  Journal of the American Statistical Association, 101(474):578–590, June 2006.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1198/016214505000001230.
  * [21] Bette Loef, Albert Wong, Nicole A. H. Janssen, Maciek Strak, Jurriaan Hoekstra, H. Susan J. Picavet, H. C. Hendriek Boshuizen, W. M. Monique Verschuren, and Gerrie-Cor M. Herber.  Using random forest to identify longitudinal predictors of health in a 30-year cohort study.  Scientific Reports, 12(1):10372, June 2022.  Number: 1 Publisher: Nature Publishing Group.
  * [22] Shawn McGuire and Charles Delahunt.  Predicting United States Policy Outcomes with Random Forests, 2020.
  * [23] Jaouad Mourtada, Stéphane Gaïffas, and Erwan Scornet.  Minimax optimal rates for Mondrian trees and forests, April 2019.  arXiv:1803.05784 [math, stat].
  * [24] Erwan Scornet.  Learning with random forests.  PhD thesis, Université Pierre et Marie Curie - Paris VI, November 2015.
  * [25] Erwan Scornet.  On the asymptotics of random forests.  Journal of Multivariate Analysis, 146:72–83, April 2016.
  * [26] Erwan Scornet.  Random Forests and Kernel Methods.  IEEE Transactions on Information Theory, 62(3):1485–1500, March 2016.  Conference Name: IEEE Transactions on Information Theory.
  * [27] Erwan Scornet, Gérard Biau, and Jean-Philippe Vert.  Consistency of random forests.  The Annals of Statistics, 43(4):1716–1741, August 2015.  Publisher: Institute of Mathematical Statistics.
  * [28] P. J. Green Silverman, Bernard W.  Nonparametric Regression and Generalized Linear Models: A roughness penalty approach.  Chapman and Hall/CRC, New York, May 1993.
  * [29] Zhiqiang Tan and Cun-Hui Zhang.  Doubly penalized estimation in additive regression with high-dimensional data.  The Annals of Statistics, 47(5):2567–2600, October 2019.  Publisher: Institute of Mathematical Statistics.
  * [30] Gerard J. Van Den Berg, Max Kunaschk, Julia Lang, Gesine Stephan, and Arne Uhlendorff.  Predicting Re-Employment: Machine Learning Versus Assessments by Unemployed Workers and by Their Caseworkers.  SSRN Electronic Journal, 2023.
  * [31] Hal R. Varian.  Big Data: New Tricks for Econometrics.  Journal of Economic Perspectives, 28(2):3–28, May 2014.
  * [32] Zhi-Hua Zhou and Ji Feng.  Deep Forest, July 2020.  arXiv:1702.08835 [cs, stat].

## Appendix

### Alternate form of the CART criterion

$\displaystyle L_{n}(j,z;A)$
$\displaystyle=\hat{\mathbb{E}}_{n}\left[(Y(X)-\bar{Y}_{A})^{2}|X\in
A\right]-\frac{N(A_{L})}{N(A)}\hat{\mathbb{E}}_{n}\left[(Y(X)-\bar{Y}_{A_{L}})^{2}|X\in
A_{L}\right]$
$\displaystyle\>\>\>-\frac{N(A_{R})}{N(A)}\hat{\mathbb{E}}_{n}\left[(Y(X)-\bar{Y}_{A_{R}})^{2}|X\in
A_{R}\right]$ $\displaystyle=\hat{\mathbb{E}}_{n}\left[Y(X)^{2}|X\in
A\right]-\bar{Y}_{A}^{2}-\frac{N(A_{L})}{N(A)}\left(\hat{\mathbb{E}}_{n}\left[Y(X)^{2}|X\in
A_{L}\right]-\bar{Y}_{A_{L}}^{2}\right)$
$\displaystyle\>\>\>-\frac{N(A_{R})}{N(A)}\left(\hat{\mathbb{E}}_{n}\left[Y(X)^{2}|X\in
A_{R}\right]-\bar{Y}_{A_{R}}^{2}\right)$
$\displaystyle=\hat{\mathbb{E}}_{n}\left[Y(X)^{2}|X\in
A\right]-\left(\frac{N(A_{L})}{N(A)}\bar{Y}_{A_{L}}+\frac{N(A_{R})}{N(A)}\bar{Y}_{A_{R}}\right)^{2}$
$\displaystyle\>\>\>-\frac{N(A_{L})}{N(A)}\left(\hat{\mathbb{E}}_{n}\left[Y(X)^{2}|X\in
A_{L}\right]-\bar{Y}_{A_{L}}^{2}\right)$
$\displaystyle\>\>\>-\frac{N(A_{R})}{N(A)}\left(\hat{\mathbb{E}}_{n}\left[Y(X)^{2}|X\in
A_{R}\right]-\bar{Y}_{A_{R}}^{2}\right)$
$\displaystyle=-\left(\frac{N(A_{L})}{N(A)}\bar{Y}_{A_{L}}+\frac{N(A_{R})}{N(A)}\bar{Y}_{A_{R}}\right)^{2}+\frac{N(A_{L})}{N(A)}\bar{Y}_{A_{L}}^{2}+\frac{N(A_{R})}{N(A)}\bar{Y}_{A_{R}}^{2}$
$\displaystyle=\left(1-\frac{N(A_{L})}{N(A)}\right)\frac{N(A_{L})}{N(A)}\bar{Y}_{A_{L}}^{2}-2\frac{N(A_{L})}{N(A)}\frac{N(A_{R})}{N(A)}\bar{Y}_{A_{L}}\bar{Y}_{A_{R}}+\left(1-\frac{N(A_{R})}{N(A)}\right)\frac{N(A_{R})}{N(A)}\bar{Y}_{A_{R}}^{2}$
$\displaystyle=\frac{N(A_{L})}{N(A)}\frac{N(A_{R})}{N(A)}\left(\bar{Y}_{A_{L}}-\bar{Y}_{A_{R}}\right)^{2}\text{.}$

Using

$\hat{\mathbb{E}}_{n}\left[F(X)|X\in
A\right]=\frac{N(A_{L})}{N(A)}\hat{\mathbb{E}}_{n}\left[F(X)|X\in
A_{L}\right]+\frac{N(A_{R})}{N(A)}\hat{\mathbb{E}}_{n}\left[F(X)|X\in
A_{R}\right]$

and

$N(A)=N(A_{L})+N(A_{R})\text{.}$

### Proofs

Theorem 1: We follow a classic bias-variance decomposition, as is usual in
non-parametric theory. Define
$\bar{m}_{n}(x)=\mathbb{E}\left[\hat{m}_{n}(x)|X_{1},X_{2},...,X_{n}\right]$,
i.e. the “denoised” version of the estimator. In particular, note here that

$\bar{m}_{n}(x)=\sum_{i}w_{i}(x)m(X_{i})$

Then note that for the $L^{2}(\mathbb{P})$ error we have

$\mathbb{E}\left[\hat{m}_{n}(X)-m(X)\right]^{2}=\mathbb{E}\left[m(X)-\bar{m}_{n}(X)\right]^{2}+\mathbb{E}\left[\hat{m}_{n}(X)-\bar{m}_{n}(X)\right]^{2}$

The first part we may refer to as the estimation error, the second as the
training error. We know from intuitive reasoning, that the training error
should decay with $\frac{\sigma^{2}}{q_{n}}$. We now consider these two
components in turn:

(1) Estimation error:

$\mathbb{E}\left[m(X)-\bar{m}_{n}(X)\right]^{2}=\mathbb{E}\left[\mathbb{E}_{\theta}\left[\sum_{i}W_{i}(X)[m(X)-m(X_{i})]\right]\right]^{2}$

Consider $\sum_{i}W_{i}(x)[m(X)-m(X_{i})]$ assuming that we are given a
realisation of $\theta$. We have

$\displaystyle\sum_{i}W_{i}(X)[m(X)-m(X_{i})]$ $\displaystyle=\sum_{i\in
I(A(X))}W_{i}(X)[m(X)-m(X_{i})]$ $\displaystyle\leq\sum_{i\in
I(A(X))}W_{i}(X)\sum_{j=1}^{p}{L_{j}|X^{(j)}_{i}-X^{(j)}|}$
$\displaystyle\leq\sum_{i\in I(A(X))}W_{i}(X)\sum_{j=1}^{p}{L_{j}\sup_{x,z\in
A(X)}{|y^{(j)}-z^{(j)}|}}$
$\displaystyle=\sum_{j=1}^{p}{L_{j}\ell_{j}(A(X))}\text{.}$

Here, $I(A(X))$ is a set of indices for the samples falling into $A(X)$. Then,
note that

$\displaystyle\mathbb{E}\left[m(X)-\bar{m}_{n}(X)\right]^{2}$
$\displaystyle=\mathbb{E}\left[\mathbb{E}_{\theta}\left[\sum_{i}W_{i}(X)[m(X)-m(X_{i})]\right]\right]^{2}$
$\displaystyle\leq\mathbb{E}\left[\mathbb{E}_{\theta}\left[\sum_{j=1}^{p}{L_{j}\ell_{j}(A(X))}\right]\right]^{2}$
$\displaystyle\leq
p\sum_{j=1}^{p}L_{j}^{2}\mathbb{E}\left[{\ell_{j}(A(X))}\right]^{2}$
$\displaystyle\leq
p\sum_{j=1}^{p}L_{j}^{2}\mathbb{E}\left[\ell_{j}(A(X))^{2}\right]\text{,}$

using the Cauchy-Schwarz inequality and Jensen’s inequality. Note that in
$\mathbb{E}\left[\ell_{j}(A(X))^{2}\right]$ we take an expectation over all
randomness, the partitioning and and random sampling.

(2) Training error:

$\displaystyle\mathbb{E}\left[\hat{m}_{n}(X)-\bar{m}_{n}(X)\right]^{2}$
$\displaystyle=\mathbb{E}\left[\sum_{i}w_{i}(X)[Y_{i}-m(X_{i})]\right]^{2}$
$\displaystyle=\mathbb{E}\left[\sum_{i}w_{i}(X)[\varepsilon_{i}]\right]^{2}$
$\displaystyle=\mathbb{E}\left[\sum_{i\in
I(A(X))}w_{i}(X)\varepsilon_{i}\right]^{2}$
$\displaystyle\leq\mathbb{E}\left[\sum_{i\in
I(A(X))}w_{i}(X)^{2}\varepsilon_{i}^{2}\right]$
$\displaystyle\leq\sigma^{2}\mathbb{E}\left[\sum_{i\in
I(A(X))}w_{i}(X)^{2}\right]$
$\displaystyle\leq\sigma^{2}\mathbb{E}\left[\sum_{i\in
I(A(X))}w_{i}(X)\max_{i}w_{i}(X)\right]$
$\displaystyle\leq\sigma^{2}\mathbb{E}\left[\max_{i}w_{i}(X)\right]\text{.}$

Now by definition,
$w_{i}(X)=\mathbb{E}_{\theta}\left[\frac{s_{i}\mathbf{1}\\{X_{i}\in
A(X)\\}}{\sum_{j}s_{j}\mathbf{1}\\{X_{j}\in A(X)\\}}\right]$.
We know $\sum_{j}s_{j}\mathbf{1}\\{X_{j}\in A(X)\\}\geq q_{n}$ and
$s_{i}\mathbf{1}\\{X_{i}\in A(X)\\}\leq s_{i}$. Together, this gives
$w_{i}(X)\leq\mathbb{E}_{\theta}\left[s_{i}\right]/\lceil q_{n}\rceil$. Using
this, then

$\displaystyle\mathbb{E}\left[\max_{i}w_{i}(X)\right]$
$\displaystyle\leq\frac{1}{\lceil
q_{n}\rceil}\max_{i}\mathbb{E}_{\theta}\left[s_{i}\right]$
$\displaystyle\leq\tilde{C}\times 1/q_{n}\times a_{n}/n\text{.}$

For which we have used, that the sampling probability is $a_{n}/n$. Thus, we
have proved the theorem. $\square$

Theorem 2: To prove the consistency of the Grafted Tree, we exploit regularity
in the centered “scion” splits. This is similar to [13] with the significant
difference that we need to account for the CART step and leaf sample size.

Consider the successive recursive partition sets containing $X$ as we traverse
the tree. They are hyperrectangles of decreasing side-length. Consider the
distribution of the side length as we traverse the tree. Now, the first splits
use the CART criterion. Let the intermediate side-length in the $j$th feature
after the CART step terminates be distributed as $\mathcal{L}_{CART}^{j}$,
which is some distribution that is a fully determined function of the data,
$\theta$ and $\alpha_{n}q_{n}$. We know that $\mathcal{L}_{CART}^{j}\leq 1$.

After the CART splits, we have independent splits at the median of some
randomly drawn feature. Say, we know exactly, that $n_{k}$ samples are in the
$k$th node of the centered step containing $X$. In particular,
$n_{0}\geq\alpha_{n}q_{n}$ is the leaf sample size. The node has some length
$\mathcal{L}_{k-1}$ and the split is at the median, so if the split occurs
again in this feature, the length is distributed as
$\mathcal{L}_{k-1}\beta(n_{k}+1,n_{k-1}-n_{k})$, where $\beta(a,b)$ is the
beta distribution. Using this, we know that

$\ell_{j}(A(X))\sim\mathcal{L}_{CART}^{j}\times\Pi_{k>\kappa_{n}}\beta(n_{k}+1,n_{k-1}-n_{k})^{\mathbf{1}\\{j^{*}_{k}=j\\}}\text{,}$

where $j^{*}_{k}$ is the feature selected for the $k$th split. Now, as
$n_{k-1}/n_{k}\geq 2$ the expectation clearly vanishes:

$\displaystyle\mathbb{E}\left[\ell_{j}(A(X))\right]$
$\displaystyle=\mathbb{E}\left[\mathcal{L}_{CART}^{j}\right]\times\Pi_{k}\mathbb{E}\left[\beta(n_{k}+1,n_{k-1}-n_{k})^{\mathbf{1}\\{j^{*}_{k}=j\\}}\right]$
$\displaystyle\leq\Pi_{k}\left[\frac{p-1}{p}+\frac{1}{p}\mathbb{E}\left[\beta(n_{k}+1,n_{k-1}-n_{k})\right]\right]$
$\displaystyle=\Pi_{k}\left[\frac{p-1}{p}+\frac{1}{p}\frac{n_{k}+1}{n_{k-1}+1}\right]$
$\displaystyle\leq\Pi_{k}\left[\frac{p-1}{p}+\frac{1}{2p}\right]\xrightarrow[]{k\to\infty}0\text{.}$

We need

$\displaystyle\mathbb{E}\left[\ell_{j}(A(X))^{2}\right]$
$\displaystyle=\mathbb{E}\left[(\mathcal{L}_{CART}^{j})^{2}\right]\times\Pi_{k}\mathbb{E}\left[\beta(n_{k}+1,n_{k-1}-n_{k})^{2\times\mathbf{1}\\{j^{*}_{k}=j\\}}\right]$
$\displaystyle\leq\Pi_{k}\left[\frac{p-1}{p}+\frac{1}{p}\mathbb{E}\left[\beta(n_{k}+1,n_{k-1}-n_{k})^{2}\right]\right]$
$\displaystyle=\Pi_{k}\left[\frac{p-1}{p}+\frac{1}{p}\frac{(n_{k}+1)(n_{k}+2)}{(n_{k-1}+1)(n_{k-1}+2)}\right]$
$\displaystyle\leq\Pi_{k}\left[\frac{p-1}{p}+\frac{1}{4p}\frac{n_{k-1}+4}{n_{k-1}+1}\right]$
$\displaystyle\leq\Pi_{k}\left[\frac{p-1}{p}+\frac{5}{8p}\right]\text{.}$

Clearly, $\left[\frac{p-1}{p}+\frac{5}{8p}\right]<0$. Then,
$\Pi_{k}\left[\frac{p-1}{p}+\frac{5}{8p}\right]\leq\left[1-\frac{5}{8p}\right]^{\lfloor\log_{2}{\alpha_{n}}\rfloor}$.
Together with Theorem 1, we then get the required result.

Theorem 3: Consider a single tree learner. Let $\mathcal{P}$ be the partition
induced in the CART step, we have

$\mathbb{E}\left[t_{n}(X)-m(X)\right]^{2}=\mathbb{E}\left[\mathbb{E}\left[(t_{n}(X)-m(X))^{2}|\mathcal{P}\right]\right]\text{.}$

This partition is completely determined by the data and $\theta$. Then,

$\displaystyle\mathbb{E}\left[(t_{n}(X)-m(X))^{2}|\mathcal{P}\right]$
$\displaystyle=\sum_{A\in\mathcal{P}}\lambda(A)\mathbb{E}\left[(\hat{\mu}_{n}(X;A)-m(X))^{2}|X\in
A\right]$
$\displaystyle\leq\sum_{A\in\mathcal{P}}\lambda(A)\mathbb{E}\left[(\hat{\mu}_{n}(X;A)-m(X))^{2}|X\in\bar{A}\right]$
$\displaystyle\leq r(\alpha_{n}q_{n})\text{,}$

where $\hat{\mu}(x;A)$ is the estimator on the leaf sample of $A$, $\bar{A}$
is the closure of $A$, which is compact.

By extension, then, the average over the $M$ trees is also consistent.