# Multidimensional Projection Filters via Automatic Differentiation and
Sparse-Grid Integration

Muhammad Fuady Emzir<EMAIL_ADDRESS>Zheng Zhao  Simo Särkkä
Control and Instrumentation Engineering Department, King Fahd University of
Petroleum and Minerals, Dhahran, Saudi Arabia  Interdisciplinary Center of
Smart Mobility and Logistics, King Fahd University of Petroleum and Minerals,
Dhahran, Saudi Arabia  Uppsala University, Uppsala, Sweden  Aalto University,
Espoo, Finland

###### Abstract

The projection filter is a technique for approximating the solutions of
optimal filtering problems. In projection filters, the Kushner–Stratonovich
stochastic partial differential equation that governs the propagation of the
optimal filtering density is projected to a manifold of parametric densities,
resulting in a finite-dimensional stochastic differential equation. Despite
the fact that projection filters are capable of representing complicated
probability densities, their current implementations are limited to Gaussian
family or unidimensional filtering applications. This work considers a
combination of numerical integration and automatic differentiation to
construct projection filter algorithms for more generic problems.
Specifically, we provide a detailed exposition of this combination for the
manifold of the exponential family, and show how to apply the projection
filter to multidimensional cases. We demonstrate numerically that based on
comparison to a finite-difference solution to the Kushner–Stratonovich
equation and a bootstrap particle filter with systematic resampling, the
proposed algorithm retains an accurate approximation of the filtering density
while requiring a comparatively low number of quadrature points. Due to the
sparse-grid integration and automatic differentiation used to calculate the
expected values of the natural statistics and the Fisher metric, the proposed
filtering algorithms are highly scalable. They therefore are suitable to many
applications in which the number of dimensions exceeds the practical limit of
particle filters, but where the Gaussian-approximations are deemed
unsatisfactory.

###### keywords:

Projection filter , nonlinear filter , automatic differentiation , sparse-grid
integration

We use automatic differentiation and sparse-grid integration to automate the
construction of the projection filter.

We present methods for constructing projection filters for multidimensional
filtering problems using a non-Gaussian parametric density.

We show that the practical performance of the filter is comparable to the
particle filter and finite difference based solutions to the
Kushner–Stratonovich equation.

An open-source implementation of the method is available.

##  1 Introduction

Obtaining the statistically best estimate of an unobserved stochastic process
based on a stochastic observation process generated by it is known as the
optimal filtering problem (Jazwinski, 1970; Liptser and Shiryaev, 2001; Bain
and Crisan, 2009; Särkkä, 2013). The optimal filtering problem is equivalent
to determining the evolution of the conditional probability density of the
state of the unobserved process given the observation process up to the
current time. When the stochastic processes are modeled as Itô type of
stochastic differentiation equations (SDEs), under suitable regularity
conditions, the evolution of the conditional density is described by a
stochastic partial differential equation (SPDE) known as the
Kushner–Stratonovich equation (Jazwinski, 1970; Liptser and Shiryaev, 2010;
Kushner, 1967a, c; Wonham, 1963). The Kushner–Stratonovich equation is a
nonlinear SPDE that has a complicated structure which makes it hard to solve
(Ceci and Colaneri, 2013). Zakai (1969) introduced an alternative
representation to the Kushner–Stratonovich equation, resulting in a linear
SPDE for an unnormalized conditional density of the state. The solutions of
the Zakai equation are, nevertheless, still hard to obtain as they are
generally infinite-dimensional except when the Lie algebra naturally
associated with it is finite-dimensional (Maurel and Michel, 1984; Brigo et
al., 1999). In most cases, we can only hope for approximations of the
solutions.

Among approximations to the solution of the Kushner–Stratonovich equation is
the projection filter which has been proposed in the late 80s (Hanzon and Hut,
1991; Brigo et al., 1995, 1998, 1999). Essentially, it projects the dynamics
of the conditional probability density given by the Kushner–Stratonovich
equation onto a finite-dimensional manifold of parametric densities, resulting
in a finite-dimensional stochastic differential equation (SDE). So far, in
practice, the use of the projection filter has been very limited; outside of
the Gaussian density family, projection filters have mainly been developed for
unidimensional dynamical systems (Brigo et al., 1995, 1999; Azimi-Sadjadi and
Krishnaprasad, 2005; Armstrong and Brigo, 2016; Koyama, 2018; Tronarp and
Särkkä, 2019). One of the main issues in projection filter is that the filter
equation is model-specific, and depends on the chosen parametric family and
the natural statistics used. This means that in order to construct the filter,
one needs to derive the filtering equation analytically. This manual procedure
is prone to mistakes (Armstrong and Brigo, 2016). Another challenge is that in
order to propagate the parameters of SDEs, the log-partition function and
several expected statistics need to be computed at every step via numerical
integration. These numerical integrals are difficult to calculate since one
needs to take care of different integrands with different supports
simultaneously. To reduce the computational burden, so far, projection
filter’s numerical implementation has relied upon a recursive procedure to
compute the expected values of the natural statistics and the Fisher
information matrix (Brigo et al., 1995). This recursion is, in general, only
feasible for unidimensional problems.

Previously, Armstrong and Brigo (2016) proposed a functional ring as a
symbolic layer for abstraction of the operations in a function space including
addition, multiplication, scalar multiplication, and differentiation. Whilst
this abstraction provides some degree of automation to derivation of the
filter equation, numerical integrations for computing the expected values of
the natural statistics are still necessary. These integrations can be time-
consuming and repetitive, as the expected values of the natural statistics
need to be calculated by taking their different supports into account.
Moreover, the implementation in Armstrong and Brigo (2016) is restricted to
unidimensional problems.

Although the projection filter has not recently been used in many
applications, there are recent theoretical results for it. For example,
Armstrong and Brigo (2019, 2018); Armstrong et al. (2018) study alternative
projections than the standard Stratonovich projection. Different distance
measures than the Hellinger distance have also been recently investigated in
Armstrong and Brigo (2016). Apart from the exponential family, the mixture
family has been considered in Armstrong and Brigo (2016), and it was proven
that if the mixture components are chosen to be the Galerkin basis functions,
then the projection filter under the $L^{2}$ metric is equivalent to a
Galerkin method applied to the Kushner–Stratonovich SPDE. The projection
filter is also gaining popularity as an alternative approach to solve finite-
dimensional quantum filtering problems (van Handel and Mabuchi, 2005; Gao et
al., 2019, 2020b, 2020a). In these formulations, the original projection
filter formulation is modified using a different type of projection from Amari
and Nagaoka (2000, Chapter 7). For a recent review of these topics, we refer
the reader to Brigo (2022).

There is also another line of research that offers methods to approximate the
filtering densities via random samples. Among these methods is particle
filtering which is also known as sequential Monte Carlo (Doucet, 2001).
Although these methods can approximate the optimal filters for almost all
types of filtering problems, to avoid particle weight collapse, the number of
particles in these methods often scales exponentially with the dimension of
the problem (Snyder et al., 2008; Poterjoy, 2015; Beskos et al., 2017).
Another class of filters are so-called assumed density filters (Kushner,
1967b) which include sigma-point filters such as unscented Kalman filter
(Julier and Uhlmann, 1997) and other Gaussian approximations of the optimal
filter (Ito and Xiong, 2000; Särkkä and Sarmavuori, 2013). These filters can
be seen as special cases of projection filters (Brigo et al., 1999). For
comprehensive surveys of optimal filtering, we refer the reader to Bain and
Crisan (2009); Crisan and Rozovskii (2011); Särkkä (2013); Candy (2016).

The contribution of this paper is to introduce an alternative way to implement
projection filters via polynomial expansions and automatic differentiation.
Automatic differentiation is a method for determining the exact derivatives of
a function (to machine precision) by tracing basic mathematical operations
involved in the function (Rall, 1981; Griewank and Walther, 2008). It has a
computing cost that is proportional to the cost of function evaluation,
depending on whether forward, backward, or hybrid accumulations are utilized,
as well as on the dimensionality of the Jacobian. In particular, we use a
Chebyshev polynomial expansions for unidimensional problems and sparse-grid
integrations for higher dimensions in order to efficiently compute the log-
partition function. The numerical integration is only used for computing the
log-partition function, and the rest of the computations in the filter are
computed via automatic differentiation, which significantly simplifies the
filter algorithms and eliminates the repetitive numerical integrations for
calculating the expected values of the statistics of the exponential family at
hand. In order to do so, we use a smooth bijection from the state sample space
$\mathcal{X}$ to the canonical open hypercube $(-1,1)^{\otimes d}$
($\mathcal{X}$ is assumed to be fixed, see Section 3). This setting ensures
that the integrations can be performed in a fixed domain. To demonstrate the
use of the projection filter, we focus on a class of models, where the
coefficients of the dynamic model SDE can be represented by polynomials, the
natural statistics are monomials of the states, and the drift coefficient of
the observation model SDE is in the span of the natural statistics. In all,
our method allows for an automatic implementation of projection filters. That
is to say, one only needs to specify the model (symbolically) and a bijection
function to carry out the filter computations – the manual derivation for the
filter equation is no more needed.

Using sparse-grid techniques in optimal filtering is not entirely novel. For
instance, there are Gaussian approximation-based filters that use sparse-grid
quadratures (Winshcel and Krätzig, 2010; Jia et al., 2012; Baek and Bang,
2013; Radhakrishnan et al., 2016; Singh et al., 2018). Furthermore, Bao et al.
(2014) employ an adaptive sparse-grid and a finite difference approach to
solve the Zakai equation. We propose an alternative approximation for
nonlinear optimal filtering problems compared to Bao et al. (2014).
Specifically, we employ a sparse-grid method to calculate the log-partition
function that is used in the projection filter approximation. Compared to Bao
et al. (2014), our approach requires much less variables to be propagated at
each time step. This is due to that we only propagate the natural parameters
of the exponential family and not the quadrature points. Moreover, our
approach requires the quadrature points to be calculated only once at the
beginning of the computations, whereas the approach of Bao et al. (2014)
requires the quadrature points to be calculated at each time step at various
grid-refinement phases. Despite the fact that adaptive sparse grids may have
fewer nodes while retaining excellent accuracy, the grid-refinement technique
is problem-dependent and requires additional computing steps to generate new
quadrature points. These procedures involve repeated evaluations of the
integrand (Bungartz and Griebel, 2004; Ma and Zabaras, 2009). Introducing
these additional calculations for every time step incurs a substantial
computational expense. Similarly to our approach, Bao et al. (2014) employ a
hypercube as the sparse-grid domain. Nevertheless, the placement of the
hypercube is adjusted each time using a heuristic technique based on a finite
number of samples derived from calculating the quadrature point trajectory
using the state equations. Our approach is free of these extra computations.

Even if we can compute the moments and Fisher metric perfectly, the projection
filter still has local projection errors. Fortunately, the estimates for these
local projection errors are well-established for certain classes of state-
space models and parametric families (with different natural statistics)
(Brigo et al., 1999, Propositions 6.2 and 6.4). We believe that our
approximation will be useful in problems where obtaining approximate solutions
to Kushner–Stratonovich equations with particle filters is deemed too costly,
numerically unstable, or time-consuming, and where the Gaussian approximation
is not adequate.

The article is organized as follows. In Section 2, we present the problem
formulation. In particular, we define the class of state-space models (i.e.,
the SDEs of the dynamic and observation models) for which the proposed
projection filter can be constructed. Section 3 describes how to approximate
the log-partition function via polynomial interpolation and automatic
differentiation for unidimensional problems. Section 4 describes the use of
sparse-grid methods to compute the log-partition function in multidimensional
cases. Section 5 describes methods to obtain the expected values of the
natural statistics and the Fisher metric via automatic differentiation in
multidimensional problems. In Section 6, we conduct numerical experiments to
show the efficiency of the proposed projection filter in both unidimensional
and multidimensional filtering problems. Section 7 concludes our findings111We
have made the code for the approach described in this work accessible for
download. The source code is available at https://to_be_revealed..

##  2 Projection filter construction for an exponential family

In this paper, we consider optimal filtering problems consisting of the
following dynamic and observation models:

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

where $x_{t}\in\mathbb{R}^{d_{x}}$ and $y_{t}\in\mathbb{R}^{d_{y}}$. The
processes $\left\\{W_{t}\right\\}_{t\geq 0},\left\\{V_{t}\right\\}_{t\geq 0}$
are Brownian motions 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. We also define
$\alpha(x_{t},t)\coloneqq\sigma(x_{t},t)Q_{t}\sigma(x_{t},t)^{\top}$. For the
sake of exposition, from now on, we assume that $R_{t}=I$ for all $t\geq 0$.
For any vector $v_{t}\in\mathbb{R}^{n}$ that is a function of time $t$, we
denote the $k$-th element of $v_{t}$ by $v_{t,k}$. Similarly, for any matrix
$X_{t}\in\mathbb{R}^{m\times n}$ that is a function of time $t$, we denote the
$i,j$ entry of $X_{t}$ as $X_{t,(i,j)}$.

Let the $m$-dimensional exponential family EM$(c)$ be the set of probability
densities of the form
$p_{\theta}(x)=\exp(\sum_{i=1}^{m}c_{i}(x)\theta_{i}-\psi(\theta)),x\in\mathcal{X}$,
where $\theta\in\Theta\subset\mathbb{R}^{m}$ and
$\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$ are the natural parameters and the
sample space of the state $x$, respectively. We also assume that the natural
statistics $\left\\{c_{i}\right\\}_{i=1,\ldots,m}$ are linearly independent.
The log-partition function $\psi$ is given by

$\displaystyle\psi(\theta)$
$\displaystyle=\log\int_{\mathcal{X}}\exp\left(\sum_{i=1}^{m}c_{i}(x)\theta_{i}\right)dx.$
(2)

The inclusion of the log-partition function $\psi$ inside the exponential
ensures that the density $p_{\theta}$ is normalized. Moreover, by definition
of the exponential family, $\Theta\subset\mathbb{R}^{m}$ is selected such that
for any $\theta\in\Theta$, $\psi<\infty$ (Calin and Udrişte, 2014, Section
1.4).

Under suitable regularity conditions to guarantee the well-definedness of the
system (i.e., the solutions to the SDEs in Equation (1) exist, and the
projection from the space of probability densities to the finite-dimensional
manifold exists), the projection filter using Stratonovich projection to
finite-dimensional manifold is given by (Brigo et al., 1999):

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

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

In what follows, we restrict the class of the dynamics in (1) with the
following assumptions.

###### Assumption 2.1

The functions $f,h,$ and $\sigma$ are polynomial functions of $x_{t}$.

###### Assumption 2.2

The natural statistics $\left\\{c_{i}\right\\}$ are selected as monomials of
the elements of $x_{t}$.

###### Assumption 2.3

Each component of $h_{t}$, that is, $h_{t,k}$ is in the span of
$\left\\{1,c_{1},\ldots,c_{m}\right\\}$. that is, there exists collection of
real numbers $\big{\\{}\lambda^{(i)}_{k}\colon
i=1,\ldots,m,\,k=1,\ldots,d_{x}\big{\\}}$ such that

$\displaystyle h_{t,k}$
$\displaystyle=\lambda_{k}^{(0)}+\sum_{i=1}^{m}\lambda_{k}^{(i)}c_{i},$

for all $k=1,\ldots,d_{x}$.

Under Assumption 2.3, the resulting exponential family of probability
densities is known as the EM$(c^{\ast})$ family defined in Brigo et al.
(1999). For this specific family EM$(c^{\ast})$, the projection filter
equation reduces to (Brigo et al., 1999, Theorem 6.3)

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

Furthermore, (4) can be simplified further with Assumptions 2.1 and 2.2. The
idea is to express the terms in Equation (4) as products of matrices and
expectations of statistics under these assumptions. Specifically, we obtain
$\mathcal{L}_{t}[c]-\frac{1}{2}|h_{t}|^{2}c=a_{0}+A_{0}\tilde{c}$ and
$\frac{1}{2}|h_{t}|^{2}=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_{x}}\rightarrow\mathbb{R}^{m_{h}}$ is the vector of the
remaining monomials of $x_{t}$ that are not included in $c$. Therefore, we can
express (4) as

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

where the matrix of the polynomial coefficients $M$ is given by

$\displaystyle M$ $\displaystyle=A_{0}+\eta(\theta_{t})b_{h}^{\top}.$ (6)

Denoting $\tilde{\eta}:=\mathbb{E}_{\theta_{t}}[\tilde{c}]$ we arrive at

$\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}.$ (7)

Equation (7) can also be used with other exponential families where the
natural statistics are not necessarily monomials. For instance, (7) can be
used in the case where the natural statistics are allowed to have negative
powers of $x_{t}$. However, the automation of implementation for the
polynomials is much easier because the polynomials are closed under addition,
multiplication, differentiation, and integration. During the implementation,
we use SymPy (Meurer et al., 2017), which is a symbolic mathematical package
for Python programming language, to automatically calculate the matrices and
vectors $a_{0}$, $b_{0}$, and $A_{0}$ in (6) and (7).

In order to solve Equation (7), we need to compute the expectation
$\tilde{\eta}$ and the Fisher metric $g$. These quantities can be computed
from the log-partition function. Specifically, Amari (1985); Calin and Udrişte
(2014); Brigo et al. (1999) show that

$\displaystyle\eta_{i}=\partial_{i}\psi\quad\text{and}\quad
g_{ij}(\theta)=\partial_{i}\partial_{j}\psi$ (8a)

hold for all $i,j=1,\ldots,d_{x}$. In the equation above, we denote
$\partial_{i}\psi:=\partialderivative{\psi}{\theta_{i}}$ which is a commonly
used shorthand in information geometry. Provided that the log-partition
function $\psi$ is evaluated numerically, evaluating $\eta$ and $g$ is
straightforward with any numerical package that supports automatic
differentiation. We will show in the subsequent sections that the
approximation error arising from these computations is predictable.
Furthermore, for the remaining statistics $c_{h}$, a simple trick using
automatic differentiation can be used to calculate their expectations. We will
describe this in Section 5.

##  3 Calculating log-partition function in unidimensional case

The previous discussion shows that under Assumptions 2.1, 2.2, and 2.3, after
re-arranging the filter equation to (7), to construct the projection filter,
it remains to calculate $\psi$. In this section, we will describe how to
compute $\psi$ efficiently in the unidimensional case ($d_{x}=1$), while the
next section covers the multidimensional case. We will start by defining a
bijection that will be used in the numerical integration of the log-partition
function (2). In the most general setting, the support of the parametric
probability $p(\cdot,\theta)$ depends on $\theta$. However, Amari (1985)
pointed out that this condition poses a substantial difficulty to analysis.
Therefore, following Amari (1985, Section 2.1), we assume that the support of
$p(\cdot,\theta)$ is uniform for all $\theta\in\Theta$. For notational
simplicity let us set $\Omega=(-1,1)$, and denote the class of bounded
variation functions on $\Omega$ as $BV(\bar{\Omega})$, where $\bar{\Omega}$ is
the closure of $\Omega$.

Equation (2) can be written as

$\displaystyle\psi(\theta)$
$\displaystyle=\log\int_{\Omega}z(\tilde{x})d\tilde{x},$ (9a) $\displaystyle
z(\tilde{x})$
$\displaystyle=\phi^{\prime}(\tilde{x})\exp\left(\sum_{i}c_{i}(\phi(\tilde{x}))\theta_{i}\right).$
(9b)

Let $L^{2}_{w}(\Omega)$ be the set of square integrable functions on $\Omega$
with respect to a positive weight function
$w(\tilde{x})=(1-\tilde{x})^{-1/2}$. Let us then define inner product by
$\expectationvalue{f,g}:=\int_{\Omega}w(\tilde{x})f(\tilde{x})g(\tilde{x})d\tilde{x}$.
The space $L^{2}_{w}(\Omega)$ equipped with the inner product is a Hilbert
space, and the Chebyshev polynomials are a complete orthogonal basis of
$L_{w}^{2}(\Omega)$. Assuming that $z\in L^{2}_{w}(\Omega)$, we can expand
$z(\tilde{x})$ on the Chebyshev polynomial basis:

$\displaystyle z(\tilde{x})$
$\displaystyle=\sum_{i=0}^{\infty}\hat{z}_{i}(\theta)T_{i}(\tilde{x}),$
$\displaystyle\hat{z}_{i}(\theta)$ $\displaystyle=\expectationvalue{T_{i},z}.$
(10)

In this equation, $T_{i}$ is the $i$-th Chebyshev polynomial of the first
kind, and $\hat{z}_{i}(\theta)$ is the corresponding Chebyshev coefficient.
Since the Chebysev series is related to the Fourier series by substituting
$\tilde{x}=\cos(\omega)$, $\hat{z}_{i}(\theta)$ can be numerically computed
using the fast Fourier Transform (FFT). Let
$S_{n}^{T}f:=\sum_{j=0}^{n}\expectationvalue{T_{j},f}T_{j}$ be the truncated
Chebyshev sum of $f$ up to $n$ first terms. Moreover, let $\psi(\theta)^{(n)}$
be the approximation of $\psi$ where $z$ is replaced by $S_{n}^{T}z$:

$\displaystyle\psi(\theta)^{(n)}$
$\displaystyle:=\log\left(\sum_{j=0}^{n}\int_{\Omega}\hat{z}_{j}(\theta)T_{j}(\tilde{x})d\tilde{x}\right),$
(11)
$\displaystyle=\log\left(\sum_{j=0}^{n}\hat{z}_{j}(\theta)\bar{T}_{j}\right),$
$\displaystyle\bar{T}_{j}$
$\displaystyle:=\int_{\Omega}T_{j}(\tilde{x})d\tilde{x}=\begin{cases}\frac{-1^{j}+1}{1-j^{2}}&,j\neq
0,\\\ 0&,j=0.\end{cases}$ (12)

In the following proposition, we show that the approximated log-partition
function $\psi(\theta)^{(n)}$ converges to $\psi$.

###### Proposition 3.1

Suppose that the natural statistics $\left\\{c_{i}\right\\}$ and the bijection
$\phi$ are chosen so that the $z$ in (9b) belongs to $L^{2}_{w}(\Omega)$.
Then, the approximated log-partition $\psi(\theta)^{(n)}$ converges to
$\psi(\theta)$ as $n$ approaches infinity.

###### Proof 1

From the definition of $\psi$ and $\psi^{(n)}$ we obtain:

$\displaystyle\absolutevalue{\psi(\theta)-\psi(\theta)^{(n)}}$
$\displaystyle=\absolutevalue{\log\left[\dfrac{\int_{\Omega}zd\tilde{x}}{\int_{\Omega}S_{n}^{T}zd\tilde{x}}\right]}$
$\displaystyle=\absolutevalue{\log\left[\dfrac{\int_{\Omega}z-S_{n}^{T}zd\tilde{x}+\int_{\Omega}S_{n}^{T}zd\tilde{x}}{\int_{\Omega}S_{n}^{T}zd\tilde{x}}\right]}$
$\displaystyle\leq\absolutevalue{\log\left[\dfrac{\int_{\Omega}\absolutevalue{z-S_{n}^{T}z}d\tilde{x}+\int_{\Omega}S_{n}^{T}zd\tilde{x}}{\int_{\Omega}S_{n}^{T}zd\tilde{x}}\right]}$
$\displaystyle\leq\absolutevalue{\log\left[\dfrac{\left(\int_{\Omega}1^{2}d\tilde{x}\right)^{1/2}\left(\int_{\Omega}\absolutevalue{z-S_{n}^{T}z}^{2}d\tilde{x}\right)^{1/2}+\int_{\Omega}S_{n}^{T}zd\tilde{x}}{\int_{\Omega}S_{n}^{T}zd\tilde{x}}\right]}$
$\displaystyle\leq\absolutevalue{\log\left[\dfrac{\sqrt{2}\left(\int_{\Omega}\absolutevalue{z-S_{n}^{T}z}^{2}w(\tilde{x})d\tilde{x}\right)^{1/2}+\int_{\Omega}S_{n}^{T}zd\tilde{x}}{\int_{\Omega}S_{n}^{T}zd\tilde{x}}\right]}$
$\displaystyle=\absolutevalue{\log\left[\dfrac{\sqrt{2}\norm{z-S^{T}_{n}z}_{2}+\int_{\Omega}S_{n}^{T}zd\tilde{x}}{\int_{\Omega}S_{n}^{T}zd\tilde{x}}\right]},$

where $\norm{\cdot}_{2}$ is $L^{2}_{w}\left(\Omega\right)$ norm. The result
follows immediately.

###### Remark 3.1

Proposition 3.1 shows that approximation (2) converges. Indeed, since $z$ in
(9b) is likely to be several times continuously differentiable or even
analytic, even a stronger convergence result can be expected since the
Chebyshev approximation error $\absolutevalue{f-S_{n}^{T}f}$ has
$\mathcal{O}(r^{-n})$ rate when $f$ is analytic for some $r>1$ (Mason and
Handscomb, 2003, Theorem 5.16). In this case
$\absolutevalue{\psi(\theta)-\psi(\theta)^{(n)}}$ also decays exponentially as
$\mathcal{O}(r^{-n})$.

In order to use the Chebyshev sum for approximating $\psi$, it is inevitable
to compute its inner products. This can be done either by calculating the
inner products or the FFT, at the cost of complicated derivations and
demanding computations. Fortunately, by using the orthogonal interpolation
approach (Mason and Handscomb, 2003) for approximating $z$, computing the
inner products is no longer needed. This is can be done follows. An
interpolant of a function $f$ with polynomial of order $n-1$ is a polynomial
$p_{n}(x)=J_{n-1}f(x)$, such that when it is evaluated at the $n$ distinct
interpolation points $\left\\{x_{i}\right\\}$ it has exactly the same value as
$f(x_{i})$. That is, for any $i=1,\ldots,n$:

$\displaystyle p_{n}(x_{i})=c_{0}+c_{1}x_{i}+\cdots+c_{n}x_{i}^{n-1}$
$\displaystyle=f(x_{i}).$

The resulting polynomial can be written in terms of the Lagrange polynomials
(Mason and Handscomb, 2003, Lemma 6.3) as

$\displaystyle J_{n-1}f(x)$ $\displaystyle=\sum_{i=1}^{n}\ell_{i}(x)f(x_{i}).$

When the interpolation points are selected to be the zeros of the Chebyshev
polynomial of the first kind, $\ell_{i}$ is given by (Mason and Handscomb,
2003, Corollary 6.4 A)

$\displaystyle\ell_{i}(x)$
$\displaystyle=\dfrac{T_{n}(x)}{(n+1)(x-x_{i})U_{n-1}(x_{i})},$ (13)

where $U_{n}$ is the Chebyshev polynomial of the second kind of order $n$.
When approximating $f$ by its $n$-th order Chebyshev interpolant $J_{n-1}f$,
we can approximate the definite integral $\int_{\Omega}fwdx$ by what is known
as Gauss–Chebyshev quadrature (of the first kind) (Mason and Handscomb, 2003,
Theorem 8.4):

$\displaystyle\int_{\Omega}f(\tilde{x})w(\tilde{x})d\tilde{x}$
$\displaystyle\approx Q^{1}_{n}f,$ $\displaystyle Q^{1}_{n}f$
$\displaystyle:=\int_{\Omega}J_{n-1}f(\tilde{x})w(\tilde{x})d\tilde{x}=\sum_{i=1}^{n}\dfrac{\pi}{n}f(\tilde{x}_{i}),$
$\displaystyle w(\tilde{x})$ $\displaystyle=(1-\tilde{x}^{2})^{-1/2},$
$\displaystyle\tilde{x}_{i}$
$\displaystyle=\cos\left(\dfrac{(i-\frac{1}{2})\pi}{n}\right).$ (14)

This approximation is based on the continuous orthogonality of the Chebyshev
polynomials (that is if $i\neq j$, $\langle T_{i},T_{j}\rangle=0$) and is
exact if $f$ is a polynomial of degree $2n-1$ or less (Mason and Handscomb,
2003, Theorem 8.2). It is known that for a function $f\in BV(\bar{\Omega})$,
the integration error satisfies $E_{n}f\leq\pi V_{\bar{\Omega}}(f)/(2n)$ where

$\displaystyle E_{n}f$
$\displaystyle:=\absolutevalue{\int_{\Omega}f(\tilde{x})w(\tilde{x})d\tilde{x}-Q^{1}_{n}f},$
(15)

and $V_{\bar{\Omega}}(f)$ stands for the total variation of $f$ on
$\bar{\Omega}$ (Riess, 1971; Chui, 1972).

To compute the integral of $z(\tilde{x})$ given in (9b), we will apply the
quadrature rule (14). We must be cautious, though, because not every smooth
bijection $\phi$ used to compute the integral of $z$ on $\Omega$ will be
permissible. In fact, for some $\theta>0$, let
$p_{\theta}=\frac{1}{\sqrt{2\pi}}\exp(-(x-\theta)^{2}/2)$. If we choose
$\phi=\sqrt{2}\erf(\tilde{x})^{-1}$, then
$z(\tilde{x})=\frac{1}{2}\exp(\frac{2\phi(\tilde{x})\theta-\theta^{2}}{2})$,
which tends to infinity as $\tilde{x}\to 1$, meaning that $w^{-1}z$ may not be
in $BV(\bar{\Omega})$, resulting in numerical integration failure. We offer
conditions in the following lemma that may be used to verify that the
numerical quadrature (14) produces a consistent integration result, that is,
the error decreases with increasing $n$. Furthermore, because the domain of
the bijection $\phi$ is $\Omega$, we assign $z(1)=\lim_{\tilde{x}\uparrow
1}z(\tilde{x})$, and $z(-1)=\lim_{\tilde{x}\downarrow-1}z(\tilde{x})$.

In the following, we present the main results of this section that give upper-
bound estimates of the approximation errors of the log-partition functions,
and the first and second moments of $c(x)$.

###### Lemma 3.1

Let $\phi:\Omega\rightarrow\mathbb{R}$ be a smooth bijection,
$w(\tilde{x})=(1-\tilde{x}^{2})^{-1/2}$, and Assumption 2.2 be satisfied. Let
also $\zeta:=\phi^{-1}:\mathbb{R}\rightarrow\Omega$. If for any
$\theta\in\Theta$,

$\displaystyle\mathbb{E}_{\theta}\left[\dfrac{1}{\zeta^{\prime}(x)^{2}}\right]$
$\displaystyle<\infty,$ (16a)
$\displaystyle\mathbb{E}_{\theta}\left[\dfrac{\zeta^{\prime\prime}(x)^{2}}{\zeta^{\prime}(x)^{4}}\right]$
$\displaystyle<\infty,$ (16b)

then $r$,$s_{k}$, and $t_{l,k}$, defined by

$\displaystyle r(\tilde{x})$ $\displaystyle:=z(\tilde{x})w(\tilde{x})^{-1},$
$\displaystyle s_{k}(\tilde{x})$
$\displaystyle:=c_{k}(\phi(\tilde{x}))r(\tilde{x}),$ $\displaystyle
t_{l,k}(\tilde{x})$ $\displaystyle:=c_{l}(\phi(\tilde{x}))s_{k}(\tilde{x}),$

belong to $BV(\bar{\Omega})$.

###### Proof 2

First we will prove that $r\in BV(\bar{\Omega})$. To see this, it is enough to
show that $z\in BV(\bar{\Omega})$ since $w^{-1}\in BV(\bar{\Omega})$ and
$BV(\bar{\Omega})$ is an algebra.

Before showing that $z\in BV(\bar{\Omega})$, we claim that for
$\theta\in\Theta$, $\mathbb{E}_{\theta}[\absolutevalue{x^{k}}]<\infty,\forall
k\in\mathbb{N}$ when Assumption 2.2 holds. To see this, we observe that the
moment generating function $\varphi$ with respect to the natural statistics
$\left\\{c_{i}\right\\}$ is given by

$\displaystyle\varphi$
$\displaystyle:=\mathbb{E}_{\theta}\left[\exp(v^{\top}c(x))\right]$
$\displaystyle=\frac{1}{\exp(\psi)}\int_{\mathbb{R}}\exp(c^{\top}(x)\theta)\exp(v^{\top}c(x))dx$
$\displaystyle=\exp(\psi(\theta+v)-\psi(\theta)).$

Therefore,
$\mathbb{E}_{\theta}[c_{i}^{k}(x)]=\left.\partialderivative[k]{\varphi}{v_{i}}\right|_{v=0}$
since $\psi$ is infinitely differentiable on $\Theta$ (Brigo et al., 1995,
Lemma 2.3). Furthermore, since $c_{i}$ is a monomial, we have proven the
claim.

Since $z$ is differentiable on $\Omega$, then the total variation of $z$ is
given by $\int_{\Omega}\absolutevalue{z^{\prime}}d\tilde{x}$. Let
$\beta(x)=\exp(\sum_{i=1}^{m}c_{i}(x)\theta_{i})$. Without loss of generality,
assume that $\phi$ is monotonically increasing, otherwise, choose the
bijection as $-\phi$. Since $\phi$ is monotonically increasing, the total
variation of $z$ is given by

$\displaystyle V_{\bar{\Omega}}(z)$
$\displaystyle=\int_{\Omega}\absolutevalue{z^{\prime}}d\tilde{x}=\int_{\Omega}\absolutevalue{\left(\phi^{\prime}(\tilde{x})\beta(\phi(\tilde{x}))\right)^{\prime}}d\tilde{x}$
$\displaystyle\leq\int_{\Omega}\absolutevalue{\dfrac{\phi^{\prime\prime}(\tilde{x})}{\phi^{\prime}(\tilde{x})}}z(\tilde{x})d\tilde{x}+\int_{\Omega}\sum_{i=1}^{m}\absolutevalue{\phi^{\prime}(\tilde{x})c^{\prime}_{i}(x)\theta_{i}}z(\tilde{x})d\tilde{x}.$

Using $x=\phi(\tilde{x})$ and

$\displaystyle\zeta^{\prime}(x)$
$\displaystyle=\frac{1}{\phi^{\prime}(\tilde{x})},$
$\displaystyle\dfrac{\phi^{\prime\prime}(\tilde{x})}{\phi^{\prime}(\tilde{x})}$
$\displaystyle=-\dfrac{\zeta^{\prime\prime}(x)}{\zeta^{\prime}(x)^{2}},$

and by introducing a positive measure on $\mathbb{R}$,
$d\mu(x)=\beta(x)dx=z(\tilde{x})d\tilde{x}$, we can further simplify

$\displaystyle\int_{\Omega}\absolutevalue{\dfrac{\phi^{\prime\prime}(\tilde{x})}{\phi^{\prime}(\tilde{x})}}z(\tilde{x})d\tilde{x}$
$\displaystyle+\int_{\Omega}\sum_{i=1}^{m}\absolutevalue{\phi^{\prime}(\tilde{x})c^{\prime}_{i}(x)\theta_{i}}z(\tilde{x})d\tilde{x}$
$\displaystyle=\int_{\mathbb{R}}\absolutevalue{\dfrac{\zeta^{\prime\prime}(x)}{\zeta^{\prime}(x)^{2}}}d\mu+\int_{\mathbb{R}}\sum_{i=1}^{m}\absolutevalue{\frac{1}{\zeta^{\prime}(x)}c^{\prime}_{i}(x)\theta_{i}}d\mu.$

By Jensen inequality and (16b), the first part of the righthand side is
finite,

$\displaystyle\int_{\mathbb{R}}\absolutevalue{\dfrac{\zeta^{\prime\prime}(x)}{\zeta^{\prime}(x)^{2}}}d\mu=\exp(\psi)\mathbb{E}_{\theta}\left[\absolutevalue{\dfrac{\zeta^{\prime\prime}(x)}{\zeta^{\prime}(x)^{2}}}\right]<\infty.$

The second part is also finite, since

$\displaystyle\int_{\mathbb{R}}\sum_{i=1}^{m}\absolutevalue{\frac{1}{\zeta^{\prime}(x)}c^{\prime}_{i}(x)\theta_{i}}d\mu$
$\displaystyle=\exp(\psi)\mathbb{E}_{\theta}\left[\sum_{i=1}^{m}\absolutevalue{\frac{1}{\zeta^{\prime}(x)}c^{\prime}_{i}(x)\theta_{i}}\right]$
$\displaystyle\leq\exp(\psi)\sum_{i=1}^{m}\mathbb{E}_{\theta}\absolutevalue{\frac{1}{\zeta^{\prime}(x)}c^{\prime}_{i}(x)\theta_{i}}$
$\displaystyle\leq\exp(\psi)\sum_{i=1}^{m}\left[\mathbb{E}_{\theta}\left[\frac{1}{\zeta^{\prime}(x)^{2}}\right]\right]^{1/2}\left[\mathbb{E}_{\theta}\left[c^{\prime}_{i}(x)^{2}\right]\right]^{1/2}\absolutevalue{\theta_{i}}$

which is less than infinity by (16a), Assumption 2.2, and since any moments of
$x$ are finite. To show that $s_{k}\in BV(\bar{\Omega})$ it is enough to check
that $c_{k}(\phi)z$ is of bounded variation on $\Omega$. As before, we get

$\displaystyle V_{\bar{\Omega}}(c_{k}(\phi)z)$
$\displaystyle=\int_{\Omega}\absolutevalue{\left(c_{k}(\phi(\tilde{x}))z(\tilde{x})\right)^{\prime}}d\tilde{x}$
$\displaystyle\leq\int_{\mathbb{R}}\absolutevalue{\dfrac{c_{k}^{\prime}(x)}{\zeta^{\prime}(x)}}d\mu+\int_{\mathbb{R}}\absolutevalue{\dfrac{c_{k}(x)\zeta^{\prime\prime}(x)}{\zeta^{\prime}(x)^{2}}}d\mu+\int_{\mathbb{R}}\sum_{i=1}^{m}\absolutevalue{\dfrac{c_{k}(x)c_{i}^{\prime}(x)}{\zeta^{\prime}(x)}\theta_{i}}d\mu.$

The first and last terms of the right hand side of the last equation are
finite because of condition (16a) and since any moment of $x$ is finite. The
second part follows from Hölder’s inequality and from condition (16b).

By Assumption 2.2, there exists $q\in\mathbb{N}$ such that
$c_{l}(x)c_{k}(x)=x^{q}$. Therefore, the result for $t_{l,k}\in
BV(\bar{\Omega})$ follows similarly as $s_{k}\in BV(\bar{\Omega})$, which
completes the proof.

Lemma 3.1 shows that it is sufficient to have the conditions (16a) and (16b)
hold in order to have $r,s_{k},t_{l,k}\in BV(\bar{\Omega})$. Consider two
bijections $\phi=\tan(\frac{\pi}{2}\tilde{x})$ and
$\phi=\frac{\tilde{x}}{1-\tilde{x}^{2}}$. Since their corresponding
$\absolutevalue{\frac{1}{\zeta^{\prime}(x)^{2}}}$ and
$\absolutevalue{\frac{\zeta^{\prime\prime}(x)}{\zeta^{\prime}(x)^{2}}}$ are
less than absolute values of some polynomials, the conditions (16a) and (16b)
hold for any density from the exponential families with monomial natural
statistics. On the other hand, the bijection $\phi=\tanh^{-1}(\tilde{x})$
under $p_{\theta}=\frac{1}{\sqrt{2\pi}}\exp(-(x-\theta)^{2}/2)$ satisfies
these conditions for $\theta\in(-\pi,\pi)$. For related discussion on the
effect of bijection in the context of sparse-grid integration, see Griebel and
Oettershagen (2014).

As a direct consequence of Lemma 3.1 and results from Chui (1972), we have the
following upper bounds on the integration errors of $r,s_{i}$ and $t_{l,k}$.

###### Corollary 3.1

The integration errors of $r,s_{i}$ and $t_{l,k}$ for $i,l,k=1,\ldots,m$
satisfy

$\displaystyle E_{n}r$ $\displaystyle\leq\pi V_{\bar{\Omega}}(r)/(2n),$ (17a)
$\displaystyle E_{n}s_{i}$ $\displaystyle\leq\pi
V_{\bar{\Omega}}(s_{i})/(2n),$ (17b) $\displaystyle E_{n}t_{l,k}$
$\displaystyle\leq\pi V_{\bar{\Omega}}(t_{l,k})/(2n),$ (17c)

respectively.

We now present the main result for error bounds of $\psi(\theta)^{(n)}$,
$\partial_{i}\psi(\theta)^{(n)}$, and
$\partial_{i}\partial_{j}\psi(\theta)^{(n)}$.

###### Theorem 3.1

Let Assumption 2.2 be satisfied. Let $\phi$ be a smooth bijection satisfying
the conditions in Lemma 3.1, and
$r,\left\\{s_{k}\right\\},\left\\{t_{l,k}\right\\}$ be functions defined as in
Lemma 3.1 whose integrations errors are given by (17). Denote by
$\psi(\theta)^{(n)}$ the approximation of log-partition function $\psi$ where
$\int_{\Omega}zd\tilde{x}$ in (9a) is replaced by $Q_{1}^{n}r$. The log-
partition function approximation error satisfies,

$\displaystyle\absolutevalue{\psi-\psi^{(n)}}$
$\displaystyle\leq\dfrac{1}{Q^{1}_{n}r}E_{n}r,$

where $Q^{1}_{n}r>0$. Furthermore, the approximation errors of the first and
the second derivatives of $\psi$ have the following upper bounds:

$\displaystyle\absolutevalue{\partial_{i}\psi-\partial_{i}\psi^{(n)}}$
$\displaystyle\leq\dfrac{1}{\exp(\psi)}E_{n}s_{i}+\absolutevalue{\dfrac{Q^{1}_{n}s_{i}}{\exp(\psi)Q^{1}_{n}r}}E_{n}r,$
$\displaystyle\absolutevalue{\partial_{j}\partial_{i}\psi(\theta)-\partial_{j}\partial_{i}\psi(\theta)^{(n)}}$
$\displaystyle\leq\dfrac{1}{\exp(\psi)}E_{n}t_{j,i}+\absolutevalue{\dfrac{Q^{1}_{n}t_{j,i}}{\exp(\psi)Q^{1}_{n}r}}E_{n}r$
$\displaystyle+\dfrac{1}{\exp(\psi)}E_{n}s_{j}+\absolutevalue{\dfrac{Q^{1}_{n}s_{j}}{\exp(\psi)Q^{1}_{n}r}}E_{n}r.$

###### Proof 3

Without loss of generality, we can again assume that $\phi$ is monotonically
increasing (otherwise choose $-\phi$). Observe that on $\Omega$,
$r(\tilde{x})>0$ since the bijection $\phi$ is monotonically increasing.
Therefore, by (14), we have $Q^{1}_{n}r>0$. Furthermore, as in Proposition
3.1, we have

$\displaystyle\absolutevalue{\psi(\theta)-\psi(\theta)^{(n)}}$
$\displaystyle=\absolutevalue{\log\left[\dfrac{\int_{\Omega}z(\tilde{x})d\tilde{x}}{Q^{1}_{n}r}\right]}=\absolutevalue{\log\left[\dfrac{\int_{\Omega}z(\tilde{x})d\tilde{x}-Q^{1}_{n}r}{Q^{1}_{n}r}+1\right]}$
$\displaystyle\leq\log\left[\dfrac{E_{n}r}{Q^{1}_{n}r}+1\right]\leq\dfrac{E_{n}r}{Q^{1}_{n}r}.$

For the error $\absolutevalue{\partial_{i}\psi-\partial_{i}\psi^{(n)}}$, we
have

$\displaystyle\absolutevalue{\partial_{i}\psi(\theta)-\partial_{i}\psi(\theta)^{(n)}}$
$\displaystyle=\absolutevalue{\partial_{i}\left[\log\int_{\Omega}z(\tilde{x})d\tilde{x}-\log
Q^{1}_{n}r\right]}$
$\displaystyle=\absolutevalue{\dfrac{\int_{\Omega}c_{i}(\phi(\tilde{x}))z(\tilde{x})d\tilde{x}}{\int_{\Omega}z(\tilde{x})d\tilde{x}}-\dfrac{Q^{1}_{n}s_{i}}{Q^{1}_{n}r}}$
$\displaystyle\leq\absolutevalue{\dfrac{E_{n}s_{i}}{\int_{\Omega}z(\tilde{x})d\tilde{x}}+\dfrac{Q^{1}_{n}s_{i}}{\int_{\Omega}z(\tilde{x})d\tilde{x}}-\dfrac{Q^{1}_{n}s_{i}}{Q^{1}_{n}r}}$
$\displaystyle\leq\absolutevalue{\dfrac{E_{n}s_{i}}{\int_{\Omega}zd\tilde{x}}+Q^{1}_{n}s_{i}\left[\dfrac{Q^{1}_{n}r-\int_{\Omega}z(\tilde{x})d\tilde{x}}{Q^{1}_{n}r\int_{\Omega}zd\tilde{x}}\right]}$
$\displaystyle\leq\dfrac{1}{\exp(\psi)}E_{n}s_{i}+\absolutevalue{\dfrac{Q^{1}_{n}s_{i}}{\exp(\psi)Q^{1}_{n}r}}E_{n}r.$

Due to the fact that

$\displaystyle\partial_{j}\partial_{i}\psi$
$\displaystyle=\dfrac{\int_{\Omega}c_{i}(\phi(\tilde{x}))c_{j}(\phi(\tilde{x}))z(\tilde{x})d\tilde{x}}{\exp(\psi)}-\dfrac{\int_{\Omega}c_{j}(\phi(\tilde{x}))z(\tilde{x})d\tilde{x}}{\exp(\psi)^{2}},$

using a similar technique as for the second part above, we get

$\displaystyle\absolutevalue{\partial_{j}\partial_{i}\psi(\theta)-\partial_{j}\partial_{i}\psi(\theta)^{(n)}}$
$\displaystyle\leq\dfrac{1}{\exp(\psi)}E_{n}t_{j,i}+\absolutevalue{\dfrac{Q^{1}_{n}t_{j,i}}{\exp(\psi)Q^{1}_{n}r}}E_{n}r$
$\displaystyle+\dfrac{1}{\exp(\psi)}E_{n}s_{j}+\absolutevalue{\dfrac{Q^{1}_{n}s_{j}}{\exp(\psi)Q^{1}_{n}r}}E_{n}r.$

In essence, Theorem 3.1 guarantees that the approximation errors of the log-
partition function as well as the first and the second derivatives of $\psi$
have upper bounds that are linear with respect to the approximation error of
$\int_{\Omega}z(\tilde{x})d\tilde{x}$. This is reasonable since the partial
derivative with respect to the natural parameter $\theta_{i}$ is a linear
operator.

###### Remark 3.2

It is possible to have the error bounds above satisfy even stronger decaying
inequalities if $r,s_{k},t_{l,k}$ are differentiable or twice differentiable
on $\bar{\Omega}$, at the cost of imposing stronger conditions on the
bijection $\phi$ (Chui, 1972).

##  4 Calculating log-partition function in multidimensional case

In this section, we will introduce a sparse-grid integration method to
calculate the log-partition function for multidimensional problems. For
detailed exposition of sparse-grid integration methods, we refer the reader to
Gerstner and Griebel (1998).

As we mentioned earlier in Section 3, to compute the log-partition function,
we can use polynomial interpolation of order $n$, which will give exact
integration results if the integrand is a polynomial of degree $2n-1$. One of
the main issues that arises in the Gauss–Chebyshev quadrature is that the
quadrature nodes of the polynomial of order $n-1$, $\Gamma_{n-1}^{1}$, are not
a subset of the quadrature nodes of the polynomials of order $n$,
$\Gamma_{n}^{1}$, that is, the nodes of Gauss–Chebyshev quadrature are not
nested. This means that if we need to increase the accuracy of the
integration, we cannot reuse $\Gamma_{n-1}^{1}$ and add the remaining missing
nodes, but instead, the whole set of nodes needs to be computed from scratch.
In multidimensional quadratures, having a quadrature rule with nested
integration nodes is critical for computational efficiency.

Gauss–Chebyshev quadrature relies on the continuous orthogonality of the
Chebyshev polynomials. There is also a discrete orthogonality property of the
Chebyshev polynomials (Mason and Handscomb, 2003, Section 4.6). The quadrature
rule that takes benefit of this property is known as the Clenshaw–Curtis
quadrature. This quadrature rule has a polynomial degree of exactness of $n-1$
(Gerstner and Griebel, 1998) and its nodes are nested. If the inner product
weight is set to identity, the quadrature nodes become the zeros of the
Legendre polynomials and the weights can be calculated by integrating the
associated polynomials. This Gauss–Legendre quadrature is not nested. Kronrod
(1966) proposed an extension of $n$-point Gauss quadrature by other $n+1$
points so that the quadrature nodes are nested. This quadrature has the
polynomial degree of exactness equal to $3n+1$ (Laurie, 1997; Gerstner and
Griebel, 1998). Solving the quadrature nodes for this scheme is, however,
difficult, and Patterson (1968) proposed an iterated version of the Kronrod
scheme.

For the case $d_{x}>1$, we can use extensions of aforementioned quadratures to
$d$-dimensional hypercube integration (Novak and Ritter, 1996), see also
Barthelmann et al. (2000); Bungartz and Griebel (2004); Gerstner and Griebel
(1998); Judd and Skrainka (2011). This is done via what is known as Smolyak’s
construction (Smolyak, 1963). In Smolyak’s approach, it is sufficient to have
the means for solving the integration problem for unidimensional problem
efficiently, since the algorithm for multidimensional case is fully determined
in terms of the algorithm of the unidimensional one (Wasilkowski and
Woźniakowski, 1995). Consider a one-dimensional quadrature formula, where the
univariate integrand $f$ is interpolated with an $n$ degree polynomial

$Q_{n}^{1}f:=\sum_{i}^{N_{n}^{1}}w_{i,n}f(x_{i,n}).$

The differences of the quadrature formulas of different degrees are given by

$\Delta_{n}^{1}f:=(Q_{n}^{1}-Q_{n-1}^{1})f,$

where we define $Q_{0}^{1}f=0$. The differences are quadrature formulas
defined on the union of the quadrature nodes
$\Gamma_{n}^{1}\bigcup\Gamma_{n-1}^{1}$, where for the case of the nested
nodes, $\Gamma_{n}^{1}=\Gamma_{n}^{1}\bigcup\Gamma_{n-1}^{1}$. The Smolyak’s
construction for $d$-dimensional functions $f$ for $n\in\mathbb{N}$ and
$\mathbf{k}\in\mathbb{N}^{d}$ is given by (Gerstner and Griebel, 1998; Delvos,
1982)

$\displaystyle Q^{d}_{n}f$ $\displaystyle:=\sum_{\norm{\mathbf{k}}_{1}\leq
n+d-1}\left(\Delta_{k_{1}}^{1}\otimes\cdots\otimes\Delta_{k_{d}}^{1}\right)f.$

Let us denote the number of quadrature points in the unidimensional polynomial
interpolation of order $n$ as $N_{n}^{1}$. An important feature of Smolyak’s
construction is that instead of having ${(N_{n}^{1})}^{d}$ quadrature points
as the case with the tensor product of unidimensional grids, it has far fewer
quadrature points to be evaluated. If $N_{n}^{1}=\mathcal{O}(2^{n})$, under
Smolyak’s construction, the order of $N_{n}^{d}$ is
$\mathcal{O}(2^{n}n^{d-1})$, which is substantially lower than
$\mathcal{O}(2^{nd})$ resulting from the product rules (Novak and Ritter,
1996; Gerstner and Griebel, 1998).

As for the error bound, for a function $f\in\mathcal{C}^{r}$, we have
$E_{n}^{1}f=\mathcal{O}((N_{n}^{1})^{-r})$, where
$E_{n}^{d}f:=\absolutevalue{\int_{\Omega}fdx-Q^{d}_{n}f}$. This bound holds
for all interpolation-based quadrature formulas with positive weights, such as
the Clenshaw–Curtis, Gauss–Patterson, and Gauss–Legendre schemes. When we use
these quadrature formulas as the unidimensional basis, if
$f\in\mathcal{W}^{r}_{d}$, where

$\mathcal{W}^{r}_{d}:=\left\\{f:\Omega\rightarrow\mathbb{R}:\norm{\dfrac{\partial^{\norm{\mathbf{s}}_{1}}f}{\partial
x_{1}^{s_{1}}\ldots\partial
x_{d}^{s_{d}}}}_{\infty}<\infty,\norm{\mathbf{s}}_{1}\leq r\right\\},$

then we have $E_{n}^{d}f=\mathcal{O}(2^{-nr}n^{(d-1)(r+1)})$ (Wasilkowski and
Woźniakowski, 1995; Gerstner and Griebel, 1998). We may establish the error
estimate for the multidimensional case using a similar approach used in the
proof of Proposition 5.2, but with caution, because the equivalent of $z$ (9b)
for the multidimensional case may likewise have an unbounded total variation.

Specifically, for the integration problem (9a), the density $p_{\theta}(x)$ at
the points on the boundary of the canonical hypercube is almost zero since
$\exp(\sum c_{i}(x)\theta_{i})$ in (2) goes to zero as $\norm{x}$ goes to
infinity. This implies that the grid points that lie on the boundary of the
canonical cube can be neglected in the quadrature calculation. This also
reduces the required number of grid points.

##  5 Extending the projection filter to $d_{x}>1$

The numerical implementation of the projection filter for the exponential
family relies on computing the expected values of the natural statistics and
the Fisher matrix. In the previous sections, we have shown how to calculate
the log-partition function for both the unidimensional and multidimensional
cases. The repetitive uses of numerical integrations for computing the Fisher
information matrix and expectations of natural statistics can be avoided by
using automatic differentiation and relations (8a). However, it still remains
to show how to obtain the expectations of the remaining statistics in
$\tilde{c}$.

In unidimensional dynamics, it is well known that the Fisher information
matrix and expectations of the natural statistics can be obtained efficiently
via a recursion, given that we have calculated the first $m$ of them (Brigo et
al., 1995). Unfortunately, this recursion does not hold for multidimensional
problems. This can be concluded from Proposition 5.1.

Therefore, in order to solve the projection filtering problems efficiently
when $d_{x}\geq 2$, we show that it is possible to obtain the expectations of
the expanded natural statistics $\tilde{c}$ via computing certain partial
derivatives. This procedure is particularly efficient by using modern
automatic differentiation tools. The main result is given in Proposition 5.2.

Let us set some notations as follows. Let
$c:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{m}$, where $m$ is the dimension of
the parameter space $\Theta$. Define a multi-index
$\underline{i\mkern-2.0mu}\mkern 2.0mu\in\mathbb{N}^{d_{x}}$ and

$x^{\underline{i\mkern-2.0mu}\mkern
2.0mu}:=x_{1}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(1)}x_{2}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(2)}\cdots
x_{d_{x}}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(d_{x})}.$ (18)

We can identify the elements of $c$ by the multi-index
$\underline{i\mkern-2.0mu}\mkern 2.0mu$. As before, let
$\eta_{\underline{i\mkern-2.0mu}\mkern
2.0mu}:=\mathbb{E}_{\theta}[c_{\underline{i\mkern-2.0mu}\mkern 2.0mu}]$. For
two multi-indexes $\underline{i\mkern-2.0mu}\mkern
2.0mu,\underline{j\mkern-2.0mu}\mkern 2.0mu\in\mathbb{N}^{d_{x}}$, we define
an addition operation element-wise. For simplicity, we also write
$\underline{j\mkern-2.0mu}\mkern 2.0mu=\underline{i\mkern-2.0mu}\mkern
2.0mu+a$ where $\underline{j\mkern-2.0mu}\mkern
2.0mu(k)=\underline{i\mkern-2.0mu}\mkern 2.0mu(k)+a$, for $k=1,\ldots,d_{x}$
and $a\in\mathbb{N}$. Then we have the following result as an analog of Lemma
3.3 from Brigo et al. (1995).

###### Proposition 5.1

Let a set of exponential parametric densities be given by

$\displaystyle S$
$\displaystyle=\left\\{p(\cdot,\theta),\theta\in\Theta\right\\},$
$\displaystyle p(x,\theta)=\exp(c(x)^{\top}\theta-\psi),$

where $\Theta\subset\mathbb{R}^{m}$ is open. The log-partition function $\psi$
is infinitely differentiable on $\Theta$. Furthermore, the first, second, and
$k$-th moments of the natural expectations are given by

$\displaystyle\mathbb{E}_{\theta}\matrixquantity[c_{\underline{i\mkern-2.0mu}\mkern
2.0mu}]$ $\displaystyle=\partial_{\underline{i\mkern-2.0mu}\mkern
2.0mu}\psi(\theta)=:\eta_{\underline{i\mkern-2.0mu}\mkern 2.0mu}(\theta),$
(19a)
$\displaystyle\mathbb{E}_{\theta}\matrixquantity[c_{\underline{i\mkern-2.0mu}\mkern
2.0mu}c_{\underline{j\mkern-2.0mu}\mkern 2.0mu}]$
$\displaystyle=\partial^{2}_{\underline{i\mkern-2.0mu}\mkern
2.0mu,\underline{j\mkern-2.0mu}\mkern
2.0mu}\psi(\theta)+\eta_{\underline{i\mkern-2.0mu}\mkern
2.0mu}(\theta)\eta_{\underline{j\mkern-2.0mu}\mkern 2.0mu}(\theta),$ (19b)
$\displaystyle\mathbb{E}_{\theta}\matrixquantity[c_{\underline{i\mkern-2.0mu}\mkern
2.0mu_{1}}\cdots c_{\underline{i\mkern-2.0mu}\mkern 2.0mu_{k}}]$
$\displaystyle=\exp(-\psi(\theta))\partial^{k}_{\underline{i\mkern-2.0mu}\mkern
2.0mu_{1}\cdots\underline{i\mkern-2.0mu}\mkern 2.0mu_{k}}\exp(\psi(\theta)).$
(19c)

The Fisher information matrix satisfies

$\displaystyle g_{\underline{i\mkern-2.0mu}\mkern
2.0mu\underline{j\mkern-2.0mu}\mkern
2.0mu}=\partial^{2}_{\underline{i\mkern-2.0mu}\mkern
2.0mu\underline{j\mkern-2.0mu}\mkern 2.0mu}\psi(\theta).$ (20)

In particular, if the natural statistics are a collection of monomials of $x$,

$\displaystyle c=\matrixquantity[x^{\underline{i\mkern-2.0mu}\mkern
2.0mu_{1}},\ldots,x^{\underline{i\mkern-2.0mu}\mkern 2.0mu_{m}}],$

the Fisher information matrix satisfies

$\displaystyle g_{\underline{i\mkern-2.0mu}\mkern
2.0mu\underline{j\mkern-2.0mu}\mkern
2.0mu}=\eta_{\underline{i\mkern-2.0mu}\mkern
2.0mu+\underline{j\mkern-2.0mu}\mkern
2.0mu}(\theta)-\eta_{\underline{i\mkern-2.0mu}\mkern
2.0mu}(\theta)\eta_{\underline{j\mkern-2.0mu}\mkern 2.0mu}(\theta),$ (21)

where $\eta_{\underline{i\mkern-2.0mu}\mkern
2.0mu+\underline{j\mkern-2.0mu}\mkern
2.0mu}(\theta):=\mathbb{E}_{\theta}[x^{\underline{i\mkern-2.0mu}\mkern
2.0mu+\underline{j\mkern-2.0mu}\mkern 2.0mu}]$ and for any
$\underline{i\mkern-2.0mu}\mkern 2.0mu\in\mathbb{N}^{d_{x}}$

$\displaystyle\mathbb{E}_{\theta}\matrixquantity[x^{\underline{i\mkern-2.0mu}\mkern
2.0mu}]$
$\displaystyle=(-1)^{d_{x}}\int_{\mathbb{R}^{d_{x}}}\dfrac{x^{\underline{i\mkern-2.0mu}\mkern
2.0mu+1}}{\pi(\underline{i\mkern-2.0mu}\mkern
2.0mu+1)}\dfrac{\partial^{d_{x}}p(x,\theta)}{\partial x_{1}\cdots\partial
x_{d_{x}}}dx,$ (22)

where

$\displaystyle\pi(\underline{i\mkern-2.0mu}\mkern 2.0mu)$
$\displaystyle=\underline{i\mkern-2.0mu}\mkern
2.0mu(1)\times\cdots\times\underline{i\mkern-2.0mu}\mkern 2.0mu({d_{x}}).$

###### Proof 4

All results except (21) and (22) follow directly from Brigo et al. (1995,
Lemma 2.3) and Amari (1985) using multi-index notations. By definition,
$\mathbb{E}_{\theta}[c_{\underline{i\mkern-2.0mu}\mkern
2.0mu}c_{\underline{j\mkern-2.0mu}\mkern
2.0mu}]=\eta_{\underline{i\mkern-2.0mu}\mkern
2.0mu+\underline{j\mkern-2.0mu}\mkern 2.0mu}(\theta)$, hence using (19b),
(19a), and (20), we obtain (21). For (22), observe that

$\displaystyle\mathbb{E}_{\theta}\matrixquantity[x^{\underline{i\mkern-2.0mu}\mkern
2.0mu}]$
$\displaystyle=\idotsint_{\mathbb{R}^{d_{x}}}x_{1}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(1)}x_{2}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(2)}\cdots
x_{d}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(d_{x})}p(x,\theta)dx_{1}\cdots
dx_{d_{x}}$
$\displaystyle=\idotsint_{\mathbb{R}^{d_{x}-1}}x_{1}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(1)}x_{2}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(2)}\cdots
x_{d_{x}-1}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(d_{x}-1)}\left[\int_{-\infty}^{\infty}x_{d_{x}}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(d_{x})}p(x,\theta)dx_{d_{x}}\right]dx_{1}\cdots dx_{d_{x}-1}$
$\displaystyle=\idotsint_{\mathbb{R}^{d_{x}-1}}x_{1}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(1)}x_{2}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(2)}\cdots
x_{d_{x}-1}^{\underline{i\mkern-2.0mu}\mkern 2.0mu(d_{x}-1)}$
$\displaystyle\times\left[\left.\dfrac{x_{d_{x}}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(d_{x})+1}}{\underline{i\mkern-2.0mu}\mkern
2.0mu(d_{x})+1}p(x,\theta)\right|^{\infty}_{-\infty}-\int_{-\infty}^{\infty}\dfrac{x_{d_{x}}^{\underline{i\mkern-2.0mu}\mkern
2.0mu(d_{x})+1}}{\underline{i\mkern-2.0mu}\mkern
2.0mu(d_{x})+1}\partialderivative{p(x,\theta)}{x_{d_{x}}}dx_{d_{x}}\right]dx_{1}\cdots
dx_{d_{x}-1}.$

Doing a similar integration for the remaining integration variables
$x_{1},\ldots,x_{d_{x}-1}$ and since $p(x,\theta)$ vanishes as
$\norm{x}\rightarrow\infty$ faster than any polynomial, we end up with (22).

In unidimensional problems (22) admits a recursion, where we can express
$\eta_{m+i}$ as a linear combination of $\eta_{i},\ldots,\eta_{m+i-1}$ (Brigo
et al., 1995). However, the recursion does not hold for $d_{x}>1$. This can be
verified by evaluating (5.1) for the case $d_{x}=2$ and natural statistics
$\left\\{x^{\underline{i\mkern-2.0mu}\mkern 2.0mu}\right\\}$, where
$\underline{i\mkern-2.0mu}\mkern
2.0mu\in\left\\{(0,1),(0,2),(1,0),(1,1),(2,0)\right\\}$. Instead of relying on
a recursion to compute the higher moments, we can use a simple procedure which
is given in the following proposition.

###### Proposition 5.2

Let $s(x):\mathbb{R}^{d}\rightarrow\mathbb{R}$ be another statistic, linearly
independent of the natural statistics $c(x)$ with the corresponding natural
parameters $\theta$ (see Proposition 5.1). 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 the expectation of $s(x)$ under
$p_{\theta}$ is given by

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

###### Proof 5

Define a density as
$\bar{p}(x,\tilde{\theta}):=\exp(\tilde{c}^{\top}\tilde{\theta}-\tilde{\psi}(\tilde{\theta}))$,
where $\bar{p}(x,\theta_{\ast})=p(x,\theta)$. Then, by Proposition 5.1:

$\displaystyle\mathbb{E}_{\theta}[s(x)]$
$\displaystyle=\int_{\mathcal{X}}p(x,\theta)s(x)dx=\int_{\mathcal{X}}\bar{p}(x,\theta_{\ast})s(x)dx=\left.\partialderivative{\tilde{\psi}(\tilde{\theta})}{\tilde{\theta}_{m+1}}\right|_{\tilde{\theta}=\theta_{\ast}}.$

The above proposition shows that by using an expanded exponential family, the
higher moments can be computed via partial differentiation. This removes the
difficulty of obtaining the higher moments which previously could only be
cheaply obtained for unidimensional problems via recursion (Brigo et al.,
1995). The linear independence condition in this assumption is required since
otherwise, the expectation of the statistic $s(x)$ can be obtained by a linear
combination of the expectations of the elements of $c(x)$.

##  6 Numerical experiments

(a) density evolution $\phi=\text{arctanh}$ (c) Hellinger distance
$\phi=\text{arctanh}$ (d) Hellinger distance
$\phi=\frac{\tilde{x}}{1-\tilde{x}^{2}}$ Figure 1: Figures 1a and 1a show the
evolution of densities solved by approximating the Kushner–Stratonovich
equation via a finite difference scheme (KS-FDiff), and the projection filter
using different bijections and numbers of Chebyshev nodes (blue, green and red
for 12, 48, and 96 quadrature nodes respectively), while Figures 1c and 1d
show their respective Hellinger distances. It can be seen that the Hellinger
distances of the projection filter and the solution of the finite difference
scheme using 48 and 96 quadratures are practically indistinguishable.

In this section, we demonstrate the application of the proposed implementation
of the projection filter to selected optimal filtering problems. We compare
the projection filter densities against the densities obtained by solving the
Kushner–Stratonovich equation via the central finite difference approach (Bain
and Crisan, 2009) and kernel density estimate from a bootstrap particle filter
with systematic resampling (Chopin, 2020). All numerical simulations in this
section are implemented using JAX (Bradbury et al., 2018) where the automatic
differentiation is supported by default. The sparse-grid implementation is
taken from the TASMANIAN library (Stoyanov, 2015).

For the unidimensional problem below, we use a simple central finite-
difference scheme to approximate the solution of the Kushner–Stratonovich
equation. This numerical scheme is stable enough for the simulation period
that we are interested in. For the nonlinear two-dimensional problem, we use
the Crank–Nicolson scheme since for this problem the central finite-difference
scheme is unstable in the grid space and simulation time that we are
interested in. We set the ratio $\Delta t/(\Delta x^{2})$ to be significantly
below $0.5$ according to the von Neumann stability analysis (Smith, 1985).

###  6.1 Unidimensional example

The first example is the following scalar system dynamic with a nonlinear
measurement model:

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

with independent standard Brownian motions $\\{W_{t},t\geq 0\\}$ and
$\\{V_{t},t\geq 0\\}$, and where $\sigma=0.4,\beta=0.8$ are constants. We
generate one measurement trajectory realization from the model with $x_{0}=1$.
The simulation time step is set to be $10^{-4}$.

We use the exponential manifold with $c_{i}\in\\{x,x^{2},x^{3},x^{4}\\}$ and
choose the initial condition to be $\theta_{0}=[0,1,0,-1]$. We assume that
this initial condition is exactly the initial density of the dynamical system
which corresponds to a double-mode non-Gaussian density with peaks at $-1$ and
$1$. We compare two choices of the bijections, the first is
$\text{arctanh}(\tilde{x})$ and the second is
$\frac{\tilde{x}}{1-\tilde{x}^{2}}$. Furthermore, in the Chebyshev quadrature
(14) we compare three choices for the number of Chebyshev nodes, where we set
them to be $12$, $48$, and $96$. We also solve the Kushner–Stratonovich SPDE
with finite differences on an equidistant grid $(-5,5)$ with one thousand
points.

The evolution of the densities provided by the Kushner–Stratonovich equation
approach and projection filter using the different settings can be seen in
Figures 1a and LABEL:fig:boyd_wireframe. The Hellinger distances between the
projection filter densities and the numerical solutions of the
Kushner–Stratonovich equation can be seen in Figures 1c and 1d. These figures
show that for this particular filtering problem when using a proper selection
of a bijection function, the projection filter solution using a few nodes is
very close to the one using a very high number of nodes. In particular, when
using arctanh as the bijection, with the node number as low as $12$, the
Hellinger distance of solution of the projection filter to the
Kushner–Stratonovich solution until $t=1.4$ is below $10^{-3}$. This is not
the case when we use $\frac{\tilde{x}}{1-\tilde{x}^{2}}$, as its Hellinger
distance is significantly larger with lower than with higher number of nodes.
The interpolation using the latter bijection with $12$ nodes is less accurate
since the initial condition as can be seen in the Figure
LABEL:fig:boyd_wireframe. This inaccuracy leads to a false double-peaked
density that appears from time $1.0$

We can also see that the projection filter solution for the case of $48$ nodes
is almost indistinguishable from that of $96$, regardless of the choice of the
bijection, with the Hellinger distance to the Kushner–Stratonovich solution
being kept below $10^{-4}$ for most of the time until $1.4$. This shows one
numerical benefit of this approach for this particular example: the filter
equation can be propagated efficiently using a very low number of integration
nodes.

###  6.2 Two-dimensional example

Here, we show that the projection filter solved using the proposed method can
perform well in many interesting problems. In these experiments we use the
Gauss–Patterson scheme with the canonical domain $(-1,1)^{\otimes d_{x}}$.

####  6.2.1 Comparison to Kalman–Bucy filter

In this experiment, we consider the following stochastic estimation problem on
a two-dimensional linear system:

$\displaystyle d\matrixquantity[x_{1,t}\\\ x_{2,t}]$
$\displaystyle=-\matrixquantity[x_{1,t}\\\
x_{2,t}]dt+\sigma_{w}\matrixquantity[dW_{1,t}\\\ dW_{2,t}],$ (24a)
$\displaystyle d\matrixquantity[y_{1,t}\\\ y_{2,t}]$
$\displaystyle=-\matrixquantity[x_{1,t}\\\
x_{2,t}]dt+\sigma_{v}\matrixquantity[dV_{1,t}\\\ dV_{2,t}].$ (24b)

We set $dt=10^{-3}$, $\sigma_{v}=10^{-1}$, and $\sigma_{w}=1$a, and simulate
the SDE with $n_{t}=1000$ Euler–Maruyama steps and record the measurement
trajectory. To obtain a reference result, we implement the Kalman–Bucy filter
using the Euler scheme. For the projection filter, we compare the two
integration procedures: the sparse-grid methods and the quasi Monte Carlo
(qMC) (Leobacher and Pillichshammer, 2014). We choose the arctanh function as
the bijection function based on our previous examination on the unidimensional
problem. We use the Gauss–Patterson scheme for unidimensional sparse grid
where we compare various integration levels. For the qMC, we choose the Halton
low discrepancy sequence. The number of quadrature points for qMC is chosen to
be equal to the selected sparse grid quadrature points. We use the Gaussian
family in our projection filter. This means that we set the natural statistics
to be $\left\\{x^{\underline{i\mkern-2.0mu}\mkern 2.0mu}\right\\}$, where
$\underline{i\mkern-2.0mu}\mkern
2.0mu\in\left\\{(0,1),(0,2),(1,0),(1,1),(2,0)\right\\}$.

Figures 6.2.1 and 6.2.1 illustrate the Hellinger distances between the
approximated solution of the projection filter and the solution of the
Kalman–Bucy filter, whereas Figures 6.2.1 and 6.2.1 show the estimation
outcomes for both states. Using 49 quadrature points, which is comparable to
level 3 in sparse-grid integration, both integration methods provide
estimation results that are quite similar. The sparse-grid based projection
filter approximates the standard Kalman–Bucy filter solution more accurately
than qMC. As more quadrature nodes are added, the Hellinger distances rapidly
decrease, as can be seen in Figures 6.2.1 and 6.2.1. With the sparse-grid
level set to 5, the sparse-grid based projection filter generates a
practically negligible Hellinger distance to the Kalman-Bucy solution after
$t=0.4$. In contrast, the Hellinger distances of the densities determined from
qMC integration continue to oscillate around $10^{-4}$.

(a) SPG-level 3 (c) SPG-level 4 (d) qMC-129 Figure 2: The Hellinger distances
between the densities obtained by solving the Kalman–Bucy equation and those
obtained by projection filter, solved with different integrators. The number
of points in the quasi Monte Carlo integrators are set to be equal to the
respective number of sparse grid points.  (a) SPG-level 5 (c) SPG-level 6 (d)
qMC-769 Figure 3: Continuation of Figure 6.2.1 for different levels of sparse
grid integrators and quasi Monte Carlo points. (a) SPG-level 3 (c) SPG-level 6
(d) qMC-769 Figure 4: Comparison of mean and one standard deviation confidence
intervals for $x_{1,t}$ obtained with the Kalman–Bucy filter and the
projection filter solved via different integration setups.  (a) SPG-level 3
(c) SPG-level 6 (d) qMC-769 Figure 5: Comparison of mean and one standard
deviation confidence intervals for $x_{2,t}$ obtained with the Kalman–Bucy
filter and the projection filter solved via different integration setups.

####  6.2.2 Van der Pol oscillator

In this section, we compare the projection filter, the numerical solution of
the Kushner–Stratonovich equation obtained via a Crank–Nicolson scheme, and a
bootstrap particle filter with systematic resampling (Chopin, 2020). The
filtering model considered is the partially observed the Van der Pol
oscillator:

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

For simulation, we set $\mu=0.3$ and $\sigma_{v}=\sigma_{w}=1$. We also set
$dt=2.5\times 10^{-4}$. We discretize the dynamic model (25) using the
Euler–Maruyama scheme for both the particle filter and the simulation. For the
particle filter, we use $4\times 10^{4}$ samples at each time. The time step
of the finite-difference scheme is set to be $dt=1\times 10^{-3}$. The two-
dimensional grid is set to be $(-6,6)^{2}$ with $150^{2}$ equidistant points.
In the projection filter we use a sparse-grid integration scheme where we set
the level to 8, which corresponds to 4097 quadrature points, and we set the
bijection to arctanh. The natural statistics in the projection filter are set
to be $\left\\{x^{\underline{i\mkern-2.0mu}\mkern 2.0mu}\right\\}$, where
$\underline{i\mkern-2.0mu}\mkern
2.0mu\in\\{(0,1),(0,2),(0,3),(0,4),(1,0),(1,1),$
$(1,2),(1,3),(2,0),(2,1),(2,2),(3,0),(3,1)$ $,(4,0)\\}$. Further, the initial
density is set to be a Gaussian density, with variance equal to $5\times
10^{1}\sigma_{w}^{2}dtI$ and mean at the origin.

Figure 6: Hellinger distance between the projection filter densities for the
Van der Pol oscillator and the densities obtained from the Crank–Nicolson
solution to the Kushner–Stratonovich equation (solid black) along with the
Hellinger distances between the Gaussian kernel density estimates computed
from the particle filter samples and the densities obtained by the
Crank–Nicolson solution (dashed black).

Figure 6 shows time of evolution of the Hellinger distance between the
Crank–Nicolson solution to the Kushner–Stratonovich SPDE and the projection
filter, and the Hellinger distance between the Crank–Nicolson solution and the
particle filter (estimated using a Gaussian kernel density estimate). It can
be seen that the Hellinger distances of the projection filter are
significantly smaller than those of the particle filter. Figures 7 show the
evolution of the approximate densities obtained by the Crank–Nicolson method,
particle filter, and the projection filter. There figures confirm that the
projection filter approximates the Kushner–Stratonovich solution better than
the particle filter. We note though that the Crank–Nicolson solution to the
Kushner–Stratonovich equation is only stable up to around 1.75s. After that, a
high frequency oscillation starts to appear (not shown). Thus we believe that
the Crank–Nicolson solution of the Kushner–Stratonovich equation is no longer
accurate after 1.75s. This numerical stability issue of the Crank–Nicolson
scheme for the Kushner–Stratonovich equation might also contribute to the
Hellinger distance between the finite-difference approximation of the
Kushner–Stratonovich equation and the density of the projection filter in
Figure 6.

A comparison of the different moment approximations is shown in Figures 8, and
9. The figures show the absolute difference between these expected values in
log scale, where $\mathbb{E}_{a},\mathbb{E}_{b},$ and $\mathbb{E}_{c}$ are the
expectation operator where the densities are obtained using Crank–Nicolson
scheme, projection filter, and the particle filter, respectively. Figures 8,
and 9 clearly demonstrate that in terms of expectations of the natural
statistics and their respective absolute errors, in the present setup, the
projection filter offers a much better approximation to the problem than the
particle filter.

Figure 7: Comparison of solution of the Crank–Nicolson solution to the
Kushner–Stratonovich equation, the densities of the projection filter for Van
der Pol dynamics (25), and Gaussian kernel density estimates from the solution
of the particle filter. The Gaussian kernel density estimate smoothing
coefficient is set to be proportional to the Cholesky factor of the sample
covariance matrix.

Figure 8: Absolute differences of the approximate expected values of the
natural statistics. In these figures $\mathbb{E}_{a}$ is the expectation
operator where the densities are obtained using Crank–Nicolson scheme,
$\mathbb{E}_{b}$ uses densities obtained using projection filter, and
$\mathbb{E}_{c}$ is the expectation operator where the densities are the
empirical densities of the particle filter, respectively.

Figure 9: Continuation of Figure 8 for different sets of natural statistics.

##  7 Conclusions and discussion

We have proposed a novel implementation of the projection filter for a class
of dynamical systems. In this new approach, we have successfully removed some
difficulties in projection filter implementation by calculating all the
moments and the Fisher information matrix using automatic differentiation
knowing that they can be obtained as partial derivatives of the log partition
function. We calculate the log-partition function via Gauss–Chebyshev
quadrature in unidimensional problems and via sparse-grid integration in
multidimensional problems. We also replace a recursion procedure typically
used to obtain the higher moments – which is valid only for unidimensional
problems – with a simple partial derivative evaluation of the log-partition
function that works for multidimensional problems. We have shown
mathematically that the resulting approximation error of the log-partition
function and its first and second derivatives have linear relations to the
integration error of partition functions and converge to zero as the number of
integration grid points approaches infinity. This result enables us to focus
only on choosing integration grid such that the integration error of the
partition function is acceptable.

We have shown the effectiveness of the proposed method by comparing it to a
finite-difference approximation of the Kushner–Stratonovich equation and a
bootstrap particle filter. In particular, for a unidimensional problem we have
shown that our projection filter implementation can be constructed using as
low as $12$ integration nodes, while still being very close to the finite-
difference solution to the Kushner–Stratonovich equation.

We also compared sparse-grid and quasi Monte Carlo (qMC) implementations of
the projection filter to the exact Kalman–Bucy filter solution in a two
dimensional linear dynamic model. We found that the sparse-grid based
projection filter approximates the Kalman–Bucy filter solution more accurately
than qMC. With the sparse-grid level set to 5, the sparse-grid based
projection filter has a practically negligible Hellinger distance to the
Kalman-Bucy solution after $t=0.4$, while the Hellinger distances of the
densities obtained with qMC integration continue to oscillate around
$10^{-4}$.

We have also applied the projection filter to non-linear/non-Gaussian
partially observed two-dimensional Van der Pol oscillator model. For the Van
der Pol oscillator filtering problem, we have shown that the projection filter
consructed with a level 8 sparse grid gives posterior density approximations
that are closer to the Crank–Nicolson solution to the Kushner–Stratonovich
equation than a particle filter with $4\times 10^{4}$ particles. In addition
to being lightweight, our projection filter implementation does not suffer
from numerical stability issues as the finite-difference or Crank–Nicolson
approximations to the Kushner–Stratonovich equation.

Thanks to the sparse-grid integration, our projection filter algorithm is
suitable for high-dimensional problems. In addition, if an increased accuracy
is needed, it would be possible to employ adaptive sparse-grid techniques for
more precise integration of the log partition function. In the future, we will
investigate the possibility of using parametric bijections from the canonical
cube to $\mathbb{R}^{d}$ to further reduce the number of quadrature nodes
while maintaining the same accuracy level.

##  8 Acknowledgements

Muhammad Emzir would like to express his gratitude to the KFUPM Dean of
Research Oversight and Coordination for the SR211015 grant.

## References

  * Amari (1985) S. Amari.  _Differential-Geometrical Methods in Statistics_.  Springer New York, 1985.
  * Amari and Nagaoka (2000) S. Amari and H. Nagaoka.  _Methods of Information Geometry_.  Number 191 in Translations of Mathematical Monographs. American Mathematical Society, 2000.
  * Armstrong and Brigo (2016) J. Armstrong and D. Brigo.  Nonlinear filtering via stochastic PDE projection on mixture manifolds in $\mathcal{L}^{2}$ direct metric.  _Mathematics of Control, Signals, and Systems_ , 28(1):5, 2016.
  * Armstrong and Brigo (2018) J. Armstrong and D. Brigo.  Intrinsic stochastic differential equations as jets.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 474(2210):20170559, 2018.
  * Armstrong and Brigo (2019) J. Armstrong and D. Brigo.  Optimal approximation of sdes on submanifolds: the Itô-vector and Itô-jet projections.  _Proceedings of the London Mathematical Society_ , 119(1):176–213, 2019.
  * Armstrong et al. (2018) J. Armstrong, D. Brigo, and E. R. Ferrucci.  Optimal approximation of SDEs on submanifolds: the Itô-vector and Itô-jet projections.  _Proceedings of the London Mathematical Society_ , 119(1):176–213, 2018.
  * Azimi-Sadjadi and Krishnaprasad (2005) B. Azimi-Sadjadi and P. Krishnaprasad.  Approximate nonlinear filtering and its application in navigation.  _Automatica_ , 41(6):945–956, 2005.
  * Baek and Bang (2013) K. Baek and H. Bang.  Adaptive sparse grid quadrature filter for spacecraft relative navigation.  _Acta Astronautica_ , 87:96–106, 2013.
  * Bain and Crisan (2009) A. Bain and D. Crisan.  _Fundamentals of Stochastic Filtering_.  Springer New York, 2009.
  * Bao et al. (2014) F. Bao, Y. Cao, C. Webster, and G. Zhang.  A hybrid sparse-grid approach for nonlinear filtering problems based on adaptive-domain of the Zakai equation approximations.  _SIAM/ASA Journal on Uncertainty Quantification_ , 2(1):784–804, 2014.
  * Barthelmann et al. (2000) V. Barthelmann, E. Novak, and K. Ritter.  High dimensional polynomial interpolation on sparse grids.  _Advances in Computational Mathematics_ , 12(4):273–288, 2000.
  * Beskos et al. (2017) A. Beskos, D. Crisan, A. Jasra, K. Kamatani, and Y. Zhou.  A stable particle filter for a class of high-dimensional state-space models.  _Advances in Applied Probability_ , 49(1):24–48, 2017.
  * Bradbury et al. (2018) J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.
  * Brigo (2022) D. Brigo.  Optimal projection filters, 2022.  arXiv:2205.01594.
  * Brigo et al. (1995) D. Brigo, B. Hanzon, and F. L. Gland.  A differential geometric approach to nonlinear filtering: the projection filter.  Research Report 2598, INRIA Rennes - Bretagne Atlantique, 1995\.
  * Brigo et al. (1998) D. Brigo, B. Hanzon, and F. L. Gland.  A differential geometric approach to nonlinear filtering: the projection filter.  _IEEE Transactions on Automatic Control_ , 43(2):247–252, 1998.
  * Brigo et al. (1999) D. Brigo, B. Hanzon, and F. L. Gland.  Approximate nonlinear filtering by projection on exponential manifolds of densities.  _Bernoulli_ , 5(3):495, 1999.
  * Bungartz and Griebel (2004) H.-J. Bungartz and M. Griebel.  Sparse grids.  _Acta Numerica_ , 13:147–269, may 2004.
  * Calin and Udrişte (2014) O. Calin and C. Udrişte.  _Geometric Modeling in Probability and Statistics_.  Springer International Publishing, 2014.
  * Candy (2016) J. V. Candy.  _Bayesian Signal Processing: Classical, Modern, and Particle Filtering Methods_.  IEEE Press / Wiley & Sons, Inc, 2nd edition, 2016.
  * Ceci and Colaneri (2013) C. Ceci and K. Colaneri.  The Zakai equation of nonlinear filtering for jump-diffusion observations: Existence and uniqueness.  _Applied Mathematics & Optimization_, 69(1):47–82, 2013.
  * Chopin (2020) N. Chopin.  _An Introduction to Sequential Monte Carlo_.  Springer, 2020.
  * Chui (1972) C. K. Chui.  Concerning Gaussian–Chebyshev quadrature errors.  _SIAM Journal on Numerical Analysis_ , 9(2):237–240, 1972.
  * Crisan and Rozovskii (2011) D. Crisan and B. Rozovskii.  _The Oxford Handbook of Nonlinear Filtering_.  Oxford University Press, 2011.
  * Delvos (1982) F.-J. Delvos.  d-Variate boolean interpolation.  _Journal of Approximation Theory_ , 34(2):99–114, 1982.
  * Doucet (2001) A. Doucet.  _Sequential Monte Carlo Methods in Practice_.  Springer New York, 2001.
  * Gao et al. (2019) Q. Gao, G. Zhang, and I. R. Petersen.  An exponential quantum projection filter for open quantum systems.  _Automatica_ , 99:59–68, 2019.
  * Gao et al. (2020a) Q. Gao, D. Dong, I. R. Petersen, and S. X. Ding.  Design of a Quantum Projection Filter.  _IEEE Transactions on Automatic Control_ , 65(8):3693–3700, 2020a.
  * Gao et al. (2020b) Q. Gao, G. Zhang, and I. R. Petersen.  An improved quantum projection filter.  _Automatica_ , 112:108716, 2020b.
  * Gerstner and Griebel (1998) T. Gerstner and M. Griebel.  Numerical integration using sparse grids.  _Numerical Algorithms_ , 18(3/4):209–232, 1998\.
  * Griebel and Oettershagen (2014) M. Griebel and J. Oettershagen.  Dimension-adaptive sparse grid quadrature for integrals with boundary singularities.  In J. Garcke and D. Pflüger, editors, _Sparse Grids and Applications - Munich 2012_ , pages 109–136. Springer International Publishing, Cham, 2014.
  * Griewank and Walther (2008) A. Griewank and A. Walther.  _Evaluating Derivatives : Principles and Techniques of Algorithmic Differentiation_.  SIAM, 2nd edition, 2008.
  * Hanzon and Hut (1991) B. Hanzon and R. Hut.  New results on the projection filter.  In _Proceedings of the first European Control Conference_ , volume 1, pages 623–628, 1991.
  * Ito and Xiong (2000) K. Ito and K. Xiong.  Gaussian filters for nonlinear filtering problems.  _IEEE Transactions on Automatic Control_ , 45(5):910–927, may 2000.
  * Jazwinski (1970) A. Jazwinski.  _Stochastic Processes and Filtering Theory_.  Academic Press, New York, 1970.
  * Jia et al. (2012) B. Jia, M. Xin, and Y. Cheng.  Sparse-grid quadrature nonlinear filtering.  _Automatica_ , 48(2):327–341, 2012.
  * Judd and Skrainka (2011) K. L. Judd and B. Skrainka.  High performance quadrature rules: how numerical integration affects a popular model of product differentiation.  Technical report, Centre for Microdata Methods and Practice, 2011.
  * Julier and Uhlmann (1997) S. J. Julier and J. K. Uhlmann.  New extension of the Kalman filter to nonlinear systems.  In I. Kadar, editor, _Signal Processing, Sensor Fusion, and Target Recognition VI_ , 1997.
  * Koyama (2018) S. Koyama.  Projection smoothing for continuous and continuous-discrete stochastic dynamic systems.  _Signal Processing_ , 144:333–340, 2018.
  * Kronrod (1966) A. S. Kronrod.  Nodes and weights of quadrature formulas.  _Mathematics of Computation_ , 20(93):184, 1966\.
  * Kushner (1967a) H. Kushner.  Nonlinear filtering: The exact dynamical equations satisfied by the conditional mode.  _IEEE Transactions on Automatic Control_ , 12(3):262–267, 1967a.
  * Kushner (1967b) H. Kushner.  Approximations to optimal nonlinear filters.  _IEEE Transactions on Automatic Control_ , 12(5):546–556, 1967b.
  * Kushner (1967c) H. J. Kushner.  Dynamical equations for optimal nonlinear filtering.  _Journal of Differential Equations_ , 1967c.
  * Laurie (1997) D. P. Laurie.  Calculation of Gauss-Kronrod quadrature rules.  _Mathematics of Computation_ , 66(219):1133–1146, 1997.
  * Leobacher and Pillichshammer (2014) G. Leobacher and F. Pillichshammer.  _Introduction to Quasi-Monte Carlo Integration and Applications_.  Birkhäuser, 2014.
  * Liptser and Shiryaev (2001) R. Liptser and A. Shiryaev.  _Statistics of Random Processes_.  Springer Berlin Heidelberg, 2001.
  * Liptser and Shiryaev (2010) R. Liptser and A. Shiryaev.  _Statistics of Random Processes II_.  Springer Berlin Heidelberg, 2010.
  * Ma and Zabaras (2009) X. Ma and N. Zabaras.  An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations.  _Journal of Computational Physics_ , 228(8):3084–3113, 2009.
  * Mason and Handscomb (2003) J. Mason and D. Handscomb.  _Chebyshev Polynomials_.  Chapman & Hall/CRC, 2003.
  * Maurel and Michel (1984) M. C. Maurel and D. Michel.  Des resultats de non existence de filtre de dimension finie.  _Stochastics_ , 13(1-2):83–102, 1984.
  * Meurer et al. (2017) A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, Š. Roučka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman, and A. Scopatz.  SymPy: Symbolic computing in Python.  _PeerJ Computer Science_ , 3:e103, 2017.  ISSN 2376-5992.
  * Novak and Ritter (1996) E. Novak and K. Ritter.  High dimensional integration of smooth functions over cubes.  _Numerische Mathematik_ , 75(1):79–97, 1996.
  * Patterson (1968) T. N. L. Patterson.  The optimum addition of points to quadrature formulae.  _Mathematics of Computation_ , 22(104):847–847, 1968.
  * Poterjoy (2015) J. Poterjoy.  A localized particle filter for high-dimensional nonlinear systems.  _Monthly Weather Review_ , 144(1):59–76, 2015\.
  * Radhakrishnan et al. (2016) R. Radhakrishnan, A. K. Singh, S. Bhaumik, and N. K. Tomar.  Multiple sparse-grid Gauss-Hermite filtering.  _Applied Mathematical Modelling_ , 40(7-8):4441–4450, 2016.
  * Rall (1981) L. B. Rall.  _Automatic Differentiation: Techniques and Applications_.  Number 120 in Lecture Notes in Computer Science. Springer, 1981.
  * Riess (1971) R. D. Riess.  A note on error bounds for Gauss–Chebyshev quadrature.  _SIAM Journal on Numerical Analysis_ , 8(3):509–511, sep 1971.
  * Särkkä (2013) S. Särkkä.  _Bayesian Filtering and Smoothing_.  Cambridge University Press, 2013.
  * Särkkä and Sarmavuori (2013) S. Särkkä and J. Sarmavuori.  Gaussian filtering and smoothing for continuous-discrete dynamic systems.  _Signal Processing_ , 93(2):500–510, feb 2013\.
  * Singh et al. (2018) A. K. Singh, R. Radhakrishnan, S. Bhaumik, and P. Date.  Adaptive sparse-grid Gauss-Hermite filter.  _Journal of Computational and Applied Mathematics_ , 342:305–316, 2018.
  * Smith (1985) G. D. Smith.  _Numerical Solution of Partial Differential Equations : Finite Difference Methods_.  Clarendon Press Oxford University Press, Oxford Oxfordshire New York, 1985\.
  * Smolyak (1963) S. A. Smolyak.  Quadrature and interpolation formulas for tensor products of certain classes of functions.  _Dokl. Akad. Nauk SSSR_ , 148:1042–1045, 1963.
  * Snyder et al. (2008) C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson.  Obstacles to high-dimensional particle filtering.  _Monthly Weather Review_ , 136(12):4629–4640, 2008.
  * Stoyanov (2015) M. Stoyanov.  User manual: TASMANIAN sparse grids.  Technical Report ORNL/TM-2015/596, Oak Ridge National Laboratory, One Bethel Valley Road, Oak Ridge, TN, 2015.
  * Tronarp and Särkkä (2019) F. Tronarp and S. Särkkä.  Updates in Bayesian filtering by continuous projections on a manifold of densities.  In _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019.
  * van Handel and Mabuchi (2005) R. van Handel and H. Mabuchi.  Quantum projection filter for a highly nonlinear model in cavity QED.  _J. Opt. B: Quantum Semiclass. Opt._ , 7(10):S226–S236, 2005.  ISSN 1464-4266.
  * Wasilkowski and Woźniakowski (1995) G. Wasilkowski and H. Woźniakowski.  Explicit cost bounds of algorithms for multivariate tensor product problems.  _Journal of Complexity_ , 11(1):1–56, 1995.
  * Winshcel and Krätzig (2010) V. Winshcel and M. Krätzig.  Solving, estimating, and selecting nonlinear dynamic models without the curse of dimensionality.  _Econometrica_ , 78(2):803–821, 2010.
  * Wonham (1963) W. M. Wonham.  Stochastic problems in optimal control.  _Proceedings of the IEEE_ , 51(3):530–530, 1963\.
  * Zakai (1969) M. Zakai.  On the optimal filtering of diffusion processes.  _Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete_ , 11(3):230–243, 1969.