The curvature effect in Gaussian random fields

In science, many natural phenomena are modeled by a stochastic complex system in which collective properties emerge from the non-linear interactions among its parts in different levels of scale [1, 2]. Typically, in these systems, the whole is more than the sum of the parts. In this context, the characterization of the dynamics of stochastic systems using a mathematical tool is crucial for understanding the underlying processes that govern the emergence of complex behavior [3]. For example, being able to predict when a phase transition is going to happen is a relevant and challenging problem in complex systems analysis.

The dynamics of stochastic complex systems has recently garnered a lot of interest in the physics literature, directly contributing to the solution of several problems in social [4], biological [5] and economic [6] sciences. A relevant aspect concerns the prediction of phase transitions in a quantitative way by means of an objective mathematical criterion [7, 8]. In this paper, we compute intrinsic geometric properties from the underlying parametric space of random fields composed by Gaussian variables [9] and study how these quantities change along phase transitions.

Geometrodynamics is a research field whose main goal is to characterize and describe physical phenomena completely in terms of geometry [10], in an attempt to unify the fundamental forces and reformulate general relativity in terms of metric tensors of Riemannian manifolds [11, 12]. These issues have been investigated by several physicists and remain an active field in the 21st century, as a mathematical tool for the unification of gravitation [13], quantum mechanics [14] and in the study if particle systems [15, 16]. Besides, modern quantum physics studies argue that the most fundamental constituents of the universe are likely to be fields and not properly particles [17]. Information geometry is a research field that combines information theory and Riemannian geometry to study intrinsic geometric properties of parametric spaces associated with random variables [18–20]. Inspired by these ideas, our investigation proposes an experiment to measure how the variation of the curvature in the underlying manifold (parametric space) of random field models is related to variations in the system's entropy. In this context, our study can be thought as an attempt to apply information geometry to completely characterize the dynamics of Gaussian random fields in terms of intrinsic properties of their parametric spaces.

In particular, we assume some simplifying hypothesis: first, the random field is Markovian in the sense that the probability of a given variable in the field depends only on the variables belonging to a local neighborhood system around the variable [21]. Second, the model is isotropic in the sense that the inverse temperature parameter, which control the spatial dependence structure, is spatially invariant and constant for all orientations in space [22]. Last, but not least, we deal with a pairwise interaction model, which means that we allow only binary relationships. In summary, we consider a pairwise isotropic Gaussian–Markov random field (GMRF) to model the interaction between spatially dependent Gaussian random variables organized in a 2D lattice. The degree of interaction is quantified by a coupling parameter, also known as the inverse temperature. It has been observed that an increase in the inverse temperature parameter beyond a critical value induces the emergence of phase transitions in the system [23].

In the analysis of complex systems, one question that naturally arises is: how can we predict that the system is approaching a phase transition in a quantitative way? Besides, how to determine whether the system is moving towards order or randomness? Often, researchers extract several statistical and physical measures such as entropy, free energy, average magnetization, etc from both systems and use an extrinsic distance function (Euclidean distance) to measure a degree of similarity between the feature vectors. The proposed method allows the computation of intrinsic geometric properties from the parametric space, such as the first and second fundamental forms (Fisher information matrices) of the underlying manifold and its shape operator.

To summarize the big picture, we briefly describe the computational experiments conducted in this investigation: at the beginning of the dynamics, the inverse temperature parameter is zero, and the random field model degenerates to a regular Gaussian distribution (independent random variables). In this situation, the parametric space exhibit constant negative curvature (hyperbolic geometry) [24]. The idea is to analyze how the emergence of a spatial dependence structure along time leads to irreversible geometric transformations in the parametric space. The obtained results show the existence of the curvature effect, which can be described as: the variations of the curvature when the system moves towards higher entropy states is different from the variations of the curvature when the system moves towards lower entropy states, which induces the emergence of an intrinsic arrow of time as a natural one-way direction of evolution [25, 26]. The main objective of this scientific investigation is to propose an information-geometric framework to understand and characterize the dynamics of random fields defined on 2D lattices. Little is known about this process in situations where the random variables are not independent. The key idea is to analyze how the emergence of a spatial dependence structure along time leads to geometric transformations in the parametric space.

The remaining of the paper is organized as follows: in section 2, we define Gaussian random fields, showing that it is a curved statistical model. In section 3, we present important differential geometry concepts, such as surface, tangent planes, fundamental forms, Gaussian, mean and principal curvatures. Section 4 defines the first and second-order Fisher information matrices, their relation with the fundamental forms and presents the complete derivation of the components of these operators in pairwise isotropic GMRFs. Section 5 shows the results obtained by computational experiments using Markov chain Monte Carlo (MCMC) simulations. Finally, section 6 presents the conclusions and final remarks.

Random fields are mathematical models used in the study of stochastic and non-linear complex systems [27]. Among all the random field models categorized in the literature, Gaussian random fields are remarkably important [28]. Part of this relevance is due to the central limit theorem [29], which states that the summation of independent and identically distributed random variables tend to be normally distributed. Several random field applications can be found in the literature, from image processing techniques in computer vision [30] to statistical and cosmological physics [31, 32].

Pairwise isotropic GMRFs are mathematical structures particularly suitable to study spatially dependent continuous random variables by means of non-linear interactions between neighboring particles on a lattice. The main advantage of this model is related to the mathematical tractability. First, in this model, all the parameters are scalars, making the parametric space a regular manifold. Moreover, we can derive closed-form expressions for several expected values, which makes it possible the exact computation of information-theoretic measures, such as Fisher information. In other models, we have to approximate these quantities using MCMC simulations, drastically increasing the computational burden. Furthermore, by invoking the Hammersley–Clifford theorem [33], which states the equivalence between Gibbs random fields (global models) and Markov random fields (local models), we characterize a pairwise isotropic GMRF by the set of local conditional density functions:

$$\begin{equation}p\left({x}_{i}\vert {\eta }_{i},\vec{\theta }\right)=\frac{1}{\sqrt{2\pi {\sigma }^{2}}}\,\mathrm{exp}\left\{-\frac{1}{2{\sigma }^{2}}{\left[{x}_{i}-\mu -\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{j}-\mu \right)\right]}^{2}\right\}\end{equation} \tag{ 1 }$$

where η_i denotes the second-order neighborhood system comprised by the eight nearest neighbors of x_i , $\vec{\theta }=(\mu ,{\sigma }^{2},\beta )$ denotes the vector of model parameters, with μ and σ² being, the expected value (mean) and the variance of the random variables in the lattice, and β being the inverse temperature, which encodes the spatial dependence between the variables in the field. Figure 1 shows the first, second and third order neighborhood systems defined on a 2D lattice.

Zoom In Zoom Out Reset image size

Note that if β = 0, the model degenerates to a regular Gaussian distribution, used to model independent random variables. The main advantage in using the local model is that, by avoiding the joint Gibbs distribution, which has a parametrisation in terms of vectors and matrices, we considerably simplify the parametric space.

2.1. A curved statistical model

From mathematical statistics, it is possible to express a likelihood function $p(\mathbf{X}\vert \vec{\theta })$ on its natural parameters as:

$$\begin{equation}p(\mathbf{X}\vert \vec{\theta })=\mathrm{exp}\left\{\sum\limits _{j=1}^{K}{c}_{j}(\vec{\theta }){T}_{j}\left(\mathbf{X}\right)+d(\vec{\theta })+S\left(\mathbf{X}\right)\right\}\end{equation} \tag{ 2 }$$

where $S\left(\mathbf{X}\right)$ is a function of the observations only, $d(\vec{\theta })$ is a function of the parameters only, $\vec{T}=\left({T}_{1}\left(\mathbf{X}\right),{T}_{2}\left(\mathbf{X}\right),\dots ,{T}_{k}\left(\mathbf{X}\right)\right)$ denotes the vector of sufficient statistics and $\vec{c}=({c}_{1}(\vec{\theta }),{c}_{2}(\vec{\theta }),\dots ,{c}_{k}(\vec{\theta }))$ denotes the vector of natural parameters. Let X = {x₁, x₂, ..., x_n } be a sample of a pairwise isotropic GMRF model where Δ = 8 represents the number of neighbors. As the number of natural parameters is greater than the number of parameters, isotropic pairwise GMRF's are curved models. Then, the pseudo-likelihood function can be expressed as:

$$\begin{align}\hfill F\left(\mathbf{X}\vert {\eta }_{i},\vec{\theta }\right)& ={\left(2\pi {\sigma }^{2}\right)}^{-n/2}\,\mathrm{exp}\left\{-\frac{1}{2{\sigma }^{2}}\sum\limits _{i=1}^{n}{\left[\left({x}_{i}-\mu \right)-\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{j}-\mu \right)\right]}^{2}\right\}\hfill \\ \hfill & =\mathrm{exp}\left\{-\frac{n}{2}\left[\mathrm{log}(2\pi {\sigma }^{2})+\frac{{\mu }^{2}}{{\sigma }^{2}}\right]+\frac{\beta {\Delta}{\mu }^{2}n}{{\sigma }^{2}}\left[1-\frac{\beta {\Delta}}{2}\right]\right\}\hfill \\ \hfill & \quad \times \mathrm{exp}\left\{\left[\frac{\mu }{{\sigma }^{2}}\left(1-\beta {\Delta}\right)\right]\sum\limits _{i=1}^{n}{x}_{i}-\frac{1}{2{\sigma }^{2}}\sum\limits _{i=1}^{n}{x}_{i}^{2}+\frac{\beta }{{\sigma }^{2}}\sum\limits _{i=1}^{n}\sum\limits _{j\in {\eta }_{i}}{x}_{i}{x}_{j}\right.\hfill \\ \hfill & \quad \left.-\left[\frac{\beta \mu }{{\sigma }^{2}}(1-\beta {\Delta})\right]\sum\limits _{i=1}^{n}\sum\limits _{j\in {\eta }_{i}}{x}_{j}-\frac{\beta }{2{\sigma }^{2}}\sum\limits _{i=1}^{n}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{x}_{j}{x}_{k}\right\}\hfill \end{align} \tag{ 3 }$$

Through direct observation, we can identify the correspondence:

$$\begin{align}\hfill \vec{c}& =\left(\left[\frac{\mu }{{\sigma }^{2}}\left(1-\beta {\Delta}\right)\right],-\frac{1}{2{\sigma }^{2}},\frac{\beta }{{\sigma }^{2}},-\left[\frac{\beta \mu }{{\sigma }^{2}}(1-\beta {\Delta})\right],-\frac{\beta }{2{\sigma }^{2}}\right)\hfill \\ \hfill \vec{T}& =\left(\sum\limits _{i=1}^{n}{x}_{i},\sum\limits _{i=1}^{n}{x}_{i}^{2},\sum\limits _{i=1}^{n}\sum\limits _{j\in {\eta }_{i}}{x}_{i}{x}_{j},\sum\limits _{i=1}^{n}\sum\limits _{j\in {\eta }_{i}}{x}_{j},\sum\limits _{i=1}^{n}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{x}_{j}{x}_{k}\right)\hfill \end{align} \tag{ 4 }$$

with S(X) = 0 and

$$\begin{equation}d(\vec{\theta })=-\frac{n}{2}\left[\mathrm{log}(2\pi {\sigma }^{2})+\frac{{\mu }^{2}}{{\sigma }^{2}}\right]+\frac{\beta {\Delta}{\mu }^{2}n}{{\sigma }^{2}}\left[1-\frac{\beta {\Delta}}{2}\right]\end{equation} \tag{ 5 }$$

Furthermore, it is clear that if β = 0, the pseudo-likelihood function is simplified to the usual likelihood function of the Gaussian model, in which the number of parameters is equal to the number of natural parameters:

$$\begin{equation}F\left(\mathbf{X}\vert \vec{\theta }\right)=\mathrm{exp}\left\{\frac{\mu }{{\sigma }^{2}}\sum\limits _{i=1}^{n}{x}_{i}-\frac{1}{2{\sigma }^{2}}\sum\limits _{i=1}^{n}{x}_{i}^{2}-\frac{n}{2}\left[\mathrm{log}(2\pi {\sigma }^{2})+\frac{{\mu }^{2}}{{\sigma }^{2}}\right]\right\}\end{equation} \tag{ 6 }$$

where S(X) = 0 and:

$$\begin{equation}\vec{c}=\left(\frac{\mu }{{\sigma }^{2}},-\frac{1}{2{\sigma }^{2}}\right)\quad \vec{T}=\left(\sum\limits _{i=1}^{n}{x}_{i},\sum\limits _{i=1}^{n}{x}_{i}^{2}\right)\quad d(\vec{\theta })=-\frac{n}{2}\left[\mathrm{log}(2\pi {\sigma }^{2})+\frac{{\mu }^{2}}{{\sigma }^{2}}\right]\end{equation} \tag{ 7 }$$

The previous equations say that the inverse temperature parameter is responsible for making the model curved, as the number of natural parameters becomes greater than the number of parameters. Our goal is to measure how the emergence of the inverse temperature parameter into the model geometrically transforms the underlying parametric space. When the inverse temperature parameter is zero, the parametric space is a surface with constant negative curvature. However, when the inverse temperature increases up to a critical value, the geometry of parametric space suffers significant changes as the system undergoes phase transitions.

Differential geometry is a mature research field that is the mathematical background for several areas of physical sciences, remarkably to general relativity, in which it is employed to characterize the geometry of space-time [34]. But, what are the objects of study in differential geometry? Classical differential geometry is the study of local properties of curves and surfaces, where local refers to properties which depend on the behavior of the curve or surface in the neighborhood of a point [35, 36]. Before we proceed, it is interesting to give an intuition behind the concept of surface: a surface in R³ can be constructed by taking pieces of a plane, deforming them and arranging them in a way that the resulting shape has no sharp points, edges or self-intersections so that we can assign a tangent plane to every point of the surface [35–38]. In other words, a surface in R³ looks like an open subset of R² which has been smoothly deformed. An example of surface is the parametric space of a Gaussian random variable with a probability density function $p(x;\vec{\theta })$, in which $\vec{\theta }=(\mu ,{\sigma }^{2})$, where μ is the mean and σ² is the variance.

Definition 1 (Surface). A surface in R³ is a subset M ⊂ R³ such that for each point p ∈ M there exists a proper patch in M whose image contains a neighborhood of p in M, where a proper patch is a one-to-one mapping of an open set D ⊂ R² into R³.

In order to adapt the calculus of the plane R² to an arbitrary surface, a fundamental concept is the tangent plane. The intuition is that as lines are the simplest curves, the simplest surfaces are planes, so it is possible to approximate very complicated surfaces by simple planes.

Definition 2 (Tangent plane). Let M ⊂ R³ be a regular surface and p ∈ M an arbitrary point. A vector $\vec{v}$ is tangent to M at p provided there exists an > 0 and a smooth parametrized curve α : (−, ) → M with α(0) = p and ${\alpha }^{\prime }(0)=\vec{v}$. The set of all tangent vectors at p defines the tangent plane in p, denoted by T_p M.

The first thing a surface inhabitant needs to do in geometric terms is to measure the distance between two points and the angle between two vectors belonging to a tangent space. This is different from the distance between these points as measured by someone that lives in the ambient space, since, often, the straight line between the points in R³ will not be contained in the surface. The mathematical object that allows one to compute lengths on a surface, and also angles and areas, is the first fundamental form of the surface [38, 39].

Definition 3 (First fundamental form). Let p be an arbitrary point of a surface M. The first fundamental form of M at p associate to tangent vectors $\vec{w},\vec{z}\in {T}_{p}M$ the scalar:

$$\begin{equation}{\langle \vec{w},\vec{z}\rangle }_{p,M}=\vec{w}\cdot \vec{z}\end{equation} \tag{ 8 }$$

In other words, the first fundamental form allows us to compute dot products in the tangent plane.

Suppose x(u, v) defines a local parametrisation of a surface M. Then, any tangent vector to M at a point p can be expressed uniquely as a linear combination of ${\vec{x}}_{u}$ and ${\vec{x}}_{v}$, obtained by differentiating x(u, v) with respect to u and v, respectively. Hence, we can write:

$$\begin{equation}\vec{w}={\lambda }_{1}{\vec{x}}_{u}+{\lambda }_{2}{\vec{x}}_{v}\end{equation} \tag{ 9 }$$

$$\begin{equation}\vec{z}={\mu }_{1}{\vec{x}}_{u}+{\mu }_{2}{\vec{x}}_{v}\end{equation} \tag{ 10 }$$

and the inner product becomes:

$$\begin{align}\hfill {\langle \vec{w},\vec{z}\rangle }_{p,M}& =({\lambda }_{1}{\vec{x}}_{u}+{\lambda }_{2}{\vec{x}}_{v})\cdot ({\mu }_{1}{\vec{x}}_{u}+{\mu }_{2}{\vec{x}}_{v})\hfill \\ \hfill & ={\lambda }_{1}{\mu }_{1}E+({\lambda }_{1}{\mu }_{2}+{\lambda }_{2}{\mu }_{1})F+{\lambda }_{2}{\mu }_{2}G\hfill \end{align} \tag{ 11 }$$

where

$$\begin{equation}E={\vec{x}}_{u}\cdot {\vec{x}}_{u}\hspace{25.0pt}F={\vec{x}}_{u}\cdot {\vec{x}}_{v}\hspace{25.0pt}G={\vec{x}}_{v}\cdot {\vec{x}}_{v}\end{equation} \tag{ 12 }$$

The functions E, F and G are called the coefficients of the first fundamental form of the surface. This structure is important because it enables the computation of arc lengths of curves on the surface and the areas of regions on the surface. Note that an infinitesimal displacement in the surface ds may be expressed in terms of the coefficients of the first fundamental form as:

$$\begin{equation}\mathrm{d}{s}^{2}=E\,\mathrm{d}{u}^{2}+2F\,\mathrm{d}u\,\mathrm{d}v+G\,\mathrm{d}{v}^{2}\end{equation} \tag{ 13 }$$

In many applications, the first fundamental form is often written in the modern notation of the metric tensor:

$$\begin{equation}{\langle \vec{w},\vec{z}\rangle }_{p,M}={\vec{w}}^{T}\left[\begin{matrix}\hfill E\hfill & \hfill F\hfill \\ \hfill F\hfill & \hfill G\hfill \end{matrix}\right]\vec{z}\end{equation} \tag{ 14 }$$

where F = 0 implies in orthogonal coordinate systems. Note that if we have E = G = 1 and F = 0, we achieve the Euclidean geometry. It is usual to denote the first fundamental form by $\mathbb{I}$.

The first fundamental form relates to the study of intrinsic properties of surfaces. However, there are several extrinsic properties that are relevant for the complete characterization of a surface, such as the Gaussian curvature and the mean curvature. In summary, to measure how a surface is curving at a point p, we need to measure the rate of change of the unit normal vector at p. This discussion involves the definition of the second fundamental form of a surface [40].

Definition 4 (Second fundamental form). Let x(u, v) be a local parametrisation of a surface M with standard unit normal $\vec{N}$. As the local coordinates (u, v) change to (u + Δu, v + Δv) the surface moves away from the tangent plane by the distance d defined by:

$$\begin{equation}d=(x(u+{\Delta}u,v+{\Delta}v)-x(u,v))\cdot \vec{N}\end{equation} \tag{ 15 }$$

as illustrates figure 2.

Zoom In Zoom Out Reset image size

By a Taylor expansion, it is possible to express x(u + Δu, v + Δv) − x(u, v) as:

$$\begin{equation}{x}_{u}{\Delta}u+{x}_{v}{\Delta}v+\frac{1}{2}\left({x}_{uu}{({\Delta}u)}^{2}+2{x}_{uv}{\Delta}u{\Delta}v+{x}_{vv}{({\Delta}v)}^{2}\right)+r({\Delta}u,{\Delta}v)\end{equation} \tag{ 16 }$$

where the remainder function r becomes negligible when (Δu)² + (Δv)² tends to zero (infinitesimal displacement). Note that as both x_u and x_v belong to the tangent place, they are orthogonal to the normal vector $\vec{n}$, which means that the deviation d can be expressed as:

$$\begin{equation}d=\frac{1}{2}\left(L{({\Delta}u)}^{2}+2M{\Delta}u{\Delta}v+N{({\Delta}v)}^{2}\right)\end{equation} \tag{ 17 }$$

where

$$\begin{equation}L={x}_{uu}\cdot \vec{n}\hspace{25.0pt}M={x}_{uv}\cdot \vec{n}\hspace{25.0pt}N={x}_{vv}\cdot \vec{n}\end{equation} \tag{ 18 }$$

When Δu → 0 and Δv → 0, the resulting expression:

$$\begin{equation}L\,\mathrm{d}{u}^{2}+2M\,\mathrm{d}u\,\mathrm{d}v+N\,\mathrm{d}{v}^{2}\end{equation} \tag{ 19 }$$

defines the second fundamental form of the surface, which, often, is denoted by $\mathbb{II}$ and organized as a square matrix:

$$\begin{equation}\mathbb{II}=\left[\begin{matrix}\hfill L\hfill & \hfill M\hfill \\ \hfill M\hfill & \hfill N\hfill \end{matrix}\right]\end{equation} \tag{ 20 }$$

In order to compute the surface curvatures, it is necessary to obtain the expressions of the second fundamental form and of the differential of the Gauss map in a coordinate system. In computational terms, this can be done through the definition of the shape operator [35].

Definition 5 (Shape operator). Let M be a surface with first fundamental form $\mathbb{I}$ and second fundamental form $\mathbb{II}$. Then, the shape operator P can be computed by:

$$\begin{equation}P=-(\mathbb{II}){(\mathbb{I})}^{-1}=-\left[\begin{matrix}\hfill L\hfill & \hfill M\hfill \\ \hfill M\hfill & \hfill N\hfill \end{matrix}\right]{\left[\begin{matrix}\hfill E\hfill & \hfill F\hfill \\ \hfill F\hfill & \hfill G\hfill \end{matrix}\right]}^{-1}\end{equation} \tag{ 21 }$$

The shape operator encodes relevant information about the curvature of surfaces, being a powerful mathematical tool for geometric analysis. From the shape operator, we can obtain the Gaussian, mean and principal curvatures.

Definition 6 (Principal, mean and Gaussian curvatures). Let M be a surface and P its shape operator. It can be shown that [35]:

• (a)  
  The principal curvatures are the eigenvalues of the shape operator P.
• (b)  
  The mean curvature, K_M, is the trace of the shape operator P.
• (c)  
  The Gaussian curvature, K_G, is the determinant of the shape operator P.

In the next sections, we will discuss how to derive the first and second fundamental forms of Gaussian random field manifolds by means of the first-order and second-order Fisher information matrices. Undoubtedly, one of the most important findings of Carl Friedrich Gauss is that the Gaussian curvature is unchanged when the surface is bent without stretching, which is formally stated by the remarkable theorema egregium [38]. In other words, this result shows that the curvatures of a surface can be computed using only the first fundamental form. The implication of this theorem to our study is related to the amount of stretching the parametric space of Gaussian random fields suffers when it bends depends whether the inverse temperature is increasing or decreasing.

Information geometry is a recent research field that combines differential geometry and information theory to study intrinsic properties of the parametric spaces of several statistical models [18, 20]. According to this theory, the first-order Fisher information matrix is the metric tensor that equips the underlying parametric space of a statistical manifold (first fundamental form), whereas the second-order Fisher information matrix is the second fundamental form of the parametric space. It allows the application of the same mathematical tools used in the study deterministic physical phenomena, such as general relativity, into the macroscopic description of random systems [41].

In practice, the metric tensor allows us to express the square of an infinitesimal displacement in the manifold, ds², as a function of an infinitesimal displacement in the tangent space. Assuming a surface M and a matrix notation, we have:

$$\begin{equation}\mathrm{d}{s}^{2}=\left[\begin{matrix}\hfill \mathrm{d}u\hfill & \hfill \mathrm{d}v\hfill \end{matrix}\right]\left[\begin{matrix}\hfill E\hfill & \hfill F\hfill \\ \hfill F\hfill & \hfill G\hfill \end{matrix}\right]\left[\begin{matrix}\hfill \mathrm{d}u\hfill \\ \hfill \mathrm{d}v\hfill \end{matrix}\right]=E\,\mathrm{d}{u}^{2}+2F\,\mathrm{d}u\,\mathrm{d}v+G\,\mathrm{d}{v}^{2}\end{equation} \tag{ 22 }$$

where the matrix of coefficients E, F, e G is the first fundamental form. If this matrix is positive definite, we have a Riemannian manifold. Note that in the Euclidean case, the metric tensor is the identity matrix (since the space is flat), and we have the known Pythagorean relation ds² = du² + dv².

The second-order Fisher information matrix (second fundamental form) encodes information about the curvature of the manifold and it is used in the definition of the shape operator, whose determinant, trace and eigenvalues correspond to the Gaussian, mean and principal curvatures. Our goal is to measure the variations in the geometrical structure of the parametric space of pairwise isotropic GMRFs along the evolution of these systems under phase transitions.

Definition 7 (First-order Fisher information matrix). Let $p(X;\vec{\theta })$ be a probability density function where $\vec{\theta }=({\theta }_{1},\dots ,{\theta }_{n})\in {\Theta}$ is the vector of parameters. The first-order Fisher information matrix, which is the first fundamental form of underlying parametric space, is defined as:

$$\begin{equation}{\left\{\mathbb{I}(\vec{\theta })\right\}}_{ij}=\mathbb{E}\left[\left(\frac{\partial }{\partial {\theta }_{i}}\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial {\theta }_{j}}\,\mathrm{log}\,p(X;\vec{\theta })\right)\right],\quad \text{for}\enspace i,j=1,\dots ,n\end{equation} \tag{ 23 }$$

It has been shown that, under certain regularity conditions, the information equality is valid, which means that both first and second-order Fisher information matrices are identical [42–44].

Definition 8 (Second-order Fisher information matrix). Let $p(X;\vec{\theta })$ be a probability density function where $\vec{\theta }=({\theta }_{1},\dots ,{\theta }_{n})\in {\Theta}$ is the vector of parameters. The second-order Fisher information matrix, which is the second fundamental form of underlying parametric space, is defined as:

$$\begin{equation}{\left\{\mathbb{II}(\vec{\theta })\right\}}_{ij}=-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}\,\mathrm{log}\,p(X;\vec{\theta })\right],\quad \text{for}\enspace i,j=1,\dots ,n\end{equation} \tag{ 24 }$$

4.1. Considerations about the information equality

In the following, we provide a brief discussion about the information equality condition, which is a valid property for several probability density function belonging to the exponential family when the observations are independent [45]. Let X be a random variable with a probability density function $p(X;\vec{\theta })$. First, note that:

$$\begin{align}\hfill \mathbb{E}\left[\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}\,\mathrm{log}\,p(x;\vec{\theta })\right]& =\mathbb{E}\left[\frac{\partial }{\partial {\theta }_{j}}\left(\frac{\partial }{\partial {\theta }_{i}}\,\mathrm{log}\,p(x;\vec{\theta })\right)\right]\hfill \\ \hfill & =\mathbb{E}\left[\frac{\partial }{\partial {\theta }_{j}}\left(\frac{1}{p(x;\vec{\theta })}\frac{\partial }{\partial {\theta }_{i}}p(x;\vec{\theta })\right)\right]\hfill \end{align} \tag{ 25 }$$

By direct application of the product rule, we have:

$$\begin{align}\hfill \mathbb{E}\left[\frac{\partial }{\partial {\theta }_{j}}\left(\frac{1}{p(x;\vec{\theta })}\frac{\partial }{\partial {\theta }_{i}}p(x;\vec{\theta })\right)\right]& =\mathbb{E}\left[-\frac{1}{p{(x;\vec{\theta })}^{2}}\frac{\partial }{\partial {\theta }_{i}}p(x;\vec{\theta })\frac{\partial }{\partial {\theta }_{j}}p(x;\vec{\theta })+\frac{1}{p(x;\vec{\theta })}\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}p(x;\vec{\theta })\right]\hfill \\ \hfill & =-\mathbb{E}\left[\left(\frac{1}{p(x;\vec{\theta })}\frac{\partial }{\partial {\theta }_{i}}p(x;\vec{\theta })\right)\left(\frac{1}{p(x;\vec{\theta })}\frac{\partial }{\partial {\theta }_{i}}p(x;\vec{\theta })\right)\right]\hfill \\ \hfill & \quad +\mathbb{E}\left[\frac{1}{p(x;\vec{\theta })}\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}p(x;\vec{\theta })\right]\hfill \end{align} \tag{ 26 }$$

Using the definition of the expectation operator, the second term of equation (26) can be expressed as:

$$\begin{align}\hfill \mathbb{E}\left[\frac{1}{p(x;\vec{\theta })}\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}p(x;\vec{\theta })\right]& =\int p(x;\vec{\theta })\frac{1}{p(x;\vec{\theta })}\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}p(x;\vec{\theta })\mathrm{d}x\hfill \\ \hfill & =\int \frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}p(x;\vec{\theta })\mathrm{d}x\hfill \end{align} \tag{ 27 }$$

Under certain regularity conditions it is possible to interchange the integration and differentiation operators:

$$\begin{equation}\int \frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}p(x;\vec{\theta })\mathrm{d}x=\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}\int p(x;\vec{\theta })\mathrm{d}x=\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}1=0\end{equation} \tag{ 28 }$$

Note also that the arguments in the expectation in the first term of equation (26) can be rewritten as:

$$\begin{equation}\frac{1}{p(x;\vec{\theta })}\frac{\partial }{\partial {\theta }_{i}}p(x;\vec{\theta })=\frac{\partial }{\partial {\theta }_{i}}\,\mathrm{log}\,p(x;\vec{\theta })\end{equation} \tag{ 29 }$$

which finally leads to the equality:

$$\begin{equation}\mathbb{E}\left[\left(\frac{\partial }{\partial {\theta }_{i}}\,\mathrm{log}\,p(x;\vec{\theta })\right)\left(\frac{\partial }{\partial {\theta }_{j}}\,\mathrm{log}\,p(x;\vec{\theta })\right)\right]=-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial {\theta }_{i}\partial {\theta }_{j}}\,\mathrm{log}\,p(x;\vec{\theta })\right]\end{equation} \tag{ 30 }$$

In geometric terms, it means that the elements of the first fundamental form $(\mathbb{I})$ are equal to the elements of the second fundamental form $(\mathbb{II})$, which means constant Gaussian curvature. We will see that this is not the case in pairwise isotropic GMRFs when the inverse temperature deviates from zero.

4.2. Derivation of the first fundamental form

In this section, we provide the mathematical derivation of the components of the first fundamental form. As the parametric space is a 3D manifold, the first-order Fisher information matrix has the shape:

$$\begin{equation}\mathbb{I}(\vec{\theta })=\left(\begin{matrix}\hfill A\hfill & \hfill B\hfill & \hfill C\hfill \\ \hfill D\hfill & \hfill E\hfill & \hfill F\hfill \\ \hfill G\hfill & \hfill H\hfill & \hfill I\hfill \end{matrix}\right)\end{equation} \tag{ 31 }$$

We will begin with the first component of the matrix, which involves the derivatives with respect to the μ parameter. Note that the first component of $I(\vec{\theta })$ is:

$$\begin{equation}A={\mathbb{I}}_{\mu \mu }(\vec{\theta })=\mathbb{E}\left[\left(\frac{\partial }{\partial \mu }\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial \mu }\,\mathrm{log}\,p(X;\vec{\theta })\right)\right]\end{equation} \tag{ 32 }$$

where $p(X;\vec{\theta })$ is replaced by the local conditional density function of the Gaussian random field, given by equation (1). The computation of the derivatives leads to:

$$\begin{equation}A=\mathbb{E}\left\{\frac{1}{{\sigma }^{2}}{\left(1-\beta {\Delta}\right)}^{2}\frac{1}{{\sigma }^{2}}{\left[\left({x}_{i}-\mu \right)-\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{j}-\mu \right)\right]}^{2}\right\}\end{equation} \tag{ 33 }$$

Expanding the square, we have:

$$\begin{align}\hfill A& =\frac{1}{{\sigma }^{2}}{\left(1-\beta {\Delta}\right)}^{2}\mathbb{E}\left\{\frac{1}{{\sigma }^{2}}\left[{\left({x}_{i}-\mu \right)}^{2}-2\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{i}-\mu \right)\left({x}_{j}-\mu \right)+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\left({x}_{j}-\mu \right)\left({x}_{k}-\mu \right)\right]\right\}\hfill \end{align} \tag{ 34 }$$

And simplifying the expected values, we reach:

$$\begin{equation}A=\frac{{\left(1-\beta {\Delta}\right)}^{2}}{{\sigma }^{2}}\left[1-\frac{1}{{\sigma }^{2}}\left(2\beta \sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right)\right]\end{equation} \tag{ 35 }$$

where Δ is the cardinality of the neighborhood system (Δ = 8 in a second-order system), σ_ij is the covariance between the central variable x_i and one of its neighbors x_j ∈ η_i and σ_jk is the covariance between two variables x_j and x_k belonging to the neighborhood η_i . The second component of the first fundamental form is:

$$\begin{equation}B={\mathbb{I}}_{\mu {\sigma }^{2}}(\vec{\theta })=\mathbb{E}\left[\left(\frac{\partial }{\partial \mu }\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial {\sigma }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right)\right]\end{equation} \tag{ 36 }$$

which leads to:

$$\begin{equation}B=\frac{(1-\beta {\Delta})}{2{\sigma }^{6}}\mathbb{E}\left\{{\left[\left({x}_{i}-\mu \right)-\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{j}-\mu \right)\right]}^{3}\right\}-\frac{(1-\beta {\Delta})}{2{\sigma }^{4}}\mathbb{E}\left\{\left({x}_{i}-\mu \right)-\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{j}-\mu \right)\right\}\end{equation} \tag{ 37 }$$

Note that second term of equation (37) is zero, since:

$$\begin{equation}\mathbb{E}\left[{x}_{i}-\mu \right]-\beta \sum\limits _{j\in {\eta }_{i}}\mathbb{E}\left[{x}_{j}-\mu \right]=0-0=0\end{equation} \tag{ 38 }$$

and the expansion of the first term of equation (37) leads to:

$$\begin{align}\hfill \mathbb{E}\left\{{\left[\left({x}_{i}-\mu \right)-\beta \sum\limits _{j\in {\eta }_{i}}\left({x}_{j}-\mu \right)\right]}^{3}\right\}& =\mathbb{E}\left[{\left({x}_{i}-\mu \right)}^{3}\right]-3\beta \sum\limits _{j\in {\eta }_{i}}\mathbb{E}\left[({x}_{i}-\mu )({x}_{i}-\mu )({x}_{j}-\mu )\right]\hfill \\ \hfill & \quad +3{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\mathbb{E}\left[({x}_{i}-\mu )({x}_{j}-\mu )({x}_{k}-\mu )\right]\hfill \\ \hfill & \quad -{\beta }^{3}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\mathbb{E}\left(({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )\right]\hfill \end{align} \tag{ 39 }$$

The first term of (39) is zero for Gaussian random variables, since every central moment of odd order is null. According to the Isserlis' theorem [46], it is trivial to see that in fact all the third order cross terms are null, therefore, B = 0. The third component of the first fundamental form is:

$$\begin{equation}C={\mathbb{I}}_{\mu \beta }(\vec{\theta })=\mathbb{E}\left[\left(\frac{\partial }{\partial \mu }\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial \beta }\,\mathrm{log}\,p(X;\vec{\theta })\right)\right]\end{equation} \tag{ 40 }$$

Plugging the local conditional density function and doing some basic algebra, we reach:

$$\begin{align}\hfill C& =\frac{(1-\beta {\Delta})}{{\sigma }^{4}}\left\{\mathbb{E}\left[({x}_{i}-\mu )({x}_{i}-\mu )({x}_{j}-\mu )\right]-2\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\mathbb{E}\left[({x}_{i}-\mu )({x}_{j}-\mu )({x}_{k}-\mu )\right]\right.\hfill \\ \hfill & \quad \left.+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\mathbb{E}\left[({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )\right]\right\}\hfill \end{align} \tag{ 41 }$$

Once again, all the third order moments are zero by the Isserlis' theorem, resulting in C = 0. For the next component, by the symmetry of the metric tensor, D = B = 0. In order to calculate fifth component of the first fundamental form, we have to compute:

$$\begin{equation}E={\mathbb{I}}_{{\sigma }^{2}{\sigma }^{2}}(\vec{\theta })=\mathbb{E}\left[\left(\frac{\partial }{\partial {\sigma }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial {\sigma }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right)\right]\end{equation} \tag{ 42 }$$

which is given by:

$$\begin{align}\hfill E& =\mathbb{E}\left\{{\left[-\frac{1}{2{\sigma }^{2}}+\frac{1}{2{\sigma }^{4}}\left({x}_{i}-\mu -\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\right]}^{2}\right\}\hfill \\ \hfill & =\frac{1}{4{\sigma }^{4}}-\frac{1}{2{\sigma }^{6}}\mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{2}\right\}+\frac{1}{4{\sigma }^{8}}\mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{4}\right\}\hfill \end{align} \tag{ 43 }$$

Note that the first expectation leads to the following equality:

$$\begin{equation}\mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{2}\right\}={\sigma }^{2}-2\beta \sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\end{equation} \tag{ 44 }$$

For the second expectation, we have:

$$\begin{equation}\mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{4}\right\}=\mathbb{E}\left[{({x}_{i}-\mu )}^{4}\right]-4\beta \sum\limits _{j\in {\eta }_{i}}\mathbb{E}\left[{({x}_{i}-\mu )}^{3}({x}_{j}-\mu )\right]\end{equation} \tag{ 45 }$$

$$\begin{align}\hfill & +6{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\mathbb{E}\left[{({x}_{i}-\mu )}^{2}({x}_{j}-\mu )({x}_{k}-\mu )\right]\hfill \\\ & -4{\beta }^{3}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\mathbb{E}\left[({x}_{i}-\mu )({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )\right]\hfill \\ & +{\beta }^{4}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}\mathbb{E}\left[({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )({x}_{m}-\mu )\right]\hfill \end{align} \tag{ 46 }$$

leading to five different expectation terms. We invoke the Isserlis' theorem for Gaussian random variables to express higher order moments in terms of second-order moments. Hence, after some algebraic manipulations, we have:

$$\begin{align}\hfill E& =\frac{1}{2{\sigma }^{4}}-\frac{1}{{\sigma }^{6}}\left[2\beta \sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right]\hfill \\ \hfill & \quad +\frac{1}{{\sigma }^{8}}\left[3{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{ik}-{\beta }^{3}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\left({\sigma }_{ij}{\sigma }_{kl}+{\sigma }_{ik}{\sigma }_{jl}+{\sigma }_{il}{\sigma }_{jk}\right)\right.\hfill \\ \hfill & \quad \left.+{\beta }^{4}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}\left({\sigma }_{jk}{\sigma }_{lm}+{\sigma }_{jl}{\sigma }_{km}+{\sigma }_{jm}{\sigma }_{kl}\right)\right]\hfill \end{align} \tag{ 47 }$$

The sixth component of the first fundamental form is given by:

$$\begin{equation}F={\mathbb{I}}_{{\sigma }^{2}\beta }(\vec{\theta })=\mathbb{E}\left[\left(\frac{\partial }{\partial {\sigma }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial \beta }\,\mathrm{log}\,p(X;\vec{\theta })\right)\right]\end{equation} \tag{ 48 }$$

which can be computed as:

$$\begin{align}\hfill F=& \,\mathbb{E}\left\{\left[-\frac{1}{2{\sigma }^{2}}+\frac{1}{2{\sigma }^{4}}{\left(({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)}^{2}\right]\right.\hfill \\ \hfill & \left.\times \left[\frac{1}{{\sigma }^{2}}\left(({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\left(\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\right]\right\}\hfill \\ \hfill & =-\frac{1}{2{\sigma }^{4}}\mathbb{E}\left\{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right\}\hfill \\ \hfill & +\frac{1}{2{\sigma }^{6}}\mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{3}\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right\}\hfill \end{align} \tag{ 49 }$$

By computing the first expectation, we have:

$$\begin{equation}\mathbb{E}\left\{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right\}=\sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\end{equation} \tag{ 50 }$$

The expansion of the second expectation leads to:

$$\begin{align}\hfill \mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{3}\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right\}& =\mathbb{E}\left\{\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right.\hfill \\ \hfill & \quad \times \left.\left[{({x}_{i}-\mu )}^{3}-3\beta \sum\limits _{j\in {\eta }_{i}}{({x}_{i}-\mu )}^{2}({x}_{j}-\mu )\right.\right.\hfill \\ \hfill & \quad \left.\left.+3{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}({x}_{i}-\mu )({x}_{j}-\mu )({x}_{k}-\mu )\right.\right.\hfill \\ \hfill & \quad \left.\left.-{\beta }^{3}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )\right]\right\}\hfill \end{align} \tag{ 51 }$$

Again, by direct application of the Isserlis' theorem to express higher-order cross moments in terms of second-order moments and after some simplifications, we have:

$$\begin{align}\hfill F& =\frac{1}{{\sigma }^{4}}\left[\sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right]\hfill \\ \hfill & \quad -\frac{1}{2{\sigma }^{6}}\left[6\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{ik}-3{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\left({\sigma }_{ij}{\sigma }_{kl}+{\sigma }_{ik}{\sigma }_{jl}+{\sigma }_{il}{\sigma }_{jk}\right)\right.\hfill \\ \hfill & \quad \left.+{\beta }^{3}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}\left({\sigma }_{jk}{\sigma }_{lm}+{\sigma }_{jl}{\sigma }_{km}+{\sigma }_{jm}{\sigma }_{kl}\right)\right]\hfill \end{align} \tag{ 52 }$$

It is straightforward to see that C = G = 0 and F = H, since the first fundamental form is symmetric. Finally, the last component is defined as:

$$\begin{equation}I={\mathbb{I}}_{\beta \beta }(\vec{\theta })=\mathbb{E}\left[\left(\frac{\partial }{\partial \beta }\,\mathrm{log}\,p(X;\vec{\theta })\right)\left(\frac{\partial }{\partial \beta }\,\mathrm{log}\,p(X;\vec{\theta })\right)\right]\end{equation} \tag{ 53 }$$

which is given by:

$$\begin{align}\hfill I& =\frac{1}{{\sigma }^{4}}\mathbb{E}\left\{{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{2}{\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{2}\right\}\hfill \\ \hfill & =\frac{1}{{\sigma }^{4}}\mathbb{E}\left\{\left[{({x}_{i}-\mu )}^{2}-2\beta \sum\limits _{j\in {\eta }_{i}}({x}_{i}-\mu )({x}_{j}-\mu )+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}({x}_{j}-\mu )({x}_{k}-\mu )\right]\right.\hfill \\ \hfill & \quad \left.\times \left[\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}({x}_{j}-\mu )({x}_{k}-\mu )\right]\right\}\hfill \\ \hfill & =\frac{1}{{\sigma }^{4}}\mathbb{E}\left\{\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}({x}_{i}-\mu )({x}_{i}-\mu )({x}_{j}-\mu )({x}_{k}-\mu )-2\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}({x}_{i}-\mu )({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )\right.\hfill \\ \hfill & \quad \left.+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}({x}_{j}-\mu )({x}_{k}-\mu )({x}_{l}-\mu )({x}_{m}-\mu )\right\}\hfill \end{align} \tag{ 54 }$$

Once again, by using the Isserlis' formula and some algebra, we have:

$$\begin{align}\hfill \mathbb{I}& =\frac{1}{{\sigma }^{2}}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}+\frac{1}{{\sigma }^{4}}\left[2\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{ik}-2\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\left({\sigma }_{ij}{\sigma }_{kl}+{\sigma }_{ik}{\sigma }_{jl}+{\sigma }_{il}{\sigma }_{jk}\right)\right.\hfill \\ \hfill & \quad \left.+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}\left({\sigma }_{jk}{\sigma }_{lm}+{\sigma }_{jl}{\sigma }_{km}+{\sigma }_{jm}{\sigma }_{kl}\right)\right]\hfill \end{align} \tag{ 55 }$$

concluding that the first fundamental form has the following structure:

$$\begin{equation}\mathbb{I}(\vec{\theta })=\left(\begin{matrix}\hfill A\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E\hfill & \hfill F\hfill \\ \hfill 0\hfill & \hfill F\hfill & \hfill I\hfill \end{matrix}\right)\end{equation} \tag{ 56 }$$

where the non-zero elements are used to define how we compute an infinitesimal displacement in the manifold (parametric space) around the point $\vec{p}=(\mu ,{\sigma }^{2},\beta )$:

$$\begin{align}\hfill {(\mathrm{d}s)}^{2}& =\left[\begin{matrix}\hfill \mathrm{d}\mu \hfill & \hfill \mathrm{d}{\sigma }^{2}\hfill & \hfill \mathrm{d}\beta \hfill \end{matrix}\right]\left[\begin{matrix}\hfill A\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E\hfill & \hfill F\hfill \\ \hfill 0\hfill & \hfill F\hfill & \hfill I\hfill \end{matrix}\right]\left[\begin{matrix}\hfill \mathrm{d}\mu \hfill \\ \hfill \mathrm{d}{\sigma }^{2}\hfill \\ \hfill \mathrm{d}\beta \hfill \end{matrix}\right]\hfill \\ \hfill & =A{(\mathrm{d}\mu )}^{2}+E{(\mathrm{d}{\sigma }^{2})}^{2}+I{(\mathrm{d}\beta )}^{2}+2F(\mathrm{d}\beta )(\mathrm{d}{\sigma }^{2})\hfill \end{align} \tag{ 57 }$$

4.3. Derivation of the second fundamental form

In the following, we proceed with the derivation of the second fundamental form of the model, by computing the elements of the second-order Fisher information matrix. As the parametric space is a 3D manifold, the second fundamental form has the shape:

$$\begin{equation}\mathbb{II}(\vec{\theta })=\left(\begin{matrix}\hfill L\hfill & \hfill M\hfill & \hfill N\hfill \\ \hfill O\hfill & \hfill P\hfill & \hfill Q\hfill \\ \hfill R\hfill & \hfill S\hfill & \hfill T\hfill \end{matrix}\right)\end{equation} \tag{ 58 }$$

The first component of the second fundamental form $\mathbb{II}(\vec{\theta })$ is:

$$\begin{equation}L={\mathbb{II}}_{\mu \mu }(\vec{\theta })=-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial {\mu }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right]\end{equation} \tag{ 59 }$$

which is given by:

$$\begin{equation}L=-\frac{\left(1-\beta {\Delta}\right)}{{\sigma }^{2}}\mathbb{E}\left\{\frac{\partial }{\partial \mu }\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right\}=\frac{1}{{\sigma }^{2}}{\left(1-\beta {\Delta}\right)}^{2}\end{equation} \tag{ 60 }$$

where Δ = 8 is the size of the neighborhood system. The second component is defined by:

$$\begin{equation}M={\mathbb{II}}_{\mu {\sigma }^{2}}(\vec{\theta })=-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial \mu \partial {\sigma }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right]\end{equation} \tag{ 61 }$$

resulting in

$$\begin{equation}M=\frac{1}{{\sigma }^{4}}(1-\beta {\Delta})\mathbb{E}\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]=\frac{1}{{\sigma }^{4}}(1-\beta {\Delta})\left[0-0\right]=0\end{equation} \tag{ 62 }$$

The third component of the second fundamental form is also zero:

$$\begin{align}\hfill N={\mathbb{II}}_{\mu \beta }(\vec{\theta })& =-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial \mu \partial \beta }\,\mathrm{log}\,p(X;\vec{\theta })\right]\hfill \\ \hfill & =\frac{1}{{\sigma }^{2}}\mathbb{E}\left\{{\Delta}\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]+(1-\beta {\Delta})\left[\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]\right\}\hfill \\ \hfill & =0+0=0\hfill \end{align} \tag{ 63 }$$

Moving forward to the fourth component, note that $O={\mathbb{II}}_{{\sigma }^{2}\mu }(\vec{\theta })=0$, since a change in the order of the differentiation operators does not affect the result. Thus, we proceed directly to the fifth component, given by:

$$\begin{align}\hfill P={\mathbb{II}}_{{\sigma }^{2}{\sigma }^{2}}(\vec{\theta })& =-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial {({\sigma }^{2})}^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right]\hfill \\ \hfill & =-\mathbb{E}\left\{\frac{\partial }{\partial {\sigma }^{2}}\left[-\frac{1}{2{\sigma }^{2}}+\frac{1}{2{\sigma }^{4}}{\left({x}_{i}-\mu -\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)}^{2}\right]\right\}\hfill \\ \hfill & =-\mathbb{E}\left\{\frac{1}{2{\sigma }^{4}}-\frac{1}{{\sigma }^{6}}{\left[({x}_{i}-\mu )-\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right]}^{2}\right\}\hfill \\ \hfill & =\frac{1}{2{\sigma }^{4}}-\frac{1}{{\sigma }^{6}}\left[2\beta \sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right]\hfill \end{align} \tag{ 64 }$$

The next component of the second fundamental form is given by:

$$\begin{align}\hfill Q={\mathbb{II}}_{{\sigma }^{2}\beta }(\vec{\theta })& =-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial {\sigma }^{2}\partial \beta }\,\mathrm{log}\,p(X;\vec{\theta })\right]\hfill \\ \hfill & =-\mathbb{E}\left\{\frac{\partial }{\partial {\sigma }^{2}}\left[\frac{1}{{\sigma }^{2}}\left({x}_{i}-\mu -\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\left(\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\right]\right\}\hfill \\ \hfill & =\frac{1}{{\sigma }^{4}}\left[\sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right]\hfill \end{align} \tag{ 65 }$$

Note that the next two components of the second fundamental form are identical to their symmetric counterparts, that is, $R={\mathbb{II}}_{\beta \mu }(\vec{\theta })={\mathbb{II}}_{\mu \beta }(\vec{\theta })=N=0$ and $S={\mathbb{II}}_{\beta {\sigma }^{2}}(\vec{\theta })={\mathbb{II}}_{{\sigma }^{2}\beta }(\vec{\theta })=Q$. Finally, the last component of the second fundamental form is given by:

$$\begin{equation}T={\mathbb{II}}_{\beta \beta }(\vec{\theta })=-\mathbb{E}\left[\frac{{\partial }^{2}}{\partial {\beta }^{2}}\,\mathrm{log}\,p(X;\vec{\theta })\right]\end{equation} \tag{ 66 }$$

which can be computed as:

$$\begin{align}\hfill T& =-\frac{1}{{\sigma }^{2}}\mathbb{E}\left\{\frac{\partial }{\partial \beta }\left[\left({x}_{i}-\mu -\beta \sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\left(\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\right]\right\}\hfill \\ \hfill & =\frac{1}{{\sigma }^{2}}\mathbb{E}\left[\left(\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\left(\sum\limits _{j\in {\eta }_{i}}({x}_{j}-\mu )\right)\right]\hfill \\ \hfill & =\frac{1}{{\sigma }^{2}}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\hfill \end{align} \tag{ 67 }$$

concluding that the second fundamental form has the following structure:

$$\begin{equation}\mathbb{II}(\vec{\theta })=\left(\begin{matrix}\hfill L\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill P\hfill & \hfill Q\hfill \\ \hfill 0\hfill & \hfill Q\hfill & \hfill T\hfill \end{matrix}\right)\end{equation} \tag{ 68 }$$

Finally, note also that when the inverse temperature parameter is fixed at zero, both first and second fundamental forms converge to:

$$\begin{equation}g(\vec{\theta })=\left(\begin{matrix}\hfill \frac{1}{{\sigma }^{2}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{1}{2{\sigma }^{4}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\end{equation} \tag{ 69 }$$

which is exactly the metric tensor (Fisher information matrix) of the parametric space of a single Gaussian random variable.

4.4. Derivation of the Gaussian curvature

Given the fact that both the first and second fundamental forms are have been calculated, in this subsection we explicitly derive a closed-form expression for the Gaussian curvature, showing how it depends on the inverse temperature. First, we have to calculate the inverse of the first fundamental form, that is, $\mathbb{I}{(\vec{\theta })}^{-1}$. To find the elements of the matrix $\mathbb{I}{(\vec{\theta })}^{-1}$, note that:

$$\begin{equation}\left(\begin{matrix}\hfill A\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill E\hfill & \hfill F\hfill \\ \hfill 0\hfill & \hfill F\hfill & \hfill I\hfill \end{matrix}\right)\left(\begin{matrix}\hfill a\hfill & \hfill b\hfill & \hfill c\hfill \\ \hfill d\hfill & \hfill e\hfill & \hfill f\hfill \\ \hfill g\hfill & \hfill h\hfill & \hfill i\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)\end{equation} \tag{ 70 }$$

By solving this linear system, it is not difficult to see that the inverse of the first fundamental form is given by:

$$\begin{equation}\mathbb{I}{(\vec{\theta })}^{-1}=\left(\begin{matrix}\hfill \frac{1}{A}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{I}{EI-{F}^{2}}\hfill & \hfill -\frac{F}{EI-{F}^{2}}\hfill \\ \hfill 0\hfill & \hfill -\frac{F}{EI-{F}^{2}}\hfill & \hfill \frac{E}{EI-{F}^{2}}\hfill \end{matrix}\right)\end{equation} \tag{ 71 }$$

The shape operator can be calculated by the direct multiplication of the second fundamental form and the inverse of the first fundamental form, which leads to:

$$\begin{equation}P(\vec{\theta })=\left(\begin{matrix}\hfill L\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill P\hfill & \hfill Q\hfill \\ \hfill 0\hfill & \hfill Q\hfill & \hfill T\hfill \end{matrix}\right)\left(\begin{matrix}\hfill \frac{1}{A}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{I}{EI-{F}^{2}}\hfill & \hfill -\frac{F}{EI-{F}^{2}}\hfill \\ \hfill 0\hfill & \hfill -\frac{F}{EI-{F}^{2}}\hfill & \hfill \frac{E}{EI-{F}^{2}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill \frac{L}{A}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{PI-QF}{EI-{F}^{2}}\hfill & \hfill -\frac{QE-PF}{EI-{F}^{2}}\hfill \\ \hfill 0\hfill & \hfill \frac{QI-TF}{EI-{F}^{2}}\hfill & \hfill \frac{ET-QF}{EI-{F}^{2}}\hfill \end{matrix}\right)\end{equation} \tag{ 72 }$$

The Gaussian curvature, K_G, is the determinant of the shape operator. Using the rule of Sarrus, we have:

$$\begin{align}\hfill \vert P(\vec{\theta })\vert ={K}_{\text{G}}& =\frac{L}{A}\frac{(PI-QF)(ET-QF)}{{(EI-{F}^{2})}^{2}}-\frac{L}{A}\frac{(QE-PF)(QI-TF)}{{(EI-{F}^{2})}^{2}}=\hfill \\ \hfill & =\frac{L}{A}\frac{1}{{(EI-{F}^{2})}^{2}}\left(PIET+{Q}^{2}{F}^{2}-{Q}^{2}EI-P{F}^{2}T\right)\hfill \end{align} \tag{ 73 }$$

Note that the above expression is equivalent to:

$$\begin{equation}{K}_{\text{G}}=\frac{L}{A}\frac{1}{{(EI-{F}^{2})}^{2}}\left(PT(EI-{F}^{2})-{Q}^{2}(EI-{F}^{2})\right)\end{equation} \tag{ 74 }$$

which finally leads to:

$$\begin{equation}{K}_{\text{G}}=\frac{L}{A}\left(PT-{Q}^{2}\right)\end{equation} \tag{ 75 }$$

We will see in the next sections that the quantity Q in the second fundamental form is highly asymmetric with respect to variations in the inverse temperature. Knowing that we have the individual expressions for P and T, we can express their product as:

$$\begin{align}\hfill PT& =\left[\frac{1}{2{\sigma }^{4}}-\frac{1}{{\sigma }^{6}}\left(2\beta \sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right)\right]\left[\frac{1}{{\sigma }^{2}}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right]\hfill \\ \hfill & =\frac{1}{2{\sigma }^{6}}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}-\frac{1}{{\sigma }^{8}}\left[2\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{kl}-{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}{\sigma }_{jk}{\sigma }_{lm}\right]\hfill \end{align} \tag{ 76 }$$

Calculating the square of Q, we have:

$$\begin{equation}{Q}^{2}=\frac{1}{{\sigma }^{8}}\left[\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{ik}-2\beta \sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{kl}+{\beta }^{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}\sum\limits _{l\in {\eta }_{i}}\sum\limits _{m\in {\eta }_{i}}{\sigma }_{jk}{\sigma }_{lm}\right]\end{equation} \tag{ 77 }$$

Therefore, we can express the Gaussian curvature as:

$$\begin{equation}{K}_{\text{G}}=\frac{\frac{1}{2{\sigma }^{6}}{\sum }_{j\in {\eta }_{i}}{\sum }_{k\in {\eta }_{i}}{\sigma }_{jk}-\frac{1}{{\sigma }^{8}}{\sum }_{j\in {\eta }_{i}}{\sum }_{k\in {\eta }_{i}}{\sigma }_{ij}{\sigma }_{ik}}{1-\frac{1}{{\sigma }^{2}}\left(2\beta {\sum }_{j\in {\eta }_{i}}{\sigma }_{ij}-{\beta }^{2}{\sum }_{j\in {\eta }_{i}}{\sum }_{k\in {\eta }_{i}}{\sigma }_{jk}\right)}\end{equation} \tag{ 78 }$$

It is possible to express the Gaussian curvature in terms of the Pearson correlation coefficients ρ_ij . Recall that:

$$\begin{equation}{\rho }_{ij}=\frac{{\sigma }_{ij}}{{\sigma }_{i}{\sigma }_{j}}\end{equation} \tag{ 79 }$$

where σ_i denotes the standard deviation of the random variable x_i . If the variance is invariant along the random field, which is the case in our model, the Gaussian curvature can be expressed by:

$$\begin{equation}{K}_{\text{G}}=\frac{1}{{\sigma }^{4}}\frac{\left(\frac{1}{2}{\sum }_{j\in {\eta }_{i}}{\sum }_{k\in {\eta }_{i}}{\rho }_{jk}-{\sum }_{j\in {\eta }_{i}}{\sum }_{k\in {\eta }_{i}}{\rho }_{ij}{\rho }_{ik}\right)}{\left(1-2\beta {\sum }_{j\in {\eta }_{i}}{\rho }_{ij}+{\beta }^{2}{\sum }_{j\in {\eta }_{i}}{\sum }_{k\in {\eta }_{i}}{\rho }_{jk}\right)}\end{equation} \tag{ 80 }$$

4.5. Fundamental forms in tensorial notation

The direct computation of triple and quadruple summations has time complexity of O(n³) and O(n⁴), respectively, which is quite large. To reduce the computational cost of these operations in our simulations, we express the components of the first and second fundamental forms using Kronecker products (tensor products). First, note that each 3 × 3 matrix that defines a neighborhood patch composed by x_i ∪ η_i can be transformed into a vector p_i of nine elements by piling its rows. Hence, it is possible to compute the covariance matrix of these vectors, denoted by Σ_p . From this covariance matrix, we extract two components: (1) a vector of eight elements, $\vec{\rho }$, formed by the components of the central row of Σ_p , excluding the middle one, which denotes the variance of x_i (only the covariances between x_i and x_j , for j ≠ i); and (2) an 8 × 8 sub-matrix, ${{\Sigma}}_{p}^{-}$, obtained by removing the central row and central column from Σ_p (only the covariances between x_j ∈ η_i and x_k ∈ η_i ). Figure 3 shows a diagram illustrating the decomposition of the covariance matrix Σ_p into the sub-matrix ${{\Sigma}}_{p}^{-}$ and the vector $\vec{\rho }$. Therefore, making use of Kronecker products, it is possible to rewrite the first fundamental form (metric tensor) in a tensorial notation, providing a computationally efficient way to compute the elements of $I(\vec{\theta })$:

$$\begin{equation}A={\mathbb{I}}_{\mu \mu }(\vec{\theta })=\frac{1}{{\sigma }^{2}}{\left(1-\beta {\Delta}\right)}^{2}\left[1-\frac{1}{{\sigma }^{2}}\left(2\beta {{\Vert}\vec{\rho }{\Vert}}_{+}-{\beta }^{2}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right)\right]\end{equation} \tag{ 81 }$$

$$\begin{align}\hfill E& ={\mathbb{I}}_{{\sigma }^{2}{\sigma }^{2}}(\vec{\theta })=\frac{1}{2{\sigma }^{4}}-\frac{1}{{\sigma }^{6}}\left[2\beta {{\Vert}\vec{\rho }{\Vert}}_{+}-{\beta }^{2}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\hfill \\ \hfill & +\frac{1}{{\sigma }^{8}}\left[3{\beta }^{2}{{\Vert}\vec{\rho }\otimes \vec{\rho }{\Vert}}_{+}-3{\beta }^{3}{{\Vert}\vec{\rho }\otimes {{\Sigma}}_{p}^{-}{\Vert}}_{+}+3{\beta }^{4}{{\Vert}{{\Sigma}}_{p}^{-}\otimes {{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\hfill \end{align} \tag{ 82 }$$

$$\begin{align}\hfill F& ={\mathbb{I}}_{{\sigma }^{2}\beta }(\vec{\theta })={\mathbb{I}}_{\beta {\sigma }^{2}}(\vec{\theta })=\frac{1}{{\sigma }^{4}}\left[{{\Vert}\vec{\rho }{\Vert}}_{+}-\beta {{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\hfill \\ \hfill & -\frac{1}{2{\sigma }^{6}}\left[6\beta {{\Vert}\vec{\rho }\otimes \vec{\rho }{\Vert}}_{+}-9{\beta }^{2}{{\Vert}\vec{\rho }\otimes {{\Sigma}}_{p}^{-}{\Vert}}_{+}+3{\beta }^{3}{{\Vert}{{\Sigma}}_{p}^{-}\otimes {{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\hfill \end{align} \tag{ 83 }$$

$$\begin{equation}I={\mathbb{I}}_{\beta \beta }(\vec{\theta })=\frac{1}{{\sigma }^{2}}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}+\frac{1}{{\sigma }^{4}}\left[2{{\Vert}\vec{\rho }\otimes \vec{\rho }{\Vert}}_{+}-6\beta {{\Vert}\vec{\rho }\otimes {{\Sigma}}_{p}^{-}{\Vert}}_{+}+3{\beta }^{2}{{\Vert}{{\Sigma}}_{p}^{-}\otimes {{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\end{equation} \tag{ 84 }$$

where ${{\Vert}A{\Vert}}_{+}$ represents the summation of all the entries of the vector/matrix A and ⊗ denotes the Kronecker (tensor) product. Similarly, by employing the same reasoning, the elements of the second fundamental form $\mathbb{II}(\vec{\theta })$ can be expressed as:

$$\begin{equation}L={\mathbb{II}}_{\mu \mu }(\vec{\theta })=\frac{1}{{\sigma }^{2}}{\left(1-\beta {\Delta}\right)}^{2}\end{equation} \tag{ 85 }$$

$$\begin{equation}P={\mathbb{II}}_{{\sigma }^{2}{\sigma }^{2}}(\vec{\theta })=\frac{1}{2{\sigma }^{4}}-\frac{1}{{\sigma }^{6}}\left[2\beta {{\Vert}\vec{\rho }{\Vert}}_{+}-{\beta }^{2}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\end{equation} \tag{ 86 }$$

$$\begin{equation}Q={\mathbb{II}}_{{\sigma }^{2}\beta }(\vec{\theta })=I{I}_{\beta {\sigma }^{2}}(\vec{\theta })=\frac{1}{{\sigma }^{4}}\left[{{\Vert}\vec{\rho }{\Vert}}_{+}-\beta {{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right]\end{equation} \tag{ 87 }$$

$$\begin{equation}T={\mathbb{II}}_{\beta \beta }(\vec{\theta })=\frac{1}{{\sigma }^{2}}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\end{equation} \tag{ 88 }$$

Zoom In Zoom Out Reset image size

4.6. Relation between entropy and the second fundamental form

A relevant quantity in the study of stochastic complex systems is entropy, one of the most ubiquitous concepts in science, with applications in a large number of research fields. For instance, in information theory, entropy is related to the degree of uncertainty about a source of information [47]. In statistical physics, entropy plays an important role, being a central piece in the second law of thermodynamics, which states that the entropy of isolated systems left to spontaneous evolution cannot decrease, as they always arrive at a state of thermodynamic equilibrium, where the entropy is highest [48].

To compute the entropy in a pairwise isotropic GMRF, recall that it can be defined as the expected value of self-information, which leads to:

$$\begin{align}\hfill {H}_{\beta }(\vec{\theta })& =-E\left[\mathrm{log}\,p\left({x}_{i}\vert {\eta }_{i},\vec{\theta }\right)\right]\hfill \\ \hfill & =\frac{1}{2}\left[\mathrm{log}\left(2\pi {\sigma }^{2}\right)+1\right]-\frac{1}{{\sigma }^{2}}\left[\beta \sum\limits _{j\in {\eta }_{i}}{\sigma }_{ij}-\frac{{\beta }^{2}}{2}\sum\limits _{j\in {\eta }_{i}}\sum\limits _{k\in {\eta }_{i}}{\sigma }_{jk}\right]\hfill \\ \hfill & ={H}_{\text{G}}(\vec{\theta })-\frac{1}{{\sigma }^{2}}\left(\beta {{\Vert}\vec{\rho }{\Vert}}_{+}+\frac{{\beta }^{2}}{2}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right)\hfill \end{align} \tag{ 89 }$$

where ${H}_{\text{G}}(\vec{\theta })$ denotes the entropy of a Gaussian random variable. Note that the entropy is a quadratic function of the inverse temperature parameter β. Besides, for β = 0, we have ${H}_{\beta }(\vec{\theta })={H}_{\text{G}}(\vec{\theta })$, as expected. Using the Kronecker product and rewriting the summations in terms of the components Q and T from the second fundamental form, we have:

$$\begin{equation}{H}_{\beta }(\vec{\theta })={H}_{\text{G}}(\vec{\theta })-\beta {\sigma }^{2}Q-{\beta }^{2}\frac{T}{2}\end{equation} \tag{ 90 }$$

Note that all the components of the second fundamental form can be expressed in terms of the variance of the random variables and the system's entropy. First, note that:

$$\begin{equation}{\Delta}{H}_{\beta }={H}_{\beta }(\vec{\theta })-{H}_{G}(\vec{\theta })=-\frac{1}{{\sigma }^{2}}\left(\beta {{\Vert}\vec{\rho }{\Vert}}_{+}-\frac{{\beta }^{2}}{2}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}\right)\end{equation} \tag{ 91 }$$

which leads to:

$$\begin{equation}2\beta {{\Vert}\vec{\rho }{\Vert}}_{+}-{\beta }^{2}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}=-2{\sigma }^{2}{\Delta}{H}_{\beta }\end{equation} \tag{ 92 }$$

and therefore:

$$\begin{equation}P=\frac{1}{2{\sigma }^{4}}+\frac{2}{{\sigma }^{4}}{\Delta}{H}_{\beta }=\frac{1}{{\sigma }^{4}}\left(\frac{1}{2}+2{\Delta}{H}_{\beta }\right)\end{equation} \tag{ 93 }$$

Then, note that:

$$\begin{equation}{{\Vert}\vec{\rho }{\Vert}}_{+}-\beta {{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}=-{\sigma }^{2}\frac{\partial }{\partial \beta }{H}_{\beta }\end{equation} \tag{ 94 }$$

leading to:

$$\begin{equation}Q=-\frac{1}{{\sigma }^{2}}\frac{\partial }{\partial \beta }{H}_{\beta }\end{equation} \tag{ 95 }$$

Finally, note that:

$$\begin{equation}{{\Vert}{{\Sigma}}_{p}^{-}{\Vert}}_{+}={\sigma }^{2}\frac{{\partial }^{2}}{\partial {\beta }^{2}}{H}_{\beta }\end{equation} \tag{ 96 }$$

which leads to:

$$\begin{equation}T=\frac{{\partial }^{2}}{\partial {\beta }^{2}}{H}_{\beta }\end{equation} \tag{ 97 }$$

In this section we present and discuss the obtained results from our computational simulations with Gaussian random fields dynamics.

5.1. Geometric analysis of Gaussian random field dynamics

In order to simulate the dynamics of the Gaussian random field, MCMC simulation was employed to generate outcomes of the model using the Metropolis–Hastings algorithm [49]. For each simulation, arbitrary initial values for the parameters μ and σ² are selected, but the initial inverse temperature is set to zero, making all random variables independent. At the end of each iteration, we perform a small and positive small displacement Δβ until a maximum value β_MAX is reached, defining the first half of an information cycle (entropy is increasing). After that, to complete a cycle, at each iteration, the same negative small displacement −Δβ is performed, until the inverse temperature reaches zero once again (entropy is decreased). The other parameters of the random field (mean and variance) are estimated by the sample mean and sample variance at each iteration. In our simulations, a complete full cycle takes 1000 iterations, where each random field outcome is represented by a 512 × 512 matrix of real numbers. Figure 4 shows some samples of the random field during the evolution of the system.

Zoom In Zoom Out Reset image size

To move forward, it is necessary to explain how the results were generated. First, recall that a full cycle of the MCMC simulation is composed by 1000 iterations of the Metropolis–Hastings algorithm: initially, the inverse temperature parameter β is set to zero, and at the end of each iteration, β is incremented by Δβ = 0.0006, up to the 500th iteration. In the second half of the cycle, at the end of each iteration, β is decremented by Δβ = 0.0006, until it reaches zero once again. Figure 5 shows the variation of the inverse temperature and entropy along a full cycle.

Zoom In Zoom Out Reset image size

As the entropy of a Gaussian random field depends directly on the components Q and T of the second fundamental form, which are respectively the negative of the first derivative scaled by the inverse of the variance and second derivative of the entropy itself, we also plot the variations of these quantities along a full cycle of the MCMC simulation. Figure 6 illustrates how these components change over time. Note that while the plot of T (second derivative) is smooth and almost symmetric, the plot of Q (negative of first derivative) non-smooth and highly asymmetric.

Zoom In Zoom Out Reset image size

The behavior of the Gaussian curvature along a full cycle of the numerical simulations shows an interesting pattern. In the beginning, when the inverse temperature is zero, the parametric space has constant negative Gaussian curvature (K = −1), which means hyperbolic geometry. During the first half, when the system moves towards higher entropy states, the sign of the Gaussian curvature changes from negative to positive, whereas in the second half, when the system moves towards lower entropy states, the sign of the Gaussian curvature changes from positive to negative. However, the amount of curvature necessary to bend and stretch/shrink the parametric space when moving towards lower entropy states is significantly higher than that necessary to bend and stretch/shrink the parametric space when moving towards higher entropy states. We call this asymmetric pattern of evolution of Gaussian random fields as the curvature effect, which can be described as: the variations of the Gaussian curvature when the system moves towards higher entropy states is different from the variations of the Gaussian curvature when the system moves towards lower entropy states. Figure 7 illustrates the variation of the Gaussian curvature along the MCMC dynamics.

Zoom In Zoom Out Reset image size

An interesting question about the results obtained with the MCMC simulation is: when exactly the Gaussian curvature becomes positive? Our numerical computations show that the changes in the sign of the Gaussian curvature coincide with the moments of abrupt changes in the behavior of the system's entropy. Figure 8 shows the global configurations of the system in the points where the Gaussian curvature change its sign. It is possible to visualize that change in the sign of the Gaussian curvature is directly related to moments in which the random field is undergoing phase transitions. Therefore, our results suggest that the Gaussian curvature is a good measure to determine whether the system is close to a critical point, that is, if it is approaching a phase transition. However, the intrinsic geometry of the parametric space (curvature) during the first phase transition (from randomness to clusters) is quite distinct from the intrinsic geometry of the parametric space during the second phase transition (from clusters to randomness), leading to this strange pattern that we called the curvature effect.

Zoom In Zoom Out Reset image size

By inspecting the mean and the principal curvatures, we also note a highly asymmetric pattern of evolution during the dynamics. The difference in comparison with the Gaussian curvature is that the mean curvature is always negative due to the principal curvatures. The second principal curvature is the only one to become positive along the MCMC simulation. From differential geometry, we know that the mean curvature is the summation of the principal curvatures and the Gaussian curvature is the product of the principal curvatures. Figure 9 illustrates the variation of the mean and principal curvatures along the dynamics.

Zoom In Zoom Out Reset image size

In figure 10, we visualize the Gaussian curvature, entropy and the components of the second fundamental form T and Q used in the computation of the Gaussian curvature as a function of the inverse temperature. What happens with these plots is that, due to the hysteresis phenomenon, which is a characteristic of complex systems with memory, these functions does not behave like proper functions of the inverse temperature parameter along the numerical simulations, because, since we are dealing with a stochastic complex system, for values of beta belonging to a certain range, there are two distinct values, depending whether the system is moving towards higher entropy states (increasing beta) or moving towards lower entropy states (decreasing beta).

Zoom In Zoom Out Reset image size

In order to analyze how the system's entropy change as a function of the curvatures, we build a 2D representation to visualize the entropy as a function of the curvatures. The idea is to visualize what happens to the entropy as curvature changes along the MCMC simulation. In summary, we observe that this curvature effect indicates that the amount of curvature needed to bend the parametric space during an increase in the system's entropy (blue line) is smaller than the amount of curvature required to bend the parametric space during the same reduction in the system's entropy (red curve). Figure 11 shows the plots of the curvature loops for both Gaussian and mean curvatures. The resulting behavior resembles a mathematical model of hysteresis, which underlies a large number of phase transitions in physical models [50, 51].

Zoom In Zoom Out Reset image size

Note that this behavior induces a natural orientation to the process of, from a low entropy state, bringing the random field to a high entropy state and back, which is an intrinsic notion of arrow of time [52]. It is interesting to observe that in the Gaussian curvature cycle (left image), the natural orientation is clockwise (blue curve first, red curve after), while in the mean curvature cycle (right image), the natural orientation is anti-clockwise (blue curve first, red curve after). The Python source code with the computational implementation of the numerical simulations are available at: https://github.com/alexandrelevada/Curvature_GMRF.

Stochastic complex systems are composed of several random variables that interact in a non-linear way along time. Depending on how the parameter models change, a series of complex behavior can emerge from this dynamics. In this study, we addressed the problem of deriving closed-form expressions for the first and second fundamental forms of the underlying manifolds of Gaussian random fields in order to compute the principal, mean and Gaussian curvatures. Information geometry shows that the fundamental forms of parametric spaces of random variables are given by the first and second-order Fisher information matrices of the statistical model. Mathematical expressions for the components of these matrices were derived, allowing the computation of intrinsic geometric properties.

The intrinsic geometric structure of statistical models of independent random variables has been extensively studied in information geometry. For instance, it has been shown that the parametric space of a Gaussian random variable has constant negative Gaussian curvature, which means that the geometry is hyperbolic. However, little is known about the geometry of random field models, where the inverse temperature parameter induces a spatial dependence structure among the variables. In this paper, we investigated how the variation of the system's entropy is related to the variation in the principal, mean and Gaussian curvatures.

Using MCMC simulations via the Metropolis–Hastings algorithm, we performed several cycles of evolution composed by two phase transitions: a first one in which the system's entropy is increased, and a second one in which the system's entropy is decreased to the minimum value. The obtained results show that the variations in principal, mean and Gaussian curvatures are highly asymmetric, suggesting that the parametric space suffers a series of irreversible geometric deformations. Our geometric analysis has shown an unreported phenomenon: the curvature effect, which suggests that the deformations in the parametric space are more prominent during a decrease of the inverse temperature than during an increase of the inverse temperature, indicating the emergence of an arrow of time in the evolution of the random field.

Our choice for the lattice geometry was motivated by the fact that, according to discrete mathematics, lattices are partially ordered sets (POSETS), and therefore, it is possible to define neighbourhood systems as POSETS too. The main advantage is that we can organize neighborhood systems as vectors by using the lexicographic ordering, making it possible the computation of the covariance matrices of these patterns, which significantly reduces the computational cost of the numerical simulations. The choice of regular networks, that is, networks in which all the nodes have the same degree can be further investigated. However, for more complex topologies we believe that another kind of mathematical expressions are required, since, despite the possibility of replacing Δ by the average degree, the computation of the covariances between neighboring nodes would require some adjustments because some neighbourhoods would be smaller/greater than others. Anyway, this is an excellent suggestion for future works. We will try to investigate how to generalize the proposed equations to regular graphs and then to arbitrary random networks.

Future works may include a deeper study about the relationship between curvature and the geodesic distances between two random fields operating in different regimes, as a way to provide an intrinsic similarity measure. An analysis on how the components of the first and second fundamental form change as the inverse temperature varies may reveal relevant information about the underlying geometric structure of the parametric space. Besides, we intend to investigate techniques for the estimation of the inverse temperature parameter in order to simulate a situation in which we do not have direct access to the real inverse temperature value. Another idea consists in computing information-theoretic divergences, such as the KL-divergence, between pairs of random field models. Furthermore, we intend to study the feasibility of the application of the proposed information-geometric analysis in other random field models, such as the classic binary Ising model and the q-state Potts models, in which each variable assumes one of q different discrete states. Moreover, we intend to derive the Riemann curvature tensor for Gaussian random field manifolds as a way to completely characterize the curvature in terms of intrinsic measures: the Christoffel symbols of the metric tensor (first fundamental form). Finally, we intend to compute approximations for the geodesic distances between GMRFs and develop machine learning applications to perform unsupervised metric learning in classification problems.

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES)-Finance Code 001.

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/alexandrelevada/Curvature_GMRF.