Measuring Spatial Clustering via Metropolis-Hastings Diffusion Distance Thomas Weighill1 and Chidinma Williams1 1
Department of Mathematics and Statistics, University of North Carolina at Greensboro, Greensboro, NC 27402 July 17, 2026
arXiv:2607.14880v1 [math.ST] 16 Jul 2026
Abstract We propose a novel measure of the discrepancy between two probability distributions f and g on a graph – which we call the diffusion distance – that measures the rate of convergence of f to g under a graph-constrained Markov chain with stationary distribution g. As a default choice for this Markov chain, we use the Metropolis-Hastings transition matrix targeting g with proposals given by a random walk on the graph. Our primary case of interest is when the second distribution g is uniform, in which case the diffusion distance becomes a measure of spatial clustering in f . Used in this way, (MetropolisHastings) diffusion distance to uniformity extends Moran’s I-type measures of spatial autocorrelation by incorporating global graph geometry rather than just local patterns. Indeed, Moran’s I, the most well-known measure of spatial autocorrelation, can be viewed as a one-step heuristic for diffusion distance, so long as specific spatial weights are used. We establish theoretical bounds and a stability result for our measure, connecting it to graph spectra and optimal transport. We then turn our attention to outlining a statistical test for spatial clustering using diffusion distance. Under permutation null models, we derive high-probability bounds on diffusion distance underpinned by exact spectral formulas for convergence of distributions, enabling an efficient statistical test for spatial clustering on large datasets. We empirically compare diffusion distance to Moran’s I both as a numerical measure and as a statistical test. We show that diffusion distance exhibits higher power on synthetic data using a stochastic block model. Empirical analysis of Black population distributions for 100 U.S. cities shows that diffusion distance detects subtle differences in urban segregation patterns that Moran’s I does not.
1
Introduction
A fundamental phenomenon in the analysis of spatial data is that data values can be influenced by their location in relation to other data values. This phenomenon pervades ecology, geography, economics, public health, and urban planning, among other fields. Ignoring this influence and assuming spatial independence can inflate Type I error rates, yield biased regression estimates, and lead to erroneous interpretations of variable relationships [24, 25]. The most common form of spatial dependence studied in the statistics literature is the dependency of data between adjacent (or very nearby) locations. Spatial autocorrelation quantifies this dependency by measuring correlations between neighboring locations [2]. Positive autocorrelation occurs when similar values cluster together, while negative autocorrelation is associated with checkerboard-like arrangements [8, 16]. Quantifying spatial autocorrelation enables more accurate analyses and better-informed policy decisions across diverse applications [2, 34]. Scalar indices like Moran’s I [30] are a common way to quantify or test for spatial autocorrelation. Focusing on the relationship between neighboring data values, however, restricts attention only to local patterns. Large scale patterns, such as the interplay between the global connectivity of the underlying geography and the data, can go undetected unless spatial weights are specifically calibrated to incorporate the correct scale of interactions. In this paper, we address this by introducing a new method to detect spatial dependence in data. Our jumping off point is the observation in [11, Theorem 7.1] that Moran’s I with 1
a particular choice of spatial weights corresponds to how close the data is to uniformity after one step of a Markov chain. This immediately suggests an expansion of Moran’s I – namely, measuring convergence to uniformity over longer runs of a suitable Markov chain. In other words, we model spatial attributes as probability distributions evolving under natural Markov chain dynamics on graphs, and use convergence speed as a measure of spatial clustering. Beginning in the most general context, we introduce the diffusion distance diff p,ϵ (f → g), which measures the minimal steps required for a source distribution f on a graph G to converge to within ℓp tolerance ϵ of a target distribution g under a Markov transition matrix designed to target g as stationary distribution. A natural choice for this transition matrix, which we use throughout, is the Metropolis-Hastings transition matrix whose proposal probabilities are given by a random walk on G. Our main case of interest is when g is uniform (or captures uniformity in some way). In this setting, slow mixing from f reflects strong spatial clustering in f , while fast convergence indicates near-uniformity of f . Notably, diffusion distance depends on the global geometry of the underlying graph, rather than just local correlations. We begin with the theoretical study of diffusion distance and its connections to other important concepts in spatial data analysis. We establish bounds in terms of graph spectra, and in terms of Wasserstein distance on graphs. We prove a stability result controlling the effect of small changes in f on diff p,ϵ (f → g). We also make the connection to Moran’s I mentioned above explicit. In order to use diffusion distance (to uniformity) for statistical tests, we require expectations under various null models similar to those for Moran’s I. Under a random permutation null model, we derive the exact expectation and variance of the n-step reduction in squared deviation, resulting in a formula for bounding the p-value of an observed diffusion distance in terms of the spectrum of the underlying graph. This allows us to avoid sampling a prohibitive number of random permutations to obtain a p-value for large datasets. We empirically demonstrate the power of the diffusion distance test on synthetic graphs generated using a stochastic block model (SBM). Beyond statistical tests, we also demonstrate diffusion distance as a way to quantify spatial clustering patterns. We compute diffusion distances for 100 U.S. cities using Black population distributions as our data of interest; in this context, spatial clustering corresponds to spatial sorting and/or segregation. We compare to Moran’s I as a benchmark, and highlight specific instances where the two measures agree and diverge respectively.
1.1
Related work
Other measures of spatial autocorrelation. The dominant tools for quantifying spatial dependence are Moran’s I, Geary’s C, and the family of local indicators of spatial association (LISA) [31, 14, 3]. These statistics enjoy clear interpretations and closed-form inference under various null models, making them the standard benchmarks in spatial analysis. Moran’s I serves as our benchmark in this paper, and can be interpreted as a one-step heuristic for diffusion distance (Theorem 3.6). Our approach departs from these classical statistics by measuring diffusion under global graph dynamics rather than local correlations. A broad literature enriches Moran’s I-type scalar summaries by capturing spatial heterogeneity through eigenvector spatial filtering, geographically weighted regression, Bayesian spatial models, and network-based methods [8, 9, 38, 12, 15, 17, 7, 37, 13]. Our method adds to this literature by introducing a measure which is readily interpretable in terms of graph dynamics and which has closed-form theoretical results under random permutation models. Spectral graph theory. The study of graphs via the spectra of graph Laplacians (and other operators) is well-established and has been deployed in many data science applications such as spectral clustering, diffusion maps [10], and graph diffusion distance [18, 41] (a measure of similarity between two graphs, and not a distance between distributions on a common graph as defined in this paper). Graph spectra are also intimately related via classical results to mixing times for random walks. We exploit this connection to prove results about convergence from a specific initial distribution f rather than worst-case mixing bounds, so that our measure captures not only information about a graph, but also about how the given distribution f interacts with the global geometry of that graph (see Theorem 3.3).
2
Segregation measurement. Quantifying racial and socioeconomic segregation in cities has motivated a rich methodological literature [27, 35]. Classical indices such as the dissimilarity index and isolation index capture global separation but are insensitive to the spatial arrangement of groups [43]. Spatially explicit segregation measures attempt to incorporate geographic contiguity and local mixing [35, 43, 22]. Our framework contributes to this literature by providing a mixing-time interpretation of segregation: a city is more segregated insofar as its demographic distribution takes longer to equilibrate under local graph diffusion. This reinterpretation connects segregation measurement directly to the spectral geometry of the city’s adjacency graph. Our work is similar in spirit to measures of segregation that attempt to measure the propensity of members of one class to encounter another class [39, 36], but there are qualitative differences. For example, locally heterogenous population distributions with a global tilt result in many encounters, which may hide the overall trend. This might be less of an issue when measuring segregation, but can be a problem for statistical inference. Optimal transport and geometry-aware comparison. Optimal transport provides a geometry-aware way to compare probability distributions by quantifying the cost of redistributing mass [33, 42]. Recent spatial applications include spatio-temporal transport [19], and geometry-aware segregation measures [32]. Our approach is complementary: rather than asking how costly it is to optimally move mass from f to τ , we ask how many steps of graph-constrained diffusion are required.
1.2
Limitations
Our target application is geospatial data, and like almost all geospatial data analysis techniques, our method is vulnerable to the Modifiable Areal Unit Problem (MAUP), wherein reaggregation of data by different geographic units results in different outputs. The fact that our measure includes a diffusion process suggests we are less vulnerable to the “checkerboard problem” (a reaggregation of an checkerboard pattern can suddenly become uniform). Another limitation is that we do not attempt to model actual population dynamics, instead relying on random walks on graphs as an approximation, leaving aside many factors such as public transport availability or traffic patterns. Notably, our measure depends heavily on the geometry of geographic adjacency graphs, which is understudied (but see the very recent [1]).
2
Preliminaries
P We denote by ∆N the probability simplex ∆N = {v ∈ RN | ∀i vi ≥ 0, i vi = 1}. Since all our state spaces are finite, probability distributions are elements of ∆N . Vectors will be row vectors throughout this paper, and we denote by 1 the all-ones vector of the appropriate size. We denote the uniform distribution N1 1 by τ throughout. A matrix M is row-stochastic if M 1⊤ = 1⊤ , and bi-stochastic if in addition 1M = 1.
2.1
Spatial autocorrelation
Moran’s I is a classical global measure of spatial autocorrelation that quantifies whether nearby P observations N tend to have similar values [30, 3]. Let V = {v1 , . . . , vN } denote observations with mean v̄ = N1 i=1 vi , and let W = [wij ] be a nonnegative spatial weight matrix. Then Moran’s I is I = PN
i=1
N PN
j=1 wij
PN PN ·
i=1
j=1 wij (vi − v̄)(vj − v̄) . PN 2 i=1 (vi − v̄)
(1)
Positive values indicate positive spatial autocorrelation, negative values indicate dispersion, and values near zero indicate weak global structure [3]. Under the randomization null and traceless spatial weights matrix, E[I] = − N 1−1 , providing the classical baseline for spatial randomness [9, 29]. A common default choice for the spatial weights, which we use throughout the paper, are row-standardized weights based on an adjacency graph: ( 1 i and j are adjacent and i ̸= j wij = deg i 0 otherwise 3
2.2
Markov chains and Metropolis-Hastings
Let P ∈ RN ×N be a row-stochastic transition matrix on a finite state space with N states. For a distribution f ∈ ∆N modeled as a row vector, the distribution after n steps is f P n . A distribution τ is stationary for P if τ = τ P . When P is irreducible and aperiodic on a finite state space, it admits a unique stationary distribution τ , and f P n → τ as n → ∞ for any initial distribution f . Given a transition matrix Q (referred to as a proposal chain) and a target distribution g, we can define the Metropolis-Hastings Markov chain P whose transition probability from i to j ̸= i is given by Qji gj . Pij = min Qij , gi Practically speaking, this amounts to proposing a new state using Q and accepting with probability proportional to the proposed change in g. By design, g is stationary for P . Metropolis-Hastings Markov chains were originally designed to sample from complicated distributions known only up to a constant factor. We borrow this construction for a different purpose, namely to construct a canonical Markov chain with specified equilibrium distribution. Note that by [4, Theorem 1], we can view the Metropolis-Hastings Markov chain as an optimal approximation to the proposal chain Q with specified stationary distribution.
2.3
Optimal transport and total variation distance
Pn Pm Given probability measures µ = i=1 ai δvi and ν = j=1 bj δuj on a finite metric space (X, d), the pWasserstein distance measures the cost to move µ to ν, and is given by Wp (µ, ν) = min
π∈Π(µ,ν)
n X m X
1/p d(vi , uj )p πij
,
i=1 j=1
where Π(µ, ν) denotes the set of couplings with marginals µ and ν [20, 40]. Recall that the total variation PN distance between µ and ν is , ∥µ − ν∥TV = 12 i=1 |µ(vi ) − ν(vi )|. On spaces of diameter D, these quantities satisfy the following: W1 (µ, ν) ≤ 2D ∥µ − ν∥TV (see [42]).
3
(Metropolis-Hastings) diffusion distance
We begin with the general concept of the diffusion distance from one probability distribution to another on a common state space. Definition 3.1 (Diffusion distance in ℓp -norm). Let f, g ∈ ∆N be probability distributions on a finite state space, let P be a Markov transition matrix, let 1 ≤ p < ∞, and let ϵ > 0. The ℓp -diffusion distance from f to g is n diff P p,ϵ (f → g) := min n ∈ N0 : ∥f P − g∥p < ϵ . Despite the name, diffusion distance is not a metric. It is assymetric in general, as suggested by the notation f → g, and for certain choices of g and P can take the value ∞. The parameter ϵ is a tolerance which must be chosen ahead of time (see Remark 3.2 for some recommendations). For practical purposes we are mostly interested in the following special case (in which case we suppress the P in diff P ): • f is a probability distribution on a graph G and g = τ is uniform, • P is the Metropolis-Hastings Markov chain targeting τ whose proposal chain is a random walk on G. Note that P is symmetric and bistochastic in this special case. We refer to this setup as Metropolis-Hastings diffusion distance to uniformity. Note that when G is connected, P is ergodic as soon as G is non-bipartite (e.g. contains a triangle), and thus we assume throughout that this is the case.
4
Remark 3.2 (Heuristics for choosing p and ϵ). Metropolis-Hastings diffusion distance to uniformity depends on the choice of p and ϵ. While there may be case-specific reasons for a particular choice, we give some general recommendations here. For statistical tests, the underlying graph is typically fixed, and the best theoretical results and approximations are available for p = 2 (Theorem 4.4). A reasonable way to pick ϵ for a graph with N nodes is to note that E[||f − τ ||2 ] ≈ √1N for large N where f is drawn uniformly from the simplex ∆N , and set ϵ to be about a tenth of this value. This ensures that even spatially independent data will typically have positive diffusion distance, creating a strong right-tailed test. For measuring clustering directly, the graph size may vary, and p = 1 has better scale invariance. In particular, note that a distribution f which is zero on a fraction x ∈ [0, 1] of the nodes of a graph and uniform on the rest satisfies ||f − τ ||1 = 2x. We can thus visualize ϵ as the ℓ1 deviation from uniform achieved when the fraction ϵ/2 of the graph is uniformly covered, regardless of the size of the graph, which can help us calibrate ϵ.
3.1
Theoretical bounds
We now record some theoretical properties of diffusion distance. Each is stated in a convenient level of generality, but note that all results in this section apply to Metropolis-Hasting diffusion distance to uniformity. Theorem 3.3 (Spectral upper bound for ergodic chains targeting uniformity). Let P be the transition matrix of an irreducible, aperiodic Markov chain on N states with stationary distribution τ , let 1 ≤ p < ∞, and assume P is diagonalizable over C. Let 1 = λ1 , λ2 , . . . , λN be the eigenvalues of P , with right eigenvectors ⊤ N {ui }N i=1 and left eigenvectors {vi }i=1 satisfying P ui = λi ui ,
vi⊤ P = λi vi⊤ ,
v1⊤ = τ,
u1 = 1,
and the biorthogonality relations vi⊤ uj = δij , where δij = 1 if i = j and 0 otherwise. For i ≥ 2, normalize so that ∥vi⊤ ∥p = 1, and set λ := max2≤i≤N |λi | < 1. Write any initial distribution f with f 1 = 1 as PN f = i=1 ci vi⊤ with ci = f ui . If f = τ then diff P p,ϵ (f → τ ) = 0; otherwise, & diff P p,ϵ (f → τ )
≤
log
PN
i=2 |ci | / ϵ
'
log(1/λ)
.
Proof. Since P is irreducible and aperiodic, it is primitive: some power of P has all strictly positive entries [28]. The Perron–Frobenius theorem then guarantees that the eigenvalue 1 is simple and that |λi | < 1 for every i ≥ 2 [28]; in particular λ < 1. Because {vi⊤ }N i=1 is a basis of left eigenvectors, every row vector f with f 1 = 1 has a unique expansion PN f = i=1 ci vi⊤ . Multiplying on the right by uj and applying biorthogonality gives cj = f uj for each j. In particular c1 = f 1 = 1, so since v1⊤ = τ , N X f −τ = ci vi⊤ . i=2
Using vi⊤ P n = λni vi⊤ (which follows by induction from the left eigenvector relation vi⊤ P = λi vi⊤ ) together with τ P n = τ , N X f P n − τ = (f − τ )P n = ci λni vi⊤ . i=2
By the triangle inequality and ∥vi⊤ ∥p = 1, ∥f P n − τ ∥p ≤
N X
|ci | |λi |n ≤ λn
i=2
N X i=2
where the last step uses |λi | ≤ λ for all i ≥ 2.
5
|ci |,
PN If f = τ , then ci = τ ui = v1⊤ ui = δ1i for all i by biorthogonality, so i=2 |ci | = 0 and f P n = τ for all PN n n n ≥ 0, giving diff P p,ϵ (f → τ ) = 0. Otherwise, set S := i=2 |ci | > 0. The bound ∥f P − τ ∥p ≤ λ S falls below ϵ as soon as λn S < ϵ, which is equivalent to log(S/ϵ) . log(1/λ)
n>
The smallest such integer n is ⌈log(S/ϵ)/ log(1/λ)⌉, which gives the stated upper bound on diff P p,ϵ (f → τ ). Proposition 3.4 (Lower bound in terms of Wasserstein distance). Let G = (V, E) be a graph with shortestpath metric d, diameter D, and unit edge lengths. Let P be a Markov transition matrix with Pij > 0 only if (i, j) ∈ E, with stationary distribution g. Then for a probability measure f on V , W1 (f, g) − Dϵ ≤ diff P 1,ϵ (f → g). Proof. Let
n := diff P 1,ϵ (f → g).
By definition of diffusion distance, Using (2.3),
∥f P n − g∥1 < ϵ.
W1 (f P n , g) ≤ 2D ∥f P n − g∥TV = D ∥f P n − g∥1 < Dϵ.
By the triangle inequality,
W1 (f, g) ≤ W1 (f, f P n ) + W1 (f P n , g).
It therefore remains to bound W1 (f, f P n ). Consider the matrix
Γ := diag(f ) P n .
Its (i, j) entry is
Γij = fi (P n )ij .
Since P n is row-stochastic, the row sums of Γ satisfy X X Γij = fi (P n )ij = fi , j
j
and the column sums satisfy X
Γij =
i
X
fi (P n )ij = (f P n )j .
i
Thus Γ is a coupling of f and f P n , and therefore W1 (f, f P n ) ≤
X
d(i, j) Γij .
i,j
Now, if (P n )ij > 0, then there exists a path from i to j of length at most n using edges along which P assigns positive transition probability. Since Pij > 0 only if (i, j) ∈ E and each edge has unit length, it follows that d(i, j) ≤ n whenever Γij > 0. Hence X X W1 (f, f P n ) ≤ n Γij = n Γij = n. i,j
i,j
Combining the bounds gives W1 (f, g) < n + Dϵ. Equivalently,
W1 (f, g) − Dϵ < n = diff P 1,ϵ (f → g).
Since diff P 1,ϵ (f → g) is an integer, this implies W1 (f, g) − Dϵ ≤ diff P 1,ϵ (f → g), as claimed. 6
We close with a stability result controlling the effect of perturbations of the source distribution. In particular, we prove that the diffusion distance from a function g which is close to f is bounded above and below by the diffusion distance of g with shifted tolerances. Proposition 3.5 (Stability from source distribution to uniform). Let P be a bistochastic Markov transition matrix on N states with uniform stationary distribution τ , and let 1 ≤ p < ∞. Let f, g ∈ ∆N be probability distributions such that ||f − g||p < δ. Then for any ϵ > δ, P P diff P p,ϵ+δ (f → τ ) ≤ diff p,ϵ (g → τ ) ≤ diff p,ϵ−δ (f → τ )
Proof. Only the left inequality requires proof; the right follows by symmetry, after replacing ϵ with ϵ − δ. Let n := diff P p,ϵ (g → τ ). By definition of diffusion distance,
∥gP n − τ ∥p < ϵ.
Using the triangle inequality, ∥f P n − τ ∥p = ∥f P n − gP n + gP n − τ ∥p ≤ ∥(f − g)P n ∥p + ∥gP n − τ ∥p . Since P is bistochastic, the Birkhoff–von Neumann theorem [5] expresses P as a convex combination of permutation matrices: K K X X P = αk Πk , αk ≥ 0, αk = 1, k=1
k=1
where each Πk is a permutation matrix. For any x ∈ R , each permutation preserves the ℓp norm, so N
∥xΠk ∥p = ∥x∥p . Hence ∥xP ∥p =
K X k=1
Applying this inductively gives Therefore,
≤
αk xΠk
K X
αk ∥xΠk ∥p =
K X
αk ∥x∥p = ∥x∥p .
k=1
k=1
∥xP n ∥p ≤ ∥x∥p
for all n ≥ 0.
p
∥(f − g)P n ∥p ≤ ∥f − g∥p < δ.
Combining the bounds,
∥f P n − τ ∥p < δ + ϵ.
By definition of diffusion distance, this implies P diff P p,ϵ+δ (f → τ ) ≤ n = diff p,ϵ (g → τ ),
as claimed.
3.2
Connection to Moran’s I
To study one-step behavior, we use the squared deviation from a reference distribution g, devg (f ) :=
1 ∥f − g∥22 . N
(2)
The one-step squared deviation devg (f P ) provides the link between our framework to Moran-type spatial autocorrelation statistics. 7
(a) I = 0.89, diff 0.01 = 118
(b) I = 0.89, diff 0.01 = 164
(c) I = 0.89, diff 0.01 = 301
Figure 1: Three graphs with uniform distribtion f supported on the yellow nodes. The graphs have equal diameters and equal Moran’s I for f , but different diffusion distances to uniform with p = 2, ϵ = 0.01 due to the change in geometry. Theorem 3.6 (Moran’s I as normalized one-step squared deviation [11]). Let P be a bistochastic Markov matrix on N states with uniform stationary distribution τ , let f ∈ ∆N be a probability vector with devτ (f ) > 0, and set M := P P ⊤ . Then devτ (f P ) Iτ (f, M ) = . devτ (f ) Proof. Since P is bistochastic, M 1 = P P ⊤ 1 = P 1 = 1, so every row of M sums to 1 and N X
S0 =
Mij = N.
i,j=1
Substituting S0 = N into the definition of Iτ (f, M ), Iτ (f, M ) =
(f − τ )⊤ M (f − τ ) . ∥f − τ ∥22
Dividing the definitions of devτ (f P ) and devτ (f ), 1 (f − τ )⊤ M (f − τ ) devτ (f P ) (f − τ )⊤ M (f − τ ) = = N = Iτ (f, M ). 1 2 devτ (f ) ∥f − τ ∥22 N ∥f − τ ∥2
Diffusion distance therefore expands Moran’s I in two ways. Firstly, the usual spatial weights are replaced by the Markov chain-inspired weights P P ⊤ with P a transition matrix, as first proposed in [11]. And secondly, multi-step convergence is measured instead of just one-step convergence. This hints at a limitation of Moran’s I: as a one-step statistic, it cannot detect structural differences that manifest only across multiple steps. Two systems may exhibit identical one-step contraction yet very different long-run convergence behavior, which diffusion distance is designed to capture. Figure 1 shows an illustrative example with grid (sub)graphs.
4
Inference under permutation null models
In addition to quantifying the extent of spatial clustering, Metropolis-Hastings diffusion distance to uniformity can be used as a statistical test for the presence of spatial clustering. To enable this, we develop inference under permutation null models. Given a probability distribution f ∈ ∆N , the permutation null model is the uniform distribution on all permutations of the values of f . If π is a permutation of {1, . . . , N }, then we denote by fπ the vector f with values rearranged by π. Permuting spatial labels disrupts location dependence while preserving the empirical mass values, and thus provides a natural baseline for assessing whether an observed spatial pattern is unusually clustered relative to random relabelings on the same graph. It is also a common null model used in applications of Moran’s I [44], since it does not assume any particular form of the distribution at each location. We analyze one-step and n-step squared deviations, for which null moments can be computed exactly. 8
4.1
Squared deviation under the null
Let π be a uniformly random permutation of {1, . . . , N }. The following theorem collects the first and second moments of the permutation null distribution of dev(fπ P ). Theorem 4.1 (Permutation moments of squared deviation, uniform stationary distribution). Assume N ≥ 4. Let P be an N × N transition matrix with uniform stationary distribution τ = N1 1, and set M := P P ⊤ . Let π be uniformly distributed over all permutations of {1, . . . , N }, let f be a probability vector, and define Za := fa − N1 ,
dev(f ) := N1
N X
Za2 ,
S4 (Z) :=
a=1
N X
Za4 ,
∆ :=
a=1
N X
Mii2 .
i=1
Writing D := (N − 1)(N − 2)(N − 3), the first moment is tr(M ) − 1 , Eπ dev(fπ P ) = dev(f ) N −1 and the second moment is dev(f )2 h 2 (N − 3N + 3) tr(M )2 + 2tr(M 2 ) − N 2 Eπ dev(fπ P )2 = ND − 2(N 2 − 6N + 6)tr(M ) − 3N (N − 1)∆ +
i
S4 (Z) h −(N − 1) tr(M )2 + 2tr(M 2 ) + 2N N 2D i − 4tr(M ) + N (N + 1)∆ ,
and consequently 2 Vπ dev(fπ P ) = Eπ dev(fπ P )2 − Eπ dev(fπ P ) . Under the permutation null, we see that the expected one-step squared deviation factors into the initial deviation dev(f ) and a graph-dependent term controlled by the diagonal of M = P P ⊤ , which encodes the self-similarity of the chain’s one-step transitions. The second-moment formula involves the graph-structural quantities tr(M ), tr(M 2 ), and ∆. The full proof of Theorem 4.1 is outlined in Appendix B. Observation 4.2. Note that when P is symmetric, we can express the key quantities in Theorem 4.1 in terms of the eigenvalues λi and orthonormal eigenvectors vi of P . The following two require only eigenvalues X X tr(M ) = tr(P 2 ) = λ2i , tr(M 2 ) = tr(P 4 ) = λ4i i
i
While the third also requires eigenvectors: !2 ∆=
X i
2 2
Pii
=
X X i
2 λ2k vik
k
=
XX k
ℓ
λ2k λ2ℓ
X
2 2 vik viℓ
i
2 2 2 Note that since the vk are unit vectors, we have i vik viℓ ≤ i vik = 1 for all k, ℓ. The eigenvector dependence disappears when the diagonal of M is constant: if Mii = c for all i, then tr(M ) = N c and hence P ∆ = i Mii2 = N c2 = tr(M )2 /N . This holds whenever the graph is D-regular, where Pij = 1/D on each edge and so Mii = 1/D for every i. In that case, all three quantities depend on the eigenvalues of P alone.
P
P
Remark 4.3. Given that one-step squared deviation is equivalent to Moran’s I by Theorem 3.6, one might assume that we can reuse classical moment calculations for Moran’s I. However, these classical results always assume traceless spatial weights matrix (see e.g. [9]), which is very different to our case where M = P P ⊤ typically has positive trace, requiring us to derive new results. Note that if tr(M ) = 0 and dev(f ) = 1 then Theorem 4.1 recovers the classical expectation of − N 1−1 . 9
4.2
High-probability control of diffusion distance
For the n-step null and the resulting concentration bound, we restrict to normal transition matrices P , for which the spectral structure of Mn = P n (P n )⊤ is fully determined by the eigenvalues of P P ⊤ . In particular, this is the case for a Metropolis-Hastings chain targeting uniformity. The proof of the following Theorem uses Cantelli’s inequality and the expectations in Theorem 4.1. Theorem 4.4 (High-probability bound on diffusion distance under permutation). Assume N ≥ 4. Let P be an N × N normal transition matrix with uniform stationary distribution τ = N1 1, let f ∈ ∆N be a probability distribution, and let π be a uniformly random permutation of {1, . . . , N }. For n ≥ 0, define En := Eπ dev(fπ P n ) , Vn := Vπ dev(fπ P n ) , ρn (ϵ) := ϵ2 − N En where En and Vn are obtained from Theorem 4.1 by replacing P with P n . For n large enough that ρn (ϵ) > 0, we have N 2 Vn . Pπ diff P 2,ϵ (fπ → τ ) > n ≤ 2 N Vn + ρn (ϵ)2 Proof. Since P is bistochastic with uniform stationary distribution τ = N1 1, for every probability vector x, ∥x − τ ∥22 =
N X i=1
Applying this to x = fπ P ,
xi −
N 1 2 1 X 1 2 =N· = N dev(x). xi − N N i=1 N
n
∥fπ P n − τ ∥22 = N dev(fπ P n ).
Therefore n diff P 2,ϵ (fπ → τ ) > n ⇐⇒ ∥fπ P − τ ∥2 ≥ ϵ
⇐⇒ ∥fπ P n − τ ∥22 ≥ ϵ2 ⇐⇒ N dev(fπ P n ) ≥ ϵ2 . Since P is normal, P n is also normal, and P n (P n )⊤ = U diag(µn1 , . . . , µnN ) U ⊤ , where P P ⊤ = U diag(µ1 , . . . , µN ) U ⊤ is an orthogonal diagonalization with µ1 = 1. In particular, P n is bistochastic, since P n 1 = 1 and 1⊤ P n = 1⊤ . Therefore Mn := P n (P n )⊤ satisfies Mn 1 = 1, and Theorem 4.1 applies with P replaced by P n . This gives tr(Mn ) − 1 , En := Eπ dev(fπ P n ) = dev(f ) N −1 Vn := Vπ dev(fπ P n ) = Eπ dev(fπ P n )2 − En2 , where Eπ [dev(fπ P n )2 ] is given explicitly by Theorem 4.1 with M , ∆ replaced by Mn , ∆n . Since ∥fπ P n − τ ∥22 = N dev(fπ P n ), it follows by linearity and scaling that Eπ ∥fπ P n − τ ∥22 = N En , Vπ ∥fπ P n − τ ∥22 = N 2 Vn . Cantelli’s inequality. Let X := ∥fπ P n − τ ∥22 , so that Eπ [X] = N En and Vπ (X) = N 2 Vn . We wish to bound Pπ (X ≥ ϵ2 ). Assuming ρn (ϵ) := ϵ2 − N En > 0, we write Pπ (X ≥ ϵ2 ) = Pπ X − N En ≥ ϵ2 − N En = Pπ X − N En ≥ ρn (ϵ) . Cantelli’s inequality states that for any random variable X with finite variance and any t > 0, P(X − E[X] ≥ t) ≤ 10
V(X) . V(X) + t2
Applying this with t = ρn (ϵ) > 0 and Vπ (X) = N 2 Vn , Pπ X − N En ≥ ρn (ϵ) ≤ Since
N 2 Vn . 2 N Vn + ρn (ϵ)2
2 Pπ diff P 2,ϵ (fπ → τ ) > n = Pπ (X ≥ ϵ ),
the stated bound follows.
4.3
Efficient bounds for permutation p-tests
Bringing together Theorems 4.1 and 4.4, we obtain a strategy for bounding the p-value under a permutation null model when P is a symmetric matrix (for example, the Metropolis-Hastings transition matrix). This strategy avoids having to sample many permutations. Given that the number of permutations grows like N !, this is very useful for large N . Suppose we are given P and an observed probability distribution f . We perform the following steps. 1. Compute n = diff P 2,ϵ (f → τ ), the observed diffusion distance. 2. Obtain eigenvalues λk and eigenvectors vk for P . 3. Mn := P n (P n )⊤ = P 2n since P is symmetric, and P 2n has eigenvalues λ2n k and the same eigenvectors vk . Observation 4.2 applied to P n therefore gives the three structural quantities directly from the spectrum of P : 2 X X XX 2 tr(Mn ) = λ2n tr(Mn2 ) = λ4n ∆n = λ2n k , k , k vik . k
k
i
k
Note that the eigendecomposition of P is computed once for a given graph; each subsequent computation then requires only raising the eigenvalues to a power. Substituting these into Theorem 4.1 yields En := Eπ dev(fπ P n ) and Vn := Vπ dev(fπ P n ) . 4. Provided ϵ2 − N Eπ dev(fπ P n ) > 0, we use Theorem 4.4 to obtain an upper bound for the required p-value Pπ diff P 2,ϵ (fπ → τ ) > n . If the inequality is violated, then the significance is too marginal to be determined by moment methods, and direct sampling is required. For very large N , we can compute only the k largest eigenvalues and their eigenvectors, and treat Step 3 as an approximation. The quality of the approximation depends on the decay of the eigenvalues of P , which one might hope to be fast for geospatial graphs of fine granularity.
5
Numerical experiments
5.1
Measuring power on stochastic block models
We test the power of the diffusion distance permutation test on randomly generated graphs and distributions. We use a stochastic block model (SBM) to generate a graph G with 40 nodes divided evenly into two classes A and B. An edge appears between nodes of the same class with probability 0.5 and between nodes of different classes with probability 0.05. For each node, a random scalar value is drawn uniformly from [0, 1]. Values for nodes in class B are uniformly shifted by a parameter α ≥ 0; see Figure 2b. All values are normalized to sum to 1 before computing diffusion distance. As α increases, the spatial dependence in the data increases and should be detectable by diffusion distance and Moran’s I tests. We run 1, 000 trials with 1, 000 permutations each, and reject the null hypothesis of spatial independence for p-values less than 0.05. We use p = 2 and ϵ = 0.01 throughout for diffusion distance based on the heuristics in Remark 3.2. Figure 2a shows that when α = 0, there is a false positive rate of about 0.05 for both methods as expected. As α increases, diffusion distance outperforms Moran’s I in terms of power, likely because the data is locally noisy but globally biased. As an additional precaution, we rerun the experiment with ϵ = 0.001 and get almost identical results (see the Appendix). 11
(a) Rejection rates with 95% Wilson confidence interval
(b) Sample dataset with α = 0.2
Figure 2: Power analysis of diffusion distance and Moran’s I permutation tests on SBM graphs.
5.2
U.S. city demographics
We study geospatial clustering using a dataset of 100 U.S. cities with boundaries from the CDC’s 500 city project [6], and Census demographic data obtained from NHGIS [26]. For each city, we use the dual graph of the tract map, with vertices representing census tracts and edges connecting geographically adjacent tracts; all computations are restricted to the largest connected component. Figure 3 illustrates the dual graph construction for Chicago. For tract i, let Bi and Ti denote the Black and total populations, respectively, P and set si = Bi /Ti , with si = 0 when Ti = 0. Our initial distribution f is defined by fi = si / i si . For each city, we compute the Metropolis Hastings diffusion distance to uniformity using f . We choose ϵ = 0.25 and p = 1 based on Remark 3.2. We then compare diffusion distance with Moran’s I, computed on the same census-tract graph for each city. Figure 4 shows a clear positive association between the two measures when diffusion distance is plotted on a log-scale (as expected, since it measures a number of iterations of roughly exponential decay), with a rank correlation of 0.728 over 100 cities. Diffusion distance and Moran’s I therefore capture overlapping but nonidentical aspects of spatial structure. To more closely observe the types of urban segregation treated differently by the two methods, we focus on cities with 20% or higher Black population percentage citywide. Among these, the six cities labeled in Figure 4 stand out as being ranked significantly higher by Moran’s I than by diffusion distance. Roughly speaking, this means they are more locally clustered than they are globally clustered. Looking at the city graphs in Figure 5, we see at least one common pattern: non-Black areas near the city center totally or partially surrounded by high Black population areas (this pattern is much less pronounced in Toledo than in the other five, however). This could be the trace left by actions such as “slum clearance” in the 20th century [23] which displaced Black residents from specific areas. The low diffusion distance rank indicates that while these cities are segregated, this segregation is not stable under natural graph dynamics. The boundary between Black and non-Black neighborhoods is not a natural bottleneck and could be an artifact of specific historical events.
12
Figure 3: Chicago tract map and corresponding dual graph, both colored by tract-level Black population share.
Figure 4: Comparing Metropolis-Hastings diffusion distance to uniformity and Moran’s I on Black population data for 100 U.S. cities using scatter plots for the raw values with a log x-axis (left), and for the ranks among the 100 cities (right).
13
(a) Baltimore MD
(d) Memphis TN
(b) Birmingham AL
(e) Rochester NY
(c) Detroit MI
(f) Toledo OH
Figure 5: Dual graphs of six cities with substantial Black populations whose rank by Moran’s I (local clustering) significantly exceeds their rank by diffusion distance (global clustering). See the red points in Figure 4. Color scale shows the Black population share by Census tract.
14
References [1] Sara Anderson, Sarah Cannon, Brooke Feinberg, and Anne Friedman. Census dual graphs: Properties and random graph models. arXiv preprint arXiv:2604.04960, 2026. [2] Luc Anselin. Spatial econometrics: methods and models. Springer Science & Business Media, 1988. [3] Luc Anselin. Local indicators of spatial association—lisa. Geographical analysis, 27(2):93–115, 1995. [4] Louis J Billera and Persi Diaconis. A geometric interpretation of the metropolis-hastings algorithm. Statistical Science, pages 335–339, 2001. [5] Garrett Birkhoff. Tres observaciones sobre el álgebra lineal. Universidad Nacional de Tucumán Revista, Serie A, 5:147–151, 1946. [6] Centers for Disease Control and Prevention. 500 Cities: City Boundaries. https://chronicdata.cdc. gov/500-Cities-Places/500-Cities-City-Boundaries/n44h-hy2j, 2024. Accessed: 28-May-2024. [7] Yongwan Chun and Daniel A Griffith. Spatial statistics: A theory and methods review. International Encyclopedia of Geography: People, the Earth, Environment and Technology, pages 1–14, 2016. [8] Andrew D Cliff and John K Ord. Spatial autocorrelation. Pion London, 1973. [9] Andrew D Cliff and John K Ord. Spatial Processes: Models and Applications. Pion Ltd, 1981. [10] Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5–30, 2006. [11] Moon Duchin, James M. Murphy, and Thomas Weighill. Measuring segregation via analysis on graphs. SIAM Journal on Matrix Analysis and Applications, 44(1):80–105, 2023. [12] A. S. Fotheringham, Chris Brunsdon, and Martin Charlton. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. John Wiley & Sons, Chichester, UK, 2002. [13] A. S. Fotheringham, Wenbai Yang, and Wei Kang. Multiscale geographically weighted regression. Annals of the American Association of Geographers, 107(6):1247–1265, 2017. [14] Robert C Geary. The contiguity ratio and statistical mapping. The Incorporated Statistician, 5(3):115–146, 1954. [15] Alan E Gelfand, Hyon-Jung Kim, CF Sirmans, and Sudipto Banerjee. Spatial modeling with spatially varying coefficient processes. Journal of the American Statistical Association, 98(462):387–396, 2003. [16] Arthur Getis. A critique of spatial autocorrelation. Geographical Analysis, 40(3):290–300, 2008. [17] Daniel A. Griffith. Spatial Autocorrelation and Spatial Filtering: Gaining Understanding Through Theory and Scientific Visualization. Springer, Berlin, Heidelberg, 2003. [18] David K. Hammond, Yaniv Gur, and Chris R. Johnson. Graph diffusion distance: A difference measure for weighted graphs based on the graph Laplacian exponential kernel. In Proceedings of the IEEE Global Conference on Signal and Information Processing (GlobalSIP), pages 419–422, 2013. [19] Hicham Janati, Marco Cuturi, and Alexandre Gramfort. Spatio-temporal alignments: Optimal transport through space and time. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, pages 1695–1704. PMLR, 2020. [20] Leonid V. Kantorovich. Mathematical methods of organizing and planning production. Management Science, 6(4):366–422, 1939. [21] Jakini Auset Kauba and Thomas Weighill. Topological analysis of us city demographics. La Matematica, 3(4):1400–1425, 2024. 15
[22] Eli Knaap, Levi John Wolf, Rey Sergio, Su Han Chan, Taylor Oshan, and Stefanie Lumnitz. Segregation: A Python package for measuring residential segregation. Journal of Open Source Software, 5(47):1951, 2020. [23] Jessica LaVoice. The long-run implications of slum clearance: A neighborhood analysis. Journal of Public Economics, 236:105153, 2024. [24] Pierre Legendre. Spatial autocorrelation: trouble or new paradigm? Ecology, 74(6):1659–1673, 1993. [25] J.P. LeSage. Spatial econometrics. Department of Economics, University of Toledo, 1998. [26] Steven Manson, Jonathan Schroeder, David Van Riper, Tracy Kugler, and Steven Ruggles. Ipums national historical geographic information system: Version 16.0, 2021. [27] Douglas S. Massey and Nancy A. Denton. The dimensions of residential segregation. Social Forces, 67(2):281–315, 1988. [28] Carl D. Meyer. Matrix Analysis and Applied Linear Algebra. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000. [29] Alexander M. Mood, Franklin A. Graybill, and Duane C. Boes. Introduction to the Theory of Statistics. McGraw-Hill, 3rd edition, 1974. [30] Patrick Alfred Pierce Moran. The interpretation of statistical maps. Journal of the Royal Statistical Society. Series B (Methodological), 1948. [31] Patrick Alfred Pierce Moran. Notes on continuous stochastic phenomena. Biometrika, 37(1/2):17–23, 1950. [32] Daniel Peng and James M. Murphy. Optimal transport as a metric for measuring segregation. In 2025 International Conference on Sampling Theory and Applications (SampTA), pages 1–5, 2025. [33] Gabriel Peyré, Marco Cuturi, and Justin Solomon. Gromov-wasserstein averaging of kernel and distance matrices. International Conference on Machine Learning, pages 2664–2672, 2016. [34] Dirk Pfeiffer, Timothy P Robinson, Mark Stevenson, Kim B Stevens, David J Rogers, and Archie CA Clements. Spatial analysis in epidemiology. Oxford University Press, 2008. [35] Sean F. Reardon and Glenn Firebaugh. Measures of multigroup segregation. Sociological Methodology, 32(1):33–67, 2002. [36] Sandro Sousa and Vincenzo Nicosia. Quantifying ethnic segregation in cities through random walks. Nature Communications, 13:5809, 2022. [37] Jonathan B. Thayn. Eigenvector spatial filtering and spatial autoregression. In Shashi Shekhar, Hui Xiong, and Xun Zhou, editors, Encyclopedia of GIS. Springer, Cham, 2017. [38] Michael Tiefelsdorf and Barry Boots. Modelling spatial processes: The identification and analysis of spatial relationships in regression residuals by means of moran’s i. Journal of Geographical Systems, 5(1):25–45, 1998. [39] Wietske Van der Laan Bouma-Doff, Coen Van Eijl, Henk Friesema, et al. Random walk-based segregation measures. Review of Economics and Statistics, 96(3):383–398, 2014. [40] L. N. Vaserstein. Markov processes over denumerable products of spaces, describing large systems of automata. Problemy Peredači Informacii, 5(3):64–72, 1969. [41] Titouan Vayer, Laetitia Chapel, Rémi Flamary, Romain Tavenard, and Nicolas Courty. Optimal transport for structured data with application on graphs. In Proceedings of the 36th International Conference on Machine Learning (ICML), volume 97 of Proceedings of Machine Learning Research, pages 6275–6284, 2019. 16
[42] Cédric Villani. Optimal transport: old and new, volume 338. Springer Science & Business Media, 2008. [43] David W. S. Wong. Implementing spatial segregation measures in GIS. Computers, Environment and Urban Systems, 27(1):53–70, 2003. [44] Tian Zhang, Ge Lin, and Arthur Getis. A decomposition of moran’s i for clustering detection. Computational Statistics & Data Analysis, 51(12):6123–6137, 2007.
A
Computational and data details
The experiments in Section 5 were run on laptops/MacBooks using Python. The largest city dual graph is New York NY with 1848 nodes. Computing the Moran’s I and diffusion distance values in Figure 4 is almost instantaneous since the sparse Metropolis-Hastings transition matrix is easy to compute and exponentiate as needed. The power analysis in Figure 2 requires multiple permutations and multiple trials and so takes longer. It takes just under a second to perform a single diffusion distance permutation test on one SBM graph with 1,000 permutations; this was on a MacBook with a 2.3 GHz 8-Core Intel Core i9 processor and 64 GB of RAM. The entire experiment for Figure 2 therefore takes about 10 minutes per α value. For much larger graphs, we would expect to use the approximations in Section 4. A version of Figure 2 with ϵ = 0.001 is shown in Figure 6, with almost identical results.
Figure 6: Rejection rates for ϵ = 0.001 with the same setup as in Section 5.1. Census data obtained from NHGIS is used in accordance with their licence; see https://www.nhgis.org/ citation-and-use-nhgis-data. We use the preprocessed city graphs published by the authors of [21] and used in their analysis in that paper.
B
Proof of Theorem 4.1
Preliminaries. Since f is a probability vector, N X
Zi =
i=1
PN
N X
i=1 fi = 1, and therefore
fi − N1
= 1 − 1 = 0.
i=1
Since P is bistochastic, M := P P ⊤ is symmetric and satisfies M 1 = P P ⊤ 1 = P 1 = 1, so N X
Mij = 1⊤ M 1 = 1⊤ 1 = N,
i,j=1 N X i,j,k,ℓ=1
Mij Mkℓ =
N X i,j=1
17
Mij
2
= N 2.
Also, since M is symmetric, tr(M 2 ) =
PN
1 2 i,j=1 Mij . Because τ = N 1 is stationary and P 1 = 1,
fπ P − τ = fπ − τ P = Zπ P, and therefore 1 1 ∥Zπ P ∥22 = Zπ P P ⊤ Zπ⊤ N N N 1 1 X = Zπ M Zπ⊤ = Mij Zπ(i) Zπ(j) . N N i,j=1
dev(fπ P ) =
Permutation moments. Because π is uniformly distributed over all permutations of {1, . . . , N }, for any m distinct positions i1 , . . . , im , the vector (π(i1 ), . . . , π(im )) is uniformly distributed over all ordered m-tuples of distinct elements of {1, . . . , N }. We compute the following eight moments. (a) Since π(i) is uniform on {1, . . . , N }, N
Eπ [Zπ(i) ] =
1 X Za = 0. N a=1
(b) Similarly, N
2 Eπ [Zπ(i) ]=
1 X 2 Z = dev(f ). N a=1 a
(c) For i ̸= j, Eπ [Zπ(i) Zπ(j) ] =
X 1 Za Zb . N (N − 1) a̸=b
Now X
Za Zb =
N X
Za
2
−
a=1
a̸=b
N X
Za2 = 0 − N dev(f ) = −N dev(f ),
a=1
so Eπ [Zπ(i) Zπ(j) ] = −
dev(f ) . N −1
(d) N
4 Eπ [Zπ(i) ]=
1 X 4 S4 (Z) Z = . N a=1 a N
(e) For i ̸= j, 3 Eπ [Zπ(i) Zπ(j) ] =
X 1 Za3 Zb . N (N − 1) a̸=b
Using
PN
b=1 Zb = 0, so that
b̸=a Zb = −Za ,
P
X a̸=b
Za3 Zb =
N X a=1
Za3
X
Zb =
N X a=1
b̸=a
and therefore 3 Eπ [Zπ(i) Zπ(j) ] = −
(f ) For i ̸= j, 2 2 Eπ [Zπ(i) Zπ(j) ]=
Za3 (−Za ) = −S4 (Z),
S4 (Z) . N (N − 1)
X 1 Za2 Zb2 . N (N − 1) a̸=b
18
Since
PN
2 a=1 Za = N dev(f ),
X
Za2 Zb2 =
N X
Za2
2
−
a=1
a̸=b
N X
Za4 = N 2 dev(f )2 − S4 (Z),
a=1
so 2 2 Eπ [Zπ(i) Zπ(j) ]=
N 2 dev(f )2 − S4 (Z) . N (N − 1)
(g) For pairwise distinct i, j, k, 2 Eπ [Zπ(i) Zπ(j) Zπ(k) ] =
X 1 Za2 Zb Zc . N (N − 1)(N − 2) a,b,c distinct
For fixed a, since
b̸=a Zb = −Za ,
P
X
Zb Zc =
X
b,c̸=a b̸=c
Zb
b̸=a
2
−
X
Zb2
b̸=a
= Za2 − N dev(f ) − Za2
= 2Za2 − N dev(f ). Therefore X
Za2 Zb Zc =
N X
Za2 2Za2 − N dev(f )
a=1
a,b,c distinct
= 2S4 (Z) − N 2 dev(f )2 , so 2 Eπ [Zπ(i) Zπ(j) Zπ(k) ] =
2S4 (Z) − N 2 dev(f )2 . N (N − 1)(N − 2)
(h) For pairwise distinct i, j, k, ℓ, Eπ [Zπ(i) Zπ(j) Zπ(k) Zπ(ℓ) ] =
X 1 Za Zb Zc Zd . N (N − 1)(N − 2)(N − 3) a,b,c,d distinct
Expanding 0 =
PN
4
and grouping by multiplicity pattern, X X X X 0 = S4 (Z) + 4 Za3 Zb + 6 Za2 Zb2 + 12 Za2 Zb Zc + Za Zb Zc Zd . a=1 Za
a̸=b
a̸=b
a,b,c distinct
a,b,c,d distinct
Substituting the sums from (e), (f), and (g), X Za Zb Zc Zd = −S4 (Z) + 4S4 (Z) − 6 N 2 dev(f )2 − S4 (Z) − 12 2S4 (Z) − N 2 dev(f )2 a,b,c,d distinct
= 3N 2 dev(f )2 − 6S4 (Z), and therefore Eπ [Zπ(i) Zπ(j) Zπ(k) Zπ(ℓ) ] =
3N 2 dev(f )2 − 6S4 (Z) . N (N − 1)(N − 2)(N − 3)
First moment. Since Eπ [Zπ(i) ] = 0 for all i, the standard identity for expected quadratic forms gives Eπ [Zπ M Zπ⊤ ] = tr M Covπ (Zπ ) . 19
From moments (b) and (c),
dev(f ) Covπ (Zπ ) ij = dev(f ) − N −1
which can be written as Covπ (Zπ ) =
i = j, i ̸= j,
N dev(f ) 1 I− J , N −1 N
where J = 11⊤ . Therefore 1 tr M Covπ (Zπ ) N dev(f ) 1 = tr M I − J N −1 N dev(f ) tr(M J) = tr(M ) − . N −1 N
Eπ [dev(fπ P )] =
Since M 1 = 1, we have tr(M J) = 1⊤ M 1 = N , and so tr(M ) − 1 . Eπ dev(fπ P ) = dev(f ) N −1 Second moment. Squaring the expression for dev(fπ P ) and taking expectations, N 2 Eπ dev(fπ P )2 =
N X
Mij Mkℓ Eπ Zπ(i) Zπ(j) Zπ(k) Zπ(ℓ) .
(3)
i,j,k,ℓ=1
The expected value on the right depends only on the multiplicity pattern of (i, j, k, ℓ). We partition [N ]4 into five sets: Γ4 := {(i, j, k, ℓ) : i = j = k = ℓ}, Γ31 := {(i, j, k, ℓ) : exactly three indices are equal}, Γ22 := {(i, j, k, ℓ) : indices form two distinct equal pairs}, Γ211 := {(i, j, k, ℓ) : exactly one equal pair, other two distinct}, Γ1111 := {(i, j, k, ℓ) : all four indices pairwise distinct}, and write Tα :=
4 (i,j,k,ℓ)∈Γα Mij Mkℓ . Since these five sets partition [N ] ,
P
T4 + T31 + T22 + T211 + T1111 = N 2 . Computing T4 . T4 =
N X
2 Maa = ∆.
a=1
Computing T31 . For fixed a ̸= b, the four tuples in Γ31 with repeated value a andP singleton b contribute 4Maa Mab (using symmetry of M ). Summing over all ordered pairs a ̸= b and using b̸=a Mab = 1 − Maa , T31 = 4
N X
Maa
a=1
=4
N X
X
Mab
b̸=a
Maa (1 − Maa )
a=1
= 4 tr(M ) − ∆ . 20
Computing T22 . For each unordered pair {a, b} with a ̸= b, the six tuples in Γ22 with values {a, a, b, b} 2 contribute 2Maa Mbb + 4Mab . Summing over all such pairs, X X 2 T22 = 2 Maa Mbb + 4 Mab . a<b
a<b
Since 2
X
Maa Mbb = tr(M )2 − ∆,
a<b
4
X
2 Mab = 2 tr(M 2 ) − ∆ ,
a<b
we obtain
T22 = tr(M )2 + 2tr(M 2 ) − 3∆.
Computing T211 . Fix pairwise distinct a, b, c with a the repeated value. The six corresponding tuples in Γ211 contribute 2Maa Mbc + 4Mab Mac . Thus T211 = 2
N X a=1
For the first sum, using
b,c Mbc = N ,
P
X
X
Maa
Mbc + 4
Mab Mac .
a=1 b,c̸=a b̸=c
b,c̸=a b̸=c
P
N X X
b Mab = 1,
c Mca = 1, and
P
b̸=a Mbb = tr(M ) − Maa ,
P
Mbc = N − 1 − 1 + Maa − tr(M ) − Maa
b,c̸=a b̸=c
= N − 2 − tr(M ) + 2Maa . Therefore 2
N X
Maa
a=1
X
Mbc = 2
N X
Maa N − 2 − tr(M ) + 2Maa
a=1
b,c̸=a b̸=c
= 2(N − 2)tr(M ) − 2tr(M )2 + 4∆. For the second sum, since X
b̸=a Mab = 1 − Maa ,
P
Mab Mac =
b,c̸=a b̸=c
X
Mab
2
b̸=a
−
X
2 Mab = (1 − Maa )2 −
b̸=a
X
2 Mab .
b̸=a
Summing over a, N X X a=1 b,c̸=a b̸=c
Mab Mac =
N X
N X X 2 2 1 − 2Maa + Maa − Mab
a=1
a=1 b̸=a
= N − 2tr(M ) + ∆ − tr(M 2 ) − ∆
= N − 2tr(M ) + 2∆ − tr(M 2 ). Combining both parts, T211 = 4N + (2N − 12)tr(M ) − 2tr(M )2 − 4tr(M 2 ) + 12∆.
21
Computing T1111 . By the partition identity, T1111 = N 2 − T4 − T31 − T22 − T211 = N 2 − 4N − (2N − 8)tr(M ) + tr(M )2 + 2tr(M 2 ) − 6∆. Combining. Substituting moments (a)–(h) and the Tα formulas into (3), S4 (Z) S4 (Z) N 2 dev(f )2 − S4 (Z) N 2 Eπ dev(fπ P )2 = T4 − T31 + T22 N N (N − 1) N (N − 1) 2S4 (Z) − N 2 dev(f )2 3N 2 dev(f )2 − 6S4 (Z) + T211 + T1111 . N (N − 1)(N − 2) N (N − 1)(N − 2)(N − 3) We collect coefficients separately using common denominator D = (N − 1)(N − 2)(N − 3). Coefficient of dev(f )2 . The dev(f )2 contributions arise from T22 , T211 , and T1111 : (N − 2)(N − 3) T22 − (N − 3) T211 + 3 T1111 = (N 2 − 3N + 3) tr(M )2 + 2tr(M 2 ) − N 2 − 2(N 2 − 6N + 6)tr(M ) − 3N (N − 1)∆, giving a contribution of i N dev(f )2 h 2 (N − 3N + 3) tr(M )2 + 2tr(M 2 ) − N 2 − 2(N 2 − 6N + 6)tr(M ) − 3N (N − 1)∆ . D Coefficient of S4 (Z). The S4 (Z) contributions arise from all five patterns: (N − 1)(N − 2)(N − 3) T4 − (N − 2)(N − 3) T31 − (N − 2)(N − 3) T22 + 2(N − 3) T211 − 6 T1111 = −(N − 1) tr(M )2 + 2tr(M 2 ) + 2N − 4tr(M ) + N (N + 1)∆, giving a contribution of i S4 (Z) h −(N − 1) tr(M )2 + 2tr(M 2 ) + 2N − 4tr(M ) + N (N + 1)∆ . ND Dividing (4) by N 2 yields the stated formula for Eπ [dev(fπ P )2 ].
22
(4)