Information geometry for multiparameter models: new perspectives on the origin of simplicity

To judge whether different parameter values are similar, or points on the model manifold are nearby, we need some notion of distance. The idea of information geometry [3–6] is that this must come from the model itself, and the noise in our measurements: if two different values predict data which is easily distinguishable, then we should not regard them as being near to each other. Conversely, if two different parameter values are close together, they must predict nearly indistinguishable data.

For nonlinear least-squares models, the distance between the predictions of two parameters θ and $\tilde\theta$ is given by the χ² distance of equation (3), $\sum_{i = 1}^M (y_{\theta}(t_i) - y_{\tilde\theta}(t_i))^2/2\sigma_i^2$, which is basically the Euclidean distance between the two vectors y_θ and $y_{\tilde\theta}$. For more general models, distinguishability is measured between the probability distributions $p(x\vert\theta)$ predicted for a set of possible outcomes x (spin configurations for the Ising model [40], temperature maps for the cosmic background radiation [41], etc). We shall discuss three such measures, valid for general probability distributions, in section 7: the Hellinger distance [42] (equation (17)), the Bhattacharyya divergence [43] (equation (19)), and the Jeffreys [44, p 158] or symmetrized Kullback–Leibler [45, equation 2.5] divergence (equation (21)).

All of these measures are zero when $\theta = \tilde\theta$, are symmetric, and can be calculated for any two points θ and $\tilde\theta$; they are thus a global measure of distance. All of these measures reduce to the χ² Euclidean distance for the vector predictions of least-squares models. The latter two measures do not obey the triangle inequality, but each corresponds to distances in a Minkowski-like space where some directions are 'time-like' (section 7). These three, together with a wide family of other divergences [46, 47], share the same local measure, that is, the distance over an infinitesimal change $\tilde\theta = \theta + d\theta$:

$$\begin{align} ds^2 &= \sum_{\mu,\nu = 1}^D g_{\mu\nu} d\theta^\mu d\theta^\nu,\nonumber\\ g_{\mu\nu}(\theta) &= \sum_{x} \frac{\partial \log p(x \vert \theta)}{\partial \theta^\mu} \frac{\partial \log p(x \vert \theta)}{\partial \theta^\nu} p{({x} \vert {\theta})}. \end{align} \tag{ 5 }$$

This is the Fisher metric, or $g_{\mu\nu}$ for $\mu,\nu = 1,2,\ldots,D$ the Fisher information matrix (FIM) [48]. It may be thought of as measuring distance in parameter space in units of standard deviations, using the width of the distribution $p(x\vert\theta)$ on the space of possible data. (Indeed, the famous Cramér-Rao bound [49–51] uses the Fisher information metric to bound the ability to distinguish the predictions between two nearby models—the ability to know that the parameters $\tilde \theta$ are incorrect given data generated by parameters θ.) We shall study these more general model manifolds for probabilistic models by generating isometric embeddings. By preserving one of the distance measures $d^2(\theta,\tilde\theta)$ between every pair of models θ and $\tilde\theta$, we generate a geometrical representation of the model manifold that both preserves the FIM and (in the cases studied so far) exhibits the same sloppy hyperribbon structure as we have found in least-squares models.

While individual components of the FIM tell us how sensitive the model is to particular parameters, we are more interested in its sensitivity to parameter combinations. In the extreme case, many different parameters may control the same mechanism, making the model sensitive to changes in each of them in isolation. But if they all have equivalent effects, then the model is effectively one-dimensional, and it will be possible to make correlated changes which cancel, and do not affect the results. In this case, we would expect the FIM to have one large eigenvalue, and many small or zero ones. The one important eigenvector we will call the stiff direction in parameter space, and the remaining unimportant ones sloppy. The stiff direction is both important for predictions, and tightly constrained by observed data; the sloppy directions are neither. Figure 2 draws this in the context of an embedding $\mathcal{Y} \subset \mathbb{R}^3$.

Zoom In Zoom Out Reset image size

In more realistic models with many parameters, the eigenvalues of the FIM often have a characteristic spectrum: they are spread over many orders of magnitude, and are roughly evenly spaced on a log scale. Figure 3 shows some examples of this, from explicit computation of the FIM for models from a wide variety of sciences. It is a striking empirical fact that so many of them display this structure. These eigenvalues tell us about the local geometry of the model manifold, that is, the effect of small changes in the parameters. While these are extremely suggestive, they have the flaw that they depend on the particular parameterization chosen. If we re-scale some parameters, such as by changing their units, then the FIM and hence the spectrum of eigenvalues will change ⁸ .

Zoom In Zoom Out Reset image size

We can avoid this by measuring the global geometry of the model manifold [63, 64] using geodesics along the thick and thin directions. Let $\theta(\tau)$ be a path in parameter space, which then is mapped by the model onto a path on the model manifold. The length of this path, as defined by the Fisher information metric:

$$\begin{align} L = \int \sqrt{ds^2} = \int_0^1 d\tau \sqrt{\sum_{\mu\nu} g_{\mu\nu}(\theta(\tau)) \frac{\partial\theta^\mu}{\partial\tau} \frac{\partial\theta^\nu}{\partial\tau}}, \end{align} \tag{ 6 }$$

is independent of the way we define the coordinate system L on the model manifold. A geodesic is a curve connecting two points which minimizes the length between them. For models with Gaussian noise $x = y + \mathcal{N}(0, \sigma^2)$ as in equation (4) above, this geodesic distance agrees with the Euclidean path length in the prediction space, measured by walking along $\mathcal Y$.

We measure the widths of the model manifold by starting from an initial point (say, the best fit to the experimental data) and launching geodesics along one of the eigendirections in parameter space, shooting both ways until it hits the boundaries. Figure 4 shows the resulting geodesic widths across the model manifold for a number of different models. As advertised, they too show a roughly geometric progression over many orders of magnitude, with the stiff parameter directions yielding the widest directions and with the many sloppy directions corresponding to incredibly thin directions on the model manifold. We will call this emergent hierarchical structure of the model manifold a hyperribbon, after the idea that a ribbon is much longer than it is wide, and much wider than it is thick.

Zoom In Zoom Out Reset image size

We call a model sloppy, and the model manifold a hyperribbon, if the parameter-space eigenvalues and prediction-space widths are hierarchical: roughly equally spaced in log, spanning several decades. A sloppy model will have an effective low-dimensional description if many of its widths are small compared to the precision of the predictions demanded. This need not be the case—the ΛCDM model of the Universe cosmic microwave background radiation (column 3 in figure 3 and column 6 in figure 4) is sloppy (especially considering how few parameters it has) because the experiments are so precise that all parameters are well determined from the data. An impressive use of optimal experimental design [65] has shown that the individual parameters in our cell signaling network [1, 2] (column 1 in figure 3) could in principle be determined, albeit with a rather large number of experiments [66]. Alternatively, one could measure all the individual parameters separately in a model whose collective behavior is sloppy [67], but the resulting collective behavior will usually be predictable only if every one of the parameters is well specified [68]. Finally, most scientific models of complex systems like biological signaling are likely approximations to the microscopic behavior (Michaelis–Menten equations replacing enzyme reaction kinetics), and the resulting parameters are 'renormalized' to incorporate coarse-grained ignored complexity. Systematically designing experiments to measure these coarse-grained model parameters may just uncover flaws in the microscopic description that are irrelevant to the emergent behavior. We shall systematically study these useful coarse-grained models in section 4.1.

We can understand sloppiness, and give bounds for the widths of the model manifold that explain the observed hierarchies, through approximation theory. Specifically, we study polynomial approximations of model predictions. For simplicity, here we consider models whose predictions are a one-dimensional function $y_\theta(t)$ over the interval $t \in (-T,T)$, such as in figure 1 where we choose T = 0.9. For models where this is not the case, we can simply rescale t so that it fits in this region. The bounds on model predictions shall rely on the smoothness of y as the experimental control t is varied.

In the simplest case, assume y_θ is analytic in t, with Taylor series $\sum a_k t^{\,k}$ and coefficients $|a_k| \leqslant 1$ for all k. This is a sufficient condition to give the series radius of convergence equal to one, larger than the radius T of the time range being studied. Consider approximating y(t) by truncating its Taylor series. The approximating polynomial, $f_{N-1}(t) = a_0 + a_1 t + \cdots + a_{N-1} t^{\,N-1}$, differs from y(t) by at most:

$$\begin{align} \epsilon_N = \sum_{N}^\infty |a_k| t^{\,k} \leqslant \sum_{N}^\infty T^{\,k} \leqslant T^N/(1-T), \end{align} \tag{ 7 }$$

since $|a_k|\leqslant 1$ and $|t|\lt T$. Recall from equation (1) that the model manifold $\mathcal{Y}$ is defined as the space of all possible model predictions y_θ for varying model parameters θ. Understanding the constraints on $f_{N-1}$ for each y_θ will therefore allow us to understand what constraints are set on the model manifold $\mathcal{Y}$. Due to the construction of $f_{N-1}$, we have the following linear equation:

$$\begin{align*} \begin{bmatrix} f_{N-1}(t_0) \\ f_{N-1}(t_1) \\ \ldots \\ f_{N-1}(t_{N-1}) \end{bmatrix} &= \begin{bmatrix} a_0 + a_1 t_0 + \cdots + a_{N-1} t_0^{\,N-1}\\ a_0 + a_1 t_1 + \cdots + a_{N-1} t_1^{\,N-1}\\ \ldots \\ a_0 + a_1 t_{N-1} + \cdots + a_{N-1} t_{N-1}^{\,N-1} \end{bmatrix}\nonumber\\& = \begin{bmatrix} 1 & t_0 & \cdots & t_0^{\,N-1}\\ 1 & t_1 & \cdots & t_1^{\,N-1}\\ \cdots & \cdots & \cdots & \cdots\\ 1 & t_{N-1} & \cdots & t_{N-1}^{\,N-1} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \cdots\\ a_{N-1} \end{bmatrix}, \end{align*}$$

or, in particular:

$$\begin{align} f_{N-1}(\mathbf{t}) = V(\mathbf{t}) \mathbf{a}, \end{align} \tag{ 8 }$$

where V is the well-studied Vandermonde matrix $V_{jk} = (t_j)^k$. The best-known property of the Vandermonde matrix is its determinant $\det(V) = \prod_{0\leqslant i \leqslant j \lt N} (t_i-t_j)$. Notice that, since $t\in (-T,T)$, there are ${\sim} N^2/2$ factors $|t_i - t_j|$ ranging in size from 2 T down to $~1/N$: the determinant (the product of the eigenvalues) will be very small.

Because the model predictions are characterized by the coefficients in their Taylor series as well as where they are evaluated (i.e. the a_i and $V(\mathbf{t})$ in equation (8) respectively), understanding what constrains both a_i and V(t) will provide an understanding of what bounds the model.

For instance, for our models where the a_i are limited to be between −1 and 1, the Vandermonde matrix takes the hypercube of possible values of a and maps it to a hyper-parallelepiped $\mathcal{F}$ in prediction space that is incredibly skewed, as shown in figure 5; its volume is proportional to the incredibly small determinant. The model manifold $\mathcal{Y}$ is constrained to be within $\epsilon_N = T^N/(1-T)$ (equation (7)) of this hyper-parallelepiped because each $y_\theta\in\mathcal{Y}$ has a truncated power series $f_{N-1}$ which differs from it by at most ε_N . The highly skewed nature of this hyper-parallelepiped, whose widths are determined by the Vandermonde matrix, forces a similar hierarchical structure on the model manifold $\mathcal{Y}$.

Zoom In Zoom Out Reset image size

To expand the understanding of why a hyperribbon structure is observed for non-analytic models, in [30] it is proven for a relaxed smoothness condition where the coefficients a_i are bounded not by a hypercube but by an n-sphere of radius r, and not for all coefficients but only those up to a certain order. Bounds were derived both using a Chebyshev approximation and by using the Taylor series considered above. Following a very similar line of reasoning, it is possible to show that the model manifold $\mathcal{Y}$ is constrained by a hyper-ellipsoid whose $j^\mathrm{th}$ largest principal axis length is bounded by a power law ρ^−j . To illustrate this, consider three nonlinear models: (a) a sum of exponentials, such as in figure 1, (b) reaction velocities of an enzyme-catalyzed reaction, and (c) the infected population in an SIR epidemiology model. In all three cases, the models have varying model parameters which produce different predictions, shown in figure 6(a). The model manifolds are all bound by the same hyperellipsoid, shown in figure 6(b), and form a hyperribbon structure as the successive widths of the hyperellipsoid follow a geometric decay.

Zoom In Zoom Out Reset image size

Hyperribbons have also been observed for models with multiple control variables, for instance those which vary for both time and temperature [30, supplemental material]. Here, a two-dimensional Vandermonde matrix is constructed to understand the bounds in prediction space, producing a similar decay in the widths of the hyperellipsoid that bounds the model manifold with additional degeneracies based on the number of additional control variables.

We can get a different intuitive picture for why smoothness in the control parameters implies a hyperribbon structure in the predictions, by considering interpolation theory [63, 64]. Consider studying successive cross sections of the model manifold generated by experimental measurements at times $t_0,t_1,\dots$, here for analytic functions whose radius of convergence $|y_\theta^{(n)}(t)|/n! \lt1$ for all points t in the interval [53, 63]. The value of $y_\theta(t)$ at each time $t\in (-T,T)$ is an experimental prediction, and therefore a coordinate axis of our model manifold. Because there are an infinite number of points on the interval, and each point lends itself to a dimension on the model manifold, the resulting manifold is infinite-dimensional. Because of the bounded derivatives ⁹ , measuring $y_\theta(t)$ at a new time further constrains the range of possible values at other times, as shown in figure 7. Each experimental measurement $y_\theta(t_j)$ fixes a point, and has the effect of slicing the model manifold along axis t_j at y_j —keeping only the parameters θ whose predictions agree with the experiments.

Figure 7 shows a function $y_\theta(t)\in \mathcal{Y}$ taking values $y_\theta(t_j) = y_j$, and rigorous bounds on any function $z(t)\in \mathcal{Y}$ that agree with $y_\theta(t)$ at the constrained time points. The range of predictions $z(t^{^{\prime}})$ in figure 7 gives the width of the model manifold cross section after slicing along certain fixed points. Note that this width shrinks with each new experimental measurement, as the effect of additional constraints only further limits the range of allowed functions. Fixing the value of $y_\theta(t_j) = y_j$ shrinks the range of all other possible predictions, just as a ribbon snipped in two has a cut edge that typically reflects the width and thickness of the original ribbon. The convergence of a polynomial interpolation between the points $y(t_j)$ is equivalent to a hyperribbon structure of the model manifold $\mathcal{Y}$.

Zoom In Zoom Out Reset image size

Simplified models are found on the boundaries of the model manifold, and these boundaries exhibit a universal structure. Returning to the motivating example of two exponentials introduced in equation (2) and figure 1, observe that the boundary has a triangle-like structure: a two-dimensional manifold with three one-dimensional edges that meet at three zero-dimensional cusps. Although the size and shape of the model manifold depend on the observation time, this triangle-like boundary structure is independent of these details and depends only on the form of the observation function, $y_\theta(t)$. We have noted that each of these boundary elements is a simplified, approximate model. Formally, this boundary structure is a hierarchical cell complex that is a natural topological decomposition of the boundary.

Similar boundary complexes are found in models from very diverse contexts [78, 79]. That is to say, a typical N-dimensional model is bounded by a collection of N − 1-dimensional faces, which join at N − 2-dimensional edges, and so forth. The adjacency relationship among these facets of different dimensions defines a partially ordered set (often called a POSet), described as graded because each facet has a dimension associated with it. These relations can be conveniently visualized as a directed graph known as a Hasse diagram, as we have done in figure 8 for cGh theories and repeated in figure 9 for the exponential example.

Zoom In Zoom Out Reset image size

Each of the boundary cells of a model manifold corresponds to a simplifying approximation that corresponds to an extreme, idealized behavior of the model. The Hasse diagram summarizes how these approximations can be composed to yield progressively simpler models. That is to say, the Hasse diagram of the boundary complex is a road map from a complicated multi-parameter model through various approximations to the distinct behavioral regimes the model enables. When the model manifold is endowed with an information-theoretic metric, it, along with its boundary structure, is compressed in sloppy directions and stretched in stiff ones. Useful approximations, i.e. those that retain key explanatory variables, correspond to those boundary cells aligned with the long-axes of the manifold.

We translate this observation into an algorithm, the Manifold Boundary Approximation Method (MBAM) that uses information geometry to identify simplified models as boundary cells of a model manifold that is appropriate for a particular data set. The basic idea is to construct a geodesic shortest-path through the model's parameter space until it intersects a face on the boundary complex. By choosing a sloppy direction in parameter space (searching both directions along an eigenvector of the FIM with small eigenvalue), the geodesic will find a nearby boundary point with nearly identical predictions. By considering parameter values along the geodesic path, one deduces the limiting approximation that corresponds to the appropriate boundary cell. Often, this limiting approximation may involve one or more parameters going to extreme values, such as infinite or zero, in a coordinated way. We re-parameterize the model so that participating parameters are grouped together into a single combination that goes to zero at the boundary of the manifold. Finally, we analytically apply the approximation to the model to construct a simpler model with one fewer parameter. We iterate this process until the model is sufficiently simple.

MBAM has a formal similarity to using scale separation to perform model reduction. Scale separation is one of the oldest and most well-established methods of model reduction. It introduces a small parameter in which to perturb; the model reduces to a simplified form in the limit that the parameter becomes zero. Similarly, MBAM determines a combination of parameters, the perturbation parameter, that goes to zero at the boundary of the model manifold. However, sloppiness adds a new wrinkle to how we use and interpret this procedure. We do not necessarily require the true value of the perturbation parameter to be small. Rather, the true parameter value may be large provided the data to be explained do not constrain its value. Indeed, we observed above that although a model with separated scales will be sloppy, separated-scales are not necessary for sloppiness. The full model is statistically indistinguishable from the simpler separation-of-scales model. We advocate for using simple models because they reflect the informativity of the data, not because the true system to be modeled has well-separated scales.

Because the problems related to model complexity transcend scientific disciplines, there are many well-established reduction techniques that have grown-up around specific application areas. Most of these exploit the specific structure of the problem at hand, usually in conjunction with expert intuition. For example, dynamical systems are often simplified using singular perturbation theory in which the fast dynamics (often identified by an expert) are replaced by the slow manifold. Since, MBAM is agnostic to the mathematical form of the model, it can be compared to domain-specific reduction methods. In many cases, MBAM is formally equivalent to these techniques. However, because it is a data-driven method, it need not require expert insight and can guide insights into the modeling of complex systems.

For example, consider the problem of modeling transient dynamics of a power system. Power systems are large, engineered systems with many interconnected components. Time scales vary from speed-of-light electrical responses to load variations over hours and days. Models incorporate multi-physics in the form of heterogeneous components that describe not just power delivery, but generation from multiple sources, including renewables, and interactions with economic demand.

Of particular interest are models of synchronous generators, a power systems component that plays an important role in large simulations. The power community has developed a hierarchy of models describing synchronous generators that are formally connected through a series of singular perturbations. A typical inference problem involves estimating machine parameters from short disturbances that may last a few seconds. This problem is challenging because the relevant time scales are intermediate between the fastest and slowest time scales of the system, so it is difficult to judge which time scales are relevant. Indeed, the true dynamics are not necessarily faster or slower than the observation time scales, but parameters are unidentifiable because of sloppiness.

Applying MBAM to slowly-developing disturbances leads to the expected sequence of singularly perturbed reduced models advocated by theorists. However, when applied to disturbances at intermediate time scales, the results are more interesting. Recall that singularly perturbed systems have both an inner and an outer solution corresponding to the fast and slow dynamic respectively. Rather than uniformly removing fast dynamics, MBAM selects a combination of inner and outer solutions of the singularly perturbed system to remove parameters not constrained by the data. The remaining parameters are those that could be learned from the data.

Next, we turn to the theory of linear control that has produced many formal insights into the approximation of dynamical systems. We first discuss how control-theoretic quantities relate to sloppiness, and then demonstrate how MBAM relates to another standard approximation method known as balanced truncation.

A linear, time invariant (LTI) system takes the form:

$$\begin{align} \dot{x} & = Ax + Bu\nonumber\\ y & = Cx + Du, \end{align} \tag{ 9 }$$

where $x(t) \in \mathbb{R}^n$ are state space dynamical variables, $u(t) \in \mathbb{R}^p$ are inputs/control variables, and $y(t) \in \mathbb{R}^m$ are observed variables. A, B, C, and D are constant matrices of the appropriate shape. This model arises as the linearization of a generic dynamical system near a stable equilibrium, and theorems proved on LTI systems generalize to nonlinear dynamical systems under certain regularity conditions. Two important concepts to emerge in the study of LTI systems are observability and controllability. Roughly, observability is how well the internal state of the system can be inferred from its outputs while controllability is how easily a system can be driven to a particular state through a control protocol. The concepts of observability and controllability generalize to nonlinear systems, though they are most easily quantified in the context of LTI systems.

The model reduction problem as formulated in controls is not immediately analogous to the parameter reduction problem addressed by MBAM. For controls, 'model reduction' refers to reducing the dynamic order of a system. The goal is to approximate a system in the form of equation (9), with another LTI system with a smaller state space (i.e. $\tilde{x}(t) \in \mathbb{R}^{n^{^{\prime}}}$ where $n^{^{\prime}} \lt n$).

Balanced truncation uses the concepts of observability and controllability to identify accurate approximations. For LTI systems, observability and controllability are each quantified by a Gramian matrix, W_O and W_C respectively. The Gramians are symmetric, positive definite matrices that act as metrics quantifying the observability and controllability of vectors in the state space. A fundamental theorem for LTI systems shows that for an appropriate change of basis in the state space, the resulting observability and controllability Gramians are equal and diagonal. In this so-called 'balanced representation,' the vectors in the state space are equally controllable and observable. Balanced truncation is a reduction strategy in which an LTI system is first transformed into its balanced representation and the least observable/controllable states are truncated (removed) from the system.

To relate state-space reduction methods of control theory to MBAM, we parameterize the family of LTI systems of fixed dynamical order n. We generated data by choosing many potential inputs u(t) and evaluating the resulting time series y(t). The resulting model manifold is sloppy, and the sloppy parameters are those components in the A, B, and C matrices that determine the dynamics of the least observable/controllable states. When we apply MBAM to this system, we find the manifold boundary corresponding to the limit that the unobservable and uncontrollable states are constant decouple from the observed variables, y(t). This limit is equivalent to the balanced truncation method advocated by the controls community.

Balanced truncation has two major limitations. First, it is difficult to interpret because physically meaningful states are mixed together in the balanced representation. It is also difficult to generalize to nonlinear systems. However, the information geometry techniques underlying sloppy model analysis naturally accommodate these regimes. Thus, MBAM is a natural generalization of principles of balanced truncation to structure-preserving and nonlinear model reduction of interest to the control community.

As presented here, MBAM has some limitations and raises some fundamental questions about the geometric and topological properties of multi-parameter models. Here we discuss the relationship between MBAM and more general statistical theory. We also highlight some interesting formal considerations related to the parameter singularities that are beyond the scope of this review.

In statistics, simplified models that are special cases of a more general model are called nested models. Nested models are important in many statistical applications, especially hypothesis testing. The boundary models that we describe here are a specific type of nested models. But more generally, sometimes useful reduced models may also be found at non extreme values. For example, in a spin system with parameters quantifying couplings (e.g. nearest neighbor, next nearest neighbor, etc) that can be either positive or negative, a natural reduced model is to fix these couplings to zero rather than to extreme values. Future work could adapt the MBAM to identify other types of useful reduced models that are not on the boundary. Whether there are geometric ways to identify these other nested models that do not make explicit reference to their having one less parameter remains an open question.

Taking parameters to extreme values is a common reduction strategy in science because the reduced model can often be evaluated in a different way. For example, in a singularly perturbed dynamical system, an ordinary differential equation transitions into a differential-algebraic equation [80]. In general, the singular nature of these limits makes them especially relevant for inter-theoretic reasoning in mathematical physics as in figure 8. More broadly, however, model manifolds may have singular models that are not strictly boundaries, such as a the tip of a cone. These singularities often correspond to useful reduced models. More work is needed to understand the geometric and topological properties of these generalized singularities and develop reduction algorithms to efficiently find them.

The natural measure of how much information we learn from an experiment is the average reduction in the entropy over parameters θ upon seeing data x, which is the mutual information [83]:

$$\begin{align*} I(X;\Theta) = S(\Theta)-S(\Theta|X) = \sum_x \int\negthickspace d\theta\, p(x|\theta)\,p(\theta)\log\frac{p(x|\theta)}{p(x)}. \end{align*}$$

This information depends on the prior. Maximizing it may be thought of as choosing a minimally informative prior, and this choice is one way to formalize what the notion of 'uninformative' means ¹¹ . The idea of choosing a prior by maximizing this information was studied by Bernando, but always in the limit of $M\to\infty$, infinitely many repetitions or infinitely much data, where the resulting priors are called reference priors [84]. In this limit they approach Jeffreys prior (under mild conditions), and this idea may be seen as an alternative motivation for using Jeffreys prior [85]. But away from this limit, the resulting priors behave very differently, without losing reparameterization invariance. And, in particular, they have model selection properties which have been overlooked in the statistics literature.

Priors which maximize mutual information were investigated in [36], which showed that they have a preference for edges, and hence tend to select simpler subspaces. In fact they are almost always discrete, with weight on K atoms (points in parameter space) [86–89]:

$$\begin{align} p_{\star}(\theta) = \mathop{\mathrm{argmax}}_{p(\theta)}I(X;\Theta) = \sum_{a = 1}^{K}\lambda_{a}\delta(\theta-\theta_{a}). \end{align} \tag{ 11 }$$

Along a short dimension (say ${\lesssim}1$ standard deviations) the atoms will lie only at the boundaries. But along a long dimension (${\gg}1$ standard deviations) there will be many atoms, approaching a continuum. This turns out to be precisely the behavior we want on a hyperribbon.

To illustrate this behavior, consider the sum-of-exponentials model equation (2), used in figure 1, except with more parameters (and normalized to one at t = 0):

$$\begin{align} &x_{t}\sim y_{t}(\theta)+\mathcal{N}(0,\sigma^{2}),\qquad y_{t}(\theta) = \sum_{\mu = 1}^{D}\frac{1}{D}e^{-tk_{\mu}(\theta)},\nonumber\\&\quad 0\leqslant k_{\mu}(\theta)\leqslant\infty. \end{align} \tag{ 12 }$$

Figure 10 shows $p_{\star}(\theta)$ for a model with D = 2 parameters at various values of σ, drawn in y_t space for times $t = 1,3$. (This choice of co-ordinates gives an isometric embedding: Euclidean distance on the paper is proportional to the Fisher metric.) At very large σ there are just two atoms at the extreme values (i.e the zero-dimensional boundaries). At very small σ, the data-rich regime, there are many points almost tiling the space. And in the limit σ → 0 (equivalent to $M\to\infty$ repetitions) this approaches Jeffreys prior, which is flat in y. But it is the intermediate values which represent the more typical hyper-ribbon situation, in which some dimensions are data-rich, and some are data-starved.

Zoom In Zoom Out Reset image size

The idea of a prior which depends on the number of repetitions (or equivalently on the noise level σ) breaks the dictum 'yesterday's posterior is tomorrow's prior' for Bayesian updating. However, any prior which depends on the experiment to be performed (including Jeffreys) already breaks this rule. If tomorrow's experiment is upgraded to have less noise in some components, even if it otherwise explores the same model, then the combination requires a new prior, thus the update from today to tomorrow is different from a simple application of Bayes rule. For Jeffreys prior there is an exception for repetitions of the same experiment.

But, as we demonstrate below, this invariance to repetition means that Jeffreys performs very badly for hyper-ribbon model manifolds, sloppy models, in which we cannot fix all parameters [37]. Indeed, the distinction between a parameter we can fix and one we cannot is precisely that we do not yet have enough data to fix the latter.

The idea of accepting a discrete prior is also initially surprising to classical intuition. What if the truth were between two points of weight, would we eventually have poor accuracy? Three replies are possible. First, 'eventually' here implies taking much more data, and further updating the posterior, which is ruled out by the prior's dependence on the number of repetitions: the appropriate prior is the one for all the data. Second, along a relevant parameter direction, discretization is arguably no worse than truncating some measured θ to a sensible number of digits, for given error bars ¹² . And third, along an irrelevant direction, placing all weight at an extreme value is precisely what effective theories do. For example, QED assumes that the coupling constants of all forces except electromagnetism are exactly zero, which we know to be false in reality. But this false assumption is irrelevant to the experiments being described, and allows predictions of very high accuracy [91].

Jeffreys prior is, as we mentioned, the limit $M\to\infty$ of the maximum-information prior $p_{\star}(\theta)$. So according to our arguments, it is the correct thing to use in the regime when all parameters are well-constrained, when all are in the data-rich regime. How badly wrong is it when we are not in this regime? The eigenvalues of the metric tensor in sloppy, multiparameter models span an enormous range, so the square root of the determinant giving Jeffreys prior (equation (10)) is often enormously small. Worse, it varies often by many orders of magnitude even for small changes in parameters. Uninformative priors are supposed to be broadly agnostic about parameter choices, and Jeffreys prior fails in this regard.

Thanks to progress in numerical techniques for finding $p_{\star}(\theta)$ [37, 92] we can now illustrate this extreme weighting. Figure 11 shows the same sum-of-exponentials model, but now comparing a model with D = 3 parameters to one with D = 26 parameters. Both describe the same data, measuring decays at M = 26 different times in $1\leqslant t_i \leqslant5$, each with noise σ = 0.1. The discrete priors for these two models are fairly similar, but the concentration of weight in Jeffreys prior is much stronger in the higher-dimensional model. This is most clearly seen in the posterior, after observing a particular x as shown. In D = 3 Jeffreys performs well, but in D = 26 its posterior is centered 20 standard deviations away from the data. This is not a subtle error.

How many repetitions would be required in order for Jeffreys prior to be justified here? The smallest eigenvalue of the FIM is about 10⁻⁵⁰ here. To make this order 1, we need more repetitions than there are seconds in the age of the Universe [37]. While this is obviously a toy model, tens of irrelevant parameters, and astronomical timescales to measure them all, are also found in more realistic models (figures 3 and 4).

Another way to compare these priors is by how much information we will learn. The scores (expressed in bits) for the priors in figure 11 are shown in figure 12. Notice that adding additional parameters almost always increases $p_{\star}(\theta)$, but we observe that the extreme concentration of weight gives Jeffreys prior a much lower score in the high-dimensional model. Even without being able to find $p_{\star}(\theta)$ exactly, this measure offers a tool for model selection. But unlike traditional approaches to model selection, this takes place before seeing data x. It depends only on the resolution of the experiment, not its outcome.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finding $p_{\star}(\theta)$ is not always easy, as $I(X;\Theta)$ involves both an integral over X and a sum over the atoms: The K points of weight all interact with each other, via the data they generate. We present here another prior which shares many of the virtues (and avoids Jeffreys' vices), while being easier to draw samples from.

To derive this, consider first the following distribution on data space X, called a normalized maximum likelihood ¹³ :

$$\begin{align} p_{\mathrm{NML}}(x) = \frac{\max_{\hat{\theta}}p(x|\hat{\theta})}{Z},\qquad Z = \int dx\:\max_{\hat{\theta}}p(x|\hat{\theta}). \end{align} \tag{ 13 }$$

For models with Gaussian noise, this describes a cloud of points filling and surrounding the embedded model manifold $\mathcal{Y} \subset \mathbb{R}^M$ (figures 10 and 13), with a range given by the measurement error in the predictions. We then define a prior distribution on parameter space Θ, by simply projecting the weight of $p_{\mathrm{NML}}(x)$ at each x to the corresponding maximum likelihood point $\hat{\theta}$. We may write this as:

$$\begin{align} p_\sharp(\theta) = \int dx\;p_{\mathrm{NML}}(x)\:\delta\big(\theta-\mathop{\mathrm{argmax}}_{\hat{\theta}}p(x|\hat{\theta})\big). \end{align} \tag{ 14 }$$

This prior has full support, but places extra weight on points on the boundaries of the model manifold, and more weight on lower-dimensional boundaries (such as fold lines, and corners), thus generalizing the idea of a slab-and-spike prior [95] to many dimensions and curved manifolds. How much weight is on boundaries is chosen automatically, depending on the level of noise in the likelihood function.

In practice, this prior can be computed efficiently by sampling $p_{\mathrm{NML}}(x)$ using Monte Carlo methods, and then recording for each accepted point x the corresponding $\hat{\theta}$ to give a sampling distribution $p_\sharp(\theta)$ in parameter space. The Metropolis weight does not require the normalization factor Z.

Figures 10 and 13 show the result, in the same sum of exponentials models as above. For all models $x\sim y(\theta)+\mathcal{N}(0,\sigma^{2})$, the embedding of the model manifold is smeared out by a distance σ to make $p_{\mathrm{NML}}(x)$. This means that, unlike Jeffreys prior, the difference between places where the manifold is thin compared to σ, and places where it is very very thin, will not matter at all. The smeared weight is pulled back to contribute to $p_\sharp(\theta)$ exactly at the closest edge. And lower-dimensional edges have more weight; zero-dimensional points usually have finite weight. Thus $p_\sharp(\theta)$ is as a generalized slab-and-spike prior, with weight both in the volume and on edges.

Zoom In Zoom Out Reset image size

Like $p_{\star}(\theta)$, the degree of preference for edges depends on the level of noise σ, or the number of repetitions. Along irrelevant, sloppy, short, directions, almost all the weight is on the edges. Along relevant, stiff, long, directions, most of it fills the interior. As the noise is decreased, dimensions will smoothly pass from the former category into the latter. We view this as a virtue, but it does mean the prior $p_\sharp(\theta)$ is not invariant to repetition of the experiment.

While it does not maximize the information $I(X;\Theta)$, the prior $p_\sharp(\theta)$ usually does very well on this measure. Figure 12 above shows an estimate of this, compared to $p_{\star}(\theta)$ as well as Jeffreys prior, and a log-normal prior.

Because many models are microscopically defined, but are mostly observed through their emergent properties we sought to understand how the model manifold of microscopic parameters contracts when observables are coarse-grained. Rather than coarse-graining the parameters to directly build an effective description as is done in a usual RG treatment, we instead wanted to understand how coarse-graining changes one's ability to measure microscopic parameters [38, 39]. To do this we consider a continuous coarse-graining whereas the smallest length scale $l = l_0\exp{b}$ changes the parameters change according to $d\theta^\mu / db = \beta^{\mu}$ where $\vec{\beta}$ are the beta functions that define the flow. Infinitesimal distances after such a flow can be calculated by flowing the parameters under the usual RG, and then measuring their distance in the coarse grained model, where the system size is smaller. This change in the microscopic parameter space metric can be summarized by the following equation [39]:

$$\begin{align} \frac{d g_{\mu \nu}}{db} = \underbrace{\beta^{\alpha} \partial_\alpha g_{\mu \nu}}_{g^{^{\prime}}~\textrm{along flow}} + \underbrace{g_{\alpha \mu} \partial_{\nu} \beta^{\alpha} + g_{\alpha \nu} \partial_{\mu} \beta^{\alpha}}_{\textrm{contraction of}~\Delta \theta}- \underbrace{D g_{\mu \nu}}_\textrm{Vol. contraction}. \end{align} \tag{ 15 }$$

The first two terms constitute the Lie derivative of the metric under flow β. The third term arises because the metric is extensive, proportional to the system size L^D , where D is the dimensionality of the system.

While this equation summarizes the local change in the metric, and as such the change in infinitesimal distances during coarse-graining, it remains to analyze how these distances change under typical RG flows. For this analysis it is useful to use the fact that the metric can be written as the second derivative of the free energy (in units where β = 1), $g_{\mu\nu} = \frac{\partial}{\partial \theta^\mu} \frac{\partial}{\partial \theta^\nu} F$. In renormalization group analysis it is useful to separate the free energy into a singular and an analytic piece, $F = F_a+F_s$, where the singular part is non-analytic as a function of its variables near the RG fixed point [97]. Near the RG fixed point the metric, like the free energy, can be also divided into singular and analytic pieces, $g = g^s+g^a$. By assumption, the singular contributions of the free energy are those that are preserved under RG flow—the three terms in equation (15) must precisely cancel. This means that along relevant directions—which are best measured through their influence on emergent properties, the microscopic metric does not contract ¹⁴ . The preservation of the metric under flow can be used to constrain the form of the metric along relevant directions [98]. However, these constraints only apply to the lower dimensional subset of parameters which are relevant in the usual RG treatment.

Irrelevant directions are instead measured through their influence on microscopic details and their metric is dominated by analytic contributions. For these directions the first term is expected to be small and infinitesimal distances ds contract according to:

$$\begin{align} \frac{d\log (ds^2)}{db} = -D+2y_i, \end{align} \tag{ 16 }$$

where y_i is the corresponding (negative) RG exponent.

We have verified this prediction numerically for a two dimensional Ising model, where we have added couplings between non-nearest-neighbor spins (see figure 14 and [38]). We coarse-grain the system using simple decimation, only observing the red sites on a checkerboard, and so on recursively. In the Ising model these results have an intuitive interpretation; the magnetic field and temperature can be measured easily without seeing individual spins as they control the net magnetization and the correlation length. However, other parameter combinations can only be measured by examining microscopic configurations.

Zoom In Zoom Out Reset image size

While this metric centered approach provides information about how the local structure of the manifold contracts under coarsening, it does not necessarily provide information about the global structure of the manifold. While we have fewer analytical approaches, we can numerically investigate how an entire model manifold contracts under coarse-graining. We can numerically investigate the diffusion equation, where a simple renormalization group treatment yields a single relevant parameter (the drift, v) a single marginal parameter (the diffusion constant D) and an infinite number of irrelevant parameters which control higher cumulants of the distribution [38, 39]. Numerically we investigate a hopping model, where parameters θ^µ describe the probability that a single particle moves a distance µ in a discrete time step (see figure 15(A)). We consider a seven parameter model when particles are allowed to move up to three spaces to the left and right; particle number is not necessarily conserved, but there is a six dimensional subspace which conserves particle number (the subspace where $\sum_\mu \theta^\mu = 1$).

Zoom In Zoom Out Reset image size

The metric can be measured at a single point in parameter space by observing the particle density after some number of time steps. Its eigenvalue distribution is shown in figure 15(A) for the parameters where $\theta^\mu = 1/7$. In the upper panel of figure 15(A), the particle distribution is shown after t time steps. Initially the distribution is given by the parameters, $\rho_i \propto \theta^i$, and each parameter can be measured independently, yielding a metric proportional to the identity, which is not sloppy. After sufficient coarse-graining steps, the distribution approaches a Gaussian as required by the central limit theorem, characterized by a velocity and diffusion constant, and, if particle number is not conserved, by a rate of particle creation/annihilation (R). In this limit, the spectrum of the Fisher Information becomes sloppy, with the stiffest eigenvector changing R, the next stiffest changing V, then D, and with a hierarchy of irrelevant parameters which manifest only through higher cumulants in the distribution.

The non-local contraction of distances can be seen qualitatively by looking at two (two-dimensional) diffusion kernels, one a triangle, one a square, whose drift is both zero and which have the same diffusion constant (red and blue points, respectively in figure 15(B)). After a single time step these two kernels yield very clearly distinguishable particle distributions. However, after multiple time steps the two distributions become similar, approaching the same Gaussian. We numerically investigated the seven parameter model, by taking kernels uniformly from the six parameter space of particle conserving kernels, and plotting them in the first three principle components of the space of particle densities. After a single time step the points fall in a uniform three dimensional blob—the corresponding model manifold does not have a hyper-ribbon structure. The red and blue points, which each have zero drift and the same diffusion constant, are far apart in this model manifold. However, as time evolves the three dimensional blob collapses onto a two dimensional structure parameterized by a drift and diffusion constant. As this limit is approached the red and blue points approach each other rapidly.

In renormalizable models, after coarse-graining the model manifold contracts so that every point is close to a smaller dimensional surface which can often be described by a reduced model with just relevant parameters kept. This gives an information theoretic perspective for one of the central practical findings of the renormalization group that a simple 'toy model' will often suffice to describe a microscopically complicated system. Many systems are not renormalizable; however, they still often show a sloppy spectrum in their Fisher Information's eigenvalues and in the hyper-ribbon structure of their model manifolds. These systems would often not be sloppy if they were observed at some suitable microscopic scale. However, for macroscopic behavior, only a subset of their microscopic properties remain important, leading to local sloppiness and a global hyper-ribbon structure. In that sense the sloppiness that arises in renormalizable systems may be qualitatively similar to that arising in other systems. One key difference is that renormalizable models have a self-similarity during coarse-graining, allowing metric flow to be written compactly in terms of parameter space point flow (as in equation (15) from [39], see [99] for an interesting perspective).

To understand why it is that embedding in a pseudo-Riemannian space is more natural for manifolds of probabilistic models, consider the failures of a natural embedding method, using the Hellinger distance:

$$\begin{align} d^2_\mathrm{Hel}[\theta,\tilde\theta] = 4 \sum_x \left(\sqrt{p(x|\theta)}-\sqrt{p(x,\tilde\theta)}\right)^2. \end{align} \tag{ 17 }$$

Since probabilities are normalized, $1 = \sum_x p(x|\theta) = d^2_\mathrm{Hel}[\theta,\theta]/4$, the vectors $\mathbf{n}_x = 2 \sqrt{p(x|\theta)}$ are points on a hypersphere of radius two in a space whose dimension is the cardinality of the observations. The Hellinger distance is the Euclidean distance between these points. As θ is varied, the points form a manifold:

$$\begin{align} y_x(\theta) = 2\sqrt{p(\mathbf{x}_i|\theta)}, \end{align} \tag{ 18 }$$

which can be visualized with PCA. The metric between nearby points on this manifold is given by the Fisher information ¹⁵ , it is isometric. The fact that it is isometric makes it a natural embedding, in the sense that distance between points represents the distinguishability of the different probability distributions (e.g. if two distributions have no overlap, then they are orthogonal from each other in this representation, since their dot product is zero).

The curse of dimensionality that hinders Hellinger can be seen easily in this framework (see [40]). As the dimension of the data x increases, the distributions become more distinguishable from each other. They therefore become increasingly orthogonal, and all distances between non-identical distributions converge to a constant value as they become maximally apart on the hypersphere. As an illustration of this, consider the finite, two-dimensional Ising model in a field with nearest-neighbor coupling (figure 16). For small system sizes, two similar field and coupling strengths yield similar spin configurations and can be difficult to distinguish from one another. However, as the system size increases (and therefore the dimension of the observed data also increases) the increase in information allows one to more easily evaluate the underlying parameters. The probability distributions representing different parameters are more distinguishable from each other, and are therefore further apart in the hypersphere embedding. In the limit of infinite data, we can completely disentangle all possible parameter combinations: in other words, all distributions become orthogonal and equidistant from one another, forming a hypertetrahedron in probability space. The geometry of a regular hypertetrahedron with N vertices clearly cannot be captured with fewer than a N − 1-dimensional embedding in a Euclidean space. But this does not remain true if we allow for both space and time-like coordinates.

If instead, we look at the information density rather than the total amount of information in a distribution, we can avoid this orthogonality catastrophe by reducing rather than increasing the dimensionality. To accomplish this, we can make use of replica theory. We can simulate the curse of dimensionality by considering replicas of the original system, i.e. multiple, independent samples of the original distributions. As an illustration, consider the two-dimensional Ising model with field and nearest neighbor interactions. A $2\times 2$ spin system, visualized in this hypersphere embedding, has 2⁴ dimensions, one for each possible spin configuration. By considering a range of possible field and coupling strengths, we sweep out a region of the model manifold. Figure 16(a) shows the first three components of a classic PCA of this embedding, with the curvature of the hypersphere very evident. As the system size increases to $4\times 4$ spins, the dimensionality increases and the orthogonality problem becomes more evident. This effect is simulated by looking at four replicas of the $2\times 2$ system, shown in figure 16(c).

Rather than increasing the system size with many replicas we decrease it by considering the limit of zero replicas. In this way, we obtain an intensive embedding. The distance between points becomes characterized not by the Hellinger distance but the Bhattacharyya ¹⁶ divergence [43]:

$$\begin{align} d^2_\mathrm{Bhat}[\theta,\tilde\theta] = -8\log\left(\sum_x \sqrt{p(x|\theta) p(x,\tilde\theta)}\right). \end{align} \tag{ 19 }$$

By working through the derivation of PCA using the Bhattacharyya divergence, which does not respect the triangle inequality, we obtain an embedding which is pseudo-Riemannian, called InPCA [41].

Empirical results from probabilistic models show that they share the hyperribbon structures we find in least squares models. In figure 4, for example, we see geometrically decreasing widths for the two-dimensional Ising model, variable-width Gaussians, CMB skymaps, and neural networks [40, 41]. Figure 18 shows the model manifold for a $2\times2$ Ising model, exhibiting projections along the first five InPCA components, clearly showing this hyperribbon decrease of widths.

The Hellinger distance and Bhattacharyya divergence belong to a larger class of divergences [46, 47] whose metric is given by the Fisher Information. The most commonly used distance measure for probability distributions, $p(x|\theta) \equiv p$ and $p(x|\tilde{\theta}) \equiv \tilde p$, is the Kullback–Leibler divergence [45] or relative entropy:

$$\begin{align} D_\mathrm{KL}[\theta, \tilde{\theta}] = \sum_x p(x|\theta)\log\left(\frac{p(x|\theta)}{p(x|\tilde{\theta})}\right). \end{align} \tag{ 20 }$$

Because it is asymmetric, often researchers symmetrize [44, 45] the Kullback–Leibler divergence by considering the sum of the two permutations, which we write as the square of a distance $d_\mathrm{sKL}$:

$$\begin{align} d^2_\mathrm{sKL} &= D_\mathrm{KL}\big[\theta,\tilde\theta\big] + D_\mathrm{KL}\big[\tilde{p} \big\Vert p\big] \nonumber\\ &= \sum_x \Big(p(x|\theta) - p(x|\tilde{\theta})\Big)\log\left(\frac{p(x|\theta)}{p(x|\tilde{\theta})}\right). \end{align} \tag{ 21 }$$

By combining the method of projecting data along orthogonal direction from MDS [41] (section 7.4) with the symmetric KL divergence, we can derive the isKL embedding [40]. An astonishing feature of this method is that it can map exponential families into a finite dimensional embedding space, with dimension twice the number of model parameters.

Exponential families [111] have probability distributions that can be expressed in the form familiar to physicists, as a partition function described by a Boltzmann-like distribution:

$$\begin{align} p(x|\theta) = h(x)\exp\left(\sum_i \eta_i(\theta) \phi_i(x) -A(\theta)\right), \end{align} \tag{ 22 }$$

where $\eta_i(\theta)$ is a function of model parameters and is known as the ith natural parameter, $\phi(x)$ is a function only of x and known as the ith sufficient statistic, h(x) is a normalizing function, and $A(\theta)$ is the log partition function. By combining equations (21) and (22), we see that:

$$\begin{align} d^2_\mathrm{sKL} = \sum_i \left( \eta_i(\theta) - \eta_i(\tilde{\theta})\right)\big( \left\langle\phi_i\right\rangle_\theta - \left\langle\phi_i\right\rangle_{\tilde{\theta}} \big), \end{align} \tag{ 23 }$$

where $\left\langle\cdot\right\rangle_\theta$ reflects taking an average over x using $p(x|\theta)$. Rearranging the terms in the above equation, we can see that:

$$\begin{align} d^2_\mathrm{sKL} & = \frac{1}{4}\sum_i\left( [\eta_i(\theta) + \left\langle\phi_i\right\rangle_\theta] - [\eta_i(\tilde{\theta}) + \left\langle\phi_i\right\rangle_{\tilde{\theta}}]\right)^2\nonumber\\ &\quad - \left( [\eta_i(\theta) - \left\langle\phi_i\right\rangle_\theta] - [\eta_i(\tilde{\theta}) - \left\langle\phi_i\right\rangle_{\tilde{\theta}}]\right)^2 \nonumber\\ & = [\mathcal{S}_i(\theta) - \mathcal{S}_i(\tilde{\theta})]^2 - [\mathcal{T}_i(\theta) - \mathcal{T}_i (\tilde{\theta})]^2, \end{align} \tag{ 24 }$$

where $\mathcal{S}_i$ and $\mathcal{T}_i$ represent space-like and time-like components, and can be directly calculated from the model:

$$\begin{align} \mathcal{S}_i(\theta) & = \frac{1}{2}\big(\eta_i(\theta) + \left\langle\phi_i\right\rangle_\theta\big),\nonumber \\ \mathcal{T}_i(\theta) & = \frac{1}{2}\big(\eta_i(\theta) - \left\langle\phi_i\right\rangle_\theta\big). \end{align} \tag{ 25 }$$

There are a finite number of components, determined by the number of parameters, regardless of the dimension of x. This provides a systematic, computationally efficient way of visualizing model manifolds, as opposed to having to calculate the distance between all points and computing the projections [40]. For example, calculating the entire probability distribution for the Ising model (needed for the Hellinger and Bhattacharyya embeddings) becomes infeasible for all but tiny system sizes (figures 16–18), but traditional Monte-Carlo methods can be used to calculate the magnetization and interaction energy needed for the isKL embedding of the Ising model for any reasonable system size (figure 19).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Moreover, isKL provides a valuable check for data without a known model. By calculating the $d^2_\mathrm{sKL}$ between data, and visualizing both the real and imaginary components extracted using MDS, we can check to see if the resulting embedding is finite or infinite. If finite, then one is led to believe that the data results from a process described by an exponential family (basically a statistical sampling described by a free energy).

Here we discuss the methods used to generate the visualizations produced in this section. Classical PCA and MDS [41, 109] take data sets and infer features which are linearly uncorrelated. PCA allows researchers to examine data in a new basis by centering it and then extracting components of maximal variance, or diagonalize its 'moment of inertia', then extract projections in a hierarchical manner which minimizes the error in distances between points when restricted to a smaller dimensional subspace. Because PCA optimizes distances, one can extend it to take as inputs not the raw data but rather the distances themselves, and in this way obtain MDS. This approach can be generalized to intensive embeddings (InPCA and isKL) in the following way. First, given a set of N probability distributions x_i , calculate the squared distance, i.e. the Bhattacharyya divergence or Symmetrized Kullback–Leibler divergence between them, $X_{i,j} = D(\mathbf{x}_i \mid \mid \mathbf{x}_j)$ where D represents the choice of divergence. Next, generate the intensive cross-covariance matrix [41], $W = - PXP/2$, where $P_{ij} = \delta_{ij} - 1/N$ serves to mean-shift the data, so that it is centered on its center of mass.

Next, compute the eigenvalue decomposition for W, $W = U\Sigma U^T$, and note that there may be both positive and negative eigenvalues. It is important that they be ranked by magnitude, i.e. $|\Sigma_{00}| \geqslant |\Sigma_{11}| \geqslant \dots$ . The coordinate projections for the intensive embedding are determined by:

$$\begin{align} T = U\sqrt{\Sigma}. \end{align} \tag{ 26 }$$

One then plots the projections using the columns of T. Coordinates corresponding to positive eigenvalues represent space-like components, in that they have a positive-squared contribution to the distance between points. Conversely, coordinates corresponding to negative eigenvalues represent time-like components, in that they have negative-squared contributions to the distance between points. The resulting figures are interpreted much like a spacetime diagram, where points with zero distance are light-like separated. Note that, because all points are a positive-squared distance apart (since their divergence is non-negative) their overall distance will be a space-like separation, although they may be visualized using both space-like and time-like components.

The sloppiness of multiparameter models suggests that extracting parameters from data will be challenging or impossible. Model reduction methods such as MBAM show that one can describe data in complex systems with emergent models for which many parameters and model components are simply ignored (set to zero, infinity, or replaced by other parameters), even when those components of the model are microscopically known to be important and relevant. While we may find the emergent models more comprehensible, are not they wrong? The rates in a system biology model do not go to infinity, the enzymes are not in the limit of short supply—the simpler model works by shifting other parameters to compensate for these simplifying assumptions.

The progress of science often exhibits this in reverse—one develops a simple theory which captures the essential features of a problem, and then adds complexity to address new experimental facts. Can we use information geometry to guide the design of new experiments to pin down models in an optimal way? Optimal experimental design [113] is a well-developed branch of statistics, but usually is done in a linearized approximation (but see [114])—can we use hyperribbon embeddings and differential geometry to turn experimental design into a more global tool? As a particular example, the InPCA methods were developed in order to visualize the space of possible Universes as reflected by current measurements of the microwave background radiation. The next generation of instruments is designed to measure and test new physical phenomena—gravitational waves from the epoch of inflation, evidence of light relic particles in the early Universe, and signatures of dark matter in the cosmic microwave background radiation. Would it be useful to visualize the predictions of extensions of the ΛCDM model that include these phenomena, by mapping out those portions of the model manifold that are consistent with current experimental data? Can we use such visualizations to help guide the choice of the next generation of instruments?

Our information geometry gives us a theoretical justification for trusting the predictions of nonlinear models with many undetermined parameters. Climate science, for example, incorporates not only well-controlled theories of heat, hydrodynamics, and radiation, but much more speculative models of clouds, glaciers, and forest ecologies (not to mention human behavior). It is reassuring to hear that equally complex models in systems biology with even more undetermined parameters are quite useful in guiding drug design (see below). But just because complex models fitted to data can make correct predictions does not mean they must do so. Statistics provides many tests of model validity and many methods to compare the predictive power of rival models, for cases where the parameters are well determined. Can we use information geometry methods to create similar tests and comparison tools for multiparameter models where the uncertainties in model parameters are typically not only outside the linear regime, but often unbounded?

Emergence in physics relies on small parameters. Continuum limits work (elasticity, fluids, superfluids, etc) when frequencies and wavelengths are low; the renormalization group works (emergent scale invariance, avalanches, crackling noise, etc) when in addition there is a proximity to a continuous transition between states. Sections 6 and 7.3 discussed how information geometry applied to these models form hyperribbons. On the other hand, we used interpolation theory to explain the low-dimensional emergent behavior and reduced models that arise in nonlinear least-squares models fit to data (section 3). Emergence in physics and emergence in least-squares models are both characterized by sloppiness (only a few combinations of microscopic parameters dominate the behavior), and both form hyperribbons, but our explanations for the two seem completely different. For more general probabilistic models, intensive embeddings (section 7) usually lead to hyperribbons and emergent theories, but we do not yet have an explanation for this in general, except for the special (but very useful) case of exponential families, where the model manifold not only is thinner and thinner along successive cross sections, but becomes zero thickness after a finite number of embedding dimensions.

Can we find a common theme between these disparate results—a unified theory for the emergence of simplicity from complexity? One possible clue is that sloppiness is a property both of the scientific model and the predictions of interest. The Ising model is not sloppy, and has no beautiful emergent theory, unless one only cares about long length and time scales. And a systems biology model is not sloppy if one designs experiments to measure one parameter at a time. Multiparameter least squares models have emergent simplicity when one focuses on collective behavior of the system as a whole, and statistical mechanics has emergent simplicity when one focuses on collective behavior of large numbers of particles. Perhaps it is caring about the answers to the right questions that makes science possible?

In our work visualizing the model manifold for probabilistic models (section 7), we grappled with the curse of dimensionality, commonplace in big data applications. We identified the problem of having too much information: data points in too high a dimensional space naturally becoming orthogonal. We cured the problem first using replica theory and then using an analytic method for models that form 'exponential families'. Our cures work because of two features. (a) We found a way of embedding our model in a Minkowski-like space, using MDS ([41], section 7.4), that violates the triangle inequality. (b) We use measures of squared distance between data points (Bhattacharyya and symmetrized Kullback–Leibler) that diverge logarithmically with the overlap between the probability distributions. (This ensures that doubling the amount of data simply rescales all distances by a factor of 2).

Can we apply these basic ideas to other big data problems, that are not framed as probabilistic models? We know of at least one machine learning method, PHATE [115] (due to Kevin Moon) that has applied these two principles to visualize biological data (in particular single-cell RNA sequencing data sets). Here let us consider two possible applications to machine learning, as represented by the millions-of-parameter deep learning neural network models.

Deep neural networks are 'big' in the number of neurons per layer (hundreds to thousands), the number of neural connections between layers (millions of weight parameters), the number of dimensions in the initial data used to train the network (say, one per pixel in the images being recognized), and the number of data points used to train the network (millions of images). They have been remarkably successful at extracting meaningful information from these large data sets, suggesting that there is some hidden, low-dimensional structures [116] in the data, that is disentangled and flattened via the nonlinear transformations through succeeding neural layers.

With Pratik Chaudhari, some of us have been exploring [117] how one might use intensive embeddings to visualize both the training of these neural networks (as in [41]) and potentially explore low-dimensional structures discovered after the training process. Consider a network that categorizes, taking N images x and returns a preference probability $\boldsymbol{\rho}(\mathbf{x})$ whose largest value is the predicted category. The vector $\mathbf{R} = (1/N) \{\boldsymbol{\rho}(x_1), \boldsymbol{\rho}(x_2), \dots\}$ can be viewed as a probability distribution over the whole set, and one could use InPCA or isKL to visualize how the network trains with time, averaging over the random initial states and training steps. Visualization of the resulting structures for training and test data could lend geometrical insight into how and why different training protocols differ in their generalization errors. Alternatively, one could examine the neural activations $h(\mathbf{x})$ of the different layers of a single neural network—each encoding distilled information about the initial data x gleaned from the previous layer that is most relevant to the task at hand. Using PCA on the vectors h(x) is likely doomed for the standard reasons, but conceivably using MDS with an intensive metric could allow one to see the emergence of a low dimensional (perhaps hyperribbon-like) representation of the data as one passes through successive trained layers.

The manifold boundary approximation method is a principled way to select effective theories for sloppy systems. MBAM is objective, data-driven, and does not require expert intuition. Through sequences of simplifying approximations, it gives identifiable models that retain essential elements, such as dominant feedback loops, while remaining expressed in terms of microscopic physics. We anticipate that simplified models will be increasingly important by enabling clearer reasoning about big data and complex processes.

Recent result from systems biology modeling the Wnt pathway illustrate how this may be done. Wnt signaling is the primary mechanism by which short-range extra-cellular signals are relayed to the nucleus, and it plays a major role in normal embryonic development as well as cancer tumerogenesis. While the real system is notoriously complex, it regularly exhibits simple, comprehensible behavior. For example, the steady state conditions of the pathway may encode for the level of extra-cellular Wnt [118] while at other times, a negative feedback drives sustained oscillations that regulate somitogenesis [119]. In each case, the behavior is explicable in terms of simple models with only a few basic mechanisms. However, by reasoning about the relationships between these reduced models within the larger space of candidate approximations, more is possible. The supremal model, i.e. simplest model reducible to each effective behavior includes the relevant mechanisms for transitioning between states [120]. This model makes predictions that are qualitatively different from the data to which it was trained and is a starting point for controller design.

Can simple, effective models facilitate the transition of big data into practical knowledge in spite of being incomplete? Quantitative systems pharmacology uses complex models to design and engineer drugs [121], but reasoning about causality within complex models is difficult. In our conversations with the FDA and pharmaceutical companies, challenges in justifying new drugs and human trials are closely related to model complexity [122]. If it is unreasonable to identify all of the parameters in realistic biological models, how can drug developers meet standards of government approval agencies [123]? Simple models that are rigorously derived from first principles can be more easily and cheaply validated. Furthermore, simple models that contain only the relevant mechanism are more readily communicated and trusted in ways that complicated simulations are not.

Can similar reasoning be applied beyond biology? Empirical potentials in molecular modeling are designed to mimic quantum effects within a classical framework. But which combination of quantum effects are responsible for specific macroscopic material properties? The combinatorial large space of stoichiometry and material properties make the problem particularly challenging. Large databases of quantum calculations [124, 125] archive data while projects such as OpenKIM [126] satisfy the dual problem of documenting and archiving useful models. But leveraging these data and models to design new materials remains an open problem. Can simple models that are successful in a limited context help clarify what must be accurately modeled to guide material by design?

We have speculated about the relationship between sloppiness and large machine learning models. Although model size seems to play an important role in learning capacity and training rate, large models cannot be deployed on small devices without a reduction step. Furthermore, large black-box models also do not lead to explicit, declarative knowledge in the same way as minimal, parsimonious representation [127]. The nascent fields of explainable AI [128] and physics-guided machine learning [129–131] may also benefit from simplification in a context-adaptive way.

We have also mentioned climate as a field whose predictions are the result of multi-physics simulations. Uncertainties in specific components (e.g. models of clouds or forest ecology), as well as their role in the larger simulation, make it difficult to reason about causality in such simulations. Simplified models that automatically identify the essential feedback loops and effective control knobs may guide how modelers understand how these uncertain components affect conclusions. They may also aid science communication and education for policy makers and the general public.

The situation in climate is reminiscent of that in power systems that similarly use multi-physics to model large engineered systems. When complex models fail to match data (even presumably identified models), for example in well-instrumented blackouts [132], it is bewildering to identify the source of the problem. Simplified models make it easier to identify and correct errors in the model. Our model reduction methods have made inroads towards this end in power systems [32]; can they similarly find use in these other areas?

Power systems application also face the challenge that they operate at the intersection of industry and government where issues of individual privacy, proprietary information, and public security are relevant. Model reduction is a means to limit and regulate how information is shared or redacted among vested parties. Models that contain sensitive infrastructure information can be systematically and objectively coarse-grained to remove classified pieces while simultaneously meeting the demands of public disclosure and oversight.

The computational cost of solving the geodesic equation limit the scope of current algorithms to models with at most a few hundred parameters, so alternative strategies are needed. It may be possible to guide optimization algorithms toward the boundary; merging the fitting and reduction steps. Alternative approaches may focus on targeted subsets of parameters. The topological relationship among potential models remains largely unexplored, but such a theory could better inform model selection for target applications.

Efficiently reasoning with big models and big data will require new organizational paradigms [133, 134]. The physical theories of the twentieth century, interconnected through a web of classical limits as in figure 8, suggest a way forward. Effective theories organize our understanding and focus our attention to the relevant physical principles within specific regimes of validity. The hyperribbons presented here reflect a similar information hierarchy, but sloppiness adds a new wrinkle to the traditional narrative by justifying simple theories without resorting to classical limits, coarse-graining, or renormalization. Visualization tools and model reduction methods extract simple explanations in a context-specific way, and our theories of information topology organize them into a network of interrelated models, each applicable at an appropriate scale and for specific applications. While science has traditionally progressed by discovering effective theories before the underlying mechanism, the future of science may well operate in the opposite direction, cataloging and organizing effective theories for the Universe's myriad complexities. Sloppy models abetted with the formalism of information geometry explain why and how this can be done.