We can assign each point in the graph, one of a finite set S of labels. Such assignment, it is often called a configuration. We can then assign probability measures to the set Ω of all possible configurations ω. Hence, ω A represents the configuration ω restricted to the subset A of V. We may think of ω A as a configuration on the subgraph G A restricting V to points of A.

We can define local characteristics on MRFs. The local characteristics of a probability measure ℙ defined on Ω are the conditional probabilities: ℙ ( ω t   |   ω T ∖ t ) = ℙ ( ω t   |   ω N t ) (1)

This represents the probability that the point t is assigned the value ω t , given the values at all other points of the graph. Let us re-write Eq. 1. Since the probability measure will define an MRF if the local characteristics depend only on the outcomes at neighboring points, i.e., if for every ω ℙ ( ω X i   |   ω G ∖ X i ) = ℙ ( ω X i   |   ω N X i ) (2)

Given an arbitrary graph, we may refer to a set of points C, as a clique, if every pair of points in C are neighbors. This includes the empty set as a clique. A clique is then a set whose induced subgraph is complete. Cliques are also called complete induced subgraphs or maximal subgraphs.

A potential η is an assignment of a number η A ( ω ) to every subconfiguration ω A of a configuration ω in the graph G. A given η, induces an energy U ( ω ) on the set of all configurations ω as follows: U ( ω ) = ∑ A η A ( ω ) (3)

Here, for fixed ω, the sum is taken over all subsets A ⊆ V including the empty set. It is possible to define a probability measure, called the Gibbs measure induced by U as ℙ ( ω ) = e − U ( ω ) Z (4)

Z (taken from the German word zustanssumme or sum over states) is a normalization constant called the partition function. As it is known, explicit computation of the partition function is in many cases a very challenging endeavor. There is a great deal of work in the development of methods and approaches to overcome some (but not all) challenges in this regard. Some of these approximations will be discussed later on. Z = ∑ ω e − U ( ω ) (5)

The term potential is often used in connection with potential energies. In this context η A is commonly termed a potential energy in physics applications. ϕ A = e − η A is then called a potential.

Equations 4, 5 can be thus rewritten as: ℙ ( ω ) = ∏ A ϕ A ( ω ) Z (6) Z = ∑ ω ∏ A ϕ A ( ω ) (7)

Since this latter use is more common in probability and graph theory, and it is also used in theoretical physics, we will refer to Eqs. 6, 7 as the definitions of Gibbs measure and partition function (respectively) unless otherwise stated. This will also be justified given that Eq. 6 is a form of probability factorization (in this case a clique factorization) [11].

A potential is termed a nearest neighbor Gibbs potential if ϕ A ( ω ) = 1 whenever A is not a clique. We often call a Gibbs measure to any regular measure induced by a nearest neighbor Gibbs potential. However, we may define more general Gibbs measures by considering different classes of potentials.

The inclusion of all cliques in the calculation of the Gibbs measure is needed to establish the equivalence between Gibbs random fields and Markov random fields. A nearest neighbor Gibbs measure on a graph determines an MRF as follows [22]:

Let ℙ ( ω ) be a probability measure determined on Ω by a nearest neighbor Gibbs potential ϕ: ℙ ( ω ) = ∏ C ϕ C ( ω ) Z (8)

With the product taken over all cliques C on the graph G. Then, ℙ ( ω X i | ω G ∖ X i ) = ℙ ( ω ) ∑ ω ′ ℙ ( ω ′ ) (9)

Here ω ′ is any configuration which agrees with ω at all points except X i . ℙ ( ω X i | ω G ∖ X i ) = ∏ C ϕ C ( ω ) ∑ ω ′ ∏ C ϕ C ( ω ′ ) (10)

For any clique C that does not contain X i , ϕ C ( ω ) = ϕ C ( ω ′ ) , So that all the terms that correspond to the cliques that do not contain the point X i cancel both from the numerator and the denominator in Eq. 10, therefore this probability depends only on the values x i at X i and its neighbors. ℙ defines thus an MRF. A more general proof of this equivalence was given by Hammersley-Clifford theorem (see for instance [11]).

In essence, we can state that among the general class of random fields, Markov random fields are defined by obeying the Markov neighborhood law. Gibbs fields are usually understood as Markov fields with strictly positive probability measures (in particular, a strictly positive joint probability density). These Markov-Gibbs fields are thus defined by the Markov property and the positive definite probabilities and are the ones that follow the Hammersley-Clifford theorem. More general Gibbs fields can be defined by other neighborhood laws than the Markov property [23], but these will not be addressed in the present work.

To discuss the conditional independence structure induced by MRFs, let us consider the following: An adjacency matrix A i j represents the neighborhood law (as given by the Markov property) on the graph G. Every non-zero entry in this matrix represents a statistical dependency relation between two elements on X. The conditional dependence structure on MRFs is related not only to the local statistical independence conditions, but also to the dependency structure of the whole graph [11, 24].

A definition of conditional independence (CI) for the set of random variables can be given as follows: ( X i ⊥ ⊥ X j ) | X l ⇔ F X i , X j | X l = X l * ( X i * , X j * ) = F X i | X l = X l * ( X i * ) ⋅ F X j | X l = X l * ( X j * ) (11) ∀   X i , X j , X l ∈ X

Here ⊥ ⊥ refers to conditional independence between two random variables. F X i , X j | X l = X l * ( X i * , X j * ) = P r ( X i ≤ X i * , X j ≤ X j * | X l = X l * ) is the joint conditional cumulative distribution of X i and X j given X l . X i * , X j * and X l * are realizations of the corresponding random variables.

In the case of MRFs, CI is defined by means of graph separation: Hence X i ⊥ ⊥ G X j | X l iff X l separates X i from X j in G. This means that if we remove node X l there are no undirected paths from X i to X j in G.

Conditional independence in random fields can be considered in terms of subsets of V. Let A, B and C be subsets of V. The statement X A ⊥ ⊥ G ^ X B | X C , which holds only iff C separates A from B in G, means that if we remove all vertices in C there will be no paths connecting any vertex in A to any vertex in B. This is customarily called the global Markov property of TMFs [11, 24].

The smallest set of vertices that renders a vertex X i conditionally independent of all other vertices in the graph is called its Markov blanket, denoted m b ( X i ) . If we define the closure of a node X i as C ( X i ) then X i ⊥ ⊥ G ∖C ( X i ) | m b ( X i ) .

In an MRF, the Markov blanket of a vertex is its set of first neighbors. This statement is the so-called undirected local Markov property. Starting from the local Markov property, it is possible to show that two vertices X i and X j are conditionally independent given the rest if there is no direct edge between them. This is the pairwise Markov property.

If we denote by G X i → X j the set of undirected paths in the graph G connecting vertices X i and X j , then the pairwise Markov property of an MRF is given by: X i ⊥ ⊥ X j | G ∖ { X i , X j }   ⇔   G X i → X j = ∅ (12)

Hence the global Markov property implies the local Markov property which, in turn, implies the pairwise Markov property. For systems with positive definite probability densities, it has been proved that pairwise Markov actually implied global Markov (See [11] p. 119 for a proof). This is important for applications since it is easier to assess pairwise conditional independence statements.

Due to their intrinsic simplicity and generality, MRFs have attracted the attention of mathematical physicists and probability theorists looking to extend their associated theoretical foundations. Important work has been done, for instance, to incorporate geometrical properties and generalized embeddings to the theory of random fields. Extremely relevant in this regard is the monumental work presented in the monograph by Adler and Taylor [31]. There, the authors expand on the consideration of a random field as a stochastic process in a metric space (discrete, Euclidean, etc.) to consider random fields as stochastic mappings over manifolds. This extension is given via writing down differential geometry characterizations of the fields based on a measure-theoretic definition of probability. Though this work may seem quite abstract, it was indeed born out of an idea for an application of random fields to neuroscience. Nurturing from similar ideas, recent work by Ganchev [32] has expanded the notion of locality of MRFs and assimilate it to the geometric features present in lattice quantum gauge theories, to generate a gauge theory of Markov-Gibbs fields. Again, even if the setting seems to be quite theoretical, an application to the modeling of trading networks in finance is given.

Other mathematical extensions of Markov random fields are related to the nature of the graphical model considered. In general, probabilistic graphical models may belong to one of two quite general classes: Markov networks (such as MRFs) which are undirected graphs or Bayesian networks which are directed graphs. The difference between undirected and directed graphical models impose consequences in the kind of fundamental mathematical objects of the theory: joint probabilities or conditional probabilities, loopy graphs or trees–directed acyclic graphs–, clique factorization vs. conditional probability factorization via the chain rule, etc. Whether the model is undirected or directed also has modeling and computational consequences. To be fair, both models have pros and cons.

Trying to overcome the limitations of both general approaches, Freno and Trentin [33] developed a more general approach to random fields termed Hybrid random fields (HRFs). The purpose of HRFs is to allow the systems to present a wider variety of conditional independence structures. As we will discuss later, allowing for a systematic incorporation of more general classes of conditional independence structures in indeed one of the current hot topics in computational intelligence and machine learning. Actually, even when HRFs are theoretical constructs (much alike MRFs) they were designed to be learning machines, i.e., to be supplemented with training algorithms to deal with high dimensional data. HRFs were developed for logical inference in the presence of partial information or noise. As in the case of MRFs and of their gauge extensions just mentioned, HRFs were developed to rely on a principle of locality which is an extension of the Markov property that allows for sparse stochastic matrix representations amenable for the computation on actual applications. Once a (graph) structure has been given (or inferred) HRFs are able (as is the case of MRFs) to learn the local (conditional or joint-partial) probability distributions from empirical data, a task commonly known in statistics as parameter learning [34]. Hence HRFs are theoretically founded, but developed thinking in applications. The scope of applicability of MRFs has also become broader by expanding its applicability to model tensor valued quantities [35], giving rise to the so-called multilayer graphical models, also called multilayer networks [36–39].

Aside from expanding the fundamental structure of MRFs, mathematical physics applications of Gibbs random fields are abundant. In particular, the so-called Random Field Ising model (RFIM) has gained a lot of attention in the recent years. By using the monotonicity properties of the associated stochastic field, Aizenmann and Peled [40] were able to prove that there is a power law upper bound on the correlations on a two-dimensional Ising model, supplemented with a quenched random magnetic field. The fact that by combining random fields (the intrinsic Ising field and the quenched magnetic field), the nature of the phase transitions may drastically change has made the RFIM a current topic of discussion in mathematical statistical mechanics. The consequences of the induction of long range order in the RFIM, leading to the emergence of the so-called Imry-Ma phase or Imry-Ma states (named so since Imry and Ma were actually behind the first proposal of the RFIM [41]) have been the object of intense study recently. Berzin and co-workers [42] used MRFs to analyze the dynamic fluctuations of the order parameter in the Imry-Ma RFIM and its coupling with the static fluctuations of the structural random field (accounting for the defects). Interestingly, anisotropic coupling arises from two non-absolutely overlapping local fields [43]. The effects of the non-overlapping fields in anisotropy and disorder has been studied since several decades ago [44], but the actual relationship with non-locality was established relatively recently. For instance, it was until 2018 that Chatterjee was able to quantitatively describe the decay of correlations of the 2D RFIM [45] in a relevant paper that led Aizenmann to re-analyze his former, mostly qualitative proposal [40, 46].

Local stochastic phenomena in non-homogeneous and disordered media in the context of the RFIM has also attracted attention in relation to critical exponents and scaling. Trying to expand on the origins of long range order from local interactions, Fytas and coworkers have studied the 4D RFIM and its hyperscaling coefficients [47]. This is particularly interesting since it has been shown, via perturbative renormalization group calculations, that the critical exponents of the RFIM in D dimensions are the same as the exponents of the pure Ising model in D − 2 dimensions [48]. Related work has been carried out by Tarjus and Tissier, but they instead resort to the use of the so-called functional renormalization group approach in the multi-copy formalism setting [49]. Their work has extended the predictive capabilities of MRFs by incorporating ideas from symmetry breaking allowing to characterize not just long-range order (LRO) but also intermediate states characterized by quasi-long range order (QLRO). The fact that QLRO may be attained from purely Markov statistics (localized interactions) is in itself appealing for statistical physics. The fact that local dependencies may suffice to account for LRO and QLRO under certain conditions that do not violate the Markov property of the MRFs, will have relevant consequences for the applications of MRFs outside physics, such as in the case of image reconstruction and pattern recognition in machine learning. We will come back to these ideas later on.

Locality as depicted in MRFs can also have important consequences for the theory of fluctuations in fields of interacting particles. Reconstructing Boltzmann statistics from local Gibbs fields (that as we have repeatedly stated are formally equivalent to MRFs, provided strictly positive probability measures) imply that under central limit scales the fluctuation field of local functions can be represented instead as a function of the density fluctuation field, in what is known as the Boltzmann-Gibbs principle (BGP). It has been shown that the BGP induces a duality whose origins are purely probabilistic, i.e., is independent of the nature of the interactions provided their compliance with the tenets of MRFs [50].

It is worth noticing that these contemporary developments in the formal theory of MRFs are actually founded on seminal work by probability theorists and mathematical physicists such as Dobrushin, Ruelle, Gudder, Kessler and others. For instance, Dobrushin laid out the essential conditions of regularity that allow to make explicit the conditional probabilities in MRF models [8]. This work, further developed by Lanford and Ruelle [9] gives rise to the so called Dobrushin-Lanford-Ruelle (DLR) equations that established, in a formal way, the properties of general Gibbs measures. Later on, Dobrushin expanded on these ideas by applying perturbation methods to generalize Gibbs measures to even wider classes of interactions (i.e., to include other families of potentials) [51]. An application of these ideas in quantum field theory can be found in [52] within the context of (truncated) generalized Gibbs ensembles.

Aside from measure-theoretical and algebraic foundations of MRFs, important developments were made by considering explicit dependency structures. In particular, the introduction of strong independence properties led to the formal definition of Gaussian random fields by Gudder [53]. Much of this earlier work has been summarized in the monograph by Kindermann and Laurie Snell [22]. The fact that MRFs are characterized by Gibbs measures even for many-body interactions (under special conditions), and not only for paired-potentials, was already envisioned by Sherman [54], though it remained an unfinished task for decades. Many body effects have actually been reported in the context of localization in the random field Heisenberg chain [55]. One step ahead toward generalizing MRFs consisted in exploring the equivalence of some properties of random fields in terms of sample functions. In this regard, Starodubov [56] proved that there are random fields stochastically equivalent to an MRF, but defined on another probability triple whose sample functions belong to a map associated with the original MRF. The existence of such mappings has relevant implications for applications, in particular in cases in which explicit computation of the partition function is intractable.

Discrete and continuous versions of random fields have been applied to model systems in condensed matter physics and materials science (CMP/MS). The relevance of MRFs and its extensions relies on their suitability to describe the onset of spatio-temporal phenomena from localized interactions. Acar and Sundararaghavan [57] have used MRFs to model the spatio-temporal evolution of microstructures, such as grain growth in polychrystalline microstructures as captured by videomicroscopy experiments. Experimental data is the foundation for explicit calculations of the (empirical) conditional probability distributions.

Gaussian random fields have been used to model quenched random potentials in fluids via mode-coupling by Konincks and Krakoviack [58], and to model beta-distributed material properties by Liu and coworkers [59]. These and other extensions in CMP/MS made use of continuous, piecewise continuous or lattice fluid extensions of Gibbs random fields. Such is also the case of the work of Chen and coworkers [60] who introduced stochastic harmonic potentials in random fields to account for the effects of local interactions on the properties of structured materials; of the work by Singh and Adhilari [61] on Brownian motion in confined active colloids and of the work of Yamazaki [62] on stochastic Hall magnetohydrodynamics. A semi-continuous approach (called smoothed particle hydrodynamics, SPH), using discrete MRFs and extension theorems, was used by Ullah and collaborators [63] in their density dependent hydrodynamic model for crowd coherency detection in active matter.

Extending the ideas of the classic RFIM, Tadic and collaborators [64] were able to describe critical Barkhausen avalanches in quasi-2D ferromagnets with an open boundary. The use of MRFs with disordered field components has also allowed to characterize embedded inhomogeneities in the spectral properties of Rayleigh waves with application to the study of the Earth’s microseismic field [65]. Geoacustic measurements and its MRF modeling allowed these researchers to estimate the mechanical and structural properties of the Earth’s crust and upper mantle. Accurate estimates of these properties are foundational to develop seismic-resistant devices and structures.

MRFs have also been applied in other areas of physics aside from statistical mechanics and condensed matter. MRFs were applied for instance, in geophysical models of marine climate patterns [66], to study reservoir lithology [67] and subsurface soil patterns [68] from remote sensing data. Aside from geophysics, optics and acoustics have also incorporated MRF applications. In acoustics, for instance, an MRF formalism can be used for the isolation of selected signals [69]; or for the segmentation of sonar pulses [70]. In chemical physics, MRFs are applied for the analysis of molecular structures [71], and in the implementation of quantum information algorithms for molecular physics modeling [72].

Disparate as the applications of MRF in the physical sciences just presented may be, these are neither a comprehensive nor even a representative list. However, we expect that some of the essential aspects of its wide range of applicability and the large room for theoretical development still available for these types of models were captured in the previous discussion. Moving on to applications and developments in other disciplines, such as Biology/Biomedicine and the Data Sciences, we will try to convey, not just the usefulness of a quintessential model in statistical physics in other realms–which is huge, indeed–. We also intend to show how some of the implementations and theoretical improvements in other disciplines, can be exported back to physics and may help to solve some of the many remaining conundrums of the theory and applications of random fields in the physical sciences.

One somehow natural application of MRFs is imaging de-noizing or segmentation. This is a quite general problem in which one wishes to discern patterns from a blurred image. In particular an MRF is built to discern which points in imaging space (pixels, voxels) are locally correlated with each other, pointing out to their membership to the same object in the image. The Markov neighborhood structure of the MRF is hence used to un-blur patterns and being able to accurately interpret the images. Often MRFs (or its associated conditional Random fields) are used in conjunction with inference machines such as Convolutional Neural Networks (CNNs). This is the case of the work by Li and Ping [73] who used a neural conditional random field (NCRF) for metastasis detection from lymph node slide images. Their NCRF approach infers the spatial correlations among neighboring patches via a fully connected conditional MRF incorporated on top of a CNN feature extractor. Their modeling approach used a conditional distribution of an MRF with a Gibbs distribution. As is often the case the energy function (i.e., the Hamiltonian) consists of two terms, one summarizing the contributions from unary potentials characteristic for each patch, and the other one summing the pairwise potentials measuring the cost of jointly assigning two neighboring patches (i.e., the interaction potentials).

As is common in physics, estimating the marginals is an intractable problem. Li and Ping resorted to using a mean-field approach and then conditioning their results on this mean field calculations. In order to do this, they trained a CNN with the empirical data. CNN-MRF approaches have also been recently applied to successfully discern computerized tomography imaging (CT scans) [74] for prostate and other pelvic organs at risk. After processing the data with an encoder/decoder scheme, the output of CNN was used as the unary potential of the MRF. Then via a MRF block model based on local convolution layers, a global convolution layer, and a 3D max-pooling layer the authors were able to calculate the pairwise potential. The maximum likelihood optimization problem was then solved via an adaptive loss function.

A similar approach was followed by Fu and collaborators [75] to solve the retinal vessel segmentation problem, fundamental in the diagnostics and surgery of ophthalmological diseases, and, until quite recently manually performed by an ocular pathologist. The authors also used a two term energy function within a mean field approach. To minimize the energy function subject to empirical constraints they used a recurrent neural network based on Gaussian kernels on the feature vectors applying standard gradient descent methods. Blood vessel segmentation was also studied using conditional MRFs by Orlando and coworkers [76]. However, instead of using a mean-field approach and inferring the marginals using neural networks, these authors chose to perform Maximum a Posteriori (MAP) labeling with likelihood functions optimized via Support Vector Machines (SVMs). Imaging segmentation via MRFs can be applied not only at the tisular level, but also on cellular (and even supramolecular) scales. Several blood diseases, for instance, are diagnosed by discerning the quantity, morphology and other aspects of leukocytes as well as their nuclear and cytoplasmic structure. To this end, Reta and coworkers used unsupervised binary MRFs (i.e., classical Ising-like fields) to study leukocyte segmentation [77]. A Markov neighborhood and clique potential approach was followed. This classic approach has been enough since from their high quality colored imaging data, it was possible to define an energy function based on a priori Gaussian-distributed probabilities, then applying a maximum likelihood approach to calculate the posterior probability. Related ideas were used to study microvasculature disorders in glioblastomas by the group of Kurz [78].

Computational biology and bioinformatics are also disciplines that have widely adopted the random field formalism as a relevant component of their toolkits. There are several instances in which MRFs can be adapted to solve problems in these domains: from structural biology problems in which the spatio-temporal locality is naturally mapped onto random fields, to molecular regulatory networks in which the graph structure of the MRFs mimic the underlying connectivity of the networks, to semantic and linguistic segmentation problems in genomic sequences or biomedical texts.

Regarding computational models in structural biology, Rosenberg-Johansen and his group [82] used a combination of deep neural networks and conditional random fields to improve predictions on the secondary structure of proteins (i.e., the three dimensional conformation of local protein segments, the formation of alpha helices, beta sheets and so on). The CRF approach was quite useful in this case (in general non-computationally tractable), since in protein secondary structure, there is a high degree of crosstalk between neighboring elements (residues), then the local dependency structure greatly shrinks the search space. Previously, Yanover and Fromer [83] applied an MRF formalism for the prediction of low energy, protein side configurations, a relevant problem fro several aspects of structural biology such as de novo protein folding, homology modeling and protein-protein docking. The different types of local interactions among amino acid residues: hydrophobic, hydrophilic, charged, polar, etc.) modeled as pairwise potentials let to semi-empirical expressions for the potential energies used in the MRF formalism. Once explicit expressions for the field have been written, the authors resort to a belief-propagation algorithm to find the optimal solution to the MRF problem given the constraints. Several improvements were actually applied to the message-passing algorithm that allow the authors to find a method to obtain the lowest energy amino acid chain configurations. This kind of approach may also be relevant to improve solving methods of random fields in statistical physics problems since it led to approximate explicit forms of the partition function.

Improving methods to discern the structural properties of proteins are also quite used in the context of protein homology, i.e., to investigate on the functions of proteins related to their structural similarity to other proteins, perhaps in different organisms. Local homology relationships can also be investigated by means of Markov random field methods. Xu and collaborators developed a method (or better, a family of methods) called MRFalign for protein homology detection based on the alignment of MRFs [84, 85]. Aside from purely Ising approaches, other methods of random fields of statistical mechanics have been adopted in the computational biology community. One of them is the Potts model. Recently, Wilburn and Eddy used a Potts model with latent variables for the prediction of remote protein homology (involving changes such as insertions and deletions) [86] importance sampling from extensive databases was used to perform MAP optimization as commonly done in computational biology and computer science.

A topic related to homology, but also involving space-dependent electrostatic interactions (protein-protein interactions, in particular) is protein function prediction. Networked models of protein prediction have been developed: primitive models can be used to associate a function to a given protein given the functions of proteins in their interaction neighborhood and probabilistic models may do this by weighting interactions with an associated probability. Gehrman and collaborator devised a CRF method fro protein function prediction based on these premises [87]. To solve the CRF, they resort to a factor graph approach [88] to write down explicit contributions to the cliques [89] and then using an approximate Gibbs measure calculated from this clique factorization. The approximation is based on other relevant feature of Markov random fields, which we will discuss later in the context of statistics and computer science: the use of the so-called Gibbs sampler or Gibbs sampling algorithm [90]. The Gibbs sampler is a Markov chain Monte Carlo (MCMC) method used to obtain a sequence of observations–approximated from a specified multivariate probability distribution–, in those cases for which direct sampling is difficult or even impossible (e.g., NP-hard or super-combinatorial problems).

Perhaps not so well known as a relevant structural biology problem until recently, is the determination of three dimensional chromosome structure inside the cell’s nucleus. Long range chromosomal interactions are believed to be ultimately related to fundamental issues on global and local gene regulation phenomena. A recently devised experimental method for global chromosome conformation capture is known as Hi-C. Nuclear DNA is subject to formaldehyde treatment to enhance covalent interactions glueing chromosome segments that are three dimensionally adjacent. Then a battery of restriction enzymes is used to cut DNA into pieces. Such pieces are sequenced and the identity of the spatially adjacent regions are then discovered. The data is noisy and often incomplete. For these reasons, a team lead by Yun Li developed a hidden Markov random field method to analyze Hi-C data to detect long range chromosomal interactions [91]. This method combines ideas from MRFs, Bayesian networks and Hidden Markov models. In a nutshell, they assumed a mixture of negative binomials as an Ising prior [22] and supplemented it with Bayesian inference to calculate the joint probabilities via a Metropolis-Hastings pseudo-likelihood approach.

Other applications of random fields in biology include demography and selection to study weakly deleterious genetic variants in complex demographic environments [111] and for species clustering [112], in population genetics. MRFs have also been applied to understand species distribution patterns and endemism and to unveil [113] interactions between co-occurring species in processes governing community assembly [114]; as well as for spatially explicit community occupancy [115] in ecology.

Another group of disciplines in which MRFs have flourished is comprised of Data Science, Computer Science and Modern statistics. The next section will be devoted to presenting and discussing some developments of random fields in that setting.

As we already mentioned in the context of applications of random field to biomedical imaging, segmentation and pattern identification to enhance the resolution of spatial and/or spatio-temporal maps is a common use of MRFs. From the many applications in the field of computerized image processing, we will discuss some that present peculiarities or distinctive features that may be of more general interest. For instance, to face the challenge of capturing three dimensional structure from two-dimensional images, the so-called depth perception, Kozik used an MRF-based methodology [116] in which the energy function was modeled via a polynomial regression model and a depth estimation algorithm with correlated uncertainties (a sort of twofold autoregressive model). By using these entries Kozik then solved an MAP problem to obtain the maximum likelihood solution to the MRF.

In the context of AI to enhance low-resolution images (the super-resolution problem), Stephenson and Chen devised an adaptive MRF method [117] based on passing-message optimization by a loopy propagation algorithm. Also in the context of AI approaches to image processing Li and Wand developed a combination of MRFs as generative models and deep CNNs to discriminate two-dimensional images to try to solve the so-called image synthesis problem, a relevant problem in computer vision with applications both to photo-editing and neuroscience [118]. A problem related to image synthesis is image classification, in which certain features of images are discerned and used to cluster images by similitudes in these feature spaces. Applications in image recognition in security, forensics and scientific microscopy and imaging among others abound. To improve the accuracy of image classification algorithms, Wen and coworkers developed a CRF method in which machine-learned feature functions took the place of the unary and binary terms in the potential energy [119], as in previous cases Gaussian priors and loopy belief propagation algorithms were used to solve the random field.

Geostatistics and geographical information systems are also quite amenable to be modeled within the MRF paradigm due to their natural spatio-temporal dependency structures. In the context of prediction of environmental risks and the effects of limited sampling, Bohorquez and colleagues developed an approach based on multivariate functional random fields for the spatial prediction of functional features at unsampled locations by resorting to covariates [120]. As in the case of random field hydrodynamics (mentioned in the physics section), an empirical approach based on continuous field estimators was chosen. Continuous spatio-temporal correlation structures via so-called Kriging methods extending the ideas of discrete random fields are commonly used in environmental analysis and risk assessment [121, 122].

Geological modeling is another field at the intersection of geostatistics and geophysics which has adopted the MRF formalism to deal with their problems. A segmentation approach was used for stochastic geological modeling with the use of hidden MRFs [123]. Using a methodological approximation similar to the one used in computer vision and biomedical imaging, latent variable MRFs are used to perform three-dimensional segmentation. The model is supplemented with finite Gaussian mixture models for the parameter calculations and a Gibbs sampling inference framework, following a similar approach to the one developed by the group of Li [124], based on the methods of Rue and Held [125] and by Solberg et al [126] and further developed by Toftaker and Tjelmeland [127]. More refined geostatistical methods have been based on a clever combination of several developments of Markov random field theory. Along these lines, the work by Reuschen, Xu and Nowak [128] is noteworthy, since they used Bayesian inversion (based on Markov conditional independence) to develop a random field approach to hierarchical geostatistical models and used Gibbs sampling MCMC to solve them.

The combined use of ideas from Markov and Gibbs random fields in statistical learning and other approaches in modern statistics has indeed become a fruitful line of research with important theoretical developments and a multitude of applications [24, 34, 129]. The use of MRFs and CRFs as tools for statistical learning has been used in a multitude of settings in both generative and discriminative models [33]. Aside Ising models and MRFs, perhaps the most widely used applications of the random fields are Gibbs sampling and Markov chain Monte Carlo methods that we already mentioned. Due to the generality and the relatively low computational complexity of these sampling/simulation methods, several methods have been developed based on them.

Gibbs sampling is a form of Markov chain Monte Carlo (MCMC) algorithm. MCMC methods are used to obtain a sequence of observations of a random experiment by an approximation from a given (specified) multivariate probability distribution when direct sampling is challenging (computationally or otherwise). The essence of the method is building a Markov chain whose equilibrium distribution is precisely the specified multivariate distribution. Then, a sample of such distribution is just a sequence of states of the Markov chain. The use of the Markov property of an MRF allows to use Gibbs sampling as an MCMC method, when the joint probability distribution is not known (or is very complex) but the conditional distributions are known (or easier). Due to this, by using the pairwise Markov property, Gibbs sampling is particularly fit to sampling the posterior distribution of Bayesian networks (understood as a collection of conditional distributions), a quite relevant problem in both, statistical learning and in large computer simulation problems.

Aside from these basic issues, Gibbs sampling has been extensively enhanced over the years. One important improvement has been the incorporation of adaptive rejection sampling [130, 131], particularly useful for situations in which evaluation of the density distribution function is computationally expensive (e.g., non-conjugated Bayesian models). Adaptive rejection sampling can be even applied to modeling via non-linear mixed models [131]. To further minimize the computational burden of Gibbs sampling, Meyer and collaborators [132] developed an algorithm which samples via Lagrange interpolation polynomials, instead of exponential distributions. Convergence can be also improved by double-adaptive independent rejection sampling [133] which is based on a scheme of minimizing the correlation among samples. Gibbs sampling approaches also allow for the determination of dense distribution simulated sampling from sparse sampled data [134], even in high dimensional latent fields over large datasets [135].

Gao and Gormley implemented a Gibbs sampling scheme based on CRFs weighted via neural scoring factors (implemented as parameters in factor graphs) with applications to Natural Language Processing (NLP) [136]. MCMC has also been used, in the context of Gibbs random fields in data pre-processing, to reduce the computational burden of data intensive signal processing [137, 138]. Gibbs sampling can also be applied in parallel within the context of Gaussian MRFs on large grids or lattice models [139]. Parallel Gibbs sampling methods can also be developed in the context of sampling acceleration for structured graphs [140].

Markov random fields and its associated Gibbs measures can also be used to advance statistical methods in large deviation theory [141] and to develop methods of joint probability decomposition based on product measures [142]. Exact factorizability of joint probability distributions is a most relevant question in modern probability [143–146] with important applications in data analytics [147], applied mathematics [148], computational biology [149] and network science [150], among other fields. MRFs also have been applied to embed filtrations on high dimensional hyperparameter spaces. The main idea is using random fields as hierarchical models projecting the relevant hyper-parameter space to a lower dimensional filtration [135]. This general problem is closely related with the feature selection problem in computer science and data analytics. We will discuss applications of the MRF formalism in that context in the next subsection.

Feature selection (FS) refers to a quite general class of problems in computer science, data analysis and AI. Feature selection aims to find the minimum number of maximal relevant features to characterize a high dimensional data set. One outstanding family of methods of feature selection is regression methods in which a set of regression variables is used to predict one (or a few) dependent variables via functional relationships (commonly linear combinations with a distribution of weights). A subset of the whole set of regression variables is considered statistically significant, in that context those are the selected features. FS is a more general problem than linear, multivariate or even non-linear regression. MRF can be used to generalize regression procedures to more complex situations. One notable method was developed by Stoehr, Marin and Pudio [151] who used hidden Gibbs random fields to implement model selection via an information theoretical optimization criterion known as Block likelihood information. Cilla and coworkers [152] developed a FS method to be used in sequence classification based on hidden CRFs supplemented with a generalized Lasso group regularization method that instead of the colinearity condition employs L1-norm optimization of the parameters. The authors showed that FS outcomes with this method outperforms standard conditional random field approaches.

Feature selection efficacy of MRFs is closely related to the actual structure of the underlying adjacency matrices. Especially relevant is the issue of separability. Although non-trivial separability does not preclude the use of MRFs in large datasets, as long as the positive definite nature of the measures is ensured; there may be computational complexity limitations for practical uses. Recently, Sain and Furrer [153] discussed on some general properties of random fields (in particular for multivariate Gaussian MRFs) that need to be taken into account in the design of computationally efficient modeling strategies with such random fields. By designing FS schemes with MRFs based on the optimization of parameter estimation, for instance via structured learning it is possible to improve substantially on the computational complexity of such algorithms [154–158]. The graph structure of MRFs can also be optimized to enhance the FS capabilities of the algorithms [159–163]. More information along these lines can be found in the comprehensive review by Adams and Beling [164] and in the one by Vergara and Estevez [165].

As already mentioned the structure of MRF may result advantageous to solve segmentation problems or delimitation of statistical dependencies. These are problems that are extremely relevant in the context of computational linguistics and natural language processing applications. We will discuss these in the following subsection.

Automated textual identification and meaning discernment are extremely complex (and very useful) tasks in current artificial intelligence research and applications. The ability to detect text patches with semantic similarity is one of the founding steps in the ability to process natural language by a computer. By combing a deep learning approach (a convolutional neural network) with MRF models, Liu and collaborators [166] devised an effective algorithm for semantic segmentation [167], which they called a Deep Parsing Network (DPN). Within the DPN scheme, a CNN is used to calculate the unary terms of a two-term energy function, while the pairwise terms were approximated with a mean-field model. The mean field contributions were iteratively optimized using a back-propagation algorithm able to generalize to higher order perturbative contributions. Although the original application of semantic segmentation has been applied to image segmentation, its applications to NLP are somehow straight forward [168, 169].

A similar method was developed earlier by Mai, Wu and Cui and applied to improve word segmentation disambiguation in the Chinese language [170]. Main and colleagues, however, decided to use a CRF on top of a bidirectional maximum matching algorithm. Parameter estimation for the CRF was performed via maximum likelihood estimates. These ideas were further advanced by Qiu, et al [171] who used CRFs for clinical entity recognition in Chinese. Speech tagging from voice recordings was performed using a CRF devised by Khan and collaborators [172]. Even computer assisted fake news detection [173] and headline prediction [174] can be achieved using CNNs and MRFs.

Social network analysis, including online social networks, other forms of interpersonal interaction networks and even some social networks in non-human creatures, have become a relevant field of research in recent times (though the subject has been relevant in the contexts of sociology and animal behavior for decades) [175]. The analysis of social network via MRFs is becoming more and more common also. As an example, Jia and collaborators have used MRFs to infer attributes in online social network data [176]. Their model used the social network structure itself to develop a pairwise MRF. From empirical training data, the authors used the individual behaviors to learn a probability that each user has a given attribute. Then used that as an a priori probability, compute the posterior probabilities by a loopy belief propagation algorithm over the MRF, to, finally, optimizing the belief propagation algorithm by a second neighbor criteria that sparsifies the adjacency matrix. Further optimization of similar ideas was obtained by using graph convolutional networks, i.e., CNNs over CRFs [177]. Attribute inference in social network data via MRFs can also be used to improve cybersecurity algorithms [178], to learn consumer intentions [179], to study the epidemiology of depression [180] among other issues. Social networks as well as some classes of molecular interaction and ecological networks are also relevant to the development and improvement of MRF and CRF learning algorithms. This is so since often a sketch (sometimes a detailed one) of the network dependency structure is known a priori [181, 182]. This is yet another instance in which applications may nurture back the formal theory of random fields.

Graph signal theory, also called graph signal processing (GSP) is a field of signal analytics that deals with signals whose domain (as identified by a graph) is irregular [183–185]. In the context of GSP, the vertices or nodes represent probes in which the signal has been evaluated or sensed and the edges are relationships between these vertices. Data processing of the signals exploits the structure of the associated graph. GSP is often seen as an intermediate step between single channel signal processing and spatio-temporal signal analysis. The nature of the edges is determined by the relationship (spatial, contextual, relational, etc.) between the vertices. Whenever edges are defined via a statistical dependence structure, GSP can be mapped to either an MRF or a CRF, thus allowing the use of all the tools of random field theory to perform GSP [186, 187]. The networked nature of the domain of signals embedded in a graph, allows the use of spectral graph theoretical methods for signal processing [188–190]. Conversely, correlations between features on the signals are also useful to identify the structure of the underlying graph [191, 192].

GSP has a number of relevant applications, from spatio-temporal analysis of brain data [193]; to analyze vulnerabilities in power grid data [194]; to topological data analysis [195], chemoinformatics [196] and single cell transcriptomic analysis [197], to mention but a few examples. Statistical learning techniques have also being founded on a combination of MRFs and GSP [198, 199], taking advantage of both the networked structure, the statistical dependence relationships and the temporal correlations of the signals [200–202]. Random field approaches to GSP have also been applied in the context of deep convolutional networks [203, 204], often invoking features of the underlying joint conditional probability distributions such as ergodicity [205] and stationarity [206].