# Microwave fluorescence detection of spin echoes

E. Billaud1, L. Balembois1, M. Le Dantec1, M. Rančić1, E. Albertinale1, S.
Bertaina2, T. Chanelière3, P. Goldner4, D. Estève1, D. Vion1, P. Bertet1∗, E.
Flurin1<EMAIL_ADDRESS><EMAIL_ADDRESS>1Université Paris-
Saclay, CEA, CNRS, SPEC, 91191 Gif-sur-Yvette Cedex, France
2CNRS, Aix-Marseille Université, IM2NP (UMR 7334), Institut Matériaux
Microélectronique et Nanosciences de Provence, Marseille, France
3Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble,
France
4Chimie ParisTech, PSL University, CNRS, Institut de Recherche de Chimie
Paris, 75005 Paris, France

###### Abstract

Counting the microwave photons emitted by an ensemble of electron spins when
they relax radiatively has recently been proposed as a sensitive method for
electron paramagnetic resonance (EPR) spectroscopy, enabled by the development
of operational Single Microwave Photon Detectors (SMPD) at millikelvin
temperature. Here, we report the detection of spin echoes in the spin
fluorescence signal. The echo manifests itself as a coherent modulation of the
number of photons spontaneously emitted after a
$\pi/2_{X}-\tau-\pi_{Y}-\tau-\pi/2_{\Phi}$ sequence, dependent on the relative
phase $\Phi$. We demonstrate experimentally this detection method using an
ensemble of $\mathrm{Er}^{3+}$ ion spins in a scheelite crystal of
$\mathrm{CaWO}_{4}$. We use fluorescence-detected echoes to measure the erbium
spin coherence time, as well as the echo envelope modulation due to the
coupling to the ${}^{183}\mathrm{W}$ nuclear spins surrounding each ion. We
finally compare the signal-to-noise ratio of inductively-detected and
fluorescence-detected echoes, and show that it is larger with the fluorescence
method.

EPR spectroscopy is the method of choice to characterize the concentration and
properties of paramagnetic impurities in a sample [1]. For that goal,
sequences of microwave pulses are applied to the spins, the most widely used
being the Hahn echo ($\pi/2_{X}-\tau-\pi_{Y}-\tau-\mathrm{echo}$) [2]. Applied
to an ensemble of $N$ electron spins $S=1/2$ with an inhomogeneously broadened
Larmor frequency distribution, the Hahn echo sequence generates a transient
build-up of macroscopic transverse magnetization $\langle S_{Y}\rangle$ at a
time $\tau$ after the second pulse, with an amplitude that depends on multiple
factors such as spin density, coherence time, nuclear spin environment, …
Therefore, spin echoes contain important information for EPR spectroscopy.

Echoes are usually detected by coupling the spins to a microwave resonator at
frequency $\omega_{0}$; the Larmor precession of the echo transverse
magnetization induces the emission into an output waveguide of a transient
microwave pulse, which is then amplified and detected. Its amplitude $\langle
Y_{e}\rangle=pN\sqrt{T_{e}\Gamma_{r}}$ depends on the average spin
polarization $p$, the rate at which a spin spontaneously emits a photon into
the waveguide $\Gamma_{r}$, and the echo duration $T_{e}$ [3, 4]. Note that
because the energy stored in the spins is $pN\hbar\omega_{0}/2$, the number of
echo photons $\langle Y_{e}\rangle^{2}$ cannot be larger than $pN/2$, and it
is often much smaller. Enhancing the signal-to-noise ratio (SNR) can be
achieved by using small-mode-volume and high-quality-factor resonators to
increase $\Gamma_{r}$ [5, 6, 7, 8, 9, 10], by cooling the sample far below
$\hbar\omega_{0}/k_{B}$ to polarize the spins ($p=1$) and suppress thermal
noise, and by using quantum-limited amplifiers [3, 11, 12, 13]. Nevertheless,
vacuum fluctuations in the echo detection mode, with standard deviation
$\delta Y_{e}=1/2$, ultimately impose an upper limit to the SNR achievable in
this Inductive-Detection (ID) method.

In this work, we demonstrate an alternative echo detection method, which
avoids this limit and can thus reach a higher SNR than ID. We rotate the
transverse magnetization generated at the echo time into a longitudinal one by
adding a $\pi/2$ pulse applied with a phase $\Phi$ with respect to the first
$\pi/2_{X}$ pulse (called the restoring pulse in the following), similar to
the optical detection of spin echoes [14, 15]. For ideal pulses, the
longitudinal polarization becomes $\langle S_{Z}\rangle=N/2\cos\Phi$. We
detect this coherent modulation by counting the microwave photons emitted
after the sequence when the spins relax radiatively, using a Single Microwave
Photon Detector (SMPD) as demonstrated in [4]. Because the number of echo
photons $N/2$ is much larger than $\langle Y_{e}\rangle^{2}$, and because an
ideal SMPD can detect them noiselessly, the SNR of Fluorescence-Detected (FD)
echo can be much larger than in ID echo. Here, we report the observation of
this FD spin-echo signal, and use it to measure the coherence time of rare-
earth-ions in a crystal as well as the electron-spin-echo envelope modulation
(ESEEM) due to surrounding nuclear spins. We finally compare the echo SNR
reached in FD and ID in the same experimental conditions, and find it is
slightly larger with FD despite the non-idealities of our experiment.

Figure 1:  Experiment principle and spin fluorescence signal. a, Schematic of
the experiment. The resonator of frequency $\omega_{0}/2\pi$ is coupled to
spin $i$ with a strength $g_{0,i}$. It is damped via the transmission line at
a rate $\kappa_{c}$ and via internal losses at a rate $\kappa_{int}$. The
transmission line is used both for exciting the spins (black square pulses)
and for routing the fluorescence signal (red decay) towards the SMPD, which is
read out via another port. b, Sketch of the niobium planar resonator
fabricated on top of the $\mathrm{CaWO}_{4}$ sample in the $ab$ plane. The
wire is along the x axis ($\varphi_{w}=51^{\circ}$) and $B_{0}$ is in the $ab$
plane ($\varphi=31^{\circ}$). The cross section of the sample below the $2\
\mu m$-wide wire (shown as a red rectangle) displays the spatial distribution
of the coupling $g_{0}$. c, Spin detection principle: a pulse of amplitude
$\beta$ and duration $\Delta t$ excites the spins that then relax by emitting
fluorescence photons (red curve), detected and integrated over a period
$T_{int}$. Clicks (top panel-vertical bars) occur with a higher probability at
short times after the excitation, hence the decay of the measured
instantaneous count rate $\langle\dot{C}\rangle$ (bottom panel-circles) down
to the dark count rate $\alpha=5.8\cdot 10^{2}\ \mathrm{s}^{-1}$ (dashed
line). The solid line is the simulation result. d, EPR spectrum of
$\mathrm{Er}^{3+}\mathrm{:CaWO}_{4}$. The average number of counts $\langle
C\rangle$ (open blue circles) is shown as a function of $B_{0}$ (bottom axis),
also converted into frequency detuning (top axis). The solid line is a fit to
a Lorentzian with FWHM $1.6\ \mathrm{MHz}$. Data were obtained with SMPD1.

For our demonstration, we use an ensemble of erbium $\mathrm{Er}^{3+}$ ions in
a scheelite crystal of $\mathrm{CaWO}_{4}$ [16, 17, 18, 19, 20]. The
$\mathrm{Er}^{3+}$ ground state is a degenerate doublet, and behaves as an
effective electron spin $S=1/2$ in the presence of a magnetic field $B_{0}$The
main components of its anisotropic gyromagnetic tensor
$\gamma_{\parallel}/2\pi=17.45$ GHz/T and $\gamma_{\perp}/2\pi=117.3$ GHz/T
depend on whether $B_{0}$ is applied parallel or perpendicular to the crystal
$c$-axis. Here we consider only the $I=0$ erbium isotopes. $\mathrm{CaWO}_{4}$
has a low nuclear-magnetic-moment density, since only the tungsten atoms have
a stable spin-1/2 isotope (${}^{183}\mathrm{W}$, present at $14\%$ natural
abundance), and their gyromagnetic ratio $\gamma_{n}/2\pi=1.8$ MHz/T is low.
The spin properties of paramagnetic rare-earth-ion-doped crystals were
recently studied in the millikelvin regime, and shown to reach coherence times
well above the millisecond [21, 20, 22, 23]. In this work we use the same
sample as in Ref. [20], a pure scheelite crystal with residual
$\mathrm{Er}^{3+}$ concentration $7\,\mu\mathrm{m}^{-3}$ [24], enabling a
quantitative comparison between results obtained by FD and by ID.

As shown in Fig. 1a, the $\mathrm{Er}^{3+}$ ions are inductively coupled to a
planar superconducting LC resonator deposited on top of the $ab$ oriented
crystal surface. The resonator geometry schematically depicted in Fig. 1b
consists of a $2\mu\mathrm{m}$-wide wire that acts as an inductance, in
parallel with an inter-digitated capacitor leading to a frequency
$\omega_{0}/2\pi=6.999$ GHz (resonator 1 in [20]). Because the magnetic field
vacuum fluctuations $\mathbf{\delta B_{1}}(\mathbf{r})$ depend on the position
$\mathbf{r}=(y,z)$ with respect to the inductance, the spin-photon coupling
constant $g_{0}(\mathbf{r})=\mathbf{\delta
B_{1}}(\mathbf{r})\cdot\mathbf{\gamma}\cdot\langle\downarrow|\mathbf{S}|\uparrow\rangle$
is spatially inhomogeneous, reaching its largest value close to the inductance
(see Fig. 1b). The resonator has internal residual losses
($\kappa_{int}=3.6\times 10^{5}\ s^{-1}$) and is capacitively coupled to the
line ($\kappa_{c}=2\times 10^{6}\ s^{-1}$), which yields a total energy
damping rate $\kappa=\kappa_{c}+\kappa_{int}$. We send through a heavily
attenuated input line square-shaped excitation pulses of duration $\Delta t$,
pulsation $\omega_{0}$ and amplitude at the resonator input $\beta$. Each spin
at location $\mathbf{r}$ are driven at the Rabi frequency
$4\beta\sqrt{\kappa_{c}}g_{0}(\mathbf{r})/\kappa$. The reflected pulse,
together with the subsequent spin fluorescence signal, is routed via a
circulator towards a SMPD based on a superconducting transmon qubit [25, 4].
The sample and SMPD are cooled at 10mK in a dilution refrigerator. The SMPD is
operated cyclically, each cycle giving a click ($c=1$) or not ($c=0$). Two
different SMPD devices were used in this work, with different cycle duration
and dark count rates $\alpha$. For SMPD1 (described in [26]), $\alpha=5.8\cdot
10^{2}\ \mathrm{s}^{-1}$ and the cycle duration is $10.2\ \mathrm{\mu s}$. For
SMPD2 (described in [4]), $\alpha=2\cdot 10^{3}\ \mathrm{s}^{-1}$ and the
cycle duration is $8.6\ \mathrm{\mu s}$.

A typical spin fluorescence trace following an excitation pulse applied at
$t=0$ is shown in Fig. 1c. The output $c(t_{j})$ of the SMPD cycle $j$ as a
function of the cycle time $t_{j}$ shows an excess of photons following the
pulse, which decays over a time scale of seconds. By repeating the traces
$\sim 10^{2}$ times, we obtain the average count rate
$\langle\dot{C}(t)\rangle$, shown in Fig. 1c, $C(t)=\sum_{0\leq t_{j}\leq
t}c(t_{j})$ being the number of clicks between $0$ and $t$. Under the
excitation pulse, each spin $i$ undergoes a Rabi rotation leading to a
longitudinal magnetization $S_{Z,i}$. It then emits a photon with probability
$(\Gamma_{r,i}/\Gamma_{1,i})(2S_{Z,i}+1)/2$ at a rate $\Gamma_{1,i}$. Here,
$\Gamma_{1,i}=\Gamma_{sl}+\Gamma_{r,i}$ is the total spin relaxation rate, sum
of the spin-lattice rate $\Gamma_{sl}$ and of the radiative rate
$\Gamma_{r,i}=\kappa g_{0,i}^{2}/(\delta_{i}^{2}+\kappa^{2}/4)$, $g_{0,i}$
being spin $i$ coupling constant and $\delta_{i}$ its detuning to the
resonator [10]. Therefore,
$\langle\dot{C}(t)\rangle=\eta\sum_{i}\Gamma_{r,i}\left[2\langle
S_{Z,i}\rangle+1\right]e^{-\Gamma_{1,i}t}/2$, where $0\leq\eta\leq 1$ is the
probability that the radiative de-excitation of an erbium ion gives rise to a
SMPD click. Because of the inhomogeneity of $g_{0,i}$ (due to the spatial
profile of the resonator mode) and of $\delta_{i}$ (due to the erbium ions
inhomogeneous broadening), both the longitudinal polarization $\langle
S_{Z,i}\rangle$ and the radiative relaxation rate $\Gamma_{r,i}$ are also
widely distributed, explaining why the fluorescence curve in Fig. 1c is non-
exponential. This can be modeled using the known spin-lattice relaxation rate
$\Gamma_{sl}=0.2\leavevmode\nobreak\ \mathrm{s}^{-1}$ [20], the
$\mathrm{Er}^{3+}$ concentration, and the $g_{0,i}$ and $\delta_{i}$
distributions. Simulation results [4, 27] quantitatively reproduce the
observed fluorescence curve as seen in Fig. 1c and yield an efficiency
$\eta=0.12\pm 0.01$. This efficiency corresponds to the product of the
contributions from the measured spin resonator internal losses
($\kappa_{c}/\kappa=0.85$), SMPD duty cycle ($\eta_{d}=0.78$), SMPD efficiency
($\eta_{SMPD}=0.45$), and of the losses between the spin device and the SMPD
(inferred to be $\approx 0.5$).

Figure 2:  Fluorescence detection of spin echo. a, Number of counts $\langle
C\rangle$ as a function of the restoring pulse phase $\Phi$. Blue dots (red
crosses) are experimental (simulated) points, with $\tau=100\ \mu\mathrm{s}$,
$T_{int}=1.52\ \mathrm{s}$, $\beta\approx 10\ \mathrm{ns}^{-1/2}$, $\Delta
t=5\ \mathrm{\mu s}$, and $\eta=0.12$ (from the simulation). The solid line is
a fit yielding $\langle C_{\mathrm{echo}}\rangle=44$ counts and $\langle
C_{\mathrm{bg}}\rangle=730$ counts. b, Measured $\langle
C_{\mathrm{echo}}\rangle$ (open blue dots) as a function of delay $\Delta\tau$
between the restoring pulse and the expected echo time. The solid line is a
Gaussian fit yielding $T_{e}=6.1\mu\mathrm{s}$. Error bars are $1\sigma$
statistical. Data were obtained with SMPD1.

In the limit $\Gamma_{r,i}\gg\Gamma_{sl}$ and $t\gg\Gamma_{r,i}^{-1}$, the
total number of counts is simply $\langle C(t)\rangle=\eta(2\langle
S_{Z}\rangle+1)/2$, where $S_{Z}=\sum S_{Z,i}$ is the spin ensemble
longitudinal magnetization, and is therefore the quantity of interest for
spectroscopy. In practice, an integration window $T_{int}$ is defined, and the
number of counts $\langle C\rangle=\langle\sum_{0\leq t_{i}\leq
T_{\mathrm{int}}}c(t_{i})\rangle-\alpha T_{\mathrm{int}}$ is computed with the
dark count background subtracted. In Fig. 1d, we show $\langle
C\rangle(B_{0})$ around $B_{0}=59$ mT, corresponding to the erbium resonance,
with the field applied at an angle $\varphi=31^{\circ}$ with respect to the
$a$-axis in the $ab$ plane, an angle which minimizes the erbium linewidth
[28]. An approximately Lorentzian lineshape is obtained, with FWHM
$\Gamma/2\pi=1.6$ MHz, in agreement with the results obtained by ID-EPR [20].

We now turn to the fluorescence detection of spin echoes. The sequence
consists of one pulse of amplitude $\beta/2$ around the $X$ axis, followed
after a delay $\tau$ by a second pulse of amplitude $\beta$ around the $Y$
axis, and after another delay $\tau$ by a third pulse of amplitude $\beta/2$
around an axis making a variable angle $\Phi$ with $X$. In Fig. 2a, the number
of counts $\langle C\rangle$ is shown as a function of $\Phi$. The data are
well fitted by $\langle C(\Phi)\rangle=\langle C_{\mathrm{bg}}\rangle+\langle
C_{\mathrm{echo}}\rangle\cos\Phi$, thus displaying the expected FD-echo
coherent modulation. The FD echo amplitude is measured by successively
measuring the number of counts for $\Phi=0$ and $\Phi=\pi$, yielding
$C_{\mathrm{echo}}\equiv[C(0)-C(\pi)]/2$. Varying the delay $\tau+\Delta\tau$
between the second and third pulses, $\langle C_{\mathrm{echo}}\rangle$ shows
a clear echo shape with $T_{e}=6.1\ \mu\mathrm{s}$ (see Fig. 2b).

Contrary to the simple model with ideal pulses, we find $\langle
C_{\mathrm{echo}}\rangle\ll\langle C_{\mathrm{bg}}\rangle$, which calls for an
explanation. We therefore simulate the pulse sequence, and extract $\langle
C\rangle$ by performing the same analysis as with the data. The simulated
$\langle C\rangle(\Phi)$, rescaled by the known efficiency $\eta$, shows an
oscillation with the same phase as the data, an amplitude $\langle
C_{\mathrm{echo}}\rangle$ approximately twice larger (see below), and a
constant background value close to the measured one (see Fig. 2a), indicating
that the latter is due to the spread in Rabi angles among the spin ensemble,
which reduces the echo amplitude as is well-known in pulsed ID-EPR. Such a
spread in Rabi angle could be mitigated in future work using rapid adiabatic
pulses [29, 30, 31, 32].

Figure 3:  Fluorescence detection of the spin coherence time a, Measured
$\langle C_{\mathrm{echo}}\rangle$ (blue open circles) as a function of the
echo delay $2\tau$, with $T_{int}=3.7\ \mathrm{s}$, $\beta\approx 100\
\mathrm{ns}^{-1/2}$ and $\Delta t=5\ \mathrm{\mu s}$. The line is a fit
yielding $T_{2,q}=3.6\,\mathrm{ms}$ and exponent $x_{q}=1.7$. b, Sketched
dependence of the number of counts $\langle C\rangle$ as a function of the
restoring pulse angle $\Phi$ at 3 different times $2\tau$. The orange line
represents $\langle C\rangle$ and the green width the standard deviation at
each angle. The diamonds show the values of $\Phi$ used in our experiment to
compute the standard deviation $\sigma$, also schematically represented here.
c, Measured (open blue diamond) standard deviation $\langle\sigma\rangle$ as a
function of the delay $2\tau$ (see main text). Here, $T_{int}=0.58\mathrm{s}$,
$\beta\approx 10\ \mathrm{ns}^{-1/2}$ and $\Delta t=5\ \mathrm{\mu s}$. The
line is a fit yielding $T_{2}=19.8\,\mathrm{ms}$ and exponent $x=1.2$ d,
Measured $C_{\mathrm{echo}}$ as a function of delay $2\tau$, showing ESEEM.
Here, $T_{int}=0.75\mathrm{s}$, $\beta\approx 10\ \mathrm{ns}^{-1/2}$ and
$\Delta t=5\ \mathrm{\mu s}$. The blue and red open circles are the results
from two successive runs. Error bars are $1\sigma$ statistical. Data were
obtained with SMPD1.

Using the FD echo, we can measure the $\mathrm{Er}^{3+}$ spin coherence time.
First, we measure $\langle C_{\mathrm{echo}}\rangle$ as a function of the
inter-pulse delay $\tau$ (see Fig. 3a). The data are fitted by a stretched
exponential decay $\langle
C_{\mathrm{echo}}\rangle=A\mathrm{e}^{-(2\tau/T_{2,q})^{x_{q}}}$, yielding a
time constant $T_{2,q}=3.6$ ms and exponent $x_{q}=1.7$. This is consistent
with the values measured in [20] using ID-EPR in quadrature-averaged mode.
However, in [20] was found that the spin coherence time was in fact longer
than $T_{2,q}$. Indeed, global magnetic field noise causes the echo transverse
magnetization to acquire an extra phase $\delta\Phi$, which can be modeled as
being Gaussian-distributed with a standard deviation of $2\tau/T_{2,q}$.
Writing the number of counts obtained for a restoring pulse phase $\Phi$ as
$C(\Phi,\delta\Phi,\tau)=C_{\mathrm{bg}}+C_{\mathrm{echo}}\cos(\Phi+\delta\Phi)\mathrm{e}^{-(2\tau/T_{2})^{x}}$,
$T_{2}$ being the spin coherence time and $x$ a streching exponent, we see
indeed that $\langle C_{\mathrm{echo}}\rangle(\tau)=\langle
C_{\mathrm{echo}}\rangle\mathrm{e}^{-(2\tau/T_{2})^{x}}\mathrm{e}^{-(2\tau/T_{2,q})^{2}}$,
thus decaying in a time $T_{2,q}\ll T_{2}$.

On the other hand, it is readily shown that
$\sigma\equiv\sqrt{\mathrm{Var}(C(\Phi,\delta\Phi,\tau))}\approx\sqrt{\mathrm{Var}(C_{\mathrm{bg}})}(1+\frac{C^{2}_{\mathrm{echo}}}{8\mathrm{Var}(C_{\mathrm{bg}})}\mathrm{e}^{-2(2\tau/T_{2})^{x}})$,
where the variance is taken over $\delta\Phi$ and $\Phi$ (see Fig. 3b for a
schematic explanation). Experimentally, we measure the total count
$C(\Phi_{k},\tau)$ as a function of $\tau$ for $\Phi_{k}=k\pi/2$ ($k=0,1,2,3$)
and compute the standard deviation at a given $\tau$ over the 4 angles
$\Phi_{k}$. The data averaged over 950 repetitions (shown in Fig. 3c) is
fitted and yields the spin coherence time $T_{2}=19.8$ ms and stretching
exponent $x=1.2$. The value of $T_{2}$ is in good agreement with the one
obtained by ID-detected echo and magnitude-averaging [33, 34] as well as with
expectations from spectral diffusion caused by the 183W nuclear spin bath
dynamics [20]; the fitted stretching exponent is lower, for unknown reasons.

An important application of spin echoes is to probe the local environment of
paramagnetic species through pulsed hyperfine spectroscopy [1]. We now show
that this is possible also with FD, by measuring the modulation of the
$\mathrm{Er}^{3+}$ echo signal caused by the proximal nuclear spin environment
[16, 35, 36]. In Fig. 3d, $\langle C_{\mathrm{echo}}\rangle$ is shown as a
function of the inter-pulse delay $\tau$ with a $1$ $\mu\mathrm{s}$ step size.
We observe a reproducible modulation, likely due to the hyperfine coupling of
the erbium electron spin to the proximal ${}^{183}\mathrm{W}$ nuclear spins
[37]. We note that the echo modulation also likely explains the 50% reduction
in echo amplitude of the measurements in Fig. 2a with respect to the
simulations (which do not take ESEEM into account), as the chosen inter-pulse
delay $\tau=100\ \mu\mathrm{s}$ in Fig. 2a is close to an ESEEM minimum, as
seen in Fig. 3d.

Figure 4:  SNR comparison between ID and FD detection methods a, Setup 2 used
for this SNR measurement. A JPA is installed behind the SMPD2 for ID echo
measurement. $B_{0}$ is applied with an angle $\varphi=47^{\circ}$. b, ID and
FD echo detection following the exact same two-pulse sequence with
$\beta\approx 5\ \mathrm{ns}^{-1/2}$ and $\Delta t=4\ \mathrm{\mu s}$ to
calibrate the echo photon number. c, Histograms of echo detection with FD
(red, 4000 iterations) and ID (blue, 9000 iterations) taken in the same
conditions, with $\beta\approx 2.5\ \mathrm{ns}^{-1/2}$ and $\Delta t=4\
\mathrm{\mu s}$, and their gaussian fits. Black lines show the distributions
mean. Graphs are aligned on 0 in horizontal axis and are 6 standard deviation
large.

We finally compare the SNR of FD and ID echo. For that, the setup is modified
as shown in Fig. 4a: a non-degenerate Josephson Parametric Amplifier (JPA) is
included in the measurement chain behind the SMPD2, which makes it possible to
measure ID and FD echo under the exact same conditions. We first use the SMPD
to calibrate the photon number in an ID echo, and therefore the echo amplitude
$Y_{e}$. For that, a two-pulse Hahn echo sequence is applied without a
restoring pulse. A time trace of the averaged ID echo shows a Gaussian echo
shape $u(t)$ (see Fig. 4b). Detecting the same echo with the SMPD yields the
corresponding photon number $\langle Y_{e}\rangle^{2}=0.4$, and therefore the
factor needed to convert the integrated heterodyne voltage amplitude $\int
V(t)u(t)dt$ into the dimensionless echo amplitude $Y_{e}$.

We then compare the SNR of ID two-pulse and FD restored-echo sequences, using
identical first-two-pulse parameters, as well as sequence repetition time. To
avoid JPA saturation, we use pulse powers much lower than in Fig. 1
measurements, leading to a larger relaxation rate $\Gamma_{r}\sim
10\mathrm{s}^{-1}$ (data not shown). The resulting echo histograms are shown
in Fig. 4c. We first note that $\langle Y_{e}\rangle^{2}=0.08$, whereas
$\langle C_{echo}\rangle=3$, demonstrating that a larger echo signal can
indeed be obtained with FD than with ID. The ratio of ID to FD echo photon
numbers is $\langle Y_{e}\rangle^{2}/(\langle
C_{echo}\rangle/\eta_{d}\eta_{SMPD})=8\cdot 10^{-3}$, not too far off from the
simple scaling $2N\Gamma_{r}T_{e}=2\cdot 10^{-3}$ predicted for ideal pulses
and no losses, given that $N=2\langle
C_{echo}\rangle/(\eta_{d}\eta_{SMPD})=20$, and $T_{e}=6.1\mu\mathrm{s}$.

The standard deviation of the ID echo is found to be $\delta Y_{e}=1.7$. This
is close to, although larger than, the value $1/\sqrt{2}$ expected from the
sole contribution of vacuum fluctuations and added noise from the quantum-
limited amplifier, likely due to contributions from other amplifiers in the
detection chain, as well as possible over-estimation of the conversion factor.
The standard deviation of the FD echo $\delta C_{\mathrm{echo}}=12$ counts is,
on the other hand, entirely dominated by the dark count contribution. Even
then, the measured FD echo signal-to-noise ratio
$\mathrm{SNR}_{\mathrm{FD}}=\langle C_{\mathrm{echo}}\rangle/2\delta
C_{\mathrm{echo}}=0.12$ is slightly larger than the ID echo
$\mathrm{SNR}_{\mathrm{ID}}=\langle Y_{e}\rangle/2\delta Y_{e}=0.075$. This
result confirms the value of fluorescence detection in terms of sensitivity
when applied to echo detection, even with the present generation of SMPD
devices. Future SMPDs with lower dark count rates would yield further
enhancements. For instance, SMPD1 in the same operating conditions would yield
double the signal-to-noise ratio, given the $1/\sqrt{\alpha}$ SNR scaling [4].
Note that this improvement was not measured with SMPD1 due to the absence of a
parametric amplifier in the spin detection chain.

The detection of spin echoes using fluorescence detection completes the proof-
of-principle results in [4] and establishes FD-EPR as an operational
alternative to ID-EPR. To reach its full potential, FD-EPR requires spins in
the Purcell regime, implying that they should be located within the $\sim pL$
resonator mode volume. Therefore, FD-EPR may prove useful for samples with
small volumes, such as 2D materials. Finally, with larger radiative rate and
lower dark count SMPD devices, FD-EPR has the potential to reach single-spin
sensitivity in the near future.

## Acknowledgements

We acknowledge technical support from P. Sénat, P.-F. Orfila and S. Delprat,
and are grateful for fruitful discussions within the Quantronics group. We
acknowledge IARPA and Lincoln Labs for providing a JTWPA used in the
measurements. We acknowledge support from the Agence Nationale de la Recherche
(ANR) through projects MIRESPIN (ANR-19-CE47-0011), DARKWADOR
(ANR-19-CE47-0004), MARS (ANR-20-CE92-0041), and through the Chaire
Industrielle NASNIQ (ANR-17-CHIN-0001) cofunded by Atos, and from the Région
Ile-De-France through the DIM SIRTEQ (project REIMIC). This project has
received funding from the European Union’s Horizon 2020 research and
innovation program under Marie Sklodowska-Curie grant agreement no. 792727
(SMERC)

## References

  * Schweiger and Jeschke [2001] A. Schweiger and G. Jeschke,  _Principles of pulse electron paramagnetic resonance_ (Oxford University Press, 2001).
  * Hahn [1950] E. L. Hahn, Spin Echoes, Physical Review 80, 580 (1950).
  * Bienfait  _et al._ [2016a] A. Bienfait, J. J. Pla, Y. Kubo, M. Stern, X. Zhou, C. C. Lo, C. D. Weis, T. Schenkel, M. L. W. Thewalt, D. Vion, D. Esteve, B. Julsgaard, K. Moelmer, J. J. L. Morton, and P. Bertet, Reaching the quantum limit of sensitivity in electron spin resonance, Nature Nanotechnology 11, 253 (2016a).
  * Albertinale  _et al._ [2021] E. Albertinale, L. Balembois, E. Billaud, V. Ranjan, D. Flanigan, T. Schenkel, D. Estève, D. Vion, P. Bertet, and E. Flurin, Detecting spins by their fluorescence with a microwave photon counter, Nature 600, 434 (2021).
  * Narkowicz  _et al._ [2008] R. Narkowicz, D. Suter, and I. Niemeyer, Scaling of sensitivity and efficiency in planar microresonators for electron spin resonance, Review of Scientific Instruments 79, 084702 (2008).
  * Shtirberg  _et al._ [2011] L. Shtirberg, Y. Twig, E. Dikarov, R. Halevy, M. Levit, and A. Blank, High-sensitivity Q-band electron spin resonance imaging system with submicron resolution, Review of Scientific Instruments 82, 043708 (2011).
  * Malissa  _et al._ [2013] H. Malissa, D. I. Schuster, A. M. Tyryshkin, A. A. Houck, and S. A. Lyon, Superconducting coplanar waveguide resonators for low temperature pulsed electron spin resonance spectroscopy, Review of Scientific Instruments 84, 025116 (2013).
  * Sigillito  _et al._ [2014] A. J. Sigillito, H. Malissa, A. M. Tyryshkin, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, M. L. W. Thewalt, K. M. Itoh, J. J. L. Morton, A. A. Houck, D. I. Schuster, and S. A. Lyon, Fast, low-power manipulation of spin ensembles in superconducting microresonators, Applied Physics Letters 104, (2014).
  * Artzi  _et al._ [2015] Y. Artzi, Y. Twig, and A. Blank, Induction-detection electron spin resonance with spin sensitivity of a few tens of spins, Applied Physics Letters 106, 084104 (2015).
  * Bienfait  _et al._ [2016b] A. Bienfait, J. Pla, Y. Kubo, X. Zhou, M. Stern, C.-C. Lo, C. Weis, T. Schenkel, D. Vion, D. Esteve, J. Morton, and P. Bertet, Controlling Spin Relaxation with a Cavity, Nature 531, 74 (2016b).
  * Eichler  _et al._ [2017] C. Eichler, A. Sigillito, S. Lyon, and J. Petta, Electron Spin Resonance at the Level of $1{0}^{4}$ Spins Using Low Impedance Superconducting Resonators, Physical Review Letters 118, 037701 (2017).
  * Probst  _et al._ [2017] S. Probst, A. Bienfait, P. Campagne-Ibarcq, J. J. Pla, B. Albanese, J. F. D. S. Barbosa, T. Schenkel, D. Vion, D. Esteve, K. Moelmer, J. J. L. Morton, R. Heeres, and P. Bertet, Inductive-detection electron-spin resonance spectroscopy with 65 spins/Hz^(1/2) sensitivity, Applied Physics Letters 111, 202604 (2017).
  * Ranjan  _et al._ [2020] V. Ranjan, S. Probst, B. Albanese, T. Schenkel, D. Vion, D. Esteve, J. J. L. Morton, and P. Bertet, Electron spin resonance spectroscopy with femtoliter detection volume, Applied Physics Letters 116, 184002 (2020), publisher: American Institute of Physics.
  * Breiland  _et al._ [1973] W. G. Breiland, C. B. Harris, and A. Pines, Optically Detected Electron Spin Echoes and Free Precession in Molecular Excited States, Physical Review Letters 30, 158 (1973), publisher: American Physical Society.
  * Oort  _et al._ [1988] E. v. Oort, N. B. Manson, and M. Glasbeek, Optically detected spin coherence of the diamond N-V centre in its triplet ground state, Journal of Physics C: Solid State Physics 21, 4385 (1988), publisher: IOP Publishing.
  * Mims  _et al._ [1961] W. B. Mims, K. Nassau, and J. D. McGee, Spectral Diffusion in Electron Resonance Lines, Physical Review 123, 2059 (1961).
  * Antipin  _et al._ [1968] A. Antipin, A. Katyshev, I. Kurkin, and L. Shekun, Paramagnetic resonance and spin-lattice relaxation of Er3+ and Tb3+ ions in CaWO4 crystal lattice, Sov. Phys. Solid State 10, 468 (1968).
  * Enrique [1971] B. G. Enrique, Optical Spectrum and Magnetic Properties of Er3+ in CaWO4, The Journal of Chemical Physics 55, 2538 (1971), publisher: American Institute of Physics.
  * Bertaina  _et al._ [2007] S. Bertaina, S. Gambarelli, A. Tkachuk, I. N. Kurkin, B. Malkin, A. Stepanov, and B. Barbara, Rare-earth solid-state qubits, Nature Nanotechnology 2, 39 (2007).
  * [20] M. Le Dantec, M. Rančić, S. Lin, E. Billaud, V. Ranjan, D. Flanigan, S. Bertaina, T. Chanelière, P. Goldner, A. Erb, R. B. Liu, D. Estève, D. Vion, E. Flurin, and P. Bertet, Twenty-three–millisecond electron spin coherence of erbium ions in a natural-abundance crystal, Science Advances 7, eabj9786.
  * Li  _et al._ [2020] P.-Y. Li, C. Liu, Z.-Q. Zhou, X. Liu, T. Tu, T.-S. Yang, Z.-F. Li, Y. Ma, J. Hu, P.-J. Liang, X. Li, J.-Y. Huang, T.-X. Zhu, C.-F. Li, and G.-C. Guo, Hyperfine Structure and Coherent Dynamics of Rare-Earth Spins Explored with Electron-Nuclear Double Resonance at Subkelvin Temperatures, Physical Review Applied 13, 024080 (2020), publisher: American Physical Society.
  * Rančić  _et al._ [2022] M. Rančić, M. L. Dantec, S. Lin, S. Bertaina, T. Chanelière, D. Serrano, P. Goldner, R. B. Liu, E. Flurin, D. Estève, D. Vion, and P. Bertet,  _Electron-spin spectral diffusion in an erbium doped crystal at millikelvin temperatures_, Tech. Rep. arXiv:2203.15012 (arXiv, 2022) arXiv:2203.15012 [cond-mat, physics:quant-ph] type: article.
  * Alexander  _et al._ [2022] J. Alexander, G. Dold, O. W. Kennedy, M. Šimėnas, J. O’Sullivan, C. W. Zollitsch, S. Welinski, A. Ferrier, E. Lafitte-Houssat, T. Lindström, P. Goldner, and J. J. L. Morton,  _Coherent spin dynamics of rare-earth doped crystals in the high-cooperativity regime_, Tech. Rep. arXiv:2206.04027 (arXiv, 2022) arXiv:2206.04027 [quant-ph] type: article.
  * Erb and Lanfranchi [2013] A. Erb and J.-C. Lanfranchi, Growth of high-purity scintillating CaWO4 single crystals for the low-temperature direct dark matter search experiments CRESST-II and EURECA, CrystEngComm 15, 2301 (2013).
  * Lescanne  _et al._ [2020] R. Lescanne, S. Deléglise, E. Albertinale, U. Réglade, T. Capelle, E. Ivanov, T. Jacqmin, Z. Leghtas, and E. Flurin, Irreversible Qubit-Photon Coupling for the Detection of Itinerant Microwave Photons, Physical Review X 10, 021038 (2020).
  * [26] L. Balembois, E. Billaud, and P. Bertet, A new generation of SMPD with a very high sensitivity, To be published .
  * [27] E. Billaud, L. Balembois, and P. Bertet, Fluorescence based spectroscopy of rare earth ions in scheelite, To be published .
  * Mims and Gillen [1966] W. B. Mims and R. Gillen, Broadening of Paramagnetic-Resonance Lines by Internal Electric Fields, Physical Review 148, 438 (1966).
  * Doll and Jeschke [2017] A. Doll and G. Jeschke, Wideband frequency-swept excitation in pulsed EPR spectroscopy, Journal of Magnetic Resonance 280, 46 (2017).
  * Doll  _et al._ [2013] A. Doll, S. Pribitzer, R. Tschaggelar, and G. Jeschke, Adiabatic and fast passage ultra-wideband inversion in pulsed EPR, Journal of Magnetic Resonance 230, 27 (2013).
  * Spindler  _et al._ [2017] P. E. Spindler, P. Schöps, W. Kallies, S. J. Glaser, and T. F. Prisner, Perspectives of shaped pulses for EPR spectroscopy, Journal of Magnetic Resonance Special Issue on Methodological advances in EPR spectroscopy and imaging, 280, 30 (2017).
  * O’Sullivan  _et al._ [2021] J. O’Sullivan, O. W. Kennedy, K. Debnath, J. Alexander, C. W. Zollitsch, M. Šimėnas, A. Hashim, C. N. Thomas, S. Withington, I. Siddiqi, K. Mølmer, and J. J. L. Morton, Random-access quantum memory using chirped pulse phase encoding, arXiv:2103.11697 [quant-ph] (2021), arXiv: 2103.11697.
  * Tyryshkin  _et al._ [2003] A. M. Tyryshkin, S. A. Lyon, A. V. Astashkin, and A. M. Raitsimring, Electron spin relaxation times of phosphorus donors in silicon, Physical Review B 68, 193207 (2003).
  * Tyryshkin  _et al._ [2012] A. M. Tyryshkin, S. Tojo, J. J. L. Morton, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, T. Schenkel, M. L. W. Thewalt, K. M. Itoh, and S. A. Lyon, Electron spin coherence exceeding seconds in high-purity silicon, Nature Materials 11, 143 (2012).
  * Mims [1972] W. B. Mims, Envelope Modulation in Spin-Echo Experiments, Physical Review B 5, 2409 (1972).
  * Mims  _et al._ [1990] W. B. Mims, J. L. Davis, and J. Peisach, The exchange of hydrogen ions and of water molecules near the active site of cytochrome c, Journal of Magnetic Resonance (1969) 86, 273 (1990).
  * Probst  _et al._ [2020] S. Probst, G. Zhang, M. Rančić, V. Ranjan, M. Le Dantec, Z. Zhang, B. Albanese, A. Doll, R. B. Liu, J. Morton, T. Chanelière, P. Goldner, D. Vion, D. Esteve, and P. Bertet, Hyperfine spectroscopy in a quantum-limited spectrometer, Magnetic Resonance 1, 315 (2020).


# GBDT-MO: Gradient Boosted Decision Trees
for Multiple Outputs

Zhendong Zhang and Cheolkon Jung This work was supported by the National
Natural Science Foundation of China (No. 61872280) and the International S&T
Cooperation Program of China (No. 2014DFG12780).Z. Zhang and C. Jung
(corresponding author) are with the School of Electronic Engineering, Xidian
University, Xi’an 710071, China e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Gradient boosted decision trees (GBDTs) are widely used in machine learning,
and the output of current GBDT implementations is a single variable. When
there are multiple outputs, GBDT constructs multiple trees corresponding to
the output variables. The correlations between variables are ignored by such a
strategy causing redundancy of the learned tree structures. In this paper, we
propose a general method to learn GBDT for multiple outputs, called GBDT-MO.
Each leaf of GBDT-MO constructs predictions of all variables or a subset of
automatically selected variables. This is achieved by considering the
summation of objective gains over all output variables. Moreover, we extend
histogram approximation into multiple output case to speed up the training
process. Various experiments on synthetic and real-world datasets verify that
GBDT-MO achieves outstanding performance in terms of both accuracy and
training speed. Our codes are available on-line.

###### Index Terms:

gradient boosting, decision tree, multiple outputs, variable correlations,
indirect regularization.

##  I Introduction

Machine learning and data-driven approaches have achieved great success in
recent years. Gradient boosted decision tree (GBDT) [1] [2] is a powerful
machine learning tool widely used in many applications including multi-class
classification [3], flocculation process modeling [4], learning to rank [5]
and click prediction [6]. It also produces state-of-the-art results for many
data mining competitions such as the Netflix prize [7]. GBDT uses decision
trees as the base learner and sums the predictions of a series of trees. At
each step, a new decision tree is trained to fit the residual between ground
truth and current prediction. GBDT is popular due to its accuracy, efficiency
and interpretability. Many improvements have been proposed after [1]. XGBoost
[8] used the second order gradient to guide the boosting process and improve
the accuracy. LightGBM [9] aggregated gradient information in histograms and
significantly improved the training efficiency. CatBoost [10] proposed a novel
strategy to deal with categorical features.

A limitation of current GBDT implementations is that the output of each
decision tree is a single variable. This is because each leaf of a decision
tree produces a single variable. However, multiple outputs are required for
many machine learning problems including but not limited to multi-class
classification, multi-label classification [11] and multi-output regression
[12]. Other machine learning methods, such as neural networks [13], can adapt
to any dimension of outputs straightforwardly by changing the number of
neurons in the last layer. The flexibility for the output dimension may be one
of the reasons why neural networks are popular. However, it is somewhat
strange to handle multiple outputs by current GBDT implementations. At each
step, they construct multiple decision trees each of which corresponds to an
individual output variable, then concatenates the predictions of all trees to
obtain multiple outputs. This strategy is used in the most popular open-
sourced GBDT libraries: XGBoost [8], LightGBM [9], and CatBoost [10].

The major drawback of the abovementioned strategy is that correlations between
variables are ignored during the training process because those variables are
treated in isolation and they are learned independently. However, correlations
more or less exist between output variables. For example, there are
correlations between classes for multi-class classification. It is verified in
[14] that such correlations improve the generalization ability of neural
networks. Ignoring variable correlations also leads to redundancy of the
learned tree structures. Thus, it is necessary to learn GBDT for multiple
outputs via better strategies. Up to now, a few works have explored it. Geurts
et al. [15] transformed the multiple output problem into a single output
problem by kernelizing the output space. However, this method was not scalable
because the space complexity of its kernel matrix was $n^{2}$ where $n$ is the
number of training samples. Si et al. [16] proposed GBDT for sparse output
(GBDT-sparse). They mainly focused on extreme multi-label classification
problems, and the outputs were represented in sparse format. A sparse split
finding algorithm was designed for square hinge loss. [15] and [16] worked for
specific loss and they did not employ the second order gradient and histogram
approximation.

In this paper, we propose a novel and general method to learn GBDT for
multiple outputs, which is scalable and efficient, named GBDT-MO. Unlike
previous works, we employ the second order gradient and histogram
approximation to improve GBDT-MO. The learning mechanism is designed based on
them to jointly fit all variables in a single tree. Each leaf of a decision
tree constructs multiple outputs at once. This is achieved by maximizing the
summation of objective gains over all output variables. Sometimes, only a
subset of the output variables is correlated. It is expected that the proposed
method automatically selects those variables and constructs predictions for
them at a leaf. We achieve this by adding $L_{0}$ constraint to the objective
function. Since the learning mechanism of GBDT-MO enforces the learned trees
to capture variable correlations, it plays a role in indirect regularization.
Experiments on both synthesis and real-world datasets show that GBDT-MO
achieves better generalization ability than the standard GBDT. Moreover, GBDT-
MO achieves a fast training speed, especially when the number of outputs is
large.

(a) Round 1, variable 1 (b) Round 1, variable 2 (c) Round 1 (d) Round 2,
variable 1

(e) Round 2, variable 2 (f) Round 2

Figure 1: Illustration of the learned trees of GBDT-SO and GBDT-MO on Yeast
with 10 output variables. The first and second columns show GBDT-SO for the
first and second output variables, respectively. The last column shows GBDT-MO
for corresponding variables. We limit tree depth to 3 for illustration.

Compared with existing methods, main contributions of this paper are as
follows:

  * •

We formulate the problem of learning multiple outputs for GBDT, and propose a
split finding algorithm by deriving a general approximate objective for this
problem.

  * •

To learn a subset of outputs, we add a sparse constraint to the objective.
Based on it, we develop two sparse split finding algorithms.

  * •

We extend histogram approximation [17] into multiple output case to speed up
the training process.

The rest of this paper is organized as follows. First, we review GBDT for
single output and introduce basic definitions in Section II. Then, we describe
the details of GBDT-MO in Section III. We address related work in Section IV.
Finally, we perform experiments and conclude in Sections V and VI,
respectively.

##  II GBDT for Single Output

In this section, we review GBDT for single output. First, we describe the work
flow of GBDT. Then show how to derive the objective of GBDT based on the
second order Taylor expansion of the loss, which is used in XGBoost. The
objective to multiple variable cases will be generalized in Section III.
Finally, we explain the split finding algorithms which exactly or
approximately minimize the objective.

###  II-A Work Flow

Denote $\mathcal{D}=\\{(\mathbf{x}_{i},y_{i})_{i=1}^{n}\\}$ as a dataset with
$n$ samples, where $\mathbf{x}\in\mathbb{R}^{m}$ is an $m$ dimension input.
Denote $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ as the function of a decision
tree which maps $\mathbf{x}$ into a scalar. Since GBDT integrates $t$ decision
trees in an additive manner, the prediction of GBDT is
$\hat{y}_{i}=\sum_{k=1}^{t}f_{k}(\mathbf{x}_{i})$, where $f_{k}$ is the
function of $k$-th decision tree. GBDT aims at constructing a series of trees
given datasets. It first calculates gradient based on current prediction at
each boosting round, and then constructs a new tree guided by gradient.
Finally, it updates the prediction using the new tree. The most important part
of GBDT is to construct trees based on gradient.

###  II-B Objective

Based on the construction mechanism of decision trees, $f$ can be further
expressed as follows:

$f(\mathbf{x})=\mathbf{w}_{q(\mathbf{x})},\quad
q:\mathbb{R}^{m}\rightarrow[1,L],\quad\mathbf{w}\in\mathbb{R}^{L}$ (1)

where $L$ is the number of leaves of a decision tree, $q$ is a function which
selects a leaf given $\mathbf{x}$ and $\mathbf{w}_{i}$ is the value of $i$-th
leaf. That is, once a decision tree is constructed, it first maps input into a
leaf, then returns the value of that leaf.

Now, we consider the objective of $(t+1)$-th decision tree given prediction
$\hat{y}$ of the first $t$ tree.

$\sum_{i=1}^{n}l(\hat{y}_{i}+f(\mathbf{x}_{i}),y_{i})+\lambda\mathcal{R}(f)$
(2)

where the first term is fidelity term, $\mathcal{R}$ is regularization term of
$f$. In this work, $\lambda$ is a positive number to control the trade-off
between fidelity term and regularization term. We suppose $l$ is a second
order differentiable loss. Based on the space of $f$, i.e. a constant value
for each leaf, the fidelity term of (2) is separable w.r.t. each leaf. Then,
(2) is rewritten as follows:

$\sum_{j=1}^{L}\left\\{\sum_{i\in
leaf_{j}}l(\hat{y}_{i}+\mathbf{w}_{j},y_{i})\right\\}+\lambda\mathcal{R}(\mathbf{w})$
(3)

Although there are many choices of $\mathcal{R}$, we set
$\mathcal{R}(\mathbf{w})=\frac{1}{2}\lVert\mathbf{w}\rVert_{2}^{2}$, which is
commonly used. Because (3) is _separable_ w.r.t. each leaf, we only consider
the objective of a single leaf as follows:

$\mathcal{L}=\sum_{i}l(\hat{y}_{i}+w,y_{i})+\frac{\lambda}{2}w^{2}$ (4)

where $w$ is the value of a leaf and $i$ is enumerated over the samples
belonging to that leaf. $l(\hat{y}_{i}+w,y_{i})$ can be approximated by the
second order Taylor expansion of $l(\hat{y}_{i},y_{i})$. Then, we have

$\mathcal{L}=\sum_{i}\left\\{l(\hat{y}_{i},y_{i})+g_{i}w+\frac{1}{2}h_{i}w^{2}\right\\}+\frac{\lambda}{2}w^{2}$
(5)

where $g_{i}$ and $h_{i}$ are the first and second order derivatives of
$l(\hat{y}_{i},y)$ w.r.t. $\hat{y}_{i}$. By setting
$\frac{\partial\mathcal{L}}{\partial w}$ to $0$, we obtain the optimal value
of $w$ as follows:

$w^{*}=-\frac{\sum_{i}g_{i}}{\sum_{i}h_{i}+\lambda}$ (6)

Substituting (6) into (5), we get the optimal objective as follows:

$\mathcal{L}^{*}=-\frac{1}{2}\frac{\left(\sum_{i}g_{i}\right)^{2}}{\sum_{i}h_{i}+\lambda}$
(7)

We ignore $l(\hat{y},y)$ since it is a constant term given $\hat{y}$.

###  II-C Split Finding

One of the most important problems in decision tree learning is to find the
best split given a set of samples. Specifically, samples are divided into left
and right parts based on the following rule:

$\mathbf{x}_{i}\in\left\\{\begin{array}[]{cc}left,&\mathbf{x}_{ij}\leq T\\\
right,&\mathbf{x}_{ij}>T\end{array}\right.$ (8)

where $\mathbf{x}_{ij}$ is $j$-th element of $\mathbf{x}_{i}$ and $T$ is the
threshold. The goal of split finding algorithms is to find the best column $j$
and threshold $T$ such that the gain between the optimal objectives before
split and after split is maximized. The optimal objective after split is
defined as the sum of the optimal objectives of left and right parts.

$gain=\mathcal{L}^{*}-(\mathcal{L}^{*}_{left}+\mathcal{L}^{*}_{right})$ (9)

where maximizing $gain$ is equivalent of minimizing
$\mathcal{L}^{*}_{left}+\mathcal{L}^{*}_{right}$ because $\mathcal{L}^{*}$ is
fixed for a given set of samples. $gain$ is used to determine whether a tree
is grown. If $gain$ is smaller than a threshold, we stop the growth to avoid
over-fitting.

Exact and approximate split finding algorithms have been developed. The exact
algorithm enumerates over all possible splits on all columns. For efficiency,
it first sorts the samples according to the values of each column and then
visit the samples in the sorted order to accumulate $g$ and $h$ which are used
to compute the optimal objective. The exact split finding algorithm is
accurate. However, when the number of samples is large, it is time-consuming
to enumerate over all of the possible splits. Approximate algorithms are
necessary as the number of samples increases. The key idea of approximate
algorithms is that they divide samples into buckets and enumerate over these
buckets instead of individual samples. The gradient statistics are accumulated
within each bucket. The complexity of the enumeration process for split is
independent of the number of samples. In literature, there are two strategies
for bucketing: quantile-based and histogram-based. Since the latter is
significantly faster than the former [9], we focus on the latter in this work.
For histogram-based bucketing, a bucket is called a bin. Samples are divided
into $b$ bins by $b$ adjacent intervals: $(s_{0},s_{1},s_{2},\dots,s_{b})$
where $s_{0}$ and $s_{b}$ are usually set to $-\infty$ and $+\infty$
respectively. These intervals are constructed based on the distribution of
$j$th input column of the whole dataset. Once constructed, they are keeping
unchanged. Given a sample $\mathbf{x}_{i}$, it belongs to $k$th bin if and
only if $s_{k-1}<\mathbf{x}_{ij}<=s_{k}$. The bin value of $\mathbf{x}_{ij}$
is obtained by binary-search with complexity $\mathcal{O}(\log b)$. Given bin
values, the histogram of $j$-th input column is constructed by a single fast
scanning. When samples are divided into two parts, it may be unnecessary to
construct the histograms for both parts. One can store the histogram of their
parent node in memory and construct the histogram of one part. Then, the
histogram of another part is obtained by subtracting the constructed part from
the parent histogram. This trick reduces the running time of histogram
construction by at least half.

##  III GBDT for Multiple Outputs

In this section, we describe GBDT-MO in detail. We first formulate the general
problem of learning GBDT for multiple outputs. Specifically, we derive the
objective for learning multiple outputs based on the second order Taylor
expansion of loss. We approximate this objective and connect it with the
objective for single output. We also formulate the problem of learning a
subset of variables and derive its objective. This is achieved by adding
$L_{0}$ constraints. Then, we propose split finding algorithms that minimize
the corresponding objectives. Finally, we discuss our implementations and
analyze the complexity of our proposed split finding algorithms. In this work,
we denote $\mathbf{X}_{i}$ is $i$-th row of a matrix $\mathbf{X}$,
$\mathbf{X}_{.j}$ is $j$-th column and $\mathbf{X}_{ij}$ is its element of
$i$-th row and $j$-th column.

###  III-A Objective

We derive the objective of GBDT-MO. Each leaf of a decision tree constructs
multiple outputs. Denote
$\mathcal{D}=\\{(\mathbf{x}_{i},\mathbf{y}_{i})_{i=1}^{n}\\}$ as a dataset
with $n$ samples, where $\mathbf{x}\in\mathbb{R}^{m}$ is an $m$ dimensional
input and $\mathbf{y}\in\mathbb{R}^{d}$ is a $d$ dimension output instead of a
scalar. Denote $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$ as the function of
a decision tree which maps $\mathbf{x}$ into the output space. Based on the
construction mechanism of decision trees, $f$ can be further expressed as
follows:

$f(\mathbf{x})=\mathbf{W}_{q(\mathbf{x})},\quad
q:\mathbb{R}^{m}\rightarrow[1,L],\quad\mathbf{W}\in\mathbb{R}^{L\times d}$
(10)

where $L$ is the number of leaves of a decision tree, $q$ is a function which
selects a leaf given $\mathbf{x}$ and $\mathbf{W}_{i}\in\mathbb{R}^{d}$ is the
values of $i$-th leaf. That is, once a decision tree is constructed, it first
maps an input into a leaf, then returns the $d$ dimension vector of that leaf.
Next, the prediction of the first $t$ trees is
$\hat{\mathbf{y}}_{i}=\sum_{k=1}^{t}f(\mathbf{x}_{i})$.

We consider the objective of the $(t+1)$-th tree given $\hat{\mathbf{y}}$.
Because it is separable w.r.t. each leaf (see (2) and (3)), we only consider
the objective of a single leaf as follows:

$\mathcal{L}=\sum_{i}l(\hat{\mathbf{y}}_{i}+\mathbf{w},\mathbf{y}_{i})+\lambda\mathcal{R}(\mathbf{w})$
(11)

We highlight that $\mathbf{w}\in\mathbb{R}^{d}$ is a vector with $d$ elements
which belongs to a leaf. Again, we suppose $l$ is a second order
differentiable function. $l(\hat{\mathbf{y}}_{i}+\mathbf{w},\mathbf{y}_{i})$
can be approximated by the second order Taylor expansion of
$l(\hat{\mathbf{y}}_{i},\mathbf{y}_{i})$. Set
$\mathcal{R}(\mathbf{w})=\frac{1}{2}\lVert\mathbf{w}\rVert_{2}^{2}$, we have:

$\mathcal{L}=\sum_{i}\left\\{l(\hat{\mathbf{y}}_{i},\mathbf{y}_{i})+(\mathbf{g})_{i}^{T}\mathbf{w}+\frac{1}{2}\mathbf{w}^{T}(\mathbf{H})_{i}\mathbf{w}\right\\}+\frac{\lambda}{2}\lVert\mathbf{w}\rVert_{2}^{2}$
(12)

where $(\mathbf{g})_{i}=\frac{\partial l}{\partial\hat{\mathbf{y}}_{i}}$ and
$(\mathbf{H})_{i}=\frac{\partial^{2}l}{\partial{\hat{\mathbf{y}}_{i}}^{2}}$.
To avoid notation conflicts with the subscript of vectors or matrices, we use
$(\cdot)_{i}$ to indicate that an object belongs to $i$-th sample. This
notation is omitted when there is no ambiguity. By setting
$\frac{\partial\mathcal{L}}{\partial\mathbf{w}}=\mathbf{0}$ for (12), we
obtain the optimal leaf values:

$\mathbf{w}^{*}=-\left(\sum_{i}(\mathbf{H})_{i}+\lambda\mathbf{I}\right)^{-1}\left(\sum_{i}(\mathbf{g})_{i}\right)$
(13)

where $\mathbf{I}$ is an identity matrix. By substituting $\mathbf{w}^{*}$
into (12) and ignoring the constant term
$l(\hat{\mathbf{y}}_{i},\mathbf{y}_{i})$, we get the optimal objective as
follows:

$\mathcal{L}^{*}=-\frac{1}{2}\left(\sum_{i}(\mathbf{g})_{i}\right)^{T}\left(\sum_{i}(\mathbf{H})_{i}+\lambda\mathbf{I}\right)^{-1}\left(\sum_{i}(\mathbf{g})_{i}\right)$
(14)

We have derived the optimal leaf values and the optimal objective for multiple
outputs. Comparing (13) with (6) and (14) with (7), it is easy to see that
this is a natural generalization of the single output case. In fact, when the
loss function $l$ is separable w.r.t. different output dimensions, or
equivalently, when its hessian matrix $\mathbf{H}$ is diagonal, each element
of $\mathbf{w}^{*}$ is obtained by the same way as in (6).

$\widetilde{\mathbf{w}}^{*}_{j}=-\frac{\sum_{i}(\mathbf{g}_{j})_{i}}{\sum_{i}(\mathbf{h}_{j})_{i}+\lambda}$
(15)

where $\mathbf{h}\in\mathbb{R}^{d}$ is the diagonal elements of $\mathbf{H}$.
And the optimal objective in (14) can be expressed as the sum of objectives
over all output dimensions.

$\widetilde{\mathcal{L}}^{*}=-\frac{1}{2}\sum_{j=1}^{d}\left\\{\frac{\left(\sum_{i}(\mathbf{g}_{j})_{i}\right)^{2}}{\sum_{i}(\mathbf{h}_{j})_{i}+\lambda}\right\\}$
(16)

GBDT-MO and GBDT are different even when $\mathbf{H}$ is diagonal because
GBDT-MO considers the objectives of all output variables at the same time.

However, it is problematic when $l$ is not separable or equivalently
$\mathbf{H}$ is non-diagonal. First, it is difficult to store $\mathbf{H}$ for
every sample when the output dimension $d$ is large. Second, to get the
optimal objective for each possible split, it is required to compute the
inverse of a $d\times d$ matrix, which is time-consuming. Thus, it is
impractical to learn GBDT-MO using the exact objective in (14) and the exact
leaf values in (13). Fortunately, it is shown in [18] that
$\widetilde{\mathbf{w}}^{*}$ and $\widetilde{\mathcal{L}}^{*}$ are good
approximations of the exact leaf values and the exact objective when the
diagonal elements of $\mathbf{H}$ are dominated. $\widetilde{\mathbf{w}}^{*}$
and $\widetilde{\mathcal{L}}^{*}$ are derived from an upper bound of
$l(\hat{\mathbf{y}},\mathbf{y})$. See appendix A for details and further
discussions. Although it is possible to derive better approximations for some
specific loss functions, we use the above diagonal approximation in this work.
We leave better approximations as our future work.

###  III-B Sparse Objective

In Section III-A, we define the objective as the sum of all output variables
because there are correlations among variables. However, only a subset of
variables are correlated in practice, not all variables. Thus, it is required
to learn the values of a suitable subset of $\mathbf{w}$ in a leaf. Although
only a subset of variables is covered in a leaf, all variables can be covered
because there are many leaves. Moreover, sparse objective is able to reduce
the number of parameters. We obtain sparse objective by adding $L_{0}$
constraint to the non-sparse one. Based on the diagonal approximation, the
optimal sparse leaf values are as follows:

$\displaystyle\mathbf{w}^{*}_{sp}=$
$\displaystyle\arg\min_{\mathbf{w}}G^{T}\mathbf{w}+\frac{1}{2}\mathbf{w}^{T}\left(diag(H)+\lambda\mathbf{I}\right)\mathbf{w}$
(17) $\displaystyle s.t.$ $\displaystyle\quad\lVert\mathbf{w}\rVert_{0}<=k$

where $k$ is the maximum non-zero elements of $\mathbf{w}$. We denote
$G=\sum_{i}(\mathbf{g})_{i}$ and $H=\sum_{i}(\mathbf{h})_{i}$ to simplify the
notation. Note the difference between $H$ and $\mathbf{H}$.

For $j$-th element of $\mathbf{w}_{sp}^{*}$, its value is either
$-\frac{G_{j}}{H_{j}+\lambda}$ or 0. Accordingly, the objective contributed by
$j$-th column is either $-\frac{1}{2}\frac{G_{j}^{2}}{H_{j}+\lambda}$ or 0.
Thus, to minimize (17), we select $k$ columns with largest
$v_{j}=\frac{G_{j}^{2}}{H_{j}+\lambda}$. Let $\pi$ be the sorted order of $v$
such that:

$v_{\pi^{-1}(1)}\geq v_{\pi^{-1}(2)}\geq\dots\geq v_{\pi^{-1}(d)}$ (18)

Then, the solution of (17) is as follows:

$(\mathbf{w}_{sp}^{*})_{j}=\left\\{\begin{array}[]{cc}-\frac{G_{j}}{H_{j}+\lambda},&\pi(j)\leq
k\\\ 0,&\pi(j)>k\end{array}\right.$ (19)

That is, columns with $k$ largest $v$ keep their values while others set to
$0$. The corresponding optimal objective is as follows:

$\mathcal{L}^{*}_{sp}=-\frac{1}{2}\sum_{j:\pi(j)\leq
k}\frac{G_{j}^{2}}{H_{j}+\lambda}$ (20)

Our sparse objective is similar to the objective in [16]. Since GBDT-sparse
only uses the first order derivative, the second order derivative is set to a
constant (1.0) for GBDT-sparse. Thus, $H_{j}$ in (20) is the number of samples
belonging to the leaf. The final solution of the sparse objective only
contains a subset of variables. However, to get the optimal subset, all
variables are involved in comparisons. That is, all variables are indirectly
used.

Algorithm 1 Histogram for Multiple Outputs

Input: the set of used samples $\mathcal{S}$, column index $k$, output
dimension $d$ and number of bins $b$

Output: histogram of $k$th column

Initialize Hist.$\mathbf{g}$ $\in\mathbb{R}^{b\times d}$, Hist.$\mathbf{h}$
$\in\mathbb{R}^{b\times d}$ and Hist.$\mathbf{cnt}$ $\in\mathbb{R}^{b}$

for $i$ in $\mathcal{S}$ do

$bin\leftarrow$ bin value of $\mathbf{x}_{ik}$

Hist.$\mathbf{cnt}[bin]$ $\leftarrow$ Hist.$\mathbf{cnt}[bin]$ \+ $1$

for $j=1$ to $d$ do

Hist.$\mathbf{g}[bin][j]$ $\leftarrow$ Hist.$\mathbf{g}[bin][j]$ \+
$(\mathbf{g}_{j})_{i}$

Hist.$\mathbf{h}[bin][j]$ $\leftarrow$ Hist.$\mathbf{h}[bin][j]$ \+
$(\mathbf{h}_{j})_{i}$

end for

end for

Algorithm 2 Approximate Split Finding for Multiple Outputs

Input: histograms of current node, input dimension $m$ and output dimension
$d$

Output: split with maximum $gain$

$gain\leftarrow 0$

for $k=1$ to $m$ do

Hist $\leftarrow$ histogram of $k$th column

$b$ $\leftarrow$ number of bins of Hist

$G\leftarrow\sum_{i=1}^{b}$ Hist.$\mathbf{g}[i]$, $H\leftarrow\sum_{i=1}^{b}$
Hist.$\mathbf{h}[i]$

$G^{l}\leftarrow\mathbf{0}$, $H^{l}\leftarrow\mathbf{0}$

for $i=1$ to $b$ do

$G^{l}\leftarrow G^{l}+$ Hist.$\mathbf{g}[i]$, $H^{l}\leftarrow H^{l}+$
Hist.$\mathbf{h}[i]$

$G^{r}\leftarrow G-G^{l}$, $H^{r}\leftarrow H-H^{l}$

$score\leftarrow\sum_{j=1}^{d}\left\\{\frac{(G_{j}^{l})^{2}}{H_{j}^{l}+\lambda}+\frac{(G_{j}^{r})^{2}}{H_{j}^{r}+\lambda}-\frac{(G_{j})^{2}}{H_{j}+\lambda}\right\\}$,
(16)

$gain\leftarrow$ $\max(gain,score)$

end for

end for

###  III-C Split Finding

Split finding algorithms for single output maximize the gain of objective
before and after split. When dealing with multiple outputs, the objective is
defined over all output variables. To find the maximum gain, it is required to
scan all columns of outputs. The exact algorithm is inefficient because it
enumerates all possible splits. In this work, we use the histogram
approximation based one to speed up the training process. To deal with
multiple outputs, one should extend the histogram construction algorithm from
a single variable case into a multiple variable case. Such an extension is
straightforward. Denote $b$ as the number of bins of a histogram. Then, the
gradient information is stored in a $b\times d$ matrix and its $j$-th column
corresponds to the gradient information of $j$-th variable of outputs. We
describe the histogram construction algorithm for multiple outputs in
Algorithm 1. Once the histogram is constructed, we scan its bins to find the
best split. We describe this histogram based split finding algorithm in
Algorithm 2. Compared with the single output one, the objective gain is the
sum of gains over all outputs.

Algorithm 3 Gain for Sparse Split Finding

Input: gradient statistics of current split, output dimension $d$ and sparse
constraint $k$.

Output: gain of current split

$Q_{l},Q_{r}\leftarrow$ top-k priority queue

for $j=1$ to $d$ do

append $\frac{(G^{l}_{j})^{2}}{H^{l}_{j}+\lambda}$ to $Q_{l}$, append
$\frac{(G^{r}_{j})^{2}}{H^{r}_{j}+\lambda}$ to $Q_{r}$

end for

$score=\sum_{v_{l}\in Q_{l}}v_{l}+\sum_{v_{r}\in Q_{r}}v_{r}$, (21)

Algorithm 4 Gain for Restricted Sparse Split Finding

Input: gradient statistics of current split, output dimension $d$ and sparse
constraint $k$.

Output: gain of current split

$Q\leftarrow$ top-k priority queue

for $j=1$ to $d$ do

append
$\frac{(G^{l}_{j})^{2}}{H^{l}_{j}+\lambda}+\frac{(G^{r}_{j})^{2}}{H^{r}_{j}+\lambda}$
to $Q$

end for

$score=\sum_{v\in Q}v$, (22)

(a) Non-sparse split finding (b) Sparse split finding (c) Restricted sparse
split finding

Figure 2: Examples of different split algorithms for $4$ samples in rows and
$3$ outputs in columns. The sparse constraint is set to $2$. The dashed line
divides samples into two parts. Numbers denote the objectives for each column.
Selected columns of left part are marked in blue and selected columns of right
part are marked in red. Best viewed on the screen.

###  III-D Sparse Split Finding

We propose histogram approximation based sparse split finding algorithms.
Compared with the non-sparse one, the key difference is to compute their
objective gain given a possible split as follows:

$\frac{1}{2}\left\\{\sum_{i:\pi^{l}(i)\leq
k}\frac{(G^{l}_{i})^{2}}{H^{l}_{i}+\lambda}+\sum_{j:\pi^{r}(j)\leq
k}\frac{(G^{r}_{j})^{2}}{H^{r}_{j}+\lambda}\right\\}-const$ (21)

where $\pi^{l}$ is the sorted order of $\frac{(G^{l})^{2}}{H^{l}+\lambda}$ and
$\pi^{r}$ is the sorted order of $\frac{(G^{r})^{2}}{H^{r}+\lambda}$
respectively. $const$ means that the objective before split is fixed for every
possible split. When we scan over columns, we maintain the top-k columns for
both parts whose $\frac{(G)^{2}}{H+\lambda}$ is the largest. This can be
achieved by a top-k priority queue. We describe the algorithm in Algorithm 3.

In Algorithm 3, the sets of the selected columns of left and right parts are
not completely overlapping. We restrict those two sets to be completely
overlapping. In other word, the selected columns of two parts are shared.
Then, the objective gain in such a case becomes:

$\frac{1}{2}\sum_{i:\pi(i)\leq
k}\left\\{\frac{(G^{l}_{i})^{2}}{H^{l}_{i}+\lambda}+\frac{(G^{r}_{i})^{2}}{H^{r}_{i}+\lambda}\right\\}-const$
(22)

where $\pi$ is the sorted order of
$\frac{(G^{l})^{2}}{H^{l}+\lambda}+\frac{(G^{r})^{2}}{H^{r}+\lambda}$. We call
it the restricted sparse split finding algorithm. We describe the gain
computing for it in Algorithm 4. There are two advantages of the restricted
one:

  * •

it has lower computational complexity because it only maintains a single top-k
priority queue.

  * •

it introduces smoothness prior into the function space because it makes two
child nodes with the same parent more similar.

We provide an example to show the differences between non-sparse split
finding, sparse split finding and restricted sparse split finding in Fig. 2.

Algorithm 5 Tree Growth

Input: set of samples $\mathcal{S}$, gradient statistics, and hyper-parameters
for tree learning.

Output: a decision tree

$Hist$ $\leftarrow$ Histograms of $\mathcal{S}$

$Split$ $\leftarrow$ $SplitFinding(Hist)$

$Q\leftarrow$ priority queue sorted by $Split.gain$

$Q.push(Split,Hist,\mathcal{S})$

repeat

$Split,Hist,\mathcal{S}\leftarrow Q.pop()$

$Hist^{l},Hist^{r},\mathcal{S}^{l},\mathcal{S}^{r}\leftarrow
ApplySplit(Split,Hist,\mathcal{S})$

if stop condition is not meet then

$Split^{l}\leftarrow SplitFinding(Hist^{l})$

$Q.push(Split^{l},Hist^{l},\mathcal{S}^{l})$

end if

if stop condition is not meet then

$Split^{r}\leftarrow SplitFinding(Hist^{r})$

$Q.push(Split^{r},Hist^{r},\mathcal{S}^{r})$

end if

until $Q$ is empty

###  III-E Implementation Details

We implement GBDT-MO from scratch by C++. First, we implement the core
function of LightGBM for ourselves, called GBDT-SO. Then, integrate the
learning mechanism for multiple outputs into it, i.e. GBDT-MO111Our codes are
available in https://github.com/zzd1992/GBDTMO.. We also provide a Python
interface. We speed up GBDT-MO using multi-core parallelism, implemented with
OpenMP. A decision tree grows up in the best-first manner. Specifically, we
store the information of nodes that has not been divided in memory. At each
time when we need to add a node, we select the node whose objective gain is
the maximum from all stored nodes. In practice, we store up to $48$ nodes in
memory to reduce the memory cost. We describe the algorithm for growth of a
tree in Algorithm 5.

###  III-F Complexity Analysis

We analyze the training and inference complexity for GBDT-SO and GBDT-MO. For
the training complexity, we focus on the complexity of split finding
algorithms. Recall that the input dimension is $m$ and the output dimension is
$d$, sparse constraint is $k$, and $b$ is the number of bins. Suppose that $b$
is fixed for all input dimensions, $t$ is the number of boosting rounds and
$h$ is the maximum tree depth. Table I shows the training complexity. The
complexity of non-sparse split finding is $\mathcal{O}(bmd)$ because it
enumerates on histogram bins, input dimensions and output dimension. For the
sparse case, it is multipled by $\log k$ because inserting an element into a
top-k priority queue requires $\mathcal{O}(\log k)$ comparisons (see Algorithm
4). When the exact hessian is used (see (14)), it becomes
$\mathcal{O}(bmd^{3})$ because calculating the inverse of a $d\times d$ matrix
requires $\mathcal{O}(d^{3})$ operations. GBDT-SO has the same complexity as
GBDT-MO for split finding. However, it does not mean that GBDT-SO is as fast
as GBDT-MO. Beyond split finding, GBDT-SO has more overhead that slows its
training speed down. For example, samples are divided into two parts, i.e. $d$
times for GBDT-SO, but once for GBDT-MO. Thus, the training speed of GBDT-MO
is faster than GBDT-SO in practice. Table II shows the inference complexity.
It can be observed that GBDT-SO needs more comparisons than GBDT-MO.

TABLE I: Complexity of split finding. Methods | Complexity
---|---
GBDT-SO | $\mathcal{O}(bmd)$
GBDT-MO | $\mathcal{O}(bmd)$
GBDT-MO (sparse) | $\mathcal{O}(bmd\log k)$
GBDT-MO (exact hessian) | $\mathcal{O}(bmd^{3})$
TABLE II: Complexity of inference. Methods | Comparisons | Additions
---|---|---
GBDT-SO | $\mathcal{O}(tdh)$ | $\mathcal{O}(td)$
GBDT-MO | $\mathcal{O}(th)$ | $\mathcal{O}(td)$
GBDT-MO (sparse) | $\mathcal{O}(th)$ | $\mathcal{O}(tk)$

##  IV Related Work

Gradient boosted decision tree (GBDT) proposed in [1] has received much
attention due to its accuracy, efficiency and interpretability. GBDT has two
characteristics: it uses decision trees as the base learner and its boosting
process is guided by the gradient of some loss function. Many variants of [1]
have been proposed. Instead of the first order gradient, XGBoost [8] also uses
the second order gradient to guide its boost process and derives the
corresponding objective for split finding. The histogram approximation of
split finding is proposed in [17] which is used as the base algorithm in
LightGBM [9]. Because the second order gradient improves the accuracy and
histogram approximation improves the training efficiency, those two
improvements are also used in the proposed GBDT-MO.

Since machine learning problems with multiple outputs become common, many tree
based or boosting based methods have been proposed to deal with multiple
outputs. [19] [20] generalize the impurity measures defined for binary
classification and ranking tasks to a multi-label scenario for splitting a
node. However, they are random forest based methods. That is, new trees are
not constructed in a boosting manner. Several works extend adaptive boost
(AdaBoost) into multi-label cases such as AdaBoost.MH [21] and AdaBoost.LC
[22]. The spirits of AdaBoost.MH and AdaBoost.LC are different from GBDT. At
each step, a new base learner is trained from scratch on the re-weighted
samples. Moreover, AdaBoost.MH only works for Hamming loss, while AdaBoost.LC
only works for the covering loss [22].

They do not belong to GBDT families. Two works which belong to GBDT families
have been proposed for learning multiple outputs [15] [16]. [15] transforms
the multiple output problem into the single output problem by kernelizing the
output space. To achieve this, an $n\times n$ kernel matrix should be
constructed where $n$ is the number of training samples. Thus, this method is
not scalable. Moreover, it works only for square loss. GBDT for sparse output
(GBDT-sparse) is proposed in [16]. The outputs are represented in sparse
format. A sparse split finding algorithm is designed by adding $L_{0}$
constraint to the objective. The sparse split finding algorithms of GBDT-MO
are inspired by this work. There are several differences between GBDT-MO and
GBDT-sparse:

  * •

GBDT-sparse focuses on extreme multi-label classification problems, whereas
GBDT-MO focuses on general multiple output problems. GBDT-sparse requires the
loss is separable over output variables. It also requires its gradient is
sparse, i.e. $\frac{\partial l(\hat{y},y)}{\partial\hat{y}}=0$ if $\hat{y}=y$.
During training, it introduces the clipping operator into the loss to maintain
the sparsity of gradient. The facts limit the types of loss.

  * •

GBDT-sparse does not employ the second order gradient. Its objective for split
finding is derived based on the first order Taylor expansion of the loss as in
[1].

  * •

GBDT-sparse does not employ histogram approximation to speed up the training
process. It is worthwhile because it is not clear how to construct sparse
histograms, especially when combining with the second order gradient.

  * •

The main motivation of GBDT-sparse is to reduce the space complexity of
training and the size of models, whereas the main motivation of GBDT-MO is to
improve its generalization ability by capturing correlations between output
variables.

It may be hard to store the outputs in memory without sparse format when
output dimension is very large. In such a situation, GBDT-sparse is a better
choice.

##  V Experiments

We evaluate GBDT-MO on problems of multi-output regression, multi-class
classification and multi-label classification. First, we show the benefits of
GBDT-MO using two synthetic problems. Then, we evaluate GBDT-MO on six real-
world datasets. We further evaluate our sparse split finding algorithms.
Finally, we analyze the impact of diagonal approximation on the performance.
Recall that GBDT-SO is our own implementation of GBDT for single output.
Except for the split finding algorithm, all implementation details are the
same as GBDT-MO for a fair comparison. The purpose of our experiments is not
pushing state-of-the-art results on specific datasets. Instead, we would show
that GBDT-MO has better generalization ability than GBDT-SO. On synthetic
datasets, we compare GBDT-MO with GBDT-SO. On real-world datasets, we compare
GBDT-MO with GBDT-SO, XGBoost, LightGBM and GBDT-sparse222GBDT-sparse is
implemented for ourselves, and we do not report its training time due to the
slow speed.. XGBoost and LightGBM are state-of-the-art GBDT implementations
for single output, while GBDT-sparse is the most relevant work to GBDT-MO. All
experiments are conducted on a workstation with Intel Xeon CPU E5-2698 v4. We
use 4 threads on synthetic datasets, while we use 8 threads on real-world
datasets. We provide hyper-parameter settings in Appendix B. The training
process is terminated when the performance does not improve within 25 rounds.

TABLE III: RMSE on synthetic datasets | friedman1 | random projection
---|---|---
GBDT-SO | 0.1540 | 0.0204
GBDT-MO | 0.1429 | 0.0180

(a) friedman1 (b) random projection

Figure 3: RMSE curves on friedman1 and random project. RMSE curves for
training samples are drawn in solid lines, while RMSE curves for test samples
are drawn in dashed lines. Best viewed on the screen.

(a) MNIST (b) Caltech101

Figure 4: Cross entropy curves and accuracy curves on MNIST and Caltech101.
Cross entropy curves for training samples are drawn in solid lines, while
accuracy curves for test samples are drawn in dashed lines. Best viewed on the
screen.

###  V-A Synthetic Datasets

The first dataset is derived from the friedman1 regression problem [23]. Its
target $y$ is generated by:

$f(\mathbf{x})=\sin(\pi\mathbf{x}_{1}\mathbf{x}_{2})+2(\mathbf{x}_{3}-0.5)^{2}+\mathbf{x}_{4}+0.5\mathbf{x}_{5}$
(23) $y=f(\mathbf{x})+0.1\varepsilon$ (24)

where $\mathbf{x}\in\mathbb{R}^{10}$ and $\varepsilon\sim\mathcal{N}(0;1)$.
Each element of $\mathbf{x}$ is sampled from $\mathcal{U}(-1,1)$. The last
five elements of $\mathbf{x}$ is irrelevant to the target. We extend this
problem into multiple output case by adding independent noise to
$f(\mathbf{x})$. That is $\mathbf{y}_{i}=f(\mathbf{x})+0.1\varepsilon$ where
$\mathbf{y}\in\mathbb{R}^{5}$.

The second dataset is generated by random projection.

$\mathbf{y}=\mathbf{w}^{T}\mathbf{x}$ (25)

where $\mathbf{x}\in\mathbb{R}^{4}$, $\mathbf{w}\in\mathbb{R}^{4\times 8}$ and
$\mathbf{y}\in\mathbb{R}^{8}$. Each element of $\mathbf{x}$ and $\mathbf{w}$
is independently sampled from $\mathcal{U}(-1,1)$.

For both datasets, we generate $10,000$ samples for training and $10,000$
samples for test. We train them via mean square error (MSE) and evaluate the
performance on test samples via root mean square error (RMSE). We repeat the
experiments 5 times with different seeds and average the RMSE. As shown in
Table III, GBDT-MO is better than GBDT-SO. We provide the training curves in
Fig. 3. For fairness, curves of different methods are plotted with the same
learning rate and maximum tree depth. It can be observed that GBDT-MO has
better generalization ability on both datasets. This is because its RMSE on
test samples is lower and its performance gap is smaller. Here, performance
gap means the differences of performance between training samples and unseen
samples which is usually used to measure the generalization ability of machine
learning algorithms. The output variables of friedman1 are correlated because
they are observations of the same underlying variable corrupted by Gaussian
noise. The output variables of random projection are also correlated because
this projection is over-complete. This supports our claim that GBDT-MO has
better generalization abilities because its learning mechanism encourages it
to capture variable correlations. However, GBDT-MO suffers from slow
convergence speed.

TABLE IV: Dataset statistics. * means the features are pre-processed. Dataset | # Training samples | # Test samples | # Features | # Outputs | Problem type
---|---|---|---|---|---
MNIST | 50,000 | 10,000 | 784 | 10 | classification
Yeast | 1,038 | 446 | 8 | 10 | classification
Caltech101 | 6,073 | 2,604 | 324* | 101 | classification
MNIST-inpainting | 50,000 | 10,000 | 200* | 24 | regression
Student-por | 454 | 159 | 41* | 3 | regression
NUS-WIDE | 161,789 | 107,859 | 128* | 81 | multi-label
TABLE V: Performance on real-world datasets | MNIST | Yeast | Caltech101 | NUS-WIDE | MNIST-inpaining | Student-por
---|---|---|---|---|---|---
accuracy | accuracy | accuracy | top-1 accuracy | RMSE | RMSE
XGBoost | 97.86 | 62.94 | 56.52 | 43.72 | 0.26088 | 0.24623
LightGBM | 98.03 | 61.97 | 55.94 | 43.99 | 0.26090 | 0.24466
GBDT-sparse | 96.41 | 62.83 | 43.93 | 44.05 | - | -
GBDT-SO | 98.08 | 61.97 | 56.62 | 44.10 | 0.26157 | 0.24408
GBDT-MO | 98.30 | 62.29 | 57.49 | 44.21 | 0.26025 | 0.24392
Figure 5: Average training time in second for each boost round. Four relatively larger datasets are used. Y axis is plotted in log scale. TABLE VI: Training time for sparse split finding Restricted/Unrestricted | MNIST | MNIST-inpaining | Caltech101 | NUS-WIDE
---|---|---|---|---
$k=2$ | 0.619/0.599 | - | - | -
$k=4$ |  0.737/0.758 |  0.149/0.179 | - | -
$k=8$ |  0.721/0.778 |  0.161/0.195 |  0.578/0.663 |  1.248/1.266
$k=16$ | - |  0.153/0.208 |  0.620/0.738 |  1.237/1.365
$k=32$ | - | - |  0.628/0.866 |  1.256/1.453
$k=64$ | - | - |  0.650/0.937 |  1.263/1.493
TABLE VII: Performance for sparse split finding Restricted/Unrestricted | MNIST | MNIST-inpaining | Caltech101 | NUS-WIDE
---|---|---|---|---
accuracy | RMSE | accuracy | top-1 accuracy
$k=2$ | 97.54/97.79 | - | - | -
$k=4$ |  98.19/98.18 | 0.26424/0.26368 | - | -
$k=8$ |  98.07/98.06 |  0.26245/0.26391 | 53.57/54.11 | 44.20/44.20
$k=16$ | - |  0.26232/0.26248 |  55.98/54.69 | 44.22/44.27
$k=32$ | - | - |  56.71/56.64 |  44.07/44.04
$k=64$ | - | - |  58.09/57.58 |  44.23/44.22

###  V-B Real-world Datasets

In this subsection, we evaluate GBDT-MO on six real-world datasets and compare
with related methods in terms of test performance and training speed.

MNIST333yann.lecun.com/exdb/mnist/ is widely used for classification whose
samples are $28\times 28$ gray images of handwritten digits from 0 to 9. Each
sample is converted into a vector with 784 elements.

Yeast444archive.ics.uci.edu/ml/datasets/Yeast has 8 input attributions each of
which is a measurement of the protein sequence. The goal is to predict protein
localization sites with 10 possible choices.

Caltech101555www.vision.caltech.edu/Image_Datasets/Caltech101/ contains images
of objects belonging to 101 categories. Most categories have about 50 samples.
The size of each image is roughly $300\times 200$. To obtain fixed length
features, we resize each image into $64\times 64$ and compute the HOG
descriptor [24] of the resized image. Each sample is finally converted to a
vector with 324 elements.

MNIST-inpainting is a regression task based on MNIST. We crop the central
$20\times 20$ patch from the original $28\times 28$ image because most
boundary pixels are 0. Then, the cropped image is divided into upper and lower
halves each of which is a $10\times 20$ patch. The upper half is used as the
input. We further crop a $4\times 6$ small patch at the top center of the
lower half. This small patch is used as the target. The task is to predict the
pixels of this $4\times 6$ patch given upper half image.

Student-por666archive.ics.uci.edu/ml/datasets/Student+Performance predicts the
Portuguese language scores in three different grades of students based on
their demographic, social and school related features. The original scores
range from 0 to 20. We linearly transform them into $[-1,1]$. We use one-hot
coding to deal with the categorical features on this dataset.

NUS-WIDE777mulan.sourceforge.net/datasets-mlc.html is a dataset for real-world
web image retrieval [25]. Tsoumakas et.al [11] selects a subset label of this
dataset and uses it for multi-label classification. Images are represented
using 128-D cVLAD+ features described in [26].

We summarize the statistics of the above datasets in Table IV. Those datasets
are diverse in terms of scale and complexity. For MNIST, MNIST-inpainting and
NUS-WIDE, we use the official training-test split. For others, we randomly
select $70\%$ samples for training and the rest for test. We repeat this
strategy 10 times with different seeds and report the average results. For
multi-class classification, we use cross-entropy loss. For regression and
multi-label classification, we use MSE loss. For regression, the performance
is measured by RMSE. For multi-class classification and multi-label
classification, the performance is measured by top-1 accuracy.

RMSE and accuracy on real-world datasets are shown in Table V. The overall
performance of GBDT-MO is better than others. We provide statistical tests in
Appendix C to show the performance of GBDT-MO on datasets with random training
and testing split. We compare the loss curves and the accuracy curves of GBDT-
MO and GBDT-SO on MNIST and Caltech101 in Fig. 4. We conclude that GBDT-MO has
better generalization ability than GBDT-SO. Because its training loss is
higher while its test performance is better.

We also compare the training speed. Specifically, we run 10 boost rounds and
record the average training time for each boost round. We repeat this process
three times and report the average training time in seconds as shown in Fig.
5. The training speed for Yeast and Student-por is not reported here because
those two datasets are too small. GBDT-MO is remarkably faster than GBDT-SO
and XGBoost, especially when the number of outputs is large. It is slightly
faster than LightGBM. Since we eliminate the interference from outer factors,
the comparisons between GBDT-SO and GBDT-MO show the effects of the proposed
method objectively.

###  V-C Sparse Split Finding

All experiments of GBDT-MO performed in Sections V-A and V-B use non-sparse
split finding algorithm. In this subsection, We evaluate our unrestricted
sparse split finding algorithm (Algortithm 3) and restricted sparse split
finding algorithm (Algorithm 4). We compare their performance and training
speed on MNIST, MNIST-inpainting, Caltech101 and NUS-WIDE with different
sparse factor $k$. The hyper-parameters used here are the same as the
corresponding non-sparse one’s.

Their training speed in seconds is shown in Table VI. The restricted one is
faster than the unrestricted one. Their speed differences are remarkable when
$k$ is large. The performance on test samples is shown in Table VII. The
restricted one is slightly better than the unrestricted one. Note that the
performance of our sparse split algorithms is sometimes better than the non-
sparse one’s with a proper $k$ (for example, $k=64$ on Caltech101).

###  V-D Objective with Exact Hessian

GBDT-MO replaces the exact hessian with the diagonal hessian to reduce
computational complexity. We analyze the impact of the diagonal approximation
on MNIST and Yeast because their output dimensions are relatively small. Table
VIII shows accuracy and training time for a single round. Two methods have
similar accuracy but the diagonal approximation is much faster. Thus, the
diagonal approximation of hessian achieves good trade-off between accuracy and
speed.

TABLE VIII: Impact of diagonal hessian. Accuracy/Time | MNIST | Yeast
---|---|---
diagonal hessian | 98.30/0.79s | 62.29/0.0011s
exact hessian | 98.17/3.96s | 62.89/0.0094s

##  VI Conclusions

In this paper, we have proposed a general method to learn GBDT for multiple
outputs. The motivation of GBDT-MO is to capture the correlations between
output variables and reduce redundancy of tree structures. We have derived the
approximated learning objective for both non-sparse case and sparse case based
on the second order Taylor expansion of loss. For sparse case, we have
proposed the restricted algorithm which restricts the subsets of left and
right parts to be the same and the unrestricted algorithm which has no such a
restriction. We have extended the histogram approximation into multiple output
case to speed up the training process. We have evaluated that GBDT-MO is
remarkably and consistently better in generalization ability and faster in
training speed compared with GBDT for single output. We have also evaluated
that the restricted sparse split finding algorithm is slightly better than the
unrestricted one by considering both performance and training speed. However,
GBDT-MO suffers from slower convergence, especially at the beginning of
training. In our future work, we would improve its convergence speed.

## References

  * [1] Jerome H Friedman,  “Greedy function approximation: A gradient boosting machine.,”  Annals of Statistics, vol. 29, no. 5, pp. 1189–1232, 2001.
  * [2] Jerome H Friedman,  “Stochastic gradient boosting,”  Computational Statistics and Data Analysis, vol. 38, no. 4, pp. 367–378, 2002.
  * [3] Ping Li,  “Robust logitboost and adaptive base class (abc) logitboost,”  in Proceedings of the Twenty-Sixth Conference on Uncertainty in Artificial Intelligence, 2010, pp. 302–311.
  * [4] Chongchong Qi, Andy Fourie, Qiusong Chen, Xiaolin Tang, Qinli Zhang, and Rugao Gao,  “Data-driven modelling of the flocculation process on mineral processing tailings treatment,”  Journal of Cleaner Production, vol. 196, pp. 505–516, 2018.
  * [5] Christopher JC Burges,  “From ranknet to lambdarank to lambdamart: An overview,”  Learning, vol. 11, no. 23-581, pp. 81, 2010.
  * [6] Matthew Richardson, Ewa Dominowska, and Robert J Ragno,  “Predicting clicks: estimating the click-through rate for new ads,”  in Proceedings of the International Conference on World Wide Web, 2007, pp. 521–530.
  * [7] James Bennett, Stan Lanning, et al.,  “The netflix prize,”  in Proceedings of KDD Cup and Workshop. New York, NY, USA., 2007, vol. 2007, p. 35.
  * [8] Tianqi Chen and Carlos Guestrin,  “Xgboost: A scalable tree boosting system,”  Knowledge Discovery and Data Mining, pp. 785–794, 2016.
  * [9] Guolin Ke, Qi Meng, Thomas William Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tieyan Liu,  “Lightgbm: a highly efficient gradient boosting decision tree,”  in Advances in Neural Information Processing Systems, 2017, pp. 3149–3157.
  * [10] Liudmila Ostroumova Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, Anna Veronika Dorogush, and Andrey Gulin,  “Catboost: unbiased boosting with categorical features,”  in Advances in Neural Information Processing Systems, 2018, pp. 6638–6648.
  * [11] Grigorios Tsoumakas, Eleftherios Spyromitros-Xioufis, Jozef Vilcek, and Ioannis Vlahavas,  “Mulan: A java library for multi-label learning,”  Journal of Machine Learning Research, vol. 12, pp. 2411–2414, 2011\.
  * [12] Hanen Borchani, Gherardo Varando, Concha Bielza, and Pedro Larranaga,  “A survey on multi-output regression,”  Data Mining and Knowledge Discovery, vol. 5, no. 5, pp. 216–233, 2015.
  * [13] Ian Goodfellow, Yoshua Bengio, and Aaron Courville,  Deep Learning,  MIT Press, 2016,  http://www.deeplearningbook.org.
  * [14] Geoffrey E Hinton, Oriol Vinyals, and Jeffrey Dean,  “Distilling the knowledge in a neural network,”  arXiv: Machine Learning, 2015.
  * [15] Pierre Geurts, Louis Wehenkel, and Florence Dalchebuc,  “Gradient boosting for kernelized output spaces,”  in Proceedings of the International Conferenece on Machine Learning, 2007, pp. 289–296.
  * [16] Si Si, Huan Zhang, S Sathiya Keerthi, Dhruv Mahajan, Inderjit S Dhillon, and Chojui Hsieh,  “Gradient boosted decision trees for high dimensional sparse output,”  in Proceedings of the International Conferenece on Machine Learning, 2017, pp. 3182–3190.
  * [17] Stephen Tyree, Kilian Q Weinberger, Kunal Agrawal, and Jennifer Paykin,  “Parallel boosted regression trees for web search ranking,”  in Proceedings of the International Conference on World Wide Web, 2011, pp. 387–396.
  * [18] Tianqi Chen, Sameer Singh, Ben Taskar, and Carlos Guestrin,  “Efficient second-order gradient boosting for conditional random fields,”  in Proceedings of the International Conference on Artificial Intelligence and Statistics, 2015, pp. 147–155.
  * [19] Rahul Agrawal, Archit Gupta, Yashoteja Prabhu, and Manik Varma,  “Multi-label learning with millions of labels: recommending advertiser bid phrases for web pages,”  in Proceedings of the International Conference on World Wide Web, 2013, pp. 13–24.
  * [20] Yashoteja Prabhu and Manik Varma,  “Fastxml: a fast, accurate and stable tree-classifier for extreme multi-label learning,”  in Proceedings of the International Conference on Knowledge Discovery and Data Mining, 2014, pp. 263–272.
  * [21] Robert E Schapire and Yoram Singer,  “Improved boosting algorithms using confidence-rated predictions,”  in Proceedings of the International Conference on Learning Theory, 1998, vol. 37, pp. 80–91.
  * [22] Yonatan Amit, Ofer Dekel, and Yoram Singer,  “A boosting algorithm for label covering in multilabel problems,”  in Proceedings of the International Conference on Artificial Intelligence and Statistics, 2007, pp. 27–34.
  * [23] Jerome H Friedman,  “Multivariate adaptive regression splines,”  Annals of Statistics, vol. 19, no. 1, pp. 1–67, 1991.
  * [24] Navneet Dalal and Bill Triggs,  “Histograms of oriented gradients for human detection,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2005, vol. 1, pp. 886–893.
  * [25] Tat-Seng Chua, Jinhui Tang, Richang Hong, Haojie Li, Zhiping Luo, and Yan-Tao Zheng,  “Nus-wide: A real-world web image database from national university of singapore,”  in Proceedings of the ACM Conference on Image and Video Retrieval, Santorini, Greece, July 8-10, 2009.
  * [26] Eleftherios Spyromitrosxioufis, Symeon Papadopoulos, Ioannis Kompatsiaris, Grigorios Tsoumakas, and Ioannis P Vlahavas,  “A comprehensive study over vlad and product quantization in large-scale image retrieval,”  IEEE Transactions on Multimedia, vol. 16, no. 6, pp. 1713–1728, 2014.

TABLE IX: Grid searched learning rate and maximum depth Learning rate/Maximum depth | MNIST | Yeast | Caltech101 | NUS-WIDE | MNIST-inpaining | Student-por
---|---|---|---|---|---|---
XGBoost | 0.10/8 | 0.10/4 | 0.10/8 | 0.10/8 | 0.10/8 | 0.10/4
LightGBM | 0.25/6 | 0.10/4 | 0.10/9 | 0.05/8 | 0.10/8 | 0.05/4
GBDT-sparse | 0.05/8 | 0.10/5 | 0.10/9 | 0.10/8 | - | -
GBDT-SO | 0.10/6 | 0.10/4 | 0.05/8 | 0.05/8 | 0.10/8 | 0.05/4
GBDT-MO | 0.10/8 | 0.10/5 | 0.10/10 | 0.10/8 | 0.10/7 | 0.10/4
TABLE X: Number of trees of learned models | MNIST | Yeast | Caltech101 | NUS-WIDE | MNIST-inpaining | Student-por
---|---|---|---|---|---|---
XGBoost | 1700 | 299.0 | 2234.0 | 5921 | 4880 | 133.4
LightGBM | 1760 | 334.0 | 2692.0 | 11636 | 5618 | 219.6
GBDT-sparse | 645 | 21.4 | 237.4 | 494 | - | -
GBDT-SO | 2180 | 183.0 | 3697.0 | 12051 | 6373 | 373.7
GBDT-MO | 615 | 85.1 | 517.7 | 3546 | 2576 | 125.0
TABLE XI: Preliminarily searched hyper-parameters | MNIST | Yeast | Caltech101 | NUS-WIDE | MNIST-inpaining | Student-por
---|---|---|---|---|---|---
Gain threshold | 1e-3 | 1e-3 | 1e-3 | 1e-6 | 1e-6 | 1e-6
Min samples in leaf | 16 | 16 | 16 | 4 | 4 | 4
Max bins | 8 | 32 | 32 | 64 | 16 | 8

##  Appendix A Approximated Objective

We suppose there exists $\gamma>0$ such that

$\displaystyle\gamma\mathbf{H}_{ii}\geq\sum_{j}|\mathbf{H}_{ij}|,\quad\forall
i$ (26)

That is, $\mathbf{H}$ is dominated by its diagonal elements. Then, we have:

$\displaystyle\mathbf{w}^{T}\mathbf{Hw}$
$\displaystyle=\sum_{i}\sum_{j}\mathbf{H}_{ij}\mathbf{w}_{i}\mathbf{w}_{j}$
(27)
$\displaystyle\leq\frac{1}{2}\sum_{i}\sum_{j}|\mathbf{H}_{ij}|(\mathbf{w}_{i}^{2}+\mathbf{w}_{j}^{2})$
$\displaystyle=\sum_{i}\sum_{j}|\mathbf{H}_{ij}|\mathbf{w}_{i}^{2}$
$\displaystyle\leq\gamma\sum_{i}\mathbf{H}_{ii}\mathbf{w}_{i}^{2}$

The last inequality holds based on the assumption in (26). Substituting this
result into the Taylor expansion of loss

$\displaystyle l(\mathbf{\hat{y}}+\mathbf{w},\mathbf{y})$
$\displaystyle\approx
l(\mathbf{\hat{y}},\mathbf{y})+\mathbf{g}^{T}\mathbf{w}+\frac{1}{2}\mathbf{w}^{T}\mathbf{Hw}$
(28) $\displaystyle\leq
l(\mathbf{\hat{y}},\mathbf{y})+\sum_{i}\mathbf{g}_{i}\mathbf{w}_{i}+\frac{\gamma}{2}\sum_{i}\mathbf{H}_{ii}\mathbf{w}_{i}^{2}$

We ignore the remainder term of the Taylor expansion [8]. Substituting this
upper bound of $l(\mathbf{\hat{y}}+\mathbf{w},\mathbf{y})$ into (12), we get
the optimal value of $\mathbf{w}$.

$\mathbf{w}^{*}_{j}=-\frac{\sum_{i}(\mathbf{g}_{j})_{i}}{\gamma\sum_{i}(\mathbf{h}_{j})_{i}+\lambda}$
(29)

Recall that $\mathbf{h}$ is the diagonal elements of $\mathbf{H}$. This
solution is the same as (15), except the coefficient $\gamma$. In practice,
the actual leaf weight is $\alpha\mathbf{w}^{*}_{j}$ where $\alpha$ is the so-
called learning rate. Then, the effect of $\gamma$ can be canceled out by
adjusting $\alpha$ and $\lambda$. Thus, we do not need to consider the exact
value of $\gamma$ in practice.

##  Appendix B Hyper-parameter Settings

We discuss our hyper-parameter settings. We find the best maximum depth $d$
and learning rate via a grid search. $d$ is searched from
$\\{4,5,6,7,8,9,10\\}$ and learning rate is searched from
$\\{0.05,0.1,0.25,0.5\\}$. We provide the selected $d$ and learning rate in
Table IX. We set $L_{2}$ regularization $\lambda$ to 1.0 and the maximum
leaves to $0.75\times 2^{d}$ in all experiments. We preliminarily search other
hyper-parameters and fix them for all methods. We list them on real-world
datasets in Table XI. Note that the hyper-parameters in Table XI may not be
optimal. When comparing with the gain threshold, we use the average gain over
the involved output variables. We provide the number of trees of the best
models in Table X. For single variable GBDT models, the number of trees is the
summation of the number of trees for each output variable.

##  Appendix C Statistical Test of Results

There are no standard training and testing splits on Yeast, Caltech101 and
Student-por. We obtain the results in Table V by averaging the results of 10
random trials. Thus, statistical tests on how likely GBDT-MO is better than
others are necessary. Denote $X$ is a random variable which means the
performance differences between GBDT-MO and the others. $X$ is commonly
supposed to be Student’s T-distribution because its true variance is unknown.
We estimate the mean and standard deviation of $X$ from the 10 random trials.
Then, the confidence in whether GBDT-MO performs better than the others is
defined as $P(X>0)$ for classifications and $P(X<0)$ for regressions, as
reported in Table XII. Most confidence scores are significantly higher than
0.5. Thus, it can be expected that GBDT-MO performs better than others.

TABLE XII: Confidence of the superiority of GBDT-MO Confidence | Yeast | Caltech101 | Student-por
---|---|---|---
XGBoost | 0.067 | 0.930 | 0.888
LightGBM | 0.745 | 0.988 | 0.685
GBDT-SO | 0.812 | 0.955 | 0.561