Gaussian belief propagation
Least squares at maximal elaboration
November 25, 2014 — September 1, 2022
\[\renewcommand{\argmin}{\operatorname{arg min}} \renewcommand{\KL}{\operatorname{KL}} \renewcommand{\var}{\operatorname{Var}} \renewcommand{\cov}{\operatorname{Cov}} \renewcommand{\dd}{\mathrm{d}} \renewcommand{\vv}[1]{\boldsymbol{#1}} \renewcommand{\rv}[1]{\mathsf{#1}} \renewcommand{\vrv}[1]{\vv{\rv{#1}}} \renewcommand{\disteq}{\stackrel{d}{=}} \renewcommand{\gvn}{\mid} \renewcommand{\Ex}{\mathbb{E}} \renewcommand{\Pr}{\mathbb{P}} \renewcommand{\one}{\unicode{x1D7D9}}\]
A particularly tractable model assumption for message-passing inference generalizes classic Gaussian Kalman filters and Bayesian linear regression with Gaussian prior, or, in a frequentist setting, least squares regression. Essentially, we regard the various nodes in our system as jointly Gaussian RVs with given prior mean and covariance (i.e. we do not allow the variances themselves to be random, as a Gaussian is not a valid prior for a variance). This works well because the Gaussian distribution is a natural exponential family, and moreover, self-conjugate under linear transformations, and thus the message-passing updates are simple linear algebra operations.
How does this relate to least squares? Bickson (2009) puts it elegantly: \[\begin{array}{c} \mathbf{A} \mathbf{x}=\mathbf{b}\\ \mathbb{\Downarrow}\\ \mathbf{A x}-\mathbf{b}=0\\ \mathbb{\Downarrow}\\ \min _{\mathbf{x}}\left(\frac{1}{2} \mathbf{x}^{T} \mathbf{A} \mathbf{x}-\mathbf{b}^{T} \mathbf{x}\right)\\ \mathbb{\Downarrow}\\ \max _{\mathbf{x}}\left(-\frac{1}{2} \mathbf{x}^{T} \mathbf{A} \mathbf{x}+\mathbf{b}^{T} \mathbf{x}\right)\\ \mathbb{\Downarrow}\\ \max _{\mathbf{x}} \exp \left(-\frac{1}{2} \mathbf{x}^{T} \mathbf{A} \mathbf{x}+\mathbf{b}^{T} \mathbf{x}\right) \end{array}\]
GBP comes with various conveniences. For example, we can implement it mechanically without knowing the structure of the graphical model in advance. We can execute it concurrently since the model is robust against random execution schedules. This seems to be a happy coincidence of several features of multivariate Gaussians, e.g. that the Gaussian is self-conjugate and thus passes through the factor and variable nodes in commensurable form. Also, all the required operations can be represented by simple linear algebra.
NB: notation is varied and thus possibly confusing. In the original paper (Weiss and Freeman 2001) variables are presented as pairwise linear models rather than a classic graphical model, which is a good insight but possibly confusing.
A famous interactive tutorial is gaussianbp.github.io, reprising Davison and Ortiz (2019). Loeliger et al. (2007) is also nice; it includes explicit accounting for Fourier transforms and arbitrary matrix products.
1 Parameterization
Gaussian BP is often conducted using the canonical parameterization, in which a particular node prior \(m\) may be written as: \[ p_{m}\left(\mathbf{x}_{m}\right)\propto e^{-\frac{1}{2}\left[\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{m}\right)^{\top} \boldsymbol{\Lambda}_{m}\left(\mathbf{x}_{m}-\boldsymbol{\mu}_{m}\right)\right]} \] where \(\boldsymbol{\mu}_{m}\) is the mean of the distribution and \(\Lambda_{m}\) is its precision (i.e. inverse covariance.)
Davison and Ortiz (2019) recommend the “information parameterization” of Eustice, Singh, and Leonard (2006): \[ p_{m}\left(\mathbf{x}_{m}\right)\propto e^{\left[-\frac{1}{2} \mathbf{x}_{m}^{\top} \Lambda_{m} \mathbf{x}_{m}+\boldsymbol{\eta}_{m}^{\top} \mathbf{x}_{m}\right]} \] the value of either. The information vector \(\boldsymbol{\eta}_{m}\) relates to the mean vector as \[ \boldsymbol{\eta}_{m}=\Lambda_{m} \boldsymbol{\mu}_{m}. \] The information form is convenient as multiplication of distributions (when we update) is handled simply by adding the information vectors and precision matrices.
Let us recap that in Gaussian tricks.
2 Linearization
See GP noise transforms for a typical way to do it. In the information parameterization the rules look different. Suppose the expected value of the measurement is given by non-linear map \(h,\) in the sense that we observe Gaussian noise centred on \(h(\mathbf{x}_s)\). We define the associated Gaussian factor as follows: \[ f_{s}\left(\mathbf{x}_{s}\right)\propto e^{-\frac{1}{2}\left[\left(\mathbf{z}_{s}-\mathbf{h}_{s}\left(\mathbf{x}_{s}\right)\right)^{\top} \Lambda_{s}\left(\mathbf{z}_{s}-\mathbf{h}_{s}\left(\mathbf{x}_{s}\right)\right)\right]} \] We apply the first-order Taylor series expansion of non-linear measurement function \(\mathbf{h}_{s}\) to find its approximate value for state values \(\mathbf{x}_{s}\) close to \(\mathbf{x}_{0}\) : \[ \mathbf{h}_{s}\left(\mathbf{x}_{s}\right) \approx \mathbf{h}_{s}\left(\mathbf{x}_{0}\right)+\mathrm{J}_{s}\left(\mathbf{x}_{s}-\mathbf{x}_{0}\right). \] Here \(\mathrm{J}_{s}\) is the Jacobian \(\left.\frac{\partial \mathbf{h}_{s}}{\partial \mathbf{x}_{s}}\right|_{\mathbf{x}_{s}=\mathbf{x}_{0}} .\) With a bit of rearranging quadratic forms, we can massage \(\mathbf{h}_{s}\left(\mathbf{x}_{s}\right)\) around state variables \(\mathbf{x}_{0}\), turning it into a Gaussian factor expressed in terms of \(\mathbf{x}_{s}\). We use the linear factor represented in information form by information vector \(\boldsymbol{\eta}_{s}\) and precision matrix \(\Lambda_{s}^{\prime}\) calculated as follows: \[ \begin{aligned} \boldsymbol{\eta}_{s} &=\mathrm{J}_{s}^{\top} \Lambda_{s}\left[\mathrm{~J}_{s} \mathbf{x}_{0}+\mathbf{z}_{s}-\mathbf{h}_{s}\left(\mathbf{x}_{0}\right)\right] \\ \Lambda_{s}^{\prime} &=\mathrm{J}_{s}^{\top} \Lambda_{s} \mathrm{~J}_{s}. \end{aligned} \]
3 Variational inference
Intuition building. Noting that the distance between two Gaussians in KL is given by \[ \KL(\mathcal{N}(\mu_1,\Sigma_1)\parallel \mathcal{N}(\mu_2,\Sigma_2)) ={\frac {1}{2}}\left(\operatorname {tr} \left(\Sigma _{2}^{-1}\Sigma _{1}\right)+(\mu_{2}-\mu_{1})^{\mathsf {T}}\Sigma _{2}^{-1}(\mu_{2}-\mu_{1})-k+\ln \left({\frac {\det \Sigma _{2}}{\det \Sigma _{1}}}\right)\right)\] We can think about belief propagation in terms of KL-similar updates.
Suppose I have a Gaussian RV with some prior distribution \(\vrv{x}\sim \mathcal{N}(\vv{m}_{\vrv{x}}, \Sigma_{\vrv{x}})\) and I observe some datum \(\vrv{y}\sim f(\vrv{x} + \epsilon)\) where \(\epsilon \sim \mathcal{N}(0, \Sigma_{\vrv{y}})\). Define \(\hat{\vrv{x}}:=\vrv{x}\gvn \vrv{y}.\) Assert that we can approximate \(\vrv{x}\sim \mathcal{N}(\vv{m}_{\hat{\vrv{x}}}, \Sigma_{\hat{\vrv{x}}})\). What is the best approximation to the posterior distribution of \(\vrv{x}\) given the observation? I.e. how do we choose \(\vv{m}_{\hat{\vrv{x}}}, \Sigma_{\hat{\vrv{x}}}\)? One way is to choose them to minimise a KL divergence, \[\begin{aligned} \vv{m}_{\hat{\vrv{x}}}, \Sigma_{\hat{\vrv{x}}} &=\argmin_{\vv{m}, \Sigma} \KL\left(\mathcal{N}(\vv{m}, \Sigma)\parallel\vrv{x}\gvn \vrv{y}\right)\\ &=\argmin_{\vv{m}, \Sigma} \KL\left(\mathcal{N}(\vv{m}, \Sigma)\parallel \mathcal{N}(\vv{m}_{\vrv{x}}, \Sigma_{\vrv{x}}\right). \end{aligned}\] We still need to choose a way of estimating this posterior which could exploit stochastic gradient estimation or a GP noise transforms. TBC.
4 Parallelisation
El-Kurdi’s thesis (Y. M. El-Kurdi 2014) exploits some assortment of tricks to attain high-speed convergence in pde-type models.
5 Robust
There are generalised belief propagations based on elliptical laws (Agarwal et al. 2013; Davison and Ortiz 2019), making this into a kind of elliptical belief propagation. It is usually referred to as a robust factor or as dynamic covariance scaling. Usually this is done with robust losses; I’m not sure why other elliptical distributions are unpopular.
6 Use in PDEs
Finite Element Models can be expressed as locally-linear constraints and thus expressed using GBP (Y. El-Kurdi et al. 2016; Y. M. El-Kurdi 2014; Y. El-Kurdi et al. 2015).
7 Tools
7.1 jaxfg
7.2 gradslam
Differentiable GPB solver (Jatavallabhula, Iyer, and Paull 2020).
7.3 ceres solver
ceres-solver, (C++), the google least squares solver, seems to solve this kind of problem. I am not sure where the covariance matrices go in. I occasionally see mention of “CUDA” in the source repo so maybe it exploits GPUs these days.