# Demonstrating Quantum Microscopic Reversibility Using Coherent States of
Light

Marco Bellini<EMAIL_ADDRESS>Istituto Nazionale di Ottica (CNR-
INO), Largo E. Fermi 6, 50125 Florence, Italy  LENS and Department of Physics
and Astronomy, University of Firenze, 50019 Sesto Fiorentino, Florence, Italy
Hyukjoon Kwon<EMAIL_ADDRESS>Korea Institute for Advanced Study, Seoul
02455, South Korea  Nicola Biagi  Istituto Nazionale di Ottica (CNR-INO),
Largo E. Fermi 6, 50125 Florence, Italy  LENS and Department of Physics and
Astronomy, University of Firenze, 50019 Sesto Fiorentino, Florence, Italy
Saverio Francesconi  Istituto Nazionale di Ottica (CNR-INO), Largo E. Fermi 6,
50125 Florence, Italy  LENS and Department of Physics and Astronomy,
University of Firenze, 50019 Sesto Fiorentino, Florence, Italy  Alessandro
Zavatta  Istituto Nazionale di Ottica (CNR-INO), Largo E. Fermi 6, 50125
Florence, Italy  LENS and Department of Physics and Astronomy, University of
Firenze, 50019 Sesto Fiorentino, Florence, Italy  M. S. Kim
<EMAIL_ADDRESS>QOLS, Blackett Laboratory, Imperial College London,
London SW7 2AZ, United Kingdom

###### Abstract

The principle of microscopic reversibility lies at the core of fluctuation
theorems, which have extended our understanding of the second law of
thermodynamics to the statistical level. In the quantum regime, however, this
elementary principle should be amended as the system energy cannot be sharply
determined at a given quantum phase space point. In this Letter, we propose
and experimentally test a quantum generalization of the microscopic
reversibility when a quantum system interacts with a heat bath through energy-
preserving unitary dynamics. Quantum effects can be identified by noting that
the backward process is less likely to happen in the existence of quantum
coherence between the system’s energy eigenstates. The experimental
demonstration has been realized by mixing coherent and thermal states in a
beam-splitter, followed by heterodyne detection in an optical setup. We verify
that the quantum modification for the principle of microscopic reversibility
is critical in the low-temperature limit, while the quantum-to-classical
transition is observed as the temperature of the thermal field gets higher.

Fluctuation theorems Gallavotti95 ; Jarzynski97 ; Crooks99  set the
cornerstone for describing the statistical properties of non-equilibrium
thermodynamics. Throughout the development of fluctuation theorems, the
principle of microscopic reversibility has served as the main pillar, by
providing a universal symmetry relation between a phase space trajectory and
its time-reversed version. The second law of thermodynamics naturally arises
from averaging fluctuating quantities over all possible trajectories, which
leads to macroscopic irreversibility. Fluctuation theorems have been
experimentally verified for various types of systems Liphardt02 ; Collin05 ;
Blickle06 ; Saira12  in the classical regime, where the thermodynamic
elements, such as work, energy, and heat, are determined from the phase space
trajectories by well-established classical definitions.

On the other hand, there exists a fundamental difference between quantum and
classical phase space structures, as the uncertainty principle prohibits the
position and momentum of a quantum system from being measured precisely at the
same time. Furthermore, energy is not necessarily a unique basis to
investigate the thermodynamic quantities of a quantum state and its evolution.
A quantum state can be prepared in a superposition of energy eigenstates, and
measurement can be performed for an observable that does not commute with the
system Hamiltonian. Considerable efforts from both thermodynamics Esposito09 ;
Campisi11 ; Deffner11 ; Deffner13 ; Hanggi15 ; Jarzynski15 ; Funo18 ;
Taranto18  and quantum information Goold16 ; Alhambra16 ; Binder19 ; Landi21
perspectives have been put into an investigation of how quantum coherence
changes the fluctuation theorems, as well as how we generalize thermodynamic
quantities, such as work and heat, consistent with quantum mechanics Talkner07
; Bera19 . Throughout these efforts, quantum optics has played an important
role by providing a promising platform to simulate quantum processes in the
microscopic regime Scully03 ; Pikovski12 ; Vidrighin16 .

Recent developments of quantum fluctuation theorems Albash13 ; Batalhao14 ;
Aberg18 ; Manzano18 ; Holmes19 ; Kwon19 ; Micadei20 ; Micadei21 ; Goold21 ;
Yada22  have proposed that the symmetry relation between the forward and
backward quantum thermodynamic processes should be modified properly to fully
describe the quantum trajectories including the contributions from coherence.
In this Letter, we establish a method to explore the difference between
quantum and classical theory of microscopic reversibility on a very basic
level and demonstrate it using a quantum optical setup. In the experiment, an
initial coherent state of light interacts with a thermal field, and the
resulting mixed state is then measured on a coherent basis. The measured
statistics clearly follow the quantum version of microscopic reversibility,
and we identify that quantum coherence plays a crucial role in the low-
temperature limit. The experimental data agree with the theoretically
predicted behaviors of the quantum modification with respect to the degree of
coherence and the temperature of the thermal field.

Microscopic reversibility in classical dynamics.— Suppose that a physical
system is interacting with a thermal bath at temperature $T$. In classical
statistical mechanics, a phase space trajectory of the dynamics is
characterized by the system’s position and momentum at each time $t$ as
$z(t)=(x(t),p(t))$. Consequently, the time-reversed trajectory is defined as
$\bar{z}(-t)=(\bar{x}(-t)),\bar{p}(-t))=(x(-t),-p(-t))$ by the transformation
$\bar{x}=x$ and $\bar{p}=-p$. The principle of microscopic reversibility tells
us that there exists a symmetry relation between the probabilities of the
trajectory $z(t)$ following the time-forward dynamics and the trajectory
$\bar{z}(-t)$ following the time-reversed dynamics. When considering only the
initial and final phase space points $z_{i}=(x_{i},p_{i})$ and
$z_{f}=(x_{f},p_{f})$, the symmetry relation between the forward transition
probability ${\cal P}_{\rightarrow}\left(z_{f}|z_{i}\right)$ from $z_{i}$ to
$z_{f}$ and the backward transition probability ${\cal
P}_{\leftarrow}\left(\bar{z}_{i}|\bar{z}_{f}\right)$ from $\bar{z}_{f}$ to
$\bar{z}_{i}$ can be simplified as

$\frac{{\cal P}_{\rightarrow}\left(z_{f}|z_{i}\right)}{{\cal
P}_{\leftarrow}\left(\bar{z}_{i}|\bar{z}_{f}\right)}=e^{-\beta Q},$ (1)

where $Q$ is the heat flow from the bath to the system during the forward
process and $\beta=1/(k_{B}T)$ with the Boltzmann constant $k_{B}$. If the
system-bath is isolated and no work has been done to the system, the amount of
heat is equal to the system’s energy change $Q=\Delta
E=H_{f}(z_{f})-H_{i}(z_{i})$, where $H_{i}(z_{i})$ and $H_{f}(z_{f})$ denote
the system’s Hamiltonians at initial and final points. From the microscopic
reversibility condition in Eq. (1), various types of classical fluctuation
theorems, including the Jarzynski equality Jarzynski97  and Crooks fluctuation
theorem Crooks99  can be derived.

Figure 1: (a) Forward and (b) backward quantum transition probabilities when a
quantum state is interacting with the bath $\hat{\rho}_{\rm th}$ and (c) the
rescaling of time-reversed coherent states.

Microscopic reversibility in quantum dynamics.— In quantum mechanics, however,
phase space should be defined differently as position and momentum are
promoted to operators $\hat{x}$ and $\hat{p}$, which do not commute to each
other, i.e., $[\hat{x},\hat{p}]=i\hbar$. By taking into account the non-
classical nature of quantum mechanics, we can construct a phase space
structure for a quantum harmonic oscillator system with the Hamiltonian
$\hat{H}(\hat{x},\hat{p})=\frac{\hbar\omega_{0}}{2}\left(\hat{x}^{2}+\hat{p}^{2}\right)=\hbar\omega_{0}\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right)$,
where $\hat{a}=(\hat{x}+i\hat{p})/\sqrt{2}$ is a bosonic annihilation
operator. In particular, coherent states, satisfying
$\hat{a}|\alpha\rangle=\alpha|\alpha\rangle$, compose a family of the most
classical states among pure quantum states Mandel86 ; Zurek93 . The coherent
state saturates the minimum uncertainty relation, and its distribution is
centered at $z=(x,p)=({\rm Re}[\alpha],{\rm Im}[\alpha])/\sqrt{2}$ in quantum
phase space. While the coherent state resembles a classical state at a certain
phase space point, an important difference arises from the fact that the state
contains quantum coherence between different energy states $|n\rangle$;
$|\alpha\rangle=e^{-|\alpha|^{2}/2}\sum_{n=0}^{\infty}\left(\alpha^{n}/\sqrt{n!}\right)|n\rangle$.

Now let us generalize the classical scenario to the quantum harmonic
oscillator system by interacting the system with a bath in thermal
equilibrium, $\hat{\rho}_{\rm
th}=e^{-\beta\hat{H}_{B}}/\mbox{Tr}[e^{-\beta\hat{H}_{B}}]$, where
$\hat{H}_{B}$ is a bath Hamiltonian. By denoting the unitary operator
generated by an interaction Hamiltonian $\hat{H}_{I}$ as
$\hat{U}=e^{-i\tau\hat{H}_{I}}$, the transition probability from the initial
state $|\psi_{i}\rangle$ to the final state $|\psi_{f}\rangle$ is defined as

${\cal
P}_{\rightarrow}(\psi_{f}|\psi_{i})={\mbox{Tr}}\left[\hat{U}(|\psi_{i}\rangle\langle\psi_{i}|\otimes\hat{\rho}_{\rm
th})\hat{U}^{\dagger}\left(|\psi_{f}\rangle\langle\psi_{f}|\otimes\mathbb{\hat{1}}_{B}\right)\right].$
(2)

The quantum mechanical extension of the microscopic reversibility condition
can be established by finding the backward process that corresponds to a
generalized time-reversal Aberg18 ; Kwon19 . In particular, when the system’s
equilibrium state $\hat{\rho}_{\rm eq}=e^{-\beta\hat{H}}/{\rm
Tr}[e^{-\beta\hat{H}}]$ remains unchanged by the system-bath interaction, the
following relation holds Kwon19  between the forward and backward
probabilities:

$\frac{{\cal P}_{\rightarrow}(\psi_{f}|\psi_{i})}{{\cal
P}_{\leftarrow}(\tilde{\psi}_{i}|\tilde{\psi}_{f})}=\langle\psi_{i}|e^{\beta\hat{H}}|\psi_{i}\rangle\langle\psi_{f}|e^{-\beta\hat{H}}|\psi_{f}\rangle,$
(3)

where $|\tilde{\psi}_{i}\rangle\propto
e^{\beta\hat{H}/2}|\bar{\psi}_{i}\rangle$ and $|\tilde{\psi}_{f}\rangle\propto
e^{-\beta\hat{H}/2}|\bar{\psi}_{f}\rangle$ with the normalization
$\langle\tilde{\psi}_{i(f)}|\tilde{\psi}_{i(f)}\rangle=1$. Here,
$|\bar{\psi}\rangle$ denotes the time-reversed state footnote  of
$|\psi\rangle$.

However, the explicit construction of the backward process in Eq. (3) can be
very complicated in general. To consider the simplest scenario, we assume that
the system and bath Hamiltonians are time-independent and invariant under
time-reversal, i.e., $\hat{H}(\hat{x},\hat{p})=\hat{H}(\hat{x},-\hat{p})$ and
$\hat{H}_{B}(\hat{x}_{B},\hat{p}_{B})=\hat{H}_{B}(\hat{x}_{B},-\hat{p}_{B})$,
and the interaction Hamiltonian satisfies the energy conservation condition,

$[\hat{H}_{I},\hat{H}+\hat{H}_{B}]=0.$ (4)

We note that the equilibrium state $\hat{\rho}_{\rm eq}$ becomes a fixed point
under the dynamics as $\hat{U}(\hat{\rho}_{\rm eq}\otimes\hat{\rho}_{\rm
th})\hat{U}^{\dagger}=\hat{\rho}_{\rm eq}\otimes\hat{\rho}_{\rm th}$. Under
these assumptions, the backward dynamics is simply described by the time-
reversed unitary operator $\hat{U}_{R}=e^{-i\tau\hat{H}_{I}^{R}}$, where
$\hat{H}_{I}^{R}(\hat{x},\hat{p},\hat{x}_{B},\hat{p}_{B})=\hat{H}_{I}(\hat{x},-\hat{p},\hat{x}_{B},-\hat{p}_{B})$
Aberg18 . Consequently, the backward probability of the quantum trajectory
becomes

${\cal
P}_{\leftarrow}(\bar{\psi}_{i}|\bar{\psi}_{f})={\mbox{Tr}}\left[\hat{U}_{R}(|\bar{\psi}_{f}\rangle\langle\bar{\psi}_{f}|\otimes\hat{\rho}_{\rm
th})\hat{U}_{R}^{\dagger}\left(|\bar{\psi}_{i}\rangle\langle\bar{\psi}_{i}|\otimes\mathbb{\hat{1}}_{B}\right)\right]$
(5)

(see Fig. 1(a) and (b)). For such a process, often referred to as a thermal
operation Janzing00 ; Morodecki13 ; Lostaglio18 , the heat can be solely
characterized by the system’s energy change $Q=\Delta E$ so that no work is
done on the system, just as in the previously discussed classical scenario.

One may then ask whether this relation reduces to the classical microscopic
reversibility condition in Eq. (1) by taking
$|\psi_{i(f)}\rangle=|\alpha_{i(f)}\rangle$ for the transitions between two
coherent states. We note that this is not the case because coherent states are
not the eigenstates of $\hat{H}$. Instead, from Eq. (3), we find the
microscopic reversibility condition between the coherent states,

$\frac{{\cal P}_{\rightarrow}(\alpha_{f}|\alpha_{i})}{{\cal
P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})}=\exp\left[-\frac{|\alpha_{f}|^{2}}{n_{\rm
th}+1}+\frac{|\alpha_{i}|^{2}}{n_{\rm th}}\right],$ (6)

where $\tilde{\alpha}_{i}=\bar{\alpha}_{i}e^{\beta\hbar\omega_{0}/2}$ and
$\tilde{\alpha}_{f}=\bar{\alpha}_{f}e^{-\beta\hbar\omega_{0}/2}$ with
$\bar{\alpha}$ being the complex conjugate of $\alpha$, and $n_{\rm
th}=e^{-\beta\hbar\omega_{0}}/(1-e^{-\beta\hbar\omega_{0}})$ is the mean
excitation number when the system equilibrates to the bath. We note that the
ratio of the forward-to-backward transition probabilities is solely determined
by the initial and final coherent state amplitudes.

We observe that two kinds of modifications are made for the microscopic
reversibility condition in the quantum regime. The first comes from the
rescaling of the reverse trajectory with respect to the bath temperature such
that $|\tilde{\alpha}_{i}|\geq|\alpha_{i}|$ and
$|\tilde{\alpha}_{f}|\leq|\alpha_{f}|$ (see Fig. 1(c)). This is because a
quantum state can contain coherence between energy eigenstates, where each of
them should be weighted differently depending on the energy and the bath
temperature Aberg18 ; Kwon19 . While the rescaled states provide a quadratic
correction in $\beta\hbar\omega_{0}$ by noting that $\langle\alpha|\alpha
e^{\pm\beta\hbar\omega_{0}/2}\rangle\approx
1-(|\alpha|^{2}/8)(\beta\hbar\omega_{0})^{2}$, this so-called Gibbs rescaling
is essential to establish the microscopic reversibility, such that the ratio
between the forward and backward trajectories does not depend on the detailed
form of the interaction unitary $\hat{U}$.

The second quantum modification can be captured by the following quantity,

$\Upsilon=e^{\beta Q}\left[\frac{{\cal
P}_{\rightarrow}(\psi_{f}|\psi_{i})}{{\cal
P}_{\leftarrow}(\tilde{\psi}_{i}|\tilde{\psi}_{f})}\right]=\frac{\langle\psi_{i}|e^{\beta\hat{H}}|\psi_{i}\rangle}{e^{\beta\langle\psi_{i}|\hat{H}|\psi_{i}\rangle}}\frac{\langle\psi_{f}|e^{-\beta\hat{H}}|\psi_{f}\rangle}{e^{-\beta\langle\psi_{f}|\hat{H}|\psi_{f}\rangle}},$
(7)

where $Q=\Delta
E=\langle\psi_{f}|\hat{H}|\psi_{f}\rangle-\langle\psi_{i}|\hat{H}|\psi_{i}\rangle$.
The quantity $\Upsilon$ indicates how far the right-hand side of Eq. (6)
deviates from $e^{-\beta Q}$. The last expression in Eq. (7) shows that the
correction comes from the fact that
$\langle\psi|e^{\pm\beta\hat{H}}|\psi\rangle$ is always larger than
$e^{\pm\beta\langle\psi|\hat{H}|\psi\rangle}$ when a quantum state
$|\psi\rangle$ contains superposition between energy eigenstates. Thus, the
quantum effect always gives a correction in the way that the backward-to-
forward process ratio is smaller than the classical prediction, i.e.,
$\Upsilon\geq 1$, which means that the Gibbs-rescaled backward process is less
likely to happen if quantum coherence is involved. We find the explicit form
of the quantum modification for coherent states,

$\displaystyle\log\Upsilon$ $\displaystyle=A|\alpha|^{2}_{\rm
tot}-B\left(\Delta|\alpha|^{2}\right),$ (8)

where $|\alpha|^{2}_{\rm tot}=|\alpha_{i}|^{2}+|\alpha_{f}|^{2}$ and
$\Delta|\alpha|^{2}=|\alpha_{f}|^{2}-|\alpha_{i}|^{2}$ with coefficients
$A=\cosh(\beta\hbar\omega_{0})-1$ and
$B=\sinh(\beta\hbar\omega_{0})-\beta\hbar\omega_{0}$. The leading order of the
quantum modification in $\Upsilon$ scales quadratically with respect to
$\beta\hbar\omega_{0}$ and linearly with respect to $|\alpha|^{2}_{\rm tot}$.
The larger the $|\alpha|^{2}$, the more different energy terms are involved in
$|\alpha\rangle$ so the larger the coherence. More precisely, contribution of
different energy terms in a coherent state $|\alpha\rangle$ can be quantified
by the quantum Fisher information Kwon18 , which is proportional to the
variance of energy, ${\rm
Var}_{|\alpha\rangle}(\hat{H})=\langle\alpha|\hat{H}^{2}|\alpha\rangle-\langle\alpha|\hat{H}|\alpha\rangle^{2}=(\hbar\omega_{0})^{2}|\alpha|^{2}$.
For a general transition from the state $|\psi_{i}\rangle$ to
$|\psi_{f}\rangle$, the lowest order correction is given as
$\log\Upsilon=\frac{\beta^{2}}{2}\left[{\rm
Var}_{|\psi_{i}\rangle}(\hat{H})+{\rm
Var}_{|\psi_{f}\rangle}(\hat{H})\right]+{\cal O}((\beta\hbar\omega_{0})^{3})$.
This implies that the quantum correction becomes even larger for the
transition between non-classical states with high energy variances, such as
squeezed states.

Experimental test of microscopic reversibility.—

Figure 2: Simplified experimental setup for testing the quantum fluctuation
theorem. Cubes represent 50$\%$ beam-splitters, and HT-BS are high-
transmittivity beam-splitters. All other symbols are defined in the text.
Figure 3: Experimental data for the quantum microscopic reversibility
described in Eq. (6). We plot the logarithmic ratio of forward and backward
transition probabilities using different initial state amplitudes
($|\alpha_{i}|$ in the range from 1.46 to 3.36) and different mean photon
number thermal states ($n_{\rm th}$ in the range from 1.22 to 3.57). The
reflectivity of the beam-splitter combining the coherent and thermal states is
varied as well (85$\%$ for data with $n_{\rm th}=1.62$, 70$\%$ for the
others).

We experimentally demonstrate the coherence effects on the microscopic
reversibility using quantum optics. We take an optical mode $\hat{a}$ as a
system and another optical mode $\hat{b}$ as a bath state, whose Hamiltonians
are expressed as
$\hat{H}=\hbar\omega_{0}\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right)$ and
$\hat{H}_{B}=\hbar\omega_{0}\left(\hat{b}^{\dagger}\hat{b}+\frac{1}{2}\right)=\frac{\hbar\omega_{0}}{2}(\hat{x}_{B}^{2}+\hat{p}_{B}^{2})$,
respectively. The energy conservation condition in Eq. (4) of the interaction
unitary can be achieved by interacting the system and bath through passive
linear optical elements, which preserves the total photon number
$\hat{N}=\hat{a}^{\dagger}\hat{a}+\hat{b}^{\dagger}\hat{b}$. In the
experiment, we adopt the simplest model of linear system-bath interaction
given by a beam-splitter
$\hat{U}(\theta)=e^{-i\theta(\hat{a}\hat{b}^{\dagger}+\hat{a}^{\dagger}\hat{b})}=e^{-i\theta(\hat{x}\hat{x}_{B}+\hat{p}\hat{p}_{B})}$
with real values of $\theta$. By noting that the interaction Hamiltonian
$\hat{H}_{I}(\hat{x},\hat{p},\hat{x}_{B},\hat{p}_{B})=\hat{x}\hat{x}_{B}+\hat{p}\hat{p}_{B}$
is invariant under the time-reversal
$(\hat{x},\hat{p},\hat{x}_{B},\hat{p}_{B})\rightarrow(\hat{x},-\hat{p},\hat{x}_{B},-\hat{p}_{B})$,
the backward unitary becomes the same as the forward unitary, i.e.,
$U_{R}(\theta)=U(\theta)$, which simplifies the experimental procedure.

The transition probability of the forward (backward) trajectory is obtained by
preparing the coherent state $|\alpha_{i}\rangle$
$(|\tilde{\alpha}_{f}\rangle)$ and interacting with the bath using the beam-
splitter, followed by heterodyne measurements of the system state. The
distribution of the heterodyne detection outcomes becomes the Husimi
$Q$-function Richter98 , and the probability density at the point
$|\alpha_{f}\rangle$ $(|\tilde{\alpha}_{i}\rangle)$ indicates the transition
probability of the forward (backward) trajectory by noting that the
$Q$-function of a quantum state $\hat{\rho}$ is given by
$Q_{\hat{\rho}}(\alpha):=\frac{1}{\pi}\langle\alpha|\hat{\rho}|\alpha\rangle$.

Figure 2 shows a simplified scheme of the experiment to test the microscopic
reversibility condition in Eq. (6). A pulsed Ti:sapphire laser is used as the
main source to produce coherent states $|\alpha\rangle$ of adjustable
amplitude, and thermal states $\hat{\rho}_{\rm th}$ of different mean photon
number $n_{\rm th}$ by phase/amplitude randomization via a rotating glass disk
(RD) Arecchi1965 ; Zavatta2007 ; Parigi2009 . After mixing them in a variable-
ratio beam-splitter (BS), the output state undergoes heterodyne detection,
using another portion of the laser emission as the local oscillator (LO), for
the simultaneous measurement of its orthogonal field quadratures $\hat{x}$ and
$\hat{p}$ and the reconstruction of its $Q$-function (more details on the
experimental setup are available in the Supplemental Material Suppl ). The
measurement procedure starts by fixing a mean photon number $n_{\rm th}$ of
the thermal state, which is mixed with an initial coherent state
$\lvert\alpha_{i}\rangle$ in the beam-splitter BS set at the appropriate value
of reflectivity. Then, one measures the $Q$-function of the output state and,
by evaluating it in $\alpha_{f}$, retrieves the transition probability of the
forward process, ${\cal P}_{\rightarrow}(\alpha_{f}|\alpha_{i})$. For the
time-reversed process, we inject a coherent state with a Gibbs rescaled
amplitude $\tilde{\alpha}_{f}=\bar{\alpha}_{f}e^{-\beta\hbar\omega_{0}/2}$
into the BS, set at the same reflectivity value and with the same thermal
state at the other input. The $Q$-function of the output state is then finally
evaluated in the Gibbs rescaled amplitude
$\tilde{\alpha}_{i}=\bar{\alpha}_{i}e^{\beta\hbar\omega_{0}/2}$ to obtain
${\cal P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})$.

Experimental results presented in Fig. 3 clearly demonstrate that the ratio
between the forward and backward transition probabilities well matches the
quantum prediction of the microscopic reversibility in Eq. (6) for a wide
range of values of coherent state amplitudes and bath temperatures. Note that
we used both real and complex values for the coherent state amplitudes
$\alpha_{i}$. In the first case, the rescaling for the reverse process simply
implies a change in the absolute value of the coherent state amplitudes,
whereas in the second case it also involves a phase reversal due to complex
conjugation. We also successfully verify that the microscopic reversibility
condition holds regardless of the values of the beam-splitter reflectivity.

Figure 4: Quantum modification of the microscopic reversibility given by
$\log\Upsilon$ with respect to the overall coherence $|\alpha|^{2}_{\rm tot}$
and for different values of the mean thermal photon number $n_{\rm th}$. The
experimental data verify that the modification $\log\Upsilon$ scales linearly
by increasing the coherence as predicted in Eq. (8).

Figure 4 shows experimental data and theoretical predictions for the linear
scaling of $\log\Upsilon$ with respect to the overall coherence
$|\alpha|^{2}_{\rm tot}$ for different mean numbers of thermal photons.

Figure 5:  Temperature dependence of the quantum correction. The experimental
data show the quadratic scaling of $\log\Upsilon/|\alpha|^{2}_{\rm tot}$ with
the inverse temperature $\beta\hbar\omega_{0}$, regardless of the initial and
final coherent state amplitudes. Quantum to classical transition is observed
in the high temperature limit $\beta\hbar\omega_{0}\ll 1$.

The significant difference between the classical and quantum descriptions is
investigated by studying the modification factor $\Upsilon$ in Eq. (7) as a
function of $|\alpha|^{2}_{\rm tot}$. The lowest order approximation of the
quantum correction in Eq. (8) is given by
$\log\Upsilon\approx\frac{(\beta\hbar\omega_{0})^{2}}{2}|\alpha|^{2}_{\rm
tot}$, while the second order bound is obtained as
$\left|\log\Upsilon-\frac{(\beta\hbar\omega_{0})^{2}}{2}|\alpha|^{2}_{\rm
tot}\right|\leq\frac{(\beta\hbar\omega_{0})^{3}}{6}|\alpha|^{2}_{\rm
tot}+{\cal O}((\beta\hbar\omega_{0})^{4})$ by noting that
$\left|\Delta|\alpha|^{2}\right|\leq|\alpha|^{2}_{\rm tot}$.

In the low-temperature limit, the quantum modification in the experiment is
observed to reach $\Upsilon>200$ for the transition from $\alpha_{i}=2.4$ to
$\alpha_{f}=5.2$ with $n_{\rm th}=1.26$ ($\beta\hbar\omega_{0}\approx 0.58$).
This can be understood from the fact that the Gibbs rescaling for the reverse
process severely deforms quantum states having coherences between the energy
eigenstates when $\beta\hbar\omega_{0}\sim 1$. Instead, the quantum-to-
classical transition is seen in the high-temperature limit
($\beta\hbar\omega_{0}\ll 1$), where the quantum fluctuation relation
resembles the classical fluctuation relation in Eq. (1). For example, by
taking $n_{\rm th}=3.57$ (corresponding to $\beta\hbar\omega_{0}\approx
0.25$), the modification measured in the experiment is $\Upsilon\approx 1.57$
for $\alpha_{i}=3.14$ and $\alpha_{f}=2.17$. Figure 5 clearly shows the
quadratic scaling of the quantum correction with respect to the inverse
temperature $\beta$, where the leading coefficient after normalization is
given by $\frac{\log\Upsilon}{|\alpha|^{2}_{\rm
tot}}\approx\frac{(\beta\hbar\omega_{0})^{2}}{2}$.

Remarks.— We demonstrate the principle of microscopic reversibility in the
quantum regime via quantum optics experiments. Our method provides a general
extension of the classical thermodynamic interactions to a quantum scenario,
by replacing classical phase space points with coherent states in quantum
optics and by modeling the system-bath interaction using an energy-preserving
unitary operation. The experimental results clearly show that the ratio
between the forward and backward quantum thermodynamic processes obeys the
quantum description of microscopic reversibility rather than the classical
prediction. In particular, we observe that quantum coherences involved in both
initial and final states contribute to the modification of microscopic
reversibility, especially in the low-temperature limit. Our results show that
the quantum modification can aid the symmetry relation between the forward and
backward dynamics to understand classical and quantum fluctuation theorems in
a single framework, which can be applied to a wide range of physical
processes, for example, exploring the Markovian master equation derived by a
set of beam-splitters Kim95 .

###### Acknowledgements.

N.B., S.F., M.B., and A.Z. gratefully acknowledge the support of the EU under
the ERA-NET QuantERA project “ShoQC” (Grant No. 731473) and the FET Flagship
on Quantum Technologies project “Qombs” (Grant no. 820419). H.K. is supported
by the KIAS Individual Grant No. CG085301 at Korea Institute for Advanced
Study. M.S.K. acknowledges the KIST Open Research Program.

## References

  * (1) G. Gallavotti and E. G. D. Cohen, Dynamical Ensembles in Nonequilibrium Statistical Mechanics, Phys. Rev. Lett. 74, 2694 (1995).
  * (2) C. Jarzynski, Nonequilibrium Equality for Free Energy Differences, Phys. Rev. Lett. 78, 2690 (1997).
  * (3) G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Phys. Rev. E 60, 2721 (1999).
  * (4) J. Liphardt, S. Dumont, S. B. Smith, I. Tinoco, Jr., and C. Bustamante, Equilibrium Information from Nonequilibrium Measurements in an Experimental Test of Jarzynski’s Equality, Science 296, 1832 (2002).
  * (5) D. Collin, F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco, Jr., and C. Bustamante, Verification of the Crooks Fluctuation Theorem and Recovery of RNA Folding Free Energies, Nature (London) 437, 231 (2005).
  * (6) V. Blickle, T. Speck, L. Helden, U. Seifert, and C. Bechinger, Thermodynamics of a Colloidal Particle in a Time-Dependent Nonharmonic Potential, Phys. Rev. Lett. 96, 070603 (2006).
  * (7) O.-P. Saira, Y. Yoon, T. Tanttu, M. Möttönen, D. V. Averin, and J. P. Pekola, Test of the Jarzynski and Crooks Fluctuation Relations in an Electronic System, Phys. Rev. Lett. 109, 180601 (2012).
  * (8) M. Esposito, U. Harbola, and S. Mukamel, Nonequilibrium Fluctuations, Fluctuation Theorems, and Counting Statistics in Quantum Systems, Rev. Mod. Phys. 81, 1665 (2009).
  * (9) M. Campisi, P. Hänggi, and P. Talkner, Colloquium: Quantum Fluctuation Relations: Foundations and Applications, Rev. Mod. Phys. 83, 771 (2011).
  * (10) S. Deffner, M. Brunner and E. Lutz, Quantum fluctuation theorems in the strong damping limit, Europhys. Lett. 94, 30001 (2011).
  * (11) S. Deffner, Quantum Entropy Production in Phase Space, Europhys. Lett. 103, 30001 (2013).
  * (12) P. Hänggi and P. Talkner, The Other QFT, Nat. Phys. 11, 108 (2015).
  * (13) C. Jarzynski, H. T. Quan, and S. Rahav, Quantum-Classical Correspondence Principle for Work Distributions, Phys. Rev. X 5, 031038 (2015).
  * (14) K. Funo, M. Ueda, and T. Sagawa, Quantum Fluctuation Theorems in Thermodynamics in the Quantum Regime (Springer, 2018) pp. 249–273.
  * (15) P. Taranto, K. Modi, and F. A. Pollock, Emergence of a fluctuation relation for heat in nonequilibrium Landauer processes, Phys. Rev. E 97, 052111 (2018).
  * (16) J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, The Role of Quantum Information in Thermodynamics: A Topical Review, J. Phys. A 49, 143001 (2016).
  * (17) Á. M. Alhambra, L. Masanes, J. Oppenheim, and C. Perry, Fluctuating Work: From Quantum Thermodynamical Identities to a Second Law Equality, Phys. Rev. X 6, 041017 (2016).
  * (18) F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso, eds., Thermodynamics in the Quantum Regime (Springer International, 2019).
  * (19) G. T. Landi and M. Paternostro, Irreversible entropy production: From classical to quantum, Rev. Mod. Phys. 93, 035008 (2021).
  * (20) P. Talkner, E. Lutz, and P. Hänggi, Fluctuation Theorems: Work Is Not an Observable, Phys. Rev. E 75, 050102(R) (2007).
  * (21) M. N. Bera, A. Riera, M. Lewenstein, Z. B. Khanian and A. Winter, Thermodynamics as a Consequence of Information Conservation, Quantum 3, 121 (2019).
  * (22) M. O. Scully, M. S. Zubairy, G. S. Agarwal, and H. Walther, Extracting Work from a Single Heat Bath via Vanishing Quantum Coherence, Science 299, 862 (2003).
  * (23) I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. S. Kim, and Č. Brukner , Probing Planck-scale physics with quantum optics Nat. Phys. 8, 393 (2012).
  * (24) M. D. Vidrighin, O. Dahlsten, M. Barbieri, M. S. Kim, V. Vedral, and I. A. Walmsley, Photonic Maxwell’s Demon, Phys. Rev. Lett. 116, 050401 (2016).
  * (25) T. Albash, D. A. Lidar, M. Marvian, and P. Zanardi, Fluctuation theorems for quantum processes, Phys. Rev. E 88, 032146 (2013).
  * (26) T. B. Batalhão, A. M. Souza, L. Mazzola, R. Auccaise, R. S. Sarthour, I. S. Oliveira, J. Goold, G. De Chiara, M. Paternostro, and R. M. Serra, Experimental Reconstruction of Work Distribution and Study of Fluctuation Relations in a Closed Quantum System, Phys. Rev. Lett. 113, 140601 (2014).
  * (27) J. Åberg, Fully Quantum Fluctuation Theorems, Phys. Rev. X 8, 011019 (2018).
  * (28) G. Manzano, J. M. Horowitz, and J. M. R. Parrondo, Quantum Fluctuation Theorems for Arbitrary Environments: Adiabatic and Nonadiabatic Entropy Production, Phys. Rev. X 8, 031037 (2018).
  * (29) Z. Holmes, S. Weidt, D. Jennings, J. Anders, and F. Mintert, Coherent Fluctuation Relations: From the Abstract to the Concrete, Quantum 3, 124 (2019).
  * (30) H. Kwon and M. S. Kim, Fluctuation theorems for a quantum channel, Phys. Rev. X 9, 031029 (2019).
  * (31) K. Micadei, G. T. Landi, and E. Lutz, Quantum Fluctuation Theorems beyond Two-Point Measurements, Phys. Rev. Lett. 124, 090602 (2020).
  * (32) K. Micadei, J. P. S. Peterson, A. M. Souza, R. S. Sarthour, I. S. Oliveira, G. T. Landi, R. M. Serra, and E. Lutz, Experimental Validation of Fully Quantum Fluctuation Theorems Using Dynamic Bayesian Networks Phys. Rev. Lett. 127, 180603 (2021).
  * (33) J. Goold and K. Modi, Fluctuation theorem for nonunital dynamics AVS Quantum Sci. 3, 045001 (2021).
  * (34) T. Yada, N. Yoshioka, and T. Sagawa, Quantum Fluctuation Theorem under Quantum Jumps with Continuous Measurement and Feedback, Phys. Rev. Lett. 128, 170601 (2022).
  * (35) W. H. Zurek, S. Habib, and J. P. Paz, Coherent states via decoherence, Phys. Rev. Lett. 70, 1187 (1993).
  * (36) L. Mandel, Non-Classical States of the Electromagnetic Field, Phys. Scr. T12, 34 (1986).
  * (37) In the basis of position eigenstate $|x\rangle$ such that $\hat{x}|x\rangle=x|x\rangle$, the time-reversal of the state $|\psi\rangle=\int_{-\infty}^{\infty}dx\psi(x)|x\rangle$ becomes $|\bar{\psi}\rangle=\int_{-\infty}^{\infty}dx\psi^{*}(x)|x\rangle$.
  * (38) D. Janzing, P. Wocjan, R. Zeier, R. Geiss, and T. Beth, Thermodynamic Cost of Reliability and Low Temperatures: Tightening Landauer’s Principle and the Second Law, Int. J. Theor. Phys. 39, 2717 (2000).
  * (39) M. Horodecki and J. Oppenheim, Fundamental limitations for quantum and nanoscale thermodynamics, Nat. Commun. 4, 2059 (2013).
  * (40) M. Lostaglio, Á. M. Alhambra, and C. Perry, Elementary Thermal Operations, Quantum 2, 52 (2018).
  * (41) H. Kwon, H. Jeong, D. Jennings, B. Yadin, and M. S. Kim, Clock–Work Trade-Off Relation for Coherence in Quantum Thermodynamics, Phys. Rev. Lett. 120, 150602 (2018).
  * (42) Th. Richter, Determination of photon statistics and density matrix from double homodyne detection measurements, Journal of Modern Optics, 45:8, 1735-1749 (1998).
  * (43) F.T. Arecchi, Measurement of the Statistical Distribution of Gaussian and Laser Sources, Phys. Rev. Lett. 15, 912 (1965).
  * (44) A. Zavatta, V. Parigi, and M. Bellini, Experimental nonclassicality of single-photon-added thermal light states, Phys. Rev. A 75, 052106 (2007).
  * (45) V. Parigi, A. Zavatta, and M. Bellini, Implementation of single-photon creation and annihilation operators: experimental issues in their application to thermal states of light, J. Phys. B: At. Mol. Opt. Phys 42, 114005 (2009).
  * (46) Supplemental Material for detailed derivation and experimental setup.
  * (47) M. S. Kim and N. Imoto, Phase-sensitive reservoir modeled by beam splitters, Phys. Rev. A 52, 2401 (1995).

##  I Supplementary Material

###  I.1 Direct derivation of the quantum microscopic reversibility condition

The symmetry relation between the forward and backward transition
probabilities of quantum thermodynamic processes has been discussed in Refs.
Aberg18 ; Kwon19 . Here, we directly derive an explicit form of the quantum
microscopic reversibility condition for the transition between coherent
states,

$\frac{{\cal P}_{\rightarrow}(\alpha_{f}|\alpha_{i})}{{\cal
P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})}=\exp\left[-\frac{|\alpha_{f}|^{2}}{n_{\rm
th}+1}+\frac{|\alpha_{i}|^{2}}{n_{\rm th}}\right],$ (9)

where $\tilde{\alpha}_{i}=\bar{\alpha}_{i}e^{\beta\hbar\omega_{0}/2}$ and
$\tilde{\alpha}_{f}=\bar{\alpha}_{f}e^{-\beta\hbar\omega_{0}/2}$.

###### Proof.

We first note that a Gibbs-rescaled coherent state can be expressed as the
following form:

$|\alpha e^{\pm\beta\hbar\omega_{0}/2}\rangle={\cal
N}_{\pm}(\alpha)e^{\pm\beta\hat{H}/2}|\alpha\rangle,$

where
$\hat{H}=\hbar\omega_{0}\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right)$ and
${\cal
N}_{\pm}(\alpha)=\left(\langle\alpha|e^{\pm\beta\hat{H}}|\alpha\rangle\right)^{-\frac{1}{2}}$.
Using the Gibbs-rescaled coherent states, we express the backward transition
probability as

$\displaystyle{\cal P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})$
$\displaystyle={\mbox{Tr}}\left[\hat{U}_{R}\left(|\tilde{\alpha}_{f}\rangle\langle\tilde{\alpha}_{f}|\otimes\hat{\rho}_{\rm
th}\right)\hat{U}_{R}^{\dagger}\left(|\tilde{\alpha}_{i}\rangle\langle\tilde{\alpha}_{i}|\otimes\mathbb{\hat{1}}_{B}\right)\right]$
(10) $\displaystyle={\cal N}_{+}^{2}(\alpha_{i}){\cal
N}_{-}^{2}(\alpha_{f}){\mbox{Tr}}\left[\hat{U}_{R}\left(e^{-\beta\hat{H}/2}|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|e^{-\beta\hat{H}/2}\otimes\hat{\rho}_{\rm
th}\right)\hat{U}_{R}^{\dagger}\left(e^{\beta\hat{H}/2}|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|e^{\beta\hat{H}/2}\otimes\mathbb{\hat{1}}_{B}\right)\right].$

From the energy conservation condition $[\hat{H}_{I},\hat{H}+\hat{H}_{B}]=0$
and time-reversal symmetry of the system and bath Hamiltonians,
$H(\hat{x},\hat{p})=H(\hat{x},-\hat{p})$ and
$H_{B}(\hat{x}_{B},\hat{p}_{B})=H_{B}(\hat{x}_{B},-\hat{p}_{B})$, we observe
that
$[\hat{H}_{I}^{R},\hat{H}+\hat{H}_{B}]=0\Rightarrow[\hat{U}_{R},\hat{H}+\hat{H}_{B}]=0$.
We then express the quantum state following the reverse dynamics as

$\displaystyle\hat{U}_{R}\left(e^{-\beta\hat{H}/2}|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|e^{-\beta\hat{H}/2}\otimes\hat{\rho}_{\rm
th}\right)\hat{U}_{R}^{\dagger}$
$\displaystyle=\frac{1}{Z_{B}}\hat{U}_{R}\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\left(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B}\right)\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\hat{U}_{R}^{\dagger}$
$\displaystyle=\frac{1}{Z_{B}}\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\hat{U}_{R}\left(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B}\right)\hat{U}_{R}^{\dagger}\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right),$

where $\hat{\rho}_{\rm th}=\frac{e^{-\beta\hat{H}_{B}}}{Z_{B}}$ with
$Z_{B}={\rm Tr}\left[e^{-\beta\hat{H}_{B}}\right]$. By substituting this
formula to Eq. (10), we obtain

$\displaystyle{\cal P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})$
(11) $\displaystyle\quad=\frac{{\cal N}_{+}^{2}(\alpha_{i}){\cal
N}_{-}^{2}(\alpha_{f})}{Z_{B}}{\rm Tr}\left[\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\hat{U}_{R}\left(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B}\right)\hat{U}_{R}^{\dagger}\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\left(e^{\beta\hat{H}/2}|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|e^{\beta\hat{H}/2}\otimes\mathbb{\hat{1}}_{B}\right)\right]$
$\displaystyle\quad=\frac{{\cal N}_{+}^{2}(\alpha_{i}){\cal
N}_{-}^{2}(\alpha_{f})}{Z_{B}}{\rm
Tr}\left[\hat{U}_{R}^{\dagger}\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\left(e^{\beta\hat{H}/2}|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|e^{\beta\hat{H}/2}\otimes\mathbb{\hat{1}}_{B}\right)\left(e^{-\beta\hat{H}/2}\otimes
e^{-\beta\hat{H}_{B}/2}\right)\hat{U}_{R}(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B})\right]$
$\displaystyle\quad=\frac{{\cal N}_{+}^{2}(\alpha_{i}){\cal
N}_{-}^{2}(\alpha_{f})}{Z_{B}}{\rm
Tr}\left[\hat{U}_{R}^{\dagger}\left(|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|\otimes
e^{-\beta
H_{B}}\right)\hat{U}_{R}(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B})\right]$
$\displaystyle\quad={\cal N}_{+}^{2}(\alpha_{i}){\cal
N}_{-}^{2}(\alpha_{f}){\rm
Tr}\left[\hat{U}_{R}^{\dagger}\left(|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|\otimes\hat{\rho}_{\rm
th}\right)\hat{U}_{R}\left(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B}\right)\right].$

As both the thermal state and the identity operator are invariant under time-
reversal, $(\hat{x}_{B},\hat{p}_{B})\rightarrow(\hat{x}_{B},-\hat{p}_{B})$,
there exist the representations $\hat{\rho}_{\rm
th}=\sum_{\lambda}P_{\hat{\rho}_{\rm
th}}(\lambda)|\lambda\rangle\langle\lambda|=\sum_{\lambda}P_{\hat{\rho}_{\rm
th}}(\lambda)|\bar{\lambda}\rangle\langle\bar{\lambda}|$ and
$\mathbb{1}_{B}=\sum_{\lambda}P_{\mathbb{1}_{B}}(\lambda)|\bar{\lambda}\rangle\langle\bar{\lambda}|$,
where $\bar{\lambda}$ denotes time-reversal of the bath coordinates. This
observation leads to

$\displaystyle{\rm
Tr}\left[\hat{U}_{R}^{\dagger}\left(|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|\otimes\hat{\rho}_{\rm
th}\right)\hat{U}_{R}\left(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes\mathbb{1}_{B}\right)\right]$
$\displaystyle=\sum_{\lambda,\lambda^{\prime}}P_{\hat{\rho}_{\rm
th}}(\lambda)P_{\mathbb{1}_{B}}(\lambda^{\prime}){\rm
Tr}\left[\hat{U}_{R}^{\dagger}\left(|\bar{\alpha}_{i}\rangle\langle\bar{\alpha}_{i}|\otimes|\bar{\lambda}\rangle\langle\bar{\lambda}|\right)\hat{U}_{R}\left(|\bar{\alpha}_{f}\rangle\langle\bar{\alpha}_{f}|\otimes|\bar{\lambda}^{\prime}\rangle\langle\bar{\lambda}^{\prime}|\right)\right]$
(12) $\displaystyle=\sum_{\lambda,\lambda^{\prime}}P_{\hat{\rho}_{\rm
th}}(\lambda)P_{\mathbb{1}_{B}}(\lambda^{\prime})\langle\bar{\alpha}_{f},\bar{\lambda}^{\prime}|\hat{U}_{R}^{\dagger}|\bar{\alpha}_{i},\bar{\lambda}\rangle\langle\bar{\alpha}_{i},\bar{\lambda}|\hat{U}_{R}|\bar{\alpha}_{f},\bar{\lambda}^{\prime}\rangle$
$\displaystyle=\sum_{\lambda,\lambda^{\prime}}P_{\hat{\rho}_{\rm
th}}(\lambda)P_{\mathbb{1}_{B}}(\lambda^{\prime})\langle\alpha_{i},\lambda|\hat{U}^{\dagger}|\alpha_{f},\lambda^{\prime}\rangle\langle\alpha_{f},\lambda^{\prime}|\hat{U}|\alpha_{i},\lambda\rangle$
$\displaystyle={\rm
Tr}\left[\hat{U}\left(|\alpha_{i}\rangle\langle\alpha_{i}|\otimes\hat{\rho}_{\rm
th}\right)\hat{U}^{\dagger}\left(|\alpha_{f}\rangle\langle\alpha_{f}|\otimes\mathbb{1}_{B}\right)\right]$
$\displaystyle={\cal P}_{\rightarrow}(\alpha_{f}|\alpha_{i}).$

Here, the third equality is obtained from the time-reversal symmetry of the
quantum dynamics
$\left|\langle\bar{\alpha}_{i},\bar{\lambda}|\hat{U}_{R}|\bar{\alpha}_{f},\bar{\lambda}^{\prime}\rangle\right|^{2}=\left|\langle\alpha_{f},\lambda^{\prime}|\hat{U}|\alpha_{i},\lambda\rangle\right|^{2}$.
Finally, by rearranging the expressions in Eqs. (11) and (12), we obtain

$\frac{{\cal P}_{\rightarrow}(\alpha_{f}|\alpha_{i})}{{\cal
P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})}=\frac{1}{{\cal
N}_{+}^{2}(\alpha_{i}){\cal
N}_{-}^{2}(\alpha_{f})}=\langle\alpha_{i}|e^{\beta\hat{H}}|\alpha_{i}\rangle\langle\alpha_{f}|e^{-\beta\hat{H}}|\alpha_{f}\rangle.$
(13)

An explicit form of the right-hand side of the equation can be evaluated as

$\displaystyle\langle\alpha_{i}|e^{\beta\hat{H}}|\alpha_{i}\rangle$
$\displaystyle=e^{-|\alpha_{i}|^{2}(1-\exp[\beta\hbar\omega_{0}])}=e^{\frac{|\alpha_{i}|^{2}}{n_{\rm
th}}}$ $\displaystyle\langle\alpha_{f}|e^{-\beta\hat{H}}|\alpha_{f}\rangle$
$\displaystyle=e^{-|\alpha_{f}|^{2}(1-\exp[-\beta\hbar\omega_{0}])}=e^{-\frac{|\alpha_{f}|^{2}}{n_{\rm
th}+1}},$

where $\frac{1}{n_{\rm th}}=e^{\beta\hbar\omega_{0}}-1$ and $\frac{1}{n_{\rm
th}+1}=1-e^{-\beta\hbar\omega_{0}}$ by defining $n_{\rm
th}=\frac{e^{-\beta\hbar\omega_{0}}}{1-e^{-\beta\hbar\omega_{0}}}$. By noting
that $Q=\Delta
E=\langle\alpha_{f}|\hat{H}|\alpha_{f}\rangle-\langle\alpha_{i}|\hat{H}|\alpha_{i}\rangle$,
we derive a quantum correction,

$\frac{{\cal P}_{\rightarrow}(\alpha_{f}|\alpha_{i})}{{\cal
P}_{\leftarrow}(\tilde{\alpha}_{i}|\tilde{\alpha}_{f})}=\Upsilon e^{-\beta
Q},$

with

$\displaystyle\Upsilon$
$\displaystyle=\frac{\langle\alpha_{i}|e^{\beta\hat{H}}|\alpha_{i}\rangle}{e^{\beta\langle\alpha_{i}|\hat{H}|\alpha_{i}\rangle}}\frac{\langle\alpha_{f}|e^{-\beta\hat{H}}|\alpha_{f}\rangle}{e^{-\beta\langle\alpha_{f}|\hat{H}|\alpha_{f}\rangle}}$
$\displaystyle=\exp\left[{-|\alpha_{i}|^{2}(1-e^{\beta\hbar\omega_{0}}+\beta\hbar\omega_{0})-|\alpha_{f}|^{2}(1-e^{-\beta\hbar\omega_{0}}-\beta\hbar\omega_{0})}\right]$
$\displaystyle=\exp[(\cosh(\beta\hbar\omega_{0})-1)(|\alpha_{i}|^{2}+|\alpha_{f}|^{2})-(\sinh(\beta\hbar\omega_{0})-\beta\hbar\omega_{0})\left(|\alpha_{f}|^{2}-|\alpha_{i}|^{2}\right)].$

∎

###  I.2 Experimental setup

The experimental setup, shown in Figure 6, is based on a mode-locked Ti:Sa
laser emitting 1.5 ps-long pulses at 786 nm with a repetition rate of 81 MHz.
Its output is split into three parts. The first one provides the local
oscillator (LO) pulses for heterodyne detection. The second part is employed
to produce a pseudo-thermal state of light. A lens (L) focuses the beam onto a
rotating ground glass disk (RD) Arecchi1965 . Collecting a small portion of
the light scattered by the disk into a single-mode fiber, we obtain a thermal
state in a single spatial mode Zavatta2007 ; Parigi2009 . The third beam gives
the coherent states. Two variable attenuators (not shown in the figure),
placed along the thermal and coherent state paths, are used to control the
thermal mean photon number $\bar{n}_{\rm th}$ and the coherent state amplitude
$\alpha$.

Figure 6: Scheme of the experimental setup. The box explains the symbols
employed to indicate the polarization direction. Red arrows refer to LO, blue
arrows refer to signal. Abbreviations: L: lens, RD: rotating disk,
$\lambda/2$: half-wave plate, PBS: polarizing beam-splitter, P: polarizer, BC:
Babinet compensator, HD: homodyne detector.

Thermal and coherent states are mixed in a polarization-based beam-splitter
with variable reflectivity composed as follows. The thermal states are
vertically polarized while the coherent states are horizontally polarized;
their spatial modes are mixed via a polarizing beam-splitter (PBS). Then, a
half-wave plate (HWP) rotates their polarization before entering into a
polarizer (P), diagonally oriented. The combination of these three devices
plays the role of a beam-splitter with adjustable reflectivity, which allows
us to control the ratio of coherent and thermal states simulating the
thermodynamic process.

Throughout the experiment, the reflectivity of the beam-splitter mixing
coherent and thermal states has been set to two different values: 85% and 70%,
depending on the coherent state amplitudes and mean thermal photon numbers
used. This is due to the need of a sufficiently large mean field amplitude of
the states reaching the two homodyne detectors, for a reliable lock of the
relative signal-LO phase and in order to remain in the detectors’ dynamic
range. Experimental results confirm that the ${\cal P}_{\rightarrow}/{\cal
P}_{\leftarrow}$ ratio is independent on the beam-splitter reflectivity.

The detection is performed via a polarization-based heterodyne detector. The
signal, composed by the mix of thermal and coherent states, and the LO enter
the central PBS with diagonal polarizations. Therefore, they are equally split
between the two arms of the heterodyne. In each arm a HWP, which rotates the
polarization back to diagonal, and a second PBS play the role of the balanced
beam-splitter of the homodyne detector (HD). A crucial parameter of the
heterodyne detector is the phase relationship between the quadrature measured
by HD1 and the one measured by HD2, which must be exactly $\pi/2$. For this
reason we inserted a Babinet compensator (BC) in the path to HD1, which allows
us to fine tune the delay between signal and LO.

We first use coherent states with real values of $\alpha$. This simplifies the
procedure since the Gibbs rescaled $\tilde{\alpha}$ remains real and thus all
the $Q$-functions are centred on the $\hat{x}$ quadrature axis. Then, we also
test the use of coherent states with complex amplitudes (data points
corresponding to a mean thermal photon number $n_{th}=1.22$ in Fig. 3 of the
main text). In such a case, the Gibbs rescaled values are complex conjugated
and thus the $Q$-functions of the forward and time-reversed process are
located in two different quadrants of the phase space, as shown in Fig. 3 of
the main text. For the sake of simplicity, when starting from an initial
coherent state with complex amplitude $\alpha_{i}=|\alpha|\exp(i\phi)$, we
evaluate the $Q$-function at values of $\alpha_{f}$ having the same phase of
$\alpha_{i}$ ($\theta_{i}=\theta_{f}=\phi$). This simplifies the experimental
procedure since the LO phase has to be locked only to two values:
$\theta_{i}=\phi$ for the forward process, and $\tilde{\theta}_{f}=-\phi$ for
the time-reversed one.

$\hat{x}$ and $\hat{p}$ quadrature data are simultaneously collected from the
two HDs, with the global LO phase locked via a piezo-mounted mirror. These
data directly provide the Husimi $Q$-function Richter98 , whose distribution
is estimated via a maximum likelihood algorithm after 50000 $\hat{x}$ and
$\hat{p}$ quadrature measurements. The confidence intervals of the
distribution parameters are calculated using a bootstrap method: 1000 random
values are sampled from the distribution and are employed to evaluate again
the distribution parameters; the procedure is repeated 1000 times and the
confidence intervals are retrieved from the standard deviation of the
parameters.


# Mixed-Integer Linear Optimization for Cardinality-Constrained Random Forests

Jan Pablo Burgard, Maria Eduarda Pinheiro, Martin Schmidt  Trier University,
Department of Economic and Social Statistics, Universitätsring 15, 54296
Trier, Germany<EMAIL_ADDRESS>Trier University, Department of
Mathematics, Universitätsring 15, 54296 Trier, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Date: August 31, 2024)

###### Abstract.

Random forests are among the most famous algorithms for solving classification
problems, in particular for large-scale data sets. Considering a set of
labeled points and several decision trees, the method takes the majority vote
to classify a new given point. In some scenarios, however, labels are only
accessible for a proper subset of the given points. Moreover, this subset can
be non-representative, e.g., due to collection bias. Semi-supervised learning
considers the setting of labeled and unlabeled data and often improves the
reliability of the results. In addition, it can be possible to obtain
additional information about class sizes from undisclosed sources. We propose
a mixed-integer linear optimization model for computing a semi-supervised
random forest that covers the setting of labeled and unlabeled data points as
well as the overall number of points in each class for a binary
classification. Since the solution time rapidly grows as the number of
variables increases, we present some problem-tailored preprocessing techniques
and an intuitive branching rule. Our numerical results show that our approach
leads to a better accuracy and a better Matthews correlation coefficient for
biased samples compared to random forests by majority vote, even if only few
labeled points are available.

###### Key words and phrases:

Random forests, Semi-supervised learning, Mixed-integer linear optimization,
Preprocessing, Cardinality constraints

###### 2020 Mathematics Subject Classification:

90C11,90C90,90-08,68T99

##  1\. Introduction

Random forests are one of the most famous approaches in supervised learning
[3]. It has been applied to various fields such as the prediction of diseases
[23, 12], 3D object recognition [25] and Fraude and accident detection [27,
10]. The main reasons why random forests are popular are that they prevent
over-fitting [13], that they have only a few parameters to tune, and that they
can be used directly for high-dimensional problems [9, 2]. The core idea is,
given labeled data, to combine the prediction of different trees, in general,
using the majority vote to classify new points.

Nevertheless, acquiring labels for every unit of interest can be costly—in
particular when classic surveys are used to obtain the labels. In this
situation, it would be beneficial to train the random forest with only partly
labeled data. This yields a semi-supervised learning setting [30]. Algorithms
for semi-supervised learning have already been proposed for neural networks
[16, 21, 20], logistic regression [1, 7], support vector machines [8, 19], and
decision trees [15, 29, 14].

In the case of random forests, [17] propose an iterative and deterministic
annealing-like training algorithm that maximizes the multi-class margin of
labeled and unlabeled samples. Furthermore, [18] extend the co-training
paradigm to random forests, determining how certain the model is about its
predictions for unlabeled data. Moreover, [28] combine active learning and
semi-supervised learning to improve the final classification performance of
random forests by utilizing supervised clustering to categorize the unlabeled
data.

However, in many applications, it is possible to know the total amount of
elements in each class within a population, e.g., when external sources
provide this information. For instance, a company might only know the overall
number of successful transactions, but might not be able to identify which
specific customer’s transactions were successful. An intuitive example is an
online retailer that may track the total number of good customer reviews but
does not have access to individual ratings due to anonymity practices. Another
example is from healthcare, where it is possible to know how many patients
have a disease but, due to data privacy reasons, one does not know which
specific person is affected or not. [4] propose aggregating this extra
information for logistic regression. They develop a cardinality-constrained
multinomial logit model. For support vector machines, [6] present a mixed-
integer quadratic optimization model and iterative clustering techniques to
tackle cardinality constraints for each class. Moreover, for the case of
decision trees, [5] propose a mixed-integer linear optimization model for
computing semi-supervised optimal classification trees that serve the same
purpose.

Our contribution here is to propose a random forest model that imposes a
cardinality constraint on the classification of the unlabeled data. We develop
a big-$M$-based mixed-integer linear programming (MILP) model to solve the
cardinality-constrained random forest (C2RF) problem that includes the
cardinality constraint for the unlabeled data. The cardinality constraint
helps to account for biased samples since the number of predictions in each
class on the population is bounded by the constraint. In particular, our
numerical results show that our approach leads to a better accuracy and a
better Matthews correlation coefficient for biased samples compared to random
forests by majority vote, even if only few labeled points are available. The
computation time for this MILP grows with the number of variables—especially
for an increasing number of integer variables. To account for this, we present
theoretical results that lead to preprocessing techniques that significantly
reduce the computation time.

This paper is organized as follows. In Section 2 we present the optimization
model and prove the correctness of the used big-$M$ parameter. Afterward, the
preprocessing techniques are discussed in Section 3 and an intuitive branching
rule is presented in Section 4. There, we also present our algorithm that
combines the mentioned techniques and the MILP formulation. In Section 5 we
report and discuss numerical results. Finally, we conclude in Section 6.

##  2\. An MILP Formulation for Cardinality-Constrained Random Forests

Let $X\in{p\times N}=[X_{u},X_{l}]$ be the data matrix with unlabeled data
$X_{u}=[x^{1},\dotsc,x^{m}]$ and labeled data $X_{l}=[x^{m+1},\dotsc,x^{N}]$.
Hence, we are given points $x^{i}\in p$ for all $i\leavevmode\nobreak\
\in\leavevmode\nobreak\ [1,N]\mathrel{{\mathop{:}}{=}}\\{1,\dots,N\\}$. We set
$n\mathrel{{\mathop{:}}{=}}N-m$ and $y\in\\{-1,1\\}^{n}$ as the vector of
class labels for the labeled data. Let $t$ be the number of given decision
trees and let $A^{j}\in{p\times d}$ be a subset of the labeled data with size
$d$ for $j\in[1,t]$. For each $j\in[1,t]$, based on each column of $A^{j}$ and
its label, the $j$th tree generates a vector $r^{j}\in\\{-1,1\\}^{m}$ that
classifies the unlabeled data $X_{u}$. Thus, for each unlabeled point $x^{i}$
we observe a vector of classification
$r_{i}=[r^{1}_{i},\dots,r^{t}_{i}]\in\\{-1,1\\}^{t}$. Hence,
$R=[r_{1},\dotsc,r_{m}]\in\\{-1,1\\}^{t\times m}$ and $r_{i}^{j}$ is the
classification of $x^{i}\in X_{u}$ given by the tree $j$. In a random forest,
the prediction for a point $x^{i}\in X_{u}$ is the dominant class chosen by
the individual $t$ trees, i.e., the majority vote.

In many applications, aggregated information on the labels is available, e.g.,
from census data. For what follows, we assume to know the total number
$\lambda\in
ofunlabeledpointsthatbelongtothepositiveclassandproposeamodelsuchthatwecanusealinearcombinationofthetreeclassificationsaswellas$λ$asanadditionalinformation.Ourgoalistofindoptimalparameters$α^*
∈t,η^* ∈, and $z^{*}\in\\{0,1\\}^{m}$ that solve the optimization problem

$\displaystyle\min_{\alpha,\eta,z}\quad$ $\displaystyle\eta$ (0a) s.t.
$\displaystyle\alpha^{\top}r_{i}\leq-1+z_{i}M,\quad i\in[1,m],$ (0b)
$\displaystyle\alpha^{\top}r_{i}\geq 1-(1-z_{i})M,\quad i\in[1,m],$ (0c)
$\displaystyle\lambda-\eta\leq\sum_{i=1}^{m}z_{i}\leq\lambda+\eta,$ (0d)
$\displaystyle\ell\leq\alpha_{j}\leq u,\quad j\in[1,t],$ (0e) $\displaystyle
0\leq\eta\leq\bar{\eta},$ (0f) $\displaystyle z_{i}\in\\{0,1\\},\quad
i\in[1,m],$ (0g)

where $M$ needs to be chosen sufficiently large, $u>\ell>0$ holds, and we set

$\bar{\eta}\mathrel{{\mathop{:}}{=}}\max\left\\{\lambda,m-\lambda\right\\}.$
(1)

Note that the objective function in (0a) minimizes the classification error
for the unlabeled data. As $z_{i}$ is binary, Constraints (0b) and (0c) lead
to

$\displaystyle\alpha^{\top}r_{i}\geq 1\implies z_{i}=1,\quad i\in[1,m],$
$\displaystyle\alpha^{\top}r_{i}\leq-1\implies z_{i}=0,\quad i\in[1,m].$

Constraint (0d) ensures that the number of unlabeled data points classified as
positive is as close to $\lambda$ as possible. Constraint (0e) bounds the
weight of each tree’s decision for the final classification. This means that
for $j\in[1,t]$, as $\alpha_{j}$ gets closer to $u$, the $j$th tree gets
greater influence on the final classification, and as $\alpha_{j}$ gets closer
to $\ell$, the $j$th tree has less influence on the final classification.
Observe that since $\alpha_{j}\geq\ell>0$ holds for all $j\in[1,t]$, all trees
contribute to the final classification. Moreover, if $\alpha_{j}$ has the same
value for all $j\in[1,t]$, all trees contribute equally to the final
classification and we are in the standard random forest setup with majority
vote. Note that the upper bound $u$ is not necessary for the correctness of
the model but will serve as a big-$M$-type parameter as can be seen in
Proposition 1 below. The upper bound $\bar{\eta}$ in Constraint (0f) is also
not necessary for the correctness of Model (2). Nevertheless, as can be seen
in Proposition 2, this upper bound does not cut off any solution. Hence, we
include it in our implementation because we expect that the solution process
benefits from tight bounds. Problem (2) is an MILP. We refer to this problem
as C2RF (Cardinality-Constrained Random Forest). As usual for big-$M$
formulations, the choice of $M$ is crucial. If $M$ is too small, the problem
can become infeasible or optimal solutions could be cut off. If $M$ is chosen
too large, the respective continuous relaxations usually lead to bad lower
bounds and solvers may encounter numerical troubles. The choice of $M$ is
discussed in the following proposition.

######  Proposition 1.

A valid big-$M$ for Problem (2) is given by $M=ut+1$, i.e., $M$ is linear in
the number of trees in the forest.

###### Proof.

For all $i\in[1,m]$ we have $r_{i}\in\\{-1,1\\}^{t}$. Moreover, Constraint
(0e) ensures that $\alpha_{j}\leq u$ holds for all $j\in[1,t]$. Hence,

$\alpha^{\top}r_{i}\leq\sum_{j=1}^{t}\alpha_{j}\leq ut$

and

$\alpha^{\top}r_{i}\geq-\sum_{j=1}^{t}\alpha_{j}\geq-ut$

hold for all $i\in[1,m]$ and $M=ut+1$ does not cut any feasible solution. ∎

The following proposition makes a statement about the upper bound $\bar{\eta}$
in Constraint (0f).

######  Proposition 2.

Consider Problem (2) in which Constraint (0f) is replaced by $\eta\geq 0$.
Then, for every $\eta^{*}$ as being part of an optimal solution, it holds
$\eta^{*}\leq\bar{\eta}$ for $\bar{\eta}$ as defined in (1).

###### Proof.

Observe that since $z_{i}\in\\{0,1\\}$ for all $i\in[1,m]$,

$0\leq\sum_{i=1}^{m}z_{i}\leq m$

holds. Moreover, the maximum value occurs if all points are classified as
positive. If this happens, from Constraint (0d) we obtain

$\eta\geq\sum_{i=1}^{m}z_{i}-\lambda=m-\lambda.$

On the other hand, the minimum value of $\sum_{i=1}^{m}z_{i}$ occurs if all
points are classified as negative. If this happens, from Constraint (0d) we
obtain

$\eta\geq\sum_{i=1}^{m}z_{i}+\lambda=\lambda.$

Since Problem (2) is a minimization Problem, $\eta\leq\bar{\eta}$ holds. Thus,
the upper bound $\bar{\eta}$ in Constraint (0f) does not cut off any optimal
point. ∎

##  3\. Preprocessing

In this section, we present different preprocessing techniques for Problem (2)
that can be used to reduce the number of binary as well as the number of
continuous variables.

The first insight is that, if all trees have the same classification for some
unlabeled points, these points must have the same final classification and,
therefore, the respective binary variables always have the same values.

######  Proposition 3.

Let $k\in[1,m]$ and consider
$\mathcal{K}\mathrel{{\mathop{:}}{=}}\\{i\in[1,m]\colon r_{i}=r_{k}\\}$. Then,
$(\alpha,\eta,z)$ is a feasible point of Problem (2) if and only if there
exists a vector $\bar{z}\in\\{0,1\\}^{m+1-|\mathcal{K}|}$ such that
$(\alpha,\eta,\bar{z})$ is a feasible point of the problem

$\displaystyle\min_{\alpha,\eta,z}\quad$ $\displaystyle\eta$ (1a) s.t.
$\displaystyle\alpha^{\top}r_{i}\leq-1+z_{i}M,\quad
i\in\\{k\\}\cup[1,m]\setminus\mathcal{K},$ (1b)
$\displaystyle\alpha^{\top}r_{i}\geq 1-(1-z_{i})M,\quad
i\in\\{k\\}\cup[1,m]\setminus\mathcal{K},$ (1c)
$\displaystyle\lambda-\eta\leq\sum_{i\in[1,m]\setminus\mathcal{K}}z_{i}+|\mathcal{K}|z_{k}\leq\lambda+\eta,$
(1d) $\displaystyle\eqref{eq:alpha},\eqref{eq:eta2}$ (1e) $\displaystyle
z_{i}\in\\{0,1\\},\quad i\in\\{k\\}\cup[1,m]\setminus\mathcal{K}.$ (1f)

###### Proof.

Consider $(\alpha,\eta,z)$ a feasible point of (2) and

$\bar{z}_{i}=z_{i},\quad i\in\\{k\\}\cup[1,m]\setminus\mathcal{K}.$

We now prove that $(\alpha,\eta,\bar{z})$ is a feasible point of (3). Because
(0b), (0c), and (0g) hold, (1b), (1c), and (1f) are satisfied. Moreover, since
for all $i\in\mathcal{K}$ it holds $r_{i}=r_{k}$, we obtain that
$\alpha^{\top}r_{k}=\alpha^{\top}r_{i}$ holds for all $i\in\mathcal{K}$.
Hence, by Constraint (0b) and (0c), we obtain that $z_{i}=z_{k}$ also holds
for all $i\in\mathcal{K}$. This together with $\bar{z}_{k}=z_{k}$ implies that
$\lvert\mathcal{K}\rvert\bar{z}_{k}=\sum_{i\in\mathcal{K}}z_{i}$ is satisfied.
Hence,

$\sum_{i\in[1,m]\setminus\mathcal{K}}\bar{z}_{i}+|\mathcal{K}|\bar{z}_{k}=\sum_{i\in[1,m]\setminus\mathcal{K}}z_{i}+\sum_{i\in\mathcal{K}}z_{i}=\sum_{i=1}^{m}z_{i}$
(2)

is also satisfied, and, by Constraint (0d), we obtain that Constraint (1d)
holds. Therefore, $(\alpha,\eta,\bar{z})$ is a feasible point of Problem (3).

On the other hand, let $(\alpha,\eta,\bar{z})$ be a feasible point of Problem
(3) and set

$z_{i}=\begin{cases}\bar{z}_{i},\text{ if }i\notin\mathcal{K},\\\
\bar{z}_{k},\text{ otherwise}.\end{cases}$ (3)

Since (1b), (1c) and (1f) are satisfied, (0b) and (0c) hold for each
$i\notin\mathcal{K}$ and (0g) holds for all $i\in[1,m]$. Further, because each
$i\in\mathcal{K}$ satisfies $r_{i}=r_{k}$,
$\alpha^{\top}r_{i}=\alpha^{\top}r_{k}$ holds for all $i\in\mathcal{K}$.
Hence, by Constraints (1b) and (1c) we obtain that

$1-(1-z_{i})M=1-(1-\bar{z}_{k})M\leq\alpha^{\top}r_{i}\leq-1+\bar{z}_{k}M=-1+z_{i}M$

is satisfied for all $i\in\mathcal{K}$. Besides that, Expression (3) implies
that $|\mathcal{K}|\bar{z}_{k}=\sum_{i\in\mathcal{K}}z_{i}$ and, therefore,
Expression (2) also holds. Hence, by Constraint (1d), we obtain that
Constraint (0d) is satisfied. Therefore, $(\alpha,\eta,z)$ is a feasible point
of Problem (2). ∎

The following proposition shows that if one or more trees classify all points
exactly as another tree, some continuous variables of Problem (2) can be
eliminated.

######  Proposition 4.

Given $g\in[1,t]$, consider
$\mathcal{G}\mathrel{{\mathop{:}}{=}}\\{j\in[1,t]\colon r^{g}=r^{j}\\}$. Then,
$(\alpha^{*},\eta^{*},z^{*})$ is a solution to Problem (2) if and only if
there exists a vector $\bar{\alpha}\in\mathbb{R}^{t+1-|\mathcal{G}|}$ such
that $(\bar{\alpha},\eta^{*},z^{*})$ is a solution to the problem

$\displaystyle\min_{\alpha,\eta,z}\quad$ $\displaystyle\eta\ $ (3a) s.t.
$\displaystyle\sum_{j\in[1,t]\setminus\mathcal{G}}\alpha_{j}r_{i}^{j}+|\mathcal{G}|\alpha_{g}r_{i}^{g}\leq-1+z_{i}M,\quad
i\in[1,m],$ (3b)
$\displaystyle\sum_{j\in[1,t]\setminus\mathcal{G}}\alpha_{j}r_{i}^{j}+|\mathcal{G}|\alpha_{g}r_{i}^{g}\geq
1-(1-z_{i})M,\quad i\in[1,m],$ (3c) $\displaystyle\eqref{unlabeled part2},$
(3d) $\displaystyle\ell\leq\alpha_{j}\leq u,\quad
j\in\\{g\\}\cup[1,t]\setminus\mathcal{G},$ (3e)
$\displaystyle\eqref{eq:eta2},\eqref{eq:binaryz}.$ (3f)

###### Proof.

Let $(\alpha^{*},\eta^{*},z^{*})$ be a solution to Problem (2) and

$\bar{\alpha}_{j}=\begin{cases}\alpha_{j}^{*},&\text{if
}j\notin\mathcal{G},\\\
\sum_{j\in\mathcal{G}}\alpha_{j}^{*}/|\mathcal{G}|,&\text{otherwise.}\end{cases}$

Since $\ell\leq\alpha_{j}^{*}\leq u$ holds for all $j\in[1,t]$, we obtain

$\ell=\frac{\ell}{|\mathcal{G}|}(|\mathcal{G}|)\leq\bar{\alpha}_{g}\leq\frac{u}{|\mathcal{G}|}(|\mathcal{G}|)=u.$

Moreover, because $r^{g}=r^{j}$ is satisfied for all $j\in\mathcal{G}$,

$\sum_{j\in\mathcal{G}}\alpha_{j}^{*}r_{i}^{j}=r_{i}^{g}\sum_{j\in\mathcal{G}}\alpha_{j}^{*}=|\mathcal{G}|\bar{\alpha}_{g}r_{i}^{g}$
(4)

holds for all $i\in[1,m]$. Hence, for all $i\in[1,m]$,

$\sum_{j=[1,t]\setminus\mathcal{G}}\bar{\alpha}_{j}r_{i}^{j}+|\mathcal{G}|\bar{\alpha}_{g}r_{i}^{g}=\sum_{j=[1,t]\setminus\mathcal{G}}\alpha_{j}^{*}r_{i}^{j}+\sum_{j\in\mathcal{G}}\alpha_{j}^{*}r_{i}^{j}=(\alpha^{*})^{\top}r_{i}$
(5)

is satisfied and, consequently,

$1-(1-z_{i}^{*})M\leq\sum_{j=[1,t]\setminus\mathcal{G}}\bar{\alpha}_{j}r_{i}^{j}+|\mathcal{G}|\bar{\alpha}_{g}r_{i}^{g}\leq-1+z_{i}^{*}M$

holds for $i\in[1,m]$. Therefore, $(\bar{\alpha},\eta^{*},z^{*})$ is a
solution to Problem (4).

On the other hand, let $(\bar{\alpha},\eta^{*},z^{*})$ be a solution to
Problem (4) and set

$\alpha^{*}_{j}=\begin{cases}\bar{\alpha}_{j},&\text{if
}j\notin\mathcal{G},\\\ \bar{\alpha}_{g},&\text{otherwise}.\end{cases}$

Since (3e) holds, (0e) is satisfied for all $j\in[1,t]$. Besides that, since
$r^{g}=r^{j}$ is satisfied for all $j\in\mathcal{G}$, (4) and (5) also hold
for all $i\in[1,m]$. Hence, for all $i\in[1,m]$, we have

$1-(1-z_{i}^{*})M\leq(\alpha^{*})^{\top}r_{i}\leq-1+z_{i}^{*}M$

and $(\alpha^{*},\eta^{*},z^{*})$ is a solution to Problem (2). ∎

Finally, the following proposition allows to fix some binary variables $z_{i}$
of Problem (2).

######  Proposition 5.

For each $i\in[1,m]$, consider $\mathcal{A}_{i}=\\{j\in[1,t]\colon
r_{i}^{j}=-1\\}$ and $\mathcal{B}_{i}=\\{j\in[1,t]:r_{i}^{j}=1\\}$. If for
some $i\in[1,m]$,

$\varphi_{i}\mathrel{{\mathop{:}}{=}}-u\lvert\mathcal{A}_{i}\rvert+\ell\lvert\mathcal{B}_{i}\rvert\geq
1$ (6)

holds, then every feasible point $(\alpha,\eta,z)$ of Problem (2) satisfies
$z_{i}=1$. If, on the other hand,

$\phi_{i}\mathrel{{\mathop{:}}{=}}-\ell\lvert\mathcal{A}_{i}\rvert+u\lvert\mathcal{B}_{i}\rvert\leq-1$
(7)

is satisfied for some $i\in[1,m],$ then any feasible point $(\alpha,\eta,z)$
of Problem (2) satisfies $z_{i}=0$.

###### Proof.

Since $\ell\leq\alpha_{j}\leq u$ is satisfied for all $j\in[1,t]$, if for some
$i\in[1,m]$,

$-u|\mathcal{A}_{i}|+\ell|\mathcal{B}_{i}|\geq 1$

holds, we obtain

$\alpha^{\top}r_{i}=-\sum_{j\in\mathcal{A}_{i}}\alpha_{j}+\sum_{j\in\mathcal{B}_{i}}\alpha_{j}\geq-u|\mathcal{A}_{i}|+\ell|\mathcal{B}_{i}|\geq
1,$

and by Constraint (0b), $z_{i}=1$. On the other hand, if for some $i\in[1,m]$,
it holds

$-\ell|\mathcal{A}_{i}|+u|\mathcal{B}_{i}|\leq-1,$

we get

$\alpha^{\top}r_{i}=-\sum_{j\in\mathcal{A}_{i}}\alpha_{j}+\sum_{j\in\mathcal{B}_{i}}\alpha_{j}\leq-\ell|\mathcal{A}_{i}|+u|\mathcal{B}_{i}|\leq-1,$

and by Constraint (0c), $z_{i}=0$. ∎

Consider now

$\mathcal{P}\mathrel{{\mathop{:}}{=}}\left\\{i\in[1,m]\colon\varphi_{i}\geq
1\right\\}\quad\text{and}\quad\mathcal{N}\mathrel{{\mathop{:}}{=}}\left\\{i\in[1,m]\colon\phi_{i}\leq-1\right\\}.$

Then, $\lvert\mathcal{P}\rvert+\lvert\mathcal{N}\rvert$ binary variables can
be fixed. Moreover, $\lvert\mathcal{P}\rvert$ points then are already
classified as positive. If $|{\mathcal{P}}|\geq\lambda$, due to cardinality
constraint, all remaining points $x^{i}\in
X_{u}\setminus(\mathcal{P}\cup\mathcal{N})$ must be classified as negative,
and $\lambda$ can be set to $0$. On the other hand, if
$\lvert\mathcal{P}\rvert<\lambda$, only $\lambda-\lvert\mathcal{P}\rvert$
points in $X_{u}\setminus(\mathcal{P}\cup\mathcal{N})$ must be classified as
positive. This update is present in Step 1 in Algorithm 1 below.

##  4\. Branching Priorities

One aspect that can significantly affect the performance of MILP solvers is
the applied branching rule. In this brief section, we present a problem-
tailored rule for helping the MILP code to solve Problem (2). To this end, let
us consider binary variables $z_{i},z_{k}\in\\{0,1\\}$, $i,k\in[1,m]$, and
positive integer values $\xi_{i}$ and $\xi_{k}$ so that $\xi_{i}>\xi_{k}$
implies that the solver should branch on $z_{i}$ before $z_{k}$. In our
context, a point for which the percentage of trees that classify the point as
positive (or negative) is larger than for another point seems to be “easier”
to classify. Hence, we want to branch on the respective binary variable first.
Based on that, we establish a criterion for a branching strategy. We set
$\theta_{i}=\lvert\text{mean}(r_{i})\rvert$ for each $x^{i}\in
X_{u}\setminus(\mathcal{P}\cup\mathcal{N})$. Observe that the higher the value
of $\theta_{i}$, the more trees classify the point $x^{i}$ in one specific
class. Hence, we consider $\xi_{i}$ the position of $\theta_{i}$ in the vector
of the increasingly sorted values of $\theta$. Thus, the higher the value of
$\theta_{i}$, the higher the value of $\xi_{i}$, and hence, the higher the
branching priority for the binary variable $z_{i}$.

Motivated by the preprocessing techniques presented in Section 3 and the
branching priority discussed in this section, we obtain Algorithm 1 to solve
Problem (2).

Input : $R\in\\{-1,1\\}^{t\times m},u>\ell>0,$ $\lambda\in\mathbb{N}$,
$\mathcal{K}=\emptyset$, $\beta=0$, and $\gamma=0$.

1 Compute $M=ut+1$ and
$\bar{R}=[\bar{r}_{1},\dotsc,\bar{r}_{h}]\in\\{-1,1\\}^{t\times h}$ being the
set of all different $r_{i}\in R.$

2 for  _$k\in\\{1,\dotsc,h\\}$ _ do

3 Compute $w_{k}=\left|\\{i\in[1,m]:r_{i}=\bar{r}_{k}\\}\right|$,
$\varphi_{k}$ as described in (6), and $\phi_{k}$ as described in (7).

4 if  _$\varphi_{k}\geq 1$_ then

5 Set $\mathcal{K}\leftarrow\mathcal{K}\cup\\{k\\}$ and
$\beta\leftarrow\beta+w_{k}$.

6 else if  _$\phi_{k}\leq-1$_ then

7 Set $\mathcal{K}\leftarrow\mathcal{K}\cup\\{k\\}$ and
$\gamma\leftarrow\gamma+w_{k}$.

8 end if

9

10 end for

11Compute $S=[s^{1},\dotsc,s^{q}]^{\top}\in\\{-1,1\\}^{q\times h}$ being the
set of all different $\bar{r}^{j}\in\bar{R}$.

12 for  _$g\in\\{1,\dotsc,q\\}$_ do

13 Compute $v_{g}=\left|\\{j\in[1,t]:r^{j}=\bar{r}^{g}\\}\right|$ and set
$s^{g}\leftarrow v_{g}s^{g}$.

14 end for

15for  _$i\in\\{1,\dotsc,h\\}\backslash\mathcal{K}$_ do

16 Compute $\theta_{i}=|\text{mean}(s_{i})|$.

17 end for

18for  _$i\in\\{1,\dotsc,h\\}\backslash\mathcal{K}$_ do

19 Compute $\xi_{i}$, i.e., the position of $\theta_{i}$ in the vector of the
increasingly sorted values of $\theta$.

20 end for

Compute $\bar{\lambda}=\min\\{0,\lambda-\beta\\}$ and
$\bar{\eta}=\max\\{\bar{\lambda},m-\beta-\gamma-\bar{\lambda}\\}$ and solve
$\displaystyle\min_{\alpha,\eta,z}\quad$ $\displaystyle\eta$ s.t.
$\displaystyle\alpha^{\top}s_{i}\leq-1+z_{i}M,\quad
i\in[1,h]\setminus\mathcal{K},$ $\displaystyle\alpha^{\top}s_{i}\geq
1-(1-z_{i})M,\quad i\in[1,h]\setminus\mathcal{K},$
$\displaystyle\bar{\lambda}-\eta\leq\sum_{i\in[1,h]\setminus\mathcal{K}}w_{i}z_{i}\leq\bar{\lambda}+\eta,$
$\displaystyle\ell\leq\alpha_{j}\leq u,\quad j\in[1,q],$ $\displaystyle
0\leq\eta\leq\bar{\eta},$ $\displaystyle z_{i}\in\\{0,1\\},\quad
i\in[1,h]\setminus\mathcal{K}.$ with branching priorities $\xi$ to compute
$\alpha^{*},\eta^{*},z^{*}$.

Algorithm 1 p-C2RF: Preprocessing and Solving C2RF

##  5\. Numerical Results

In this section, we present and discuss our computational results that
demonstrate the impact of considering the total amount of points in each class
and of using the preprocessing techniques as well as the branching rule to
speed up the solution process.

We illustrate this on different test sets from the literature. The test sets
are discussed in Section 5.1, while the computational setup is described in
Section 5.2. The evaluation criteria are depicted in Section 5.3. Finally, the
numerical results are discussed in Section 5.4 and 5.5.

###  5.1. Test Sets

For the computational analysis of the proposed approaches, we consider the
subset of instances presented by [22] that are suitable for classification
problems and that have at most three classes and at least 5000 points.
Repeated instances are removed and instances with missing information are
reduced to the observations without missing information. If three classes are
given in an instance, we transform them into two classes such that the class
with label 1 represents the positive class and the other two classes represent
the negative class. This results in a final test set of 13 instances, as
listed in Table 1. To avoid numerical instabilities, all data sets are scaled
as follows. For each coordinate $j\in[1,p]$, we compute

$l_{j}=\min_{i\in[1,N]}\left\\{x^{i}_{j}\right\\},\quad
u_{j}=\max_{i\in[1,N]}\left\\{x^{i}_{j}\right\\},\quad
m_{j}=0.5\left(l_{j}+u_{j}\right)$

and shift each coordinate $j$ of all data points $x^{i}$ via
$\bar{x}^{i}_{j}=x^{i}_{j}-m_{j}$. Furthermore, if a coordinate $j$ of the re-
scaled points is still large, i.e., if $\tilde{l}_{j}=l_{j}-m_{j}<-10^{2}$ or
$\tilde{u}_{j}=u_{j}-m_{j}>10^{2}$ holds, it is re-scaled via

$\tilde{x}^{i}_{j}=(\bar{v}-\underline{v})\frac{\bar{x}^{i}_{j}-\tilde{l}_{j}}{\tilde{u}_{j}-\tilde{l}_{j}}+\bar{v}$

with $\bar{v}=10^{2}$ and $\underline{v}=-10^{2}$. The corresponding 4
instances that we re-scale are marked with an asterisk in Table 1.

Table 1. The entire test set with the number of points ($N$) and the dimension ($p$) ID | Instance | $N$ | $p$
---|---|---|---
1 | phoneme | 5349 | 5
2 | magic | 18 905 | 10
$3^{*}$ | adult | 48 790 | 14
$4^{*}$ | churn | 5000 | 20
$5^{*}$ | ring | 7400 | 20
6 | twonorm | 7400 | 20
7 | waveform_21 | 5000 | 21
8 | ann_thyroid | 7129 | 21
9 | agaricus_lepiota | 8124 | 22
10 | waveform_40 | 5000 | 40
11 | connect_4 | 67 557 | 42
12 | coil2000 | 8380 | 85
$13^{*}$ | clean2 | 6598 | 168

In our computational study, we focus on emphasizing the statistical importance
of cardinality constraints—mainly in the case of non-representative biased
samples. Biased samples are highly recurrent in non-probability surveys, which
are surveys with an inclusion process that is not tracked and, hence, the
inclusion probabilities are unknown. This means that correction methods such
as inverse inclusion probability weighting cannot be applicable. For a primer
on inverse inclusion probability weighting, we refer to [26] and the
references therein.

To reproduce such a scenario, we create 5 biased samples with
$1\text{\,}\mathrm{\char 37\relax}$ of the data being labeled for each
instance. Differently from a simple random sample, where each point has an
equal probability of being chosen as labeled data, in these biased samples the
labeled data is chosen with probability $85\text{\,}\mathrm{\char 37\relax}$
for belonging to the positive class. Moreover, we consider $t=20$ trees and
for each $j\in[1,t]$, the size of the training subset $A^{j}$ is
$20\text{\,}\mathrm{\char 37\relax}$ of the labeled data. For each training
subset we use the Decision Tree package [24] to generate $r^{j}$.

In addition, in Appendix A, we provide the results under simple random
sampling, which produces unbiased samples. In this scenario, the results of
the proposed methods are similar to the random forest. Hence, there is no
downside to using the proposed method in case of an unknown sampling process.

###  5.2. Computational Setup

For each one of the 65 instances described in Section 5.1, we compare the
following approaches.

  1. (a)

RF: Random Forest by majority vote.

  2. (b)

C2RF as given in Problem (2) with $\bar{\eta}$ as defined in (1) and $M$ from
Proposition 1.

  3. (c)

p-C2RF as described in Algorithm 1.

  4. (d)

only PP: Algorithm 1 without the branching rule described in Step 1.

  5. (e)

only BR: C2RF as given in Problem (2) with the branching rule as described in
Section 4 but without our problem-tailored preprocessing.

Our comparison has been implemented in Julia 1.8.5 and we use Gurobi 11.5 and
JuMP [11] to solve C2RF as well as the MILP in Algorithm 1. All computations
were executed on the high-performance cluster “Elwetritsch”, which is part of
the “Alliance of High-Performance Computing Rheinland-Pfalz” (AHRP). We use a
single Intel XEON SP 6126 core with $2.6\text{\,}\mathrm{GHz}$ and
$64\text{\,}\mathrm{GB}$ RAM as well as a time limit of
$7200\text{\,}\mathrm{s}$.

Based on our preliminary experiments, for C2RF and p-C2RF we set the bounds to
$\ell=1$ and $u=100$. Moreover, we set the MIPFocus parameter of Gurobi to 3.

###  5.3. Evaluation Criteria

The first evaluation criterion is the run time of the different methods. To
compare run times, we use empirical cumulative distribution functions (ECDFs).
Specifically, for $S$ being a set of solvers (or approaches as above) and for
$P$ being a set of problems, we denote by $t_{p,s}\geq 0$ the run time of the
approach $s\in S$ applied to the problem $p\in P$ in seconds. If
$t_{p,s}>7200$, we consider problem $p$ as not being solved by approach $s$.
With these notations, the performance profile of approach $s$ is the graph of
the function $\gamma_{s}:[0,\infty)\to[0,1]$ given by

$\gamma_{s}(\sigma)=\frac{1}{|P|}\lvert\left\\{p\in P\colon
t_{p,s}\leq\sigma\right\\}\rvert.$

Moreover, knowing the true label of all points, we categorize them into four
distinct categories: true positive (TP) or true negative (TN) if the point is
classified correctly in the positive or negative class, respectively, as well
as false positive (FP) if the point is misclassified in the positive class and
as false negative (FN) if the point is misclassified in the negative class.
Motivated by this, we compute two classification metrics, for which a higher
value indicates a better classification. The first one is accuracy (AC). It
measures the proportion of correctly classified points and is given by

$\text{AC}\mathrel{{\mathop{:}}{=}}\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}\in[0,1].$
(8)

The second metric is Matthews correlation coefficient (MCC). It measures the
correlation coefficient between the observed and predicted classifications and
is computed by

$\text{MCC}\mathrel{{\mathop{:}}{=}}\frac{\text{TP}\times\text{TN}-\text{FP}\times\text{FN}}{\sqrt{(\text{TP}+\text{FP})(\text{TP}+\text{FN})(\text{TN}+\text{FP})(\text{TN}+\text{FN})}}\in[-1,1].$
(9)

The main statistical question is the following: For a specific instance, does
using the cardinality constraint as additional information increase the
accuracy and the MCC? Since C2RF p-C2RF, only PP and only BR solve the same
problem, we only compare the difference of the accuracy and MCC according to

$\overline{\text{AC}}\mathrel{{\mathop{:}}{=}}\text{AC}_{\textsf{p-C${}^{2}$RF}}-\text{AC}_{\textsf{RF}},\quad\overline{\text{MCC}}\mathrel{{\mathop{:}}{=}}\text{MCC}_{\textsf{p-C${}^{2}$RF}}-\text{MCC}_{{\textsf{RF}}},$
(10)

where $\text{AC}_{\textsf{RF}}$ and $\text{AC}_{\textsf{p-C${}^{2}$RF}}$ are
computed as in (8), and $\text{MCC}_{\textsf{RF}}$ and
$\text{MCC}_{\textsf{p-C${}^{2}$RF}}$ as in (9).

###  5.4. Discussion of Run Times

Figure 1 shows the ECDFs for the measured run times. As expected, RF is the
fastest algorithm because it does not include any binary variable related to
the unlabeled points as the newly proposed models do. It can be seen that
p-C2RF significantly outperforms C2RF. C2RF solves only
$58\text{\,}\mathrm{\char 37\relax}$ of the instances within the time limit,
while p-C2RF solves $94\text{\,}\mathrm{\char 37\relax}$. This shows that the
preprocessing techniques and the branching rule significantly decrease the run
times. However, by comparing the two lines for “only PP” and “only BR”, we see
that most of the speed-up is obtained by the preprocessing techniques while
the branching rule only helps to improve the performance for a small amount of
instances. Since the branching rule is not harming and sometimes helps, we
decide to include it in what follows.

Figure 1. ECDFs for run times (in seconds) Table 2. Median of run times (in seconds) ID | RF | C2RF | p-C2RF | only PP | only BR
---|---|---|---|---|---
1 | 0.042 | 3859.72 | 2.939 | 3.000 | 1249.15
2 | 0.261 | — | 109.603 | 127.63 | —
3 | 0.513 | — | 2.865 | 2.865 | —
4 | 0.107 | 68.255 | 12.304 | 12.495 | 76.813
5 | 0.074 | — | 999.798 | 1018.56 | —
6 | 0.069 | — | — | — | —
7 | 0.172 | 1834.94 | 4.054 | 4.149 | —
8 | 0.058 | 3.791 | 0.168 | 0.191 | 3.585
9 | 0.087 | 899.891 | 1.599 | 1.568 | 1018.80
10 | 0.066 | 883.18 | 144.127 | 182.252 | —
11 | 1.129 | — | 204.337 | 147.684 | —
12 | 0.194 | 70.102 | 14.443 | 12.619 | 75.937
13 | 0.229 | 0.986 | 0.275 | 0.367 | 1.124

In Table 2 we present the median run times of the 5 biased samples for each
instance. A “—” means that the approach did not solve at least $3$ of the
samples of the instance within the time limit. Once can see that RF almost
always takes less than $1\text{\,}\mathrm{s}$ to solve the problem. When
comparing the two novel approaches and only the instances that C2RF solves at
least one sample, Table 2 shows that our techniques decrease the time
computation by $89\text{\,}\mathrm{\char 37\relax}$ on average.

###  5.5. Discussion of Accuracy and MCC

Observe that for both metrics $\overline{\text{AC}}$ and
$\overline{\text{MCC}}$, a value greater than zero indicates that p-C2RF had a
better result than RF. Besides that, the box in the boxplot depicts the range
of the medium $50\text{\,}\mathrm{\char 37\relax}$ of the values;
$25\text{\,}\mathrm{\char 37\relax}$ of the values are below and
$25\text{\,}\mathrm{\char 37\relax}$ are above the box. Figure 2 shows that
the $\overline{\text{AC}}$ values are greater than zero in
$75\text{\,}\mathrm{\char 37\relax}$ of the results (left plot). Hence, our
proposed method has better accuracy than the conventional random forest. It
can also be seen in Figure 2 that the $\overline{\text{MCC}}$ values are
greater than zero in most cases (right plot). Therefore, our method has a
better MCC than RF.

Figure 2. Comparison of $\overline{\text{AC}}$ (left) and $\overline{\text{MCC}}$ (right); see (10) Table 3. Median of AC and MCC (in percentage) ID | Accuracy | MCC
---|---|---
| RF | p-C2RF | RF | p-C2RF
1 | 62.16 | 72.51 | 69.28 | 69.20
2 | 65.14 | 75.03 | 70.35 | 73.13
3 | 76.28 | 77.78 | 55.66 | 67.16
4 | 61.98 | 79.64 | 55.17 | 57.09
5 | 50.80 | 60.20 | 54.20 | 61.45
6 | 58.90 | 66.93 | 65.76 | 67.05
7 | 75.88 | 78.59 | 78.55 | 76.47
8 | 98.58 | 98.74 | 84.03 | 84.73
9 | 81.30 | 87.59 | 83.95 | 87.57
10 | 61.35 | 71.13 | 71.14 | 70.00
11 | 24.79 | 56.69 | 52.19 | 55.98
12 | 89.10 | 88.72 | 54.61 | 53.57
13 | 99.43 | 100 | 98.89 | 100

When comparing each instance, Table 3 shows the median of AC and MCC of the 5
biased samples for RF and p-C2RF. It can be seen that, in the majority of
instances, our approach has a greater value of accuracy and MCC than RF.
Especially in terms of accuracy, we obtained a better median value in 12 of
the 13 instances. Regarding MCC, our approach has a better median value than
RF in 8 of the 13 instances. When RF has better MCC than p-C2RF, it is never
better than $2.5\text{\,}\mathrm{\char 37\relax}$.

Figure 2 and Table 3 show that using the cardinality constraint of each class
as additional information allows to correctly classify the points with higher
accuracy and better MCC than with the RF by majority vote.

##  6\. Conclusion

For several classification problems, it can be expensive to acquire labels for
the entire population of interest. Nevertheless, external sources can, in some
cases, offer additional information on how many points are in each class. For
the case of binary classification, we proposed a semi-supervised random forest
that can be modeled using a big-$M$-based MILP formulation. We also presented
problem-tailored preprocessing techniques and a branching rule to reduce the
computational cost of solving the MILP model.

Under the condition of simple random sampling, our proposed semi-supervised
method has very similar accuracy and MCC as a standard random forest. In many
applications, however, the available data come from non-probability samples.
In this case, the data collection mechanism is largely unknown and there is
the risk of obtaining biased samples. Our numerical results show that our
model has better accuracy and MCC than the conventional random forest even
with a small number of labeled points and biased samples.

## Acknowledgements

The authors thank the DFG for their support within RTG 2126 “Algorithmic
Optimization”.

## References

  * [1] Massih-Reza Amini and Patrick Gallinari  “Semi-Supervised Logistic Regression”  In _Proceedings of the 15th European Conference on Artificial Intelligence_ , ECAI’02  Lyon, France: IOS Press, 2002, pp. 390–394
  * [2] Gérard Biau and Erwan Scornet  “A random forest guided tour”  In _TEST: An Official Journal of the Spanish Society of Statistics and Operations Research_ 25.2, 2016, pp. 197–227  DOI: 0.1007/s11749-016-0481-7
  * [3] Leo Breiman  “Random Forests”  In _Machine Learning_ 45.1 Kluwer Academic Publishers, 2001, pp. 5–32  DOI: 10.1023/A:1010933404324
  * [4] Jan Pablo Burgard, Joscha Krause and Simon Schmaus  “Estimation of regional transition probabilities for spatial dynamic microsimulations from survey data lacking in regional detail”  In _Computational Statistics & Data Analysis_ 154, 2021, pp. 107048  DOI: 10.1016/j.csda.2020.107048
  * [5] Jan Pablo Burgard, Maria Eduarda Pinheiro and Martin Schmidt  “Mixed-Integer Linear Optimization for Semi-Supervised Optimal Classification Trees”, 2024  arXiv:2401.09848 [math.OC]
  * [6] Jan Pablo Burgard, Maria Eduarda Pinheiro and Martin Schmidt  “Mixed-integer quadratic optimization and iterative clustering techniques for semi-supervised support vector machines” To appear  In _TOP_ , 2024  DOI: 10.1007/s11750-024-00668-w
  * [7] Danilo Bzdok et al.  “Semi-Supervised Factored Logistic Regression for High-Dimensional Neuroimaging Data”  In _Advances in Neural Information Processing Systems_ 28 Cambridge, MA, USA: MIT Press, 2015, pp. 3348–3356
  * [8] Olivier Chapelle, Mingmin Chi and Alexander Zien  “A Continuation Method for Semi-Supervised SVMs”  In _Proceedings of the 23rd International Conference on Machine Learning_ , ICML ’06  New York, NY, USA: Association for Computing Machinery, 2006, pp. 185–192  DOI: 10.1145/1143844.1143868
  * [9] Adele Cutler, D. Richard Cutler and John R. Stevens  “Random Forests”  In _Ensemble Machine Learning: Methods and Applications_ New York, NY: Springer New York, 2012, pp. 157–175  DOI: 10.1007/978-1-4419-9326-7_5
  * [10] Nejdet Dogru and Abdulhamit Subasi  “Traffic accident detection using random forest classifier”  In _2018 15th learning and technology conference (L &T)_, 2018, pp. 40–45  IEEE  DOI: 10.1109/LT.2018.8368509
  * [11] Iain Dunning, Joey Huchette and Miles Lubin  “JuMP: A Modeling Language for Mathematical Optimization”  In _SIAM Review_ 59.2, 2017, pp. 295–320  DOI: 10.1137/15M1020575
  * [12] Vishan Kumar Gupta, Avdhesh Gupta, Dinesh Kumar and Anjali Sardana  “Prediction of COVID-19 confirmed, death, and cured cases in India using random forest model”  In _Big Data Mining and Analytics_ 4.2, 2021, pp. 116–123  DOI: 10.26599/BDMA.2020.9020016
  * [13] Trevor Hastie, Robert Tibshirani and Jerome Friedman  “The elements of statistical learning: data mining, inference and prediction”  Springer, 2009  DOI: 10.1007/978-0-387-84858-7
  * [14] Kyoungok Kim  “A hybrid classification algorithm by subspace partitioning through semi-supervised decision tree”  In _Pattern Recognition_ 60, 2016, pp. 157–163  DOI: 10.1016/j.patcog.2016.04.016
  * [15] Michelangelo Cćeci Dragi Kocev, Jurica Levatić and Sašo Džeroski  “Semi-supervised classification trees”  In _Journal of Intelligent Information Systems_ 49, 2017, pp. 461–486  DOI: 10.1007/s10844-017-0457-4
  * [16] Dong-Hyun Lee  “Pseudo-Label : The Simple and Efficient Semi-Supervised Learning Method for Deep Neural Networks”  In _ICML 2013 Workshop : Challenges in Representation Learning (WREPL)_ , 2013
  * [17] Christian Leistner, Amir Saffari, Jakob Santner and Horst Bischof  “Semi-Supervised Random Forests”  In _2009 IEEE 12th International Conference on Computer Vision_ , 2009, pp. 506–513  DOI: 10.1109/ICCV.2009.5459198
  * [18] Ming Li and Zhi-Hua Zhou  “Improve Computer-Aided Diagnosis With Machine Learning Techniques Using Undiagnosed Samples”  In _IEEE Transactions on Systems, Man, and Cybernetics - Part A: Systems and Humans_ 37.6, 2007, pp. 1088–1098  DOI: 10.1109/TSMCA.2007.904745
  * [19] Stefano Melacci and Mikhail Belkin  “Laplacian Support Vector Machines Trained in the Primal”  In _Journal of Machine Learning Research_ 12, 2009  DOI: 10.48550/ARXIV.0909.5422
  * [20] Thao N. N. Nguyen, Bharadwaj Veeravalli and Xuanyao Fong  “A Semi-Supervised Learning Method for Spiking Neural Networks Based on Pseudo-Labeling”  In _2023 International Joint Conference on Neural Networks (IJCNN)_ , 2023, pp. 1–7  DOI: 10.1109/IJCNN54540.2023.10191317
  * [21] Avital Oliver et al.  “Realistic Evaluation of Deep Semi-Supervised Learning Algorithms”  In _Advances in Neural Information Processing Systems_ 31 Curran Associates, Inc., 2018  DOI: 10.48550/arXiv.1804.09170
  * [22] Randal S. Olson et al.  “PMLB: a large benchmark suite for machine learning evaluation and comparison”  In _BioData Mining_ 10.36, 2017, pp. 1–13  DOI: 10.1186/s13040-017-0154-4
  * [23] Madhumita Pal and Smita Parija  “Prediction of Heart Diseases using Random Forest”  In _Journal of Physics: Conference Series_ 1817.1 IOP Publishing, 2021, pp. 012009  DOI: 10.1088/1742-6596/1817/1/012009
  * [24] Ben Sadeghi et al.  “DecisionTree.jl - A Julia implementation of the CART Decision Tree and Random Forest algorithms”  Zenodo, 2022  DOI: 10.5281/zenodo.7359268
  * [25] Jamie Shotton et al.  “Real-time human pose recognition in parts from single depth images”  In _CVPR 2011_ , 2011, pp. 1297–1304  DOI: 10.1109/CVPR.2011.5995316
  * [26] C. J. Skinner and D’arrigo  “Inverse probability weighting for clustered nonresponse”  In _Biometrika_ 98.4, 2011, pp. 953–966  DOI: 10.1093/biomet/asr058
  * [27] Shiyang Xuan et al.  “Random forest for credit card fraud detection”  In _2018 IEEE 15th International Conference on Networking, Sensing and Control (ICNSC)_ , 2018, pp. 1–6  DOI: 10.1109/ICNSC.2018.8361343
  * [28] Youqiang Zhang et al.  “Active Semi-Supervised Random Forest for Hyperspectral Image Classification”  In _Remote Sensing_ 11.24, 2019  DOI: 10.3390/rs11242974
  * [29] Arman Zharmagambetov and Miguel A. Carreira-Perpinan  “Semi-Supervised Learning with Decision Trees: Graph Laplacian Tree Alternating Optimization”  In _Advances in Neural Information Processing Systems_ 35, 2022, pp. 2392–2405
  * [30] Xiaojin Zhu and Andrew B. Goldberg  “Introduction to Semi-Supervised Learning”  In _Introduction to Semi-Supervised Learning_ , Synthesis Lectures on Artificial Intelligence and Machine Learning  Morgan & Claypool Publishers, 2009  DOI: 10.2200/S00196ED1V01Y200906AIM006

##  Appendix A Numerical Results for Simple random samples

In Section 5 we present a computational study on non-representative biased
samples. To complement our numerical results, we also present the results for
simple random sampling. For simple random sampling, each element in the data
set has the same probability $(n/N)$ to be included in the sample of labeled
data of size $n$. The instances are the same as described in Section 5.1. The
computational setup follows the description in Section 5.2. As before, the
used evaluation criteria are $\overline{\text{AC}}$ and
$\overline{\text{MCC}}$ as in (10).

Figure 3. Comparison of $\overline{\text{AC}}$ (left) and
$\overline{\text{MCC}}$ (right); see (10)

It can be seen in Figure 3 that $75\text{\,}\mathrm{\char 37\relax}$ of the
values of $\overline{\text{AC}}$ are between $-0.05$ and $0.05$ (left plot).
Figure 3 (right plot) also shows that $\overline{\text{MCC}}$ has a value
greater than 0 and lower than 0 in $50\text{\,}\mathrm{\char 37\relax}$ of the
cases.

Table 4 shows the median of AC and MCC for each instance for p-C2RF and RF. In
the majority of instances, our approach has a better or a very similar
accuracy and MCC compared to the conventional random forest. Especially in
terms of MCC, this is the case for all 13 instances. From Figure 3 and Table 4
we can conclude that the accuracy and MCC of our proposed method p-C2RF and
the standard random forest are very similar in the context of simple random
sampling. This is expected because the cardinality constraint aims to balance
the class distribution and since the sample is not biased, this constraint
does not introduce additional meaningful information to the problem.

Table 4. Median of AC and MCC (in percentage) ID | Accuracy | MCC
---|---|---
| RF | p-C2RF | RF | p-C2RF
1 | 76.68 | 76.32 | 71.73 | 71.34
2 | 78.28 | 78.14 | 75.82 | 75.86
3 | 81.32 | 80.85 | 70.79 | 73.70
4 | 85.86 | 76.57 | 50.0 | 51.63
5 | 71.81 | 67.35 | 72.61 | 67.35
6 | 84.32 | 85.09 | 85.32 | 85.10
7 | 76.75 | 77.10 | 71.97 | 74.10
8 | 97.68 | 98.61 | 50.0 | 84.43
9 | 88.15 | 87.17 | 88.17 | 87.15
10 | 78.57 | 79.84 | 75.13 | 77.23
11 | 75.36 | 67.71 | 50.0 | 55.89
12 | 93.23 | 87.01 | 50.0 | 52.62
13 | 97.64 | 100 | 95.37 | 100