Continuous Latent Variables

Many high-dimensional data sets lie close to a surface of much lower dimension. This chapter models that structure with a continuous hidden variable: principal component analysis first as a projection, then as the maximum likelihood fit of a Gaussian latent variable model (with EM, missing data, Bayesian model selection and factor analysis), and finally its kernel, non-Gaussian and nonlinear generalisations.

TrungTin Nguyen · Following Bishop (2006), Pattern Recognition and Machine Learning, Chapter 12 · Course home · Mathematical toolkit · Python code for the figures

Ch. 1 Introduction Ch. 2 Probability distributions Ch. 3 Linear regression Ch. 4 Linear classification Ch. 5 Neural networks Ch. 6 Kernel methods Ch. 7 Sparse kernel machines Ch. 8 Graphical models Ch. 9 Mixture models and EM Ch. 10 Approximate inference Ch. 11 Sampling methods Ch. 12 Continuous latent variables Ch. 13 Sequential data Ch. 14 Combining models

Introduction: data near a low-dimensional manifold

In Chapter 9 the hidden variable was a discrete label (which mixture component produced a point). In this chapter it is continuous. The motivating observation (Bishop, pp. 559–560) is that the points of many data sets lie close to a manifold whose dimension is much smaller than that of the data space.

Bishop's example: paste one image of a handwritten digit into a blank 100 × 100 image, at a random position and with a random rotation (Figure 12.1). Every such image is a point of \(\R^{10\,000}\), yet the whole collection is described by three numbers: the horizontal shift, the vertical shift and the angle. The images therefore lie on a three-dimensional manifold. The manifold is curved: as the digit slides across a pixel, that pixel goes from white to black and back to white, which is not a linear function of the shift. The shifts and the angle are latent variables: they produced the images but are not recorded with them. Real handwriting adds a few more degrees of freedom (size, slant, writing style), still far fewer than the number of pixels. The oil flow data (Appendix A) are a second example: each observation has 12 measurements, but within one flow configuration only two quantities vary (the oil and water fractions; the gas fraction is what remains), so the points lie near a noisy two-dimensional manifold made of several pieces, one or more per flow regime.

The departures of real data from such a manifold are treated as noise. This suggests a generative description: draw a point \(\mathbf{z}\) from a distribution on a low-dimensional latent space, map it into the data space, and add noise. The simplest version makes the latent distribution and the noise Gaussian and the map linear (a linear-Gaussian model, Bishop Section 8.1.4). It leads to probabilistic PCA and to factor analysis. The chapter first treats classical, non-probabilistic PCA (Section 12.1), then shows that it is the maximum likelihood solution of a linear-Gaussian latent variable model (Section 12.2), and ends with nonlinear and non-Gaussian generalisations (Sections 12.3 and 12.4).

12.1 Principal Component Analysis

Principal component analysis (PCA), also called the Karhunen–Loève transform, is used for dimensionality reduction, lossy data compression, feature extraction and visualisation (Jolliffe, 2002). It has two classical definitions that lead to the same algorithm. Hotelling (1933) asks for the orthogonal projection onto an \(M\)-dimensional subspace that keeps as much variance of the data as possible; Pearson (1901) asks for the \(M\)-dimensional (affine) subspace that is closest to the data in mean squared distance. Section 12.1.1 solves the first problem, Section 12.1.2 the second, and Result 12.3 shows that the two answers are the same subspace.

12.1.1 Maximum variance formulation

Fix \(M\lt D\) (choosing \(M\) is discussed in Section 12.2.3) and start with \(M=1\). A one-dimensional projection is described by a direction, which we take to be a unit vector \(\mathbf{u}_1\) (only the direction matters). The point \(\mathbf{x}_n\) is projected to the number \(\mathbf{u}_1^\T\mathbf{x}_n\). The projected data have mean \(\mathbf{u}_1^\T\bar{\mathbf{x}}\), where

\[ \bar{\mathbf{x}}=\frac1N\sum_{n=1}^N\mathbf{x}_n , \tag{12.1} \]

and variance

\[ \frac1N\sum_{n=1}^N\big\{\mathbf{u}_1^\T\mathbf{x}_n-\mathbf{u}_1^\T\bar{\mathbf{x}}\big\}^2=\mathbf{u}_1^\T\mathbf{S}\mathbf{u}_1 , \tag{12.2} \] \[ \mathbf{S}=\frac1N\sum_{n=1}^N(\mathbf{x}_n-\bar{\mathbf{x}})(\mathbf{x}_n-\bar{\mathbf{x}})^\T . \tag{12.3} \]

Maximising \(\mathbf{u}_1^\T\mathbf{S}\mathbf{u}_1\) without a constraint is pointless (it grows like \(\lVert\mathbf{u}_1\rVert^2\)), so we impose \(\mathbf{u}_1^\T\mathbf{u}_1=1\) with a Lagrange multiplier \(\lambda_1\) and look for stationary points of

\[ \mathbf{u}_1^\T\mathbf{S}\mathbf{u}_1+\lambda_1\big(1-\mathbf{u}_1^\T\mathbf{u}_1\big). \tag{12.4} \]

They satisfy

\[ \mathbf{S}\mathbf{u}_1=\lambda_1\mathbf{u}_1 , \tag{12.5} \] \[ \mathbf{u}_1^\T\mathbf{S}\mathbf{u}_1=\lambda_1 . \tag{12.6} \]

Further directions are defined one at a time: each new direction maximises the projected variance among all unit vectors orthogonal to the directions already chosen. The book states that the resulting \(M\) directions are the eigenvectors of the \(M\) largest eigenvalues, "easily shown using proof by induction" (Exercise 12.1). Two statements are hidden in this sentence, and both are proved below: the one-at-a-time construction picks the eigenvectors in order, and the \(M\)-dimensional subspace so obtained maximises the total projected variance among all \(M\)-dimensional projections, not only among those built greedily.

Proof

(a), with Lagrange multipliers as the exercise asks. Suppose, as induction hypothesis, that the first \(k\) directions are the eigenvectors \(\mathbf{u}_1,\dots,\mathbf{u}_k\) (true for \(k=1\) by Result 12.1). The next direction \(\mathbf{u}\) maximises \(\mathbf{u}^\T\mathbf{S}\mathbf{u}\) subject to \(\mathbf{u}^\T\mathbf{u}=1\) and \(\mathbf{u}^\T\mathbf{u}_j=0\) for \(j=1,\dots,k\). With multipliers \(\lambda\) and \(\eta_1,\dots,\eta_k\), the Lagrangian is

\[ \mathbf{u}^\T\mathbf{S}\mathbf{u}+\lambda\big(1-\mathbf{u}^\T\mathbf{u}\big)+\sum_{j=1}^k\eta_j\,\mathbf{u}^\T\mathbf{u}_j , \]

and its gradient vanishes when \(2\mathbf{S}\mathbf{u}-2\lambda\mathbf{u}+\sum_j\eta_j\mathbf{u}_j\)\({}=\mathbf{0}\). (The constraint gradients \(2\mathbf{u},\mathbf{u}_1,\dots,\mathbf{u}_k\) are orthogonal and non-zero at every feasible point, hence linearly independent, so the multiplier rule applies.) Multiply on the left by \(\mathbf{u}_m^\T\) for some \(m\le k\). Because \(\mathbf{S}\) is symmetric and \(\mathbf{u}_m\) is an eigenvector, \(\mathbf{u}_m^\T\mathbf{S}\mathbf{u}\)\({}=(\mathbf{S}\mathbf{u}_m)^\T\mathbf{u}\)\({}=\lambda_m\mathbf{u}_m^\T\mathbf{u}\)\({}=0\); also \(\mathbf{u}_m^\T\mathbf{u}=0\) and \(\mathbf{u}_m^\T\mathbf{u}_j\)\({}=\delta_{mj}\). What is left is \(\eta_m=0\). Hence \(\mathbf{S}\mathbf{u}\)\({}=\lambda\mathbf{u}\): the new direction is an eigenvector orthogonal to \(\mathbf{u}_1,\dots,\mathbf{u}_k\), and its projected variance is \(\mathbf{u}^\T\mathbf{S}\mathbf{u}\)\({}=\lambda\). The admissible eigenvalues are \(\lambda_{k+1},\dots,\lambda_D\) (an eigenvector orthogonal to \(\mathbf{u}_1,\dots,\mathbf{u}_k\) lies in \(\operatorname{span}\{\mathbf{u}_{k+1},\dots,\mathbf{u}_D\}\)), so the maximum is \(\lambda_{k+1}\), attained at \(\mathbf{u}_{k+1}\). Directly: writing \(\mathbf{u}\)\({}=\sum_{i\gt k}c_i\mathbf{u}_i\) gives \(\mathbf{u}^\T\mathbf{S}\mathbf{u}\)\({}=\sum_{i\gt k}\lambda_ic_i^2\)\({}\le\lambda_{k+1}\). This completes the induction.

(b). The identity on the left is the computation of (12.2) applied to each column of \(\mathbf{V}\). Write \(\mathbf{S}\)\({}=\sum_i\lambda_i\mathbf{u}_i\mathbf{u}_i^\T\), so that

\[ \Tr(\mathbf{V}^\T\mathbf{S}\mathbf{V})=\sum_{i=1}^D\lambda_i\,\mathbf{u}_i^\T\mathbf{V}\mathbf{V}^\T\mathbf{u}_i=\sum_{i=1}^D\lambda_ip_i,\qquad p_i=\lVert\mathbf{V}^\T\mathbf{u}_i\rVert^2 . \]

\(\mathbf{P}=\mathbf{V}\mathbf{V}^\T\) is the orthogonal projection onto the column space of \(\mathbf{V}\), so \(0\)\({}\le p_i\)\({}=\mathbf{u}_i^\T\mathbf{P}\mathbf{u}_i\)\({}\le\lVert\mathbf{u}_i\rVert^2\)\({}=1\), and \(\sum_ip_i\)\({}=\Tr(\mathbf{P}\sum_i\mathbf{u}_i\mathbf{u}_i^\T)\)\({}=\Tr(\mathbf{P})\)\({}=\Tr(\mathbf{V}^\T\mathbf{V})\)\({}=M\). Using \(\lambda_i\ge\lambda_M\) for \(i\le M\), \(\lambda_i\le\lambda_M\) for \(i\gt M\), and \(p_i\le1\),

\[ \sum_{i=1}^D\lambda_ip_i-\sum_{i=1}^M\lambda_i=\sum_{i\le M}\lambda_i(p_i-1)+\sum_{i\gt M}\lambda_ip_i\le\lambda_M\sum_{i\le M}(p_i-1)+\lambda_M\sum_{i\gt M}p_i=\lambda_M\Big(\sum_ip_i-M\Big)=0 . \]

If the columns of \(\mathbf{V}\) span \(\operatorname{span}\{\mathbf{u}_1,\dots,\mathbf{u}_M\}\), then \(p_i=1\) for \(i\le M\) and \(p_i=0\) otherwise, and equality holds. Conversely, if \(\lambda_M\gt\lambda_{M+1}\), equality in the second step forces \(p_i=0\) for every \(i\gt M\) (and then \(p_i=1\) for \(i\le M\)), that is, the column space of \(\mathbf{V}\) is \(\operatorname{span}\{\mathbf{u}_1,\dots,\mathbf{u}_M\}\). Since the greedy directions of (a) attain the bound, they are also optimal among all \(M\)-dimensional projections. This is the "Ky Fan" form of the Rayleigh–Ritz theorem. ∎

In summary, PCA computes \(\bar{\mathbf{x}}\) and \(\mathbf{S}\) and finds the eigenvectors of the \(M\) largest eigenvalues. A full eigendecomposition of a \(D\times D\) matrix costs \(O(D^3)\); when only \(M\) eigenpairs are needed, iterative methods such as the power method cost \(O(MD^2)\) (Golub and Van Loan, 1996), and the EM algorithm of Section 12.2.2 avoids forming \(\mathbf{S}\) altogether. The power method itself is a one-line consequence of the eigen-expansion: \(\mathbf{S}^t\mathbf{v}\)\({}=\sum_i\lambda_i^tc_i\mathbf{u}_i\), so the normalised iterates \(\mathbf{S}^t\mathbf{v}/\lVert\mathbf{S}^t\mathbf{v}\rVert\) converge to \(\pm\mathbf{u}_1\) at the rate \((\lambda_2/\lambda_1)^t\) whenever \(c_1=\mathbf{u}_1^\T\mathbf{v}\ne0\).

12.1.2 Minimum-error formulation

Now take any complete orthonormal basis \(\{\mathbf{u}_i\}_{i=1}^D\) of \(\R^D\),

\[ \mathbf{u}_i^\T\mathbf{u}_j=\delta_{ij} . \tag{12.7} \]

Every data point has an exact expansion

\[ \mathbf{x}_n=\sum_{i=1}^D\alpha_{ni}\mathbf{u}_i , \tag{12.8} \]

and taking the inner product with \(\mathbf{u}_j\) gives \(\alpha_{nj}\)\({}=\mathbf{x}_n^\T\mathbf{u}_j\), so that

\[ \mathbf{x}_n=\sum_{i=1}^D\big(\mathbf{x}_n^\T\mathbf{u}_i\big)\mathbf{u}_i . \tag{12.9} \]

This is just a rotation of the coordinate system. To reduce the dimension we approximate \(\mathbf{x}_n\) using only \(M\) coefficients that depend on \(n\), the remaining \(D-M\) coefficients being constants \(b_i\) shared by all points:

\[ \widetilde{\mathbf{x}}_n=\sum_{i=1}^Mz_{ni}\mathbf{u}_i+\sum_{i=M+1}^Db_i\mathbf{u}_i . \tag{12.10} \]

(This is the general approximation by a point of an \(M\)-dimensional affine subspace: choose the basis so that \(\mathbf{u}_1,\dots,\mathbf{u}_M\) span the direction of the subspace; the offset is \(\sum_{i\gt M}b_i\mathbf{u}_i\).) The distortion is the mean squared error

\[ J=\frac1N\sum_{n=1}^N\lVert\mathbf{x}_n-\widetilde{\mathbf{x}}_n\rVert^2 . \tag{12.11} \]

Minimising over the coefficients and the offsets gives (Bishop, p. 564; derivation below)

\[ z_{nj}=\mathbf{x}_n^\T\mathbf{u}_j\quad(j=1,\dots,M), \tag{12.12} \] \[ b_j=\bar{\mathbf{x}}^\T\mathbf{u}_j\quad(j=M+1,\dots,D), \tag{12.13} \] \[ \mathbf{x}_n-\widetilde{\mathbf{x}}_n=\sum_{i=M+1}^D\big\{(\mathbf{x}_n-\bar{\mathbf{x}})^\T\mathbf{u}_i\big\}\mathbf{u}_i , \tag{12.14} \] \[ J=\frac1N\sum_{n=1}^N\sum_{i=M+1}^D\big(\mathbf{x}_n^\T\mathbf{u}_i-\bar{\mathbf{x}}^\T\mathbf{u}_i\big)^2=\sum_{i=M+1}^D\mathbf{u}_i^\T\mathbf{S}\mathbf{u}_i . \tag{12.15} \]

It remains to choose the basis. For \(D=2\), \(M=1\) we must choose a unit vector \(\mathbf{u}_2\) minimising \(J\)\({}=\mathbf{u}_2^\T\mathbf{S}\mathbf{u}_2\); with a multiplier \(\lambda_2\) the Lagrangian

\[ \widetilde{J}=\mathbf{u}_2^\T\mathbf{S}\mathbf{u}_2+\lambda_2\big(1-\mathbf{u}_2^\T\mathbf{u}_2\big) \tag{12.16} \]

is stationary when \(\mathbf{S}\mathbf{u}_2\)\({}=\lambda_2\mathbf{u}_2\), and then \(J=\lambda_2\). The minimum takes \(\mathbf{u}_2\) to be the eigenvector of the smaller eigenvalue, so the retained direction \(\mathbf{u}_1\) is the eigenvector of the larger one: the same answer as maximum variance. In general (Exercise 12.2) the solution is given by eigenvectors,

\[ \mathbf{S}\mathbf{u}_i=\lambda_i\mathbf{u}_i\quad(i=1,\dots,D), \tag{12.17} \]

with distortion equal to the sum of the discarded eigenvalues,

\[ J=\sum_{i=M+1}^D\lambda_i . \tag{12.18} \]

Proof (Exercise 12.2, and a shortcut)

The identity. \(\sum_{i=1}^D\mathbf{u}_i^\T\mathbf{S}\mathbf{u}_i\)\({}=\Tr(\mathbf{S}\sum_i\mathbf{u}_i\mathbf{u}_i^\T)\)\({}=\Tr(\mathbf{S}\mathbf{I})\), since \(\sum_i\mathbf{u}_i\mathbf{u}_i^\T\)\({}=\mathbf{I}\) for a complete orthonormal basis. Minimising \(J\) is therefore the same as maximising \(\Tr(\mathbf{V}^\T\mathbf{S}\mathbf{V})\) with \(\mathbf{V}\)\({}=(\mathbf{u}_1,\dots,\mathbf{u}_M)\), and Result 12.2(b) gives the answer at once. The exercise asks for a direct argument with a matrix of multipliers, which also shows what all the stationary points look like.

Lagrangian with a matrix of multipliers. Put the discarded vectors into the \(D\times(D-M)\) matrix \(\widehat{\mathbf{U}}\)\({}=(\mathbf{u}_{M+1},\dots,\mathbf{u}_D)\). Then \(J\)\({}=\Tr(\widehat{\mathbf{U}}^\T\mathbf{S}\widehat{\mathbf{U}})\) and the constraints are \(\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\)\({}=\mathbf{I}\) (orthogonality to the retained vectors does not affect \(J\); the retained vectors can be chosen afterwards to complete the basis). With one multiplier per constraint, arranged in a matrix \(\mathbf{H}\),

\[ \widetilde{J}=\Tr\big\{\widehat{\mathbf{U}}^\T\mathbf{S}\widehat{\mathbf{U}}\big\}+\Tr\big\{\mathbf{H}\big(\mathbf{I}-\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\big)\big\}. \tag{12.93} \]

\(\mathbf{H}\) may be taken symmetric. The matrix \(\mathbf{B}\)\({}=\mathbf{I}-\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\) is symmetric, so \(\Tr(\mathbf{H}\mathbf{B})\)\({}=\Tr((\mathbf{H}\mathbf{B})^\T)\)\({}=\Tr(\mathbf{B}\mathbf{H}^\T)\)\({}=\Tr(\mathbf{H}^\T\mathbf{B})\), and \(\Tr(\mathbf{H}\mathbf{B})\)\({}=\Tr(\tfrac12(\mathbf{H}+\mathbf{H}^\T)\mathbf{B})\): only the symmetric part of \(\mathbf{H}\) enters (the \((i,j)\) and \((j,i)\) constraints are the same equation).

Stationarity. With \(\partial\Tr(\mathbf{A}^\T\mathbf{S}\mathbf{A})/\partial\mathbf{A}\)\({}=2\mathbf{S}\mathbf{A}\) and \(\partial\Tr(\mathbf{H}\mathbf{A}^\T\mathbf{A})/\partial\mathbf{A}\)\({}=\mathbf{A}(\mathbf{H}+\mathbf{H}^\T)\)\({}=2\mathbf{A}\mathbf{H}\) (toolkit), the derivative of (12.93) vanishes when

\[ \mathbf{S}\widehat{\mathbf{U}}=\widehat{\mathbf{U}}\mathbf{H}. \]

All stationary points. Diagonalise the symmetric \(\mathbf{H}\)\({}=\bm{\Psi}\widehat{\bm{\Lambda}}\bm{\Psi}^\T\) with \(\bm{\Psi}\) orthogonal and \(\widehat{\bm{\Lambda}}\) diagonal. The matrix \(\widehat{\mathbf{V}}\)\({}=\widehat{\mathbf{U}}\bm{\Psi}\) still has orthonormal columns and \(\mathbf{S}\widehat{\mathbf{V}}\)\({}=\widehat{\mathbf{U}}\mathbf{H}\bm{\Psi}\)\({}=\widehat{\mathbf{V}}\widehat{\bm{\Lambda}}\): its columns are \(D-M\) orthonormal eigenvectors of \(\mathbf{S}\), with eigenvalues on the diagonal of \(\widehat{\bm{\Lambda}}\). So at a stationary point the columns of \(\widehat{\mathbf{U}}\) are an arbitrary orthonormal basis (\(\widehat{\mathbf{U}}\)\({}=\widehat{\mathbf{V}}\bm{\Psi}^\T\)) of a subspace spanned by \(D-M\) eigenvectors, and

\[ J=\Tr(\widehat{\mathbf{U}}^\T\mathbf{S}\widehat{\mathbf{U}})=\Tr(\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\mathbf{H})=\Tr\mathbf{H}=\Tr\widehat{\bm{\Lambda}}, \]

the sum of the eigenvalues of those \(D-M\) eigenvectors. The particular solution with \(\bm{\Psi}=\mathbf{I}\) (\(\mathbf{H}\) diagonal) has the eigenvectors themselves as columns; every other stationary point with the same subspace gives the same \(J\), which is why the book may "choose the eigenvector solution".

The minimum. The feasible set \(\{\widehat{\mathbf{U}}:\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\)\({}=\mathbf{I}\}\) is closed and bounded, so a minimiser exists; it satisfies the multiplier rule because the constraint map has full rank at every feasible point (for any symmetric \(\mathbf{B}\), moving \(\widehat{\mathbf{U}}\) in the direction \(\tfrac12\widehat{\mathbf{U}}\mathbf{B}\) changes \(\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\) by \(\mathbf{B}\) to first order). Hence the minimum is the smallest possible sum of \(D-M\) eigenvalues, \(\sum_{i\gt M}\lambda_i\), which is (12.18), and the principal subspace is spanned by the eigenvectors of the \(M\) largest eigenvalues. ∎

If \(M=D\) nothing is discarded and PCA is simply a rotation of the axes onto the eigenvectors. A related linear method, canonical correlation analysis (CCA; Hotelling, 1936; Bach and Jordan, 2002), works with two random vectors \(\mathbf{x}\) and \(\mathbf{y}\) and looks for pairs of directions \(\mathbf{a},\mathbf{b}\) whose projections are maximally correlated. Maximising \(\mathbf{a}^\T\mathbf{S}_{xy}\mathbf{b}\) subject to \(\mathbf{a}^\T\mathbf{S}_{xx}\mathbf{a}\)\({}=\mathbf{b}^\T\mathbf{S}_{yy}\mathbf{b}\)\({}=1\) with two multipliers gives \(\mathbf{S}_{xy}\mathbf{b}\)\({}=\rho\,\mathbf{S}_{xx}\mathbf{a}\) and \(\mathbf{S}_{yx}\mathbf{a}\)\({}=\rho\,\mathbf{S}_{yy}\mathbf{b}\) (multiplying by \(\mathbf{a}^\T\) and \(\mathbf{b}^\T\) shows the two multipliers both equal the correlation \(\rho\)), and eliminating \(\mathbf{b}\) yields the generalised eigenproblem \(\mathbf{S}_{xy}\mathbf{S}_{yy}^{-1}\mathbf{S}_{yx}\mathbf{a}\)\({}=\rho^2\mathbf{S}_{xx}\mathbf{a}\) mentioned by the book.

12.1.3 Applications of PCA

Compression of digit images

Take the 1010 images of the digit three in the MNIST test set (the book's "off-line digits" are MNIST, Appendix A); each is a vector of \(D=28\times28=784\) grey levels. Every eigenvector of \(\mathbf{S}\) is itself a vector in \(\R^{784}\) and can be shown as an image (Figure 12.3). Figure 12.4 shows the eigenvalue spectrum and the distortion \(J\) as a function of \(M\), which by (12.18) is the tail sum of the spectrum. With grey levels 0 to 255, as in the book, the first four eigenvalues are \(3.42\), \(2.80\), \(2.38\) and \(1.61\times10^5\), matching the book's \(3.4\), \(2.8\), \(2.4\) and \(1.6\times10^5\) (the 6131 threes of the training set give \(3.64\), \(2.71\), \(2.32\) and \(1.60\times10^5\), so the test set is the better match). In our computation, keeping \(M=1\), 10, 50 and 250 components retains 12.1%, 54.7%, 85.4% and 99.0% of the total variance \(\Tr\mathbf{S}\). Since 254 of the 784 pixels are blank in every one of these images, at least 254 eigenvalues are exactly zero. (Errata: the captions of Figures 12.3 and 12.4 should say "digit three", and the text on p. 566 shows the first four, not five, eigenvectors.)

Substituting (12.12) and (12.13) into (12.10), the PCA approximation of \(\mathbf{x}_n\) is

\[ \begin{aligned}\widetilde{\mathbf{x}}_n&=\sum_{i=1}^M\big(\mathbf{x}_n^\T\mathbf{u}_i\big)\mathbf{u}_i+\sum_{i=M+1}^D\big(\bar{\mathbf{x}}^\T\mathbf{u}_i\big)\mathbf{u}_i &(12.19)\\ &=\bar{\mathbf{x}}+\sum_{i=1}^M\big(\mathbf{x}_n^\T\mathbf{u}_i-\bar{\mathbf{x}}^\T\mathbf{u}_i\big)\mathbf{u}_i , &(12.20)\end{aligned} \]

where the second line replaces \(\sum_{i\gt M}(\bar{\mathbf{x}}^\T\mathbf{u}_i)\mathbf{u}_i\) by \(\bar{\mathbf{x}}-\sum_{i\le M}(\bar{\mathbf{x}}^\T\mathbf{u}_i)\mathbf{u}_i\), using the expansion of the mean in the complete basis,

\[ \bar{\mathbf{x}}=\sum_{i=1}^D\big(\bar{\mathbf{x}}^\T\mathbf{u}_i\big)\mathbf{u}_i . \tag{12.21} \]

This is a compression: once \(\bar{\mathbf{x}}\) and \(\mathbf{u}_1,\dots,\mathbf{u}_M\) are stored (\(D(M+1)\) numbers, shared by the whole data set), each image is described by the \(M\) scores \(\mathbf{x}_n^\T\mathbf{u}_i-\bar{\mathbf{x}}^\T\mathbf{u}_i\) instead of its \(D\) pixels. Figure 12.5 shows reconstructions for several \(M\).

Pre-processing: standardisation and whitening

PCA is also used to put data into a standard form before another algorithm is applied, for example when variables are measured in different units. In the Old Faithful data the waiting time is about an order of magnitude larger than the eruption duration, which is why Section 9.1 first standardised each variable to zero mean and unit variance. The covariance matrix of the standardised data is the correlation matrix

\[ \rho_{ij}=\frac1N\sum_{n=1}^N\frac{(x_{ni}-\bar{x}_i)}{\sigma_i}\frac{(x_{nj}-\bar{x}_j)}{\sigma_j} , \tag{12.22} \]

where \(\sigma_i\) is the standard deviation of variable \(i\) (computed with divisor \(N\); erratum: the first printing says "variance"). It has \(\rho_{ii}=1\) and \(\lvert\rho_{ij}\rvert\le1\) (Cauchy–Schwarz), with \(\rho_{ij}=\pm1\) for perfectly (linearly) correlated variables and 0 for uncorrelated ones. Standardisation rescales the variables but does not remove their correlation. PCA can do more: write the eigenvector equations (12.17) together as

\[ \mathbf{S}\mathbf{U}=\mathbf{U}\mathbf{L} , \tag{12.23} \]

with \(\mathbf{L}\)\({}=\diag(\lambda_1,\dots,\lambda_D)\) and \(\mathbf{U}\) the orthogonal matrix with columns \(\mathbf{u}_i\), and define

\[ \mathbf{y}_n=\mathbf{L}^{-1/2}\mathbf{U}^\T(\mathbf{x}_n-\bar{\mathbf{x}}) . \tag{12.24} \]

Proof

Mean. The map (12.24) is affine, so \(\frac1N\sum_n\mathbf{y}_n\)\({}=\mathbf{L}^{-1/2}\mathbf{U}^\T\big(\frac1N\sum_n\mathbf{x}_n-\bar{\mathbf{x}}\big)\)\({}=\mathbf{0}\); hence the sample covariance of the \(\mathbf{y}_n\) is \(\frac1N\sum_n\mathbf{y}_n\mathbf{y}_n^\T\), the left side of (12.25).

Covariance. The first equality in (12.25) is the definition of \(\mathbf{y}_n\) (note \((\mathbf{L}^{-1/2})^\T\)\({}=\mathbf{L}^{-1/2}\) because it is diagonal); the second takes the constant matrices out of the sum and uses (12.3). For the third, multiply (12.23) on the left by \(\mathbf{U}^\T\) and use \(\mathbf{U}^\T\mathbf{U}=\mathbf{I}\): \(\mathbf{U}^\T\mathbf{S}\mathbf{U}\)\({}=\mathbf{L}\). Diagonal matrices commute, so \(\mathbf{L}^{-1/2}\mathbf{L}\mathbf{L}^{-1/2}\)\({}=\mathbf{I}\). The inverse square root exists because every \(\lambda_i\gt0\).

Distances. \(\mathbf{y}_n-\mathbf{y}_m\)\({}=\mathbf{L}^{-1/2}\mathbf{U}^\T(\mathbf{x}_n-\mathbf{x}_m)\), so \(\lVert\mathbf{y}_n-\mathbf{y}_m\rVert^2\)\({}=(\mathbf{x}_n-\mathbf{x}_m)^\T\mathbf{U}\mathbf{L}^{-1}\mathbf{U}^\T(\mathbf{x}_n-\mathbf{x}_m)\), and \(\mathbf{U}\mathbf{L}^{-1}\mathbf{U}^\T\)\({}=\mathbf{S}^{-1}\). ∎

Two remarks. Whitening is not unique: for any orthogonal \(\mathbf{Q}\), the points \(\mathbf{Q}\mathbf{y}_n\) are white too, since their covariance is \(\mathbf{Q}\mathbf{I}\mathbf{Q}^\T\)\({}=\mathbf{I}\) (the symmetric choice \(\mathbf{Q}=\mathbf{U}\), i.e. \(\mathbf{S}^{-1/2}(\mathbf{x}_n-\bar{\mathbf{x}})\), changes the data least). This rotational freedom reappears in probabilistic PCA (12.39) and is the reason why independent component analysis needs non-Gaussian sources (Result 12.16). And if some \(\lambda_i=0\) the data lie in a lower-dimensional subspace, and one whitens only the coordinates with \(\lambda_i\gt0\).

PCA versus Fisher's linear discriminant

Both PCA and Fisher's linear discriminant (Section 4.1.4) reduce the dimension linearly, but PCA is unsupervised: it sees only the \(\mathbf{x}_n\). Fisher's discriminant also uses the class labels, and the two can disagree badly (Figure 12.7).

Visualisation

For visualisation one takes \(M=2\) and plots each point at the coordinates \((\mathbf{x}_n^\T\mathbf{u}_1,\mathbf{x}_n^\T\mathbf{u}_2)\) (or the centred scores, which only shifts the picture). Figure 12.8 does this for the oil flow data.

12.1.4 PCA for high-dimensional data

Sometimes \(N\lt D\): a few hundred images, each with millions of pixels. Then \(N\) points span an affine subspace of dimension at most \(N-1\), so at least \(D-N+1\) eigenvalues of \(\mathbf{S}\) are zero and components beyond \(N-1\) are meaningless; worse, an \(O(D^3)\) eigendecomposition is out of reach. The remedy works with the centred \(N\times D\) data matrix \(\mathbf{X}\) (rows \((\mathbf{x}_n-\bar{\mathbf{x}})^\T\)), for which \(\mathbf{S}\)\({}=N^{-1}\mathbf{X}^\T\mathbf{X}\):

\[ \frac1N\mathbf{X}^\T\mathbf{X}\mathbf{u}_i=\lambda_i\mathbf{u}_i . \tag{12.26} \]

Multiplying by \(\mathbf{X}\),

\[ \frac1N\mathbf{X}\mathbf{X}^\T(\mathbf{X}\mathbf{u}_i)=\lambda_i(\mathbf{X}\mathbf{u}_i) , \tag{12.27} \]

so \(\mathbf{v}_i\)\({}=\mathbf{X}\mathbf{u}_i\) solves the \(N\times N\) eigenproblem

\[ \frac1N\mathbf{X}\mathbf{X}^\T\mathbf{v}_i=\lambda_i\mathbf{v}_i . \tag{12.28} \]

Conversely, multiplying (12.28) by \(\mathbf{X}^\T\),

\[ \frac1N\mathbf{X}^\T\mathbf{X}(\mathbf{X}^\T\mathbf{v}_i)=\lambda_i(\mathbf{X}^\T\mathbf{v}_i) , \tag{12.29} \]

so \(\mathbf{X}^\T\mathbf{v}_i\) is an eigenvector of \(\mathbf{S}\), and if \(\lVert\mathbf{v}_i\rVert=1\) the unit eigenvector is

\[ \mathbf{u}_i=\frac{1}{(N\lambda_i)^{1/2}}\mathbf{X}^\T\mathbf{v}_i . \tag{12.30} \]

Proof

(a) The rows of \(\mathbf{X}\) sum to zero, \(\mathbf{1}^\T\mathbf{X}\)\({}=\sum_n(\mathbf{x}_n-\bar{\mathbf{x}})^\T\)\({}=\mathbf{0}^\T\), so they are linearly dependent and the row space has dimension at most \(N-1\). Hence \(\operatorname{rank}\mathbf{S}\)\({}=\operatorname{rank}\mathbf{X}\)\({}\le N-1\) and the null space of \(\mathbf{S}\) has dimension at least \(D-N+1\). Also \(\mathbf{X}^\T\mathbf{1}=\mathbf{0}\) gives \(\mathbf{X}\mathbf{X}^\T\mathbf{1}\)\({}=\mathbf{0}\). (When the points are in general position the rank is exactly \(N-1\), which is the book's "same \(N-1\) eigenvalues".)

(b) Let \(\lambda\ne0\) be an eigenvalue of \(N^{-1}\mathbf{X}^\T\mathbf{X}\) with eigenspace \(E\). For \(\mathbf{u}\in E\setminus\{\mathbf{0}\}\), (12.27) shows that \(\mathbf{X}\mathbf{u}\) is an eigenvector of \(N^{-1}\mathbf{X}\mathbf{X}^\T\) with eigenvalue \(\lambda\), and \(\mathbf{X}\mathbf{u}\ne\mathbf{0}\) because \(\mathbf{X}^\T(\mathbf{X}\mathbf{u})\)\({}=N\lambda\mathbf{u}\ne\mathbf{0}\). So \(\mathbf{u}\mapsto\mathbf{X}\mathbf{u}\) is an injective linear map from \(E\) into the \(\lambda\)-eigenspace of \(N^{-1}\mathbf{X}\mathbf{X}^\T\), and the multiplicity of \(\lambda\) can only grow. The symmetric argument with \(\mathbf{X}^\T\) and (12.29) shows it can only shrink in the other direction, so the multiplicities are equal.

(c) By (12.29) each \(\mathbf{u}_i\) is an eigenvector of \(\mathbf{S}\) with eigenvalue \(\lambda_i\). For the inner products, use (12.28) in the form \(\mathbf{X}\mathbf{X}^\T\mathbf{v}_j\)\({}=N\lambda_j\mathbf{v}_j\):

\[ \mathbf{u}_i^\T\mathbf{u}_j=\frac{\mathbf{v}_i^\T\mathbf{X}\mathbf{X}^\T\mathbf{v}_j}{N(\lambda_i\lambda_j)^{1/2}}=\frac{N\lambda_j\,\mathbf{v}_i^\T\mathbf{v}_j}{N(\lambda_i\lambda_j)^{1/2}}=\delta_{ij}. \]

In particular \(\lVert\mathbf{X}^\T\mathbf{v}_i\rVert^2\)\({}=N\lambda_i\), which is where the factor \((N\lambda_i)^{-1/2}\) in (12.30) comes from. In terms of the thin singular value decomposition \(\mathbf{X}\)\({}=\mathbf{V}\bm{\Sigma}\mathbf{U}^\T\) (toolkit), the \(\mathbf{v}_i\) and \(\mathbf{u}_i\) are the left and right singular vectors and \(\lambda_i=\sigma_i^2/N\); (12.30) is the relation \(\mathbf{X}^\T\mathbf{v}_i\)\({}=\sigma_i\mathbf{u}_i\). ∎

This "snapshot" trick (Sirovich, 1987) reappears in kernel PCA (Section 12.3), where the feature space may be infinite-dimensional and only the \(N\times N\) matrix of inner products is available.

12.2 Probabilistic PCA

So far PCA is a projection with no probability model behind it. This section shows that it is the maximum likelihood solution of a simple latent variable model, probabilistic PCA (PPCA), proposed independently by Tipping and Bishop (1997, 1999b) and by Roweis (1998), and closely related to factor analysis (Section 12.2.4). Having a model brings several things that plain PCA lacks (Bishop, pp. 570–571):

• a Gaussian density whose number of free parameters can be controlled while it still captures the dominant correlations (Result 12.10);
• EM algorithms that need only a few leading eigenvectors, never form \(\mathbf{S}\), cope with missing values and extend to mixtures of PCA models (Sections 12.2.2 and 12.4.3);
• a likelihood, and therefore Bayesian methods that choose the dimension \(M\) automatically (Section 12.2.3), and comparison with other density models: PCA gives a small reconstruction error to any point near the principal subspace, even one far away from all the data, whereas the PPCA density is small there;
• class-conditional densities for classifiers, and samples from the fitted distribution.

Writing \(\widetilde{\mathbf{w}}_d^\T\) for the \(d\)th row of \(\mathbf{W}\), (12.32) factorises as \(p(\mathbf{x}\given\mathbf{z})\)\({}=\prod_{d=1}^D\N(x_d\given\widetilde{\mathbf{w}}_d^\T\mathbf{z}+\mu_d,\sigma^2)\): the observed coordinates are independent given \(\mathbf{z}\), which is the structure of the naive Bayes model (Exercise 12.6). The standard Gaussian prior loses nothing: a general prior \(\N(\mathbf{z}\given\mathbf{m},\bm{\Sigma})\) gives the same family of distributions for \(\mathbf{x}\) after redefining \(\mathbf{W}\) and \(\bm{\mu}\) (Exercise 12.4). The columns of \(\mathbf{W}\) will turn out to span the principal subspace. Note that the model maps the latent space into the data space, the reverse of the projection view of Section 12.1; the projection is recovered later from the posterior \(p(\mathbf{z}\given\mathbf{x})\).

To fit the model by maximum likelihood we need the marginal distribution of \(\mathbf{x}\),

\[ p(\mathbf{x})=\int p(\mathbf{x}\given\mathbf{z})\,p(\mathbf{z})\,\d\mathbf{z}. \tag{12.34} \]

Because the model is linear-Gaussian, this is Gaussian (Exercise 12.7):

\[ p(\mathbf{x})=\N(\mathbf{x}\given\bm{\mu},\mathbf{C}), \tag{12.35} \] \[ \mathbf{C}=\mathbf{W}\mathbf{W}^\T+\sigma^2\mathbf{I}. \tag{12.36} \]

Given that \(\mathbf{x}\) is Gaussian, its mean and covariance follow from (12.33) and the independence of \(\mathbf{z}\) and \(\bm{\epsilon}\) (so \(\E[\mathbf{z}\bm{\epsilon}^\T]\)\({}=\E[\mathbf{z}]\E[\bm{\epsilon}]^\T\)\({}=\mathbf{0}\)):

\[ \E[\mathbf{x}]=\E[\mathbf{W}\mathbf{z}+\bm{\mu}+\bm{\epsilon}]=\bm{\mu}, \tag{12.37} \] \[ \cov[\mathbf{x}]=\E\big[(\mathbf{W}\mathbf{z}+\bm{\epsilon})(\mathbf{W}\mathbf{z}+\bm{\epsilon})^\T\big]=\E\big[\mathbf{W}\mathbf{z}\mathbf{z}^\T\mathbf{W}^\T\big]+\E[\bm{\epsilon}\bm{\epsilon}^\T]=\mathbf{W}\mathbf{W}^\T+\sigma^2\mathbf{I}. \tag{12.38} \]

(It is Gaussian because \((\mathbf{z},\bm{\epsilon})\) is jointly Gaussian, being a pair of independent Gaussian vectors, and \(\mathbf{x}\) is an affine function of it; toolkit.) Bishop's picture: an isotropic Gaussian "spray can" moved across the principal subspace, spraying ink with density weighted by the prior, leaves a pancake-shaped density in data space.

The parametrisation is redundant. If \(\widetilde{\mathbf{W}}\)\({}=\mathbf{W}\mathbf{R}\) with \(\mathbf{R}\) an \(M\times M\) orthogonal matrix (\(\mathbf{R}\mathbf{R}^\T=\mathbf{I}\)), then

\[ \widetilde{\mathbf{W}}\widetilde{\mathbf{W}}^\T=\mathbf{W}\mathbf{R}\mathbf{R}^\T\mathbf{W}^\T=\mathbf{W}\mathbf{W}^\T , \tag{12.39} \]

so \(\mathbf{W}\) and \(\mathbf{W}\mathbf{R}\) define the same distribution of \(\mathbf{x}\): rotating the latent coordinates (\(\mathbf{z}\mapsto\mathbf{R}^\T\mathbf{z}\) leaves \(\N(\mathbf{0},\mathbf{I})\) unchanged) and counter-rotating \(\mathbf{W}\) changes nothing that can be observed. Evaluating the density requires \(\mathbf{C}^{-1}\), which the Woodbury identity reduces to an \(M\times M\) inversion:

\[ \mathbf{C}^{-1}=\sigma^{-2}\mathbf{I}-\sigma^{-2}\mathbf{W}\mathbf{M}^{-1}\mathbf{W}^\T , \tag{12.40} \] \[ \mathbf{M}=\mathbf{W}^\T\mathbf{W}+\sigma^2\mathbf{I}. \tag{12.41} \]

(Erratum: the first printing has \(\sigma^{-1}\mathbf{I}\) as the first term of (12.40).) The posterior distribution of the latent variable is

\[ p(\mathbf{z}\given\mathbf{x})=\N\big(\mathbf{z}\given\mathbf{M}^{-1}\mathbf{W}^\T(\mathbf{x}-\bm{\mu}),\,\sigma^2\mathbf{M}^{-1}\big) \tag{12.42} \]

(erratum: the first printing has \(\sigma^{-2}\mathbf{M}\) as the covariance). Its mean depends on \(\mathbf{x}\), its covariance does not.

Proof (all four facts from one completion of the square)

The book obtains (12.35) from the general linear-Gaussian results (2.115) and (12.42) from (2.116) (Exercises 12.7 and 12.8; see the toolkit). Here everything is derived directly. The joint density is \(p(\mathbf{x}\given\mathbf{z})p(\mathbf{z})\)\({}=(2\pi)^{-(D+M)/2}\sigma^{-D}\exp(-\tfrac12E)\) with

\[ E=\frac{1}{\sigma^2}\lVert\mathbf{x}-\bm{\mu}-\mathbf{W}\mathbf{z}\rVert^2+\mathbf{z}^\T\mathbf{z}=\mathbf{z}^\T\Big(\mathbf{I}+\frac{\mathbf{W}^\T\mathbf{W}}{\sigma^2}\Big)\mathbf{z}-\frac{2}{\sigma^2}\mathbf{z}^\T\mathbf{W}^\T\mathbf{y}+\frac{\lVert\mathbf{y}\rVert^2}{\sigma^2},\qquad\mathbf{y}=\mathbf{x}-\bm{\mu}. \]

The matrix of the quadratic form is \(\mathbf{I}+\sigma^{-2}\mathbf{W}^\T\mathbf{W}\)\({}=\sigma^{-2}\mathbf{M}\). Completing the square (toolkit) with \(\mathbf{m}\)\({}=\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{y}\),

\[ E=\frac{1}{\sigma^2}(\mathbf{z}-\mathbf{m})^\T\mathbf{M}(\mathbf{z}-\mathbf{m})+\frac{1}{\sigma^2}\Big(\lVert\mathbf{y}\rVert^2-\mathbf{y}^\T\mathbf{W}\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{y}\Big). \]

(Expand the first term: \(\sigma^{-2}\{\mathbf{z}^\T\mathbf{M}\mathbf{z}-2\mathbf{z}^\T\mathbf{W}^\T\mathbf{y}+\mathbf{y}^\T\mathbf{W}\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{y}\}\), so the last term of the second bracket cancels the extra piece.) (c): as a function of \(\mathbf{z}\), the joint density is proportional to \(\exp\{-\tfrac{1}{2\sigma^2}(\mathbf{z}-\mathbf{m})^\T\mathbf{M}(\mathbf{z}-\mathbf{m})\}\), a Gaussian with mean \(\mathbf{m}\) and covariance \(\sigma^2\mathbf{M}^{-1}\), which is (12.42). (b), first part: the second bracket is \(\mathbf{y}^\T\mathbf{A}\mathbf{y}\) with \(\mathbf{A}\)\({}=\sigma^{-2}(\mathbf{I}-\mathbf{W}\mathbf{M}^{-1}\mathbf{W}^\T)\), and \(\mathbf{A}=\mathbf{C}^{-1}\):

\[ \mathbf{C}\mathbf{A}=\frac{1}{\sigma^2}\big(\mathbf{W}\mathbf{W}^\T+\sigma^2\mathbf{I}-\mathbf{W}(\mathbf{W}^\T\mathbf{W}+\sigma^2\mathbf{I})\mathbf{M}^{-1}\mathbf{W}^\T\big)=\frac{1}{\sigma^2}\big(\mathbf{W}\mathbf{W}^\T+\sigma^2\mathbf{I}-\mathbf{W}\mathbf{W}^\T\big)=\mathbf{I}. \]

(This is the Woodbury identity with \(\mathbf{A}=\sigma^2\mathbf{I}\), \(\mathbf{U}=\mathbf{W}\), \(\mathbf{V}=\mathbf{W}^\T\).) (a) and the determinant: integrating the Gaussian factor over \(\mathbf{z}\) gives \((2\pi)^{M/2}\lvert\sigma^2\mathbf{M}^{-1}\rvert^{1/2}\)\({}=(2\pi)^{M/2}\sigma^M\lvert\mathbf{M}\rvert^{-1/2}\), so

\[ p(\mathbf{x})=(2\pi)^{-D/2}\,\sigma^{-(D-M)}\lvert\mathbf{M}\rvert^{-1/2}\exp\big\{-\tfrac12(\mathbf{x}-\bm{\mu})^\T\mathbf{C}^{-1}(\mathbf{x}-\bm{\mu})\big\}. \]

This is a Gaussian with mean \(\bm{\mu}\) and covariance \(\mathbf{C}\), provided the constant equals \((2\pi)^{-D/2}\lvert\mathbf{C}\rvert^{-1/2}\), that is \(\lvert\mathbf{C}\rvert\)\({}=\sigma^{2(D-M)}\lvert\mathbf{M}\rvert\). That holds by the matrix determinant lemma \(\lvert\mathbf{I}_D+\mathbf{P}\mathbf{Q}\rvert\)\({}=\lvert\mathbf{I}_M+\mathbf{Q}\mathbf{P}\rvert\): \(\lvert\sigma^2\mathbf{I}_D+\mathbf{W}\mathbf{W}^\T\rvert\)\({}=\sigma^{2D}\lvert\mathbf{I}_M+\sigma^{-2}\mathbf{W}^\T\mathbf{W}\rvert\)\({}=\sigma^{2D}\sigma^{-2M}\lvert\mathbf{M}\rvert\). (In fact no lemma is needed: a density integrates to one, so the constant is forced.) (d): \(\mathbf{C}\mathbf{W}\)\({}=\mathbf{W}\mathbf{W}^\T\mathbf{W}+\sigma^2\mathbf{W}\)\({}=\mathbf{W}\mathbf{M}\); multiply by \(\mathbf{C}^{-1}\) on the left and \(\mathbf{M}^{-1}\) on the right. The cost statement follows from (b): \(\ln\lvert\mathbf{C}\rvert\)\({}=(D-M)\ln\sigma^2+\ln\lvert\mathbf{M}\rvert\) and the quadratic form needs only \(\mathbf{b}=\mathbf{W}^\T\mathbf{y}\). ∎

12.2.1 Maximum likelihood PCA

For independent observations \(\mathbf{X}=\{\mathbf{x}_n\}\), (12.35) gives the log likelihood

\[ \ln p(\mathbf{X}\given\bm{\mu},\mathbf{W},\sigma^2)=\sum_{n=1}^N\ln p(\mathbf{x}_n\given\mathbf{W},\bm{\mu},\sigma^2)=-\frac{ND}{2}\ln(2\pi)-\frac N2\ln\lvert\mathbf{C}\rvert-\frac12\sum_{n=1}^N(\mathbf{x}_n-\bm{\mu})^\T\mathbf{C}^{-1}(\mathbf{x}_n-\bm{\mu}). \tag{12.43} \]

Maximising over \(\bm{\mu}\) gives \(\bm{\mu}_{\rm ML}=\bar{\mathbf{x}}\), and substituting back,

\[ \ln p(\mathbf{X}\given\mathbf{W},\bm{\mu},\sigma^2)=-\frac N2\Big\{D\ln(2\pi)+\ln\lvert\mathbf{C}\rvert+\Tr\big(\mathbf{C}^{-1}\mathbf{S}\big)\Big\}. \tag{12.44} \]

Proof

Only the last term of (12.43) involves \(\bm{\mu}\). Its gradient is \(\sum_n\mathbf{C}^{-1}(\mathbf{x}_n-\bm{\mu})\)\({}=\mathbf{C}^{-1}(N\bar{\mathbf{x}}-N\bm{\mu})\) (using the symmetry of \(\mathbf{C}^{-1}\)), which vanishes only at \(\bm{\mu}=\bar{\mathbf{x}}\) since \(\mathbf{C}^{-1}\) is invertible. The Hessian is \(-N\mathbf{C}^{-1}\), negative definite because \(\mathbf{C}\) is positive definite: \(\mathbf{v}^\T\mathbf{C}\mathbf{v}\)\({}=\lVert\mathbf{W}^\T\mathbf{v}\rVert^2+\sigma^2\lVert\mathbf{v}\rVert^2\)\({}\gt0\) for \(\mathbf{v}\ne\mathbf{0}\), and the inverse of a positive definite matrix is positive definite (its eigenvalues are the reciprocals). A strictly concave quadratic has its stationary point as unique global maximiser. For (12.44), use the trace trick \(\mathbf{a}^\T\mathbf{B}\mathbf{a}\)\({}=\Tr(\mathbf{B}\mathbf{a}\mathbf{a}^\T)\) (toolkit):

\[ \sum_n(\mathbf{x}_n-\bar{\mathbf{x}})^\T\mathbf{C}^{-1}(\mathbf{x}_n-\bar{\mathbf{x}})=\Tr\Big(\mathbf{C}^{-1}\sum_n(\mathbf{x}_n-\bar{\mathbf{x}})(\mathbf{x}_n-\bar{\mathbf{x}})^\T\Big)=N\Tr(\mathbf{C}^{-1}\mathbf{S}). \]

Collecting the factor \(-N/2\) gives (12.44). ∎

Maximisation over \(\mathbf{W}\) and \(\sigma^2\) is harder but still has a closed form. The book quotes it from Tipping and Bishop (1999b): every stationary point has the form

\[ \mathbf{W}_{\rm ML}=\mathbf{U}_M(\mathbf{L}_M-\sigma^2\mathbf{I})^{1/2}\mathbf{R}, \tag{12.45} \]

where the columns of \(\mathbf{U}_M\) are any \(M\) eigenvectors of \(\mathbf{S}\), \(\mathbf{L}_M\) is the diagonal matrix of their eigenvalues and \(\mathbf{R}\) is an arbitrary \(M\times M\) orthogonal matrix; the maximum is attained when the eigenvectors are those of the \(M\) largest eigenvalues (all other choices being saddle points), and then

\[ \sigma^2_{\rm ML}=\frac{1}{D-M}\sum_{i=M+1}^D\lambda_i , \tag{12.46} \]

the average variance of the discarded directions. The book gives no proof; the next result supplies it.

Proof (the full Tipping–Bishop argument, with the existence of the maximum added)

Write \(\mathcal{L}(\mathbf{W},\sigma^2)\) for (12.44).

Step 1: the gradient and the stationarity condition. With \(\d\ln\lvert\mathbf{C}\rvert\)\({}=\Tr(\mathbf{C}^{-1}\d\mathbf{C})\), \(\d(\mathbf{C}^{-1})\)\({}=-\mathbf{C}^{-1}(\d\mathbf{C})\mathbf{C}^{-1}\) (toolkit) and \(\d\mathbf{C}\)\({}=(\d\mathbf{W})\mathbf{W}^\T+\mathbf{W}(\d\mathbf{W})^\T\), and using the cyclic property and symmetry of the trace,

\[ \d\ln\lvert\mathbf{C}\rvert=2\Tr\big(\mathbf{W}^\T\mathbf{C}^{-1}\d\mathbf{W}\big),\qquad \d\Tr(\mathbf{C}^{-1}\mathbf{S})=-\Tr\big(\mathbf{C}^{-1}\d\mathbf{C}\,\mathbf{C}^{-1}\mathbf{S}\big)=-2\Tr\big(\mathbf{W}^\T\mathbf{C}^{-1}\mathbf{S}\mathbf{C}^{-1}\d\mathbf{W}\big). \]

Since \(\d f\)\({}=\Tr(\mathbf{G}^\T\d\mathbf{W})\) identifies the gradient \(\mathbf{G}\),

\[ \frac{\partial\mathcal{L}}{\partial\mathbf{W}}=N\big(\mathbf{C}^{-1}\mathbf{S}\mathbf{C}^{-1}\mathbf{W}-\mathbf{C}^{-1}\mathbf{W}\big),\qquad \frac{\partial\mathcal{L}}{\partial\sigma^2}=-\frac N2\Tr\big(\mathbf{C}^{-1}-\mathbf{C}^{-1}\mathbf{S}\mathbf{C}^{-1}\big). \]

Multiplying the first by \(\mathbf{C}/N\), the stationarity condition in \(\mathbf{W}\) is

\[ \mathbf{S}\mathbf{C}^{-1}\mathbf{W}=\mathbf{W}. \]

Tipping and Bishop distinguish three kinds of solution: (i) \(\mathbf{W}=\mathbf{0}\), which satisfies it trivially; (ii) \(\mathbf{C}=\mathbf{S}\), the case of an exact covariance model, possible only when \(\mathbf{S}-\sigma^2\mathbf{I}\) is positive semi-definite of rank at most \(M\), i.e. when the \(D-M\) smallest eigenvalues of \(\mathbf{S}\) all equal \(\sigma^2\); (iii) \(\mathbf{W}\ne\mathbf{0}\) and \(\mathbf{C}\ne\mathbf{S}\), the generic case. The singular value decomposition handles all three at once.

Step 2: the singular value decomposition of \(\mathbf{W}\). Write \(\mathbf{W}\)\({}=\mathbf{U}\mathbf{L}\mathbf{V}^\T\) with \(\mathbf{U}\) (\(D\times M\)) having orthonormal columns \(\mathbf{u}_j'\), \(\mathbf{L}=\diag(l_1,\dots,l_M)\) with \(l_j\ge0\), and \(\mathbf{V}\) orthogonal (toolkit; if \(\operatorname{rank}\mathbf{W}\lt M\), some \(l_j=0\) and the corresponding \(\mathbf{u}_j'\) complete the orthonormal set arbitrarily). By Result 12.6(d), \(\mathbf{C}^{-1}\mathbf{W}\)\({}=\mathbf{W}\mathbf{M}^{-1}\), and \(\mathbf{W}^\T\mathbf{W}\)\({}=\mathbf{V}\mathbf{L}^2\mathbf{V}^\T\) gives \(\mathbf{M}^{-1}\)\({}=\mathbf{V}(\mathbf{L}^2+\sigma^2\mathbf{I})^{-1}\mathbf{V}^\T\), so

\[ \mathbf{C}^{-1}\mathbf{W}=\mathbf{U}\mathbf{L}(\mathbf{L}^2+\sigma^2\mathbf{I})^{-1}\mathbf{V}^\T . \]

The stationarity condition becomes \(\mathbf{S}\mathbf{U}\mathbf{L}(\mathbf{L}^2+\sigma^2\mathbf{I})^{-1}\mathbf{V}^\T\)\({}=\mathbf{U}\mathbf{L}\mathbf{V}^\T\); multiplying on the right by \(\mathbf{V}(\mathbf{L}^2+\sigma^2\mathbf{I})\) (diagonal matrices commute),

\[ \mathbf{S}\mathbf{U}\mathbf{L}=\mathbf{U}\mathbf{L}(\mathbf{L}^2+\sigma^2\mathbf{I}),\qquad\text{i.e.}\qquad l_j\,\mathbf{S}\mathbf{u}_j'=l_j(\sigma^2+l_j^2)\,\mathbf{u}_j'\quad(j=1,\dots,M). \]

So either \(l_j=0\) (the \(j\)th column of \(\mathbf{W}\mathbf{V}\) vanishes and \(\mathbf{u}_j'\) is unconstrained), or \(\mathbf{u}_j'\) is an eigenvector of \(\mathbf{S}\) whose eigenvalue \(\lambda\) satisfies \(\lambda=\sigma^2+l_j^2\), i.e. \(l_j=(\lambda-\sigma^2)^{1/2}\), which requires \(\lambda\gt\sigma^2\). Setting \(k_j=\sigma^2+l_j^2\) and \(\mathbf{R}=\mathbf{V}^\T\) gives \(\mathbf{W}\)\({}=\mathbf{U}(\mathbf{K}-\sigma^2\mathbf{I})^{1/2}\mathbf{R}\), as in (a); conversely every such \(\mathbf{W}\) satisfies the condition (run the steps backwards). The columns of \(\mathbf{U}\) are orthonormal, so no eigenvector is used twice; within a repeated eigenvalue any orthonormal choice is allowed. Case (i) is all \(l_j=0\); case (ii) is the special case where the unused eigenvalues equal \(\sigma^2\).

Step 3: the noise variance. Fix a stationary \(\mathbf{W}\) whose \(q\) non-zero columns use the eigenvectors \(\mathbf{u}_j\), \(j\in\mathcal{I}\) (\(\lvert\mathcal{I}\rvert=q\le M\)). Then \(\mathbf{C}\)\({}=\sigma^2\mathbf{I}+\sum_{j\in\mathcal{I}}(\lambda_j-\sigma^2)\mathbf{u}_j\mathbf{u}_j^\T\), so \(\mathbf{C}\) has the same eigenvectors as \(\mathbf{S}\), with eigenvalue \(\lambda_j\) on \(\mathbf{u}_j\) for \(j\in\mathcal{I}\) and \(\sigma^2\) on the other \(D-q\) eigenvectors. Therefore

\[ \ln\lvert\mathbf{C}\rvert=\sum_{j\in\mathcal{I}}\ln\lambda_j+(D-q)\ln\sigma^2,\qquad \Tr(\mathbf{C}^{-1}\mathbf{S})=\sum_{j\in\mathcal{I}}\frac{\lambda_j}{\lambda_j}+\sum_{i\notin\mathcal{I}}\frac{\lambda_i}{\sigma^2}=q+\frac{1}{\sigma^2}\sum_{i\notin\mathcal{I}}\lambda_i , \] \[ \mathcal{L}=-\frac N2\Big\{D\ln(2\pi)+\sum_{j\in\mathcal{I}}\ln\lambda_j+(D-q)\ln\sigma^2+q+\frac{1}{\sigma^2}\sum_{i\notin\mathcal{I}}\lambda_i\Big\}. \]

As a function of \(\sigma^2\) this has the form \(-\tfrac N2\{a\ln\sigma^2+b/\sigma^2\}+\const\) with \(a=D-q\gt0\), \(b\)\({}=\sum_{i\notin\mathcal{I}}\lambda_i\)\({}\gt0\), whose derivative \(-\tfrac N2\{a/\sigma^2-b/\sigma^4\}\) vanishes only at \(\sigma^2=b/a\), a maximum. So

\[ \sigma^2=\frac{1}{D-q}\sum_{i\notin\mathcal{I}}\lambda_i , \]

the average of the eigenvalues not used, which for \(q=M\) and \(\mathcal{I}=\{1,\dots,M\}\) is (12.46). (The same equation comes from \(\partial\mathcal{L}/\partial\sigma^2\)\({}=0\): \(\Tr\mathbf{C}^{-1}\)\({}=\sum_{\mathcal{I}}\lambda_j^{-1}+(D-q)\sigma^{-2}\) and \(\Tr(\mathbf{C}^{-1}\mathbf{S}\mathbf{C}^{-1})\)\({}=\sum_{\mathcal{I}}\lambda_j^{-1}+\sigma^{-4}\sum_{i\notin\mathcal{I}}\lambda_i\).) The retained eigenvalues must exceed this average: a choice of \(\mathcal{I}\) that violates this gives no stationary point at all. Substituting \(\sigma^2\) back, \(q+\sigma^{-2}\sum_{i\notin\mathcal{I}}\lambda_i\)\({}=q+(D-q)\)\({}=D\), and

\[ \mathcal{L}(\mathcal{I})=-\frac N2\Big\{D\big(1+\ln(2\pi)\big)+\sum_{j\in\mathcal{I}}\ln\lambda_j+(D-q)\ln\Big(\frac{1}{D-q}\sum_{i\notin\mathcal{I}}\lambda_i\Big)\Big\}, \]

which with \(\mathcal{I}=\{1,\dots,M\}\) is the maximal value stated in (b).

Step 4: zero columns are not maxima. Suppose \(q\lt M\), so \(\mathbf{W}\) (after the rotation \(\mathbf{R}\)) has a zero column, and let \(\mathbf{u}_i\), \(i\notin\mathcal{I}\), be an unused eigenvector. Put \(\varepsilon\mathbf{u}_i\) into the zero column. Then \(\mathbf{C}_\varepsilon\)\({}=\mathbf{C}+\varepsilon^2\mathbf{u}_i\mathbf{u}_i^\T\), and since \(\mathbf{u}_i\) is an eigenvector of \(\mathbf{C}\) with eigenvalue \(\sigma^2\), \(\lvert\mathbf{C}_\varepsilon\rvert\)\({}=\lvert\mathbf{C}\rvert(1+t)\) and \(\Tr(\mathbf{C}_\varepsilon^{-1}\mathbf{S})\)\({}=\Tr(\mathbf{C}^{-1}\mathbf{S})-\lambda_i\big(\tfrac{1}{\sigma^2}-\tfrac{1}{\sigma^2+\varepsilon^2}\big)\), with \(t=\varepsilon^2/\sigma^2\). Hence

\[ \mathcal{L}(\varepsilon)-\mathcal{L}(0)=\frac N2\Big\{\frac{\lambda_i}{\sigma^2}\,\frac{t}{1+t}-\ln(1+t)\Big\}=\frac N2\Big(\frac{\lambda_i}{\sigma^2}-1\Big)t+O(t^2), \]

which is positive for small \(\varepsilon\) whenever \(\lambda_i\gt\sigma^2\). Because \(\sigma^2\) is the average of the unused eigenvalues, one of them exceeds it unless they are all equal. So such a point is not a local maximum; in particular \(\mathbf{W}=\mathbf{0}\) (with \(\sigma^2=\Tr\mathbf{S}/D\)) is not a maximum unless \(\mathbf{S}\propto\mathbf{I}\).

Step 5: non-principal eigenvectors give saddle points. Suppose that at a stationary point (any \(q\)) a used eigenvector \(\mathbf{u}_j\) (\(j\in\mathcal{I}\)) has a smaller eigenvalue than an unused one \(\mathbf{u}_i\) (\(i\notin\mathcal{I}\)): \(\lambda_i\gt\lambda_j\). Rotate in the plane of these two vectors: let \(\mathbf{Q}_\theta\) map \(\mathbf{u}_j\mapsto\cos\theta\,\mathbf{u}_j+\sin\theta\,\mathbf{u}_i\), \(\mathbf{u}_i\mapsto\cos\theta\,\mathbf{u}_i-\sin\theta\,\mathbf{u}_j\) and fix every other eigenvector, and move to \(\mathbf{W}(\theta)\)\({}=\mathbf{Q}_\theta\mathbf{W}\). Then \(\mathbf{C}(\theta)\)\({}=\mathbf{Q}_\theta\mathbf{C}\mathbf{Q}_\theta^\T\), so \(\lvert\mathbf{C}(\theta)\rvert\)\({}=\lvert\mathbf{C}\rvert\), and with \(c_k\) the eigenvalue of \(\mathbf{C}\) on \(\mathbf{u}_k\),

\[ \Tr\big(\mathbf{C}(\theta)^{-1}\mathbf{S}\big)=\Tr\big(\mathbf{C}^{-1}\mathbf{Q}_\theta^\T\mathbf{S}\mathbf{Q}_\theta\big)=\sum_k\frac{(\mathbf{Q}_\theta\mathbf{u}_k)^\T\mathbf{S}(\mathbf{Q}_\theta\mathbf{u}_k)}{c_k}. \]

Only \(k=i,j\) change: \((\mathbf{Q}_\theta\mathbf{u}_j)^\T\mathbf{S}(\mathbf{Q}_\theta\mathbf{u}_j)\)\({}=\lambda_j\cos^2\theta+\lambda_i\sin^2\theta\) and \((\mathbf{Q}_\theta\mathbf{u}_i)^\T\mathbf{S}(\mathbf{Q}_\theta\mathbf{u}_i)\)\({}=\lambda_i\cos^2\theta+\lambda_j\sin^2\theta\), with \(c_j=\lambda_j\) and \(c_i=\sigma^2\). Therefore

\[ \mathcal{L}(\theta)-\mathcal{L}(0)=\frac N2\sin^2\theta\,(\lambda_i-\lambda_j)\Big(\frac{1}{\sigma^2}-\frac{1}{\lambda_j}\Big)\gt0\qquad(\theta\ne0), \]

because \(\lambda_i\gt\lambda_j\gt\sigma^2\). In the opposite direction, changing only the length of column \(j\) changes only the eigenvalue \(c=\sigma^2+l_j^2\) of \(\mathbf{C}\) along \(\mathbf{u}_j\), and \(-\tfrac2N\mathcal{L}\) changes by \(\ln c+\lambda_j/c\), which is minimal at \(c=\lambda_j\): the likelihood decreases. A point that is a maximum along one direction and a minimum along another is a saddle point. (The rotation is also "pure": the likelihood is flat along the \(M(M-1)/2\) directions \(\mathbf{W}\mapsto\mathbf{W}\mathbf{R}\).)

Step 6: the maximum exists. For any \(\mathbf{W}\), \(\mathbf{C}\) has eigenvalue \(\sigma^2\) on the orthogonal complement of the column space of \(\mathbf{W}\); let \(\mathbf{P}\) be the orthogonal projection onto that complement (rank at least \(D-M\)). Then \(\mathbf{C}^{-1}-\sigma^{-2}\mathbf{P}\) is positive semi-definite, so \(\Tr(\mathbf{C}^{-1}\mathbf{S})\)\({}\ge\sigma^{-2}\Tr(\mathbf{P}\mathbf{S})\)\({}\ge\sigma^{-2}\sum_{i\gt M}\lambda_i\), the last step by Result 12.2(b) applied to the complement (the trace of \(\mathbf{S}\) over a subspace of dimension \(r\) is at least the sum of its \(r\) smallest eigenvalues). Also \(\ln\lvert\mathbf{C}\rvert\)\({}\ge D\ln\sigma^2\) and \(\ln\lvert\mathbf{C}\rvert\)\({}\ge\ln(\sigma^2+\lVert\mathbf{W}\rVert_2^2)+(D-1)\ln\sigma^2\). Hence

\[ -\frac2N\mathcal{L}\ \ge\ D\ln(2\pi)+D\ln\sigma^2+\frac{(D-M)\,\sigma^2_{\rm ML}}{\sigma^2}, \]

with \(\sigma^2_{\rm ML}\gt0\) by the assumption on the rank. The right side tends to \(+\infty\) as \(\sigma^2\to0\) or \(\sigma^2\to\infty\), and for \(\sigma^2\) in a closed interval inside \((0,\infty)\) the second lower bound on \(\ln\lvert\mathbf{C}\rvert\) tends to \(+\infty\) as \(\lVert\mathbf{W}\rVert\to\infty\). So the sets \(\{\mathcal{L}\ge c\}\) are closed and bounded, \(\mathcal{L}\) attains its maximum, and the maximiser is a stationary point.

Conclusion. By Step 5 the maximiser has no unused eigenvalue larger than a used one, so it uses the eigenvectors of the \(q\) largest eigenvalues. If \(q\lt M\), Step 4 would force \(\lambda_{q+1}=\dots=\lambda_D\), hence \(\lambda_M=\lambda_D\), contrary to the assumption; so \(q=M\), the eigenvectors are \(\mathbf{u}_1,\dots,\mathbf{u}_M\) (and \(\lambda_M\gt\sigma^2_{\rm ML}\), as required for real \(l_j\), again because \(\lambda_M\gt\lambda_D\)), and by Step 3 the noise variance is (12.46). All such points have the same likelihood: they differ only by \(\mathbf{R}\) and, if \(\lambda_M=\lambda_{M+1}\), by the choice of eigenvectors in the tied eigenspace. ∎

A complementary view compares the stationary values directly. If the top \(q\) eigenvectors are used, Step 3 gives \(-\tfrac2N\mathcal{L}\)\({}=D(1+\ln2\pi)+\psi_q\) with \(\psi_q\)\({}=\sum_{i\le q}\ln\lambda_i+(D-q)\ln\bar\lambda_q\), where \(\bar\lambda_q\) is the average of \(\lambda_{q+1},\dots,\lambda_D\). Since \(\bar\lambda_q\)\({}=\frac{1}{D-q}\lambda_{q+1}+\frac{D-q-1}{D-q}\bar\lambda_{q+1}\), concavity of the logarithm (Jensen) gives \((D-q)\ln\bar\lambda_q\)\({}\ge\ln\lambda_{q+1}+(D-q-1)\ln\bar\lambda_{q+1}\), that is \(\psi_q\ge\psi_{q+1}\): every additional principal direction increases the maximised likelihood (strictly, unless \(\lambda_{q+1},\dots,\lambda_D\) are all equal).

Interpretation. Take \(\mathbf{R}=\mathbf{I}\): the columns of \(\mathbf{W}_{\rm ML}\) are the principal eigenvectors scaled by \((\lambda_i-\sigma^2)^{1/2}\). (Erratum: the first printing omits "square root of" in this sentence on p. 575.) Variances of independent contributions add, so along \(\mathbf{u}_i\) the data variance \(\lambda_i\) is split into \(\lambda_i-\sigma^2\) from the latent variable (pushed through the \(i\)th column of \(\mathbf{W}\)) and \(\sigma^2\) from the noise. For a unit vector \(\mathbf{v}\), the model variance is \(\mathbf{v}^\T\mathbf{C}\mathbf{v}\): if \(\mathbf{v}\) is orthogonal to the principal subspace (\(\mathbf{v}^\T\mathbf{U}_M\)\({}=\mathbf{0}\)) it is \(\sigma^2\), the average discarded eigenvalue; if \(\mathbf{v}=\mathbf{u}_i\) with \(i\le M\) it is \((\lambda_i-\sigma^2)+\sigma^2\)\({}=\lambda_i\). So the fitted model reproduces the data variance exactly along the principal axes and replaces the variance in every other direction by a single average; in particular \(\Tr\mathbf{C}\)\({}=\sum_{i\le M}\lambda_i+(D-M)\sigma^2_{\rm ML}\)\({}=\Tr\mathbf{S}\) and \(\Tr(\mathbf{C}^{-1}\mathbf{S})=D\).

If the maximum likelihood solution is found by numerical optimisation (for example conjugate gradients) or by EM, \(\mathbf{R}\) is essentially arbitrary and the columns of \(\mathbf{W}\) need not be orthogonal. The principal axes are easy to recover: \(\mathbf{W}^\T\mathbf{W}\)\({}=\mathbf{R}^\T(\mathbf{L}_M-\sigma^2\mathbf{I})\mathbf{R}\) is an eigendecomposition, so if \(\mathbf{W}^\T\mathbf{W}\)\({}=\mathbf{Q}\bm{\Delta}\mathbf{Q}^\T\) then the columns of \(\mathbf{W}\mathbf{Q}\bm{\Delta}^{-1/2}\) are the principal eigenvectors (up to sign and order, when the eigenvalues are distinct) and \(\bm{\Delta}+\sigma^2\mathbf{I}\) holds their eigenvalues. Ahn and Oh (2003) modify EM to produce ordered orthonormal directions directly. The rotational freedom is a continuous non-identifiability, analogous to (but unlike) the discrete relabelling symmetry of mixture models.

With \(M=D\) there is no reduction: \(\mathbf{U}_M=\mathbf{U}\), \(\mathbf{L}_M=\mathbf{L}\), and \(\mathbf{U}\mathbf{U}^\T\)\({}=\mathbf{R}\mathbf{R}^\T\)\({}=\mathbf{I}\) give

\[ \mathbf{C}=\mathbf{U}(\mathbf{L}-\sigma^2\mathbf{I})^{1/2}\mathbf{R}\mathbf{R}^\T(\mathbf{L}-\sigma^2\mathbf{I})^{1/2}\mathbf{U}^\T+\sigma^2\mathbf{I}=\mathbf{U}\mathbf{L}\mathbf{U}^\T=\mathbf{S}, \tag{12.47} \]

the unrestricted Gaussian maximum likelihood estimate (any \(\sigma^2\le\lambda_D\) then gives the same \(\mathbf{C}\), so \(\sigma^2\) is no longer determined).

From the data space back to the latent space

For visualisation and compression we reverse the mapping with Bayes' theorem and summarise a point \(\mathbf{x}\) by its posterior. From (12.42), the posterior mean is

\[ \E[\mathbf{z}\given\mathbf{x}]=\mathbf{M}^{-1}\mathbf{W}_{\rm ML}^\T(\mathbf{x}-\bar{\mathbf{x}}) , \tag{12.48} \]

which corresponds to the point

\[ \mathbf{W}\E[\mathbf{z}\given\mathbf{x}]+\bm{\mu} \tag{12.49} \]

in data space. As the book notes, this has the form of regularised linear regression (Section 3.3.1): \(\E[\mathbf{z}\given\mathbf{x}]\) minimises \(\lVert\mathbf{x}-\bm{\mu}-\mathbf{W}\mathbf{z}\rVert^2+\sigma^2\lVert\mathbf{z}\rVert^2\), whose normal equations are \((\mathbf{W}^\T\mathbf{W}+\sigma^2\mathbf{I})\mathbf{z}\)\({}=\mathbf{W}^\T(\mathbf{x}-\bm{\mu})\). The posterior covariance \(\sigma^2\mathbf{M}^{-1}\) does not depend on \(\mathbf{x}\). As \(\sigma^2\to0\) the posterior mean becomes

\[ (\mathbf{W}_{\rm ML}^\T\mathbf{W}_{\rm ML})^{-1}\mathbf{W}_{\rm ML}^\T(\mathbf{x}-\bar{\mathbf{x}}) , \tag{12.50} \]

an orthogonal projection onto the principal subspace, so standard PCA is recovered (Exercise 12.11); the posterior covariance tends to zero and the density becomes singular. For \(\sigma^2\gt0\) the projection is pulled towards the origin (Exercise 12.12).

Proof

(a) \(\mathbf{M}\)\({}=\mathbf{Q}(\bm{\Delta}+\sigma^2\mathbf{I})\mathbf{Q}^\T\), so with \(\mathbf{c}\)\({}=\mathbf{Q}^\T\mathbf{W}^\T(\mathbf{x}-\bm{\mu})\), \(\mathbf{Q}^\T\E[\mathbf{z}\given\mathbf{x}]\)\({}=(\bm{\Delta}+\sigma^2\mathbf{I})^{-1}\mathbf{c}\), while the least-squares coordinates \(\mathbf{z}_\perp\)\({}=(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T(\mathbf{x}-\bm{\mu})\) satisfy \(\mathbf{Q}^\T\mathbf{z}_\perp\)\({}=\bm{\Delta}^{-1}\mathbf{c}\). Componentwise the ratio is \(\delta_k/(\delta_k+\sigma^2)\). As \(\sigma^2\to0\), \(\mathbf{M}\)\({}\to\mathbf{W}^\T\mathbf{W}\), which is invertible, and (12.48) tends to (12.50). The point \(\mathbf{W}\mathbf{z}_\perp+\bm{\mu}\) is the orthogonal projection of \(\mathbf{x}\) onto \(\bm{\mu}+\operatorname{col}\mathbf{W}\), because \(\mathbf{W}(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T\) is the orthogonal projection matrix onto the column space.

(b) With \(\mathbf{R}=\mathbf{I}\), \(\mathbf{W}_{\rm ML}^\T\mathbf{W}_{\rm ML}\)\({}=\mathbf{L}_M-\sigma^2\mathbf{I}\) and \(\mathbf{M}=\mathbf{L}_M\), so \(\E[\mathbf{z}\given\mathbf{x}]\)\({}=\mathbf{L}_M^{-1}(\mathbf{L}_M-\sigma^2\mathbf{I})^{1/2}\mathbf{U}_M^\T(\mathbf{x}-\bar{\mathbf{x}})\) and \(\mathbf{W}_{\rm ML}\E[\mathbf{z}\given\mathbf{x}]\)\({}=\mathbf{U}_M(\mathbf{L}_M-\sigma^2\mathbf{I})\mathbf{L}_M^{-1}\mathbf{U}_M^\T(\mathbf{x}-\bar{\mathbf{x}})\), which is the stated sum. The orthogonal projection onto the principal subspace is \(\mathbf{U}_M\mathbf{U}_M^\T(\mathbf{x}-\bar{\mathbf{x}})\). With \(\sigma^2=0\) the posterior mean is \(\mathbf{L}_M^{-1/2}\mathbf{U}_M^\T(\mathbf{x}-\bar{\mathbf{x}})\), the whitened PCA scores of (12.24).

(c) By (12.48), \(\mathbf{M}\E[\mathbf{z}\given\mathbf{x}]\)\({}=\mathbf{W}_{\rm ML}^\T(\mathbf{x}-\bar{\mathbf{x}})\); substituting into the orthogonal projection \(\bar{\mathbf{x}}+\mathbf{W}(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T(\mathbf{x}-\bar{\mathbf{x}})\) gives the formula. This is the optimal reconstruction in the least-squares sense of Section 12.1.2, since the closest point of an affine subspace is the orthogonal projection. ∎

How many parameters?

PPCA defines a Gaussian whose number of free parameters is controlled by \(M\). A full covariance has \(D(D+1)/2\) parameters (plus \(D\) for the mean), quadratic in \(D\); a diagonal one has \(D\) but cannot express correlations. PPCA keeps the \(M\) dominant correlations with a count that grows only linearly in \(D\). The matrix \(\mathbf{W}\) and \(\sigma^2\) have \(DM+1\) entries, but the rotation (12.39) makes some of them redundant. An \(M\times M\) orthogonal matrix has \(M-1\) free parameters in its first column (a unit vector), \(M-2\) in the second (unit and orthogonal to the first), and so on, \(M(M-1)/2\) in total. The number of degrees of freedom in \(\mathbf{C}\) is therefore

\[ DM+1-M(M-1)/2 . \tag{12.51} \]

Proof

If \(\mathbf{W}_1\mathbf{W}_1^\T\)\({}=\mathbf{W}_2\mathbf{W}_2^\T\) then the column spaces agree (the column space of \(\mathbf{W}\mathbf{W}^\T\) is that of \(\mathbf{W}\)). As \(\mathbf{P}\)\({}=\mathbf{W}_1(\mathbf{W}_1^\T\mathbf{W}_1)^{-1}\mathbf{W}_1^\T\) projects onto it, \(\mathbf{W}_2\)\({}=\mathbf{P}\mathbf{W}_2\)\({}=\mathbf{W}_1\mathbf{R}\) with \(\mathbf{R}\)\({}=(\mathbf{W}_1^\T\mathbf{W}_1)^{-1}\mathbf{W}_1^\T\mathbf{W}_2\). Then \(\mathbf{W}_1\mathbf{R}\mathbf{R}^\T\mathbf{W}_1^\T\)\({}=\mathbf{W}_1\mathbf{W}_1^\T\); multiplying on the left by \((\mathbf{W}_1^\T\mathbf{W}_1)^{-1}\mathbf{W}_1^\T\) and on the right by \(\mathbf{W}_1(\mathbf{W}_1^\T\mathbf{W}_1)^{-1}\) gives \(\mathbf{R}\mathbf{R}^\T=\mathbf{I}\). So the parameters giving one covariance form a copy of the orthogonal group, of dimension \(M(M-1)/2\) (the column count above), and \(\sigma^2\) is determined by \(\mathbf{C}\) when \(M\lt D\) (it is the smallest eigenvalue of \(\mathbf{C}\), of multiplicity at least \(D-M\)). This gives (12.51). For \(M=D-1\): \(D(D-1)+1-\tfrac{(D-1)(D-2)}{2}\)\({}=\tfrac{2D^2-2D+2-D^2+3D-2}{2}\)\({}=\tfrac{D(D+1)}{2}\). For \(M=0\): \(0+1-0=1\). (For \(M=D\) the formula over-counts by one, because then \(\sigma^2\) can be traded against \(\mathbf{W}\mathbf{W}^\T\).) The count was also checked numerically as the rank of the Jacobian of \((\mathbf{W},\sigma^2)\mapsto\mathbf{C}\). ∎

With \(M=D-1\), the columns of \(\mathbf{W}\) control the variance along \(D-1\) independent directions and \(\sigma^2\) the variance along the remaining one, so PPCA interpolates between the isotropic Gaussian (\(M=0\)) and the full-covariance Gaussian (\(M=D-1\)).

12.2.2 EM algorithm for PCA

PPCA has a latent \(\mathbf{z}_n\) for every \(\mathbf{x}_n\), so the general EM algorithm of Section 9.4 applies. With a closed-form maximum already available (Result 12.8) this may look pointless, but EM has real advantages: in high dimensions each iteration is cheaper than an eigendecomposition and never forms \(\mathbf{S}\); the same algorithm extends to factor analysis, which has no closed form (Section 12.2.4); and it handles missing data. Since the data points are independent, the complete-data log likelihood is

\[ \ln p\big(\mathbf{X},\mathbf{Z}\given\bm{\mu},\mathbf{W},\sigma^2\big)=\sum_{n=1}^N\big\{\ln p(\mathbf{x}_n\given\mathbf{z}_n)+\ln p(\mathbf{z}_n)\big\}, \tag{12.52} \]

where the \(n\)th row of \(\mathbf{Z}\) is \(\mathbf{z}_n^\T\). We set \(\bm{\mu}=\bar{\mathbf{x}}\), its exact maximum likelihood value (Result 12.7). Taking the expectation with respect to the posterior of the latent variables gives

\[ \begin{aligned}\E\big[\ln p(\mathbf{X},\mathbf{Z}\given\bm{\mu},\mathbf{W},\sigma^2)\big]=-\sum_{n=1}^N\Big\{&\frac D2\ln(2\pi\sigma^2)+\frac12\Tr\big(\E[\mathbf{z}_n\mathbf{z}_n^\T]\big)+\frac{1}{2\sigma^2}\lVert\mathbf{x}_n-\bm{\mu}\rVert^2\\ &-\frac{1}{\sigma^2}\E[\mathbf{z}_n]^\T\mathbf{W}^\T(\mathbf{x}_n-\bm{\mu})+\frac{1}{2\sigma^2}\Tr\big(\E[\mathbf{z}_n\mathbf{z}_n^\T]\mathbf{W}^\T\mathbf{W}\big)+\frac M2\ln(2\pi)\Big\}.\end{aligned} \tag{12.53} \]

(Erratum: the first printing omits the term \(\frac M2\ln(2\pi)\); it does not affect the maximisation.) The posterior enters only through two expected sufficient statistics. In the E step they are computed from (12.42) with the current parameters:

\[ \E[\mathbf{z}_n]=\mathbf{M}^{-1}\mathbf{W}^\T(\mathbf{x}_n-\bar{\mathbf{x}}), \tag{12.54} \] \[ \E[\mathbf{z}_n\mathbf{z}_n^\T]=\sigma^2\mathbf{M}^{-1}+\E[\mathbf{z}_n]\E[\mathbf{z}_n]^\T , \tag{12.55} \]

using \(\E[\mathbf{z}\mathbf{z}^\T]\)\({}=\cov[\mathbf{z}]+\E[\mathbf{z}]\E[\mathbf{z}]^\T\). The M step maximises (12.53) with these statistics fixed (Exercise 12.15):

\[ \mathbf{W}_{\rm new}=\Big[\sum_{n=1}^N(\mathbf{x}_n-\bar{\mathbf{x}})\E[\mathbf{z}_n]^\T\Big]\Big[\sum_{n=1}^N\E[\mathbf{z}_n\mathbf{z}_n^\T]\Big]^{-1}, \tag{12.56} \] \[ \sigma^2_{\rm new}=\frac{1}{ND}\sum_{n=1}^N\Big\{\lVert\mathbf{x}_n-\bar{\mathbf{x}}\rVert^2-2\E[\mathbf{z}_n]^\T\mathbf{W}_{\rm new}^\T(\mathbf{x}_n-\bar{\mathbf{x}})+\Tr\big(\E[\mathbf{z}_n\mathbf{z}_n^\T]\mathbf{W}_{\rm new}^\T\mathbf{W}_{\rm new}\big)\Big\}. \tag{12.57} \]

Proof

Write \(\mathbf{A}\)\({}=\sum_n(\mathbf{x}_n-\bar{\mathbf{x}})\E[\mathbf{z}_n]^\T\) and \(\mathbf{B}\)\({}=\sum_n\E[\mathbf{z}_n\mathbf{z}_n^\T]\). By (12.54), \(\mathbf{A}\)\({}=\sum_n(\mathbf{x}_n-\bar{\mathbf{x}})(\mathbf{x}_n-\bar{\mathbf{x}})^\T\mathbf{W}\mathbf{M}^{-1}\)\({}=N\mathbf{S}\mathbf{W}\mathbf{M}^{-1}\), and by (12.55), \(\mathbf{B}\)\({}=N\sigma^2\mathbf{M}^{-1}+\mathbf{M}^{-1}\mathbf{W}^\T(N\mathbf{S})\mathbf{W}\mathbf{M}^{-1}\). Hence \(\mathbf{B}\)\({}=N(\sigma^2\mathbf{I}+\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{S}\mathbf{W})\mathbf{M}^{-1}\), so \(\mathbf{B}^{-1}\)\({}=N^{-1}\mathbf{M}(\sigma^2\mathbf{I}+\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{S}\mathbf{W})^{-1}\) and \(\mathbf{W}_{\rm new}\)\({}=\mathbf{A}\mathbf{B}^{-1}\)\({}=\mathbf{S}\mathbf{W}\mathbf{M}^{-1}\mathbf{M}(\sigma^2\mathbf{I}+\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{S}\mathbf{W})^{-1}\), as stated. For \(\sigma^2\): since \(\mathbf{W}_{\rm new}\mathbf{B}\)\({}=\mathbf{A}\), the last term of (12.57) summed over \(n\) is \(\Tr(\mathbf{W}_{\rm new}^\T\mathbf{W}_{\rm new}\mathbf{B})\)\({}=\Tr(\mathbf{W}_{\rm new}^\T\mathbf{A})\), and the middle term is \(-2\Tr(\mathbf{W}_{\rm new}^\T\mathbf{A})\). So the braces sum to \(\sum_n\lVert\mathbf{x}_n-\bar{\mathbf{x}}\rVert^2-\Tr(\mathbf{W}_{\rm new}^\T\mathbf{A})\)\({}=N\Tr\mathbf{S}-N\Tr(\mathbf{W}_{\rm new}^\T\mathbf{S}\mathbf{W}\mathbf{M}^{-1})\), which is both stated forms. (The second needs only running sums, so the algorithm can be run in a single pass over the data.)

Fixed points. Suppose \(\mathbf{S}\mathbf{C}^{-1}\mathbf{W}\)\({}=\mathbf{W}\), i.e. (Result 12.6(d)) \(\mathbf{S}\mathbf{W}\)\({}=\mathbf{W}\mathbf{M}\). Then \(\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{S}\mathbf{W}\)\({}=\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{W}\mathbf{M}\)\({}=\mathbf{M}^{-1}(\mathbf{M}-\sigma^2\mathbf{I})\mathbf{M}\)\({}=\mathbf{M}-\sigma^2\mathbf{I}\), so \(\mathbf{W}_{\rm new}\)\({}=\mathbf{W}\mathbf{M}\mathbf{M}^{-1}\)\({}=\mathbf{W}\). Conversely, let \(\mathbf{W}_{\rm new}=\mathbf{W}\) and put \(\mathbf{P}\)\({}=\mathbf{M}^{-1}\mathbf{W}^\T\mathbf{S}\mathbf{W}\); then \(\mathbf{S}\mathbf{W}\)\({}=\mathbf{W}(\sigma^2\mathbf{I}+\mathbf{P})\), so \(\mathbf{W}^\T\mathbf{S}\mathbf{W}\)\({}=(\mathbf{M}-\sigma^2\mathbf{I})(\sigma^2\mathbf{I}+\mathbf{P})\) and \(\mathbf{P}\)\({}=(\mathbf{I}-\sigma^2\mathbf{M}^{-1})(\sigma^2\mathbf{I}+\mathbf{P})\). Expanding, \(\mathbf{0}\)\({}=\sigma^2\mathbf{I}-\sigma^4\mathbf{M}^{-1}-\sigma^2\mathbf{M}^{-1}\mathbf{P}\), i.e. \(\mathbf{M}\)\({}=\sigma^2\mathbf{I}+\mathbf{P}\), so \(\mathbf{S}\mathbf{W}\)\({}=\mathbf{W}\mathbf{M}\): the stationarity condition. At such a point \(\mathbf{S}\mathbf{W}\mathbf{M}^{-1}\mathbf{W}^\T\)\({}=\mathbf{W}\mathbf{W}^\T\), so the \(\sigma^2\) fixed-point equation reads \(D\sigma^2\)\({}=\Tr\mathbf{S}-\Tr(\mathbf{W}\mathbf{W}^\T)\)\({}=\Tr\mathbf{S}-\sum_{j\in\mathcal{I}}(\lambda_j-\sigma^2)\), i.e. \((D-q)\sigma^2\)\({}=\sum_{i\notin\mathcal{I}}\lambda_i\), as in Step 3 of Result 12.8.

Monotonicity. This is the general property of EM (Section 9.4): the E step makes the lower bound touch the log likelihood, and the M step (which here maximises (12.53) exactly) cannot decrease the bound. ∎

In our numerical check (8-dimensional synthetic data, \(M=3\)), EM from a random start reached the closed-form maximum of Result 12.8 to eight significant figures, with \(\sigma^2\) equal to the average discarded eigenvalue; the saddle points are unstable fixed points, so EM started at random does not stop there.

Cost. A full eigendecomposition of \(\mathbf{S}\) costs \(O(D^3)\) and the leading \(M\) eigenpairs \(O(MD^2)\), but forming \(\mathbf{S}\) already costs \(O(ND^2)\). The snapshot method of Section 12.1.4 (Sirovich, 1987) avoids \(\mathbf{S}\) but costs \(O(N^3)\). An EM iteration is dominated by the sums over the data, \(O(NDM)\), which for \(M\ll D\) can more than compensate for needing several iterations. Because the E step works one point at a time and the M step needs only the running sums \(\mathbf{A}\), \(\mathbf{B}\) and \(\sum_n\lVert\mathbf{x}_n-\bar{\mathbf{x}}\rVert^2\) (Result 12.11), the algorithm can be run on-line, reading each \(D\)-dimensional point once.

Missing data

Because PPCA is a probability model, values that are missing at random can be handled by marginalising over them, again with EM (Exercise 12.16). Figure 12.11 shows the result on the oil flow data with 30% of the values deleted.

The zero-noise limit: EM for PCA

As \(\sigma^2\to0\) PPCA becomes standard PCA, and EM still gives a valid algorithm (Roweis, 1998). By (12.55) only \(\E[\mathbf{z}_n]\) is needed, and \(\mathbf{M}=\mathbf{W}^\T\mathbf{W}\). Let \(\widetilde{\mathbf{X}}\) be the \(N\times D\) centred data matrix and \(\bm{\Omega}\) the \(M\times N\) matrix whose \(n\)th column is \(\E[\mathbf{z}_n]\) (erratum: the first printing says "\(D\times M\) whose \(n\)th row"). The E step (12.54) and the M step (12.56) become

\[ \bm{\Omega}=(\mathbf{W}_{\rm old}^\T\mathbf{W}_{\rm old})^{-1}\mathbf{W}_{\rm old}^\T\widetilde{\mathbf{X}}^\T , \tag{12.58} \] \[ \mathbf{W}_{\rm new}=\widetilde{\mathbf{X}}^\T\bm{\Omega}^\T(\bm{\Omega}\bm{\Omega}^\T)^{-1} \tag{12.59} \]

(erratum: the first printing has \(\widetilde{\mathbf{X}}\) instead of \(\widetilde{\mathbf{X}}^\T\) in (12.58)). The E step projects the data orthogonally onto the current subspace; the M step re-fits the subspace by least squares with the projections held fixed (Exercise 12.17). Both are on-line friendly.

Proof

(a) As \(\sigma^2\to0\), \(\mathbf{M}\)\({}\to\mathbf{W}^\T\mathbf{W}\) (invertible for full-rank \(\mathbf{W}\)), so \(\E[\mathbf{z}_n]\)\({}\to(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T(\mathbf{x}_n-\bar{\mathbf{x}})\), whose stacking as columns is (12.58), and \(\E[\mathbf{z}_n\mathbf{z}_n^\T]\)\({}\to\E[\mathbf{z}_n]\E[\mathbf{z}_n]^\T\). Hence \(\sum_n\E[\mathbf{z}_n\mathbf{z}_n^\T]\)\({}\to\bm{\Omega}\bm{\Omega}^\T\) and \(\sum_n(\mathbf{x}_n-\bar{\mathbf{x}})\E[\mathbf{z}_n]^\T\)\({}=\widetilde{\mathbf{X}}^\T\bm{\Omega}^\T\), and (12.56) becomes (12.59).

(b) Substitute (12.58) into (12.59), using \(\widetilde{\mathbf{X}}^\T\widetilde{\mathbf{X}}\)\({}=N\mathbf{S}\) and \(\mathbf{G}=\mathbf{W}^\T\mathbf{W}\): \(\widetilde{\mathbf{X}}^\T\bm{\Omega}^\T\)\({}=N\mathbf{S}\mathbf{W}\mathbf{G}^{-1}\) and \(\bm{\Omega}\bm{\Omega}^\T\)\({}=N\mathbf{G}^{-1}\mathbf{W}^\T\mathbf{S}\mathbf{W}\mathbf{G}^{-1}\), so \(\mathbf{W}_{\rm new}\)\({}=\mathbf{S}\mathbf{W}\mathbf{G}^{-1}\mathbf{G}(\mathbf{W}^\T\mathbf{S}\mathbf{W})^{-1}\mathbf{G}\), as stated. The right factor is invertible, so the column space of \(\mathbf{W}_{\rm new}\) is that of \(\mathbf{S}\mathbf{W}\). For convergence write \(\mathbf{W}_0\)\({}=\mathbf{U}_M\mathbf{A}+\mathbf{U}_\perp\mathbf{B}\) with \(\mathbf{A}\)\({}=\mathbf{U}_M^\T\mathbf{W}_0\) invertible and \(\mathbf{U}_\perp\) the remaining eigenvectors. Then \(\mathbf{S}^t\mathbf{W}_0\)\({}=\big(\mathbf{U}_M+\mathbf{U}_\perp\mathbf{E}_t\big)\mathbf{L}_M^t\mathbf{A}\) with \(\mathbf{E}_t\)\({}=\mathbf{L}_\perp^t\mathbf{B}\mathbf{A}^{-1}\mathbf{L}_M^{-t}\), whose norm is at most \((\lambda_{M+1}/\lambda_M)^t\lVert\mathbf{B}\mathbf{A}^{-1}\rVert\)\({}\to0\); the column space is that of \(\mathbf{U}_M+\mathbf{U}_\perp\mathbf{E}_t\), which tends to that of \(\mathbf{U}_M\). This is orthogonal (subspace) iteration, the block version of the power method (Golub and Van Loan, 1996).

(c) \(\lVert\widetilde{\mathbf{X}}-\bm{\Omega}^\T\mathbf{W}^\T\rVert_F^2\)\({}=\sum_n\lVert\mathbf{x}_n-\bar{\mathbf{x}}-\mathbf{W}\bm{\omega}_n\rVert^2\). For fixed \(\mathbf{W}\) each \(\bm{\omega}_n\) solves a least-squares problem with normal equations \(\mathbf{W}^\T\mathbf{W}\bm{\omega}_n\)\({}=\mathbf{W}^\T(\mathbf{x}_n-\bar{\mathbf{x}})\), which is (12.58). For fixed \(\bm{\Omega}\) the problem in \(\mathbf{W}^\T\) is a multivariate regression of \(\widetilde{\mathbf{X}}\) on \(\bm{\Omega}^\T\), with normal equations \(\bm{\Omega}\bm{\Omega}^\T\mathbf{W}^\T\)\({}=\bm{\Omega}\widetilde{\mathbf{X}}\), which is (12.59). Each step minimises the cost over its block, so the cost cannot increase (see also Exercise 12.17). ∎

Bishop's physical picture (\(D=2\), \(M=1\)): the subspace is a rod and each point is tied to it by a spring obeying Hooke's law, with energy proportional to the squared length. In the E step the rod is held fixed and the attachment points slide along it to minimise the energy, ending at the orthogonal projections. In the M step the attachment points are frozen and the rod is released to its minimum-energy position. Figure 12.12 shows the iterations.

12.2.3 Bayesian PCA

So far \(M\) was given. For visualisation one takes \(M=2\), but in general the right value is unclear. One can look for a gap in the eigenvalue spectrum (Figure 12.4), but often there is none. Cross-validation on the likelihood is possible (PPCA has one) but costly, especially for a mixture of PPCA models where each component needs its own dimension (Tipping and Bishop, 1999a). A Bayesian approach would integrate out \(\bm{\mu}\), \(\mathbf{W}\) and \(\sigma^2\) and compare the marginal likelihoods of different \(M\); this is intractable but can be approximated variationally (Bishop, 1999b). The book presents a simpler method based on the evidence approximation (Bishop, 1999a), appropriate when \(N\) is large and the posterior is sharply peaked. It uses automatic relevance determination (ARD, Section 7.2.2): each column \(\mathbf{w}_i\) of \(\mathbf{W}\) gets its own Gaussian prior with precision \(\alpha_i\),

\[ p(\mathbf{W}\given\bm{\alpha})=\prod_{i=1}^M\Big(\frac{\alpha_i}{2\pi}\Big)^{D/2}\exp\Big\{-\frac12\alpha_i\mathbf{w}_i^\T\mathbf{w}_i\Big\}. \tag{12.60} \]

The hyperparameters are chosen by maximising the marginal likelihood

\[ p(\mathbf{X}\given\bm{\alpha},\bm{\mu},\sigma^2)=\int p(\mathbf{X}\given\mathbf{W},\bm{\mu},\sigma^2)\,p(\mathbf{W}\given\bm{\alpha})\,\d\mathbf{W}, \tag{12.61} \]

with \(\bm{\mu}\) and \(\sigma^2\) treated as parameters (no priors), and \(\ln p(\mathbf{X}\given\mathbf{W},\bm{\mu},\sigma^2)\) given by (12.43). Some \(\alpha_i\) go to infinity, forcing \(\mathbf{w}_i\to\mathbf{0}\) (its posterior collapses onto the origin): a sparse solution. The effective dimension of the principal subspace is the number of finite \(\alpha_i\), and the surviving columns are the "relevant" ones. The Bayesian criterion trades the improved fit from more columns against the complexity they add, as in the relevance vector machine (Section 7.2). The integral is intractable, and with a Laplace approximation and a sharply peaked posterior the re-estimation equation becomes

\[ \alpha_i^{\rm new}=\frac{D}{\mathbf{w}_i^\T\mathbf{w}_i} , \tag{12.62} \]

the analogue of (3.98) with \(D\) parameters per column. These updates are interleaved with EM for \(\mathbf{W}\) and \(\sigma^2\): the E step is still (12.54)–(12.55), the \(\sigma^2\) update is still (12.57), and only the \(\mathbf{W}\) update changes, to

\[ \mathbf{W}_{\rm new}=\Big[\sum_{n=1}^N(\mathbf{x}_n-\bar{\mathbf{x}})\E[\mathbf{z}_n]^\T\Big]\Big[\sum_{n=1}^N\E[\mathbf{z}_n\mathbf{z}_n^\T]+\sigma^2\mathbf{A}\Big]^{-1},\qquad\mathbf{A}=\diag(\alpha_i). \tag{12.63} \]

With \(M=D-1\) the model can represent everything from a full-covariance Gaussian (all \(\alpha_i\) finite) to an isotropic one (all \(\alpha_i\to\infty\)); a smaller \(M\) saves computation but caps the dimension that can be found. A caution about what the iteration optimises: substituting (12.62) into \(F\) of Gap 12.6 gives \(F\)\({}=\ln p(\mathbf{X}\given\mathbf{W})-D\sum_i\ln\lVert\mathbf{w}_i\rVert+\const\), and the second term tends to \(+\infty\) as a column shrinks to zero, so the joint objective has no finite maximum and the scheme should be read as a local fixed-point iteration. Columns strongly supported by the data settle at a local maximum away from zero; weakly supported ones are pulled to zero, which is the pruning. The fully Bayesian variational treatment (Bishop, 1999b) avoids this degeneracy.

Gibbs sampling for Bayesian PCA

Bayesian PCA also illustrates Gibbs sampling (Section 11.3). Figure 12.15 follows \(\ln\alpha_i\) for data in \(D=4\) dimensions with \(M=3\), generated from a PPCA model with one direction of high variance and low-variance noise elsewhere. The posterior has three modes: at any time one \(\alpha_i\) is small (its column carries the strong direction) and the other two are large (their columns are switched off), and the sampler jumps abruptly between the modes. The three modes are images of each other under permutation of the columns of \(\mathbf{W}\), which leaves the likelihood and the prior unchanged; this is the continuous-latent analogue of label switching in mixtures.

Derivation of the conditionals

Each conditional is proportional to the joint density

\[ \prod_n\N(\mathbf{x}_n\given\mathbf{W}\mathbf{z}_n+\bm{\mu},\tau^{-1}\mathbf{I})\,\N(\mathbf{z}_n\given\mathbf{0},\mathbf{I})\times p(\mathbf{W}\given\bm{\alpha})\prod_ip(\alpha_i)\,p(\tau)\,p(\bm{\mu}), \]

viewed as a function of one block. (1) The \(\mathbf{z}_n\) are conditionally independent and each has the PPCA posterior (Result 12.6). (2) The prior (12.60) makes all entries \(W_{di}\) independent with variance \(\alpha_i^{-1}\), so each row has prior \(\N(\widetilde{\mathbf{w}}_d\given\mathbf{0},\mathbf{A}^{-1})\), and given \(\mathbf{Z}\) the likelihood factorises over rows: \(\prod_n\N(x_{nd}\given\widetilde{\mathbf{w}}_d^\T\mathbf{z}_n+\mu_d,\tau^{-1})\). Row \(d\) is therefore a Bayesian linear regression of \(x_{nd}-\mu_d\) on \(\mathbf{z}_n\), and completing the square gives the stated Gaussian (as in (3.53)–(3.54)). (3) As a function of \(\alpha_i\), \(\alpha_i^{D/2}e^{-\alpha_i\lVert\mathbf{w}_i\rVert^2/2}\times\alpha_i^{a_\alpha-1}e^{-b_\alpha\alpha_i}\) is a Gamma kernel (toolkit). (4) Likewise \(\tau^{ND/2}e^{-\tau Q/2}\tau^{a_\tau-1}e^{-b_\tau\tau}\) with \(Q\)\({}=\sum_n\lVert\mathbf{x}_n-\mathbf{W}\mathbf{z}_n-\bm{\mu}\rVert^2\). (5) As a function of \(\bm{\mu}\), the exponent is \(-\frac{\tau}{2}\sum_n\lVert\mathbf{x}_n-\mathbf{W}\mathbf{z}_n-\bm{\mu}\rVert^2-\frac\beta2\lVert\bm{\mu}\rVert^2\), a Gaussian with precision \((\beta+N\tau)\mathbf{I}\) and mean \(\tau\sum_n(\mathbf{x}_n-\mathbf{W}\mathbf{z}_n)/(\beta+N\tau)\). ∎

The model here puts a prior on \(\mathbf{W}\) only; a fully Bayesian treatment with priors on \(\bm{\mu}\), \(\sigma^2\) and \(\bm{\alpha}\), fitted variationally, is given by Bishop (1999b), and Minka (2001c) discusses Bayesian choices of the dimension of PCA.

12.2.4 Factor analysis

The only difference from PPCA is a diagonal instead of an isotropic noise covariance. As in PPCA, the observed variables are independent given \(\mathbf{z}\): the model explains the covariance between variables through \(\mathbf{W}\) and gives each variable its own independent noise variance in \(\bm{\Psi}\). Factor analysis is as old as PCA (Everitt, 1984; Bartholomew, 1987; Basilevsky, 1994). Lawley (1953) and Anderson (1963) showed that at the stationary points of the likelihood of a factor model with \(\bm{\Psi}=\sigma^2\mathbf{I}\) the columns of \(\mathbf{W}\) are scaled eigenvectors of \(\mathbf{S}\) and \(\sigma^2\) is the average discarded eigenvalue; Tipping and Bishop (1999b) showed that the maximum uses the principal eigenvectors (Result 12.8).

Factor analysis has the same rotational invariance as PPCA (Exercise 12.19): with \(\widetilde{\mathbf{W}}\)\({}=\mathbf{W}\mathbf{R}\), \(\mathbf{R}\) orthogonal, \(\widetilde{\mathbf{W}}\widetilde{\mathbf{W}}^\T+\bm{\Psi}\)\({}=\mathbf{W}\mathbf{W}^\T+\bm{\Psi}\), so the likelihood cannot tell \(\mathbf{W}\) from \(\mathbf{W}\mathbf{R}\). This is why attempts to interpret individual factors (the coordinates of \(\mathbf{z}\)) have been controversial. Here the model is a density model, in which the latent subspace matters but its coordinates do not; removing the rotational degeneracy requires non-Gaussian latent variables, which leads to ICA (Section 12.4.1).

The maximum likelihood estimate of \(\bm{\mu}\) is again \(\bar{\mathbf{x}}\) (Exercise 12.20; the argument of Result 12.7 applies verbatim because \(\mathbf{C}\)\({}=\mathbf{W}\mathbf{W}^\T+\bm{\Psi}\) is positive definite). Unlike PPCA there is no closed form for \(\mathbf{W}\) and \(\bm{\Psi}\), and EM is used (Rubin and Thayer, 1982). The E step is

\[ \E[\mathbf{z}_n]=\mathbf{G}\mathbf{W}^\T\bm{\Psi}^{-1}(\mathbf{x}_n-\bar{\mathbf{x}}), \tag{12.66} \] \[ \E[\mathbf{z}_n\mathbf{z}_n^\T]=\mathbf{G}+\E[\mathbf{z}_n]\E[\mathbf{z}_n]^\T , \tag{12.67} \] \[ \mathbf{G}=(\mathbf{I}+\mathbf{W}^\T\bm{\Psi}^{-1}\mathbf{W})^{-1}, \tag{12.68} \]

which needs only \(M\times M\) inversions (and the trivial inverse of the diagonal \(\bm{\Psi}\)). The M step is

\[ \mathbf{W}_{\rm new}=\Big[\sum_{n=1}^N(\mathbf{x}_n-\bar{\mathbf{x}})\E[\mathbf{z}_n]^\T\Big]\Big[\sum_{n=1}^N\E[\mathbf{z}_n\mathbf{z}_n^\T]\Big]^{-1}, \tag{12.69} \] \[ \bm{\Psi}_{\rm new}=\diag\Big\{\mathbf{S}-\mathbf{W}_{\rm new}\frac1N\sum_{n=1}^N\E[\mathbf{z}_n](\mathbf{x}_n-\bar{\mathbf{x}})^\T\Big\}, \tag{12.70} \]

where \(\diag\{\cdot\}\) sets the off-diagonal elements to zero.

A Bayesian factor analysis follows from the techniques of the book (for example ARD priors on the loadings, as in Section 12.2.3). The number of free parameters of the factor analysis covariance is \(DM+D-M(M-1)/2\) (Exercise 12.18): the single \(\sigma^2\) of PPCA is replaced by \(D\) uniquenesses, and the rotations are still redundant.

The two models differ in how they behave when the data are transformed (Exercise 12.25). Fitting PCA or PPCA to rotated data gives the same fit with \(\mathbf{W}\) rotated accordingly, but rescaling the variables (changing units) changes the principal directions. Factor analysis has the opposite property: a component-wise rescaling of the data is absorbed by rescaling the rows of \(\mathbf{W}\) and the entries of \(\bm{\Psi}\).

Proof

By the change-of-variables formula for a Gaussian, \(\N(\mathbf{A}\mathbf{x}\given\bm{\mu}',\mathbf{C}')\)\({}=\lvert\det\mathbf{A}\rvert^{-1}\N(\mathbf{x}\given\mathbf{A}^{-1}\bm{\mu}',\mathbf{A}^{-1}\mathbf{C}'\mathbf{A}^{-\T})\). Hence the log likelihood of the transformed data at parameters \((\bm{\mu}',\mathbf{W}',\bm{\Phi}')\) equals the log likelihood of the original data at \((\mathbf{A}^{-1}\bm{\mu}',\mathbf{A}^{-1}\mathbf{W}',\mathbf{A}^{-1}\bm{\Phi}'\mathbf{A}^{-\T})\), minus the constant \(N\ln\lvert\det\mathbf{A}\rvert\) (note \(\mathbf{A}^{-1}(\mathbf{W}'\mathbf{W}'^\T+\bm{\Phi}')\mathbf{A}^{-\T}\)\({}=(\mathbf{A}^{-1}\mathbf{W}')(\mathbf{A}^{-1}\mathbf{W}')^\T+\mathbf{A}^{-1}\bm{\Phi}'\mathbf{A}^{-\T}\)). If the class of \(\bm{\Phi}\) is preserved by \(\bm{\Phi}\mapsto\mathbf{A}\bm{\Phi}\mathbf{A}^\T\) and by its inverse, the map \((\bm{\mu},\mathbf{W},\bm{\Phi})\mapsto(\mathbf{A}\bm{\mu},\mathbf{A}\mathbf{W},\mathbf{A}\bm{\Phi}\mathbf{A}^\T)\) is a bijection of the parameter space that turns one likelihood into the other up to a constant, so it maps maximisers to maximisers. (i) If \(\mathbf{A}\) and \(\bm{\Phi}\) are diagonal, so is \(\mathbf{A}\bm{\Phi}\mathbf{A}^\T\) (entries \(a_d^2\psi_d\)). (ii) If \(\mathbf{A}\mathbf{A}^\T\)\({}=c\,\mathbf{I}\), then \(\mathbf{A}\sigma^2\mathbf{I}\mathbf{A}^\T\)\({}=c\sigma^2\mathbf{I}\). In general neither property holds for the other model: a diagonal rescaling of isotropic noise is not isotropic, and a rotation of diagonal noise is not diagonal. (The Svensén solution to this exercise writes \(\mathbf{A}\bm{\Phi}^{-1}\mathbf{A}^\T\); the correct transformed noise covariance is \(\mathbf{A}\bm{\Phi}\mathbf{A}^\T\).) ∎

12.3 Kernel PCA

Chapter 6 showed that an algorithm that uses the data only through inner products \(\mathbf{x}^\T\mathbf{x}'\) can be made nonlinear by replacing them with a kernel \(k(\mathbf{x},\mathbf{x}')\). Applied to PCA this gives kernel PCA (Schölkopf, Smola and Müller, 1998). Assume for now that the data are centred, \(\sum_n\mathbf{x}_n=\mathbf{0}\). The principal components solve

\[ \mathbf{S}\mathbf{u}_i=\lambda_i\mathbf{u}_i , \tag{12.71} \] \[ \mathbf{S}=\frac1N\sum_{n=1}^N\mathbf{x}_n\mathbf{x}_n^\T , \tag{12.72} \]

with \(\mathbf{u}_i^\T\mathbf{u}_i=1\). Now map each point into a feature space, \(\mathbf{x}_n\mapsto\bm{\phi}(\mathbf{x}_n)\), and do PCA there; linear projections in feature space become nonlinear functions of \(\mathbf{x}\) (Figure 12.16 of the book). Notation warning: in this section \(M\) is the dimension of the feature space (possibly infinite) and \(\mathbf{C}\) is the feature-space covariance, not the PPCA quantities of Section 12.2. Assuming for the moment that the features are centred too, \(\sum_n\bm{\phi}(\mathbf{x}_n)\)\({}=\mathbf{0}\),

\[ \mathbf{C}=\frac1N\sum_{n=1}^N\bm{\phi}(\mathbf{x}_n)\bm{\phi}(\mathbf{x}_n)^\T , \tag{12.73} \] \[ \mathbf{C}\mathbf{v}_i=\lambda_i\mathbf{v}_i\quad(i=1,\dots,M), \tag{12.74} \]

that is,

\[ \frac1N\sum_{n=1}^N\bm{\phi}(\mathbf{x}_n)\big\{\bm{\phi}(\mathbf{x}_n)^\T\mathbf{v}_i\big\}=\lambda_i\mathbf{v}_i . \tag{12.75} \]

If \(\lambda_i\gt0\), (12.75) expresses \(\mathbf{v}_i\) as a combination of the feature vectors,

\[ \mathbf{v}_i=\sum_{n=1}^Na_{in}\bm{\phi}(\mathbf{x}_n),\qquad a_{in}=\frac{\bm{\phi}(\mathbf{x}_n)^\T\mathbf{v}_i}{N\lambda_i} . \tag{12.76} \]

Substituting back,

\[ \frac1N\sum_{n=1}^N\bm{\phi}(\mathbf{x}_n)\bm{\phi}(\mathbf{x}_n)^\T\sum_{m=1}^Na_{im}\bm{\phi}(\mathbf{x}_m)=\lambda_i\sum_{n=1}^Na_{in}\bm{\phi}(\mathbf{x}_n), \tag{12.77} \]

and multiplying by \(\bm{\phi}(\mathbf{x}_l)^\T\) turns every feature vector into a kernel value \(k(\mathbf{x}_n,\mathbf{x}_m)\)\({}=\bm{\phi}(\mathbf{x}_n)^\T\bm{\phi}(\mathbf{x}_m)\):

\[ \frac1N\sum_{n=1}^Nk(\mathbf{x}_l,\mathbf{x}_n)\sum_{m=1}^Na_{im}k(\mathbf{x}_n,\mathbf{x}_m)=\lambda_i\sum_{n=1}^Na_{in}k(\mathbf{x}_l,\mathbf{x}_n) \tag{12.78} \]

(erratum: the first printing has the upper limit \(m\) instead of \(N\) in the inner sum). With the Gram matrix \(\mathbf{K}\), \(K_{nm}\)\({}=k(\mathbf{x}_n,\mathbf{x}_m)\), and \(\mathbf{a}_i\)\({}=(a_{i1},\dots,a_{iN})^\T\) (erratum: "\(a_{ni}\)" should be "\(a_{in}\)"), this reads

\[ \mathbf{K}^2\mathbf{a}_i=\lambda_iN\mathbf{K}\mathbf{a}_i , \tag{12.79} \]

and it suffices to solve

\[ \mathbf{K}\mathbf{a}_i=\lambda_iN\mathbf{a}_i . \tag{12.80} \]

The coefficients are normalised by requiring unit feature-space eigenvectors,

\[ 1=\mathbf{v}_i^\T\mathbf{v}_i=\sum_{n=1}^N\sum_{m=1}^Na_{in}a_{im}\bm{\phi}(\mathbf{x}_n)^\T\bm{\phi}(\mathbf{x}_m)=\mathbf{a}_i^\T\mathbf{K}\mathbf{a}_i=\lambda_iN\mathbf{a}_i^\T\mathbf{a}_i , \tag{12.81} \]

and the projection of a point \(\mathbf{x}\) onto the \(i\)th component needs only kernel evaluations:

\[ y_i(\mathbf{x})=\bm{\phi}(\mathbf{x})^\T\mathbf{v}_i=\sum_{n=1}^Na_{in}\bm{\phi}(\mathbf{x})^\T\bm{\phi}(\mathbf{x}_n)=\sum_{n=1}^Na_{in}k(\mathbf{x},\mathbf{x}_n). \tag{12.82} \]

Proof

(a) \(\mathbf{C}\mathbf{v}\)\({}=N^{-1}\bm{\Phi}^\T\bm{\Phi}\bm{\Phi}^\T\mathbf{a}\)\({}=N^{-1}\bm{\Phi}^\T\mathbf{K}\mathbf{a}\)\({}=(\kappa/N)\bm{\Phi}^\T\mathbf{a}\)\({}=(\kappa/N)\mathbf{v}\), and \(\mathbf{v}\ne\mathbf{0}\) because \(\lVert\mathbf{v}\rVert^2\)\({}=\mathbf{a}^\T\mathbf{K}\mathbf{a}\)\({}=\kappa\lVert\mathbf{a}\rVert^2\)\({}\gt0\). Conversely, if \(\mathbf{C}\mathbf{v}\)\({}=\lambda\mathbf{v}\) with \(\lambda\gt0\), then \(\mathbf{v}=\bm{\Phi}^\T\mathbf{a}\) with \(\mathbf{a}\)\({}=\bm{\Phi}\mathbf{v}/(N\lambda)\), which is (12.76); multiplying \(\mathbf{C}\mathbf{v}\)\({}=\lambda\mathbf{v}\) on the left by \(\bm{\Phi}\) gives \(N^{-1}\mathbf{K}\bm{\Phi}\mathbf{v}\)\({}=\lambda\bm{\Phi}\mathbf{v}\), i.e. \(\mathbf{K}\mathbf{a}\)\({}=N\lambda\mathbf{a}\) (with \(\bm{\Phi}\mathbf{v}\ne\mathbf{0}\) because \(\bm{\Phi}^\T\bm{\Phi}\mathbf{v}\)\({}=N\lambda\mathbf{v}\ne\mathbf{0}\)). This is the argument of Result 12.5 with \(\bm{\Phi}\) in place of \(\mathbf{X}\).

(b) \(\lVert\mathbf{v}\rVert^2\)\({}=\mathbf{a}^\T\mathbf{K}\mathbf{a}\)\({}=\kappa\mathbf{a}^\T\mathbf{a}\), which is (12.81).

(c) Any solution of (12.80) solves (12.79) (multiply by \(\mathbf{K}\)). Conversely, \(\mathbf{K}\) is symmetric, so \(\R^N\) is the orthogonal sum of its range and its null space; write \(\mathbf{a}\)\({}=\mathbf{a}_r+\mathbf{a}_0\) accordingly. (12.79) says \(\mathbf{K}(\mathbf{K}\mathbf{a}-\lambda N\mathbf{a})\)\({}=\mathbf{0}\), i.e. \(\mathbf{K}\mathbf{a}-\lambda N\mathbf{a}\)\({}=(\mathbf{K}\mathbf{a}_r-\lambda N\mathbf{a}_r)-\lambda N\mathbf{a}_0\) lies in the null space. The first bracket lies in the range, so it must vanish: \(\mathbf{a}_r\) solves (12.80), and \(\mathbf{a}_0\) is arbitrary. Since \(\lVert\bm{\Phi}^\T\mathbf{a}_0\rVert^2\)\({}=\mathbf{a}_0^\T\mathbf{K}\mathbf{a}_0\)\({}=0\), \(\bm{\Phi}^\T\mathbf{a}_0\)\({}=\mathbf{0}\): the feature vector \(\mathbf{v}=\bm{\Phi}^\T\mathbf{a}\) and hence every projection \(y(\mathbf{x})\)\({}=\bm{\phi}(\mathbf{x})^\T\mathbf{v}\) are unchanged.

(d) The rank of \(\mathbf{C}\)\({}=N^{-1}\bm{\Phi}^\T\bm{\Phi}\) equals the rank of \(\mathbf{K}\), at most \(N\), however large \(M\) is; but \(M\) can exceed \(D\), so the number of non-zero components can exceed the \(D\) of linear PCA. ∎

So far the features were assumed centred. In general they are not, and the mean \(\frac1N\sum_l\bm{\phi}(\mathbf{x}_l)\) cannot be subtracted explicitly when the feature space is not available. The centred features

\[ \widetilde{\bm{\phi}}(\mathbf{x}_n)=\bm{\phi}(\mathbf{x}_n)-\frac1N\sum_{l=1}^N\bm{\phi}(\mathbf{x}_l) \tag{12.83} \]

have Gram matrix

\[ \begin{aligned}\widetilde{K}_{nm}&=\widetilde{\bm{\phi}}(\mathbf{x}_n)^\T\widetilde{\bm{\phi}}(\mathbf{x}_m)\\ &=\bm{\phi}(\mathbf{x}_n)^\T\bm{\phi}(\mathbf{x}_m)-\frac1N\sum_{l=1}^N\bm{\phi}(\mathbf{x}_n)^\T\bm{\phi}(\mathbf{x}_l)\\ &\qquad-\frac1N\sum_{l=1}^N\bm{\phi}(\mathbf{x}_l)^\T\bm{\phi}(\mathbf{x}_m)+\frac{1}{N^2}\sum_{j=1}^N\sum_{l=1}^N\bm{\phi}(\mathbf{x}_j)^\T\bm{\phi}(\mathbf{x}_l)\\ &=k(\mathbf{x}_n,\mathbf{x}_m)-\frac1N\sum_{l=1}^Nk(\mathbf{x}_l,\mathbf{x}_m)\\ &\qquad-\frac1N\sum_{l=1}^Nk(\mathbf{x}_n,\mathbf{x}_l)+\frac{1}{N^2}\sum_{j=1}^N\sum_{l=1}^Nk(\mathbf{x}_j,\mathbf{x}_l),\end{aligned} \tag{12.84} \]

or in matrix form

\[ \widetilde{\mathbf{K}}=\mathbf{K}-\mathbf{1}_N\mathbf{K}-\mathbf{K}\mathbf{1}_N+\mathbf{1}_N\mathbf{K}\mathbf{1}_N , \tag{12.85} \]

where \(\mathbf{1}_N\) is the \(N\times N\) matrix with every entry \(1/N\).

Proof

In matrix form (12.83) is \(\widetilde{\bm{\Phi}}\)\({}=\mathbf{H}\bm{\Phi}\) (subtract the column means), so \(\widetilde{\mathbf{K}}\)\({}=\widetilde{\bm{\Phi}}\widetilde{\bm{\Phi}}^\T\)\({}=\mathbf{H}\mathbf{K}\mathbf{H}\) since \(\mathbf{H}\) is symmetric; expanding \((\mathbf{I}-\mathbf{1}_N)\mathbf{K}(\mathbf{I}-\mathbf{1}_N)\) gives (12.85), whose \((n,m)\) entry is (12.84) because \((\mathbf{1}_N\mathbf{K})_{nm}\)\({}=\frac1N\sum_lK_{lm}\), \((\mathbf{K}\mathbf{1}_N)_{nm}\)\({}=\frac1N\sum_lK_{nl}\) and \((\mathbf{1}_N\mathbf{K}\mathbf{1}_N)_{nm}\)\({}=\frac{1}{N^2}\sum_{j,l}K_{jl}\). Result 12.14 applied to the centred features gives the eigenproblem and normalisation. For a point \(\mathbf{x}\), the projection onto \(\mathbf{v}_i\)\({}=\sum_na_{in}\widetilde{\bm{\phi}}(\mathbf{x}_n)\) of the centred feature \(\widetilde{\bm{\phi}}(\mathbf{x})\)\({}=\bm{\phi}(\mathbf{x})-\frac1N\sum_l\bm{\phi}(\mathbf{x}_l)\) is \(\sum_na_{in}\widetilde{\bm{\phi}}(\mathbf{x})^\T\widetilde{\bm{\phi}}(\mathbf{x}_n)\), and expanding the inner product exactly as in (12.84) gives \(\widetilde{k}\). Finally \(\mathbf{H}\mathbf{1}=\mathbf{0}\) gives \(\widetilde{\mathbf{K}}\mathbf{1}\)\({}=\mathbf{0}\). (Checked numerically with an explicit polynomial feature map: the kernel projections equal the feature-space projections.) ∎

With the linear kernel \(k(\mathbf{x},\mathbf{x}')\)\({}=\mathbf{x}^\T\mathbf{x}'\), kernel PCA is ordinary PCA (Exercise 12.27): then \(\widetilde{\mathbf{K}}\)\({}=\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^\T\), (12.80) is the \(N\times N\) problem (12.28), the normalised \(\mathbf{a}_i\) are \(\mathbf{v}_i/(N\lambda_i)^{1/2}\), and \(\sum_na_{in}(\mathbf{x}_n-\bar{\mathbf{x}})\)\({}=\widetilde{\mathbf{X}}^\T\mathbf{a}_i\)\({}=\mathbf{u}_i\) by (12.30), so \(y_i(\mathbf{x})\)\({}=\mathbf{u}_i^\T(\mathbf{x}-\bar{\mathbf{x}})\).

Figure 12.17 applies kernel PCA with the "Gaussian" kernel

\[ k(\mathbf{x},\mathbf{x}')=\exp\big(-\lVert\mathbf{x}-\mathbf{x}'\rVert^2/0.1\big) \tag{12.86} \]

to a synthetic two-dimensional data set, drawing contours along which the projection

\[ \bm{\phi}(\mathbf{x})^\T\mathbf{v}_i=\sum_{n=1}^Na_{in}k(\mathbf{x},\mathbf{x}_n) \tag{12.87} \]

is constant (in practice with the centred kernel \(\widetilde{k}\) of Result 12.15).

Limitations. Kernel PCA needs the eigenvectors of the \(N\times N\) matrix \(\widetilde{\mathbf{K}}\) rather than the \(D\times D\) matrix \(\mathbf{S}\), so for large data sets approximations are used. And while linear PCA can approximate a point by its projection onto the principal subspace,

\[ \widehat{\mathbf{x}}_n=\sum_{i=1}^L\big(\mathbf{x}_n^\T\mathbf{u}_i\big)\mathbf{u}_i , \tag{12.88} \]

kernel PCA in general cannot: \(\bm{\phi}\) maps the \(D\)-dimensional data space onto a \(D\)-dimensional curved surface in feature space, and the projection of \(\bm{\phi}(\mathbf{x})\) onto the kernel principal subspace usually does not lie on that surface, so it has no pre-image \(\mathbf{x}\). Approximate pre-images can be computed (Bakir, Weston and Schölkopf, 2004).

12.4 Nonlinear Latent Variable Models

The linear-Gaussian models of this chapter are important, easy to analyse and fit, and useful as building blocks. This section sketches models that are nonlinear, non-Gaussian, or both. The two generalisations are related: any density can be obtained from a fixed simple density by a nonlinear change of variables (Exercise 12.28). In one dimension, if \(q\) and \(p\) are densities with distribution functions \(F_q\) and \(F_p\) and \(q\gt0\) everywhere, the increasing map \(y=f(x)=F_p^{-1}(F_q(x))\) turns \(x\sim q\) into \(y\sim p\), because it satisfies the differential equation

\[ p\big(f(x)\big)\,f'(x)=q(x), \]

which is the change-of-variables formula (1.27) read backwards (see Exercise 12.28). Several practical latent variable models are built on this idea.

12.4.1 Independent component analysis

Keep the observed variables linear in the latent ones but let the latent distribution be non-Gaussian. Independent component analysis (ICA) takes a factorised latent distribution,

\[ p(\mathbf{z})=\prod_{j=1}^Mp(z_j). \tag{12.89} \]

The motivating example is the "cocktail party": two people speak at once and two microphones record linear mixtures of the two voices (ignoring delays and echoes). If the constant mixing coefficients can be inferred from the recordings alone, the mixing can be inverted to recover each voice separately. This is blind source separation: neither the sources nor the mixing coefficients are observed (Cardoso, 1998). Following MacKay (2003), ignore the time ordering and treat the samples as i.i.d.; there are as many latent variables (voices) as observed ones (microphones), \(M=D\), and no noise term is needed, because a square nonsingular mixing does not produce a singular distribution.

Derivation of the gradient

\(\d\ln\lvert\det\mathbf{B}\rvert\)\({}=\Tr(\mathbf{B}^{-1}\d\mathbf{B})\), so its gradient is \(\mathbf{B}^{-\T}\) (toolkit). For one term, \(\partial\ln p_j(\mathbf{b}_j^\T\mathbf{x})/\partial\mathbf{b}_j\)\({}=g_j(\mathbf{b}_j^\T\mathbf{x})\,\mathbf{x}\), which is row \(j\) of \(\mathbf{g}(\mathbf{B}\mathbf{x})\mathbf{x}^\T\). For the heavy-tailed choice (12.90) below, \(\ln p(z)=-\ln\pi-\ln\cosh z\) and \(g(z)=-\tanh z\), so the gradient is \(N\mathbf{B}^{-\T}-\sum_n\tanh(\mathbf{B}\mathbf{x}_n)\mathbf{x}_n^\T\) (checked by finite differences). Multiplying the gradient on the right by \(\mathbf{B}^\T\mathbf{B}\) gives the "natural gradient" direction \(\big[N\mathbf{I}+\sum_n\mathbf{g}(\mathbf{y}_n)\mathbf{y}_n^\T\big]\mathbf{B}\), \(\mathbf{y}_n\)\({}=\mathbf{B}\mathbf{x}_n\), which avoids the matrix inverse (Amari, Cichocki and Yang, 1996). ∎

The method needs non-Gaussian sources. In PPCA and factor analysis the latent distribution is a standard Gaussian, and the likelihood cannot distinguish \(\mathbf{W}\) from \(\mathbf{W}\mathbf{R}\) (12.39): the latent coordinates are determined only up to a rotation. A general Gaussian latent distribution does not help, because it is equivalent to the standard one (Exercise 12.4). Put differently, PCA rotates the data so that the coordinates are uncorrelated, and zero correlation is necessary for independence but not sufficient (Exercise 12.29). A common choice is the heavy-tailed density

\[ p(z_j)=\frac{1}{\pi\cosh(z_j)}=\frac{2}{\pi\big(e^{z_j}+e^{-z_j}\big)} \tag{12.90} \]

(erratum: the first printing has numerator 1 in the rightmost expression). It integrates to one (\(\int\operatorname{sech}z\,\d z=\pi\)), has variance \(\pi^2/4\) and heavier tails than a Gaussian (excess kurtosis 2), as many real signals do.

Proof of the Gaussian part

\(\mathbf{x}=\mathbf{A}\mathbf{z}\) is Gaussian with mean \(\mathbf{0}\) and covariance \(\mathbf{A}\mathbf{A}^\T\) (an affine image of a Gaussian), and \((\mathbf{A}\mathbf{R})(\mathbf{A}\mathbf{R})^\T\)\({}=\mathbf{A}\mathbf{R}\mathbf{R}^\T\mathbf{A}^\T\)\({}=\mathbf{A}\mathbf{A}^\T\), so the two distributions, and hence the two likelihoods, coincide: this is (12.39) with \(\sigma^2=0\). Unmixing with \((\mathbf{A}\mathbf{R})^{-1}\) gives \(\mathbf{R}^\T\mathbf{z}\sim\N(\mathbf{0},\mathbf{R}^\T\mathbf{R})\)\({}=\N(\mathbf{0},\mathbf{I})\), whose density factorises: its components are independent, so every rotation produces "independent components". Equivalently, after whitening (Result 12.4) any further rotation keeps the coordinates uncorrelated, and for jointly Gaussian variables uncorrelated means independent. For non-Gaussian sources the higher-order structure breaks the rotational symmetry. ∎

The original ICA algorithm (Bell and Sejnowski, 1995) maximised an information-theoretic objective (infomax). A probabilistic formulation suggests generalisations, for example independent factor analysis (Attias, 1999a), with noisy observations, possibly different numbers of latent and observed variables, and source densities modelled by Gaussian mixtures, fitted by EM with a variational approximation for the posterior over the sources. The ICA literature is large (Jutten and Hérault, 1991; Comon et al., 1991; Amari et al., 1996; Hyvärinen and Oja, 1997, 2000; Stone, 2004, and others listed in the book).

12.4.2 Autoassociative neural networks

Neural networks (Chapter 5) can also reduce dimension without supervision. An autoassociative network has \(D\) inputs, \(D\) outputs and \(M\lt D\) hidden units (Figure 12.18 of the book), and is trained to reproduce its input, by minimising

\[ E(\mathbf{w})=\frac12\sum_{n=1}^N\lVert\mathbf{y}(\mathbf{x}_n,\mathbf{w})-\mathbf{x}_n\rVert^2 . \tag{12.91} \]

Because \(M\lt D\), perfect reconstruction is impossible in general. The book states that with linear hidden units the error has a unique global minimum at which the network projects onto the principal subspace (Bourlard and Kamp, 1988; Baldi and Hornik, 1989), and that nonlinear hidden units do not help. Both claims follow from Section 12.1.2.

Proof

Whatever \(\mathbf{g}\) is, every output lies in the set \(\mathcal{A}\)\({}=\{\mathbf{W}_2\mathbf{h}+\mathbf{b}_2:\mathbf{h}\in\R^M\}\), an affine subspace of dimension at most \(M\). Hence \(\lVert\mathbf{y}(\mathbf{x}_n)-\mathbf{x}_n\rVert^2\)\({}\ge\operatorname{dist}(\mathbf{x}_n,\mathcal{A})^2\), and \(\frac2NE\) is at least the mean squared distance from the data to an affine subspace of dimension at most \(M\) (enlarge it to dimension exactly \(M\) if necessary; distances can only decrease). By Result 12.3 this is at least \(\sum_{i\gt M}\lambda_i\). With linear hidden units, \(\mathbf{W}_1=\mathbf{U}_M^\T\), \(\mathbf{b}_1\)\({}=-\mathbf{U}_M^\T\bar{\mathbf{x}}\), \(\mathbf{W}_2=\mathbf{U}_M\), \(\mathbf{b}_2=\bar{\mathbf{x}}\) give exactly the PCA reconstruction and attain the bound; conversely, equality requires every output to be the orthogonal projection of its input onto the principal subspace (Result 12.3), which for an affine network means \(\mathbf{W}_2\mathbf{W}_1\)\({}=\mathbf{U}_M\mathbf{U}_M^\T\). That matrix has rank \(M\), so \(\mathbf{W}_1\) (an \(M\times D\) matrix) has rank \(M\) and its row space contains, hence equals, \(\operatorname{span}\{\mathbf{u}_1,\dots,\mathbf{u}_M\}\): \(\mathbf{W}_1\)\({}=\mathbf{G}\mathbf{U}_M^\T\) with \(\mathbf{G}\) invertible, and then \(\mathbf{W}_2\mathbf{G}\)\({}=\mathbf{U}_M\). (The global minimum is unique as a map, not as weights; Baldi and Hornik (1989) show that the other critical points are saddles.) A sigmoid is nearly linear near the origin, so small input weights followed by compensating large output weights approach the linear solution. In our numerical check a tanh hidden layer reached an error 0.3% above the PCA bound, never below it. ∎

So two-layer networks bring nothing over standard PCA, which (via the singular value decomposition) is guaranteed to give the correct answer in finite time, with ordered eigenvalues and orthonormal eigenvectors. The situation changes with more layers (Figures 12.19 and 12.20 of the book): a network with a nonlinear first hidden layer, an \(M\)-unit second hidden layer (linear, say), a nonlinear third hidden layer and linear outputs computes two successive mappings, \(F_1:\R^D\to\R^M\) (a general nonlinear projection onto a subspace \(\mathcal{S}\)) and \(F_2:\R^M\to\R^D\) (a general embedding of \(\mathcal{S}\) back into the data space, which may be curved). Such a network performs nonlinear PCA and includes linear PCA as a special case, but the error (12.91) is no longer quadratic in the weights: training is a nonlinear optimisation with possible local minima, and \(M\) must be chosen in advance.

12.4.3 Modelling nonlinear manifolds

Many natural data lie near low-dimensional, noisy, curved manifolds, and capturing this explicitly can give better density models than general-purpose methods. The book surveys several approaches.

Mixtures of linear models

A curved manifold can be approximated piecewise linearly: cluster the data (for example with K-means) and apply PCA within each cluster. It is better to assign points by reconstruction error, so that one cost is optimised throughout (Kambhatla and Leen, 1997; Hinton, Dayan and Revow, 1997), but these methods still lack an overall density model. With PPCA one obtains a fully probabilistic mixture of probabilistic PCA models (Tipping and Bishop, 1999a),

\[ p(\mathbf{x})=\sum_{k=1}^K\pi_k\,\N\big(\mathbf{x}\given\bm{\mu}_k,\mathbf{W}_k\mathbf{W}_k^\T+\sigma_k^2\mathbf{I}\big), \]

which has a discrete latent variable (the component) and a continuous one (the position within the component's subspace), and is fitted by EM.

Why the M step is "PPCA on the weighted covariance"

With only the component label treated as latent (the continuous \(\mathbf{z}\) integrated out analytically), the expected complete-data log likelihood is \(\sum_n\sum_kR_{nk}\{\ln\pi_k+\ln\N(\mathbf{x}_n\given\bm{\mu}_k,\mathbf{C}_k)\}\). The \(\pi_k\) update is that of the Gaussian mixture (9.22). For component \(k\), the gradient in \(\bm{\mu}_k\) is \(\mathbf{C}_k^{-1}\sum_nR_{nk}(\mathbf{x}_n-\bm{\mu}_k)\), zero at the weighted mean whatever \(\mathbf{C}_k\) is. Substituting it, the trace trick gives

\[ \sum_nR_{nk}\ln\N(\mathbf{x}_n\given\bm{\mu}_k,\mathbf{C}_k)=-\frac{N_k}{2}\Big\{D\ln(2\pi)+\ln\lvert\mathbf{C}_k\rvert+\Tr\big(\mathbf{C}_k^{-1}\mathbf{S}_k\big)\Big\}, \]

which is (12.44) with \(N_k\) and \(\mathbf{S}_k\) in place of \(N\) and \(\mathbf{S}\). Result 12.8 therefore gives its exact maximiser. (Tipping and Bishop, 1999a, also describe a variant that keeps \(\mathbf{z}\) as a latent variable and avoids the eigendecomposition.) ∎

A variational Bayesian treatment infers both the number of components and their dimensions (Bishop and Winn, 2000). Parameters such as \(\mathbf{W}\) or the noise variances can be tied across components, and replacing isotropic by diagonal noise gives a mixture of factor analysers (Ghahramani and Hinton, 1996a; Ghahramani and Beal, 2000). A hierarchical mixture of PPCA models supports interactive visualisation (Bishop and Tipping, 1998). The graphical models are drawn in Exercise 12.23.

Principal curves

Instead of many linear pieces, one can fit a single nonlinear model. PCA finds a linear subspace close to the data in least squares; principal curves (Hastie and Stuetzle, 1989) extend this to one-dimensional curves. A curve is a vector function \(\mathbf{f}(\lambda)\) of a scalar \(\lambda\), naturally parametrised by arc length. For a point \(\widehat{\mathbf{x}}\) let \(\lambda\)\({}=g_{\mathbf{f}}(\widehat{\mathbf{x}})\) be the parameter of the closest point of the curve. For a continuous density \(p(\mathbf{x})\), \(\mathbf{f}\) is a principal curve if every point of the curve is the mean of the points that project onto it,

\[ \E\big[\mathbf{x}\given g_{\mathbf{f}}(\mathbf{x})=\lambda\big]=\mathbf{f}(\lambda). \tag{12.92} \]

A density can have many principal curves, and with finite data one restricts attention to smooth curves. Hastie and Stuetzle alternate, in a way reminiscent of EM for PCA, a projection step (each point gets the \(\lambda\) of its closest point on the curve) and a re-estimation step (each point of the curve becomes a weighted average of the data points that project nearby, the nearest ones weighted most), starting from the first principal component. If the curve is constrained to be a straight line, the procedure converges to the first principal component and is the power method: with \(\mathbf{f}(\lambda)\)\({}=\bar{\mathbf{x}}+\lambda\mathbf{u}\), \(\lVert\mathbf{u}\rVert=1\), the projection step gives \(\lambda_n\)\({}=\mathbf{u}^\T(\mathbf{x}_n-\bar{\mathbf{x}})\); the best line through \(\bar{\mathbf{x}}\) with these \(\lambda_n\) has direction \(\sum_n\lambda_n(\mathbf{x}_n-\bar{\mathbf{x}})/\sum_n\lambda_n^2\)\({}=\mathbf{S}\mathbf{u}/(\mathbf{u}^\T\mathbf{S}\mathbf{u})\), and re-parametrising by arc length normalises it: \(\mathbf{u}_{\rm new}\)\({}=\mathbf{S}\mathbf{u}/\lVert\mathbf{S}\mathbf{u}\rVert\). (This is Result 12.12 with \(M=1\).) Principal surfaces generalise the idea to more dimensions but are little used, because smoothing data in several dimensions is hard.

Multidimensional scaling, LLE and Isomap

Multidimensional scaling (MDS; Torgerson, 1952; Cox and Cox, 2000) seeks a low-dimensional configuration of points whose pairwise distances match given distances as closely as possible; with Euclidean distances the classical version is equivalent to PCA. Using general dissimilarities gives nonmetric MDS.

Proof

Distances do not change under translation, so with \(\widetilde{\mathbf{x}}_n\)\({}=\mathbf{x}_n-\bar{\mathbf{x}}\), \(\Delta_{nm}\)\({}=\lVert\widetilde{\mathbf{x}}_n\rVert^2+\lVert\widetilde{\mathbf{x}}_m\rVert^2-2\widetilde{\mathbf{x}}_n^\T\widetilde{\mathbf{x}}_m\). With \(\mathbf{G}\)\({}=\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^\T\) and \(\mathbf{g}\) its diagonal, \(\bm{\Delta}\)\({}=\mathbf{g}\mathbf{1}^\T+\mathbf{1}\mathbf{g}^\T-2\mathbf{G}\). Since \(\mathbf{H}\mathbf{1}=\mathbf{0}\), the first two terms are annihilated by the double centring, and \(\mathbf{H}\widetilde{\mathbf{X}}\)\({}=\widetilde{\mathbf{X}}\) (the columns of \(\widetilde{\mathbf{X}}\) already sum to zero) gives \(\mathbf{H}\mathbf{G}\mathbf{H}\)\({}=\mathbf{G}\). Hence \(\mathbf{B}\)\({}=-\frac12(-2\mathbf{G})\)\({}=\mathbf{G}\). By Result 12.5, the eigenvectors of \(\mathbf{G}\)\({}=N\cdot N^{-1}\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^\T\) are the \(\mathbf{v}_i\) of (12.28) with \(\kappa_i=N\lambda_i\), and the PCA scores are \(\widetilde{\mathbf{X}}\mathbf{u}_i\)\({}=\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^\T\mathbf{v}_i/(N\lambda_i)^{1/2}\)\({}=N\lambda_i\mathbf{v}_i/(N\lambda_i)^{1/2}\)\({}=\sqrt{\kappa_i}\,\mathbf{v}_i\), using (12.30). So MDS needs only the distances, which is what lets it be applied to other dissimilarities. ∎

Locally linear embedding (LLE; Roweis and Saul, 2000) first finds, for each point, the weights that best reconstruct it from its \(K\) nearest neighbours, \(\min\lVert\mathbf{x}_n-\sum_mW_{nm}\mathbf{x}_m\rVert^2\) subject to \(\sum_mW_{nm}=1\) (with \(W_{nm}=0\) for non-neighbours). Because the weights sum to one, the residual equals \(\lVert\sum_mW_{nm}(\mathbf{x}_n-\mathbf{x}_m)\rVert^2\)\({}=\mathbf{w}_n^\T\mathbf{G}_n\mathbf{w}_n\) with the local Gram matrix \((\mathbf{G}_n)_{jk}\)\({}=(\mathbf{x}_n-\mathbf{x}_j)^\T(\mathbf{x}_n-\mathbf{x}_k)\), and a Lagrange multiplier for the constraint gives \(\mathbf{w}_n\)\({}=\mathbf{G}_n^{-1}\mathbf{1}/(\mathbf{1}^\T\mathbf{G}_n^{-1}\mathbf{1})\) (regularised if \(\mathbf{G}_n\) is singular). The weights are unchanged by translating, rotating or rescaling a point together with its neighbours, so they describe the local geometry. LLE then finds low-dimensional points \(\mathbf{y}_n\) that are reconstructed by the same weights, minimising \(\sum_n\lVert\mathbf{y}_n-\sum_mW_{nm}\mathbf{y}_m\rVert^2\)\({}=\Tr\{\mathbf{Y}^\T(\mathbf{I}-\mathbf{W})^\T(\mathbf{I}-\mathbf{W})\mathbf{Y}\}\) subject to \(\sum_n\mathbf{y}_n=\mathbf{0}\) and \(N^{-1}\mathbf{Y}^\T\mathbf{Y}\)\({}=\mathbf{I}\). As in Result 12.2 (with minimum instead of maximum) the solution uses the eigenvectors of the smallest eigenvalues of \((\mathbf{I}-\mathbf{W})^\T(\mathbf{I}-\mathbf{W})\), discarding the constant vector (eigenvalue 0, because \(\mathbf{W}\mathbf{1}=\mathbf{1}\)). Although the embedding is nonlinear, both steps are solved exactly, with no local minima.

Isomap (Tenenbaum, de Silva and Langford, 2000) applies classical MDS to geodesic distances measured along the manifold (for two points on a circle, the arc length rather than the chord). It (i) links each point to its \(K\) nearest neighbours, or to all points within a radius \(\epsilon\), with edges weighted by Euclidean distance; (ii) approximates geodesic distances by shortest paths in this graph (Dijkstra's or Floyd's algorithm); (iii) runs classical MDS (Result 12.18) on the matrix of graph distances.

Latent trait models and density networks

With continuous latent variables and discrete observed variables one obtains latent trait models (Bartholomew, 1987); for binary data, for example, \(p(x_d=1\given\mathbf{z})\)\({}=\sigma(\mathbf{w}_d^\T\mathbf{z}+b_d)\). The integral over \(\mathbf{z}\) can no longer be done analytically, even with a linear dependence on \(\mathbf{z}\). Tipping (1999) used the variational bound of Section 10.6 to visualise binary data with a two-dimensional latent space, as PCA visualises continuous data. This is the dual of Bayesian logistic regression (Section 4.5): there a single uncertain weight vector multiplies \(N\) known feature vectors; here each of the \(N\) data points has its own uncertain latent vector, which multiplies a shared set of weights, so the roles of the two kinds of vector are exchanged. Exponential-family generalisations are described by Collins, Dasgupta and Schapire (2002).

A density network (MacKay, 1995; MacKay and Gibbs, 1999) takes \(\mathbf{z}\sim\N(\mathbf{0},\mathbf{I})\) and maps it through a multilayer neural network (Chapter 5), which with enough hidden units can approximate any nonlinear map, so very general densities result. The price is that the marginal likelihood is intractable. It is approximated by Monte Carlo (Chapter 11), \(p(\mathbf{x})\)\({}\approx\frac1L\sum_{l=1}^Lp(\mathbf{x}\given\mathbf{z}^{(l)})\) with \(\mathbf{z}^{(l)}\) drawn from the prior, but many samples may be needed for an accurate approximation, which is costly.

The generative topographic mapping

Restricting the nonlinear map and choosing a convenient latent distribution gives a model that is both nonlinear and cheap to train: the generative topographic mapping (GTM; Bishop, Svensén and Williams, 1996, 1997a, 1998b).

Integrating over the latent space is a sum over the grid, and because \(\mathbf{y}\) is linear in \(\mathbf{W}\), the likelihood is available in closed form and can be maximised by EM. The curse of dimensionality that limits basis-function regression (Section 1.4) does not bite, because the basis functions live on the two-dimensional latent space whatever the dimension of the data.

Derivation of the M step

The E step is Bayes' theorem for a mixture with equal weights \(1/K\) and a common covariance, so the normalising constants cancel. The expected complete-data log likelihood is

\[ \mathcal{Q}(\mathbf{W},\beta)=\sum_{n=1}^N\sum_{k=1}^KR_{kn}\Big\{\frac D2\ln\beta-\frac\beta2\lVert\mathbf{x}_n-\mathbf{W}\bm{\phi}_k\rVert^2\Big\}+\const,\qquad\bm{\phi}_k=\bm{\phi}(\mathbf{z}_k). \]

Its gradient in \(\mathbf{W}\) is \(\beta\sum_n\sum_kR_{kn}(\mathbf{x}_n-\mathbf{W}\bm{\phi}_k)\bm{\phi}_k^\T\). Setting it to zero, \(\mathbf{W}\sum_kG_{kk}\bm{\phi}_k\bm{\phi}_k^\T\)\({}=\sum_k\sum_nR_{kn}\mathbf{x}_n\bm{\phi}_k^\T\), i.e. \(\mathbf{W}\bm{\Phi}^\T\mathbf{G}\bm{\Phi}\)\({}=\mathbf{X}^\T\mathbf{R}^\T\bm{\Phi}\); transposing gives the stated linear system. It is a weighted least-squares problem (concave in \(\mathbf{W}\)), and a Gaussian prior \(-\frac\lambda2\lVert\mathbf{W}\rVert_F^2\) adds \(\lambda\mathbf{W}\) to the normal equations, i.e. \(\frac\lambda\beta\mathbf{I}\) after dividing by \(\beta\). The derivative in \(\beta\), \(\frac{ND}{2\beta}-\frac12\sum_{n,k}R_{kn}\lVert\mathbf{x}_n-\mathbf{W}\bm{\phi}_k\rVert^2\) (using \(\sum_kR_{kn}=1\)), vanishes at the stated \(\beta_{\rm new}\). (Our numerical check confirmed that these updates never decrease the log likelihood.) ∎

GTM and the self-organising map. GTM can be seen as a probabilistic version of the self-organising map (SOM; Kohonen, 1982, 1995), which also represents a two-dimensional manifold by a grid of prototype vectors \(\mathbf{m}_k\). The SOM resembles K-means: each data point is assigned to its nearest prototype \(c\)\({}=\argmin_k\lVert\mathbf{x}_n-\mathbf{m}_k\rVert\), and the prototypes are moved towards it, \(\mathbf{m}_k\leftarrow\mathbf{m}_k+\eta\,h_{ck}(\mathbf{x}_n-\mathbf{m}_k)\), with a neighbourhood function \(h_{ck}\) on the grid that makes neighbouring prototypes move together, so the initially random prototypes "self-organise" into a smooth sheet. But the SOM does not optimise any well-defined cost function (Erwin, Obermayer and Schulten, 1992), which makes its parameters hard to set and its convergence hard to assess, and organisation is not guaranteed. GTM instead maximises a likelihood that EM increases monotonically and defines a density in data space. The smoothness that SOM imposes through \(h_{ck}\) is built into GTM through the smooth map \(\mathbf{W}\bm{\phi}(\mathbf{z})\), and the hard winner-take-all assignment is replaced by soft responsibilities, which become hard as \(\beta\to\infty\) (as K-means is the hard limit of EM for Gaussian mixtures, Section 9.3.2). The probabilistic foundation makes extensions straightforward (Bishop, Svensén and Williams, 1998a): Bayesian treatments, missing values, discrete variables, Gaussian-process manifolds (Section 6.4), hierarchical GTM (Tiño and Nabney, 2002). And because the manifold is a continuous surface, not just a set of prototypes, one can compute its local stretching, the magnification factor \(\det(\mathbf{J}^\T\mathbf{J})^{1/2}\) with \(\mathbf{J}\)\({}=\mathbf{W}\,\partial\bm{\phi}/\partial\mathbf{z}^\T\) (Bishop, Svensén and Williams, 1997b), and its directional curvatures (Tiño, Nabney and Sun, 2001), and display them with the projected data.

Exercises

All 29 exercises of Chapter 12 (Bishop, pp. 599–603), with full solutions (hidden). Exercises marked www have official solutions (Svensén and Bishop); the solutions below are self-contained and were checked independently, and two slips in the published solution manuals are pointed out where they occur. Where the main text already contains the complete derivation, the solution summarises it and links to it. Difficulty stars are the book's.

1. Exercise 12.1 (★★, www). Prove by induction that the \(M\)-dimensional projection with maximal variance is given by the eigenvectors of the \(M\) largest eigenvalues of \(\mathbf{S}\): assuming the result for \(M\), maximise the variance along a new unit direction \(\mathbf{u}_{M+1}\) orthogonal to \(\mathbf{u}_1,\dots,\mathbf{u}_M\), using Lagrange multipliers for all the constraints.

   Solution

   The case \(M=1\) is Result 12.1. Assume \(\mathbf{u}_1,\dots,\mathbf{u}_M\) are the eigenvectors of the \(M\) largest eigenvalues. Maximise \(\mathbf{u}^\T\mathbf{S}\mathbf{u}\) subject to \(\mathbf{u}^\T\mathbf{u}=1\) and \(\mathbf{u}^\T\mathbf{u}_j=0\) (\(j\le M\)) with the Lagrangian \(\mathbf{u}^\T\mathbf{S}\mathbf{u}+\lambda(1-\mathbf{u}^\T\mathbf{u})+\sum_{j\le M}\eta_j\mathbf{u}^\T\mathbf{u}_j\). Stationarity: \(2\mathbf{S}\mathbf{u}-2\lambda\mathbf{u}+\sum_j\eta_j\mathbf{u}_j\)\({}=\mathbf{0}\). Left-multiplying by \(\mathbf{u}_m^\T\) and using \(\mathbf{u}_m^\T\mathbf{S}\mathbf{u}\)\({}=\lambda_m\mathbf{u}_m^\T\mathbf{u}\)\({}=0\), \(\mathbf{u}_m^\T\mathbf{u}=0\), \(\mathbf{u}_m^\T\mathbf{u}_j\)\({}=\delta_{mj}\) gives \(\eta_m=0\). Hence \(\mathbf{S}\mathbf{u}\)\({}=\lambda\mathbf{u}\), the variance along \(\mathbf{u}\) is \(\mathbf{u}^\T\mathbf{S}\mathbf{u}\)\({}=\lambda\), and \(\mathbf{u}\) must be an eigenvector orthogonal to the first \(M\), so \(\lambda\in\{\lambda_{M+1},\dots,\lambda_D\}\); the maximum is \(\lambda_{M+1}\) at \(\mathbf{u}_{M+1}\). The total variance of the \((M+1)\)-dimensional projection is then \(\lambda_1+\dots+\lambda_{M+1}\), and Result 12.2(b) shows that no \((M+1)\)-dimensional projection does better, so the greedy construction is also globally optimal. Details: Result 12.2.

2. Exercise 12.2 (★★). Show that the minimum of the distortion (12.15) subject to (12.7) is obtained when the \(\mathbf{u}_i\) are eigenvectors of \(\mathbf{S}\), using the matrix of multipliers \(\mathbf{H}\) in (12.93); show that \(\mathbf{H}\) may be taken symmetric and that the general solution of \(\mathbf{S}\widehat{\mathbf{U}}\)\({}=\widehat{\mathbf{U}}\mathbf{H}\) gives the same \(J\) as the eigenvector solution.

   Solution

   Only the symmetric part of \(\mathbf{H}\) enters (12.93), because \(\mathbf{I}-\widehat{\mathbf{U}}^\T\widehat{\mathbf{U}}\) is symmetric. Differentiating, \(2\mathbf{S}\widehat{\mathbf{U}}-2\widehat{\mathbf{U}}\mathbf{H}\)\({}=\mathbf{0}\). With \(\mathbf{H}\)\({}=\bm{\Psi}\widehat{\bm{\Lambda}}\bm{\Psi}^\T\), the columns of \(\widehat{\mathbf{U}}\bm{\Psi}\) are orthonormal eigenvectors of \(\mathbf{S}\), so \(\widehat{\mathbf{U}}\) is any orthonormal basis of a span of \(D-M\) eigenvectors, and \(J\)\({}=\Tr(\widehat{\mathbf{U}}^\T\mathbf{S}\widehat{\mathbf{U}})\)\({}=\Tr\mathbf{H}\)\({}=\Tr\widehat{\bm{\Lambda}}\), the same value as for the eigenvector solution \(\widehat{\mathbf{U}}\bm{\Psi}\). The minimum (which exists and satisfies the multiplier rule) is the sum of the \(D-M\) smallest eigenvalues. Full details: Result 12.3.

3. Exercise 12.3 (★). Verify that the vectors (12.30) have unit length when the \(\mathbf{v}_i\) do.

   Solution

   \(\lVert\mathbf{X}^\T\mathbf{v}_i\rVert^2\)\({}=\mathbf{v}_i^\T\mathbf{X}\mathbf{X}^\T\mathbf{v}_i\)\({}=\mathbf{v}_i^\T(N\lambda_i\mathbf{v}_i)\)\({}=N\lambda_i\) by (12.28), so \(\lVert\mathbf{u}_i\rVert^2\)\({}=N\lambda_i/(N\lambda_i)\)\({}=1\). The same computation with \(\mathbf{v}_i^\T\mathbf{X}\mathbf{X}^\T\mathbf{v}_j\)\({}=N\lambda_j\mathbf{v}_i^\T\mathbf{v}_j\) shows that the \(\mathbf{u}_i\) are mutually orthogonal (Result 12.5).

4. Exercise 12.4 (★, www). Replace the prior (12.31) by \(\N(\mathbf{z}\given\mathbf{m},\bm{\Sigma})\). Show that, after redefining the parameters, the marginal \(p(\mathbf{x})\) is the same as for the standard model, for any \(\mathbf{m}\) and \(\bm{\Sigma}\).

   Solution

   With \(\mathbf{x}\)\({}=\mathbf{W}\mathbf{z}+\bm{\mu}+\bm{\epsilon}\) and \(\mathbf{z}\sim\N(\mathbf{m},\bm{\Sigma})\) independent of \(\bm{\epsilon}\), \(\mathbf{x}\) is Gaussian with mean \(\mathbf{W}\mathbf{m}+\bm{\mu}\) and covariance \(\mathbf{W}\bm{\Sigma}\mathbf{W}^\T+\sigma^2\mathbf{I}\). Factor \(\bm{\Sigma}\)\({}=\mathbf{L}\mathbf{L}^\T\) (Cholesky or symmetric square root) and set \(\widetilde{\mathbf{W}}\)\({}=\mathbf{W}\mathbf{L}\), \(\widetilde{\bm{\mu}}\)\({}=\bm{\mu}+\mathbf{W}\mathbf{m}\). Then \(p(\mathbf{x})\)\({}=\N(\mathbf{x}\given\widetilde{\bm{\mu}},\widetilde{\mathbf{W}}\widetilde{\mathbf{W}}^\T+\sigma^2\mathbf{I})\), the marginal of the standard model with parameters \(\widetilde{\mathbf{W}},\widetilde{\bm{\mu}},\sigma^2\). Equivalently, \(\mathbf{z}\)\({}=\mathbf{m}+\mathbf{L}\mathbf{z}'\) with \(\mathbf{z}'\sim\N(\mathbf{0},\mathbf{I})\) and \(\mathbf{W}\mathbf{z}+\bm{\mu}\)\({}=\widetilde{\mathbf{W}}\mathbf{z}'+\widetilde{\bm{\mu}}\). Conversely, the standard model is the case \(\mathbf{m}=\mathbf{0}\), \(\bm{\Sigma}=\mathbf{I}\), so the two families of distributions coincide and the extra parameters are redundant. (The solution manuals write the covariance as \(\sigma^2\mathbf{I}+\mathbf{W}^\T\bm{\Sigma}^{-1}\mathbf{W}\) or \(\mathbf{W}\bm{\Sigma}^{-1}\mathbf{W}^\T\); the correct form is \(\mathbf{W}\bm{\Sigma}\mathbf{W}^\T\).)

5. Exercise 12.5 (★★). If \(\mathbf{x}\sim\N(\bm{\mu},\bm{\Sigma})\) in \(\R^D\) and \(\mathbf{y}\)\({}=\mathbf{A}\mathbf{x}+\mathbf{b}\) with \(\mathbf{A}\) an \(M\times D\) matrix, show that \(\mathbf{y}\) is Gaussian, find its mean and covariance, and discuss \(M\lt D\), \(M=D\) and \(M\gt D\).

   Solution

   The moment generating function is \(\E[e^{\mathbf{t}^\T\mathbf{y}}]\)\({}=e^{\mathbf{t}^\T\mathbf{b}}\E[e^{(\mathbf{A}^\T\mathbf{t})^\T\mathbf{x}}]\)\({}=\exp\{\mathbf{t}^\T(\mathbf{A}\bm{\mu}+\mathbf{b})+\tfrac12\mathbf{t}^\T\mathbf{A}\bm{\Sigma}\mathbf{A}^\T\mathbf{t}\}\), using the Gaussian MGF \(\E[e^{\mathbf{s}^\T\mathbf{x}}]\)\({}=\exp(\mathbf{s}^\T\bm{\mu}+\tfrac12\mathbf{s}^\T\bm{\Sigma}\mathbf{s})\). This is the MGF of \(\N(\mathbf{A}\bm{\mu}+\mathbf{b},\mathbf{A}\bm{\Sigma}\mathbf{A}^\T)\), and the MGF determines the distribution. Assume \(\bm{\Sigma}\) positive definite. \(M\lt D\): \(\mathbf{A}\bm{\Sigma}\mathbf{A}^\T\) is positive definite iff \(\mathbf{A}\) has full row rank \(M\) (then \(\mathbf{v}^\T\mathbf{A}\bm{\Sigma}\mathbf{A}^\T\mathbf{v}\)\({}=0\) forces \(\mathbf{A}^\T\mathbf{v}=\mathbf{0}\), i.e. \(\mathbf{v}=\mathbf{0}\)), and \(\mathbf{y}\) has a proper density on \(\R^M\) (marginals of Gaussians are the special case where \(\mathbf{A}\) selects coordinates). \(M=D\): if \(\mathbf{A}\) is invertible, \(p(\mathbf{y})\)\({}=p_{\mathbf{x}}(\mathbf{A}^{-1}(\mathbf{y}-\mathbf{b}))/\lvert\det\mathbf{A}\rvert\), a nondegenerate Gaussian. \(M\gt D\): \(\operatorname{rank}(\mathbf{A}\bm{\Sigma}\mathbf{A}^\T)\)\({}\le D\)\({}\lt M\), so the covariance is singular: \(\mathbf{y}\) lives on the affine subspace \(\mathbf{b}+\operatorname{col}\mathbf{A}\) of dimension at most \(D\) and has no density on \(\R^M\). This is PPCA with \(\sigma^2=0\) (the map \(\mathbf{z}\mapsto\mathbf{W}\mathbf{z}+\bm{\mu}\) from \(\R^M\) to \(\R^D\), \(D\gt M\)), and it is why the noise term is needed for a proper density.

6. Exercise 12.6 (★, www). Draw the directed graph of PPCA with the components of \(\mathbf{x}\) as separate nodes and verify that it has the independence structure of the naive Bayes model (Section 8.2.2).

   Solution

   Because \(\sigma^2\mathbf{I}\) is diagonal, \(p(\mathbf{x}\given\mathbf{z})\)\({}=\prod_d\N(x_d\given\widetilde{\mathbf{w}}_d^\T\mathbf{z}+\mu_d,\sigma^2)\): each \(x_d\) has the single parent \(\mathbf{z}\) (the parameters \(\widetilde{\mathbf{w}}_d,\mu_d,\sigma^2\) are omitted) and there are no edges between the \(x_d\). This is the naive Bayes graph of Figure 8.24, with the class variable replaced by the continuous \(\mathbf{z}\). Every path between \(x_i\) and \(x_j\) passes tail-to-tail through \(\mathbf{z}\), so by d-separation \(x_i\perp\!\!\!\perp x_j\given\mathbf{z}\), exactly the naive Bayes assumption. Marginally the \(x_d\) are dependent (the path is open when \(\mathbf{z}\) is unobserved), which is how the model represents correlations: \(\cov[x_i,x_j]\)\({}=\widetilde{\mathbf{w}}_i^\T\widetilde{\mathbf{w}}_j\) for \(i\ne j\).

7. Exercise 12.7 (★★). Using the results (2.270)–(2.271) for the mean and covariance of a general distribution, derive the marginal (12.35).

   Solution

   (2.270)–(2.271) are the laws of total expectation and total covariance: \(\E[\mathbf{x}]\)\({}=\E_{\mathbf{z}}[\E[\mathbf{x}\given\mathbf{z}]]\) and \(\cov[\mathbf{x}]\)\({}=\E_{\mathbf{z}}[\cov[\mathbf{x}\given\mathbf{z}]]+\cov_{\mathbf{z}}[\E[\mathbf{x}\given\mathbf{z}]]\). (The second follows by writing \(\mathbf{x}-\E\mathbf{x}\)\({}=(\mathbf{x}-\E[\mathbf{x}\given\mathbf{z}])+(\E[\mathbf{x}\given\mathbf{z}]-\E\mathbf{x})\); the cross term has zero expectation because the first bracket has zero conditional mean.) Here \(\E[\mathbf{x}\given\mathbf{z}]\)\({}=\mathbf{W}\mathbf{z}+\bm{\mu}\) and \(\cov[\mathbf{x}\given\mathbf{z}]\)\({}=\sigma^2\mathbf{I}\), so \(\E[\mathbf{x}]\)\({}=\mathbf{W}\E[\mathbf{z}]+\bm{\mu}\)\({}=\bm{\mu}\) and \(\cov[\mathbf{x}]\)\({}=\sigma^2\mathbf{I}+\mathbf{W}\cov[\mathbf{z}]\mathbf{W}^\T\)\({}=\sigma^2\mathbf{I}+\mathbf{W}\mathbf{W}^\T\), which are (12.37)–(12.38). The distribution is Gaussian because \(\mathbf{x}\) is an affine function of the Gaussian vector \((\mathbf{z},\bm{\epsilon})\) (Exercise 12.5), or by the linear-Gaussian marginal (2.115); a direct computation is in Result 12.6.

8. Exercise 12.8 (★★, www). Using (2.116), show that the posterior \(p(\mathbf{z}\given\mathbf{x})\) is (12.42).

   Solution

   Match \(p(\mathbf{z})\)\({}=\N(\mathbf{z}\given\mathbf{0},\mathbf{I})\) with \(\N(\mathbf{z}\given\mathbf{m},\bm{\Lambda}^{-1})\) (so \(\mathbf{m}=\mathbf{0}\), \(\bm{\Lambda}=\mathbf{I}\)) and \(p(\mathbf{x}\given\mathbf{z})\) with \(\N(\mathbf{x}\given\mathbf{A}\mathbf{z}+\mathbf{b},\mathbf{L}^{-1})\) (so \(\mathbf{A}=\mathbf{W}\), \(\mathbf{b}=\bm{\mu}\), \(\mathbf{L}=\sigma^{-2}\mathbf{I}\)). By (2.116)–(2.117) the posterior covariance is \((\bm{\Lambda}+\mathbf{A}^\T\mathbf{L}\mathbf{A})^{-1}\)\({}=(\mathbf{I}+\sigma^{-2}\mathbf{W}^\T\mathbf{W})^{-1}\)\({}=\sigma^2(\sigma^2\mathbf{I}+\mathbf{W}^\T\mathbf{W})^{-1}\)\({}=\sigma^2\mathbf{M}^{-1}\), and the mean is \(\sigma^2\mathbf{M}^{-1}\{\mathbf{W}^\T\sigma^{-2}(\mathbf{x}-\bm{\mu})+\mathbf{0}\}\)\({}=\mathbf{M}^{-1}\mathbf{W}^\T(\mathbf{x}-\bm{\mu})\). This is (12.42); the first printing of the book has the covariance wrong (see the erratum noted at (12.42)). A self-contained derivation is in Result 12.6.

9. Exercise 12.9 (★). Show that maximising (12.43) over \(\bm{\mu}\) gives \(\bm{\mu}_{\rm ML}=\bar{\mathbf{x}}\).

   Solution

   \(\partial/\partial\bm{\mu}\) of (12.43) is \(\mathbf{C}^{-1}\sum_n(\mathbf{x}_n-\bm{\mu})\)\({}=N\mathbf{C}^{-1}(\bar{\mathbf{x}}-\bm{\mu})\), and \(\mathbf{C}^{-1}\) is invertible, so the only stationary point is \(\bm{\mu}=\bar{\mathbf{x}}\) (Result 12.7).

10. Exercise 12.10 (★★). Using second derivatives, show that \(\bm{\mu}_{\rm ML}=\bar{\mathbf{x}}\) is the unique maximum.

    Solution

    The Hessian with respect to \(\bm{\mu}\) is \(-N\mathbf{C}^{-1}\), constant and negative definite (\(\mathbf{C}\succ0\) because \(\mathbf{v}^\T\mathbf{C}\mathbf{v}\)\({}=\lVert\mathbf{W}^\T\mathbf{v}\rVert^2+\sigma^2\lVert\mathbf{v}\rVert^2\)). So (12.43) is a strictly concave quadratic in \(\bm{\mu}\); exactly,

    \[ -\frac12\sum_n(\mathbf{x}_n-\bm{\mu})^\T\mathbf{C}^{-1}(\mathbf{x}_n-\bm{\mu})=-\frac N2(\bm{\mu}-\bar{\mathbf{x}})^\T\mathbf{C}^{-1}(\bm{\mu}-\bar{\mathbf{x}})+\text{(terms free of }\bm{\mu}), \]

    which is uniquely maximised at \(\bar{\mathbf{x}}\) (Result 12.7).

11. Exercise 12.11 (★★, www). Show that as \(\sigma^2\to0\) the posterior mean becomes an orthogonal projection onto the principal subspace, as in conventional PCA.

    Solution

    As \(\sigma^2\to0\), \(\mathbf{M}\)\({}\to\mathbf{W}^\T\mathbf{W}\) and (12.48) tends to (12.50). With \(\mathbf{W}\)\({}=\mathbf{U}_M(\mathbf{L}_M-\sigma^2\mathbf{I})^{1/2}\) (taking \(\mathbf{R}=\mathbf{I}\)) and \(\sigma^2\to0\), \((\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T(\mathbf{x}-\bar{\mathbf{x}})\)\({}=\mathbf{L}_M^{-1/2}\mathbf{U}_M^\T(\mathbf{x}-\bar{\mathbf{x}})\): the (whitened) PCA scores (12.24). Mapped back into data space, \(\mathbf{W}\E[\mathbf{z}\given\mathbf{x}]\)\({}\to\mathbf{U}_M\mathbf{U}_M^\T(\mathbf{x}-\bar{\mathbf{x}})\), the orthogonal projection onto the principal subspace. For general \(\mathbf{W}\), \(\mathbf{W}(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T\) is the orthogonal projection onto \(\operatorname{col}\mathbf{W}\) (Result 12.9).

12. Exercise 12.12 (★★). For \(\sigma^2\gt0\), show that the posterior mean is shifted towards the origin relative to the orthogonal projection.

    Solution

    Diagonalise \(\mathbf{W}^\T\mathbf{W}\)\({}=\mathbf{Q}\bm{\Delta}\mathbf{Q}^\T\). In the coordinates \(\mathbf{Q}^\T\mathbf{z}\), the posterior mean is \((\bm{\Delta}+\sigma^2\mathbf{I})^{-1}\mathbf{c}\) and the orthogonal-projection coordinates are \(\bm{\Delta}^{-1}\mathbf{c}\), with \(\mathbf{c}\)\({}=\mathbf{Q}^\T\mathbf{W}^\T(\mathbf{x}-\bm{\mu})\). Each component is multiplied by \(\delta_k/(\delta_k+\sigma^2)\in(0,1)\), so the posterior mean is strictly shorter (closer to the origin) unless it is zero. At the maximum likelihood solution the factors are \((\lambda_i-\sigma^2)/\lambda_i\) along the principal directions (Result 12.9).

13. Exercise 12.13 (★★). Show that the optimal least-squares reconstruction of a data point from its PPCA posterior mean is (12.94), \(\widetilde{\mathbf{x}}\)\({}=\mathbf{W}_{\rm ML}(\mathbf{W}_{\rm ML}^\T\mathbf{W}_{\rm ML})^{-1}\mathbf{M}\,\E[\mathbf{z}\given\mathbf{x}]\).

    Solution

    The reconstruction must lie in the principal subspace \(\bar{\mathbf{x}}+\operatorname{col}\mathbf{W}\), and the point of an affine subspace closest to \(\mathbf{x}\) is the orthogonal projection \(\bar{\mathbf{x}}+\mathbf{W}(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T(\mathbf{x}-\bar{\mathbf{x}})\). By (12.48), \(\mathbf{W}^\T(\mathbf{x}-\bar{\mathbf{x}})\)\({}=\mathbf{M}\,\E[\mathbf{z}\given\mathbf{x}]\), which gives (12.94) for centred data (add \(\bar{\mathbf{x}}\) otherwise). The naive reconstruction \(\mathbf{W}\E[\mathbf{z}\given\mathbf{x}]+\bar{\mathbf{x}}\) is shrunk towards the mean (Exercise 12.12), and the factor \((\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{M}\) undoes the shrinkage (Result 12.9).

14. Exercise 12.14 (★). Verify that (12.51) gives the count of a general covariance when \(M=D-1\) and of an isotropic covariance when \(M=0\).

    Solution

    \(M=D-1\): \(D(D-1)+1-\frac{(D-1)(D-2)}{2}\)\({}=\frac{2D^2-2D+2-(D^2-3D+2)}{2}\)\({}=\frac{D^2+D}{2}\)\({}=\frac{D(D+1)}{2}\), the number of free entries of a symmetric \(D\times D\) matrix. \(M=0\): \(0+1-0=1\), the single parameter \(\sigma^2\) of \(\sigma^2\mathbf{I}\). Why (12.51) is a genuine dimension count is shown in Result 12.10.

15. Exercise 12.15 (★★, www). Derive the M-step equations (12.56) and (12.57) by maximising (12.53).

    Solution

    The two derivatives are

    \[ \begin{aligned}\frac{\partial\E[\ln p]}{\partial\mathbf{W}}&=\frac{1}{\sigma^2}\Big\{\sum_n(\mathbf{x}_n-\bar{\mathbf{x}})\E[\mathbf{z}_n]^\T-\mathbf{W}\sum_n\E[\mathbf{z}_n\mathbf{z}_n^\T]\Big\},\\ \frac{\partial\E[\ln p]}{\partial\sigma^2}&=-\frac{ND}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_n\Big\{\lVert\mathbf{x}_n-\bar{\mathbf{x}}\rVert^2-2\E[\mathbf{z}_n]^\T\mathbf{W}^\T(\mathbf{x}_n-\bar{\mathbf{x}})\\ &\qquad+\Tr(\E[\mathbf{z}_n\mathbf{z}_n^\T]\mathbf{W}^\T\mathbf{W})\Big\}.\end{aligned} \]

    The first vanishes at (12.56), a maximum because (12.53) is a concave quadratic in \(\mathbf{W}\); the second vanishes at (12.57) with \(\mathbf{W}=\mathbf{W}_{\rm new}\). The \(\frac M2\ln(2\pi)\) term missing from the first printing of (12.53) is constant and does not matter. Full derivation, including why the stationary point is the maximum: Gap 12.4.

16. Exercise 12.16 (★★★). Derive the EM algorithm for PPCA when some data values are missing at random, treating both the \(\mathbf{z}_n\) and the missing values as latent; show that it reduces to Section 12.2.2 when nothing is missing.

    Solution

    Under MAR the missingness mechanism factors out of the likelihood, which becomes \(\prod_n\N(\mathbf{x}_n^o\given\bm{\mu}_o,\mathbf{W}_o\mathbf{W}_o^\T+\sigma^2\mathbf{I})\). The posterior over \((\mathbf{z}_n,\mathbf{x}_n^m)\) factorises as \(p(\mathbf{z}_n\given\mathbf{x}_n^o)p(\mathbf{x}_n^m\given\mathbf{z}_n)\), the first factor being (12.42) computed with the observed rows \(\mathbf{W}_o\). The M step is one least-squares problem per row \((\widetilde{\mathbf{w}}_d^\T,\mu_d)\), with missing entries contributing \(\E[\widetilde{\mathbf{z}}_n\widetilde{\mathbf{z}}_n^\T]\) times the old row, followed by a \(\sigma^2\) update that includes the noise variance of the imputed entries. With no missing values and \(\bm{\mu}=\bar{\mathbf{x}}\), the cross moments separate and the updates become (12.56)–(12.57). All formulas and the reduction are in Gap 12.5.

17. Exercise 12.17 (★★, www). For \(J\)\({}=\sum_n\lVert\mathbf{x}_n-\bm{\mu}-\mathbf{W}\mathbf{z}_n\rVert^2\) (12.95), show that minimising over \(\bm{\mu}\) reduces the problem to centred variables, that minimising over \(\mathbf{z}_n\) with \(\mathbf{W}\) fixed gives the E step (12.58), and that minimising over \(\mathbf{W}\) with \(\{\mathbf{z}_n\}\) fixed gives the M step (12.59).

    Solution

    \(\bm{\mu}\). \(\partial J/\partial\bm{\mu}\)\({}=-2\sum_n(\mathbf{x}_n-\bm{\mu}-\mathbf{W}\mathbf{z}_n)\)\({}=\mathbf{0}\) gives \(\bm{\mu}\)\({}=\bar{\mathbf{x}}-\mathbf{W}\bar{\mathbf{z}}\), and substituting, \(J\)\({}=\sum_n\lVert(\mathbf{x}_n-\bar{\mathbf{x}})-\mathbf{W}(\mathbf{z}_n-\bar{\mathbf{z}})\rVert^2\): the same cost in centred variables. With \(\widetilde{\mathbf{X}}\) the \(N\times D\) matrix of rows \((\mathbf{x}_n-\bar{\mathbf{x}})^\T\) and \(\mathbf{Z}\) the \(N\times M\) matrix of rows \((\mathbf{z}_n-\bar{\mathbf{z}})^\T\) (the Svensén solution says \(D\times M\), a slip), \(J\)\({}=\Tr\{(\widetilde{\mathbf{X}}-\mathbf{Z}\mathbf{W}^\T)(\widetilde{\mathbf{X}}-\mathbf{Z}\mathbf{W}^\T)^\T\}\).

    E step. \(\partial J/\partial\mathbf{Z}\)\({}=-2(\widetilde{\mathbf{X}}-\mathbf{Z}\mathbf{W}^\T)\mathbf{W}\)\({}=\mathbf{0}\) gives \(\mathbf{Z}\)\({}=\widetilde{\mathbf{X}}\mathbf{W}(\mathbf{W}^\T\mathbf{W})^{-1}\), i.e. \(\bm{\Omega}\)\({}=\mathbf{Z}^\T\)\({}=(\mathbf{W}^\T\mathbf{W})^{-1}\mathbf{W}^\T\widetilde{\mathbf{X}}^\T\), which is (12.58) (with \(\bm{\Omega}\) the \(M\times N\) matrix whose columns are the \(\mathbf{z}_n\); these automatically have zero mean because \(\widetilde{\mathbf{X}}\) does).

    M step. \(\partial J/\partial\mathbf{W}\)\({}=-2(\widetilde{\mathbf{X}}-\mathbf{Z}\mathbf{W}^\T)^\T\mathbf{Z}\)\({}=\mathbf{0}\) gives \(\mathbf{W}\)\({}=\widetilde{\mathbf{X}}^\T\mathbf{Z}(\mathbf{Z}^\T\mathbf{Z})^{-1}\)\({}=\widetilde{\mathbf{X}}^\T\bm{\Omega}^\T(\bm{\Omega}\bm{\Omega}^\T)^{-1}\), which is (12.59). In each case \(J\) is a convex quadratic in the block being optimised (a sum of squares of affine functions), so the stationary point is the minimum and \(J\) never increases (Result 12.12).

18. Exercise 12.18 (★). Find the number of independent parameters of the factor analysis model.

    Solution

    The covariance \(\mathbf{W}\mathbf{W}^\T+\bm{\Psi}\) has \(DM\) entries in \(\mathbf{W}\) and \(D\) in \(\bm{\Psi}\), and, exactly as in Result 12.10, \(\mathbf{W}\) and \(\mathbf{W}\mathbf{R}\) give the same covariance, removing \(M(M-1)/2\) parameters: \(DM+D-M(M-1)/2\) for the covariance (plus \(D\) for the mean). It differs from (12.51) by replacing the 1 of \(\sigma^2\) with \(D\) uniquenesses. (We confirmed this count as the rank of the Jacobian of \((\mathbf{W},\bm{\Psi})\mapsto\mathbf{C}\) at random parameter values.) The model is a genuine restriction of the full covariance only if this does not exceed \(D(D+1)/2\); the difference is \(\frac12\{(D-M)^2-(D+M)\}\), so one needs \((D-M)^2\ge D+M\) (the Ledermann bound): for example at most \(M=2\) factors for \(D=5\) variables.

19. Exercise 12.19 (★★, www). Show that factor analysis is invariant under rotations of the latent space coordinates.

    Solution

    Let \(\mathbf{R}\) be orthogonal, \(\widetilde{\mathbf{z}}\)\({}=\mathbf{R}^\T\mathbf{z}\) and \(\widetilde{\mathbf{W}}\)\({}=\mathbf{W}\mathbf{R}\). Then \(\widetilde{\mathbf{z}}\sim\N(\mathbf{0},\mathbf{R}^\T\mathbf{R})\)\({}=\N(\mathbf{0},\mathbf{I})\), the same prior, and \(\widetilde{\mathbf{W}}\widetilde{\mathbf{z}}\)\({}=\mathbf{W}\mathbf{R}\mathbf{R}^\T\mathbf{z}\)\({}=\mathbf{W}\mathbf{z}\), so \(p(\mathbf{x}\given\widetilde{\mathbf{z}})\) with loadings \(\widetilde{\mathbf{W}}\) is the same as \(p(\mathbf{x}\given\mathbf{z})\) with \(\mathbf{W}\): the joint distribution of the observable \(\mathbf{x}\) and the (relabelled) latent variable is unchanged. In particular \(p(\mathbf{x})\)\({}=\N(\bm{\mu},\widetilde{\mathbf{W}}\widetilde{\mathbf{W}}^\T+\bm{\Psi})\)\({}=\N(\bm{\mu},\mathbf{W}\mathbf{W}^\T+\bm{\Psi})\), so the likelihood cannot distinguish \(\mathbf{W}\) from \(\mathbf{W}\mathbf{R}\).

20. Exercise 12.20 (★★). Show that the only stationary point of the factor analysis likelihood with respect to \(\bm{\mu}\) is \(\bar{\mathbf{x}}\), and that it is a maximum.

    Solution

    The log likelihood is (12.43) with \(\mathbf{C}\)\({}=\mathbf{W}\mathbf{W}^\T+\bm{\Psi}\). This \(\mathbf{C}\) is positive definite because \(\mathbf{v}^\T\mathbf{C}\mathbf{v}\)\({}=\lVert\mathbf{W}^\T\mathbf{v}\rVert^2+\sum_d\psi_dv_d^2\)\({}\gt0\) for \(\mathbf{v}\ne\mathbf{0}\). The gradient \(N\mathbf{C}^{-1}(\bar{\mathbf{x}}-\bm{\mu})\) vanishes only at \(\bar{\mathbf{x}}\), and the Hessian \(-N\mathbf{C}^{-1}\) is negative definite, so this is the unique maximum (the argument of Result 12.7).

21. Exercise 12.21 (★★). Derive the E-step formulae (12.66) and (12.67) for factor analysis.

    Solution

    The posterior is linear-Gaussian with prior precision \(\mathbf{I}\) and noise precision \(\bm{\Psi}^{-1}\): its covariance is \((\mathbf{I}+\mathbf{W}^\T\bm{\Psi}^{-1}\mathbf{W})^{-1}\)\({}=\mathbf{G}\) (12.68) and its mean \(\mathbf{G}\mathbf{W}^\T\bm{\Psi}^{-1}(\mathbf{x}_n-\bar{\mathbf{x}})\) (12.66), with \(\bm{\mu}\) replaced by \(\bar{\mathbf{x}}\) thanks to Exercise 12.20. Then \(\E[\mathbf{z}_n\mathbf{z}_n^\T]\)\({}=\mathbf{G}+\E[\mathbf{z}_n]\E[\mathbf{z}_n]^\T\) (12.67). See Gap 12.7, which also gives the equivalent forms \(\mathbf{W}^\T\mathbf{C}^{-1}\) and \(\mathbf{I}-\mathbf{W}^\T\mathbf{C}^{-1}\mathbf{W}\).

22. Exercise 12.22 (★★). Write down the expected complete-data log likelihood for factor analysis and derive the M step (12.69)–(12.70).

    Solution

    With \(\mathbf{y}_n\)\({}=\mathbf{x}_n-\bar{\mathbf{x}}\),

    \[ \mathcal{Q}=-\frac N2\ln\lvert\bm{\Psi}\rvert-\frac12\sum_n\Tr\Big(\bm{\Psi}^{-1}\E\big[(\mathbf{y}_n-\mathbf{W}\mathbf{z}_n)(\mathbf{y}_n-\mathbf{W}\mathbf{z}_n)^\T\big]\Big)-\frac12\sum_n\Tr\E[\mathbf{z}_n\mathbf{z}_n^\T]+\const . \]

    The \(\mathbf{W}\) gradient vanishes at (12.69), independently of \(\bm{\Psi}\); maximising \(-\frac N2\ln\psi_d-\frac{T_{dd}}{2\psi_d}\) over each \(\psi_d\) gives \(\psi_d=T_{dd}/N\), and simplifying \(\mathbf{T}\) at \(\mathbf{W}_{\rm new}\) gives (12.70). Details in Gap 12.7.

23. Exercise 12.23 (★, www). Draw the graphical model of a mixture of PPCA models with component-specific \(\mathbf{W}\), \(\bm{\mu}\), \(\sigma^2\), and the modified graph in which these parameters are shared.

    Solution

    Left: each point has a discrete component indicator \(s_n\) (one-of-\(K\), with prior \(\pi\)) and a continuous \(\mathbf{z}_n\sim\N(\mathbf{0},\mathbf{I})\); \(p(\mathbf{x}_n\given s_n=k,\mathbf{z}_n)\)\({}=\N(\mathbf{x}_n\given\mathbf{W}_k\mathbf{z}_n+\bm{\mu}_k,\sigma_k^2\mathbf{I})\), with the component parameters in a plate over \(k\). Right: \(\mathbf{W},\bm{\mu},\sigma^2\) shared by all components. Then \(p(\mathbf{x}_n\given s_n,\mathbf{z}_n)\) does not depend on \(s_n\), the indicator is irrelevant, and the model is a single PPCA: sharing everything is useless. Useful intermediate models share only some parameters (for example \(\sigma^2\), or \(\bm{\mu}\)).

24. Exercise 12.24 (★★★). Using the representation of Student's t-distribution as an infinite mixture of Gaussians (Section 2.3.7), derive an EM algorithm for the maximum likelihood fit of a multivariate t-distribution.

    Solution

    Latent variable representation (Section 2.3.7; toolkit):

    \[ \St(\mathbf{x}\given\bm{\mu},\bm{\Lambda},\nu)=\int_0^\infty\N\big(\mathbf{x}\given\bm{\mu},(\eta\bm{\Lambda})^{-1}\big)\Gam(\eta\given\nu/2,\nu/2)\,\d\eta . \]

    Give each point a latent precision scale \(\eta_n\). With \(\Delta_n^2\)\({}=(\mathbf{x}_n-\bm{\mu})^\T\bm{\Lambda}(\mathbf{x}_n-\bm{\mu})\), the complete-data log likelihood is

    \[ \sum_n\Big\{\frac D2\ln\eta_n+\frac12\ln\lvert\bm{\Lambda}\rvert-\frac{\eta_n}{2}\Delta_n^2+\frac\nu2\ln\frac\nu2-\ln\Gamma\Big(\frac\nu2\Big)+\Big(\frac\nu2-1\Big)\ln\eta_n-\frac\nu2\eta_n\Big\}+\const . \]

    E step. As a function of \(\eta_n\) the posterior is \(\propto\eta_n^{D/2}e^{-\eta_n\Delta_n^2/2}\eta_n^{\nu/2-1}e^{-\nu\eta_n/2}\), so \(\eta_n\given\mathbf{x}_n\sim\Gam\big(\frac{\nu+D}{2},\frac{\nu+\Delta_n^2}{2}\big)\) and

    \[ \E[\eta_n]=\frac{\nu+D}{\nu+\Delta_n^2},\qquad\E[\ln\eta_n]=\psi\Big(\frac{\nu+D}{2}\Big)-\ln\frac{\nu+\Delta_n^2}{2}, \]

    with \(\psi\) the digamma function (toolkit), all evaluated at the old parameters. Points far from the centre (large \(\Delta_n^2\)) get small weights \(\E[\eta_n]\): this is the robustness of the t-distribution.

    M step. The expected complete-data log likelihood is a weighted Gaussian log likelihood plus a term in \(\nu\). Setting derivatives to zero: \(\sum_n\E[\eta_n]\bm{\Lambda}(\mathbf{x}_n-\bm{\mu})\)\({}=\mathbf{0}\) and \(\frac N2\bm{\Lambda}^{-1}-\frac12\sum_n\E[\eta_n](\mathbf{x}_n-\bm{\mu})(\mathbf{x}_n-\bm{\mu})^\T\)\({}=\mathbf{0}\) give

    \[ \bm{\mu}_{\rm new}=\frac{\sum_n\E[\eta_n]\mathbf{x}_n}{\sum_n\E[\eta_n]},\qquad\bm{\Lambda}_{\rm new}^{-1}=\frac1N\sum_n\E[\eta_n](\mathbf{x}_n-\bm{\mu}_{\rm new})(\mathbf{x}_n-\bm{\mu}_{\rm new})^\T , \]

    and for \(\nu\), \(\frac{\partial}{\partial\nu}\) gives the one-dimensional equation

    \[ 1+\ln\frac\nu2-\psi\Big(\frac\nu2\Big)+\frac1N\sum_n\big\{\E[\ln\eta_n]-\E[\eta_n]\big\}=0 , \]

    solved numerically. It has exactly one root: \(\ln(\nu/2)-\psi(\nu/2)\) decreases from \(+\infty\) to 0 as \(\nu\) goes from 0 to \(\infty\), and the last term is below \(-1\), because Jensen's inequality and \(\ln t\le t-1\) give \(\E[\ln\eta_n]\)\({}\lt\ln\E[\eta_n]\)\({}\le\E[\eta_n]-1\) (strictly, since the Gamma posterior is not degenerate). This is the algorithm of Liu and Rubin (1995); our implementation increased the log likelihood at every iteration.

25. Exercise 12.25 (★★, www). For the linear-Gaussian model with noise covariance \(\bm{\Phi}\), show that the maximum likelihood solution for the transformed data \(\mathbf{A}\mathbf{x}_n\) is \((\mathbf{A}\bm{\mu}_{\rm ML},\mathbf{A}\mathbf{W}_{\rm ML},\mathbf{A}\bm{\Phi}_{\rm ML}\mathbf{A}^\T)\), and that the form of the model is preserved for (i) diagonal \(\mathbf{A}\) and diagonal \(\bm{\Phi}\) (factor analysis) and (ii) orthogonal \(\mathbf{A}\) and \(\bm{\Phi}=\sigma^2\mathbf{I}\) (PPCA).

    Solution

    The likelihood of \(\{\mathbf{A}\mathbf{x}_n\}\) at \((\bm{\mu}',\mathbf{W}',\bm{\Phi}')\) equals the likelihood of \(\{\mathbf{x}_n\}\) at \((\mathbf{A}^{-1}\bm{\mu}',\mathbf{A}^{-1}\mathbf{W}',\mathbf{A}^{-1}\bm{\Phi}'\mathbf{A}^{-\T})\) times \(\lvert\det\mathbf{A}\rvert^{-N}\), so maximisers correspond under \((\bm{\mu},\mathbf{W},\bm{\Phi})\mapsto(\mathbf{A}\bm{\mu},\mathbf{A}\mathbf{W},\mathbf{A}\bm{\Phi}\mathbf{A}^\T)\), as long as this map preserves the class of allowed \(\bm{\Phi}\): diagonal is preserved by diagonal \(\mathbf{A}\), and \(\sigma^2\mathbf{I}\) by orthogonal \(\mathbf{A}\) (giving \(\sigma^2\mathbf{A}\mathbf{A}^\T\)\({}=\sigma^2\mathbf{I}\)). Hence factor analysis is covariant under rescaling of the variables and PPCA (like PCA) under rotations, but in general not conversely. Full proof: Result 12.13. (The official solution writes \(\bm{\Phi}_A\)\({}=\mathbf{A}\bm{\Phi}^{-1}\mathbf{A}^\T\); it should be \(\mathbf{A}\bm{\Phi}\mathbf{A}^\T\).)

26. Exercise 12.26 (★★). Show that every solution of (12.80) solves (12.79); that adding to a solution of (12.80) any null vector of \(\mathbf{K}\) gives a solution of (12.79) with the same eigenvalue; and that this does not change the projections (12.82).

    Solution

    Multiplying (12.80) by \(\mathbf{K}\) gives (12.79). If \(\mathbf{K}\mathbf{a}\)\({}=\lambda N\mathbf{a}\) and \(\mathbf{K}\mathbf{a}_0=\mathbf{0}\), then \(\mathbf{K}^2(\mathbf{a}+\mathbf{a}_0)\)\({}=\mathbf{K}^2\mathbf{a}\)\({}=\lambda N\mathbf{K}\mathbf{a}\)\({}=\lambda N\mathbf{K}(\mathbf{a}+\mathbf{a}_0)\). Conversely every solution of (12.79) has this form (Result 12.14(c)). Finally \(\lVert\sum_na_{0n}\bm{\phi}(\mathbf{x}_n)\rVert^2\)\({}=\mathbf{a}_0^\T\mathbf{K}\mathbf{a}_0\)\({}=0\), so \(\sum_na_{0n}\bm{\phi}(\mathbf{x}_n)\)\({}=\mathbf{0}\) and \(\sum_na_{0n}k(\mathbf{x},\mathbf{x}_n)\)\({}=\bm{\phi}(\mathbf{x})^\T\mathbf{0}\)\({}=0\) for every \(\mathbf{x}\): the projection (12.82) is unchanged.

27. Exercise 12.27 (★★). Show that kernel PCA with the linear kernel \(k(\mathbf{x},\mathbf{x}')\)\({}=\mathbf{x}^\T\mathbf{x}'\) is conventional PCA.

    Solution

    With \(\bm{\phi}(\mathbf{x})=\mathbf{x}\), the feature space is the data space and the feature covariance (12.73) of the centred features is \(\mathbf{S}\), so the eigenvectors \(\mathbf{v}_i\) of (12.74) are the principal components \(\mathbf{u}_i\). In kernel form: \(\mathbf{K}=\mathbf{X}\mathbf{X}^\T\) and, by Result 12.15, \(\widetilde{\mathbf{K}}\)\({}=\mathbf{H}\mathbf{X}\mathbf{X}^\T\mathbf{H}\)\({}=\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^\T\). Then (12.80) reads \(\frac1N\widetilde{\mathbf{X}}\widetilde{\mathbf{X}}^\T\mathbf{a}_i\)\({}=\lambda_i\mathbf{a}_i\), which is (12.28), so \(\mathbf{a}_i\propto\mathbf{v}_i\); the normalisation (12.81) gives \(\mathbf{a}_i\)\({}=\mathbf{v}_i/(N\lambda_i)^{1/2}\); and the feature-space eigenvector \(\sum_na_{in}(\mathbf{x}_n-\bar{\mathbf{x}})\)\({}=\widetilde{\mathbf{X}}^\T\mathbf{a}_i\)\({}=\widetilde{\mathbf{X}}^\T\mathbf{v}_i/(N\lambda_i)^{1/2}\)\({}=\mathbf{u}_i\) by (12.30). The projections are \(y_i(\mathbf{x})\)\({}=\mathbf{u}_i^\T(\mathbf{x}-\bar{\mathbf{x}})\), the PCA scores. (Checked numerically.)

28. Exercise 12.28 (★★, www). Using (1.27), show that any density \(p(y)\) can be obtained from a fixed density \(q(x)\) that is everywhere non-zero by a monotonic change of variable \(y=f(x)\); write down the differential equation for \(f\) and sketch the transformation.

    Solution

    For strictly increasing differentiable \(f\), \(P(Y\le f(x))=P(X\le x)\), i.e. \(F_p(f(x))=F_q(x)\), and differentiating gives the differential equation \(p(f(x))\,f'(x)=q(x)\); equivalently \(p(y)\)\({}=q(f^{-1}(y))\,\frac{\d f^{-1}}{\d y}\), which is (1.27). Its solution is \(f=F_p^{-1}\circ F_q\). Because \(q\gt0\) everywhere, \(F_q\) is continuous and strictly increasing, so \(f\) is well defined and increasing for any target \(p\) (where \(p\) vanishes on an interval, \(F_p^{-1}\) and hence \(f\) jump across it; \(f'=q/p\) is finite wherever \(p(f(x))\gt0\)). The sketch: plot \(q\) along the \(x\) axis and \(p\) along the \(y\) axis; the curve \(y=f(x)\) maps each interval \([x,x+\d x]\) to \([f(x),f(x)+f'(x)\d x]\), carrying the same probability mass, so \(f\) is steep where \(q\) is large and \(p\) small, and flat where the reverse holds. Combined with a Gaussian latent variable, this is how nonlinear latent variable models produce non-Gaussian densities (Section 12.4).

29. Exercise 12.29 (★★, www). Show that independent variables have a diagonal covariance matrix; then, for \(y_1\) symmetric about 0 and \(y_2=y_1^2\), write down \(p(y_2\given y_1)\), observe that the variables are dependent, and show that their covariance matrix is nevertheless diagonal. (Erratum: the first printing states the second part incorrectly.)

    Solution

    If \(p(z_1,z_2)=p(z_1)p(z_2)\), then \(\cov[z_1,z_2]\)\({}=\iint(z_1-\bar z_1)(z_2-\bar z_2)p(z_1)p(z_2)\,\d z_1\d z_2\)\({}=\int(z_1-\bar z_1)p(z_1)\d z_1\int(z_2-\bar z_2)p(z_2)\d z_2\)\({}=0\cdot0\). For the second part, \(p(y_2\given y_1)=\delta(y_2-y_1^2)\), a point mass whose location depends on \(y_1\), so the variables are dependent. By symmetry \(\E[y_1]=0\) and \(\E[y_1^3]=0\) (odd moments vanish; assume \(\E\lvert y_1\rvert^3\lt\infty\)), so \(\cov[y_1,y_2]\)\({}=\E[y_1y_2]-\E[y_1]\E[y_2]\)\({}=\E[y_1^3]-0\)\({}=0\). Hence the covariance matrix is diagonal although the variables are dependent: zero correlation does not imply independence. This is why decorrelating the data (PCA, whitening) cannot by itself recover independent sources, and ICA needs more than second-order statistics (Section 12.4.1).

Demonstrations for this chapter

Each demonstration runs in the browser and is tied to the results it illustrates; the ▶ Demo boxes above open them already set to the relevant case. Random data use a fixed seed, so a class sees the same pictures.

All figures of the chapter are reproduced by code/ch12_continuous_latent.py (one function per figure); numerical checks of the chapter's identities are in code/checks_ch12.py.

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 12, pp. 559–603.
2. Svensén, M. and Bishop, C. M. (2011). Pattern Recognition and Machine Learning: Errata and Additional Comments; and Svensén, M. and Bishop, C. M. (2009). Pattern Recognition and Machine Learning: Solutions to the Exercises (Web edition). Microsoft Research.
3. Ahn, J.-H. and Oh, J.-H. (2003). A constrained EM algorithm for principal component analysis. Neural Computation, 15(1), 57–65. doi:10.1162/089976603321043694.
4. Amari, S., Cichocki, A. and Yang, H. H. (1996). A new learning algorithm for blind signal separation. Advances in Neural Information Processing Systems, 8, 757–763.
5. Anderson, T. W. (1963). Asymptotic theory for principal component analysis. Annals of Mathematical Statistics, 34(1), 122–148. doi:10.1214/aoms/1177704248.
6. Attias, H. (1999). Independent factor analysis. Neural Computation, 11(4), 803–851. doi:10.1162/089976699300016458.
7. Bach, F. R. and Jordan, M. I. (2002). Kernel independent component analysis. Journal of Machine Learning Research, 3, 1–48.
8. Bakir, G. H., Weston, J. and Schölkopf, B. (2004). Learning to find pre-images. Advances in Neural Information Processing Systems, 16, 449–456.
9. Baldi, P. and Hornik, K. (1989). Neural networks and principal component analysis: learning from examples without local minima. Neural Networks, 2(1), 53–58. doi:10.1016/0893-6080(89)90014-2.
10. Bartholomew, D. J. (1987). Latent Variable Models and Factor Analysis. Charles Griffin, London.
11. Basilevsky, A. (1994). Statistical Factor Analysis and Related Methods: Theory and Applications. Wiley, New York. doi:10.1002/9780470316894.
12. Bell, A. J. and Sejnowski, T. J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6), 1129–1159. doi:10.1162/neco.1995.7.6.1129.
13. Bishop, C. M. (1999a). Bayesian PCA. Advances in Neural Information Processing Systems, 11, 382–388.
14. Bishop, C. M. (1999b). Variational principal components. Proceedings of the Ninth International Conference on Artificial Neural Networks (ICANN 99), IEE, 509–514. doi:10.1049/cp:19991160.
15. Bishop, C. M. and James, G. D. (1993). Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A, 327(2–3), 580–593. doi:10.1016/0168-9002(93)90728-Z.
16. Bishop, C. M., Svensén, M. and Williams, C. K. I. (1997b). Magnification factors for the GTM algorithm. Proceedings of the Fifth IEE International Conference on Artificial Neural Networks, 64–69. doi:10.1049/cp:19970703.
17. Bishop, C. M., Svensén, M. and Williams, C. K. I. (1998a). Developments of the generative topographic mapping. Neurocomputing, 21(1–3), 203–224. doi:10.1016/S0925-2312(98)00043-5.
18. Bishop, C. M., Svensén, M. and Williams, C. K. I. (1998b). GTM: the generative topographic mapping. Neural Computation, 10(1), 215–234. doi:10.1162/089976698300017953.
19. Bishop, C. M. and Tipping, M. E. (1998). A hierarchical latent variable model for data visualization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 20(3), 281–293. doi:10.1109/34.667885.
20. Bishop, C. M. and Winn, J. M. (2000). Non-linear Bayesian image modelling. Computer Vision, ECCV 2000, Lecture Notes in Computer Science 1842, 3–17. doi:10.1007/3-540-45054-8_1.
21. Bourlard, H. and Kamp, Y. (1988). Auto-association by multilayer perceptrons and singular value decomposition. Biological Cybernetics, 59(4–5), 291–294. doi:10.1007/BF00332918.
22. Cardoso, J.-F. (1998). Blind signal separation: statistical principles. Proceedings of the IEEE, 86(10), 2009–2025. doi:10.1109/5.720250.
23. Collins, M., Dasgupta, S. and Schapire, R. E. (2002). A generalization of principal components analysis to the exponential family. Advances in Neural Information Processing Systems, 14, 617–624.
24. Comon, P. (1994). Independent component analysis, a new concept? Signal Processing, 36(3), 287–314. doi:10.1016/0165-1684(94)90029-9.
25. Comon, P., Jutten, C. and Hérault, J. (1991). Blind separation of sources, part II: problems statement. Signal Processing, 24(1), 11–20. doi:10.1016/0165-1684(91)90080-3.
26. Cox, T. F. and Cox, M. A. A. (2000). Multidimensional Scaling, 2nd edition. Chapman and Hall/CRC, Boca Raton. doi:10.1201/9781420036121.
27. Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1), 1–38. doi:10.1111/j.2517-6161.1977.tb01600.x.
28. Erwin, E., Obermayer, K. and Schulten, K. (1992). Self-organizing maps: ordering, convergence properties and energy functions. Biological Cybernetics, 67(1), 47–55. doi:10.1007/BF00201801.
29. Everitt, B. S. (1984). An Introduction to Latent Variable Models. Chapman and Hall, London. doi:10.1007/978-94-009-5564-6.
30. Fan, K. (1949). On a theorem of Weyl concerning eigenvalues of linear transformations I. Proceedings of the National Academy of Sciences, 35(11), 652–655. doi:10.1073/pnas.35.11.652.
31. Ghahramani, Z. and Beal, M. J. (2000). Variational inference for Bayesian mixtures of factor analysers. Advances in Neural Information Processing Systems, 12, 449–455.
32. Ghahramani, Z. and Hinton, G. E. (1996a). The EM algorithm for mixtures of factor analyzers. Technical Report CRG-TR-96-1, University of Toronto.
33. Golub, G. H. and Van Loan, C. F. (1996). Matrix Computations, 3rd edition. Johns Hopkins University Press, Baltimore.
34. Hastie, T. and Stuetzle, W. (1989). Principal curves. Journal of the American Statistical Association, 84(406), 502–516. doi:10.1080/01621459.1989.10478797.
35. Hinton, G. E., Dayan, P. and Revow, M. (1997). Modeling the manifolds of images of handwritten digits. IEEE Transactions on Neural Networks, 8(1), 65–74. doi:10.1109/72.554192.
36. Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24(6), 417–441. doi:10.1037/h0071325.
37. Hotelling, H. (1936). Relations between two sets of variates. Biometrika, 28(3–4), 321–377. doi:10.1093/biomet/28.3-4.321.
38. Hyvärinen, A. and Oja, E. (1997). A fast fixed-point algorithm for independent component analysis. Neural Computation, 9(7), 1483–1492. doi:10.1162/neco.1997.9.7.1483.
39. Hyvärinen, A. and Oja, E. (2000). Independent component analysis: algorithms and applications. Neural Networks, 13(4–5), 411–430. doi:10.1016/S0893-6080(00)00026-5.
40. Ilin, A. and Raiko, T. (2010). Practical approaches to principal component analysis in the presence of missing values. Journal of Machine Learning Research, 11, 1957–2000.
41. Jolliffe, I. T. (2002). Principal Component Analysis, 2nd edition. Springer, New York. doi:10.1007/b98835.
42. Jutten, C. and Hérault, J. (1991). Blind separation of sources, part I: an adaptive algorithm based on neuromimetic architecture. Signal Processing, 24(1), 1–10. doi:10.1016/0165-1684(91)90079-X.
43. Kambhatla, N. and Leen, T. K. (1997). Dimension reduction by local principal component analysis. Neural Computation, 9(7), 1493–1516. doi:10.1162/neco.1997.9.7.1493.
44. Kohonen, T. (1982). Self-organized formation of topologically correct feature maps. Biological Cybernetics, 43(1), 59–69. doi:10.1007/BF00337288.
45. Kohonen, T. (1995). Self-Organizing Maps. Springer, Berlin. doi:10.1007/978-3-642-97610-0.
46. Lawley, D. N. (1953). A modified method of estimation in factor analysis and some large sample results. Uppsala Symposium on Psychological Factor Analysis, Nordisk Psykologi Monograph Series 3, 35–42.
47. LeCun, Y., Bottou, L., Bengio, Y. and Haffner, P. (1998). Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11), 2278–2324. doi:10.1109/5.726791. (The MNIST data.)
48. Little, R. J. A. and Rubin, D. B. (2019). Statistical Analysis with Missing Data, 3rd edition. Wiley, Hoboken. doi:10.1002/9781119482260.
49. Liu, C. and Rubin, D. B. (1995). ML estimation of the t distribution using EM and its extensions, ECM and ECME. Statistica Sinica, 5(1), 19–39.
50. MacKay, D. J. C. (1995). Bayesian neural networks and density networks. Nuclear Instruments and Methods in Physics Research A, 354(1), 73–80. doi:10.1016/0168-9002(94)00931-7.
51. MacKay, D. J. C. (2003). Information Theory, Inference, and Learning Algorithms. Cambridge University Press.
52. MacKay, D. J. C. and Gibbs, M. N. (1999). Density networks. In J. W. Kay and D. M. Titterington (eds), Statistics and Neural Networks: Advances at the Interface, Oxford University Press, 129–146.
53. Minka, T. P. (2001c). Automatic choice of dimensionality for PCA. Advances in Neural Information Processing Systems, 13, 598–604.
54. Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559–572. doi:10.1080/14786440109462720.
55. Roweis, S. (1998). EM algorithms for PCA and SPCA. Advances in Neural Information Processing Systems, 10, 626–632.
56. Roweis, S. T. and Saul, L. K. (2000). Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500), 2323–2326. doi:10.1126/science.290.5500.2323.
57. Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581–592. doi:10.1093/biomet/63.3.581.
58. Rubin, D. B. and Thayer, D. T. (1982). EM algorithms for ML factor analysis. Psychometrika, 47(1), 69–76. doi:10.1007/BF02293851.
59. Schölkopf, B., Smola, A. and Müller, K.-R. (1998). Nonlinear component analysis as a kernel eigenvalue problem. Neural Computation, 10(5), 1299–1319. doi:10.1162/089976698300017467.
60. Sirovich, L. (1987). Turbulence and the dynamics of coherent structures. I. Coherent structures. Quarterly of Applied Mathematics, 45(3), 561–571. doi:10.1090/qam/910462.
61. Stone, J. V. (2004). Independent Component Analysis: A Tutorial Introduction. MIT Press, Cambridge, MA. doi:10.7551/mitpress/3717.001.0001.
62. Tenenbaum, J. B., de Silva, V. and Langford, J. C. (2000). A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500), 2319–2323. doi:10.1126/science.290.5500.2319.
63. Tiño, P. and Nabney, I. T. (2002). Hierarchical GTM: constructing localized nonlinear projection manifolds in a principled way. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(5), 639–656. doi:10.1109/34.1000238.
64. Tiño, P., Nabney, I. T. and Sun, Y. (2001). Using directional curvatures to visualize folding patterns of the GTM projection manifolds. Artificial Neural Networks, ICANN 2001, Lecture Notes in Computer Science 2130, 421–428. doi:10.1007/3-540-44668-0_59.
65. Tipping, M. E. (1999). Probabilistic visualisation of high-dimensional binary data. Advances in Neural Information Processing Systems, 11, 592–598.
66. Tipping, M. E. and Bishop, C. M. (1999a). Mixtures of probabilistic principal component analyzers. Neural Computation, 11(2), 443–482. doi:10.1162/089976699300016728.
67. Tipping, M. E. and Bishop, C. M. (1999b). Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B, 61(3), 611–622. doi:10.1111/1467-9868.00196.
68. Torgerson, W. S. (1952). Multidimensional scaling: I. Theory and method. Psychometrika, 17(4), 401–419. doi:10.1007/BF02288916.

Numbers quoted from "our" computations (the MNIST eigenvalues, the EM, missing-data, Bayesian PCA, saddle-point, autoencoder and other checks) were produced with numpy/scipy while writing these notes; every displayed identity was also checked numerically.