Diffusion maps

We begin with the basic definitions required to construct a diffusion map (Table 1 below summarizes our notation). This requires a set of data points, X, and a metric on those points, d. The data may be composed of many observables, making each x_k ∈ X a vector in a space of potentially very high dimension. Essentially we employ the algorithm in [25], with α = 1 and ϵ = 1. This involves

(i)constructing a similarity measure, K(x_i, x_j), between pairs of data points,

(ii)normalizing K to create a Laplacian operator $$ \widehat{A}$$ on the data,

(iii)computing eigenvector-eigenvalue pairs {(ϕ_k, λ_k)} for the Laplacian, and

(iv)constructing the diffusion map $$ \Phi (t,x)≔{\left({\lambda }_k^t{\varphi }_k\left(x\right)\right)}_{k=1}^n$$.

Table 1

List of symbols used.

Symbol Reference
$$	data set
n	number of data points
[k]	the set of integers 1, 2, …, k
$$	micro-scale observable
d(⋅, ⋅)	metric on X
D	pairwise distance matrix
A	affinity matrix
$$	heat operator
$$	diffusion distance on X
{ϕk}	eigenvectors of $$
{λk}	eigenvalues of $$
Φ(t, ⋅)	diffusion map at scale t
{ψk(t, ⋅)}	diffusion coordinates of Φ(t, ⋅)
ξi,j	inner product of $$ and $$

The choice of similarity measure in (i) is highly system-dependent, for agent-wise interactions may manifest themselves in many different ways and so require different metrics. In (ii), the choice of normalization in constructing $$ \widehat{A}$$ is subtler, and based on a choice of Laplace operator; this is not sensitive to the data under observation, but rather reflects an assumption of Manifold Learning, that the data points of X are drawn from a manifold embedded in euclidean space. For rigorous results on when such a manifold exists for a given dataset, see [26].

Given n points X = {x₁, …, x_n}, we construct an n × n matrix of pairwise distances, D = (D_i,j) ≔ (d(x_i, x_j)), associated to a metric d. Writing k(x) for the k^th nearest neighbor (kNN) of x ∈ X, let K(x, y) be the gaussian kNN (GkNN) kernel defined by

This kernel could be replaced with any other, as another function may better reflect the user’s impression of the similarity between two agents a distance d(x, y) apart. This particular choice is studied in [27], and serves to allow pairwise distances to scale according to their local geometries. One could consider this the standard gaussian kernel applied to the modified distance function $$ d(x,y)/\sqrt{d(x,k(x\left)\right)d(y,k(y\left)\right)}$$. The more sparse the point cloud (X, d) is near x or y, the smaller this distance becomes. This kernel also has the desirable property that it is scale-free, in that multiplying all distances by a fixed constant does not change the kernel K.

Let σ denote a row sum operator: if A is an m × n matrix, let σ(A) be the m-vector whose i^th component is $$ \sigma {\left(A\right)}_i≔{\Sigma }_{j=1}^n{A}_{i,j}$$. We next calculate an affinity matrix A = (A_i,j) ≔ (K(x_i, x_j)), its row sums σ(A), and construct the matrix operator A′ via

Fig 2

Conceptual illustration of diffusion maps applied to data on fish schooling.

A) For the same time period, different affinity operators can be calculated. These affinities could be correlations in simply spatial proximity (Ai) or correlations in velocities (Aii) for example. The operator is used to describe a network of agents. B) Once calculated, the heat diffusion operator allows us to quantify the geometry of the network and with this in hand C) diffusion map coordinates may be used to efficiently embed the fish in a low-dimensional space. Doing this for different metrics / affinities allows us to quantify the alignment between various diffusion map coordinate systems, which reveals changes in macro-scale behavior.

We remark that though increasing t decreases the resolution of small scale features in the embedding, in practice we take t = 1 throughout and define ψ_k(x) = λ_k ϕ_k, so that

Map alignment

In this section we formulate a statistic to measure the degree to which the macro-scale organization of two observables measured on the micro-scale level are related. Let X = {x₁, …, x_n} be the collection of data for the n interacting agents making up the system, and let f₁ and f₂ be two vector-valued observables, such as position or velocity, defined on each agent,

For each i = 1, 2 the diffusion coordinates $$ \left\{{\psi }_j^{\left(i\right)}\right\}$$ are ordered by decreasing contribution to the global organization of the network X under the heat operator $$ {\widehat{A}}^{\left(i\right)}$$. Therefore if there is a coherent large-scale organization on X manifest in the observables f_i, then the leading diffusion coordinates of Φ⁽¹⁾ should be determined in some manner by the leading coordinates of Φ⁽²⁾. To confirm this, we expand the unit vectors $$ {\varphi }_j^{\left(1\right)}$$ in the subspace generated by the top k diffusion coordinates of $$ {\widehat{A}}^{\left(2\right)}$$, i.e. we project the eigenvectors into the subspace spanned by $$ \{{\varphi }_1^{\left(2\right)},\dots ,{\varphi }_k^{\left(2\right)}\}$$ via the linear operator

Then $$ \Vert {\Pi }_k^{\left(2\right)}\left({\varphi }_j^{\left(1\right)}\right)\Vert $$ is a measure of how well the coordinate $$ {\varphi }_j^{\left(1\right)}$$ associated to f₁ is correlated with the large-scale network structure induced on X by the variable f₂. One attempts to choose k in a manner that ensures that the subspace spanned by $$ \{{\varphi }_1^{\left(2\right)},\dots ,{\varphi }_k^{\left(2\right)}\}$$ encompasses the relevant macro-scale coordinates induced by f₂ on X, but without k becoming so large that we begin to recover a significant portion of ϕ simply due to the large dimension of the image space. For example, if we take k = n the projection $$ {\Pi }_n^{\left(2\right)}$$ is the identity, simply expressing ϕ in the basis $$ \left\{{\varphi }_j^{\left(2\right)}\right\}$$. However, taking k ≪ n we recover only the portion of the vector which can be expressed as a linear combination of the top k diffusion coordinates of Φ⁽²⁾. We return to this issue below.

Now we derive the expected projection size under a null hypothesis that the two variables f₁ and f₂ are entirely unrelated. Given the pair of diffusion bases, we fix $$ \varphi \in \left\{{\varphi }_j^{\left(1\right)}\right\}\subset {\mathbb{R}}^n$$ and calculate the distribution of the random variable $$ \Vert {\Pi }_k^{\left(2\right)}\varphi \Vert $$ assuming that the k-dimensional range, $$ \text{Span}({\varphi }_1^{\left(2\right)},\dots ,{\varphi }_k^{\left(2\right)})$$, is chosen uniformly at random from the set of all k-dimensional subspaces. One can show, and with a little thought it is clear, that this is equivalent to asking for the magnitude of a random unit vector of $$ {\mathbb{R}}^n$$ projected onto its leading k coordinates. When the vector is sampled uniformly at random from $$ \{\psi \in {\mathbb{R}}^n:\Vert \psi {\Vert }_2=1\}$$, the quantity $$ {\Vert {\Pi }_k^{\left(2\right)}\varphi \Vert }^2$$ is a beta distributed random variable with parameters (n −k)/2 and k/2. Recognizing this, one can calculate the expectation and variance to be k/n and 2(n − k)k/(n³ + 2n²), respectively. Normalizing the projection size accordingly, we arrive at the statistic

Again, the leading diffusion coordinates of Φ⁽ⁱ⁾ are data-driven macro-scale variables, expressing the large-scale organization of the network as viewed through $$ {\widehat{A}}^{\left(i\right)}$$. If the data obeys a (possibly nonlinear and multi-valued) relationship between f₁ and f₂ on the macro-scale, then it should be evidenced by a relationship between the bases $$ \left\{{\varphi }_j^{\left(1\right)}\right\}$$ and $$ \left\{{\varphi }_j^{\left(2\right)}\right\}$$. Note that the observables f_i are used only in providing a metric on X. One can perform an identical analysis on the data given only the metrics d₁ and d₂. In the context of this paper we will assume, however, that the metrics come from some micro-scale variable defined on X, as this is the approach most useful in analyzing a new complex system.

Let $$ u=({\lambda }_1^{\left(1\right)},\dots ,{\lambda }_k^{\left(1\right)})$$ be the vector of decreasing eigenvalues of $$ {\widehat{A}}_1$$, and take $$ \widehat{u}≔u/{\Vert u\Vert }_2$$. Then the entries of $$ \widehat{u}$$ give a normalized notion of the energy of each $$ {\varphi }_j^{\left(1\right)}$$, and so gives a measure of how important each mode $$ {\varphi }_j^{\left(1\right)}$$ is to the macrostructure derived from f₁. Finally, we define the map alignment statistic to be the sum of the $$ {\mathcal{P}}_k^{\left(2\right)}\left({\varphi }_j^{\left(1\right)}\right)$$, weighted by the contribution of $$ {\varphi }_j^{\left(1\right)}$$ to the large-scale organization of f₁:

Informative choice of k

Here we briefly discuss how to choose k so that Z_k gives an accurate measure of the degree to which Φ⁽¹⁾ is determined by $$ {\Phi }_k^{\left(2\right)}$$. Recall that the head of the spectrum of $$ \widehat{A}$$ is sensitive to large-scale changes, but the bulk of the spectrum is associated to high frequency eigenfunctions, capturing small-scale organization often attributable to noise or local idiosyncrasies in the data. As a result, given training data or simulated data one may simply estimate a value of k after which the variance explained by the $$ j^{\text{th}}$$ eigenspace, $$ {\lambda }_j/{\sum }_{i=1}^n{\lambda }_i$$, is approximately independent of the macrostate of the system. That is, choose k so that for j > k, the value λ_j makes up a fixed fraction of the sum of the eigenvalues.

Alternatively, in using empirical data we found the bulk of the spectrum often follows an approximate power-law decay. As a result, one may take k to be the smallest integer such that for j > k, λ_j ≈ Cj^−β for positive constants C, β.

These are necessarily rough guidelines, and in practice one often finds a range of plausible choices for k. However, when increasing k the additional summands of $$ {\Pi }_k^{\left(2\right)}$$ must eventually become essentially independent ($$ \langle {\varphi }_i^{\left(1\right)},{\varphi }_k^{\left(2\right)}\rangle \approx 0$$) of the leading elements of $$ \left\{{\varphi }_i^{\left(1\right)}\right\}$$. It follows that after a point, a larger value of k has little influence on the computed MAS, Z_k(f₁, f₂). Therefore an overestimate of k has very little effect on the results, so it is advised that in practice one err on the side of larger choices of k.

Unsupervised macrostate classification

If the system under study exhibits multiple emergent behaviors relating the two variables f₁ and f₂, then coherence alone is not enough to identify system behavior. In passing to the norm of the vector $$ {\Pi }_k^{\left(2\right)}{\varphi }_j^{\left(1\right)}$$, we clearly give up a significant amount of information. In order to distinguish multiple coherent states, we inspect the individual dot products $$ \left\{\langle {\varphi }_i^{\left(1\right)},{\varphi }_j^{\left(2\right)}\rangle \right\}$$ that comprise the map alignment statistic.

However, to remove the effects of the trailing diffusion coordinates, which correspond to smaller scale variations in the geometry of $$ (X,{\mathcal{D}}^{\left(i\right)})$$, we consider the sequence of inner products given by

While the diffusion coordinates provide a natural organization of the data, the individual components of the map are eigenvectors, and so are computed via an iterated maximization problem: given the linear operator A, the i^th eigenvector is given by

Therefore, for each collection {ξ_i,j: 1 ≤ i, j ≤ k} we define the map σ: [k] × [k] → [k²] so that the components of (|ξ_{σ(i, j)}|) are ordered by decreasing magnitude. (Here [k] denotes the first k positive integers). This removes the ambiguity in both index and sign for these products. We call the resulting vector the f₁ f₂-covariance vector (at time t), abbreviated as ξ (resp. ξ(t)). One may then use these vectors as a fingerprint for the macro-scale relationship exhibited between f₁ and f₂ at a given time.

Fig 3, panel A, gives an example of a toy flocking system comprised of 50 agents, exhibiting a circling behavior at time t₁ and a (mostly) aligned state at time t₂; its top four diffusion coordinates associated to both Φ⁽¹⁾ and Φ⁽²⁾ are displayed for each time. (The choice of affinities $$ {\widehat{A}}^{\left(1\right)}$$ and $$ {\widehat{A}}^{\left(2\right)}$$ here are the same as in the flocking model below, being based on spatial proximity and velocity, respectively). In panel B we see the absolute values of the inner products of the diffusion coordinates, where the matrix (represented by a color plot) has its (i, j)-entry given by $$ \langle {\psi }_i^{\left(1\right)},{\psi }_j^{\left(2\right)}\rangle $$. Finally, in panel C we see the covariance vectors at times t₁ and t₂.

Fig 3

A toy system of 50 agents; the left and right sides of the figure represent observations of the system at two different times when the system is in two different organized states.

(A) Agents in each subpanel are colored by the indicated diffusion coordinate, (B) A color plot of the matrix (ξ_i,j) with yellow indicating larger values, blue indicating smaller values, (C) The values of the above inner product matrix, ordered by absolute magnitude; these are the vectors used to classify the system’s regimes.

Covariance vectors allow us to embed the activity of the agents in a shared space, $$ {\mathbb{R}}^{k^2}$$, and perform clustering analyses to group similar macro-scale relationships between f₁ and f₂ over time. In the next section we will perform a boilerplate clustering analysis to distinguish between distinct coherent group structures found in a data set of golden shiner schooling behavior, a task for which the map alignment statistic is insufficient on its own.

Analysis of bird flocking simulations

As a first test of our method we explored a well studied agent-based model of flocking behavior [28] (see S2 and S3 Movies for a visual reference). The model consists of a set of N agents, {x_j}, moving on the torus T¹ = S¹ × S¹, each with position p_i and velocity determined by the angle θ_i and the magnitude v_i. In this model each agent moves with a fixed speed, v_i = v_j for all i and j, so we simply write v for this speed.

Each agent’s direction is updated at each time step by the average local direction plus a noise term, a normal random variable distributed as N(0, ϵ²). Here local refers to the agents within a certain distance r of the focal agent. With all other parameters of the model fixed, the coupling parameter r governs the degree to which the system self-organizes. In the hydrodynamic limit, long-range order emerges for r above some critical value r_c(v, ϵ) > 0, meaning that the mean velocity of the flock is nonzero: $$ {\displaystyle \underset{N\to \infty }{\text{lim}}\Vert \frac{1}{N}\sum _{i=1}^N(\text{cos}{\theta }_i,\text{sin}{\theta }_i)\Vert \ne 0}$$. The direction of group travel is chosen randomly, and breaks the symmetry of the system. See S1 Text for model details.

In the discrete time, finite-agent setting we observe the mean velocity of the flock, 〈v_j〉 = (v cos〈θ_j〉, v sin〈θ_j〉), taking this as our ‘ground truth’ measurement of the level of emergent, coherent behavior in the system. As r increases past r_c(v, ϵ), the mean velocity (in the large N limit) transitions from zero to a nonzero value. However, this shift occurs continuously; as r grows we find the steady-state mean velocity of the flock increases towards its maximum value. Even with r = 0 the global mean velocity will have magnitude $$ O(1/\sqrt{N})$$ and suffer fluctuations, so determining whether a preferred flock-wide direction of motion has been chosen in the finite case makes little sense. Instead, below we compare the value of Z_k(p, v) to 〈v_j〉 to determine how well the MAS detects large-scale organization.

We chose two metrics to define a pair of time-evolving graphs as follows: Since we are interested in flocking behavior, we let d⁽¹⁾(x_i, x_j; t) be the euclidean distance between the $$ i^{\text{th}}$$ and $$ j^{\text{th}}$$ agents at time t, while d⁽²⁾(x_i, x_j; t) denotes the kernel-smoothed L² norm of the two agents’ relative speed profiles. (See (S3), (S4) and (S5) for the explicit definitions of the metrics used). In particular, d⁽¹⁾ captures the spatial information in the flock, while d⁽²⁾ measures the degree that agents are moving in concert over a short (on the order of 10 time steps) recent period.

We then applied diffusion maps to the resulting sequence of graphs and calculated the map alignment between the two variables, p and v. The choice of subspace size to use, k = 10, was empirical. We observe that at various levels of organization the variance explained by the first five modes (i.e. the sum of the magnitudes of the first four components of $$ \widehat{u}$$) is highly variable, while the variance explained by the higher modes ($$ {\varphi }_6^{\left(i\right)},{\varphi }_7^{\left(i\right)},\dots ,{\varphi }_{500}^{\left(i\right)}$$) shows little dependence on system organization, though we decided to increase k to 10 to be conservative. This suggests that the macro-scale information is well captured by $$ {\Phi }_{10}^{\left(i\right)}$$, while components 11 through 500 are largely accounting for noise and small-scale features of $$ {\widehat{A}}^{\left(i\right)}$$. Note, however, that one may take k larger than 50 without any qualitative change in the analysis.

In Fig 4 below we plot the map alignment statistic Z₁₀(p, v) (blue) calculated from two sets of 100 independent simulations of the flocking model with 500 agents and fixed parameters (ϵ = π/5, v = 1/320, see SI). In the left panels the coupling parameter r is held at a low value (r = r₀ = 0.0065) for 800 time steps before increasing linearly to a high value over 400 steps and remaining at that value (r = r₁ = 0.06) for the rest of the trial; see the lower panel of Fig 4A. The mean velocity 〈v_j〉 of the flock, averaged over the 100 trials, is plotted in orange in the upper panel.

Fig 4

Top: Map alignment statistic Z₁₀(p, v) averaged over 100 flocking simulations of 2000 time steps (blue). The grey region denotes the 90% empirical confidence interval of Z₁₀(p, v). The mean velocity 〈v_j〉 of the flock is averaged over the 100 trials and plotted in orange. Bottom: The model’s coupling constant r plotted over time. Notice that in (B) the MAS does not suffer the same hysteresis that the flock’s speed does.

One can see that the statistic Z₁₀(p, v), measuring the dependence of position diffusion coordinates on the velocity diffusion coordinate system, has a strong correspondence with the velocity correlation statistic. That is, large-scale organization was clearly detected and quantified using only the two micro-scale inter-agent distance time series.

In a second trial we explored how quickly map alignment and long-range correlations vary in response to change in the micro-scale dynamics of the agents. This is done by allowing r to vary sinusoidally between r₀ and r₁, with increasing frequency. The black line plotted in the lower panel of Fig 4B shows the coupling parameter as a function of time; in the upper panel the orange line shows the magnitude of the mean velocity of the flock, and the blue plot gives the mean value of Z₁₀ across trials, with the filled area indicating the 90% confidence interval for Z₁₀.

Clearly Z₁₀(p, v) indicates strongly the presence or absence of correlation. The map alignment statistic shows both a faster response to changing r and suffers from less attenuation as the speed of parameter change increases. Times t = 800 to 1000 especially highlight how the drop in Z₁₀ is a precursor to the value of the macrovariable 〈v_j〉 falling. In this case and the previous case, we find we are able to track the transition from unorganized motion to flocking by the alignment between the velocity diffusion coordinates, $$ {\Phi }_{10}^{\left(1\right)}$$, and the position diffusion coordinates, $$ {\Phi }_{10}^{\left(2\right)}$$.

Analysis of empirical fish schooling data

As an empirical application we use position and velocity data for a school of golden shiners (notemigonus crysoleucas); details of how the data were collected can be found in S1 Text (see also S1 Movie for a visual reference to these data). These fish were studied in [29, 30]. In [30] the authors found that the school, constrained to an essentially two dimensional environment, would generally be found in one of three macro-states defined according to the school’s collective polarization, O_p, and rotation, O_r: milling (circling), swarming (disordered, stationary), or polarized (translational). See S1 Fig for example frames of these behaviors. Writing v(x_i) for the velocity of fish i and p(x_i) for its position, the chosen macrovariables are defined as

Fig 5

Statistics calculated from 5000 frames of the golden shiner dataset.

A): Map alignment Z₁₅(p, v) calculated over 5000 video frames. The plots are colored according to our classification of macro-scale state; milling = red, swarming = yellow, α-polarized = blue, β-polarized = green. B): The group polarization macro-variable O_p calculated from the same 5000 frames. The coloring shows the states according to [29]; milling = red, swarming = yellow, polarized = teal, unclassified/transitional = grey.

For the study of fish schooling, as with the simulated data above, we built a pair of distance matrices at each time step. The first was calculated from the euclidean distances between agents (individual fish), d⁽¹⁾, defined as before. The second, d⁽²⁾, using the velocity data only. In this case, our goal was to capture any interaction between agents, either instantaneous or at a time lag due to one fish leading another. Therefore we define d⁽²⁾(x_i, x_j;t) as in the previous section (a weighted L² norm of the differences in the agents’ recent velocity profile), but now allow one of the two agent’s velocity profiles to begin up to two seconds earlier. This amounts to starting the trajectory at an earlier time step, t − ℓ, with 0 ≤ ℓ ≤ 60. We take d⁽²⁾(x_i, x_j; t) to be the minimum of this collection of distances, as this choice of time translation gives the best alignment of the two fishes’ trajectories.

As with the flocking model, the first diffusion map’s leading coordinates identify the large groups of fish who form spatially proximate clusters; the leading diffusion coordinates associated to d⁽²⁾ highlight groups whose velocity profiles are (up to a 2s, or 60 frame, translation in time) quite similar. The map alignment between these two systems, Z₁₅(p, v), is plotted in the upper panel of Fig 5. When the fish are swarming (see S1 Fig, center column) we see that these two coordinate systems are unrelated, and that Z₁₅(p, v) is correspondingly low.

Here again, the choice of k = 15 is driven by observing that outside of the top fifteen entries of $$ \widehat{u}=u/\Vert u\Vert $$ for $$ u=({\lambda }_1^{\left(i\right)},\dots ,{\lambda }_n^{\left(i\right)})$$, i ∈ {1, 2}, the spectral decay appears to follow a power law for most frames (though this observation does break down for some frames when the eigenvalues become small enough). The results below prove to be robust to varying k from 10 to 25. Therefore the span of $$ \{{\varphi }_1^{\left(1\right)},\dots ,{\varphi }_{15}^{\left(1\right)}\}$$ contains the majority of the structure within the data, while k = 15 is not so large that the signal in Z_k becomes diluted.

The two coherent behaviors we would like to separate are the polarized, or linear, motion (right column in S1 Fig) and milling, or circling, motion (left column in S1 Fig). Both states exhibit high levels of organization, and we must use a sharper tool to distinguish them.

To do this we employ the covariance vectors ξ_i,j, as defined in the Unsupervised Macrostate Classification section above, over the course of 5000 frames of video. First, those states with Z₁₅(x, v) < 3 are considered to exhibit incoherent behavior and to have no macro-scale dynamic, so we classify them as swarming frames and remove them from the set of frames to be classified. (NB: the summands $$ {\mathcal{P}}_{15}^{\left(v\right)}({\varphi }_j^{\left(x\right)})$$ for 1 ≤ j ≤ 15 are centered, unit variance random variables if we assume the two coordinate systems are unrelated. For this reason the threshold of 3 can be considered a reasonably large deviation from the mean for Z₁₅(x, v). However, the choice is ad hoc, and users may modify this threshold in practice if it enables them to better subdivide their system into different behavioral regimes). Then for each time t we compute the covariance vector ξ(t) = (|ξ_σ(i,j)(t)|). We cluster the resulting vectors using the k-means algorithm, with k = 3. This leaves us with three groups representing distinct coherent behaviors, labelled G1, G2, and G3, as well as the incoherent group, which we label N. These groups define the coloration of the top panel of Fig 5. Comparing these group memberships to those given by calculating the group polarization and rotation, O_p and O_r, and classifying behavior as in [29] we find that group G1 contains 85.86% of the milling frames, G2 and G3 together comprise 87.47% of the polarized frames, while N holds 65.85% of the swarming frames (another 32.01% are relegated to G2). The full results of the classification appear in Table 2.

Table 2

The percentage of the frames of a given state (milling, polarized, etc.) which is found in each of our unsupervised groups (G1, G2, G3, N) are displayed below the counts.

Counts of classifications by our map alignment-based statistics (rows) and by the macro-scale variable classification (columns) for the 5000 frames of fish movement data.

	Milling	Swarming	Polarized	Transitional
G1	1463	5	73	348
	(85.9%)	(1.5%)	(5.1%)	(22.5%)
G2	0	105	738	454
	(0%)	(32.0%)	(51.9%)	(29.4%)
G3	241	2	505	352
	(14.1%)	(0.6%)	(35.5%)	(22.8%)
N	0	216	105	391
	(0%)	(65.9%)	(7.4%)	(25.3%)

That the polarized frames are divided into two separate clusters suggests that there are two distinct modes of polarized behavior exhibited by the fish. Upon inspecting ξ(t) for the frames from each cluster, we find that the two groups are essentially distinguished by their relative levels of organization; for example, the mean value of Z₁₅ for cluster G2 was 6.3186, while for G3 the mean map alignment statistic is 12.0878. Practically, the frames of cluster G2 typically consist of polarized group motion with multiple subgroups of fish following distinct paths; we call this β-polarized behavior. This creates a macro-scale division in the network structure of $$ {\widehat{A}}^{\left(v\right)}$$ that is not present in the network structure of $$ {\widehat{A}}^{\left(x\right)}$$, which leads to a lower alignment between their top eigenvectors; here 〈ξ_1,1(t)〉 = 0.176. On the other hand, the frames of G3 feature a single, unified group motion, leading to a high alignment (〈ξ_1,1(t)〉 = 0.331) between the top diffusion coordinates of each network. We refer to this as α-polarized schooling behavior. It is possible that there is a correspondence between the empirically defined β-polarized activity and the dynamic parallel group behavior described in the family of models of [31].

One may approach the problem of classifying the covariance vectors {ξ(t)} in numerous ways; we chose the k-means algorithm for its simplicity. However, it requires the user to fix the number of desired clusters beforehand. As the prior work of [29] only defined two separate coherent behaviors, milling and polarized, one may wonder what results we would find if we instead restricted the classification to two clusters. We did this and, labelling the groups G1* and G2*, found that G1* successfully captures all the milling frames, and G2* is predominantly composed of polarized frames (see S1 Table). However, roughly 1 in 4 polarized frames are relegated to G1*. S2 Table shows how the two k-means clusterings apportion the frames. In particular, the α-polarized frames tend to be classified as milling. By increasing the number of subgroups to k = 3, we allow the algorithm to perform a more sensitive clustering, which improved the correlation between our groupings and those of [29], as well as distinguished between two forms of polarized behavior that were grouped together in that work. (Note that while one could increase the number of clusters k to four or five, one experiences diminishing returns in the sense that the differences in schooling behavior will become increasingly subtle).