# Quantum Induced Coherence Light Detection and Ranging

Gewei Qian1, Xingqi Xu1,∗, Shun-An Zhu1, Chenran Xu1, Fei Gao2, V. V.
Yakovlev3, Xu Liu1, Shi-Yao Zhu1,4 and Da-Wei Wang1,4,5,†

###### Abstract

Quantum illumination has been proposed and demonstrated to improve the signal-
to-noise ratio (SNR) in light detection and ranging (LiDAR). When relying on
coincidence detection, such a quantum LiDAR is limited by the response time of
the detector and suffers from jamming noise. Inspired by the Zou-Wang-Mandel
experiment, we design, construct and validate a quantum induced coherence
(QuIC) LiDAR which is inherently immune to ambient and jamming noises. In
traditional LiDAR the direct detection of the reflected probe photons suffers
from deteriorating SNR for increasing background noise. In QuIC LiDAR we
circumvent this obstacle by only detecting the entangled reference photons,
whose single-photon interference fringes are used to obtain the distance of
the object, while the reflected probe photons are used to erase path
information of the reference photons. In consequence, the noise accompanying
the reflected probe light has no effect on the detected signal. We demonstrate
such noise resilience with both LED and laser light to mimic the background
noise and jamming attack. The proposed method paves a new way of battling
noise in precise quantum electromagnetic sensing and ranging.

††preprint: APS/123-QED

Quantum properties of light such as entanglement and squeezing can be used to
enhance the sensitivity and signal-to-noise ratio (SNR) in light detection and
ranging (LiDAR) Giovannetti et al. (2001, 2002); Lloyd (2008); Maccone and Ren
(2020). Such a quantum LiDAR can be accomplished by using quantum illumination
(QI) of entangled photons combined with correlation detection Lloyd (2008);
Maccone and Ren (2020); Tan et al. (2008); Shapiro (2020); Nair and Gu (2020);
Lanzagorta (2011); Wilde et al. (2017); Sanz et al. (2017); Zhang et al.
(2013); Zhuang and Shapiro (2022), where one beam of entangled photons from
spontaneous parametric down conversion (SPDC) serves as a reference while the
other beam is sent to probe the presence of an object. The extra information
provided by the reference photons in correlation detection significantly
improves the SNR in detecting the object in noisy and lossy environment
Lopaeva et al. (2013); Pirandola et al. (2018); Zhang et al. (2015); Zhuang et
al. (2017) (see Fig. 1 (a) and (b)). The spatial and temporal correlation of
entangled photons can be used to extract the image Gregory et al. (2020) and
distance Zhao et al. (2022) of the object based on joint measurement, which
has been experimentally extended to the microwave regime Barzanjeh et al.
(2020); Chang et al. (2019); Barzanjeh et al. (2020). However, in the above
described approach the noise from the object still enters the detector, which
can be jammed by saturation attacks when the noise is much stronger than the
signal.

Quantum induced coherence, first demonstrated in the Zou-Wang-Mandel (ZWM)
experiment Zou et al. (1991) without employing coincidence detection, has been
used to image objects with undetected photons Lemos et al. (2014); Paterova et
al. (2020a); Kviatkovsky et al. (2020); Paterova et al. (2020b); Hochrainer et
al. (2022, 2017). In such a setup entangled photons are generated in the SPDC
of two nonlinear crystals. The information on the origin of the entangled
photon pairs plays a key role in the single-photon interference of the SPDC
light. The interference of signal photons depends on whether we in principle
know which crystal generates the photon pairs. By overlapping the idler modes
of the two crystals, which-way information of the signal photons is erased
Zajonc et al. (1991); Scully and Drühl (1982); Scully et al. (1991), the
single-photon interference of the signal photons can be observed. On the other
hand, if an object is placed between the two crystals to block the idler
photons, we can identify which crystal generates the entangled photon pairs,
such that the single-photon interference pattern of the signal photons
disappears. This mechanism has also been widely used in spectroscopy Kutas et
al. (2021); Kalashnikov et al. (2016); Paterova and Krivitsky (2020), optical
coherence tomography Paterova et al. (2018); Vanselow et al. (2020); Mukamel
et al. (2020); Gilaberte Basset et al. (2021) and holography Töpfer et al.
(2022). In these applications visible photons are usually used to trigger the
detector while their entangled infrared partners are used to interact with the
object, such that infrared imaging and spectroscopy are achieved with visible
light detectors.

Figure 1:  Physical principles for conventional, QI and QuIC LiDARs. (a)
Conventional-illumination-based LiDAR. (b) Quantum illumination LiDAR, where
the SNR is enhanced by joint measurement. (c) QuIC LiDAR, instead of the probe
light which is mixed with noise, a locally generated reference light, which
benefits from quantum induced coherence, is detected. (d) The setup consists
of three modules (enclosed in dashed lines), the pump laser, the local
scanning module and the detection module. A collimated 532 nm laser (red line)
is used to pump a PPLN crystal to generate SPDC entangled photon pairs whose
frequencies can be tuned by adjusting the temperature of the crystal. Here the
wavelengths of the reference (blue line) and probe (yellow line) lights are
centered at 893 and 1316 nm. The two light beams are separated by a dichroic
mirror DM2. The probe beam is sent to the object while the reference beam is
sent to a scanning mirror Mr. The reflected probe, reference and pump beams go
through the crystal a second time. We detect the reference beam for QuIC LiDAR
with a Si-based CMOS camera. The probe light is also detected with an InGaAs
detector to compare the performance of QuIC LiDAR and traditional OCT.

In this Letter, we design and implement a quantum induced coherence (QuIC)
LiDAR based on the ZWM interference. It has a substantial advantage with
respect to previous QI LiDAR proposals since QuIC LiDAR detects the image and
distance of objects from the visibility of single-photon interference fringes
rather than two-photon coincidence counting, significantly enhancing the
detection efficiency, shortening the collection time and improving the ranging
resolution (see Fig.1 (b) and (c)). More importantly, since we detect the
reference photons rather than the probe photons that are used to interact with
the object, such a scheme is intrinsically immune to background noise and
intentional jamming. We demonstrate the ranging and imaging ability of QuIC
LiDAR and show its robustness in LED and laser backgrouond noise. This proof-
of-principle demonstration obtains a 20 dB enhancement in noise resilience of
the QuIC LiDAR, in comparison with that using classical coherence. QuIC LiDAR
can be readily extended to mid-infrared and terahertz regime, where signal
detection and background noise are challenging problems.

Figure 2:  Position acquisition of QuIC LiDAR. (a) Real-time images captured
by the CMOS camera during the scanning of Mr. The two images carrying
interference patterns correspond to the two surfaces of the plate. The
exposure time of each image is 50 ms. The scanning step is 60 nm with a total
scanning length 3.5 mm. (b) Interference fringes along the white lines in (a).
The inset shows the interference near the position marked by the asterisk. (c)
Interference visibility obtained from STFT with an integration window 100
$\mu$m. The step in the horizontal axis is 1 $\mu$m. Inset: the Fourier
transform of the interference in the red frame in (b). The maximum peak value
in the spatial frequency range from 2.0 to 2.4 $\mu$m-1 is extracted as the
visibility. The two peaks in the main frame with the weights being marked by
red dashed lines indicate the positions of the two surfaces. (d) Images on the
two surfaces obtained from QuIC LiDAR.  Figure 3:  Noise robustness of QuIC
LiDAR. The images in the upper (lower) row are extracted from the reference
(probe) light cameras. (a) and (b) are the images without background noise.
(c) and (d) are the images with background LED light. (e), (f), (g) and (h)
are images with a jamming laser shining on the sample. In (e) and (f) the
direction of the laser satisfies the phase matching condition for stimulated
PDC in the PPLN crystal, while in (g) and (h) there is a phase mismatch.
Except for a phase matched jamming attack with a laser source, the QuIC LiDAR
images from the reference light are unaffected by the noises. In comparison,
the probe light images are prone to all the noises. The powers of the
reflected probe beam, the LED light and the laser before entering the probe
camera are 72 nW, 44 nW and 2.6 $\mu$W. The exposure time is 6 ms for (a),
(c), (e) and (g), 15ms for (b) and (d), and 20 $\mu$s for (f) and (h). Here
(f) and (h) are unprocessed images from the probe camera. The reference image
has lower resolution compared with the one from the probe beam, due to the
imperfection of the transverse phase matching in the crystal, which is
restricted by the crystal aperture and the pump beam waist size Fuenzalida et
al. (2022).

Three light beams are used to build the QuIC LiDAR. A pump laser beam pumps a
nonlinear crystal, which generates a probe (idler) beam and a reference
(signal) beam in SPDC. Both the pump and reference beams are kept as local
while the probe beam is sent out to detect the object. The probe beam
reflected from the object is collected to go through the nonlinear crystal a
second time together with the locally reflected pump and reference beams in a
Michalson configuration Chekhova and Ou (2016) of the ZWM experiment Zou et
al. (1991) (see Fig. 1(d)). We use a single-frequency 532 nm (10 kHz
linewidth) laser to pump a type-0 periodically poled lithium niobate (PPLN)
crystal (length 2 cm). The narrow linewidth (long coherence length) of the
pump light guarantees that the coherence in the two SPDC processes is
maintained. We observe interference when the light paths of the reference and
probe beams (from leaving to re-entering the nonlinear crystal) are equal
within the coherence length of the SPDC modes. The light intensity of the
reference beam at the detector (see Supplementary Materials),

$\displaystyle
I_{r}\propto\eta^{2}[1+\gamma(\tau)|r_{p}|\cos(\phi_{p}+\phi_{r}-\phi_{0})],$
(1)

where $\eta$ is the probability amplitude of generating an entangled photon
pair in a single pass through the nonlinear crystal, $r_{p}$ is the
reflectivity amplitude of the probe beam from the object, $\phi_{p}$,
$\phi_{r}$ and $\phi_{0}$ are the accumulated phases of the probe, reference
and pump beams in the round loops from leaving to returning to the PPLN
crystal. Here we assume $\gamma(\tau)=\exp(-\tau^{2}/2\sigma^{2})$ be the
overlapping function between the returned probe and reference light modes with
$\tau$ being their traveling time difference. The interference fringes are
significant only when $\tau$ is smaller than $1/\Delta\omega$, where
$\Delta\omega$ is the frequency bandwidth of the SPDC light. We obtain the
distance of the object by scanning the light path of the reference beam to
achieve maxima in the interference visibility,

$V=|r_{p}|\gamma(\tau).$ (2)

Compared to QI LiDAR which relies on coincidence counting with signals
proportional to $|r_{p}|^{2}$, the visibility in QuIC LiDAR is proportional to
$|r_{p}|$, i.e., linear with the reflectivity amplitude of the object Wang et
al. (1991). Such a linear scaling makes QuIC LiDAR highly advantageous for
detecting surfaces with low reflectivity. This sensitivity enhancement is the
result of the mixing between the returned probe light and the local reference
light before being detected, similar to the mixing between the signal and the
local oscillator in homodyne detection.

The performance of QuIC LiDAR is tested with two silica plates with hollowed
symbols (Fig 1 (d)). We scan the reference light path until observing
interference patterns on the CMOS camera, which indicates that an object
reflects the probe light from a distance the same as the one between Mr and
the PPLN crystal. We find two such positions corresponding to the two surfaces
of the silica plates (see Fig. 2 (a) and (b)). The envelopes of the
interference fringes during the scanning have a spatial width 0.4 mm, which is
determined by the coherence length of the entangled photons and sets a limit
on the ranging resolution of QuIC LiDAR. This is in the same spirit of
improving resolution through mode gating Ansari et al. (2021), where the
gating mode is provided by the reflected probe photons. In order to filter the
noise in determining the distance of the object, we perform short-time Fourier
transform (STFT) with a scanning window 100 $\mu$m and obtain the visibility
as a function of the scanning distance (Fig. 2 (c)). The weights of the two
peaks are regarded as the expected distances of the two surfaces, such that a
3D reconstruction of objects is achieved with a lateral resolution 412 $\mu$m
and a ranging accuracy 5.1 $\mu$m (see Fig. 2 (d) and Supplementary
Materials), which can be considered as a quantum version of the optical
coherence tomography (OCT) Paterova et al. (2018); Vallés et al. (2018);
Vanselow et al. (2020). The ranging accuracy is higher than that of the LiDAR
relying on the measurement of time of flight (including QI LiDAR), which is
limited by the jitter time (tens of picoseconds resulting in milimeter
uncertainty) of electronic devices. The visibility of the interference
patterns varies across the beam due to the finite lateral coherence of the
SPDC light.

While QI can significantly enhance the SNR when the noise is much lower than
the signal Lloyd (2008); Lopaeva et al. (2013), its performance in LiDAR
quickly deteriorates when the background noise increases to the same level of
the signal. It is also prone to jamming attack since all photons from the
object are detected (see Fig.1(b)). In comparison, QuIC LiDAR is superior in
noise resilience, originating from the essence of ZWM experiment, i.e., the
position of the object is obtained by detecting the local reference light
instead of the reflected probe light (see Fig.1(c)). The ambient noise that
has the same wavelength as the probe light, although can go into the local
module, is filtered by the dichroic mirror DM4 before the reference light
detector. Such noise neither has any effect on the interference pattern of the
reference light, since they cannot alter its which-way-information. On the
other hand, the noise that has the same polarization and wavelength as the
detected reference light is filtered by another dichiroic mirror, DM2, before
entering the local module.

We test such noise resilience of QuIC LiDAR by shining the sample from behind
with a LED light and a laser. In order to make a comparison, we also extract
the images from the probe light following the same procedure as that of the
reference light to mimick the performance of ordinary optical coherence
tomography. Without background noise, the two images extracted from the
reference and probe light cameras are similar except for a slight difference
in resolution Fuenzalida et al. (2022) (Fig. 3(a) and (b)). We use LED light
with a similar wavelength as that of the probe light to simulate background
noise. The spectrum power density of the LED is about 7 mW/m2/nm, much lower
than that of the sun (300 mW/m2/nm). As a result, the noise enters the probe
light camera and substantially obscures the image (Fig. 3(d)), while no
difference is detected by the reference light camera (Fig. 3(c)), which
demonstrates the robustness of QuIC LiDAR against the background noise.

Figure 4:  The SNR of QuIC LiDAR in LED (a) and laser (b) noises. In (a), we
compare the SNR in QuIC LiDAR (blue triangles) and in direct detection of the
probe light (red circles). For noise level below 23.4 dB, the LED noise is
weaker than the systematic noise, such that the SNR of the probe light remains
unchanged. For higher LED noise, the SNR of the probe light decreases, while
the one of the QuIC LiDAR remains unchanged. In (b), we simulate saturation
jamming attacks by sending a laser back to the crystal, which induces
stimulated PDC noise. The red circles and blue triangles are the SNRs of the
reference light in the areas satisfying and violating the phasing matching
condition (of the jamming laser), respectively. The powers of probe and noise
light are measured by an InGaAs power meter with an iris covering the region
of interest on the camera.

The only practical way to jam QuIC LiDAR is to send a laser that has the same
frequency and polarization as the probe light back to the local module. The
stimulated parametric down conversion (PDC) light can jam the detector of the
reference light (see the black spot in Fig. 3 (e)). However, in order to jam
the LiDAR through the PDC process, the pump and jamming lasers have to overlap
in space, be colinear, and the wavelength and propagating direction of the
jamming laser should satisfy the phase matching conditions for a given crystal
temperature. A slight deviation of the laser from the phase matching angle
substantially reduces the stimulated PDC noise. As a result, the jamming is
only efficient in specific directions covering a few pixels, leaving other
pixels undisturbed (see Fig. 3 (e)). Such a strict requirement on phase-
matching makes it difficult to jam QuIC LiDAR, and also allows us to counter
fight such an attack by adjusting the orientation or the temperature of the
PPLN crystal (Fig. 3 (g)). In comparison, the probe light camera can be easily
jammed by the laser (see Fig. 3 (f) and (h)).

To quantitatively characterize the noise robustness of QuIC LiDAR, we measure
the SNR, $\mathcal{F}_{s}/\mathcal{F}_{n}$, where $\mathcal{F}_{s}$ and
$\mathcal{F}_{n}$ are the signal peak and noise baseline in the Fourier
transform spectrum (see Fig.S4 in the Supplementary Materials), at different
noise levels defined in decibel, $10{\rm log}_{10}(P_{n}/P_{p})$, where
$P_{n}$ and $P_{p}$ are the powers of the noise and reflected probe light,
respectively (see Fig.4). We compare the SNR of the interference of the probe
and reference light to quantify the advantage of QuIC LiDAR in noise
resilience. For the LED noise, the SNR of the probe light starts dropping at a
noise level 23.4 dB until going below the detection line at 36 dB, while the
SNR of the reference light is unaffected (see Fig. 4(a)). While only 43.3 dB
of noise power can be imposed with our LED, QuIC LiDAR can resist much
stronger noise. For saturation attacks from a laser source, the stimulated PDC
light at the phase matching direction can reduce the SNR of the reference
light to 1 with a laser noise 43.7 dB. However, in pixels where the phase
matching condition cannot be satisfied, the SNR remains unaffected, as shown
by the blue triangles in Fig. 4(b).

In conclusion, we demonstrate that quantum induced coherence can be used to
build a quantum LiDAR that has superior noise resilience compared with QI
LiDAR which relies on coincidence counting. Although a similar setup has been
used in imaging, spectroscopy and optical coherence tomography, featured by
visible light detection with infrared light probing, its ability in ranging
and battling noise remains largely unnoticed. With the freedom in selecting
wavelengths in SPDC Hochrainer et al. (2022), we can build mid-IR and
terahertz QuIC LiDAR by using proper crystals (e.g., AgGaS2 for mid-infrared
up to 13um) Kumar et al. (2021); Kutas et al. (2021); Mukai et al. (2022);
Paterova et al. (2022). In such wavelengths the LiDAR detection is better
concealed and difficult to jam because efficient detectors and lasers are
rare. QuIC LiDAR can be used to meet the challenges of background noise and
efficient detection in those wavelengths. For applications in microwave regime
where background noise brings a huge challenge, we can integrate the
microwave-to-optics conversion Chang et al. (2019); Andrews et al. (2014); Tu
et al. (2022) to the current design or use microwave-optical entangled light
sources Rueda et al. (2019). The weak intensity of SPDC light makes it
difficult to detect objects with low reflection. A solution to this problem is
to use an SU(1,1) interferometer by adding a seed to our current setup to
improve the signal-to-noise ratio Heuer et al. (2015). The QuIC LiDAR
described here can be integrated with on-chip entangled photon sources and
delay lines Ono et al. (2019) to make it small and compact to be used in
modern transportation vehicles.

This work was supported by the National Natural Science Foundation of China
(Grants No. 11934011), Zhejiang Province Key Research and Development Program
(Grant No. 2020C01019), the Strategic Priority Research Program of Chinese
Academy of Sciences (Grant No. XDB28000000), and the Fundamental Research
Funds for the Central Universities. VVY acknowledges partial support from the
NSF (CMMI-1826078), AFOSR (FA9550-20-1-0366), and NIH (1R01GM127696,
1R21GM142107, 1R21CA269099).

<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
## References

  * Giovannetti et al. (2001) V. Giovannetti, S. Lloyd, and L. Maccone, Nature 412, 417 (2001).
  * Giovannetti et al. (2002) V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. A 65, 022309 (2002).
  * Lloyd (2008) S. Lloyd, Science 321, 1463 (2008).
  * Maccone and Ren (2020) L. Maccone and C. Ren, Phys. Rev. Lett. 124, 200503 (2020).
  * Tan et al. (2008) S.-H. Tan, B. I. Erkmen, V. Giovannetti, S. Guha, S. Lloyd, L. Maccone, S. Pirandola, and J. H. Shapiro, Phys. Rev. Lett. 101, 253601 (2008).
  * Shapiro (2020) J. H. Shapiro, IEEE Aerospace and Electronic Systems Magazine 35, 8 (2020).
  * Nair and Gu (2020) R. Nair and M. Gu, Optica 7, 771 (2020).
  * Lanzagorta (2011) M. Lanzagorta, Synth. Lect. Quantum Comput. 3, 1 (2011).
  * Wilde et al. (2017) M. M. Wilde, M. Tomamichel, S. Lloyd, and M. Berta, Phys. Rev. Lett. 119, 120501 (2017).
  * Sanz et al. (2017) M. Sanz, U. Las Heras, J. J. García-Ripoll, E. Solano, and R. Di Candia, Phys. Rev. Lett. 118, 070803 (2017).
  * Zhang et al. (2013) Z. Zhang, M. Tengner, T. Zhong, F. N. Wong, and J. H. Shapiro, Phys. Rev. Lett. 111, 010501 (2013).
  * Zhuang and Shapiro (2022) Q. Zhuang and J. H. Shapiro, Phys. Rev. Lett. 128, 010501 (2022).
  * Lopaeva et al. (2013) E. Lopaeva, I. R. Berchera, I. P. Degiovanni, S. Olivares, G. Brida, and M. Genovese, Phys. Rev. Lett. 110, 153603 (2013).
  * Pirandola et al. (2018) S. Pirandola, B. R. Bardhan, T. Gehring, C. Weedbrook, and S. Lloyd, Nat. Photonics 12, 724 (2018).
  * Zhang et al. (2015) Z. Zhang, S. Mouradian, F. N. Wong, and J. H. Shapiro, Phys. Rev. Lett. 114, 110506 (2015).
  * Zhuang et al. (2017) Q. Zhuang, Z. Zhang, and J. H. Shapiro, Phys. Rev. Lett. 118, 040801 (2017).
  * Gregory et al. (2020) T. Gregory, P.-A. Moreau, E. Toninelli, and M. J. Padgett, Sci. Adv. 6, eaay2652 (2020).
  * Zhao et al. (2022) J. Zhao, A. Lyons, A. C. Ulku, H. Defienne, D. Faccio, and E. Charbon, Opt. Express 30, 3675 (2022).
  * Barzanjeh et al. (2020) S. Barzanjeh, S. Pirandola, D. Vitali, and J. M. Fink, Sci. Adv. 6, eabb0451 (2020).
  * Chang et al. (2019) C. S. Chang, A. Vadiraj, J. Bourassa, B. Balaji, and C. Wilson, Appl. Phys. Lett. 114, 112601 (2019).
  * Zou et al. (1991) X. Zou, L. J. Wang, and L. Mandel, Phys. Rev. Lett. 67, 318 (1991).
  * Lemos et al. (2014) G. B. Lemos, V. Borish, G. D. Cole, S. Ramelow, R. Lapkiewicz, and A. Zeilinger, Nature 512, 409 (2014).
  * Paterova et al. (2020a) A. V. Paterova, S. M. Maniam, H. Yang, G. Grenci, and L. A. Krivitsky, Sci. Adv. 6, eabd0460 (2020a).
  * Kviatkovsky et al. (2020) I. Kviatkovsky, H. M. Chrzanowski, E. G. Avery, H. Bartolomaeus, and S. Ramelow, Sci. Adv. 6, eabd0264 (2020).
  * Paterova et al. (2020b) A. V. Paterova, H. Yang, Z. S. Toa, and L. A. Krivitsky, Appl. Phys. Lett. 117, 054004 (2020b).
  * Hochrainer et al. (2022) A. Hochrainer, M. Lahiri, M. Erhard, M. Krenn, and A. Zeilinger, Rev. Mod. Phys. 94, 025007 (2022).
  * Hochrainer et al. (2017) A. Hochrainer, M. Lahiri, R. Lapkiewicz, G. B. Lemos, and A. Zeilinger, Optica 4, 341 (2017).
  * Zajonc et al. (1991) A. Zajonc, L. Wang, X. Zou, and L. Mandel, Nature 353, 507 (1991).
  * Scully and Drühl (1982) M. O. Scully and K. Drühl, Phys. Rev. A 25, 2208 (1982).
  * Scully et al. (1991) M. O. Scully, B.-G. Englert, and H. Walther, Nature 351, 111 (1991).
  * Kutas et al. (2021) M. Kutas, B. Haase, J. Klier, D. Molter, and G. von Freymann, Optica 8, 438 (2021).
  * Kalashnikov et al. (2016) D. A. Kalashnikov, A. V. Paterova, S. P. Kulik, and L. A. Krivitsky, Nat. Photonics 10, 98 (2016).
  * Paterova and Krivitsky (2020) A. V. Paterova and L. A. Krivitsky, Light Sci. Appl. 9, 1 (2020).
  * Paterova et al. (2018) A. V. Paterova, H. Yang, C. An, D. A. Kalashnikov, and L. A. Krivitsky, Quantum Sci. Technol. 3, 025008 (2018).
  * Vanselow et al. (2020) A. Vanselow, P. Kaufmann, I. Zorin, B. Heise, H. M. Chrzanowski, and S. Ramelow, Optica 7, 1729 (2020).
  * Mukamel et al. (2020) S. Mukamel, M. Freyberger, W. Schleich, M. Bellini, A. Zavatta, G. Leuchs, C. Silberhorn, R. W. Boyd, L. L. Sánchez-Soto, A. Stefanov, et al., J. Phys. B 53, 072002 (2020).
  * Gilaberte Basset et al. (2021) M. Gilaberte Basset, A. Hochrainer, S. Töpfer, F. Riexinger, P. Bickert, J. R. León-Torres, F. Steinlechner, and M. Gräfe, Laser Photonics Rev. 15, 2000327 (2021).
  * Töpfer et al. (2022) S. Töpfer, M. Gilaberte Basset, J. Fuenzalida, F. Steinlechner, J. P. Torres, and M. Gräfe, Sci. Adv. 8, eabl4301 (2022).
  * Fuenzalida et al. (2022) J. Fuenzalida, A. Hochrainer, G. B. Lemos, E. A. Ortega, R. Lapkiewicz, M. Lahiri, and A. Zeilinger, Quantum 6, 646 (2022).
  * Chekhova and Ou (2016) M. Chekhova and Z. Ou, Adv. Opt. Photonics 8, 104 (2016).
  * Wang et al. (1991) L. Wang, X. Zou, and L. Mandel, Phys. Rev. A 44, 4614 (1991).
  * Ansari et al. (2021) V. Ansari, B. Brecht, J. Gil-Lopez, J. M. Donohue, J. Řeháček, Z. c. v. Hradil, L. L. Sánchez-Soto, and C. Silberhorn, PRX Quantum 2, 010301 (2021).
  * Vallés et al. (2018) A. Vallés, G. Jiménez, L. J. Salazar-Serrano, and J. P. Torres, Phys. Rev. A 97, 023824 (2018).
  * Kumar et al. (2021) M. Kumar, P. Kumar, A. Vega, M. A. Weissflog, T. Pertsch, and F. Setzpfandt, Appl. Phys. Lett. 119, 244001 (2021).
  * Mukai et al. (2022) Y. Mukai, R. Okamoto, and S. Takeuchi, Opt. Express 30, 22624 (2022).
  * Paterova et al. (2022) A. V. Paterova, Z. S. Toa, H. Yang, and L. A. Krivitsky, ACS Photonics 9, 2151 (2022).
  * Andrews et al. (2014) R. W. Andrews, R. W. Peterson, T. P. Purdy, K. Cicak, R. W. Simmonds, C. A. Regal, and K. W. Lehnert, Nat. Phys. 10, 321 (2014).
  * Tu et al. (2022) H. Tu, K. Liao, Z. Zhang, X. Liu, S. Zheng, S. Yang, X. Zhang, H. Yan, and S. Zhu, Nat. Photonics 16, 291 (2022).
  * Rueda et al. (2019) A. Rueda, W. Hease, S. Barzanjeh, and J. M. Fink, npj Quantum Information 5, 108 (2019).
  * Heuer et al. (2015) A. Heuer, R. Menzel, and P. Milonni, Phys. Rev. A 92, 033834 (2015).
  * Ono et al. (2019) T. Ono, G. F. Sinclair, D. Bonneau, M. G. Thompson, J. C. Matthews, and J. G. Rarity, Opt. Lett. 44, 1277 (2019).


# Sparse Infinite Random Feature Latent Variable Modeling

Michael Minyi Zhang
Department of Statistics and Actuarial Science
University of Hong Kong
Hong Kong SAR
<EMAIL_ADDRESS>
###### Abstract

We propose a non-linear, Bayesian non-parametric latent variable model where
the latent space is assumed to be sparse and infinite dimensional a priori
using an Indian buffet process prior. A posteriori, the number of instantiated
dimensions in the latent space is guaranteed to be finite. The purpose of
placing the Indian buffet process on the latent variables is to: 1.)
Automatically and probabilistically select the number of latent dimensions.
2.) Impose sparsity in the latent space, where the Indian buffet process will
select which elements are exactly zero. Our proposed model allows for sparse,
non-linear latent variable modeling where the number of latent dimensions is
selected automatically. Inference is made tractable using the random Fourier
approximation and we can easily implement posterior inference through Markov
chain Monte Carlo sampling. This approach is amenable to many observation
models beyond the Gaussian setting. We demonstrate the utility of our method
on a variety of synthetic, biological and text datasets and show that we can
obtain superior test set performance compared to previous latent variable
models.

##  1 Introduction

Latent variable modeling (LVM) is a popular statistical method for
representing high-dimensional observations with a lower-dimensional structure.
Popular examples of continuous latent variable modeling include methods like
factor analysis, where the observed data are assumed to decompose into a
linear combination of underlying factors, which are lower-dimensional random
variables which partially explain the variability of the observed data, and
the factor loadings, which represents the extent to which an observation
expresses the latent factor. LVMs can be interpreted as matrix factorization
techniques, where the observed data can be modeled as a linear combination of
two low-rank matrices. However, these linear models are limited in their
capacity to capture the underlying structure in the data.

Non-linear extensions of latent variable models have proven to be quite
powerful for modeling high-dimensional data. The Gaussian process latent
variable model (GPVLM) is one notable example of a Bayesian, non-linear,
latent variable model where the functional relationship between the latent
space and the observed data is assumed to be Gaussian process distributed.
However, in the Bayesian setting, inference is quite challenging as
computation is generally tractable under Gaussian observation models. One
method to make posterior inference tractable outside the Gaussian setting is
to approximate the kernel function in the Gaussian process with the random
Fourier feature approximation. Using the random projections from this
approximation, we can compute simple Gibbs sampling updates in a Markov chain
Monte Carlo (MCMC) sampler for GPLVM-based models outside of Gaussian
observation models.

Moreover, selecting the number of latent dimensions to use in these models is
an important variable but generally unknown a priori. Typical strategies for
selecting the latent dimensionality include cross validation, or selecting an
large number of latent dimensions and shrinking the effect of the excess
dimensions towards zero. These techniques are computationally demanding and
ignore the underlying uncertainty of the latent dimensions.

Additionally, the latent space can be difficult to interpret and may contain
spurious latent dimensions that do not offer any additional information about
the data. It can also be difficult to disentangle the what contributions a
latent dimension might provide to modeling some underlying attribute in the
observed data. Thus, enforcing sparsity in the latent variables results in a
concise and interpretable latent space that furthermore can produce a
regularization effect for superior predictive performance compared to non-
sparse methods.

In this paper, we present a novel Bayesian method for sparse, non-linear,
latent variable modeling. By combining a random Fourier features-based latent
variable model with a beta-Bernoulli prior on the latent space, we can
simultaneously automatically explore the number of latent dimensions, learn a
sparse representation of the latent variables, and tractability perform
posterior inference for Gaussian and non-Gaussian observations. We will
proceed with an introduction to previous work in latent variable modeling in
Section 2. Then we will formally introduce our model and the inference
procedure in Section 3. Next, we will examine the performance of our model on
a wide variety of data sets in Section 4. Lastly, we conclude our paper with a
discussion of possible future work in Section 5.

##  2 Background

###  2.1 Latent Variable Models

The basic, linear latent variable model is expressed as a low-rank matrix
decomposition: $\mathbf{Y}=\mathbf{X}\mathbf{W}$. We denote the observed data
with $\mathbf{Y}=[\mathbf{y}_{1},\ldots,\mathbf{y}_{N}]^{T}$, which is an
$N\times J$ matrix, for $N$ observations and $J$ observed dimensions. The
latent variables are denoted as
$\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{N}]^{T}$ which is an $N\times
D$ matrix, for $D\ll J$ latent dimensions. We denote
$\mathbf{W}=[\mathbf{w}_{1},\ldots\mathbf{w}_{J}]$ to be a $D\times J$
projection matrix. If we set $\mathbf{X}$ to be the product of the left
singular vectors and the diagonal matrix of singular values,
$\mathbf{U}\boldsymbol{\Sigma}$ and $\mathbf{W}$ to be the right singular
vectors corresponding to the $D$ largest eigenvalues of
$\mathbf{Y}\mathbf{Y}^{T}$ , then we obtain the optimal solution to principal
components analysis (Pearson,, 1901).

We can place probabilistic assumptions on the data generating process for the
linear LVM to obtain a Bayesian formulation of the basic factor analysis model
(Press and Shigemasu,, 1989):

$\mathbf{y}_{i}\sim\mathcal{N}_{J}(\mathbf{x}_{i}\mathbf{W},\sigma^{2}\mathbf{I}),\;\mathbf{x}_{i}\sim\mathcal{N}_{D}(0,\mathbf{I}),\;\mathbf{w}_{j}\sim\mathcal{N}_{D}(\boldsymbol{\mu}_{\mathbf{w}},\boldsymbol{\Sigma}_{\mathbf{w}}).$
(1)

In this context, $\mathbf{X}$ are the latent factors and $\mathbf{W}$ are the
factor loadings matrix. By marginalizing out the latent variables,
$\mathbf{X}$, and estimating $\mathbf{W}$ via maximum likelihood, we obtain
the probabilistic variant of principal components analysis (Tipping and
Bishop,, 1999). Thus, we can see the connection between factor analysis,
dimensionality reduction and low-rank matrix factorization and we can unify
these methods as latent variable modeling (Engelhardt and Stephens,, 2010).
LVMs are commonly used in applied settings, where we observe different
recordings corresponding to a single datum, particularly in scientific
applications where we are interested in measuring traits that are not directly
unobservable. For example, Conti et al., (2014) analyze the 1970 British
Cohort Study data set using a Bayesian factor analysis model and discover the
contribution of latent, unobservable factors like depression and anxiety play
in a person’s life outcomes.

Outside the Gaussian setting, we can model positive integer-valued data sets
using LVMs as well. Non-negative matrix factorization can model such data sets
by constraining $\mathbf{X}$ and $\mathbf{W}$ in the matrix factorization
setting to be positive valued (Lee and Seung,, 1999). In the probabilistic
setting, we can assume the data in the generating process is Poisson
distributed using Poisson factor analysis (PFA, Zhou et al., 2012a, ):

$\displaystyle\begin{split}y_{ij}&\sim\mbox{Poisson}(\mathbf{x}_{i}\mathbf{w}_{j}^{T}),\;x_{id}\sim\mbox{Gamma}\left(r_{d},\frac{1-p_{d}}{p_{d}}\right),\;\mathbf{w}_{j}\sim\mbox{Dirichlet}(w_{0},\ldots,w_{0})\\\
r_{d}&\sim\mbox{Gamma}(a_{0},b_{0}),\;p_{d}\sim\mbox{Beta}(c_{0}\epsilon_{0},c_{0}(1-\epsilon_{0})).\end{split}$
(2)

We can use model admixtures in using latent Dirichlet allocation (LDA), where
we can discover the underlying topical content of a corpus of documents (Blei
et al.,, 2003; Griffiths and Steyvers,, 2004):

$\displaystyle\begin{split}\mathbf{y}_{i}\sim\mbox{Mult}(n_{i};\mathbf{x}_{i}\mathbf{w}_{1}^{T},\ldots,\mathbf{x}_{i}\mathbf{w}_{J}^{T}),\;\mathbf{x}_{i}\sim\mbox{Dirichlet}(x_{0},\ldots,x_{0}),\;\mathbf{w}_{j}\sim\mbox{Dirichlet}(w_{0},\ldots,w_{0}).\end{split}$
(3)

where $N_{i}$ is the total number of words in document $i$:
$n_{i}=\sum_{j=1}^{J}y_{ij}$. This basic admixture model can decompose each
document into a probability distribution over topics, which represents the
semantical content of the document, and model a topic as a probability
distribution over the vocabulary.

###  2.2 Sparse Latent Variable Models

The basic latent variable model is a powerful tool for discovering the
underlying structure in a data set. However, often times we may assume that an
observation is associated with only a small subset of latent factors. The
sparsity constraint in a latent variable model acts as a regularizer which
enforces a concise latent representation of the data and can enhance the
interpretability of the latent factors (Olshausen and Field, 1996a, ;
Olshausen and Field, 1996b, ). The latent variables in the dense setting are
unidentifiable up to an arbitrary rotation. This poses a computational problem
as posterior inference algorithms for the latent variables can easily be
trapped in poor local optima. Sparsity can partially solve this indeterminacy
problem in latent variable models because the imposition of sparsity will
eliminate many equivalent representations of the latent variables under
rotational ambiguity.

A popular sparsity-inducing prior in Bayesian methods is the “spike-and-slab”
prior that selects which latent variables are exactly zero and which latent
variables are activated (Mitchell and Beauchamp,, 1988; West,, 2003):

$x_{id}\sim(1-\pi_{d})\cdot\delta_{0}+\pi_{d}\cdot\mathcal{N}(0,1),$ (4)

where the parameter $\pi_{d}\in(0,1)$ controls the sparsity level of the
latent dimension $d$.

###  2.3 Beta-Bernoulli and Indian Buffet Processes

One drawback in the aforementioned latent variable models is that the number
of latent dimensions, $D$ is assumed to be known and fixed a priori. In
actuality, this quantity is rarely known and typical techniques for selecting
the latent dimensionality like cross validation or reversible jump MCMC are
extraordinarily expensive to implement in practice. A more elegant approach
for selecting the latent dimensionality is to use Bayesian non-parametric
(BNP) methods. The BNP perspective for this problem is to assume the number of
latent dimensions is infinite a priori. This assumption is realistic for
latent variable models as we assume the number of latent dimensions grows to
infinity as we observe an infinite number of observations.

One BNP variant of spike-and-slab variable selection is using the beta-
Bernoulli process, which is a distribution over infinite dimensional binary
matrices. Suppose $\mathbf{Z}=[\mathbf{z}_{1},\ldots\mathbf{z}_{N}]^{T}$ is an
$N\times D$ binary-valued matrix and let the elements of $\mathbf{Z}$,
$z_{id}$, be distributed as:

$z_{id}|\pi_{d}\sim\mbox{Bernoulli}(\pi_{d}),\;\pi_{d}\sim\mbox{Beta}(\alpha/D,1).$
(5)

By letting $D\rightarrow\infty$ we obtain the beta-Bernoulli process and by
integrating out $\pi_{d}$ we obtain the Indian buffet process (IBP) which
results in the following stochastic process for generating binary matrices
(Thibaux and Jordan,, 2007; Griffiths and Ghahramani,, 2011): The first
observation occupies a $\mbox{Poisson}(\alpha)$ number of latent dimensions.
Each subsequent observation $i=2,\ldots,N$ selects the previously instantiated
latent dimensions with probability proportional to its popularity and then
samples $\mbox{Poisson}(\alpha/i)$ new dimensions. $\alpha$ has plays a
similar role to $\eta$ in the Dirichlet process as the concentration parameter
for a beta-Bernoulli process, where this parameter controls the number of new
features to sample.

Knowles and Ghahramani, (2011) uses the Indian buffet process for sparse
latent factor modeling:

$\displaystyle\begin{split}\mathbf{Z}\sim\mathcal{IBP}(\alpha),\;\mathbf{x}_{i}\sim\mathcal{N}_{D}(0,\mathbf{I}),\mathbf{w}_{j}\sim\mathcal{N}_{D}(\boldsymbol{\mu}_{\mathbf{w}},\boldsymbol{\Sigma}_{\mathbf{w}})\\\
\alpha\sim\mbox{Gamma}(\alpha_{0},\beta_{0}),\;\mathbf{y}_{i}\sim\mathcal{N}_{J}((\mathbf{x}_{i}\odot\mathbf{z}_{i})\mathbf{W},\sigma^{2}\mathbf{I}).\end{split}$
(6)

This sparse factor model is a Bayesian non-parametric variant of the factor
analysis model in Equation 1 with an infinite dimensional spike-and-slab prior
to control the sparsity of the latent space.

###  2.4 Non-linear Latent Variable Models and Gaussian Processes

The linear latent variable model is limited in its expressivity and cannot
capture non-linear relationships between the latent variables and the observed
data. One common prior over real-valued functions to model non-linear
relationships is the Gaussian process (Williams and Rasmussen,, 2006):
$f\sim\mathcal{GP}(\mu(\cdot),\mathbf{K}(\cdot,\cdot))$, given a mean function
$\mu(\cdot)$ and a covariance function $\mathbf{K}(\cdot,\cdot)$, which has
the property that a GP-distributed function evaluated at a finite set of
points is multivariate normal:
$f(\mathbf{X})\sim\mathcal{N}(\mu(\mathbf{X}),\mathbf{K}(\mathbf{x},\mathbf{x}^{\prime}))$.

Gaussian processes are have been widely used in non-linear latent variable
modeling (Lawrence,, 2005). The basic GPLVM model is

$\mathbf{y}_{j}\sim\mathcal{N}_{N}(f(\mathbf{X}),\sigma^{2}_{j}\mathbf{I}),\;f(\mathbf{X})\sim\mathcal{N}(\mu(\mathbf{X}),\mathbf{K}(\mathbf{x},\mathbf{x}^{\prime})),\;\mathbf{x}_{i}\sim\mathcal{N}_{D}(0,\mathbf{I}).$
(7)

Due to the conjugate relationship between the likelihood and prior, we can
integrate out the function $f$ in closed form–this property in crucial for
computation in GPLVMs:

$\mathbf{y}_{j}\sim\mathcal{N}_{N}(0,\mathbf{K}(\mathbf{x},\mathbf{x}^{\prime})+\sigma^{2}_{j}\mathbf{I}),\;\mathbf{x}_{i}\sim\mathcal{N}_{D}(0,\mathbf{I}).$
(8)

The GPLVM has proven to be a very popular method for non-linear latent
variable modeling. However, posterior inference for the GPLVM is
computationally tractable when the likelihood is Gaussian. For non-Gaussian
likelihoods, posterior inference becomes difficult because the mapping
function $f$ cannot be marginalized out in closed form. Without
marginalization, we cannot evaluate closed form expressions of the posterior
gradient with respect to the latent variables and therefore cannot use
sampling algorithms like Hamiltonian Monte Carlo (Duane et al.,, 1987) or
optimization techniques like evaluating the maximum a posteriori (MAP)
estimate.

###  2.5 Random Fourier Features and Kernel Learning

To achieve computational tractability in general settings, we can use random
Fourier features (RFF) to approximate the kernel function in GPs (Rahimi and
Recht,, 2008). Although RFFs are typically used to reduce the computational
complexity of kernel methods from $\mathcal{O}(N^{3})$ to
$\mathcal{O}(NM^{2})$, we can take advantage of the random feature mapping
induced by the RFF approximated GP-model that makes computation tractable for
non-Gaussian likelihoods.

According to Mercer’s theorem, we can equivalently represent a positive
definite kernel function as an inner product of a feature mapping (Mercer,,
1909):

$\mathbf{K}(\mathbf{x},\mathbf{x}^{\prime})=\left\langle\varphi(\mathbf{x}),\varphi(\mathbf{x}^{\prime})\right\rangle,\quad\mathbf{x},\mathbf{x}^{\prime}\in\mathbb{R}^{D}.$
(9)

Bochner’s theorem states that any continuous, stationary, kernel function,
$\mathbf{K}(\mathbf{x},\mathbf{x}^{\prime})=\mathbf{K}(\mathbf{x}-\mathbf{x}^{\prime})$,
is positive definite if and only if $\mathbf{K}(\cdot,\cdot)$ is the Fourier
transform of a non-negative measure $p(\mathbf{w})$, which is guaranteed to be
a density with a properly scaled kernel,
$\mathbf{K}(\mathbf{x},\mathbf{x})=\mathbf{K}(0)=1$ (Bochner,, 1959). Let
$\varphi(\mathbf{x})=\exp(i\mathbf{w}^{\top}(\mathbf{x}-\mathbf{x}^{\prime}))$,
then:

$\mathbf{K}(\mathbf{x}-\mathbf{x}^{\prime})=\int_{\mathbb{R}^{D}}p(\mathbf{w})\exp(i\mathbf{w}^{\top}(\mathbf{x}-\mathbf{x}^{\prime}))\text{d}\mathbf{w}=\mathbb{E}_{p(\mathbf{w})}[\varphi(\mathbf{x})\varphi(\mathbf{x}^{\prime})^{*}].$
(10)

Using these two theorems, we can obtain an unbiased Monte Carlo estimate of
the kernel function:

$\mathbf{K}(\mathbf{x}-\mathbf{x}^{\prime})\approx\frac{1}{M}\sum_{m=1}^{M}\varphi(\mathbf{x})\varphi(\mathbf{x}^{\prime})^{*},$
(11)

for $M$ Monte Carlo samples. If we have a real-valued kernel, then we can drop
the imaginary portion of $\varphi(\mathbf{x})$, so that

$\displaystyle\varphi(\mathbf{x})=\frac{1}{\sqrt{M}}\left[\sin(\mathbf{w}_{1}^{\top}\mathbf{x}),\;\cos(\mathbf{w}_{1}^{\top}\mathbf{x})\cdots\sin(\mathbf{w}_{M}^{\top}\mathbf{x}),\cos(\mathbf{w}_{M}^{\top}\mathbf{x})\right]^{T},\quad\mathbf{w}_{m}\sim
p(\mathbf{w}).$ (12)

Using Mercer’s theorem, we can equivalently represent a kernel method as a
linear model with respect to the basis function projection,
$\varphi(\mathbf{x})\beta$, given some regression weights
$\beta\in\mathbb{R}^{M}$. This random basis function representation is vital
to make posterior inference tractable in non-Gaussian likelihoods for our
method.

###  2.6 Random Feature Latent Variable Modeling

Lázaro-Gredilla et al., (2010) first used RFFs to approximate a GP-regression
model, but only estimated the random frequencies, $\mathbf{W}$, and kernel
hyperparameters using MAP estimation. Oliva et al., (2016) introduced a fully
Bayesian variant of a random Fourier features, called Bayesian non-parametric
kernel learning (BaNK), and additionally placed Dirichlet process (DP) mixture
of Gaussian-Inverse Wishart distributions on $\mathbf{W}$.

$\displaystyle\mathbf{y}_{j}$
$\displaystyle\sim\mathcal{N}_{N}\big{(}\varphi_{w}(\mathbf{X})\boldsymbol{\beta}_{j},\sigma^{2}_{j}\mathbf{I}\big{)},\boldsymbol{\beta}_{j}\sim\mathcal{N}_{M}(\boldsymbol{\beta}_{0},\mathbf{B}_{0}),\mathbf{w}_{m}\sim\mathcal{N}_{D}(\boldsymbol{\mu}_{\zeta_{m}},\boldsymbol{\Sigma}_{\zeta_{m}}),$
(13) $\displaystyle(\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k})$
$\displaystyle\sim\mathcal{NIW}(\boldsymbol{\mu}_{0},\kappa_{0},\lambda_{0},\boldsymbol{\Psi}_{0}),\zeta_{m}\sim\mathcal{CRP}(\eta),\eta\sim\mbox{Gamma}(a_{\eta},b_{\eta}),\sigma^{-2}_{j}\sim\mbox{Gamma}(a_{0},b_{0}).$

BaNK assigns each $\mathbf{w}_{m}$ in
$\smash{\mathbf{W}=[\mathbf{w}_{1}\dots\mathbf{w}_{M}]^{\top}}$ to a mixture
component using variable $\zeta_{m}$, which is distributed according to a
Chinese restaurant process (CRP, Aldous,, 1985), the posterior predictive
distribution of the Dirichlet process. Using this DP mixture, we can explore
the posterior space of stationary kernels by sampling $\mathbf{w}_{m}$, as a
random Fourier feature variant of the spectral mixture kernel (Wilson and
Adams,, 2013). The parameter $\eta$ is the concentration parameter of the
Dirichlet process, which controls the number of new mixtures the DP mixture
model will allocate.

$\boldsymbol{\beta}$ acts as the regression weights for the random basis
functions of $\varphi_{w}(\mathbf{X})$. Placing a Gaussian prior on
$\boldsymbol{\beta}$ allows and an inverse gamma prior on $\sigma^{2}_{j}$ for
$\boldsymbol{\beta}$ and $\sigma^{2}_{j}$ to be marginalized out in closed
form for the Gaussian model. For count data likelihoods, like the Bernoulli
(Polson et al.,, 2013), negative binomial (Zhou et al., 2012b, ), and
multinomial distributions (Chen et al.,, 2013) we can easily Gibbs sample
$\boldsymbol{\beta}$ in closed form using the Pólya-gamma augmentation.

Gundersen et al., (2021) extended BaNK for approximating GPLVMs and showed
that inference was made tractable for non-Gaussian likelihoods. The resulting
random feature latent variable modeling (RFLVM) extents the BaNK model in
Equation 13 with:

$\displaystyle\mathbf{y}_{j}$
$\displaystyle\sim\mathcal{L}\big{(}g\big{(}\varphi_{w}(\mathbf{X})\boldsymbol{\beta}_{j}\big{)},\boldsymbol{\theta}\big{)},$
$\displaystyle\boldsymbol{\theta}$ $\displaystyle\sim p(\boldsymbol{\theta}),$
$\displaystyle\mathbf{x}_{i}$
$\displaystyle\sim\mathcal{N}_{D}(\mathbf{0},\mathbf{I}).$ (14)

where $\mathcal{L}(\cdot)$ is a likelihood function, $g(\cdot)$ is a link
function that maps the real-valued
$\varphi_{w}(\mathbf{X})\boldsymbol{\beta}_{j}$ to the support of the
likelihood, and $\boldsymbol{\theta}$ are other parameters to $\mathcal{L}$,
if they exist.

##  3 Method

In our method, we extend the contributions of BaNK in Equation 13 and RFLVM in
Equation 14 with a infinite, sparsity inducing prior on $\mathbf{X}$.
Specifically, we place the sparse latent factor model from Knowles and
Ghahramani, (2011) on the RFLVM from Gundersen et al., (2021):

$\displaystyle\mathbf{y}_{j}\sim\mathcal{L}\big{(}g\big{(}\varphi_{w}(\mathbf{X}\odot\mathbf{Z})\boldsymbol{\beta}_{j},\boldsymbol{\theta}\big{)},\;\mathbf{x}_{i}\sim\mathcal{GP}(0,\mathbf{I}),\;\mathbf{w}_{m}\sim\mathcal{GP}(\boldsymbol{\mu}_{\zeta_{m}},\boldsymbol{\Sigma}_{\zeta_{m}}),\;\zeta_{m}\sim\mathcal{CRP}(\eta),$
(15)
$\displaystyle\boldsymbol{\mu}_{k}\sim\mathcal{GP}\left(\boldsymbol{\mu}_{0},\lambda_{0}^{-1}\boldsymbol{\Sigma}_{k}\right),\;\boldsymbol{\Sigma}_{k}^{-1}\sim\mathcal{WP}(\boldsymbol{\Psi}_{0},\nu_{0}),\;\mathbf{Z}\sim\mathcal{IBP}(\alpha),\alpha\sim\mbox{Gamma}(\alpha_{0},\beta_{0}).$

Now that the latent dimension is infinite a priori, the $D$-dimensional priors
are all assumed to be infinite dimension. Thus, the prior on the latent
variables $\mathbf{X}$, the random frequencies $\mathbf{W}$, and the mixture
locations $(\mu_{k},\boldsymbol{\Sigma}_{k}^{-1})$ are now Gaussian process
and Wishart process distributed (Wilson and Ghahramani,, 2011).

###  3.1 Posterior Inference

We adopt a fully Bayesian approach for our proposed model and take samples
from the posterior using MCMC. To sample the sparsity indicators,
$\mathbf{Z}$, we adopt the stick-breaking representation of the IBP (Teh et
al.,, 2007), where we explicitly instantiate the Beta distributed weights of
the latent dimensions according to a generative process for the weights that
exhibit a decreasing ordering:

$\nu_{d}\sim\mbox{Beta}(\alpha,1),\;\pi_{d}=\prod_{d^{\prime}=1}^{d}\nu_{d^{\prime}},\;z_{id}\sim\mbox{Bernoulli}(\pi_{d})$
(16)

This stick-breaking construction is equivalent to the IBP, but allows for an
efficient slice-sampling algorithm for sampling the posterior of $\pi_{d}$ as
opposed to the fully collapsed sampler typically used in the IBP (Damien et
al.,, 1999). We sample the sparsity indicator, $z_{id}$ from the full
conditional:

$P(z_{id}=1|-)\propto\frac{\pi_{d}}{\pi^{\ast}}P(\mathbf{y}_{i}|\mathbf{w}_{d},\mathbf{x}_{i},\mathbf{z}_{i,-d},z_{id}=1,\pi_{d},\boldsymbol{\beta}),\;\pi^{\ast}=\min\left\\{1,\min_{d:\exists
z_{id}=1}\pi_{d}\right\\}.$ (17)

The $\pi^{\ast}$ term in the denominator is necessary for singleton values
where changing the value of $z_{id}$ may result in a different value of the
last active dimension weight. Taking advantage of the conjugate relationship
between the Bernoulli and beta distributions, we can update the
$d=1,\ldots,D^{+}$ instantiated latent dimensions:

$\pi_{d}\sim\mbox{Beta}\left(n_{d},1+N-n_{d}\right),\;n_{d}=\sum_{i=1}^{N}z_{id}$
(18)

To sample new dimensions, we first sample $s\sim\mathcal{U}(0,\pi^{\ast})$, a
slice variable which controls the instantiation of the dimension weights. Then
we draw the new dimension weights, $\pi_{d^{\ast}}\in(0,\pi_{d^{\ast}-1})$ for
$d^{\ast}=1,\ldots,D^{\ast}$, while $\pi_{d^{\ast}}>s$ from the distribution:

$P(\pi_{d^{\ast}})\propto\exp\left(\alpha\sum_{i=1}^{N}\frac{1}{i}\left(1-\pi_{d^{\ast}}\right)^{i}\right)\pi_{d^{\ast}}^{\alpha-1}(1-\pi_{d^{\ast}})^{N},$
(19)

When we instantiate new dimension weights we also instantiate the new
parameters corresponding to the new dimensions,
$(\mathbf{x}_{d^{\ast}},\mathbf{w}_{d^{\ast}},\boldsymbol{\mu}_{d^{\ast}},\boldsymbol{\Sigma}_{d^{\ast}})$,
by sampling from their respective priors.

The full conditional for the IBP concentration parameter $\alpha$ can be
sampled in closed form, if the prior on $\alpha$ is Gamma distributed
(Griffiths and Ghahramani,, 2011):

$P(\alpha|\mathbf{Z})\sim\mbox{Gamma}\left(\alpha_{0}+D^{+},\beta_{0}+\sum_{i=1}^{N}\frac{1}{i}\right)$
(20)

To take posterior samples of $\mathbf{X}$, we use the elliptical slice sampler
(ESS) (Murray et al.,, 2010). This is an efficient MCMC algorithm to take
posterior samples of non-conjugate models with Gaussian priors. We choose to
use this sampling algorithm because it is easy to implement, the is ESS
guaranteed to transition to a new state every iteration, and there are no
tuning parameters to set.

We take posterior draws of $\mathbf{W}$ using a Metropolis–Hastings sampler,
where our proposal distribution is set to the prior distribution on
$\mathbf{w}_{m}$. This simplifies the acceptance probability to the ratio of
likelihoods:

$\displaystyle\mathbf{w}_{m}^{\star}\sim
p(\mathbf{w}_{m}\mid\zeta,\boldsymbol{\mu},\boldsymbol{\Sigma}),\;\rho_{{MH}}=\min\Bigg{\\{}1,\frac{p(\mathbf{Y}\mid\mathbf{X},\mathbf{w}_{m}^{\star},\boldsymbol{\theta})}{p(\mathbf{Y}\mid\mathbf{X},\mathbf{w}_{m},\boldsymbol{\theta})}\Bigg{\\}},$
(21) $\displaystyle
P(\mathbf{w}_{m}|-)\sim\rho_{{MH}}\cdot\delta_{\mathbf{w}_{m}^{\star}}+(1-\rho_{{MH}})\cdot\delta_{\mathbf{w}_{m}}.$

and we sample the latent indicators of the DP mixture,
$\mathbf{\zeta}=[\zeta_{1}\dots\zeta_{M}]$, using the standard “Algorithm 8”
for sampling for Dirichlet process mixture models, with the mixture locations
integrated out (Neal,, 2000):

$\displaystyle
p(\zeta_{m}=k\mid\boldsymbol{\mu},\boldsymbol{\Sigma},\mathbf{W},\alpha)\sim\begin{cases}\frac{n_{k}^{-m}}{M-1+\alpha}\mathcal{T}(\boldsymbol{\mu}_{k}^{\prime},\boldsymbol{\Sigma}_{k}^{\prime},\kappa_{k}^{\prime})&n_{k}^{-m}>0\\\
\frac{\alpha}{M-1+\alpha}\mathcal{T}(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0},\kappa_{0})&n_{k}^{-m}=0.\end{cases}$
(22)

where $\mathcal{T}$ refers to the multivariate-$t$ distribution, parameterized
by the mean, $\boldsymbol{\mu}_{k}^{\prime}$, covariance,
$\boldsymbol{\Sigma}_{k}^{\prime}$, and the degrees of freedom,
$\kappa_{k}^{\prime}$:

$\displaystyle n_{k}$
$\displaystyle=\sum_{m=1}^{M}I(\zeta_{m}=k),\;\bar{\mathbf{w}}^{(k)}=\frac{1}{n_{k}}\sum_{m:\zeta_{m}=k}\mathbf{w}_{m},$
$\displaystyle\kappa_{k}^{\prime}$ $\displaystyle=\kappa_{0}+n_{k}-D^{\ast}+1$
(23) $\displaystyle\boldsymbol{S}_{\bar{\mathbf{w}}_{k}}$
$\displaystyle=\sum_{m:z_{m}=k}(\mathbf{w}_{m}-\bar{\mathbf{w}}^{(k)})(\mathbf{w}_{m}-\bar{\mathbf{w}}^{(k)})^{\top},$
$\displaystyle\boldsymbol{S}_{k}$
$\displaystyle=\sum_{m:z_{m}=k}\frac{\lambda_{0}n_{k}}{\lambda_{0}+n_{k}}(\mathbf{w}_{m}-\boldsymbol{\mu}_{0})(\mathbf{w}_{m}-\boldsymbol{\mu}_{0})^{\top}$
$\displaystyle\boldsymbol{\Sigma}_{k}^{\prime}$
$\displaystyle=\boldsymbol{\Psi}_{0}+\boldsymbol{S}_{\bar{\mathbf{w}}_{k}}+\boldsymbol{S}_{k},$
$\displaystyle\boldsymbol{\mu}_{k}^{\prime}$
$\displaystyle=\frac{\lambda_{0}\boldsymbol{\mu}_{0}+n_{k}\bar{\mathbf{w}}^{k}}{\lambda_{0}+n_{k}},$

For the purpose of convenience, we re-instantiate the parameter locations for
the mixture model on $\mathbf{w}$. We can Gibbs sample the scale and location
parameters using the typical inverse Wishart-Gaussian conjugacy relationship
with a Gaussian likelihood.

$\displaystyle\boldsymbol{\Sigma}_{k}$
$\displaystyle\sim\mathcal{W}^{-1}(\boldsymbol{\Sigma}_{k}^{\prime},\kappa_{0}+n_{k}),\quad\boldsymbol{\mu}_{k}\sim\mathcal{N}(\boldsymbol{\mu}_{k}^{\prime},(\lambda_{0}+n_{k})^{-1}\boldsymbol{\Sigma}_{k}).$
(24)

We sample the DPMM concentration parameter $\alpha$ by using an augmentation
scheme to make sampling $\alpha$ conditionally conjugate with a Gamma prior
(Escobar and West,, 1995):

$\displaystyle\rho$
$\displaystyle\sim\text{Beta}(\alpha+1,M),\frac{\pi_{\rho}}{1-\pi_{\rho}}=\frac{a_{\alpha}+K^{+}-1}{M(b_{\alpha}-\log(\rho))},\;\;K^{+}=\left|\\{k:n_{k}>0\\}\right|,$
(25) $\displaystyle\alpha$
$\displaystyle\sim\pi_{\rho}\cdot\mbox{Gamma}(a_{\alpha}+K^{+},b_{\alpha}-\log(h))+(1-\pi_{\rho})\cdot\mbox{Gamma}(a_{\alpha}+K^{+}-1,b_{\alpha}-\log(\rho)).$

Sampling the regression weights $\boldsymbol{\beta}$ and the other likelihood
specific parameters (if they exist) depend on the likelihood chosen for the
data. For the Gaussian model we integrate out $\boldsymbol{\beta}_{j}$ and
$\sigma^{2}_{j}$, due to the conjugacy between the Gaussian and inverse Gamma
priors with the Gaussian likelihood. For likelihoods that can be written as:

$P(y_{ij}|\mathbf{x}_{i},\beta_{j},\mathbf{w}_{j})\propto
c_{ij}\frac{(\exp(\varphi_{w}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}))^{a_{ij}}}{(1+\exp(\varphi_{w}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}))^{b_{ij}}}$
(26)

which include the Bernoulli, negative binomial and multinomial likelihoods,
then we can Gibbs sample $\boldsymbol{\beta}_{j}$ in closed form using the
Pólya-Gamma augmentation:

$\displaystyle\omega_{ij}$
$\displaystyle\mid\boldsymbol{\beta}_{j}\sim\text{PG}(b_{ij},\varphi_{w}(\mathbf{x}_{i})^{\top}\boldsymbol{\beta}_{j}),$
$\displaystyle\textbf{V}_{\boldsymbol{\omega}_{j}}$
$\displaystyle=(\varphi_{w}(\mathbf{X})^{\top}\boldsymbol{\Omega}_{j}\varphi_{w}(\mathbf{X})+\mathbf{B}_{0}^{-1})^{-1},$
(27) $\displaystyle\boldsymbol{\beta}_{j}$
$\displaystyle\mid\boldsymbol{\Omega}_{j}\sim\mathcal{N}(\textbf{m}_{\boldsymbol{\omega}_{j}},\textbf{V}_{\boldsymbol{\omega}_{j}}),$
$\displaystyle\textbf{m}_{\boldsymbol{\omega}_{j}}$
$\displaystyle=\textbf{V}_{\boldsymbol{\omega}_{j}}(\varphi_{w}(\mathbf{X})^{\top}\boldsymbol{\tau}_{j}+\mathbf{B}_{0}^{-1}\boldsymbol{\beta}_{0}),$

where $\tau_{ij}=a_{ij}-b_{ij}/2$ and
$\boldsymbol{\tau}_{j}=[\tau_{1j}\dots\tau_{Nj}]^{\top}$. The full details of
the Pólya-gamma augmentation are available in the Appendix. The remaining
inference step is sampling the posterior of the other likelihood specific
parameters, $\boldsymbol{\theta}$. In some special cases, we have a Gibbs
sampling update for the parameter, like the dispersion parameter of a negative
binomial distribution (Zhou and Carin,, 2012), but outside of these cases we
require a tailored solution for posterior sampling.

##  4 Experiments

We apply our sparse, infinite RFLVM (IBP RFLVM) model on a wide variety of
datasets. In these experiments we evaluate the performance of our model with
competing models in terms of predictive test set log-likelihood for the real-
valued data sets and test set perplexity on the count-valued data sets. Our
competing models in the real-valued case are the RFLVM (Gundersen et al.,,
2021), the Bayesian inducing point GPLVM (Damianou et al.,, 2016) implemented
in GPy (GPy,, 2012), and the infinite linear factor model (IBP LFM) (Knowles
and Ghahramani,, 2011). In the count-valued case, we compare the sparse and
dense Poisson, negative binomial and multinomial RFLVMs, with Poisson factor
analysis (Zhou et al., 2012a, ), and latent Dirichlet allocation (Griffiths
and Steyvers,, 2004).

In the following experiments, we run each model for $100$ MCMC iterations. We
set the number of random features to be $M=50$ and the number of inducing
points for GPLVM to be equal to $M$. The hyperpriors in the model are set to
be vague. We initialize the number of latent dimensions to be two for the
infinite models and we set the number of latent dimensions for the finite
models to be the number of instantiated dimensions of the comparable infinite
models (ex: we set the number of dimensions for Multinomial RFLVM and LDA to
be the number of instantiated dimensions for IBP Multinomial RFLVM). For test
set evaluation, we randomly hold out $20\%$ of the $y_{ij}$ in $\mathbf{Y}$
and evaluate each method over five trials.

###  4.1 Synthetic Data

We generate a synthetic data set based off the canonical “Cambridge bars” data
set used in many Indian buffet process papers. In our synthetic example we
assume there are four latent dimensions, corresponding to each shape. We
sample the parameters of the IBP RFLVM from the prior data generating process
and set the random frequencies $W$ to be the shapes of the Cambridge data set.
In this data set, we know that the true data generating process is sparse in
the latent space, $X$. The predictive log likelihood results for the held out
test set data in this example show that the IBP RFLVM performs better than the
RFLVM, GPLVM, and the IBP LFM if there is sparsity in the latent space of the
data generating process as well as a non-linear process between the latent
space and observations (see Table 1).

###  4.2 Empirical Data

We evaluate our method on a wide array of empirical data sets. The real-valued
data sets we look at are the Wisconsin breast cancer data set, a protein
expression data set for mice with trisomy and a control group, gene expression
data from a study of ovarian cancer (Martoglio et al.,, 2002), the Olivetti
faces data set. The test set predictive log likelihood is reported Table 1 for
the aforementioned data sets, along with the Cambridge data set. We can see
that the IBP RFLVM exhibits superior performance compared to the other
competing methods in this experiment.

One possible reason why the IBP RFLVM performs so well in these examples is
that sparsity will provide some level of regularization compared to dense
models and produce stronger predictive performance. Another reason is that the
latent space for the data generating processes of these examples indeed
exhibit sparsity, so the prior of the IBP RFLVM is more appropriate than the
dense models which was evidenced in the Cambridge example. The IBP RFLVM and
dense RFLVM generally perform better than the Bayesian GPLVM in terms of
predictive log likelihood because the RFLVM models are exact Bayesian methods
and can better capture the posterior uncertainty compared to the variational
Bayesian GPLVM which will only provide a good point estimate but the mean-
field posterior approximation will poorly capture the posterior uncertainty.
Moreover, the expressivity of the non-linear mapping of IBP RFLVM, RFLVM and
GPLVM produces better predictive performance than the linear IBP LFM, which
performs the worst in this experimental setting.

Next, we evaluate our model on a collection of integer-valued data sets. We
look at the following data sets: A data set that counts the number of bicycles
on different roads in Montreal, a data set counting the words said by members
of the U.S. congress in speeches, a collection of spam and non-spam SMS
messages, and the 20 Newsgroups data set. In Table 2, we can see the test set
perplexity results of the competing methods on these data sets. Except for the
Congress data set, the NB IBP RLFVM produces the best perplexity results.
Furthermore, the IBP variant models generally perform better than their dense
counterparts and the RFLVM models generally perform better than the linear
models. The reason for the superior performance could be for the
aforementioned reasons. Additionally, negative binomial models can model over-
dispersed counts where the variance of the data exceeds the mean, which the
Poisson and multinomial models cannot do.

###  4.3 Qualitative Assessment

Lastly, we can look at the latent space for the Wisconsin breast cancer data
set and the Congress data set and observe a clear delineation between the two
classes in both of the data sets. The label information is not used in
training the model, but we can uncover the two separate classes by inspecting
$\mathbf{X}$. Figures 1 and 2 shows the latent space of these two data sets
for each of the models investigated in this experimental section. The
observations in each of these figures are sorted by their label. For the IBP
RFLVM and RFLVM results on the cancer data set, we can see that there is a
clear delineation between cancerous and non-cancerous samples particularly
with the first two latent dimensions. In the NB IBP RFLVM and NB RFLVM for the
Congress data set, we can see that the second and third latent dimensions are
particularly informative for discriminating between Democratic and Republican
members. This structure is less obvious to see in the other models, which
suggests that the best performing models in our experiments in terms of
predictive performance can actually discover a clear latent structure for some
of our data sets.

(a) IBP RFLVM.

(b) RFLVM.

(c) GPLVM.

(d) IBP LFM.

Figure 1: Latent space for Wisconsin breast cancer data set. The first 212
observations are malignant samples, the remaining are benign.

(a) NB IBP.

(b) NB.

(c) Pois. IBP.

(d) Pois.

(e) Mult. IBP.

(f) Mult.

(g) PFA.

(h) LDA.

Figure 2: Latent space for Congress data set. The first 285 observations are
Democratic members of Congress, the remaining are Republicans. Table 1: Test
set log likelihood for Gaussian models on real-valued datasets over five
trials. One standard error is reported in parentheses.

| Cancer | Mice | Ovarian | Olivetti | Cambridge
---|---|---|---|---|---
IBP RFLVM | 6.68e+03 (114.16) | 1.39e+04 (148.76) | 1.81e+03 (94.02) | 1.43e+05 (726.41) | 1.09e+04 (1212.76)
RFLVM | 6.19e+03 (119.32) | 1.36e+04 (163.09) | 1.73e+03 (93.84) | 1.42e+05 (957.49) | 9.91e+03 (896.54)
GPLVM | -4.79e+03 (26.31) | -1.22e+04 (63.53) | -6.96e+02 (92.65) | -1.16e+05 (168.80) | -2.32e+04 (134.68)
IBP LFM | -2.46e+04 (33.90) | -6.10e+04 (96.61) | -3.99e+03 (21.14) | -5.84e+05 (694.63) | -2.20e+05 (6020.21)

Table 2: Test set perplexity for count data models on integer-valued datasets over five trials. One standard error is reported in parentheses. | Montreal | Congress | Newsgroups | Spam
---|---|---|---|---
NB IBP RFLVM | 16.25 (0.10) | 719.05 (9.59) | 19.14 (0.27) | 158.41 (2.20)
NB RFLVM | 16.69 (0.10) | 727.74 (6.86) | 21.18 (0.12) | 193.26 (2.42)
Poisson IBP RFLVM | 24.76 (0.89) | 2029.75 (38.47) | 23.94 (0.74) | 389.70 (13.02)
Poisson RFLVM | 39.34 (1.58) | 2278.67 (13.82) | 26.66 (0.74) | 410.75 (18.81)
Multinomial IBP RFLVM | 21.58 (0.37) | 1005.36 (12.71) | 20.78 (0.39) | 418.42 (7.09)
Multinomial RFLVM | 22.49 (0.24) | 940.92 (19.86) | 26.53 (1.04) | 370.99 (13.41)
PFA | 27.11 (2.03) | 678.55 (1.94) | 27.43 (1.08) | 186.97 (0.57)
LDA | 117.03 (1.86) | 938.74 (3.44) | 19.16 (0.03) | 218.32 (0.04)

##  5 Conclusion

In this paper we propose a Bayesian model for sparse non-linear latent
variable modeling that can explore the number of latent dimension a posteriori
using the Indian buffet process. Our experimental results show that our method
is generally superior to the dense variant of our model and other competing
linear and non-linear latent variable models. The IBP RFLVM can easily be
generalized to non-Gaussian likelihoods, which is typically difficult for
GPLVMs due to a lack of computational tractability in posterior inference in
the non-Gaussian setting.

In future work, we would like to explore variants of our model using the two-
parameter Indian buffet process. The one-parameter IBP used in this paper
controls the overall sparsity and latent dimensionality with $\alpha$. This
constraint may be too restrictive, so we may want to use the two-parameter IBP
where these two behaviors are governed by their own parameters. Lastly, we are
interested in using other Bayesian non-parametric priors on the sparsity
indicator, $\mathbf{Z}$, like the gamma-Poisson process (Titsias,, 2007) or
the negative binomial process (Zhou and Carin,, 2013) and compare their
performance with our IBP model. Lastly, we are also interested in applying the
IBP RFLVM to scientific settings. For example, GPLVMs are commonly used for
modeling single cell data where the observed data are high-dimensional but we
are interested in modeling a lower dimensional representation of the data. In
such applications, sparse methods are appealing but there are few examples of
sparse, non-linear, non-Gaussian LVMs. We fill this gap with our proposed
model.

## References

  * Aldous, (1985) Aldous, D. J. (1985).  Exchangeability and related topics.  In École d’Été de Probabilités de Saint-Flour XIII—1983, pages 1–198. Springer.
  * Blei et al., (2003) Blei, D. M., Ng, A. Y., and Jordan, M. I. (2003).  Latent Dirichlet allocation.  Journal of Machine Learning Research, 3(Jan):993–1022.
  * Bochner, (1959) Bochner, S. (1959).  Lectures on Fourier integrals, volume 42.  Princeton University Press.
  * Chen et al., (2013) Chen, J., Zhu, J., Wang, Z., Zheng, X., and Zhang, B. (2013).  Scalable inference for logistic-normal topic models.  In Neural Information Processing Systems, pages 2445–2453.
  * Conti et al., (2014) Conti, G., Frühwirth-Schnatter, S., Heckman, J. J., and Piatek, R. (2014).  Bayesian exploratory factor analysis.  Journal of Econometrics, 183(1):31–57.
  * Damianou et al., (2016) Damianou, A. C., Titsias, M. K., and Lawrence, N. D. (2016).  Variational inference for latent variables and uncertain inputs in Gaussian processes.  Journal of Machine Learning Research, 17(1):1425–1486.
  * Damien et al., (1999) Damien, P., Wakefield, J., and Walker, S. G. (1999).  Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(2):331–344.
  * Duane et al., (1987) Duane, S., Kennedy, A., Pendleton, B. J., and Roweth, D. (1987).  Hybrid Monte Carlo.  Physics Letters B, 195(2):216–222.
  * Engelhardt and Stephens, (2010) Engelhardt, B. E. and Stephens, M. (2010).  Analysis of population structure: a unifying framework and novel methods based on sparse factor analysis.  PLoS Genetics, 6(9):e1001117.
  * Escobar and West, (1995) Escobar, M. D. and West, M. (1995).  Bayesian density estimation and inference using mixtures.  Journal of the American Statistical Association, 90(430):577–588.
  * GPy, (2012) GPy (2012).  GPy: A Gaussian process framework in Python.  http://github.com/SheffieldML/GPy. BSD 3-clause license.
  * Griffiths and Ghahramani, (2011) Griffiths, T. L. and Ghahramani, Z. (2011).  The Indian buffet process: An introduction and review.  Journal of Machine Learning Research, 12(4).
  * Griffiths and Steyvers, (2004) Griffiths, T. L. and Steyvers, M. (2004).  Finding scientific topics.  Proceedings of the National Academy of Sciences, 101(1):5228–5235.
  * Gundersen et al., (2021) Gundersen, G. W., Zhang, M. M., and Engelhardt, B. E. (2021).  Latent variable modeling with random features.  Artificial Intelligence and Statistics, 130:1333–1341.
  * Holmes et al., (2006) Holmes, C. C., Held, L., et al. (2006).  Bayesian auxiliary variable models for binary and multinomial regression.  Bayesian Analysis, 1(1):145–168.
  * Knowles and Ghahramani, (2011) Knowles, D. and Ghahramani, Z. (2011).  Nonparametric Bayesian sparse factor models with application to gene expression modeling.  Annals of Applied Statistics, 5(2B):1534–1552.
  * Lawrence, (2005) Lawrence, N. D. (2005).  Probabilistic non-linear principal component analysis with Gaussian process latent variable models.  Journal of Machine Learning Research, 6(Nov):1783–1816.
  * Lázaro-Gredilla et al., (2010) Lázaro-Gredilla, M., Quiñonero-Candela, J., Rasmussen, C. E., and Figueiras-Vidal, A. R. (2010).  Sparse spectrum Gaussian process regression.  Journal of Machine Learning Research, 11:1865–1881.
  * Lee and Seung, (1999) Lee, D. D. and Seung, H. S. (1999).  Learning the parts of objects by non-negative matrix factorization.  Nature, 401(6755):788–791.
  * Martoglio et al., (2002) Martoglio, A.-M., Miskin, J. W., Smith, S. K., and MacKay, D. J. (2002).  A decomposition model to track gene expression signatures: Preview on observer-independent classification of ovarian cancer.  Bioinformatics, 18(12):1617–1624.
  * Mercer, (1909) Mercer, J. (1909).  Functions of positive and negative type, and their connection the theory of integral equations.  Philosophical Transactions of the Royal Society of London, Series A, 209(441-458):415–446.
  * Mitchell and Beauchamp, (1988) Mitchell, T. J. and Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  Journal of the American Statistical Association, 83(404):1023–1032.
  * Murray et al., (2010) Murray, I., Adams, R., and MacKay, D. (2010).  Elliptical slice sampling.  In Artificial Intelligence and Statistics, pages 541–548. PMLR.
  * Neal, (2000) Neal, R. M. (2000).  Markov chain sampling methods for Dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 9(2):249–265.
  * Oliva et al., (2016) Oliva, J. B., Dubey, A., Wilson, A. G., Póczos, B., Schneider, J., and Xing, E. P. (2016).  Bayesian nonparametric kernel-learning.  In Artificial Intelligence and Statistics, pages 1078–1086.
  * (26) Olshausen, B. A. and Field, D. J. (1996a).  Emergence of simple-cell receptive field properties by learning a sparse code for natural images.  Nature, 381(6583):607–609.
  * (27) Olshausen, B. A. and Field, D. J. (1996b).  Natural image statistics and efficient coding.  Network: Computation in Neural Systems, 7(2):333.
  * Pearson, (1901) Pearson, K. (1901).  On lines and planes of closest fit to systems of points in space.  London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):559–572.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using Pólya–gamma latent variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Press and Shigemasu, (1989) Press, S. J. and Shigemasu, K. (1989).  Bayesian inference in factor analysis.  In Contributions to Probability and Statistics, pages 271–287. Springer.
  * Rahimi and Recht, (2008) Rahimi, A. and Recht, B. (2008).  Random features for large-scale kernel machines.  In Neural Information Processing Systems, pages 1177–1184.
  * Teh et al., (2007) Teh, Y. W., Gorür, D., and Ghahramani, Z. (2007).  Stick-breaking construction for the Indian buffet process.  In Artificial Intelligence and Statistics, pages 556–563. PMLR.
  * Thibaux and Jordan, (2007) Thibaux, R. and Jordan, M. I. (2007).  Hierarchical beta processes and the Indian buffet process.  In Artificial Intelligence and Statistics, pages 564–571. PMLR.
  * Tipping and Bishop, (1999) Tipping, M. E. and Bishop, C. M. (1999).  Probabilistic principal component analysis.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3):611–622.
  * Titsias, (2007) Titsias, M. (2007).  The infinite gamma-Poisson feature model.  Neural Information Processing Systems, 20.
  * West, (2003) West, M. (2003).  Bayesian factor regression models in the “large $p$, small $n$” paradigm.  Bayesian Statistics, 7:733–742.
  * Williams and Rasmussen, (2006) Williams, C. K. I. and Rasmussen, C. E. (2006).  Gaussian processes for machine learning, volume 2.  MIT Press.
  * Wilson and Adams, (2013) Wilson, A. G. and Adams, R. P. (2013).  Gaussian process kernels for pattern discovery and extrapolation.  In International Conference on Machine Learning, pages 1067–1075.
  * Wilson and Ghahramani, (2011) Wilson, A. G. and Ghahramani, Z. (2011).  Generalised Wishart processes.  In Uncertainty in Artificial Intelligence, pages 736–744.
  * Zhou and Carin, (2012) Zhou, M. and Carin, L. (2012).  Augment-and-conquer negative binomial processes.  In Neural Information Processing Systems, pages 2546–2554.
  * Zhou and Carin, (2013) Zhou, M. and Carin, L. (2013).  Negative binomial process count and mixture modeling.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(2):307–320.
  * (42) Zhou, M., Hannah, L., Dunson, D., and Carin, L. (2012a).  Beta-negative binomial process and Poisson factor analysis.  In Artificial Intelligence and Statistics, pages 1462–1471. PMLR.
  * (43) Zhou, M., Li, L., Dunson, D., and Carin, L. (2012b).  Lognormal and gamma mixed negative binomial regression.  In International Conference on Machine Learning, volume 2012, page 1343.

##  Appendix A Appendix

###  A.1 Experimental Details

The Wisconsin breast cancer, the mouse trisomy, and spam data sets are from
the UCI Machine Learning Repository. The Olivetti and 20 Newsgroups data sets
are from the scikit-learn package. The Congress data set is available at
https://github.com/jgscott/STA380/blob/master/data/congress109.csv. The
Montreal data set is available at
https://www.kaggle.com/pablomonleon/montreal-bike-lanes. The ovarian cancer
data set is available at: https://tinyurl.com/3chw6znr. For the real-valued
data sets, we standardize the data to have zero mean and unit variance. For
the count data sets, we take the square-root transformation of the data as a
normalization method for variance stabalization.

###  A.2 Pólya-Gamma Augmentation Derivations

Refer to Equation 26, which allows us to rewrite the likelihood as
proportional to a Gaussian:

$P(y_{ij}|\mathbf{x}_{i},\beta_{j},\mathbf{w}_{j})\propto
c_{ij}\frac{(\exp(\varphi_{w}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}))^{a_{ij}}}{(1+\exp(\varphi_{w}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}))^{b_{ij}}}=2^{-b_{ij}}e^{\tau_{ij}\varphi_{w}(\mathbf{x}_{i})_{j}}\int_{0}^{\infty}e^{-\omega\varphi_{w}(\mathbf{x}_{i})_{j}^{2}/2}p(\omega)\text{d}\omega.$
(28)

A random variable $\omega$ is Pólya-gamma distributed with parameters $b>0$
and $c\in\mathbb{R}$, denoted $\omega\sim\text{PG}(b,c)$, if

$\omega\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{2\pi^{2}}\sum_{k=1}^{\infty}\frac{g_{k}}{(k-1/2)^{2}+c^{2}/(4\pi^{2})},$
(29)

where $\stackrel{{\scriptstyle d}}{{=}}$ denotes equality in distribution and
$g_{k}\sim\mbox{Gamma}(b,1)$ are independent gamma random variables. We can
sample $\omega$ conditioned on $\psi_{ij}$ as
$p(\omega\mid\psi_{ij})\sim\text{PG}(b_{ij},\psi_{ij})$. If we set
$a_{ij}=y_{ij}$ and $b_{ij}=1$, then we have a sampler for Bernoulli
observations (Polson et al.,, 2013). If we set $a_{ij}=y_{ij}$ and
$b_{ij}=y_{ij}+r_{j}$, then we have a sampler for negative binomial
observations for a dispersion parameter $r_{j}$ (Zhou et al., 2012b, ).

The derivation for the Pólya-gamma augmentation with the multinomial
likelihood depends on the reparameterization of the likelihood (Holmes et
al.,, 2006; Chen et al.,, 2013). We rewrite the likelihood as

$\displaystyle\begin{split}p(\mathbf{Y}|-)&=\prod_{i=1}^{N}\frac{\Gamma\left(\sum_{j=1}^{J}y_{ij}+1\right)}{\prod_{j=1}^{J}\Gamma\left(y_{ij}+1\right)}\prod_{j=1}^{J}\left(\frac{\exp\left\\{\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}\right\\}}{\sum_{j=1}^{J}\exp{\left\\{\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}\right\\}}}\right)^{y_{ij}}\\\
&\propto\prod_{i=1}^{N}\prod_{j=1}^{J}\frac{\left(\exp\left\\{\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\xi_{ij}\right\\}\right)^{y_{ij}}}{\left(1+\exp\left\\{\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\xi_{ij}\right\\}\right)^{y_{ij}+\sum_{j=1}^{J}y_{ij}}}\end{split}$
(30)

Where $\xi_{ij}=\log\sum_{j^{\prime}\neq
j}\exp\\{\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j^{\prime}}\\}$.
By convention and for identifiability purposes, we set
$\boldsymbol{\beta}_{J}=0$. We let $\tau_{ij}=y_{ij}-\sum_{j=1}^{J}y_{ij}/2$.
The likelihood is proportional to:

$\displaystyle\begin{split}p(\mathbf{Y}|-)&\propto\prod_{i=1}^{N}\prod_{j=1}^{J}\exp\Big{\\{}\tau_{ij}\Big{(}\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\xi_{ij}\Big{)}-\frac{\omega_{nj}}{2}\Big{(}\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\xi_{ij}\Big{)}^{2}\Big{\\}}.\end{split}$
(31)

So this gives us a posterior w.r.t. $\boldsymbol{\beta}_{j}$ as

$p(\boldsymbol{\beta}_{j}\mid\mathbf{y}_{j},\mathbf{X})\propto
p(\boldsymbol{\beta}_{j})\prod_{n=1}^{N}\exp\Big{\\{}\boldsymbol{\kappa}_{i}\Big{(}\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\boldsymbol{\xi}_{j}\Big{)}-\frac{1}{2}\Big{(}\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\boldsymbol{\xi}_{j}\Big{)}^{T}\boldsymbol{\Omega}_{n}\Big{(}\varphi_{\mathbf{W}}(\mathbf{x}_{i})\boldsymbol{\beta}_{j}-\boldsymbol{\xi}_{j}\Big{)}\Big{\\}}$
(32)

which we can rewrite into a closed form update as

$\displaystyle\boldsymbol{\beta}_{j}\mid\boldsymbol{\omega}_{j}$
$\displaystyle\sim\mathcal{N}(\textbf{m}_{\boldsymbol{\omega}_{j}},\textbf{V}_{\boldsymbol{\omega}_{j}})$
(33)

where

$\displaystyle\boldsymbol{\Omega}_{j}$
$\displaystyle=\text{diag}([\boldsymbol{\omega}_{1j},\dots,\boldsymbol{\omega}_{Nj}])$
(34) $\displaystyle\mathbf{V}_{\boldsymbol{\omega}_{j}}$
$\displaystyle=(\varphi_{w}(\mathbf{X})^{\top}\boldsymbol{\Omega}_{j}\varphi_{w}(\mathbf{X})+\mathbf{B}_{0}^{-1})^{-1},$
$\displaystyle\textbf{m}_{\boldsymbol{\omega}_{j}}$
$\displaystyle=\textbf{V}_{\boldsymbol{\omega}_{j}}(\varphi_{w}(\mathbf{X})^{\top}(\boldsymbol{\tau}_{j}+\boldsymbol{\xi}_{j}^{T}\boldsymbol{\Omega}_{j})+\mathbf{B}_{0}^{-1}\boldsymbol{\beta}_{0}),$
$\displaystyle\boldsymbol{\tau}_{j}$
$\displaystyle=\mathbf{y}_{j}-\frac{1}{2}\sum_{j=1}^{J}y_{ij}$

and we sample $\boldsymbol{\Omega}_{j}$ with

$\displaystyle\boldsymbol{\omega}_{j}\mid\boldsymbol{\beta}_{j}$
$\displaystyle\sim\text{PG}\left(\sum_{j=1}^{J}y_{ij},\varphi_{w}(\mathbf{X})\boldsymbol{\beta}_{j}-\boldsymbol{\xi}_{j}\right)$
(35)


Generate $I\sim\textrm{uniform}(\\{1,\cdots,n+i-1\\})$ $\triangleright$
Branching individual

$I_{\textrm{div}}[i]\leftarrow I$

Generate $v,q\sim\textrm{uniform}(0,1)$

append $F^{-1}_{\textrm{div}}(u)$ to $L_{\textrm{div}}$ $\triangleright$
Fraction $\lambda$ at birth

if $E=i$ and $p>F^{-1}_{\textrm{div}}(u)$ then $E\leftarrow n+i$
$\triangleright$ New spinal individual position

end for

end function

From the list $I_{\textrm{div}}$, it is possible to retrieve the lineage of
all individuals using any chosen labeling method. In our case we used the
Ulam-Harris-Neveu notations and we propose the following algorithm to generate
the list $U$ of labels from the list $I$ previously obtained and the initial
size of the population $n$.

function $U$ = Labels($I_{\textrm{div}},n$)

$N_{\textrm{div}}=\textrm{length}(I_{\textrm{div}})$

$U=[1,\cdots,n,0,\cdots,0]$ $\triangleright$ $U$ is of size
$n+N_{\textrm{div}}$

for $i\in\\{1,\cdots,N_{\textrm{div}}\\}$ do

$(U[I_{\textrm{div}}[i]],U[n+i]])\leftarrow(U[I_{\textrm{div}}[i]]1,U[I_{\textrm{div}}[i]]2)$
$\triangleright$ $U[I_{\textrm{div}}[i]]\in\mathcal{U}$

end for

$U\leftarrow\textrm{QuickSort}(U)$ $\triangleright$ sorting algorithm in
$O((N_{\textrm{div}}+n)\log(N_{\textrm{div}}+n))$

end function

The last line of this algorithm lists the labels of individuals in increasing
order, grouping siblings together. The label of the spinal individual is thus
$U[E]$. To compute a trajectory of the spinal process on the interval $[0,T]$,
we need to establish the time of jumps and their outcomes. Using the
deterministic evolution of the traits and knowing the branching events, it is
easy to recover the traits at all time for the individuals. Furthermore,
depending on the statistics that one want to evaluate on the system, it might
not be necessary to compute the traits of all individuals. For that reason we
propose an algorithmic construction of a trajectory of the spinal process,
returning the list of jumps times, events and labels of the individuals living
at time $T$. This algorithm also distinguishes the spinal individual in the
population by returning the label of the spinal individual in the living ones
at time $T$. To simplify we will denote
$\textrm{PPP}_{[t_{1},t_{2}]}(c(\cdot))$ the list of times given by a Poisson
point process of intensity $c(\cdot)$ on the interval $[t_{1},t_{2}]$,
computed using Lewis’ thinning algorithm [39].

Algorithm 1 Simulation algorithm of the jumps times until time $T$

model parameters:
$d,F^{-1}_{\textrm{div}},F^{-1}_{\textrm{loss}},F^{*-1}_{\textrm{loss}},T$,
initial condition: $z=[x^{1},\cdots,x^{n}]$.

$(T_{\textrm{div}},I_{\textrm{div}},E,L_{\textrm{div}},T_{\textrm{loss}},T_{\textrm{loss}}^{*},I_{\textrm{loss}},L_{\textrm{loss}},L_{\textrm{loss}}^{*})$

% Initializing the spine

Generate $u\sim$uniform$(0,1)$

$i_{0}=1$

while $\sum_{i=1}^{i_{0}}x^{i}/(\sum_{i=1}^{n}x^{i})<u$ do

$i_{0}\leftarrow i_{0}+1$

end while

% Generating binary spinal tree on division events

$(T_{\textrm{div}},I_{\textrm{div}},E,L_{\textrm{div}})$ =
Tree($0,T,n,i_{0},F^{-1}_{\textrm{div}}$)

% Generating loss events for the spine and individuals outside the spine

$\left(T_{\textrm{loss}},T_{\textrm{loss}}^{*},I_{\textrm{loss}}\right)=\left([\cdot],[\cdot],[\cdot]\right)$

for $i\in\\{1,\cdots,\textrm{length}(T_{\textrm{div}})\\}$ do

% For the individuals outside the spine

Generate
$P\sim\textrm{PPP}_{\big{[}T_{\textrm{div}}[i-1],T_{\textrm{div}}[i]\big{]}}(K_{\textrm{loss}}(n+i-1)^{2}d(\cdot))$$\triangleright$
Times of loss events

append $P$ to $T_{\textrm{loss}}$

Generate
$I_{u}\sim\left(\textrm{uniform}\left(\\{1,\cdots,n+i-1\\}\right)^{\textrm{length}(T_{\textrm{loss}}[i])}\right)$
$\triangleright$ Distributing loss events

append $I_{u}$ to $I_{\textrm{loss}}$

% For the spine

Generate
$P^{*}\sim\textrm{PPP}_{\big{[}T_{\textrm{div}}[i-1],T_{\textrm{div}}[i]\big{]}}((n+i-1)d(\cdot))$
$\triangleright$ Times of loss events

append $P^{*}$ to $T^{*}_{\textrm{loss}}$

end for

Generate
$L_{u}\sim\left(\textrm{uniform}(0,1)^{\textrm{length}(T_{\textrm{loss}})}\right)$
$\triangleright$ Generating fractions lost

$L_{\textrm{loss}}=F^{-1}_{\textrm{loss}}(L_{u})$

Generate
$L_{u}^{*}\sim\left(\textrm{uniform}(0,1)^{\textrm{length}(T_{\textrm{loss}}^{*})}\right)$
$\triangleright$ Generating fractions lost

$L_{\textrm{loss}}^{*}=F^{*-1}_{\textrm{loss}}(L_{u}^{*})$

The output tuple
$(T_{\textrm{div}},I_{\textrm{div}},E,L_{\textrm{div}},T_{\textrm{loss}},T_{\textrm{loss}}^{*},I_{\textrm{loss}},L_{\textrm{loss}},L_{\textrm{loss}}^{*})$
is the minimal information needed to construct the spinal process introduced
in Section 4. The list $U$ of labels of the individuals can be computed with
the function Labels($I_{\textrm{div}},n$). Notice that a lot of operations can
be parallelized in this algorithm, unlike in the classical Lewis’ algorithm.
Moreover, depending on the statistic that one wants to compute on this
process, the algorithm can be further simplified. For example, if one takes
interest in the total biomass of the population, the indexes of the
individuals that branched are not necessary.

## Acknowledgments

I express my gratitude to L. Coquille, A. Marguet, and C. Smadi for their
guidance and valuable feedback throughout this work. I would like to extend
special thanks to A. Marguet and C. Smadi for their numerous corrections made
to this article. I thank V. Bansaye for stimulating discussions on the subject
of this paper and S. Billiard for fruitful discussions on the Yule model. I
acknowledge partial support by the Chair ”Modélisation Mathématique et
Biodiversité” of VEOLIA-Ecole Polytechnique-MNHN-F.X. This work is supported
by the French National Research Agency in the framework of the ”France 2030”
program (ANR-15-IDEX-0002) and by the LabEx PERSYVAL-Lab
(ANR-11-LABX-0025-01).

## References

  * [1] David J. Aldous.  Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today.  Statistical Science, 16(1):23–34, 2001.
  * [2] Krishna B. Athreya.  Change of measures for Markov chains and the LlogL theorem for branching processes.  Bernoulli, 6(2):323–338, 2000.
  * [3] Vincent Bansaye.  Spine for interacting populations and sampling.  arXiv preprint, 2105.03185, 2022.
  * [4] Vincent Bansaye, Jean-François Delmas, Laurence Marsalle, and Viet Chi Tran.  Limit theorems for Markov processes indexed by continuous time Galton-Watson trees.  The Annals of Applied Probability, 21(6):2263–2314, 2011.
  * [5] Vincent Bansaye and Sylvie Méléard.  Stochastic models for structured populations, volume 16.  Springer, 2015.
  * [6] Jean Bertoin and Bastien Mallein.  Biggins’ martingale convergence for branching Lévy processes.  Electronic Communications in Probability, 23:1 – 12, 2018.
  * [7] J. D. Biggins and A. E. Kyprianou.  Measure change in multitype branching.  Advances in Applied Probability, 36(2):544–581, 2004.
  * [8] John D Biggins.  Martingale convergence in the branching random walk.  Journal of Applied Probability, 14(1):25–37, 1977.
  * [9] K.J Brown and S.S Lin.  On the existence of positive eigenfunctions for an eigenvalue problem with indefinite weight function.  Journal of Mathematical Analysis and Applications, 75(1):112–120, 1980.
  * [10] Nicolas Champagnat and Sylvie Méléard.  Invasion and adaptive evolution for individual-based spatially structured populations.  Journal of Mathematical Biology, 55:147–188, 2007.
  * [11] Brigitte Chauvin and Alain Rouault.  K-P-P equation and supercritical branching brownian motion in the subcritical speed area. application to spatial trees.  Probability theory and related fields, 80(2):299–314, 1988.
  * [12] Brigitte Chauvin, Alain Rouault, and Anton Wakolbinger.  Growing conditioned trees.  Stochastic Processes and their Applications, 39(1):117–130, 1991\.
  * [13] David Cheek and Samuel GG Johnston.  Ancestral reproductive bias in branching processes.  Journal of Mathematical Biology, 86(5):70, 2023.
  * [14] Raphaël Chetrite and Hugo Touchette.  Nonequilibrium Markov processes conditioned on large deviations.  Annales Henri Poincaré, 16:2005–2057, 2015.
  * [15] Bertrand Cloez.  Limit theorems for some branching measure-valued processes.  Advances in Applied Probability, 49(2):549–580, 2017.
  * [16] Joel E Cohen.  Population growth and earth’s human carrying capacity.  Science, 269(5222):341–346, 1995.
  * [17] Jérôme Coville.  On a simple criterion for the existence of a principal eigenfunction of some nonlocal operators.  Journal of Differential Equations, 249(11):2921–2953, 2010.
  * [18] Pierre Del Moral.  Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications.  Springer New York, New York, NY, 2004.
  * [19] Luc Devroye.  General Principles in Random Variate Generation, pages 27–82.  Springer New York, New York, NY, 1986.
  * [20] Rick Durrett.  Probability: theory and examples, volume 49.  Cambridge university press, 2019.
  * [21] Maren Eckhoff, Andreas E. Kyprianou, and Matthias Winkel.  Spines, skeletons and the strong law of large numbers for superdiffusions.  The Annals of Probability, 43(5):2545–2610, 2015.
  * [22] János Engländer.  Branching diffusions, superdiffusions and random media.  Probability Surveys, 4:303 – 364, 2007.
  * [23] János Engländer, Simon C. Harris, and Andreas E. Kyprianou.  Strong Law of Large Numbers for branching diffusions.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 46(1):279 – 298, 2010.
  * [24] Nicolas Fournier and Sylvie Méléard.  A microscopic probabilistic description of a locally regulated population and macroscopic approximations.  The Annals of Applied Probability, 14(4):1880 – 1919, 2004.
  * [25] Coralie Fritsch, Jérôme Harmand, and Fabien Campillo.  A modeling approach of the chemostat.  Ecological Modelling, 299:1–13, 2015.
  * [26] Hans-Otto Georgii and Ellen Baake.  Supercritical multitype branching processes: The ancestral types of typical individuals.  Advances in Applied Probability, 35(4):1090–1110, 2003.
  * [27] Y. Git, J. W. Harris, and S. C. Harris.  Exponential growth rates in a typed branching diffusion.  The Annals of Applied Probability, 17(2):609–653, 2007.
  * [28] Robert Hardy and Simon C. Harris.  A spine approach to branching diffusions with applications to Lp-convergence of martingales.  Séminaire de Probabilités XLII, 1979:281–330, 2009.
  * [29] Simon C. Harris, Samuel G. G. Johnston, and Matthew I. Roberts.  The coalescent structure of continuous-time Galton-Watson trees.  arXiv preprint, 1703.00299, 2019.
  * [30] Simon C. Harris and Matthew I. Roberts.  The many-to-few lemma and multiple spines.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 53(1):226 – 242, 2017.
  * [31] Simon C Harris and David Williams.  Large-deviations and martingales for a typed branching diffusion, 1.  Asterique, 236:133–154, 1996.
  * [32] Theodore Edward Harris et al.  The theory of branching processes, volume 6.  Springer Berlin, 1963.
  * [33] Benoît Henry, Sylvie Méléard, and Viet Chi Tran.  Time reversal of spinal processes for linear and non-linear branching processes near stationarity.  Electronic Journal of Probability, 28:1 – 27, 2023.
  * [34] Samuel G. G. Johnston and Amaury Lambert.  The coalescent structure of uniform and Poisson samples from multitype branching processes.  arXiv preprint, 1912.00198, 2021.
  * [35] Olav Kallenberg.  Stability of critical cluster fields.  Mathematische Nachrichten, 77(1):7–43, 1977.
  * [36] H. Kesten and B. P. Stigum.  A limit theorem for multidimensional Galton-Watson processes.  The Annals of Mathematical Statistics, 37(5):1211–1223, 1966.
  * [37] Thomas Kurtz, Russell Lyons, Robin Pemantle, and Yuval Peres.  A conceptual proof of the Kesten-Stigum theorem for multi-type branching processes.  Classical and Modern Branching Processes, 84:181–185, 1997.
  * [38] A.E. Kyprianou.  Travelling wave solutions to the K-P-P equation: alternatives to Simon Harris’ probabilistic analysis.  Annales de l’Institut Henri Poincare (B) Probability and Statistics, 40(1):53–72, 2004.
  * [39] P. A. W Lewis and G. S. Shedler.  Simulation of nonhomogeneous Poisson processes by thinning.  Naval Research Logistics Quarterly, 26(3):403–413, 1979.
  * [40] Rong-Li Liu, Yan-Xia Ren, and Renming Song.  LlogL criterion for a class of superdiffusions.  Journal of Applied Probability, 46(2):479–496, 2009.
  * [41] Rong-Li Liu, Yan-Xia Ren, and Renming Song.  LlogL condition for supercritical branching Hunt processes.  Journal of Theoretical Probability, 24(1):170–193, 2011.
  * [42] Rong-Li Liu, Yan-Xia Ren, and Renming Song.  Strong Law of Large Numbers for a class of superdiffusions.  Acta applicandae mathematicae, 123(1):73–97, 2013.
  * [43] Russell Lyons.  A simple path to Biggins’ martingale convergence for branching random walk.  Classical and Modern Branching Processes, 84:217–221, 1997.
  * [44] Russell Lyons, Robin Pemantle, and Yuval Peres.  Conceptual proofs of LLogL criteria for mean behavior of branching processes.  The Annals of Probability, 23(3):1125–1138, 1995.
  * [45] Aline Marguet.  Uniform sampling in a structured branching population.  Bernoulli, 25(4A):2649 – 2695, 2019.
  * [46] R.E. Miles.  On the homogeneous planar Poisson point process.  Mathematical Biosciences, 6:85–127, 1970.
  * [47] J. Neveu.  Arbres et processus de Galton-Watson.  Annales de l’I.H.P. Probabilités et statistiques, 22(2):199–207, 1986.
  * [48] Y. Ogata.  On Lewis’ simulation method for point processes.  IEEE Transactions on Information Theory, 27(1):23–31, 1981.
  * [49] Paul A. Rees.  Chapter 5 - Population Growth, pages 121–137.  Academic Press, 2018.
  * [50] Yan-Xia Ren, Renming Song, and Ting Yang.  Spine decomposition and LlogL criterion for superprocesses with non-local branching mechanisms.  arXiv preprint, 1609.02257, 2020.