Learning Random Geometric Graphs Drawn in Probabilistic Metric Spaces

generic

Learning Random Geometric Graphs Drawn in Probabilistic Metric Spaces

arXiv:2608.19082v1 [stat.ML] 19 Aug 2026

Dalia Chakrabarty† , Kangrui Wang∗ , Chuqiao Zhang† , Ye Liu‡ , ∗ Alan Turing Institute

London NW1 2DB U.K. kangrui.wang@turing.ac.uk † Department of Mathematics

University of York York YO1 7DD U.K. dalia.chakrabarty@york.ac.uk joe.zhang@brunel.ac.uk ‡ Howbery Park, Wallingford

OX10 8BA U.K. y.liu@hwallingford.com Abstract: We present a new data-driven learning of a Random Geometric Graph (RGG) of a multivariate dataset, where the graph is drawn in a probabilistic metric space. This graph learning works for generic datasets, irrespective of the type of the observables; their probability distributions; or size of the data. We identify a metric of the space that the graph is drawn in, as a probability distribution of a random variable that we introduce, namely, a variable that represents the disparity between the connectedness of two vertices of the graph, and the correlation between the two random variables that are attached to the respective vertex. It is the closed-form cdf of this disparity variable that we advance as the distance function of the host space of the learnt RGG, such that the edge exists between any two nodes, if this inter-nodal distance falls short of a chosen cutoff probability. Drawing the RGG in this probabilistic space leads to the graph being an Soft RGG, such that any edge - if it exists - exists with an identified probability. We forward a simple Rejection Sampling-based technique for learning the probability of any edge. The expected degree distribution of a vertex of this RGG is identified as local, and dependent on the inter-observable correlation matrix. If said correlation matrix is not known, it can be learnt given the data, using its closed-form posterior probability density function, that we forward. We illustrate our graph learning method by learning multiple RGGs of highly multivariate real datasets. Primary: Random graphs 60-XX; Secondary: Distance in graphs 05C12; Measures of association (correlation, canonical correlation, etc.) 62H20; Probabilistic metric spaces 54E70. Keywords and phrases: Probabilistic geometric graphs, Probabilistic metric spaces, Distance function on graphs, Bayesian inference using sampling.

1

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

2

1. Introduction Graphs of complex multivariate datasets manifest intuitive illustrations of correlation structures of data, and are of pan-disciplinary interest (Airoldi, 2007; Anandkumar et al., 2012; Bandyopadhyay and Canale, 2016; Benner et al., 2014; Carvalho and West, 2007; Whittaker, 2008). We present a new method for learning a probabilistic graph of a given multivariate dataset, using the correlation structure of this data, where the graph variable in our consideration is a Random Geometric Graph or RGG (Giles et al., 2016; Penrose, 2003, 2016), that is drawn in probabilistic metric space, (Menger, 1942; Schweizer and Sklar, 1983), s.t. the distance between a pair of nodes, is a probability distribution over positive support. We will show that such an RGG is in fact, a Soft Random Geometric Graph or SRGG (Dettmann, 2018; Penrose, 2016). We recall that the probability for an edge to exist in an SRGG is the “connection function”, that takes as its input, the inter-nodal distance, while in an RGG, an edge exists if the distance between the nodes - that straddle the edge - falls short of a cutoff. Said cutoff is a probability in the RGG drawn in a probabilistic metric space. The inter-nodal distance of this RGG is advanced as a metric in the probabilistic metric space that the graph is drawn in, and this distance is used to compute the expected degree distribution of the learnt graph, as local, and dependent on the inter-observable correlation matrix of the given multivariate dataset that comprises values of such observables. We provide a closed-form posterior probability density of this inter-observable correlation matrix. We also provide a closed-form probability of the RGG variable, computed given this correlation matrix. Our graph learning is fast, and is not constrained by the nature of the observables or that of the correlation structure of the data at hand. The size of the data is also not a constraint on our graph learning - we perform an empirical illustration of learning an RGG with a large (∼ 7000) number of nodes. We start with a preliminary recollection of the basic definitions in Section 2. In Section 3, we present the closed-form density of the inter-observable correlation matrix of the given dataset. Thereafter, we develop the closed-form probability of an edge of the RGG variable (in Section 2.6), and identify the inter-nodal distance of this graph, as a distance function in the probabilistic metric space in which the RGG is drawn. Thus, this distance function is shown to be a probability distribution over a positive support (Section 2.5). This RGG is shown to be an SRGG in Section 2.10, where the connection function of this SRGG is shown to be a function of our identified inter-nodal distance. We compute the expected degree distribution of our learnt RGG in Section 2.8 and subsequently, declare the host probabilistic metric space in which we draw the RGG. The technique for learning the RGG is delineated in Section 2.6. Then in Section 4 we present the Bayesian inference that is undertaken to learn the graph alone - when the correlation of the given dataset is known; else, the graph and this correlation are learnt simultaneously with a Metropolis-with-a2-block-update. Learning of the RGG of a cuboidally-shaped data is discussed

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

3

in Section 5. Thereafter, we present an illustration on real data - to learn the large network of diseases that afflict humans, correlated by phenotypes (in Section 6). Empirical illustration on simulated data is included in the Supplement, which includes demonstration of the effect of noise in the data on learning the RGG, 2. Preliminaries In this section we provide the background on concepts that we use in the learning of the graph of a multivariate data set. 2.1. RGGs and SRGGs We seek a graph variable defined on vertex set V = {1, 2, . . . , p}, with the random variable Xi attached to the i-th node, ∀i ∈ V and ∀Xi ∈ X ⊆ R. Edges form independently in this graph and there are no self-loops. Let Gi,j = gi,j ∈ {0, 1} be the edge variable between the i-th and jth nodes, i, j ∈ [M ], where [M ] := {i, j : i < j; i, j ∈ V }. Definition 2.1. In an RGG, Gi,j ∼Poisson Point Process with its realisation gi,j = 1 ⇐⇒ D(Xi , Xj ) < τ , for a chosen cutoff τ ≥ 0, where distance between the i-th and j-th nodes is D(Xi , Xj ), with D : X × X −→ R≥0 . We note that Gi,j is the binary random variable that we assign a value of 1 to, if the i-th and j-th nodes are connected by an edge in the graph; if the edge does not exist, then we set Gi,j = 0. Definition 2.2. An SRGG defined on vertex set V is s.t. Pr(Gi,j = 1) is given by the connection function ϕ(d), where distance function between the i-th and jth nodes is d. Thus, ϕ : R≥0 −→ [0, 1] is the connection function of this SRGG, i.e. edge exists between the i-th and j-th nodes in this SRGG, with probability ϕ(d). 2.2. Partial correlation matrix Let the given multivariate dataset D, comprise n observations of each of p ob(i) (i) (i) servables X1 , X2 , . . . , Xp , such that the i-th row of D is: x1 , x2 , . . . , xp , where i ∈ {1, 2, . . . , n}. We standardise data on the j-th observable, by the estimated (1) (2) (n) mean and standard deviation of the sample {xj , xj , . . . , xj }, s.t. the dataset comprising n values of each of the standardised observables, is DS . (Later in Section 5 we will discuss the general case of X j a d-dimensional random variable, ∀j ∈ {1, . . . , p}). Here j ∈ {1, . . . , p}. Then the inter-observable correlation (S) matrix of the data DS is ΣC = [corr(Xi , Xj )].

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

4

 −1 (S) Definition 2.3. For the precision matrix Ψ = [ψi,j ] := ΣC , the matrix of absolute partial correlations is R = [ρi,j ], with ψi,j ρi,j = − p , ψii ψjj

i ̸= j, and ρii = 1 for i = j.

(2.1)

Here we want to learn the graph of dataset D, using the known absolute partial correlation ρi,j between Xi and Xj ∀i, j ∈ [M ]. If however, we wish to learn the graph given the known (absolute) correlation |corr(Xi , Xj )| between this variable pair - we only have to replace ρi,j in the discussion below, with |corr(Xi , Xj )|, ∀i, j ∈ [M ]. Keeping this in mind, in the discussions below, we will include the possibilities of graph learning given the correlation or the partial correlation, by referring to “(partial) correlation” throughout, unless otherwise stated. Below, the developed methodology considers that the (partial) correlation matrix of the considered dataset is known. 2.3. RGG drawn in probabilistic metric space: an SRGG We construct an RGG variable in a probabilistic metric space. This renders any edge of this graph a probability distribution (over positive support). Then the inter-nodal distance is anticipated to be a distribution with positive support, and this distribution will be shown to be conditional on the observation that is the absolute partial correlation ρi,j between Xi and Xj . To reflect this conditioning, we update the notation for distance between i-th and j-th nodes to D(Xi , Xj ; ρi,j ), which we will sometimes abbreviate to Di,j . We learn edge Gi,j given the value ρi,j of the absolute partial correlation variable Ri,j , between Xi and Xj , i, j ∈ [M ]. Thus, Ri,j = ρi,j ∈ [0, 1]. In this section we discuss the learning of the graph of dataset D, the (partial) correlation matrix of which we treat as known. As a supplement to this current agenda, later in Section 3, we will discuss how we can learn the correlation matrix of the given dataset by developing the closed-form posterior probability density of the correlation matrix. In Section 2.10, we will show that by constructing the random graph as an RGG drawn in a probability metric space, the graph is rendered an SRGG. We will prove this using the identified distance function in the space in which we draw the graph, and learn using the (partial) correlation structure of the given multivariate dataset. 2.4. Disparity between two variables We embed the graph variable in a probabilistic metric space (Menger, 1942), where in such a space, to any pair of points, we assign a probability distribution over positive support. Thus, to the i, j-th pair of nodes in the graph, we assign a probability distribution FS(Xi ,Xj ) (s) of the disparity variable S(Xi , Xj ) that is introduced in Definition 2.4.

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

5

Definition 2.4. The disparity random variable is defined as S(Xi , Xj ) ≡ Si,j := |Gi,j − Ri,j |. In our work, disparity Si,j measures the difference between the known (partial) correlation between Xi , Xj , and the connectedness of the i-th and j-th nodes of the graph, where the absolute partial correlation variable is Ri,j , and the binary edge variable Gi,j represents the connectedness between the i-th and j-th nodes. This holds ∀i, j ∈ [M ]. Then the disparity Si,j = s ∈ [0, 1], ∀i, j ∈ [M ]. Given the sparse information that we have on the behaviour of the disparity variable Si,j , we recall (using its value s := |gi,j − ρi,j |) that — if the given absolute value ρi,j of the (partial) correlation between Xi and Xj is high, (i.e. close to 1), then the probability for the edge to exist between the i-th and j-th nodes is high; probability for Gi,j to attain the value 0 then, is low. Thus, the probability for Si,j to attain a low value is higher than for disparity to attain a higher value. — again, if the given ρi,j is close to 0, the probability for Gi,j to be 0 is high, while that for Gi,j = 1 is low. Thus, again, we find that probability of Si,j at a low s is higher than that at a higher s. Using this intuition, we model the probability density of Si,j as a decreasing function of s ∈ [0, 1]. However, the scale with which this density decreases with s, is an unknown, i.e. this scale is a real random variable. We denote such a squared scale, (or a variance variable) to be υi,j . There is no available observation to further constrain the model of the joint density of Si,j and υi,j , than what we have motivated above. Hence we model this joint density of the disparity variable Si,j and the variance variable υi,j , at Si,j = s and υi,j = υ as:  2 s 1 fSi,j ,υi,j (s, υ) = K √ exp − , if s ∈ [0, 1]; v ∈ R>0 2υ 2πυ = 0 otherwise, ∀i, j ∈ [M ], where K > 0 is a constant that we compute below. Also, lack of further information motivates our assumption that the prior on the variance variable υi,j is Uniform[0,1]. The above holds ∀i, j ∈ [M ]. For other models of this joint density, Bayesian inference on the edge variable is likely to converge in mean. Keeping in mind our ulterior aim of learning the edge given the known interobservable correlation, we will need to marginalise out the variance variable from this modelled joint density of variance and disparity, (where the latter disparity bears information on the sought edge). Since we assume υ ∼ U nif orm[0, 1], the marginal density of the disparity

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

6

variable Si,j is Z1 fSi,j (s) =

√

K

 2 1 s exp − dυ, s ∈ [0, 1], 2υ 2πυ

0

= 0 otherwise, ∀i, j ∈ [M ], where K > 0 is a global scale. Then given that s ∈ [0, 1], this marginalisation over the variance parameter gives   2 KsΓ − 12 , s2 √ fSi,j (s) = 2 π "r #   2  2 s s = K −s , (2.2) exp − + s erf √ π 2 2 for s ∈ [0, 1]; else fSi,j (s) = 0. This holds ∀i, j ∈ [M ]. Since fS (s) is a density, the constant K is s.t "Z ! #−1 r  2   1 2 s s exp − K = + s erf √ − s ds π 2 2 0  −1   1 1 1 − +√ = erf √ 2 2 2πe

(2.3)

We will show below that when the absolute partial correlation between Xi and Xj is known to be ρi,j , the distance between the i-th and j-th points in this host space is D(Xi , Xj ; ρi,j ) = FS(Xi ,Xj ) (s), which is the cdf of the disparity computed at a value s := |gi,j − ρi,j |. We will compute said cdf by invoking the aforementioned pdf of the disparity presented in Equation 2.2. Then, given the multivariate dataset (that leads to the matrix R = [ρi,j ] of absolute inter-observable (partial) correlations), the learnt random graph variable is GS,V (R, τ ) that is defined on vertex set V s.t. the edge exists between the i-th and j-th nodes ⇐⇒ D(Xi , Xj ; ρi,j ) = FS(Xi ,Xj ) (s) ≤ τ , for a chosen cutoff probability τ . 2.5. cdf of Si,j and distance D(Xi , Xj ; ρi,j ) Definition 2.5. The density fSi,j (s) of Si,j (in Equation 2.2) is used to write the cdf FSi,j (s) of Si,j , computed at the disparity value s: "Z ! # r  2   s 2 u u FSi,j (s) = K exp − + u erf √ − u du π 2 2 0    2  s2 (s + 1) s s 2 − (2.4) =K erf √ + √ exp(−s /2) , 2 2 2 2π where K is defined in Equation 2.3. This holds ∀i, j ∈ [M ].

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

7

Remark 2.1. We see from Equation 2.4 that FSi,j (0) = 0 and FSi,j (1) = 1, given the definition of K in Equation 2.3. Also, FSi,j (s) is monotonically increasing with s ∈[0,1]; continuous in [0,1]; and one-to-one in [0,1]. This holds ∀i, j, ∈ [M ]. Next we check if D(Xi , Xj ; ρi,j ) = FSi,j (s) is a distance function. Proposition 2.1. D(Xi , Xj ; ρi,j ) = FSi,j (s) for s = |gi,j − ρi,j |, is the distance between the i-th and j-th nodes, given that the absolute partial√correlation be
tween Xi and Xj is ρi,j . Here FSi,j (s) = K ((s2 + 1)/2) erf s/ 2 − s2 /2 √  +(s/ 2π) exp(−s2 /2) Proof. We condense the notation D(Xi , Xj ; ρi,j ) to Di,j . We prove that Di,j given by FSi,j (s) that is stated in Equation 2.4 is a distance function by proving that it is non-negative; symmetric; is 0 ⇐⇒ Xi = Xj ; and obeys the triangle rule. — Non-negative: Di,j is non-negative. This follows since Di,j = FSi,j (s) is a cdf, ∀s ∈ [0, 1], ∀i, j ∈ [M ]. — Symmetric: Di,j = Dj,i by definition of Si,j := |Gi,j − ρi,j |. — Di,j = FSi,j (s) = 0 if s = 0, from Equation 2.4, ∀i, j ∈ [M ]. Again, if D(s) = FSi,j (s) = 0, it implies s = 0 since FSi,j (·) is one-to-one and monotonic increasing in [0,1]. Therefore, D(Si,j = s) = 0 ⇐⇒ s = 0, ∀i, j ∈ [M ]. — {Si,j }i≤j;i,j∈V is a set comprising identically distributed variables, each distributed as FSi,j (·) as given in Equation 2.4, ∀i, j ∈ [M ]. Then for Di,j = FSi,j (s), the triangle rule trivially holds ∀s ∈ [0, 1]: FSi,j (s) + FSj,k (s) ≥ FSi,k (s) i.e. Di,j + Dj,k ≥ Di,k , ∀i, j ∈ [M ], j < k; j, k ∈ {1, . . . , p}; i < k; i, k ∈ {1, . . . , p}.

Thus, in the RGG drawn in a probabilistic metric space, the distance function D(·, ·; ·) is the probability distribution of the disparity variable, as defined in Definition 2.5. 2.6. Learning the RGG at the given (partial) correlation matrix & probability of the graph Since we develop the graph learning in the information paradigm that R = [ρi,j ] is known, the disparity variable between Xi and Xj reduces to S(Xi , Xj ) = |Gi,j − ρi,j |, ∀i, j ∈ [M ]. At the given ρi,j , the distance Di,j between the i-th and j-th nodes is FSi,j (s), with s = |gi,j − ρi,j |, for gi,j ∈ {0, 1}, ∀i, j ∈ [M ]. Remark 2.2. In the RGG, we expect that Di,j < τ , would imply Gi,j = 1. However, this assignment of a value to Gi,j then raises a concern: there appears a circularity in our assignment of the value gi,j to the edge variable Gi,j , using the gi,j -dependent inter-nodal distance Di,j .

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

8

We avoid any such circularity, by stating that if at the given ρi,j , Di,j = FSi,j (s = |gi,j − ρi,j |) < τ =⇒ Gi,j = gi,j ; else Gi,j = 1 − gi,j in the learnt RGG, ∀i, j ∈ [M ]. Thus, in the RGG GS,V (R, τ ), we will compute Di,j = FSi,j (s = 1 − ρi,j ) and check if this computed distance falls short of a chosen τ , or not. If it does, then the edge between the i-th and j-th nodes exists in the RGG; otherwise it does not. At the same time, if this edge exists, we learn the probability with which this edge exists. This then reflects the SRGG signature of the random graph that we learn, in addition to the graph being an RGG (drawn in a probabilstic space). In order to identify the probability with which an edge exists, the following is undertaken. At the given ρi,j , we use Rejection Sampling to learn the edge variable that is straddled by the i-th and j-th nodes of the sought RGG. Such sampling will need to be done from the probability of the edge variable Gi,j given the known value ρi,j of the variable Ri,j that represents the absolute partial correlation between Xi and Xj . Here, Z 1 fSi,j (s = |gi,j − ρi,j |) = Pr(Gi,j = gi,j |ρi,j )fRi,j (ρi,j )/

fRi,j (u)du, u=0

which reduces to fSi,j (s = |gi,j − ρi,j |) = Ci,j Pr(Gi,j = gi,j |ρi,j ), given that the R1 density of Ri,j is known, and ρi,j is known, s.t. Ci,j := fRi,j (ρi,j )/[ u=0 fRi,j (u)du] is a known constant. Then summing over all values of Gi,j , at the known absolute (partial) correlation p Xi and Xj , the constant Ci,j = K[(1 √
− √
ρi,j between ρi,j ) erf (1 − ρi,j )/ 2 −(1−ρi,j )+ 2/π exp(−(1−ρi,j )2 /2)+ρi,j erf ρi,j / 2 − p ρi,j + 2/π exp(−(ρi,j )2 /2)] where K is defined in Equation 2.3. Thus, Ci,j is known, and hence, Pr(Gi,j = gi,j |ρi,j ) is known in a closed-form way, ∀i, j ∈ [M ]. We will draw ni,j number of samples of the edge variable Gi,j from its probability mass function conditional on the known pi,j = gi,j |ρi,j ) = √
ρi,j . But Pr(G fSi,j (s)/Ci,j = K[|gi,j −ρi,j | erf |gi,j − ρi,j |/ 2 −|gi,j −ρi,j |+ 2/π exp(−(gi,j − ρi,j )2 /2)]/Ci,j . It follows that Pr(Gi,j = gi,j |ρi,j ) = fSi,j (s)/Ci,j ≤ 1 =⇒ fSi,j (s)/K ≤ Ci,j /K, for the positive K. In fact, definitions of fSi,j (·), K and Ci,j show that fSi,j (s)/K < Ci,j /K ∀ρi,j ∈ [0, 1], ∀i, j ∈ [M ]. Hence, in our implementation of Rejection Sampling, we draw ni,j samples of Gi,j from the proposal density that we choose to be the Bernoulli(0.5) density, and scale this by the constant Ci,j√ /K, s.t. Ci,j /K envelopes the target p
fSi,j (s)/K = (|gi,j − ρi,j | erf |gi,j − ρi,j |/ 2 − |gi,j − ρi,j | + 2/π exp(−(gi,j − ρi,j )2 /2))/K. (⋆) (⋆) We accept the r-th proposed sample gi,j if fSi,j (|gi,j −ρi,j |)/K ≥ (Ci,j /K)u,   p √ (⋆) (⋆) (⋆) (⋆) i.e. if |gi,j −ρi,j | erf |gi,j − ρi,j |/ 2 −|gi,j −ρi,j |+ 2/π exp(−(gi,j −ρi,j )2 /2) ≥ Ci,j u, where u is the value of U ∼ U nif orm[0, 1]. Then the r-th sample is

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces (r)

(⋆)

9

(⋆)

denoted gi,j = gi,j . If on the other hand, gi,j is not accepted, then we set (r)

(⋆)

gi,j = 1 − gi,j . (R)

(1)

(2)

(n

)

We collate these samples into the set P i,j := {gi,j , gi,j , . . . , gi,ji,j }, to define (R)

the relative frequency of Gi,j =1 in this set P i,j as: nP i,j

(k)

gi,j

wi,j := k=1 ni,j

.

Typically ni,j = n, ∀i, j ∈ [M ] and we have used n ∼ 103 to 104 in undertaken applications. Thus, the sampling of gi,j at the known ρi,j is undertaken using Rejection Sampling. The algorithm used to undertake the sampling is presented in Section 4. Remark 2.3. We note that the relative frequency wi,j of achieving Gi,j = 1 in our sample generated via Rejection Sampling using fSi,j (s), is an approximation of Pr(Gi,j = 1|ρi,j ), i.e. wi,j ≈ Pr(Gi,j = 1|ρi,j ). Then at the chosen τ ∈ [0, 1], if FSi,j (s = 1 − ρi,j ) < τ , Gi,j = 1, and this edge between the i-th and j-th nodes exists in GS,V (R, τ ), with probability wi,j . If FSi,j (s = 1 − ρi,j ) ≥ τ , then the edge between the i-th and j-th nodes exists with probability 0, i.e. Gi,j = 0 then. This holds ∀i, j ∈ [M ]. Thus, we notice that it is possible for us to learn a random graph with the known (partial) correlation matrix R of the dataset D, without imposing any thresholding on any edge, i.e. at τ = 1. At the same time, our graph learning can be undertaken while acknowledging the geometric-ness of the graph. In particular, we identify a data-driven optimal τ in an application given a real dataset (Section 7). To compute the adjacency matrix of the RGG GS,V (R, τ ), we first compute the indicator function 1H (τ ) (FSi,j (1 − ρi,j ))), where H (τ ) := {δ : δ < τ, δ ∈ [0, 1]},

(2.5)

at a chosen τ ∈ [0, 1]. Then the adjacency matrix of the learnt RGG is A = [wi,j 1H (τ ) (FSi,j (1 − ρi,j ))].

(2.6)

We define the posterior probability of the RGG of the given dataset D, using the probability of the edges that are accepted at a chosen τ . Definition 2.6. Recalling that edges form independently in the RGG GS,V (R, τ ) that we seek of the dataset D, and that there are no self-loops in this RGG, the posterior probability of this RGG variable, given the known partial correlation

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

10

matrix R = [ρi,j ] of this dataset, is: Y π(GS,V (R, τ )|R) = Pr(Gi,j = 1|ρi,j )1H (τ ) (FSi,j (1 − ρi,j )) (2.7) i,j∈[M ]

≈

Y

wi,j 1H (τ ) (FSi,j (1 − ρi,j ))

i,j∈[M ]

where the set H (τ ) is defined in Equation 2.5. 2.7. Probability metric space that the RGG is drawn in Definition 2.7. In our work, the triple {X , D(·), ∆} is the probabilistic metric space that the RGG GS,V (R, τ ) is drawn in. Here the distance function D(Xi , Xj ; ρi,j ) ≡ Di,j = FSi,j (s) is defined in Proposition 2.1 for the cdf FSi,j (·) ∈ F+ , ∀Xi , Xj ∈ X , where F+ is the set of distributions over positive support. Lastly, for a given s ∈ [0, 1], a triangle function ∆ can be defined as the binary operation: ∆(FSj,k (s), FSi,k (s)) := FSj,k (s) + FSi,k (s) ≥ FSk,j (s), since {Sa,b }a,b∈[M ] is a set of identically distributed random variables. Hence, ∆(·, ·) is commutative, associative and takes FSi,i (·) = 0 as its identity. Here Si,j ≡ S(Xi , Xj ) and the above holds ∀Xi , Xj , Xk ∈ X . i, j, k ∈ V . 2.8. Degree distribution Vertex set of RGG GS,R (V , τ ) is V = {1, . . . , p}. We assume that the density at point x is f (x). For any i ∈ V , define (open or close, but bound) Borel measurable ball Bi,a , centred at the i-th node, with radius a, where a ∈ [0, 1], given that radius of a ball in the space of the RGG, is a probability. Let random variable N (Bi,a ) be the number of elements of vertex set V that lie inside Bi,a , as connected to the i-th node at the centre of Bi,a . We√recall
that if Di,j = FSi,j (s = 1 − ρi,j ) = K[(((1 − ρi,j )2 + 1)/2) erf (1 − ρi,j )/ 2 − √ (1 − ρi,j )2 /2+ ((1 − ρi,j )/ 2π) exp(−(1 − ρi,j )2 /2)] < τ , Gi,j = 1 in our RGG, at the chosen cutoff τ . At the given partial correlation matrix R (of the dataset for which this RGG is learnt), the number density of edges that will form between the i-th node at the centre of the ball Bi,a , and a point within distance x ≤ a within this ball, is X (R) ni (x) := f (x)H(τ − Di,j ), where j∈C i (x)

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

11

C i (x) := {j : j ̸= i, j ∈ V , Di,j ≤ x}. Here H(·, ·) is a Heaviside function. Then the expectation of the degree of the i-th node in the leant RGG is: Z a (R) (R) E[Ni ] = dx 2πxni (x) (2.8) x=0

≡

2 λ(i) a,τ (R)πa

Z a

where 2 λ(i) τ (R) :=

(R)

dx xni

(x)

x=0

a2

,

s.t. 2 E[N (Bi,a )] = λ(i) τ (R)πa

for the given a. If τ < a, there is no contribution to the last integral from x ∈ [τ, a]; else, the expectation is contributed to by all x ∈ [0, a]. Thus, the number density of edges within a bound region of radius a, centred (i) at point Xi in X , is given as the location-dependent λτ (R) for a given partial correlation matrix R, at a chosen threshold cutoff probability τ . Even if the density f (·) is homogeneous, the RGG variable is generated by an inhomogeneous point process. We note that the local intensity of the generative process is a function of the threshold τ , and the correlation matrix of the dataset for which the RGG is learnt. For a given Xi , if Xj and Xj ′ are s.t. disparity Si,j = s is less than Si,j ′ = s′ , ′ ′ (s ). Then Di,j < Di,j ′ , s.t. the i-th node is more cdf FSi,j (s) is lower than FSi,j likely to be connected to the j-th node, than to the j ′ -th node, at a given τ . Thus, variables with low values of mutual disparity, could form clusters within the RGG of a the given dataset. Remark 2.4. |C i (x)| is approximately the same for all points that lie in the bulk of such a cluster - as is f (x). Thus, nodes that lie in the bulk of a cluster will share similar values of their degrees. On the other hand, for the i-th node that lies near the edge of a cluster, |C i (x)| is lower than if this node were in the bulk. Then, in general, the degree of a node near the edge of a cluster is lower than that in its bulk. Within a “small” distance inwards from the edge of a cluster, the degree distribution is likely to show an increasing trend, but how quickly this trend shows up, i.e. quantification of the aforesaid “small”ness, is affected by the choice of the cutoff τ used in the learning of the RGG of a given dataset. 2.9. Choosing the cutoff probability τ Since the learnt graph GS,V (R, τ ) is a random graph-valued variable, we can define its posterior probability π(GS,V (R, τ )|R) given the partial inter-observable

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

12

correlation matrix R of the given dataset. This then allows for an organic identification of the optimal cutoff, at which the most-robust RGG is realised, where said robustness is to changes in τ : the τ at which the rate of change of the RGG posterior is minimised, i.e. d(π(GS,V (R, τ )|R))/dτ is minimised, is the optimal cutoff that produces the most-robust RGG. Since we typically compute logarithm of the posterior probability of the RGG, the optimal τ is implemented as the cutoff, at which the rate of change (with τ ) of the logarithm of this posterior probability is minimised. 2.10. RGG in probabilistic metric space is an SRGG We draw an RGG in the probabilistic metric space, and an edge is accepted if the inter-nodal distance falls below a chosen cutoff probability τ , where the accepted edge exists with an identified probability. Proposition 2.2. The RGG GS,V (R, τ ), drawn in a probabilistic metric space is an SRGG. Proof. At the known partial correlation matrix R = [ρi,j ], and for s = |gi,j − ρi,j |, sampling Gi,j from fSi,j (s), ni,j times, produces the set (R

(t) n

i,j P i,j = {gi,j }t=1 ,

Pni,j (t) with wi,j := t=1 gi,j /ni,j . Then edge probability Pr(Gi,j = 1|ρi,j ) ≈ wi,j , ∀i, j, ∈ [M ]. We have seen above that inter-nodal distance Di,j ≡ FSi,j (s = 1 − ρi,j ) < τ =⇒ Gi,j = 1; else Gi,j = 0. Then we define the function ϕ : [0, 1] −→ [0, 1] s.t. at a chosen τ ∈ [0, 1], the edge between the i-th and j-th nodes exists with probability ϕ(Di,j ) := wi,j 1H (τ ) (Di,j ), where H (τ ) := {δ : δ < τ, δ ∈ [0, 1]}. This holds ∀i, j ∈ [M ]. Thus, in this RGG variable, we have identified the probability with which an edge exists, where said probability is the function ϕ(·) of the distance between the nodes that this edge connects. An SRGG is a random graph in which the probability with which an edge exists, is a function of the inter-nodal distance. Hence this RGG drawn in a probabilistic metric space, is an SRGG. 3. Posterior density of the inter-observable correlation matrix In this section we present a closed-form posterior probability density of the inter-observable correlation matrix, where n observations of each of p observables X1 , X2 , . . . , Xp comprise the dataset. Later in Section 5, we discuss learning the RGG for vector-valued observables. We standardise the data on the

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

13

variable Xi using the mean and standard deviation of the sample comprising (1) (2) (n) the n observations, i.e. the sample {xi , xi , . . . , xi }. Such standardisation is undertaken for ∀i ∈ {1, . . . , p}. While values of X1 , . . . , Xp comprise the dataset D, data on each of p standardised observables Z1 , . . . , Zp , comprise the dataset Ds . The observable X m with standardised components, is denoted Z m , s.t. we (1) (2) (p) can define vector Z m = (Zm , Zm , . . . , Zm )T , ∀m ∈ {1, . . . , n}. Following earlier work, (Gruber and West, 2016; Ni et al., 2017; Wang and West, 2009), we model the vector-valued observable output at any observational instance, with a vector-valued Gaussian Process. This implies that the joint density of the n realisations of this observable, namely Z 1 , . . . , Z n , is a matrix Nor(S) mal density MN (0, Σ(S) m , ΣC ). Here, the undertaken standardisation implies that the mean matrix of this density is a null matrix, while its two covariance matrices are rendered correlation matrices. These two correlation matrices are: ′ the d × d-dimensional inter-observable matrix Σ(S) m , the m, m -th element of ′ which is the correlation between Z m and Z m′ , ∀m, m ∈ {1, . . . , n}; and the (S) p × p-dimensional inter-variable correlation matrix ΣC , the i, j-th element of which is the correlation between Zi and Zj , ∀i, j ∈ {1, . . . , p}. (S)

Theorem 3.1. When the prior on ΣC is Uniform, the joint posterior proba(S) bility density of the correlation matrices ΣC and Σ(S) m , given the standardised (S) data DS , can be marginalised over all Σm , to yield the marginal posterior pdf of the inter-observable correlation matrix given data DS , as: (S)

(S)

π(ΣC |DS ) ∝ ΣC

−p/2

DS (ΣC )−1 (DS )T (S)

− m+1 2

,

α d (S) (S) , with α = − − 1, where prior on Σ(S) m is non-informative: π0 (Σm ) = Σm 2 (S) and ΣC is assumed invertible.

The proof follows from marginalising over all Σ(S) m , using the uninformative prior on this matrix, as per the theorem statement, and recalling that the poste(S) (S) rior pdf π(ΣC |DS ) is obtained for a Uniform prior on ΣC , s.t. it is proportional (S) (S) (S) to the likelihood of ΣC in DS , i.e. π(ΣC |DS ) ∝ L (ΣC ; DS ). Towards the −(d+1) computation of this likelihood, we use the result that d(Σ(S) dY m ) = |Y | (S) −1 (Mathai and Pederzoli, 1997), where Y := (Σm ) . Then invoking invertibility (S) of DS (ΣC )−1 (DS )T , the result follows. (S) The posterior probability density of the correlation matrix ΣC , given the data DS , as stated in Theorem 3.1, is reminiscent of the density of the matrixvalued t-distribution, but differences exist. 4. Inference Remark 4.1. When the inter-observable correlation matrix ΣC of the given multivariate data D is known, then the probability of the RGG variable com-

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

14

puted with the known partial correlation matrix R = [ρi,j ] - which is itself computed using the known ΣC - is closed-form. We then learn Gi,j (given ρi,j ) by undertaking Rejecton Sampling from the density of the disparity variable Si,j = |Gi,j − ρi,j |, ∀i, j ∈ [M ]. In Section 4.1, we present the steps for undertaking Rejection Sampling-based learning of the RGG, given the known correlation matrix of the given data. This approach is referred to as Algorithm 1. Remark 4.2. When we wish to learn ΣC of data D, and then learn the RGG given this learnt correlation matrix, we undertake MCMC-based Bayesian inference on the RGG and the partial correlation matrix R simultaneously, within a Metropolis-with-2-block-update algorithm. We update R (by updating ΣC ) in the first block of an iteration of the undertaken MCMC chain, and then in the second block of the iteration, update the RGG variable using the updated R. However, there is no feedback from this second block to the first block of the subsequent iteration. Thus, this inferential scheme differs from Metropolis-within-Gibbs. Hence we refer to it as Metropolis-with-2-block-update. The algorithm for the learning of the correlation matrix and the RGG using MCMC, is shown in Algorithm 2. A schematic flow of the undertaken inference is in Figure 1. In the applications we undertake, we do not possess strong prior information on the unknowns; in such situations, the priors used are retained as weak. This is manifest in high variances used for the prior probability density on any unknown. A graph edge parameter has a prior of Bernoulli(0.5). In applications that are blessed with more information, stronger priors can be used. 4.1. Algorithm 1: implementation of Rejection Sampling At a chosen τ ∈ [0, 1], and knowing that the (partial) correlation matrix of the data D is R = [ρi,j ], the following algorithm is used to perform Rejection Sampling to learn the RGG of D. 1. In the k-th trial of Rejection Sampling, we propose the edge between the (⋆,k) i-th and j-th nodes as Gi,j = gi,j , from a proposal distribution p(gi,j ) = Bernoulli(0.5), and scale it with the constant Ci,j /K that envelopes the (⋆,k) target function fSi,j (s = |gi,j −ρi,j |)/K. (See section 2.6 for details). This is equivalent to performing Rejection Sampling from Pr(Gi,j = gi,j |ρi,j ) = fSi,j (s)/Ci,j , using the proposal p(gi,j ) = Bernoulli(0.5), ∀ρi,j ∈ [0, 1]. We undertake this ∀i, j ∈ [M ]. (⋆,k) 2. Then (fSi,j (s = |gi,j −ρi,j |)/K)/((Ci,j /K))) ≥ u, implies that we accept (k)

(⋆,k)

(k)

(⋆,k)

the k-th sample as gi,j = gi,j ; else we set gi,j = 1 − gi,j . Here, U = u, where U ∼ Uniform[0, 1]. We undertake the Rejection Sampling till we (R) have generated ni,j samples that we use to populate the set P i,j as:

Undertake Rejection Samof edge Gi,j from Chakrabarty et. al/Learning RGGs in Probabilistic pling Spaces

Need to learn RGG variable only - using Rejection Sample

Compute absolute partial correlation matrix R = [ρi,j ] given Σc , if learning RGG using R; else skip this step, and use ci,j = corr(Xi , Xj ) in place of ρi,j .

yes

Inter-observable correlation matrix Σc = [ci,j ] - where ci,j = corr(Xi , Xj ) is known, or already learnt/estimated given data D

(t)

no

In the 1st block of the t(⋆,t) th iteration, propose ci,j (t−1)

Need to learn both ΣC and RGG, using MCMCby-2-block-update

from TN (ci,j , vi,j , −1, 1), for chosen variance vi,j . (⋆,t) Accept/reject ci,j ∀i, j ∈ [M ] using posterior pdf of Σc , computed at current & proposed values of this matrix

15

Pr(Gi,j = gi,j |Ri,j ρi,j ), which is ∝ fSi,j (s) given that density of Ri,j is known and edge is learnt at known value ρi,j of Ri,j . Here the pdf fSi,j (s) of Si,j = |Gi,j − ρi,j | is computed at s = |gi,j − ρi,j |. Generate Ni,j samples of Gi,j s.t. t-th sam(t) pled edge is gi,j . Then Pr(Gi,j = 1|ρi,j ) ≈ wi,j = Pni,j (t) t=1 gi,j /ni,j . If internodal distance Di,j = FSi,j (s = 1 − ρi,j ) < τ , edge between i-th and j-th nodes exists in the RGG with probability wi,j ; else Gi,j = 0. This holds ∀i, j ∈ [M ]

(t)

Compute ρi,j using ci,j . In the 2nd block of the t-th iteration, propose (⋆,t) (t) gi,j ∼ Bernoulli(ρi,j ). (⋆,t)

Accept/Reject gi,j using the edge posterior probability s.t. sampled edge (t) is gi,j , ∀i, j ∈ [M ]. Post MCMC, edge between ith and j-th nodes is found to exist with probability Pni,j (t) ≈ wi,j = t=1 gi,j /ni,j , if FSi,j (s = 1 − ρi,j ) < τ ; else Gi,j = 0

Fig 1. Flow depicting inference undertaken, when the correlation matrix of the given dataset is known or to-be-learnt.

(R)

(1)

(n

)

P i,j = {gi,j , . . . , gi,ji,j }. The above is undertaken for ∀i, j ∈ [M ], with ni,j = n ∼ 103 or 104 . P(ni,j (k) 3. Compute wi,j = k=1 gi,j /ni,j , ∀i, j ∈ [M ], where Pr(Gi,j = 1|ρi,j ) ≈ wi,j . 4. Compute FSi,j (s = 1 − ρi,j ) (using Equation 2.4). FSi,j (s = 1 − ρi,j ) < τ =⇒ edge between the i-th and j-th nodes exists with probability Pr(Gi,j = 1|ρi,j ) ≈ wi,j . On the other hand, FSi,j (s = 1 − ρi,j ) ≥ τ =⇒ Gi,j = 0. This is undertaken ∀i, j ∈ [M ]. 5. Thus, learning the value of Gi,j ∀i, j, ∈ [M ], will define the realisation of the RGG variable GS,V (R, τ ) at a chosen τ , given the known partial

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

16

correlation matrix R (of a given dataset D), for which the edge set is; E = {gi,j : H(τ − FSi,j (1 − ρi,j )), i, j ∈ [M ]}, where if the edge exists between the i-th and j-th nodes, it will do so, with probability approximated by wi,j , ∀i, j ∈ [M ]. 4.2. Algorithm 2: implementation of Metropolis-with-2-block-update Updating of the correlation matrix is performed in the first block, and each edge (S) variable is updated in the second block. For the correlation matrix ΣC = [ci,j ] of the given data D, in the t-th iteration, Cij is proposed from a Truncated Normal density (TN ) that is left truncated at -1 and right truncated at 1, with (t−1) an experimentally chosen variance (vi,j ), and the current value ci,j as the proposal mean. When we learn the graph with partial correlations - instead of (t) correlations - we update the partial correlation between Xi and Xj to ρi,j using the value of Ci,j that is updated in the first block of the t-th iteration. Then Gi,j is proposed in the second block of this iteration, from a Bernoulli pmf (t) with parameter ρi,j . We discuss the algorithm below. 1. In the first block of t-th iteration of the MCMC chain, the i, j-th element of the correlation matrix is proposed as: (t−1)

Ci,j = c⋆,t i,j ∼ TN (ci,j

, vi,j , −1, 1), ∀i, j ∈ [M ].

Let the proposed Ci,j value be the i, j-th element of the proposed inter(⋆,t) observable correlation matrix ΣC , i.e. (⋆,t)

ΣC

= [c⋆,t i,j ]. (t)

2. This proposed covariance matrix may or may not be accepted, (as ΣC = (t) [ci,j ]), depending on whether the acceptance ratio leads to acceptance or rejection. Acceptance occurs if (⋆,t)

π(ΣC

Q

|DS ) ×

(t−1)

pdf of TN (ci,j

, vi,j , −1, 1)

i,j∈[M ] (t−1)

π(ΣC

|DS ) ×

Q

(⋆,t)

pdf of TN (ci,j , vi,j , −1, 1)

≥ u,

i,j∈[M ] (t−1)

where U = u, with U ∼ U nif orm[0, 1]. Here ΣC

(t−1)

= [ci,j

] and poste-

(·) rior probability density of a correlation matrix is π(ΣC |DS ) as given in (t) (⋆,t) (t) (t−1) Theorem 3.1. If accepted, we set ΣC = ΣC ; else, ΣC = ΣC . Here,

DS is the dataset that results from standardisation of D. (S) 3. Theorem 3.1 gives the posterior of ΣC for Uniform (U nif orm[−1, 1]) (S) prior on ci,j , ∀i, j ∈ {1, . . . , p}. Thus the likelihood of ΣC in the data, is

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

17

given by this theorem. If in an application, a more informative prior on ci,j - than Uniform - is available, then such a prior is multiplied with the likelihood, to provide the posterior, that is in turn used in the acceptance ratio. (S) 4. Then using this updated ΣC matrix, the updated partial correlation (t) matrix R(t) = [ρi,j ] is populated using the updated absolute partial cor(t)

relations {ρi,j }i,j∈V that are computed using Equation 2.1. 5. Subsequently, in the 2nd block of the t-th iteration, the RGG is updated, using the recently updated absolute partial correlation matrix R(t) . The proposed edge variable connecting the i-th to j-th vertex is (⋆,t)

Gi,j = gi,j

(t)

∼ Bernoulli(ρi,j ). (t)

We compute the density of the disparity variable Si,j = |Gi,j − ρi,j |, at the (⋆,t)

proposed edge, where this density fSi,j (s) computed at s = |gi,j (⋆,t)

(t)

− ρi,j |,

(t)

is proprtional to Pr(Gi,j = gi,j |ρi,j ). This holds ∀i, j ∈ [M ]. 6. We accept the proposed edge if the acceptance ratio exceeds or equals u that is the value of U ∼ U nif orm[0, 1], where acceptance ratio relevant (t) to the 2nd block of the t-th iteration is computed at Ri,j = ρi,j as: Q

(⋆,t)

fSi,j (s = |gi,j

(t)

− ρi,j |)

i,j∈[M ]

Q

(t−1)

fSi,j (s = |gi,j

(t−1)

− ρi,j

)

×

i,j∈[M ]

Q

(t−1)

pmf of Bernoulli(ρi,j

)

i,j∈[M ]

Q

(t)

pmf of Bernoulli(ρi,j )

i,j∈[M ] (t) (⋆,t) (t) (t−1) If accepted, we set gi,j = gi,j ; else, gi,j = gi,j , ∀i, j ∈ [M ]. Here the

density of the disparity variable is given in Equation 2.2. 7. Prior used on Gi,j is Bernoulli(0.5). Values of the edge variable Gi,j current in the n post-burnin iterations are (R) collated into the set P i,j , ∀i, j ∈ [M ]. Then out of the n samples collected in (R)

P i,j , we compute the fraction wi,j of samples that have a value 1, where wi,j estimates the sample mean of Pr(Gi,j = 1|ρij ). The edge between the i-th and j-th nodes then exists with probability Pr(Gi,j = 1|ρi,j )1H (τ ) (FSi,j (1 − ρi,j )), where the set H (τ ) = {δ : δ < τ, δ ∈ [0, 1]}. Thus, the RGG GS,V (R, τ ) is learnt, with the probability of the existent edges estimated. 4.3. When learning large networks When our interest is in learning a graphical model on a vertex set with |V | = p ≳20 - as in a large network - the computational cost of making MCMC-based

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

18

inference on the off-diagonal elements in the upper (or lower) triangle of a p × pdimensional correlation matrix, is prohibitive, thereby prohibiting computation (S) of R. Then we use the plugin estimate ĉi,j of ci,j , s.t. correlation matrix ΣC is estimated as [ĉi,j ]. Then the partial correlation matrix is estimated as R̂, using this estimated correlation matrix, as long as p ≲ 1000. Thus, such a network is learnt as the RGG GS,V (R̂, τ ) if p ≲ 1000. However, if p > 1000, (S) the cost of inverting the estimated ΣC - to achieve the corresponding partial correlation matrix - might be desired to be avoided. In that case, we will learn the RGG using the estimated correlation matrix, instead of the partial correlation (S) matrix, i.e. as GS,V (ΣC , τ ). Then we also reduce the number Niter of samples drawn from the marginal of Gi,j , depending on the value of p, when undertaking Rejection Sampling, ∀i, j ∈ [M ]. We have learnt such a network bearing > 5000 nodes, each of which is a human disease, (Section 6). In this network, any pair of such diseases is correlated by the rank of the extent of overlap of corresponding phenotypes. 5. Generalisation to learning an RGG of data on vector-valued observables (k)

(k)

(k)

Let us consider learning the graph of the dataset D = {x1 , x2 , . . . , xp }nk=1 , (·) where the observable X i = xi ∈ X ⊆ Rd , ∀i ∈ {1, . . . , p}. Here we define X i = (Xi,1 , . . . , Xi,d )T . It is also possible that the dataset D comprises a varying number of observations of X i , distinguished from that of X j , for i, j ∈ {1, 2, . . . , p}, but here we develop the learning of the graph using the same number (n) of observations of each of the p observables. Consequently, the dataset D is cuboidal in shape. In fact, we standardise the data D to the dataset Ds , as suggested in Section 3. Motivated by our RGG learning given scalar observables, we draw the RGG of Ds in the probabilistic metric space {X , D(·), ∆}, with the distance D(·) between a pair of distinct nodes, given by the cdf of the disparity between the observables attached to these nodes, conditional on the absolute (partial) correlation between these variables. The triangle function ∆ is defined in Definition 2.7. This RGG is defined on the vertex set V ′ = {(1, 1), (1, 2), . . . , (1, d), (2, 1), . . . , (2, d), . . . , (p, 1), . . . , (p, 1), . . . , (p, d)}, where the observable Xi,m is attached to the (i, m)-th node, ∀i ∈ {1, . . . , p} and ∀m ∈ {1, . . . , d}. The learnt RGG is a flat graph, s.t. edges can exist between nodes (i, m) and (i, m′ ), as well as between nodes (i, m) and (j, m′ ) for m < m′ ; m, m′ ∈ {1, . . . , d} and i, j ∈ [M ], though self-loops are disallowed in our learnt RGGs. Definition 5.1. In this exercise of learning the RGG with vector-valued vari(m,m′ ) ables attached to respective nodes, we invoke the disparity Si,j between Xi,m and Xj,m′ , i.e. the disparity between the m-th component of X i and the m′ -th component of X j , ∀m, m′ ∈ [M ] and i < j; i, j{1, . . . .d}. We invoke the defined (m,m′ ) (m,m′ ) disparity between two random variables from above, as: Si,j = |Gi,j −

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces (m,m′ )

ρi,j

(m,m′ )

|, and Si,j

19

attains values ∈ [0, 1]. Here, the edge variable that con(m,m′ )

(m,m′ )

nects the (i, m)-th node to the (j, m′ )-th node is Gi,j , and ρi,j is the known, absolute (partial) correlation between the variables attached to these two nodes. Then the distance between √ is D(Xi,m , Xj,m′ ) = √ these two nodes FS (m,m′ ) (s) = K[(s2 + 1)/2 erf(s/ 2) − s2 /2 + (s/ 2π) exp(−s2 /2)], as we i,j

have seen in Equation 2.4, with K defined in Equation 2.3. In the RGG of dataset Ds , if D(Xi,m , Xj,m′ ) < τ , the edge exists between the (i, m)-th node (to which Xi,m is attached), and the (j, m′ )-th node (to which Xj,m′ is attached). If D(Xi,m , Xj,m′ ) ≥ τ , the edge between these two nodes is absent in the RGG. This way of learning the RGG for vector-valued observables is then subject to knowledge of the correlation between variable Xi,m and variable Xj,m′ . Such correlation learning sets the learning of the RGG in this case, different from the learning of the RGG given data on scalar-valued observables. 5.1. Learning the correlation between Xi,m and Xj,m′ We now discuss learning the inter-variable given dataset Ds , in order to learn its RGG. Here Ds is built with n observations of Xi,m , for i ∈ {1, . . . , p} and m ∈ {1, . . . , d}. Hence, the pd × pd-dimensional ΣC ⊗ Σm = [corr(Xi,m , Xj,m′ )], where — the p × p-dimensional matrix of inter-observable correlations is ΣC = [corr(X i , X j )]; — the d × d-dimensional matrix of inter-component correlations is Σm = ′ [corr(W (m) , W (m ) )], — with W defined as: W (m) = (X1,m , X2,m , . . . , Xp,m )T . Thus, corr(Xi,m , Xj,m′ ) is the d(i−1)+m, d(j −1)+m′ -th element of ΣC ⊗Σm . We need to know corr(Xi,m , Xj,m′ ) to define the absolute (partial) correlation (m,m′ ) ρi,j that allows us to compute the cdf FS (m,m′ ) (s) of the disparity between i,j

Xi,m and Xj,m′ , (Equation 2.4), and thereby learn if the edge exists between the nodes that these two variables are attached to, at a chosen τ . 5.1.1. Inter-observable correlation matrix The inter-observable correlation ΣC matrix can be estimated before we start learning the RGG, or could be inferred upon using the posterior probability density that was reported earlier in Section 3. This posterior density is π(ΣC |DS ) ∝ ΣC

−p/2

DS (ΣC )−1 (DS )T

− n+1 2

.

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

20

5.1.2. Inter-component correlation matrix The m, m′ -th element of the inter-component correlation Σm can be computed by transforming a statistical distance or divergence - such as the Hellinger distance or Kullbeck-Leibler divergence - between the probability of the RGG variable learnt given the data on W (m) and the probability of the RGG given ob′ servations of W (m ) . Here, n values of W (m) comprise the dataset DW (m) := (i) (i) {(x1,m , . . . , xp,m )}ni=1 , ∀m ∈ {1, . . . , d}. (m)

(m)

Let Gi,j = gi,j ∈ {0, 1} be the edge between the (i, m)-th node that hosts the m-th component of the i-th observable, and the (j, m)-th node, (to which (m) the m-th component of the j-th observable is attached). Again, let ρi,j be the absolute (partial) correlation between these two observables. To compute the probability π(GS,V (R(m) , τ )|R(m) ) of the RGG variable GS,V (R(m) , τ ), given DW (m) , we recall from Equation 2.8 that we first need to compute the prob(m) ability with which the edge Gi,j exists in this RGG. We use Rejection Sampling to draw n samples from the unnormalised density of the disparity be(m) (m) tween variables Xi,m and Xj,m , where said density is fSi,j (m) (|gi,j −ρi,j |)/K = q   (m) (m) (m) (m) √ (m) (m) (m) (m) |gi,j −ρi,j | erf |gi,j − ρi,j |/ 2 −|gi,j −ρi,j |+ π2 exp(−(gi,j −ρi,j )2 /2), (m)

where K is defined in Equation 2.3. The fraction wi,j of the samples that is a (m)

(m)

1, is the approximation of the probability Pr(Gi,j = 1|ρi,j ). To then acknowledge the effect of the thresholding, we compute the distance D(Xi,m , Xj,m ) (m) between the nodes to which Xi,m and Xj,m are attached - at Gi,j = 1. This is the cdf of the disparity Si,j (m) between these two random variables, computed (m)

at 1 − ρi,j , and is given in Equation 2.4. If this cdf falls short of the chosen τ , (m)

(m)

then edge Gi,j = 1, with probability wi,j ; else, there is no edge between the nodes to which Xi,m and Xj,m are attached. Then Y (m) (m) π(GS,V (R(m) , τ )|R(m) ) = wi,j 1H (τ ) (FS (m) (1 − ρi,j )) (5.1) i,j

i,j∈[M ]

Similarly, the probability of the RGG variable of the data DW (m′ ) is computed. The Hellinger distance (or Kullbeck Leibler divergence): ′

δ(GS,V (R(m) , τ ), GS,V (R(m ) , τ )) ′

′

is computed between π(GS,V (R(m) , τ )|R(m) ) and π(GS,V (R(m ) , τ )|R(m ) ). ′

Definition 5.2. δ(GS,V (R(m) , τ ), GS,V (R(m ) , τ )), is computed using the set of probability values that are computed using samples (generated in ni,j trials, with Rejection Sampling) of the edge variables relevant to each RGG, where we normalise each computed probability by the maximal - out of all computed - probability values. Additionally, the inter-observable correlation matrices of datasets DW (m) and DW (m′ ) can be learnt using the closed-form posterior density of the

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

21

inter-observable correlation matrix discussed in Section 3. Then the m, m′ -th element of the inter-component correlation Σm is modelled as: ′

Σm := [exp(−δ(GS,V (R(m) , τ ), GS,V (R(m ) , τ )))].

(5.2)

5.1.3. Correlation between any two components of any two observables Knowing the d × d-dimensional Σm and the p × p-dimensional ΣC , we compute the correlation between the m-th component of the i-th observable and the m′ -th component of the j-th observable, ∀i, j ∈ V , and ∀m, m′ ∈ {1, . . . , d}: ΣC ⊗ Σm = [corr(Xi,m , Xj,m′ )], 5.2. Learning RGG given inter-observable and inter-component correlation matrices The RGG that we learn given the cuboidally-shaped dataset Ds , is GS,V ′ (ΣC ⊗ (m,m′ )

Σm , τ ). Then for (i, m), (j, m′ ) ∈ V ′ , the edge Gi,j between the (i, m)-th node and the (j, m′ )-th node exists, if the inter-nodal distance D(Xi,m , Xj,m′ ) = p √ (m,m′ ) 2 (m,m′ ) (m,m′ ) (m,m′ ) [ 2/π exp(−|1−ρi,j | /2)+|1−ρi,j | erf((|1−ρi,j |)/ 2)−|1−ρi,j |] < τ ; else this edge does not exist. As delineated above, we can again estimate the probability with which this edge exists - if at all - using the relative frequency of an edge sample to be 1, amongst samples drawn using Rejection Sampling, from the pdf of the disparity between the variables Xi,m and Xj,m′ . Thus, the RGG is (m,m′ ) learnt given the absolute correlation ρi,j between variables Xi,m and Xj,m′ , (m,m′ )

where ρi,j is the d(i−1)+m, d(j −1)+m′ -th element of ΣC ⊗Σm . We might know the inter-observable matrix ΣC , or learn it using its closed-form posterior, as delineated in Section 3. The inter-component matrix Σm might again be known, or it can be computed as a transformation - as stated in Equation 5.2 - of the distance between the probability of the RGG of the data on the m-th component of all observables, and that of the data on the m′ -th component of all observables. 6. Empirical illustration: learning the human disease-phenotype network The human disease-phenotype network was learnt by Hoehndorf et al. (2015) (HSG hereon), by considering the similarity parameter for each pair of diseases that is an element of an identified set of diseases in the Human Disease Ontology (DO), that contains information about rare and common diseases, and spans heritable, developmental, infectious and environmental diseases. Here, the “similarity parameter” between one disease and another, is computed using the ranked vectors of “normalised pointwise mutual information” (NMPI) parameters for the two diseases, where the NMPI parameter describes the relevance of a

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

Fig 2. RGG G

(St)

S,ΣC

22

(V , 0.1), representing the human disease-phenotype network that we

learn using the computed Spearman rank correlation between the rank vectors of a list of phenotypes relevant to a disease, where the phenotype ranking reflects semantic relevance of a phenotype to the disease in question (quantified by HSG as the NPMI parameter in the DDP h dataset). In our learnt RGG, τ = 0.1. Here the vertex set V has 8676 elements, but nodes with no edges are discarded from this visualised graph, resulting in 6052 diseases (nodes) and 145210 edges that are shown this figure. Diseases identified by HSG, to belong to one of the 19 given disease class, are presented above in the same colour; the colour key identifying these classes, is attached.

phenotype, to the disease in question. HSG define the NMPI parameter semantically, as the normalised number of co-occurrences of a given phenotype and a

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

23

disease in the titles and abstracts of 5 million articles in Medline. The diseasedisease pairwise semantic similarity parameters – computed using the degree of overlap in the relevance ranks of phenotypes associated with each disease – result in a similarity matrix, which HSG turn into an inter-disease network based on phenotypes. They choose from the top-ranking 0.5% of inter-disease similarity values. Phenotypes associated with diseases, and corresponding scoring functions (such as the NPMI), exist in the file “doid2hpo-fulltext.txt.gz” at http://aber-owl.net/aber-owl/diseasephe notypes. In fact, at the site http://aber-owl.net/aber- owl/diseasephenotypes/data/, HSG have uploaded all the data that they have used. The file ”doid2hpo-fulltext.txt.gz” available at this site, contains information about Ndis diseases, and the semantic relevance of each of the Npheno phenotypes to each disease, as quantified by NPMI parameter values, in addition to other scores such as t-scores and z-scores. In this file, Ndis is 8676 and Npheno is 19323. In the phenotypic similarity network between diseases that HSG report, diseases are the nodes, and the edge between two nodes exists in this undirected graph, if the similarity between the nodes (diseases) is in the highest-ranking 0.5% of the 38,688,400 similarity values. They remove all self-loops and nodes with a degree of 0. Their network is presented in http://aber-owl.net/aberowl/diseasephenotypes/network/. The “Group Selector” function on their visualisation kit, allows for the identification of 19 clusters in their disease-disease network, with each cluster corresponding to a disease-class. Total number of nodes over their identified 19 clusters, is 5059; number of edges is 65,795; average node degree≈26.2. HSG’s network then manifests a similarity-structure that is computed using available NPMI parameter values. Our interest is in learning the disease-disease network as an RGG, with each edge existent at a learnt probability. We perform such learning using the NPMI semantic-relevance data that is made available for each of the Ndis number of diseases, by HSG; so Ndis is the norm of the vertex set V of our sought RGG. We refer to this human disease-phenotype data as DDP h . Using DDP h , we first compute the correlation Sij between the i-th and j-th diseases in V , for each of which, information on the ranked (semantic) relevance of each of the Npheno phenotypes exist, in this given dataset. Upon computation of pairwise correlations, the RGG for the data DDP h is learnt. We compute the correlation between the i-th and j-th diseases in the DDP h data, (i, j ∈ {1, . . . , Ndis }, i < j), in the following way. We rank the NPMI parameter values - that indicate association between the i-th disease and each of the Npheno phenotypes - with phenotypes of highest semantic relevance to the i-th disease, assigned a rank 1. Let the rank vector of phenotypes, by semantic relevance to the i-th disease take the value ri and similarly, that for the j-th (rank) disease is rj . We compute the Spearman rank correlation sij , between vectors ri and rj ∀ i, j ∈ {1, . . . , Ndis }; i < j. Spearman rank correlation is preferred to the correlation between vectors of normalised NPMI values, since we intend to correlate the i-th disease with the j-th disease, depending on how relevant

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

24

Fig 3. Left: comparison of the relative number of nodes (diseases) that we recover in each of the 19 disease classes (in filled circles joined by solid lines), with the relative class-membership reported by HSG overplotted as open circles threaded by broken lines, Right: our computed ratios of the averaged intra-class to inter-class variance for each of the 19 classes, shown in filled circles; the ROC Area Under Curve values reported by HSG for each class, is overplotted as open circles joined by broken lines. The disease class indices are assigned values 1 to 19; these are the following disease classes respectively: cellular proliferation diseases, integumentary diseases, diseases of the nervous system, genetic diseases, diseases of metabolism, diseases by infectious agents, diseases of mental health, physical disorders, diseases of the reproductive system, of the immune system, of the respiratory system, of the musculoskeletal system, syndromes, gastrointestinal diseases,cardiovascular diseases, urinary diseases, viral infections, thoracic diseases, diseases of the endocrine system.

a given list of phenotypes is, to each disease, i.e. on the ranked relevance of phenotypes. We learn the network given this correlation, that is itself computed using data DDP h (see Section 4.3 on learning large networks). Remark 6.1. The RGG visualised in Figure 2 is a subnet of the full network learnt as the RGG GS,Σ(St) (V , 0.1) where |V | = Ndis , and the inter-observable C

(St)

(rank)

(Spearman rank) correlation matrix of data DDP h is ΣC = [sij ], i < j, i, j ∈ {1, . . . , Ndis }, s.t. this visualised graph has 6052 number of nodes, each with a non-zero degree, and 145210 edges, so that the average node degree is ≈24. This RGG represents our learning of the disease-phenotype network (Figure 2). 6.1. Comparing against earlier work We use Figure 3 to present comparison of our results to HSG’s, including a comparison between the relative number of nodes i.e. diseases, in each of the 19 disease classes that HSG classify their reported network into and our results (shown in the left panel of this figure). Our learnt RGG tallies very well with the earlier result. The right panel of Figure 3 displays the ratio of intra-class to inter-class variance of each disease-class that we identify; value of the area under the Receiver Operating Characteristic curve (ROCAUC) for each cluster identified by HSG is overplotted, where the ROCAUC value for the i-th cluster

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

25

can be interpreted as probability that a randomly chosen node is ranked as more likely to be in the i-th class than in the j-th; i ̸= j; i, j = 1, . . . , 19. Thus, our method of learning a large network allows for the learning of the clustering distribution of the large dataset, for which this network has been learnt. 7. Empirical illustration: identifying optimal cut-off In this application we consider data on the in-channel water level at 893 river cross-sections, within the wider Humber region in northern England. These water levels are modelled response of the river system to a collection of 132 simulated storm events, and the Environment Agency, U.K. holds proprietorship rights on this data. In this dataset, all storms are characterised by five hydraulic boundary conditions (or parameters), namely the upstream fluvial flow from Aire, Don, Ouse, Trent and the downstream water level at the mouth of the Humber Estuary. By construction, the dataset includes 132 blocks - each block corresponding to a simulated storm - with a block comprising 893 rows and six columns. s.t. the first five columns of the i-th block are populated by the values of the storm parameters relevant to the i-th storm, while the last column holds values of the in-channel water level, (referred hereafter as “water level”), at each of the 893 river cross-section locations, (referred hereon as “location”), that are considered; i ∈ {1, . . . , 132}. Then the observations in the first five columns of any block are the same across all 893 rows of this block, though the observed water level values vary from one row to another, within this block. We define the random variable Xi as the spatially-local water level variable within the Humber region, when the i-th storm strikes the region. We denote the j-th observation of variable Xi as xi,j ; this is the water level in the j-th location within the region, due to the i-th storm striking, for j ∈ {1, . . . , 893}. We also define the random variable Wj that represents the water level in the j-th location, triggered by a storm event. Then the observed water level in the j-th location, due to the i-th storm is xi,j , i.e. the i-th observed value of Wj in the available data is xi,j . We learn one RGG on the vertex set V storm := {1, . . . , 132}, in which, the random variable Xi is attached to the i-th node, ∀i ∈ V storm . Thus, the dataset that is used to learn this RGG is Dstorm that comprises the observed values of each of X1 , . . . , X132 . The 132 × 132-dimensional correlation function Σstorm = [|corr(Xi , Xi′ )|] is used to learn this RGG. We refer to this RGG learnt at the cutoff τ as GS,Σstorm (V storm , τ ). We learn a second RGG GS,Σloc (V loc , τ ) on the vertex set V loc := {1, . . . , 893}, in which the variable Wj is attached to the j-th node ∀j ∈ V loc , where the dataset used to learn this RGG is Dloc that holds observations of W1 , . . . , W893 , s.t. the 893 × 893-dimensional correlation matrix Σloc = [|corr(Wj , Wj ′ )|]. We in fact identify the optimal cutoff τ ⋆ in each dataset, by identifying the τ at which the rate of change of the logarithm of the posterior probability of the RGG variable is minimised, s.t. the resulting RGG is most resilient to changes

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

in τ . Thus,

 ⋆

τ := arg min

26

 d log π(GS,Σstorm (V storm , τ )|Dstorm ) . dτ

We find τ ⋆ = 0.07808 in learning given dataset Dstorm . Again, the identification of τ ⋆ in dataset Dloc is yields τ ⋆ = 0.2643.

Fig 4. Block of two left-most panels: panel on the bottom left of this block displays the logarithm of the posterior of the RGG variable of dataset Dloc , learnt at varying values of the cut-off probability. The slope of this graph posterior (with respectto τ ) is plotted against τ on the bottom right. The vertical line in black marks the τ value of 0.2643, at which this slope appears to be minimised. That this τ value is a minima of the slope, is confirmed by plotting the parametric fit to the numerically-computed first derivative of the slope function with respect to τ (in the top left), while the fit to the second derivative of the slope function is plotted on the top right of this block. Numerical derivative computation is done via differencing in Python, and fits to the computed data are identified using NumPy’s polynomial fitting, (in Python). The vertical lines in these panels respectively confirm the approximately zero value of the first derivative of the slope, as well as the positivity of the second derivative - at τ = 0.2643. Block of two right-most panels: the same as the two left panels discussed above, except here results relevant to the RGG variable of daatset Dstorm are presented. The optimal τ value for the RGG of this data is identified as τ ⋆ = 0.07808, at which the slope of the RGG variable with respect to τ attains a minima (shown by vertical line in the bottom right panel of this block); where first derivative of the slope is approximately 0, (shown in the top left); and its second derivative is positive, as depicted by the vertical line at τ ≈ 0.07808 in the top right panel.

Log posterior of the RGG variables computed at τ ∈ {0.001, 0.002, . . . , 0.999}, learnt given the correlation between Xi and Xi′ , ∀i, i′ ∈ {1, . . . , 132}, are displayed in Figure 4. To illustrate the minimisation of the function that we refer to as slope(τ ) := π(GS,Σstorm (V storm , τ )|Σstorm )/dτ , we check the closeness of the (fit to the data on the) derivative of slope(τ ) (with respect to τ ) to 0, and check if the second derivative is positive. Under these checks, τ ⋆ for dataset Dstorm is noted to be about 0.07808. The two right-most panels of Figure 4 show the same for dataset Dloc - for which the corresponding RGG is defined on an 893-sized vertex set - to suggest an optimal τ of about 0.2643. As stated above, in the computation of the derivatives, we use differencing, and fit polynomials (using

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

27

NumPy’s polynomial fitting) to the generated data. The usage of such fitting defines our algorithm of identifying τ ⋆ . Once the optimal values of τ are learnt for Dstorm and that for Dloc , we learn the RGG of each of these datasets, at these identified τ ⋆ . The learnt RGGs are visualised in Figure 5. Each of the 893 locations that we consider in Dloc , is within the channel of one of four rivers (Aire, Don, Ouse and Trent) that flow in the Humber region. Thus, the RGG in which a node is attached to the local water level variable, has four clusters. This is noted in the RGG learnt given data Dloc .

Fig 5. Left: the learnt RGG GS,Σloc (V loc , τ ⋆ ), where τ ⋆ = 0.2643. The RGG displays four clusters, where the nodes of each cluster are represented in distinct colours. Right: GS,Σstorm (V storm , 0.07808) learnt at the optimal τ , given dataset Dstorm . Nodes that belong to the distinct clusters of this graph, are presented in distinct colours. Membership in a cluster is identified using the clustering function in the Network package in Python.

8. Empirical illustration: degree distribution at varying τ We learn RGGs of a dataset comprising binary information, to demonstrate (1) capacity in the RGG-learning to address data on categorical variables, (as well an variables of mixed type); (2) degree distribution, and he effect of the cutoff τ on this. The purpose of this empirical illustration is to showcase such capcity, while our treatment of the available dataset does not offer any realistic interpretation. The data Duser we use, is available in “118379821279745746467.feat”, at: https://snap.stanford.edu/data/ego-Gplus.html (Mcauley and Leskovec, 2014), on the p = 500-dimensional feature vector of n = 699 users of Google+ Circles which was a core facility of the defunct Google+ Social Network, used to categorise “friends’“ circles. The presence of any of the p features for any user, is marked by a “1” in this dataset, while the absence of the same is marked with a “0”. We interpret the p-dimensional vector of 0s and 1s as a vector comprising p observed values of a binary parameter X that characterises a user, where each observation is made at each of p distinct instances. X can only attain values of 0 and 1. Thus, a variable Xi ∈ {0, 1} is attached to the i-th node of the RGG of dataset Duser , ∀i ∈ V user where the RGG is defined on the

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

28

vertex set V user = {1, . . . , n}. Since Xi is categorical, we use the Cramer’s V measure (Sheskin, 2003) to compute corr(Xi , Xj ) ∀i, j ∈ V user . Then we learn the n × n-dimensional correlation matrix Σuser = [corr(Xi , Xj )] and learn the RGG GS,Σuser (V user , τ ), at different τ values, including at the optimal τ for this dataset, which is identified as discussed in the previous section, to be about 0.363. Figure 6 presents four RGGs learnt at τ values of 0.1752, 0.3630, 0.6036, 0.7968; the degree distribution of each of the learnt RGGs is presented below the corresponding RGG. In the RGG learnt at the best cut-off of about 0.363, the two highest values of the degree are the uniqiue two degree values of all points marked respectively in green and blue, in the corresponding RGG. We note this RGG to bear two core clusters, with the remaining points - i.e. noncolourised points bearing low degrees - contributing to the rest of the RGG. The colourised points in the RGG learnt at the optimal τ , are the respective members of the two clusters that distinguish this learnt graph. Thus, almost all members nodes of each cluster bear the same degree, which are also the highest values of the degree amongst all nodes in the graph. At higher values of τ , we can identify the emergence of further structures in the clustering distribution of the learnt RGG, compared to the cleaner, bimodal clustering distribution that we spot for the RGG learnt at the optimal τ . At these high τ values, a narrow interval of degree values is identified, s.t. these degrees correspond to nodes close to the edges of the individual clusters. The loglog frequency distribution of degrees in this narrow interval is akin to a straight line with a steep slope. At lower degrees, the log-log degree distribution is noted to be similar to a power-law at τ values similar to the optimal τ ; however, as τ increases, increasingly more nodes get pulled into the clusters that manifest increasingly more structure, and the degree distribution of the nodes with low degree values then turn Poisson, than power-law. 9. Conclusion In this paper, we have presented a new learning of a random graph of a given multivariate dataset, as an RGG drawn in a probabilistic metric space, to result in an SRGG. We forward the closed-form probability of such an undirected inhomogeneous graph, conditioned on the correlation structure of the dataset for which the graph is learnt. We give the metric of the space that the RGG is drawn in, as the cdf of the disparity variable that we introduce here. Disparity is defined as the absolute difference between the “connectedness” of two nodes, (where said connectedness is given by the mutual edge), and the absolute (partial) correlation between the random variables that are respectively attached to these nodes. Then at the known inter-observable correlation/partial correlation matrix, the edge between these nodes exists in the RGG, if the inter-nodal distance falls short of a chosen cutoff probability τ , i.e. if τ exceeds the cdf of the disparity computed at edge=1. To learn the RGG of a given dataset, we use Rejection Sampling to construct a set of samples drawn randomly from the closed-form probability of the edge

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

29

Fig 6. Top row: RGGs learnt at τ values of 0.1752, 0.3630, 0.6036, 0.7968, from left to right. Nodes that have a degree of 0, are coloured red in these RGGs. The RGG in the second from left panel, represents the RGG learnt at the optimal τ , given the dataset Duser ; nodes in the two modal clusters of this graph are coloured green and blue. Bottom row: figure representing the degree distribution of a learnt RGG, is placed below the panel that presents the learnt RGG.

between a pair of nodes, (conditioned on the correlation between the pair of observables that are attached respectively to the nodal pair). In this set of edge samples, the relative frequency for the edge variable to attain the value of 1, is then proportional to this conditional edge probability. Thus, if the cdf of the disparity computed at edge=1 falls short of a chosen cutoff probability, the edge exists in the RGG, with probability that is (approximately) given by this computed relative frequency. Learning such an RGG in a probabilistic metric space allows for an easy way of illustrating the correlation structure of a dataset, including data comprising information from vector-valued observables, irrespective of whether the observables are numerical or not. We realise that the optimal τ at which the RGG is learnt for a given dataset (Zhang et al., 2025) is also organically realised. The τ value at which the slope of the probability of the RGG is minimised, is identified as the threshold that produces the graph that is the most robust to changes in the threshold. Thus, the optimal τ can be identified because we have advanced the closed-form probability of the RGG, given the inter-observable correlation matrix of the given dataset. In fact, the posterior density of this correlation matrix is closed-form, as we have discussed in Section 3. While the correlation matrix can be updated using its closed-form posterior given the dataset at hand, each edge of the graph can be sampled from the closed-form edge probability, computed using the known inter-observable correlation matrix of the dataset. Inference in such a situation is undertaken using a Metropolis-with-2-block-update scheme. Alternatively, if the correlation structure of the given dataset is known, then we advocate learning the edges by undertaking Rejection Sampling from the closed-form edge probability, (given

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

30

the known correlation). Although omitted from the paper, and included in the Supplement, our approach potentially allows for acknowledgement of measurement errors of the observables, in learning the RGG. It can be demonstrated that the effect of ignoring existent measurement errors, is to distort the learnt RGG, even in the case of a simple, low-dimensional dataset - as demonstrated in Section 1 of the Supplement. A useful fallout of our RGG learning is the formulation on an inter-graph distance function in the space of such RGGs. Subsequent to the learning of the graph variable given the correlation structure of each of two multivariate datasets, we can define an inter-graph distance as a statistical distance or divergence. Indeed, it is possible to compute the Hellinger distance (Banerjee et al., 2015; Matusita, 1953) between the probabilities of the corresponding pair of random graph variables, given the respective dataset. We have indicated the usage of such an inter-graph distance in Equation 5.2, and real-world applications of such a distance function have been undertaken (Zhang et al., 2025, 2024). In our implementation of the correlation between the m-th and m′ -th components of the i-th and j-th observables, we have used the Hellinger distance between the RGG of the data on the m-th component of all observables, and that of the data on the m′ -th component. This inter-graph distance in our work is a generic distance that is informed by the probability distribution of the random graph variable, conditional on the inter-variable correlation matrix of the dataset comprising (noisy) data available on mixed variables in general. The inter-graph distance that we propose, is unlike distances based on counts of node-wise differences (Dekker and Colbert, 2004; Imrich, 2000), or distances defined for certain types of graphs (Bang et al., 2015; Godsil, 1988; Mass, 1987). Appendix A: Proof of Theorem 3.1 (S)

Proof. Likelihood of correlation matrices ΣC and Σ(S) m , given data DS is matrix Normal, (see Section 3). By invoking the matrix Normal density, we write the joint posterior of these correlation matrices using this likelihood, and the priors stated in the statement of the theorem. The joint posterior probability (S) density of ΣC , Σ(S) m , given data DS : h i  h i (S) (S) (S) ΣC , Σ(S) ∝ L Σ(S) ΣC , Σ(S) , i.e. m |DS m , ΣC ; D S m h i 1 (S) ΣC , Σ(S) ∝ p d × m |DS dp 2 (S) 2 Σ(S) (2π) 2 ΣC m  o −d 1 n 2 −1 (S) −1 (S) −1 T exp − tr (Σ(S) Σ , ) (D )(Σ ) (D ) S S m m C 2 (9.1)

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces (S) (S) —using prior on Σ(S) m to be π0 (Σm ) = Σm (S)

α

31

d where α = − − 1, and 2

—using prior on ΣC to be uniform. Then marginalising this joint posterior over Σ(S) m , we get:   (S) π ΣC |DS ∝ Z

1 (S)

ΣC

p 2

× M

Σ(S) m

−d 2 −1

Σ(S) m

d 2



o 1 n (S) −1 (S) −1 T exp − tr (Σm ) DS (ΣC ) (DS ) d(Σ(S) m ) 2 (9.2)

(d×d) Here Σ(S) . Now, m ∈M ⊆R −1 −(d+1) – let Y := (Σ(S) . Then d(Σ(S) dY (Mathai and Pederzoli, m ) m ) = |Y | 1997), h i

−1 – let V −1 := DS (ΣC )−1 (DS )T , =⇒ tr (Σ(S) DS (ΣC )−1 (DS )T m )  −1  tr V Y (using commutativeness of trace), (S)

(S)

≡

so that in Equation 9.2, we get   Z   d d 1  1 (S) |Y | 2 |Y | 2 +1 × exp − tr V −1 Y |Y |−(d+1) dY . π ΣC |DS ∝ p 2 (S) 2 ΣC M (9.3) The integral in the RHS of Equation 9.3 represents the unnormalised Wishart pdf Wd (V , q), over all values of the random matrix Y , where the scale matrix and degrees of freedom of this pdf are V and q = d+1 respectively, i.e. q > d−1. Thus, integral in the RHS of Equation 9.3 is theintegral of the unnormalised  −1  (S) , pdf of Y ∼ Wd (V , q), over the full support of Y ≡ ΣR i.e. the integral in the RHS of Equation 9.3 is the normalisation of this pdf : q qd q 2 2 Γd |V | 2 ≡ 2    −1 d+1 (d+1)(d) d+1 2 (S) 2 2 Γd DS (ΣC )−1 (DS )T , 2  −1 d+1 2 (S) i.e. integral on RHS of Equation 9.3 is proportional to DS (ΣC )−1 (DS )T ,   (S) i.e. π ΣC |DS ∝



1 (S)

ΣC

p 2

DS (ΣC )−1 (DS )T (S)

−1 d+1 2

Now, if DS (ΣC )−1 (DS )T is invertible, −1 ·  −· (S) (S) DS (ΣC )−1 (DS )T = DS (ΣC )−1 (DS )T . (S)

(9.4)

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

32

 −1 (S) (S) – It is given that ΣC is invertible, i.e. ΣC exists. – The original dataset is examined to discard rows that are linear transformations of each other, leading to data matrix DS , no two rows of which are linear transformations of each other =⇒ DS (ΣC )−1 (DS )T is positive definite, i.e. DS (ΣC )−1 (DS )T is invertible,  −1 (d+1)/2 −(d+1)/2 (S) (S) =⇒ DS (ΣC )−1 (DS )T = DS (ΣC )−1 (DS )T . Using this in Equation 9.4: (S)

(S)

  −p/2 −(d+1)/2 (S) (S) (S) π ΣC |DS ∝ ΣC DS (ΣC )−1 (DS )T .

(9.5)

References

Airoldi, E. M. (2007), “Getting Started in Probabilistic Graphical Models,” PLoS Computational Biology, 3(12), e252. Anandkumar, A., Tan, V., Huang, F., and Willsky, A. (2012), “HighDimensional Gaussian Graphical Model Selection: Walk Summability and Local Separation Criterion,” JMLR, 13, 2293–2337. Bandyopadhyay, D., and Canale, A. (2016), “Sparse Multi-Dimensional Graphical Models: A Unified Bayesian Framework,” Journal of Royal Statistical Society Series C, 65(4), 619–640. Banerjee, S., Basu, A., Bhattacharya, S., Bose, S., Chakrabarty, D., and Mukherjee, S. (2015), “Minimum distance estimation of Milky Way model parameters and related inference,” SIAM/ASA Journal on Uncertainty Quantification, 3(1), 91–115. Bang, S., Dubickas, A., Koolen, J. H., and Moulton, V. (2015), “There are only finitely many distance-regular graphs of fixed valency greater than two,” Advances in Mathematics, 269(Supplement C), 1–55. Benner, P., Findeisen, R., Flockerzi, D., Reichl, U., and Sundmacher, K. (2014), Large-Scale Networks in Engineering and Life Sciences, Modeling and Simulation in Science, Engineering and Technology, Switzerland: Springer. Carvalho, C. M., and West, M. (2007), “Dynamic matrix-variate graphical models,” Bayesian Analysis, 2(1), 69–97. URL: https://doi.org/10.1214/07-BA204 Dekker, A. H., and Colbert, B. D. (2004), Network robustness and graph topology„ in Proceedings of the 27th Australasian conference on Computer science, Volume 26, ACSC ’04, Darlinghurst, Australia, Australia: Australian Computer Society, Inc.,, pp. 359–368. Dettmann, C. P. (2018), “Isolation and Connectivity in Random Geometric Graphs with Self-similar Intensity Measures,” Journal of Statistical Physics, 672, 679–700.

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

33

Giles, A. P., Georgiou, O., , and Dettmann, C. P. (2016), “Connectivity of Soft Random Geometric Graphs,” Journal of Statistical Physics, 162(4), 1068– 1083. Godsil, C. D. (1988), “There are only finitely many distance-regular graphs of fixed valency greater than two,” Combinatorica, 8(4), 333–343. Gruber, L., and West, M. (2016), “GPU-Accelerated Bayesian Learning and Forecasting in Simultaneous Graphical Dynamic Linear Models,” Bayesian Analysis, 11(1), 125–149. Hoehndorf, R., Schofield, P. N., and Gkoutos, G. V. (2015), “Analysis of the human diseasome using phenotype similarity between common, genetic, and infectious diseases,” Scientific Reports, 5(10888). URL: http://dx.doi.org/10.1038/srep10888 Imrich, Wilfried; Klavar, S. (2000), Hamming graphs, Product graphs, WileyInterscience, New York: Wiley-Interscience Series in Discrete Mathematics and Optimization. Mass, C. (1987), “Transportation in graphs and the admittance spectrum,” Discrete Applied Mathematics, 16(1), 31–49. Mathai, A. M., and Pederzoli, G. (1997), “Some Properties of Matrix-Variate Laplace Transforms and Matrix-Variate Whittaker Functions,” Linear Algebra and its Applications, 253, 209–226. Matusita, K. (1953), “On the estimation by the minimum distance method,” Annals of the Institute of Statistical Mathematics, 5(1), 59–65. Mcauley, J., and Leskovec, J. (2014), “Discovering social circles in ego networks,” ACM Trans. Knowl. Discov. Data, 8(1). URL: https://doi.org/10.1145/2556612 Menger, K. (1942), “Statistical metrics.,” Proc. Nat. Acad. Sci. USA, 28 (12), 535–537. Ni, Y., Stingo, F. C., and Baladandayuthapani, V. (2017), “Sparse MultiDimensional Graphical Models: A Unified Bayesian Framework,” Journal of the American Statistical Association, 112(518), 779–793. Penrose, M. (2003), Random Geometric Graphs, Oxford: Oxford Studies in Probability, OUP. Penrose, M. D. (2016), “Connectivity of Soft Random Geometric Graphs,” Annals of Applied Probability, 26, 986–1028. Schweizer, B., and Sklar, A. (1983), Probabilistic Metric Spaces, : North-Holland. Sheskin, D. (2003), Handbook of Parametric and Nonparametric Statistical Procedures: Third Edition, Boca Raton: Taylor and Francis. URL: https://books.google.co.uk/books?id=oX5GnwEACAAJ Wang, H., and West, M. (2009), “Bayesian analysis of matrix normal graphical models,” Biometrika, 96, 821–834. Whittaker, J. (2008), Graphical Models in Applied Multivariate Statistics, Switzerland: Wiley. Zhang, C., Dantu, S., Mitra, D., and Chakrabarty, D. (2025), “Identifying critical residues of a protein using meaningfully-thresholded Random Geometric Graphs,”. URL: https://arxiv.org/abs/2506.10015

Chakrabarty et. al/Learning RGGs in Probabilistic Spaces

34

Zhang, C., Grosan, C., and Chakrabarty, D. (2024), “Individualised recovery trajectories of patients with impeded mobility, using distance between probability distributions of learnt graphs,” Artificial Intelligence in Medicine, .

Related documents

• Faster algorithms for approximating combinatorial and geometric data#213219
• Information geometry of absorbing Markov-chain and discriminative random walks#195023
• Probabilistic Methods in Combinatorics#67362
• Parallel Radical and Screen Distributions on Lightlike Hypersurfaces of Indefinite Statistical Manifolds#1011643
• Entropy of Soft Random Geometric Graphs in General Geometries#662394