8.1 Gaussian recognition distributions and sparse coding

We can now write a more specific expression for the joint relative entropy that will allow us to see how to carry out the M step, as well. Let us simplify even further and treat $\lambda$ and the parameters of the prior as fixed. Then substituting the cross entropy of the joint, Eq. 8.7, into the joint relative entropy for Gaussian recognition distributions, Eq. 8.1, yields

(All “constant” terms have been lumped into a single scalar $c$.) The only remaining model parameter to be optimized is the emission matrix, $\mathbf{{C}}$. We can compute the gradient of the joint relative entropy, Eq. 8.8, with respect to $\mathbf{{C}}$ using some results from the appendix (Section B.1):

Thus, the emission matrix is once again given by some version of the normal equations. In practice, however, if each $\bm{y}$ is (e.g.) an image, and the model is overcomplete, then the matrix inversions will be computationally challenging. Instead, the gradient in Eq. 8.9 (or the natural gradient [29, 31]) can be descended.

To carry out the E step, we apply Eqs. 8.3 and 8.4 to the energy of the generative model’s joint or posterior distribution (they are equivalent up to constant terms), to wit, the bracketed terms in Eq. 8.7. Let us also commit to the source energies specified by Eq. 8.6, i.e., the energy of the Laplace distribution. Then the resulting joint energy,

is quadratic in $\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}$ everywhere except at points where $\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\hat{x}_{k}$ is zero (for any $k$). Therefore the energy-gradient function (of $\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{x}}$) and the expectation in Eq. 8.3 “almost” commute with each other, and the gradient of the joint relative entropy becomes

Thus¹¹ 1 Technically the final line is only a sufficient, not a necessary, condition for the gradient to be zero, but we are assuming the Hessian is locally positive definite., the optimal choice for the mean is the solution to the $L^{1}$-penalized least-squares problem. For “fat” matrices $\mathbf{{C}}$, this is known as basis-pursuit denoising,²² 2 This is also sometimes referred to as the LASSO problem , although typically that term is reserved for “tall” output matrices $\mathbf{{C}}$. a convex optimization for which efficient solvers exist.

That leaves the covariance matrix, given by Eq. 8.4. Again we exploit the fact that the Hessian of the energy is a sub-cubic—indeed, constant—function of $\bm{\hat{x}}$, except at points where any $\hat{x}_{k}$ is zero. Therefore,

The approximation in the final line follows from the approximation to the Laplace energy in Eq. 8.6.

We have arrived at a sparse-coding algorithm very similar to that of Lewicki and colleagues [29, 31], but by a somewhat different route. In those papers, the matrix of basis functions, $\mathbf{{C}}$, is fit by directly descending along the gradient of the marginal relative entropy (MRE), ${\text{H}_{p\hat{p}}{\mathopen{}\mathclose{{}\left[{\bm{Y}}{};\mathbf{{C}}}\right]}}$, rather than descending via the gradient of its upper bound, the joint relative entropy (JRE). The price, of course, is the intractability of the MRE. In the classical papers, this price is paid with Laplace’s method, approximating the model marginal, ${\hat{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{y}}{};\mathbf{{C}}}\right)}$, with the integral $\int_{\bm{\hat{x}}{}}{\hat{p}\mathopen{}\mathclose{{}\left(\bm{\hat{x}}}\right)}{\hat{p}\mathopen{}\mathclose{{}\left(\bm{\hat{x}}\middle|\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\hat{y}}{};\mathbf{{C}}}\right)}\mathop{}\!\mathrm{d}{\bm{\hat{x}}{}}$ in the vicinity of one of the modes, $\bm{\hat{x}}^{0}$, of the model joint (Eq. 2.34). This in some sense removes the latent variables from the problem, although the modes still need to be estimated. Consequently, there is no E step, and no fixed parameters ${\bm{\nu}_{\bm{y}}}$ and ${\mathbf{{P}}_{\bm{y}}}$ (for any $\bm{y}$).

Nevertheless, the Hessian of the joint (or posterior) energy still shows up in the classical derivation’s loss, this time as a byproduct of Laplace’s method, as a measure of the volume of the posterior. This is no coincidence: we have allowed the E step to collapse into a Laplace approximation [36] by using an energy function that is sub-cubic almost everywhere. However, in the classical method, this precision matrix is not fixed during optimization of $\mathbf{{C}}$, but rather is interpreted as a function of $\mathbf{{C}}$ and $\bm{\hat{x}}^{0}$: $f(\mathbf{{C}},\bm{\hat{x}}^{0})\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.} }}=\log\mathopen{}\mathclose{{}\left\lvert{\mathbf{{P}}_{\bm{y}}}\mathopen{}\mathclose{{}\left(\mathbf{{C}},\bm{\hat{x}}^{0}}\right)}\right\rvert$. On the other hand, the trace term in Eq. 8.8, $\lambda\text{tr}\mathopen{}\mathclose{{}\left[\mathbf{{C}}{\mathbf{{P}}^{-1}_{\bm{y}}}\mathbf{{C}}^{\text{T}}}\right]$, never materializes, because the squared reconstruction error is simply evaluated at the mode, rather than being averaged under a recognition distribution.

Conceptually, these two terms encode the contribution of the recognition precision (${\mathbf{{P}}_{\bm{y}}}$) to the two losses, ours and the classical. Evidently, increasing the precision has opposite effects on the two losses, lowering the JRE (via $\lambda\text{tr}\mathopen{}\mathclose{{}\left[\mathbf{{C}}{\mathbf{{P}}^{-1}_{\bm{y}}}\mathbf{{C}}^{\text{T}}}\right]$) but raising the approximate MRE (via $\log\mathopen{}\mathclose{{}\left\lvert{\mathbf{{P}}_{\bm{y}}}\mathopen{}\mathclose{{}\left(\mathbf{{C}},\bm{\hat{x}}^{0}}\right)}\right\rvert$). But increasing the magnitude of the emission matrix $\mathbf{{C}}$ has the same effect on these two terms: It increases the contribution of recognition uncertainty to the average reconstruction error ($\lambda\text{tr}\mathopen{}\mathclose{{}\left[\mathbf{{C}}{\mathbf{{P}}^{-1}_{\bm{y}}}\mathbf{{C}}^{\text{T}}}\right]$), and it increases the recognition precision itself (and therefore $\log\mathopen{}\mathclose{{}\left\lvert{\mathbf{{P}}_{\bm{y}}}\mathopen{}\mathclose{{}\left(\mathbf{{C}},\bm{\hat{x}}^{0}}\right)}\right\rvert$) via Eq. 8.11. Indeed, the partial derivative of $f(\mathbf{{C}},\bm{\hat{x}}^{0})\mathrel{\vbox{\hbox{\scriptsize.}\hbox{\scriptsize.} }}=\log\mathopen{}\mathclose{{}\left\lvert{\mathbf{{P}}_{\bm{y}}}\mathopen{}\mathclose{{}\left(\mathbf{{C}},\bm{\hat{x}}^{0}}\right)}\right\rvert$ with respect to $\mathbf{{C}}$ is precisely the same as the derivative of the trace term $\lambda\text{tr}\mathopen{}\mathclose{{}\left[\mathbf{{C}}{\mathbf{{P}}^{-1}_{\bm{y}}}\mathbf{{C}}^{\text{T}}}\right]$ with respect to $\mathbf{{C}}$ (namely, $2\lambda\mathbf{{C}}{\mathbf{{P}}^{-1}_{\bm{y}}}$).

There is one last discrepancy between the loss gradients. In the classical derivation, the mode $\bm{\hat{x}}^{0}$ is not fixed, either, and in turns depends on $\mathbf{{C}}$. So the total derivative of $\log\mathopen{}\mathclose{{}\left\lvert{\mathbf{{P}}_{\bm{y}}}\mathopen{}\mathclose{{}\left(\mathbf{{C}},\bm{\hat{x}}^{0}(\mathbf{{C}})}\right)}\right\rvert$ generates extra terms. On the other hand, in the original papers, these are derived and then promptly dropped, under the conjecture that curvature unrelated to volume is irrelevant to the optimization.