# Self-regulated ligand-metal charge transfer upon lithium ion de-
intercalation process from LiCoO2 to CoO2

Roberto Fantin  Ambroise van Roekeghem  Anass Benayad  Université Grenoble
Alpes, CEA-LITEN, 17 rue des Martyrs, 38054 Grenoble, France

###### Abstract

Understanding the role of metal and oxygen in the redox process of layered 3d
transition metal oxides is crucial to build high density and stable next
generation Li-ion batteries. We combine hard X-ray photoelectron spectroscopy
and ab-initio-based cluster model simulations to study the electronic
structure of prototypical end-members LiCoO2 and CoO2. The role of cobalt and
oxygen in the redox process is analyzed by optimizing the values of d-d
electron repulsion and ligand-metal p-d charge transfer to the Co 2p spectra.
We clarify the nature of oxidized cobalt ions by highlighting the transition
from positive to negative ligand-to-metal charge transfer upon Li+ de-
intercalation.

††preprint: APS/123-QED

The widespread success of layered lithium transition metal oxides as positive
electrode materials in Li-ion batteries is based on their ability to
reversibly intercalate Li+ ions and exchange electrons while preserving
crystal integrity [1]. The archetype for these materials is LiCoO2, introduced
three decades ago and still one of the most used cathode materials [2, 3].
Although LiCoO2 is now regarded as a conventional material in the battery
community, the fundamental electron transfer mechanism associated with Li+ de-
intercalation and Li-ion cell operation is not yet fully understood,
preventing the increase of usable capacity of this material [4]. Considering
the crystal field splitting of the Co 3d states due to the distorted
octahedral coordination, Li+ de-intercalation from LiCoO2 should be
compensated by cobalt oxidation from t${}_{2g}^{6}$e${}_{g}^{0}$ (Co3+) to
t2g5e${}_{g}^{0}$ (Co4+), both in the low-spin configuration. In practice,
only about half of lithium ions are typically de-intercalated to avoid fast
degradation [4]. Such a limit was explained by the Co 3d band pinning to the
top of the O 2p one, using qualitative band diagram models [5, 6].

The role of oxygen in the redox process was highlighted in the nineties by
density functional theory (DFT) calculations showing that a significant part
of the electron transfer happens at the O sites [7, 8], explaining the O-O
interlayer shrinking observed by in-situ X-ray diffraction (XRD) [9].
Including electronic correlations via dynamical mean field theory (DMFT) does
not change qualitatively this picture: the average total occupation of the Co
3d shell does not change significantly along LiCoO2, Li0.5CoO2, and CoO2 [10].
While oxygen participation to the redox process of LiCoO2 is nowadays accepted
in the literature [11, 12], the nature of the co-participating cobalt is
doubtful, in particular with respect to the commonly-referred Co3+ to Co4+
reaction. Its understanding would help develop new strategies to increase Li-
ion battery energy density by exploiting anionic redox mechanism in transition
metal oxides [4] and [13, 14].

From the structural point of view, the distinction between Co3+ and Co4+ in
LixCoO2 cannot be easily established. Upon delithiation, the crystal structure
of LixCoO2 is overall preserved even through various phase transitions. These
include gliding of the CoO2 layers, distortion of the CoO6 octahedra, and
specific Li+ ordering, but only a gradual contraction of the Co-O bond as
observed by in situ XRD [9, 15, 16]. Even in the case of Li0.5CoO2 it is
unclear if the low-temperature charge ordered phase displays a complete
Co3+/Co4+ separation [17, 18].

X-ray photoelectron spectroscopy (XPS) is a suitable technique to unveil the
redox state of cobalt and oxygen in LixCoO2 since it directly probes the local
electronic structure of ions. Previous XPS studies on the deintercalation
process of LixCoO2 suggested that both cobalt and oxygen participate to the
redox process, as deduced by the analysis of O 1s and Co 2p XPS spectra based
on cluster theory assumptions but without supporting simulation [19, 20]. The
Co 2p core level spectra present satellite structures on the high binding
energy side, a signature of correlation effects that can be exploited to get
deeper insight on the electronic structure of transition metal oxides [21, 22,
23]. The local electronic structure of cobalt and oxygen sites was also
investigated by X-ray absorption spectroscopy (XAS) [24, 25, 26]. Mizokawa et
al. interpreted the O K-edge XAS spectra changes supported by unrestricted
Hartee-Fock density of states calculations for different Co3+/Co4+ mixtures in
the CoO2 triangular lattice. They observed a larger O 2p hole concentration
around the Co4+ ions [26]. While this interpretation is appealing, it does not
explain the satellite structure of XPS spectra. In fact, since the XPS final
state is ionized, this technique is more sensitive to the charge transfer
satellite structures, revealing the complex interplay within the metal-ligand
framework [23, 27].

To our knowledge, a direct comparison between experimental XPS and theoretical
simulations of de-intercalated LixCoO2, including ligand-metal charge transfer
and d-d correlations, is not yet present in the literature. Such aspects are
nonetheless critical to characterize the electronic structure of 3d transition
metal compounds and quantify the d-d electron repulsion (Udd) and the ligand-
metal p-d charge transfer energy ($\Delta$) [28]. In this Letter, we bring new
insight on the charge transfer mechanism of LixCoO2 by combining valence band
XPS and core-level Co 2p hard X-ray photoelectron spectroscopy (HAXPES)
measurements on thin films electrodes with DFT-based single cluster model
calculations for the end-members LiCoO2 and CoO2. We find that delithiation
drives the compound from the mixed valence regime in LiCoO2 to the negative
charge transfer regime in CoO2, leaving the net number of electrons in the Co
3d shell nearly constant. Yet, we observe a reorganization of the electronic
structure: the de-lithiation process is compensated by an electron withdrawal
from the t2g states while there is an electron density backflow from O 2p to
the eg states.

To set the basis for our cluster model Hamiltonian, we first obtained a p-d
tight binding model from the wannierization of converged paramagnetic DFT
calculations of LiCoO2 and CoO2 electronic structures. Therefore, we reduced
at minimum the dependence of our model to semi-empirical parameters, in
particular the on-site energies and hopping parameters, to emphasize the role
of Udd and $\Delta$ on the electronic structure of LiCoO2 and CoO2. The
hopping terms for the CoO6 cluster were extracted from the wannier tight
binding model. The Hamiltonian was then augmented with the Coulomb interaction
and spin-orbit coupling of the 3d shell, while the charge transfer was treated
by means of configuration interaction model [29].

Figure 1:  Ground state character by cluster model calculations of (left panels) LiCoO2 and (right panels) CoO2: (a,e) total Co 3d, (b,f) t2g, and (c,g) eg electronic occupations and (d,h) total spin S as a function of Udd and $\Delta$. The dashed lines indicate the values obtained by exact diagonalization of the p-d tight binding model only. The gray areas indicate the range where we observed best agreement with our photoelectron spectroscopy measurements. Table 1:  Hubbard Udd and charge transfer $\Delta$ energies for LiCoO2 evaluated by (a) semi-empirical cluster model calculations and (b) linear response approach, and (c) constrained random phase approximation, compared to our results. The values in parenthesis are for CoO2 [29]. Udd (eV) |  $\Delta$ (eV) | Reference
---|---|---
3.5 | 4 |  [30] (a)
6.5 | 1 |  [26, 31] (a)
5.5 | -0.5 |  [32] (a)
4.91 (5.37) | – |  [33] (b)
4.40 (3.68) | – |  [34] (c)
4.5 $\pm$ 1.0 | 5.0 $\pm$ 1.5 | This work
(4.0 $\pm$ 0.5) | (-2 $\pm$ 0.5) |

The ground state for the cluster model Hamiltonian was obtained by exact
diagonalization method as implemented in Quanty [35], leaving only Udd and
$\Delta$ as empirical parameters. In all our calculations, we restricted the
Hamiltonian to evaluate configurations between dn and d8 with n = 6 for LiCoO2
and n = 5 for CoO2. Figure 1 shows the expected values of the number of
electrons in the Co 3d, t2g, and eg shells and the total spin S as a function
of Udd and $\Delta$. In each panel, the dashed line shows the value obtained
by the solution of the tight binding Hamiltonian only. For LiCoO2, the
occupation of the Co 3d shell estimated by this method is larger than the
nominal value of six, as expected for the Co3+ oxidation state (Fig. 1a).
Specifically, while the t2g states are almost filled (Fig. 1b), an occupation
of about 1 is observed for the eg orbitals (Fig. 1c). This is a clear effect
of the strong covalence between Co 3d and O 2p orbitals, recalling that the eg
orbitals point towards the ligand O 2p ones in the (distorted) octahedral
crystal field. By introducing electronic correlations in our many-electrons
cluster model, we observe a net decrease of the Co 3d total occupation.
Nevertheless, the cluster calculation is more influenced by $\Delta$ than by
Udd, in particular regarding the eg occupations (Fig. 1c). With low $\Delta$,
the displacement of electrons from O 2p to Co 3d is energetically favored; in
the case of a negative $\Delta$, this process dominates against the intra-
atomic Coulombic repulsion.

To find out which parameters represent better the electronic structure of
LiCoO2, we compared the simulated and HAXPES Co 2p spectra [29]. The gray area
in the panels indicate the range of values that better satisfied these
conditions and we compare our result with other values of Udd and $\Delta$ in
Table 1. In the literature, no agreement is found regarding the electronic
structure of LiCoO2, which was defined as intermediate [30], charge transfer
[31] or even negative charge transfer insulator [32], while from DFT
calculations LiCoO2 is a band insulator, due to the complete filling of the
t2g band [10]. Our simulations highlight that the dominant aspect in the
electronic structure of LiCoO2 is the Co 3d - O 2p covalence, already
described by DFT and subsequently corrected to better describe the effects of
correlations. This allowed us to classify LiCoO2 as a mixed-valence phase
following the modern classification for high-valence transition metal oxides
[36], in line with van Elp et al. [30]. To get a general guideline for Udd,
the screened repulsion values evaluated by first-principles approaches are
also reported in the Table and in line with our findings.

Regarding the electronic structure of CoO2, the values Udd $\approx$ $\Delta$
$\approx$ 4.5 eV obtained for LiCoO2 gave a reasonable starting point from
which we followed the trends identified by Bocquet et al. [22]: with
increasing oxidation state, Udd slightly increases while $\Delta$ abruptly
decreases, because of orbital shrinking and larger electronegativity. The DFT-
based tight binding solution of CoO2 shows a slightly lower Co 3d occupation
than in LiCoO2, although still larger than six electrons (Fig. 1e). While the
t2g shell lost about one electron, the occupation of the eg orbitals even
increased (Fig. 1f,g). The system was found to be in low spin configuration,
in agreement with experimental data [37]. In the cluster calculations, with
decreasing $\Delta$, we observed a transition at $\Delta$ $\approx$ -2 eV from
cobalt high-spin (t${}_{2g}^{3}$e${}_{g}^{2}$, S = 3/2, HS) to low spin
(t${}_{2g}^{5}$e${}_{g}^{1}$, S = 1/2, LS) configuration. The best agreement
with the experimental spectra was obtained for values just below such
transition (Tab. 1). As a direct consequence of the clear negative charge
transfer nature of this compound, the electronic structure in the low spin
region is not in the t${}_{2g}^{5}$e${}_{g}^{0}$ configuration usually
referred to Co4+ oxidation state in the literature. This result gives another
perspective on the role of oxygen in the charge transfer mechanism. The O 2p
participation directly results from the increased hybridization between the eg
and O 2p states while one net electron is extracted from the t2g band. The
oxygen atoms coordinated around cobalt are all equally involved in this
process, with an average hole concentration increasing from 0.09 in LiCoO2 to
0.23 in CoO2 per oxygen atom in the cluster. Moreover, for CoO2, the
electronic occupations obtained by our cluster model are similar to our DFT
predictions and to DFT+DMFT results [10].

Figure 2:  Comparison of LiCoO2 and CoO2 (a,d) XPS valence bands, (b,e) DFT
PDOS, and (c,f) electron removal (PES) and addition (IPES) spectral function.
The PDOS for LiCoO2 in panel b was shifted to match with the main experimental
peak. The O 2p spectral functions were obtained by averaging the curves
obtained from all O atoms in the cluster and scaling them to the actual
composition and to the photoelectron cross section $\sigma$
($\sigma_{Co_{3d}}$ / $\sigma_{O_{2p}}$ $\approx$ 10 [38]). The Udd and
$\Delta$ parameters for LiCoO2 (LS CoO2) are 4.5 and 4.5 eV (4.5 and -2 eV).
For the HS CoO2 phase, $\Delta$ = 4.5 eV. Figure 3:  (a,b) Co 2p core level
photoemission spectra computed within a cluster model and compared to in-lab
HAXPES experimental data with h$\nu$ = 5.4 keV. The depth sensitivity was
estimated with the TPP2M method as $\sim$ 15 nm [39, 40]. The experimental
data is shown as empty circles after background subtraction of an iterated
Shirley function. All simulated spectra were normalized and shifted to the
experimental Co 2p3/2 peak maximum.

Figure 2 shows the experimental XPS valence band of the (a) pristine LiCoO2
and (b) cycled Li0.12CoO2 thin films [29]. The valence band of LiCoO2 presents
a narrow peak at 1-2 eV (A) followed by a large band between 3 to 8 eV (B) and
a smaller one at 11-12 eV (C), in accordance with the literature [19, 20, 41].
For the cycled LixCoO2, we note that the contribution of F 2p, and O 2p from
surface species deposed after Li+ cycling cannot be neglected [42], although
the main valence band contribution is to be assigned to LixCoO2.

The experimental valence spectra are compared to the DFT partial density of
states (PDOS) and the Co 3d and O 2p electron removal spectra for LiCoO2 and
CoO2. For LiCoO2, the DFT PDOS fit well with features A and B, but do not
reproduce the small feature C at higher binding energies (Fig. 2b). In
contrast, the cluster calculation (Fig. 2c) allows to reproduce this peak,
assigned to a charge transfer satellite, while preserving the overall
structure of the PDOS. Features A and C are therefore related to the screened
and unscreened t2g photoelectrons, respectively. Their relative intensity is
proportional to $\alpha$ and $\beta$ in the 1-electron removal wave function,
expressed in the configuration interaction framework of our model as
$\Psi=\alpha\ket{d^{n}}+\beta\ket{d^{n+1}\bm{L}}+\gamma\ket{d^{n+2}\bm{L^{2}}}$,
where $\bm{L}$ denotes a hole in the ligand shells ($\gamma$ being
negligible). Feature B is related to O 2p and eg mixed orbitals, as evident
from both the PDOS and the electron removal spectra. For the CoO2 valence
band, we observe an increasing mixing between the O 2p and Co 3d states and a
closing of the gap between the two bands (Figs. 2e,f), which follows the
experimental observation. Finally, we note that the simulated Co 3d spectra
for the HS phase does not match with the experimental valence band, indicating
that the cobalt ions are indeed in LS state.

To simulate the core photoemission spectra, the ground state Hamiltonian was
augmented with Co 2p core level including core-valence multiplet interactions
and spin-orbit coupling [29]. The resulting spectra obtained with
representative values for Udd and $\Delta$ are compared to experimental hard
X-ray photoelectron spectroscopy (HAXPES) Co 2p spectra in Fig. 3. These were
obtained using a Cr K$\alpha$ X-ray source (5.4 keV) to reduce the
contribution of the uppermost surface contamination (pristine) or degraded
(cycled) surface layers and overcome the overlapping with Co LMM Auger
transitions [40]. As shown in Figs. 3a,b the Co 2p3/2 and Co 2p1/2 spin orbit
components are separated by about 15 eV and are each constituted by a mainline
($\sim$ 780, 795 eV) and a satellite ($\sim$ 790, 805 eV) peak. This
structure, typical for 2p core photoemission in transition metal oxides, has
been explained in the framework of cluster model theory, assigning the
mainline and satellite peaks to the locally screened
($\ket{2p^{5}d^{n+1}\bm{L}}$) and unscreened ($\ket{2p^{5}d^{n}}$) states,
respectively [23, 43]. Our simulated Co 2p spectrum of LiCoO2 (Fig. 3a)
matches with this interpretation, so far only assumed in the literature [19,
20, 41]. To get more insight into the core spectra, we also computed partial
spectra by applying restrictions to the configurations considered by the
transition operator. These highlight that d6 and d7 are the main contributions
(56% and 38%, respectively), in agreement with the ground state electronic
occupations (Figs. 1a-e). The spectral distributions for these configurations
agree with the interpretation given in literature: while a mixed contribution
is observed in the main line, the satellite peak is uniquely present in the d6
partial spectrum. The simulation fits overall well with the experimental data
except for the asymmetry of the main line. This has been referred to non-local
charge transfer screening processes, which indeed are not included in our
single cluster model but can be obtained by either multicluster or DMFT
calculations, as shown for other transition metal oxides [19, 44, 45, 46]. The
larger experimental intensity between the asymmetric main line and the
satellite ($\sim$ 785 eV) can be referred instead to Co2+ ions formed by
surface degradation typically associated to battery cycling and eventually
visible even with HAXPES measurements [40].

In the literature, the formation of Co4+ was related to the broadening of the
main peak [19, 20]. The absence of a new satellite for Co4+ is a consequence
of the low weight of a d5 electronic configuration. The simulated Co 2p
spectra (fig. 3b) of CoO2 assigned to t${}_{2g}^{5}$e${}_{g}^{1}$ electronic
structure (fig. 1), present similarities to LiCoO2, without significant change
of the satellite position. The ground state of CoO2 has a mixed character with
a majority weight for the d6 configuration (10%, 45% and 38% for d5, d6 and d7
configurations, respectively), due to the increase in Co 3d - O 2p
hybridization. Because of such distribution, the d6 satellite peak is still
present while the contribution from d5 configuration is negligible. Note that
the partial spectra in Fig. 3b refers to the total electron number in the 3d
shell, however we observe a reorganization between t2g and eg in our ground
state calculations (Figs. 1f,g). The asymmetry of the main line can be related
to the non-local screening channel due to the metallic character of CoO2. For
the sake of comparison, in Fig. 3b we show the simulated spectra for a HS
configuration (red line), in which the d5 contribution dominates the spectra.

Finally, to understand the correlation between crystal structure change from
O3 to O1 (using Delmas’ notation [47]) and the negative charge transfer
transition, we performed a calculation using the LiCoO2 cluster model with one
extra hole and the same parameters obtained for CoO2. The resulting Co 2p
spectrum (blue line), agrees well with the one obtained by the formal O1 CoO2
structure, suggesting that the electronic structure reorganization is not
mainly driven by the small structural changes observed upon de-lithiation.

In summary, we found a reorganization in the local electronic structure of
CoO2 driven by the decrease of the charge transfer energy towards negative
values. As a negative charge transfer material, the electronic structure of
CoO2 is better described as 3d${}^{6}\bm{L}$, in which the charge extracted
from the t2g states is balanced by a backflow from O 2p to eg orbitals. The
decrease of both O 2p and t2g electron occupations finds good agreement with
the published O K-edge XAS as well as resonant inelastic X-ray spectroscopy
(RIXS) studies [12, 26], which however did not highlight the non-zero
occupation of the eg states.

This combined HAXPES and ab-initio-based cluster model simulations study shows
the importance of considering electronic correlations and charge transfer
theory to understand the redox process of layered 3d transition metal oxides.
As for LixCoO2, the anionic contribution in the redox process was proposed for
Li-rich and Ni-rich oxides at high state of charge, the mechanism behind it
being nowadays under discussion in the literature [48]. Recent studies on the
parent material LixNiO2 proposed a central role of negative charge transfer in
the charge compensation mechanism [49, 50]. This example using the
prototypical LixCoO2 material can therefore lead towards better understanding
of the redox mechanism for these materials, targeted as next generation
positive electrode materials.

We thank J. Ast and S. Motellier for their support on XRD and ICP-Ms
investigations. We acknowledge C. Secouard for providing the LiCoO2 thin films
used in this study. This work was supported by the “Recherches Technologiques
de Base” program of the French National Research Agency (ANR) and by CEA
FOCUS-Battery Program. The work was carried out at the platform of nano-
characterization (PFNC).

## References

  * Liu  _et al._ [2016] C. Liu, Z. G. Neale, and G. Cao, “Understanding electrochemical potentials of cathode materials in rechargeable batteries,” Materials Today 19, 109–123 (2016).
  * Manthiram and Goodenough [2021] A. Manthiram and J. B. Goodenough, “Layered lithium cobalt oxide cathodes,” Nat. Energy 6, 323–323 (2021).
  * Gent  _et al._ [2022] W. E. Gent, G. M. Busse, and K. Z. House, “The predicted persistence of cobalt in lithium-ion batteries,” Nat. Energy 7, 1132–1143 (2022).
  * Lyu  _et al._ [2021] Y. Lyu, X. Wu, K. Wang, Z. Feng, T. Cheng, Y. Liu, M. Wang, R. Chen, L. Xu, J. Zhou, Y. Lu, and B. Guo, “An Overview on the Advances of LiCoO${}_{\textrm{2}}$ Cathodes for Lithium‐Ion Batteries,” Adv. Energy Mater. 11, 2000982 (2021).
  * Goodenough and Park [2013] J. B. Goodenough and K.-S. Park, “The Li-Ion Rechargeable Battery: A Perspective,” J. Am. Chem. Soc. 135, 1167–1176 (2013).
  * Chebiam  _et al._ [2001] R. V. Chebiam, F. Prado, and A. Manthiram, “Soft Chemistry Synthesis and Characterization of Layered $\textrm{Li}_{\textrm{1-x}}\textrm{Ni}_{\textrm{1-y}}\textrm{Co}_{\textrm{y}}\textrm{O}_{\textrm{2-}\delta}$ (0 $\leq$ x $\leq$ 1 and 0 $\leq$ y $\leq$ 1),” Chem. Mater. 13, 2951–2957 (2001).
  * Carlier  _et al._ [2003] D. Carlier, A. Van der Ven, C. Delmas, and G. Ceder, “First-Principles Investigation of Phase Stability in the O2-LiCoO${}_{\textrm{2}}$ System,” Chem. Mater. 15, 2651–2660 (2003).
  * Van der Ven  _et al._ [1998] A. Van der Ven, M. K. Aydinol, G. Ceder, G. Kresse, and J. Hafner, “First-principles investigation of phase stability in Li${}_{\textrm{x}}$CoO${}_{\textrm{2}}$,” Phys. Rev. B 58, 2975–2987 (1998).
  * Tarascon  _et al._ [1999] J. Tarascon, G. Vaughan, Y. Chabre, L. Seguin, M. Anne, P. Strobel, and G. Amatucci, “In Situ Structural and Electrochemical Study of Ni${}_{\textrm{1-x}}$Co${}_{\textrm{x}}$O${}_{\textrm{2}}$ Metastable Oxides Prepared by Soft Chemistry,” Journal of Solid State Chemistry 147, 410–420 (1999).
  * Isaacs and Marianetti [2020] E. B. Isaacs and C. A. Marianetti, “Compositional phase stability of correlated electron materials within DFT + DMFT,” Phys. Rev. B 102, 045146 (2020).
  * Assat and Tarascon [2018] G. Assat and J.-M. Tarascon, “Fundamental understanding and practical challenges of anionic redox activity in Li-ion batteries,” Nat. Energy 3, 373–386 (2018).
  * Hu  _et al._ [2021] E. Hu, Q. Li, X. Wang, F. Meng, J. Liu, J.-N. Zhang, K. Page, W. Xu, L. Gu, R. Xiao, H. Li, X. Huang, L. Chen, W. Yang, X. Yu, and X.-Q. Yang, “Oxygen-redox reactions in LiCoO${}_{\textrm{2}}$ cathode without O-O bonding during charge-discharge,” Joule 5, 720–736 (2021).
  * Chen  _et al._ [2021] J. Chen, W. Deng, X. Gao, S. Yin, L. Yang, H. Liu, G. Zou, H. Hou, and X. Ji, “Demystifying the Lattice Oxygen Redox in Layered Oxide Cathode Materials of Lithium-Ion Batteries,” ACS Nano 15, 6061–6104 (2021).
  * Saubanère  _et al._ [2016] M. Saubanère, E. McCalla, J.-M. Tarascon, and M.-L. Doublet, “The intriguing question of anionic redox in high-energy density cathodes for Li-ion batteries,” Energy Environ. Sci. 9, 984–991 (2016).
  * Chen and Dahn [2004] Z. Chen and J. Dahn, “Methods to obtain excellent capacity retention in LiCoO${}_{\textrm{2}}$ cycled to 4.5 V,” Electroch. Acta 49, 1079–1090 (2004).
  * Mukai  _et al._ [2020] K. Mukai, T. Uyama, and T. Nonaka, “Revisiting LiCoO${}_{\textrm{2}}$ Using a State-of-the-Art In Operando Technique,” Inorg. Chem. 59, 11113–11121 (2020).
  * Takahashi  _et al._ [2007] Y. Takahashi, N. Kijima, K. Tokiwa, T. Watanabe, and J. Akimoto, “Single-crystal synthesis, structure refinement and electrical properties of Li${}_{\textrm{0.5}}$CoO${}_{\textrm{2}}$,” J. Phys.: Condens. Matter 19, 436202 (2007).
  * Motohashi  _et al._ [2009] T. Motohashi, T. Ono, Y. Sugimoto, Y. Masubuchi, S. Kikkawa, R. Kanno, M. Karppinen, and H. Yamauchi, “Electronic phase diagram of the layered cobalt oxide system Li${}_{\textrm{x}}$CoO${}_{\textrm{2}}$ (0.0 $\leq$ x $\leq$ 1.0 ),” Phys. Rev. B 80, 165114 (2009).
  * Dahéron  _et al._ [2008] L. Dahéron, R. Dedryvère, H. Martinez, M. Ménétrier, C. Denage, C. Delmas, and D. Gonbeau, “Electron Transfer Mechanisms upon Lithium Deintercalation from LiCoO${}_{\textrm{2}}$ to CoO${}_{\textrm{2}}$ Investigated by XPS,” Chem. Mater. 20, 583–590 (2008).
  * Ensling  _et al._ [2014] D. Ensling, G. Cherkashinin, S. Schmid, S. Bhuvaneswari, A. Thissen, and W. Jaegermann, “Nonrigid Band Behavior of the Electronic Structure of LiCoO${}_{\textrm{2}}$ Thin Film during Electrochemical Li Deintercalation,” Chem. Mater. 26, 3948–3956 (2014).
  * Zaanen  _et al._ [1986] J. Zaanen, C. Westra, and G. A. Sawatzky, “Determination of the electronic structure of transition-metal compounds: 2p x-ray photoemission spectroscopy of the nickel dihalides,” Phys. Rev. B 33, 8060–8073 (1986).
  * Bocquet  _et al._ [1996] A. E. Bocquet, T. Mizokawa, K. Morikawa, A. Fujimori, S. R. Barman, K. Maiti, D. D. Sarma, Y. Tokura, and M. Onoda, “Electronic structure of early 3d-transition-metal oxides by analysis of the 2p core-level photoemission spectra,” Phys. Rev. B 53, 1161–1170 (1996).
  * van Veenendaal [2006] M. van Veenendaal, “Competition between screening channels in core-level x-ray photoemission as a probe of changes in the ground-state properties of transition-metal compounds,” Phys. Rev. B 74, 085118 (2006).
  * Montoro  _et al._ [2000] L. A. Montoro, M. Abbate, and J. M. Rosolen, “Changes in the Electronic Structure of Chemically Deintercalated LiCoO${}_{\textrm{2}}$,” Electrochem. Solid-State Lett. (2000).
  * Yoon  _et al._ [2002] W.-S. Yoon, K.-B. Kim, M.-G. Kim, M.-K. Lee, H.-J. Shin, J.-M. Lee, J.-S. Lee, and C.-H. Yo, “Oxygen Contribution on Li-Ion Intercalation-Deintercalation in LiCoO${}_{\textrm{2}}$ Investigated by O K-Edge and Co L-Edge X-ray Absorption Spectroscopy,” J. Phys. Chem. B 106, 2526–2532 (2002).
  * Mizokawa  _et al._ [2013] T. Mizokawa, Y. Wakisaka, T. Sudayama, C. Iwai, K. Miyoshi, J. Takeuchi, H. Wadati, D. G. Hawthorn, T. Z. Regier, and G. A. Sawatzky, “Role of Oxygen Holes in Li${}_{\textrm{x}}$CoO${}_{\textrm{2}}$ Revealed by Soft X-Ray Spectroscopy,” Phys. Rev. Lett. 111, 056404 (2013).
  * de Groot  _et al._ [2021] F. M. de Groot, H. Elnaggar, F. Frati, R.-p. Wang, M. U. Delgado-Jaime, M. van Veenendaal, J. Fernandez-Rodriguez, M. W. Haverkort, R. J. Green, G. van der Laan, Y. Kvashnin, A. Hariki, H. Ikeno, H. Ramanantoanina, C. Daul, B. Delley, M. Odelius, M. Lundberg, O. Kuhn, S. I. Bokarev, E. Shirley, J. Vinson, K. Gilmore, M. Stener, G. Fronzoni, P. Decleva, P. Kruger, M. Retegan, Y. Joly, C. Vorwerk, C. Draxl, J. Rehr, and A. Tanaka, “2p x-ray absorption spectroscopy of 3d transition metal systems,” J. of Electron Spectros. Relat. Phenomena 249, 147061 (2021).
  * Zaanen  _et al._ [1985] J. Zaanen, G. A. Sawatzky, and J. W. Allen, “Band gaps and electronic structure of transition-metal compounds,” Phys. Rev. Lett. 55, 418–421 (1985).
  * [29] See Supplemental Material at [URL] for details on delithiation of LiCoO2, XPS and HAXPES characterization, DFT calculations, wannierization, and cluster model calculations.
  * van Elp  _et al._ [1991] J. van Elp, J. L. Wieland, H. Eskes, P. Kuiper, G. A. Sawatzky, F. M. F. de Groot, and T. S. Turner, “Electronic structure of CoO, Li-doped CoO, and LiCoO${}_{\textrm{2}}$,” Phys. Rev. B 44, 6090–6103 (1991).
  * Ikedo  _et al._ [2010] K. Ikedo, Y. Wakisaka, T. Mizokawa, C. Iwai, K. Miyoshi, and J. Takeuchi, “Electronic structure of Li${}_{\textrm{x}}$CoO${}_{\textrm{2}}$ studied by photoemission spectroscopy and unrestricted Hartree-Fock calculations,” Phys. Rev. B 82, 075126 (2010).
  * Mesilov  _et al._ [2013] V. V. Mesilov, V. R. Galakhov, B. A. Gizhevskii, A. S. Semenova, D. G. Kellerman, M. Raekers, and M. Neumann, “Charge states of cobalt ions in nanostructured lithium cobaltite: X-ray absorption and photoelectron spectra,” Phys. Solid State 55, 943–948 (2013).
  * Zhou  _et al._ [2004] F. Zhou, M. Cococcioni, C. A. Marianetti, D. Morgan, and G. Ceder, “First-principles prediction of redox potentials in transition-metal compounds with LDA + U,” Phys. Rev. B 70, 235121 (2004).
  * Kim  _et al._ [2021] B. Kim, K. Kim, and S. Kim, “Quantification of Coulomb interactions in layered lithium and sodium battery cathode materials,” Phys. Rev. Materials 5, 035404 (2021).
  * Haverkort  _et al._ [2012] M. W. Haverkort, M. Zwierzycki, and O. K. Andersen, “Multiplet ligand-field theory using Wannier orbitals,” Phys. Rev. B 85, 165113 (2012).
  * Sawatzky and Green [2016] G. A. Sawatzky and R. J. Green, in  _Quantum Materials: Experiments and Theory_ (Forschungszentrum, Zentralbibliothek, Jülich, 2016).
  * Motohashi  _et al._ [2007] T. Motohashi, Y. Katsumata, T. Ono, R. Kanno, M. Karppinen, and H. Yamauchi, “Synthesis and Properties of CoO${}_{\textrm{2}}$, the x = 0 End Member of the Li${}_{\textrm{x}}$CoO${}_{\textrm{2}}$ and Na${}_{\textrm{x}}$CoO${}_{\textrm{2}}$ Systems,” Chem. Mater. 19, 5063–5066 (2007).
  * Leckey [1976] R. C. G. Leckey, “Subshell photoionization cross sections of the elements for Al K$\alpha$ radiation,” Phys. Rev. A 13, 1043–1051 (1976).
  * Tanuma  _et al._ [1994] S. Tanuma, C. J. Powell, and D. R. Penn, “Calculations of electron inelastic mean free paths. V. Data for 14 organic compounds over the 50-2000 eV range,” Surf. Interface Anal. 21, 165–176 (1994).
  * Fantin  _et al._ [2022] R. Fantin, A. Van Roekeghem, and A. Benayad, “Revisiting Co 2p core‐level photoemission in LiCoO${}_{\textrm{2}}$ by in‐lab soft and hard X‐ray photoelectron spectroscopy: A depth‐dependent study of cobalt electronic structure,” Surface & Interface Analysis , sia.7167 (2022).
  * Galakhov  _et al._ [2002] V. R. Galakhov, V. V. Karelina, D. G. Kellerman, V. S. Gorshkov, N. A. Ovechkina, and M. Neumann, “Electronic Structure, X-ray Spectra, and Magnetic Properties of the LiCoO${}_{\textrm{2}-\delta}$ and Na${}_{\textrm{x}}$CoO${}_{\textrm{2}}$ Nonstoichiometric Oxides,” Phys. Solid State 44 (2002).
  * Dedryvère  _et al._ [2005] R. Dedryvère, S. Laruelle, S. Grugeon, L. Gireaud, J.-M. Tarascon, and D. Gonbeau, “XPS Identification of the Organic and Inorganic Components of the Electrode/Electrolyte Interface Formed on a Metallic Cathode,” J. Electrochem. Soc. 152, A689 (2005).
  * Hufner [2003] S. Hufner,  _Photoelectron Spectroscopy_ , third edition ed. (Springer, 2003).
  * Ghiasi  _et al._ [2019] M. Ghiasi, A. Hariki, M. Winder, J. Kuneš, A. Regoutz, T.-L. Lee, Y. Hu, J.-P. Rueff, and F. M. F. de Groot, “Charge-transfer effect in hard x-ray 1 s and 2 p photoemission spectra: LDA + DMFT and cluster-model analysis,” Phys. Rev. B 100, 075146 (2019).
  * van Veenendaal and Sawatzky [1993] M. A. van Veenendaal and G. A. Sawatzky, “Nonlocal screening effects in 2 p x-ray photoemission spectroscopy core-level line shapes of transition metal compounds,” Phys. Rev. Lett. 70, 2459–2462 (1993).
  * Taguchi  _et al._ [2008] M. Taguchi, M. Matsunami, Y. Ishida, R. Eguchi, A. Chainani, Y. Takata, M. Yabashi, K. Tamasaku, Y. Nishino, T. Ishikawa, Y. Senba, H. Ohashi, and S. Shin, “Revisiting the Valence-Band and Core-Level Photoemission Spectra of NiO,” Phys. Rev. Lett. 100, 206401 (2008).
  * Delmas  _et al._ [1980] C. Delmas, C. Fouassier, and P. Hagenmuller, “Structural classification and properties of the layered oxides,” Physica B+C 99, 81–85 (1980).
  * Zhang  _et al._ [2022] M. Zhang, D. A. Kitchaev, Z. Lebens-Higgins, J. Vinckeviciute, M. Zuba, P. J. Reeves, C. P. Grey, M. S. Whittingham, L. F. J. Piper, A. Van der Ven, and Y. S. Meng, “Pushing the limit of 3d transition metal-based layered oxides that use both cation and anion redox for energy storage,” Nat; Rev. Mater. 7, 522–540 (2022).
  * Foyevtsova  _et al._ [2019] K. Foyevtsova, I. Elfimov, J. Rottler, and G. A. Sawatzky, “LiNiO${}_{\textrm{2}}$ as a high-entropy charge- and bond-disproportionated glass,” Phys. Rev. B 100, 165104 (2019).
  * Menon  _et al._ [2023] A. Menon, B. Johnston, S. Booth, L. Zhang, K. Kress, B. Murdock, G. Paez Fajardo, N. Anthonisamy, N. Tapia-Ruiz, S. Agrestini, M. Garcia-Fernandez, K. Zhou, P. Thakur, T. Lee, A. Nedoma, S. Cussen, and L. Piper, “Oxygen-Redox Activity in Non-Lithium-Excess Tungsten-Doped LiNiO${}_{\textrm{2}}$ Cathode,” PRX Energy 2, 013005 (2023).


# Gaussian-Based Parametric Bijections For Automatic Projection Filters

Muhammad F. Emzir  Zheng Zhao  Lahouari Cheded  Simo Särkkä

###### Abstract

The automatic projection filter is a recently developed numerical method for
projection filtering that makes use of sparse-grid integration and automatic
differentiation. However, its accuracy is highly dependent on the accuracy of
the cumulant-generating function calculation done via sparse-grid integration,
which in turn depends on the selection of the bijection from the canonical
hypercube to the state space. In this paper, we propose two new adaptive
parametric bijections for the automatic projection filter. The first bijection
relies on the minimization of Kullback–Leibler divergence, whereas the second
method employs the sparse-grid Gauss–Hermite quadrature formula. These
bijections let the sparse-grid nodes move adaptively with the high-density
region of the state space, resulting in a substantially improved approximation
while using only a small number of quadrature nodes. The practical
applicability of the methodology is illustrated in three simulated nonlinear
filtering problems.

projection filter, adaptive bijection, numerical quadrature, automatic
differentiation, sparse-grid integration

##  1 Introduction

The projection filter [1, 2] is an approximate optimal filtering method that
is based on projections of probability densities on manifolds, and it
generalizes, for example, assumed density filtering. The projection filter
projects the Kushner–Stratonovich equation [3] of the optimal filtering
solution onto a finite-dimensional manifold of parametric densities, producing
a finite-dimensional stochastic differential equation (SDE). Because the
projection filter is a Hilbert-space projection from a set of probability
densities to a manifold of parametric densities, it converges, under well-
defined conditions, to the solution of Kushner–Stratonovich SPDE as the number
of parameters approaches infinity. When the manifold is the manifold formed by
mixtures of probability densities and the $L^{2}$ metric is used, then the
projection filter is equivalent to a Galerkin method for solving the
Kushner–Stratonovich equation (see [4, Theorem 5.1]).

The projection filter has so far been used in just a few scenarios. Outside of
the Gaussian family, the most prevalent applications have been in univariate
dynamical systems [2, 5, 6, 7]. Recently, Emzir et al. [8] introduced an
automatic projection filter for a wider class of filtering problems. Although
the method is applicable to a large class of multivariate dynamical systems,
its computational demand can be high due to the need for a large number of
quadrature nodes [8]. The large number of nodes is needed because of the fixed
(rather than an adaptive) bijection between the canonical hypercube and the
integration domain. Consequently, the fixed bijection also practically limits
the numerical accuracy of the filter when the computational budget is limited.

In this study, we propose two adaptive bijections for the automatic projection
filter. These bijections are constructed using parametric Gaussian probability
densities. In the first bijection, the Kullback–Leibler divergence between the
current projection density and the Gaussian density is minimized. The second
bijection is constructed such that the quadrature nodes correspond to the
Gauss–Hermite formula for sparse integration. When the automatic projection
filter algorithm is combined with these bijections, the part of the state
space which has high filtering density can be automatically tracked. This
requires much less quadrature nodes than would be needed with a fixed
bijection method. We illustrate the effectiveness of the proposed bijections
using three numerical examples where it is challenging to confine conditional
densities inside a fixed domain. We also extend the class of filtering
problems in comparison to [8] by only requiring that the drift and diffusion
functions belong to a vector space which is closed under partial
differentiation.

The paper is organized as follows. We first review the projection filter for
the exponential family in Section 2, and then explain how to propagate the SDE
parameters using numerical integration and automatic differentiation. Section
3 contains the main contributions of the paper, which is are the two new
bijections from the canonical hypercube to the state space of the projection
filter. Section 4 shows the practical applicability of the bijections in the
automatic projection filter in three simulated nonlinear filtering problems.
Finally, Section 5 summarizes the results and concludes the paper.

##  2 Automatic Projection Filter

In this section, the aim is to review the relevant theoretical results that
constitute the foundation of the automatic projection filter [8] and also to
present the extension of the model considered in this paper. We consider
optimal filtering problems on the following state-space model consisting of
continuous-time stochastic dynamic and observation models:

$\displaystyle dx_{t}$ $\displaystyle=f(x_{t})\,dt+\varrho(x_{t})\,dW_{t},$
(1a) $\displaystyle dy_{t}$ $\displaystyle=h(x_{t})\,dt+dV_{t},$ (1b)

where
$x_{t}\in\mathcal{X}=\mathbb{R}^{d},y_{t}\in\mathcal{Y}=\mathbb{R}^{d_{y}}$.
The processes $\\{W_{t},t\geq 0\\}$ and $\\{V_{t},t\geq 0\\}$ are independent
Wiener processes taking values in $\mathbb{R}^{d_{w}}$ and
$\mathbb{R}^{d_{y}}$ with invertible spectral density matrices $Q_{t}$ and
$R_{t}$ for all $t\geq 0$, respectively. For the sake of exposition, and
without a loss of generality, in the rest of the paper, we assume that
$R_{t}=I$ for all $t\geq 0$.

The conditional probability density of the state at time $t$, $x_{t}$, given a
history of measurements $y_{\tau}$, $0\leq\tau\leq t$, satisfies the
Kushner–Stratonovich equation [3]. Let us define a class of probability
densities $\mathcal{P}$ with respect to the Lebesque measure on $\mathcal{X}$
as $\mathcal{P}=\\{p\in L^{1}:\int_{\mathcal{X}}p(x)\,dx=1,p(x)\geq 0,\forall
x\in\mathcal{X}\\}$. In particular, let us consider the exponential family

$\displaystyle\mathrm{EM}(c):=\left\\{p\in\mathcal{P}:p(x)=\exp(c(x)^{\top}\theta-\psi(\theta))\right\\},$
(2)

where $\theta\in\Theta\subset\mathbb{R}^{m}$ is the natural parameter and
$c:\mathbb{R}^{d}\to\mathbb{R}^{m}$ is a vector of natural statistics that are
assumed to be linearly independent. The natural parameter space $\Theta$ is
defined as

$\displaystyle\Theta:=\left\\{\theta\in\mathbb{R}^{m}:\int_{\mathcal{X}}\exp(c(x)^{\top}\theta)dx<\infty\right\\}.$
(3)

An exponential family is said to be regular if $\Theta$ is an open subset of
$\mathbb{R}^{m}$. In this work, we specifically work with regular exponential
families. In the development of the proposed improved projection filter, we
extensively use the cumulant-generating function (i.e., the log Laplace
transform or log partition function [9, 8]) which is defined as

$\displaystyle\psi(\theta)=\log\left[\int_{\mathcal{X}}\exp(c(x)^{\top}\theta)dx\right],\quad\theta\in\Theta.$
(4)

Because the exponential family is assumed to be regular and the natural
statistics are linearly independent, the exponential family is minimal [10,
11]. We recall the following standard result for a minimal regular exponential
family [10, Theorems 2.2.1 and 2.2.5].

######  Theorem 1.

In a regular exponential family, the set $\Theta$ as defined in (3) is convex.
The cumulant-generating function $\psi(\theta)$ is strictly convex on $\Theta$
and it is differentiable up to an arbitrary order. The moments of the natural
statistics $c_{i}(x)$, $i=1,\ldots,m$ exist for any order, and the
expectations of $c_{i}$ and the corresponding Fisher information matrix $g$
are, respectively, given by,

$\displaystyle\mathbb{E}_{\theta}\left[c_{i}\right]$
$\displaystyle=\partialderivative{\psi(\theta)}{\theta_{i}},$ $\displaystyle
g_{i,j}(\theta)$
$\displaystyle=\partialderivative{\psi(\theta)}{\theta_{i}}{\theta_{j}}.$ (5)

If the representation is minimal, then $g$ is positive definite.

Let us denote by
$\mathcal{S}=\left\\{p_{\theta}:\theta\in\Theta\subseteq\mathbb{R}^{m}\right\\}$
a class of parametric densities which does not have to be the exponential
family. There are a few approaches that can be used to project the
Kushner–Stratonovich equation onto the manifold of parametric densities [12].
The standard way is to leverage the property that the square root of the
density $\sqrt{p_{t}}$ belongs to $L^{2}$. Therefore, by requiring that the
stochastic differential $d\sqrt{p_{t}}$ is an element of $L^{2}$ and by
expressing the Kushner–Stratonovich equation in its Stratonovich form, we can
project $d\sqrt{p_{t}}$ onto the tangent space
$T_{\sqrt{p_{\theta}}}\mathcal{S}^{1/2}$. This is elucidated in the following
lemma [2, Lemma 2.1].

######  Lemma 1.

Let $p_{\theta}\in\mathcal{S}$ and $u$ be a function such that
$\mathbb{E}_{p_{\theta}}[\absolutevalue{u}^{2}]<\infty$. Then
$v:=u\sqrt{p_{\theta}}\in L^{2}$, and its projection onto the tangent space
$T_{\sqrt{p_{\theta}}}\mathcal{S}^{1/2}$ is given by

$\displaystyle\Pi_{\theta}v$
$\displaystyle=\sum_{i=1}^{m}\sum_{j=1}^{m}4g^{-1}_{ij}\expectationvalue{v,\frac{1}{2\sqrt{p_{\theta}}}\partialderivative{p_{\theta}}{\theta_{j}}}\frac{1}{2\sqrt{p_{\theta}}}\partialderivative{p_{\theta}}{\theta_{i}},$
(6)

where the inner product is defined as $\expectationvalue{u,v}=\int
u(x)\,v(x)\,dx$.

For the exponential family $\text{EM}(c)$, the projection filter using
Stratonovich projection is given by [2]:

$\displaystyle d\theta_{t}$
$\displaystyle=g(\theta_{t})^{-1}\mathbb{E}_{\theta_{t}}\left[\mathcal{L}\left[c\right]-\frac{1}{2}h^{\top}h\left[c-\eta(\theta_{t})\right]\right]dt$
$\displaystyle\;+g(\theta_{t})^{-1}\sum_{k=1}^{d_{y}}\mathbb{E}_{\theta_{t}}\left[h_{k}\left[c-\eta(\theta_{t})\right]\right]\circ
dy_{t,k}.$ (7)

In the equation above, $\mathbb{E}_{\theta}$ is the expectation with respect
to the parametric probability density $p_{\theta}$, $h_{k}$ is the $k$-th
element of $h$,
$\theta_{t}\mapsto\eta(\theta_{t}):=\mathbb{E}_{\theta_{t}}[c]$, $\mathcal{L}$
is the backward diffusion operator, and $\circ$ denotes the Stratonovich
multiplication.

We can now extend the class of state-space models from what was considered in
[8]. Let $\mathscr{P}$ be a set of smooth functions
$\varphi:\mathbb{R}^{d}\to\mathbb{R}$ such that $\mathscr{P}$ is a finite-
dimensional vector space that is closed under partial differentiation with
respect to $x_{1},\ldots,x_{d}$. Examples of sets $\mathscr{P}$ are
polynomials on $x_{1},\ldots,x_{d}$ with a total order less or equal to some
$k\in\mathbb{N}$ and functions of the form
$f(x)=\sum_{j=1}^{n_{k}}a_{j}\text{cos}(k_{j}^{\top}x)+b_{j}\text{sin}(\ell_{j}^{\top}x),$
with $a_{j},b_{j}\in\mathbb{R},k_{j},\ell_{j}\in\mathbb{Z}^{d},$ and
$n_{k}\in\mathbb{N}$. We impose the following assumptions on the model (1).

######  Assumption 1.

Elements of functions $f,h,$ and $\varrho\varrho^{\top}$ belong to
$\mathscr{P}$.

######  Assumption 2.

The natural statistics $\left\\{c_{i}\right\\}$ are selected as linearly
independent elements of $\mathscr{P}$.

######  Assumption 3.

Each element of $h$ is in the span of $\left\\{1,c_{1},\ldots,c_{m}\right\\}$.
This means there exists $\lambda_{k}\in\mathbb{R}^{d_{y}}$ for $k=0,\ldots,m$,
such that $h=\lambda_{0}+\sum_{k=1}^{d_{y}}\lambda_{k}c_{k}$.

Under Assumption 3, the resulting exponential family of probability densities
is known as the EM$(c^{\ast})$ family [2]. For this specific family, the
projection filter equation reduces to [2, Theorem 6.3]

$\displaystyle d\theta_{t}$
$\displaystyle=g(\theta_{t})^{-1}\mathbb{E}_{\theta_{t}}\left[\mathcal{L}\left[c\right]-\frac{1}{2}(h^{\top}h)\left[c-\eta(\theta_{t})\right]\right]dt+\sum_{k=1}^{d_{y}}\lambda_{k}dy_{t,k},$
(8)

Assumptions 1 and 2 ensure that every element of
$\\{\mathcal{L}[c],h^{\top}h,h^{\top}hc\\}$ belongs to a larger vector space
$\tilde{\mathscr{P}}$ spanned by the basis of $\mathscr{P}$ and a finite
number of multiplications of basis functions of $\mathscr{P}$. Explicitly, we
can write $\mathcal{L}[c]-\frac{1}{2}(h^{\top}h)c=a_{0}+A_{0}\tilde{c}$ and
$\frac{1}{2}(h^{\top}h)=b_{0}+b_{h}^{\top}\tilde{c}$, where
$a_{0}\in\mathbb{R}^{m}$, $A_{0}\in\mathbb{R}^{m\times(m+m_{h})}$,
$b_{0}\in\mathbb{R}$, $b_{h}\in\mathbb{R}^{(m+m_{h})}$, and
$\tilde{c}^{\top}=[c^{\top},c_{h}^{\top}]$. Note that
$c_{h}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m_{h}}$ is the vector of the
remaining statistics of $x$ that are linearly independent of the elements of
$c$. Therefore, under these assumptions, equation (8) can be expressed as

$\displaystyle
d\theta_{t}=g(\theta_{t})^{-1}\left[a_{0}+b_{0}\eta(\theta_{t})+M(\theta_{t})\tilde{\eta}(\theta_{t})\right]dt+\lambda
dy_{t},$ (9)

where $\lambda=[\lambda_{1},\ldots,\lambda_{m}]^{\top}$,
$M(\theta_{t})=A_{0}+\eta(\theta_{t})b_{h}^{\top}$, and
$\tilde{\eta}(\theta_{t}):=\mathbb{E}_{\theta_{t}}[\tilde{c}]$.

In order to solve (9), we need to compute the expectation $\tilde{\eta}$ and
the Fisher metric $g$. In the automatic projection filter [8], the Fisher
metric $g$ is obtained by automatic differentiation of the approximated
cumulant-generating function via Theorem 1. The cumulant-generating function
can approximated via the Gauss–Chebyshev quadrature method in the univariate
case or via the sparse-grid integration method in the multivariate case. For
the expected values of the extended statistics $\tilde{\eta}$, the automatic
projection filter uses the following lemma to calculate the remaining
expectations in (9) [8, Proposition 5.2].

######  Lemma 2.

Let $s(x):\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a statistic, linearly
independent of the natural statistics $c(x)$. If there exists an open
neighborhood of zero $\Theta_{0}$ such that
$\tilde{\psi}(\tilde{\theta}):=\log(\int_{\mathcal{X}}\exp(\tilde{c}^{\top}\tilde{\theta})dx)<\infty$
for any $\tilde{\theta}\in\Theta\times\Theta_{0}$ where
$\tilde{c}=[c^{\top}(x)\;s(x)]^{\top}$, then for any $\theta\in\Theta$

$\displaystyle\mathbb{E}_{\theta}[s]$
$\displaystyle=\left.\partialderivative{\tilde{\psi}(\tilde{\theta})}{\tilde{\theta}_{m+1}}\right|_{\tilde{\theta}=\theta_{\ast}},$
$\displaystyle\theta_{\ast}=\matrixquantity[\theta\\\ 0].$

##  3 Adaptive Parametric Bijections

In this section, we will present the key contributions of this paper. Before
doing so, we outline briefly the difficulties associated with adopting a
static bijection, as was done in [8].

###  3.1 Challenges of static bijection

Although the automatic projection filter has been shown to perform well with
multivariate dynamics [8], it unfortunately requires a large number of
quadrature nodes. This is due to the need for the accurate computation of the
exponential of the cumulant-generating function (4), which is done via
numerical integration. To compute the cumulant-generating function for a
parametric density $p_{\theta}$ a fixed smooth bijection
$\phi:\mathcal{D}\rightarrow\mathcal{X}$ is used, where
$\mathcal{D}:=(-1,1)^{d}$ is the canonical hypercube. Because during
filtering, the high-density region of the filtering density moves within
$\mathbb{R}^{d}$, fixing the bijection can lead to numerical instabilities. To
illustrate this, let us denote by $\tilde{X}=\\{\tilde{x}_{i}\\}_{i=1}^{N}$
the set of quadrature nodes in the canonical hypercube and consider the
univariate case with
$p_{\theta}(x)=\frac{1}{\sqrt{2\pi}}\exp(-(x-\theta)^{2}/2)$ and
$\phi=\tanh^{-1}(\tilde{x})$. In the domain $\mathcal{D}$, when a substantial
part of the bijected nodes $\phi(\tilde{X})$ lies outside the high-density
interval, the function $\phi^{\prime}(\tilde{x})\,p_{\theta}(\phi(\tilde{x}))$
moves to edge of $\mathcal{D}$, and is almost zero at the other edge. Fig. 1
illustrates this effect with $\theta$ being either $0$ or $\pi/2$, and by
using only 16 quadrature nodes, $\left\\{\phi(\tilde{x}_{i})\right\\}$ does
not cover the region of high $p_{\theta}$ value as $\theta$ moves away from
zero. Even worse, when $\theta$ is far from zero, increasing the number of
quadrature nodes does not lead to a significant decrease in the integration
error. This serious shortcoming is addressed next.

(a) $\theta=0$ (b) $\theta=0$

(c) $\theta=\pi/2$ (d) $\theta=\pi/2$

Figure 1: Illustration of 16 bijected Gauss–Chebyshev’s quadrature nodes,
parametric density $p_{\theta}$, and interpolated approximation using
Chebyshev’s nodes with the bijection $\phi=\tanh^{-1}$. Figs. (a) and (c) are
in the original sample space $\mathcal{X}$ and Figs. (b) and (d) are in the
bijected sample space $\phi(\mathcal{X})$. When $\theta$ is shifted quite far
from zero, the quadrature nodes do not entirely cover most of the high-density
region of $p_{\theta}$.

###  3.2 First Gaussian-Based Parametric Bijection

The aforementioned challenges lead us to propose a new parametric bijection
$\phi_{\xi}$, where the bijection parameter
$\xi\in\mathbb{R}^{n_{\xi}},n_{\xi}\in\mathbb{N}$ can be used to reduce the
numerical integration error. Let us first consider the univariate case and
introduce a parametric probability density $q_{\xi}$ with support
$\mathbb{R}$. Then $\zeta_{\xi}(x):=2\int_{-\infty}^{x}q_{\xi}(y)dy-1$ is a
valid bijection from $\mathbb{R}\rightarrow\mathcal{D}$.

Now let us choose an exponential family manifold with $m$ natural parameters.
Choosing $\phi_{\xi}=\zeta_{\xi}^{-1}$, we can rewrite the exponential of the
cumulant-generating function as

$\displaystyle\exp(\psi(\theta))$
$\displaystyle=\int_{\mathcal{D}}\exp(c(\phi(\tilde{x}))^{\top}\theta)\dfrac{1}{2q_{\xi}(\phi(\tilde{x}))}d\tilde{x}.$
(10)

Instead of directly numerically computing this function using its definition
(4), it is beneficial to use this transformed integral representation instead.
The accuracy of the quadrature approximation to this integral will vary with
the selection of $q_{\xi}$. Recall that a Gaussian quadrature rule with $N$
quadrature nodes is exact when the integrand is chosen to be a polynomial of
order $2N-1$ or less [13]. Therefore, in order to get an accurate integration
result, we can attempt to choose $q_{\xi}$ such that the integrand is close to
a polynomial of order $2N-1$ or less. In fact, if we chose $q_{\xi}$ to be
exactly $p_{\theta}$, then we would have the integrand equal to
$\frac{1}{2}\exp(\psi(\theta))$. This would mean that the integration result
would be exact, thus avoiding any integration errors. However, in the
following proposition we show that there are no possible bijections that can
be evaluated explicitly and which make the integrand in (10) assume a
polynomial form.

######  Proposition 1.

Let $\phi$ be a bijection from $\mathcal{D}$ to $\mathbb{R}$ and let
$\mathrm{EM}(c)$ be the exponential family given by (2) with parameters
$\theta\in\Theta$, where $\Theta$ is given by (3). If there is no explicit
antiderivative for $\beta(x):=\exp(c(x)^{\top}\theta)$ as a function of $x$,
then there is no explicit bijection $\phi$ such that
$\beta(\phi(\tilde{x}))\,\phi(\tilde{x})^{\prime}$ is a polynomial of
$\tilde{x}$.

###### Proof.

Suppose
$\beta(\phi(\tilde{x}))\,\phi(\tilde{x})^{\prime}=\rho_{n}(\tilde{x})$, where
$\rho_{n}$ is a polynomial of order $n$ in $\tilde{x}$. Integrating both sides
gives

$\displaystyle\lim_{s_{\ast}\downarrow-1}\int_{s=s_{\ast}}^{\tilde{x}}\beta(\phi(s))\derivative{\phi(s)}{s}\,ds$
$\displaystyle=\lim_{s_{\ast}\downarrow-1}\int_{s=s_{\ast}}^{\tilde{x}}\rho_{n}(s)\,ds,$
$\displaystyle\lim_{s_{\ast}\downarrow-1}\int_{\phi=\phi(s_{\ast})}^{\phi(\tilde{x})}\beta(\phi)\,d\phi$
$\displaystyle=\int_{s=-1}^{\tilde{x}}\rho_{n}(s)\,ds=\bar{\rho}_{n}(\tilde{x}),$

where $\bar{\rho}_{n}$ is a polynomial because $\rho_{n}$ is. Since the
antiderivative of $\beta$ has no explicit form, there is no explicit form for
$\phi$. ∎

Even if the explicit antiderivative of $\beta(x)$ does exist, that is, we can
form $\int\beta(x)\,dx=\bar{\beta}(x)+C$, then still solving for the bijection
$\phi$ might be impossible as $\bar{\beta}(\phi)$ might not be invertible.
Therefore, instead of attempting to make the integration error exactly zero,
we focus on finding a parametric bijection, or a family of them, which leads
to a smaller integration error. For this purpose, let us use a $d$-dimensional
numerical quadrature with a positive weight function $\omega(\tilde{x})$ and
with $N$ quadrature nodes as

$\displaystyle Q^{d}_{N}g$
$\displaystyle:=\int_{\mathcal{D}}g(\tilde{x})\,\omega(\tilde{x})\,d\tilde{x}\approx\sum_{i=1}^{N}w_{i}\,g(\tilde{x}_{i}),$
(11)

where $\tilde{x}_{i}$ are the nodes and $w_{i}$ are the weights of the
quadrature rule. The multiplier function $\omega(\tilde{x})$ varies for
different quadrature schemes. It is well-known that for a univariate Gauss-
type quadrature, the quadrature nodes are distinct, and the weights
$\\{w_{i}\\}_{i=1}^{N}$ are positive [13, Theorem 1.46].

As a generalization of one-dimensional bijection $\phi_{\xi}$, where
$\partialderivative{\phi_{\xi}}{\tilde{x}}=(2q_{\xi}(x))^{-1}$, in
multivariate context, the bijection operator (also denoted as $\phi_{\xi}$) is
constructed from $q_{\xi}$, where
$\absolutevalue{\det\partialderivative{\phi_{\xi}}{\tilde{x}}}=(2^{d}q_{\xi}(x))^{-1}$.
By using the quadrature rule, we can approximate the expectation of an
arbitrary function $f(x)$ with respect to the density $p_{\theta}$ by

$\displaystyle\mathbb{E}_{\theta,N}\left[f;\xi\right]$
$\displaystyle:=Q^{d}_{N}\left[f(\phi_{\xi})u(\phi_{\xi})\omega^{-1}\right]$
(12)

In this equation, $u(x):=p_{\theta}(x)/(2^{d}q_{\xi}(x))$.

The operator $f\mapsto\mathbb{E}_{\theta,N}[f;\xi]$ is linear, but it does not
guarantee that $\mathbb{E}_{\theta,N}[1;\xi]=1$. The following lemma lists
some important properties of this operator that will be used in the subsequent
developments.

######  Lemma 3.

If the approximation of the cumulant-generating function using $N$ quadrature
nodes is given by
$\psi(\theta)^{(N)}=\log(Q^{d}_{N}\left[\exp(c(\phi_{\xi})^{\top}\theta)\absolutevalue{\det\partialderivative{\phi_{\xi}}{\tilde{x}}}\omega^{-1}\right])$,
then the following hold:

  1. 1.

$\mathbb{E}_{\theta,N}\left[1;\xi\right]=\dfrac{\exp(\psi(\theta)^{(N)})}{\exp(\psi(\theta))}$,

  2. 2.

$\dfrac{\partial\psi(\theta)^{(N)}}{\partial\theta_{j}}=\dfrac{\mathbb{E}_{\theta,N}\left[c_{j};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}$,

  3. 3.

$\dfrac{\partial^{2}\psi(\theta)^{(N)}}{\partial\theta_{j}\partial\theta_{k}}=\dfrac{\mathbb{E}_{\theta,N}\left[c_{j}c_{k};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}-\dfrac{\mathbb{E}_{\theta,N}\left[c_{j};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}\dfrac{\mathbb{E}_{\theta,N}\left[c_{k};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}$,

  4. 4.

$\dfrac{\mathbb{E}_{\theta,N}\left[f;\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}=\sum_{i=1}^{N}w_{i}f(\phi_{\xi}(\tilde{x}_{i}))\\\
\times\exp(c(\phi(\tilde{x}_{i})^{\top}\theta-\psi(\theta)^{(N)}))\absolutevalue{\det\partialderivative{\phi_{\xi}}{\tilde{x}}}\omega(\tilde{x}_{i})^{-1}$.

###### Proof.

For the first equality, using the definition (12), we get:

$\displaystyle\mathbb{E}_{\theta,N}\left[1;\xi\right]$
$\displaystyle=Q^{d}_{N}\left[u(\phi_{\xi})\omega^{-1}\right],$
$\displaystyle=Q^{d}_{N}\left[\exp(c(\phi_{\xi})^{\top}\theta-\psi(\theta))\absolutevalue{\det\partialderivative{\phi_{\xi}}{\tilde{x}}}\omega^{-1}\right]$
$\displaystyle=\frac{\exp(\psi(\theta)^{(N)})}{\exp(\psi(\theta))}.$

The second and third equalities follow directly by taking derivatives of
$\psi(\theta)^{(N)}$ with respect to $\theta_{j}$, and with respect to
$\theta_{j}$ and $\theta_{k}$. The fourth equality follows from the first. ∎

A natural way to choose the parameter $\xi$ is by minimizing the square of
numerical integration error of (10), that is,
$(\exp(\psi(\theta))-\exp(\psi(\theta)^{(N)}))^{2}$. Let us define the
integration error $E_{N}[f(x);\xi]$ as follows:

$\displaystyle E_{N}[f;\xi]$
$\displaystyle:=\mathbb{E}_{\theta}[f]-\mathbb{E}_{\theta,N}[f;\xi].$ (13)

Then because
$\exp(\psi(\theta))-\exp(\psi(\theta)^{(N)})=\exp(\psi(\theta))(1-\mathbb{E}_{\theta,N}\left[1;\xi\right])$,
minimizing the square of the numerical integration error of (10) is equivalent
to minimizing $E_{N}[1;\xi]^{2}$. The following proposition gives a necessary
and sufficient condition such that the squared error $E_{N}[1;\xi]^{2}$ is a
convex function. We denote the partial ordering of two squared matrices as
$A\succeq B$ if $A-B$ is a positive semidefinite matrix.

######  Proposition 2.

Let $u(x)=\frac{p_{\theta}(x)}{2^{d}q_{\xi}(x)}$, where
$q_{\xi}=\exp(\tilde{c}(x)^{\top}\xi-\psi(\xi))$, $c(x)$ and
$\tilde{c}(x)\in\mathcal{C}^{2}(\mathbb{R})$ are the natural statistics of
$p_{\theta}$ and $q_{\xi}$, respectively, $\psi(\xi)$ is the cumulant-
generating function corresponding to $q_{\xi}$, and $\phi_{\xi}$ is a smooth
bijection from $\mathcal{D}$ to $\mathbb{R}^{d}$ constructed from $q_{\xi}$.
Define $\Xi:=\left\\{\xi\in\mathbb{R}^{n_{\xi}}:\psi(\xi)<\infty\right\\}$.
Then $E_{N}[1;\xi]^{2}$ is convex on an open convex set $\Xi_{0}\subseteq\Xi$,
if and only if for any $\xi\in\Xi_{0}$, the following condition is satisfied:

$\displaystyle\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative{u}{\xi};\xi\right]\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative{u}{\xi};\xi\right]^{\top}\succeq
E_{N}[1;\xi]\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative[2]{u}{\xi};\xi\right]$
(14)

###### Proof.

Using the definitions of $E_{N}[1;\xi]$ and
$\mathbb{E}_{\theta,N}\left[\cdot\right]$ from (13) and (12) respectively, we
can write

$\displaystyle\frac{1}{2}\partialderivative{E_{N}[1;\xi]^{2}}{\xi}=-E_{N}[1;\xi]\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative{u}{\xi};\xi\right]$
(15)

and

$\displaystyle\frac{1}{2}\partialderivative[2]{E_{N}[1;\xi]^{2}}{\xi}$
$\displaystyle=\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative{u}{\xi};\xi\right]\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative{u}{\xi};\xi\right]^{\top}$
$\displaystyle-
E_{N}[1;\xi]\mathbb{E}_{\theta,N}\left[\frac{1}{u}\derivative[2]{u}{\xi};\xi\right].$
(16)

Therefore, if for any $\xi\in\Xi_{0}$, condition (14) is satisfied, then
$\partialderivative[2]{E_{N}[1;\xi]^{2}}{\xi}$ is positive semidefinite. By
[14, Theorem 4.5] $E_{N}[1;\xi]^{2}$ is a convex function on $\Xi_{0}$. The
necessary part can be obtained using a similar argument to the proof of [14,
Theorem 4.5]. ∎

We can now highlight two difficulties of using gradient-based method to obtain
an optimal $\xi$ that will minimize $E_{N}[1,\xi]^{2}$. First, by parts 1) and
4) of Lemma 3, the Jacobian and Hessian of $E_{N}[1;\xi]^{2}$ with respect to
$\xi$, which are given respectively by (15) and (16), cannot be explicitly
calculated unless the cumulant-generating function $\psi(\theta)$ is known in
a closed form. Secondly, it is impossible to ensure the criterion that
$E_{N}[1,\xi]^{2}$ is locally convex even on some bounded interval since its
Hessian cannot be evaluated. As a consequence, using gradient methods or
(quasi) Newton methods for finding a locally optimal $\xi$ is not feasible.

With that being said, a more promising approach is to select $\xi$ such that
$q_{\xi}$ covers the high-density area of $p_{\theta}$ as tightly as possible
by optimizing another criterion as explained below, and the resulting
bijection $\phi_{\xi}$ should then be amenable to direct computation. This
turns out to be a viable approach to take, and the thus-constructed bijection
could be shown to work in practice as well as be optimal with respect to the
squared integration error, under some technical conditions. Our starting point
is the following lemma [15].

######  Lemma 4.

Let $q_{\xi}$ be a density from an exponential family with natural statistics
given by $\tilde{c}(x)$. For any distribution $p$, the distribution
$q_{\xi}^{\ast}$ that minimizes $\text{KL}(p||q_{\xi})$, satisfies

$\displaystyle\mathbb{E}_{q_{\xi}^{\ast}}\left[\tilde{c}\right]=\mathbb{E}_{p}\left[\tilde{c}\right].$
(17)

Essentially, the exponential density $q_{\xi}$ that minimizes the
$\text{KL}(p||q_{\xi})$ distance satisfies the moment-matching equality (17).
If we choose $q_{\xi}$ to be a Gaussian density, then the mean $\mu$ and
variance $\Sigma$ of $q_{\xi}$ should satisfy

$\displaystyle\mu_{i}$ $\displaystyle=\mathbb{E}_{\theta}\left[x_{i}\right],$
$\displaystyle i=1,\ldots,d,$ (18a) $\displaystyle\Sigma_{ij}$
$\displaystyle=\mathbb{E}_{\theta}\left[(x_{i}-\mu_{i})(x_{j}-\mu_{j})\right],$
$\displaystyle i,j=1,\ldots,d.$ (18b)

In practical implementation, we replace the cumulant-generating function by
its approximation $\psi(\theta)^{(N)}$ that is obtained from numerical
integration using $N$ quadrature nodes. Since the approximation errors of the
parameters $\mu$, and $\Sigma$, as well as the cumulant-generating function
$\psi(\theta)$, depend linearly on the integration errors on
$\int\exp(c^{\top}x)dx$, then selecting $N$ such that all the approximation
errors of $\mu$, $\Sigma$, and $\psi(\theta)$ are simultaneously kept below
their acceptable upper limits, is feasible (see [8, Theorem 3.1]).

In the following proposition, we show how to construct a smooth bijection
$\zeta_{\xi}$ from $\mathbb{R}^{d}$ to $\mathcal{D}$ using a parametric
density $q_{\xi}$. We also show that its inverse $\phi_{\xi}=\zeta_{\xi}^{-1}$
can be expressed in a closed form when using Gaussian density $q_{\xi}$ with
parameters (18).

######  Proposition 3.

Let $s:\mathbb{R}^{d}\times\mathcal{P}\to\mathcal{D}$ defined by

$\displaystyle s(z,q_{\xi})$
$\displaystyle=\matrixquantity[s_{1}(z_{1},q_{\xi})\\\ \vdots\\\
s_{d}(z_{d},q_{\xi})],$ (19a) $\displaystyle s_{i}(z_{i},q_{\xi})$
$\displaystyle=2\dfrac{\int_{-\infty}^{z_{i}}q_{\xi}(z_{1},\ldots,y_{i},\ldots,z_{d})dy_{i}}{\int_{-\infty}^{\infty}q_{\xi}(z_{1},\ldots,y_{i},\ldots,z_{d})dy_{i}}-1.$
(19b)

If $q_{\xi}$ is a smooth density with support equal to $\mathbb{R}^{d}$, then
$s(\cdot,q_{\xi})$ is a smooth bijection from $\mathbb{R}^{d}$ onto
$\mathcal{D}$.

Assume that $q_{\xi}$ is a Gaussian density with mean $\mu$ and covariance
matrix $\Sigma\succ 0$, and let the eigendecomposition of $\Sigma$ be given by
$\Sigma=T^{-1}\Lambda T$, where $T$ is unitary. Let $\tilde{q}_{\xi}$ be
another Gaussian density with mean at $T\mu$ and variance $\Lambda$. Define
$\zeta_{\xi}:\mathbb{R}^{d}\to\mathcal{D}$ as
$\zeta_{\xi}(x)=s(Tx,\tilde{q}_{\xi})$. Then
$\absolutevalue{\det\partialderivative{\zeta_{\xi}}{x}}=2^{d}q_{\xi}$.
Moreover, the inverse of $\zeta_{\xi}$, denoted as
$\phi_{\xi}:\mathcal{D}\to\mathbb{R}^{d}$, is given by

$\displaystyle\phi_{\xi}(\tilde{x})=\mu+\sqrt{2}T^{-1}\Lambda^{1/2}\erf^{-1}(\tilde{x}),\quad\tilde{x}\in\mathcal{D}$
(20)

where
$\erf^{-1}(\tilde{x})=\matrixquantity[\erf^{-1}(\tilde{x}_{1}),\cdots,\erf^{-1}(\tilde{x}_{d})]^{\top}$.

###### Proof.

Each $s_{i}(z_{i},q_{\xi})$ is a smooth bijection from $\mathbb{R}$ to
$(-1,1)$. By the definition of $s$, for each $\tilde{x}\in\mathcal{D}$ there
exists $z\in\mathbb{R}^{d}$ such that $s(z,q_{\xi})=\tilde{x}$. Let
$z,y\in\mathbb{R}^{d}$. Suppose $s(z,q_{\xi})=s(y,q_{\xi})$. Then since
$s_{i}(z_{i},q_{\xi})$ is a smooth bijection, $z_{i}$ must equal to $y_{i}$,
which leads to $z=y$. The smoothness of $\zeta_{\xi}$ follows from that of
$s_{i}(z_{i},q_{\xi})$ for each $i$.

For the second part, it can be verified that with $z=Tx$ and $\mu_{z}=T\mu$ we
have
$s_{i}(z_{i},\tilde{q}_{\xi})=\erf(\frac{\Lambda_{ii}^{-1/2}(z_{i}-\mu_{z,i})}{\sqrt{2}}).$
Then since $T$ is unitary, we get
$\absolutevalue{\det\partialderivative{\zeta_{\xi}}{x}}=\absolutevalue{\det\partialderivative{s(z,\tilde{q}_{\xi})}{z}}=\prod_{i=1}^{d}2\frac{1}{\sqrt{2\pi\Lambda_{ii}}}\exp(-\frac{1}{2}(z_{i}-\mu_{z,i})^{2})=2^{d}q_{\xi}$.
The inverse of $\zeta_{\xi}$ can be obtained directly from the definition of
$s_{i}$∎

Proposition 3 essentially tells us that for any non-degenerate Gaussian
density $q_{\xi}$, there exists a bijection $\zeta_{\xi}$ such that the
transformation of an infinitesimal volume $dx$ on $\mathbb{R}^{d}$ under
$\zeta_{\xi}$ is equivalent to $2^{d}q_{\xi}$ times $d\tilde{x}$. The
bijection $\phi_{\xi}$ given by (20) can be interpreted as the following
consecutive operations. First it transforms the quadrature nodes from the
canonical hypercube to $\mathbb{R}^{d}$ by the inverse error function, then it
scales each axis by the appropriate square roots of the eigenvalues of
$\Sigma$. Afterwards, it rotates the quadrature nodes according to the unitary
matrix $T$ and lastly it shifts the quadrature nodes by $\mu$ from the origin;
see the middle columns of Figures 4 in Section 4 for an illustration of these
operations. In the numerical implementation, we will use this bijection to
project the Gauss–Chebyshev nodes (for univariate case), or the sparse
Gauss–Patterson nodes (for multivariate case), from $\mathcal{D}$ to
$\mathbb{R}^{d}$ (see [16] for details). We will refer to these numerical
integrations as the Gauss–Chebyshev quadrature (GCQ) and the Gauss–Patterson
quadrature (GPQ), respectively.

Let us turn our attention to checking if the parameters $\xi$ that are
selected via the moment-matching rule (17) also optimize $E_{N}[1;\xi]^{2}$
locally. Proposition 4 below shows that if we use an approximated version of
the moment-matching criterion (17) for a general exponential family
$\text{EM}(\tilde{c})$ (it does not have to be a Gaussian family), then the
selected parameters $\xi$ also optimize $E_{N}[1;\xi]^{2}$ under certain
constraints in the numerical expectation. Explicitly, in a special case where
both natural statistics $c$ and $\tilde{c}$ are equivalent, this constraint
requires that the numerical expectation
$\mathbb{E}_{N,\theta}\left[c;\xi\right]/\mathbb{E}_{N,\theta}\left[1;\xi\right]$
is equal to the true value of $\mathbb{E}_{\theta}\left[c\right]$.

######  Proposition 4.

Let $q_{\xi}$ in (19a) be a density from $\text{EM}(\tilde{c})$ where the
natural statistics $\tilde{c}_{i}\in\mathcal{C}^{2}(\mathbb{R}^{d})$ are
linearly independent. The parameter $\xi$ optimizes $E_{N}[1;\xi]^{2}$ if
$\partialderivative{\psi(\xi)}{\xi}-\frac{\mathbb{E}_{\theta,N}\left[\tilde{c};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}+\frac{\mathbb{E}_{\theta,N}\left[\derivative{x}{\xi}\derivative{m}{x};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}=0.$
In particular, if the following approximated moment-matching rule is used to
choose the parameter $\xi$:

$\displaystyle\dfrac{\mathbb{E}_{\theta,N}\left[\tilde{c};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}=$
$\displaystyle\partialderivative{\psi(\xi)}{\xi},$ (21)

then, the selected parameter $\xi$ is a local optimum of $E_{N}[1;\xi]^{2}$ if
$\mathbb{E}_{\theta,N}\left[\derivative{x}{\xi}\derivative{m}{x};\xi\right]=0.$
If $\tilde{c}=c$ then $\xi$ is a local optimum of $E_{N}[1;\xi]^{2}$ if
$\xi=\theta$.

###### Proof.

Let us write $x(\xi):=\phi_{\xi}(\tilde{x})$. Using (15) from the proof of
Proposition 2, and
$m(\xi):=c^{\top}(x(\xi))\theta-\tilde{c}^{\top}(x(\xi))\xi-(\psi(\theta)-\psi(\xi))$,
we can write

$\displaystyle\frac{1}{2\mathbb{E}_{\theta,N}\left[1;\xi\right]}\derivative{E_{N}[1;\xi]^{2}}{\xi}=$
$\displaystyle-[\partialderivative{\psi(\xi)}{\xi}-\frac{\mathbb{E}_{\theta,N}\left[\tilde{c};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}$
$\displaystyle+\frac{\mathbb{E}_{\theta,N}\left[\derivative{x}{\xi}\derivative{m}{x};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}]E_{N}[1;\xi].$

Hence, if the approximated moment-matching rule (21) is satisfied, and
$\mathbb{E}_{\theta,N}\left[\derivative{m}{x}\derivative{x}{\xi};\xi\right]=0$,
then $\frac{1}{2}\derivative{E_{N}[1;\xi]^{2}}{\xi}=0$, which ensures the
local optimality. The case when $\tilde{c}=c$ follows directly by
substitution. ∎

Note that the approximated moment-matching rule (21) can be implemented as an
iterative procedure, where using initial parameters $\xi_{i}$, one computes
the updated parameters $\xi_{i+1}$ via
$\frac{\mathbb{E}_{\theta,N}\left[\tilde{c};\xi_{i}\right]}{\mathbb{E}_{\theta,N}\left[1;\xi_{i}\right]}=\partialderivative{\psi(\xi_{i+1})}{\xi_{i+1}}$.
The update is repeated until the current iterate is close enough to the fixed
point. To analyze the convergence of this iterative procedure, let us denote
$\hat{\xi}:=\partialderivative{\psi(\xi)}{\xi}$. The mapping from $\xi$ to
$\hat{\xi}$ is a diffeomorphism via the Legendre transformation [10, Theorem
2.2.3]. Therefore, we can write
$F_{N}(\hat{\xi}):=\frac{\mathbb{E}_{\theta,N}\left[\tilde{c};\xi\right]}{\mathbb{E}_{\theta,N}\left[1;\xi\right]}$.
The approximated moment-matching rule can then be written as a Picard
iteration $\hat{\xi}_{i+1}=F_{N}(\hat{\xi}_{i})$. Using Banach’s fixed point
theorem, in the proposition below, we show that there exists a fixed point of
the mapping $F_{N}$ on some subset of
$\hat{\Xi}:=\left\\{\hat{\xi}:\xi\in\Xi\right\\}$.

######  Proposition 5.

Using the notations of Proposition 4, suppose that for each $\theta\in\Theta$,
$\tilde{c}_{i}(\phi_{\xi})u(\phi_{\xi})\omega^{-1}$ belongs to
$\mathcal{W}^{r}_{d}(\mathcal{D})$ for any
$i\in\left\\{1,\ldots,n_{\xi}\right\\}$ and for some $r>0\in\mathbb{N}$,
uniformly for any $\xi\in\Xi_{\theta}$, where,
$\mathcal{W}^{r}_{d}(\mathcal{D}):=\\{f:\mathcal{D}\to\mathbb{R}:\|\frac{\partial^{\absolutevalue{s}}f}{\partial\tilde{x}_{(1)}^{s_{1}}\ldots\partial\tilde{x}_{(d)}^{s_{d}}}\|_{\infty}<\infty,s_{i}<r\\}$,
with $\absolutevalue{s}=\sum_{i=1}^{d}s_{i}$ and $\Xi_{\theta}\subset\Xi$
open. Moreover, assume that $\tilde{c}_{i}(\phi_{\xi})u(\phi_{\xi})$ is
continuously differentiable $w.r.t$ $\xi$ on $\Xi_{\theta}$. Then there exists
an $N_{0}\in\mathbb{N}$ and a subset $\hat{\Xi}_{c}\subset\hat{\Xi}$ such that
$F_{N}$ is a contraction in $\hat{\Xi}_{c}$ for $N\geq N_{0}$.

###### Proof.

$F_{N}$ is continuously differentiable on $\Xi_{\theta}$ by virtue of the
hypothesis that $\tilde{c}_{i}(\phi_{\xi})u(\phi_{\xi})$ is continuously
differentiable $w.r.t$ $\xi$ on $\Xi_{\theta}$ and the definition of $F_{N}$.
Using the mean value theorem [17, §0.27] we can write
$\|F_{N}(\hat{\xi}_{1})-F_{N}(\hat{\xi}_{2})\|\leq\sup_{\hat{\xi}\in\hat{\Xi}_{c}}\|\partialderivative{F_{N}}{\hat{\xi}}\|\|\hat{\xi}_{1}-\hat{\xi}_{2}\|$.
Therefore, if
$\sup_{\hat{\xi}\in\hat{\Xi}_{c}}\|\partialderivative{F_{N}}{\hat{\xi}}\|<1$
in an open convex subset $\hat{\Xi}_{c}$ then $F_{N}$ is a contraction in the
set. Notice that using
$\partialderivative{\xi}{\hat{\xi}}=(\partialderivative[2]{\psi(\xi)}{\xi})^{-1}$
[18, p. 17], we obtain
$\partialderivative{F_{N}}{\hat{\xi}}=(\partialderivative[2]{\psi(\xi)}{\xi})^{-1}\partialderivative{F_{N}}{\xi}$.
Since $\\{\tilde{c}_{i}\\}$ are linearly independent, then
$(\partialderivative[2]{\psi(\xi)}{\xi})$ is invertible for any $\xi\in\Xi$
[11]. Therefore, finding $\hat{\Xi}_{c}$ such that
$\sup_{\hat{\xi}\in\hat{\Xi}_{c}}\|\partialderivative{F_{N}}{\hat{\xi}}\|<1$
is equivalent to finding an open convex subset $\Xi_{c}\subset\Xi_{\theta}$
such that

$\displaystyle\sup_{\xi\in\Xi_{c}}\norm{\left(\partialderivative[2]{\psi(\xi)}{\xi}\right)^{-1}\partialderivative{F_{N}}{\xi}}<1.$
(22)

Since,
$\sup_{\xi\in\Xi_{c}}\|{(\partialderivative[2]{\psi(\xi)}{\xi})^{-1}\partialderivative{F_{N}}{\xi}}\|$
is less than
$\sup_{\xi\in\Xi_{c}}\|(\partialderivative[2]{\psi(\xi)}{\xi})^{-1}\|\sup_{\xi\in\Xi_{c}}\|\partialderivative{F_{N}}{\xi}\|,$
and, for an invertible matrix $T$ we have $\frac{1}{\|T\|}<\|T^{-1}\|$,
condition in (22) is satisfied if
$\sup_{\xi\in\Xi_{c}}\|\partialderivative{F_{N}}{\xi}\|<\inf_{\xi\in\Xi_{c}}\|(\partialderivative[2]{\psi(\xi)}{\xi})^{-1}\|^{-1}<\inf_{\xi\in\Xi_{c}}\|(\partialderivative[2]{\psi(\xi)}{\xi})\|$.
Since $\Xi_{\theta}$ is open, we can select an open subset
$\Xi_{a}\subset\Xi_{\theta}$ away from the boundary of $\Xi$ such that there
exists a positive $\alpha$ satisfying $\alpha
I\prec\partialderivative[2]{\psi(\xi)}{\xi}$ for any $\xi\in\Xi_{a}$. This is
always possible since $\text{EM}(\tilde{c})$ is a regular exponential family,
which means $\Xi$ is an open convex subset of $\mathbb{R}^{n_{\xi}}$ [9].
Since for any
$i,\tilde{c}_{i}(\phi_{\xi})u(\phi_{\xi})\omega^{-1}\in\mathcal{W}^{r}_{d}(\mathcal{D})$
on $\Xi_{\theta}$, then
$\|F_{N}(\hat{\xi}(\xi))-\mathbb{E}_{\theta}\left[\tilde{c}\right]\|=\mathcal{O}(N^{-r}\log(N)^{(d-1)(r-1)})$
on $\Xi_{\theta}$ [19, §4.1.1]. As $N$ approaches infinity, the Jacobian
$\partialderivative{F_{N}}{\xi}$ at any $\xi_{a}\in\Xi_{a}$ decreases to zero.
Therefore, there exists $N_{0}$ and $r_{1}>0$ such that for any $\xi\in
B(\xi_{a},r_{1})\subset\Xi_{a}$, the requirement
$\|\partialderivative*{F_{N}}{\xi}\|<\alpha$ is satisfied for $N\geq N_{0}$ in
$\hat{\Xi}_{c}:=\left\\{\hat{\xi}:\xi\in B(\xi_{a},r_{1})\right\\}$. ∎

The conditions in Proposition 5 can be shown to be satisfied for some $r>0$
when $|p_{\theta}/(q_{\xi}^{r})|$ is bounded and goes to zero as
$\norm{x}\to\infty$ with $\Xi_{\theta}=\Xi$. This condition is valid for the
three numerical simulations considered in Section 4.

###  3.3 Second Gaussian-Based Parametric Bijection

In Section 3.2, we have used the approximated moment-matching rule (21) to
construct the bijection from the hypercube $\mathcal{D}$ to $\mathbb{R}^{d}$.
However, there are other quadrature methods that do not operate on the
hypercube, for example, the Gauss–Hermite quadrature which operates on
$\mathbb{R}^{d}$. To compute the cumulant-generating function efficiently
using these numerical integration methods, we introduce another bijection
where we specifically focus on the Gauss–Hermite quadrature (GHQ) case.

We use $N(d,l)$ to denote the number of quadrature nodes for a dimension $d$
and a level $l$, such that $N:=N(d,l)$. We obtain the sparse multivariate
Gauss–Hermite nodes and weights
$\left\\{\tilde{x}i,w_{i}\right\\}_{i=1}^{N(d,l)}$ by using the Smolyak
construction on one-dimensional Gauss–Hermite nodes (see [20] for details). We
use a similar choice of univariate quadrature nodes per level as in GPQ,
where, for a level $l\in\mathbb{N}$, the corresponding univariate nodes are
given by $N(1,l)=2^{l+1}-1$. This coincides with the choice of nodes made
previously in [21], although in their work, quadrature levels greater than
three were not used. Interestingly enough, if we select
$N(1,l)=\mathcal{O}(2^{l})$, then both GPQ and GHQ would satisfy [16]
$N(d,l)=\mathcal{O}(2^{l}l^{d-1})$.

The exponential of the cumulant-generating function is approximated as
follows:

$\displaystyle\int_{\mathbb{R}^{d}}\exp(c(x)^{\top}\theta)dx$
$\displaystyle=\int_{\mathbb{R}^{d}}\exp(c(\phi_{\xi}(\tilde{x}))^{\top}\theta)\exp(\tilde{x}^{\top}\tilde{x})$
$\displaystyle\times$
$\displaystyle\absolutevalue{\det\partialderivative{\phi_{\xi}(\tilde{x})}{\tilde{x}}}\exp(-\tilde{x}^{\top}\tilde{x})d\tilde{x}\approx\sum_{i=1}^{N(d,l)}w_{i}z(\tilde{x}_{i}),$
(23)

where,

$\displaystyle\phi_{\xi}(\tilde{x})$
$\displaystyle:=\mu+\sqrt{2}T^{-1}\Lambda^{1/2}\tilde{x},$ (24a)
$\displaystyle z(\tilde{x})$
$\displaystyle:=\exp(c(\phi_{\xi}(\tilde{x}))^{\top}\theta)\exp(\tilde{x}^{\top}\tilde{x})\absolutevalue{\det\partialderivative{\phi_{\xi}(\tilde{x})}{\tilde{x}}}.$
(24b)

In (23), $\left\\{\tilde{x}_{i}\right\\}$ is the set of the Gauss–Hermite
quadrature nodes, and $\left\\{w_{i}\right\\}$ are their corresponding
weights. We select $\mu$ and $\Sigma$ according to the moment-matching rule
(17). Therefore, the bijection (24a) also ensures that the sparse
Gauss–Hermite quadrature nodes are always placed in the high-density domain in
$\mathbb{R}^{d}$.

The GHQ scheme suffers from a very weak nesting capability since the
intersection between roots of Hermite polynomials of successive orders
contains only the origin $x=0$. This is in contrast with the GPQ scheme used
in the previous section, which is highly efficient since it has a polynomial
exactness up to order $3N+1$, in addition to being fully nested. Therefore,
more integration nodes would be required by GHQ to achieve the same accuracy
as GPQ (see [13, Section 3.1.2]). A problem with high-order Gauss–Hermite
methods is that the quadrature weights can easily lie below machine precision
[22]. In some applications, fortunately, one can ignore the Gauss–Hermite
nodes that have weights below machine precision and still obtain satisfying
integration results. As we will see in Section 4, with a similar sparse
integration level to that of GPQ, the GHQ scheme combined with the adaptive
bijection (24a) might offer a competitive advantage compared to GPQ combined
with (20) as their quadrature nodes spread wider than those of latter; see
Figure 4.

###  3.4 Integrating the parametric bijection in the automatic projection
filter

We present Algorithm 1 that combines the parametric bijection with the
automatic projection filter algorithm. The bijections (20) and (24a) have
identical parameters $\mu$ and $\Sigma$, and thus can be implemented
similarly.

1:procedure AutomaticProjectionFilter($p_{\theta},\theta,\xi,t,dy_{t}$)

2: $g\leftarrow$ FisherMetric($p_{\theta};\xi$) $\triangleright$ Calculate
Fisher Metric of $p_{\theta}$

3: $\tilde{\eta}\leftarrow$ ExtendedExpectation($p_{\theta};\xi$)
$\triangleright$ Calculate expectation of the extended natural statistics
$\tilde{c}$

4: $\eta\leftarrow\tilde{\eta}[:m]$ $\triangleright$ The expectation of
natural statistics $c$ is the first $m$ of $\tilde{\eta}$

5: $d\theta\leftarrow$ FilterEquation($g,\xi,\tilde{\eta},\eta,t,dy_{t}$)
$\triangleright$ Eq. (9)

6: $\theta\leftarrow\theta+d\theta$ $\triangleright$ Update the parameter

7: $\xi\leftarrow$ UpdateBijectionParameters($\eta,g$) $\triangleright$ Eq.
(21)

8: return $\theta,\xi$$\triangleright$ Return the updated value of the natural
parameter $\theta$ and bijection parameter $\xi$

9:end procedure

Algorithm 1 Single step from the automatic projection filter using parametric
bijection.

##  4 Numerical Examples

In this section, we will consider three numerical examples to assess the
effectiveness of our proposed bijections. For these examples, we choose
$\mathscr{P}$ as the set of polynomials in $x_{t}$ with order less than or
equal to some $n_{p}\in\mathbb{N}$.

###  4.1 Univariate example

The first example is a scalar dynamical system with a nonlinear measurement
model [1, 2]:

$\displaystyle dx_{t}$ $\displaystyle=\kappa dt+\sigma dW_{t},$ (25a)
$\displaystyle dy_{t}$ $\displaystyle=\beta x^{3}_{t}dt+dV_{t},$ (25b)

with two independent standard Wiener processes $\\{W_{t},t\geq 0\\}$ and
$\\{V_{t},t\geq 0\\}$, and three constants given by the process noise constant
$\sigma=0.4$, the drift constant $\kappa=0.25$, and the nonlinear measurement
scale $\beta=0.8$. It is well known that the optimal filter for this problem
is infinite dimensional [23]. We use an exponential manifold with
$c_{i}\in\\{x,x^{2},x^{3},x^{4}\\}$. We remark that the drift term in (25a)
makes the filtering problem significantly harder than that of [8]. We choose
as initial parameters of the projection filter $\theta_{0}=[0,2,0,-1]^{\top}$.
Further, we assume that this initial condition vector reflects exactly the
initial density of the dynamical system which corresponds to a double-mode
non-Gaussian density with peaks at $-1$ and $1$. We solve the
Kushner–Stratonovich stochastic PDE on a uniform grid. The simulation time
step is set to be $10^{-4}$. We generate one measurement realization with
$x_{0}=1$. We compare three different bijections: The first one is a static
bijection $\tanh(\tilde{x})^{-1}$, the second one is the Gaussian-based
bijection (20), and the third one is the Gauss–-Hermite bijection (24a).
Furthermore, we use $9$ and $18$ integration nodes for the static bijection
and $9$ such nodes for both parametrized bijections.

The simulation results, shown in Figure 2, depict the evolution of the
densities obtained with the finite difference approach and with the projection
filters, along with their corresponding Hellinger distances. Even though both
projection filters based on parametric bijections (20) and (24a) require a
smaller number of integration nodes than does the static bijection, they do
however lead to sufficiently reduced Hellinger distances when compared to the
finite difference approximation, as shown in Figure 2. The projection filter
based on the bijection (20) with GCQ produces better approximated densities
than does the GHQ scheme across the entire simulation time. This is also shown
in Figure 2, where the Hellinger distances associated with the GCQ scheme (FD-
vs-Proj-GCQ-9) only oscillate between one and ten times the lowest Hellinger
distance obtained from the static bijection with 18 quadrature nodes (FD-vs-
static-18). In contrast to this and at the end of the simulation, the
Hellinger distance between the finite difference solution and the projection
filter approximation using the static bijection with nine quadrature nodes
(FD-vs-static-9) is about one hundred times greater than that from the lowest
Hellinger distance (FD-vs-static-18). Thus, this simulation clearly shows that
the parametric bijection (20) is superior to both the static and Gauss–Hermite
bijections (24a) for this illustrative example which portrays a hard
stochastic nonlinear filtering problem to solve.

Figure 2: The Hellinger distances from the finite difference solutions to
densities solved by projection filters with static bijections with $9$ and
$18$ nodes (FD-vs-static-9 and FD-vs-static-18), or with parametric bijections
(20) and (24a) with $9$ nodes, (FD-vs-Proj-GCQ-9 and FD-vs-Proj-GHQ-9),
respectively.

###  4.2 Modified Van-der-Pol Oscillator

In this section, we compare the projection filter with a bootstrap particle
filter with systematic resampling [24]. The dynamic model considered here is a
modified Van-der-Pol oscillator, the standard form of which is widely used as
a model for oscillatory processes in physics, electronics, biology, neurology,
sociology and economics:

$\displaystyle d\matrixquantity[x_{1,t}\\\ x_{2,t}]$
$\displaystyle=\matrixquantity[\kappa x_{1,t}+x_{2,t}\\\ -x_{1,t}+\kappa
x_{2,t}+\mu(1-x_{1,t}^{2})x_{2,t}]dt+\matrixquantity[0\\\ \sigma_{w}]dW_{t},$
(26a) $\displaystyle dy$ $\displaystyle=x_{1,t}dt+\sigma_{v}dV_{t}.$ (26b)

In this simulation, we set $\mu=0.3$, $\kappa=1.25$ and
$\sigma_{v}=\sigma_{w}=1$. We also set the simulation time step to be
$2.5\times 10^{-3}$. Unlike the case with $\kappa=0$, where the probability
density evolution can be easily contained using a compact support [8], the
probability densities corresponding to (26) with $\kappa>0$ expand quickly in
time. We use here both GPQ and GHQ with their sparse-grid integration schemes
where we set the level to four. For the GPQ scheme, the number of nodes used
is 129 while for the GHQ scheme, it is 189 after ignoring all nodes with
weights less than $10^{-9}$. For the particle filter, we use $9.6\times
10^{6}$ samples in our simulation. As a reference and for comparison purposes,
4097 quadrature nodes were needed to solve a related Van-der-Pol oscillator
filtering problem in [8] with $\kappa=0$. We discretize the dynamic model (26)
using the Euler-Maruyama scheme for both the particle filter and the
measurement process. For a multi index $\underline{i\mkern-2.0mu}\mkern
2.0mu\in\mathbb{N}$, define $x^{\underline{i\mkern-2.0mu}\mkern
2.0mu}=x_{1}^{i(1)}x_{2}^{i(2)}$. The natural statistics are set to be
${x^{\underline{i\mkern-2.0mu}\mkern 2.0mu}}$, where
$1\leq\absolutevalue{\underline{i\mkern-2.0mu}\mkern 2.0mu}\leq 4$. Further,
the initial density is set to be the standard Gaussian density.

To show the performance of the projection filter obtained by both sparse
integration schemes, we calculate the empirical densities of the particle
filter on a fixed grid. The comparison of the empirical densities from the
particle filters with those from the projection filters can be seen in Figure
4, while the Hellinger distances between the empirical densities and those
from the projection filters can be seen in Figure 3. Figure 4 shows that the
bijected quadrature nodes of both sparse GPQ and GHQ schemes are
systematically adapting to the shapes of the densities as they evolve in time.
The shapes of the projection densities resemble those of the empirical
densities from the particle filter, except that the projection filter’s
densities cannot capture some sharp notch-like shapes with low values as in
Figure 4. The Hellinger distances from the projection filter’s densities
obtained using the GHQ scheme to the empirical densities are slightly lower
compared to those obtained using the GPQ scheme.

Figure 3: Hellinger distance from the empirical densities to the projection
filter’s densities solved using both GPQ and GHQ.

Figure 4: Comparison of empirical densities from particle filter (left) and
densities from the projection filters solved using GPG (center) and GHQ
(right) at $t=0.995$. The grey dots represents the position of the bijected
quadrature nodes of the sparse Gauss–Patterson or Gauss–Hermite schemes,
respectively.

###  4.3 Stochastic Epidemiology Application

In this section, we consider an application of the projection filter to the
stochastic suspected-infected-recovery (SIR) nonlinear filtering problem,
widely used for example in the study of the spread of infectious diseases.
Consider the following SDE [25]:

$\displaystyle d\matrixquantity[x_{1,t}\\\ x_{2,t}]$
$\displaystyle=\matrixquantity[-\beta x_{1,t}x_{2,t}-\mu x_{1,t}+\mu\\\ \beta
x_{1,t}x_{2,t}-\left(\lambda+\mu\right)x_{2,t}]dt+\matrixquantity[-\sigma
x_{1,t}x_{2,t}\\\ \sigma x_{1,t}x_{2,t}]dW_{t},$ (27a) $\displaystyle
dy_{1,t}$ $\displaystyle=x_{2,t}dt+kdV_{t}.$ (27b)

In this equation, $x_{1,t}$ is the fraction of the suspected population,
$x_{2,t}$ is the fraction of the infected population, assuming that the
population size is constant. The fraction of the recovered population is given
by $1-(x_{1,t}+x_{2,t})$, and its dynamic is non-stochastic which can be
excluded from the SDE. The constants $\beta,\mu,\lambda$ correspond to the
average number of contacts per infection per day, the birth rate, and the
recovery rate of the infected people, respectively. The constants $\sigma$,
and $k$ are positive constants associated with the process and measurement
noises, respectively.

If $0<\beta<\min(\lambda+\mu-\frac{\sigma^{2}}{2},2\mu)$ then the disease-free
equilibrium $x_{\ast}=(1,0)$ is globally asymptotically stable [25, Theorem
2.1]. Therefore, the stationary probability measure has no density with
respect to the Lebesgue measure. The high-density domain of the conditional
probability density of SIR dynamics will be shrinking in time. In our
simulation, we chose $\mu=0.2,\beta=0.126,\lambda=0.1,\sigma=0.2,k=10^{-4}$
and set the initial density to be a Gaussian density with mean
$[0.95,0.02]^{\top}$ and variance $\text{diag}[0.95,0.02]\times 10^{-3}$.
Using a similar exponential family to that of the previous section, we compare
here the performance of the projection filter achieved by using the bijection
(20) and GPQ level 5 with that of the particle filter. The Hellinger distance
between the two densities at different times can be seen in Figure 5. This,
therefore, clearly shows that the bijection (20) can also perform very well in
accurately tracking probability densities that shrink with time, as is the
case here with SIR dynamics, and, by the same token, also demonstrates the
adaptive capability of the proposed bijection (20).

Figure 5: Hellinger distance from the empirical densities to the projection
filter’s densities solved using GPQ level 5.

##  5 Conclusions

In this work, we have introduced two new parametric bijections to be used in
conjunction with the automatic projection filter that was recently proposed in
[8]. The first bijection was constructed by selecting a Gaussian density
$q_{\xi}$, whose parameters were obtained by minimizing the KL divergence
between the projected density $p_{\theta}$ and a Gaussian density $q_{\xi}$,
which is equivalent to evaluating moment-matching conditions. We have shown
that this bijection could also optimize the squared integration error under
some sufficient theoretical conditions. The second bijection was also
constructed via the same moment-matching conditions, but it was tailored for
the GHQ scheme. We then applied these bijections to three practically-
motivated numerical examples and found that they all achieved far superior
performance in terms of the Hellinger distances to the ground truth than did
the static bijection, using far fewer quadrature nodes and thus achieving both
higher accuracy and reduced computational time.

## References

  * [1] B. Hanzon and R. Hut, “New results on the projection filter,” in _European Control Conference_ , Grenoble, Jul. 1991, p. 9.
  * [2] D. Brigo, B. Hanzon, and F. L. Gland, “Approximate nonlinear filtering by projection on exponential manifolds of densities,” _Bernoulli. Official Journal of the Bernoulli Society for Mathematical Statistics and Probability_ , vol. 5, no. 3, p. 495, Jun. 1999.
  * [3] H. J. Kushner, “On the differential equations satisfied by conditional probability densities of Markov processes, with applications,” _Journal of the Society for Industrial and Applied Mathematics, Series A: Control_ , vol. 2, no. 1, pp. 106–119, 1964.
  * [4] J. Armstrong and D. Brigo, “Nonlinear filtering via stochastic PDE projection on mixture manifolds in L2 direct metric,” _Mathematics of Control, Signals, and Systems_ , vol. 28, no. 1, p. 5, Dec. 2016.
  * [5] ——, “Optimal approximation of SDEs on submanifolds: The Itô-vector and Itô-jet projections,” Oct. 2016.
  * [6] S. Koyama, “Projection smoothing for continuous and continuous-discrete stochastic dynamic systems,” _Signal Processing_ , vol. 144, pp. 333–340, Mar. 2018.
  * [7] A. Kutschireiter, L. Rast, and J. Drugowitsch, “Projection Filtering with Observed State Increments with Applications in Continuous-Time Circular Filtering,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 686–700, 2022.
  * [8] M. F. Emzir, Z. Zhao, and S. Särkkä, “Multidimensional projection filters via automatic differentiation and sparse-grid integration,” _Signal Processing_ , vol. 204, p. 108832, Mar. 2023.
  * [9] L. D. Brown, “Fundamentals of Statistical Exponential Families with Applications in Statistical Decision Theory,” _Lecture Notes-Monograph Series_ , vol. 9, pp. i–279, 1986.
  * [10] R. E. Kass and P. W. Vos, _Geometrical Foundations of Asymptotic Inference_ , ser. Wiley Series in Probability and Statistics. New York: Wiley, 1997, ”A Wiley Interscience publication.”.
  * [11] O. Calin and C. Udrişte, _Geometric Modeling in Probability and Statistics_. Springer International Publishing, 2014.
  * [12] D. Brigo, “Optimal Projection Filters,” May 2022, comment: arXiv admin note: text overlap with arXiv:1610.03887.
  * [13] W. Gautschi, _Orthogonal Polynomials: Computation and Approximation_ , ser. Numerical Mathematics and Scientific Computation. Oxford New York: Oxford university press, 2004.
  * [14] R. T. Rockafellar, _Convex Analysis_ , ser. Princeton Mathematical Series. Princeton, N.J: Princeton University Press, 1970, no. 28.
  * [15] R. Herbrich, “Minimising the Kullback–Leibler Divergence,” Microsoft, Tech. Rep., 2005.
  * [16] T. Gerstner and M. Griebel, “Numerical integration using sparse grids,” _Numerical Algorithms_ , vol. 18, no. 3/4, pp. 209–232, 1998.
  * [17] T. Aubin, _A Course in Differential Geometry_ , ser. Graduate Studies in Mathematics. Providence, R.I: American Mathematical Society, 2001, no. v. 27.
  * [18] S.-i. Amari, _Information Geometry and Its Applications_ , ser. Applied Mathematical Sciences. Tokyo: Springer Japan, 2016, vol. 194.
  * [19] M. Holtz, _Sparse Grid Quadrature in High Dimensions with Applications in Finance and Insurance_ , ser. Lecture Notes in Computational Science and Engineering. Heidelberg New York: Springer, 2011, no. 77.
  * [20] H.-J. Bungartz and M. Griebel, “Sparse grids,” _Acta Numerica_ , vol. 13, pp. 147–269, May 2004.
  * [21] B. Jia, M. Xin, and Y. Cheng, “Sparse Gauss-Hermite Quadrature Filter with Application to Spacecraft Attitude Estimation,” _Journal of Guidance, Control, and Dynamics_ , vol. 34, no. 2, pp. 367–379, Mar. 2011.
  * [22] L. N. Trefethen, “Exactness of Quadrature Formulas,” _SIAM Review_ , vol. 64, no. 1, pp. 132–150, Feb. 2022.
  * [23] M. Hazewinkel, S. I. Marcus, and H. J. Sussmann, “Nonexistence of finite-dimensional filters for conditional statistics of the cubic sensor problem,” _Systems & Control Letters_, vol. 3, no. 6, pp. 331–340, Dec. 1983.
  * [24] N. Chopin, _An Introduction to Sequential Monte–Carlo_. Cham, Switzerland: Springer, 2020.
  * [25] E. Tornatore, S. Maria Buccellato, and P. Vetro, “Stability of a stochastic SIR system,” _Physica A: Statistical Mechanics and its Applications_ , vol. 354, pp. 111–126, Aug. 2005.