Note: Derivatives of Gaussian KL-Divergence for some parameterizations of the posterior covariance for variational Gaussian-process inference

These notes provide the derivatives of the KL-divergence $D_{\text{KL}}[Q(\mathbf{z})\Vert P(\mathbf{z})]$ between two multivariate Gaussian distributions $Q(\mathbf{z})$ and $P(\mathbf{z})$ with respect to a few parameterizations $\theta$ of the covariance matrix $\mathbf{\Sigma }(\theta )$ of $Q$. This is useful for variational Gaussian process inference, where clever parameterizations of the posterior covariance are required to make the problem tractable. Tables for differentiating matrix-valued functions can be found in The Matrix Cookbook .

Consider two multivariate Gaussian distributions $Q(\mathbf{z})=\mathcal{N}({\mathit{\mu }}_q,\mathbf{\Sigma }(\theta ))$ and $P(\mathbf{z})=\mathcal{N}({\mathit{\mu }}_0,{\mathbf{\Sigma }}_0={\mathbf{\Lambda }}^{-1})$ with dimension $L$. The KL divergence $D_{\text{KL}}[Q(\mathbf{z})\Vert P(\mathbf{z})]$ has the closed form

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathcal{D}:=& D_{\text{KL}}[Q(\mathbf{z})\Vert Pr(\mathbf{z})]\\ =& \frac{1}{2}\{({\mathit{\mu }}_0-{\mathit{\mu }}_q{)}^{\mathrm{\top }}\mathbf{\Lambda }({\mathit{\mu }}_0-{\mathit{\mu }}_q)\\ & +\mathrm{tr}\left(\mathbf{\Lambda }\mathbf{\Sigma }\right)-\ln |\mathbf{\Sigma }|-\ln |\mathbf{\Lambda }|\}+\text{constant.}\end{array}\end{array}$$

In variational Bayesian inference, we minimize $\mathcal{D}$ while maximizing the expected log-probability of some observations with respect to $Q(\mathbf{z})$. Closed-form derivatives of $\mathcal{D}$ in terms of the parameters of $Q$ are useful for manually optimizing code for larger problems. The derivatives of $\mathcal{D}$ in terms of ${\mathit{\mu }}_q$ are straightforward: ${\partial }_{{\mathit{\mu }}_q}\mathcal{D}=\mathbf{\Lambda }({\mathit{\mu }}_q-{\mathit{\mu }}_z)$ and ${\mathrm{H}}_{{\mathit{\mu }}_q}\mathcal{D}=\mathbf{\Lambda }$. In these notes, we explore derivatives of $\mathcal{D}$ with respect to a few different parameterizations ("$\theta$") of $\mathbf{\Sigma }(\theta )$.

We evaluate the following parameterizations for $\mathbf{\Sigma }$:

1. Optimizing the full $\mathbf{\Sigma }$ directly
2. $\mathbf{\Sigma }\approx \mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$
3. $\mathbf{\Sigma }\approx {\mathbf{A}}^{\mathrm{\top }}\mathrm{diag}[\mathbf{v}]\mathbf{A}$
4. $\mathbf{\Sigma }\approx [\mathbf{\Lambda }+\mathrm{diag}[\mathbf{p}]{]}^{-1}$
5. ${\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{\mathbf{Q}}^{\mathrm{\top }}\mathbf{F}$, where $\mathbf{Q}\in {\mathbb{R}}^{K\times K}$, $K<L$ and $\mathbf{F}\in {\mathbb{R}}^{K\times L}$, $\mathbf{F}{\mathbf{F}}^{\mathrm{\top }}=\mathbf{I}$

---

Optimizing $\mathbf{\Sigma }$ directly

We first obtain gradients of $\mathcal{D}$ in $\mathbf{\Sigma }$ (assuming $\mathbf{\Sigma }$ is full-rank). These can be used to derive gradients in $\theta$ for some parameterizations $\mathbf{\Sigma }(\theta )$ using the chain rule. The gradient of $\mathcal{D}$ in $\mathbf{\Sigma }$ can be obtained using identities (57) and (100) in The Matrix Cookbook:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\partial }_{\mathbf{\Sigma }}\mathcal{D}& ={\partial }_{\mathbf{\Sigma }}\{\mathrm{tr}\left(\mathbf{\Lambda }\mathbf{\Sigma }\right)-\ln |\mathbf{\Sigma }|\}\\ & =\frac{1}{2}(\mathbf{\Lambda }-{\mathbf{\Sigma }}^{-1}).\end{array}\end{array}$$

The Hessian in $\mathbf{\Sigma }$ is a fourth-order tensor. It's simpler to express the Hessian in terms of a Hessian-vector product, which can be used with Krylov subspace solvers to efficiently compute the update in Newton's method. Considering an $L\times L$ matrix $\mathbf{M}$, the Hessian-vector product is given by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\left[{\mathbf{H}}_{\mathbf{\Sigma }}\mathcal{D}\right]\mathbf{M}& ={\partial }_{\mathbf{\Sigma }}⟨{\partial }_{\mathbf{\Sigma }}\mathcal{D},\mathbf{M}⟩={\partial }_{\mathbf{\Sigma }}\mathrm{tr}[({\partial }_{\mathbf{\Sigma }}\mathcal{D}{)}^{\mathrm{\top }}\mathbf{M}],\end{array}\end{array}$$

where $⟨\cdot ,\cdot ⟩$ denotes the scalar (Frobenius) product. This is given by identity (124) in the Matrix Cookbook:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\partial }_{\mathbf{\Sigma }}\mathrm{tr}\left[\frac{1}{2}{(\mathbf{\Lambda }-{\mathbf{\Sigma }}^{-1})}^{\mathrm{\top }}\mathbf{M}\right]& =-\frac{1}{2}{\partial }_{\mathbf{\Sigma }}\mathrm{tr}\left[{\mathbf{\Sigma }}^{-1}\mathbf{M}\right]=\frac{1}{2}{\mathbf{\Sigma }}^{-1}{\mathbf{M}}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}.\end{array}\end{array}$$ 

 

---

Optimizing $\mathbf{\Sigma }\approx \mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$

We consider an approximate posterior covariance of the form

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathbf{\Sigma }& \approx \mathbf{X}{\mathbf{X}}^{\mathrm{\top }},\mathbf{X}\in {\mathbb{R}}^{L\times K}\end{array}\end{array}$$

where $\mathbf{X}$ is a rank-$K<L$ matrix with $L$ rows and $K$ columns.

Since $\mathbf{X}$ is not full rank, the log-determinant $\ln |\mathbf{\Sigma }|=\ln |\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}|$ in $\text{(1)}$ diverges, due to the zero eigenvalues in the null space of $\mathbf{X}$. However, since this null-space is not being optimized, it does not affect our gradient. It is sufficient to replace the log-determinant with that of the reduced-rank representation, $\ln |{\mathbf{X}}^{\mathrm{\top }}\mathbf{X}|$. Identity (55) in The Matrix Cookbook provides the derivative of this, ${\partial }_{\mathbf{X}}\ln |{\mathbf{X}}^{\mathrm{\top }}\mathbf{X}|=2{{\mathbf{X}}^{+}}^{\mathrm{\top }}$, where $(\cdot {)}^{+}$ is the pseudoinverse. Combined with identity (112), this gives the following gradient of $\mathcal{D}(\mathbf{X})$:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\partial }_{\mathbf{X}}\mathcal{D}& ={\partial }_{\mathbf{X}}\frac{1}{2}\{\mathrm{tr}\left[\mathbf{\Lambda }\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}\right]-\ln |{\mathbf{X}}^{\mathrm{\top }}\mathbf{X}|\}=\mathbf{\Lambda }\mathbf{X}-{{\mathbf{X}}^{+}}^{\mathrm{\top }}.\end{array}\end{array}$$

The Hessian-vector product requires the derivative of ${\partial }_{\mathbf{X}}\mathrm{tr}\left[{\mathbf{X}}^{+}\mathbf{M}\right]$:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\partial }_{\mathbf{X}}⟨{\partial }_{\mathbf{\Sigma }}\mathcal{D},\mathbf{M}⟩& ={\partial }_{\mathbf{X}}\mathrm{tr}\left[{(\mathbf{\Lambda }\mathbf{X}-{{\mathbf{X}}^{+}}^{\mathrm{\top }})}^{\mathrm{\top }}\mathbf{M}\right]={\partial }_{\mathbf{X}}\mathrm{tr}\left[\mathbf{\Lambda }\mathbf{X}\mathbf{M}\right]-{\partial }_{\mathbf{X}}\mathrm{tr}\left[{\mathbf{X}}^{+}\mathbf{M}\right].\end{array}\end{array}$$

Goulob and Pereya (1972) Eq. 4.12 gives the derivative of a fixed-rank pseudoinverse:

$$\begin{array}{}\text{(8)}& \begin{array}{r}\partial {\mathbf{X}}^{+}=-{\mathbf{X}}^{+}(\partial \mathbf{X}){\mathbf{X}}^{+}+{\mathbf{X}}^{+}\mathbf{X}{{}^{+}}^{\mathrm{\top }}(\partial \mathbf{X}{)}^{\mathrm{\top }}(1-\mathbf{X}{\mathbf{X}}^{+})+(1-{\mathbf{X}}^{+}\mathbf{X})(\partial \mathbf{X}{)}^{\mathrm{\top }}\mathbf{X}{{}^{+}}^{\mathrm{\top }}{\mathbf{X}}^{+}\end{array}\end{array}$$

Since $\mathbf{X}$ is $N\times K$ with rank $K$, ${\mathbf{X}}^{+}\mathbf{X}$ is full-rank. Therefore ${\mathbf{X}}^{+}\mathbf{X}={\mathbf{I}}_k$ and the final term in $\text{(8)}$ vanishes. The derivative of the pseudoinverse can now be written as:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\partial {\mathbf{X}}^{+}& =-{\mathbf{X}}^{+}(\partial \mathbf{X}){\mathbf{X}}^{+}+{\mathbf{X}}^{+}\mathbf{X}{{}^{+}}^{\mathrm{\top }}(\partial \mathbf{X}{)}^{\mathrm{\top }}({\mathbf{I}}_n-\mathbf{X}{\mathbf{X}}^{+})\end{array}\end{array}$$

Since the derivative of a trace of a matrix-valued function is just the (transpose) of the scalar derivative,

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\partial }_{\mathbf{X}}\mathrm{tr}\left[{\mathbf{X}}^{+}\mathbf{M}\right]& ={\{-{\mathbf{X}}^{+}\mathbf{M}{\mathbf{X}}^{+}+{\mathbf{X}}^{+}\mathbf{X}{{}^{+}}^{\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}({\mathbf{I}}_n-\mathbf{X}{\mathbf{X}}^{+})\}}^{\mathrm{\top }}\\ & =-{{\mathbf{X}}^{+}}^{\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{{\mathbf{X}}^{+}}^{\mathrm{\top }}+(\mathbf{I}-{{\mathbf{X}}^{+}}^{\mathrm{\top }}{\mathbf{X}}^{\mathrm{\top }})\mathbf{M}{\mathbf{X}}^{+}{{\mathbf{X}}^{+}}^{\mathrm{\top }}.\end{array}\end{array}$$

Overall, we obtain the following Hessian-vector product:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\partial }_{\mathbf{X}}⟨{\partial }_{\mathbf{\Sigma }}\mathcal{D},\mathbf{M}⟩& =\mathbf{\Lambda }{\mathbf{M}}^{\mathrm{\top }}+{{\mathbf{X}}^{+}}^{\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{{\mathbf{X}}^{+}}^{\mathrm{\top }}-(\mathbf{I}-{{\mathbf{X}}^{+}}^{\mathrm{\top }}{\mathbf{X}}^{\mathrm{\top }})\mathbf{M}{\mathbf{X}}^{+}{{\mathbf{X}}^{+}}^{\mathrm{\top }}\end{array}\end{array}$$

---

Optimizing $\mathbf{\Sigma }=\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$ when $\mathbf{X}$ is full-rank

Equations $\text{(6)}$ and $\text{(11)}$ are also valid if $\mathbf{X}$ is a rank-$L$ triangular (Choleskey) factorization of $\mathbf{\Sigma }$. In this case the pseudoinverse can be replaced by the full inverse, and various terms simplify:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\partial }_{\mathbf{X}}\mathcal{D}& =\mathbf{\Lambda }\mathbf{X}-{\mathbf{X}}^{-\mathrm{\top }}\\ {\partial }_{\mathbf{X}}⟨{\partial }_{\mathbf{x}}\mathcal{D},\mathbf{M}⟩& =\mathbf{\Lambda }{\mathbf{M}}^{\mathrm{\top }}+{\mathbf{X}}^{-\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{\mathbf{X}}^{-\mathrm{\top }}\end{array}\end{array}$$

---

Optimizing $\mathbf{\Sigma }={\mathbf{A}}^{\mathrm{\top }}\mathrm{diag}[\mathbf{v}]\mathbf{A}$

Let $\mathbf{\Sigma }={\mathbf{A}}^{\mathrm{\top }}\mathrm{diag}[\mathbf{v}]\mathbf{A}$, where $\mathbf{A}$ is fixed and $\mathbf{v}\in {\mathbb{R}}^L$ are free parameters. Define $\mathrm{diag}[\cdot ]$ as an operator that constructs a diagonal matrix from a vector, or extracts the main diagonal from a matrix if its argument is a matrix. The gradient of $\mathcal{D}$ in $\mathbf{v}$ is:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\partial }_{\mathbf{X}}\mathcal{D}& ={\partial }_{\mathbf{X}}\frac{1}{2}\{\mathrm{tr}[\mathbf{\Lambda }{\mathbf{A}}^{\mathrm{\top }}\mathrm{diag}[\mathbf{v}]\mathbf{A}]-\ln |{\mathbf{A}}^{\mathrm{\top }}\mathrm{diag}[\mathbf{v}]\mathbf{A}|\}\\ & =\frac{1}{2}\{\mathrm{diag}[\mathbf{A}\mathbf{\Lambda }{\mathbf{A}}^{\mathrm{\top }}]-\frac{1}{\mathbf{v}}\}\end{array}\end{array}$$

The hessian in $\mathbf{v}$ is a matrix in this case:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathrm{H}}_{\mathbf{v}}\mathcal{D}& =\frac{1}{2}\mathrm{diag}\left[\frac{1}{{\mathbf{v}}^2}\right].\end{array}\end{array}$$

This parameterization is useful for spatiotemporal inference problems, where the matrix $\mathbf{A}$ represents a fixed convolution which can be evaluated using the Fast Fourier Transform (FFT).

---

Inverse-diagonal approximation

Let ${\mathbf{\Sigma }}^{-1}=\mathbf{\Lambda }+\mathrm{diag}\left[\mathbf{p}\right]$. To obtain the gradient in $\mathbf{p}$, combine the derivatives ${\partial }_{\mathbf{\Sigma }}\mathcal{D}$ (Eq. $\text{(2)}$) and ${\partial }_{\mathbf{p}}\mathbf{\Sigma }$ using the chain rule. If $\mathcal{f}(\mathbf{\Sigma })$ is a function of $\mathbf{\Sigma }$, and $\mathbf{\Sigma }({\theta }_i)$ is a function of a parameter ${\theta }_i$, then the chain rule is (The Matrix Cookbook; Eq. 136):

$$\begin{array}{}\text{(15)}& \begin{array}{r}{\partial }_{{\theta }_i}\mathcal{f}=⟨{\partial }_{\mathbf{\Sigma }}\mathcal{f},{\partial }_{{\theta }_i}\mathbf{\Sigma }⟩=\sum _{kl}({\partial }_{{\mathbf{\Sigma }}_{kl}}\mathcal{f})({\partial }_{{\theta }_i}{\mathbf{\Sigma }}_{kl})\end{array}\end{array}$$

From $\text{(2)}$ we have ${\partial }_{\mathbf{\Sigma }}\mathcal{D}=\frac{1}{2}(\mathbf{\Lambda }-{\mathbf{\Sigma }}^{-1})$; Since ${\mathbf{\Sigma }}^{-1}=\mathbf{\Lambda }+\mathrm{diag}\left[\mathbf{p}\right]$, this simplifies to:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\partial }_{\mathbf{\Sigma }}\mathcal{D}& =\frac{1}{2}(\mathbf{\Lambda }-{\mathbf{\Sigma }}^{-1})\\ & =\frac{1}{2}(\mathbf{\Lambda }-\mathbf{\Lambda }-\mathrm{diag}\left[\mathbf{p}\right])\\ & =-\frac{1}{2}\mathrm{diag}\left[\mathbf{p}\right]\end{array}\end{array}$$

We also need ${\partial }_{{\mathbf{p}}_i}\mathbf{\Sigma }$. Let $\mathbf{Y}={\mathbf{\Sigma }}^{-1}$. The derivative $\partial {\mathbf{Y}}^{-1}$ is given as identity (59) in The Matrix Cookbook as $\partial {\mathbf{Y}}^{-1}=-{\mathbf{Y}}^{-1}(\partial \mathbf{Y}){\mathbf{Y}}^{-1}$. Using this, we can obtain ${\partial }_{{\mathbf{p}}_i}\mathbf{\Sigma }$:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\partial }_{{\mathbf{p}}_i}\mathbf{\Sigma }& ={\partial }_{{\mathbf{p}}_i}{\mathbf{Y}}^{-1}=-{\mathbf{Y}}^{-1}\left({\partial }_{{\mathbf{p}}_i}\mathbf{Y}\right){\mathbf{Y}}^{-1}=-\mathbf{\Sigma }\left({\partial }_{{\mathbf{p}}_i}{\mathbf{\Sigma }}^{-1}\right)\mathbf{\Sigma }\\ & =-\mathbf{\Sigma }{\partial }_{{\mathbf{p}}_i}[\mathbf{\Lambda }+\mathrm{diag}[{\mathbf{p}}_i]]\mathbf{\Sigma }=-\mathbf{\Sigma }{\mathbf{J}}_{ii}\mathbf{\Sigma }\\ & =-{\mathit{\sigma }}_i{\mathit{\sigma }}_i^{\mathrm{\top }}\end{array}\end{array}$$

where ${\mathit{\sigma }}_i$ is the $i^{\text{th}}$ row of $\mathbf{\Sigma }$ and ${\mathbf{J}}_{ii}$ is a matrix which is zero everywhere, except for at index $(i,i)$, where it is $1$.

Applying $\text{(15)}$ to $\text{(16)}$ and $\text{(17)}$ for a particular element ${\mathbf{p}}_i$ gives:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\partial }_{{\mathbf{p}}_i}\mathcal{D}& =\sum _{kl}[{\partial }_{{\mathbf{\Sigma }}_{kl}}\mathcal{D}][{\partial }_{{\mathbf{p}}_i}{\mathbf{\Sigma }}_{kl}]=\sum _{kl}{\{-\frac{1}{2}\mathrm{diag}\left[\mathbf{p}\right]\}}_{kl}{\{-{\mathit{\sigma }}_i{\mathit{\sigma }}_i^{\mathrm{\top }}\}}_{kl}\\ & =\frac{1}{2}\sum _{kl}{\delta }_{k=l}{\mathbf{p}}_k{\mathit{\sigma }}_{ik}{\mathit{\sigma }}_{il}=\frac{1}{2}\sum _k{\mathbf{p}}_k{\mathit{\sigma }}_{ik}{\mathit{\sigma }}_{ik}=\frac{1}{2}\sum _k{\mathbf{p}}_k{\mathit{\sigma }}_{ik}^2\\ & =\frac{1}{2}\mathbf{p}{\mathit{\sigma }}_i^{\circ 2}\end{array}\end{array}$$

where $(\cdot {)}^{\circ 2}$ denotes the element-wise square of a vector or matrix. In matrix notation, this is:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\partial }_{\mathbf{p}}\mathcal{D}& =\frac{1}{2}\mathbf{p}{\mathbf{\Sigma }}^{\circ 2}=\frac{1}{2}\mathrm{diag}\left[\mathbf{\Sigma }\mathrm{diag}\left[\mathbf{p}\right]\mathbf{\Sigma }\right],\end{array}\end{array}$$

The Hessian-vector product is cumbersome, since each term in the expression $\mathbf{\Sigma }\left(\mathrm{diag}\left[\mathbf{p}\right]\right)\mathbf{\Sigma }$ depends on $\mathbf{p}$. In the case of the log-linear Poisson GLM, the gradient $\text{(???)}$ simplifies further and optimization becomes tractable. We will explore this further in later notes.

This parameterization resembles the closed-form covariance update for a linear, Gaussian model, where $1/\mathbf{p}$ is a vector of measurement noise variances. It is also a useful parameterization for variational Bayesian solutions for non-conjugate Generalized Linear Models (GLMs), where $\mathbf{p}$ becomes a free parameter to be estimated.

---

Optimizing $\mathbf{\Sigma }={\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{\mathbf{Q}}^{\mathrm{\top }}\mathbf{F}$

Let $\mathbf{\Sigma }={\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{\mathbf{Q}}^{\mathrm{\top }}\mathbf{F}$, where $\mathbf{Q}\in {\mathbb{R}}^{K\times K}$; $K<L$ is the free parameter and $\mathbf{F}\in {\mathbb{R}}^{K\times L}$ is a fixed transformation. If $\mathbf{Q}$ is a lower-triangular matrix, then this approximation involves optimizing $K(K+1)/2$ parameters.

Since the trace is invariant under cyclic permutation, $\mathrm{tr}\left[\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{\mathbf{Q}}^{\mathrm{\top }}\mathbf{F}\right]=\mathrm{tr}\left[\mathbf{F}\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{\mathbf{Q}}^{\mathrm{\top }}\right]$. The derivatives have the same form as $\text{(12)}$ with $\tilde{\mathbf{\Lambda }}=\mathbf{F}\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}$:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\partial }_{\mathbf{Q}}\mathcal{D}& =\tilde{\mathbf{\Lambda }}\mathbf{Q}-{\mathbf{Q}}^{-\mathrm{\top }}\\ & =\mathbf{F}\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}-{\mathbf{Q}}^{-\mathrm{\top }}\\ {\partial }_{\mathbf{Q}}⟨{\partial }_{\mathbf{Q}}\mathcal{D},\mathbf{M}⟩& =\tilde{\mathbf{\Lambda }}{\mathbf{M}}^{\mathrm{\top }}+{\mathbf{Q}}^{-\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{\mathbf{Q}}^{-\mathrm{\top }}\\ & =\mathbf{F}\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}+{\mathbf{Q}}^{-\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{\mathbf{Q}}^{-\mathrm{\top }}\end{array}\end{array}$$

This form is convenient for spatiotemporal inference problems that are sparse in frequency space. In this application, $\mathbf{F}$ corresponds a (unitary) Fourier transform with all by $K$ of the resulting frequency components discarded. The product of $\mathbf{F}$ with a vector $\mathbf{v}$ can be computed in $\mathcal{O}[L\log (L)]$ time using the Fast Fourier Transform (FFT). Alternatively, if $K\le \mathcal{O}(\log (L))$, it is faster to simply multiply $\mathbf{F}\mathbf{v}$ directly. Furthermore, if $\mathbf{F}$ is semi-orthogonal ($\mathbf{F}{\mathbf{F}}^{\mathrm{\top }}=\mathbf{I}$), then calculation of ${\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}$ can be re-used (for example $\mathrm{diag}[\mathbf{\Sigma }]=[({\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{)}^{\circ 2}{]}^{\mathrm{\top }}\mathbf{1}$).

---

Conclusion

These notes provide the gradients and Hessian-vector products for four simplified parameterizations of the posterior covariance matrix for variational Gaussian process inference. If combined with the gradients and Hessian-vector products for the expected log-likelihood, these expressions can be used with Krylov-subspace solvers to compute the Newton-Raphson update to optimize $\mathbf{\Sigma }$.

We evaluated the following parameterizations for $\mathbf{\Sigma }$:

• $\mathbf{\Sigma }$:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}\partial & =\frac{1}{2}(\mathbf{\Lambda }-{\mathbf{\Sigma }}^{-1})\\ \partial ⟨\partial ,\mathbf{M}⟩& =\frac{1}{2}{\mathbf{\Sigma }}^{-1}{\mathbf{M}}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\end{array}\end{array}$$

• $\mathbf{\Sigma }\approx \mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\partial & =\mathbf{\Lambda }\mathbf{X}-{{\mathbf{X}}^{+}}^{\mathrm{\top }}.\\ \partial ⟨\partial ,\mathbf{M}⟩& =\mathbf{\Lambda }{\mathbf{M}}^{\mathrm{\top }}+{{\mathbf{X}}^{+}}^{\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{{\mathbf{X}}^{+}}^{\mathrm{\top }}-(\mathbf{I}-{{\mathbf{X}}^{+}}^{\mathrm{\top }}{\mathbf{X}}^{\mathrm{\top }})\mathbf{M}{\mathbf{X}}^{+}{{\mathbf{X}}^{+}}^{\mathrm{\top }}\end{array}\end{array}$$

• $\mathbf{\Sigma }\approx {\mathbf{A}}^{\mathrm{\top }}\mathrm{diag}[\mathbf{v}]\mathbf{A}$:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}\partial & =\frac{1}{2}\{\mathrm{diag}[\mathbf{A}\mathbf{\Lambda }{\mathbf{A}}^{\mathrm{\top }}]-\frac{1}{\mathbf{v}}\}\\ \partial ⟨\partial ,\mathbf{u}⟩& =\frac{1}{2}{\left[\frac{1}{{\mathbf{v}}^2}\right]}^{\mathrm{\top }}\mathbf{u}\end{array}\end{array}$$

• $\mathbf{\Sigma }\approx [\mathbf{\Lambda }+\mathrm{diag}[\mathbf{p}]{]}^{-1}$:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\partial & =\frac{1}{2}\mathbf{p}{\mathbf{\Sigma }}^{\circ 2}=\frac{1}{2}\mathrm{diag}\left[\mathbf{\Sigma }\mathrm{diag}\left[\mathbf{p}\right]\mathbf{\Sigma }\right],\end{array}\end{array}$$

• ${\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}{\mathbf{Q}}^{\mathrm{\top }}\mathbf{F}$:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}\partial & =\mathbf{F}\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}\mathbf{Q}-{\mathbf{Q}}^{-\mathrm{\top }}\\ \partial ⟨\partial ,\mathbf{M}⟩& =\mathbf{F}\mathbf{\Lambda }{\mathbf{F}}^{\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}+{\mathbf{Q}}^{-\mathrm{\top }}{\mathbf{M}}^{\mathrm{\top }}{\mathbf{Q}}^{-\mathrm{\top }}\end{array}\end{array}$$

In future notes, we will consider the full derivatives required for variational latent Gaussian-process inference for the Poisson and probit generalized linear models.