MITP/23-050

Endpoint divergences in inclusive $\bar{B}\to X_{s}\gamma$111Based on talks
given by R.S. at the 20th annual workshop on Soft-Collinear Effective Theory
(SCET 2023), Berkeley, March 27-30, 2023 and by T.H. at 21st Conference on
Flavor Physics and CP Violation (FPCP 2023), Lyon, May 29 - June 2, 2023 and
on Ref. [17].

Tobias Hurtha and Robert Szafron,b

aPRISMA+ Cluster of Excellence and Institute of Physics (THEP),
Johannes Gutenberg University, D-55099 Mainz, Germany
bDepartment of Physics, Brookhaven National Laboratory, Upton, N.Y., 11973,
U.S.A.

###### Abstract

The subleading so-called resolved contributions represent the largest
uncertainty in the inclusive decay mode $\bar{B}\to X_{s}\gamma$. However
there had been no complete proof of factorization of these subleading
contributions. This failure of factorisation can be traced back to endpoint
divergences and cured by recently proposed refactorisation techniques.

##  1 Introduction

The inclusive decay modes $\bar{B}\to X_{s,d}\gamma$ and $\bar{B}\to
X_{s,d}\ell^{+}\ell^{-}$ are golden modes of the indirect search for new
physics at the Belle-II experiment [1, 2, 3]. The Belle II experiment will
collect two orders of magnitude larger data samples than the $B$ factories [4]
and will measure the inclusive decays with the highest precision. Therefore,
subleading effects have become now relevant on the theory side.

Within the heavy mass expansion (HME), these inclusive so-called penguin modes
are dominated by the partonic contributions, which can be calculated
perturbatively, and sub-leading contributions start at the quadratic level,
$(\Lambda/m_{b})^{2}$ only. However, it is well known that this operator
product expansion breaks down in these inclusive modes if one considers
operators beyond the leading ones. This breakdown manifests in nonlocal power
corrections, also called resolved contributions. They are characterised by
containing subprocesses in which the photon couples to light partons instead
of connecting directly to the effective weak-interaction vertex [5].

These resolved contributions can be systematically calculated using soft-
collinear effective theory (SCET). In the case of the inclusive $\bar{B}\to
X_{s}\gamma$ decay, all resolved contributions to $O(1/m_{b})$ have been
calculated some time ago [6]. Also, the analogous contributions to the
inclusive $\bar{B}\to X_{s,d}\ell^{+}\ell^{-}$ decays have been analysed to
$O(1/m_{b})$ [7, 8].

Recently, the uncertainty due to the resolved contribution was reduced with
the help of a new hadronic input [9, 10]. But these resolved contributions
still represent the largest uncertainty of order $4-5\%$ in the prediction of
the inclusive decay rate $\bar{B}\to X_{s}\gamma$ [12] and of the low-$q^{2}$
observables of $\bar{B}\to X_{s,d}\ell^{+}\ell^{-}$ [13, 14, 15].

Moreover, a large scale dependence and also a large charm mass dependence were
identified in the lowest order result of the resolved contribution, which
calls for a systematic calculation of $\alpha_{s}$ corrections and
renormalisation group (RG) summation [10, 16]. A mandatory prerequisite for
this task is an all-order in the strong coupling constant $\alpha_{s}$
factorisation formula for the subleading power corrections.

One finds a factorisation formula for the various contributions to the
inclusive penguin decays [6], where the symbol $\otimes$ denotes the
convolution of the soft and jet functions.

$\displaystyle d\Gamma(\bar{B}\to
X_{s}\,\gamma,\ell^{+}\ell^{-})=\sum_{n=0}^{\infty}\,\frac{1}{m_{b}^{n}}\,\sum_{i}\,H_{i}^{(n)}J_{i}^{(n)}\otimes
S_{i}^{(n)}$
$\displaystyle\qquad+\sum_{n=1}^{\infty}\,\frac{1}{m_{b}^{n}}\,\bigg{[}\sum_{i}\,H_{i}^{(n)}J_{i}^{(n)}\otimes
S_{i}^{(n)}\otimes\bar{J}_{i}^{(n)}+\sum_{i}\,H_{i}^{(n)}J_{i}^{(n)}\otimes
S_{i}^{(n)}\otimes\bar{J}_{i}^{(n)}\otimes\bar{J}_{i}^{(n)}\bigg{]}\,.$ (1)

The first line describes the so-called direct contributions, while the second
line contains the resolved contributions. The latter appear first only at the
order $1/m_{b}$ in the heavy-quark expansion. Here hard functions
$H_{i}^{(n)}$ describe physics at the high scale $m_{b}$. $J_{i}^{(n)}$ are
the so-called jet functions which represent the physics of the hadronic final
state $X_{s}$ at the intermediate hard-collinear scale
$\sqrt{m_{b}\Lambda_{\rm QCD}}$. The soft functions $S_{i}^{(n)}$, the so-
called shape functions, parameterise the hadronic physics at the scale
$\Lambda_{\text{QCD}}$. Within the resolved contributions, we have a new
ingredient in the factorisation formula, the so-called anti-hardcollinear jet
functions $\bar{J}_{i}^{(n)}$ due to the coupling of virtual photons with
virtualities of order $\sqrt{m_{b}\Lambda_{\text{QCD}}}$ to light partons
instead of the weak vertex directly. They are not represented by cut
propagators as the usual jet functions but as full propagator functions
dressed by Wilson lines.

However, the specific resolved ${O}_{8g}-{O}_{8g}$ contribution does not
factorise because the convolution integral is ill-defined. The authors of Ref.
[6] claimed that there is an essential difference between divergent
convolution integrals in power-suppressed contributions of exclusive $B$
decays and the divergent convolution integral in the present case, while the
former were of IR origin, the latter divergence were of UV nature.
Nevertheless, using a hard cut-off in the resolved contribution, the sum of
direct and resolved ${O}_{8g}-{O}_{8g}$ contributions was shown to be scale
and scheme independent at the lowest order. However the failure of
factorisation did not allow for a consistent resummation of large logarithms.
In a recent paper, the divergences in the resolved and in the direct
contributions were identified as endpoint divergences. It was shown that also
the divergence in the direct contribution can be traced back to a divergent
convolution integral [17].

Recently new techniques [18, 19, 20, 21] were presented in specific collider
applications, which allow for an operator-level reshuffling of terms within
the factorisation formula so that all endpoint divergences cancel out. This
idea of refactorisation was now implemented in this flavour example of the
resolved contributions, which includes nonperturbative soft functions, and the
subleading shape functions, not present in collider applications [17]. This is
the first QCD application of these new refactorsiation techniques in flavour
physics.

In the following, we present the various steps of this analysis. This analysis
leads to a renormalised factorisation theorem on the operator level for these
resolved contributions to all orders in the strong coupling constant. This
result establishes the validity of the general factorisation theorem, given in
Eq. 1, - also for the ${O}_{8g}-{O}_{8g}$ contributions. This theorem now
allows for higher-order calculations of the resolved contributions and
consistent summation of large logarithms [17].

##  2 Matching on SCET and degeneracies in the EFT

The first step in the derivation of a factorisation theorem is hard matching.
We have to match the electroweak operator onto SCET. The operator under
consideration, ${O}_{8}$ 222For the various conventions we refer the reader to
Ref. [17]:

${O}_{8g}=-\frac{g_{s}}{8\pi^{2}}\,m_{b}\,\bar{s}\sigma_{\mu\nu}(1+\gamma_{5})G^{\mu\nu}b\,,$
(2)

matches onto two possible SCET operators. The direct contribution is
represented by a next-to-leading power (NLP) $B$-type current in SCET, i.e. a
power-suppressed current composed of two collinear building blocks (see left
figure of Figure 2):

$\mathcal{O}_{8g}^{B1}\left(u\right)=\int\frac{dt}{2\pi}e^{-ium_{b}t}\overline{\chi}_{hc}\left(t\bar{n}\right)\gamma_{\nu\perp}\,Q_{s}\,\mathcal{B^{\nu}}_{\overline{hc}\perp}\left(0\right)\gamma_{\mu\perp}\,\mathcal{A^{\mu}}_{hc\perp}\left(0\right)\left(1+\gamma_{5}\right)h\left(0\right),$
(3)

with $Q_{s}$ as the electric charge of the strange quark in units of $e$ and
the electromagnetic gauge-invariant transverse photon field

$\mathcal{B^{\nu}}_{\overline{hc}\perp}=e\left(A_{\perp}^{\nu}-\frac{\partial_{\perp}^{\nu}}{n\partial}nA\right).\,$
(4)

The second operator in SCET on which the electroweak operator ${O}_{8g}$
matches is a leading-power (LP) $A$-type current, i.e. a leading-power current
including one collinear building block, :

$\mathcal{O}_{8g}^{A0}\left(0\right)=\overline{\chi}_{\overline{hc}}\left(0\right)\frac{\not{n}}{2}\gamma_{\mu\perp}\mathcal{A^{\mu}}_{hc\perp}\left(0\right)\left(1+\gamma_{5}\right)h\left(0\right)\,,$
(5)

with the gauge invariant gluon field
$\mathcal{A^{\mu}}_{hc\perp}=W_{hc}^{\dagger}\left[D_{hc\perp}^{\mu}W_{hc}\right]=\mathcal{A}_{hc\perp}^{a\mu}t^{a}$.
Then the resolved contribution is represented by a time-ordered product of
this leading-power (LP) $A$-type current, with a subleading $\mathcal{L}_{\xi
q}^{\left(1\right)}$ Lagrangian, see right figure of Fig. 2.

Figure 1:  QCD diagram at LO, see text.

Figure 2:  SCET diagrams: direct (left) and resolved (right) contributions,
see text.

We can compare the different kinematics of the diagrams with $A$\- and
$B$-type currents in Figure 2. The external s-quark carries hardcollinear
momentum. Therefore the intermediate propagator is hard. This situation is
represented in SCET by the $B$-type current. When the momentum of the external
s-quark tends to zero, the propagator becomes anti-hardcollinear and cannot be
integrated out – it must be reproduced by a dynamical field in the low energy
EFT. This situation is represented in SCET by the time-order product of
subleading Lagrangian and the $A$-type current. The degeneracy in the EFT
description is the reason why the SCET develops divergencies in the
convolution integrals.

For the explicit matching of the two currents at leading order (LO) we find

${{\bf\color[rgb]{0.5078125,0.21484375,0.1484375}{C_{LO}^{B1}\left(m_{b},u\right)}}=(-1)\frac{\overline{u}}{u}\,\frac{m_{b}^{2}}{4\pi^{2}}\,\frac{G_{F}}{\sqrt{2}}\lambda_{t}\,C_{8g}}=(-1)\frac{\overline{u}}{u}{\bf\color[rgb]{0.5078125,0.21484375,0.1484375}C_{LO}^{A0}\left(m_{b}\right)}\,,$
(6)

where we use the hardcollinear momentum fraction
$u=\frac{\bar{n}p_{s}}{m_{b}}$ and
$\overline{u}=1-u=\frac{\bar{n}p_{r}}{m_{b}}$ with the hardcollinear momenta
of the strange quark and the gluon, $\bar{n}p_{s}$ and $\bar{n}p_{r}$ in the
direct contribution.

##  3 Bare factorisation theorem and endpoint divergences

The derivation of the factorisation theorem follows the standard approach [22,
23, 24]. We first perform the soft decoupling transformation [25], but we do
not use a new notation for the hardcollinear fields after decoupling. The
decay rate is obtained from the imaginary part of the current-current
correlator. The states factorise and thus allow taking matrix elements
separately in hardcollinear, anti-hardcollinear and soft sectors. The
hardcollinear matrix elements lead to jet functions. Integration of the anti-
hardcollinear fields at the amplitude level leads to antihardcollinear
(radiative) jet functions. We can then simplify the Dirac and colour
structure, and the remaining soft matrix element defines the shape functions.

For the direct contribution we find the bare factorisation theorem on the
operator level (see Figure 3 for a schematic description):

Figure 3:  Schematic description of the factorisation of the direct
contribution, see text.

$\frac{d\Gamma}{dE_{\gamma}}=\mathcal{N}_{B}\,\int_{0}^{1}du\,{\bf\color[rgb]{0.5078125,0.21484375,0.1484375}C^{B1}\left(m_{b},u\right)}\int_{0}^{1}du^{\prime}{\bf\color[rgb]{0.5078125,0.21484375,0.1484375}C^{B1*}\left(m_{b},u^{\prime}\right)}\int_{-p_{+}}^{\overline{\Lambda}}d\omega{\bf\color[rgb]{0.2421875,0.37890625,0.41015625}\,J\left(M_{B}\left(p_{+}+\omega\right),u,u^{\prime}\right)}{\bf\color[rgb]{0.421875,0.4296875,0.25390625}{\cal
S}\left(\omega\right)}$ (7)

$\mathcal{N}_{B}$ is a prefactor
333$\mathcal{N}_{B}=e^{2}Q_{s}^{2}\frac{E_{\gamma}}{2\pi}$. The hard function
${\bf\color[rgb]{0.5078125,0.21484375,0.1484375}C^{B1}}$ is given at LO in the
last section as a result of the matching of QCD on SCET. The hardcollinear jet
function is a genuine next-to-leading (NLP) object. We define it as a vacuum
matrix element of a product of hardcollinear fields:

$\displaystyle{\bf\color[rgb]{0.2421875,0.37890625,0.41015625}J\left(p^{2},u,u^{\prime}\right)}=\frac{(-1)}{2N_{c}}\,\frac{1}{2\pi}\,\int\frac{dtdt^{\prime}}{\left(2\pi\right)^{2}}\,\,d^{4}x\,\,e^{-im_{b}\left(ut-u^{\prime}t^{\prime}\right)+ipx}$
(8) $\displaystyle\text{Disc}\left[\left\langle
0\right|tr\left[\frac{1+\not{v}}{2}(1-\gamma_{5})\not{\mathcal{A}}_{hc\perp}\left(x\right)\gamma_{\perp}^{\nu}\chi_{hc}\left(t^{\prime}\bar{n}+x\right)\overline{\chi}_{hc}\left(t\bar{n}\right)\gamma_{\nu\perp}\not{\mathcal{A}}_{hc\perp}\left(0\right)\left(1+\gamma_{5}\right)\right]\left|0\right]\rangle\right]\,.$

The soft function - the leading power (LP) shape function - is defined as [26]

$\displaystyle{{\bf\color[rgb]{0.421875,0.4296875,0.25390625}{\cal
S}\left(\omega\right)}}=\frac{1}{2m_{B}}\int\frac{dt}{2\pi}e^{-i\omega
t}\left\langle
B\right|h\left(tn\right)S_{n}\left(tn\right)S_{n}^{\dagger}\left(0\right)h\left(0\right)\left|B\right\rangle\,.$
(9)

with soft Wilson lines denoted by $S$.

There is an endpoint divergence in the convolution of the hard matching
coefficients and the jet function within the direct contribution. At LO, one
finds explicitly with the hard matching coefficient
${{\bf\color[rgb]{0.5078125,0.21484375,0.1484375}{C_{LO}^{B1}\left(m_{b},u\right)}}}$
given in Eq. 6 and the LO jet function given by 444$A(\epsilon)$ is a function
which is finite for $\epsilon\to 0$.

$\displaystyle{\bf\color[rgb]{0.2421875,0.37890625,0.41015625}J\left(p^{2},u,u^{\prime}\right)}=C_{F}\frac{\alpha_{s}}{4\pi\,m_{b}}\theta(p^{2})\,A(\epsilon)\,\delta(u-u^{\prime})u^{1-\epsilon}(1-u)^{-\epsilon}\left(\frac{p^{2}}{\mu^{2}}\right)^{-\epsilon}$
(10)

that the convolution of both quantities. diverges in the $u\to 0$ limit:

$\int_{0}^{1}\frac{du}{u}\,\int_{0}^{1}\frac{du^{\prime}}{u^{\prime}}\delta(u-u^{\prime})u^{1-\epsilon}\approx\int_{0}^{1}\frac{du}{u^{1+\epsilon}}\,.$
(11)

Consequently, due to endpoint divergences, the bare factorisation formula of
the direct contribution is already invalid for the $d\to 4$ limit at LO.

For the resolved contribution we find the following factorisation theorem on
the operator level (for a schematic description of the factorisation of the
resolved contribution see Figure 4):

$\frac{d\Gamma}{dE_{\gamma}}=\mathcal{N}_{A}\left|{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{A0}\left(m_{b}\right)}\right|^{2}{\int_{-p_{+}}^{\overline{\Lambda}}}d\omega\,{{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J_{g}\left({m_{b}}\left(p_{+}+\omega\right)\right)}}\int d\omega_{1}\int
d\omega_{2}\,{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}\left(\omega_{1}\right)\,\overline{J}^{*}\left(\omega_{2}\right)}\,{\cal
S}\left(\omega,\omega_{1},\omega_{2}\right),$ (12)

with the prefactor $\mathcal{N}_{A}=\mathcal{N}_{B}\equiv\mathcal{N}$. The
hard function ${\bf\color[rgb]{0.5078125,0.21484375,0.1484375}C^{A}}$ is given
at LO in the last section as a result of the matching of QCD on SCET again. In
the hardcollinear sector, there are only gluon fields, so the standard LO
gluon jet function ${\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J_{g}(p^{2})}$ appears. The anti-hardcollinear jet function is defined on the
amplitude level because there is no energetic particle emitted in the anti-
hardcollinear directions besides the photons:

$\displaystyle\mathcal{O}_{T\xi q}$ $\displaystyle=i\int
d^{d}xT\left[\mathcal{L}_{\xi
q}\left(x\right),\mathcal{O}_{8g}^{A0}\left(0\right)\right]$
$\displaystyle=\int
d\omega\int\frac{dt}{2\pi}e^{-it\omega}\left[\overline{q_{s}}\right]_{\alpha}\left(tn\right){\color[rgb]{0.8984375,0.5,0.20703125}\bf\left[\overline{J}\left(\omega\right)\right]_{\alpha\beta}^{a\,\nu\mu}}Q_{s}\,\mathcal{B^{\nu}}_{\overline{hc}\perp}\left(0\right)\mathcal{A}_{hc\perp}^{\mu\,a}\left(0\right)\left[h\left(0\right)\right]_{\beta}.$
(13)

One can prove that the antihardcollinear jet function has a unique structure
to all orders in $\alpha_{s}$:

${{\color[rgb]{0.8984375,0.5,0.20703125}\bf\left[\overline{J}\left(\omega\right)\right]_{\alpha\beta}^{a\,\nu\mu}}}={{\color[rgb]{.75,.5,.25}\overline{J}\left(\omega\right)\,}}t^{a}\,\left[\gamma_{\perp}^{\nu}\gamma_{\perp}^{\mu}\frac{\not{\bar{n}}\not{n}}{4}\right]_{\alpha\beta},$
(14)

Finally, the soft function is given by

$\displaystyle{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}\left(u,t,t^{\prime}\right)}$
$\displaystyle={(d-2)^{2}}g_{s}^{2}\left\langle
B\right|\overline{h}\left(un\right)\left(1-\gamma_{5}\right)\left[S_{n}\left(un\right)t^{a}S_{n}^{\dagger}\left(un\right)\right]S_{\bar{n}}\left(un\right)S_{\bar{n}}^{\dagger}\left(t^{\prime}\bar{n}+un\right)$
(15)
$\displaystyle\frac{\not{n}\not{\bar{n}}}{4}q_{s}\left(t^{\prime}\bar{n}+un\right)\overline{q}_{s}\left(t\bar{n}\right)\frac{\not{\bar{n}}\not{n}}{4}S_{\bar{n}}\left(t\bar{n}\right)S_{\bar{n}}^{\dagger}\left(0\right)\left[S_{n}\left(0\right)t^{a}S_{n}^{\dagger}\left(0\right)\right]\left(1+\gamma_{5}\right)h\left(0\right)\left|B\right\rangle\,/\,(2m_{B})\,.$

The Fourier transform is defined as

${\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}\left(\omega,\omega_{1},\omega_{2}\right)}=\int\frac{du}{2\pi}e^{-iu\omega}\int\frac{dt}{2\pi}e^{-it\omega_{1}}\int\frac{dt^{\prime}}{2\pi}e^{it^{\prime}\omega_{2}}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}\left(u,t,t^{\prime}\right)}\,.$ (16)

Figure 4:  Schematic description of the factorisation of the resolved
contribution, see text.

At LO the factorization formula is well behaved as long as
$\omega_{1},\omega_{2}\sim\omega$ :

$\frac{d\Gamma}{dE_{\gamma}}=2\mathcal{N}_{A}\,\left|{{\bf\color[rgb]{0.5078125,0.21484375,0.1484375}C_{LO}^{A0}\left(m_{b}\right)}}\right|^{2}\int^{\overline{\Lambda}}_{-p_{+}}d\omega\,\,{{\bf\color[rgb]{0.2421875,0.37890625,0.41015625}\delta\left({m_{b}}\left(p_{+}+\omega\right)\right)}}\int_{-\infty}^{\infty}d\omega_{1}\int_{-\infty}^{\omega_{1}}d\omega_{2}\,{{\color[rgb]{0.8984375,0.5,0.20703125}\frac{1}{(\omega_{1}-i\epsilon)}\frac{1}{(\omega_{2}+i\epsilon)}}}{{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}\left(\omega,\omega_{1},\omega_{2}\right)}}.$ (17)

However, for $\omega_{1,2}\gg\omega$, the soft function can be shown to be
asymptotically constant, which leads to endpoint divergence in the convolution
integrals of jet and shape functions for large $\omega_{1,2}$.

In this limit, the soft strange quarks become hardcollinear and can be
decoupled from the soft gluons. This way, the soft function reduces to the LP
shape function, which also occurs in the direct contribution. This can be
shown explicitly taking into account that the $\omega_{1,2}\to\infty$ limit
corresponds to the $t,t^{\prime}\to 0$ limit in Eq. 3. In fact, the soft
function in the limit $\omega_{1,2}$, denoted by
${\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{{\cal S}}}$, can be
matched on the LP shape function with a perturbative kernel $K$:

$\displaystyle{\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{{\cal
S}}\left(\omega,\omega_{1},\omega_{2}\right)}=\int
d\omega^{\prime}K(\omega,\omega^{\prime},\omega_{1},\omega_{2}){\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\mathcal{S}}(\omega^{\prime})}\,.$
(18)

At LO, we find:

${\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{{\cal
S}}\left(\omega,\omega_{1},\omega_{2}\right)}=C_{F}A(\epsilon)\frac{\alpha_{s}}{(4\pi)}\,\,\omega_{1}^{1-\epsilon}\delta(\omega_{1}-\omega_{2})\int_{\omega}^{{\overline{\Lambda}}}d\omega^{\prime}\,{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}(\omega^{\prime})}\,\left(\frac{(\omega^{\prime}-\omega)}{\mu^{2}}\right)^{-\epsilon}\,,$
(19)

which includes the leading power shape function  ${\cal S}(\omega)$.

This makes the endpoint divergence in the convolution integral of jet and
shape function within the resolved contribution in the asymptotic limit
$\omega_{1,2}\to\infty$ manifest. In this limit, the resolved contribution at
LO can be written as:

$\displaystyle\frac{d\Gamma}{dE_{\gamma}}|^{\rm asy}_{A}$
$\displaystyle=2\mathcal{N}\,{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf|C_{LO}^{A0}(m_{b})|^{2}}\int_{-p_{+}}^{\overline{\Lambda}}d\omega{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J_{g}^{LO}(m_{b}(p_{+}+\omega))}\int_{m_{b}}^{\infty}d\omega_{1}{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}_{LO}(\omega_{1})}\int_{0}^{\omega_{1}}d\omega_{2}{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}^{*}_{LO}(\omega_{2})}\,\,{\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{{\cal
S}}(\omega,\omega_{1},\omega_{2})}$
$\displaystyle=\mathcal{N}{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf|C_{LO}^{A0}\left(m_{b}\right)|^{2}}\,\frac{\alpha_{s}C_{F}}{(4\pi)\,m_{b}}\,\frac{1}{\epsilon}A(\epsilon)\,\int_{-p_{+}}^{\overline{\Lambda}}d\omega\,{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}_{LO}(\omega^{\prime})}\left(\frac{m_{b}(\omega+p_{+})}{\mu^{2}}\right)^{-\epsilon}\,.$
(20)

##  4 Refactorisation to all orders in $\alpha_{s}$

Evaluating the direct contribution in the asymptotic limit $u,u^{\prime}\to 0$
(before integration)

$\displaystyle\frac{d\Gamma}{dE_{\gamma}}|^{u,u^{\prime}\to 0}_{B}$
$\displaystyle=-\mathcal{N}\,{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf\left|C_{LO}^{A0}\left(m_{b}\right)\right|^{2}}\,\frac{\alpha_{s}C_{F}}{(4\pi)\,m_{b}}\,\frac{1}{\epsilon}\,A(\epsilon)\int_{-p_{+}}^{\overline{\Lambda}}d\omega\,{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}_{LO}(\omega)}\left(\frac{m_{b}(\omega+p_{+})}{\mu^{2}}\right)^{-\epsilon}\,,$
(21)

and comparing with the asymptotic limit ($\omega_{1},\omega_{2}\to\infty$) of
the resolved contribution given in Eq. 20 one verifies at LO

$\frac{d\Gamma}{dE_{\gamma}}|^{\rm
asy}_{A}=(-1)\frac{d\Gamma}{dE_{\gamma}}|^{u,u^{\prime}\to 0}_{B}\,.$ (22)

Thus, the sum of the two terms is finite and equal to zero, and the endpoint
divergences cancel. This reflects the fact that in the limits $u\to 0$ and
$\omega_{1}\sim\omega_{2}\gg\omega$ the two terms of the subleading ${\cal
O}_{8}-{\cal O}_{8}$ contribution have the same structure. This LO result is a
special case of the all-order relation which can be formulated at the operator
level.

One can derive three refactorisation conditions which reflect this fact that
in the limits $u\sim u^{\prime}\ll 1$ and $\omega_{1}\sim\omega_{2}\gg\omega$
the resolved and the direct contribution have the same structure to all orders
in $\alpha_{s}$. The refactorisation relations are operatorial relations that
then guarantee the cancellation of endpoint divergences between the two terms
to all orders in $\alpha_{s}$ as we will see below.

  * •

We find that in the limit $u\to 0$, the matching coefficient can be further
factorised

${{\left\llbracket{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{B1}\left(m_{b},u\right)}\right\rrbracket=(-1){\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{A0}\left(m_{b}\right)}\,m_{b}\,{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}\left(um_{b}\right)\,,}}}$
(23)

where $\left\llbracket g(u)\right\rrbracket$ only denotes the leading term of
a function $g(u)$ in the limit $u\to 0$ and without any higher power
corrections in $u\ll 1$.

  * •

We find the new soft function $\widetilde{{\cal
S}}\left(\omega,\omega_{1},\omega_{2}\right)$ which corresponds to the
function  ${\cal S}\left(\omega,\omega_{1},\omega_{2}\right)$ in the limit
$\omega_{1}\sim\omega_{2}\gg\omega$. In this limit, we can consider the light
soft quarks to be hardcollinear. In this function  $\widetilde{{\cal
S}}\left(\omega,\omega_{1},\omega_{2}\right)$ higher power corrections in
$\omega/\omega_{1,2}$ are neglected.

  * •

We find that the jet function
$J\left(m_{b}\left(p_{+}+\omega\right),u,u^{\prime}\right)$ fulfills the
following relation in the limit $u\to 0$ and $u^{\prime}\to 0$:

$\int_{-p_{+}}^{\overline{\Lambda}}d\omega\left\llbracket{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J\left(m_{b}\left(p_{+}+\omega\right),u,u^{\prime}\right)}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}(\omega)}\right\rrbracket=\int_{-p_{+}}^{\overline{\Lambda}}d\omega{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J_{g}(m_{b}(p_{+}+\omega))}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{\cal
S}(\omega,m_{b}u,m_{b}u^{\prime})}\,,$ (24)

where the brackets indicate that the $u\to 0$ and $u^{\prime}\to 0$ limits
have to be taken.

Using these all-orders refactorisation conditions we then can cast the two
asymptotic subtraction terms into the following form:

$\displaystyle 0=\>$ $\displaystyle
2\mathcal{N}\left|{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{A0}}\left(m_{b}\right)\right|^{2}\int_{-p_{+}}^{\overline{\Lambda}}d\omega{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J_{g}\left(m_{b}\left(p_{+}+\omega\right)\right)}\int_{m_{b}}^{\infty}d\omega_{1}{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}\left(\omega_{1}\right)}\,\int_{0}^{\omega_{1}}d\omega_{2}{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}^{*}\left(\omega_{2}\right)}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{{\cal
S}}\left(\omega,\omega_{1},\omega_{2}\right)}$ $\displaystyle+\>$
$\displaystyle
2\mathcal{N}\int_{0}^{1}du\left\llbracket{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{B1}\left(m_{b},u\right)}\right\rrbracket\int_{u}^{1}du^{\prime}\left\llbracket{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{B1*}\left(m_{b},u^{\prime}\right)}\right\rrbracket\int_{-p_{+}}^{\overline{\Lambda}}d\omega\left\llbracket{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J\left(m_{b}\left(p_{+}+\omega\right),u,u^{\prime}\right)}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}(\omega)}\right\rrbracket\>.$ (25)

We then subtract these asymptotic terms from the sum of the two all-order bare
factorisation theorems for resolved and direct contribution derived in the
last section (see Eqs. 7 and 12) and obtain the endpoint finite factorisation
theorem.

$\displaystyle\frac{d\Gamma}{dE_{\gamma}}|_{A+B}$
$\displaystyle=2\mathcal{N}\int_{-p_{+}}^{\overline{\Lambda}}d\omega\bigg{\\{}{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J_{g}(m_{b}(p_{+}+\omega))}\left|{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{A0}\left(m_{b}\right)}\right|^{2}$ (26)
$\displaystyle\times\int_{-\infty}^{\infty}d\omega_{1}\,\int_{-\infty}^{\omega_{1}}d\omega_{2}{\color[rgb]{0.8984375,0.5,0.20703125}\bf\overline{J}(\omega_{1})\,\overline{J}^{*}(\omega_{2})}\left[{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}\left(\omega,\omega_{1},\omega_{2}\right)}-\theta(\omega_{1}-m_{b})\theta(\omega_{2}){\color[rgb]{0.421875,0.4296875,0.25390625}\bf\widetilde{\cal
S}(\omega,\omega_{1},\omega_{2})}\right]$
$\displaystyle+\int_{0}^{1}du\int_{u}^{1}du^{\prime}\,\Big{[}{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{B1}\left(m_{b},u\right)\,C^{B1*}\left(m_{b},u^{\prime}\right)}\,{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J\left(m_{b}\left(p_{+}+\omega\right),u,u^{\prime}\right)}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{{\cal
S}}\left(\omega\right)}$
$\displaystyle-\left\llbracket{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{B1}\left(m_{b},u\right)}\right\rrbracket\left\llbracket{\color[rgb]{0.5078125,0.21484375,0.1484375}\bf
C^{B1*}\left(m_{b},u^{\prime}\right)}\right\rrbracket\left\llbracket{\color[rgb]{0.2421875,0.37890625,0.41015625}\bf
J\left(m_{b}\left(p_{+}+\omega\right),u,u^{\prime}\right)}{\color[rgb]{0.421875,0.4296875,0.25390625}\bf{\cal
S}(\omega)}\right\rrbracket\Big{]}\bigg{\\}}\>,$

At this point, the convolutions integrals in the $A$\- and $B$-type
contributions are no longer divergent, and we can renormalise the functions
entering the factorisation theorem and take the limit $d\to 4$. In a final
step, we are able to show that endpoint reshuffling and renormalisation
effectively commute [17].

Acknowledgements  TH is supported by the Cluster of Excellence “Precision
Physics, Fundamental Interactions, and Structure of Matter" (PRISMA+ EXC
2118/1) funded by the German Research Foundation (DFG) within the German
Excellence Strategy (Project ID 390831469), RS is supported by the United
States Department of Energy under Grant Contract DESC0012704. R.S. and TH
thank the organisers of SCET23 and FPCP23 for their excellent organisation and
generous hospitality.

## References

  * [1] T. Hurth, “Present status of inclusive rare B decays,” Rev. Mod. Phys. 75 (2003) 1159 [arXiv:hep-ph/0212304].
  * [2] T. Hurth and M. Nakao, “Radiative and Electroweak Penguin Decays of B Mesons,” Ann. Rev. Nucl. Part. Sci. 60 (2010) 645 [arXiv:1005.1224 [hep-ph]].
  * [3] T. Hurth and F. Mahmoudi, “New physics search with flavour in the LHC era,” Rev. Mod. Phys. 85 (2013) 795 [arXiv:1211.6453 [hep-ph]].
  * [4] E. Kou et al. [Belle-II], “The Belle II Physics Book,” PTEP 2019 (2019) no.12, 123C01 [erratum: PTEP 2020 (2020) no.2, 029201] [arXiv:1808.10567 [hep-ex]].
  * [5] S. J. Lee, M. Neubert and G. Paz, “Enhanced Non-local Power Corrections to the $\bar{B}\to X_{s}\gamma$ Decay Rate,” Phys. Rev. D 75, 114005 (2007) [arXiv:hep-ph/0609224].
  * [6] M. Benzke, S. J. Lee, M. Neubert and G. Paz, “Factorization at Subleading Power and Irreducible Uncertainties in $\bar{B}\to X_{s}\gamma$ Decay,” JHEP 1008 (2010) 099 [arXiv:1003.5012 [hep-ph]].
  * [7] T. Hurth, M. Fickinger, S. Turczyk and M. Benzke, “Resolved Power Corrections to the Inclusive Decay $\bar{B}\to X_{s}\ell^{+}\ell^{-}$,” Nucl. Part. Phys. Proc. 285-286 (2017) 57 [arXiv:1711.01162 [hep-ph]].
  * [8] M. Benzke, T. Hurth and S. Turczyk, “Subleading power factorization in $\bar{B}\to{X}_{s}{\ell}^{+}{\ell}^{-}$,” JHEP 1710 (2017) 031 [arXiv:1705.10366 [hep-ph]].
  * [9] A. Gunawardana and G. Paz, “Reevaluating uncertainties in $\overline{B}\to X_{s}\gamma$ decay,” JHEP 11, 141 (2019) [arXiv:1908.02812 [hep-ph]].
  * [10] M. Benzke and T. Hurth, “Resolved $1/m_{b}$ contributions to $\bar{B}\to X_{s,d}\ell^{+}\ell^{-}$ and $\bar{B}\to X_{s}\gamma$,” Phys. Rev. D 102, 114024 (2020) [arXiv:2006.00624 [hep-ph]].
  * [11] M. Benzke and T. Hurth, “Recent progress on the nonlocal power corrections to the inclusive penguin decays B → X s $\gamma$ andB → X s $\ell$$\ell$,” [arXiv:2303.06447 [hep-ph]].
  * [12] M. Misiak et al., “Updated NNLO QCD predictions for the weak radiative $B$-meson decays,” Phys. Rev. Lett. 114 (2015) no.22, 221801 [arXiv:1503.01789 [hep-ph]].
  * [13] T. Huber, T. Hurth and E. Lunghi, “Inclusive $\bar{B}\to{X}_{s}{\ell}^{+}{\ell}^{-}$ : complete angular analysis and a thorough study of collinear photons,” JHEP 1506 (2015) 176 [arXiv:1503.04849 [hep-ph]].
  * [14] T. Huber, T. Hurth, J. Jenkins, E. Lunghi, Q. Qin and K. K. Vos, “Long distance effects in inclusive rare $B$ decays and phenomenology of $\bar{B}\to X_{d}\ell^{+}\ell^{-}$,” JHEP 1910 (2019) 228 [arXiv:1908.07507 [hep-ph]].
  * [15] T. Huber, T. Hurth, J. Jenkins, E. Lunghi, Q. Qin and K. K. Vos, JHEP 10 (2020), 088 doi:10.1007/JHEP10(2020)088 [arXiv:2007.04191 [hep-ph]].
  * [16] M. Benzke and T. Hurth, “Recent progress on the nonlocal power corrections to the inclusive penguin decays B → X s $\gamma$ andB → X s $\ell$$\ell$,” [arXiv:2303.06447 [hep-ph]].
  * [17] T. Hurth and R. Szafron, Nucl. Phys. B 991 (2023), 116200 doi:10.1016/j.nuclphysb.2023.116200 [arXiv:2301.01739 [hep-ph]].
  * [18] Z. L. Liu and M. Neubert, “Factorization at subleading power and endpoint-divergent convolutions in $h\to\gamma\gamma$ decay,” JHEP 04, 033 (2020) [arXiv:1912.08818 [hep-ph]].
  * [19] Z. L. Liu, B. Mecaj, M. Neubert and X. Wang, “Factorization at subleading power and endpoint divergences in $h\to\gamma\gamma$ decay. Part II. Renormalization and scale evolution,” JHEP 01, 077 (2021) [arXiv:2009.06779 [hep-ph]].
  * [20] M. Beneke, M. Garny, S. Jaskiewicz, R. Szafron, L. Vernazza and J. Wang, JHEP 10 (2020), 196 doi:10.1007/JHEP10(2020)196 [arXiv:2008.04943 [hep-ph]].
  * [21] M. Beneke, M. Garny, S. Jaskiewicz, J. Strohm, R. Szafron, L. Vernazza and J. Wang, “Next-to-leading power endpoint factorization and resummation for off-diagonal “gluon” thrust,” JHEP 07, 144 (2022) [arXiv:2205.04479 [hep-ph]].
  * [22] I. Moult, I. W. Stewart and G. Vita, “Subleading Power Factorization with Radiative Functions,” JHEP 11, 153 (2019) [arXiv:1905.07411 [hep-ph]].
  * [23] G. Vita, “QCD Beyond Leading Power,” [arXiv:2008.10606 [hep-ph]].
  * [24] S. E. Jaskiewicz, “Factorization and Resummation at Subleading Powers,” (Doctoral dissertation, Technische Universität München)
  * [25] C. W. Bauer, D. Pirjol and I. W. Stewart, “Soft collinear factorization in effective field theory,” Phys. Rev. D 65, 054022 (2002) [arXiv:hep-ph/0109045 [hep-ph]].
  * [26] M. Neubert, “Analysis of the photon spectrum in inclusive B —$>$ X(s) gamma decays,” Phys. Rev. D 49 (1994), 4623-4633 [arXiv:hep-ph/9312311 [hep-ph]].


# Unified Linearization-based Nonlinear Filtering

Anton Kullberg1, Isaac Skog2, , and Gustaf Hendeby1, .

This work was partially supported by the Wallenberg AI, Autonomous Systems and
Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.
1Dept. Electrical Engineering, Linköping University, Linköping, Sweden  2Dept.
Electrical Engineering, Uppsala University, Uppsala, Sweden  Email:
{anton.kullberg<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This letter shows that the following three classes of recursive state
estimation filters: standard filters, such as the extended Kalman filter;
iterated filters, such as the iterated unscented Kalman filter; and
dynamically iterated filters, such as the dynamically iterated posterior
linearization filters; can be unified in terms of a general algorithm. The
general algorithm highlights the strong similarities between specific
filtering algorithms in the three filter classes and facilitates an in-depth
understanding of the pros and cons of the different filter classes and
algorithms. We end with a numerical example showing the estimation accuracy
differences between the three classes of filters when applied to a nonlinear
localization problem.

††publicationid: pubid:

##  I Introduction

State estimation in nonlinear dynamical systems has been extensively studied
in a wide variety of research fields. Typical approaches employ some form of
linearization-based approximate inference, which we focus on here. These
approaches linearize the nonlinear model locally (in each time instance) to
employ the Kalman filter, which is the optimal estimator in the _mean-squared
error_ (mse) sense [1]. Analytical linearization leads to the _extended Kalman
filter_ (ekf), while sigma-point filters, such as the _unscented Kalman
filter_ (ukf), _cubature Kalman filter_ (ckf), and similar can be thought of
as statistical linearization filters [2, 1, 3]. Statistical linearization
filters also includes the Gaussian particle filter [4, 5].

The estimation accuracy of linearization-based filters highly depend on the
point (distribution in the statistical case) about which the models are
linearized. Typically, the linearization point (distribution) is chosen to be
the mean (distribution) of the current state estimate. With a large error in
the state estimate, this can lead to compounding errors which, in the worst
case, may cause the filter to diverge. To alleviate this problem, several
variants of iterated filters have been developed, such as the _iterated
extended Kalman filter_ (iekf), the _iterated unscented Kalman filter_ (iukf),
and the _iterated posterior linearization filter_ (iplf) [6, 7, 8, 9, 10].
These types of filters essentially iterate the measurement update, each time
re-linearizing the measurement model with the “latest” iterate. The efforts in
iterated filtering have primarily been focused on finding a better
linearization point for the measurement model, which has been motivated by the
fact that nonlinearities in the measurement model affect the resulting state
estimate to a greater extent than nonlinearities in the transition model.

Iterated filters have also been generalized to improve the linearization point
for the transition model [11, 12]. These algorithms, which we refer to as
dynamically iterated filters, are essentially iterated one-step fixed-lag
smoothers that extract information from the measurement at time $k$ to improve
the linearization of the transition model at time $k-1$. Examples of such
algorithms are the _dynamically iterated extended Kalman filter_ (diekf),
_dynamically iterated unscented Kalman filter_ (diukf), and the _dynamically
iterated posterior linearization filter_ (diplf) [11, 12].

In this letter, we seek to provide a “general” algorithm, from which all of
the aforementioned filter algorithms can be derived as special cases. In this
way, we aim to clarify and highlight the strong similarities between different
linearization-based filtering algorithms. Thus, the contribution is a
unification of linearization-based filters in a single general algorithm,
encompassing analytically and statistically linearized, as well as iterated
and non-iterated filters. We also illustrate the performance differences
between the three kinds of filter classes on an acoustic localization problem.

##  II Background

For clarity, we here present analytical and statistical linearization within a
common framework. The well-known Kalman filter and _Rauch-Tung-Stribel_ (rts)
smoother equations are also recapitulated.

###  II-A Kalman Smoother

Assume an affine state-space model with additive Gaussian noise, of the form

$\displaystyle\mathbf{x}_{k+1}$
$\displaystyle=\mathbf{A}_{\mathbf{f}}\mathbf{x}_{k}+\mathbf{b}_{\mathbf{f}}+\tilde{\mathbf{w}}_{k},$
$\displaystyle\tilde{\mathbf{w}}_{k}$
$\displaystyle\sim\mathcal{N}(\tilde{\mathbf{w}}_{k};\mathbf{0},\mathbf{Q}\\!+\\!\bm{\Omega}_{\mathbf{f}})$
(1a) $\displaystyle\mathbf{y}_{k}$
$\displaystyle=\mathbf{A}_{\mathbf{h}}\mathbf{x}_{k}+\mathbf{b}_{\mathbf{h}}+\tilde{\mathbf{e}}_{k},$
$\displaystyle\tilde{\mathbf{e}}_{k}$
$\displaystyle\sim\mathcal{N}(\tilde{\mathbf{e}}_{k};\mathbf{0},\mathbf{R}\\!+\\!\bm{\Omega}_{\mathbf{h}}).$
(1b)

Here, $\mathbf{x}_{k},~{}\mathbf{y}_{k},~{}\tilde{\mathbf{w}}_{k}$ and
$\tilde{\mathbf{e}}_{k}$ denote the state, the measurement, the process noise
and the measurement noise at time $k$, respectively. Lastly, assume that
$\mathbf{x}_{k}\in\mathcal{X},\forall k$, where $\mathcal{X}$ is some set,
typically $\mathbb{R}^{n_{x}}$, and that $\tilde{\mathbf{w}}_{k}$ and
$\tilde{\mathbf{e}}_{k}$ are mutually independent. Note that usually,
$\bm{\Omega}_{\mathbf{f}}=\bm{\Omega}_{\mathbf{h}}=\mathbf{0}$. For this
model, the (affine) Kalman smoother update equations are given by Alg. 1,
where subscript k|k denotes an estimate at time $k$ given measurements up
until time $k$ and $K$ is the final time [13].

###  II-B Analytical and Statistical Linearization

Given a nonlinear model

$\mathbf{z}=\mathbf{g}(\mathbf{x}),$

we wish to find an affine representation

$\mathbf{g}(\mathbf{x})\approx\mathbf{A}\mathbf{x}+\mathbf{b}+\eta,$ (2)

with $\eta\sim\mathcal{N}(\eta;\mathbf{0},\bm{\Omega})$. In this affine
representation, there are three free parameters, $\mathbf{A},\mathbf{b}$, and
$\bm{\Omega}$. Analytical linearization through first-order Taylor expansion
selects the parameters as

$\mathbf{A}=\frac{d}{d\mathbf{x}}\mathbf{g}(\mathbf{x})|_{\mathbf{x}=\bar{\mathbf{x}}},\quad\mathbf{b}=\mathbf{g}(\mathbf{x})|_{\mathbf{x}=\bar{\mathbf{x}}}-\mathbf{A}\bar{\mathbf{x}},\quad\bm{\Omega}=\mathbf{0},$
(3)

where $\bar{\mathbf{x}}$ is the point about which the function
$\mathbf{g}(\mathbf{x})$ is linearized. Note that $\bm{\Omega}=\mathbf{0}$
essentially implies that the linearization is assumed to be error free.

Statistical linearization instead linearizes w.r.t. a distribution
$p(\mathbf{x})$. Assuming that
$p(\mathbf{x})=\mathcal{N}(\mathbf{x};\hat{\mathbf{x}},\mathbf{P})$,
statistical linearization selects the affine parameters as

$\displaystyle\mathbf{A}$ $\displaystyle=\Psi^{\top}\mathbf{P}^{-1}$ (4a)
$\displaystyle\mathbf{b}$
$\displaystyle=\bar{\mathbf{z}}-\mathbf{A}\hat{\mathbf{x}}$ (4b)
$\displaystyle\bm{\Omega}$
$\displaystyle=\Phi-\mathbf{A}\mathbf{P}\mathbf{A}^{\top}$ (4c)
$\displaystyle\bar{\mathbf{z}}$
$\displaystyle=\mathbb{E}[\mathbf{g}(\mathbf{x})]$ (4d) $\displaystyle\Psi$
$\displaystyle=\mathbb{E}[(\mathbf{x}-\hat{\mathbf{x}})(\mathbf{g}(\mathbf{x})-\bar{\mathbf{z}})^{\top}]$
(4e) $\displaystyle\Phi$
$\displaystyle=\mathbb{E}[(\mathbf{g}(\mathbf{x})-\bar{\mathbf{z}})(\mathbf{g}(\mathbf{x})-\bar{\mathbf{z}})^{\top}],$
(4f)

where the expectations are taken w.r.t. $p(\mathbf{x})$. The major difference
from analytical linearization is that ${\bm{\Omega}\neq 0}$, which implies
that the error in the linearization is captured.

Typically, the expectations in 4 are not analytically tractable and thus,
practically, one often resorts to some numerical integration technique.

Alg. 1 Kalman smoother

  1. 1.

Time update (TU)

$\displaystyle\hat{\mathbf{x}}_{k+1|k}$
$\displaystyle=\mathbf{A}_{\mathbf{f}}\hat{\mathbf{x}}_{k|k}+\mathbf{b}_{\mathbf{f}}$
(5a) $\displaystyle\mathbf{P}_{k+1|k}$
$\displaystyle=\mathbf{A}_{\mathbf{f}}\mathbf{P}_{k|k}\mathbf{A}_{\mathbf{f}}^{\top}+\mathbf{Q}+\bm{\Omega}_{\mathbf{f}}$
(5b)

  2. 2.

Measurement update (MU)

$\displaystyle\hat{\mathbf{x}}_{k|k}$
$\displaystyle=\hat{\mathbf{x}}_{k|k-1}+\mathbf{K}_{k}(\mathbf{y}_{k}-\mathbf{A}_{\mathbf{h}}\hat{\mathbf{x}}_{k|k-1}-\mathbf{b}_{\mathbf{h}})$
(6a) $\displaystyle\mathbf{P}_{k|k}$
$\displaystyle=\mathbf{P}_{k|k-1}-\mathbf{K}_{k}\mathbf{A}_{\mathbf{h}}\mathbf{P}_{k|k-1}$
(6b) $\displaystyle\mathbf{K}_{k}$
$\displaystyle\triangleq\mathbf{P}_{k|k-1}\mathbf{A}_{\mathbf{h}}^{\top}(\mathbf{A}_{\mathbf{h}}\mathbf{P}_{k|k-1}\mathbf{A}_{\mathbf{h}}^{\top}+\mathbf{R}+\bm{\Omega}_{\mathbf{h}})^{-1}$
(6c)

  3. 3.

Smoothing step (S)

$\displaystyle\hat{\mathbf{x}}_{k|K}$
$\displaystyle=\hat{\mathbf{x}}_{k|k}+\mathbf{G}_{k}(\hat{\mathbf{x}}_{k+1|K}-\hat{\mathbf{x}}_{k+1|k})$
(7a) $\displaystyle\mathbf{P}_{k|K}$
$\displaystyle=\mathbf{P}_{k|k}+\mathbf{G}_{k}(\mathbf{P}_{k+1|K}-$
$\displaystyle\quad~{}\mathbf{A}_{\mathbf{f}}\mathbf{P}_{k|k}\mathbf{A}_{\mathbf{f}}^{\top}-\mathbf{Q}-\bm{\Omega}_{\mathbf{f}})\mathbf{G}_{k}^{\top}$
(7b) $\displaystyle\mathbf{G}_{k}$
$\displaystyle\triangleq\mathbf{P}_{k|k}\mathbf{A}_{\mathbf{f}}^{\top}(\mathbf{A}_{\mathbf{f}}\mathbf{P}_{k|k}\mathbf{A}_{\mathbf{f}}^{\top}+\mathbf{Q}+\bm{\Omega}_{\mathbf{f}})^{-1}$
(7c)

##  III Problem Formulation

To set the stage for the unification of the different filter algorithms, the
general state estimation problem is described here from a probabilistic
viewpoint. To that end, consider a discrete-time state-space model (omitting a
possible input $\mathbf{u}_{k}$ for notational brevity) given by

$\displaystyle\mathbf{x}_{k+1}$
$\displaystyle=\mathbf{f}(\mathbf{x}_{k})+\mathbf{w}_{k},$ $\displaystyle
p(\mathbf{w}_{k})$
$\displaystyle=\mathcal{N}(\mathbf{w}_{k};\mathbf{0},\mathbf{Q})$ (8a)
$\displaystyle\mathbf{y}_{k}$
$\displaystyle=\mathbf{h}(\mathbf{x}_{k})+\mathbf{e}_{k},$ $\displaystyle
p(\mathbf{e}_{k})$
$\displaystyle=\mathcal{N}(\mathbf{e}_{k};\mathbf{0},\mathbf{R}).$ (8b)

Note that 8a and 8b can equivalently be written as a _transition density_ and
a _measurement density_ as

$\displaystyle p(\mathbf{x}_{k+1}|\mathbf{x}_{k})$
$\displaystyle=\mathcal{N}(\mathbf{x}_{k+1};\mathbf{f}(\mathbf{x}_{k}),\mathbf{Q})$
(9a) $\displaystyle p(\mathbf{y}_{k}|\mathbf{x}_{k})$
$\displaystyle=\mathcal{N}(\mathbf{y}_{k};\mathbf{h}(\mathbf{x}_{k}),\mathbf{R}).$
(9b)

Further, the initial state distribution is assumed to be given by

$p(\mathbf{x}_{0})=\mathcal{N}(\mathbf{x}_{0};\hat{\mathbf{x}}_{0|0},\mathbf{P}_{0|0}).$
(10)

Given the transition and measurement densities and a sequence of measurements
$\mathbf{y}_{1:k}=\\{\mathbf{y}_{i}\\}_{i=1}^{k}$, the filtering problem
consists of computing the marginal posterior of the state at time $k$. This
can be done via the Bayesian recursions

$\displaystyle p(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})$
$\displaystyle=\int_{\mathcal{X}}\\!p(\mathbf{x}_{k}|\mathbf{x}_{k-1})p(\mathbf{x}_{k-1}|\mathbf{y}_{1:k-1})d\mathbf{x}_{k-1}$
(11a) $\displaystyle p(\mathbf{x}_{k}|\mathbf{y}_{1:k})$
$\displaystyle=\frac{p(\mathbf{y}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})}{\mathbf{Z}_{k}}$
(11b) $\displaystyle\mathbf{Z}_{k}$
$\displaystyle=\int_{\mathcal{X}}p(\mathbf{y}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})d\mathbf{x}_{k}.$
(11c)

In the case where $\mathbf{f}$ and $\mathbf{h}$ are linear, the (analytical)
solution is given by the Kalman filter [1].

In the general case, the marginal posteriors can not be computed analytically.
Inspecting 11, there are two integrals that require attention. We turn first
to the Chapman-Kolmogorov equation 11a. Assuming that
$p(\mathbf{x}_{k-1}|\mathbf{y}_{1:k-1})$ is Gaussian, 11a has a closed form
solution given by 5, _if_ $p(\mathbf{x}_{k}|\mathbf{x}_{k-1})$ is Gaussian and
8a is affine. Therefore, as 9a is Gaussian, we seek an affine approximation of
the transition function $\mathbf{f}$ as

$\mathbf{f}(\mathbf{x}_{k-1})\approx\mathbf{A}_{\mathbf{f}}\mathbf{x}_{k-1}+\mathbf{b}_{\mathbf{f}}+\eta_{\mathbf{f}},$
(12)

with
$p(\eta_{\mathbf{f}})=\mathcal{N}(\eta_{\mathbf{f}};\mathbf{0},\bm{\Omega}_{\mathbf{f}})$.
Hence, the transition density $p(\mathbf{x}_{k}|\mathbf{x}_{k-1})$ is
approximated by $q(\mathbf{x}_{k}|\mathbf{x}_{k-1})$ as

$q(\mathbf{x}_{k}|\mathbf{x}_{k-1})=\mathcal{N}(\mathbf{x}_{k};\mathbf{A}_{\mathbf{f}}\mathbf{x}_{k-1}+\mathbf{b}_{\mathbf{f}},\mathbf{Q}+\bm{\Omega}_{\mathbf{f}}).$
(13)

If $\mathbf{A}_{\mathbf{f}},\mathbf{b}_{\mathbf{f}}$, and
$\bm{\Omega}_{\mathbf{f}}$ are chosen to be the analytical linearization of
$\mathbf{f}$ about the mean of the posterior
$p(\mathbf{x}_{k-1}|\mathbf{y}_{1:k-1})$, the ekf time update is recovered
through 5. Similarly, statistical linearization about
$p(\mathbf{x}_{k-1}|\mathbf{y}_{1:k-1})$ recovers the sigma-point filter time
updates. This yields an approximate predictive distribution
$q(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})$, which can then be used to approximate
the second integral of interest (and subsequently, the posterior at time $k$).
Explicitly, the second integral is approximated by

$\mathbf{Z}_{k}\approx\int_{\mathcal{X}}p(\mathbf{y}_{k}|\mathbf{x}_{k})q(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})d\mathbf{x}_{k}.$
(14)

Similarly to 12, 14 has a closed form solution _if_
$p(\mathbf{y}_{k}|\mathbf{x}_{k})$ is Gaussian and 8b is affine. Thus, as 9b
is Gaussian, we seek an affine approximation of the measurement function
$\mathbf{h}$ as

$\mathbf{h}(\mathbf{x}_{k})\approx\mathbf{A}_{\mathbf{h}}\mathbf{x}_{k}+\mathbf{b}_{\mathbf{h}}+\eta_{\mathbf{h}},$
(15)

with
$p(\eta_{\mathbf{h}})=\mathcal{N}(\eta_{\mathbf{h}};\mathbf{0},\bm{\Omega}_{\mathbf{h}})$.
Hence, the measurement density $p(\mathbf{y}_{k}|\mathbf{x}_{k})$ is
approximated by $q(\mathbf{y}_{k}|\mathbf{x}_{k})$ as

$q(\mathbf{y}_{k}|\mathbf{x}_{k})=\mathcal{N}(\mathbf{y}_{k};\mathbf{A}_{\mathbf{h}}\mathbf{x}_{k}+\mathbf{b}_{\mathbf{h}},\mathbf{R}+\bm{\Omega}_{\mathbf{h}}),$
(16)

which leads to an analytically tractable integral. With 13 and 16, the
(approximate) marginal posterior 11 is now given by

$q(\mathbf{x}_{k}|\mathbf{y}_{1:k})=\frac{q(\mathbf{y}_{k}|\mathbf{x}_{k})q(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})}{\int_{\mathcal{X}}q(\mathbf{y}_{k}|\mathbf{x}_{k})q(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})d\mathbf{x}_{k}},$
(17)

which is analytically tractable and given by 6. Note that analytical
linearization of 15 about the mean of $q(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})$
recovers the ekf measurement update, whereas statistical linearization
recovers the sigma-point measurement update(s).

The quality of the approximate marginal posterior 17 directly depends on the
quality of the approximations 13 and 16. The quality of 13 and 16 in turn
directly depends on the choice of linearization points or densities, which is
typically chosen to be the approximate predictive and previous approximate
posterior distributions. This choice is of course free and iterated filters,
such as the iekf, iukf, and iplf have been proposed to improve the
approximation 16 [14, 15, 7, 8]. These filters essentially iterate the
measurement update to find an approximate posterior
$q^{i}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$, which is used to re-linearize the
function $\mathbf{h}$ to produce a new approximation
$q^{i+1}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$. Iterated filters were recently
generalized to dynamically iterated filters, which improve both the
approximation 16, as well as the approximation 13 [12, 11]. Dynamically
iterated filters are essentially one-step iterated fixed-lag smoothers that
produce both a better posterior approximation
$q^{i+1}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$ as well as a smoothed approximation
$q^{i+1}(\mathbf{x}_{k-1}|\mathbf{y}_{1:k})$.

Next, we describe a unification of all of these filters in terms of one
general algorithm, encompassing all possible variants of filters based on
either analytical or statistical linearization.

##  IV Unified Linearization-based Filtering

TUMUS$q(\mathbf{x}_{k-1}|\mathbf{y}_{1:k-1})$$q^{i}(\mathbf{x}_{k}|\mathbf{y}_{1:k-1})$$q^{i}(\mathbf{x}_{k}|\mathbf{y}_{1:k})$$q^{i}(\mathbf{x}_{k-1}|\mathbf{y}_{1:k})$Iterated
FiltersDynamically Iterated Filters Figure 1: Schematic illustration of
linearization-based filters. Iterated filters re-linearize the measurement
update (MU). Dynamically iterated filters also re-linearize the time update
(TU) through a smoothing step (S).

We propose a unified linearization-based filtering algorithm that encapsulates
a wide variety of existing algorithms. The main idea behind the unification is
that all linearization-based filters may be thought of as a single general
algorithm, reduces to various special cases depending on specific
implementation choices. All of the filters are essentially centered around the
three key steps 5, 6 and 7. They differ only in the choice of linearization
strategy, as well as in which steps of the general (approximative affine)
Kalman filter/smoother that are repeated or not. The general algorithm is
presented in Alg. 2 and encompasses standard linearization-based filters,
iterated filters, and dynamically iterated filters. For clarity, it is also
illustrated schematically in Fig. 1. Note that the unified algorithm is
purposefully restricted to algorithms that only require access to the latest
measurement $\mathbf{y}_{k}$, which, e.g., excludes the L-scan iplf [16].

The linearization choices, which are assumed to be the same for all of the
steps in the general filter algorithm, and the specific filter algorithms
these choices lead to, are summarized in Table I. In Table I, iterating either
the measurement update (MU), both the time update (TU) and MU or none (–) is
captured vertically, and the choice of particular linearization strategy
horizontally. Choosing analytical linearization inevitably leads to some form
of “extended” version, i.e., either the ekf, iekf, or diekf. Statistical
linearization is a bit more nuanced for two reasons. Firstly, it encapsulates
a wide variety of algorithms, depending on the particular chosen statistical
linearization, be it exact or approximated by, e.g., some form of cubature
such as the ckf or ukf. Note that we use ckf as a collective term for any
statistically linearized Kalman filter based on sigma points, such as the
smart sampling Kalman filter [17], the spherical simplex-radial ckf [18], or
the multiple quadrate Kalman filter [19]. Secondly, iterated versions of
statistical linearization filters fall into two distinct categories, iukf
style that “freezes” the covariance update until the last iterate [7], or iplf
style that continuously updates the covariance matrix – essentially changing
the sigma point spread each iteration [15]. In Table I, the “frozen”
statistical linearization based filters are summarized by, e.g., the ickf and
iukf, but should be read as encapsulating any imaginable version of
statistical linearization where the resulting filter has an update structure
similar to that of the iekf/diekf, i.e., with a “delayed” covariance update.
Note that the “freezing” or “delayed” behaviour of the iukf/diukf is not
explicitly defined in the algorithm Alg. 2 but amounts to setting
$\mathbf{P}_{k-1|k}^{i+1}:=\mathbf{P}_{k-1|k}^{i}$ after Alg. 5 and
$\mathbf{P}_{k|k}^{i+1}:=\mathbf{P}_{k|k}^{i}$ after Alg. 4 until the last
iteration.

TABLE I: Specific algorithms that Alg. 2 reduces to given different
assumptions

| Linearization Method
---|---
Iteration Type | Analytical | Statistical | Frozen Statistical
– |  ekf [1] |  ckf [20], ukf [2] | –
MU |  iekf [6] |  iplf [15] |  ickf, iukf [9]
TU & MU |  diekf [12, 21] |  diplf [12, 11] |  dickf, diukf [12]

Alg. 2 Linearization-based filter

1:$\hat{\mathbf{x}}_{k-1|k-1},\mathbf{P}_{k-1|k-1},\mathbf{y}_{k}$

2:Alg. 3 with
$\hat{\mathbf{x}}_{k-1|k}^{0}=\hat{\mathbf{x}}_{k-1|k-1},\mathbf{P}_{k-1|k}^{0}=\mathbf{P}_{k-1|k-1}$

3:Alg. 4 with
$\hat{\mathbf{x}}_{k|k}^{0}=\hat{\mathbf{x}}^{1}_{k|k-1},\mathbf{P}_{k|k}^{0}=\mathbf{P}^{1}_{k|k-1}$

4:if Standard Filter then

5: return $\hat{\mathbf{x}}_{k|k},\mathbf{P}_{k|k}$

6:else if Dynamically Iterated Filter then

7: Alg. 5

8:end if

9:$i\leftarrow 0$

10:while not converged do

11: if Dynamically Iterated Filter then

12: Algs. 3, 4, and 5

13: else if Iterated Filter then

14: Alg. 4

15: end if

16: $i\leftarrow i+1$

17:end while

18:return $\hat{\mathbf{x}}_{k|k}^{i},\mathbf{P}_{k|k}^{i}$ (and
$\hat{\mathbf{x}}_{k-1|k}^{i},\mathbf{P}_{k-1|k}^{i}$)

Alg. 3 Linearized time update

1:$\hat{\mathbf{x}}_{k-1|k-1},\mathbf{P}_{k-1|k-1},\hat{\mathbf{x}}_{k-1|k}^{i},\mathbf{P}_{k-1|k}^{i}$

2: Linearize $\mathbf{f}$ about
$\hat{\mathbf{x}}_{k-1|k}^{i},\mathbf{P}_{k-1|k}^{i}$through 3 or 4

3:Compute $\hat{\mathbf{x}}_{k|k-1}^{i+1},\mathbf{P}_{k|k-1}^{i+1}$ by 5

4:return $\hat{\mathbf{x}}_{k|k-1}^{i+1},\mathbf{P}_{k|k-1}^{i+1}$

Alg. 4 Linearized measurement update

1:$\hat{\mathbf{x}}_{k|k-1},\mathbf{P}_{k|k-1},\hat{\mathbf{x}}_{k|k}^{i},\mathbf{P}_{k|k}^{i},\mathbf{y}_{k}$

2: Linearize $\mathbf{h}$ about
$\hat{\mathbf{x}}_{k|k}^{i},\mathbf{P}_{k|k}^{i}$ through 3 or 4

3:Compute $\hat{\mathbf{x}}_{k|k}^{i+1},\mathbf{P}_{k|k}^{i+1}$ by 6

4:return $\hat{\mathbf{x}}_{k|k}^{i+1},\mathbf{P}_{k|k}^{i+1}$

Alg. 5 Linearized smoothing step

1:$\hat{\mathbf{x}}_{k-1|k-1},\mathbf{P}_{k-1|k-1},\hat{\mathbf{x}}_{k|k},\mathbf{P}_{k|k},\hat{\mathbf{x}}_{k-1|k}^{i},\mathbf{P}_{k-1|k}^{i}$

2: Linearize $\mathbf{f}$ about
$\hat{\mathbf{x}}_{k-1|k}^{i},\mathbf{P}_{k-1|k}^{i}$through 3 or 4

3:Compute $\hat{\mathbf{x}}_{k-1|k}^{i+1},\mathbf{P}_{k-1|k}^{i+1}$ by 7

4:return $\hat{\mathbf{x}}_{k-1|k}^{i+1},\mathbf{P}_{k-1|k}^{i+1}$

##  V Numerical Example

To demonstrate the application of the three types of filters, we consider a
localization problem modeled by a nonlinear state-space model. To keep the
results uncluttered, we only consider analytical linearization and focus our
comparison on the ekf, iekf, and diekf.

We consider a target maneuvering in a plane and describe the target state
using the state vector
$\mathbf{x}_{k}=\begin{bmatrix}p^{x}_{k}&v^{x}_{k}&p^{y}_{k}&v^{y}_{k}&\omega_{k}\end{bmatrix}^{\top}$.
Here, $p^{x}_{k},~{}p^{y}_{k},~{}v^{x}_{k}$, and $v^{y}_{k}$ are the Cartesian
coordinates and velocities of the target, respectively. Further, $\omega_{k}$
is the turn rate. The transition model is given by

$\mathbf{x}_{k+1}=\mathbf{f}(\mathbf{x}_{k})+\mathbf{w}_{k},$ (18)

where

$\mathbf{f}(\mathbf{x}_{k})=\begin{bmatrix}1&\frac{\sin(T\omega_{k})}{\omega_{k}}&0&-\frac{(1-\cos(T\omega_{k}))}{\omega_{k}}&0\\\
0&\cos(T\omega_{k})&0&-\sin(T\omega_{k})&0\\\
0&\frac{(1-\cos(T\omega_{k}))}{\omega_{k}}&1&\frac{\sin(T\omega_{k})}{\omega_{k}}&0\\\
0&\sin(T\omega_{k})&0&\cos(T\omega_{k})&0\\\
0&0&0&0&1\end{bmatrix}\mathbf{x}_{k},$

and $T$ is the sampling period. Further,
$\mathbf{w}_{k}\sim\mathcal{N}(\mathbf{w}_{k};\mathbf{0},\mathbf{Q})$ is the
process noise at time $k$, with

$\mathbf{Q}=\mathrm{blkdiag}\left(\begin{bmatrix}q_{1}\frac{T^{3}}{3}&q_{1}\frac{T^{2}}{2}\\\
q_{1}\frac{T^{2}}{2}&q_{1}T\end{bmatrix},\begin{bmatrix}q_{1}\frac{T^{3}}{3}&q_{1}\frac{T^{2}}{2}\\\
q_{1}\frac{T^{2}}{2}&q_{1}T\end{bmatrix},q_{2}\right),$

where $q_{1}$ and $q_{2}$ are tunable parameters of the model.

The target emits a known sound pulse at a rate of $T=$1.5\text{\,}\mathrm{s}$$
that is picked up by a set of four microphones. With this, we construct time-
difference-of-arrival (tdoa) observations, where each observation $i$ is
modeled as

$\mathbf{y}_{k}^{i}=r^{1}_{k}-r^{i}_{k}+\mathbf{e}_{k},\quad i=1,\dots,3$ (19)

where
$r^{i}_{k}\triangleq\lVert\begin{bmatrix}p^{x}_{k}&p^{y}_{k}\end{bmatrix}^{\top}-s^{i}\rVert$,
and $s^{i}$ denotes the 2D position of the $i$:th microphone. Further,
$\mathbf{e}_{k}\sim\mathcal{N}(\mathbf{0},\mathbf{R})$, where $\mathbf{R}$ has
been computed through a static calibration experiment.

We set $q_{1}=10^{-j},~{}q_{2}=10^{-l}$ and let
$j=-6,\dots,0,~{}l=-5,\dots,0$, and sweep over all such pairs, i.e., 42
different process noise configurations. For each configuration we compute the
rmse against a ground truth trajectory, obtained from a high-precision IR-
marker positioning system.

The positional rmse per noise configuration is presented in Fig. 2. Clearly,
the diekf performs the best overall and is non-divergent in most cases. Here,
divergence corresponds to an rmse higher than $1\text{\,}\mathrm{m}$. As the
process noise is increased, the difference between the algorithms decreases,
but the iterative procedure of the iekf and diekf is still clearly beneficial.

Figure 2: Positional rmse for the ekf as blue dots, iekf as orange crosses and
diekf as green squares. Each subplot corresponds to a different value of
$q_{1}$, indicated by the text in each subplot. An rmse higher than
approximately $1\text{\,}\mathrm{m}$ corresponds to a “divergent” filter based
on visual inspection of resulting estimate trajectories and is left out of the
plots.

##  VI Conclusion

A unifying view of linearization–based nonlinear filtering algorithms has been
presented. It facilitates a comprehensive understanding of the commonalities
and relationships between linearization–based standard, iterated, and
dynamically iterated filters. The presented algorithm is simple, easy to
implement, and encompasses a wide range of existing filtering algorithms.
Lastly, the three classes of unified filtering algorithms were compared in a
nonlinear localization problem, where the dynamically iterated filters were
shown to be more resilient to poor process noise parameter tuning.

## References

  * [1] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” _J. of Basic Eng._ , vol. 82, no. 1, pp. 35–45, Mar. 1960.
  * [2] S. Julier, J. Uhlmann, and H. Durrant-Whyte, “A New Approach for Filtering Nonlinear Systems,” in _Proc. 1995 Am. Control Conf. - ACC95_ , vol. 3, Seattle, WA, USA, 1995, pp. 1628–1632.
  * [3] T. Lefebvre, H. Bruyninckx, and J. De Schuller, “Comment on ”A new method for the nonlinear transformation of means and covariances in filters and estimators” [with authors’ reply],” _IEEE Trans. Automat. Contr._ , vol. 47, no. 8, pp. 1406–1409, Aug. 2002.
  * [4] J. Kotecha and P. Djuric, “Gaussian particle filtering,” _IEEE Trans. Signal Process._ , vol. 51, no. 10, pp. 2592–2601, Oct. 2003.
  * [5] Y. Wu, X. Hu, D. Hu, and M. Wu, “Comments on ”Gaussian particle filtering”,” _IEEE Trans. Signal Process._ , vol. 53, no. 8, pp. 3350–3351, Aug. 2005.
  * [6] W. Denham and S. Pines, “Sequential estimation when measurement function nonlinearity is comparable to measurement error,” in _Guidance Control Conference_ , ser. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics, Aug. 1965.
  * [7] M. A. Skoglund, F. Gustafsson, and G. Hendeby, “On Iterative Unscented Kalman Filter using Optimization,” in _2019 22th Int. Conf. Inf. Fusion FUSION_ , Ottawa, ON, Canada, Jul. 2019, pp. 1–8.
  * [8] R. Zhan and J. Wan, “Iterated Unscented Kalman Filter for Passive Target Tracking,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 43, no. 3, pp. 1155–1163, Jul. 2007.
  * [9] G. Sibley, G. Sukhatme, and L. Matthies, “The Iterated Sigma Point Kalman Filter with Applications to Long Range Stereo,” in _Robot. Sci. Syst. II_ , Aug. 2006.
  * [10] Á. F. García-Fernández, L. Svensson, M. R. Morelande, and S. Särkkä, “Posterior Linearization Filter: Principles and Implementation Using Sigma Points,” _IEEE Trans. Signal Process._ , vol. 63, no. 20, pp. 5561–5573, Oct. 2015.
  * [11] M. Raitoharju, R. Hostettler, and S. Sarkka, “Posterior Linearisation Filter for Non-linear State Transformation Noises,” in _2022 25th Int. Conf. Inf. Fusion FUSION_. Linköping, Sweden: IEEE, Jul. 2022, pp. 1–6.
  * [12] A. Kullberg, I. Skog, and G. Hendeby, “Iterated Filters for Nonlinear Transition Models,” Apr. 2023.
  * [13] S. Särkkä, _Bayesian Filtering and Smoothing_ , 1st ed. Cambridge University Press, Sep. 2013.
  * [14] A. H. Jazwinski, _Stochastic Processes and Filtering Theory_. Academic Press, 1970.
  * [15] Á. F. García-Fernández, L. Svensson, and M. R. Morelande, “Iterated Statistical Linear Regression for Bayesian Updates,” in _17th Int. Conf. Inf. Fusion FUSION_ , Salamanca, Spain, Jul. 2014, pp. 1–8.
  * [16] Á. F. García-Fernández, L. Svensson, and S. Särkkä, “Iterated Posterior Linearization Smoother,” _IEEE Trans. Automat. Contr._ , vol. 62, no. 4, pp. 2056–2063, Apr. 2017.
  * [17] J. Steinbring and U. D. Hanebeck, “LRKF Revisited: The Smart Sampling Kalman Filter (S2KF),” _J. Adv. Inf. Fusion_ , vol. 9, no. 2, pp. 106–123, 2014.
  * [18] S. Wang, J. Feng, and C. K. Tse, “Spherical Simplex-Radial Cubature Kalman Filter,” _IEEE Signal Process. Lett._ , vol. 21, no. 1, pp. 43–46, Jan. 2014.
  * [19] P. Closas, C. Fernandez-Prades, and J. Vila-Valls, “Multiple Quadrature Kalman Filtering,” _IEEE Trans. Signal Process._ , vol. 60, no. 12, pp. 6125–6137, Dec. 2012.
  * [20] I. Arasaratnam and S. Haykin, “Cubature Kalman Filters,” _IEEE Trans. Automat. Contr._ , vol. 54, no. 6, pp. 1254–1269, Jun. 2009.
  * [21] R. P. Wishner, J. A. Tabaczynski, and M. Athans, “On the Estimation of the State of Noisy Nonlinear Multivariable Systems,” in _IFAC Symp. Multivariable Control Syst._ , Düsseldorf, 1968.