Hierarchical models in the brain

This paper describes a general model that subsumes many parametric models for continuous data. The model comprises hidden layers of state-space or dynamic causal models, arranged so that the output of one provides input to another. The ensuing hierarchy furnishes a model for many types of data, of arbitrary complexity. Special cases range from the general linear model for static data to generalised convolution models, with system noise, for nonlinear time-series analysis. Crucially, all of these models can be inverted using exactly the same scheme, namely, dynamic expectation maximization. This means that a single model and optimisation scheme can be used to invert a wide range of models. We present the model and a brief review of its inversion to disclose the relationships among, apparently, diverse generative models of empirical data. We then show that this inversion can be formulated as a simple neural network and may provide a useful metaphor for inference and learning in the brain.

This paper describes hierarchical dynamic models (HDMs) and reviews a generic variational scheme for their inversion. We then show that the brain has evolved the necessary anatomical and physiological equipment to implement this inversion, given sensory data. These models are general in the sense that they subsume simpler variants, such as those used in independent component analysis, through to generalised nonlinear convolution models. The generality of HDMs renders the inversion scheme a useful framework that covers procedures ranging from variance component estimation, in classical linear observation models, to blind deconvolution, using exactly the same formalism and operational equations. Critically, the nature of the inversion lends itself to a relatively simple neural network implementation that shares many formal similarities with real cortical hierarchies in the brain.

Recently, we introduced a variational scheme for model inversion (i.e., inference on models and their parameters given data) that considers hidden states in generalised coordinates of motion. This enabled us to derive estimation procedures that go beyond conventional approaches to time-series analysis, like Kalman or particle filtering. We have described two versions; variational filtering [1] and dynamic expectation maximisation (DEM; [2]) that use free and fixed-form approximations to the posterior or conditional density respectively. In these papers, we used hierarchical dynamic models to illustrate how the schemes worked in practice. In this paper, we focus on the model per se and the relationships among its special cases. We will use DEM to show how their inversion relates to conventional treatments of these special cases.

A key aspect of DEM is that it was developed with neuronal implementation in mind. This constraint can be viewed as formulating a neuronally inspired estimation and inference framework or conversely, as providing heuristics that may inform our understanding of neuronal processing. The basic ideas have already been described, in the context of static models, in a series of papers [3]–[5] that entertain the notion that the brain may use empirical Bayes for inference about its sensory input, given the hierarchical organisation of cortical systems. In this paper, we generalise this idea to cover hierarchical dynamical systems and consider how neural networks could be configured to invert HDMs and deconvolve sensory causes from sensory input.

This paper comprises five sections. In the first, we introduce hierarchical dynamic models. These cover many observation or generative models encountered in the estimation and inference literature. An important aspect of these models is their formulation in generalised coordinates of motion; this lends them a hierarchal form in both structure and dynamics. These hierarchies induce empirical priors that provide structural and dynamic constraints, which can be exploited during inversion. In the second and third sections, we consider model inversion in general terms and then specifically, using dynamic expectation maximisation (DEM). This reprises the material in Friston et al. [2] with a special focus on HDMs. DEM is effectively a variational or ensemble learning scheme that optimises the conditional density on model states (D-step), parameters (E-step) and hyperparameters (M-step). It can also be regarded as a generalisation of expectation maximisation (EM), which entails the introduction of a deconvolution or D-step to estimate time-dependent states. In the fourth section, we review a series of HDMs that correspond to established models used for estimation, system identification and learning. Their inversion is illustrated with worked-examples using DEM. In the final section, we revisit the DEM steps and show how they can be formulated as a simple gradient ascent using neural networks and consider how evoked brain responses might be understood in terms of inference under hierarchical dynamic models of sensory input.

To simplify notation we will use fx: = fx(x) = ∂xf = ∂f/∂x to denote the partial derivative of the function, f, with respect to the variable x. We also use x ˙ = ∂tx for temporal derivatives. Furthermore, we will be dealing with variables in generalised coordinates of motion, which will be denoted by a tilde; x̃: = [x,x′,x″,…]T = [x [0],x [1],x [2],…]T, where x [i] denotes ith order motion. A point in generalised coordinates can be regarded as encoding the instantaneous trajectory of a variable, in the sense it prescribes its location, velocity, acceleration etc.

In this section, we cover hierarchal models for dynamic systems. We start with the basic model and how generalised motion furnishes empirical priors on the dynamics of the model's hidden states. We then consider hierarchical forms and see how these induce empirical priors in a structural sense. We will try to relate these perspectives to established treatments of empirical priors in static and state-space models.

Dynamic causal models are probabilistic generative models p(y,ϑ) based on state-space models. As such, they entail the likelihood, p(y|ϑ) of getting some data, y, given some parameters ϑ = {x,v,θ,λ} and priors on those parameters, p(ϑ). We will see that the parameters subsume different quantities, some of which change with time and some which do not. These models are causal in a control-theory sense because they are state-space models, formulated in continuous time.

A dynamic input-state-output model can be written as(1)The continuous nonlinear functions f and g of the states are parameterised by θ. The states v(t) can be deterministic, stochastic, or both. They are variously referred to as inputs, sources or causes. The states x(t) meditate the influence of the input on the output and endow the system with memory. They are often referred to as hidden states because they are seldom observed directly. We assume the stochastic terms (i.e., observation noise) z(t) are analytic, such that the covariance of z̃ = [z,z′,z″,…]T is well defined; similarly for the system or state noise, w(t), which represents random fluctuations on the motion of the hidden states. Under local linearity assumptions (i.e., ignoring high-order derivatives of the generative model functions), the generalised output or response ỹ = [y,y′,y″,…]T obtains from recursive differentiation with respect to time using the chain rule(2)Note that the derivatives are evaluated at each point in time and the linear approximation is local to the current state. The first (observer) equation show that the generalised states u = [ṽ,x̃,]T are needed to generate a generalised response that encodes a path or trajectory. The second (state) equations enforce a coupling between neighbouring orders of motion of the hidden states and confer memory on the system.

At this point, readers familiar with standard state-space models may be wondering where all the extra equations in Equation 2 come from and, in particular, what the generalised motions; w′, w″, … represent. These terms always exist but are ignored in standard treatments based on the theory of Markovian processes [6]. This is because standard Markovian (c.f., Wiener) processes have generalised motion that has infinite variance and are infinitely ‘jagged’ or rough. This means w′, w″, … and x″, x‴, … have no precision (inverse variance) and can be ignored with impunity. It is important to realise that this approximation is not appropriate for real or actual fluctuations, as noted at the inception of the standard theory; “a certain care must be taken in replacing an actual process by Markov process, since Markov processes have many special features, and, in particular, differ from the processes encountered in radio engineering by their lack of smoothness… any random process actually encountered in radio engineering is analytic, and all its derivative are finite with probability one” ([6], pp 122–124). So why have standard state-space models, and their attending inversion schemes like Kalman filtering, dominated the literature over the past half-century? Partly because it is convenient to ignore generalised motion and partly because they furnish reasonable approximations to fluctuations over time-scales that exceed the correlation time of the random processes: “Thus the results obtained by applying the techniques of Markov process theory are valuable only to the extent to which they characterise just these ‘large-scale’ fluctuations” ([6], p 123). However, standard models fail at short time-scales. This is especially relevant in this paper because the brain has to model continuous sensory signals on a fast time-scale.

Having said this, it is possible to convert the generalised state-space model in Equation 2 into a standard form by expressing the components of generalised motion in terms of a standard [uncorrelated] Markovian process, ς(t):(3)The first line encodes the autocorrelation function or spectral density of the fluctuations w(t) in term of smoothness parameters, γ 1,…γn, where n is the order of generalised motion. These parameters can be regarded as the coefficients of a polynomial expansion Pn(∂t)w = ς(t) (see [6], Equation 4.288 and below). The second line obtains by substituting Equation 2 into the first and prescribes a standard state-space model, whose states cover generalised motion; x [0],…,x [n]. When n = 0 we recover the state equation in Equation 1, namely, x ˙ = f(x,v)+ς. This corresponds to the standard Markovian approximation because the random fluctuations are uncorrelated and w = ς; from Equation 3. When n = 1⇒w+γ 1 w′ = ς, the fluctuations w(t) correspond to an exponentially correlated process, with a decay time of γ 1 ([6], p 121). However, generally n = ∞: “Therefore we cannot describe an actual process within the framework of Markov process theory, and the more accurately we wish to approximate such a process by a Markov process, the more components the latter must have.” ([6], p 165). See also [7] (pp 122–125) for a related treatment.

If there is a formal equivalence between standard and generalised state-space models, why not use the standard formulation, with a suitably high-order approximation? The answer is that we do not need to; by retaining an explicit formulation in generalised coordinates we can devise a simple inversion scheme (Equation 23) that outperforms standard Markovian techniques like Kalman filtering. This simplicity is important because we want to understand how the brain inverts dynamic models. This requires a relatively simple neuronal implementation that could have emerged through natural selection. From now on, we will reserve ‘state-space models’ (SSM) for standard n = 0 models that discount generalised motion and, implicitly, serial correlations among the random terms. This means we can treat SSMs as special cases of generalised state-space models, in which the precision of generalised motion on the states noise is zero.

Given the form of generalised state-space models we now consider what they entail as probabilistic models of observed signals. We can write Equation 2 compactly as(4)Where the predicted response g̃ = [g,g′,g″,…]T and motion f̃ = [f,f′,f″,…]T in the absence of random fluctuations areand D is a block-matrix derivative operator, whose first leading-diagonal contains identity matrices. This operator simply shifts the vectors of generalised motion so x [i] that is replaced by x [i+1].

Gaussian assumptions about the fluctuations provide the likelihood, p(ỹ|x̃,ṽ). Similarly, Gaussian assumptions about state-noise furnish empirical priors, p(x̃|ṽ) in terms of predicted motion(5)We will assume Gaussian priors on the generalised causes, with mean and covariance . The density on the hidden states p(x̃|ṽ) is part of the prior on quantities needed to evaluate the likelihood of the response or output. This prior means that low-order motion constrains high-order motion (and vice versa). These constraints are discounted in standard state-space models because the precision on the generalised motion of a standard Markovian process is zero. This means the only constraint is mediated by the prior p(x ˙|x,v). However, it is clear from Equation 5 that high-order terms contribute. In this work, we exploit these constraints by adopting more plausible models of noise, which are encoded by their covariances and (or precisions and ). These are functions of unknown hyperparameters, λ which control the amplitude and smoothness of the random fluctuations.

Figure 1 (left) shows the directed graph depicting the conditional dependencies implied by this model. Next, we consider hierarchal models that provide another form of hierarchical constraint. It is useful to note that hierarchical models are special cases of Equation 1, in the sense that they are formed by introducing conditional independencies (i.e., removing edges in Bayesian dependency graphs).

HDMs have the following form, which generalises the (m = 1) model above(6)Again, f (i) : = f(x (i),v (i)) and g (i) : = g(x (i),v (i)) are continuous nonlinear functions of the states. The processes z (i) and w (i) are conditionally independent fluctuations that enter each level of the hierarchy. These play the role of observation error or noise at the first level and induce random fluctuations in the states at higher levels. The causes v = [v (1),…,v (m)]T link levels, whereas the hidden states x = [x (1),…,x (m)]T link dynamics over time. The corresponding directed graphical model is shown in Figure 1 (right). In hierarchical form, the output of one level acts as an input to the next. When the state-equations are linear, the hierarchy performs successive convolutions of the highest level input, with random fluctuations entering at each level. However, inputs from higher levels can also enter nonlinearly into the state equations and can be regarded as changing its control parameters to produce quite complicated generalised convolutions with ‘deep’ (i.e., hierarchical) structure.

The conditional independence of the fluctuations at different hierarchical levels means that the HDM has a Markov property over levels, which simplifies attending inference schemes. See [8] for a discussion of approximate Bayesian inference in conditionally independent hierarchical models of static data. Consider the empirical prior implied by Equation 6(7)where the full prior is now restricted to the last level. Equation 7 is similar in form to the prior in Equation 5 but now factorises over levels; where higher causes place empirical priors on the dynamics of the level below. The factorisation in Equation 7 is important because one can appeal to empirical Bayes to interpret the conditional dependences. In empirical Bayes [9], factorisations of the likelihood create empirical priors that share properties of both the likelihood and priors. For example, the prediction g̃ (i) = g̃(x̃ (i),ṽ (i)) plays the role of a prior expectation on ṽ (i−1), yet it has to be estimated in terms of x̃ (i),ṽ (i). In short, a hierarchical form endows models with the ability to construct their own priors. These formal or structural priors are central to many inference and estimation procedures, ranging from mixed-effects analyses in classical covariance component analysis to automatic relevance determination in machine learning. The hierarchical form and generalised motion in HDMs furnishes them with both structural and dynamic empirical priors respectively.

In generalised coordinates, the precision, is the Kronecker tensor product of a temporal precision matrix, S(γ) and the precision over random fluctuations, which has a block diagonal form in hierarchical models; similarly for . The temporal precision encodes temporal dependencies among the random fluctuations and can be expressed as a function of their autocorrelations(8)Here is the second derivative of the autocorrelation function evaluated at zero. This is a ubiquitous measure of roughness in the theory of stochastic processes [10]. Note that when the random fluctuations are uncorrelated, the curvature (and higher derivatives) of the autocorrelation are infinite. In this instance, the precision of high-order motion falls to zero. This is the limiting case assumed by state-space models; it corresponds to the assumption that incremental fluctuations are independent (c.f., a Wiener process or random walk). Although, this is a convenient assumption that is exploited in conventional Bayesian filtering schemes and appropriate for physical systems with Brownian processes, it is less plausible for biological and other systems, where random fluctuations are themselves generated by dynamical systems ([6], p 81).

S(γ) can be evaluated for any analytic autocorrelation function. For convenience, we assume that the temporal correlations have the same Gaussian form. This gives(9)Here, γ is the precision parameter of a Gaussian autocorrelation function. Typically, γ>1, which ensures the precisions of high-order motion converge quickly. This is important because it enables us to truncate the representation of an infinite number of generalised coordinates to a relatively small number; because high-order prediction errors have a vanishingly small precision. An order of n = 6 is sufficient in most cases [1]. A typical example is shown in Figure 2, in generalised coordinates and after projection onto the time-bins (using a Taylor expansion, whose coefficients comprise the matrix Ẽ). It can be seen that the precision falls quickly with order and, in this case, we can consider just six orders of motion, with no loss of precision.

When dealing with discrete time-series it is necessary to map the trajectory implicit in the generalised motion of the response onto discrete samples, [y(t 1),…,y(tN)]T = Ẽ ỹ(t) (note that this is not necessary with continuous data such as sensory data sampled by the brain). After this projection, the precision falls quickly over time-bins (Figure 2, right). This means samples in the remote past or future do not contribute to the likelihood and the inversion of discrete time-series data can proceed using local samples around the current time bin; i.e., it can operate ‘on-line’.

We can now write down the exact form of the generative model. For dynamic models, under Gaussian assumptions about the random terms, we have a simple quadratic form (ignoring constants)(10)The auxiliary variables comprise prediction errors for the generalised response and motion of hidden states, where and are the respective predictions, whose precision is encoded by . The use of prediction errors simplifies exposition and may be used in neurobiological implementations (i.e., encoded explicitly in the brain; see last section and [4]). For hierarchical models, the prediction error on the response is supplemented with prediction errors on the causes(11)Note that the data and priors enter the prediction error at the lowest and highest level respectively. At intermediate levels the prediction errors, v (i−1)−g (i) mediate empirical priors on the causes. In the next section, we will use a variational inversion of the HDM, which entails message passing between hierarchical levels. These messages are the prediction errors and their influence rests on the derivatives of the prediction error with respect to the unknown states(12)This form highlights the role of causes in linking successive hierarchical levels (the DT matrix) and the role of hidden states in linking successive temporal derivatives (the D matrix). The DT in the upper-left block reflects the fact that that the prediction error on the causes depends on causes at that level and the lower level being predicted; ε (i)v = v (i−1)−g(x (i),v (i)). The D in the lower-right block plays a homologous role, in that the prediction error on the motion of hidden states depends on motion at that order and the higher order; ε [i]x = x [i+1]−f(x [i],v [i]).

These constraints on the structural and dynamic form of the system are specified by the functions g = [g (1),…,g (m)]T and f = [f (1),…,f (m)]T, respectively. The partial derivatives of these functions have a block diagonal form, reflecting the model's hierarchical separability(13)Note that the partial derivatives of g(x,v) have an extra row to accommodate the top level. To complete model specification we need priors on the parameters and hyperparameters. We will assume these are Gaussian, where (ignoring constants)(14)

In this section, we have introduced hierarchical dynamic models in generalised coordinates of motion. These models are about as complicated as one could imagine; they comprise causes and hidden states, whose dynamics can be coupled with arbitrary (analytic) nonlinear functions. Furthermore, these states can have random fluctuations with unknown amplitude and arbitrary (analytic) autocorrelation functions. A key aspect of the model is its hierarchical form, which induces empirical priors on the causes. These recapitulate the constraints on hidden states, furnished by the hierarchy implicit in generalised motion. We now consider how these models are inverted.

This section considers variational inversion of models under mean-field and Laplace approximations, with a special focus on HDMs. This treatment provides a heuristic summary of the material in [2]. Variational Bayes is a generic approach to model inversion that approximates the conditional density p(ϑ|y,m) on some model parameters, ϑ, given a model m and data y. This is achieved by optimising the sufficient statistics (e.g., mean and variance) of an approximate conditional density q(ϑ)with respect to a lower bound on the evidence (marginal or integrated likelihood) p(y|m) of the model itself. These two quantities are used for inference on the parameters of any given model and on the model per se. [11]–[15]. The log-evidence for any parametric model can be expressed in terms of a free-energy F(ỹ,q) and a divergence term, for any density q(ϑ) on the unknown quantities(15)The free-energy comprises the internal energy, U(y,ϑ) = ln p(y,ϑ) expected under q(ϑ) and an entropy term, which is a measure of its uncertainty. In this paper, energies are the negative of the corresponding quantities in physics; this ensures the free-energy increases with log-evidence. Equation 15 indicates that F(ỹ,q) is a lower-bound on the log-evidence because the cross-entropy or divergence term is always positive.

The objective is to optimise q(ϑ) by maximising the free-energy and then use F≈ln p(ỹ|m) as a lower-bound approximation to the log-evidence for model comparison or averaging. Maximising the free-energy minimises the divergence, rendering q(ϑ)≈p(ϑ|y,m) an approximate posterior, which is exact for simple (e.g., linear) systems. This can then be used for inference on the parameters of the model selected.

Invoking an arbitrary density, q(ϑ) converts a difficult integration problem (inherent in computing the evidence; see discussion) into an easier optimisation problem. This rests on inducing a bound that can be optimised with respect to q(ϑ). To finesse optimisation, one usually assumes q(ϑ) factorises over a partition of the parameters(16)In statistical physics this is called a mean-field approximation. This factorisation means that one assumes the dependencies between different sorts of parameters can be ignored. It is a ubiquitous assumption in statistics and machine learning. Perhaps the most common example is a partition into parameters coupling causes to responses and hyperparameters controlling the amplitude or variance of random effects. This partition greatly simplifies the calculation of things like t-tests and implies that, having seen some data, knowing their variance does not tell you anything more about their mean. Under our hierarchical dynamic model we will appeal to separation of temporal scales and assume, q(ϑ) = q(u(t))q(θ)q(λ), where u = [ṽ,x̃,]T are generalised states. This means that, in addition to the partition into parameters and hyperparameters, we assume conditional independence between quantities that change (states) and quantities that do not (parameters and hyperparameters).

In this dynamic setting q(u(t)) and the free-energy become functionals of time. By analogy with Lagrangian mechanics, this calls on the notion of action. Action is the anti-derivative or path-integral of energy. We will denote the action associated with the free energy by F̅, such that ∂tF̅ = F. We now seek q(ϑi) that maximise the action. It is fairly easy to show [2] that the solution for the states is a function of their instantaneous energy, U(t): = U(u|θ,λ) = ln p(ỹ,u|θ,λ)(17)where V(t) = ∂tV̅ u is their variational energy. The variational energy of the states is simply their instantaneous energy averaged over their Markov blanket (i.e., averaged over the conditional density of the parameters and hyperparameters). Because the states are time-varying quantities, their conditional density is a function of time-dependent energy. In contrast, the conditional density of the parameters and hyperparameters are functions of their variational action, which are fixed for a given period of observation.(18)Where Uθ = ln p(θ) and Uλ = ln p(λ) are the prior energies of the parameters and hyperparameters respectively and play the role of integration constants in the corresponding variational actions; V̅ θ and V̅ λ.

These equations provide closed-form expressions for the conditional or variational density in terms of the internal energy defined by our model; Equation 10. They are intuitively sensible, because the conditional density of the states should reflect the instantaneous energy; Equation 17. Whereas the conditional density of the parameters can only be determined after all the data have been observed; Equation 18. In other words, the variational energy involves the prior energy and an integral of time-dependent energy. In the absence of data, when the integrals are zero, the conditional density reduces to the prior density.

If the analytic forms of Equations 17 and 18 were tractable (e.g., through the use of conjugate priors), q(ϑi) could be optimised directly by iterating these self-consistent nonlinear equations. This is known as variational Bayes; see [16] for an excellent treatment of static conjugate-exponential models. However, we will take a simpler approach that does not require bespoke update equations. This is based on a fixed-form approximation to the variational density.

Under the Laplace approximation, the marginals of the conditional density assume a Gaussian form q(ϑi) = N(ϑi: µi,Ci) with sufficient statistics µi and Ci, corresponding to the conditional mean and covariance of the ith marginal. For consistency, we will use µi for the conditional means or modes and ηi for prior means. Similarly, we will use Σi and Ci for the prior and conditional covariances and Πi and Pi for the corresponding inverses (i.e., precisions).