The geodesic dispersion phenomenon in

random fields dynamics

Alexandre L. M. Levada1

Computing Department, Federal University of São Carlos, São Carlos, SP, Brazil

<EMAIL_ADDRESS>
###### Abstract

Random fields are ubiquitous mathematical structures in physics, with
applications ranging from thermodynamics and statistical physics to quantum
field theory and cosmology. Recent works on information geometry of Gaussian
random fields proposed mathematical expressions for the components of the
metric tensor of the underlying parametric space, allowing the computation of
the Gaussian curvature in each point of the manifold that represents the space
of all possible parameter values that define such mathematical model. A key
result in the dynamics of these random fields concerns the curvature effect, a
series of variations in the curvature that happens in the parametric space
when there are significant increase/decrease in the inverse temperature
parameter. In this paper, we propose a numerical algorithm for the computation
of geodesic curves in the Gaussian random fields manifold by deriving the 27
Christoffel symbols of the metric required for the definition of the Euler-
Lagrange equations. The fourth-order Runge-Kutta method is applied to solve
the Euler-Lagrange equations using an iterative approach based in Markov Chain
Monte Carlo simulation. Our results reveal that, when the system undergoes
phase trasitions, the geodesic dispersion phenomenon emerges: the geodesic
curve obtained by reversing the system of differential equations in time,
diverges from the original geodesic curve, as we move from zero curvature
configurations (Euclidean geometry) to negative curvature configurations
(hyperbolic-like geometry), and vice-versa. This phenomenon suggest that, time
irreversibility in random field dynamics can be a direct consequence of the
geometry of the underlying parametric space.

##  1 Introduction

The exploration of the fundamental principles governing the behavior of
physical systems has been a longstanding pursuit in several branches of
theoretical and applied physics [1]. In this pursuit, the incorporation of
randomness and disorder has emerged as a pivotal avenue for understanding the
intricate dynamics exhibited by a diverse array of systems [2, 3]. Among the
mathematical frameworks employed to model and analyze such inherent
uncertainties, random fields have proven to be indispensable [4]. Examples of
scenarios where random fields have profound implications in physics are
statistical physics [5], condensed matter physics [6], phase transitions [7],
critical phenomena [8], quantum systems [9], non-equilibrium statistical
mechanics [10], cosmology [11] and astrophysics [12].

Information geometry is a branch of mathematics that explores the geometric
structures associated with statistical models and probability distributions,
by modeling the parametric space of random variables as a Riemannian manifold
[13, 14, 15]. It provides a powerful framework for analyzing the relationships
and structures within these manifolds, offering deep insights into the
intrinsic geometry of random fields [16]. In particular, information geometry
provides a formal way to study random fields dynamics, offering insights into
local geometric properties of the parametric space as the system ehxibits
critical phenomena and phase transitions [17].

The concept of irreversible dynamics in random fields refers to the phenomenon
where the evolution of the field over time exhibits a lack of symmetry or
reversibility, meaning that the system does not return to its initial state
when allowed to evolve backward [18]. Irreversible dynamics often relate to an
increase in entropy over time. Entropy is a measure of disorder or randomness
in a system. As a random field evolves, the system tends to explore a wider
range of configurations, leading to increased entropy. The tendency for
entropy to increase is a key aspect of irreversibility in many physical
systems [19]. In the context of thermodynamics, irreversible dynamics can be
associated with processes that involve dissipation, such as heat transfer,
friction, or irreversible chemical reactions. These processes contribute to
the irreversibility of the overall system, and the random field may reflect
these irreversible thermodynamic behaviors. Phase transitions in random
fields, such as those occurring in condensed matter systems, can contribute to
irreversible dynamics. The system may undergo abrupt changes in its state, and
once the system has transitioned to a new phase, it might not easily revert to
its previous state [20].

Recent studies about information geomerty have proposed expressions for the
metric tensor of the Gaussian random fields manifold, which allows the
computation of the curvatures of the space along random field dynamics [16].
Computational simulations revealed the curvature effect, which means that a
highly asymmetric geometric deformation happens in the underlying parametric
space when there are significant increase/decrease in the inverse temperature
parameter, leading to phase transitions. When entropy is increasing, the
curvature transitions from negative to positive in a smoother way than in
comparison to the transition from positive to negative when entropy is
decreasing [17].

Motivated by this findings, this paper presents a computational method to
build geodesic curves in the manifold of Gaussian random fields, by
numerically solving the Euler-Lagrange equations that govern how a point moves
along this parametric space, given an initial location and an initial tangent
vector (velocity). Basically, the main contributions of this paper are: 1) to
derive the Christoffel symbols of the metric tensor of the Gaussian random
fields manifold [21, 22, 23]; 2) to build the Euler-Lagrange equations, which
defines a system of nonlinear differential equations that govern the dynamics
of a moving reference along geodesic curves in the underlying manifold [24,
25, 26]; and 3) to solve the Euler Lagrange equations using the fourth-order
Runge-Kutta [27, 28, 29] method and Markov-Chain Monte Carlo simulation [30,
31, 32]. Computational experiments reveal an interesting behavior, that we
called the geodesic dispersion phenomenon. In several cases, the geodesic
curve obtained by time reversing the simulations diverge from the original
geodesic curve, which suggests that, when curvature decreases from zero to
negative values, the emergence of a hyperbolic-like geometry makes a given
geodesic curve to deviate from its original trajectory.

The remaining of the paper is organized as follows: in Section 2, we introduce
the random field model and in Section 3, a brief overview of Markov Chain
Monte Carlo simulations is presented. In Section 4, we show the components of
the metric tensor for the random fields manifold and in Section 5, we discuss
geodesic distances in Riemannian manifolds. In Section 6, we briefly describe
the fourth-order Runge-Kutta (RK4) method and in Section 7, we derive the
Christoffel symbols of the metric. Section 8 shows the proposed iterative MCMC
based RK4 method for generating goedesic curves and in Section 9, we present
the computational experiments and results. Finally, Section 10 presents our
conclusions and final remarks.

##  2 Gaussian random fields

Gaussian random fields (GRFs) are stochastic processes that exhibit Gaussian
(normal) distribution at every spatial location. These fields are
characterized by the property that any finite collection of random variables,
corresponding to different locations in the field, follows a multivariate
normal distribution. Gaussian random fields are extensively used in various
scientific disciplines, including physics, due to their mathematical
tractability and versatility in modeling complex spatial structures and
uncertainties.

GRFs are models in which the set of possible states for a given variable is
continuous and innumerable. In practice, each variable in the field can assume
an infinite number of states [33]. We assume some simplifying hypothesis due
to mathematical tractability. First, we model only binary interactions in a 2D
lattice, constraining the analysis to pairwise interaction random field
models. Additionally, the inverse temperature parameter, which controls the
spatial dependence structure, is invariant and isotropic along the random
field. Finally, by invoking the Hammersley-Clifford theorem [34], which states
the equivalence between Gibbs random fields (global models) and Markov random
fields (local models), we characterize our model by the set of local
conditional density functions:

$p\left(x_{i}|\eta_{i},\vec{\theta}\right)=\frac{1}{\sqrt{2\pi\sigma^{2}}}exp\left\\{-\frac{1}{2\sigma^{2}}\left[x_{i}-\mu-\beta\sum_{j\in\eta_{i}}\left(x_{j}-\mu\right)\right]^{2}\right\\}$
(1)

where $\eta_{i}$ denotes the the second-order neighborhood system, comprised
by the 8 nearest neighbors of $x_{i}$ in the 2D lattice,
$\vec{\theta}=(\mu,\sigma^{2},\beta)$ denotes the parameters vector, with
$\mu$ and $\sigma^{2}$ being, the expected value (mean) and the variance of
the random variables in the lattice, and $\beta$ the inverse temperature. Note
that if $\beta=0$, the model degenerates to a regular Gaussian density and the
variables in the field become independent. Figure 1 shows the first, second
and third order neighborhood systems defined on a 2D lattice.

Fig. 1: First, second and third order neighborhood systems on a 2D lattice.

In summary, Gaussian random fields find applications in various branches of
physics due to their ability to model spatial variability and uncertainty in
physical processes.

##  3 Markov Chain Monte Carlo simulation

Markov Chain Monte Carlo (MCMC) algorithms are a class of computational
methods used for sampling from complex probability distributions, such as
random field models. These algorithms are particularly valuable when direct
sampling from the distribution is difficult or impractical. The fundamental
idea behind MCMC is to construct a Markov chain that has the desired
probability distribution as its equilibrium distribution. Once the Markov
chain reaches its equilibrium, samples from it can be considered as samples
from the target distribution [35, 36, 37].

A Markov chain is a mathematical model that describes a sequence of events,
where the probability of transitioning from one state to another depends only
on the current state, not on the sequence of events that preceded it. In the
context of MCMC, the states represent different points in the parameter space
of the distribution being sampled. Monte Carlo methods involve the use of
random sampling to obtain numerical results. In the case of MCMC, Monte Carlo
methods are employed to generate a sequence of samples from the target
distribution [38, 39, 40].

MCMC algorithms construct a Markov chain that satisfies the detailed balance
condition, ensuring that the desired distribution is its equilibrium
distribution. This involves defining a transition probability from one state
to another in such a way that the resulting stationary distribution is the
target distribution. One of the most commonly used MCMC algorithms is the
Metropolis-Hastings algorithm [30]. In each iteration of this algorithm, a
proposed new state is generated, and a decision is made to accept or reject
this proposal based on a acceptance probability. The acceptance probability is
designed to satisfy the detailed balance condition, ensuring convergence to
the target distribution. Algorithm 1 shows the pseudocode for the Metropolis-
Hastings method, used in all experiments along this paper.

Algorithm 1 Metropolis-Hastings algorithm for Markov Chain Monte Carlo
simulation

1:function Metropolis($n$, $m$, $N$, $\mu$, $\sigma^{2}$, $\beta$)

2: Generate a $n\times m$ matrix $X$ with samples from a Gaussian density with
mean $\mu$ and variance $\sigma^{2}$.

3: while $i\leq N$ do

4: for $i=1;i\leq n;i++$ do

5: for $i=j;j\leq m;j++$ do

6: $y_{ij}=Normal(\mu,\sigma^{2})$

7: Let $\displaystyle
P=min\left\\{1,\frac{p(y_{ij}|\eta_{i},\vec{\theta})}{p(x_{ij}|\eta_{i},\vec{\theta})}\right\\}$

8: Replace $x_{ij}$ by $y_{ij}$ with probability $p$

9: end for

10: end for

11: end while

12:end function

Recent works in applied probability and computational statistics link known
results about the diffusion limits of Markov chain Monte Carlo (MCMC)
algorithms to the concept of algorithm complexity in computer science. It has
been shown that, under appropriate assumptions, the random-walk in Metropolis-
adjusted algorithms takes $O(d^{1/3})$ iterations to converge to stationarity
in $d$ dimensions [41].

##  4 The metric tensor of the Gaussian random fields manifold

It has been shown by information geometry that the Fisher information matrix
is the metric tensor that equips the underlying parametric space of a
statistical manifold [42, 43]. In practical terms, the metric tensor makes it
possible to express the square of an infinitesimal displacement in the
manifold, $ds^{2}$, as a function of an infinitesimal displacement in the
tangent space. In order to define the Euler-Lagrange equations and compute
geodesic distances in the Gaussian random fields manifold, we need to compute
the metric tensor of the underlying parametric space. It has been shown in
previous works that the metric tensor of the Gaussian random fields manifold
has the following shape [16, 17]:

$g(\vec{\theta})=\left(\begin{array}[]{ccc}g_{11}&g_{12}&g_{13}\\\
g_{21}&g_{22}&g_{23}\\\ g_{31}&g_{32}&g_{33}\end{array}\right)$ (2)

where $g_{ij}$ is the Fisher information regarding parameters $\theta_{i}$ and
$\theta_{j}$:

$g_{ij}=E\left[\left(\frac{\partial}{\partial\theta_{i}}log~{}p(x;\vec{\theta})\right)\left(\frac{\partial}{\partial\theta_{j}}log~{}p(x;\vec{\theta})\right)\right]$
(3)

in which $p(x;\vec{\theta})$ denotes the local conditional density function of
the random field model. The closed-form expressions for all components of the
metric tensor are listed below [16, 17]:

$\displaystyle g_{11}$
$\displaystyle=\frac{\left(1-\beta\Delta\right)^{2}}{\sigma^{2}}\left[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)\right]$
(4) $\displaystyle g_{12}$ $\displaystyle=g_{13}=g_{21}=0$ (5) $\displaystyle
g_{22}$
$\displaystyle=\frac{1}{2\sigma^{4}}-\frac{1}{\sigma^{6}}\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]$
(6)
$\displaystyle+\frac{1}{\sigma^{8}}\left[3\beta^{2}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sigma_{ij}\sigma_{ik}-\beta^{3}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sum_{l\in\eta_{i}}\left(\sigma_{ij}\sigma_{kl}+\sigma_{ik}\sigma_{jl}+\sigma_{il}\sigma_{jk}\right)\right.$
$\displaystyle\hskip
28.45274pt\left.+\beta^{4}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sum_{l\in\eta_{i}}\sum_{m\in\eta_{i}}\left(\sigma_{jk}\sigma_{lm}+\sigma_{jl}\sigma_{km}+\sigma_{jm}\sigma_{kl}\right)\right]$
$\displaystyle g_{23}$
$\displaystyle=g_{32}=\frac{1}{\sigma^{4}}\left[\sum_{j\in\eta_{i}}\sigma_{ij}-\beta\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sigma_{jk}\right]$
(7)
$\displaystyle-\frac{1}{2\sigma^{6}}\left[6\beta\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sigma_{ij}\sigma_{ik}-3\beta^{2}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sum_{l\in\eta_{i}}\left(\sigma_{ij}\sigma_{kl}+\sigma_{ik}\sigma_{jl}+\sigma_{il}\sigma_{jk}\right)\right.$
$\displaystyle\hskip
85.35826pt\left.+\beta^{3}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sum_{l\in\eta_{i}}\sum_{m\in\eta_{i}}\left(\sigma_{jk}\sigma_{lm}+\sigma_{jl}\sigma_{km}+\sigma_{jm}\sigma_{kl}\right)\right]$
$\displaystyle g_{31}$ $\displaystyle=0$ (8) $\displaystyle g_{33}$
$\displaystyle=\frac{1}{\sigma^{2}}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sigma_{jk}+\frac{1}{\sigma^{4}}\left[2\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sigma_{ij}\sigma_{ik}-2\beta\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sum_{l\in\eta_{i}}\left(\sigma_{ij}\sigma_{kl}+\sigma_{ik}\sigma_{jl}+\sigma_{il}\sigma_{jk}\right)\right.$
(9)
$\displaystyle\left.+\beta^{2}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sum_{l\in\eta_{i}}\sum_{m\in\eta_{i}}\left(\sigma_{jk}\sigma_{lm}+\sigma_{jl}\sigma_{km}+\sigma_{jm}\sigma_{kl}\right)\right]$

where $\Delta$ is the cardinality of the neighborhood system ($\Delta=8$ in a
second-order system), $\sigma_{ij}$ is the covariance between the central
variable $x_{i}$ and one of its neighbors $x_{j}\in\eta_{i}$ and $\sigma_{jk}$
is the covariance between two variables $x_{j}$ and $x_{k}$ belonging to the
neighborhood $\eta_{i}$.

In order to reduce the computational cost in the numerical simulations, it is
possible to express the components of the metric tensor using Kronecker
products (tensor products). First, note that we can convert each $3\times 3$
neighborhood patch formed by $x_{i}\cup\eta_{i}$ into a vector $p_{i}$ of 9
elements by piling its rows. Then, we compute the covariance matrix of these
vectors, for $i=1,2,...,n$ denoted by $\Sigma_{p}$. From this covariance
matrix, we extract two main components: 1) a vector of size 8, $\vec{\rho}$,
composed by the the elements of the central row of $\Sigma_{p}$, excluding the
middle one, which denotes the variance of $x_{i}$ (we want only the
covariances between $x_{i}$ and $x_{j}$, for $j\neq i$; and 2) a sub-matrix of
dimensions $8\times 8$, $\Sigma_{p}^{-}$, obtained by removing the central row
and central column from $\Sigma_{p}$ (we want only the covariances between
$x_{j}\in\eta_{i}$ and $x_{k}\in\eta_{i}$). Figure 2 shows the decomposition
of the covariance matrix $\Sigma_{p}$ into the sub-matrix $\Sigma_{p}^{-}$ and
the vector $\vec{\rho}$. By employing Kronecker products, we rewrite the first
fundamental form (metric tensor) in a tensorial notation, providing a
computationally efficient way to compute the elements of $g(\vec{\theta})$:

Fig. 2: Decomposition of $\Sigma_{p}$ into $\Sigma_{p}^{-}$ and $\vec{\rho}$
on a second-order neighborhood system ($\Delta=8$). By rewriting the
components of the metric tensor in terms of Kronocker products, we can make
numerical simulations faster.

$g_{11}=\frac{1}{\sigma^{2}}\left(1-\beta\Delta\right)^{2}\left[1-\frac{1}{\sigma^{2}}\left(2\beta\left\|\vec{\rho}\right\|_{+}-\beta^{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right)\right]$
(10)

$\displaystyle g_{22}=\frac{1}{2\sigma^{4}}$
$\displaystyle-\frac{1}{\sigma^{6}}\left[2\beta\left\|\vec{\rho}\right\|_{+}-\beta^{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right]+\frac{1}{\sigma^{8}}\left[3\beta^{2}\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-3\beta^{3}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{4}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right]$
(11)

$\displaystyle
g_{23}=g_{32}=\frac{1}{\sigma^{4}}\left[\left\|\vec{\rho}\right\|_{+}-\beta\left\|\Sigma_{p}^{-}\right\|_{+}\right]-\frac{1}{2\sigma^{6}}\left[6\beta\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-9\beta^{2}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{3}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right]$
(12)

$g_{33}=\frac{1}{\sigma^{2}}\left\|\Sigma_{p}^{-}\right\|_{+}+\frac{1}{\sigma^{4}}\left[2\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-6\beta\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{2}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right]$
(13)

where $\left\|A\right\|_{+}$ represents the summation of all the entries of
the vector/matrix $A$ and $\otimes$ denotes the Kronecker (tensor) product.
The advantage of this tensorial notation using Kronecker products is that we
avoid the computation of triple and quadruple summations, which have
complexity $O(n^{3})$ and $O(n^{4})$, respectively.

###  4.1 Entropy in Gaussian random fields

The concept of entropy is crucial in many areas of science. In summary,
entropy is a thermodynamic quantity that measures the degree of molecular
freedom of a system, being associated with its number of configurations (or
microstates), that is, how many ways particles can be distributed in quantized
energy levels. It is also generally associated the randomness, dispersion of
matter and energy, and “disorder” of a complex system. In random fields, it is
possible to compute the entropy for each possible observed configuration. The
entropy of a Gaussian random field can be computed by the expected value of
self-information, which leads to:

$\displaystyle
H_{\beta}(\vec{\theta})=-E\left[log~{}p\left(x_{i}|\eta_{i},\vec{\theta}\right)\right]=H_{G}(\vec{\theta})-\frac{\beta}{\sigma^{2}}\sum_{j\in\eta_{i}}\sigma_{ij}+\frac{\beta^{2}}{2\sigma^{2}}\sum_{j\in\eta_{i}}\sum_{k\in\eta_{i}}\sigma_{jk}$
(14)

where $H_{G}(\vec{\theta})$ denotes the entropy of a Gaussian random variable
with mean $\mu$ and variance $\sigma^{2}$. Note that the entropy is a
quadratic function of the inverse temperature parameter $\beta$. Besides, for
$\beta=0$, we have $H_{\beta}(\vec{\theta})=H_{G}(\vec{\theta})$, as expected.
Using the tensorial notation, the entropy can be expressed as:

$\displaystyle
H_{\beta}(\vec{\theta})=H_{G}(\vec{\theta})-\frac{1}{\sigma^{2}}\left(\beta\left\|\vec{\rho}\right\|_{+}-\frac{\beta^{2}}{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right)$
(15)

Note that $H_{G}(\vec{\theta})$ does not depend on the inverse temperature
parameter, so differetiating $H_{\beta}(\vec{\theta})$ with respect to $\beta$
leads to:

$\displaystyle\frac{\partial}{\partial\beta}H_{\beta}$
$\displaystyle=-\frac{1}{\sigma^{2}}\left(\left\|\vec{\rho}\right\|_{+}-\beta\left\|\Sigma_{p}^{-}\right\|_{+}\right)$
(16)

which means that the second derivative of entropy with respect to $\beta$ is:

$\frac{\partial^{2}}{\partial\beta^{2}}H_{\beta}=\frac{1}{\sigma^{2}}\left\|\Sigma_{p}^{-}\right\|_{+}$
(17)

Hence, the components of the metric tensor can be expressed in terms of the
system’s entropy as:

$\displaystyle g_{11}$
$\displaystyle=\frac{1}{\sigma^{2}}\left(1-\beta\Delta\right)^{2}\left[1+2(H_{\beta}-H_{G})\right]$
(18) $\displaystyle g_{22}$
$\displaystyle=\frac{1}{2\sigma^{4}}+\frac{1}{\sigma^{4}}\left[2(H_{\beta}-H_{G})\right]+\frac{1}{\sigma^{8}}\left[3\beta^{2}\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-3\beta^{3}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{4}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right]$
(19) $\displaystyle g_{23}$
$\displaystyle=g_{32}=-\frac{1}{\sigma^{2}}\frac{\partial}{\partial\beta}H_{\beta}-\frac{1}{2\sigma^{6}}\left[6\beta\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-9\beta^{2}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{3}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right]$
$\displaystyle g_{33}$
$\displaystyle=\frac{\partial^{2}}{\partial\beta^{2}}H_{\beta}+\frac{1}{\sigma^{4}}\left[2\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-6\beta\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{2}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right]$
(20)

showing that variations in the system’s entropy directly affect the metric
tensor.

##  5 Geodesic distances in Riemannian manifolds

In a Riemannian manifold $M$ equipped with a metric tensor $g$, the length of
a curve $\gamma(t)$, defined for $t\in[a,b]$ is:

$L(\gamma)=\int_{a}^{b}\sqrt{g(\gamma^{\prime}(t),\gamma^{\prime}(t))}dt$ (21)

where $\gamma^{\prime}(t)$ denotes the first derivative (tangent vector). It
is easier to work with the energy of the curve. The energy of a curve
$\gamma(t)$ is given by:

$E(\gamma)=\frac{1}{2}\int_{a}^{b}\sqrt{g(\gamma^{\prime}(t),\gamma^{\prime}(t))}dt$
(22)

It is not difficult to show that if $M$ is a manifold equipped with a metric
tensor $g(\theta)$, then:

$L(\gamma)\leq 2(b-a)E(\gamma)$ (23)

with equality if and only if $\gamma$ is parametrized by arc length, that is,
if we move along the curve with constant velocity (the tangent vectors are
unitary). In summary, a curve $\gamma(t)$ has minimum energy if and only if it
has minimum length and arc length parametrization. A geodesic is a curve
$\gamma(t)$ that minimizes the functional $L(\gamma)$. In order to obtain an
exact solution, we need to solve the following optimization problem [44]:

$\nabla_{\gamma^{\prime}}\gamma^{\prime}=0$ (24)

where $\nabla$ is the covariant derivative. The geodesics are precisely the
“straight lines” that have constant velocity. From differential geometry, it
has been shown that a solution to this equation must satisfy the geodesic
equations, expressed by [45]:

$\frac{d^{2}\theta_{k}}{dt^{2}}+\sum_{ij}\Gamma_{ij}^{k}\frac{d\theta_{i}}{dt}\frac{d\theta_{j}}{dt}=0$
(25)

where, in our pairwise isotropic Gaussian-Markov random field,
$\theta_{1}=\mu$, $\theta_{2}=\sigma^{2}$ and $\theta_{3}=\beta$ are the model
parameters and $\Gamma_{ij}^{k}$ are the Christoffel symbols of the metric,
defined as [46]:

$\Gamma_{ij}^{k}=\frac{1}{2}\sum_{m=1}^{3}\left(\frac{\partial}{\partial\theta_{i}}g_{jm}+\frac{\partial}{\partial\theta_{j}}g_{im}-\frac{\partial}{\partial\theta_{m}}g_{ij}\right)g^{mk}$
(26)

where $g_{ij}$ a component of the metric tensor (Fisher information matrix)
and $g^{ij}$ is a component of the inverse metric tensor. As the parametric
space of pairwise isotropic Gaussian-Markov random fields are 3D manifolds, we
have a total of $3\times 3\times 3=27$ Christoffel symbols arranged in three
$3\times 3$ matrices $\Gamma^{1}$, $\Gamma^{2}$ and $\Gamma^{3}$, leading to
the following system of non-linear differential equations:

$\displaystyle\theta_{1}^{\prime\prime}(t)-\sum_{ij}\Gamma_{ij}^{1}\theta_{i}^{\prime}(t)\theta_{j}^{\prime}(t)=0$
(27)
$\displaystyle\theta_{2}^{\prime\prime}(t)-\sum_{ij}\Gamma_{ij}^{2}\theta_{i}^{\prime}(t)\theta_{j}^{\prime}(t)=0$
(28)
$\displaystyle\theta_{3}^{\prime\prime}(t)-\sum_{ij}\Gamma_{ij}^{3}\theta_{i}^{\prime}(t)\theta_{j}^{\prime}(t)=0$
(29)

where $\theta_{i}^{\prime}(t)$ and $\theta_{i}^{\prime\prime}(t)$ denotes the
first and second derivatives of the coordinates, the values of
$\theta_{1}^{(0)}=\mu_{0}$, $\theta_{2}^{(0)}=\sigma_{0}^{2}$ and
$\theta_{3}^{(0)}=\beta_{0}$ represent the starting point A, and the values of
$\theta_{1}^{\prime(0)}$, $\theta_{2}^{\prime(0)}$ and
$\theta_{3}^{\prime(0)}$ defines the initial tangent vector at A (it points to
the initial direction that we want move along the space). The fourth-order
Runge-Kutta method can be applied to obtain a numerical solution for this
system of differential equations, where the final value of the parameter $t$
(number of iterations) determines how long the geodesic curve will grow. Note
that there are no explicit values for the final parameters
$\theta_{1}^{(t)}=\mu_{t}$, $\theta_{2}^{(t)}=\sigma_{t}^{2}$ and
$\theta_{3}^{(t)}=\beta_{t}$, since the manifold is a curved space and the
geodesic is build locally, using one tangent space at a time.

##  6 The fourth-order Runge-Kutta method

The Runge-Kutta method is a numerical technique used for solving ordinary
differential equations (ODEs). Its importance lies in its ability to provide
accurate and efficient solutions to a wide range of differential equations
encountered in various scientific and engineering applications [47, 48]. The
main advantages of the Runge-Kutta method over other approaches are [49, 50,
51]:

  * •

Accuracy: The Runge-Kutta method is known for its higher order accuracy
compared to simpler methods like Euler’s method. By employing multiple stages
in its computations, such as the classic fourth-order Runge-Kutta (RK4), it
achieves a balance between accuracy and computational efficiency. This method
is particularly valuable for solving stiff differential equations where other
numerical methods might struggle due to rapid changes in the solution or
disparate time scales. The stability of Runge-Kutta methods is crucial in
ensuring reliable results across a broad spectrum of problems.

  * •

Applicability to Systems of ODEs: Many real-world problems involve systems of
interconnected ODEs. The Runge-Kutta method naturally extends to handle
systems of differential equations, making it a versatile tool for modeling
complex dynamic systems in various scientific domains. It is widely employed
in physics, biology, engineering, and other fields to simulate and analyze the
behavior of interconnected systems over time.

  * •

Adaptive Step Size: The Runge-Kutta method supports adaptive step size
control. This feature is essential when dealing with ODEs where the solution
undergoes rapid changes in certain regions and remains relatively constant in
others. Adaptive step size algorithms help optimize the trade-off between
accuracy and computational cost by adjusting the step size dynamically during
the integration process.

  * •

Numerical Stability and Robustness: Runge-Kutta methods, especially higher-
order variants, exhibit good numerical stability characteristics. This
property is crucial for preventing numerical instabilities that may arise
during the solution of differential equations, ensuring that the computed
results are reliable and physically meaningful.

  * •

Versatility in Scientific Computing: Due to its robustness and versatility,
the Runge-Kutta method has become a standard numerical tool in scientific
computing. It is employed in diverse fields, including physics, chemistry,
biology, and engineering, for simulating and understanding complex dynamical
systems.

In summary, the Runge-Kutta method is a fundamental numerical technique with
widespread applications in scientific research, especially in physics. Its
accuracy, stability, and versatility make it an indispensable tool for solving
systems of ordinary differential equations and gaining insights into the
behavior of random fields dynamics. The most popular version of Runge–Kutta
method is the fourth-order algorithm, often denoted by RK4, the classic
Runge–Kutta method, which is described in the following.

Consider a first-order ordinary differential equation in the form:

$\frac{dy}{dt}=f(t,y)$ (30)

where $y$ is the unknown function of $t$ and $f(t,y)$ is a given function
representing the rate of change of $y$ with respect to $t$. At the initial
time $t_{0}$, the corresponding $y$ value is $y_{0}$. It should be mentioned
that the function $f$ and the initial conditions $t_{0},y_{0}$ must be given
as input to the problem. The fourth-order Runge-Kutta algorithm is given by
the following steps [29]:

  1. 1.

Define $n=0$ and the step size $h>0$.

  2. 2.

Compute the four increments:

$\displaystyle k_{1}$ $\displaystyle=hf(t_{n},y_{n})$ (31) $\displaystyle
k_{2}$ $\displaystyle=hf\left(t_{n}+\frac{h}{2},y_{n}+\frac{k_{1}}{2}\right)$
(32) $\displaystyle k_{3}$
$\displaystyle=hf\left(t_{n}+\frac{h}{2},y_{n}+\frac{k_{2}}{2}\right)$ (33)
$\displaystyle k_{4}$ $\displaystyle=hf\left(t_{n}+h,y_{n}+k_{3}\right)$ (34)

  3. 3.

Update the value of the function:

$\displaystyle y_{n+1}$
$\displaystyle=y_{n}+\frac{1}{6}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)$ (35)
$\displaystyle t_{n+1}$ $\displaystyle=t_{n}+h$ (36)

  4. 4.

Repeat steps 2 and 3 for $n=1,2,3,...,N$, where $N$ is the number of
iterations.

Basically, the fourth-order Runge-Kutta algorithm works by iteratively
updating the solution of an ordinary differential equation based on weighted
averages of incremental steps, providing an accurate and efficient numerical
approximation of the solution over a specified range.

##  7 The Christoffel symbols of the metric

In order to compute the Christoffel symbols of the underlying parametric space
of pairwise isotropic Gaussian-Markov random field models, we need to find the
components of the inverse metric tensor $g^{-1}(\theta)$:

$g^{-1}(\vec{\theta})=\left(\begin{array}[]{ccc}g^{11}&0&0\\\
0&g^{22}&g^{23}\\\ 0&g^{32}&g^{33}\end{array}\right)$ (37)

It is straightforward to check that the components of the inverse metric
tensor are:

$\displaystyle g^{11}=\frac{1}{g_{11}}\qquad
g^{22}=\frac{g_{33}}{g_{22}g_{33}-g_{23}^{2}}\qquad
g^{23}=\frac{g_{23}}{g_{23}^{2}-g_{22}g_{33}}\qquad
g^{32}=\frac{g_{32}}{g_{32}^{2}-g_{22}g_{33}}\qquad
g^{33}=\frac{g_{22}}{g_{22}g_{33}-g_{32}^{2}}$ (38)

where $g^{23}=g^{32}$. In the computational experiments, we perform
regularization to avoid numerical instability by adding a small value
$\lambda=0.001$ to the main diagonal of the metric tensor. The next step
consists in the derivation of the components of the Christoffel symbols.
First, note that:

$\frac{\partial}{\partial\theta_{1}}g_{ij}=\frac{\partial}{\partial\mu}g_{ij}=0,\qquad\qquad\forall
i,j$ (39)

since none of the components of the metric tensor is a function of
$\theta_{1}=\mu$. Moving to the derivatives with respect to
$\theta_{2}=\sigma^{2}$, we have:

$\displaystyle\frac{\partial}{\partial\theta_{2}}g_{11}$
$\displaystyle=-\frac{1}{\sigma^{4}}(1-\beta\Delta)^{2}+\frac{2}{\sigma^{6}}(1-\beta\Delta)^{2}\left(2\beta\left\|\vec{\rho}\right\|_{+}-\beta^{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right)$
(40) $\displaystyle\frac{\partial}{\partial\theta_{2}}g_{22}$
$\displaystyle=-\frac{1}{\sigma^{6}}+\frac{3}{\sigma^{8}}\left(2\beta\left\|\vec{\rho}\right\|_{+}-\beta^{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right)-\frac{4}{\sigma^{10}}\left(3\beta^{2}\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-3\beta^{3}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{4}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right)$
(41) $\displaystyle\frac{\partial}{\partial\theta_{2}}g_{23}$
$\displaystyle=\frac{\partial}{\partial\theta_{2}}g_{32}=-\frac{2}{\sigma^{6}}\left(\left\|\vec{\rho}\right\|_{+}-\beta\left\|\Sigma_{p}^{-}\right\|_{+}\right)+\frac{3}{2\sigma^{8}}\left(6\beta\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-9\beta^{2}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{3}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right)$
(42) $\displaystyle\frac{\partial}{\partial\theta_{2}}g_{33}$
$\displaystyle=-\frac{1}{\sigma^{4}}\left\|\Sigma_{p}^{-}\right\|_{+}-\frac{2}{\sigma^{6}}\left(2\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-6\beta\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+3\beta^{2}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right)$
(43)

Finally, we have to differentiate the components of the metric tensor with
respect to the inverse temperature parameter, that is, $\theta_{3}=\beta$:

$\displaystyle\frac{\partial}{\partial\theta_{3}}g_{11}$
$\displaystyle=-\frac{1}{\sigma^{2}}2\Delta(1-\beta\Delta)\left[1-\frac{1}{\sigma^{2}}\left(2\beta\left\|\vec{\rho}\right\|_{+}-\beta^{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right)\right]-\frac{1}{\sigma^{4}}(1-\beta\Delta)^{2}\left(2\left\|\vec{\rho}\right\|_{+}-2\beta\left\|\Sigma_{p}^{-}\right\|_{+}\right)$
(44) $\displaystyle\frac{\partial}{\partial\theta_{3}}g_{22}$
$\displaystyle=-\frac{1}{\sigma^{6}}\left(2\left\|\vec{\rho}\right\|_{+}-2\beta\left\|\Sigma_{p}^{-}\right\|_{+}\right)+\frac{1}{\sigma^{8}}\left(6\beta\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-9\beta^{2}\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+12\beta^{3}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right)$
(45)

$\displaystyle\frac{\partial}{\partial\theta_{3}}g_{23}$
$\displaystyle=\frac{\partial}{\partial\theta_{3}}g_{32}=-\frac{1}{\sigma^{4}}\left\|\Sigma_{p}^{-}\right\|_{+}-\frac{1}{2\sigma^{6}}\left(6\left\|\vec{\rho}\otimes\vec{\rho}\right\|_{+}-18\beta\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}+9\beta^{2}\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right)$
(46) $\displaystyle\frac{\partial}{\partial\theta_{3}}g_{33}$
$\displaystyle=-\frac{1}{\sigma^{4}}\left(6\beta\left\|\vec{\rho}\otimes\Sigma_{p}^{-}\right\|_{+}-6\beta\left\|\Sigma_{p}^{-}\otimes\Sigma_{p}^{-}\right\|_{+}\right)$
(47)

Given the above, it is possible to compute the 27 Christoffel symbols, denoted
by:
$\Gamma_{11}^{1},\Gamma_{11}^{2},\Gamma_{11}^{3},...,\Gamma_{33}^{1},\Gamma_{33}^{2},\Gamma_{33}^{3}$,
organized in a tensor composed by three $3\times 3$ matrices
$\Gamma^{1},\Gamma^{2},\Gamma^{3}$. However, it can be shown that the
Christoffel symbols are symmetric with respect to the lower indices, that is,
$\Gamma_{ij}^{k}=\Gamma_{ji}^{k}$, resulting in a total of 18 different
symbols. First, note that the components
$\Gamma_{11}^{1},\Gamma_{11}^{2},\Gamma_{11}^{3}$ are:

$\displaystyle\Gamma_{11}^{1}$ $\displaystyle=0$ (48)
$\displaystyle\Gamma_{11}^{2}$ $\displaystyle=-\frac{1}{2}\left(\frac{\partial
g_{11}}{\partial\theta_{2}}g^{22}+\frac{\partial
g_{11}}{\partial\theta_{3}}g^{32}\right)$ (49) $\displaystyle\Gamma_{11}^{3}$
$\displaystyle=-\frac{1}{2}\left(\frac{\partial
g_{11}}{\partial\theta_{2}}g^{23}+\frac{\partial
g_{11}}{\partial\theta_{3}}g^{33}\right)$ (50)

Next, the components $\Gamma_{12}^{1},\Gamma_{12}^{2},\Gamma_{12}^{3}$ can be
expressed by:

$\displaystyle\Gamma_{12}^{1}=\frac{1}{2}\frac{\partial
g_{11}}{\partial\theta_{2}}g^{11}=\Gamma_{21}^{1}\qquad\Gamma_{12}^{2}=\Gamma_{21}^{2}=0\qquad\Gamma_{12}^{3}=\Gamma_{21}^{3}=0$
(51)

Similarly, the components $\Gamma_{13}^{1},\Gamma_{13}^{2},\Gamma_{13}^{3}$
also have the same structure:

$\displaystyle\Gamma_{13}^{1}=\frac{1}{2}\frac{\partial
g_{11}}{\partial\theta_{3}}g^{11}=\Gamma_{31}^{1}\qquad\Gamma_{13}^{2}=\Gamma_{31}^{2}=0\qquad\Gamma_{13}^{3}=\Gamma_{31}^{3}=0$
(52)

Moving forward to the components
$\Gamma_{22}^{1},\Gamma_{22}^{2},\Gamma_{22}^{3}$, we have:

$\displaystyle\Gamma_{22}^{1}$ $\displaystyle=0$ (53)
$\displaystyle\Gamma_{22}^{2}$ $\displaystyle=\frac{1}{2}\left[\frac{\partial
g_{22}}{\partial\theta_{2}}g^{22}+\left(2\frac{\partial
g_{23}}{\partial\theta_{2}}-\frac{\partial
g_{22}}{\partial\theta_{3}}\right)g^{32}\right]$ (54)
$\displaystyle\Gamma_{22}^{3}$ $\displaystyle=\frac{1}{2}\left[\frac{\partial
g_{22}}{\partial\theta_{2}}g^{23}+\left(2\frac{\partial
g_{23}}{\partial\theta_{2}}-\frac{\partial
g_{22}}{\partial\theta_{3}}\right)g^{33}\right]$ (55)

The components $\Gamma_{23}^{1},\Gamma_{23}^{2},\Gamma_{23}^{3}$ are given by:

$\displaystyle\Gamma_{23}^{1}$ $\displaystyle=\Gamma_{32}^{1}=0$ (56)
$\displaystyle\Gamma_{23}^{2}$ $\displaystyle=\frac{1}{2}\left(\frac{\partial
g_{22}}{\partial\theta_{3}}g^{22}+\frac{\partial
g_{33}}{\partial\theta_{2}}g^{32}\right)=\Gamma_{32}^{2}$ (57)
$\displaystyle\Gamma_{23}^{3}$ $\displaystyle=\frac{1}{2}\left(\frac{\partial
g_{22}}{\partial\theta_{3}}g^{23}+\frac{\partial
g_{33}}{\partial\theta_{2}}g^{33}\right)=\Gamma_{32}^{3}$ (58)

Finally, the components $\Gamma_{33}^{1},\Gamma_{33}^{2},\Gamma_{33}^{3}$ can
be computed by:

$\displaystyle\Gamma_{33}^{1}$ $\displaystyle=0$ (59)
$\displaystyle\Gamma_{33}^{2}$
$\displaystyle=\frac{1}{2}\left[\left(2\frac{\partial
g_{32}}{\partial\theta_{3}}-\frac{\partial
g_{33}}{\partial\theta_{2}}\right)g^{22}+\frac{\partial
g_{33}}{\partial\theta_{3}}g^{32}\right]$ (60) $\displaystyle\Gamma_{33}^{3}$
$\displaystyle=\frac{1}{2}\left[\left(2\frac{\partial
g_{32}}{\partial\theta_{3}}-\frac{\partial
g_{33}}{\partial\theta_{2}}\right)g^{23}+\frac{\partial
g_{33}}{\partial\theta_{3}}g^{33}\right]$ (61)

Therefore, from the 27 Christoffel symbols, 13 are zero and the remaining 14
non-zero symbols assume 10 different values due to symmetry properties.

###  7.1 Approximating the curvature

In differential geometry, to compute a surface curvature, 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 [46]. Let $S$ be a surface
with first fundamental form $\mathbb{I}$ and second fundamental form
$\mathbb{II}$. Then, the shape operator $P$ is:

$P=-(\mathbb{II})(\mathbb{I})^{-1}=-\begin{bmatrix}L&M\\\ M&N\\\
\end{bmatrix}\begin{bmatrix}E&F\\\ F&G\\\ \end{bmatrix}^{-1}$ (62)

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 as
the determinant, the trace and the eigenvalues of the shape operator. In this
study, we compute the shape operator of the Gaussian random fields manifold to
estimate the curvatures. The metric tensor is the generalization of the first
fundamental form, and we compute the second fundamental form as:

$\left\\{\mathbb{II}(\vec{\theta})\right\\}_{ij}=-\mathbb{E}\left[\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}log~{}p(X;\vec{\theta})\right],\text{~{}~{}~{}~{}
for }i,j=1,\ldots,n$ (63)

It has been shown that, in pairwise isotropic Gaussian random fields, the
elements of the second fundamental form are given by [16, 17]:

$\displaystyle\mathbb{II}_{\mu\mu}(\vec{\theta})$
$\displaystyle=\frac{1}{\sigma^{2}}\left(1-\beta\Delta\right)^{2}$ (64)
$\displaystyle\mathbb{II}_{\sigma^{2}\sigma^{2}}(\vec{\theta})$
$\displaystyle=\frac{1}{2\sigma^{4}}-\frac{1}{\sigma^{6}}\left[2\beta\left\|\vec{\rho}\right\|_{+}-\beta^{2}\left\|\Sigma_{p}^{-}\right\|_{+}\right]=\frac{1}{2\sigma^{4}}+\frac{1}{\sigma^{4}}\left[2(H_{\beta}-H_{G})\right]$
(65) $\displaystyle\mathbb{II}_{\sigma^{2}\beta}(\vec{\theta})$
$\displaystyle=II_{\beta\sigma^{2}}(\vec{\theta})=\frac{1}{\sigma^{4}}\left[\left\|\vec{\rho}\right\|_{+}-\beta\left\|\Sigma_{p}^{-}\right\|_{+}\right]=-\frac{1}{\sigma^{2}}\frac{\partial}{\partial\beta}H_{\beta}$
(66) $\displaystyle\mathbb{II}_{\beta\beta}(\vec{\theta})$
$\displaystyle=\frac{1}{\sigma^{2}}\left\|\Sigma_{p}^{-}\right\|_{+}=\frac{\partial^{2}}{\partial\beta^{2}}H_{\beta}$
(67)

showing that they are closely related to the system’s entropy and its first
and second derivatives.

##  8 The proposed MCMC based RK4 simulation algorithm

In order to numerically solve our system of non-linear second-order
differential equations, first we have to convert it to a first-order system.
This can be done by a simple variable substitution. Let
$\gamma_{1}(t)=\theta_{1}(t)$ and $\alpha_{1}(t)=\theta_{1}^{\prime}(t)$.
Then, we have:

$\displaystyle\gamma_{1}^{\prime}(t)$
$\displaystyle=\theta_{1}^{\prime}(t)=\alpha_{1}(t)$ (68)
$\displaystyle\alpha_{1}^{\prime}(t)$
$\displaystyle=\theta_{1}^{\prime\prime}(t)=-\sum_{ij}\Gamma_{i,j}^{1}\alpha_{i}(t)\alpha_{j}(t)$
(69)

By direct application of this procedure to the other equations, we finally
reach our system of non-linear first-order differential equations:

$\displaystyle\gamma_{1}^{\prime}(t)$ $\displaystyle=\alpha_{1}(t)$ (70)
$\displaystyle\gamma_{2}^{\prime}(t)$ $\displaystyle=\alpha_{2}(t)$ (71)
$\displaystyle\gamma_{3}^{\prime}(t)$ $\displaystyle=\alpha_{3}(t)$ (72)
$\displaystyle\alpha_{1}^{\prime}(t)$
$\displaystyle=-\sum_{ij}\Gamma_{i,j}^{1}\alpha_{i}(t)\alpha_{j}(t)$ (73)
$\displaystyle\alpha_{2}^{\prime}(t)$
$\displaystyle=-\sum_{ij}\Gamma_{i,j}^{2}\alpha_{i}(t)\alpha_{j}(t)$ (74)
$\displaystyle\alpha_{3}^{\prime}(t)$
$\displaystyle=-\sum_{ij}\Gamma_{i,j}^{3}\alpha_{i}(t)\alpha_{j}(t)$ (75)

Note that we can rewrite our system in a more compact way as:

$\displaystyle\gamma_{1}^{\prime}(t)$
$\displaystyle=F(t,\vec{\alpha})=\alpha_{1}$ (76)
$\displaystyle\gamma_{2}^{\prime}(t)$
$\displaystyle=G(t,\vec{\alpha})=\alpha_{2}$ (77)
$\displaystyle\gamma_{3}^{\prime}(t)$
$\displaystyle=H(t,\vec{\alpha})=\alpha_{3}$ (78)
$\displaystyle\alpha_{1}^{\prime}(t)$
$\displaystyle=P(t,\vec{\alpha},\Gamma_{1})=-\vec{\alpha}^{T}\Gamma^{1}\vec{\alpha}$
(79) $\displaystyle\alpha_{2}^{\prime}(t)$
$\displaystyle=Q(t,\vec{\alpha},\Gamma_{2})=-\vec{\alpha}^{T}\Gamma^{2}\vec{\alpha}$
(80) $\displaystyle\alpha_{3}^{\prime}(t)$
$\displaystyle=R(t,\vec{\alpha},\Gamma_{3})=-\vec{\alpha}^{T}\Gamma^{3}\vec{\alpha}$
(81)

where $\vec{\alpha}^{T}=[\alpha_{1},\alpha_{2},\alpha_{3}]$ is the tangent
vector and, for $k=1,2,3$:

$\Gamma^{k}=\left(\begin{array}[]{ccc}\Gamma_{11}^{k}&\Gamma_{12}^{k}&\Gamma_{13}^{k}\\\
\Gamma_{21}^{k}&\Gamma_{22}^{k}&\Gamma_{23}^{k}\\\
\Gamma_{31}^{k}&\Gamma_{32}^{k}&\Gamma_{33}^{k}\\\ \end{array}\right)$ (82)

Algorithm 2 Computing geodesic curves in random fields manifold with RK4 and
MCMC simulation

1:function GeodesicDistanceRandomField($a$, $b$, $n$, $N$,
$\vec{\gamma}^{(0)}$, $\vec{\alpha}^{(0)}$)

2: $h=(b-a)/n$

3: $t=a$

4: $L=0$ $\triangleright$ L denotes the geodesic distance

5: for $i=0$, $i<N$, $i++$ do

6: Generate an outcome of the GMRF model using $\vec{\gamma}^{(i)}$.

7: Compute the metric tensor $g(\vec{\gamma}^{(i)})$ with equations (10),
(11), (12) and (13).

8: Regulize the metric tensor adding a small $\epsilon$ in its main diagonal.

9: Compute the inverse metric tensor $g^{-1}(\vec{\gamma}^{(i)})$.

10: Compute the Christoffel symbols ($\Gamma_{1}$, $\Gamma_{2}$, $\Gamma_{3}$)
with equations (48) to (61).

11: // Perform a 4th order Runge Kutta iteration

12: $k_{0}=hF(t,\vec{\alpha})$

13: $l_{0}=hG(t,\vec{\alpha})$

14: $m_{0}=hH(t,\vec{\alpha})$

15: $x_{0}=hP(t,\vec{\alpha})$

16: $y_{0}=hQ(t,\vec{\alpha})$

17: $z_{0}=hR(t,\vec{\alpha})$

18: $\vec{\delta}=[k_{0},l_{0},m_{0}]$

19: $k_{1}=hF(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

20: $l_{1}=hG(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

21: $m_{1}=hH(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

22: $x_{1}=hP(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{1})$

23: $y_{1}=hQ(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{2})$

24: $z_{1}=hR(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{3})$

25: $\vec{\delta}=[k_{1},l_{1},m_{1}]$

26: $k_{2}=hF(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

27: $l_{2}=hG(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

28: $m_{2}=hH(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

29: $x_{2}=hP(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{1})$

30: $y_{2}=hQ(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{2})$

31: $z_{2}=hR(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{3})$

32: $\vec{\delta}=[k_{2},l_{2},m_{2}]$

33: $k_{3}=hF(t+h,\vec{\alpha}+\vec{\delta})$

34: $l_{3}=hG(t+h,\vec{\alpha}+\vec{\delta})$

35: $m_{3}=hH(t+h,\vec{\alpha}+\vec{\delta})$

36: $x_{3}=hP(t+h,\vec{\alpha}+\vec{\delta},\Gamma_{1})$

37: $y_{3}=hQ(t+h,\vec{\alpha}+\vec{\delta},\Gamma_{2})$

38: $z_{3}=hR(t+h,\vec{\alpha}+\vec{\delta},\Gamma_{3})$

39: Update the vectors $\vec{\gamma}^{(i)}$ (position) and
$\vec{\alpha}^{(i)}$ (velocity) as

40: $\gamma_{1}^{(i+1)}=\gamma_{1}^{(i)}+(1/6)(k_{0}+2k_{1}+2k_{2}+k_{3})$

41: $\gamma_{2}^{(i+1)}=\gamma_{2}^{(i)}+(1/6)(l_{0}+2l_{1}+2l_{2}+l_{3})$

42: $\gamma_{3}^{(i+1)}=\gamma_{3}^{(i)}+(1/6)(m_{0}+2m_{1}+2m_{2}+m_{3})$

43: $\alpha_{1}^{(i+1)}=\alpha_{1}^{(i)}+(1/6)(x_{0}+2x_{1}+2x_{2}+x_{3})$

44: $\alpha_{2}^{(i+1)}=\alpha_{2}^{(i)}+(1/6)(y_{0}+2y_{1}+2y_{2}+y_{3})$

45: $\alpha_{3}^{(i+1)}=\alpha_{3}^{(i)}+(1/6)(z_{0}+2z_{1}+2z_{2}+z_{3})$

46: Update the geodesic distance: $L=L+\lVert\alpha(t)\rVert h$

47: end for

48:end function

Fig. 3: The computation of the arc length can be approximated by the norm of
the difference vector (hypotenuse of a right triangle).

The basic idea is that the initial condition
$\vec{\gamma}^{(0)}=[\gamma_{1}^{(0)},\gamma_{2}^{(0)},\gamma_{3}^{(0)}]$
defines the starting point in the manifold and the initial condition
$\vec{\alpha}^{(0)}=[\alpha_{1}^{(0)},\alpha_{2}^{(0)},\alpha_{3}^{(0)}]$
defines the initial tangent vector (velocity). The arc length of the resulting
geodesic curve defines our geodesic distance. Figure 3 illustrates how the arc
length between $t=r$ and $t=s$ can be approximated by the norm of the vector
$\vec{s}$, which is the difference between $\vec{\gamma}(r)$ and
$\vec{\gamma}(s)$. If we increase the number of subdivisions of the interval
$[r,s]$, the better will be the approximation, given by:

$s(t)\approx\sum_{t=r}^{s}\sqrt{\left(x(t+\Delta
t)-x(t)\right)^{2}+\left(y(t+\Delta t)-y(t)\right)^{2}+\left(z(t+\Delta
t)-z(t)\right)^{2}}$ (83)

Note that we can rewrite $s(t)$ as:

$s(t)\approx\sum_{t=r}^{s}\sqrt{\left(\frac{x(t+\Delta t)-x(t)}{\Delta
t}\right)^{2}+\left(\frac{y(t+\Delta t)-y(t)}{\Delta
t}\right)^{2}+\left(\frac{z(t+\Delta t)-z(t)}{\Delta t}\right)^{2}}\Delta t$
(84)

From the definition of the derivative of a function, we have:

$f^{\prime}(t)=\lim\limits_{\Delta t\to 0}\frac{f(t+\Delta t)-f(t)}{\Delta t}$
(85)

Thus, in the limiting case, when $\Delta t\to 0$, the approximation becomes
exact:

$\displaystyle s(t)$
$\displaystyle=\int_{a}^{b}\sqrt{x^{\prime}(t)^{2}+y^{\prime}(t)^{2}+z^{\prime}(t)^{2}}dt$
(86)
$\displaystyle=\int_{a}^{b}\sqrt{\vec{\gamma}^{\prime}(t)^{T}\vec{\gamma}^{\prime}(t)}dt=\int_{a}^{b}\lVert\vec{\gamma}^{\prime}(t)\rVert
dt$

where the integrand $\lVert\vec{\gamma}^{\prime}(t)\rVert$ is the norm of the
tangent vector at $t$. Therefore, we compute the geodesic distance as:

$L(\gamma)=\int_{a}^{b}\lVert\gamma^{\prime}(t)\rVert
dt\approx\sum_{i=a}^{b}\lVert\gamma^{\prime}(i)\rVert h$ (87)

using a small enough step size $h=(b-a)/n$.

##  9 Computational experiments

In the first set of experiments, we use the fouth-order Runge-Kutta method to
build geodesic curves in the regular 2D Gaussian random variables manifold.
Information geometry in Gaussian random variables has been extensively studied
in the literature [52, 53]. It is not difficult to see that the metric tensor
of the Gaussian random variables manifold is:

$g(\theta)=\left(\begin{array}[]{cc}\displaystyle\frac{1}{\sigma^{2}}&\displaystyle
0\\\ \displaystyle 0&\displaystyle\frac{1}{2\sigma^{4}}\\\ \end{array}\right)$
(88)

Moreover, the derivation of the Christoffel symbols of the metric reveals
that:

$\Gamma^{1}=\left(\begin{array}[]{cc}\displaystyle
0&\displaystyle-\frac{1}{2\sigma^{2}}\\\
\displaystyle-\frac{1}{2\sigma^{2}}&\displaystyle
0\end{array}\right)\qquad\Gamma^{2}=\left(\begin{array}[]{cc}\displaystyle
1&\displaystyle 0\\\ \displaystyle
0&\displaystyle-\frac{1}{\sigma^{2}}\end{array}\right)$ (89)

Given the above, we define our system of equations as:

$\displaystyle\gamma_{1}^{\prime}(t)$
$\displaystyle=F(t,\vec{\alpha})=\alpha_{1}$ (90)
$\displaystyle\gamma_{2}^{\prime}(t)$
$\displaystyle=G(t,\vec{\alpha})=\alpha_{2}$ (91)
$\displaystyle\alpha_{1}^{\prime}(t)$
$\displaystyle=P(t,\vec{\alpha},\Gamma_{1})=-\vec{\alpha}^{T}\Gamma^{1}\vec{\alpha}$
(92) $\displaystyle\alpha_{2}^{\prime}(t)$
$\displaystyle=Q(t,\vec{\alpha},\Gamma_{2})=-\vec{\alpha}^{T}\Gamma^{2}\vec{\alpha}$
(93)

where $\vec{\alpha}^{T}=[\alpha_{1},\alpha_{2}]$ is the tangent vector. The
main difference in comparison to the random field model is that we do not need
to perform MCMC simulations to generate outcomes of the random variable.
Algorithm 3 shows the pseudocode employed to generate the results in the first
set of experiments.

Algorithm 3 Computing geodesic curves in Gaussian manifolds with RK4

1:function GeodesicDistanceGaussian($a$, $b$, $n$, $N$, $\vec{\gamma}^{(0)}$,
$\vec{\alpha}^{(0)}$)

2: $h=(b-a)/n$

3: $t=a$

4: $L=0$ $\triangleright$ L denotes the geodesic distance

5: for $i=0$, $i<N$, $i++$ do

6: Compute the metric tensor $g(\vec{\gamma}^{(i)})$ as matrix in equation
(82).

7: Compute the Christoffel symbols ($\Gamma_{1}$, $\Gamma_{2}$) as matrices in
equation (83).

8: // Perform a 4th order Runge Kutta iteration

9: $k_{0}=hF(t,\vec{\alpha})$

10: $l_{0}=hG(t,\vec{\alpha})$

11: $x_{0}=hP(t,\vec{\alpha})$

12: $y_{0}=hQ(t,\vec{\alpha})$

13: $\vec{\delta}=[k_{0},l_{0}]$

14: $k_{1}=hF(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

15: $l_{1}=hG(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

16: $x_{1}=hP(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{1})$

17: $y_{1}=hQ(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{2})$

18: $\vec{\delta}=[k_{1},l_{1}]$

19: $k_{2}=hF(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

20: $l_{2}=hG(t+0.5h,\vec{\alpha}+0.5\vec{\delta})$

21: $x_{2}=hP(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{1})$

22: $y_{2}=hQ(t+0.5h,\vec{\alpha}+0.5\vec{\delta},\Gamma_{2})$

23: $\vec{\delta}=[k_{2},l_{2}]$

24: $k_{3}=hF(t+h,\vec{\alpha}+\vec{\delta})$

25: $l_{3}=hG(t+h,\vec{\alpha}+\vec{\delta})$

26: $x_{3}=hP(t+h,\vec{\alpha}+\vec{\delta},\Gamma_{1})$

27: $y_{3}=hQ(t+h,\vec{\alpha}+\vec{\delta},\Gamma_{2})$

28: Update the vectors $\vec{\gamma}^{(i)}$ (position) and
$\vec{\alpha}^{(i)}$ (velocity) as

29: $\gamma_{1}^{(i+1)}=\gamma_{1}^{(i)}+(1/6)(k_{0}+2k_{1}+2k_{2}+k_{3})$

30: $\gamma_{2}^{(i+1)}=\gamma_{2}^{(i)}+(1/6)(l_{0}+2l_{1}+2l_{2}+l_{3})$

31: $\alpha_{1}^{(i+1)}=\alpha_{1}^{(i)}+(1/6)(x_{0}+2x_{1}+2x_{2}+x_{3})$

32: $\alpha_{2}^{(i+1)}=\alpha_{2}^{(i)}+(1/6)(y_{0}+2y_{1}+2y_{2}+y_{3})$

33: Update the geodesic distance: $L=L+\lVert\alpha(t)\rVert h$

34: end for

35:end function

Table 1 shows some obtained results using the proposed simulation method. Note
that, as expected, the geodesic distances are greater than the simple improper
2D Euclidean distance between the initial and final parameter vectors in
$R^{2}$. Additionally, the difference between the original geodesic distance
and the geodesic distances obtained by time reversing the Runge-Kutta
simulations are negligible, probably due to small numerical instabilities and
minor rounding errors. It is worth to mention that the Gaussian random
variables manifold has constant negative curvature, which means a hyperbolic-
like geometry. The parameters for all the simulations were set to $a=0$,
$b=10$ and $n=2000$, which leads to step size $h=0.005$.

Table 1: Results of the fourth-order Runge-Kutta simulation for the regular Gaussian density manifold. Init. point (A) |  Initial $\vec{\alpha}$ | Final point (B) | Final $\vec{\alpha}$ | Distances
---|---|---|---|---
$\mu$ | $\sigma^{2}$ | $\alpha_{1}$ | $\alpha_{2}$ | $\mu$ | $\sigma^{2}$ | $\alpha_{1}$ | $\alpha_{2}$ | Geo (AB) | Geo (BA) | Eucl.
0 | 1 | 0.1000 | 0.0500 | 1.0765 | 0.9576 | 0.0957 | -0.0554 | 1.1246 | 1.1250 | 1.0769
2 | 1 | 0.2000 | 0.05 | 3.4723 | 0.2818 | 0.0561 | -0.0689 | 1.7293 | 1.7251 | 1.6373
1 | 2 | 0.2500 | 0.2500 | 3.6214 | 1.1789 | 0.1471 | -0.2396 | 3.3177 | 3.3178 | 2.7461
10 | 5 | 0.5000 | 0.5000 | 13.8957 | 1.2913 | 0.1286 | -0.3741 | 6.4719 | 6.4656 | 5.3787
10 | 10 | 0.5000 | 2.0000 | 17.6233 | 11.3884 | 0.5696 | -2.0745 | 17.0900 | 17.0826 | 7.7444
100 | 100 | 1.0000 | 1.0000 | 109.0051 | 68.3735 | 0.6830 | -5.4847 | 34.8442 | 35.0553 | 32.8869

Figure 4 illustrates some of the geodesic curves obtained during our
simulations. The blue curve is built during the regular execution of the
fourth-order Runge-Kutta method and the red curve is built during the time
reversed version of the simulation. Note that the curves are almost identical,
suggesting that the dynamics modeled by the system of differential equations
is reversible.

Fig. 4: Geodesic curves in the Gaussian density manifold obtained by fourth-
order Runge-Kutta and its time reverse simulations. A denotes the initial
point and B denotes the final point in the regular simulation and vice-versa
in the time reversed version.

In the second set of experiments, we use the proposed iterative MCMC based
fouth-order Runge-Kutta method to build geodesic curves in the Gaussian random
fields manifold using the metric tensor and the Christoffel symbols derived in
the previous sections. In all simulations, due to the computational cost, the
dimensions of each outcome of the random field are $128\times 128$, the number
of iterations in the MCMC simulation is $100$, $a=0$, $b=10$ and $n=1000$,
which leads to a step size $h=0.01$.

Roughly speaking, the intuition behind the numerical simulation of the Euler-
Lagrange equations is that the initial condition
$\vec{\gamma}^{(0)}=[\gamma_{1}^{(0)},\gamma_{2}^{(0)},\gamma_{3}^{(0)}]$
defines the initial point in the manifold and the initial condition
$\vec{\alpha}^{(0)}=[\alpha_{1}^{(0)},\alpha_{2}^{(0)},\alpha_{3}^{(0)}]$
defines the initial tangent vector at the initial point (as the initial
velocity vector). Making a simple analogy, the process can be understood as
follows: suppose we place a tiny marble in the manifold defined by the
parametric space of our random field model. At each iteration of the
algorithm, the system of differential equations updates the position and the
velocity vector of the tiny marble, based in the initial conditions. With
these equations, it is possible to completely characterize the trajectory of
the moving marble along geodesics in the underlying manifold. A limitation,
however, is that, in order to compute the Christoffel symbols, we need to
estimate the covariance matrix of the patches along the random field. In order
to do so, we use Markov Chain Monte Carlo simulation to generate an outcome of
the Gaussian random field at each iteration of the Runge-Kutta method using
the parameter vector $\vec{\gamma}$, which is nothing more than the position
of the marble in the manifold. Hence, as the tiny marble moves along the
manifold, we have different outcomes for our random field model. At the end of
the simulation, the arc length of the resulting geodesic curve defines our
geodesic distance.

Table 2: Results of the proposed iterative MCMC based fourth-order Runge-Kutta simulation for the Gaussian random fields manifold. Init. point (A) |  Initial $\vec{\alpha}$ | Final point (B) | Distances
---|---|---|---
$\mu$ | $\sigma^{2}$ | $\beta$ | $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\mu$ | $\sigma^{2}$ | $\beta$ | Geo (AB) | Geo (BA) | Eucl.
0 | 1 | 0 | 0.1 | 0.1 | 0.2 | 0.7846 | 1.1079 | -0.9172 | 1.7617 | 1.8622 | 1.2126
0 | 1 | 0 | 0.01 | 0.02 | 0.08 | 1.0057 | 0.9054 | -0.4039 | 1.5918 | 1.3889 | 1.0882
0 | 1 | 0 | 0.01 | 0.02 | 0.09 | 1.1394 | 0.8927 | -0.4512 | 1.7850 | 2.2224 | 1.2305
5.0 | 10.0 | 0.5 | -0.1 | -0.1 | -0.2 | 2.0308 | 12.7482 | 1.5256 | 5.6583 | 6.2141 | 4.1742
10.0 | 5.0 | 0.0 | 2.0 | 0.05 | 0.1 | 10.9229 | 5.1097 | -1.3634 | 2.1595 | 1.8070 | 1.6397
5.0 | 10.0 | -0.5 | 0.01 | 0.01 | 0.2 | 9.3335 | 17.7629 | -1.3973 | 11.7472 | 7.3078 | 8.9358
5.0 | 10.0 | -0.5 | 0.1 | 0.1 | 0.5 | 7.6611 | 17.8300 | -2.2929 | 10.1241 | 10.844 | 8.4618
10.0 | 100.0 | 1.5 | -0.1 | -5.0 | -0.5 | 5.3494 | 28.4040 | 1.8986 | 71.759 | 70.8599 | 70.0239
0.0 | 1.0 | 0.0 | 0.5 | 0.5 | 0.5 | 0.7189 | 1.4497 | -1.6750 | 2.4432 | 2.3669 | 1.8788
0.0 | 1.0 | 0.0 | -1.0 | 0.25 | 1.0 | -0.6844 | 0.9558 | -2.0375 | 2.7383 | 2.0959 | 2.1563
0.0 | 1.0 | 0.0 | 1.0 | 1.0 | 1.0 | 0.7200 | 1.4367 | -2.3230 | 3.1187 | 2.4209 | 2.4757
3.1415 | 9.8696 | 0.0 | 0.01 | -0.02 | 0.08 | 3.6003 | 11.5537 | 0.5078 | 2.0014 | 1.8976 | 1.8178

Looking at the obtained results, it is possible to note a curious fact: the
geodesic distances from the initial points A to the final points B are
different from the reversed geodesic distances from final points B to initial
points A, obtained by time reversing the Euler-Lagrange equations. In order to
further investigate this behavior, we build the geodesic curves for the
simulations described in Table 2. Figures 5, 6, 7 and 8 show the geodesic
curves, entropies and curvatures for some numerical simulations.
Interestingly, a rather strange behavior is observed. Under certain
circunstances, the geodesic distances and their time reversed versions seem to
diverge. We call this unexpected way of acting as the geodesic dispersion
phenomenon in random fields dynamics.

Fig. 5: Geodesic curves, entropy and curvature evolution in the Gaussian
random fields manifold obtained by the proposed iterative MCMC based fourth-
order Runge-Kutta algorithm and its time reverse simulations (first row in
Table 2). The blue curve denotes the orignal geodesic and the red curve
denotes its time reversed version.

Fig. 6: Geodesic curves, entropy and curvature evolution in the Gaussian
random fields manifold obtained by the proposed iterative MCMC based fourth-
order Runge-Kutta algorithm and its time reverse simulations (third row in
Table 2). The blue curve denotes the orignal geodesic and the red curve
denotes its time reversed version.

Fig. 7: Geodesic curves, entropy and curvature evolution in the Gaussian
random fields manifold obtained by the proposed iterative MCMC based fourth-
order Runge-Kutta algorithm and its time reverse simulations (sixth row in
Table 2). The blue curve denotes the orignal geodesic and the red curve
denotes its time reversed version.

Fig. 8: Geodesic curves, entropy and curvature evolution in the Gaussian
random fields manifold obtained by the proposed iterative MCMC based fourth-
order Runge-Kutta algorithm and its time reverse simulations (eleventh row in
Table 2). The blue curve denotes the orignal geodesic and the red curve
denotes its time reversed version.

At the moment, we still do not have a clear and definitive explanation for the
emergence of this phenomenon. One possibility for the geodesic dispersion
phenomenon in random fields dynamics can be related to the abrupt changes in
the curvature along the trajectory of the geodesic curve described by the
Euler-Lagrange equations. When the original geodesic curve is built by the
numerical simulations and the system undergoes a phase transition, curvature
goes from minus one to zero, or vice-versa. The construction of the geodesic
curve from A (initial point) to B (final point) starts in a region of the
manifold that has negative curvature. It is known that negative curvature
leads to hyperbolic geometry, in which two “parallel” geodesics tend to
diverge. However, the construction of the time-reversed geodesic curve starts
in a regions of the manifold that has zero curvature, in which two “parallel”
geodesics remain parallel. When we move towards the region of negative
curvature from the region of zero curvature, some degree of divergence is
expected to happen during the process. Moreover, note that in Figures 6, 7 and
8, it is possible to see some fluctuations in the curvatures along the
numerical simulations, leading to some asymmetries in the processes. For the
interested reader, all the Python code used to generate the results presented
in this paper, can be found at
https://github.com/alexandrelevada/GeodesicDispersion.

##  10 Conclusions and final remarks

Random fields dynamics are present in many physical and complex systems due to
their ability to capture the inherent variability and complexity in many
natural phenomena as well as the stochastic nature of several physical
process. Examples of random field applications in physics arise in a variety
of fields, such as statistical physics, condensed matter physics, quantum
field theory, chaos theory, high-energy physics, fluid dynamics, cosmology and
astrophysics.

In this scientific investigation, our main motivating question was: how to
build geodesic curves in the pairwise isotropic Gaussian random fields
manifold? By using MCMC simulation and the fourth-order Runke-Kutta method, we
proposed an iterative algorithm to solve the Euler-Lagrange equations, after
the derivation of the metric tensor and the Christoffel symbols of the metric.
The analysis of the computational simulations revealed a strange fact: the
geodesic dispersion phenomenon in random fields dynamics, which is related to
how geodesic curves diverge from their time-reversed versions when the system
undergoes phase transitions, due to abrupt changes in the local curvature of
the underlying manifold. The obtained results from our numerical methods
presented a series of intriguing patterns and behaviors that significantly
contributed to our understanding of how random fields evolve in time and how
irreversible dynamics can be a direct consequence of geometric properties.
Through a meticulous exploration of the underlying mathematical framework, we
have unveiled the intricate dynamics governing the geodesic dispersion,
shedding light on its nuanced interplay with random fields.

Moreover, our research has established a formal connection between geodesic
dispersion and the intrinsic geometric properties of the underlying space. By
leveraging concepts from differential geometry, we have provided a
mathematical foundation for the observed phenomena. This not only enhances our
theoretical understanding but also opens avenues for developing more accurate
models and predictive frameworks in the context of random fields. The
implications of our findings may lead to relevant theoretical considerations
in the study of several kinds of physical systems, elucidating how the passage
of time in a system can be related to how the intrinsic geometry of its
parametric space changes as we move through the underlying manifold. In random
fields dynamics, the concepts of time and space seen to be intrinsically
related.

In conclusion, our study on the geodesic dispersion phenomenon in random
fields dynamics advances the frontier of knowledge in spatial processes. The
integration of differential geometry, correlation analysis, and practical
applications underscores the multidisciplinary significance of our findings.
We hope that our work will stimulate further research in this domain, offering
new perspectives and methodologies for investigating complex phenomena in
diverse scientific and technological applications.

Future works may include the derivation of the Ricci curvature tensor, which
measures how a shape is deformed as one moves along geodesics in the space.
Intead of using the simple Gaussian curvature, the more appropriate Ricci
curvature is a more precise definition. The Ricci curvature of a manifold is
often characterized as “geodesic dispersion”, namely, whether two parallel
geodesics shot from nearby points converge (spherical geometry or positive
geometry), remain parallel (Euclidean geometry or zero curvature), or diverge
(hyperbolic geometry or negative curvature). However, the components of the
Ricci curvature tensor in Gaussian random fields depend on the derivatives of
the Christoffel symbols, which are not simple to derive analytically.

## Acknowledgments

This study was partially funded by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior — Brasil (CAPES), Finance Code 001.

## References

  * [1] G. Temple, “The physical principles of the quantum theory,” Proc. The Royal Society A: Mathematical, Physical EngineeringSciences 138, 479–494 (1932).
  * [2] A. Nowak, R. Vallacher, A. Tesser, and W. Borkowski, “Society of self: the emergence of collective properties in self-structure.” Psychologicalreview 107.1, 39–6 (2000).
  * [3] G. Wainrib and L. C. G. D. Molino, “Optimal system size for complex dynamics in random neural networks near criticality.” Chaos 23.4, 043134 (2013).
  * [4] E. Hernández-Lemus, “Random fields in physics, biology and data science,” Frontiers Physics 9, 641859 (2021).
  * [5] M. Kardar, _Statistical Physics of Fields_ (Cambridge University Press, 2007).
  * [6] M. Kumar, V. Banerjee, S. Puri, and M. Weigel, “Critical behavior of the three-state random-field potts model in three dimensions,” Phys. Rev.Res. 4, L042041 (2022).
  * [7] A. Turkoglu and A. N. Berker, “Phase transitions of the variety of random-field potts models,” Physica A: Statistical Mechanics itsApplications 583, 126339 (2021).
  * [8] V. Ardourel and S. Bangu, “Finite-size scaling theory: Quantitative and qualitative approaches to critical phenomena,” Stud. History Philosophy Science PartA 100, 99–106 (2023).
  * [9] D. Bachtis, G. Aarts, and B. Lucini, “Quantum field theories, markov random fields and machine learning,” J. Physics: ConferenceSeries 2207, 012056 (2022).
  * [10] Y. Aoun and K. Khettabi, “On the two-point function of the ising model with infinite-range interactions,” J. StatisticalPhysics 190, 168 (2023).
  * [11] A. Malyarenko and M. Ostoja-Starzewski, “Tensor- and spinor-valued random fields with applications to continuum physics and cosmology,” ProbabilitySurveys 20, 1–86 (2023).
  * [12] E. Thélie, D. Aubert, N. Gillet, J. Hiegel, and P. Ocvirk, “Topology of reionisation times: Concepts, measurements, and comparisons to gaussian random field predictions,” Astronomy &Astrophysics 672, A184 (2023).
  * [13] S. Amari, “Information geometry,” International StatisticalReview 89, 250–273 (2021).
  * [14] F. Nielsen, “An elementary introduction to information geometry,” Entropy 22 (2020).
  * [15] F. Nielsen, “The many faces of information geometry,” Notices American MathematicalSociety 69, 36–45 (2022).
  * [16] A. L. Levada, “On the numerical computation of fisher–rao based distances and kl-divergences between gaussian random fields,” J. ComputationalScience 62, 101714 (2022).
  * [17] A. L. Levada, “The curvature effect in gaussian random fields,” J. Physics:Complexity 3, 035001 (2022).
  * [18] U. Lucia, “Irreversible and quantum thermodynamic considerations on the quantum zeno effect,” ScientificReports 13, 10763 (2023).
  * [19] A. Ben-Naim, “Is time’s asymmetry related to irreversible processes and the second law?” Entropy 25, 1297 (2023).
  * [20] R. Zakine and E. Vanden-Eijnden, “Minimum-action method for nonequilibrium phase transitions,” Phys. Rev.X 13, 041044 (2023).
  * [21] S. Weinberg, _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_ (Wiley, New York, NY, 1972).
  * [22] C. W. Misner, K. S. Thorne, , and J. A. Wheeler, _Gravitation_ (Freeman, San Francisco, CA, 1973).
  * [23] G. Arfken, _Mathematical Methods for Physicists_ (Academic Press, 1985), 3rd ed.
  * [24] C. Lanczos, _The Variational Principles of Mechanics_ (Cover, 1986), 4th ed.
  * [25] A. Islam, “Deriving the cosmological constant from the euler–lagrange equation of second-order differentiable gravitational field lagrangian,” AIPAdvances 13, 125130 (2023).
  * [26] V. E. Tarasov, “General fractional classical mechanics: Action principle, euler–lagrange equations and noether theorem,” Physica D: NonlinearPhenomena 457, 133975 (2024).
  * [27] C. Runge, “Ueber die numerische auflusung von differentialgleichungen,” MathematischeAnnalen 46, 167–178 (1895).
  * [28] M. Kutta, “Beitrag zur naherungsweisen integration totaler differentialgleichungen,” Zeitschrift fur Mathematik undPhysik 46, 435–453 (1901).
  * [29] J. C. Butcher, _Numerical Methods for Ordinary Differential Equations_ (Wiley, 2016), 3rd ed.
  * [30] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The J. ChemicalPhysics 21, 1087–1092 (1953).
  * [31] W. K. Hastings, “Monte Carlo sampling methods using Markov chains and their applications,” Biometrika 57, 97–109 (1970).
  * [32] G. L. Jones and Q. Qin, “Markov chain monte carlo in practice,” Annual Review Statistics ItsApplication 9, 557–578 (2021).
  * [33] H. Rue and L. Held, _Gaussian Markov Random Fields: Theory and Applications_ (Chapman & Hall/CRC, 2005).
  * [34] J. M. Hammersley and P. Clifford, “Markov field on finite graphs and lattices (preprint),” (1971). Unpublished.
  * [35] S. D. Hill and J. C. Spall, “Stationarity and convergence of the metropolis-hastings algorithm: Insights into theoretical aspects,” IEEE Control SystemsMagazine 39, 56–67 (2019).
  * [36] C. P. Robert and G. Casella, _Monte Carlo Statistical Methods_ (Springer, 2004), 2nd ed.
  * [37] P. Diaconis, “The markov chain monte carlo revolution,” Bulletin American MathematicalSociety 46, 179–205 (2008).
  * [38] W. Gilks, S. Richardson, and D. Spiegelhalter, _Markov Chain Monte Carlo in Practice_ , Chapman & Hall/CRC Interdisciplinary Statistics (Taylor & Francis, 1995).
  * [39] D. Landau and K. Binder, _A Guide to Monte Carlo Simulations in Statistical Physics_ (Cambridge University Press, 2021), 5th ed.
  * [40] W. Wang, “An introduction to the markov chain monte carlo method,” American J. Physics 90, 921–934 (2022).
  * [41] G. O. Roberts and J. S. Rosenthal, “Complexity bounds for markov chain monte carlo algorithms via diffusion limits,” J. Appl.Probability 53, 410–420 (2016).
  * [42] S. Amari and H. Nagaoka, _Methods of Information Geometry_ (American Mathematical Society, 2000).
  * [43] F. Nielsen, “An elementary introduction to information geometry,” Entropy 22, 1100 (2020).
  * [44] J. M. Lee, _Riemannian manifolds: an introduction to curvature_ (Springer, 1997).
  * [45] T. Needham, _Visual Differential Geometry and Forms: A Mathematical Drama in Five Acts_ (Princeton University Press, 2021).
  * [46] M. P. do Carmo, _Differential Geometry of Curves and Surfaces_ (Dover Publications Inc., 2017), 2nd ed.
  * [47] P. L. Devries and E. Hasbun, J, _A First Course in Computational Physics_ (Jones and Bartlett Publishers, 2011), 2nd ed.
  * [48] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ (Cambridge University Press, 2007), 3rd ed.
  * [49] J. C. Butcher, “On runge-kutta processes of high order,” J. Australian MathematicalSociety 4, 179–194 (1964).
  * [50] J. Butcher, “A history of runge-kutta methods,” Appl. NumericalMathematics 20, 247–260 (1996).
  * [51] E. A, “A study of runge-kutta method for n-th order differential equations,” in _2023 2nd International Conference on Advancements in Electrical, Electronics, Communication, Computing and Automation (ICAECA),_ (2023), pp. 1–4.
  * [52] S. I. Costa, S. A. Santos, and J. E. Strapasson, “Fisher information distance: A geometrical reading,” Discrete Appl.Mathematics 197, 59–69 (2015). Distance Geometry and Applications.
  * [53] W. W. S. Chen, “On computing gaussian curvature of some well known distribution,” Theoretical Mathematics &Applications 3, 85–104 (2013).