2.1 Static models

Inference in this model is a simple application of Bayes’ rule. We start by noting that the posterior is necessarily another categorical distribution (i.e., a toss of a die); the only question is what the probabilities of each category (each side of the die) are.

We begin somewhat pedantically with Bayes’s theorem as it is expressed in Eq. 2.1, but for a categorical random variable represented one-hot. To emphasize that we are making inferences from actual observations, we write this in terms of an observation from the data distribution, $\bm{y}$:

i.e., the $j^{\text{th}}$ output of the softmax function^††margin: softmax function, or more properly the soft argmax , where $\hat{z}_{k}=\log\mathopen{}\mathclose{{}\left({\hat{p}\mathopen{}\mathclose{{}\left(\bm{y}\middle|{\hat{X}}_{k}=1;\bm{\theta}}\right)}\pi_{k}}\right)$. For the vector of all possible categories, ${\bm{\hat{X}}}$, the posterior is therefore

The derivation so far is general to any mixture model. For the GMM’s Gaussian emissions, Eq. 2.4, $\hat{z}_{k}$ becomes

where the constant $c$ is irrelevant since it occurs in all the terms and therefore factors out of the softmax. This quantity, $\hat{z}_{k}$, is (up to an additive constant) the log of the posterior probability of class $k$. What does it mean?

One way to interpret Eq. 2.8 is (mentally) to fix $\hat{z}_{k}$ and see what the contours of constant (log-)probability look like. Or again, one can set $\hat{z}_{k}\stackrel{{\scriptstyle\text{set}}}{{=}}\hat{z}_{j}$ for some $j\neq k$, in which case the resulting expression (in $\bm{y}$) is the boundary between classes $k$ and $j$. In either case, these expressions are second-order polynomials in $\bm{y}$, a property inherited from the normal distribution. This is illustrated in Fig. LABEL:fig:XXX. If we sought fancier (higher-order) boundaries between classes, we would need a distribution with more “shape” parameters. ^††margin: insert figure of quadratic class boundaries here; also a subfig of linear boundaries with the projection…..

Moving toward simpler, rather than more complex, contours, consider now the special case of constant covariance matrices across classes. Then the terms $-\frac{1}{2}\log|\mathbf{{\Sigma}}|$ and $-\frac{1}{2}\bm{y}^{\text{T}}\mathbf{{\Sigma}}^{-1}\bm{y}$ could be factored out of all the elements of $\bm{\hat{z}}$, and cancelled (like $c$). That leaves

a linear function of $\bm{y}$. Setting two such expressions equal to each other for different classes $k$ and $j$ (and applying a little bit of algebra), we see that two classes are equiprobable when

Class $k$ is more probable than class $j$ when the left-hand side is larger than the right, and vice versa. (Notice, though, that in neither case is one of these classes guaranteed to be most likely, since we are here ignoring all other classes).

We can even ignore the covariance matrix if we are willing to work in the whitened space (i.e., absorb a factor of $\mathbf{{\Sigma}}^{{-1/2}}$ into $\bm{y}$ and $\bm{\mu}$). Then Eq. 2.9 becomes transparent. To decide the relative probability of a (whitened) observation belonging to class $k$ or class $j$, we first measure how far it is from the midpoint between the class centers, $\frac{1}{2}\mathopen{}\mathclose{{}\left(\bm{\mu}_{k}+\bm{\mu}_{j}}\right)$. We then project this (whitened) displacement onto the vector connecting their (whitened) centers. Here the vector runs from $j$ to $k$, so the probability of class $k$ increases with (positive) distance along the projection. But the point of equality on the projection is not zero displacement, because one class may be a priori more probable than the other. The term on the right-hand side accounts for this.

Acquiring the parameters of the Gaussian mixture model with identical covariances (Eq. 2.1.1), and then using it to classify points according to the posterior distribution (Eq. 2.7), is known as linear discriminant analysis (LDA)^††margin: linear discriminant analysis (LDA) , a nearly century-old method that goes back to Fisher. When the covariances differ, as in the more generic GMM of Eq. 2.8, and the discrimating boundary is quadratic, the method is known, appropriately, as QDA^††margin: QDA . If we generalize along a different avenue, however, allowing the emissions to be members of the same but otherwise arbitrary exponential family with identical dispersions, the boundary remains linear [21]. (We will discuss exponential families in detail in Chapter 5.) In all of these cases, estimating the parameters is conceptually straightforward: we separate the data into classes and compute the sample statistics—e.g., means and covariances for the GMM—within each. We discuss this parameter estimation more formally in Chapter 6.

As noted at the outset, our evaluation of the model often requires ${\hat{p}\mathopen{}\mathclose{{}\left(\bm{y};\bm{\theta}}\right)}$. It is easily read off the denominator in our derivation of a mixture model’s posterior distribution:

For any easily evaluable emission distributions, this is computed painlessly: the probability of the datum $\bm{y}$ under mixture component $k$, scaled by the prior probability of component $k$, and then summed over all $k$.

Comparing Eqs. 2.7 and 2.2, we see that the posterior distribution over ${\bm{\hat{X}}}$ is in the same family as the prior distribution—to wit, the family of categorical distributions. We need not even have carried out any derivation to see this: the support of ${\bm{\hat{X}}}$ is the set of length-${K}$ one-hot vectors (or alternatively, the support of ${\hat{X}}$ is the set of integers from 1 to ${K}$), any distribution over which can be described as categorical. Now recall the definition of a conjugate prior^††margin: conjugate prior : the prior distribution of a parameter is conjugate to the likelihood of that parameter when the resulting posterior distribution over the parameter is in the same family as the prior. So we interpret the “source” variable ${\bm{\hat{X}}}$ as a set of parameters, albeit with some trepidation, since the parameters of exponential families are always continuous rather than discrete (although see [47]). Still, comparing Eqs. 2.3 and 2.4, we see that the “prior” distribution mimics the “likelihood” in the usual way, with the sufficient statistics being the vector of log probabilities. Thus encouraged, we put the emission in exponential-family form,

and the prior distribution into the standard form for conjugate priors:

The lack of a log-partition function should give us pause—it says that the partition function is unity for all values of the “parameters”—but we could nevertheless interpret $\log\bm{\pi}$ as the (first) standard hyperparameter of the conjugate prior. To compute the posterior distribution, this gets updated in the standard way for conjugate priors, namely, by adding the sufficient statistics: $\log\pi_{k}+\log\mathopen{}\mathclose{{}\left({\hat{p}\mathopen{}\mathclose{{}\left(\bm{y}\middle|\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};\bm{\theta}}\right)}}\right)$ for all $k$.

In fine, the categorical distribution is in some sense “conjugate” to any likelihood function of log probability. This is one way of understanding why mixture models are amenable to Bayesian analyses….

To compute the posterior, we of course use Bayes’s theorem, Eq. 2.1, but this time we take advantage of a property of the prior and emission distributions alluded to at the start of the chapter. As this example nicely illustrates, it will not be necessary to compute the normalizer explicity. That is because we will recognize the form of the distribution—in this case, another Gaussian—even without the normalizer. This property follows from the fact that both the prior and the emission are exponentiated second-order polynomials 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}}$. Multiplying two such functions together yields a third exponentiated second-order polynomial 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}}$ (adding second-order polynomials can’t increase their order), which must be a third normal distribution.

There are circumstances under which this is really all we need to know; but what if we want the normalizer explicity? We still don’t have to take the integral, because the normalizer for a Gaussian can be expressed simply as a function³³ 3 We use $\tau=2\pi$ throughout this book [11]. of the covariance matrix: $Z=\tau^{{K}/2}|\mathbf{{\Sigma}}_{\hat{x}|\hat{y}}|^{1/2}$. (Another way to think about this is that someone had to take the integral once, but we don’t have to do it again.) Naturally, to compute this we do need an expression for the covariance matrix, so it won’t do to leave the second-order polynomial in any form whatever; we need to re-arrange it into a quadratic form, with the inverse covariance matrix in the middle. That is the only algebraic work that actually needs to be done in applying Bayes’s theorem:

where we have defined

and

These posterior cumulants should be somewhat intuitive. The posterior mean is a convex combination of the information from the emission and the prior distributions. More precisely, it is a convex combination of the prior mean, $\bm{\mu}_{\hat{x}}$, and the centered observation that has been transformed into the source space, $\mathbf{{C}}^{\dagger}(\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}}-\bm{c})$. Here $\mathbf{{C}}^{\dagger}$ is a pseudo-inverse of $\mathbf{{C}}$.⁴⁴ 4 When ${\bm{\hat{X}}}$ is higher-dimensional than ${\bm{\hat{Y}}}$, this could be any right inverse. When ${\bm{\hat{Y}}}$ is higher-dimensional than ${\bm{\hat{X}}}$, the pseudo-inverse should be intepreted in the least-squares sense. The weights are the (normalized) precisions of these two distributions, although the precision of the emission about ${\bm{\hat{X}}}$ must be transformed into the source space first with $\mathbf{{C}}^{\text{T}}$. The posterior precision, for its part, is just the sum of these two (unnormalized) precisions. As we shall see shortly, it will be convenient to have a name for the normalized precision of the emission:

So much for the posterior. What about ${\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}}{};\bm{\theta}}\right)}$? We suspect it will be yet another normal distribution, although we need to make sure. Given this suspicion, and since a normal distribution is completely characterized by its first two cumulants, we use the laws of total expectation,

and of total covariance,

The higher cumulants of ${\bm{\hat{Y}}}$ are necessarily zero because each term in the law of total cumulance invokes a higher-order (than two) cumulant of either ${\bm{\hat{X}}}$ or ${\bm{\hat{Y}}}|{\bm{\hat{X}}}$, and these are all zero. So the marginal distribution is indeed again normal:

Between this marginal distribution, Eqs. 2.18, 2.16, and 2.17, and the posterior, Eqs. 2.12, 2.13, and 2.14, we have “reversed the arrow” in Fig. 2.2A, providing the alternative parameterization in the equivalent graphical model of Fig. 2.2B.

There is a third convenient way to characterize the joint distribution, and that is to give the distribution of the vector of ${\bm{\hat{X}}}$ and ${\bm{\hat{Y}}}$ concatenated together. Unsurprisingly, this vector is also normally distributed, and indeed we have almost all the pieces for this already in hand.⁵⁵ 5 Notice that the joint distribution of ${\bm{\hat{X}}}$ and ${\bm{\hat{Y}}}$ does not have a convenient form for the Gaussian mixture model. The mean is given by the concatenation of $\bm{\mu}_{\hat{x}}$ and $\bm{\mu}_{\hat{y}}$. The covariance matrix has $\mathbf{{\Sigma}}_{\hat{x}}$ and $\mathbf{{\Sigma}}_{\hat{y}}$ in its upper-left and lower-right blocks. What remains is the upper-right (or its transpose in the lower-left) block, $\mathbf{\Sigma}_{\hat{x},\hat{y}}$, i.e. the covariance between ${\bm{\hat{X}}}$ and ${\bm{\hat{Y}}}$.

From Eq. 2.11, we can write ${\bm{\hat{Y}}}=\mathbf{{C}}{\bm{\hat{X}}}+\bm{c}+{\bm{\hat{Z}}}$, where ${\bm{\hat{Z}}}\sim\mathcal{N}\mathopen{}\mathclose{{}\left(\bm{0},\>\mathbf{{\Sigma}}_{\hat{y}|\hat{x}}}\right)$ is independent of ${\bm{\hat{X}}}$. Hence:

This covariance will turn out to be useful later when we take on the learning problem.

For the models described by Eq. 2.11, it is perhaps unsurprising that the posterior is again normal (Eq. 2.12). After all, the normal distribution is well known to be the conjugate prior for the mean of another normal distribution (it induces a normal posterior distribution over that mean). However, in the family of jointly Gaussian models we are now considering (Eq. 2.11), the mean of the emission distribution is an affine function of, rather than identical with, the source random variable. Evidently, conjugacy holds for this more generic class as well as the classical result for a normally distributed mean, which can be construed as a special case.

We consider now the special case where the underlying normal source variable ${\bm{\hat{X}}}$ is never observed—it is latent. For example, suppose that ${\bm{Y}}$ is the (random) vector of performances on a battery of intelligence tests of various kinds. Each sample ${\bm{y}}_{n}$ corresponds to individual $n$’s performance on all $L$ exams. We hypothesize that this performance can be explained by a smaller, $K<L$, number of latent intelligence “factors,” ${\bm{\hat{X}}}$, that are normally distributed across individuals. But we can’t observe them directly; we only observe test scores, which are (ideally) an affine function of these latent factors. Finally, we allow for errors in the model, errors in our measurements, unmodelled enviromental inputs, etc. by allowing our test results to have been corrupted by Gaussian noise.

This sounds like the jointly Gaussian model just discussed, but the fact that we never observe ${\bm{\hat{X}}}$ obliges us to restrict the degress of freedom of the model, as discussed at the outset of this chapter. In particular, we have just seen that the marginal distribution of observed data is normal, Eq. 2.18, and in this factor analysis, this is the only distribution we ever see. Eq. 2.16 shows that, as far as this distribution is concerned, any change in the mean of the prior distribution, $\bm{\mu}_{\hat{x}}$, could equally well be attributed to a change in the emission offset, $\bm{c}$. So without loss of generality, we can let $\bm{\mu}_{\hat{x}}=0$.

A similar consideration applies to the prior covariance. From Eq. 2.17, we see that any non-isotropic covariance could simply be absorbed into $\mathbf{{C}}$. Indeed, this makes the intepretation of the model more intuitive: we are asking for the independent factors that together explain test performances. Any covariance across different kinds of intelligence tests is due to their having similar “loading” onto these independent factors.

[[We remove one more degree of freedom. In general, $\mathbf{{\Sigma}}_{\hat{y}|\hat{x}}$ could be any positive definite matrix, but assuming it to be diagonal, $\mathbf{\Lambda}$, yields a useful interpretation…..]].

So our complete model is:

In what follows, we shall assume that one has the correct model, and needs to solve only the inference problem. We shall return to learning in Section 7.3.

We take the marginal probability first this time. Applying Eqs. 2.16, 2.17, and 2.18 to the special case of Eq. 2.20, we have

Inference in a factor analyzer is just a special case of inference in the jointly Gaussian model introduced at the beginning of this section. That is, it requires only the computation of the posterior cumulants, Eqs. 2.13 and 2.14, under the distributions that define the model, Eq. 2.20:

Anticipating the computational costs of this, and more complicated, inference problems with Gaussian random variables, we note with some dismay the matrix inversions, which cost at least $\mathcal{O}(M^{2.373})$ for $M\times M$ matrices. In factor analysis, as we have noted, the observations $\bm{y}$ are higher-dimensional than the latent variables ${\bm{\hat{X}}}$. It is fortunate, then, that the only matrix inversion carried out in the observation space is of $\mathbf{\Lambda}$, which is just $\mathcal{O}(M)$ (invert the entries on the diagonal). But what about models where the “source space” is larger than the “emission space”? What if the emission covariance is not diagonal? And how can we limit the total number of inversions when the inference step has to be applied repeatedly?

To change the space in which the posterior covariance is computed, we can apply the Woodbury matrix-inversion lemma, Eq. B.14, to Eq. 2.14:

This moves the inversion from source (${\bm{\hat{X}}}$) to emission (${\bm{\hat{Y}}}$) space. (Notice, incidentally, that the inverted quantity is precisely the marginal covariance, $\mathbf{{\Sigma}}_{\hat{y}}$.)

The expression for the posterior mean, Eq. 2.13, likewise invokes an inversion in source space ($\mathbf{{\Sigma}}^{-1}_{\hat{x}}$), as well as one in the emission space ($\mathbf{{\Sigma}}^{-1}_{\hat{y}|\hat{x}}$). In factor analysis, the first inversion was not required because the prior mean, $\bm{\mu}_{\hat{x}}$, was assumed to be zero. Here we eliminate it by substituting the “Woodburied” posterior covariance, Eq. 2.23, into Eq. 2.13. We also make use of the definition of the emission precision, Eq. 2.15:

Intuitively, the posterior mean is the prior mean, adjusted by a “correction factor.” The correction measures the difference between the emission as observed ($\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}}$) and as predicted by the prior mean ($\mathbf{{C}}\bm{\mu}_{\hat{x}}+\bm{c}$), weighted by the precision of the emission ($\mathbf{K}$).

Despite appearances, matrix inversions are still required in Eq. 2.24; they occur in the calculation of $\mathbf{K}$. Consulting its definition, Eq. 2.15, we see one inversion of the emission covariance, and another in the computation of the posterior covariance, albeit in the space of emissions via Eq. 2.23. With a little work, one of these inversions can be eliminated:

Thus in fact only one matrix need be inverted in computing the entire posterior distribution, namely $\mathbf{{C}}\mathbf{{\Sigma}}_{\hat{x}}\mathbf{{C}}^{\text{T}}+\mathbf{{\Sigma}}_{\hat{y}|\hat{x}}$—which, by Eq. 2.17, is just the marginal covariance, $\mathbf{{\Sigma}}_{\hat{y}}$. So we can also write

To make the complete calculation more perspicuous, we collect here the marginal cumulants from Eqs. 2.16 and 2.17:

along with this alternative (“Woodburied”) computation of the posterior cumulants:

where on the last line we have rewritten the posterior covariance yet again using the expression for the emission precision.

Both types of model are problematic for complex distributions, like that of natural images. Gaussian mixture models can indeed approximate any marginal density over ${\bm{Y}}$ (Eq. 2.10) to arbitrary accuracy [49]: intuitively, a distribution of any shape can be constructed by placing little blobs of probability mass at various locations. This is essential for the distribution of images, where the presence of (e.g.) lines makes the marginal distribution of (say) triplets of adjacent pixels inexplicable by their mean values or even pairwise correlations alone (Fig. LABEL:fig:pixelDistribution). But each such bit of structure must be “explained” with a separate basis vector! There is no way to explain a letter “X” as a superimposition of two lines; rather, the “X” and the two lines into which it might be decomposed would each require its own basis vector. More generally, every image (or image patch) that could not be expressed as a Gaussian disturbance of some other image would require another column in $\mathbf{{C}}$. Thus the ability to model structure comes at the price of a very large number of mixture components. It must be emphasized that we believe this price to be too high precisely because we believe the structure in images to be “componential,” explicable in terms of combinations of simpler basis elements.

A factor analyzer, in contrast, can avail itself of all basis vectors at once in attempting to explain or generate an observation, and consequently can in theory recognize an image in terms of its constituent parts. By the same token, it needs no more hidden units than there are pixels in the images: $\mathbf{{C}}$ can be square, i.e. a complete basis. However, the marginal distribution of observations under a factor analyzer is just another normal distribution (recall Eq. 2.21)—a fact that looked helpful from the perspective of calculating with the model, but is fatal for modeling natural images. This is a consequence of the higher-order structure in images alluded to above.

The existence of such higher-order structure as seen in Fig. LABEL:fig:pixelDistribution is not perhaps obvious a priori. After all, the size of the covariance matrix for 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}}{};\bm{\theta}}\right)}$, scales quadratically in the number of pixels—a lot of representational power to work with! For example, the covariance matrix of even modestly sized, (256 $\times$ 256)-pixel, images has $256^{2}(256^{2}-1)/2\approx$ $2\text{\times}{10}^{9}$ independent entries. It would seem that we could represent any distribution with this many parameters. But in fact, by taking certain symmetries into account, we reduce the number of covariance entries to a linear function of the number of pixels. The statistics of natural images are, for the most part⁶⁶ 6 Many natural images, as seen by land-dwelling animals, have sky in the upper portion and ground below, which does introduce statistical regularities that depend on absolute position. Overrepresentation of vertical and (especially) horizontal lines would likewise introduce assymmetries. , stationary: the correlation between any two pairs of pixels depends on the distance between them but not their absolute locations. In our example, there are only $256^{2}$ distances—from, say, the upper left pixel to any pixel in the image, including itself. So the covariance cannot encode so much information after all.

Covariance matrices arising from shift-invariant data will in general have this “symmetric Toeplitz” structure. More specifically, when the shift-invariant data are two-dimensional, like images, the covariance matrix will be block-Toeplitz with Toeplitz blocks [10]. Furthermore, if we assume that the covariance drops to zero at some distance within the image, and we ignore edge effects, the covariance matrix can be approximated as block-circulant with circulant blocks, in which case its eigenvectors are the (two-dimensional) discrete Fourier modes, and its eigenvalues constitute the power spectrum [10].⁷⁷ 7 Since the eigenvectors of the covariance matrix are the principal components of the data (see also Section 7.3), this likewise implies that the Fourier modes are the principal components of shift-invariant data. In natural images, the power spectrum appears to fall off with frequency $f$ at a rate of very nearly $1/f^{2}$ [44].

Thus, natural images can be decorrelated by projection onto the Fourier modes, and then “whitened” by normalizing by the amplitude spectrum (the square roots of the power spectrum). Yet it is easily seen Fig. LABEL:fig:SimoncelliOlshausenFig5 that this process leaves most of the structure intact [44]. Or again, coming from the other direction, restoring the (approximately) $1/f^{2}$ power spectrum to white noise does not produce a natural image (Fig. LABEL:fig:SimoncelliOlshausenFig5). Factor analyzers thus have no means of capturing the structure in natural images that lies beyond the mean and the power spectrum.

We have established, then, that neither factor analyzers nor GMMs are good models for natural images, albeit for roughly opposite reasons. Furthermore, the inadequacy of the multivariate Gaussian as a marginal distribution for natural images is relevant to a wider range of models than the factor analyzer. When the source variables are densely, even though non-normally, distributed, 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}}{};\bm{\theta}}\right)}$ may very well end up approximately normal anyway, for central-limit-like reasons: the marginal distribution is an average of many independent, “random” basis vectors. Averaging can wash out higher-order structure in these vectors. This is (I take it) another way of expressing Hinton’s conjecture that the most useful codes will be componential, in order to represent the properties—pose, deformation, etc.—of the objects in the image; but also non-dense, in order to accommodate multiple objects in the same image [15].

Between the “extremes” of the GMM and the factor analyzer lie models in which the “source” variables are sparsely distributed: some small subset are “on” for any given image. For example, they could be distributed according to a product over Bernoulli distributions with small means. Or again, we could interpret “on” more liberally and assign fat-tailed (leptokurtotic) distributions to the source variables, so that they usually take on values close to zero but can take on extreme probabilities more often than normally distributed variables. Common choices are the Laplace and Cauchy distributions [35, 29]. In this case it is convenient to write the prior distribution in the form of a product of Boltzmann distributions,

with the energy functions $E_{k}$ chosen so as to enforce some form of sparseness. For example,

The last is not a standard distribution, but for large $\beta$ it approximates the Laplace distribution while retaining differentiability at 0 [30]; see Fig. LABEL:fig:sparseDistributions. It will serve as our default example in the sparse-coding model that follows. The emissions, for their part, are once again normally distributed about a linear function of the source variables, and for simplicity we have given them identical variance, $1/\lambda$. Eq. 2.28 describes the complete sparse-coding model, a graphical model that is shown in Fig. 2.4A.

Since we are aiming at economy of expression (few descriptors per image), we must concede prodigality in vocabulary (many descriptors) [39]—i.e., an overcomplete set of basis vectors. This is made explicit in Fig. 2.4B, along with the independence statements implied by our choices of prior and emission distributions. A factor of 2 overcompleteness is common [29].

…

There is a reason that factor analyzers and GMMs remain popular despite their limitations: the model posterior and marginal distributions are computable in closed form. This computability is due to the prior distributions being conjugate (or pseudo-conjugate) to the emission likelihoods. When they are not, the marginal and posterior distributions can typically be computed only approximately. When the prior is Laplace, e.g., it is not clear that the posterior can be computed (although see [38]).

Perhaps the simplest approach is to make a Taylor approximation.⁸⁸ 8 Another, more powerful class of techniques for approximate inference is based on (iterative) optimization of the posterior. However, this endows it with the character of learning, and accordingly we shall defer discussing these techniques until Chapter 6. In particular, consider approximating the energy, $E(\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}})$, of a distribution near its mode, $\bm{\hat{x}}^{0}$, with the first three terms of the Taylor series. It is easy to see that, for $\bm{\hat{x}}^{0}$ to be a mode of the distribution, it must satisfy the following conditions:

i.e., at the mode, the energy’s gradient must be zero and its Hessian must be positive definite. For a generic energy, the second-order Taylor approximation at the mode is therefore

The only remaining $\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}}$-dependent term should look familiar: it is the energy of a normal distribution with mean at the mode and inverse covariance equal to the Hessian of the energy (evaluated at the mode). The first term, $E(\bm{\hat{x}}^{0})$, on the other hand, is constant in the argument $\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 can therefore be absorbed into the normalizer. Hence we can approximate the posterior as

Note that this approximation is valid at realizations of ${\bm{\hat{X}}}$ even “moderately” close to $\bm{\hat{x}}^{0}$, since the neglected higher-order terms need only be less than unity in order for the exponentiation to drive their contribution to the distribution down to much smaller values. On the other hand, it is obviously inappropriate for distributions having multiple modes, or over discrete random variables.

So much for approximating the posterior; we now move on to the model marginal. Using the definition of conditional probability with the approximate posterior distribution, Eq. 2.31, and then evaluating at the mode, we see that

This is the approximate model marginal under Laplace’s method.

Let us specialize to an isotropic Gaussian emission and a factorial prior distribution, Eq. 2.28. Start with the Hessian:

where $\text{diag}\mathopen{}\mathclose{{}\left(\bm{a}}\right)$ is a diagonal matrix with the vector $\bm{a}$ as its diagonal. For concreteness, let us specialize further to Laplace-distributed source variables, as in Eq. 2.29. Unfortunately, the curvature of this energy function is undefined at the origin, and (equally bad) zero everywhere else. Lewicki and colleagues [29, 31] therefore suggest approximating this energy—in the Hessian only—with the hyperbolic-cosine energy (see again Eq. 2.29), which approaches the Laplace energy in the limit of large $\beta$. The second derivative of this energy is

Hence, from Eq. 2.31, an approximate posterior distribution for the sparse-coding model is

where $\circ$ is the element-wise product and $\operatorname{sech}(\cdot)$ is supposed to act element-wise. Notice that although this may appear to be independent of the observations $\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{y}$, they enter through the mode, which obviously varies as a 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{y}$.

To specialize the marginal distribution to the Laplace prior, we use Eq. 2.32, substituting the Laplace energies into the joint, but the hyperbolic-cosine energies into the Hessian (to maintain differentiability):

…