Approximations of the measurement update and model likelihood for nonlinear spatiotemporal Cox processes

Epilogue: These notes were part of a project to infer unobserved ("latent") states in a neural field model, as well as the parameters of that model, from spiking observations. It has since been published. Ultimately, for speed we ended up selecting the Laplace approximation for the measurement update, solved via Newton-Raphson.

[get PDF]

We are interested in approximations to the measurement update for a spatially-extended latent-variable point-process, where the latent variables are also spatial fields that undergo dynamics similar to chemical reaction-diffusion systems. 

The latent variables or fields are concentrations, activations, or some similar physical quantity. They are therefore constrained to be non-negative, and also typically must obey conservation laws. Additionally, the observed point-process intensity field must also be constrained to be non-negative. 

Such systems arise in chemical reaction diffusion systems, epidemiological models, and neural field models, where the measurement is a point-process that is coupled indirectly to the latent spatiotemporal system.

Problem statement

We will use a multivariate Gaussian approximation to model the joint distribution of latent variables. In the continuous filed case, this is a Gaussian process. In the numerical implementation, we project this process onto a finite basis to give a finite-dimensional multivariate Gaussian distribution. The filtering update therefore requires finding a multivariate Gaussian approximation to the non-conjugate update of a Poisson observation and a multivariate Gaussian prior. 

Estimating this posterior involves an integral that becomes intractable in high-dimensions. There are numerous approximation methods to handle this, including moment-matching, variational Bayes, expected log-likelihoods, and the Laplace approximation. Often, the choice of link function furthers constrain which methods are applicable, as efficient algorithms may only be available in some classes. The update must be also constrained to prevent negative activity in the estimated latent states.

Multivariate Gaussian model

Let $\mathbf{A}(\mathbf{x})$ denote a latent vector of activity, defined over a spatial domain with coordinates $\mathbf{x}\in \mathrm{\Omega }$ on some (likely bounded) domain $\mathrm{\Omega }$.  We approximate the prior distribution over the latent activity vector $\mathbf{A}(\mathbf{x})$ in terms of its first two moments. We interpret these moments to reflect the mean and variance of a multivariate Gaussian distribution, for the purposes of both moment closure and the measurement update. This assumption is an approximation, as in general the tails of the Gaussian that extent to negative activation values are unphysical. 

$$\begin{array}{}\text{(1)}& Pr(\mathbf{A}(\mathbf{x}))\sim \mathrm{Gaussian}({\mu }_A(\mathbf{x}),{\mathrm{\Sigma }}_A(\mathbf{x}))\end{array}$$

$\mathbf{A}(x)$ is a field defined over a continuous spatial region, and so Eq. $\text{1}$ denotes a Gaussian process. In practice, we project this continuous process onto a finite basis and work with $N$ discrete activation variables corresponding to spatial regions, i.e. $\mathbf{A}=\{A_1,..,A_N\}$. In this numerical implementation, Eq. $\text{1}$ represents a finite-dimensional multivariate Gaussian distribution.

Latent field definition and discretization

This activation field is mapped to point-process intensities via a link function ${\lambda }_0=f(\mathbf{A})$. Furthermore, our observation model may include heterogeneity in terms of the density of agents or background level of activity, and we therefore incorporate a spatially inhomogeneous gain $\gamma (\mathbf{x})$ and bias $\beta (\mathbf{x})$ that adjust the observe point-process intensity. The intensity as a function of spatial coordinates $\mathbf{x}$ is then:

$$\begin{array}{}\text{(1)}& \lambda (\mathbf{x})=f[\mathbf{A}(\mathbf{x})]\cdot \gamma (\mathbf{x})+\beta (\mathbf{x}).\end{array}$$

In practice, we project this continuous, infinite-dimensional process onto a finite set of discrete basis elements $\mathbf{B}=\{b_1,..,b_N\}$, where the expected firing rate for the $n^{th}$ basis element is:

$$\begin{array}{}\text{(2)}& {\lambda }_n={\int }_{\mathbf{x}\in \mathrm{\Omega }}{b}_n(\mathbf{x})\cdot \lambda (\mathbf{x}),\end{array}$$

where ${\int }_{\mathbf{x}\in \mathrm{\Omega }}$ denotes integration over the spatial domain $\mathrm{\Omega }$ parameterized by $\mathbf{x}$. If the variations in the activity, gain, and bias, are small relative to the scale of the basis elements, we may approximate this integral as:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\lambda }_n& \approx \lambda \left({⟨\mathbf{x}⟩}_{b_n}\right)\cdot v_n\\ v_n& ={\int }_{\mathbf{x}\in \mathrm{\Omega }}{b}_n(\mathbf{x}),\end{array}\end{array}$$

where ${⟨\mathbf{x}⟩}_{b_n}$ is the centre of mass of the basis element $b_n$, and $v_n$ is the volume of said basis element. We consider an especially simple case where basis functions all have identical volume $v$, so we may write the regional intensity as

$$\begin{array}{}\text{(4)}& {\lambda }_n\approx v\cdot [{\gamma }_{n}f(A_n)+{\beta }_n],\end{array}$$

where $v$ is the (uniform) volume of the basis elements, e.g. $v={\mathrm{\Delta }}_x^2{\mathrm{\Delta }}_t$ for a process with two spatial dimensions and one time dimensions, with a fixed region size. The spatially-varying intensity is then represented as a vector of per-region intensities, $\lambda =\{{\lambda }_1,..,{\lambda }_N\}$. Since the volume parameter $v$ is redundant to the gain parameter $\gamma$, in the derivations that follow we assume that the bias and gain parameters have been premultiplied by the volume. 

We introduce the gain and bias parameters to decouple inhomogeneity in the measurements from the underlying, latent spatiotemporal process. For example, in the case of retinal waves, different regions have differing densities of retinal ganglion cells, which amounts to a spatially inhomogeneous gain. Additionally, the amount of spontaneously-active background activity varies. In order to build up mathematical solutions that are immediately useful for numerical implementation, we carry-through these bias and gain parameters in the derivations that follow. 

Count observations and the measurement posterior

Given the Gaussian prior, the Posterior estimate of the latent activations $\mathbf{A}$ is given by Bayes' rule:

$$\begin{array}{}\text{(5)}& Pr\left(\mathbf{A}\mid \mathbf{Y}\right)=\frac{Pr\left(\mathbf{A}\right)}{Pr\left(\mathbf{Y}\right)}\cdot Pr\left(\mathbf{Y}\mid \mathbf{A}\right)\end{array}$$

We observe event counts over a finite number of spatial regions, and these observed counts are independent conditioned on the per-region intensity, so we can write:

$$\begin{array}{}\text{(6)}& Pr\left(\mathbf{A}\mid \mathbf{Y}\right)=\frac{Pr\left(\mathbf{A}\right)}{Pr\left(\mathbf{Y}\right)}\prod _{n\in 1..N}Pr\left(y_n\mid \mathbf{A}\right)\end{array}$$

The dependence of the counts on the latent activation can be expanded to include the regional intensity ${\lambda }_n$ as:

$$\begin{array}{}\text{(7)}& Pr\left(y_n\mid A_n\right)=Pr\left(y_n\mid {\lambda }_n\right)\cdot Pr({\lambda }_n\mid \mathbf{A})\end{array}$$

Since the dependence of ${\lambda }_n$ on $\mathbf{A}$ is deterministic (and vice-versa), we can therefore write

$$\begin{array}{}\text{(8)}& Pr\left(y_n\mid A_n\right)=\frac{Pr\left(\mathbf{A}\right)}{Pr\left(\mathbf{Y}\right)}\prod _{n\in 1..N}Pr\left(y_n\mid {\lambda }_n\right)\end{array}$$

We observe regional counts $\mathbf{Y}=\{y_1,..,y_N\}$, which are Poisson distributed, with the observation likelihood

$$\begin{array}{}\text{(9)}& Pr(y_n\mid {\lambda }_n)=\frac{{\lambda }_n^{y_n}}{y_n!}{e}^{-{\lambda }_n}.\end{array}$$

The log-posterior

Consider the logarithmic form of the measurement update:

$$\begin{array}{}\text{(10)}& \log Pr\left(\mathbf{A}\mid \mathbf{Y}\right)=\log Pr\left(\mathbf{A}\right)-\log Pr\left(\mathbf{Y}\right)+\sum _{n\in 1..N}\log Pr\left(y_n\mid {\lambda }_n\right)\end{array}$$

The prior $\mathbf{A}$ is approximated by a Gaussian distribution $\mathcal{N}({\mu }_A,{\mathrm{\Sigma }}_A)$, and so the log-prior on $\mathbf{A}$ is:

$$\begin{array}{}\text{(11)}& \log Pr(\mathbf{A})=-\frac{1}{2}[\log |2\pi {\mathrm{\Sigma }}_A|+(\mathbf{A}-{\mu }_A{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}(\mathbf{A}-{\mu }_A)].\end{array}$$

The conditional log-likelihood is given by the Poisson observation model with regional intensity ${\lambda }_n={\gamma }_{n}f(A_n)+{\beta }_n$:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\log Pr(\mathbf{Y}\mid \mathbf{A})& =\sum _{n\in 1..N}\log Pr\left(y_n\mid {\lambda }_n\right)\\ & =\sum _{n\in 1..N}[y_n\log ({\lambda }_n)-{\lambda }_n].\\ & =\sum _{n\in 1..N}[y_n\log ({\gamma }_{n}f(A_n)+{\beta }_n)-({\gamma }_{n}f(A_n)+{\beta }_n)].\end{array}\end{array}$$

The marginal log-likelihood of the count observations $\log Pr(\mathbf{Y})$ cannot be computed, except via an intractable integral. Approximating this integral will be a major challenge for computing the model likelihood, which we will address later. However, for fixed count observations $\mathbf{Y}$, the $Pr\mathbf{Y}$ term is constant, and so:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\log Pr\left(\mathbf{A}\mid \mathbf{Y}\right)=& -\frac{1}{2}[\log |2\pi {\mathrm{\Sigma }}_A|+(\mathbf{A}-{\mu }_a{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}(\mathbf{A}-{\mu }_a)]\\ & +\sum _{n\in 1..N}[y_n\log ({\lambda }_n)-{\lambda }_n]+\text{constant}\end{array}\end{array}$$ 

Approximation methods

In this section, we explore various approaches to obtaining a Gaussian approximation to the posterior $Pr(\mathbf{A}|\mathbf{Y})\approx Q(\mathbf{A})\sim \mathcal{N}({\hat{\mu }}_A,{\hat{\mathrm{\Sigma }}}_A)$. We examine three approaches: the Laplace approximation, the variational Bayes approach, and moment-matching. 

Laplace approximation

For the Laplace approximation, we find the mode of the posterior and interpret the curvature at this mode as the inverse of the covariance matrix. Since we are dealing with spatiotemporal processes driven by physical agents (e.g. molecules, neurons, humans), we constrain the posterior mode to be non-negative. This departs slightly from the traditional Laplace approximation, in which the posterior mode is a non-extremal local maximum with zero slope. For this reason, the interpretation of the curvature at the mode as the inverse of the covariance must be treated with caution. 

$$\begin{array}{}\text{(13)}& \widehat{{\mu }_A}=\underset{\mathbf{A}}{\mathrm{argmax}}[\log Pr(\mathbf{A}\mid \mathbf{Y})]\end{array}$$

This can be solved with gradient descent or the Newton-Raphson method, which requires the gradient and Hessian of the log-posterior with respect to $\mathbf{A}$. 

We introduce some abbreviations to simplify the notation in the derivations that follow. Denote log-probabilities "$\log Pr(\cdot )$'' as ${\mathcal{L}}_{(\cdot )}$, and denote the first and second derivatives of the log-measurement likelihood with respect to individual activation variables $A\in \{A_1,..,A_N\}$ as ${\mathcal{L}}_{y\mid A}^{\prime }$ and ${\mathcal{L}}_{y\mid A}^{″}$, respectively. 

Note that  ${\lambda }_n$ is synonymous with the locally adjusted rate ${\gamma }_{n}f(A_n)+{\beta }_n$, such that ${\lambda }_n^{\prime }={\gamma }_n{f}^{\prime }(A_n)$. We also omit indexing by the basis function number $n$ when unambiguous. With these abbreviations, the gradient and Hessian of the log-posterior in $A$ are:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathrm{\nabla }}_{\mathbf{A}}{\mathcal{L}}_{\mathbf{A}\mid \mathbf{Y}}& =({\mu }_a-\mathbf{A}{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}+{\mathrm{\nabla }}_{\mathbf{A}}{\mathcal{L}}_{\mathbf{Y}\mid \mathbf{A}}\\ {\mathrm{\nabla }}_{\mathbf{A}}^2{\mathcal{L}}_{\mathbf{A}\mid \mathbf{Y}}& ={\mathrm{\Sigma }}_A^{-1}+{\mathrm{\nabla }}_{\mathbf{A}}^2{\mathcal{L}}_{\mathbf{Y}\mid \mathbf{A}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathcal{L}}_{y\mid A}^{\prime }& =\frac{\gamma }{\lambda }{f}^{\prime }(A)(y-\lambda )\\ {\mathcal{L}}_{y\mid A}^{″}& =\frac{\gamma }{\lambda }[(y-\lambda )f^{″}(A)-y\frac{\gamma }{\lambda }{f}^{\prime }(A{)}^2]\end{array}\end{array}$$\subsubsection{The identity link function}

In the case that $f(A_n)=A_n$, these gradients simplify to:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\mathcal{L}}_{y\mid A}^{\prime }& =(y-\lambda )\frac{\gamma }{\lambda }\\ {\mathcal{L}}_{y\mid A}^{″}& =-y{\left(\frac{\gamma }{\lambda }\right)}^2\end{array}\end{array}$$

The exponential link function

In the case that $f(A)=\exp (A)$:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathcal{L}}_{y\mid A}^{\prime }& =(y-\lambda )\left(\frac{\gamma e^A}{\lambda }\right)\\ {\mathcal{L}}_{y\mid A}^{″}& =\left(\frac{\gamma e^A}{\lambda }\right)[y(1-\frac{\gamma e^A}{\lambda })-\lambda ]\end{array}\end{array}$$

For more flexibility, one might add another gain parameter inside the exponentiation, i.e. $f(A_n)=\exp (\delta A_n)$, which gives:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\mathcal{L}}_{y\mid A}^{\prime }& =(y-\lambda )\left(\frac{\gamma \delta }{\lambda }{e}^{\delta A}\right)\\ {\mathcal{L}}_{y\mid A}^{″}& =\gamma {\delta }^2[-\frac{y\gamma }{{\lambda }^2}{e}^{2\delta A}+(\frac{y}{\lambda }-1)e^{\delta A}]\\ & =\delta (y-\lambda )\left(\frac{\gamma \delta }{\lambda }{e}^{\delta A}\right)-y{\left(\frac{\gamma \delta }{\lambda }{e}^{\delta A}\right)}^2\end{array}\end{array}$$

The quadratic link function

Let $\lambda =A^2\gamma +\beta$. That is, $f(A)=A^2$, $f^{\prime }(A)=2A$, and $f^{″}(A)=2$.

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\mathcal{L}}_{y\mid A}^{\prime }& =\frac{\gamma }{2}A(y-\lambda )\\ {\mathcal{L}}_{y\mid A}^{″}& =\frac{\gamma }{\lambda }[(y-\lambda )2-y\frac{\gamma }{\lambda }4A^2]\end{array}\end{array}$$

A more flexible parameterization is $\lambda =\gamma (A+b{)}^2+\beta$ gives:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\mathcal{L}}_{y\mid A}^{\prime }& =\frac{\gamma }{\lambda }(2A+b)(y-\lambda )\\ {\mathcal{L}}_{y\mid A}^{″}& =\frac{\gamma }{\lambda }[(y-\lambda )2-y\frac{\gamma }{\lambda }(2A+b{)}^2]\end{array}\end{array}$$

Logistic link 1

We might want to consider the logistic link function, which maps the range $(-\mathrm{\infty },\mathrm{\infty })$ in activation $A$ to $(0,1)$, which then may be further adjusted to span a given range using the gain/bias parameters:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}f& =\frac{1}{1+e^{-\delta A}}\\ f^{\prime }& =\delta \frac{e^{-\delta A}}{(1+e^{-\delta A}{)}^2}=\delta [1-f(A)]f(A)\\ f^{″}& =\delta [1-2f]f^{\prime }.\end{array}\end{array}$$

Logistic link 2

If the activation is bounded on $[0,\mathrm{\infty })$, it might make more sense to apply the logistic function to the log-activation, yielding the following link function:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}f& =\frac{A}{ϵ+A}\\ f^{\prime }& =\frac{ϵ}{(ϵ+A{)}^2}\\ f^{″}& =-2\frac{ϵ}{(ϵ+A{)}^3}.\end{array}\end{array}$$

where ${ϵ}_n$ is an additional free parameter that acts a like an inverse gain.

Variational approximation

In the variational approximation, we find a Gaussian distribution $Q(\mathbf{A})\sim \mathcal{N}({\hat{\mu }}_A,{\hat{\mathrm{\Sigma }}}_A)$ that approximates the true posterior by minimizing the KL divergence of the true posterior from the approximating distribution $Q$. This is conceptually equivalent to jointly maximizing the entropy of $Q$ while also maximizing the expected log-probability of the true posterior under $Q$. 

$$\begin{array}{}\text{(3)}& \underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmin}}{D}_{\mathrm{KL}}(Q\Vert P)=\underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmin}}{\int }_{\mathbf{A}}Q(\mathbf{A})\log \frac{Q(\mathbf{A})}{Pr(\mathbf{A}\mid \mathbf{Y})}=\underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmax}}[\mathrm{H}(Q)+{⟨\log Pr(\mathbf{A}\mid \mathbf{Y})⟩}_Q]\end{array}$$

To obtain a tractable form of the above, first expand $\log Pr(\mathbf{A}\mid \mathbf{Y})$ using the logarithmic form of Bayes' rule:

$$\begin{array}{}\text{(23)}& \log Pr(\mathbf{A}\mid \mathbf{Y})=\log Pr(\mathbf{Y}\mid \mathbf{A})+\log Pr\left(\mathbf{A}\right)-\log Pr\left(\mathbf{Y}\right)\end{array}$$

This gives convenient simplifications, as the prior $\log Pr\left(\mathbf{A}\right)$ is often Gaussian and has closed-form solutions, and the marginal data likelihood $\log Pr\left(\mathbf{Y}\right)$ is constant and can be dropped from the optimization. Expanding ${⟨\log Pr(\mathbf{A}\mid \mathbf{Y})⟩}_Q$ in Eq. $\text{3}$ gives:

$$\begin{array}{}\text{(24)}& \underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmax}}[\mathrm{H}(Q)+{⟨\log Pr(\mathbf{Y}\mid \mathbf{A})⟩}_Q+{⟨\log Pr\left(\mathbf{A}\right)⟩}_Q-{⟨\log Pr\left(\mathbf{Y}\right)⟩}_Q]\end{array}$$

Dropping the constant ${⟨\log Pr\left(\mathbf{Y}\right)⟩}_Q$ term and recognizing that the remaining terms reflect the KL divergence of the approximating posterior from the prior, i.e. $D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))$, we get the following optimization problem:

$$\begin{array}{}\text{(25)}& \underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmax}}[{⟨\log Pr(\mathbf{Y}\mid \mathbf{A})⟩}_Q-D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))].\end{array}$$

The objective function for variational Bayes amounts to maximizing the data likelihood ${⟨\log Pr(\mathbf{Y}\mid \mathbf{A})⟩}_Q$ under the approximation $Q$, while also minimizing the KL divergence of the prior from the approximating posterior. It can therefore be interpreted as a regularized maximum-likelihood approach. This form also connects to the objective functions often seen in variational autoencoders and in variational free energy. 

In this application both the prior and approximating posterior are Gaussian, and the KL divergence term $D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))$ has a closed-form solution reflecting the KL divergence between two multivariate Gaussian distributions:

$$\begin{array}{}\text{(26)}& D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))=\frac{1}{2}[\log \frac{|{\mathrm{\Sigma }}_A|}{|{\hat{\mathrm{\Sigma }}}_A|}-D+\text{tr}[{\mathrm{\Sigma }}_A^{-1}{\hat{\mathrm{\Sigma }}}_A]+({\mu }_A-{\hat{\mu }}_A{)}^T{\mathrm{\Sigma }}_A^{-1}({\mu }_A-{\hat{\mu }}_A)],\end{array}$$

where $D$ is the dimensionality of the multivariate Gaussian. It remains then to calculate the expected log-likelihood, ${⟨\log Pr(\mathbf{Y}\mid \mathbf{A})⟩}_Q$. As discussed in the next section, this integral is not always tractable.

Challenges for the variational approximation in this application

The variational approximation integrates over the domain for $A$, which is truncated to $[0,\mathrm{\infty })$ since negative values for $\mathbf{A}$ are unphysical. Typically, this means that efficient algorithms are challenging to derive, as closed-form solutions for the relevant integrals do not exist, or at best involve the multivariate Gaussian cumulative distribution function, its inverses, and derivatives, which are numerically expensive to compute. 

One may relax the constraint that $\mathbf{A}$ be non-negative, extending the domain of integration to $(-\mathrm{\infty },\mathrm{\infty })$, but then one must constrain optimization to return only positive means for the variational posterior. However, unless a rectifying (e.g. exponential, quadratic) link function is used, the inclusion of negative rates in the domain will make the Poisson observation likelihood undefined. For this reason, the variational update has been explored only for the exponential link function [[PARK]]. Because small changes in activation can lead to large fluctuations in rate owing to the amplification of the exponential link function, we have found that the exponential link is numerically unstable. 

(An implementation of variational optimization using the rectifying quadratic link function may be more numerically stable, and remains to be explored.)

Moment-matching

Moment matching calculates or approximates the mean and covariance of the true posterior, and uses these moments to form a multivariate Gaussian approximation. When applied as a message-passing algorithm in a graphical model, moment matching is an important step of the expectation-propagation algorithm. Moment-matching can be performed explicitly by integrating the posterior moments, but in high dimensions there is no computationally tractable way to evaluate such integrals. Since spatial correlations are essential in spatiotemporal phenomena, we cannot discard this higher dimensional structure. 

Another approach to moment-matching to note is that the the Gaussian distribution $Q$ that minimizes KL divergence from $Q$ to the true posterior will also match the moments of the posterior. We can therefore perform moment matching by minimizing this KL divergence:

$$\begin{array}{}\text{(27)}& \underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmin}}{\int }_{\mathbf{A}}Pr(\mathbf{A}\mid \mathbf{Y})\log \frac{Pr(\mathbf{A}\mid \mathbf{Y})}{Q(\mathbf{A})}=\underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmax}}\mathrm{H}[Pr(\mathbf{A}\mid \mathbf{Y})]+{⟨\log Q⟩}_{Pr(\mathbf{A}\mid \mathbf{Y})}\end{array}$$

Note that the first term is the entropy of the (true) posterior distribution. It is constant for a given update, and therefore does not affect our optimization. We can focus on the second term, and optimize:

$$\begin{array}{}\text{(28)}& \underset{{\mu }_Q,{\mathrm{\Sigma }}_Q}{\mathrm{argmax}}{⟨\log Q⟩}_{Pr(\mathbf{A}\mid \mathbf{Y})}\end{array}$$

The log-probability of a Gaussian approximation $Q$ is with mean ${\hat{\mu }}_A$ and covariance matrix ${\hat{\mathrm{\Sigma }}}_A$ is:

$$\begin{array}{}\text{(29)}& \log Q\left(A\right)=-\frac{1}{2}[\log |2\pi {\mathrm{\Sigma }}_A|+(\mathbf{A}-{\hat{\mu }}_A{)}^{\mathrm{\top }}{{\hat{\mathrm{\Sigma }}}_A}^{-1}(\mathbf{A}-{\hat{\mu }}_A)]\end{array}$$

We cannot calculate the true posterior $Pr(\mathbf{A}\mid \mathbf{Y})$, and so the integral ${⟨\log Q⟩}_{Pr(\mathbf{A}\mid \mathbf{Y})}$ cannot be computed directly. However, the normalization constant, although unknown, is constant with respect to this optimization, and is suffices to take the weighted expectation with respect to an un-normalized form of $Pr(\mathbf{A}\mid \mathbf{Y})$. 

$$\begin{array}{}\text{(30)}& \begin{array}{rl}Pr\left(\mathbf{A}\mid \mathbf{Y}\right)& ={\left|2\pi {\mathrm{\Sigma }}_A\right|}^{-\frac{1}{2}}{e}^{-\frac{1}{2}(\mathbf{A}-{\mu }_a{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}(\mathbf{A}-{\mu }_a)}\cdot \prod _{n\in 1..N}\left[\frac{({\gamma }_{n}f(A_n)+{\beta }_n{)}^{y_n}}{\exp ({\gamma }_{n}f(A_n)+{\beta }_n)}\right]\end{array}\end{array}$$

This integral, however, remains essentially intractable, due to the product of Gaussian and Poisson-like terms. In high dimension there is no (to our knowledge) computationally efficient way to estimate this integral or its derivatives. 

Expected log-likelihoods and variational Bayes

So far, we have explored three approaches to finding a Gaussian approximation to the measurement posterior: the Laplace approximation, variational Bayes, and moment matching. Moment matching is unsuitable because, to the best of our knowledge, there is no computationally tractable way to estimate the relevant moments in high-dimensions. The Laplace approximation and variational Bayes remain computationally tractable, with limitations.

The Laplace approximation suffers from errors arising from the non-negativity constraint on activity levels, and the high skewness of the distributions causes the mode to be far from the mean. Errors in estimating the covariance are especially severe, as the covariance controls the trade-off between propagating past information, and incorporating new measurements, during the filtering procedure.

Variational Bayes also has a number of challenges. First, an efficient way of evaluating the relevant integrals and their derivatives is needed. In practice, this is simplest when we can tolerate the approximation of integrating over the full domain of the prior, including the unphysical negative activations. We must also use a rectifying link function, because if the predicted point-process intensity is negative, then the Poisson likelihood for our observations is not defined. 

To our knowledge, the only rectifying link function that has been explored to-date is the exponential link function, which suffers from unacceptable numerical instability in our application. There are a few other link functions that we might explore, but we will address another approximate solution in this section based on Laplace-approximated expected log-likelihoods.

In order to minimize $D_{\text{KL}}$ in the variational Bayes approach, we must maximize the data log-likelihood under the approximating distribution $Q$, while simultaneously minimizing the KL divergence of the prior from this posterior approximation. Provided we interpret the multivariate Gaussian prior for $\mathbf{A}$ as having support over $(-\mathrm{\infty },\mathrm{\infty }{)}^D$, the KL divergence term has a closed-form solution and well-defined derivatives. The challenge, then, is to calculate expected log-likelihood term:

$$\begin{array}{}\text{(31)}& {⟨\log Pr(\mathbf{Y}\mid \mathbf{A})⟩}_Q\end{array}$$

We now derive a general second-order approximation for the expected log-likelihood for a Gaussian latent-variable process with Poisson observations, where the latent variable may be linked to the Poisson intensity via an arbitrary link function. (See Zhou and Park, plus the expected log-likehood papers, for more detail). 

In the case of intractable ${⟨\log Pr(\mathbf{Y}\mid \mathbf{A})⟩}_Q$, we approximate this integral via Laplace approximation. This yields and approximate variational inference method that is similar, but not identical, to the Laplace approximation.

Second-order approximations to the expected log-likelihood

To briefly review the notation, let $\mathbf{A}=\{A_1,..,A_N\}$ be a multivariate Gaussian latent variable reflecting our prior estimate for the distribution of latent activation, with mean ${\mu }_A$ and covariance ${\mathrm{\Sigma }}_A$. Let ${\lambda }_n^0=f(A_n)$ be an link function mapping the latent activity to a baseline intensity ${\lambda }_n^0$, which then might be further scaled and shifted due to spatially inhomogeneous gain ${\gamma }_n$ or background activity ${\beta }_n$. 

We need to compute expected log-likelihoods under the approximating posterior distribution $Q$, that is:

$$\begin{array}{}\text{(32)}& ⟨\mathcal{L}(y)n|A_n)⟩=y_n⟨\log ({\lambda }_n)⟩-⟨{\lambda }_n⟩\end{array}$$

Where $⟨\cdot ⟩$ denotes averaging over the posterior distribution $Q(\mathbf{A})$ with mean ${\hat{\mu }}_{\mathbf{A}}$ and covariance ${\hat{\mathrm{\Sigma }}}_{\mathbf{A}}$. 

In certain cases, the expectations $⟨\log ({\lambda }_n)⟩$ and $⟨{\lambda }_n⟩$ may have closed-form solutions, for example in the log-Gaussian instance (cite Rule, Zhou). Here, however, we explore a general approach based on second-order Taylor expansions, which is accurate for small variances. If $A_n$ is normally distributed with mean ${{\mu }_A}_n$ and variance ${{{\sigma }_A}_n}^2$, then out to second order:

$$\begin{array}{}\text{(33)}& \begin{array}{rl}⟨\mathcal{L}\left(y_n\mid A_n\right)⟩& \approx \mathcal{L}\left(y_n\mid {\hat{\mu }}_{A_n}\right)+\frac{{\hat{\sigma }}_{A_n}^2}{2}{\mathcal{L}}^{″}\left(y_n\mid {\hat{\mu }}_{A_n}\right)\end{array}\end{array}$$

Gradient-based methods for optimizing the expected log-likelihood require derivatives of these approximated expectations. In general, derivatives with respect to the mean are:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\frac{d^n}{d{\mu }_{A_n}^m}⟨\mathcal{L}(y_{A_n}|A_n)⟩& \approx {\mathcal{L}}^{(m)}\left(y_n\mid {\hat{\mu }}_{A_n}\right)+\frac{{\hat{\sigma }}_{A_n}^2}{2}{\mathcal{L}}^{(m+2)}\left(y_n\mid {\hat{\mu }}_{A_n}\right)\end{array}\end{array}$$

A Newton-Raphson solver for optimizing the mean ${\hat{\mu }}_{\mathbf{A}}$ requires the Hessian of the objective function. Since the approximated expectations include second derivatives, the Hessian involves derivatives out to fourth order. The chain rule for higher-order derivatives of the logarithm is too cumbersome to state for the general case. Instead, we derive the equations for three versions of $f(A)$: $A$, $A^2$, and $e^A$.

Optimizing the likelihood may also involve optimizing the variance ${\sigma }^2$ or in general, the covariance. We will address this in later sections.

For the case that $f(A)=A$

Interpreting the distribution of latent activations $\mathbf{A}$ as a multivariate Gaussian over $(-\mathrm{\infty },\mathrm{\infty }{)}^D$ allows closed-form estimation of the $D_{\text{KL}}$ contribution to the variational Bayes objective function. However, unless point-process intensities are artificially constrained to the domain $[0,\mathrm{\infty })$, the expected log-likelihood is undefined for the identity link function. This is because the Poisson measurement likelihood is not defined for negative intensities. 

In this second-order approximation, we circumvent this issue by considering a locally-quadratic approximation of the likelihood function that continues the Poisson likelihood to negative intensities. Provided variance is small, and the posterior mean is constrained to be positive, this approximation may provide an accurate estimate of the expected log-likelihood.

If $f(A)=A$ and so $\lambda =v\cdot \gamma \cdot (x+\beta /\gamma )$. Computing out to the 4th derivative.

$$\begin{array}{}\text{(35)}& \begin{array}{rl}\mathcal{L}\left(y\mid A\right)& =y\log [\lambda ]-\lambda \\ & =y[\log (v\cdot \gamma )+\log (\lambda /\gamma )]-v\cdot \gamma (x+\beta /\gamma )\\ & =y\log (\lambda /\gamma )-v\cdot \gamma A+\text{constant}\\ {\mathcal{L}}^{\prime }\left(y\mid A\right)& =y(\lambda /\gamma {)}^{-1}-v\cdot \gamma \\ {\mathcal{L}}^{″}\left(y\mid A\right)& =-y(\lambda /\gamma {)}^{-2}\\ {\mathcal{L}}^{(3)}\left(y\mid A\right)& =2y(\lambda /\gamma {)}^{-3}\\ {\mathcal{L}}^{(4)}\left(y\mid A\right)& =-6y(\lambda /\gamma {)}^{-4}\end{array}\end{array}$$

As $⟨\lambda ⟩\to 0$, the fourth derivative of the likelihood tends rapidly to infinity, which may create issues for numerical stability and accuracy. This behavior near $⟨\lambda ⟩\to 0$ is similar to the issues that plague the Laplace approximation. In particular, the distribution may become highly skewed, which means that third or higher moments may be needed, and the second-order approximation may not be accurate. However, I have reason to suspect that the issues might be less severe for the expected log-likelihood compared to the Laplace approximation. In this case, we are using a quadratic expansion about an estimate of the posterior mean, whereas the Laplace approximation seeks the posterior mode. I expect that this will have a stabilizing effect, encouraging $⟨\lambda ⟩$ toward more positive values.

For the case that $f(A)=e^A$

A closed-form solution for the expected log-likelihood exists under this link function, and closed-form expressions for the moments of log-Gaussian random variables are known (see Zhou and Park for application to log-Gaussian point processes). However, in this application the amplification of positive tails of the distribution by the exponential link function is numerically unstable and unphysical, indicating that the log-Gaussian model is inappropriate. In Rule et al. 2018, we noted that a second-order approximation to the expected likelihood was more stable and more accurate. 

$$\begin{array}{}\text{(36)}& \begin{array}{rl}\mathcal{L}\left(y\mid A\right)& =y\log [\lambda ]-\lambda \\ & =y\log [v\cdot \gamma (e^A+\beta /\gamma )]-v\cdot \gamma (e^A+\beta /\gamma )\\ & =y\log [e^A+\beta /\gamma ]-v\cdot \gamma e^A+\text{constant}\\ {\mathcal{L}}^{(n)}\left(y\mid A\right)& =y{C}^{(n-1)}-v\cdot \gamma e^A,\\ C& =\frac{e^A}{e^A+\beta /\gamma }=\frac{1}{1+\beta /\gamma e^{-A}}\\ C^{\prime }& =C(1-C)\\ C^{″}& =C^{\prime }(1-2C)\\ C^{(3)}& =C^{\prime }-6{C^{\prime }}^2\end{array}\end{array}$$

For the case that $f(A)=A^2$

A quadratic link function is rectifying, so the Poisson likelihood remains well-defined even if the domain of $\mathbf{A}$ is extended to $(-\mathrm{\infty },\mathrm{\infty }{)}^D$. However, to my knowledge there is no tractable closed-form for the expectation of logarithm of a generalized noncentral ${\chi }^2$ distributed variable, and so the second-order approximation remains useful:

$$\begin{array}{}\text{(37)}& \begin{array}{rl}\mathcal{L}\left(y\mid A\right)& =y\log [\lambda ]-\lambda \\ & =y\log [v\cdot \gamma (A^2+\beta /\gamma )]-v\cdot \gamma (A^2+\beta /\gamma )\\ & =y\log [A^2+\beta /\gamma ]-v\gamma A^2+\text{constant}\\ {\mathcal{L}}^{\prime }\left(y\mid A\right)& =yC-2v\gamma A,& C& =\frac{2A}{A^2+\beta /\gamma }\\ {\mathcal{L}}^{″}\left(y\mid A\right)& =y{C}^{\prime }-2v\gamma ,& C^{\prime }& =C(C^{-1}-C)\\ {\mathcal{L}}^{(3)}\left(y\mid A\right)& =y{C}^{″},& C^{″}& =C^{\prime }(A^{-1}-2C)-C{A}^{-2}\\ {\mathcal{L}}^{(4)}\left(y\mid A\right)& =y{C}^{(3)},& C^{(3)}& =3(C/A{)}^2-6{C^{\prime }}^2\end{array}\end{array}$$

Incorporating the prior

So far we have focused on the expected log-likelihood contribution to the variational posterior objective function. We also need to derive the gradients of the KL-divergence term. For the most part, this is identical to the derivation in Zhou and Park, so I present only some quick notes here. We are interested in the gradients (and hessians) for the KL divergence between two $k$ dimensional multivariate Gaussians, which is:

$$\begin{array}{}\text{(38)}& D_{KL}({\mathcal{N}}_0\Vert {\mathcal{N}}_1)=\frac{1}{2}[\mathrm{tr}({\mathrm{\Pi }}_1{\mathrm{\Sigma }}_0)+({\mu }_1-{\mu }_0{)}^{\mathrm{\top }}{\mathrm{\Pi }}_1({\mu }_1-{\mu }_0)-\ln \left|{\mathrm{\Pi }}_1\right|-\ln \left|{\mathrm{\Sigma }}_0\right|-k]\end{array}$$

Note that I have chosen to denote this in terms of the precision matrix ${\mathrm{\Pi }}_1={\mathrm{\Sigma }}_1^{-1}$, as it makes some of the derivations below more straightforward. 

The derivative with respect to ${\mu }_0$ is the same as the derivative with respect to ${\mu }_1$ and is:

$$\begin{array}{}\text{(39)}& {\mathrm{\nabla }}_{{\mu }_1}{D}_{KL}({\mathcal{N}}_0\Vert {\mathcal{N}}_1)={\mathrm{\nabla }}_{{\mu }_0}{D}_{KL}({\mathcal{N}}_0\Vert {\mathcal{N}}_1)=\frac{1}{2}[({\mu }_1-{\mu }_0{)}^{\mathrm{\top }}{\mathrm{\Pi }}_1]\end{array}$$

The derivative with respect to ${\mathrm{\Sigma }}_0$ is:

$$\begin{array}{}\text{(40)}& {\mathrm{\nabla }}_{{\mathrm{\Sigma }}_0}{D}_{KL}({\mathcal{N}}_0\Vert {\mathcal{N}}_1)=\frac{1}{2}[{\mathrm{\Pi }}_1-{\mathrm{\Pi }}_0]\end{array}$$

The derivative with respect to ${\mathrm{\Pi }}_1$ is:

$$\begin{array}{}\text{(41)}& {\mathrm{\nabla }}_{{\mathrm{\Pi }}_1}{D}_{KL}({\mathcal{N}}_0\Vert {\mathcal{N}}_1)=\frac{1}{2}[{\mathrm{\Sigma }}_0-{\mathrm{\Sigma }}_1+({\mu }_1-{\mu }_0)({\mu }_1-{\mu }_0{)}^{\mathrm{\top }}]\end{array}$$

For the variational interpretation we approximate the posterior $P$ with approximation $Q$ and minimize $D_{KL}(Q\Vert P)$. This involves taking the derivative with respect to ${\mathrm{\Sigma }}_0$ above. 

Optimizing the (approximate) variational approximation

We have derived approximations for the expected log-likelihood contribution to the variational Bayesian objective function, which must be optimized jointly over the posterior mean ${\hat{\mu }}_A$ and the posterior covariance ${\hat{\mathrm{\Sigma }}}_A$. The above derivations provide gradients and Hessians for optimizing ${\hat{\mu }}_A$ for a fixed ${\hat{\mathrm{\Sigma }}}_A$. In Zhou and Park, they explore the joint optimization for the (exact) objective function for a log-Gaussian variational approximation. They prove that ${\hat{\mathrm{\Sigma }}}_A$ can be optimized using a fixed-point iteration. 

(In numerical experiments, I extended this approach by interleaving one-step of the Newton-Raphson optimization for ${\hat{\mu }}_A$ with one step of the fixed-point update for ${\hat{\mathrm{\Sigma }}}_A$. In my experience this accelerated convergence. Does the fixed-point iteration for ${\hat{\mathrm{\Sigma }}}_A$ convergs for the second-oder approximated expected log-likelihood? Could a similar approach be found for the other (non-exponential) link functions explored here.)

The covariance update

To complete the variational approximation, we also need to optimize the posterior covariance. This involves the derivative of the expected log-likelihood with respect to ${\mathrm{\Sigma }}_0$. In the second-order (Laplace-approximation-like) expected log-likelihood, the dependence on the covariance enters through the second-order terms, which are:

$$\begin{array}{}\text{(42)}& \frac{1}{2}\mathrm{diag}\left({\mathrm{\Sigma }}_0\right){\mathcal{L}}^{″}\left(y\mid {\mu }_0\right)\end{array}$$

The derivative of the above with respect to ${\mathrm{\Sigma }}_0$ is:

$$\begin{array}{}\text{(43)}& \frac{1}{2}{\mathcal{L}}^{″}\left(y\mid {\mu }_0\right)\end{array}$$

along the diagonal, and 0 elsewhere. The total gradient for ${\mathrm{\Sigma }}_0$, incorporating both the $D_{\text{KL}}$ and expected log-likelihood contributions, is: 

$$\begin{array}{}\text{(44)}& {\mathrm{\nabla }}_{{\mathrm{\Sigma }}_0}\mathcal{L}(\mathbf{A}\mid \mathbf{Y})=\frac{1}{2}[{\mathrm{\Pi }}_1-{\mathrm{\Pi }}_0+{\mathcal{L}}^{″}\left(y\mid {\mu }_0\right)]\end{array}$$

Setting this gradient to zero and solving for ${\mathrm{\Sigma }}_0$ gives:

$$\begin{array}{}\text{(45)}& {\mathrm{\Pi }}_0={\mathrm{\Pi }}_1+{\mathcal{L}}^{″}\left(y\mid \mu \right)\end{array}$$

Which amounts to adding the curvature of the log likelihood ${\mathcal{L}}^{″}\left(y\mid \mu \right)$ (approximated as the curvature at the current posterior mean), to the prior precision matrix ${\mathrm{\Pi }}_1$. This is similar to the Hessian observed for the Laplace update, with the exception that we use the curvature at the estimated posterior mean, rather than posterior mode.

Computing the model likelihood

Once the Gaussian approximation is computed, how should we estimate the likelihood of the data given the model parameters $\theta$, $Pr(\mathbf{Y}|\theta )$ (or just $Pr(\mathbf{Y})$ for short)? There are numerous approximation methods available, and it remains unclear to me which is best.

Integration via Laplace approximation

This likelihood is the integral of the prior $Pr(\mathbf{A})$ times the measurement likelihood $Pr(\mathbf{Y}\mid \mathbf{A})$, and is also the normalization constant for the posterior distribution: 

$$\begin{array}{}\text{(46)}& Pr(\mathbf{Y})={\int }_{\mathbf{A}}Pr(\mathbf{Y}\mid \mathbf{A})Pr(\mathbf{A})={⟨Pr(\mathbf{Y}\mid \mathbf{A})⟩}_{Pr(\mathbf{A})}\end{array}$$

Given a Gaussian approximation $Q$ for the posterior $Pr(\mathbf{A}\mid \mathbf{Y})$, we can approximate this integral using the posterior mean ${\hat{\mu }}_A$ and covariance ${\hat{\mathrm{\Sigma }}}_A$. This Gaussian approximation can be obtained from any of the previously mentioned approximation methods, Laplace approximation, or exact or approximated variational inference. 

Under this approximation, we evaluate $Pr(\mathbf{Y}\mid \mathbf{A})Pr(\mathbf{A})$ at ${\hat{\mu }}_A$, and also compute the curvature at this point (which should match ${\hat{\mathrm{\Sigma }}}_A$ if everything has gone as planned!).

$$\begin{array}{}\text{(47)}& \log Pr(\mathbf{A}\mid \mathbf{Y})\approx -\frac{1}{2}[\log |2\pi {\hat{\mathrm{\Sigma }}}_A|+(\mathbf{A}-{\hat{\mu }}_a{)}^{\mathrm{\top }}{\hat{\mathrm{\Sigma }}}_A^{-1}(\mathbf{A}-{\hat{\mu }}_a)]\end{array}$$

We use the following logarithmic relationship derived from Bayes' rule

$$\begin{array}{}\text{(48)}& \log Pr(\mathbf{Y})=\log [Pr(\mathbf{A})Pr(\mathbf{Y}\mid \mathbf{A})]-\log Pr(\mathbf{A}\mid \mathbf{Y})\end{array}$$

Substituting the forms for the above

$$\begin{array}{}\text{(49)}& \begin{array}{rl}\log Pr(\mathbf{Y})& \approx -\frac{1}{2}[\log |2\pi {\mathrm{\Sigma }}_A|+(\mathbf{A}-{\mu }_A{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}(\mathbf{A}-{\mu }_A)]\\ & +\sum _{n\in 1..N}[y_n\log ({\gamma }_{n}f(A_n)+{\beta }_n)-({\gamma }_{n}f(A_n)+{\beta }_n)]\\ & +\frac{1}{2}[\log |2\pi {\hat{\mathrm{\Sigma }}}_A|+(\mathbf{A}-{\hat{\mu }}_A{)}^{\mathrm{\top }}{\hat{\mathrm{\Sigma }}}_A^{-1}(\mathbf{A}-{\hat{\mu }}_A)]\end{array}\end{array}$$ 

Evaluating the above at the posterior mean ${\hat{\mu }}_A$ amounts to a Laplace approximation of the integral for the likelihood.

$$\begin{array}{}\text{(50)}& \begin{array}{rl}\log Pr(\mathbf{Y})& \approx \log [Pr({\hat{\mu }}_A)Pr(\mathbf{Y}\mid {\hat{\mu }}_A)]-\log Pr({\hat{\mu }}_A\mid \mathbf{Y})\\ & =-\frac{1}{2}[\log |2\pi {\mathrm{\Sigma }}_A|+({\hat{\mu }}_A-{\mu }_A{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}({\hat{\mu }}_A-{\mu }_A)]\\ & +\sum _{n\in 1..N}[y_n\log ({\gamma }_{n}f({\hat{\mu }}_{A_n})+{\beta }_n)-({\gamma }_{n}f({\hat{\mu }}_{A_n})+{\beta }_n)]\\ & +\frac{1}{2}[\log |2\pi {\hat{\mathrm{\Sigma }}}_A|+({\hat{\mu }}_A-{\hat{\mu }}_A{)}^{\mathrm{\top }}{\hat{\mathrm{\Sigma }}}_A^{-1}({\hat{\mu }}_A-{\hat{\mu }}_A)]\end{array}\end{array}$$ 

Cleaning things up, and denoting ${\hat{\lambda }}_n={\gamma }_{n}f({\hat{\mu }}_{A_n})+{\beta }_n$, we get:

$$\begin{array}{}\text{(51)}& \begin{array}{rl}\log Pr(\mathbf{Y})& \approx -\frac{1}{2}[\log \frac{|{\mathrm{\Sigma }}_A|}{|{\hat{\mathrm{\Sigma }}}_A|}+({\hat{\mu }}_A-{\mu }_A{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_A^{-1}({\hat{\mu }}_A-{\mu }_A)]\\ & +\sum _{n\in 1..N}[y_n\log ({\hat{\lambda }}_n)-{\hat{\lambda }}_n]\end{array}\end{array}$$ 

On the similarity between the Laplace-approximated likelihood and $D_{\text{KL}}$

Note the similarity to KL divergence

$$\begin{array}{}\text{(52)}& D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))=\frac{1}{2}[\log \frac{|{\mathrm{\Sigma }}_A|}{|{\hat{\mathrm{\Sigma }}}_A|}-D+\text{tr}\left[{\mathrm{\Sigma }}_A^{-1}{\hat{\mathrm{\Sigma }}}_A\right]+({\mu }_A-{\hat{\mu }}_A{)}^T{\mathrm{\Sigma }}_A^{-1}({\mu }_A-{\hat{\mu }}_A)]\end{array}$$

Which gives the relation

$$\begin{array}{}\text{(53)}& -\frac{1}{2}[\log \frac{|{\mathrm{\Sigma }}_A|}{|{\hat{\mathrm{\Sigma }}}_A|}+({\mu }_A-{\hat{\mu }}_A{)}^T{\mathrm{\Sigma }}_A^{-1}({\mu }_A-{\hat{\mu }}_A)]=\frac{1}{2}[\text{tr}\left[{\mathrm{\Sigma }}_A^{-1}{\hat{\mathrm{\Sigma }}}_A\right]-D]-D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))\end{array}$$

Which allows the likelihood to be written as

$$\begin{array}{}\text{(54)}& \begin{array}{rl}\log Pr(\mathbf{Y})& \approx \frac{1}{2}[\text{tr}\left[{\mathrm{\Sigma }}_A^{-1}{\hat{\mathrm{\Sigma }}}_A\right]-D]-D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))\\ & +\sum _{n\in 1..N}[y_n\log ({\hat{\lambda }}_n)-{\hat{\lambda }}_n]\end{array}\end{array}$$ 

We can play with this further, bringing in the second-order expected log-likelihood:

$$\begin{array}{}\text{(55)}& \mathcal{L}\left(y_n\mid {\hat{\mu }}_{A_n}\right)=\sum _{n\in 1..N}[y_n\log ({\hat{\lambda }}_n)-{\hat{\lambda }}_n]\end{array}$$

Which, from the second-order approximation, is:

$$\begin{array}{}\text{(56)}& \begin{array}{rl}\mathcal{L}\left(y_n\mid {\hat{\mu }}_{A_n}\right)& \approx ⟨\mathcal{L}\left(y_n\mid {\hat{\mu }}_{A_n}\right)⟩-\frac{{\hat{\sigma }}_{A_n}^2}{2}{\mathcal{L}}^{″}\left(y_n\mid {\hat{\mu }}_{A_n}\right)\end{array}\end{array}$$

So:

$$\begin{array}{}\text{(57)}& \begin{array}{rl}\log Pr(\mathbf{Y})& \approx \frac{1}{2}[\text{tr}\left[{\mathrm{\Sigma }}_A^{-1}{\hat{\mathrm{\Sigma }}}_A\right]-D]-D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}))\\ & +⟨\mathcal{L}\left(y\mid {\hat{\mu }}_A\right)⟩-\sum _{n\in 1..N}\frac{{\hat{\sigma }}_{A_n}^2}{2}{\mathcal{L}}^{″}\left(y_n\mid {\hat{\mu }}_{A_n}\right)\\ & =-D_{\mathrm{KL}}(Q\Vert Pr(\mathbf{A}\mid \mathbf{Y}))+\frac{1}{2}\text{tr}\left({\mathrm{\Sigma }}_A^{-1}{\hat{\mathrm{\Sigma }}}_A\right)-\sum _{n\in 1..N}\frac{{\hat{\sigma }}_{A_n}^2}{2}{\mathcal{L}}^{″}\left(y_n\mid {\hat{\mu }}_{A_n}\right)\\ & +\text{constant}\end{array}\end{array}$$ 

Which is to say that a point estimate of the log-likelihood is very similar to the (negative) KL-divergence penalty, differing only by the trace term and the curvature correction (and the dimensionality constant $D$). 

Empirically, these terms are small and do not dominate the likelihood. As we shall see in the following sections, this similarity is not a coincidence: the Laplace approximation is connected to the Evidence Lower Bound (ELBO) in the variational Bayesian approach, especially when a second-order approximation is used to evaluate the expected log-likelihood.

The expected log-likelihood

One can consider the expected value of the log-likelihood relative to the prior distribution $Pr(\mathbf{A})$. Starting from the logarithmic form of Bayes' rule, we have:

$$\begin{array}{}\text{(4)}& \begin{array}{r}\log Pr(\mathbf{Y})=\log Pr(\mathbf{A})+\log Pr(\mathbf{Y}\mid \mathbf{A})-\log Pr(\mathbf{A}\mid \mathbf{Y})\end{array}\end{array}$$

For the true posterior $Pr(\mathbf{A}\mid \mathbf{Y})$, this equality holds for all $\mathbf{A}$, as we could recover the log-likelihood by evaluating this expression at any point. In the Laplace approximation, we evaluated this quantity at the posterior mean. 

%For the variational approach (outlined below), we consider a bound based on taking the expectation of Eq. $\text{4}$ with respect to the approximating posterior distribution $Q$.

For the expected log-likelihood approach here, we take the expectation with respect to the prior distribution $Pr(\mathbf{A})$:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{⟨\log P_Y⟩}_{P_A}& ={⟨\log P_A⟩}_{P_A}+{⟨\log P_{Y|A}⟩}_{P_A}-{⟨\log P_{A|Y}⟩}_{P_A}\\ & ={⟨\log P_{Y|A}⟩}_{P_A}+D_{\text{KL}}(P_A\Vert P_{A|Y})\end{array}\end{array}$$

We cannot compute the above exactly, because we do not have access to the true posterior $P_{A|Y}$, but we do have access to an approximating posterior $Q\approx P_{A|Y}$, which we can use to approximate the expected log-likelihood. Note that the additional $D_{\text{LK}}$ term increases the expected log likelihood, opposite of its role in the variational approach. 

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{⟨\log P_Y⟩}_{P_A}& \approx {⟨\log P_{Y|A}⟩}_{P_A}+D_{\text{KL}}(P_A\Vert Q)\end{array}\end{array}$$

To estimate the expected log-likelihood term ${⟨\log P_{Y|A}⟩}_{P_A}$, we may either use the second-order approach that we derived for variational Bayes', or simply evaluate $\log P_{Y|A}$ at the prior mean for a faster approximation.

ELBO and variational Bayes

When deriving the variational update, we omitted the (constant) data likelihood term. The full form for $D_{\text{KL}}(Q\Vert P)$ is:

$$\begin{array}{}\text{(58)}& \begin{array}{rl}{D}_{\text{KL}}(Q\Vert P)={\int }_A{Q}_A\log \frac{Q_A}{P_{A|Y}}& ={⟨\log Q⟩}_Q-{⟨\log P_{A|Y}⟩}_Q\\ & =-H_Q-{⟨\log P_{Y,A}⟩}_Q+\log P_Y\\ & =-{⟨\log P_{Y|A}⟩}_Q+D_{\text{KL}}(Q\Vert P_A)+\log P_Y,\end{array}\end{array}$$ 

which implies the following identity for the log-likelihood:

$$\begin{array}{}\text{(59)}& \begin{array}{rl}\log P_Y& =D_{\text{KL}}(Q\Vert P)+{⟨\log P_{Y|A}⟩}_Q-D_{\text{KL}}(Q\Vert P_A)\\ & =D_{\text{KL}}(Q\Vert P)+H_Q+{⟨\log P_{Y,A}⟩}_Q.\end{array}\end{array}$$ 

That is: the data likelihood is the expected log-likelihood under the variational posterior, minus the KL divergence of the posterior form the prior, plus the KL divergence of the variational posterior from the true posterior. Since this last term is always positive, the following bound holds:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\log P_Y& \ge {⟨\log P_{Y|A}⟩}_Q-D_{\text{KL}}(Q\Vert P_A)\\ & =H_Q+{⟨\log P_{Y,A}⟩}_Q\end{array}\end{array}$$ 

This inequality is sometimes called the Evidence Lower Bound (ELBO). The more accurately we can approximate the posterior as a Gaussian, the smaller $D_{\text{KL}}(Q\Vert P)$ becomes and the tight this bound becomes. 

Going forward

In practice, we use the Laplace approximation for obtaining the approximate posterior $Q$. This avoids needing to jointly optimize the posterior covariance, trading accuracy for speed. We have explore the Laplace, ELBO, and expected log-likelihood approaches to estimating the likelihood, and found them to be broadly similar. Note, however, that the accuracy of expressions involving second-order expansions is expected to break down if the variance of the latent state becomes large.