# Coherent phonon spectroscopy and interlayer modulation of charge density
wave order in the kagome metal CsV3Sb5

Noah Ratcliff  Materials Department, University of California, Santa Barbara,
California 93106, USA  Lily Hallett  Materials Department, University of
California, Santa Barbara, California 93106, USA  Brenden R. Ortiz  Materials
Department, University of California, Santa Barbara, California 93106, USA
Stephen D. Wilson  Materials Department, University of California, Santa
Barbara, California 93106, USA  John W. Harter<EMAIL_ADDRESS>Materials
Department, University of California, Santa Barbara, California 93106, USA

###### Abstract

The recent discovery of the $A$V3Sb5 ($A=$ K, Rb, Cs) material family offers
an exciting opportunity to investigate the interplay of correlations,
topology, and superconductivity in kagome metals. The low energy physics of
these materials is dominated by an unusual charge density wave phase, but
little is understood about the true nature of the order parameter. In this
work, we use a combination of ultrafast coherent phonon spectroscopy and
first-principles density functional theory calculations to investigate the
charge density wave order in CsV3Sb5. We find that the charge density wave is
the result of a simultaneous condensation of three optical phonon modes at one
$M$ and two $L$ points. This distortion can be described as tri-hexagonal
ordering with an interlayer modulation along the $c$ axis. It breaks the
$C_{6}$ rotational symmetry of the crystal and may offer a natural explanation
for reports of uniaxial order at lower temperatures in this material family.

Identifying and characterizing novel phases of matter is at the forefront of
materials science research. In this regard, kagome materials hold great
promise for realizing exotic correlated and topological ground states [1, 2,
3, 4, 5, 6, 7, 8, 9, 10]. It therefore comes as no surprise that the recently
discovered quasi-two-dimensional vanadium-based kagome metal series $A$V3Sb5
($A=$ K, Rb, Cs) has generated a flurry of interest [11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. These materials
possess a rare and coveted combination of nontrivial band topology and
superconductivity. Furthermore, the high quality and exfoliable nature of
single crystal samples makes them accessible to a wide variety of experimental
techniques, which have already uncovered a complex web of intertwined
properties. For example, the observation of a giant anomalous Hall effect [16,
17] and chiral charge order [18] suggest proximity to a time-reversal symmetry
breaking instability, even while there is no evidence of local moment
magnetism [12, 15]. In addition, thermodynamic signatures of nodal
quasiparticles [19], multiple superconducting domes [20, 21, 22], spin-triplet
supercurrents [23], and zero-bias conductance peaks inside superconducting
vortex cores [24] hint at the likelihood of unconventional—possibly
topological—superconductivity.

The low energy physics of $A$V3Sb5 is dominated by an unusual charge density
wave (CDW) phase transition at $T_{\mathrm{CDW}}=$ 78, 102, and 94 K for $A=$
K, Rb, and Cs, respectively [12, 13, 14]. This phase transition can be clearly
identified through anomalies in heat capacity, electrical resistivity, and
magnetic susceptibility, and is widely suspected to be related to a Peierls-
like nesting-driven instability at wave vectors connecting the $M$ points at
the Brillouin zone boundary. However, the true nature of the CDW order
parameter is not yet fully understood, and there are a number of open
questions. One question concerns whether or not the CDW modulates along the
$c$ axis [24, 25, 26]. A second question involves reports of uniaxial CDW
order which onsets well below $T_{\mathrm{CDW}}$ [27, 24, 28, 29]. As
superconductivity competes with and emerges from the CDW phase, understanding
the CDW order in $A$V3Sb5 is essential. In this work, we use ultrafast optical
transient reflectivity experiments in conjunction with first-principles
density functional theory (DFT) calculations to study the CDW phase transition
in CsV3Sb5. We show that the CDW is modulated along the $c$ axis and
corresponds to a simultaneous condensation of three optical phonon modes ($3Q$
order) at one $M$ and two $L$ points. This “MLL” distortion may offer an
explanation for the uniaxial order observed in CsV3Sb5 at lower temperatures.

Figure 1:  Coherent phonon spectroscopy data for CsV3Sb5. (a) Transient
reflectivity curves for a series of temperatures above and below the CDW phase
transition. Temperatures are approximately evenly spaced and scans are offset
for clarity. Above $T_{\mathrm{CDW}}$, the reflectivity increases after the
pump pulse and a single oscillation frequency is apparent. Below
$T_{\mathrm{CDW}}$, in contrast, the reflectivity decreases after the pump
pulse and the presence of multiple oscillation frequencies results in a
complex beat pattern. (b) Magnitude of the Fourier transform of the
reflectivity oscillations after subtraction of a double exponential
background. Three resonances are identified and labeled $\alpha$ (1.3 THz),
$\beta$ (3.1 THz), and $\gamma$ (4.1 THz). (c) Two-dimensional temperature-
frequency map of the Fourier magnitude of the coherent phonon oscillations.
While resonance $\gamma$ is present at all temperatures, resonance $\alpha$
only becomes Raman active below $T_{\mathrm{CDW}}$. Below $T^{*}\approx 60$ K,
the broad $\beta$ resonance appears and gradually grows in amplitude. A weak
softening of all three frequencies is apparent with increasing temperature.

Time-resolved optical reflectivity measurements were performed on freshly-
cleaved surfaces of CsV3Sb5 single crystals mounted in an optical cryostat.
Synthesis is described in Ref. 12. A non-collinear optical parametric
amplifier was used to generate $\sim 50$ fs signal (800 nm) and idler (1515
nm) pulses at a repetition rate of 500 kHz, which were used as probe and pump
beams, respectively. To avoid sample heating, both pulses had a low fluence of
$\sim 100$ $\mu$J/cm2 and were linearly polarized in-plane. A lock-in
amplifier and optical chopper were used to measure the small pump-induced
transient change in reflectivity. Phonon frequencies were calculated using DFT
within the Perdew-Burke-Ernzerhof generalized gradient approximation, as
implemented in the Vienna ab initio Simulation Package (vasp) [32]. The
projector augmented wave potentials considered 9 valence electrons for the
cesium atoms, the plane wave basis cutoff energy was 300 eV, and the zero
damping DFT-D3 van der Waals correction was employed. The unit cell structure
was relaxed using a $\Gamma$-centered $18\times 18\times 12$ k-point mesh. The
relaxed lattice parameters were $a=b=5.450$ Å and $c=9.297$ Å and the out-of-
plane antimony atoms were located at fractional height $z=0.7435$. Phonon
dispersion relations were calculated with the phonopy software package [33]
via the finite displacement method using a $3\times 3\times 2$ supercell. A
$2\times 2\times 2$ supercell was used to calculate the energies of modes
simultaneously condensing at the $M$ and $L$ points.

Figure 1(a) shows raw transient reflectivity data measured for CsV3Sb5 at
several temperatures. Above $T_{\mathrm{CDW}}$, the reflectivity increases
after the pump pulse and shows clear phonon oscillations at a single
frequency. Below $T_{\mathrm{CDW}}$, however, the reflectivity decreases and
exhibits a complex beat pattern indicative of the presence of multiple
oscillation frequencies. The stark difference in the transient optical
response across the CDW phase transition ($\Delta R/R>0$ above
$T_{\mathrm{CDW}}$, $\Delta R/R<0$ below $T_{\mathrm{CDW}}$) is likely the
result of changes in the density of states and the opening of a partial energy
gap near the Fermi level at the CDW phase transition [30]. Indeed, the sign of
$\Delta R/R$ is known to depend sensitively on several electronic parameters
influencing carrier dynamics, such as band filling and band gap
renormalization [34]; qualitative changes in the transient optical response
are therefore unsurprising at $T_{\mathrm{CDW}}$.

To better understand the nature of the coherent phonon oscillations, a double
exponential background ($A_{0}+A_{1}e^{-t/\tau_{1}}+A_{2}e^{-t/\tau_{2}}$) is
fitted and then subtracted from the reflectivity curves after $t=100$ fs. The
oscillations that remain are Fourier transformed and displayed in Fig. 1(b).
Three resonance modes are clearly present in the material: $\alpha$ at 1.3
THz, $\beta$ at 3.1 THz, and $\gamma$ at 4.1 THz. No other frequencies are
detected up to 10 THz. The striking temperature dependence of the coherent
phonon oscillations can be more clearly seen in Fig. 1(c), which displays a
two-dimensional temperature-frequency map of the oscillation spectrum. While
the $\gamma$ mode is present at all temperatures, the $\alpha$ mode shows a
sudden, intense appearance at $T\approx 92$ K. This temperature is close to
$T_{\mathrm{CDW}}=94$ K, as determined through independent heat capacity
measurements [12], and we therefore identify it with the CDW critical
temperature. The slightly lower detection temperature can be explained through
modest local heating of the sample by the laser. At even lower temperatures,
below $T^{*}\approx 60$ K, the broad $\beta$ mode appears. All three
resonances show a weak frequency softening as temperature is increased, but
there is no evidence of complete softening ($\omega\rightarrow 0$) for any of
the modes.

Figure 2:  Behavior near the CDW phase transition. (a) Transient reflectivity
curves for a finely-spaced series of temperatures across the CDW phase
transition. Temperatures are approximately evenly spaced. Curves are not
offset and finite $\Delta R/R$ values for $t<0$ are real. Inset shows
pronounced changes in the transient optical response above (red), at (green),
and below (blue) $T_{\mathrm{CDW}}$. (b) Average $\Delta R/R$ value for time
delays between $-1$ ps and $-0.2$ ps versus temperature. A metastable-like
divergence in the lifetime of the transient optical response is observed at
$T_{\mathrm{CDW}}$. Insets are cartoons of the free energy and the gray curve
is a guide to the eye. (c) Average $\left|\Delta R/R\right|$ value for
negative time delays versus pump-probe pulse period at the CDW phase
transition. The fit is described in the main text.

To investigate the critical behavior of the phase transition, transient
reflectivity data were collected within a narrow temperature range centered at
$T_{\mathrm{CDW}}$. As Fig. 2(a) illustrates, we observed a pronounced
qualitative change in the shape of the transient reflectivity curve within
$\sim 1$ K of $T_{\mathrm{CDW}}$. In particular, we measured an anomalously
long recovery behavior together with a finite $\Delta R/R$ value for negative
time delays, as shown in Fig. 2(b). This is likely the result of a pump-
induced metastable change in the material, such that recovery to equilibrium
after one pulse is still occurring by the time the next pulse arrives. To test
this hypothesis, we measured the negative time delay anomaly as a function of
the pulse period $T$. Figure 2(c) shows that a double exponential recovery,
$\Delta R(t<0)/R=A_{1}e^{-T/\tau_{1}}+A_{2}e^{-T/\tau_{2}}$, is consistent
with our data, with a least squares fit giving $\tau_{1}=0.92$ $\mu$s and
$\tau_{2}=11.7$ $\mu$s. These values represent extremely long lifetimes for
typical electronic and phononic degrees of freedom in a solid, and instead
point to metastable phase coexistence associated with a first-order phase
transition. We therefore conclude that the CDW phase transition in CsV3Sb5 is
discontinuous.

Figure 3:  First-principles phonon calculations. (a) Unit cell for CsV3Sb5. The vanadium atoms (red) form a perfect kagome net and are coordinated with in-plane and out-of-plane antimony atoms. (b) First Brillouin zone of the hexagonal lattice, with high-symmetry points labeled. (c) Calculated phonon dispersion relations. There are two unstable modes with imaginary frequencies: one at the $M$ point with irreducible representation $M_{1}^{+}$, and the other at the $L$ point with irreducible representation $L_{2}^{-}$. (d) Illustration of the fully symmetric $\Gamma_{1}^{+}$ phonon at 4.10 THz experimentally detected at all temperatures. (e) Illustration of the $L_{2}^{-}$ phonon at 1.27 THz experimentally detected below $T_{\mathrm{CDW}}$. Table 1:  Calculated phonon frequencies up to 5 THz. $\Gamma$ point |  $M$ point |  $L$ point
---|---|---
Irrep |  $f$ (THz) | Irrep |  $f$ (THz) | Irrep |  $f$ (THz)
$\Gamma_{6}^{-}$ | 1.28 | $M_{1}^{+}$ | $-2.03$ | $L_{2}^{-}$ | $-2.40$
$\Gamma_{2}^{-}$ | 1.51 | $M_{2}^{-}$ | 1.11 | $L_{2}^{-}$ | 1.27
$\Gamma_{6}^{-}$ | 1.99 | $M_{3}^{-}$ | 1.30 | $L_{3}^{-}$ | 1.30
$\Gamma_{4}^{-}$ | 2.16 | $M_{3}^{+}$ | 1.33 | $L_{4}^{-}$ | 1.43
$\Gamma_{6}^{+}$ | 2.19 | $M_{4}^{-}$ | 1.49 | $L_{4}^{-}$ | 1.52
$\Gamma_{2}^{-}$ | 2.70 | $M_{2}^{+}$ | 1.79 | $L_{1}^{-}$ | 1.78
$\Gamma_{5}^{-}$ | 3.30 | $M_{2}^{-}$ | 1.88 | $L_{1}^{+}$ | 1.84
$\Gamma_{5}^{+}$ | 3.81 | $M_{1}^{-}$ | 2.39 | $L_{2}^{+}$ | 2.39
$\Gamma_{1}^{+}$ | 4.10 | $M_{2}^{-}$ | 2.78 | $L_{1}^{+}$ | 2.76
$\Gamma_{3}^{+}$ | 4.42 | $M_{4}^{-}$ | 2.92 | $L_{3}^{+}$ | 2.90
|  | $M_{4}^{+}$ | 3.24 | $L_{3}^{-}$ | 3.27
|  | $M_{3}^{-}$ | 3.32 | $L_{4}^{+}$ | 3.29
|  | $M_{3}^{+}$ | 3.46 | $L_{4}^{-}$ | 3.52
|  | $M_{1}^{+}$ | 3.63 | $L_{2}^{-}$ | 3.67
|  | $M_{2}^{-}$ | 3.85 | $L_{1}^{+}$ | 3.84
|  | $M_{4}^{-}$ | 3.89 | $L_{3}^{+}$ | 3.87
|  | $M_{3}^{-}$ | 4.14 | $L_{4}^{+}$ | 4.12
|  | $M_{4}^{-}$ | 4.93 | $L_{3}^{+}$ | 4.92
|  |  |  | $L_{2}^{-}$ | 4.98
Figure 4:  Energy lowering by $3Q$ distortions. For MMM and MLL
configurations, positive displacements correspond to tri-hexagonal distortions
and negative displacements correspond to “Star of David” distortions. The
lowest energy occurs for the MLL configuration consisting of in-plane tri-
hexagonal distortions that are laterally shifted in neighboring planes. Curves
are least squares fits to sixth-degree polynomials.

To explain our observations, we performed DFT calculations of the phonon mode
frequencies in CsV3Sb5. At high temperatures, the material crystallizes in the
$P6/mmm$ space group, as shown in Fig. 3(a). Of particular note in the first
Brillouin zone is the $M$ point with $k_{z}=0$, which connects saddle points
in the kagome electronic band structure and is associated with an ostensible
Peierls-like nesting-driven instability [12], and the $L$ point with
$k_{z}=\pi/c$, as shown in Fig. 3(b). Table 1 lists all calculated $\Gamma$-,
$M$-, and $L$-point phonon frequencies and irreducible representations up to 5
THz, and Fig. 3(c) shows the full phonon dispersion relation. Two unstable
modes are predicted, one at the $M$ point and one at the $L$ point, with the
latter having a slightly larger imaginary frequency. A symmetry analysis
indicates that the irreducible representations of the modes are $M_{1}^{+}$and
$L_{2}^{-}$, respectively, using the notation of Miller and Love [35].

As shown previously, a simultaneous condensation of the three equivalent
$M$-point modes leads to a fully symmetric breathing mode with two possible
configurations: a “Star of David” distortion (breathing out), or a tri-
hexagonal—also called “inverse Star of David”—distortion (breathing in) [25].
These two distortions, which we label MMM, are only one possible type of $3Q$
order, however. Other possible combinations of $M$\- and $L$-point modes
exist, including MML, MLL, and LLL. To determine which of these four
possibilities represents the most likely ground state in CsV3Sb5, we
calculated their energies as a function of atomic displacement. For these
calculations, we assume a frozen phonon model and do not further relax the
unit cell (i.e. nonlinear phonon couplings are ignored). As Fig. 4 shows, the
lowest energy state corresponds to the tri-hexagonal MLL distortion and not
the MMM distortion as widely believed. This ground state distortion,
consisting of one $M$\- and two $L$-point modes, can be described as in-plane
tri-hexagonal distortions with a lateral shift of one lattice vector between
neighboring kagome planes. Further support for this conclusion comes from our
earlier experimental determination that the CDW phase transition is first-
order, which immediately eliminates the odd-parity MML and LLL distortions.

Now that the true CDW structure of CsV3Sb5 has been determined, we can revisit
our coherent phonon spectroscopy data. Firstly, the $\gamma$ mode at 4.1 THz,
which is detected at all temperatures, must be a fully symmetric mode. Sure
enough, the DFT calculations find only a single $\Gamma_{1}^{+}$ ($A_{1g}$)
mode which exactly matches the experimental frequency. This mode, illustrated
in Fig. 3(d), consists of coherent motion of the out-of-plane antimony atoms
along the $c$ axis towards/away from the kagome planes. The excellent
agreement between experiment and theory confirms the accuracy of our DFT
calculations. We next address the $\alpha$ mode detected at 1.3 THz, which
appears abruptly at $T_{\mathrm{CDW}}$ and therefore represents a phonon that
becomes Raman active in the CDW phase. This indicates that the $\alpha$ mode
must have either $M_{1}^{+}$ or $L_{2}^{-}$ symmetry. Indeed, there is an
$L_{2}^{-}$ phonon with a frequency of 1.27 THz, which is illustrated in Fig.
3(e) and involves cesium atom motion along the $c$ axis. The observation of
this mode therefore serves as strong experimental confirmation that the CDW
order in CsV3Sb5 is of the MLL type and modulated along the $c$ axis. We also
note that an $M_{1}^{+}$ phonon is predicted at 3.63 THz and a second
$L_{2}^{-}$ phonon is predicted at 3.67 THz. Weak spectral intensity is
discernible in the Fourier map below $T_{\mathrm{CDW}}$ near these
frequencies, suggesting that these modes also make an appearance in the data,
albeit with a weaker amplitude.

Finally, we examine the $\beta$ resonance, which extends over a broad spectral
range centered at 3.1 THz. This is an unusual spectral feature because it
seems to appear at $T^{*}\approx 60$ K, well below $T_{\mathrm{CDW}}$. For
this reason, we believe this resonance is related to the uniaxial order
observed below this temperature by scanning tunneling microscopy experiments
[27, 24, 28]. A $1Q$ uniaxial ordering would break the $C_{6}$ rotational
symmetry of the crystal and result in a myriad of newly Raman-active modes
excitable by a pump pulse, several with frequencies near 3.1 THz. In fact, MLL
order already breaks the $C_{6}$ rotational symmetry (whereas MMM order does
not), and we speculate that this may offer a natural explanation for the
uniaxial order reported in this material system. We further speculate that the
lower onset temperature $T^{*}$ of the $1Q$ order may be related to an order-
disorder crossover associated with $c$-axis coherence, as there are three
equivalent MLL configurations which are distinguished only by the direction of
the lateral plane-to-plane shift of the tri-hexagonal CDW pattern. While these
ideas appear promising, further experimental and theoretical work will be
needed to fully solve the $1Q$ puzzle in CsV3Sb5.

In conclusion, we used coherent phonon spectroscopy experiments in conjunction
with first-principles DFT calculations to investigate the unusual CDW order in
CsV3Sb5. We uncovered a first-order $3Q$ phase transition characterized by a
simultaneous condensation of one $M$\- and two $L$-point phonons, which we
call “MLL” CDW order. This ordering, distinct from the more commonly discussed
“MMM” ordering, involves interlayer modulation of the CDW along the $c$ axis
and may offer a natural explanation for the uniaxial order observed at lower
temperatures. We have only considered $k_{z}=\pi/c$ ($L$ point) wave vectors
in this work, but it is possible that the true ground state of the system
could have a longer wavelength modulation, as recently suggested by x-ray
diffraction experiments [31]. A $4c$ wavelength, for example, could arise from
a chiral screw-axis CDW in which the plane-to-plane shift direction
continuously rotates. While many puzzles still surround the $A$V3Sb5
materials, our results represent a large step forward in understanding the CDW
order in CsV3Sb5.

We would like to thank Sam Teicher and Hengxin Tan for helpful discussions.
This work was supported by the National Science Foundation (NSF) through
Enabling Quantum Leap: Convergent Accelerated Discovery Foundries for Quantum
Materials Science, Engineering, and Information (Q-AMASE-i): Quantum Foundry
at UC Santa Barbara (DMR-1906325). Use was made of computational facilities
purchased with funds from the NSF (CNS-1725797) and administered by the Center
for Scientific Computing (CSC). The CSC is supported by the California
NanoSystems Institute and the Materials Research Science and Engineering
Center (MRSEC; NSF DMR-1720256) at UC Santa Barbara. L.H. acknowledges support
from the Roy T. Eddleman Center for Quantum Innovation at UC Santa Barbara.
B.R.O. acknowledges support from the California NanoSystems Institute through
the Elings Fellowship program.

## References

  * [1] W.-H. Ko, P. A. Lee, X.-G. Wen. Phys. Rev. B 79, 214502 (2009).
  * [2] H.-M. Guo, N. Franz. Phys. Rev. B 80, 113102 (2009).
  * [3] A. O’Brien, F. Pollmann, P. Fulde. Phys. Rev. B 81, 235115 (2010).
  * [4] J. Wen, A. Ruegg, C.-C. J. Wang, G. A. Fiete. Phys. Rev. B 82, 075125 (2010).
  * [5] S.-L. Yu, J.-X. Li. Phys. Rev. B 85, 144402 (2012).
  * [6] W.-S. Wang, Z.-Z. Li, Y.-Y. Xiang, Q.-H. Wang. Phys. Rev. B 87, 115135 (2013).
  * [7] M. L. Kiesel, C. Platt, R. Thomale. Phys. Rev. Lett. 110, 126405 (2013).
  * [8] I. I. Mazin, H. O. Jeschke, F. Lechermann, H. Lee, M. Fink, R. Thomale, R. Valenti. Nat. Commun. 5, 4261 (2014).
  * [9] L. Ye, M. Kang, J. Liu, F. von Cube, C. R. Wicker, T. Suzuki, C. Jozwiak, A. Bostwick, E. Rotenberg, D. C. Bell, L. Fu, R. Comin, J. G. Checkelsky. Nature 555, 638 (2018).
  * [10] M. Kang, L. Ye, S. Fang, J.-S. You, A. Levitan, M. Han, J. I. Facio, C. Jozwiak, A. Bostwick, E. Rotenberg, M. K. Chan, R. D. McDonald, D. Graf, K. Kaznatcheev, E. Vescovo, D. C. Bell, E. Kaxiras, J. van den Brink, M. Richter, M. P. Ghimire, J. G. Checkelsky, R. Comin. Nat. Mater. 19, 163 (2020).
  * [11] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, E. Ertekin, T. M. McQueen, E. S. Toberer. Phys. Rev. Mater. 3, 094407 (2019).
  * [12] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, S. D. Wilson. Phys. Rev. Lett. 125, 247002 (2020).
  * [13] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. L. Teicher, R. Seshadri, S. D. Wilson. Phys. Rev. Mater. 5, 034801 (2021).
  * [14] Q. Yin, Z. Tu, C. Gong, Y. Fu, S. Yan, H. Lei. Chin. Phys. Lett. 38, 037403 (2021).
  * [15] E. M. Kenney, B. R. Ortiz, C. Wang, S. D. Wilson, M. J. Graf. J. Phys.: Condens. Matter X, XXXXX (2021).
  * [16] S.-Y. Yang, Y. Wang, B. R. Ortiz, D. Liu, J. Gayles, E. Derunova, R. Gonzalez-Hernandez, L. Smejkal, Y. Chen, S. S. P. Parkin, S. D. Wilson, E. S. Toberer, T. McQueen, M. B. Ali. Sci. Adv. 6, eabb6003 (2020).
  * [17] F. H. Yu, T. Wu, Z. Y. Wang, B. Lei, W. Z. Zhuo, J. J. Ying, X. H. Chen. arXiv:2102.10987.
  * [18] Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, J. He, X. Liu, S. S. Zhang, G. Chang, I. Belopolski, Q. Zhang, M. S. Hossain, T. A. Cochran, D. Multer, M. Litskevich, Z.-J. Cheng, X. P. Yang, Z. Guguchia, G. Xu, Z. Wang, T. Neupert, S. D. Wilson, M. Z. Hasan. arXiv:2012.15709.
  * [19] C. C. Zhao, L. S. Wang, W. Xia, Q. W. Yin, J. M. Ni, Y. Y. Huang, C. P. Tu, Z. C. Tao, Z, J. Tu, C. S. Gong, H. C. Lei, Y. F. Guo, X. F. Yang, S. Y. Li. arXiv:2102.08356.
  * [20] K. Y. Chen, N. N. Wang, Q. W. Yin, Z. J. Tu, C. S. Gong, J. P. Sun, H. C. Lei, Y. Uwatoko, J.-G. Cheng. arXiv:2102.09328.
  * [21] Z. Zhang, Z. Chen, Y. Zhou, Y. Yuan, S. Wang, L. Zhang, X. Zhu, Y. Zhou, X. Chen, J. Zhou, Z. Yang. arXiv:2103.12507.
  * [22] X. Chen, X. Zhan, X. Wang, J. Deng, X. Liu, X. Chen, J. Guo, X. Chen. arXiv:2103.13759.
  * [23] Y. Wang, S.-Y. Yang, P. K. Sivakumar, B. R. Ortiz, S. M. L. Teicher, H. Wu, A. K. Srivastava, C. Garg, D. Liu, S. S. P. Parkin, E. S. Toberer, T. McQueen, S. D. Wilson, M. N. Ali, arXiv:2012.05898.
  * [24] Z. Liang, X. Hou, W. Ma, F. Zhang, P. Wu, Z. Zhang, F. Yu, J.-J. Ying, K. Jiang, L. Shan, Z. Wang, X.-H. Chen. arXiv:2103.04760.
  * [25] H. Tan, Y. Liu, Z. Wang, B. Yan. arXiv:2103.06325.
  * [26] H. X. Li, T. T. Zhang, Y.-Y. Pai, C. Marvinney, A. Said, T. Yilmaz, Q. Yin, C. Gong, Z. Tu, E. Vescovo, R. G. Moore, S. Murakami, H. C. Lei, H. N. Lee, B. Lawrie, H. Miao. arXiv:2106.09769.
  * [27] H. Zhao, H. Li, B. R. Ortiz, S. M. L. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, I. Zeljkovic. arXiv:2103.03118.
  * [28] H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Yu, Q. Yin, C. Gong, Z. Tu, H. Lei, S. Ma, H. Zhang, S. Ni, H. Tan, C. Shen, X. Dong, B. Yan, Z. Wang, H.-J. Gao. arXiv:2103.09188.
  * [29] S. Ni, S. Ma, Y. Zhang, J. Yuan, H. Yang, Z. Lu, N. Wang, J. Sun, Z. Zhao, D. Li, S. Liu, H. Zhang, H. Chen, K. Jin, J. Cheng, L. Yu, F. Zhou, X. Dong, J. Hu, H.-J. Gao, Z. Zhao. arXiv:2104.00374.
  * [30] E. Uykur, B. R. Ortiz, S. D. Wilson, M. Dressel, A. A. Tsirlin. arXiv:2103.07912.
  * [31] B. R. Ortiz, S. M. L. Teicher, L. Kautzsch, P. M. Sarte, J. P. C. Ruff, R. Seshadri, S. D. Wilson. arXiv:2104.07230.
  * [32] G. Kresse, J. Furthmuller. Phys. Rev. B 54, 11169 (1996).
  * [33] A. Togo, I. Tanaka. Scr. Mater. 108, 1 (2015).
  * [34] N. P. Wells, P. M. Belden, J. R. Demers, W. T. Lotshaw. J. Appl. Phys. 116, 073506 (2014).
  * [35] A. P. Cracknell, B. L. Davies, S. C. Miller, W. F. Love. Kronecker Product Tables. Vol. 1. General Introduction and Tables of Irreducible Representations of Space Groups (IFI/Plenum, New York, 1979).


# Fine-Grained Element Identification in Complaint Text of Internet Fraud

Tong Liu1111Tong Liu and Siyuan Wang have equal contribution., Siyuan
Wang1111Tong Liu and Siyuan Wang have equal contribution., Jingchao Fu1, Lei
Chen1, Zhongyu Wei1222Corresponding author.,
Yaqi Liu2, Heng Ye2, Liaosa Xu2, Weiqiang Wang2, Xuanjing Huang1 1School of
Data Science, Fudan University, China; 2Ant Group, China
<EMAIL_ADDRESS>fuuuuugcn<EMAIL_ADDRESS>
(2021)

###### Abstract.

Existing system dealing with online complaint provides a final decision
without explanations. We propose to analyse the complaint text of internet
fraud in a fine-grained manner. Considering the complaint text includes
multiple clauses with various functions, we propose to identify the role of
each clause and classify them into different types of fraud element. We
construct a large labeled dataset originated from a real finance service
platform. We build an element identification model on top of BERT and propose
additional two modules to utilize the context of complaint text for better
element label classification, namely, global context encoder and label
refiner. Experimental results show the effectiveness of our model.

complaint text mining, fraud element identification, dataset

††journalyear: 2021††copyright: acmcopyright††conference: Proceedings of the
30th ACM International Conference on Information and Knowledge Management;
November 1–5, 2021; Virtual Event, QLD, Australia††booktitle: Proceedings of
the 30th ACM International Conference on Information and Knowledge Management
(CIKM ’21), November 1–5, 2021, Virtual Event, QLD, Australia††price:
15.00††doi: 10.1145/3459637.3482108††isbn:
978-1-4503-8446-9/21/11††submissionid: 2498††ccs: Applied computing Secure
online transactions

##  1\. Introduction

In the field of e-commercial business, various payment platforms provide an
easy way of capital transferring but also exposes huge threat of Internet
fraud. Every year, financial companies receive thousands of fraud complaints,
from imprudence remittance to mendacious business contract. Although complaint
text is usually full of sore and loss, it is valuable for organizations to
understand incident mechanism and avoid potential risks (Brennan et al., 2017;
Ferri, 2018). Handling fraud complaints is of necessity to resolve victims’
burning issues and recognizing fraud elements enables financial organizations
and law enforcement agencies to detect security flaws and put an end to unfair
and misleading swindle practices.

However, processing numerous complaint texts throws huge challenge for
Internet finance companies. Different from credit and trading record, fraud
complaints are composed of informal and unstructured text that takes
painstaking efforts to understand and analyze. Besides, the inspect outcome is
usually short (accept further investigation or not), which might make victims
indignant for they cannot figure out why their loss never comes back.

Figure 1.  Ideal pipelines of fraud complaints processing.

Figure 1 introduces a more well-grounded scheme for complaint text processing
and fraud judgement. When a user realizes being swindled and submits
complaints to the transfer platform, the third-party payment service provider
gets involved in and segments the complaint text into pieces so as to inspect
fraud elements inside. The text on the right briefly shows the three
fundamental elements to ascertain a fraud: (1) the fraudster deliberately
fabricate an illusory circumstance to ask for money, (2) the remittance is
material and (3) the injured party realizes the falsehood afterwards and asks
for compensation. With fraud elements automatically detected, it is more
convenient for inspectors to check whether the statement is valid and more
interpretable for users to comprehend the inspect results. Once the fraud case
is established, the third-party and legal institution will follow the hints of
fraud elements to investigate the transaction history and then present a
fairer arbitrament.

To this end, we propose a novel task to identify fraud elements in a complaint
paragraph. We design an annotation criterion, split paragraphs into clauses
and ask undergraduates to annotate each clause. Moreover, we analyze clause
distribution in complaint paragraphs and explore the connection between
successional clauses while finding that position and global coherence are
influential for identifying clause role. Therefore, we propose a hierarchical
architecture which integrate context information to obtain more accurate
prediction of clauses. Our contributions are of three-folds:

  * •

We propose a novel task and construct a dataset for fine-grained fraud element
identification on clause level to further analyse internet fraud issues.

  * •

We formulate the task as a form of sequence labeling as plentiful analysis
indicates that position and relation are significant features for complaint
clauses.

  * •

We build our model on top of BERT with a global context encoder to capture the
textual context and a label-refining mechanism to utilize the label context.
Experiment results on the dataset show the effectiveness of our proposed
model.

##  2\. Dataset

###  2.1. Annotation Framework

With the support of an Internet finance service company, we collect a
considerable Chinese complaint corpus with 7 categories of fraud elements:
content fabrication (CF), identity fabrication (IF), remittance excuse (RE),
contact platform (CP), fraud realization(FR), user demand (UD) and non-
fraudulent statement (NONE).

The annotation process consists of three steps. First, we split each complaint
paragraph into clauses according to Chinese punctuation marks (comma,
semicolon and space). Second, each clause is assigned to 2 undergraduates to
acquire a unique label of fraud element. If any two annotators disagree with
each other, the instance will be checked by a third annotator following the
majority rule. If any two cannot reach an agreement, the instance will be
discarded. Finally, we calculate the Cohen’s kappa coefficient between two
annotators to assess annotation quality and get an averaged kappa coefficient
of 78%. Overall, we construct a fraud complaint dataset containing 41,103
paragraphs and 197,878 labeled clauses.

###  2.2. Dataset Analysis

Different from existing complaint-related work (Tong et al., 2018; HaCohen-
Kerner et al., 2019; Filgueiras et al., 2019) that regards every single
complaint text as independent, our corpus contains a hierarchical relationship
between clauses and paragraphs. Therefore, except for categorical statistics,
we further explore the distribution and relation of clauses in specific
paragraph.

#### Categorical Statistics

Table 1 shows the statistics of each category. The proportion of each category
is quite uneven which makes the task challenging. Non-fraudulent statements
account for a large proportion since clauses exclusive of a specific fraud
element will be regard as non-fraudulent. The elements of content fabrication,
remittance excuse and fraud realization have large amounts which are the
critical elements for a fraud case. The number of identity fabrication
elements is relatively small, but as an important supplement of factuality
modification, they are essential for companies and legal institutions to
locate credulous people and scam artists. Although the elements of contact
platform and user demand only take up a minor space, they are indicative for
the third-party companies to determine whether they have the authority of
intervention and stimulate them to provide better customer service.

Table 1. Statistics of the fraud complaint dataset.

Statistics | CF | IF | RE | CP | FR | UD | NONE
---|---|---|---|---|---|---|---
# of clauses | 39,207 | 2,739 | 19,546 | 7,882 | 35,289 | 2,608 | 90,607
Proportion | 19.81% | 1.38% | 9.88% | 3.98% | 17.83% | 1.32% | 45.79%
Avg. length of clauses | 12.0 | 12.0 | 11.4 | 10.5 | 9.9 | 10.0 | 8.5
Vocabulary size | 12,646 | 2,337 | 6,041 | 3,597 | 7,382 | 1,595 | 17,067
Clause novelty | 0.323 | 0.853 | 0.309 | 0.456 | 0.209 | 0.612 | 0.188

Besides, we measure the semantic richness of each clause by calculating their
average length. As shown in Table 1, the elements of fabricating content and
identity are longer which confirms that they contain more abundant information
to depict the fraud action. Non-fraudulent statements are the shortest because
they are more casual and untargeted. We also compute vocabulary size for each
category and divide the size by the number of clauses to obtain clause
novelty. Since identity fabrication and user demand are barely mentioned,
their repetitive use of words is quite small. For frequently emerged elements
such as content fabrication and fraud realization, their novelty difference
reflects that fraudsters practice various deception, but most of the cheated
are suffering alike.

Figure 2.  Distribution of fraud elements in different stages of complaint
paragraph.

#### Positional Distribution

Since clauses in the same paragraph are mutually complementary, it is
worthwhile to investigate how fraud elements distribute among complaint
paragraphs. Given a paragraph, we record the serial number of each clause and
divide by the sequence length to obtain the relative position. Further we
segment the whole paragraph into 5 equidistant stages. The proportion of
different stages for each fraud element is shown in Figure 2. The elements of
content fabrication, identity fabrication and contact platform are more likely
to appear in the early stage as they come straight to describe the origin,
development and transition of the fraud case. The towering bars of fraud
realization and user demand in the rear stages illustrate victims tend to
state outcomes and make requests at the end. Non-fraudulent statements are
quite evenly dispersed in different stages for they have the effect of
lubrication to make the whole paragraph clear and coherent.

Figure 3.  Probabilities of successional fraud element pairs: (a) the original
result; (b) the balanced result.

#### Ordinal Relation

We investigate all the possible combinations of successional element pair
(Mladenic and Grobelnik, 1998). Since non-fraudulent statements lack of
explicit semantics, we filter them out. For each of other categories, we
calculate the proportion of the adjacent parent. Similarly, we remove the
effect of class imbalance by dividing the prior possibility of the previous
fraud element before computing the proportion. The original and balanced
results are shown in Figure 3. Before relieving the imbalance effect, it is
inevitable that content fabrication is more likely to appear before other
elements due to its vast existence. Even so, the user demand is less likely to
follow content fabrication which means an absence of consistency between them.
However, with a nearly same proportion, the fraud realization is less possible
to become the prerequisite except for user demand. After normalizing with
prior distribution, the elements of content fabrication, identity fabrication
and contact platform are the most possible previous statements. Interestingly,
both (a) and (b) in Figure 3 have a deep-colored diagonal indicating that
clauses frequently succeed the previous of the same type and there exists a
strong semantic transitivity between consecutive clauses.

With all the findings, we claim that the position and relation are significant
attributes for complaint clauses, and it is reasonable to treat fraud element
identification as a sequence labeling task.

##  3\. Task and Model

This section presents the task of fraud element identification in complaint
text and our classification model. The overall framework is shown in Figure 4,
which is based on a pre-trained BERT model (Devlin et al., 2018). It consists
of a Local Clause Encoder, a Global Context Encoder and a Label Refiner. The
Global Context Encoder is designed to capture the textual context across the
complaint text and Label Refiner is proposed to utilize label context.

Figure 4.  The overall architecture of our proposed model.

###  3.1. Task Definition

We first introduce some notations in our task:

  1. -

$x=(x_{1},...,x_{n})$: an complaint text with $n$ clauses, where $x_{i}$ is
the $i$-th clause;

  2. -

$x_{i}=(x_{i,1},...,x_{i,m})$: a clause with $m$ tokens, where $x_{i,j}$ is
the $j$-th token in $x_{i}$;

  3. -

$y_{i}$: the fraud element of the $i$-th clause $x_{i}$ in the complaint text
$x$.

We formulate the task as a sequence-level classification task that assigns a
element type $y_{i}$ for each clause $x_{i}$ in the complaint text $x$.

###  3.2. Our Model

In our model, the Local Clause Encoder first utilizes BERT to encode each
clause in a complaint text $x$ and take the representations of [CLS] symbol as
the clause representations $(h_{1},...,h_{n})$. Then we introduce the Global
Context Encoder and the Label Refiner.

#### Global Context Encoder

Considering a complaint text contains several fine-grained clauses, we
introduce a global context encoder (Niu et al., 2019) to capture the global
textual context information $c$ to help identify the element of clauses in the
complaint text $x$. We concatenate all clauses in the complaint text
$x^{c}=concat([x_{1},...,x_{n}])$ and take the embeddings of them as the input
of the global encoder, and use a bi-directional GRU (Huang et al., 2015) to
encode the whole text sequence. Then a mean pooling is conducted over encoding
states of all time steps to get the textual context representation as Eq.(1),
where $E_{w}$ is the word embedding matrix.

(1) $\displaystyle h^{c}_{i}$ $\displaystyle={\rm
BiGRU}(E_{w}x^{c}_{i},h^{c}_{i-1})$ $\displaystyle c$ $\displaystyle={\rm
mean}(h^{c}_{1},...,h^{c}_{i},...)$

#### Label Refiner

After encoding the clauses and getting the textual context representation of
the complaint text, we identify the fraud element of each clause. Instead of
directly feeding [CLS] representation from BERT into a output layer, we take
the clause representations $(h_{1},...,h_{n})$ in the complaint text $x$ as
the input and employ a BiGRU to learn the dependencies between the clauses and
get their hidden states as $(s^{1}_{1},s^{1}_{2},...,s^{1}_{n})$. To be aware
of the global textual context information, we concatenate the state
$s^{1}_{i}$ of clause $x_{i}$ with the textual context representation $c$
(Tanaka et al., 2019) to get the classification probability through a feed-
forward network as Eq.(2).

(2) $\displaystyle P(y_{i}|x_{i},x)$
$\displaystyle=Softmax(W_{1}[s^{1}_{i};c])$

Then we propose a _Label-Refining_ mechanism to utilize the label context
information to refine the identification results. As shown in Figure 4, after
the first BiGRU layer for classification, we take the concatenation of each
clause representation $h_{i}$ and the previous classification probability
$P(y_{i}|x_{i},x)$ as the input and employ another BiGRU to again encode the
clauses sequence in the complaint text $x$. Then we get the refined
probability $P(y^{\prime}_{i}|x_{i},x)$ as Eq.(3).

(3) $\displaystyle s^{2}_{i}$ $\displaystyle={\rm
BiGRU}([P(y_{i}|x_{i},x);h^{i}],s^{2}_{i-1})$ $\displaystyle P$
$\displaystyle(y^{\prime}_{i}|x_{i},x)=Softmax(W_{2}[s^{2}_{i};c])$

##  4\. Experiments

###  4.1. Experimental Setup

We randomly split the dataset into training, validation and test set by the
proportion of 8:1:1. The numbers of complaint paragraphs and clauses in each
data split are (32882, 4110, 4111) and (158521, 19731, 19626), respectively.

We adopt BERT-base, Chinese model as backbone. The global encoder is a two-
layer GRU while the GRUs in label refiner are both of one-layer, and the size
of hidden units in all GRUs is set as 256. Mini-batch of size 32 is taken and
the dropout rate is 0.3. We first fine-tune BERT for the classification task
as the base model for 4 epochs using AdamW optimizer with learning rate of
2e-5. Then we load the fine-tuned BERT as the local clause encoder and freeze
its parameters, and then train our model for another 10 epochs using Adam
optimizer with learning rate 2e-4.

###  4.2. Overall Performance

We compare _Our Model_ with some baseline and state-of-the-art models,
including _SVM_ (Suykens and Vandewalle, 1999), _BiGRU_ (Liu and Guo, 2019),
_BERT_ (Devlin et al., 2018), _BERT-wwm-ext_ (Cui et al., 2019), _BERT+BiGRU_
and _RoBERTa_ (Liu et al., 2019). _Our Model-GC_ and _Our Model-LR_ are
ablation tests of our model, without the global context encoder and the
_label-refining_ mechanism, respectively.

Table 2. Evaluation results for element identification in complaint text of
different models.

Model | Accuracy(%) | Precision(%) | Recall(%) | F1-score(%)
---|---|---|---|---
_SVM_ | 77.22 | 79.38 | 71.41 | 75.19
_BiGRU_ | 80.51 | 79.46 | 77.31 | 78.37
_BERT_ | 83.53 | 83.66 | 80.64 | 82.13
_BERT-wwm-ext_ | 83.27 | 82.80 | 81.63 | 82.21
_BERT+BiGRU_ | 82.70 | 82.38 | 80.44 | 81.40
_RoBERTa_ | 83.40 | 82.39 | 81.97 | 82.18
_Our Model-LR_ | 84.11 | 83.28 | 82.31 | 82.79
_Our Model-GC_ | 84.30 | 83.73 | 82.35 | 83.04
_Our Model_ | 84.47 | 84.11 | 82.30 | 83.19

We take Accuracy, Precision, Recall, and F1-score as evaluation metrics, and
results are shown in the Table 2. We have several findings as follows: (1)
_Our Model_ outperforms others in terms of all metrics except recall. This
indicates the effectiveness of our model considering both textual context and
label context information. (2) _BiGRU_ performs worse than all other BERT-
based models, which demonstrate the effectiveness of pre-trained models. (3)
_BERT_ , _BERT-wwm-ext_ and _RoBERTa_ achieve similar performance, which means
that _BERT_ is enough to capture clause representation for classification.
Thus we adopt _BERT_ as our base model. (4) The improvements of _Our Model_
compared to _Our Model-LR_ and _Our Model-GC_ respectively reveal the
effectiveness of the global context encoder and _label-refining_ mechanism for
our element identification task.

###  4.3. Further Analysis

#### Results of different fraud element

We investigate the performance of different fraud element types of our model
in terms of Precision, Recall and F1-score in Table 3. We can see that for
clauses of Fraud Realization, User Demand and Non-Fraudulent Statement, our
model achieves a good performance. From the Figure 2 which shows clauses of
Fraud Realization and User Demand concentrated in the last stage of the
complaint text, our model is easier to capture this pattern so that it
performs better for these two elements. As for Non-Fraudulent Statement, we
assume the good performance is owing to its large proportion in the dataset.

Table 3. Evaluation results for complaint text classification of different
fraud element types.

Fraud Element | Precision(%) | Recall(%) | F1-score(%)
---|---|---|---
CF | 82.68 | 82.33 | 82.50
IF | 83.90 | 82.49 | 83.19
RE | 75.85 | 63.74 | 69.27
CP | 85.90 | 83.23 | 84.54
FR | 86.07 | 89.16 | 87.59
UD | 88.43 | 86.99 | 87.70
NONE | 85.92 | 88.15 | 87.02

#### Error Analysis

We compute the confusion matrix to analyse which pairs of fraud elements are
more confused to be classified. We also compare the confusion matrices of _Our
Model-LR_ and _Our Model_ to see which elements are improved by the label
refining mechanism.

Figure 5.  The confusion matrices for classification results of _Our Model-LR_
(a) and _Our Model_ (b).

As shown in Figure 5, we can see that clauses of Non-Fraudulent Statement are
easier to be confused with all other elements, because it makes up the
majority of the imbalanced dataset. And the pairs of (Identity Fabrication,
Content Fabrication), (Remittance Excuse, Content Fabrication) and (Remittance
Excuse, Fraud Realization) are likely to be classified incorrectly. From
comparison between the diagonal elements of matrix (a) and (b), we find that
the classification of clauses in Remittance Excuse, Non-Fraudulent Statement
and Identity Fabrication are indeed improved, and the corresponding confused
numbers are decreased, for example, (Remittance Excuse, Content Fabrication)
is reduced from 142 to 126.

##  5\. Related Work

Recently many complaint classification tasks have been widely studied (Kos et
al., 2017). Researchers categorize complaints in order to figure out the
complaint reasons (Tong et al., 2018), identify downstream companies (HaCohen-
Kerner et al., 2019), and explore user demands (Joung et al., 2019).
Filgueiras et al. (2019) analyze the complaints about food safety and economic
surveillance to help government structuralize complaint letters and decide the
fine-grained department ultimately responsible for the complaint. In our work,
rather than assigning labels to the whole complaint text, we map each
component of the paragraph and list all the fraud elements to provide more
interpretable judicial outcomes.

Traditional methods for text classification utilize manually crafted features,
and employ machine learning algorithms such as Naive Bayes (McCallum et al.,
1998) and Support Vector Machines (Cortes and Vapnik, 1995) to obtain category
bounds. With the popularity of deep learning, RNN, CNN (Kim, 2014), and pre-
trained language models(Howard and Ruder, 2018; Radford et al., 2018; Devlin
et al., 2018; Liu et al., 2019) refresh the performance of text
classification. In this paper, we build our model on top of BERT and propose a
fine-grained classification schema for element identification using textual
and label context information.

##  6\. Conclusion

In this paper, we aim at analyzing complaint text of internet fraud in a fine-
grained level. We first propose an annotation scheme to distinguish various
types of fraud elements and construct a dataset as benchmark. We build a
classification model on top of BERT and use context information of complaint
text to better identify element types. Experiment results confirm the
effectiveness of our model. In the future, we will explore how to utilize
domain knowledge for complaint text modeling. We will also utilize the fine-
grained element type information for better complaint text classification.

## References

  * (1)
  * Brennan et al. (2017) Carol Brennan, Tania Sourdin, Jane Williams, Naomi Burstyner, and Chris Gill. 2017.  Consumer vulnerability and complaint handling: Challenges, opportunities and dispute system design.  _International journal of consumer studies_ 41, 6 (2017), 638–646.
  * Cortes and Vapnik (1995) Corinna Cortes and Vladimir Vapnik. 1995.  Support-vector networks.  _Machine learning_ 20, 3 (1995), 273–297.
  * Cui et al. (2019) Yiming Cui, Wanxiang Che, Ting Liu, Bing Qin, Ziqing Yang, Shijin Wang, and Guoping Hu. 2019.  Pre-training with whole word masking for chinese bert.  _arXiv preprint arXiv:1906.08101_ (2019).
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ (2018).
  * Ferri (2018) Eric Ferri. 2018\.  The Evolving Practice Of Complaint Management.  _Bloomberg Law_ (2018), 1–8.
  * Filgueiras et al. (2019) Joao Filgueiras, Luís Barbosa, Gil Rocha, Henrique Lopes Cardoso, Luís Paulo Reis, Joao Pedro Machado, and Ana Maria Oliveira. 2019\.  Complaint Analysis and Classification for Economic and Food Safety. In _Proceedings of the Second Workshop on Economics and Natural Language Processing_. 51–60.
  * HaCohen-Kerner et al. (2019) Yaakov HaCohen-Kerner, Rakefet Dilmon, Maor Hone, and Matanya Aharon Ben-Basan. 2019.  Automatic classification of complaint letters according to service provider categories.  _Information Processing & Management_ 56, 6 (2019), 102102\.
  * Howard and Ruder (2018) Jeremy Howard and Sebastian Ruder. 2018.  Universal language model fine-tuning for text classification.  _arXiv preprint arXiv:1801.06146_ (2018).
  * Huang et al. (2015) Zhiheng Huang, Wei Xu, and Kai Yu. 2015.  Bidirectional LSTM-CRF models for sequence tagging.  _arXiv preprint arXiv:1508.01991_ (2015).
  * Joung et al. (2019) Junegak Joung, Kiwook Jung, Sanghyun Ko, and Kwangsoo Kim. 2019\.  Customer complaints analysis using text mining and outcome-driven innovation method for market-oriented product development.  _Sustainability_ 11, 1 (2019), 40.
  * Kim (2014) Yoon Kim. 2014.  Convolutional neural networks for sentence classification.  _arXiv preprint arXiv:1408.5882_ (2014).
  * Kos et al. (2017) William Kos, MP Schraagen, MJS Brinkhuis, and FJ Bex. 2017\.  Classification in a Skewed Online Trade Fraud Complaint Corpus. In _Preproceedings of the 29th Benelux Conference on Artificial Intelligence November 8–9, 2017 in Groningen, The Netherlands_. 172–183.
  * Liu and Guo (2019) Gang Liu and Jiabao Guo. 2019.  Bidirectional LSTM with attention mechanism and convolutional layer for text classification.  _Neurocomputing_ 337 (2019), 325–338.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_ (2019).
  * McCallum et al. (1998) Andrew McCallum, Kamal Nigam, et al. 1998\.  A comparison of event models for naive bayes text classification. In _AAAI-98 workshop on learning for text categorization_ , Vol. 752. Citeseer, 41–48.
  * Mladenic and Grobelnik (1998) Dunja Mladenic and Marko Grobelnik. 1998.  Word sequences as features in text-learning. In _In Proceedings of the 17th Electrotechnical and Computer Science Conference (ERK98_. Citeseer.
  * Niu et al. (2019) Guocheng Niu, Hengru Xu, Bolei He, Xinyan Xiao, Hua Wu, and GAO Sheng. 2019\.  Enhancing Local Feature Extraction with Global Representation for Neural Text Classification. In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_. 496–506.
  * Radford et al. (2018) Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever. 2018.  Improving language understanding by generative pre-training.
  * Suykens and Vandewalle (1999) Johan AK Suykens and Joos Vandewalle. 1999.  Least squares support vector machine classifiers.  _Neural processing letters_ 9, 3 (1999), 293–300.
  * Tanaka et al. (2019) Koji Tanaka, Junya Takayama, and Yuki Arase. 2019\.  Dialogue-Act Prediction of Future Responses Based on Conversation History. In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics: Student Research Workshop_. 197–202.
  * Tong et al. (2018) Xuesong Tong, Bin Wu, Shuyang Wang, and Jinna Lv. 2018\.  A complaint text classification model based on character-level convolutional network. In _2018 IEEE 9th International Conference on Software Engineering and Service Science (ICSESS)_. IEEE, 507–511.


# Rao–Blackwellization in the MCMC era

Christian P. Robert1,2,3 and Gareth O. Roberts1111 The work of the first
author was partly supported in part by the French government under management
of Agence Nationale de la Recherche as part of the “Blanc SIMI 1” program,
reference ANR-18-CE40-0034 and in part by the French government under
management of Agence Nationale de la Recherche as part of the “Investissements
d’avenir” program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute). The
first author also acknowledges the support of l’Institut Universitaire de
France through two consecutive senior chairs.
1University of Warwick,
2Université Paris Dauphine PSL,
and 3CREST-ENSAE

###### Abstract

Rao–Blackwellization is a notion often occurring in the MCMC literature, with
possibly different meanings and connections with the original Rao–Blackwell
theorem (Rao, , 1945; Blackwell, , 1947), including a reduction of the
variance of the resulting Monte Carlo approximations. This survey reviews some
of the meanings of the term.

Keywords: Monte Carlo, simulation, Rao–Blackwellization, Metropolis-Hastings
algorithm, Gibbs sampler, importance sampling, mixtures, parallelisation,

> This paper is dedicated to Professor C.R. Rao in honour of his 100th
> birthday.

##  1 Introduction

The neologism Rao–Blackwellization222We will use the American English spelling
of the neologism as this version is more commonly used in the
literature.333Berkson,  (1955) may have been the first one to use (p.142) this
neologism. stems from the famous Rao–Blackwell theorem (Rao, , 1945;
Blackwell, , 1947), which states that replacing an estimator by its
conditional expectation given a sufficient statistic improves estimation under
any convex loss. This is a famous mathematical statistics result, both dreaded
and appreciated by our students for involving conditional expectation and for
producing a constructive improvement, respectively. While Monte Carlo
approximation techniques cannot really be classified as estimation, since they
operate over controlled simulations, rather than observations, with the
ability to increase the sample size if need be, and since there is rarely a
free and unknown parameter involved, hence almost never a corresponding notion
of sufficiency, seeking improvement in Monte Carlo approximation via partial
conditioning has nonetheless been named after this elegant theorem. As shown
in Figure 1, the use of the expression Rao–Blackwellization has considerably
increased in the 1990’s, once the foundational paper popularising MCMC
techniques refered to this technique to reduce Monte Carlo variability,

Figure 1:  Google assessment of the popularity of the names Rao–Blackwell,
Rao–Blackwellisation, and Rao–Blackwellization since the derivation of the
Rao–Blackwell theorem.

The concept indeed started in the 1990 foundational paper by Gelfand and Smith
(“foundational” as it launched the MCMC revolution, see Green et al., , 2015).
While this is not exactly what is proposed in the paper, as detailed in the
following section, it is now perceived444See e.g. the comment in the
Introduction Section of Liu et al.,  (1994). that the authors remarked that,
given a Gibbs sampler whose component $\theta_{1}$ is simulated from the
conditional distribution, $\pi(\theta_{1}|\theta_{2},x)$, the estimation of
the marginal $\pi(\theta_{1}|x)$ is improved by considering the average of the
(full) conditionals across iterations,

$\nicefrac{{1}}{{T}}\sum_{t=1}^{r}\pi(\theta_{1}|\theta_{2}^{(t)},x)$

which provides a parametric, unbiased and
$\text{O}(\nicefrac{{1}}{{\sqrt{T}}})$ estimator. Similarly, the approximation
to $\mathbb{E}[\theta_{1}|x]$ based on this representation

$\nicefrac{{1}}{{T}}\sum_{t=1}^{r}\mathbb{E}[\theta_{1}|\theta_{2}^{(t)},x]$

is using conditional expectations with lesser variances than the original
$\theta_{1}^{(t)}$ and may thus lead to a reduced variance for the estimator,
if correlation does not get into the way. (In that specific two-step sampler,
this is always the case Liu et al., , 1994.

We are thus facing the difficult classification task of separating what is
Rao–Blackwellization from what is not Rao–Blackwellization in simulation and
in particular MCMC settings.

The difficulty resides in setting the limits as

  * •

there is no clear notion of sufficiency in simulation and, further,
conditioning may increase the variance of the resulting estimator or slow down
convergence;

  * •

variance reduction and unbiasedness are not always relevant (as shown by the
infamous harmonic mean estimator, Neal, , 1999; Robert and Wraith, , 2009), as
for instance in infinite variance importance sampling (Chatterjee and
Diaconis, , 2018; Vehtari et al., , 2019);

  * •

there are (too) many forms of potential conditioning in simulation settings to
hope for a ranking (see, e.g., the techniques of partitioning, antithetic or
auxiliary variables, control variates as in Berg et al., , 2019, delayed
acceptance as in Banterle et al., , 2019; Beskos et al., 2006a, , adaptive
mixtures as in Owen and Zhou, , 2000; Cornuet et al., , 2012; Elvira et al., ,
2019, the later more closely connected to Rao–Blackwellization);

  * •

the large literature on the approximation of normalising constants and Bayes
factors (Robert and Marin, , 2008; Marin and Robert, , 2010, 2011) contains
many proposals that relate to Rao-Blackwellisation, as, e.g., through the
simulation of auxiliary samples from instrumental distributions as initiated
in Geyer,  (1993) and expanded into bridge sampling by Chopin and Robert,
(2010) and noise-constrastive estimation by Gutmann and Hyvärinen,  (2012);

  * •

in connection with the above, many versions of demarginalization such as slice
sampling (Roberts and Rosenthal, , 1999; Mira et al., , 2001) introduce
auxiliary variables that could be exploited towards bringing a variance
reduction;555This may however be seen as a perversion of Rao–Blackwellization
in that the dimension of the random variable used in the simulation is
increased, with the resulting estimate being obtained by the so-called Law of
the Unconscious Statistician.

  * •

there is no optimal solution in simulation as, mathematically, a quantity such
as an expectation is uniquely and exactly defined once the distribution is
known: if computation time is not accounted for, the exact value is the
optimal solution;

  * •

while standing outside a probabilistic framework, quasi-Monte Carlo techniques
(Liao, , 1998) can also be deemed to constitute an ultimate form of
Rao–Blackwellization, with the proposal of Kong et al.,  (2003) being an
intermediate solution;666As mentioned by the authors, the “group-averaged
estimator may be interpreted as Rao–Blackwellization given the orbit, so group
averaging cannot increase the variance” (p. 592)

but we will not cover any further these aspects here.

The rest of this review paper discusses Gibbs sampling in Section 2, other
MCMC settings in 3, particle filters and SMC in 5, and conclude in Section 6.

##  2 Gibbs sampling

Let us recall that a Gibbs sampler (Geman and Geman, , 1984) is a specific way
of building a Markov chain with stationary density $\pi(\cdot)$ through the
iterative generation from conditional densities associated with the joint
$\pi(\cdot)$. Its simplest version consists in partitioning the argument
$\theta$ into $\theta=(\theta_{1},\theta_{2})$ and generating alternatively
from $\pi_{1}(\theta_{1}|\theta_{2})$ and from
$\pi_{2}(\theta_{2}|\theta_{1})$. This binary version is sometimes called data
augmentation in reference to Tanner and Wong,  (1987), who implemented an
algorithm related to the Gibbs sampler for latent variable models.

When proposing this algorithm as a way to simulating from marginal densities
and (hence) posterior distributions, Gelfand and Smith (1990) explicitely
relate to the Rao–Blackwell theorem, as shown by the following quote777The
text has been retyped and may hence contains typos. The notations are those
introduced by Gelfand and Smith,  (1990) and used for a while in the
literature, see e.g. Spiegelhalter et al.,  (1995) with $[X\mid Y]$ denoting
the conditional density of $X$ given $Y$. The double indexation of the
sequence is explained below.

> …we consider the problem of calculating a final form of marginal density
> from the final sample produced by either the substitution or Gibbs sampling
> algorithms. Since for any estimated marginal the corresponding full
> conditional has been assumed available, efficient inference about the
> marginal should clearly be based on using this full conditional
> distribution. In the simplest case of two variables, this implies that
> $[X\mid Y]$ and the $y^{(i)}_{j}$’s $(j=1,\ldots,m)$ should be used to make
> inferences about $[X]$, rather than imputing $X^{(i)}_{j}$ $(j=1,\ldots,nm)$
> and basing inference on these $X^{(i)}_{j}$’s. Intuitively, this follows,
> because to estimate $[X]$ using the $x^{(i)}_{j}$’s requires a kernel
> density estimate. Such an estimate ignores the known form $[X\mid Y]$ that
> is mixed to obtain $[X]$. The formal argument is essentially based on the
> Rao–Blackwell theorem. We sketch a proof in the context of the density
> estimator itself. If $X$ is a continuous p-dimensional random variable,
> consider any kernel density estimator of $[X]$ based on the $X^{(i)}_{j}$’s
> (e.g., see Devroye and Györfi, 1985) evaluated at $x_{0}$:
> $\Delta^{(i)}_{x_{0}}=(1/h_{m}^{p})\sum_{j=1}^{m}K[(X_{0}-X^{(i)}_{j})/h_{m}]$,
> say, where $K$ is a bounded density on $\mathbb{R}^{p}$ and the sequence
> $\\{h_{m}\\}$ is such that as $m\rightarrow\infty$, $h_{m}\rightarrow 0$,
> whereas $mh_{m}\rightarrow\infty$. To simplify notation, set
> $Q_{m,x_{0}}(X)=(1/h_{m}^{p})K[(X-X^{(i)}_{j})/h_{m}]$ so that
> $\Delta^{(i)}_{x_{0}}=(1/m)\sum_{j=1}^{m}Q_{m,x_{0}}(X_{j}^{(i)})$. Define
> $\gamma_{x_{0}}^{i}=(1/m)\sum_{j=1}^{m}\mathbb{E}[Q_{m,x_{0}}(X)\mid
> Y^{(i)}_{j}]$. By our earlier theory, both $\Delta^{(i)}_{x_{0}}$ and
> $\gamma_{x_{0}}^{i}$ have the same expectation. By the Rao–Blackwell
> theorem, $\text{var}\,\mathbb{E}[Q_{m,x_{0}}(X)]\mid
> Y)\leq\text{var}\,Q_{m,x_{0}}(X)$, and hence
> $\text{MSE}(\gamma_{x_{0}}^{i})\leq\text{MSE}(\Delta^{(i)}_{x_{0}})$, where
> MSE denotes the mean squared error of the estimate of $[X]$.

This Section 2.6. of the paper calls for several precisions:

  * •

the simulations $x^{(i)}_{j}$ and $y^{(i)}_{j}$ are double-indexed because the
authors consider $m$ parallel and independent runs of the Gibbs sampler, $i$
being the number of iterations since the initial step, in continuation of
Tanner and Wong,  (1987),

  * •

the Rao–Blackwell argument is more specifically a conditional expectation
step,

  * •

as later noted by Geyer,  (1994), the conditioning argument is directed at
(better) approximating the entire density $[X]$, even though the authors
mention on the following page that the argument is “simpler for estimation of”
a posterior expectation,

  * •

they compare the mean squared errors of the expected density estimate rather
than the rates of convergence of a non-parametric kernel estimator(in
$n^{\nicefrac{{-1}}{{4+d}}}$) versus an unbiased parametric density estimator
(in $n^{\nicefrac{{-1}}{{2}}}$), which does not call for a
Rao–Blackwellization argument,

  * •

they do not (yet) mention “Rao–Blackwellization” as a technique,

  * •

and they do not envision (more) ergodic averages across iterations, possibly
fearing the potential impact of the correlation between the terms for a given
chain.

A more relevant step in the use of Rao–Blackwellization techniques for the
Gibbs sampler is found in Liu et al.,  (1994). This later article establishes
in particular that, for the two-step Gibbs sampler, Rao–Blackwellization
always produces a decrease in the variance of the empirical averages. This is
established in a most elegant manner by showing that each extra conditioning
(or further lag) decreases the correlation, which is always positive.888The
“always” qualification applies to every transform of the chain and to every
time lag. The proof relies on the associated notion of interleaving and
expresses the above correlation as the variance of a multiply conditioned
expectation:

$\text{cov}(h(\theta_{1}^{(0)}),h(\theta_{1}^{(n)})=\text{var}\,(\mathbb{E}(\ldots\mathbb{E}[\mathbb{E}\\{h(\theta_{1})\mid\theta_{2}\\}]\ldots)\\}\,,$

where the number of conditional expectations on the rhs is $n$. The authors
also warn that a “fast mixing scheme gains an extra factor in efficiency if
the mixture estimate can be easily computed” and give a counter-example when
Rao–Blackwellization increases the variance.999The function leading to the
counter-example is however a function of both $\theta_{1}$ and $\theta_{2}$,
which may be considered as less relevant in latent variable settings. This
counter-example is exploited in a contemporary paper by Geyer,
(1994)101010Geyer,  (1994) also points out that a similar Rao–Blackwellization
was proposed by Pearl,  (1987). where a necessary and sufficient but highly
theoretical condition is given for an improvement. As the author puts it in
his conclusion,

> The point of this article is not that Rao–Blackwellized estimators are a
> good thing or a bad thing. They may be better or worse than simple averaging
> of the functional of interest without conditioning. The point is that, when
> the autocorrelation structure of the Markov chain is taken into account, it
> is not a theorem that Rao–Blackwellized estimators are always better than
> simple averaging. Hence the name Rao–Blackwellized should be avoided,
> because it brings to mind optimality properties that these estimators do not
> really possess. Perhaps “averaging a conditional expectation” is a better
> name.

but his recommendation was not particularly popular, to judge from the
subsequent literature resorting to this denomination.

Another connection between Rao–Blackwellization and Gibbs sampling can be
found in Chib,  (1995), where his approximation to the marginal likelihood

$m(x)=\dfrac{\pi(\theta^{*})f(x|\theta^{*})}{\hat{\pi}(\theta^{*}|x)}$

is generaly based on an estimate111111Chib,  (1995) mentions this connection
(p.1314) but seems to restrict it to the two-stage Gibbs sampler. In the
earlier version known as “the candidate’s formula”, due to a Durham student
coming up with it, Besag,  (1989) points out the possibility of using an
approximation such as a Laplace approximation, rather than an MCMC
estimation.121212A question found on the statistics forum Cross-Validated
illustrates the difficulty with understanding demarginalisation and joint
simulation: “Chib suggests that we can insert the Gibbs sampling outputs of
$\mu$ into the summation [of the full conditionals]. But aren’t the outputs
obtained from Gibbs about the joint posterior $p(\mu,\phi|y)$? Why suddenly
can we use the results from joint distribution to replace the marginal
distribution?” of the posterior density using a latent (or auxiliary)
variable, as in Gelfand and Smith,  (1990),

$\hat{\pi}(\theta^{*}|x)=\nicefrac{{1}}{{T}}\sum_{t=1}^{T}\pi(\theta^{*}|x,z^{(t)})$

The stabilisation brought by this parametric approximation is notable when
compared with kernel estimates, even though it requires that the marginal
distribution on $z$ is correctly simulated (Neal, , 1999).

##  3 Markov chain Monte Carlo methods

In the more general setting of Markov chain Monte Carlo (MCMC) algorithms
(Robert and Casella, , 2004), further results characterise the improvement
brought by Rao–Blackwellization. Let us briefly recall that the concept behind
MCMC is to create a Markov sequence $\theta^{(n)})_{n}$ of dependent variables
that converge (in distribution) to the distribution of interest (also called
target). One of the most ubiquitous versions of an MCMC algorithm is the
Metropolis–Hastings algorithm (Metropolis et al., , 1953; Hastings, , 1970;
Green et al., , 2015)

One direct exploitation of the Rao–Blackwell theorem is found in McKeague and
Wefelmeyer,  (2000), who sho in particular that, when estimating the mean of
$\theta$ under the target distribution, a Rao–Blackwellized version based on
$\mathbb{E}[h(\theta^{(n)}|\theta^{(n-1)}]$ will improve the asymptotic
variance of the ordinary empirical estimator when the chain $(\theta^{(n)})$
is reversible. While the setting may appear quite restrictive, the authors
manage to recover data augmentation with a double conditional expectation
(when compared with Liu et al., , 1994) as well as reversible Gibbs and
Metropolis samplers of the Ising model. The difficulty in applying the method
resides in computing the conditional expectation, since a replacement with a
Monte Carlo approximation cancels its appeal.

Casella and Robert,  (1996) consider an altogether different form of
Rao–Blackwellization for both accept-reject and Metropolis–Hastings samples.
The core idea is to integrate out via a global conditional expectation the
Uniform variates used to accept or reject the proposed values.

A sample produced by the Metropolis–Hastings algorithm,
$\theta^{(1)},\ldots,\allowbreak\theta^{(T)}$, is in fact based on two
simulated samples, the sample of proposed values $\eta_{1},\ldots,\eta_{T}$
and the sample of decision variates $u_{1},\ldots,u_{T}$, with $\eta_{t}\sim
q(y|\theta^{(t-1)})$ and $u_{t}\sim\mathcal{U}([0,1])$. Since $\theta^{(t)}$
is equal to one of the earlier proposed values, an empirical average
associated with this sample can be written131313In order to avoid additional
notations, we assume a continuous model where all $\eta_{i}$’s are different
with probability one.

$\delta^{MH}=\nicefrac{{1}}{{T}}\sum_{t=1}^{T}h(\theta^{(t)})=\nicefrac{{1}}{{T}}\sum_{t=1}^{T}\sum_{i=1}^{t}\;\mathbb{I}_{\theta^{(t)}=\eta_{i}}\;h(\eta_{i})\,.$

Therefore, taking a conditional expectation of the above by integrating the
decision variates,

$\displaystyle\delta^{RB}$
$\displaystyle=\displaystyle{\nicefrac{{1}}{{T}}\sum_{i=1}^{T}\;h(\eta_{i})\;\mathbb{E}\left[\sum_{t=i}^{T}\;\mathbb{I}_{\theta^{(t)}=\eta_{i}}\bigg{|}\eta_{1},\ldots,\eta_{T}\right]}$
$\displaystyle=\displaystyle{\nicefrac{{1}}{{T}}\sum_{i=1}^{T}\;h(\eta_{i})\;\sum_{t=i}^{T}\;\mathbb{P}(\theta^{(t)}=\eta_{i}|\eta_{1},\ldots,\eta_{T})}\,,$

leads to an improvement of the empirical average, $\delta^{MH}$, under convex
losses.

While, for the independent Metropolis–Hastings algorithm, the conditional
probability can be obtained in closed form (see also Atchadé and Perron, ,
2005 and Jacob et al., , 2011), the general case, based on an arbitrary
proposal distribution $q(\cdot|\theta)$ is such that $\delta^{RB}$ is less
tractable but Casella and Robert,  (1996) derive a tractable recursive
expression for the weights of $h(\eta_{i})$ in $\delta^{RB}$, with complexity
of order $\mathcal{O}(T^{2})$. Follow-up papers are Perron,  (1999) and
Casella et al.,  (2004).

While again attempting at integrating out the extraneous uniform variates
exploited by the Metropolis–Hastings algorithm, Douc and Robert,  (2011)
derive another Rao–Blackwellized improvement over the regular
Metropolis–Hastings algorithm by following a different representation of
$\delta^{MH}$, using the accepted chain $(\xi_{i})_{i}$ instead of the
proposed sequence of the $\eta_{i}$’s as in Casella and Robert,  (1996). The
version based on accepted values is indeed rewritten as

$\delta^{MH}=\nicefrac{{1}}{{T}}\,\sum_{i=1}^{M}\mathfrak{n}_{i}h(\xi_{i})\,,$

where the $\xi_{i}$’s are the accepted $\eta_{j}$’s, $M$ is the number of
accepted $\eta_{j}$’s till iteration $T$, and $\mathfrak{n}_{i}$ is the number
of times $\xi_{i}$ appears in the sequence $(\theta^{(t)})_{t}$. This
representation is also exploited in Sahu and Zhigljavsky,  (1998); Gȧsemyr,
(2002); Sahu and Zhigljavsky,  (2003), and Malefaki and Iliopoulos,  (2008).
The Rao–Blackwellisation construct of Douc and Robert,  (2011) exploits the
following properties:

  1. 1.

$(\xi_{i},\mathfrak{n}_{i})_{i}$ is a Markov chain;

  2. 2.

$\xi_{i+1}$ and $\mathfrak{n}_{i}$ are independent given $\xi_{i}$;

  3. 3.

$\mathfrak{n}_{i}$ is distributed as a Geometric random variable with
probability parameter

$p(\xi_{i}):=\int\alpha(\xi_{i},\eta)\,q(\eta|\xi_{i})\,\mathrm{d}\eta\,;$ (1)

  4. 4.

$(\xi_{i})_{i}$ is a Markov chain with transition kernel
$\tilde{Q}(\xi,\mathrm{d}\eta)=\tilde{q}(\eta|\xi)\mathrm{d}\eta$ and
stationary distribution $\tilde{\pi}$ such that

$\tilde{q}(\cdot|\xi)\propto\alpha(\xi,\cdot)\,q(\cdot|\xi)\quad\mbox{and}\quad\tilde{\pi}(\cdot)\propto\pi(\cdot)p(\cdot)\,.$

Since the Metropolis–Hastings estimator $\delta^{MH}$ only involves the
$\xi_{i}$’s, i.e. the accepted $\eta_{t}$’s, an optimal weight for those
random variables is the importance weight $1/p(\xi_{i})$, leading to the
corresponding importance sampling estimator

$\delta^{IS}=\nicefrac{{1}}{{N}}\,\sum_{i=1}^{M}\dfrac{h(\xi_{i})}{p(\xi_{i})}\,,$

but this quantity is almost invariably unavailable in closed form and need be
estimated by an unbiased estimator. The geometric $\mathfrak{n}_{i}$ is the de
facto solution that is used in the original Metropolis-Hastings estimate, but
solutions with smaller variance also are available, based on the property that
(if $\alpha(\xi,\eta)$ denotes the Metropolis–Hastings acceptance probability)

$\hat{\zeta}_{i}=1+\sum_{j=1}^{\infty}\prod_{\ell\leq
j}\left\\{1-\alpha(\xi_{i},\eta_{\ell})\right\\}\qquad\eta_{\ell}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}q(\eta|\xi_{i})$

is an unbiased estimator of $1/p(\xi_{i})$ whose variance, conditional on
$\xi_{i}$, is lower than the conditional variance of $\mathfrak{n}_{i}$,
$\\{1-p(\xi_{i})\\}/p^{2}(\xi_{i})$. For practical implementation, in the
event $\alpha(\xi,\eta)$ is too rately equal to one, the number of terms where
the indicator funtion is replaced with its expectation
$\alpha(\xi_{i},\eta_{\ell})$ may be limited, without jeopardising the
variance domination.

##  4 Retrospective: Continuous time Monte Carlo methods

Retrospective simulation (Beskos et al., 2006a, ) is an attempt to take
advantage of the redundancy inherent in modern simulation algorithms
(particularly MCMC, rejection sampling) by subverting the traditional order of
algorithm steps. It is connected to demarginalisation and pseudo-marginal
(Andrieu and Roberts, , 2009) techniques in that it replaces a probability of
acceptance with an unbiased estimation of the said probability, hence creating
an auxiliary variable in the process. In the case of the Metropolis-Hastings
algorithm, this means substituting the ratio

$\dfrac{\pi(\theta^{\prime})q(\theta^{(t)}|\theta^{\prime})}{\pi(\theta^{(t)})q(\theta^{\prime}|^{(t)})}$

with

$\dfrac{\hat{\pi}^{\prime}q(\theta^{(t)}|\theta^{\prime})}{\hat{\pi}^{(t)}q(\theta^{\prime}|^{(t)})}$

where $\hat{\pi}^{\prime}$ is an auxiliary variable such that

$\mathbb{E}[\hat{\pi}^{\prime}|\theta^{\prime}]=\kappa\pi(\theta^{\prime})\qquad\mathbb{E}[\hat{\pi}^{(t)}|\theta^{(t)}]=\kappa\pi(\theta^{(t)})$

Retrospective simulation is most powerful in infinite dimensional contexts,
where its natural competitors are approximate and computationally expensive.
The solution advanced by Beskos et al., 2006a and Beskos et al., 2006b to
simulate diffusions in an exact manner (for a finite number of points) relies
on an auxiliary and bounding Poisson process. The selected points produced
this way actually act as a random sufficient statistic in the sense that the
stochastic process can be generated from Brownian bridges between these points
and closed form estimators conditional of these points may be available and
with a smaller variance. See also Fearnhead et al.,  (2017) for related
results on continuous-time importance sampling (CIS). This includes a
sequential importance sampling procedure with a random variable whose
expectation is equal to the importance weight.141414One difficulty with the
approach is the possible occurrence of negative importance weights (Jacob and
Thiery, , 2015).

##  5 Rao–Blackwellized particle filters

Also known as particle filtering151515An early instance, called bootstrap
filter (Gordon et al., , 1993), involved one of the authors of Gelfand and
Smith,  (1990), who thus contributed to the birth of two major advances in the
field. sequential Monte Carlo (Liu and Chen, , 1998; Doucet et al., , 1999;
Del Moral et al., , 2006) is another branch of the Monte Carlo methodology
where the concept of Rao–Blackwellisation has had an impact. We briefly recall
here that sequential Monte Carlo is used in state-space and other Bayesian
dynamic models where the magnitude of the latent variable prevents the call to
traditional Monte Carlo (and MCMC) techniques. It is also relevant for dealing
with complex static problems by creating a sequence of intermediate and
artificial models, a technique called tempering (Marinari and Parisi, , 1992).

Doucet et al.,  (2000) introduce a general version of the Rao–Blackwellized
particle filter by commenting on the inherent inefficiency of particle filters
in large dimensions, compounded by the dynamic nature of the sampling scheme.
The central filtering equation is a Bayesian update of the form

$p(z_{1:t}|y_{1:t}\propto
p(z_{1:(t-1)}|y_{1:(t-1)})p(y_{t}|z_{t})p(z_{t}|z_{t-1})$ (2)

in a state-space formulation where $(z_{t})_{t}$ (also denoted $z_{1:T}$) is
the latent Markov chain and $(y_{t})_{t}$ the observed sequence. In this
update, the conditional densities of $z_{1:t}$ and $z_{1:(t-1)}$ are usually
unavailable and need be approximated by sampling solutions.

If some marginalisation of the sampling is available for the model at hand,
this reduces the degeneracy phenomenon at the core of particle filters. The
example provided in Doucet et al.,  (2000) is one where $z_{t}(x_{t},r_{t})$,
with

$p(x_{1:t},r_{1:t}|y_{1:t})=p(x_{1:t}|y_{1:t},r_{1:t})p(r_{1:t}|y_{1:t})$

and $p(x_{1:t}|y_{1:t},r_{1:t})$ available in closed form.161616Doucet et al.,
(2000) provide a realistic illustration for a neural network where the
manageable part is obtained via a Kalman filter. This component can then be
used in the approximation of the filtering distribution (2), instead of
weighted Dirac masses, which improves its precision if only by bringing a
considerable decrease in the dimension of the particles (Doucet et al., ,
2000, Proposition 2). It is indeed sufficient to resort only to particles for
the intractable part.

See Andrieu et al.,  (2001); Johansen et al.,  (2012); Lindsten,  (2011);
Lindsten et al.,  (2011) for further extensions on this principle. In
particular, the PhD thesis of Lindsten,  (2011) contains the following and
relevant paragraph:

> Moving from [the particle estimator] to [its Rao–Blackwellized version]
> resembles a Rao–Blackwellisation of the estimator (see also Lehmann, 1983).
> In some sense, we movefrom a Monte Carlo integration to a partially
> analytical integration. However, it is not clear that the Rao–Blackwellized
> particle filter truly is a Rao-Blackwellisation of [the original], in the
> factual meaning of the concept. That is, it is not obvious that the
> conditional expectation of [the original] results in the [its
> Rao–Blackwellized version]. This is due to the nontrivial relationship
> between the normalised weights generated by the [particle filter], and those
> generated by [its Rao–Blackwellized version]. It can thus be said that [it]
> has earned its name from being inspired by the Rao–Blackwell theorem, and
> not because it is a direct application of it.

Nonetheless, any exploitation of conditional properties that does not induce a
(significant) bias is bound to bring stability and faster convergence to
particle filters.

##  6 Conclusion

The term of Rao–Blackwellisation is therefore common enough in the MCMC
literature to be considered as a component of the MCMC toolbox. As we pointed
out in the introduction, many tricks and devices introduced in the past can
fall under the hat of that term and, while a large fraction of them does not
come with a demonstrated improvement over earlier proposals, promoting the
concepts of conditioning and demarginalising as central to the field should be
seen as essential for researchers and students alike. Linking such concepts,
shared by statistics and Monte Carlo, with an elegant and historical result
like the Rao–Blackwell theorem stresses both the universality and the
resilience of the idea.

## References

  * Andrieu et al., (2001) Andrieu, C., de Freitas, N., and Doucet, A. (2001).  Rao–Blackwellised particle filtering via data augmentation.  In Advances in Neural Information Processing Systems (NIPS), pages 561–567.
  * Andrieu and Roberts, (2009) Andrieu, C. and Roberts, G. (2009).  The pseudo-marginal approach for efficient Monte Carlo computations.  Ann. Statist., 37(2):697–725.
  * Atchadé and Perron, (2005) Atchadé, Y. and Perron, F. (2005).  Improving on the independent Metropolis–Hastings algorithm.  Statistica Sinica, 15:3–18.
  * Banterle et al., (2019) Banterle, M., Grazian, C., Lee, A., and Robert, C. P. (2019).  Accelerating Metropolis–Hastings algorithms by delayed acceptance.  Foundations of Data Science, 1(2):103.
  * Berg et al., (2019) Berg, S., Zhu, J., and Clayton, M. K. (2019).  Control variates and Rao–Blackwellization for deterministic sweep Markov chains.  arXiv:1912.06926.
  * Berkson, (1955) Berkson, J. (1955).  Maximum likelihood and minimum $\chi^{2}$ estimates of the logistic function.  J. American Statist. Assoc., 50(269):130–162.
  * Besag, (1989) Besag, J. (1989).  A candidate’s formula: A curious result in Bayesian prediction.  Biometrika, 76(1):183–183.
  * (8) Beskos, A., Papaspiliopoulos, O., Roberts, G., and Fearnhead, P. (2006a).  Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discussion).  J. Royal Statist. Society Series B, 68(3):333–382.
  * (9) Beskos, A., Papaspiliopoulos, O., and Roberts, G. O. (2006b).  Retrospective exact simulation of diffusion sample paths with applications.  Bernoulli, 12(6):1077–1098.
  * Blackwell, (1947) Blackwell, D. (1947).  Conditional expectation and unbiased sequential estimation.  Ann. Statist., 18(1):105–110.
  * Casella and Robert, (1996) Casella, G. and Robert, C. (1996).  Rao-Blackwellization of sampling schemes.  Biometrika, 83:81–94.
  * Casella et al., (2004) Casella, G., Robert, C. P., and Wells, M. T. (2004).  Generalized accept-reject sampling schemes.  Lecture Notes-Monograph Series, 45:342–347.
  * Chatterjee and Diaconis, (2018) Chatterjee, S. and Diaconis, P. (2018).  The sample size required in importance sampling.  Ann. Appl. Probab., 28(2):1099–1135.
  * Chib, (1995) Chib, S. (1995).  Marginal likelihood from the Gibbs output.  J. American Statist. Assoc., 90:1313–1321.
  * Chopin and Robert, (2010) Chopin, N. and Robert, C. (2010).  Properties of nested sampling.  Biometrika, 97:741–755.
  * Cornuet et al., (2012) Cornuet, J.-M., Marin, J.-M., Mira, A., and Robert, C. (2012).  Adaptive multiple importance sampling.  Scandinavian Journal of Statistics, 39(4):798–812.
  * Del Moral et al., (2006) Del Moral, P., Doucet, A., and Jasra, A. (2006).  Sequential Monte Carlo samplers.  J. Royal Statist. Society Series B, 68(3):411–436.
  * Douc and Robert, (2011) Douc, R. and Robert, C. (2011).  A vanilla variance importance sampling via population Monte Carlo.  Ann. Statist., 39(1):261–277.
  * Doucet et al., (1999) Doucet, A., de Freitas, N., and Gordon, N. (1999).  Sequential MCMC in Practice.  Springer-Verlag.
  * Doucet et al., (2000) Doucet, A., Freitas, N. d., Murphy, K. P., and Russell, S. J. (2000).  Rao–Blackwellised particle filtering for dynamic Bayesian networks.  In Proceedings of the 16th Conference on Uncertainty in Artificial Intelligence, UAI ’00, pages 176–183, San Francisco, CA, USA. Morgan Kaufmann Publishers Inc.
  * Elvira et al., (2019) Elvira, V., Martino, L., Luengo, D., and Bugallo, M. (2019).  Generalized multiple importance sampling.  Statist. Science, 34(1):129–155.
  * Fearnhead et al., (2017) Fearnhead, P., Łatuszynski, K., Roberts, G. O., and Sermaidis, G. (2017).  Continious-time importance sampling: Monte Carlo methods which avoid time-discretisation error.  arXiv:1712.06201.
  * Gȧsemyr, (2002) Gȧsemyr, J. (2002).  Markov chain Monte Carlo algorithms with independent proposal distribution and their relation to importance sampling and rejection sampling.  Technical Report 2, Department of Statistics, Univ. of Oslo.
  * Gelfand and Smith, (1990) Gelfand, A. and Smith, A. (1990).  Sampling based approaches to calculating marginal densities.  J. American Statist. Assoc., 85:398–409.
  * Geman and Geman, (1984) Geman, S. and Geman, D. (1984).  Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images.  IEEE Trans. Pattern Anal. Mach. Intell., 6:721–741.
  * Geyer, (1993) Geyer, C. (1993).  Estimating normalizing constants and reweighting mixtures in Markov chain Monte Carlo.  Technical Report 568, School of Statistics, Univ. of Minnesota.
  * Geyer, (1994) Geyer, C. (1994).  Conditioning in Markov chain Monte Carlo.  J. Comput. Graph. Statis., 4:148–154.
  * Gordon et al., (1993) Gordon, N., Salmond, J., and Smith, A. (1993).  A novel approach to non-linear/non-Gaussian Bayesian state estimation.  IEEE Proceedings on Radar and Signal Processing, 140:107–113.
  * Green et al., (2015) Green, P. J., Łatuszyński, K., Pereyra, M., and Robert, C. P. (2015).  Bayesian computation: a summary of the current state, and samples backwards and forwards.  Statistics and Computing, 25(4):835–862.
  * Gutmann and Hyvärinen, (2012) Gutmann, M. U. and Hyvärinen, A. (2012).  Noise-contrastive estimation of unnormalized statistical models, with applications to natural image statistics.  J. Mach. Learn. Res., 13(1):307–361.
  * Hastings, (1970) Hastings, W. (1970).  Monte Carlo sampling methods using Markov chains and their application.  Biometrika, 57:97–109.
  * Jacob et al., (2011) Jacob, P., Robert, C., and Smith, M. (2011).  Using parallel computation to improve independent Metropolis–Hastings based estimation.  J. Comput. Graph. Statist., 20(3):616–635.
  * Jacob and Thiery, (2015) Jacob, P. E. and Thiery, A. H. (2015).  On nonnegative unbiased estimators.  Ann. Statist., 43(2):769–784.
  * Johansen et al., (2012) Johansen, A. M., Whiteley, N., and Doucet, A. (2012).  Exact approximation of Rao–Blackwellised particle filters.  IFAC Proceedings Volumes, 45(16):488 – 493.  16th IFAC Symposium on System Identification.
  * Kong et al., (2003) Kong, A., McCullagh, P., Meng, X. L., Nicolae, D., and Tan, Z. (2003).  A theory of statistical models for Monte-Carlo integration.  J. Royal Statist. Society Series B, 65(3):585–618.
  * Liao, (1998) Liao, J. (1998).  Variance reduction in Gibbs sampler using quasi-random numbers.  J. Comput. Graph. Statist., 7(3):253–266.
  * Lindsten, (2011) Lindsten, F. (2011).  Rao-Blackwellised particle methods for inference and identification.  PhD thesis, University of Linköping, Sweden.
  * Lindsten et al., (2011) Lindsten, F., Schön, T. B., and Olsson, J. (2011).  An explicit variance reduction expression for the Rao–Blackwellised particle filter.  IFAC Proceedings Volumes, 44(1):11979 – 11984.  18th IFAC World Congress.
  * Liu and Chen, (1998) Liu, J. and Chen, R. (1998).  Sequential Monte-Carlo methods for dynamic systems.  J. American Statist. Assoc., 93:1032–1044.
  * Liu et al., (1994) Liu, J., Wong, W., and Kong, A. (1994).  Covariance structure of the Gibbs sampler with applications to the comparisons of estimators and sampling schemes.  Biometrika, 81:27–40.
  * Malefaki and Iliopoulos, (2008) Malefaki, S. and Iliopoulos, G. (2008).  On convergence of importance sampling and other properly weighted samples to the target distribution.  J. Statist. Plann. Inference, 138:1210–1225.
  * Marin and Robert, (2010) Marin, J. and Robert, C. (2010).  On resolving the Savage–Dickey paradox.  Electron. J. Statist., 4:643–654.
  * Marin and Robert, (2011) Marin, J. and Robert, C. (2011).  Importance sampling methods for Bayesian discrimination between embedded models.  In Chen, M.-H., Dey, D., Müller, P., Sun, D., and Ye, K., editors, Frontiers of Statistical Decision Making and Bayesian Analysis, pages 513–527. Springer-Verlag, New York.
  * Marinari and Parisi, (1992) Marinari, E. and Parisi, G. (1992).  Simulated tempering: A new Monte Carlo schemes.  Europhysics letters, 19:451–458.
  * McKeague and Wefelmeyer, (2000) McKeague, I. and Wefelmeyer, W. (2000).  Markov chain Monte Carlo and Rao–Blackwellisation.  J. Statist. Plann. Inference, 85:171–182.
  * Metropolis et al., (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953).  Equations of state calculations by fast computing machines.  J. Chem. Phys., 21:1087–1092.
  * Mira et al., (2001) Mira, A., Møller, J., and Roberts, G. (2001).  Perfect slice samplers.  J. Royal Statist. Society Series B, 63:583–606.
  * Neal, (1999) Neal, R. (1999).  Erroneous results in “Marginal likelihood from the Gibbs output“.  Technical report, University of Toronto.
  * Owen and Zhou, (2000) Owen, A. and Zhou, Y. (2000).  Safe and effective importance sampling.  J. American Statistical Association, 95:135–143.
  * Pearl, (1987) Pearl, J. (1987).  Evidential reasoning using stochastic simulation in causal models.  Artificial Intelligence, 32:247–257.
  * Perron, (1999) Perron, F. (1999).  Beyond accept–reject sampling.  Biometrika, 86(4):803–813.
  * Rao, (1945) Rao, C. (1945).  Information and the accuracy attainable in the estimation of statistical parameters.  Bulletin of the Calcutta Mathematical Society, 37:81–91.
  * Robert and Casella, (2004) Robert, C. and Casella, G. (2004).  Monte Carlo Statistical Methods.  Springer-Verlag, New York, second edition.
  * Robert and Marin, (2008) Robert, C. and Marin, J.-M. (2008).  On some difficulties with a posterior probability approximation technique.  Bayesian Analysis, 3(2):427–442.
  * Robert and Wraith, (2009) Robert, C. and Wraith, D. (2009).  Computational methods for Bayesian model choice.  In Goggans, P. M. and Chan, C.-Y., editors, MaxEnt 2009 proceedings, volume 1193. AIP.
  * Roberts and Rosenthal, (1999) Roberts, G. and Rosenthal, J. (1999).  Convergence of slice sampler Markov chains.  J. Royal Statist. Society Series B, 61:643–660.
  * Sahu and Zhigljavsky, (1998) Sahu, S. and Zhigljavsky, A. (1998).  Adaptation for self regenerative MCMC.  Technical report, Univ. of Wales, Cardiff.
  * Sahu and Zhigljavsky, (2003) Sahu, S. and Zhigljavsky, A. (2003).  Self regenerative Markov chain Monte Carlo with adaptation.  Bernoulli, 9:395–422.
  * Spiegelhalter et al., (1995) Spiegelhalter, D., Thomas, A., Best, N., and Gilks, W. (1995).  BUGS: Bayesian inference using Gibbs sampling.  Technical report, Medical Research Council Biostatistics Unit, Institute of Public Health, Cambridge University.
  * Tanner and Wong, (1987) Tanner, M. and Wong, W. (1987).  The calculation of posterior distributions by data augmentation.  J. American Statist. Assoc., 82:528–550.
  * Vehtari et al., (2019) Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2019).  Pareto smoothed importance sampling.  arXiv:1507.02646.