2.1. Inner product

We endow the vector space with an inner product [Reference van den Boogaart, Egozcue and Pawlowsky-Glahn43] defined as

(8) \begin{equation} \left \langle{p_1},{p_2}\right \rangle = \frac{1}{2} \int _{\mathcal{X}} \int _{\mathcal{X}} \ln \left ( \frac{p_1(\textbf{x})}{p_1(\textbf{y})} \right ) \ln \left ( \frac{p_2(\textbf{x})}{p_2(\textbf{y})} \right ) \nu (\textbf{x}) \, \nu (\textbf{y}) \, d\textbf{x} \, d\textbf{y}, \end{equation}

where $\nu (\!\cdot\!)$ is again a density function corresponding to an appropriate measure for $\mathcal{X}$ . Notably, we see that because of the way the inner product is defined the normalization constants, $c_1$ and $c_2$ , associated with $p_1$ and $p_2$ play no role.

Because $\nu$ is a valid PDF, we can also write the inner product in (8) as

(9) \begin{equation} \left \langle{p_1,}\;{p_2}\right \rangle = \mathbb{E}_\nu \!\left [ \ln p_1 \ln p_2 \right ] - \mathbb{E}_\nu \!\left [ \ln p_1 \right ] \mathbb{E}_\nu \!\left [ \ln p_2 \right ], \end{equation}

where $\mathbb{E}_\nu [\cdot]$ is the expectation with respect to $\nu$ . To be clear, when we use expectations the argument, $f(\textbf{x})$ , and the measure (a PDF), $\nu (\textbf{x})$ , are defined over the same space, $\mathcal{X}$ :

(10) \begin{equation} \mathbb{E}_{\nu (\textbf{x})}[f(\textbf{x})] = \int _{{\mathcal{X}}} f(\textbf{x}) \nu (\textbf{x}) \, d\textbf{x}, \end{equation}

although sometimes we will abbreviate this as $\mathbb{E}_\nu [f]$ . In this work, we shall always take the measure to be a PDF (and from ${\mathcal{B}}^2$ ); however, we shall refer to it as the measure to distinguish it from the other densities involved. Following [Reference van den Boogaart, Egozcue and Pawlowsky-Glahn43], then, we can claim that ${\mathcal{B}}^2$ with inner product (8) forms a separable Hilbert space, which is referred to as a Bayesian Hilbert space. We shall sometimes briefly refer to it as a Bayesian space or Bayes space.

2.2. Information and divergence

The norm [Reference van den Boogaart, Egozcue and Pawlowsky-Glahn43] of $p \in{\mathcal{B}}^2$ can be taken as $\|p\| ={ \langle{p},{p} \rangle }^{1/2}$ . Accordingly, we can define the distance between two members of ${\mathcal{B}}^2$ , $p$ and $q$ , simply as $d(p,q) = \|p\ominus q\|$ , which induces a metric on Bayes space.

The norm of $p$ can be used to express the information content of the PDF (if it is in ${\mathcal{B}}^2$ ). In fact, we shall define

(11) \begin{equation} I(p) = \frac{1}{2}\|p\|^2 = \frac{1}{2}\left \langle{p},{p}\right \rangle \end{equation}

as the information Footnote ³ in $p$ . (The reason for the factor of $\frac{1}{2}$ will become evident.) As an example, consider $p ={\mathcal{N}}(\mu,\sigma ^2)$ (over the domain $\mathbb{R}$ ) and measure $\nu ={\mathcal{N}}(0,1)$ . The information is $I(p) = (1+2\mu ^2)/4\sigma ^4$ . The smaller the variance, the larger the information indicating that the PDF concentrates its probability mass more tightly about its mean; that is, we know better where to expect the state so we may say that we have more information about it.

We shall furthermore find it useful to define a divergence between two members of ${\mathcal{B}}^2$ , $p$ and $q$ , as

(12) \begin{equation} I(p\ominus q) = \frac{1}{2}\left \langle{p\ominus q},{p\ominus q}\right \rangle. \end{equation}

This is the information contained in the difference of $p$ and $q$ . Unlike the KL divergence, this divergence is symmetric in $p$ and $q$ and quadratic in Bayesian space. Clearly, $p = q$ if and only if $I(p\ominus q) = 0$ . Geometrically, the divergence is (half) the squared Euclidean distance between $p$ and $q$ in Bayes space.

2.3. Stochastic derivative

Accompanying this algebra is a functional calculus. Consider $p(\textbf{x}|\theta ) \in{\mathcal{B}}^2$ depending continuously on some parameter $\theta$ . We define the stochastic partial derivative of $p$ with respect to $\theta$ as [Reference Barfoot and D’Eleuterio13, Reference Egozcue, Pawlowsky-Glahn, Tolosana-Delgado, Ortego and van den Boogaart22]

(13) \begin{equation} \frac{\eth p}{\eth \theta } = \lim _{\lambda \rightarrow 0} \frac{1}{\lambda }\cdot \bigl ( p(\textbf{x}|\theta + \lambda ) \ominus p(\textbf{x}|\theta ) \bigr ). \end{equation}

Note that the result of this operation remains an element in ${\mathcal{B}}^2$ . We can also define directional derivatives and a gradient operator, but these will not be required here.

3.1. Projections

Given a subspace $\mathcal{Q}$ of ${\mathcal{B}}^2$ and $p \in{\mathcal{B}}^2$ , the $q^\star \in{\mathcal{Q}}$ that minimizes the distance to, as well as the divergence (12) from, $p$ is the projection of $p$ onto $\mathcal{Q}$ , that is,

(17) \begin{equation} q^\star = \underset{{\mathcal{Q}}}{\text{proj}}\, p. \end{equation}

As in Euclidean geometry, we can view $p$ as being decomposed into a component $p_\|$ lying in $\mathcal{Q}$ and a component $p_\perp$ perpendicular to it; therefore $q^\star = p_\|$ (see Fig. 1).

Figure 1. Projection onto a subspace, $\mathcal{Q}$ , of the Bayesian Hilbert space, ${\mathcal{B}}^2$ .

The coordinates of $q^\star$ can be calculated as

(18) \begin{equation} \boldsymbol \alpha ^\star = \left \langle{\textbf{b}},{\textbf{b}}\right \rangle ^{-1}\left \langle{\textbf{b}},{p}\right \rangle. \end{equation}

We may also write the projection as an outer product operation on $p$ , as detailed in Appendix B.

3.2. Example: one-dimensional Gaussian

To make the concept of Bayes space more tangible, consider the canonical one-dimensional Gaussian PDF defined over $x \in \mathbb{R}$ :

(19) \begin{equation} p(x) = \frac{1}{\sqrt{2 \pi \sigma ^2}} \exp \left ( -\frac{1}{2} \frac{(x-\mu )^2}{\sigma ^2} \right ), \end{equation}

where $\mu$ is the mean and $\sigma ^2$ the variance. In the language of ${\mathcal{B}}^2$ , we can write this as

(20) \begin{equation} p(x) = \underbrace{\left ( -\frac{\mu }{\sigma ^2}\right )}_{\alpha _1} \cdot \underbrace{\exp (\!-\!x)}_{b_1} \oplus \underbrace{\left (\frac{1}{\sqrt{2}\sigma ^2}\right )}_{\alpha _2} \cdot \underbrace{\exp \left ( -\frac{(x^2 - 1)}{\sqrt{2}} \right )}_{b_2} = \alpha _1 \cdot b_1 \oplus \alpha _2 \cdot b_2. \end{equation}

In other words, every Gaussian can be written as a linear combination of the two vectors, $b_1$ and $b_2$ , where the coefficients, $\alpha _1$ and $\alpha _2$ , depend on the mean and variance. Note, we can neglect the normalizing constant as equivalence is implied in the “ $=$ ” operator.

The choice of $b_1$ and $b_2$ is not arbitrary in this example. They constitute the first two basis vectors in an orthonormal basis for ${\mathcal{B}}^2$ , which can be established using the probabilist’s Hermite polynomials; Appendix C.1 provides the details of this Hermite basis. In fact, we can define a new space $\mathcal{G}$ as the span of these two basis vectors:

(21) \begin{equation}{\mathcal{G}} = \text{span}\left \{ b_1, b_2 \right \}, \end{equation}

which is a subspace of ${\mathcal{B}}^2$ . Importantly, every Gaussian PDF of the form in (19) is a member of $\mathcal{G}$ , but not every member of $\mathcal{G}$ is a valid Gaussian PDF. Only those members of $\mathcal{G}$ that have $\sigma ^2 \gt 0$ are valid Gaussian PDFs. We shall refer to $\mathcal{G}$ as the indefinite-Gaussian subspace of ${\mathcal{B}}^2$ while $\,\downarrow \!{{\mathcal{G}}} \subset{\mathcal{G}}$ will denote the positive-definite-Gaussian subset. Figure 2 shows the relationships between the various spaces and how we can view a Gaussian as a point in $\mathcal{G}$ .

3.3. Example: projecting to a Gaussian

Let us consider a simple one-dimensional, nonlinear estimation problem as a numerical example motivated by the type of inverse-distance nonlinearity found in a stereo camera model. This same experiment (with the same parameter settings) was used as a running example by [[Reference Barfoot10], Section 4]. We assume that our true state is drawn from a Gaussian prior:

(22) \begin{equation} x \sim \mathcal{N}(\mu _p, \sigma _p^2). \end{equation}

We then generate a measurement according to

(23) \begin{equation} z = \frac{fb}{x} + n, \quad n \sim \mathcal{N}(0,\sigma _r^2), \end{equation}

where $n$ is the measurement noise. The numerical values of the parameters used were

(24) \begin{equation} \begin{gathered} \mu _p = 20 \text{ [m]}, \quad \sigma _p^2 = 9 \text{ [m$^2$]}, \\[2pt] f = 400 \text{ [pixel]}, \quad b = 0.1 \text{ [m]}, \quad \sigma _r^2 = 0.09 \text{ [pixel$^2$]}. \end{gathered} \end{equation}

Figure 2. On the left is a depiction of the relationships between Bayes space, ${\mathcal{B}}^2$ , and the set of all PDFs. We see the subset of strictly positive PDFs, $\,\downarrow \!{{\mathcal{B}}^2}$ , the indefinite-Gaussian subspace, $\mathcal{G}$ , and the positive-definite-Gaussian subset, $\,\downarrow \!{{\mathcal{G}}}$ . The right image shows how a valid Gaussian PDF can be viewed as a point in a plane with coordinates that depend on its mean $\mu$ and variance $\sigma ^2$ ; only the open upper-half plane admits valid Gaussian PDFs since we must have $\sigma ^2 \gt 0$ .

Figure 3. An example of projecting a non-Gaussian posterior onto the indefinite-Gaussian subspace. The top panel shows the case where the measure associated with ${\mathcal{B}}^2$ was chosen to be the same as the (Gaussian) prior, $\nu (x) = \mathcal{N}(20,9)$ . The bottom panel does the same with a Gaussian measure selected to be closer to the posterior, $\nu (x) = \mathcal{N}(24,4)$ . We see that the Gaussian projection of the posterior is much closer to the true posterior in the bottom case.

The true posterior is given by

(25) \begin{equation} p(x|z) = \,\downarrow \!{\exp (\!-\phi (x))}, \quad \phi (x) = \underbrace{\frac{1}{2} \frac{(x - \mu _p)^2}{\sigma _p^2}}_{\text{prior}} + \underbrace{\frac{1}{2} \frac{\left (z - \frac{fb}{x}\right )^2}{\sigma _r^2}}_{\text{measurement}}. \end{equation}

This problem can also be viewed as the correction step of the Bayes filter [Reference Jazwinski26]: Start from a prior and correct it based on the latest (nonlinear) measurement.

We seek to find $q(x) \in \,\downarrow \!{{\mathcal{G}}}$ that is a good approximation to the true posterior $p(x|z)$ . To do this, we will simply project the posterior onto the indefinite-Gaussian subspace, using the method described in Section 3.1, and then normalize the result. Figure 3 shows the results of doing this for two cases that differ only in the measure $\nu$ that we associate with Bayes space. The expectations used in the projections were computed with generic numerical integration although, as discussed by [Reference Barfoot, Forbes and Yoon14], there are several other options including Gaussian quadrature. In the top case, the measure is chosen as the prior estimate, while in the bottom case it is chosen to be closer to the posterior. In both cases, we see that our projection method produces a valid PDF, but in the bottom case the result is much closer to the true posterior. This simple example provides motivation for the main point of this paper, which is that to use the tools of Bayes space effectively, we will seek to iteratively update the measure used to carry out our projections such that we can best approximate a posterior. Intuitively, this makes sense since the measure is providing a weighting to different parts of $\mathbb{R}$ so we would like to choose it to pay close attention where the posterior ends up.

4.1. Variation on the KL divergence

In variational Bayesian inference, we seek to find an approximation, $q$ , from some family of distributions constituting a subspace $\mathcal{Q}$ to the true Bayesian posterior $p \in{\mathcal{B}}^2$ . In general,

(26) \begin{equation}{\mathcal{Q}} \subseteq{\mathcal{B}}^2, \end{equation}

where equality will always ensure that $q=p$ will match the posterior exactly. But ${\mathcal{B}}^2$ is infinite-dimensional and, in practice, ${\mathcal{Q}} \subset{\mathcal{B}}^2$ is a finite-dimensional subspace.

There are many possible divergences that can be defined to characterize the “closeness” of $q$ to $p$ including the KL divergence [Reference Kullback and Leibler28], Bregman divergence [Reference Bregman19], Wasserstein divergence/Earth mover’s distance [Reference Monge35], and Rényi divergence [Reference Rényi38]. We shall focus on the KL divergence, which is defined as

(27) \begin{equation} \text{KL}(q || p) = - \int _{\mathcal{X}} q(\textbf{x}) \ln \left ( \frac{p(\textbf{x} | \textbf{z})}{q(\textbf{x})} \right ) \, d\textbf{x} = -\mathbb{E}_q [ \ln p - \ln q]. \end{equation}

Sometimes the reverse of this is used: $\text{KL}(p || q)$ . Note, we show the divergence with respect to the posterior, $p(\textbf{x} | \textbf{z})$ , but in practice during the calculations we use that $p(\textbf{x} | \textbf{z}) = p(\textbf{x},\textbf{z})/ p(\textbf{z}) = \,\downarrow \!{p(\textbf{x},\textbf{z})}$ since the joint likelihood is easy to construct and then the $p(\textbf{z})$ can be dropped for it does not result in a KL term that depends on $\textbf{x}$ . We will generically use $p$ in what follows to keep the notation clean.

4.2. KL gradient

We assume a basis $B = \{b_1,b_2\cdots b_M\}$ for $\mathcal{Q}$ and we write $q$ as

(28) \begin{equation} q = \,\downarrow \!{} \bigoplus _{m=1}^M \alpha _m\cdot b_m. \end{equation}

We desire to minimize the KL divergence with respect to the coordinates $\alpha _m$ . The gradient of $\text{KL}(q || p)$ can be computed as follows:

(29) \begin{equation} \frac{\partial \text{KL}}{\partial \alpha _n} = -\int _{{\mathcal{X}}} \left (\frac{\partial q}{\partial \alpha _n}(\ln p - \ln q) - q\frac{\partial \ln q}{\partial \alpha _n}\right )d\textbf{x}. \end{equation}

Exploiting (E4) and (E6), this reduces to

(30) \begin{equation} \frac{\partial \text{KL}}{\partial \alpha _n} = -{\mathbb{E}}_q[\ln b_n(\ln p - \ln q)] +{\mathbb{E}}_q[\ln b_n]{\mathbb{E}}_q[\ln p - \ln q] = -\left \langle{b_n},{p\ominus q}\right \rangle _q \end{equation}

or, collecting these in matrix form,

(31) \begin{equation} \frac{\partial \text{KL}}{\partial \boldsymbol \alpha ^T} = -\left \langle{\textbf{b}},{p\ominus q}\right \rangle _q, \end{equation}

where $\textbf{b} = \begin{bmatrix} b_1 & b_2 & \cdots & b_M \end{bmatrix}^T$ . Implicit in this statement is that when employed as the measure, we have normalized the current approximation, $\,\downarrow \!{q^{(i)}}$ , since we always take the measure to be a valid PDF. The necessary condition for a minimum of the KL divergence is that the gradient is zero. Newton’s method suggests the manner in which we might iteratively solve for the optimal distribution. Following the established procedure, the iteration for the coordinates is given by

(32) \begin{equation} \boldsymbol \alpha ^{(i+1)} = \boldsymbol \alpha ^{(i)} +{\textbf{H}^{(i)}}^{-1}\left \langle{\textbf{b}},{p\ominus q^{(i)}}\right \rangle _{q^{(i)}}, \end{equation}

where $\textbf{H}$ is the Hessian of the KL divergence.

4.3. KL Hessian

The $(m,n)$ entry of the Hessian is

(33) \begin{equation} \frac{\partial ^2\text{KL}}{\partial \alpha _m\partial \alpha _n} = -\frac{\partial }{\partial \alpha _m}\left \langle{b_n},{p\ominus q}\right \rangle _q. \end{equation}

This differentiation must take into account the effect of the “measure” $q$ . The product rule applies here and we can break down the differentiation as

(34) \begin{equation} \frac{\partial ^2\text{KL}}{\partial \alpha _m\partial \alpha _n} = -\left (\frac{\partial }{\partial \alpha _m}\left \langle{b_n},{p\ominus q}\right \rangle \right )_q - \left \langle{b_n},{p\ominus q}\right \rangle _{\partial q/\partial \alpha _m}, \end{equation}

the first term of which is to be read as the derivative of the inner product holding the measure fixed and the second of which deals with the derivative of the measure while holding the arguments of inner product fixed. The first term is

(35) \begin{equation} \left (\frac{\partial }{\partial \alpha _m}\left \langle{b_n},{p\ominus q}\right \rangle \right )_q = \left (\frac{\partial }{\partial \alpha _m}\left \langle{b_n},{p}\right \rangle - \frac{\partial }{\partial \alpha _m}\sum _k\alpha _k\left \langle{b_n},{b_k}\right \rangle \right )_q = -\left \langle{b_n},{b_m}\right \rangle _q = -\left \langle{b_m},{b_n}\right \rangle _q. \end{equation}

As shown in Appendix E.2, the second becomes

(36) \begin{equation} \left \langle{b_n},{p\ominus q}\right \rangle _{\partial q/\partial \alpha _m} = \left \langle{b_n},{\frac{\partial \ln q}{\partial \alpha _m}\cdot (p\ominus q)}\right \rangle _q -{\mathbb{E}}_q[\ln p - \ln q]\left \langle{b_m},{b_n}\right \rangle _q. \end{equation}

We advise that the coefficient $\partial \ln q/\partial \alpha _m$ of $p\ominus q$ is in fact a function of the state and as such cannot be transferred to the other argument of the inner product as would be possible for a scalar in the field $\mathbb{R}$ . We also recognize the factor of the last term as $\text{KL}(q||p)$ . Therefore, substituting (35) and (36) into (34) yields

(37) \begin{equation} \frac{\partial ^2\text{KL}}{\partial \alpha _m\partial \alpha _n} = \left (1 - \text{KL}(q||p)\right )\left \langle{b_m},{b_n}\right \rangle _q - \left \langle{b_n},{\frac{\partial \ln q}{\partial \alpha _m}\cdot (p\ominus q)}\right \rangle _q. \end{equation}

We observe that the second term on the right-hand side is symmetric in the indices as the substitution of (E6) will attest. In matrix form, the Hessian is

(38) \begin{equation} \textbf{H} = \frac{\partial ^2\text{KL}}{\partial \boldsymbol \alpha ^T\partial \boldsymbol \alpha } = \left (1 - \text{KL}(q||p)\right )\textbf{I}_{\boldsymbol \alpha } - \left \langle{\textbf{b}},{\frac{\partial \ln q}{\partial \boldsymbol \alpha ^T}\cdot (p\ominus q)}\right \rangle _q, \end{equation}

where $\textbf{I}_{\boldsymbol \alpha }$ is the Fisher information matrix (FIM) or Gram matrix and is described in detail in Appendix E.1. Newton’s method (32) can now be implemented. But the Hessian bears a closer look.

The Hessian can also be explicitly written as

(39) \begin{align} \frac{\partial ^2\text{KL}}{\partial \alpha _m\partial \alpha _n} &= \left \langle{b_m},{b_n}\right \rangle _q -{\mathbb{E}}_q[\ln b_m\ln b_n(\ln p - \ln q)] \nonumber \\[5pt] &+{\mathbb{E}}_q[\ln b_m\ln b_n]{\mathbb{E}}_q[\ln p - \ln q] +{\mathbb{E}}_q[\ln b_n]{\mathbb{E}}_q[\ln b_m(\ln p - \ln q)] \nonumber \\[5pt] &+{\mathbb{E}}_q[\ln b_m]{\mathbb{E}}_q[\ln b_n(\ln p - \ln q)] - 2{\mathbb{E}}_q[\ln b_m]{\mathbb{E}}_q[\ln b_n]{\mathbb{E}}_q[\ln p - \ln q], \end{align}

the terms of which can be collected as

(40) \begin{equation} \frac{\partial ^2\text{KL}}{\partial \alpha _m\partial \alpha _n} = \left \langle{b_m},{b_n}\right \rangle _q + \left \langle{-b_{mn} \oplus {\mathbb{E}}_q[\ln b_n]\cdot b_m \oplus{\mathbb{E}}_q[\ln b_m]\cdot b_n},{p\ominus q}\right \rangle _q, \end{equation}

where $b_{mn} = \exp (\ln b_m\ln b_n)$ . The symmetry in the Hessian is plainly evident in this version.

4.4. Iterative projection

In the vicinity of the optimal distribution, with a sufficiently large subspace $\mathcal{Q}$ , we may expect $p\ominus q$ to be small almost everywhere. This makes all the terms in the Hessian of first order except $\textbf{I}_{\boldsymbol \alpha }$ , which is of zeroth order. The gradient (68) is also of first order. Thus to keep Newton’s descent to this order, we may approximate the Hessian as $\textbf{H} \simeq \textbf{I}_{\boldsymbol \alpha }$ and the iterative procedure (32) becomes simply

(41) \begin{equation} \boldsymbol \alpha ^{(i+1)} = \boldsymbol \alpha ^{(i)} +{\textbf{I}_{\boldsymbol \alpha }^{(i)}}^{-1}\left \langle{\textbf{b}},{p\ominus q^{(i)}}\right \rangle _{q^{(i)}}. \end{equation}

However, as $q^{(i)} = \,\downarrow \!{\bigoplus }_m \alpha _m^{(i)}\cdot b_m$ ,

(42) \begin{equation} \left \langle{\textbf{b}},{p\ominus q^{(i)}}\right \rangle _{q^{(i)}} = \left \langle{\textbf{b}},{p}\right \rangle _{q^{(i)}} - \left \langle{\textbf{b}},{\textbf{b}}\right \rangle _{q^{(i)}}\boldsymbol \alpha ^{(i)} = \left \langle{\textbf{b}},{p}\right \rangle _{q^{(i)}} - \textbf{I}_{\boldsymbol \alpha }^{(i)}\boldsymbol \alpha ^{(i)}. \end{equation}

Hence, (41) becomes

(43) \begin{equation} \boldsymbol \alpha ^{(i+1)} ={\textbf{I}_{\boldsymbol \alpha }^{(i)}}^{-1}\left \langle{\textbf{b}},{p}\right \rangle _{q^{(i)}}. \end{equation}

The iterative update to $q$ is $q^{(i+1)} = \,\downarrow \!{\bigoplus }_m \alpha _m^{(i+1)}\cdot b_m$ . That is, the procedure can be viewed as an iterative projection,

(44) \begin{equation} q^{(i+1)} = \underset{({\mathcal{Q}}, q^{(i)})}{\,\downarrow \!{} \text{proj}} p, \end{equation}

where we explicitly indicate that we normalize the result as the output of our algorithm should be a PDF. Figure 4 depicts the scheme. The procedure is essentially the application of Newton’s method on the KL divergence with the Hessian approximated as the FIM. This is precisely the approximation made in natural gradient descent [Reference Amari4]. In our Bayesian space, the operating point of the Newton step becomes the measure for the inner product. This highlights a key aspect of using the algebra associated with a Bayesian space. It recognizes the dual role of $q$ : On the one hand, it is the approximating PDF, and on the other it serves as a measure that weights the difference between the approximation and the approximated.

Figure 4. Iterative projection onto a sequence of Bayesian Hilbert spaces, $(\mathcal{Q},q^{(i)})$ .

Convergence of iterative projection is guaranteed if the Hessian is positive definite. Provided that the subspace is large enough, we can expect convergence when we begin in a neighborhood of optimal $q$ where the first-order terms in the Hessian are sufficiently small.

It is notable that each step of the iterative projection is equivalent to the local minimization of the divergence $I(p\ominus q)$ with the measure fixed at $q^{(i)}$ because

(45) \begin{equation} I\!\left (p\ominus ( q^{(i)}\oplus \delta q)\right ) = I(p\ominus q^{(i)}) + \delta \boldsymbol \alpha ^T\left (\frac{\partial I}{\partial \boldsymbol \alpha ^T}\right )_{q^{(i)}} + \frac{1}{2}\delta \boldsymbol \alpha ^T \left (\frac{\partial ^2 I}{\partial \boldsymbol \alpha ^T\partial \boldsymbol \alpha }\right )_{q^{(i)}} \delta \boldsymbol \alpha, \end{equation}

where $\delta \boldsymbol \alpha = \boldsymbol \alpha ^{(i+1)} - \boldsymbol \alpha ^{(i)}$ and

(46) \begin{equation} \left (\frac{\partial I}{\partial \boldsymbol \alpha ^T}\right )_{q^{(i)}} = -\left \langle{\textbf{b}},{p\ominus q^{(i)}}\right \rangle _{q^{(i)}}, \quad \left (\frac{\partial ^2 I}{\partial \boldsymbol \alpha ^T\partial \boldsymbol \alpha }\right )_{q^{(i)}} = \left \langle{\textbf{b}},{\textbf{b}}\right \rangle _{q^{(i)}} \equiv \textbf{I}_{\boldsymbol \alpha }^{(i)}, \end{equation}

which are identical to the linearized forms for the KL divergence.

Throughout this section, we have assumed that the basis $B$ remains constant across iterations, but this need not be the case. We may also choose to update the basis along with the measure to maintain, for example, orthonormality. This is explored in the next example and further in Section 5 on Gaussian variational inference.

4.5. Example: iteratively projecting to a Gaussian

In the example of Section 3.3, we saw that selecting a measure that was closer to the posterior resulted in a projection that was also closer to the posterior. We now redo this example using the iterative projection concepts from this section. We will still project onto the indefinite-Gaussian subspace and employ a Gaussian measure, only now with each iteration the measure will be taken to be the (normalized) projection from the previous iteration.

We initialized the estimate to the prior, which corresponds to the first panel in Fig. 5. The next three panels show subsequent iterations of the estimate. The last panel shows the KL divergence between the estimate and the true posterior for 10 iterations. We see the estimate converged in a few iterations.

Figure 5. Example of iterative projection onto the indefinite-Gaussian subspace spanned by two Hermite basis functions, where the measure is taken to be the estimate $q^{(i)}$ at the previous iteration and the basis reorthogonalized at each iteration as described in Section 5. The estimate was initialized to the prior (first panel) and then iteratively updated (next three panels). The last panel shows the KL divergence between the estimate and the true posterior for 10 iterations, with convergence occurring at approximately 5 iterations.

Note that as the measure changes from one iteration to the next, we then have to update the basis to retain the desired orthogonality. This can be accomplished by using the reparameterization “trick” (see Appendix C.1) to adjust the basis to be orthogonal with respect to the current Gaussian measure.

4.6. Exploiting sparsity

One of the major advantages of thinking of ${\mathcal{B}}^2$ as a vector space with the definition of vector addition $\oplus$ is that Bayesian inference in general can be viewed as the addition of vectors. Consider the posterior $p(\textbf{x} | \textbf{z})$ where $\textbf{z}$ are some measurements. Bayes’ rule states that

(47) \begin{equation} p(\textbf{x} | \textbf{z}) = \frac{p(\textbf{z} | \textbf{x}) p(\textbf{x})}{p(\textbf{z})} = p(\textbf{z} | \textbf{x}) \oplus p(\textbf{x}), \end{equation}

where $p(\textbf{x})$ is a prior, $p(\textbf{z} | \textbf{x})$ is a measurement factor and, as mentioned earlier, we need not introduce the normalization constant $p(\textbf{z})$ explicitly when writing the posterior as a vector addition in Bayesian space. To be clear, addition is defined for two members of the same Bayes space and here we interpret $p(\textbf{z} | \textbf{x})$ as a function of $\textbf{x}$ since $\textbf{z}$ is a constant (e.g., a known measurement).

If we have several measurements that are statistically independent, then this can be factored as

(48) \begin{equation} p(\textbf{x} | \textbf{z}) = p(\textbf{x}) \oplus \bigoplus _{k=1}^K p(\textbf{z}_k | \textbf{x}_k), \end{equation}

where $\textbf{x}_k = \textbf{P}_k \textbf{x}$ is a subset of the variables in $\textbf{x}$ , $\textbf{P}_k$ is a projection matrix, and $\textbf{z}_k$ is the $k$ th measurement. This expresses sparsity in the state description and in the measurements. To keep the notation economical, we shall simply write

(49) \begin{equation} p = \bigoplus _{k=0}^K p_k, \end{equation}

where $p$ is the posterior and the $p_k$ comprise the prior and the measurements, corresponding to statistically independent data. In other words, the factorization becomes a summation in the Bayesian space ${\mathcal{B}}^2$ .

Now consider our projective approach to inference. As usual, given a subspace ${\mathcal{Q}} \subset{\mathcal{B}}^2$ , the optimal estimate to (49) is given by

(50) \begin{equation} q^\star = \underset{{\mathcal{Q}}}{\text{proj}} \, p = \underset{{\mathcal{Q}}}{\text{proj}} \,\bigoplus _{k=0}^K p_k = \bigoplus _{k=0}^K \underset{{\mathcal{Q}}}{\text{proj}} \, p_k. \end{equation}

That is, the projection of the sum is the sum of the projections. Each individual projection can be done separately because we are in a linear space. This is of enormous practical advantage because it means that we do not need all of $\mathcal{Q}$ to represent each projection.

We can see this more clearly by defining ${\mathcal{B}}^2_k \subset{\mathcal{B}}^2$ as the subspace corresponding to the variables $\textbf{x}_k$ . Then

(51) \begin{equation} p = \bigoplus _{k=0}^K p_k \in{\mathcal{B}}^2_0\oplus{\mathcal{B}}^2_1\oplus \cdots \oplus{\mathcal{B}}^2_K \subseteq{\mathcal{B}}^2. \end{equation}

In other words, $p$ is contained in the direct sum of the subspaces ${\mathcal{B}}^2_k$ . Each constituent part $p_k$ may be confined to a smaller subspace of ${\mathcal{B}}^2$ , depending on the variable dependencies in each term.

If we wish to project $p_k \in{\mathcal{B}}^2_k$ onto $\mathcal{Q}$ it will suffice to consider the projection on just ${\mathcal{Q}}_k ={\mathcal{B}}^2_k \cap{\mathcal{Q}}$ , that is,

(52) \begin{equation} \underset{{\mathcal{Q}}}{\text{proj}} \, p_k = \underset{{\mathcal{Q}}_k}{\text{proj}} \, p_k. \end{equation}

The subspace ${\mathcal{Q}}_k$ may, and ideally would, be smaller than $\mathcal{Q}$ . We may refer to ${\mathcal{Q}}_k$ as the marginal subspace of $\mathcal{Q}$ with respect to the subset of variables $\textbf{x}_k$ .

Therefore, the optimal estimate will be given by

(53) \begin{equation} q^\star = \underset{{\mathcal{Q}}}{\text{proj}} \, p = \bigoplus _{k=0}^K \underset{{\mathcal{Q}}_k}{\text{proj}} \, p_k. \end{equation}

Figure 6. Exploiting sparsity by projecting individual measurements onto marginal subspaces, $\mathcal{Q}_k$ , and then recombining the results.

This means that we can project the PDF associated with each measurement onto a smaller subspace and simply add up the estimates, lifting the overall estimate up into a potentially much larger space. Naturally, when employed in practice we will normalize $q^\star$ to ensure our algorithm outputs a valid PDF. The decomposition and reconstitution are illustrated in Fig. 6. Just as with the total posterior, we may describe $q^\star$ as being an element of a direct sum of the individual subspaces of $\mathcal{Q}$ , that is,

(54) \begin{equation} q^\star \in{\mathcal{Q}}_0 \oplus{\mathcal{Q}}_1 \oplus \cdots \oplus{\mathcal{Q}}_K \subseteq{\mathcal{Q}}. \end{equation}

The subspace sum may be substantially smaller than $\mathcal{Q}$ but again it will depend on the variable dependencies of each term.

This is the key result that allows most practical inference frameworks to function in a tractable way. Depending on the chosen basis for $\mathcal{Q}$ , many of the coordinates can potentially be zero and thus it will not be necessary to waste effort computing them or space storing them.

5.1. Projections

As mentioned at the end of the last section, we do not have to maintain the same basis from step to step as long as each basis spans the same subspace. This is a particularly useful maneuver when using the subspace $\mathcal{G}$ of indefinite Gaussians, which are discussed in detail for the multivariate case in Appendix D.1. Denote the mean and variance of $q^{(i)} \in{\mathcal{G}}$ as ${\boldsymbol{\mu }}^{(i)}$ and $\boldsymbol{\Sigma }^{(i)}$ and let the basis $\textbf{g}^{(i)}$ be defined as in (D3) and (D4). Note that this basis is orthonormal with respect to $q^{(i)}$ . As such, $\textbf{I}_{\boldsymbol \alpha }^{(i)} = \langle{\textbf{g}^{(i)}},{\textbf{g}^{(i)}} \rangle _{q^{(i)}} = \textbf{1}$ . Imagine the PDF to be approximated is expressed as $p = \,\downarrow \!{\exp (\!-\phi (\textbf{x}))} \in \,\downarrow \!{{\mathcal{B}}^2}$ . The coordinates resulting from the next projection are given by (D13), namely,

(55) \begin{equation} \begin{gathered} \boldsymbol \alpha _1^{(i+1)} = \left \langle{\textbf{g}_1^{(i)}},{p}\right \rangle ={\textbf{L}^{(i)}}^T{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi (\textbf{x})}{\partial \textbf{x}^T}\right ], \\[10pt] \boldsymbol \alpha _2^{(i+1)} = \left \langle{\textbf{g}_2^{(i)}},{p}\right \rangle = \sqrt{\text{$\frac{1}{2}$}\textbf{D}^T\textbf{D}}\, \text{vech}\left ({\textbf{L}^{(i)}}^T{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial ^2\phi (\textbf{x})}{\partial \textbf{x}^T\partial \textbf{x}}\right ]{\textbf{L}^{(i)}}\right ), \end{gathered} \end{equation}

where $\textbf{L}^{(i)}$ issues from the Cholesky decomposition of $\boldsymbol{\Sigma }^{(i)}$ .

The new iteration is

(56) \begin{equation} q^{(i+1)} = \underset{({\mathcal{G}}, q^{(i)})}{\,\downarrow \!{}\text{proj}} \,p = \,\downarrow \!{\exp }\left (\!-{\boldsymbol \alpha _1^{(i+1)}}^T{\boldsymbol{\gamma }}_1^{(i)} -{\boldsymbol \alpha _1^{(i+1)}}^T{\boldsymbol{\gamma }}_1^{(i)}\right ). \end{equation}

Using (D14), this becomes

(57) \begin{equation} q^{(i+1)} = \,\downarrow \!{\exp }\biggl (\!-(\textbf{x} -{\boldsymbol{\mu }}^{(i)})^T{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi (\textbf{x})}{\partial \textbf{x}^T}\right ] -\frac{1}{2}(\textbf{x} -{\boldsymbol{\mu }}^{(i)})^T{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial ^2\phi (\textbf{x})}{\partial \textbf{x}^T\partial \textbf{x}}\right ] (\textbf{x} -{\boldsymbol{\mu }}^{(i)})\biggr ), \end{equation}

which we may cast into the form,

(58) \begin{equation} q^{(i+1)} = \,\downarrow \!{\exp }\!\left (\!-\frac{1}{2}(\textbf{x} -{\boldsymbol{\mu }}^{(i+1)})^T{\boldsymbol{\Sigma }^{(i+1)}}^{-1} (\textbf{x} -{\boldsymbol{\mu }}^{(i+1)})\right ). \end{equation}

Herein

(59a) \begin{equation} {\boldsymbol{\Sigma }^{(i+1)}}^{-1} ={\mathbb{E}}_{q^{(i)}}\left [\frac{\partial ^2\phi (\textbf{x})}{\partial \textbf{x}^T\partial \textbf{x}}\right ], \end{equation}

(59b) \begin{equation} {\boldsymbol{\Sigma }^{(i+1)}}^{-1}\delta{\boldsymbol{\mu }} = -{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi (\textbf{x})}{\partial \textbf{x}^T}\right ], \end{equation}

(59c) \begin{equation} {\boldsymbol{\mu }}^{(i+1)} ={\boldsymbol{\mu }}^{(i)} + \delta{\boldsymbol{\mu }} \end{equation}

give the updates from $q^{(i)} ={\mathcal{N}}({\boldsymbol{\mu }}^{(i)},\boldsymbol{\Sigma }^{(i)})$ to $q^{(i+1)} ={\mathcal{N}}({\boldsymbol{\mu }}^{(i+1)},\boldsymbol{\Sigma }^{(i+1)})$ and these are exactly the same as those used in the iterative Gaussian variational inference approach presented by [Reference Barfoot, Forbes and Yoon14]. We have arrived at the same variational updates but have done so from the framework of a Bayesian Hilbert space, where it becomes abundantly clear that the minimization algorithm is in fact a slightly simplified version of Newton’s method. This also provides the connection back to the classic Gaussian filtering and smoothing algorithms as discussed by [Reference Barfoot, Forbes and Yoon14].

5.2. Sparsity in Gaussian inference

The effect of sparsity as it applies to iterative Gaussian inference is of particular interest. Let us consider the decomposition of a posterior $p$ in accordance to the general sparsity discussion in Section 4.6, that is,

(60) \begin{equation} p = \,\downarrow \!{\exp (\!-\phi (\textbf{x}))} = \,\downarrow \!{\exp \left (\!-\sum _{k=0}^K \phi _k( \textbf{x}_k )\right )} = \,\downarrow \!{\bigoplus_{k=0}^K} \exp \left (\!-\phi _k( \textbf{x}_k )\right ) = \,\downarrow \!{\bigoplus_{k=0}^K} p_k, \end{equation}

where $\phi _k(\textbf{x}_k)$ is the $k$ th (negative log) factor expression and $\textbf{x}_k = \textbf{P}_k\textbf{x}$ .

As in (D14), we may express the variational estimate as

(61) \begin{align} q^{(i+1)} = \underset{({\mathcal{G}},q^{(i)})}{\,\downarrow \!{}\text{proj}} \, p = \,\downarrow \!{\exp }\biggl ( -(\textbf{x} -{\boldsymbol{\mu }}^{(i)})^T\mathbb{E}_{q^{(i)}} \left [ \frac{\partial \phi (\textbf{x})}{\partial \textbf{x}^T}\right ]- \frac{1}{2} (\textbf{x} -{\boldsymbol{\mu }}^{(i)})^T \, \mathbb{E}_{q^{(i)}} \left [ \frac{\partial ^2 \phi (\textbf{x})}{\partial \textbf{x}^T \partial \textbf{x}} \right ] (\textbf{x} -{\boldsymbol{\mu }}^{(i)})\biggr ), \end{align}

using the measure $q^{(i)} = \mathcal{N}({\boldsymbol{\mu }}^{(i)}, \boldsymbol{\Sigma }^{(i)})$ . To take advantage of sparsity, we need to have it reflected in the expectations herein. The first one leads to

(62) \begin{align} {\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi (\textbf{x})}{\partial \textbf{x}^T}\right ] &= \sum _{k=0}^K{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi _k(\textbf{x}_k)}{\partial \textbf{x}^T}\right ] = \sum _{k=0}^K \left (\frac{\partial \textbf{x}_k}{\partial \textbf{x}}\right )^T{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi _k(\textbf{x}_k)}{\partial \textbf{x}_k^T}\right ]\nonumber \\[5pt] &= \sum _{k=0}^K \textbf{P}_k^T{\mathbb{E}}_{q^{(i)}}\left [\frac{\partial \phi _k(\textbf{x}_k)}{\partial \textbf{x}_k^T}\right ] = \sum _{k=0}^K \textbf{P}_k^T{\mathbb{E}}_{q^{(i)}_k}\left [\frac{\partial \phi _k(\textbf{x}_k)}{\partial \textbf{x}_k^T}\right ], \end{align}

given that $\textbf{x}_k = \textbf{P}_k\textbf{x}$ . For each factor $k$ , then, we are able to shift the differentiation from $\textbf{x}$ to $\textbf{x}_k$ . We draw attention to the last equality, where the expectation simplifies to using $q^{(i)}_k = q^{(i)}_k(\textbf{x}_k)$ , the marginal of the measure for just the variables in factor $k$ . In a similar fashion,

(63) \begin{equation} \mathbb{E}_{q^{(i)}}\left [ \frac{\partial ^2 \phi (\textbf{x})}{\partial \textbf{x}^T \partial \textbf{x}} \right ] = \sum _{k=0}^K \textbf{P}_k^T \,\mathbb{E}_{q_k^{(i)}} \left [ \frac{\partial ^2 \phi _k( \textbf{x}_k )}{\partial \textbf{x}_k^T \partial \textbf{x}_k} \right ] \textbf{P}_k \end{equation}

accounts for the second expectation in (61).

The implication of the factorization is that each factor, identified by $\phi _k(\textbf{x}_k)$ , is projected onto ${\mathcal{G}}_k$ , the marginal subspace associated with variables $\textbf{x}_k$ . The results can then be recombined for the full variational estimate as

(64) \begin{equation} q^{(i+1)} = \underset{({\mathcal{G}},q^{(i)})}{\,\downarrow \!{}\text{proj}} \, p = \,\downarrow \!{}\bigoplus _{k=0}^K \underset{({\mathcal{G}}_k,q_k^{(i)})}{\text{proj}} \, p_k = \,\downarrow \!{} \bigoplus _{k=0}^K q_k^{(i+1)}. \end{equation}

The individual projections of $p_k = \,\downarrow \!{\exp (\!-\phi _k(\textbf{x}_k) )}$ onto $({\mathcal{G}}_k,q^{(i)})$ are

(65) \begin{align} q^{(i+1)}_k & = \underset{({\mathcal{G}}_k,q_k^{(i)})}{\,\downarrow \!{}\text{proj}} \, p_k = \,\downarrow \!{\exp }\left ( -\frac{1}{2} (\textbf{x}_k -{\boldsymbol{\mu }}_k^{(i+1)} )^T \boldsymbol{\Sigma }_k^{(i+1)^{-1}} (\textbf{x}_k -{\boldsymbol{\mu }}_k^{(i+1)} ) \right ) \\[5pt] & = \,\downarrow \!{\exp }\left ( -\frac{1}{2} (\textbf{x} - \textbf{P}_k^T{\boldsymbol{\mu }}_k^{(i+1)} )^T \left ( \textbf{P}_k^T \boldsymbol{\Sigma }_k^{(i+1)^{-1}} \textbf{P}_k \right )(\textbf{x} - \textbf{P}_k^T{\boldsymbol{\mu }}_k^{(i+1)} ) \right ), \end{align}

where

(66) \begin{equation}{\boldsymbol{\mu }}_k^{(i+1)} ={\boldsymbol{\mu }}_k^{(i)} - \boldsymbol{\Sigma }_k^{(i+1)}{\mathbb{E}}_{q_k^{(i)}}\left [\frac{\partial \phi _k(\textbf{x}_k)}{\partial \textbf{x}_k^T}\right ], \quad \boldsymbol{\Sigma }_k^{(i+1)^{-1}} ={\mathbb{E}}_{q_k^{(i)}}\left [\frac{\partial ^2 \phi _k( \textbf{x}_k )}{\partial \textbf{x}_k^T \partial \textbf{x}_k}\right ]. \end{equation}

It is straightforward to show that the vector sum of $q_k$ from (65) reproduces (61). (Note that $\textbf{P}_k\textbf{P}_k^T = \textbf{1}$ as $\textbf{P}_k$ is a projection matrix and $\textbf{P}_k^T$ the corresponding dilation.)

As explained in detail by [Reference Barfoot, Forbes and Yoon14], it would be too expensive for practical problems to construct first $\boldsymbol{\Sigma }^{(i)}$ and then extract the required blocks for the marginals, $q_k^{(i)} = \mathcal{N}({\boldsymbol{\mu }}_k^{(i)}, \boldsymbol{\Sigma }_k^{(i)}) = \mathcal{N}(\textbf{P}_k{\boldsymbol{\mu }}^{(i)}, \textbf{P}_k\boldsymbol{\Sigma }^{(i)}\textbf{P}_k^T)$ . We see from the above development that we actually only require the blocks of $\boldsymbol{\Sigma }^{(i)}$ corresponding to the nonzero blocks of its inverse and the method of [Reference Takahashi, Fagan and Chen41] can be used to extract the required blocks efficiently. In Ref. [Reference Barfoot, Forbes and Yoon14], the authors provide numerical experiments showing the efficacy of this approach.

5.3. Example: SLAM

A main purpose in the current paper was to show the connection between Bayes space [Reference van den Boogaart, Egozcue and Pawlowsky-Glahn43] and Gaussian variational inference [Reference Barfoot, Forbes and Yoon14]. We see that minimizing the KL divergence between a true Bayesian posterior and an approximation can be viewed as iterative projection in Bayes space. Moreover, by exploiting the sparsity that Bayes space makes clear, the method has the potential to be applied to quite large inference problems.

Following [Reference Barfoot, Forbes and Yoon14], we consider a batch SLAM problem with a robot driving around and building a map of landmarks as depicted in Fig. 7. The robot is equipped with a laser rangefinder and wheel odometers and must estimate its own trajectory and the locations of a number of tubular landmarks. This dataset has been used previously by [Reference Barfoot, Tong and Sarkka15] to test SLAM algorithms. Groundtruth for both the robot trajectory and landmark positions (this is a unique aspect of this dataset) is provided by a Vicon motion capture system. The whole dataset is 12,000 timesteps long (approximately 20 min of real time). It was assumed that the data association (i.e., which measurement corresponds to which landmark) is known in this experiment to restrict testing to the state estimation part of the problem.

Figure 7. A special case of the approach presented in this paper was previously demonstrated [Reference Barfoot, Forbes and Yoon14] to be a useful and practical tool for robotic state estimation. The method, called exactly sparse Gaussian variational inference (ESGVI), was used to solve the simultaneous localization and mapping (SLAM) problem, outperforming the standard maximum a posteriori (MAP) estimation in certain cases. The current paper reinterprets this earlier work in a new mathematical formalism. Figure reproduced from [Reference Barfoot, Forbes and Yoon14].

The state to be estimated is

(67) \begin{equation} \textbf{x} = \begin{bmatrix} \textbf{x}_0 \\[5pt] \textbf{x}_1 \\[5pt] \vdots \\[5pt] \textbf{x}_K \\[5pt] \textbf{m}_1 \\[5pt] \vdots \\[5pt] \textbf{m}_L \end{bmatrix}, \quad \textbf{x}_k = \begin{bmatrix} x_k \\[5pt] y_k \\[5pt] \theta _k \\[5pt] \dot{x}_k \\[5pt] \dot{y}_k \\[5pt] \dot{\theta }_k \end{bmatrix}, \quad \textbf{m}_\ell = \begin{bmatrix} x_\ell \\[5pt] y_\ell \end{bmatrix}, \end{equation}

where $\textbf{x}_k$ is a robot state (position, orientation, velocity, angular velocity) and $\textbf{m}_\ell$ a landmark position.

For the (linear) prior factor on the robot states, we have

(68) \begin{equation} \varphi _k = \left \{ \begin{array}{ll} \dfrac{1}{2} (\textbf{x}_0 - \check{\textbf{x}}_0)^T \check{\textbf{P}}^{-1} (\textbf{x}_0 - \check{\textbf{x}}_0) & \quad k=0 \\[12pt] \dfrac{1}{2} (\textbf{x}_k - \textbf{A} \textbf{x}_{k-1})^T \textbf{Q}^{-1} (\textbf{x}_k - \textbf{A} \textbf{x}_{k-1}) \;\;& \quad k\gt 0 \end{array} \right ., \end{equation}

with

(69) \begin{align} \check{\textbf{P}} = \text{diag}(\sigma _x^2, \sigma _y^2,\sigma _\theta ^2, \sigma _{\dot{x}}^2, \sigma _{\dot{y}}^2, \sigma _{\dot{\theta }}^2), \;\; \textbf{A} = \begin{bmatrix} \textbf{1}\;\;\;\; & T \textbf{1} \\[12pt] \textbf{0}\;\;\;\; & \textbf{1} \end{bmatrix}, \nonumber \\[12pt] \textbf{Q} = \begin{bmatrix} \dfrac{1}{3} T^3 \textbf{Q}_C\;\;\;\; & \dfrac{1}{2} T^2 \textbf{Q}_C \\[12pt] \dfrac{1}{2} T^2 \textbf{Q}_C\;\;\;\; & T \textbf{Q}_C \end{bmatrix}, \quad \textbf{Q}_C = \text{diag}(Q_{C,1}, Q_{C,2},Q_{C,3}), \end{align}

where $T$ is the discrete-time sampling period, $Q_{C,i}$ are power spectral densities, and $\sigma _x^2$ , $\sigma _y^2$ , $\sigma _\theta ^2$ , $\sigma _{\dot{x}}^2$ , $\sigma _{\dot{y}}^2$ , $\sigma _{\dot{\theta }}^2$ are variances. The robot state prior encourages constant velocity [[Reference Barfoot10], §3, p.85].

The (nonlinear) odometry factors, derived from wheel encoder measurements, are

(70) \begin{equation} \psi _k = \frac{1}{2} \left ( \textbf{v}_k - \textbf{C}_k \textbf{x}_k \right )^T \textbf{S}^{-1} \left ( \textbf{v}_k - \textbf{C}_k \textbf{x}_k \right ), \end{equation}

where

(71) \begin{gather} \textbf{v}_k = \begin{bmatrix} u_k \\[5pt] v_k \\[5pt] \omega _k \end{bmatrix}, \quad \textbf{C}_k = \begin{bmatrix} 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & \cos \theta _k \;\;\;\;\; & \sin \theta _k \;\;\;\;\; & 0 \\[5pt] 0\;\;\;\;\; & 0\;\;\;\;\; & 0\;\;\;\;\; & -\sin \theta _k \;\;\;\;\; & \cos \theta _k \;\;\;\;\; & 0 \\[5pt] 0 \;\;\;\;\; & 0 \;\;\;\;\; & 0 \;\;\;\;\; & 0 \;\;\;\;\; & 0 \;\;\;\;\; & 1 \end{bmatrix}, \quad \textbf{S} = \text{diag}\left (\sigma ^2_u, \sigma _v^2, \sigma _\omega ^2 \right ). \end{gather}

The $\textbf{v}_k$ consists of measured forward, lateral, and rotational speeds in the robot frame, derived from wheel encoders; we set $v_k = 0$ , which enforces the nonholonomy of the wheels as a soft constraint. The $\sigma _u^2$ , $\sigma _v^2$ , and $\sigma _\omega ^2$ are measurement noise variances.

The (nonlinear) bearing measurement factors, derived from a laser rangefinder, are

(72) \begin{equation} \psi _{\ell,k} = \frac{1}{2} \frac{\left ( \beta _{\ell,k} - g(\textbf{m}_{\ell },\textbf{x}_k) \right )^2}{\sigma _r^2}, \end{equation}

with

(73) \begin{equation} g(\textbf{m}_{\ell },\textbf{x}_k) = \text{tan}^{-1}(y_\ell - y_k - d \sin \theta _k, x_\ell - x_k - d \cos \theta _k) - \theta _k, \end{equation}

where $\beta _{\ell,k}$ is a bearing measurement from the $k$ th robot pose to the $\ell$ th landmark, $d$ is the offset of the laser rangefinder from the robot center in the longitudinal direction, and $\sigma _r^2$ is the measurement noise variance. Although the dataset provides range to the landmarks as well, we chose to neglect these measurements to make the problem difficult and therefore show the benefit of taking a variational inference approach. Our setup is similar to a monocular camera situation, which is known to be a challenging SLAM problem.

Putting all the factors together, the posterior that we would like to project to a Gaussian is

(74) \begin{equation} p(\textbf{x}) = \,\downarrow \!{\exp }(\!-\phi (\textbf{x})), \quad \phi (\textbf{x}) = \sum _{k=0}^K \varphi _k + \sum _{k=0}^K \psi _k + \sum _{k=1}^K \sum _{\ell =1}^L \psi _{\ell,k} + \text{constant}, \end{equation}

where it is understood that not all $L = 17$ landmarks are actually seen at each timestep and thus we must remove the factors for unseen landmarks.

We then carry out iterative projection to a multivariate Gaussian estimate, $q^{(i)} ={\mathcal{N}}({\boldsymbol{\mu }}^{(i)},\boldsymbol{\Sigma }^{(i)})$ , according to

(75) \begin{equation} q^{(i+1)} = \underset{({\mathcal{G}},q^{(i)})}{\,\downarrow \!{}\text{proj}} \, p, \end{equation}

where we use (59) at implementation.

In a sub-sequence of the full dataset comprising $2000$ timestamps with $(x,y,\theta,\dot{x},\dot{y},\dot{\theta })$ for the robot state at each time and $17$ landmarks with $(x,y)$ position, the dimension of the state estimation problem is $N = 12034$ . Spanning the indefinite-Gaussian subspace requires $N$ basis functions for the mean and $\frac{1}{2} N(N+3)$ basis functions for the covariance, or $72,438,663$ basis functions total. To naively apply the idea of iterative projection to Bayes space would be intractable. However, by exploiting the sparsity that Bayes space affords (see Sections 4.6 and 5.2) we are able to do this in a very computationally efficient way. In Ref. [Reference Barfoot, Forbes and Yoon14], the authors provide further implementation details of this experiment.

Figure 8 shows estimation error plots for the section of the robot’s trajectory in Fig. 7. Not only does this show the concepts of Bayes space can be applied to a large problem ( $N$ in the range of thousands) but also that there are some situations where it performs slightly better than the standard MAP Gauss-Newton (GN) algorithm.

Figure 8. Error plots for a portion of the trajectory in the SLAM problem conducted by [Reference Barfoot, Forbes and Yoon14] and discussed in Section 5.3. The exactly sparse Gaussian variational inference (ESGVI) algorithm (red) is equivalent to the iterative projection approach described herein. The maximum a posteriori (MAP) Gauss-Newton (GN) algorithm (blue) is the more standard approach to solving this type of problem. Here, we see ESGVI performing slightly better than MAP GN in terms of smaller errors and more consistency (i.e., errors staying within covariance envelope). Note, in the heading error plot, the red mean line is hidden behind the blue one.