Expected value of the one-dimensional Earth Mover's Distance

TL;DR

This paper derives a method to compute the expected Earth Mover's Distance between probability distributions on a finite line, linking combinatorics and algebraic geometry, with applications in clustering analysis.

Contribution

It introduces an easy computation method for the EMD generating function on a finite line, connecting it to algebraic geometry and enabling expected value calculations.

Findings

Derived a generating function for EMD values on [n]

Computed the expected EMD in the one-dimensional case

Applied EMD in clustering analysis for a data set

Abstract

From a combinatorial point of view, we consider the Earth Mover's Distance (EMD) associated with a metric measure space. The specific case considered is deceptively simple: Let the finite set [n] = {1,...,n} be regarded as a metric space by restricting the usual Euclidean distance on the real numbers. The EMD is defined on ordered pairs of probability distributions on [n]. We provide an easy method to compute a generating function encoding the values of EMD in its coefficients, which is related to the Segre embedding from projective algebraic geometry. As an application we use the generating function to compute the expected value of EMD in this one-dimensional case. The EMD is then used in clustering analysis for a specific data set.

Equations

$${\mathcal{J}}_{\mu\nu}=\left\{J\in{\mathbb{R}}^{n\times n}:\begin{array}[]{l}J\mbox{ is a non-negative real number $n$ by $n$ matrix such that }\\ \sum_{i=1}^{n}J_{ij}=\mu_{j}\mbox{ for all }j\mbox{ and }\sum_{j=1}^{n}J_{ij}=\nu_{i}\mbox{ for all }i\end{array}\right\}.$$

$$\mathbb{EMD}(\mu ,\nu )=\underset{J\in {\mathcal{J}}_{\mu \nu }}{\inf }\sum _{i,j=1}^n|i-j|J_{ij}.$$

$${\mathcal{M}}_{p,q}=\frac{(p-1){\mathcal{M}}_{p-1,q}+(q-1){\mathcal{M}}_{p,q-1}+|p-q|}{p+q-1}$$

$$\begin{array}[]{c|ccccc}&A&B&C&D&F\\ \hline\cr X&0&19&8&2&1\\ Y&12&2&5&11&0\\ Z&2&20&2&3&3\end{array}$$

$$\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 12&2&5&0&0\\ 0&0&0&8&0\\ 0&0&0&2&0\\ 0&0&0&1&0\\ \end{array}\right]$$

$$\left[\begin{array}[]{ccccc}2&10&0&0&0\\ 0&2&0&0&0\\ 0&5&0&0&0\\ 0&3&2&3&3\\ 0&0&0&0&0\\ \end{array}\right].$$

$$\left[\begin{array}[]{ccccc}0&0&0&0&0\\ 2&17&0&0&0\\ 0&3&2&3&0\\ 0&0&0&0&2\\ 0&0&0&0&1\\ \end{array}\right].$$

$$\begin{array}[]{c|ccccc}&A&B&C&D&F\\ \hline\cr U&13&13&0&0&4\\ V&9&1&13&2&5\\ W&9&7&8&6&0\end{array}$$

$$\begin{array}[]{c|cccccc}\text{EMD}&\text{U}&\text{V}&\text{W}&\text{X}&\text{Y}&\text{Z}\\ \hline\cr\text{U}&0&24&20&24&24&18\\ \text{V}&24&0&12&26&16&22\\ \text{W}&20&12&0&16&10&16\\ \text{X}&24&26&16&0&26&10\\ \text{Y}&24&16&10&26&0&26\\ \text{Z}&18&22&16&10&26&0\\ \end{array}$$

$$⟨X,Y⟩=\text{mbox}Trace(X^{T}Y).$$

$${\omega }_j=\sum _{i=1}^j{e}_i$$

$${\mathbb{R}}_0^n=\left\{v\in {\mathbb{R}}^n\mid ⟨v,{\omega }_n⟩=0\right\}.$$

$$\begin{array}[]{cccc}\pi_{0}:&\mathbb{R}^{n}&\rightarrow&\mathbb{R}^{n}_{0}\\ &v&\mapsto&\pi_{0}(v)\end{array}$$

$${\pi }_0(v)=v-\frac{⟨v,{\omega }_n⟩}{n}{\omega }_n.$$

$$\Pi =\{{\alpha }_1,\cdots ,{\alpha }_{n-1}\}$$

$$\delta_{i,j}=\left\{\begin{array}[]{ll}1&\mbox{if $i=j$}\\ 0&\mbox{otherwise}\end{array}\right.$$

$$\mathcal{E}(v)=|v_1|+|v_1+v_2|+|v_1+v_2+v_3|+\cdots +|v_1+\cdots +v_n|$$

$$⟨v,{\omega }_j⟩=v_1+\cdots +v_j$$

$$\mathbb{EMD}(\mu ,\nu )=\mathcal{E}(\mu -\nu ).$$

$$\mathcal{C}(s,n)=\{(a_1,a_2,\cdots ,a_n)\in {\mathbb{N}}^n:a_1+\cdots +a_n=s\}.$$

$$\left(\genfrac{}{}{0ex}{}{s+n-1}{n-1}\right)\sim \frac{s^{n-1}}{(n-1)!}.$$

$$C=\left[\begin{array}[]{ccccc}0&1&2&\cdots&n-1\\ 1&0&1&\cdots&n-2\\ 2&1&0&\cdots&n-3\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ n-1&n-2&n-3&\cdots&0\end{array}\right].$$

$$\text{mbox}{EMD}_s(\mu ,\nu )=\underset{J\in {\mathcal{J}}_{\mu \nu }}{\min }⟨J,C⟩,$$

$$H_n(z,t):=\sum _{s=0}^{\infty }\left(\sum _{(\mu ,\nu )\in \mathcal{C}(s,n)\times \mathcal{C}(s,n)}{z}^{\text{mbox}{EMD}_s(\mu ,\nu )}\right)t^s,$$

$$H_1(z,t)$$

$$H_2(z,t)$$

$$H_3(z,t)$$

$$H_{p,q}(z,t):=\sum _{s=0}^{\infty }\left(\sum _{(\mu ,\nu )\in \mathcal{C}(s,p)\times \mathcal{C}(s,q)}{z}^{\text{mbox}{EMD}_s(\mu ,\nu )}\right)t^s.$$

$$H_{p,q}(z,t)=\frac{H_{p-1,q}(z,t)+H_{p,q-1}(z,t)-H_{p-1,q-1}(z,t)}{1-z^{|p-q|}t}$$

$$(i,j)⪯(i^{\prime },j^{\prime })⟺i^{\prime }-i\in \mathbb{N}\text{mbox}and{j}^{\prime }-j\in \mathbb{N}$$

$$\text{mbox}support(J):=\{(i,j):J_{ij}>0\}.$$

Full text

1. Introduction

Fix a positive integer $n$. We will denote the finite set of integers $\{1,\cdots,n\}$ by $[n]$. By a probability measure on $[n]$ we mean, as usual, a non-negative real-valued function $f$ on the set $[n]$ such that $f(1)+\cdots+f(n)=1$. By the probability simplex on $[n]$ we mean the set of all probability measures on $[n]$, denoted ${\mathcal{P}}_{n}$. We view ${\mathcal{P}}_{n}$ as embedded in ${\mathbb{R}}^{n}$. Given $\mu,\nu\in{\mathcal{P}}_{n}$ define the set of joint distribution

[TABLE]

For results concerning the geometry of ${\mathcal{J}}_{\mu\nu}$ see [DK14] and [PSU19] where they are referred to as transportation polytopes and discrete copulas respectively.

The Earth Mover’s Distance is defined as

[TABLE]

We remark that the set ${\mathcal{P}}_{n}\times{\mathcal{P}}_{n}$ is a compact subset of ${\mathbb{R}}^{2n}$ and so by continuity the infinum is actually a minimum value. Also, the “EMD” is sometimes referred to by other names, for example, in a two-dimensional setting it is called the image distance. More generally it is called the Wasserstein metric (see [Mgw]).

We recall that the set of all finite distributions, ${\mathcal{P}}_{n}$, embeds as a compact polyhedron on a hyperplane in ${\mathbb{R}}^{n}$ and inherits Lebesgue measure and has finite volume. We normalize this measure so that the total mass of ${\mathcal{P}}_{n}$ is one. We then obtain a probability measure on ${\mathcal{P}}_{n}$, which is uniform. Similarly, ${\mathcal{P}}_{n}\times{\mathcal{P}}_{n}$ may be embedded in ${\mathbb{R}}^{n}\times{\mathbb{R}}^{n}$ and can be given the (uniform) product probability measure.

From its definition, the function ${\mathbb{EMD}}$ is a metric on ${\mathcal{P}}_{n}$. The subject of this paper concerns the expected value of ${\mathbb{EMD}}$ with respect to the uniform probability measure. In this light, we define a function ${\mathcal{M}}$ on ordered pairs of non-negative integers, $(p,q)$, as

[TABLE]

with ${\mathcal{M}}_{p,q}=0$ if either $p$ or $q$ is not positive. Let ${\mathcal{M}}_{n}={\mathcal{M}}_{n,n}$ for any non-negative integer $n$.

We will prove the following theorem in Section 5.

Theorem 1**.**

Fix a positive integer $n$. Let ${\mathcal{P}}_{n}\times{\mathcal{P}}_{n}$ be given the uniform probability measure defined by Lebesgue measure from the embedding into ${\mathbb{R}}^{2n}$. The expected value of ${\mathbb{EMD}}$ on ${\mathcal{P}}_{n}\times{\mathcal{P}}_{n}$ is ${\mathcal{M}}_{n}$.

From a theoretical point of view, this paper concerns the expected value of ${\mathbb{EMD}}$. Additionally, we consider a discrete version of the EMD and compute the mean. In turn, we discuss the relationship to cluster analysis. Then, we end with a comparison of the theoretical results to grade distributions where we have noticed persistent clustering.

The above theorem is obtained as a limit of a discrete version of ${\mathbb{EMD}}$ (denoted by ${\mbox{EMD}_{s}}$, for non-negative integer $s$) which is described using a generating function. The generating function is a deformation of the Hilbert series of the Segré embedding. Some standard tools from algebraic combinatorics show up in a new way in the proofs.

From a practical point of view, we will also consider a “real world” data set with a finite number of joint probability measures derived from letter grade distributions. Specifically, we consider a network of grade distributions from the University of Wisconsin - Milwaukee campus, where two nodes are joined when the ${\mbox{EMD}_{s}}$ between them falls below a pre-specified distance threshold. Determining this threshold so that data features are revealed is a subject of research. The expected value of the EMD in the uniformly random situation helps guide this choice.

The family of networks obtained by varying the distance threshold will be used for metric hierarchical clustering. When the threshold is set so that the network is connected, the spectrum of the corresponding Laplacian matrix (see [FC]) will be computed, as it also relates to clustering. Our use of the EMD in this context should be viewed as an attempt at exploratory data analysis to identify the “communities” in this network rather than rigorous hypothesis testing.

Acknowledgments: The second author would like to thank Anthony Gamst for conversation concerning cluster analysis (and especially the reference [K]), and Ryan (Skip) Garibaldi for pointing out the relationship between the Earth Mover’s Distance and comparing grade distributions.

We would also like to thank the referees for making substantial improvements to the mathematical content and exposition of this article.

2. Non-technical preliminaries

In this section we consider some specific examples of finite distributions. A motivating situation comes from grade distributions, which in the United States are often considered with five outcomes: A, B, C, D, and F. The standard Grade Point Average (GPA) assigns 4.0 to A, 3.0 to B, 2.0 to C, 1.0 to D, and 0.0 to F. The relative distances between these five grades is computed by the absolute value of the difference of point values. That is, a B grade is three units away from an F, while one unit away from an A.

Suppose we are given three distributions in the five grade setting for classes with $s=30$ students. That is,

[TABLE]

To compare distribution X to distribution Y, one notices that if the 12 A grades in Y were moved down to B, 5 C grades moved up to B, 8 D grades moved up to C, and one D grade moved down to F, then the distributions would be identical. The matrix

[TABLE]

encodes the “conversion”. That is, if the rows and columns correspond to the grades (A,B,C,D,F) then the entry in row i and column j records how many grades to move from position i in Y to position j in X. The entries on the diagonal reflect no “Earth” movement, while the entries in the first sub and super diagonals reflect one unit of movement. The row sums return the X distribution, while the column sums return the Y distribution. In total, the value of ${\mbox{EMD}_{s}}$ is 26.

The Y and Z distributions compare as follows: move 10 B’s up one unit to a grade of A, to reflect the fact that Y had 12 grades of A. We move 5 grades down from B to C, and 3 grades from B to D. This latter change is noted as a jump across two positions which will “cost” 2 units in EMD, and since there are 3 grades to move, this makes an overall contribution of 6. Finally, 2 C grades are moved down to D, and the 3 F grades in the Z distribution would be moved to D in the distribution Y.

The joint distribution matrix for Y and Z is

[TABLE]

Interestingly, the total EMD is again 26.

Finally, the X and Z distributions are compared. The joint matrix is

[TABLE]

The total EMD is 10.

The three distributions above have the same GPA of 2.5. We note that this is a feature, which we point out to give indication that the EMD will clearly distinguish between distributions even if the GPA is constant.

To help the reader gain some intuitive feel for the EMD we augment the three distributions by:

[TABLE]

As an exercise, one can compute the 36 pairwise distances between each of the six provided distributions. We computed them with Mathematica as shown below:

[TABLE]

One can sample distributions of five grades with 30 students in many distinct ways. For each sampling method one can ask how the EMD is distributed. The sampling method could be chosen to accurately simulate synthetic data to match previously observed samples from a particular subject at a particular institution. Or, a prior distribution on grade distributions could be assumed, such as a discretization of the multivariate normal distribution.

Upon exploration of observed data one notices clear clustering of the distributions relative to the EMD. Indeed, if some distributions are encountered more frequently than others in a particular model then clustering should be expected. With this fact in mind one is led to question of sampling distributions at flat random. That is, sampling independently with each distribution being equally likely. The theoretical behavior of the uniform model can then be compared to observed data. Clustering in the uniform model can be considered “random”, while additional observed clustering in a specific data set is likely related to a causal feature. Statistics describing clustering should be understood for the uniform distribution as it has maximal entropy.

For any given distribution, one seeks a theoretical understanding of any given descriptive statistic. The present article restricts the focus to the mean of EMD. Other statistics will be considered in future work. Moreover, we focus on the uniform distribution only over the space of finite probability distributions.

Finally, we note that the results of this article imply that the mean discrete EMD on 30 student, five grade distributions is slightly larger than 26. The maximum EMD is 120 reflecting the fact that the distance between all 30 students with A grades is 120 units away from the distribution with all 30 students with F grade. Such large values of EMD are unlikely. The distribution of actual grade data, as we expect, is skewed to the right (i.e. mean larger than median).

3. Technical Preliminaries

In this section we recall some basic notation from combinatorics and linear algebra that are used throughout the paper.

3.1. Notation from linear algebra

Let ${\mathbb{M}}_{n,m}$ be the vector space of real matrices with $n$ rows and $m$ columns. Throughout, we assume that the field of scalars is the real numbers, $\mathbb{R}$. For $i$ and $j$ with $1\leq i\leq n$, $1\leq j\leq m$ we let $e_{i,j}$ denote the $n$ by $m$ matrix with 1 in the $i$-th row, $j$-th column, and 0 elsewhere. A matrix, $M\in{\mathbb{M}}_{n,m}$ is written as $M=(M_{i,j})$ where $M_{i,j}$ is the entry in the $i$-th row and $j$-th column. So, $M=\sum M_{i,j}\,e_{i,j}$. We assume standard notation for the algebra of matrices. For example, the standard inner product of $X,Y\in{\mathbb{M}}_{n,m}$ is

[TABLE]

In the case that $m=1$ we write $e_{i}=e_{i,1}$ for $1\leq i\leq n$. As usual, Let $\mathbb{R}^{n}$ denote the $n$-dimensional real vector space consisting of column vectors of length $n$. The set $\{e_{1},\cdots,e_{n}\}$ is a basis for $\mathbb{R}^{n}$. For our purposes a very useful alternative basis is given by

[TABLE]

where $1\leq j\leq n$. We call the set of $\omega_{j}$ the fundamental basis for $\mathbb{R}^{n}$. The terminology here comes from the the root system of type A in Lie theory (see [GW]).

Let the orthogonal complement, $\omega_{n}^{\perp}$, to $\omega_{n}$ be denoted by

[TABLE]

Column vectors in $\mathbb{R}^{n}_{0}$ have coordinates that sum to zero. We let $\pi_{0}$ denote the orthogonal projection from $\mathbb{R}^{n}$ onto $\mathbb{R}^{n}_{0}$,

[TABLE]

defined by the formula

[TABLE]

Note that the image of $\pi_{0}$ is $\mathbb{R}^{n}_{0}$, and the kernel contains $\omega_{n}$. For $1\leq j\leq n-1$, let $\tilde{\omega}_{j}=\pi_{0}(\omega_{j})$. Observe that $\tilde{\omega}_{1},\cdots,\tilde{\omega}_{n-1}$ span $\mathbb{R}^{n}_{0}$, and by considering dimension, are a basis for $\mathbb{R}^{n}_{0}$. We will call this set the fundamental basis for $\mathbb{R}^{n}_{0}$.

The subspace $\mathbb{R}^{n}_{0}$ has another basis that is of importance to us:

[TABLE]

where $\alpha_{j}=e_{j}-e_{j+1}$. We refer to $\Pi$ as the simple basis for $\mathbb{R}^{n}_{0}$. An essential point is that $\Pi$ is dual to the fundamental basis. That is $\langle\alpha_{i},\tilde{\omega}_{j}\rangle=\delta_{i,j}$ where

[TABLE]

Let ${\mathcal{E}}:\mathbb{R}^{n}\rightarrow{\mathbb{R}}$ be defined as

[TABLE]

for $v=\sum v_{j}e_{j}\in\mathbb{R}^{n}$.

Observe that since $\mathbb{R}^{n}_{0}\subset\mathbb{R}^{n}$, ${\mathcal{E}}$ is defined on $\mathbb{R}^{n}_{0}$ by restriction.

Let $v=\sum_{i=1}^{n-1}c_{i}\alpha_{i}\in\mathbb{R}^{n}_{0}$. Then, ${\mathcal{E}}(v)=\sum_{i=1}^{n-1}|c_{i}|$. This fact is easily seen since the fundamental basis is dual to the simple basis relative to the standard inner product, and

[TABLE]

for $1\leq j\leq n$.

In Section 5 we will prove that

Theorem 2**.**

For all $\mu,\nu\in{\mathcal{P}}_{n}$,

[TABLE]

This allows for a much more explicit combinatorial analysis of EMD. For situations in which the metric space is not a subset of the real line, the analysis is more difficult. Indeed, ${\mathbb{EMD}}$ is often computed as an optimization problem that minimizes the cost under the constraints imposed by the marginal distribution. Consequently, the computational complexity is the same as for linear programming.

3.2. Compositions and related combinatorics

Let ${\mathbb{N}}$ be the set of non-negative integers. Given $s\in{\mathbb{N}}$, and a positive integer $n$, define:

[TABLE]

An element of the set ${\mathcal{C}}(s,n)$ will be referred to as a composition of $s$ into $n$ parts. (We note that in some places of the literature these are refered to as weak compositions, since we allow zero. However, the distinction is not needed for us.)

It is an elementary fact that there are $\binom{s+n-1}{n-1}$ compositions, and therefore for fixed $n$, the number of compositions of $s$ grows as a polynomial function of $s$ with degree $n-1$. An essential fact for this paper is the asymptotic approximation

[TABLE]

As in the introduction, given compositions $\mu$ and $\nu$ of $s$ we let ${\mathcal{J}}_{\mu\nu}$ denote the set of $n$ by $n$ matrices with row sums $\mu$ and column sums $\nu$. There is a slight difference here in that we are not requiring $\mu$ and $\nu$ to be normalized to sum to one. In the same light, EMD can be extended as a metric on ${\mathcal{C}}(s,n)$.

We fix an $n$ by $n$ matrix $C$ with $i$-th row and $j$-th column entry to be $|i-j|$. That is,

[TABLE]

So, for $\mu,\nu\in{\mathcal{C}}(s,n)$ and regarding the set ${\mathcal{J}}_{\mu\nu}$ as non-negative integer matrices with prescribed row and column sums, we arrive at

[TABLE]

which is a discrete version of EMD. When we take $s\rightarrow\infty$ we recover the value referred to in the introduction.

The function EMD may be further generalized to the case where $C$ has $p$ rows and $q$ columns, with $i,j$ entry $|i-j|$. In this case $\mu$ has $p$ components and $\nu$ has $q$ components (each a composition of $s$). This generalization will be needed in an induction argument in Section 5. However, applications need only consider the $p=q$ case.

A further generalization beyond the scope of this paper is to consider more general cost matrices than $C$. This is equivalent to a variation of the metric.

4. Generating Functions

In algebraic combinatorics it is often useful to record discrete data in a formal (multivariate) power series – sometimes called a generating function. By “formal” we mean that the variables are indeterminates rather than numbers. In fact, from this point of view one can consider formal power series that only converge at zero, yet encode combinatorial data in their coefficients. Consequently, convergence is not an issue. Nonetheless, our series are all geometric series expansions of rational functions, and so will be convergent if, say, all complex variables have modulus less than 1.

Starting from the viewpoint of algebraic combinatorics we define

[TABLE]

where $t$ and $z$ are indeterminates. We see that the coefficient of $t^{s}$ in $H_{n}(z,t)$ is a polynomial in $z$ whose coefficients record the distribution of the values of ${\mbox{EMD}_{s}}$.

It is useful to see the first few values of $H$, which we compute using Mathematica and Theorem 3:

[TABLE]

As before, when considering $p$ by $q$ matrices we can analogously define $H_{p,q}(z,t)$. We also extend the definition so that $H_{p,q}=0$ if either of $p$ or $q$ is not positive. We obtain a similar series, namely

[TABLE]

Note that if $p<q$ we can regard $p$-tuples as $q$-tuples by appending zeros, so the function ${\mbox{EMD}_{s}}$ is defined.

Theorem 3**.**

For positive integers $p$ and $q$,

[TABLE]

if $(p,q)\neq(1,1)$ and $H_{1,1}=\frac{1}{1-t}$.

This proof will also be given in Section 5, after we have developed some of the consequences in the remainder of this section.

4.1. The partially ordered set $[p]\times[q]$

Recall that a partially ordered set is a set $S$ together with a relation, $\preceq$, which is required to be reflexive, antisymmetric and transitive. In particular, given positive integers $p$ and $q$, we define $S=[p]\times[q]$, and

[TABLE]

for $1\leq i,i^{\prime}\leq p$ and $1\leq j,j^{\prime}\leq q$, which is a partially ordered set.

It is important to note that not all elements are comparable with respect to this order. For example, if $p=q=2$, clearly $(1,2)\not\preceq(2,1)$ and $(2,1)\not\preceq(1,2)$. We say that $(1,2)$ and $(2,1)$ are incomparable. A subset of $S$ in which all pairs are comparable is called a chain.

Given a $p$ by $q$ matrix, $J$, we define the support as

[TABLE]

Proposition 4**.**

Let $p$ and $q$ be positive integers, and $s\in{\mathbb{N}}$. For $\mu\in{\mathcal{C}}(s,p)$, $\nu\in{\mathcal{C}}(s,q)$ and $J\in{\mathcal{J}}_{\mu\nu}$, there exists a $J^{\prime}\in{\mathcal{J}}_{\mu\nu}$ such that the support of $J^{\prime}$ is a chain in $[p]\times[q]$, and $\langle J^{\prime},C\rangle\leq\langle J,C\rangle$.

Proof.

Suppose there exist incomparable elements $(i,j),(i^{\prime},j^{\prime})$ such that $J_{ij},J_{i^{\prime}j^{\prime}}>0$. Without loss of generality assume $i^{\prime}<i$, $j<j^{\prime}$, and $0<J_{ij}\leq J_{i^{\prime}j^{\prime}}$. We construct $J^{\prime}$ as follows:

[TABLE]

For $(k,l)\not\in\{(i,j),(i^{\prime},j^{\prime}),(i^{\prime},j),(i,j^{\prime})\}$ let $J^{\prime}_{kl}=J_{kl}$. Next, let $J^{\prime}_{ij}=0$ and

[TABLE]

which is non-negative. And,

[TABLE]

[TABLE]

The result follows. ∎

The point here is that we will only need to consider matrices $J$ with support on a chain.

4.2. A special case of the Robinson-Schensted-Knuth correspondence

In this subsection we recall, in detail, a special case of the Robinson-Schensted-Knuth correspondence (RSK), see [F]. For those familiar with RSK, we consider the case where the Young diagrams have only one row. For non-experts, the exposition here does not require any knowledge of RSK.

Given $\mu\in{\mathcal{C}}(s,n)$ with $\mu=(\mu_{1},\ldots,\mu_{n})$ we define the word of $\mu$ to be a finite weakly increasing sequence of positive integers, $w=w(\mu)=w_{1}w_{2}w_{3}\cdots$ where the number of times $k$ occurs in $w$ is equal to $\mu_{k}$.

As an example, if $\mu=(3,0,2,1,0)$ then

[TABLE]

Note that the length of the word is equal to $s=\sum\mu_{i}$, and all components of $w$ are at most $n$.

Next, for $s\in{\mathbb{N}}$ and positive integers $p$ and $q$, we define

[TABLE]

Proposition 5**.**

For a given $s\in{\mathbb{N}}$ and positive integers $p$ and $q$, we have a bijection, $\Phi$, between ${\mathcal{C}}(s,p)\times{\mathcal{C}}(s,q)$ and $R(p,q;s)$.

Proof.

Given $(\mu,\nu)\in{\mathcal{C}}(s,p)\times{\mathcal{C}}(s,q)$, let $u,v$ be the words of $\mu$ and $\nu$ respectively,

[TABLE]

with $1\leq u_{i}\leq p$, and

[TABLE]

with $1\leq v_{j}\leq q$. Define a $p$ by $q$ matrix by

[TABLE]

for $1\leq i\leq p$ and $1\leq j\leq q$. Note that the support of $J$ is a chain. We define $\Phi(\mu,\nu)=J$. Given $J$, we can recover $\mu$ and $\nu$ as the row and column sums of $J$. ∎

4.3. Rank one matrices and the Segré embedding

We let ${\mathbb{D}}^{\leq k}(p,q)$ denote the set of $p$ by $q$ matrices with rank at most $k$, which is a closed affine algebraic set, called a determinantal variety. For a relatively recent expository article about the role these varieties play in algebraic geometry and representation theory see [EHP].

In this section we consider the $k=1$ case, in our context. Define $P:{\mathbb{R}}^{p}\times{\mathbb{R}}^{q}\rightarrow{\mathbb{M}}_{p,q}$ by

[TABLE]

for $v\in{\mathbb{R}}^{p}$ and $w\in{\mathbb{R}}^{q}$. Note that if $P(v,w)\neq 0$ then the rank is 1. In fact, the image of $P$ consists of those matrices with rank at most 1. So if $p,q>1$ then $P$ is not surjective.

Injectivity of $P$ fails as well since for non-zero $c\in{\mathbb{R}}$, $v$, $w$ we have $P(v,w)=P(cv,\frac{1}{c}w)$. However, if we pass to projective space we recover an injective map.

In this light, let ${\mathbb{R}}{\mathbb{P}}^{n}$ be an $n$-dimensional real projective space, that is:

[TABLE]

where ${\mathbb{R}}v$ denotes the 1-dimensional subspace spanned by non-zero $v$. We will also write ${\mathbb{R}}{\mathbb{P}}^{n}:={\mathbb{P}}({\mathbb{R}}^{n+1})$.

The Segré embedding,

[TABLE]

is defined as follows: first, we note that we can identify ${\mathbb{R}}^{pq}$ with the $p$ by $q$ matrices by choosing bases. Next, given an ordered pair of projective points (i.e. one-dimensional subspaces) we can choose non-zero vectors $v$ and $w$ respectively. The value of the Segré embedding is the one-dimensional subspace in ${\mathbb{M}}_{p,q}$ spanned by the matrix $P(v,w)$. It is easily checked that this map is well-defined and injective.

The image of the Segré embedding gives rise to a projective variety structure on the set-cartesian product of the two projective varieties. The projective coordinate algebra of the Segré embedding is intimately related to $H_{p,q}(z,t)$, which we will see next.

Let $m_{ij}$ be a choice of (algebraically independent) indeterminates. We consider the polynomial algebra

[TABLE]

Then for $1\leq i<i^{\prime}\leq p$, and $1\leq j<j^{\prime}\leq q$ define

[TABLE]

The ideal, ${\mathcal{I}}$, generated by the $\Delta(i,i^{\prime};j,j^{\prime})$ vanishes exactly on the matrices of rank 1. Conversely, any polynomial function that vanishes on the rank at most 1 matrices is in ${\mathcal{I}}$. The algebra of coordinate functions on the rank at most 1 matrices is then isomorphic to the quotient of ${\mathcal{A}}_{p,q}$ by ${\mathcal{I}}$. Define ${\mathcal{R}}(p,q):={\mathcal{A}}_{p,q}/{\mathcal{I}}$.

The point here is that monomials involving variables which have indices that are not comparable with respect to $\preceq$ may be replaced (modulo ${\mathcal{I}}$) with comparable indices. That is, $m_{ij}m_{i^{\prime}j^{\prime}}$ can be replaced with $m_{i^{\prime}j}m_{ij^{\prime}}$. This process may be thought of as “straightening” and is related to the non-negative integer matrices $J$ with support in a chain. The matrix $J$ may be thought of as the exponents in a monomial.

More generally, the situation may be put into the context of Gröbner bases. The cost matrix $C$ used here assigns a number to each pair of indices. This number can be used to scale the degree of $m_{ij}$. Using this new notion of degree, we can set up a partial order of the monomials, which can then be extended (say, lexicographically) to a well ordering of the monomials that is compatible with multiplication. That is, we can create a term order (see [CLO]). The minors generating the ideal ${\mathcal{I}}$ are indeed a Gröbner basis. The complement of the ideal of leading terms is then a vector space basis for the quotient by ${\mathcal{I}}$.

For $s\in{\mathbb{N}}$, let ${\mathcal{A}}_{p,q}^{s}$ denote the subspace of homogeneous degree $s$ polynomials, and set ${\mathcal{R}}_{p,q}^{s}={\mathcal{A}}_{p,q}^{s}/({\mathcal{A}}_{p,q}^{s}\cap{\mathcal{I}})$. Since ${\mathcal{I}}$ is generated by homogeneous polynomials, we have

[TABLE]

That is, we have an algebra gradation by polynomial degree.

The polynomial functions on ${\mathbb{R}}^{p}$ (resp. ${\mathbb{R}}^{q}$) will be denoted ${\mathcal{A}}_{p}$ (resp. ${\mathcal{A}}_{q}$). Given vectors $v\in{\mathbb{R}}^{p}$ and $w\in{\mathbb{R}}^{q}$ an element of the tensor product ${\mathcal{A}}_{p}\otimes{\mathcal{A}}_{q}$ defines a function on ${\mathbb{R}}^{p}\times{\mathbb{R}}^{q}$ with value $f(v)g(w)$. Given an element $(v,w)\in{\mathbb{R}}^{p}\times{\mathbb{R}}^{q}$, and $f\otimes g\in{\mathcal{A}}_{p}\otimes{\mathcal{A}}_{q}$, the value of $f\otimes g$ on $(v,w)$ is given by $f(v)g(w)$. Extending by linearity we obtain an algebra isomorphism from the polynomials on ${\mathbb{R}}^{p}\times{\mathbb{R}}^{q}$ to the tensor product algebra ${\mathcal{A}}_{p}\otimes{\mathcal{A}}_{q}$.

The quadratic map $P$, defined above, gives rise to an algebra homomorphism,

[TABLE]

defined such that $P^{*}(F)$ is a function on ${\mathbb{R}}^{p}\times{\mathbb{R}}^{q}$ from a function $F$ on ${\mathbb{D}}^{\leq 1}_{p,q}$. This is done via the usual adjoint map given by $[P^{*}(F)](v,w)=F(P(v,w))$.

The image of $P^{*}$ is given as the “diagonal” subalgebra:

[TABLE]

From our point of view, the significance of this structure is as follows:

• •

The dimension

[TABLE]

which is equal to the cardinality of ${\mathcal{C}}(s,p)\times{\mathcal{C}}(s,q)$. That is, a basis may be parameterized by a pair of compositions.

• •

The (finite dimensional) vector space ${\mathcal{R}}_{p,q}^{s}$ will have a dimension also equal to the above since $P^{*}$ is an isomorphism.

• •

The bijection $\Phi$ from Proposition 5 establishes that the above dimension is given by the cardinality of $R(p,q;s)$.

• •

A basis for ${\mathcal{R}}_{p,q}^{s}$ may be given by the monomials of the form $\prod m_{ij}^{P_{ij}}$ where $P\in R(p,q;s)$.

• •

These monomials correspond to the monomials in ${\mathcal{A}}_{p}^{s}\otimes{\mathcal{A}}_{q}^{s}$ with exponents $(\mu,\nu)\in{\mathcal{C}}(s,p)\times{\mathcal{C}}(s,q)$.

4.4. The specialization and a derivative

From the definition it is relatively easy to see that

[TABLE]

We next turn to the specialization $H_{p,q}(1,t)$, which turns out to be the Hilbert series of the rank at most 1 matrices. That is to say

[TABLE]

In [EW] this series was computed in Equation (6.4) as

[TABLE]

In this sense, $H(z,t)$ interpolates between the generating function for compositions and the Hilbert series of the determinantal varieties (at least in the rank one case). Moreover, for generic $z$ we have a relationship to the EMD.

Our goal is to compute the expected value of ${\mbox{EMD}_{s}}$. Therefore, it is natural to compute the partial derivative of $H_{p,q}(z,t)$ with respect to $z$, and then set $z=1$. From the definition of $H_{p,q}$,

[TABLE]

To expand this we start with

[TABLE]

the recursive relationship from Theorem 3. Then let the partial derivative of $H_{p,q}$ with respect to $z$ be denoted $H^{\prime}_{p,q}$. We find the derivative using the “quotient rule”

[TABLE]

[TABLE]

When $z=1$ this becomes

[TABLE]

Before proceeding it is useful to see some initial values:

[TABLE]

Both $H_{p,q}(1,t)$ and $H^{\prime}_{p,q}(1,t)$ are rational functions. We anticipate that their numerators are

[TABLE]

and

[TABLE]

Thus, multiplying by $(1-t)^{p+q}$ on both sides of Equation 4.2 gives:

[TABLE]

[TABLE]

or

[TABLE]

If we also note that

[TABLE]

we ultimately obtain

[TABLE]

An easy induction shows that both $W_{p,q}(t)$ and $N_{p,q}(t)$ are polynomials in $t$.

Before proceeding it is instructive to recall our goal of finding the expected value of ${\mbox{EMD}_{s}}$. In this light, define

[TABLE]

for $s\in{\mathbb{N}}$ and positive integers $p$ and $q$. The expected value of ${\mbox{EMD}_{s}}$ will be then obtained from

[TABLE]

which we will show in the proof of Theorem 1 to be ${\mathcal{M}}_{p,q}$. In the next subsection we will use Proposition 6 to find the asymptotic value as $s\rightarrow\infty$ for fixed values of $p$ and $q$. First we need more information about ${\mathcal{N}}(p,q;s)$.

Proposition 6**.**

Given positive integers $p$ and $q$,

[TABLE]

Proof.

We have seen that

[TABLE]

Differentiating the definition $H_{p,q}(z,t)$ term by term and then setting $z=1$ gives the result. ∎

Using Equation 4.3 and the formula for $W_{p,q}(t)$ one can efficiently compute $N_{p,q}(t)$ for specific values of $p$ and $q$. Following the method for generating functions, we multiply $N_{p,q}(t)$ and the series expansion of $\frac{1}{(1-t)^{p+q}}$. Because the series expansion involves only binomial coefficients we are led to an efficient method for finding the expected value of ${\mbox{EMD}_{s}}$ on ${\mathcal{C}}(s,p)\times{\mathcal{C}}(s,q)$ for any given values of $p$, $q$ and $s$. Consequently we determine the values of ${\mathcal{N}}(p,q;s)$.

Some initial data for $N_{n,n}(t)$ for $n=1,\cdots,12$ are:

[TABLE]

The coefficients of the above are non-negative integers, which we prove inductively. Furthermore, they apparently are palindromic – that is the coefficient of $t^{i}$ matches the coefficient of $t^{d-i}$ where $d$ is the polynomial degree. Lastly we note that the coefficients rise in value until the middle and then decrease – that is to say they are unimodal. Polynomials with these properties are often of interest. We conjecture:

Conjecture 1**.**

The coefficients of the polynomials $N_{n,n}(t)$ are unimodal and palindromic.

In Theorem 1, the ${\mathbb{EMD}}$ has been normalized so that $\mu$ and $\nu$ are probability distributions. In terms of pairs of compositions of $s$, this amounts to multiplying by $\frac{1}{s}$. Note that since the order of the pole at $t=1$ in $H_{p,q}(1,t)$ is one less than the order of the pole at $t=1$ in $H^{\prime}(1,t)$ we see that the (normalized) EMD is approaching a constant as $s\rightarrow\infty$. Alternatively, we could choose not to normalize and then obtain a linear growth of $s\,{\mathcal{M}}_{p,q}$.

The following is a table with approximate values of ${\mathcal{M}}_{p,q}$ for $1\leq p,q,\leq 5$,

[TABLE]

We display a 3D plot for $1\leq p,q\leq 12$ in Figure 1.

5. Further Calculations and Proofs of the theorems

This section includes the main technical points of the paper, including the proofs of the main theorems. It is more convenient to prove Theorem 1 after Theorems 2 and 3.

First, however, we make explicit two useful scaled versions of the EMD.

5.1. Unit normalized Earth Mover’s Distance

Averaging $\frac{{\mathcal{E}}(\mu-\nu)}{s}$ over $(\mu,\nu)\in{\mathcal{C}}(s,n)\times{\mathcal{C}}(s,n)$ gives rise to the expected value of the discrete EMD. Taking the limit as $s\rightarrow\infty$ gives the expected value of the normalized EMD on ${\mathcal{P}}_{n}$. Observe that the maximum value of the normalized EMD on ${\mathcal{P}}_{n}$ is $n-1$. The unit normalized EMD will be defined as the normalized EMD scaled by $\frac{1}{n-1}$. This scaling makes ${\mathcal{P}}_{n}$ into a metric space with diameter 1. When working with real data in the last section of this paper, we will always use the unit normalized Earth Mover’s Distance.

As an example, we consider the discrete case where ${\mathcal{P}}_{n}$ is replaced by ${\mathcal{C}}(s,n)$. In particular, choose $s=30$ and $n=5$ and calculate the exact histogram for the unit normalized distance. The mean of the distribution is obtained by expanding

[TABLE]

as a series around $t=0$, then taking the coefficient of $t^{30}$ and dividing by $\binom{30+5-1}{5-1}^{2}$ (the number of ordered pairs of distributions). The approximate value is 26.2938.

For the unit normalized distance we divide by $s(n-1)=30(5-1)=120$. That is, we divide by $s$ to obtain a probability distribution, and then divide by $n-1$ to scale the diameter to 1. The unit normalized mean is approximately 0.219115 as shown in Figure 2.

In the limiting case as $s\rightarrow\infty$, the mean decreases slightly from the $s=30$ case. Again, the scaling sets the diameter of the metric space ${\mathcal{P}}_{n}$ to 1. Thus, the expected value of the unit normalized EMD is

[TABLE]

Recall the notation that ${\mathcal{M}}_{n}={\mathcal{M}}_{n,n}$, for positive integer $n$.

It should be noted that all values of $\widetilde{\mathcal{M}}_{n}$ are rational numbers. We present the following approximate values of $\widetilde{\mathcal{M}}_{n}$ for $n=2,\cdots,12$,

[TABLE]

which are the limiting values as $s\rightarrow\infty$. For finite choices of $s$ and $n$ we can compute the exact mean of the unit normalized EMD by expanding

[TABLE]

in the case when $p=q=n$, and then dividing by $s(n-1)\binom{s+n-1}{n-1}^{2}$. We show the approximate values in the table below.

[TABLE]

5.2. Proof of Theorem 2

Proof.

By continuity it will suffice to show that this is true on a dense subset of ${\mathcal{P}}_{n}$. Specifically, we will consider the special case that if $\mu$ (resp. $\nu$) is of the form

[TABLE]

for some positive integer $s$ and $(a_{1},\cdots,a_{n})\in{\mathcal{C}}(s,n)$. As $s\rightarrow\infty$, such points are dense in ${\mathcal{P}}_{n}$

For $p\leq n$ (resp. $q\leq n$) we can regard ${\mathcal{C}}(s,p)$ (resp. ${\mathcal{C}}(s,q)$) as being embedded in ${\mathcal{C}}(s,n)$ by appending zeros onto the right.

By induction on $p+q$ we will show that for any non-negative integer $s$,

[TABLE]

for $\mu\in{\mathcal{C}}(s,p)$ and $\nu\in{\mathcal{C}}(s,q)$.

If $p=q=1$ the result is trivial since there is only one composition of $s$. Consider $p+q\geq 2$. Without loss of generality assume $p\leq q$.

We proceed by induction on $s$ (inside the induction on $p+q$). If $s=0$ the statement is vacuous, and so the base case is clear.

For positive integer $s$ let $J$ be a p-by-q non-negative integer matrix such that $J$ has row and column sums $\mu$ and $\nu$ respectively and $\langle J,C\rangle$ is minimal. We shall show that $\langle J,C\rangle={\mathcal{E}}(\mu-\nu)$.

If the first row (resp. column) of $J$ is zero we can delete it and reduce to the inductive hypothesis on $p+q$. Therefore, we assume that there is a positive entry in the first row (resp. column) of $J$. If $J_{11}>0$ then we can subtract $J_{11}$ from $s$ and reduce to the inductive hypothesis on $s$.

We are therefore left with $J_{11}=0$ and the existence of $i>1$, $j>1$ with $J_{1j}>0$ and $J_{i1}>0$. However, $(1,j)$ and $(i,1)$ are incomparable in the poset $[p]\times[q]$. However, by Proposition 4 we can assume that the support of $J$ is a chain.

The cost for the joint distribution, $J$, is minimized when the support is a chain. Upon inspection, one sees that the value of ${\mathcal{E}}$ agrees with the cost in the case that the support of $J$ is on a chain. The value of $J_{i,j}$ is multiply counted $|i-j|$ times – once for each of the contributing terms of ${\mathcal{E}}$. In the non-contributing terms there is a telescoping inside the absolute value.

Theorem 2 follows for chains, to which we have inductively reduced the problem. ∎

5.3. Proof of Theorem 3

Proof.

The vector space of degree $s$ homogeneous polynomial functions on the rank at most one p-by-q matrices is denoted ${\mathcal{R}}_{p,q}^{s}$. By Proposition 4, we obtain a basis for this space by considering the monomials

[TABLE]

where $J$ is a non-negative integer matrix with support on a chain. The row and column sums of $J$ are a pair of compositions of $s$ with $p$ and $q$ parts respectively. We denote these by $\mu$ and $\nu$. Note that by Proposition 5, $\mu$ and $\nu$ determine J.

If we assign each of these monomials the formal expression $z^{{\mbox{EMD}_{s}}(\mu,\nu)}t^{s}$ and sum them as formal series, we obtain the Hilbert series of ${\mathcal{R}}_{p,q}$,

[TABLE]

which we then recognize as the definition of $H_{p,q}(z,t)$.

Each monomial has a non-negative integer matrix $J$ as its exponents, with support on a chain. This chain terminates at or before $x_{p,q}^{J_{p,q}}N$. From Theorem 2 the variable $x_{p,q}$ is multiplied by $z^{|p-q|}\,t$, and contributes

[TABLE]

to all monomials. The geometric series sums to $\displaystyle\frac{1}{1-z^{|p-q|}\,t}$.

The preceding variables in the monomial may contain $x_{p,j}$ for some $1\leq j\leq q$, or $x_{i,q}$ for some $1\leq i\leq p$, but not both – since the exponent matrix has support in a chain. In the former case these monomials are in the sum $H_{p,q-1}$, while in the latter are counted in $H_{p,q-1}$.

The sum $H_{p-1,q}+H_{p,q-1}$ over counts monomials. That is to say, if a monomial has an exponent with support involving variables $x_{i,j}$ with $i<p$ and $j<q$ then it is counted once in $H_{p-1,q}$ and once in $H_{p,q-1}$. It also appears once in $H_{p-1,q-1}$. We therefore observe that such monomials are counted exactly once in the expression

[TABLE]

Finally, we see that all such monomials are counted exactly once in the product

[TABLE]

if $(p,q)\neq(1,1)$ and $H_{1,1}=\frac{1}{1-t}$. ∎

5.4. Proof of Theorem 1

Proof.

Fix positive integers $p$ and $q$. The coefficient of $t^{s}$ in $\frac{1}{(1-t)^{p+q}}$ is

[TABLE]

And, we have

[TABLE]

for some non-negative integers $c_{0},\cdots,c_{k}$. Thus the coefficient of $t^{s}$ in

[TABLE]

is therefore asymptotic to

[TABLE]

We will next find an inhomogeneous three term recursive formula for $N_{p,q}(1)$ from Equation 4.3.

First we observe that $W_{p,q}(1)=\binom{p+q-2}{p-1}$. Therefore, we have

[TABLE]

Then we divide by $(p+q-1)!$ to obtain the asymptotic. However, our goal is to obtain the expected value of ${\mbox{EMD}_{s}}$. So, in light of Proposition 6, we will need to divide by

[TABLE]

Thus, we find that the expected value is

[TABLE]

which we can rewrite as

[TABLE]

∎

6. Relation to spectral graph theory

We next turn our attention to some results from spectral graph theory and describe a connection with the expected value of the EMD.

The concept of a graph (or network) is likely familiar to the reader. We recall the terminology briefly. By a graph we mean an ordered pair, $(V,E)$, where $V$ is a finite set whose elements are called vertices and $E$ is a finite set whose elements are called edges, together with an injective mapping from $E$ to unordered pairs of distinct vertices. The elements of $E$ are said to join the corresponding pair of vertices. The number of vertices joined to a given vertex, $v\in V$, is called the degree of $v$, denoted $deg(v)$.

A sequence of distinct vertices, $v_{1},v_{2},\cdots,v_{t}$ with $v_{i}$ joined to $v_{i+1}$ for each $i$ is called a path. If a path exists between between all pairs of vertices then we say that the graph is connected.

Given vertices $v$ and $w$ in a connected graph, the distance between vertex $v$ and vertex $w$ is the length of the shortest path starting with $v$ and ending with $w$, and will be denoted $\rho(v,w)$. The function $\rho$ is a metric on $V$. For an integer $r$, the ball of radius $r$ centered at $v\in V$ will be defined as

[TABLE]

Furthermore, let $n_{r}(v)=|B_{r}(v)|-|B_{r-1}(v)|$, and set

[TABLE]

which is a finite sum giving the expected distance a vertex is from $v$.

As we shall see, this metric is related to the topic of this article. Specifically, the mean distance in a graph is defined to be:

[TABLE]

where $m=|V|$.

The above notation is from [M2], Section 3. Observe that the definition is equivalent to averaging the the distance between all two-element subsets of $V$. We also point out that if we consider all ordered pairs of vertices we have:

[TABLE]

From our point of view, we consider the graph to be on the vertex set ${\mathcal{C}}(s,n)$ where two compositions are joined when the (unnormalized) Earth Mover’s Distance is exactly 1. We call this graph the Earth Mover’s Graph, denoted $G(s,n)$. In this case the length of the shortest path between two vertices is the Earth Mover’s Distance.

The average distance between ordered pairs of vertices in $G(s,n)$ is the subject of this article. More generally, one can consider a graph $G$ whose vertices are a subset of ${\mathcal{C}}(s,n)$ and two vertices $v$ and $w$ are joined when ${\mbox{EMD}_{s}}(v,w)\leq t$ for some fixed constant $t\geq 0$. We will call $t$ the threshold. In practice, determining values of the threshold that uncover features in the data is an important research topic. Understanding the expected EMD in the uniform case is only one line of research.

The connected components of $G$ are of interest in cluster analysis. For example, the connected components of $G$ may be interpreted as clusters. If $t=0$ there are no edges in $G$, thus no connected components. For sufficiently large $t$, every pair of vertices is joined and there is only one component. As $t$ decreases the graph disconnects. Hierarchical connection of components defined by $t$ gives rise to many different types of clustering. An example of this type of analysis will be given in the next section.

There are several “off the shelf” methods for clustering analysis. See [DM] Chapter 20 for some commentary, especially about the popular “k-means” algorithm. Here we present only metric hierarchical clustering and spectral clustering as they relate to Theorem 1. However, the data that we consider in this article can and should be looked at from several points of view. In particular, unsupervised machine learning techniques are merited. See [HTF] as a general reference.

The average distance of a graph is also related to other invariants; we recommend the survey [M1]. First we recall some additional terminology. Given a graph $G$ with vertices $\{v_{1},\cdots,v_{m}\}$ and $k$ edges, one can form the Laplacian matrix $L_{G}=D_{G}-A_{G}$ where $D_{G}$ is the diagonal matrix with the degree of vertex $v_{i}$ in the i-th row and i-th column, while $A(G)$ is the adjacency matrix in which the entry in row i and column j is a one if $v_{i}$ and $v_{j}$ are joined by an edge, and zero otherwise.

The spectrum of $L_{G}$ is of interest. To begin, $L_{G}$ is a positive semidefinite matrix. The multiplicity of the 0-eigenspace is equal to the number of connected components of $G$. If the spectrum of $L_{G}$ is denoted by $0=\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{m}$, the algebraic connectivity is given by $\lambda_{2}$. Intuitively, we expect “clustering” when $\lambda_{2}$ is small relative to the rest of the spectrum.

Related to the algebraic connectivity are inequalities proved in [M2]. We recall them here because they partially describe the structure of the graph based on $\lambda_{2}$. In fact we can use $\lambda_{2}$ to calculate bounds on the average distance between vertices in a graph and another invariant to be defined next.

The discrete Cheeger inequality asserts that the isoperimetric number, $i(G)$, is closely related to the spectrum of a graph. This number is defined as

[TABLE]

where $\delta(X)$ is defined to be the boundary of a set of vertices $X$ (that is $v\in X$ iff $v$ is in $X$ but is joined to a vertex not in $X$). One result from [M2] is

[TABLE]

where $d_{max}$ is the maximum degree of a vertex in $G$. These results have their underpinnings in geometry and topology, see [Mgw], for example. Intuitively, the point here is that if $G$ has two large subgraphs that are joined only by a small set of edges then $i(G)$ is small. Unfortunately, computing $i(G)$ exactly is difficult. However, the spectrum of $G$ can be computed more easily, providing the stated bounds for $i(G)$.

A third invariant for $G$ is $\overline{\rho}(G)$, the mean distance. This value can also be bounded. The inequality presented in [M2] is

[TABLE]

where $m$ is the number of vertices in $G$.

7. Real world data

In this section we consider a real world data set coming from the Section Attrition and Grade Report published by the Office of Assessment and Institutional Research at the University of Wisconsin - Milwaukee for Fall semesters of academic years 2013-2017. We selected data for courses with enrollments greater than 1,000. Analysis is done for the 12 grades A through F (with plus/minus grading). “W” grades were not reported by the University for the entire period and are not included.

Most universities collect and analyze similar types of retention and attrition data. Analysis of grade distribution data may be an as-yet untapped portal for self-examination that reveals patterns in large, multi-section courses, or allows for insight into cross-divisional course boundaries. Cluster analysis provides a visual reinforcement of the calculated EMD.

The data for 2013-2017 comprises a total of twenty one courses. In Fall 2013 there are five courses: English 101 and 102, Math 095 and 105, and Psychology 101. Math 095 was redesigned after Fall 2013 and enrollment dropped below 1,000; the other four courses were offered every year. Input data are sorted by division and year.

[TABLE]

7.1. Histogram of EMD sample

We form the Earth Mover’s Graph by computing the unit normalized EMD for each pair of the 21 courses. A histogram of the results is presented in Figure 3. The rough structure reflects some aspects of the distribution of the theoretical case for 30 students and 5 grades depicted in Figure 2. For example, it is almost unimodal and skewed to the right. It is also interesting to note that the histogram in Figure 3 shows maximal EMD between courses at around 0.25. Additionally, for each course one can compute the average EMD to all others. This course pair-wise average has a minimum EMD of 0.067, mean 0.086, and maximum 0.129.

To better understand the relative sizes of these numbers, we can compare them with the mean in the uniform model with $n=12$, which is approximately 0.1300 when $s=1000$ and then drops to 0.1293 when $s\rightarrow\infty$. For actual grade distribution data, not all grade distributions are equally likely. Nonetheless, for this data set the maximum mean EMD corresponds to Psychology 101, Fall 2014 and is surprisingly close to the theoretical value. No doubt this is an anomaly, but we were surprised by how large the unit normalized EMD was between courses when compared to uniform sampling.

7.2. Hierarchical Clustering

Using the unit normalized EMD, we next determine a threshold, $t$, such that two courses are joined by an edge when the EMD falls below $t$. Thus, we obtain a family of graphs parameterized by $t$. We consider the connected component structure for each value of $t$. In practice, one sees a single giant component when the threshold is “large”. As a heuristic, we consider large to mean greater than the expected distance in the uniform model.

In Figure 4 we step through various distance threshold values $t$ and count the number of connected components. When $t$ is very small, the graph has 21 components. As $t$ increases, the components of the graph continue to connect, with the largest range of persistence for $t\in[0.0319,0.0477]$ with two components. It is useful to compare these values to the expected EMD in the uniform model. The endpoints of this interval are 25% and 37% of the uniform mean, respectively.

Set $t=0.0478$. We form a graph on the vertex set of these 21 courses. Two courses are joined by an edge when the unit normalized distance between them falls below $t$. Since $t$ is larger than 0.0477, there is only one connected component as depicted in Figure 5.

7.3. Spectral Analysis

A key component of spectral clustering analysis is the construction of the Laplacian Matrix $L_{G}$ as defined in Section 6. The second smallest eigenvalue of $L_{G}$ gives the algebraic connectivity. Figure 6 presents a plot of the full spectrum for $t=0.0478$.

For the 21 vertex graph in Figure 5, the algebraic connectivity is $\lambda_{2}\cong 0.1213$, and the maximum degree is $d_{max}=8$. It appears as if elements $2$ and $16$ create a tenuous bridge between two “clusters”.

The theoretical bounds on the isoperimetric number from Equation 6.1 are

[TABLE]

[TABLE]

with a computed value of $i(G)\cong 0.1$ for this data set. The relatively small isoperimetric number is consistent with the single edge of connection between elements $2$ and $16$ in Figure 5.

The theoretical bounds on the mean distance $\overline{\rho}(G)$ from Equation 6.2 are

[TABLE]

[TABLE]

and one can compute $\overline{\rho}(G)\cong 2.910$ for this data set.

The most curious part of this analysis appears in the plot of the components. At threshold $t=0.0477$ the algebraic connectivity (i.e. $\lambda_{2}$) goes to zero, and the EMD separates all of the English courses from the cluster of Math and Psychology courses. More specifically, EMD splits English (Course ID’s $1$ through $10$) off from Math (ID’s $11$ - $16$) and Psychology (ID’s $17$-$21$). Or equivalently, grade distributions in English courses are most similar to grade distributions in other English courses, and least similar to grade distributions in both Math and Psychology courses.

Although we see this clear partition of the data in Figure 7, one needs to be cautious about clustering algorithms. For example, clustering identified in one algorithm may not be the same as another. See the paper [K] for a careful treatment of this topic in general. In our more specific setting the clustering is determined by the threshold $t$. Thus, the question of where to set this value is delicate: too small and we see too many components, while too large we see very little clustering at all. Furthermore, our intention is only for exploratory data analysis and not to suggest policy regarding instructional assessment.