Model Manifolds of Probability Distributions

Any measurement obtained from an experiment with uncertainty can generally be understood as a probability distribution. For example, when some data x are observed with normally distributed noise ξ of variance σ2, under experimental conditions θj, a model is expressed asx=f(θj)+ξ where L(ξ)∽ N(0,σ2),[1]and f(θi) is a prediction given the experimental conditions. This relationship is equivalent to saying that the probability of measuring data x given some conditions θ isL(x∣θ)∽ N(f(θ),σ2).[2]More complicated noise profiles with asymmetry or correlations can be accommodated with this picture. Measurements without an underlying model can also be seen as distributions, where a measurement xi with uncertainty σ can induce a probability Lx∣xi,σ of observing new data x.

We define a probabilistic model Lx∣θ, the likelihood of observing data x given parameters θ. The model manifold is defined as the set of all possible predictions, {Lx∣θi}, which is a surface parameterized by the model parameters {θi}. The parameter directions related to the longest distances along the model manifold have been shown to predict emergent behavior (how microscopic parameters lead to macroscopic behavior) (13). We will see that InPCA orders its principal components by the length of the model manifold along their direction, highlighting global structure. The boundaries of the model manifold represent simplified models which retain predictive power (14), and the constraint of data lying near the model manifold has been used to optimize experimental design (15). In this article, we study the Ising model, which defines probabilities of spin configurations given interaction strengths; a neural network, which predicts the probability of an image representing a single handwritten digit given weights and biases; and ΛCDM, which predicts the distribution of CMB radiation given fundamental constants of nature.

Hypersphere Embedding

We promised an embedding which both is isometric and preserves global structures. We satisfy the first promise by considering the hypersphere embedding{zx(θi)}=2Lx∣θi,[3]where the normalization constraint of Lx∣θ forces zx to lie on the positive orthant of a sphere. A natural measure of distance on the hypersphere is the Euclidean distance, in this case also known as the Hellinger divergence (16)d2(θ1,θ2)=z(θ1)−z(θ2)2=81−Lx∣θ1⋅Lx∣θ22,[4]where ⋅ represents the inner product over x. Now we can see that the hypersphere embedding is isometric: The Euclidean metric of this embedding is equal to the Fisher information metric I of the model manifold (17),d2(zi,zi+dzi)=∑idzidzi=∑klIkldθkdθl.[5]The Fisher information metric (FIM) is the natural metric of the model manifold (18), so the hypersphere embedding preserves the local structure of the manifold defined by Lx∣θ.

As the dimension of the data increases, almost all features become orthogonal to each other, and most measures of distance lose their ability to discriminate between the smallest and largest distances (19). For the hypersphere embedding, we see that as the dimension of x increases, the inner product in the Hellinger distance of Eq. 4 becomes smaller as the probability is distributed over more dimensions. In the limit of large dimension, all nonidentical pairs of points become orthogonal and equidistant around the hypersphere (a constant distance 8 apart), frustrating effective dimensional reductions and visualization.

To illustrate this problem with the hypersphere embedding, consider the Ising model, which predicts the likelihood of observing a particular configuration of binary random variables (spins) on a lattice. The probability of a spin configuration is determined by the Boltzmann distribution and is a function of a local pairwise coupling and a global applied field. The dimension is determined by the number of spin configurations, 2N, where N is the number of spins. Holding temperature fixed at one, we vary h and J: external magnetic field (h∈(−1.3,1.3)) and nearest-neighbor coupling (J∈(−0.4,0.6)), using a Monte Carlo method weighted by Jeffrey’s prior to sample 12,000 distinct points. From the resulting set of parameters, we compute Xij={zi(θj)} using the Boltzmann distribution and visualize the model manifold in the N-sphere embedding of Eq. 3 by projecting the predictions onto the first three principal components of X. Fig. 1A shows this projection of the model manifold of a 2×2 Ising model which is embedded in 24 dimensions. Fig. 1B shows a larger, 4×4 Ising model, of dimension 216. As the dimension is increased from 24 to 216, we see the points starting to wrap around the hypersphere, becoming increasingly equidistant and less distinguishable.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

(A–C) Hypersphere embedding, illustrating an embedding of the 2D Ising model. Points were generated through a Monte Carlo sampling and visualized by projecting the probability distributions onto the first three principal components (28). The points are colored by magnetic field strength. As the system size increases from 2×2 to 4×4, the orthogonality problem is demonstrated by an increase in “wrapping” around the hypersphere. This effect can also be produced by instead considering four replicas of the original system, motivating the replica trick which takes the embedding dimension or number of replicas to zero.

A natural way to increase the dimensionality of a probabilistic model is by drawing multiple samples from the distribution. If D is the dimension of x, then N identical draws from the distribution will have dimension DN. The more samples drawn, the easier it is to distinguish between distributions, mimicking the curse of dimensionality for large systems. We see this demonstrated for our Ising model in Fig. 1C, where we drew four replica samples from the same model. Note that compared with the original 2 × 2 model, the model manifold of the four-replica 2 × 2 model “wraps” more around the hypersphere, just like the larger, 4×4 Ising model. High-dimensional systems have “too much information,” in the same way that large numbers of samples have too much information. In the next section, we consider the contraposition of the insight that a large number of replicas lead to the curse of dimensionality and discover an embedding which not only is isometric but also ameliorates the high-dimensional wrapping around the n-sphere.

Replica Theory and the Intensive Embedding

We saw in Fig. 1 that increasing the dimension of the data led to a saturation of the distance function Eq. 4. This problem is referred to as the loss of relative contrast or the concentration of distances (19), and to overcome it requires a non-Euclidean distance function, discussed below. In the previous section we saw the same saturation of distance could be achieved by adding replicas, increasing the embedding dimension. Fig. 2A shows this process taken to an extreme: the model manifold of the 2×2 Ising model with the number of replicas taken to infinity. All of the points cluster together, obscuring the fact that the underlying manifold is 2D. To cure the abundance of information which makes all points on the hypersphere equidistant, we seek an intensive distance, such as the distance per number of replicas observed. Next, because the limit of many replicas artificially leads to the same symptoms of the curse of dimensionality, we consider the limit of zero replicas, a procedure which is often used in the study of spin glasses and disordered systems (20). Fig. 2B shows the result of this analysis, the intensive embedding, where the distance concentration has been cured, and the inherent 2D structure of the Ising model has been recovered.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Replicated Ising model illustrating the derivation of our intensive embedding. All points are colored by magnetic field strength. (A) Large dimensions are characterized by large system sizes; here we mimic a 128×128 Ising model which is of dimension 21282. The orthogonality problem becomes manifest as all points are effectively orthogonal, producing a useless visualization with all points clustered in the cusp. (B) Using replica theory, we tune the dimensionality of the system and consider the limit as the number of replicas goes to zero. In this way, we derive our intensive embedding. Note that the z axis reflects a negative-squared distance, a property which allows violations of the triangle inequality and is discussed in the text.

To find the intensive embedding, we must first find the distance between replicated models. The likelihood for N replicas of a system is given by their productL{x1,…,xN}∣θ(N)=Lx1∣θ⋯LxN∣θ,[6]where the set {x1,…,xN} represents the observed data in the replicated systems. Writing the inner product or cosine angle between two distributions asθ1;θ2=Lx∣θ1⋅Lx∣θ2,[7]and using Eq. 4, the distance per replica dN2 between two points on the model manifold isdN2(θ1,θ2)=d2(θ1,θ2)N=−8θ1;θ2N−1N.[8]We are now poised to define the intensive distance by taking the number of replicas to zero:dI2(θ1,θ2)=limN→0dN2(θ1,θ2)=−8⁡logθ1;θ2.[9]The last equality is achieved using the standard trick in replica theory, (xN−1)/N→logx as N→∞, which is a basis trick used to solve challenging problems in statistical physics (20). The trick is most evident using the identity xN=exp(log⁡Nx)≈1+N⁡log⁡x. One can check that the intensive distance is isometric,dI2(θ,θ+δθ)=δθαδθβgαβ=δθαδθβIαβ,[10]where again I is the Fisher information metric in Eq. 5, so that we can be confident the intensive embedding distance preserves local structures.

Importantly, the intensive distance does not satisfy the triangle inequality (and is thus non-Euclidean): The distance between points on the hypersphere can go to infinity, rather than lie constrained to the finite radius of the hypersphere embedding. Because of this, the intensive embedding can overcome the loss of relative contrast (19) discussed at the beginning of this section. Distances in the intensive embedding maintain distinguishability in high dimensions, as illustrated in Fig. 2B, wherein the 2D nature of the Ising model has been recovered. We hypothesize that this process, which cures the curse of dimensionality for models with too many samples, will also cure it for models with intrinsically high dimensionality. The intensive distance obtained here is proportional to the Bhattacharyya distance (21). Considering the zero-replica limit of the Hellinger divergence, we discovered a way to derive the Bhattacharyya distance. The importance of this is discussed further in the following section.

Intensive Principal Component Analysis

Classical PCA takes a set of data examples and infers features which are linearly uncorrelated. (7). The features to be analyzed with PCA are compared via their Euclidean distance. Can we generalize this comparison to use our intensive embedding distance? Given a matrix of data examples X∈Rm×p (with features along the rows), PCA first requires the mean-shifted matrix Mij=Xij−X¯i=PX, where Pij=δij−1/p is the mean-shift projection matrix and p is the number of sampled points. The covariance and its eigenvalue decomposition are thencov(X,X)=1pMTM=XTPPX=VΣVT,[13]where the orthogonal columns of the matrix V are the natural basis onto which the rows of M are projected,MV=(UDVT)V=UD=UΣ,[14]where the columns of UΣ are called the principal components of the data X.

The principal components can also be obtained from the cross-covariance matrix, MMT, sinceMMT=PXXTP=(UDVT)(UDVT)T=UΣUT.[15]The eigenbasis U of the cross-covariance is the natural basis for the components of the data, and the eigenbasis V of the covariance is the natural basis of the data points. For us this flexibility is invaluable, as the cross-covariance is more natural for expressing the distances between distributions of different parameters.

Writing our data matrix as Xij=zi(θj) using Eq. 3 for replicated systems, the cross-covariance is(MMT)ij(N)=(PXXTP)ij=(z(θi)−z¯)⋅(z(θj)−z¯)=4θi;θjN+4p2∑k,k′=1pθk;θk′N −4p∑k=1pθi;θkN+θj;θkN,[16]where z¯ is the average over all sampled parameters, and we used the definition of z in Eq. 6. As with the intensive embedding, we can take the limit as the number of replicas goes to zero to findWij=limN→01N(MMT)ij(N).[17]Explicitly, the intensive cross-covariance matrix isWij=4⁡logθi;θj+4p2∑k,k′=1plogθk′;θk −4p∑k=1plogθi;θk+logθj;θk[18]=(PLP)ij,[19]where Li,j=4⁡logθi;θj and P is the same projection matrix as defined above. In taking the limit of zero replicas, the structure of the cross-covariance has transformedPXXTP→N→0PLP,[20]and thus the symmetric Wishart structure is lost. It is therefore possible to obtain negative eigenvalues in this decomposition, which give rise to imaginary components in the projections. Note the similarity between the form of this cross-covariance and the double-centered distance matrix used in PCA and multidimensional scalding (MDS). This arises because both InPCA and PCA/MDS rely on mean shifing the input data before finding an eigenbasis. Thus, we view InPCA as a natural generalization of PCA to probability distributions and MDS to non-Euclidean embeddings.

In summary, InPCA is achieved by the following procedure: (i) Compute the cross-covariance matrix from a set of probability samples: Compute Wij as derived in Eq. 18. (ii) Compute the eigenvalue decomposition W=UΣUT. (iii) Compute the coordinate projections, T=UΣ. (iv) Plot the projections using the columns of T.

Neural Network MNIST Digit Classifier.

To demonstrate the utility of InPCA, we use it to visualize the training of a two-layer convolution neural network (CNN), constructed using TensorFlow (22), trained on the MNIST dataset of hand-written digits (23). A set of 55,000 images was used to train the network, which was then used to predict the likelihood that an additional set of 10,000 images is classified each as a specific digit between 0 and 9. We use softmax (24) to calculate the probabilities from the category estimates supplied by the network. The CNN defines the likelihood Lx∣θ that some input image θ contains the image of a particular handwritten digit x. The InPCA projections of the CNN output in Fig. 3 visualize the clustering learned by the CNN as a function of the number of learning epochs. The initialized network’s model manifold shows no knowledge of the digits (colored dots), but as training commences, the network separates digits into separate regions of its manifold (Movie S1). InPCA can be used as a fast, interpretable, and deterministic method for qualitatively evaluating what a neural network has learned.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Stages of training a CNN. Each point in the 3D projections represents one of 10,000 test images supplied to the CNN (29). At the first epoch, the neural network is untrained and so is unable to reliably classify images, with about a 90% error rate—an effect reflected in the cloud of points. As training progresses and error rate decreases, the cloud begins to cluster as shown by InPCA at the 20th epoch. Finally, when completely trained, the clustered regions are manifest at the 2000th epoch with 10 clusters representing the 10 digits.

Properties of the Intensive Embedding and InPCA

The space characterized by our intensive embedding has two weird properties: First, it is formally 1D, yet there are multiple orthogonal directions upon which it can be projected; and second, it is Minkowski-like, in that it has negative squared distances, violating the triangle inequality. We posit that, fundamentally, this second property is what allows InPCA to cure the orthogonality problem.

We begin with a discussion of the the 1D nature of the embedding space. The embedding dimension is given by DN, where D is the original dimension of data x and N is the number of replicas. In the case of noninteger replicas the space becomes “fractional” in dimension and in the limit of zero replicas ultimately goes to one. However, it is still possible to obtain projections themselves along the dominant components of this space, by leveraging the cross-covariance instead of the covariance, summarized in step ii of our algorithm. Visualizations produced by InPCA are cross-sections of a space of the dimension equal to the number of sampled points of the model manifold p, instead of the dimension D or DN.

In the limit of zero replicas in Eq. 18, the positive-definite, Wishart structure of the cross-covariance matrix is lost. It is therefore possible to have negative squared distances. The non-Euclidean nature of the embedding (flat, but Minkowski-like) does not suffer from the concentration of distances which plagues Euclidean measures in high dimensions, thus allowing the model manifold to be “unwound” from the N-sphere and for InPCA to produce useful, low-dimensional representations.

Finally, the eigenvalues of InPCA correspond to the cross-sectional widths of the model manifold. We see this quite explicitly with the following example of a biased coin (specifically, in Fig. 4B) where the eigenvalues extracted from InPCA map directly to the manifold widths measured along the direction of the corresponding InPCA eigenvector. Therefore, we see that InPCA produces a hierarchy of directions, ordered by the global widths of the model manifold. Note that, as with classical PCA, this correspondence depends on how faithfully the model manifold was originally sampled; that is, InPCA can tell you about the structure of the manifold only from observed points.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

InPCA visualization of biased coins (30). (A) The first two InPCA components correspond to the coin bias and variance, yet the first one is real and the second one is imaginary (the aspect ratio between axes is one). The contour lines represent constant distances from a fair coin and are hyperbolas: Points can be a finite distance from a fair coin yet an infinite distance from each other. (B) The ordered eigenvalues correspond to the manifold lengths, illustrating the hierarchical nature of the components extracted from InPCA.

Comparing with t-SNE and Diffusion Maps.

We compare our manifold learning technique to two standard methods, t-SNE and the diffusion maps, by applying each one to the six-parameter ΛCDM cosmological model predictions of the CMB. The ΛCDM predicts Lx∣θ, where x represents fluctuations in the CMB, and θ are the different cosmological parameters (i.e., it predicts the angular power spectrum of temperature and polarization anisotropies in sky maps of the CMB). Observations of the CMB from telescopes on satellites, balloons, and the ground provide thousands of independent measurements from large angular scales to a few arcminutes, which are used to fit model parameters. Here we consider only CMB observations from the 2015 Planck data release (25). The ΛCDM model we consider has six parameters, the Hubble constant (H0) which we sampled in a range of 20–100 km⋅s−1⋅Mpc−1, the physical baryon density (Ωbh2) and the physical cold dark matter density (Ωch2) both sampled from 0.0009 to 0.8, the primordial fluctuation amplitude (As) sampled from 10−11 to 10−8, the scalar spectral index (η) sampled from 0 to 0.98, and the optical depth at reionization (τ) sampled from 0.001 to 0.9.

To determine the likelihood functions, we use the Code for Anisotropies in the Microwave Background software package to generate power spectra (26). We perform a Monte Carlo sampling of 50,000 points around the best-fit parameters provided by the 2015 Planck data release (25), with sample weights based on the intensive distance to the best fit.

In Fig. 5 we show the first two components of the manifold embedding for InPCA, t-SNE, and diffusion maps. To apply t-SNE and the diffusion map to probabilistic data we must provide a distance. We therefore use our intensive distance, from Eq. 9, for consistency and ease of comparison. In all three cases, the first component from each method is directly related to the primordial fluctuation amplitude As, which reflects the amplitude of density fluctuations in the early universe and is the dominant feature in real data (25). The second InPCA component predicts the Hubble constant, whereas the diffusion map predicts the scalar spectral index (a reflection of the size variance of primordial density fluctuations). In all cases, the projected components were plotted against the corresponding parameters to determine correlations, such as how one can see that As corresponds with the first component in all three cases. Detailed plots and correlation coefficients for all three methods are provided in SI Appendix.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Model manifold of the six-parameter ΛCDM cosmological model predictions of temperature and polarization power spectra in the CMB using InPCA, t-SNE, and the diffusion map. Axes reflect the true aspect ratio from extracted components in all cases. Here the model manifold is colored by the primordial fluctuation amplitude, the most prominent feature in CMB data. (A) InPCA extracts, as the first and second component, this amplitude term as well as the Hubble constant. These parameters control the two most dominant features in the Planck data and so reflect a physically meaningful hierarchy of importance. In contrast, (B) t-SNE extracts only the amplitude term and (C) the diffusion map extracts the amplitude term and a different parameter, the scalar spectral index η, which reflects the scale variance of the density fluctuations in the early universe. In all plots, the orange point represents our universe, as represented by Planck 2015 data.

Such stark differences between manifold learning methods are surprising, as all techniques aim to extract important features in the data distribution, i.e., important geometric features in the manifolds. Given the ranges of sampled parameters, one would expect the variation in the Hubble constant to relate in some way to one of the dominant components, as it does for InPCA. Figures illustrating the effect of different parameters are provided in SI Appendix, following results from ref. 27.

There are two important differences between InPCA and other methods. First, InPCA has no tunable parameters and yields a geometric object defined entirely by the model distribution. For example, t-SNE embeddings rely on parameters such as the perplexity, a learning rate, and a random seed (yielding nondeterministic results), and the diffusion maps rely on a diffusion parameter and choice of diffusion operator, all of which must be manually optimized to obtain good results. Second, t-SNE and diffusion maps embed manifolds in Euclidean spaces in a way which aims to preserve local features. However, InPCA seeks to preserve both global and local features, by embedding manifolds in a non-Euclidean space.

Summary

In this article, we introduce an unsupervised manifold learning technique, InPCA, which captures low-dimensional features of general, probabilistic models with wide-ranging applicability. We consider replicas of a probabilistic system to tune its dimensionality and consider the limit of zero replicas, deriving an intensive embedding that ameliorates the canonical orthogonality problem. Our intensive embedding provides a natural, meaningful way to characterize a symmetric distance between probabilistic data and produces a simple, deterministic algorithm to visualize the resulting manifold.