# Investigating the ‘old boy network’ using latent space models

Ian Hamilton, University of Warwick, Coventry, U.K

## Abstract

This paper investigates the nature of institutional ties between a group of
English schools, including a large proportion of private schools that might be
thought of as contributing to the ‘old boy network’. The analysis is based on
a network of bilaterally-determined school rugby union fixtures. The primary
importance of geographical proximity in the determination of these fixtures
supplies a spatial ‘ground truth’ against which the performance of models is
assessed. A Bayesian fitting of the latent position cluster model is found to
provide the best fit of the models examined. This is used to demonstrate a
variety of methods that together provide a consistent and nuanced
interpretation of the factors influencing community and edge formation in the
network. The influence of homophily in fees and the proportion of boarders is
identified as notable, with evidence that this is driven by a community of
schools having the highest proportion of boarders and charging the highest
fees, suggestive of the existence and nature of an ‘old boy network’ at an
institutional level.

Keywords— clustering; community detection; latent position cluster model;
relative importance; schools; sport

##  1 Introduction

‘Old boy network’ is an English phrase used to refer to the informal system
through which men assist other men of a similar socio-economically privileged
background, reflected in attending the same school or university. It derives
from the term, ‘old boy’, used to refer to a former pupil at a British ‘public
school’, the confusing name given to a subset of traditional, and often
especially expensive, private schools in England. Generally the old boy
network is considered to act at an individual level, but it may be reasonable
to think that it will be stronger if and where there are institutional links
between schools. However, it is not clear to what degree there are such links
and, to the degree that there are, what features of the school might define
them. In the English school system, there are a number of distinctions that
might be thought to be identifiable in the functioning of such a network, for
example private/state, boarding/day, level of fees. In this paper, latent
space models are applied to a network of school rugby union fixtures to
investigate whether such functional networks exist and, if they do, what the
nature of them may be.

The fixture data provides an interesting data set on which to apply these
approaches, for a couple of reasons. First, since one may reasonably expect
that a primary consideration in agreeing a fixture would be geographical
proximity, then there is a spatial ‘ground truth’ against which to assess the
models considered. Second, there is reason to believe that the set of fixtures
may be informative. Fixtures change from year to year but not drastically, and
the process of fixtures being scheduled and changed has been taking place for
many decades, with the first school rugby union fixture taking place over a
hundred years ago. Together these considerations suggest that there has been
time for relevant factors to exert an influence such that the observed
situation represents a steady state with respect to the schools’ current
relationships, allowing inferences to be informative.

Network latent space models position nodes in an unobserved latent space, with
the probability of an edge existing between any two nodes being related to the
proximity of the two nodes in the latent space. Modelling in this way can
allow a number of features common to networks to be captured — transitivity;
homophily by attributes; and clustering — as well as often allowing for
informative graphical representations.

In the context of networks, transitivity is the phenomenon that two nodes,
which each share an edge with the same third node, will have a higher
probability of having an edge with each other than a pair that do not.
Homophily by attributes describes the greater propensity for an edge to exist
when two nodes share observed attributes. For example in a friendship network
these could be attributes such as age, sex, geographical location, or
recreational interests. Transitivity and homophily by attributes will both
lead to clustering, but it is not uncommon to observe clustering beyond what
may be explained by these features. This may be due to homophily by unobserved
attributes, self-organisation of actors, the popularity of particular actors,
or endogenous attributes such as position in the network (Handcock et al., ,
2007).

Early latent space network models used multidimensional scaling, and while
these captured transitivity and homophily by attributes they were reliant on
the arbitrary choice of a distance measure, leading to variable
interpretations. Hoff et al.,  (2002) proposed a stochastic model in which the
latent space positions may be estimated through standard statistical
techniques. Handcock et al.,  (2007) extended that model to account for
clustering by assuming that the latent positions are drawn from a finite
mixture of multivariate normal distributions.

In the present setting there is reason to believe that there will be
clustering beyond that due to transitivity or observed homophily of
attributes. First, because there are likely to be factors not captured in
available data that account for unobserved homophilies, such as sport
orientation within school culture, and social networks between relevant staff.
Second, because the fixture information is publicly accessible, which might
increase the potential for further clustering through self-organisation; if a
school sees that a number of the schools they play have an opponent in common
who they do not play they may consider proposing a fixture.

The paper proceeds in Section 2 by discussing the data, with details of its
collection, and a brief analysis highlighting some features and inter-
relations of the school covariates considered. In Section 3 the latent space
models of Hoff et al.,  (2002) and Handcock et al.,  (2007) are fitted. In
Section 4 the influence of the different covariates is explored. Section 5
provides some concluding remarks.

##  2 Data

###  2.1 Background

The Daily Mail Trophy is an annual tournament between some of the best school
rugby teams in England, along with a single school in Wales. In order to
qualify for a ranking in the Daily Mail Trophy a school must register for the
tournament and compete against a minimum of five other teams in the
tournament. There are in total 118 schools included, playing between 3 and 37
matches over the course of the three seasons. The matches played as part of
the tournament are typically only a subset of the matches played by these
schools, with other matches taking place as friendlies with non-tournament
schools, or as part of a centrally scheduled knock-out competition. In almost
all cases however, if they have played a match against one of the other
tournament teams, outside of the centrally scheduled knock-out competition, it
would be classed as a Daily Mail Trophy match and would appear in the data
set. So the existence or absence of a match with another school is an accurate
representation of the bilaterally arranged fixtures within the set of schools
in the tournament.

The network to be investigated is based on all matches in the Daily Mail
Trophy in the 2015–16, 2016–17 and 2017–18 seasons. Schools who registered but
were unable to complete five eligible matches in any given season are
maintained in our analysis, despite being excluded from an official ranking.
The network being considered is undirected; each node is a school, and each
edge has a value equal to the number of seasons, out of the three for which
there are data, in which the two schools played each other. This network is
considered in relation to school-level data representing variables that could
contribute to some observable homophily.

###  2.2 Data collection

The fixtures and results for the Daily Mail Trophy were kindly shared by
Schoolsrugby.co.uk, the organiser of the tournament, but are also available at
the tournament website (www.schoolsrugby.co.uk/dailymailtrophy.aspx). In the
analysis that follows it is the number of seasons, out of the three
considered, that two teams play each other that are used, removing the five
instances where teams played each other twice in a season.

For each school, data were also collected on:

  1. 1.

annual fees (Fees)

  2. 2.

year of foundation (Founded)

  3. 3.

the number of boys in sixth form (6th Form boys)

  4. 4.

the proportion of pupils in the school that are boys (Percent boys)

  5. 5.

the proportion of pupils that are boarders (Percent boarders)

  6. 6.

whether the school is privately or state funded (School type)

  7. 7.

whether the school played one or two terms of rugby (Term type)

  8. 8.

performance rating of the rugby-playing strength of the school (Rating)

While it would be plausible to consider other variables, these were either not
readily available, for example the size of sports bursaries or proportion of
pupils from overseas, or were substantially accounted for by this set of
covariates, for example if a school is co-educational throughout or just in
sixth form, or not at all, is well captured in the proportion of pupils that
are boys and the number of sixth form boys. All covariate data were collected
during the first week of June 2019. This means they are not contemporaneous
with the fixtures occurring in the period 2015–2018. This is not expected to
have materially impacted any conclusions because the covariates are not
subject to large year on year changes, .

Fees information was sourced from the schools’ own websites. Annual fees were
taken to be the minimum fees for full-time education of a pupil in Upper Sixth
form, not accounting for bursaries or scholarships. Thus it is zero for state-
funded schools, standard day fees for schools that admit day pupils, and
minimum boarding fees for schools that admit only boarders.

Year of foundation was generally more readily available on the Wikipedia page
for the school than the school’s website and so the Wikipedia date was used.
This was originally included as it was suspected that it might define a
meaningful similarity between schools. But in collecting the data it became
clear that the trajectories of schools were very diverse and it was often even
difficult to be sure of a definitive foundation date with occurrences of, for
example, schools moving geographically, amalgamating, moving from private to
state or vice versa, and renaming not uncommon. All of these contributed to
subjectivity in the definition of year of foundation. It is included in the
analysis however as it provides a useful sanity check to some of the later
methods in the degree to which those methods identify year of foundation as a
non-informative covariate.

The data on the number of pupils, including the proportions of boys and
boarders, and the absolute numbers of sixth form boys, were generally not
available from schools’ websites. For privately funded schools, these data
were collected from the Independent Schools Council (ISC) website
(www.isc.co.uk). For boys and for girls, this reports the number of boarders,
the number of day pupils, and the number of sixth form pupils. From these the
proportion of day pupils, and the proportion of boys are calculated. However
it should be noted that since schools admit pupils from different ages, the
relative proportions as they pertain to a consistent age group, say 13–18, are
unlikely to be the same. Given the purpose of including these proportions —
identifying homophilies — and that no ready alternative was available, then
this seems acceptable, as they still provide evidence on the nature of the
school. Pupil numbers for state schools were taken based on the most recent
relevant government report (Department for Education, , 2019). These did not
include a delineation for numbers of sixth form students or of boarders. For
this, the most recent available Ofsted report with such data was used
(www.gov.uk/government/organisations/ofsted). Since these were from previous
years and so total numbers differed, the proportions of sixth form students or
boarders were assumed to be constant, and the absolute number of sixth form
boys was adjusted for the current total number of pupils. The total number of
pupils was never materially different from current numbers, and so one may
reasonably be confident that these numbers are accurate. However it should be
noted that these Ofsted reports were quite commonly from as much as a decade
ago. In one case, Ampleforth College, the population information was not
available on the ISC website and so the method used was the same as for the
state schools, with the most recent Independent Schools Inspectorate report
(www.isi.net) used, in place of an Ofsted report, to derive the relevant
proportions. Whether a school was state or privately funded was determined by
its classification in the most recent relevant government report (Department
for Education, , 2019).

In order to determine if the school played one or two terms of rugby, the
month of the school’s first and last fixtures of the 2018 season were
considered, including those fixtures outside the Daily Mail Trophy
competition. If these both lay in a single term, then it was interpreted as
being a single term rugby-playing school. If there were two or more matches in
the second term then it was deemed to be a two term school. If there was a
single match outside of first term then previous seasons’ fixtures as well as
any information from the school’s website was used to make the categorisation.

The performance rating used here is described in Hamilton and Firth,  (2021),
applied to an aggregation of the three seasons’ results. The method accounts
for the varying schedule strengths of participating schools in a manner
consistent with the predominant league points system used in rugby union.
Taking the projected league points per match were each team to play every
other team home and away in a round robin format provides a positive-valued
measure of a comparable scale to the other covariates as well as an intuitive
interpretation to the rating measure.

For all of these factors, when considering nodal covariates the absolute
difference is used. The binary variables of school type and term type take
value 0 if identical and 1 if different for each pair. The two percentage
variables of proportion of boys and proportion of boarders are taken as the
absolute difference in these percentages.

Postcodes were taken from the Daily Mail Trophy website and then used via a
Google Maps API in order to calculate travel times and distances between
schools using the R package googleway (Cooley, , 2017) and to project
locations onto a relevant map. Distances were calculated as at Saturday 5th
October 2019 12pm using the “best guess” methodology, assuming a journey by
road. In order to plot schools geographically, latitude and longitude for each
postcode was sourced from www.freemaptools.com.

The five continuous variables are scaled to allow for better comparability and
interpretability, specifically the following units are used:

Travel Time |  hours
---|---
Fees |  £10,000
Founded |  centuries
6th Form Boys |  100 boys
Rating |  projected league points per match

###  2.3 Exploratory Data Analysis

The two binary covariates of whether a school is privately or state funded,
and whether it plays one or two terms of rugby are compared first. Table 1
shows that there is a clear relation, with the private schools substantially
more likely to play one term and state funded schools more likely to play two.

| One term | Two terms
---|---|---
Private | 79 | 17
State | 8 | 14
Table 1: Number of schools in the tournament playing one or two terms of rugby
and private or state funded Figure 1: Scatterplots of school covariates.
Private schools in grey, state funded schools in black

The other covariates are presented in Figure 1. A number of things may be
noted. First the individual covariate distributions are largely in line with
expectations. A number of the covariates are bimodal. This is unsurprising in
the case of fees, with one of the modes at zero, or in the case of the
percentage of boys, with modes close to 50% and 100%, reflecting a
predominance of all-through co-educational and single sex schools, but the
bimodal nature to the year of foundation is less intuitive without more
historical context. The number of sixth form boys is right-skewed with a mode
around 150. The proportion of boarders has a clear mode at just above zero,
reflecting a high proportion of day schools and some schools with just a
handful of boarders. The remaining proportions are, perhaps not so
intuitively, distributed quite evenly up to full boarding status, but without
another clearly identifiable mode. Rating is left-skewed with its mode at
around 2.7 league points per match, which suggests teams would be expected to
share almost one and a half bonus point per match on average. Importantly for
our analysis, with the exception of the percentage of boarders and fees, and
the proportion of boys and the number of sixth form boys, which have Pearson
correlation of 70% and 68% respectively, the continuous covariates all have
absolute correlations of less than 25%, which is helpful for being able to
discern independent effects. The correlation between the proportion of pupils
that are boarders and the fees is not due, as one might first suspect, to
boarding fees being higher as they must also account for living expenses,
since the minimum upper sixth form fees have been used in all cases, and there
are only four pure boarding schools where those minimum fees include boarding
fees. In all other cases they represent a day fee. A correlation of 68%
between the proportion of pupils that are boys and the number of sixth form
boys is less surprising. Here shading has been used in the scatterplots to
differentiate state funded and private schools. Apart from the self-
explanatory difference in fees, these also show that within this set of
schools the state schools are more likely to have high numbers of sixth form
boys, to be single sex (or close to), and to be day schools rather than
boarding. We might therefore expect some confounding of these factors in later
analyses. These charts were also looked at with a differentiation based on the
number of terms of rugby played. This highlighted similar features, as one
might expect from Table 1, but less strongly.

##  3 Latent Space Model

###  3.1 Model specification

Given the three seasons of fixture data, a binomial latent space model is
fitted here, the general form of which is

$\displaystyle P(A\mid\boldsymbol{Z},\boldsymbol{x},\boldsymbol{\beta})$
$\displaystyle=\prod_{i<j}\binom{3}{A_{ij}}\mu_{ij}^{A_{ij}}(1-\mu_{ij})^{3-A_{ij}}$
(1) $\displaystyle\text{logit}(\mu_{ij})$
$\displaystyle=\beta_{0}+\sum_{k=1}^{p}x_{ijk}\beta_{k}-d(\boldsymbol{Z_{i}},\boldsymbol{Z_{j}}).$
(2)

Here $A=[A_{ij}]$ is the symmetric adjacency matrix of fixtures with $A_{ij}$
equal to the number of seasons out of the three in which teams $i$ and $j$
played each other, $x_{ijk}$ is the $k$th nodal covariate for teams $i$ and
$j$, $\beta_{k}$ is the coefficient for the $k$th covariate,
$\boldsymbol{Z_{i}}$ is the latent position for team $i$, and
$d(\boldsymbol{Z_{i}},\boldsymbol{Z_{j}})$ is the Euclidean distance between
teams $i$ and $j$ in the latent space. This is the binomial version of the
model proposed by Hoff et al.,  (2002).

It is worth noting that a sociality parameter is not included. In this context
a sociality parameter would describe the propensity for a particular school to
have matches. It is not included here because the fixtures represented are an
incomplete set of fixtures for the participating teams, with teams generally
playing additional matches, against teams outside of the tournament, as
friendlies or as part of other competitions. While these data are not fully
available, the schedule for a sample of teams has been inspected and generally
they have played a similar total number of matches, so observing a team to
have higher degree within the network is not a reflection of that team having
a higher propensity to play matches in general. Including such effects could
therefore, for example, misleadingly diminish the extent to which one might
infer a lower homophily from the absence of a fixture in the case of teams
with lower degree. On the other hand, the current rules of the tournament
encourage teams to play as many other tournament participants as possible,
since they are awarded bonus points merely for playing matches in the
tournament, independent of result. Based on conversations with the tournament
organisers and observations of historic results, it is likely that there will
be a difference in motivation between teams in their desire to be competitive
in the tournament, and some may actively seek to schedule more tournament
matches. However there are only a small number of teams, consistent across
years, for whom the tournament is a goal in itself. For the vast majority,
they enter simply because they can, as a by-product of their standard fixture
list, which is largely similar from year to year. So while this effect could
be argued to be a genuine sociality effect, within the context of just this
network of fixtures, it is likely to apply to only a small minority of teams.
Therefore on balance, the distorting impact of inclusion is considered to be
more of a danger than that of exclusion and so no sociality effect is
included. All models are fitted using the latentnet package in R (Krivitsky
and Handcock, , 2008).

###  3.2 Hierarchical Clustering

Initially the parameters of the model represented in equations (1)–(2) are
estimated through the method of maximum likelihood and with no covariates
included.

It was suspected that geographical proximity would be a primary driver of the
propensity for a match to occur, and therefore of model distance. Figure 2
plots the pairwise Euclidean distances, calculated using the latent space
positions, against the pairwise estimated travel times. This shows a strong
relationship, with a Pearson correlation of 74%. Travel time is used here as
this would seem to be a more relevant motivating condition for a fixture than
geographical distance, but substantially similar results are found when using
geographical distance. The comparison between latent space distance and travel
time may also be used as a means of testing the choice to use a two
dimensional latent space. When fitting with three dimensions the correlation
increased only by 1%, and goodness of fit measures based on the posterior
predictive distribution of degree and minimal geodesic distance (Krivitsky and
Handcock, , 2008) showed no clear improvement with increased dimension,
strongly suggesting that modelling in two dimensions, with the
representational benefits it brings, is a reasonable choice. As such, all
further models will be applied using a two dimensional latent space.

Figure 2: Scatterplot of latent space distance against travel time. School
pairs who do not play each other during the three seasons are in red, pairs
who play each other at least once are in blue. OLS regression line is shown in
black.

While the fit is reasonable, there seems to be a notable skew to the
residuals. Looking at the residual plot in Figure 3 it can be seen that this
is driven by two groups. The first has high travel time and considerably lower
latent space distance than the linear regression would suggest, the second low
travel time but considerably higher latent space distance than the regression
would suggest.

Figure 3: Scatterplot of residuals from OLS fit of latent space distance
against travel time

It might be supposed that the former group is likely to be due to the
requirement for geographically extreme teams to travel longer distances in
order to complete a sufficient number of matches in the tournament, so that
travel time for them has a different level of consideration than for teams
with greater geographical proximity to other schools. However Figure 4
suggests this does not account for the entire effect. While the most northerly
and southerly teams, as well as the single team in East Anglia do stand out in
the figure on the left, teams from all over the country are represented and
the third highest weighted is King Edward’s, Birmingham in the middle of the
country, suggesting there is something else at work here. The other group is
made up of geographically proximate teams with a large latent space distance.
Again it is perhaps to be expected that many of these are in the more densely
geographically clustered south east, but some of the most westerly teams also
feature strongly. Even in the case of the south east teams, it remains
unexplained as to why these particular teams have this greater latent space
distance. This will be examined further in Section 4 when the influence of the
other nodal covariates is investigated.

(a)  (b)

Figure 4: Location maps for the teams featuring in the extreme residuals. Size
of dot represents the number of times that the team appears in the relevant
set of pairs. Left hand chart includes pairs where $d<\hat{d}-4$, and right
hand where $d>\hat{d}+4$, where $d$ is the latent space distance and $\hat{d}$
the expected latent space distance based on the linear regression with travel
time.

The geographical implications of model distance may also be considered by
investigating community detection in relation to the geographical location of
the schools. Numerous methods of clustering could be applied given the latent
space distance matrix. As an example, by applying a hierarchical clustering
with complete linkage the dendrogram presented in Figure 5 is obtained.

Figure 5: Dendrogram of hierarchical clustering of schools by latent space

Inspection of the dendrogram suggests that five communities could be an
appropriate partition with a number of communities converging at a distance of
around 45. In Figure 6 communities are plotted for $G=2,3,4,5,6,7$, where $G$
is the number of groups, in order to show the geographical detection ability
with different numbers of groups.

(a)  (b)  (c)

(d)  (e)  (f)

Figure 6: Communities of size G=2,3,4,5,6,7 based on complete linkage
hierarchical clustering of latent space distances

###  3.3 Latent Position Cluster Model

Handcock et al.,  (2007) provided a model with the extended feature that the
latent positions are drawn from a finite mixture of multivariate normal
distributions. That is,

$\boldsymbol{Z_{i}}\overset{i.i.d.}{\sim}\sum_{g=1}^{G}\lambda_{g}\text{MVN}_{d}(\mu_{g},{\sigma_{g}}^{2}I_{d}).$

Two methods of estimation were proposed in Handcock et al.,  (2007). The first
method is a two-stage maximum likelihood procedure, and the second a fully
Bayesian estimation using Markov chain Monte Carlo sampling. In this paper,
this is fitted based on a burn-in period of 10,000 iterations and a sample run
of 1,000,000 iterations of which every fiftieth was sampled, giving a sample
size of 20,000.

Figures 7 and 8 present the results from the two methods. In Figures 6, 7, and
8, communities are identified such that the number of schools remaining in the
same community as in the previous clustering (as represented by a particular
colour) is maximised. The different colours representing different regions in
the three Figures is thus a result of different evolutions of the community
detection in each case. The fittings using the Gaussian clustering, both based
on the two step MLE and the MCMC, appear to show better geographic separation
than did the hierarchical clustering. The MCMC clustering arguably shows a
better separation in the London area, though there is substantial agreement in
the community membership up to $G=5$.

(a)  (b)  (c)

(d)  (e)  (f)

Figure 7: Communities of size $n=2,3,4,5,6,7$ based on two stage MLE

(a)  (b)  (c)

(d)  (e)  (f)

Figure 8: Communities of size n=2,3,4,5,6,7 based on MCMC parameter estimation

For the purposes of further investigation it is useful to select a single
number of communities, $G$, with which to work. Graphical inspection of
Figures 7 and 8 suggests that between three and five communities may fit best,
given the clear geographical separation they evidence. Handcock et al.,
(2007) suggest the use of a Bayesian Information Criterion (BIC) for the
purpose of selecting the number of communities when using the MCMC fitting.
Figure 9 shows this BIC for each value of $G$. Experimentation with different
specifications of the algorithm showed the BIC value to be somewhat unstable.
As such this is taken to be indicative rather than definitive and so taking
into account both the BIC and the geographical separation $G=4$ is chosen.

Figure 9: Bayesian Information Criterion for different values of G

In Figure 10 the travel time in minutes is plotted against the latent space
distance. This is shown side by side for the two stage MLE and the MCMC
estimation to aid comparison. The most notable feature is that the effect of
applying the simultaneous Bayesian estimator obtained from the MCMC fit is to
reduce and also compress latent space distances particularly for the large
latent space distances. The interquartile range of the residuals, for example,
is thus more than 25% lower when compared to the MLE fit. The Pearson
correlation is slightly higher at 78% compared to 74% from the MLE fit, and
the intercept is closer to zero. This suggests that the assumption of
community membership usefully constrains the fit. The community memberships
are substantially similar with 108 of the 118 schools belonging to equivalent
groups. However based on the seemingly improved latent space fit, the MCMC
estimation with $G=4$ is the model selected for further analysis.

(a)  (b)

Figure 10: Scatterplot of latent space distance against travel time for two
stage MLE and MCMC fits with $G=4$. School pairs who do not play each other
during the three seasons are in red, pairs who play each other at least once
are in blue. Linear regression line is shown in black.

A feature of the mixture model is that it allows us to view the probability of
each community membership for each node. Figure 11 provides a graphical
interpretation of these probabilistic community memberships. Perhaps the most
notable feature is the difference in delineation of the communities. What
might be referred to as the northern community is most clearly delineated,
with all but two members having greater than 75% probability of being part of
that community. On the other hand, the two central communities are
substantially more indeterminate with the southernmost (that represented in
the chart in pink) having no member with a greater than 75% probability of
being part of the community. The schools in and around London show perhaps
more delineation than one might expect given the short distances. Due to the
geographical proximity of a number of schools, some are obscured in the chart.
This does not alter the interpretation when the distributions of these schools
are known. It could also be noted that the most geographically isolated
school, St Joseph’s Ipswich, has a substantial probability of being in three
of the four different groups.

(a)  Figure 11: Probability of membership of each community.

An alternative way of considering the link between the latent and geographical
spaces is presented in Figure 12. This suggests that the two latent space
dimensions are picking up something like a polar coordinate system with a
radial coordinate eminating from a centre in the south-east (as seen on the
size scale), and an angular coordinate going from the south to the north east
(as seen in the colour scale), an observation which is consistent with the
probabilistic clustering, as well as the evolution of the community detection
seen in Figures 7 and 8.

Figure 12: Map of schools with colour and size representing first and second
latent space coordinate respectively for each school, based on MCMC estimation
of latent position cluster model with $G=4$

##  4 A consideration of covariates

###  4.1 Relative Importance

The analysis so far has shown a clear link between the travel time and the
propensity for the existence of fixtures, but this is to be expected and it is
desirable to consider how it compares to the other covariates of interest.
Grömping,  (2015) provides a thorough overview of various measures of relative
importance, and the factors to consider in selecting one. These are presented
in the context of linear regression. Some of these measures may be extended to
wider families of models, and in particular, Thomas et al.,  (2008) proposed
relative importance measures for logistic regression that are consistent with
some of the features of the measure axiomatically justified by Pratt,  (1987)
for linear regression.

In the present case, the unreliability and the computational expense of the
coefficient estimates obtained when the model is fitted with all covariates
simultaneously are considerable hindrances to applying relative importance
measures directly to the model. This is especially so since a number of the
measures require at least $p!$ computations, where $p$ is the number of
regressor variables, in order to average over the permutations of regressor
variables. One alternative is instead to use the pairwise latent space
distances estimated by the selected model. A linear regression of these
distances against the nodal covariates can then be considered. Doing so
provides an estimate for the relative importance of the covariates in the
determination of pairwise latent space distance. One reason to be cautious of
this approach is that it may be viewed as a measure of relative importance of
the covariates only in so much as latent space distance is a good measure of
the propensity for a fixture to exist. However if the conclusions drawn from
the analysis prove to be corroborated by other approaches then we may have
more confidence about the meaningfulness of the results.

Two specific measures are considered here. Both have the desirable and
intuitive properties of providing: independence from the order of the
regressors in the model; scale invariance; a decomposition of the model
variance ($R^{2}$); and a proper decomposition of the model variance for any
orthogonal regressor subgroups. The first (which we will refer to as ‘Pratt’)
is the one advocated by Pratt,  (1987), which provides an axiomatic
justification for using $b_{k}\rho_{k}$ to assess relative importance, where
$b_{k}$ is the standardised coefficient for the $k$th regressor variable,
$x_{k}$, and $\rho_{k}$ is the correlation of $x_{k}$ with the independent
variable $y$, in this case the latent space distance. The principal objection
to this approach is that it can produce a negative measure, which is not
clearly interpretable. The second (LMG, after the original authors) is the
method originally due to Lindeman et al.,  (1980), which was also
influentially and independently advocated by Kruskal,  (1987). For each
regressor, this takes the marginal increase of explained variance from the
addition of the regressor to the model, and takes the mean of these values
over all regressor order permutations of the model. Taking the notation of
Grömping,  (2015), denote the explained variance and the sequential additional
variance respectively as

$\displaystyle\text{evar}(S)$ $\displaystyle=\text{var}(y)-\text{var}(y\mid
x_{j};j\in S),$ $\displaystyle\text{svar}(M)$ $\displaystyle=\text{evar}(M\cup
S)-\text{evar}(S).$

Then we can define, without loss of generality, the measure for the first
regressor as

$\text{LMG}(1)=\frac{1}{p!}\sum_{\pi}\text{svar}(\\{1\\}\mid S_{1}(\pi)),$

where $\pi$ are the regressor permutations and $S_{1}(\pi)$ the set of
regressors preceding regressor 1 in permutation $\pi$. The model has an
$R^{2}$ of 62%, and the proportion of that coming from each covariate is shown
in Table 2.

Variable | Pratt | LMG
---|---|---
Travel time | 98.09% | 97.77%
Percent Boarder | 0.88% | 0.88%
Fees | 0.72% | 0.63%
Percent Boys | 0.09% | 0.08%
Boys | 0.08% | 0.12%
School type | 0.05% | 0.30%
Rating | 0.04% | 0.14%
Term type | 0.02% | 0.02%
Founded | 0.01% | 0.04%
Table 2: Proportion of explained variance ($R^{2}$) attributable to each
covariate in linear regression of pairwise latent space distance against
covariates. Ordered in descending order of relative importance under the Pratt
measure.

It is perhaps unsurprising based on the previous analysis that travel time is
dominant, but the degree to which this is the case is nevertheless notable.
The other covariates have negligible relative importance in comparison to
travel time, but relatively they suggest that fees and the proportion of
boarders may have greater influence. The two measures of relative importance
assessed here, Pratt and LMG, are very substantially consistent, with Pratt
giving slightly higher values for travel time and fees and lower for school
type and rating.

###  4.2 Covariate inclusion

An alternative means to investigate the influence of the different covariates
is by considering models with each covariate included individually,

$\text{logit}(\mu_{ij})=\beta_{0}+x_{ijk}\beta_{k}-d(\boldsymbol{Z_{i}},\boldsymbol{Z_{j}}),\quad\quad(k=1,\dots
p),$

where the $k$ represent the nodal covariates. Table 3 presents the parameter
estimates based on the posterior mean of the coefficient, as well as a
statistic representing $q=2\times\text{min}(P(\beta_{k}>0),P(\beta_{k}<0))$,
where $\beta_{k}$ is the coefficient under consideration, and the Bayesian
Information Criterion (BIC) for each model, which may be compared with a BIC
of 3647 for the model with no covariates included. Comparing individually in
this way thus allows us to consider the importance of a particular covariate
based on the change in BIC, and the $q$-value statistic.

|  $\hat{\beta_{k}}$ |  $q$ |  BIC
---|---|---|---
Travel Time |  $-$1.523 |  $<10^{-15}$ |  3316
Fees |  $-$0.368 |  $<10^{-15}$ |  3618
Percent boarder |  $-$0.860 |  $<10^{-15}$ |  3630
Term type |  $-$0.309 |  0.0016 |  3642
6th Form boys |  $-$0.199 |  0.0039 |  3642
School type |  $-$0.370 |  0.0017 |  3643
Founded |  $-$0.038 |  0.0473 |  3646
Percent boys |  $-$0.756 |  0.0070 |  3649
Rating |  $-$0.057 |  0.3239 |  3651
Table 3: Coefficient, $q$-value and BIC when model fitted with individual
additional covariates. Ordered in increasing BIC.

The directions of the coefficient estimates are as one would expect, so that a
smaller covariate distance implies a higher probability of a fixture for all
covariates. The BIC estimates are consistent with what was seen in the
relative importance analysis in highlighting travel time as the dominant
covariate. In that fees and the proportion of boarders are the only others
that have a notably lower BIC and also $q$-value, these are also somewhat
consistent with the previous results. With the exception of rating, the
$q$-values are all low, indicating that it is likely that these factors are
influencing the propensity for a fixture to exist. This is perhaps surprising
in the case of year of foundation, though the $q$-value there is less
conclusive.

Alternatively we may choose to fit the model with all covariates, as in
equation 4.2.

$\text{logit}(\mu_{ij})=\beta_{0}+\sum^{p}_{k=1}x_{ijk}\beta_{k}-d(\boldsymbol{Z_{i}},\boldsymbol{Z_{j}}),$

The results are presented in Table 4. The computation of this shows greater
sensitivity to the specification of the algorithm and even the order in which
covariates are presented, so we might be cautious about the results. They do
however show broad agreement with previous analyses, and where there are
discrepancies between the coefficient estimates when taken together and
individually these may reasonably be thought of as being as a result of the
confounding of covariates. For example, the correlation between the percentage
of boarders and fees provides an explanation for the smaller absolute impact
from boarding when the covariates are considered together.

In this context however it is perhaps more intuitive to interpret the impact
of the covariates individually as in the models represented by equation 4.2,
in the sense that knowing nothing else about the schools in each case one
would expect the odds of a fixture to decrease by 78% for every additional
hour of travel between schools, by 31% for every £10,000 per annum difference
in fees, and by 58% for a fully boarding school (100% boarders) playing a
fully day school (0% boarders) as compared to a match between two schools with
the same proportion of boarders.

|  $\hat{\beta_{k}}$ |  $q$
---|---|---
Travel Time |  $-$1.552 |  $<10^{-15}$
Fees |  $-$0.363 |  $<10^{-15}$
Percent boarder |  $-$0.317 |  $<10^{-15}$
Term type |  $-$0.283 |  $<10^{-15}$
6th Form boys |  $-$0.252 |  $<10^{-15}$
Percent boys |  $-$0.734 |  $<10^{-15}$
Founded |  $-$0.042 |  0.0478
School type |  $+$0.154 |  0.1541
Rating |  $-$0.002 |  0.9613
Table 4: Coefficient $q$-values when model fitted with all covariates. Ordered
in increasing $q$-value.

###  4.3 Graphical inspection

Another method of investigation that the model allows is by inspection of the
latent space graphically. Given the dominance of travel time, it is useful to
control for it in considering the other covariates. Thus the following model
is considered with G=4 based on BIC values.

$\text{logit}(\mu_{ij})=\beta_{0}+x_{ij\text{Travel Time}}\beta_{\text{Travel
Time}}-d(\boldsymbol{Z_{i}},\boldsymbol{Z_{j}}),$

The colours are used to show the fees and proportion of boarders for each
school. As was noted in Section 2 these covariates are strongly linked and
this is evident in the plots. It can be seen that while schools with lower
fees and low proportions of boarders are distributed quite evenly throughout
the latent space, those with higher fees and a high proportion of boarders are
to be found disproportionately often in a cluster at the bottom of the latent
space. This effect is perhaps clearer for the proportion of boarders. This
suggests that there is a small community of schools with a higher proportion
of boarders and higher fees, who are more likely to play each other, and that
it is this that was driving the greater relative importance of these
covariates noted in sections 4.1 and 4.2.

(a)  (b)

Figure 13: Latent space plots based on model with Travel Time controlled for
and with G=4. Colour scale representing Fees on left hand side, and % Boarders
on right hand side

##  5 Concluding Remarks

The latent space models demonstrated a number of desirable features in the
context of addressing the question of the nature of institutional ties. They
showed a good ability to detect geographical communities in the data set, and
the spatial nature of the output provided a natural means of interpreting the
result in terms of pairwise distances. It was possible to employ a variety of
methods in order to analyse the relative importance of the nodal covariates in
the propensity for a fixture to exist — by using the methods of Pratt,  (1987)
and Lindeman et al.,  (1980) to provide quantitative estimates based on latent
space distance; by inspection of
$q=2\times\text{min}(P(\beta_{k}>0),P(\beta_{k}<0))$ and BIC on fitting models
with covariates included; and by using graphical assessment based on latent
space plots as in Figure 13. These produced a consistent interpretation. The
latent position cluster model of Handcock et al.,  (2007) was found to
usefully constrain the latent space model such that it appeared to provide a
better fit. It also usefully facilitated a means of assessing the strength of
community attribution at a school level by producing the probabilities of each
school being a member of each community as seen in Figure 11. It is notable
that, while a relatively high number of iterations was chosen to be run based
on the diagnostics, most of these conclusions were substantially unchanged
when running off the default settings of burn-in of 10,000, with 40,000
sampling iterations of which every tenth was sampled giving a sample size of
4,000. The BIC values and coefficient ranges were however notably more
affected by changes in other algorithm specifications when run on a lower
number of simulations. The latent space models thus performed well in
providing quantitative and qualitative insights into the nature of communities
and the influential elements of edge formation in this network.

Perhaps unsurprisingly travel time was found to be the dominant factor in
contributing to the prevalence of fixtures. The degree of this dominance is
notable however, representing 98% of $R^{2}$ based on the two relative
importance measures employed. The proportion of boarders and fees were
identified as the next most important factors, with the dependency driven by a
community of schools with the highest proportion of boarders and the highest
fees. It should be noted however that the data represent a self-selecting
group of the top rugby-playing schools. It seems not unreasonable to expect
that were a full set of data available then the impact of some of the other
covariates would take on a greater importance. For example, over a wider data
set with more state schools and a greater range of abilities, perhaps fees,
school type or rating would take on greater relative importance. However for
the question at hand, it is still notable that rather than clustering seeming
to correlate with, for example, rating, as meritocrats might have liked to
believe, or with school type as others might have suspected, it was a school’s
status as a majority-boarding school, with a high level of fees, that produced
a detectable clustering effect. This might be seen as consistent at an
institution level with a traditional understanding of the old boy network.

## References

  * Cooley, (2017) Cooley, D. (2017).  googleway: Accesses google maps APIs to retrieve data and plot maps.  R package version, 2(0).
  * Department for Education, (2019) Department for Education (2019).  Schools, pupils and their characteristics: January 2019.  https://www.gov.uk/government/statistics/schools-pupils-and-their-characteristics-january-2019.
  * Grömping, (2015) Grömping, U. (2015).  Variable importance in regression models.  Wiley Interdisciplinary Reviews: Computational Statistics, 7(2):137–152.
  * Hamilton and Firth, (2021) Hamilton, I. and Firth, D. (2021).  Retrodictive modelling of modern rugby union: Extension of Bradley-Terry to multiple outcomes.
  * Handcock et al., (2007) Handcock, M. S., Raftery, A. E., and Tantrum, J. M. (2007).  Model-based clustering for social networks.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 170(2):301–354.
  * Hoff et al., (2002) Hoff, P. D., Raftery, A. E., and Handcock, M. S. (2002).  Latent space approaches to social network analysis.  Journal of the American Statistical Association, 97(460):1090–1098.
  * Krivitsky and Handcock, (2008) Krivitsky, P. N. and Handcock, M. S. (2008).  Fitting position latent cluster models for social networks with latentnet.  Journal of Statistical Software, 24.
  * Kruskal, (1987) Kruskal, W. (1987).  Relative importance by averaging over orderings.  The American Statistician, 41(1):6–10.
  * Lindeman et al., (1980) Lindeman, R. H., Merenda, P., and Gold, R. Z. (1980).  Introduction to bivariate and multivariate analysis.  Scott, Foresman and company.
  * Pratt, (1987) Pratt, J. W. (1987).  Dividing the indivisible: using simple symmetry to partition variance explained.  In Pukkila, T. and Puntanen, S., editors, Proceedings of the Second Tampere International Conference in Statistics, pages 245–260. University of Tampere, Finland.
  * Thomas et al., (2008) Thomas, D. R., Zhu, P., Zumbo, B. D., and Dutta, S. (2008).  On measuring the relative importance of explanatory variables in a logistic regression.  Journal of Modern Applied Statistical Methods, 7(1):4.


[a]Sam Offler

# Reconstruction of bottomonium spectral functions in thermal QCD using Kernel
Ridge Regression

Gert Aarts  Chris Allton  Benjamin Jäger  Seyong Kim  Maria Paola Lombardo
Ben Page  Sinead M. Ryan  Jon-Ivar Skullerud  Thomas Spriggs

###### Abstract

We discuss results for bottomonium at nonzero temperature obtained using NRQCD
on Fastsum Generation 2L ensembles, as part of the Fastsum collaboration’s
programme to determine the spectrum of the bottomonium system as a function of
temperature using a variety of approaches. Here we give an update on results
for spectral functions obtained using Kernel Ridge Regression. We pay in
particular attention to the generation of training data and introduce the
notion of using lattice QCD ensembles to learn how to improve the generation
of training data. A practical implementation is given.

##  1 Introduction

Since the initial proposal that suppression of $J/\psi$ could be a signal for
the formation of the quark-gluon plasma (QGP) [1], quarkonia, bound states of
heavy quark-antiquark pairs, have been intensively studied. These states serve
as important probes for the plasma created in heavy-ion collisions as they do
not reach thermal equilibrium with the medium. Initial work focused on the
charmonium system, but more recently the focus of both experimental [2] and
theoretical (see e.g. Refs. [3, 4]) studies has shifted to bottomonium.

The purpose of this contribution is to continue the investigation of
bottomonium at nonzero temperature by the Fastsum collaboration [5, 6, 7],
with the aim of determining masses and widths of ground and excited states,
provided they exist, as well as the state-dependent melting or dissociation
temperatures. In Ref. [8] Kernel Ridge Regression (KRR) was introduced to
reconstruct bottomonium spectral functions from Euclidean lattice correlators.
KRR is a machine learning method that requires training data and here we
revisit intricacies of generating this data. In particular we introduce the
notion of using lattice QCD ensembles to learn how to improve the generation
of mock data and demonstrate how this may be implemented.

##  2 Lattice details

Given that the bottom quark mass scale $M$ lies above any other scale in the
problem, including the temperature $T$ of the quark-gluon plasma, it can be
integrated out to obtain standard non-relativistic QCD (NRQCD) [9, 10], also
at nonzero temperature [11]. We use the same formulation as in Refs. [6, 7],
i.e. terms up to $\mathcal{O}(v^{4})$ are included, where $v$ is the heavy
quark velocity in the bottomonium rest frame. In the relation between the
bottomonium Euclidean correlator and the spectral function, terms involving
$M/T$ are exponentially suppressed, leading to a relation with a simplified
(“$T=0$”) kernel,

$G(\tau)=\int^{\omega_{\rm max}}_{\omega_{\rm
min}}\frac{d\omega}{2\pi}\,K(\tau,\omega)\rho(\omega),\qquad\qquad
K(\tau,\omega)=e^{-\omega\tau},$ (1)

where it is understood that an additive constant $\sim 2M$ has to be added to
the lower limit of the integral, to make contact with phenomenology in full
QCD.

$1/a_{\tau}$ [GeV] |  $a_{s}$ [fm] | $\xi=a_{s}/a_{\tau}$ | $N_{s}$ |  $m_{\pi}$ [MeV] | $m_{\pi}L$ |  $T_{\rm pc}$ [MeV]
---|---|---|---|---|---|---
$5.997(34)$ | $0.01136(6)$ | $3.453(6)$ | $32$ | $236(2)$ | $4.36$ | $164(2)$
Table 1:  Parameters relevant for the ensembles: $a_{s}$ ($a_{\tau}$) is the
spatial (temporal) lattice spacing; $\xi$ is the renormalised anisotropy;
$N_{s}$ is the number of points in the spatial direction; $m_{\pi}$ is the
mass of the pion; $T_{\rm pc}$ is the pseudocritical temperature, obtained
from the inflection point of the renormalised chiral condensate [12].

For the finite-temperature study, we use the anisotropic Fastsum Generation 2L
ensembles, with $N_{f}=2+1$ flavours of clover-improved Wilson fermions, see
Ref. [12]. The light quarks are heavier than in nature, but the strange quark
mass is at its physical value. Relevant parameters are given in Table 1. In
the fixed-scale approach, the temperature is varied by changing $N_{\tau}$,
using the relation $T=1/(N_{\tau}a_{\tau})$. The $N_{\tau}$ values and
corresponding temperatures are shown in Table 2.

$N_{\tau}$ | 128 | 64 | 56 | 48 | 40 | 36 | 32 | 28 | 24 | 20 | 16
---|---|---|---|---|---|---|---|---|---|---|---
T [MeV] | 47 | 94 | 107 | 125 | 150 | 167 | 187 | 214 | 250 | 300 | 375
Table 2: Temporal lattice sizes and temperatures for the fastsum Generation 2L
ensembles [12].

##  3 Data generation

As is the case with many machine learning methods, Kernel Ridge Regression
(KRR) requires training in order to make predictions. In the current context,
this means it is necessary to generate a set of mock spectral functions and
calculate the corresponding correlators according to Eq. (1), yielding the
training set $\\{\rho_{i}(\omega),G_{i}(\tau)\\}$ ($i=1,\ldots,N_{\rm
train}$). The importance of generating appropriate training data cannot be
understated; it is only at this stage that we can incorporate any
constraints/information we know of.

In this study we choose to construct mock spectral functions from the
combination of $N_{p}=5$ Gaussian peaks, such that a single spectral function
is written as

$\rho(\omega)=\sum^{N_{p}}_{p=1}Z_{p}\exp\left(-\frac{(\omega-
m_{p})^{2}}{\Gamma_{p}}\right).$ (2)

Here $Z_{p}$ is the amplitude, $m_{p}$ is the position and $\Gamma_{p}$ is
related to the width of the Gaussian peak. Other mock functions can be
constructed; in Ref. [8] the logarithm of the spectral function was expanded
in an orthonormal set of incomplete basis functions. Continuing with the
Gaussian peaks here, we generate a collection of peaks by sampling values for
$\\{Z_{p},m_{p},\Gamma_{p}\\}$ from a set of distributions. Masses are chosen
from the following exponential distribution,

$m_{\mathrm{GeV}}=9+\frac{1}{\zeta}e^{-\omega/\zeta},$ (3)

where $\zeta$ determines the decay rate of the distribution. These values are
then converted to lattice units using

$m_{\mathrm{lat}}=a_{\tau}\left(m_{\mathrm{GeV}}-\Delta M\right),$ (4)

where $\Delta M=7.465$ GeV is the additive constant mentioned below Eq. (1).
As in Refs. [6, 7], it is determined by equating the calculated zero-
temperature ground state mass in the $\Upsilon$ channel to its experimental
value. Widths are determined by writing $\Gamma=10^{-x}$ and selecting $x$
from a uniform distribution, $0<x<6$. The amplitudes were selected for a
uniform distribution, $0<Z<10$. We put $\rho(\omega)=0$ for
$\omega<\omega_{\rm min}$ and $\omega>\omega_{\rm max}$. Finally, the mock
spectral functions (2) are normalised, using the relation at $\tau=0$,

$G(0)=\int^{\omega_{\rm max}}_{\omega_{\rm
min}}\frac{d\omega}{2\pi}\,K(0,\omega)\rho(\omega)=\int^{\omega_{\rm
max}}_{\omega_{\rm min}}\frac{d\omega}{2\pi}\,\rho(\omega),$ (5)

where $G(0)$ is the source in the NRQCD formulation. In total we generated a
set of 20 000 mock spectral functions, each consisting of $N_{p}=5$ peaks and
characterised by $3N_{p}=15$ parameters. Only $N_{\rm train}=15\,000$
functions from this dataset were used for training the KRR model.

Figure 1: Distribution of the natural logarithm of correlators in the training
set at each timeslice for $N_{\tau}=128$ for S-wave (left) and P-wave (right)
states. Dashed lines represent the upper and lower limits of the correlators
in the training set for each value of $n_{\tau}$. Solid lines represent the
Euclidean correlators obtained in simulations. Note that the two channels in
each figure – $\eta_{b}$, $\Upsilon$ (left) and $\chi_{b_{0}}$, $\chi_{b_{1}}$
(right) – are hardly distinguishable.

To assess whether the training set is potentially relevant for the correlators
obtained from the actual lattice QCD simulations, we carry out the following
comparison. At each timeslice we determine the distribution of (the logarithm
of) the correlators in the training set, as obtained from the mock spectral
functions. Given that the mock set contains 20 000 spectral functions, the
density quickly falls below 10%. This distribution is shown in Fig. 1 via a
heat map. The dashed lines represent the upper and lower limits of the
correlators in the training set for each value of $n_{\tau}$. In addition we
show the actual correlators obtained in the simulations. Both the left and the
right panel in Fig. 1 show two correlators – $\eta_{b}$, $\Upsilon$ (left) and
$\chi_{b_{0}}$, $\chi_{b_{1}}$ (right) –, but note that these are hardly
distinguishable.

The key observation is that the training data compares poorly to the actual
correlators, i.e. there is little overlap between the most common values seen
in the training set and the real data. This effect is worse in the P-wave
channels. Since machine learning algorithms assume that training and
validation data sets are drawn from the same ensembles as (or at least are
representative of) the real data sets, this is not a desired feature. We
remark that the spread of the training data itself is not a problem; it is
beneficial for KRR to have access to outliers in principle.

In order to generate a better set of training data, we use the actual
correlators to learn which members of the training set are the most
representative. This is implemented as follows: from the original training
set, a subset of correlators are chosen by determining which correlator
$G_{i}(\tau)$ lies closest to the actual correlator $G_{\rm NRQCD}(\tau)$, at
each value of the temperature and for S-waves ($\eta_{b}$, $\Upsilon$) and
P-waves ($\chi_{b_{0}}$, $\chi_{b_{1}}$) separately. Closeness is here defined
via

$i=\mathrm{argmin}\sum^{N_{\tau}-1}_{n_{\tau}=4}\left|\mathrm{ln}\left[\frac{G_{i}(n_{\tau})}{G_{\rm{NRQCD}}(n_{\tau})}\right]\right|.$
(6)

Note that the sum starts at $n_{\tau}=4$; the first few timeslices are not
included in the summation. Given that there are ten different temperatures,
see Table 2, and two S-wave or P-waves channels, this procedure yields a
subset of up to 20 representative correlators (note that repeats are
possible).

This subset is used to construct a larger, more representative, training set,
as follows. For each correlator in the subset, the corresponding spectral
function and, more importantly, the set of $3N_{p}$ parameters
{$Z_{p},m_{p},\Gamma_{p}$} are known. From those values the mean and variance
is computed and new parameter values are sampled from Gaussian distributions
with these means and variances. Note that for the amplitude $Z$ and the
parameter $\Gamma$ this is done indirectly, using $\Gamma=10^{-x}$ and an
equivalent relation for $Z$, as above. From this point onwards the procedure
is the same as earlier, but with different distributions for the parameters
characterising the mock spectral functions.

Figure 2: As above, for the improved training set.

Fig. 2 displays the same comparison as in Fig. 1, but now for the improved
training set. As can be seen, the most common correlators in the training set
now coincide with the actual correlators, leading potentially to a more
appropriate training set. This is further demonstrated in Fig. 3, where the
distribution of correlators obtained from the original and the improved
training sets are shown at selected timeslices for the $N_{\tau}=128$ lattice.
The actual correlator values at the selected timeslices are shown with the
vertical lines. We indeed observe a better overlap between the actual
correlators and the data in the training set. The latter still has a
nonnegligible width, allowing for flexibility in the application to spectral
reconstruction.

Figure 3: Distribution of the original and the improved training set for the
correlators at $n_{\tau}$ = 4 (left), 40 (centre), and 127 (right) for
$N_{\tau}=128$. Solid vertical lines represent the value of the actual
Euclidean correlators at the specified timeslices.

##  4 Kernel Ridge Regression

Kernel Ridge Regression (see e.g. Ref. [13]) combines two techniques: the
kernel method and ridge regression. Consider first linear regression in which
$y=w^{T}\phi(x)$, where $y$ are the target data, $\phi(x)$ is a vector of
functions of the input data $x$, and $w$ is a vector of parameters to be
determined. The cost function to be minimized is
$E=\frac{1}{2}[y-w^{T}\phi(x)]^{2}$. In our application, the input data are
the training set of Euclidean correlators $G_{i}(\tau)$ ($i=1,\ldots,N_{\rm
train}$). However, rather than using these directly, they are used to generate
an $N_{\rm train}\times N_{\rm train}$ matrix $\mathbf{C}$ (kernel function),
with matrix elements

$C_{ij}=\exp\left(-\gamma\sum^{N_{\tau}-1}_{n_{\tau}=4}\left[\widetilde{G}_{i}(n_{\tau})-\widetilde{G}_{j}(n_{\tau})\right]^{2}\right),\qquad\qquad\widetilde{G}_{i}(n_{\tau})=\frac{G_{i}(n_{\tau})}{\overline{G}(n_{\tau})}.$
(7)

Note that the form of the kernel function is not unique. Here $\gamma$ is a
hyper-parameter, which sets a correlation length in the space of correlators.
The correlator data appearing is normalised with
$\overline{G}(n_{\tau})=\sum_{i}G_{i}(n_{\tau})/N_{\rm train}$, the mean of
the training data at each timeslice. This is done to account for the
difference in absolute size of the correlators as $\tau$ changes.

The target data in our application are the spectral functions. To control the
size of the target data, we use here expansion (2) in terms of Gaussian peaks,
leading to $N_{a}\equiv 3N_{p}=15$ parameters for spectral function. The
target data is then represented by a $N_{\rm train}\times N_{a}$ matrix
$\mathbf{Y}$, in which each row contains the parameters of a single spectral
function. It is now assumed that the input and target data are related as

$\mathbf{Y}=\mathbf{C}\mathbf{\alpha},$ (8)

in which the $N_{\rm train}\times N_{a}$ matrix $\mathbf{\alpha}$ is the
analogue of the vector $w$ in linear regression. The aim of the training stage
is to determine the matrix $\mathbf{\alpha}$.

To prevent overfitting, it is common to add an additional term in the cost
function, proportional to the square of the parameters (ridge regression).
Hence the cost function to be minimised reads

$E(\mathbf{Y,C,\alpha})=\frac{1}{2}(\mathbf{Y}-\mathbf{C}\mathbf{\alpha})^{2}+\frac{1}{2}\lambda\mathbf{\alpha}^{T}\mathbf{C}\mathbf{\alpha},$
(9)

where $\lambda$ is the second hyper-parameter, used to regularise the
influence of the additional term. Minimising this cost function with respect
to $\mathbf{\alpha}$ then determines the optimal parameter matrix, for given
$\sigma$ and $\lambda$, as

$\mathbf{\alpha}_{\rm
opt}=\left(\mathbf{C}+\lambda\mathbf{I}\right)^{-1}\mathbf{Y}.$ (10)

The hyper-parameters are determined via a cross-validation procedure, see
below.

After this training stage, it is possible to make predictions for a spectral
function ($\mathbf{Y^{\prime}}$) given an actual Euclidean correlator using

$\mathbf{Y^{\prime}}=\mathbf{C^{\prime}}\mathbf{\alpha}_{\rm opt},$ (11)

where $\mathbf{C^{\prime}}$ is determined from the squared rescaled difference
between the actual correlator and the training correlators, i.e. it is a
matrix of size $1\times N_{\rm train}$.

##  5 Optimization

Figure 4: Heatmap of $\overline{R}^{2}$ used for cross-validation in the plane
of hyper-parameters ($\lambda,\gamma$), for $N_{\tau}=20$ (left) and
$N_{\tau}=128$ (right). The maxima, indicated with the little squares, provide
the optimal choice of hyper-parameters. In cases where $\overline{R}^{2}$ <0,
it has manually been set equal to 0 for clarity.

To determine the optimal choice of hyper-parameters, we consider the so-called
$R^{2}$ score, the coefficient of determination, defined in general as [14]

$R^{2}(y)=1-\frac{\sum_{i}(y^{\rm true}_{i}-y^{\rm
pred}_{i})^{2}}{\sum_{i}(y^{\rm true}_{i}-\overline{y}^{\rm true})^{2}}.$ (12)

In our application, $y$ is one of the $N_{a}=15$ parameters parametrising a
spectral function and we consider the mean $R^{2}$ score, averaged over the
parameters,

$\overline{R}^{2}=\frac{1}{N_{a}}\sum_{a=1}^{N_{a}}R^{2}(y_{a}).$ (13)

Rather than using a training and a validation set, we use cross-validation
[13] in which a single set is used for both training and testing. The
optimisation procedure seeks to find the value of $\overline{R}^{2}$ closest
to 1. Results for $\overline{R}^{2}$ are shown in Fig. 4, for all the hyper-
parameters pairs tested, and for two ensembles, $N_{\tau}=20$ (left) and
$N_{\tau}=128$ (right). Though the plots are similar, there does appear to be
some thermal dependence. In both cases there is a broad region of hyper-
parameter pairs that perform to a similar standard. The most optimal choice of
hyper-parameters is shown by the little squares, with
$\left(\ln(\gamma),\ln(\lambda)\right)$ = (-1.25, -3.25) (left) and
$\left(\ln(\gamma),\ln(\lambda)\right)$ = (-1.5, -2.25) (right). It is noted
that the optimal values of $\overline{R}^{2}$ are rather low, i.e. not close
to 1. Looking at the $R^{2}$ scores for each of the $N_{a}$ parameters
individually, we note that the $R^{2}$ score for the ground state mass is
$R^{2}_{m_{0}}=0.869$, whereas the $R^{2}$ scores for the groundstate width
and amplitude are significantly lower, $0.148$ and $0.336$ respectively. We
may therefore further improve the procedure by treating the various parameters
on unequal footing.

Figure 5: Spectral functions in the $\Upsilon$ channel at 4 different
temperatures, obtained using Kernel Ridge Regression (KRR). The left pane
zooms in on the ground state peak.

As a preliminary result, Fig. 5 finally presents the predicted spectral
functions in the $\Upsilon$ channel at four different temperatures. It is
noted that the masses of both the ground and first excited state are smaller
than the expected values. This may be due to the choice of expressing the
spectral functions in terms of $N_{p}=5$ Gaussian peaks and requires further
investigation. Returning to the representation described in [8] may improve
this. For the application of alternative methods to the same data, see Refs.
[15, 16], also presented at this conference.

##  6 Summary

The main message of this contribution is that it is possible to adapt the
training set by comparing its predictions to the actual correlator data from
lattice QCD simulations, at the level of the correlators. A large discrepancy
would suggest that the training sets are not representative of the real data
sets, an undesired feature. It is noted that improving the overlap is
necessary, but not sufficient. As a general remark, we note that this method
can be iterated in principle and is applicable to other parametrisations of
training spectral functions as well. The parametrisation used here – spectral
functions are written as sums of Gaussian peaks – may well be an important
limiting factor.

## Acknowledgments

We thank Jonas Glesaaen for collaboration at the early stages of this project.
This work is supported by STFC grant ST/T000813/1. SK is supported by the
National Research Foundation of Korea under grant NRF-2021R1A2C1092701 funded
by the Korean government (MEST). BP has been supported by a Swansea University
Research Excellence Scholarship (SURES). This work used the DiRAC Extreme
Scaling service at the University of Edinburgh, operated by the Edinburgh
Parallel Computing Centre on behalf of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). This equipment was funded by BEIS capital funding via STFC
capital grant ST/R00238X/1 and STFC DiRAC Operations grant ST/R001006/1. DiRAC
is part of the National e-Infrastructure. This work was performed using PRACE
resources at Cineca via grants 2015133079 and 2018194714. We acknowledge the
support of the Supercomputing Wales project, which is part-funded by the
European Regional Development Fund (ERDF) via Welsh Government, and the
University of Southern Denmark for use of computing facilities. We are
grateful to the Hadron Spectrum Collaboration for the use of their zero
temperature ensemble.

## References

  * [1] T. Matsui and H. Satz, Phys. Lett. B 178 (1986), 416-422.
  * [2] S. Chatrchyan et al. [CMS], Phys. Rev. Lett. 109 (2012), 222301 [arXiv:1208.2826 [nucl-ex]].
  * [3] R. Rapp, D. Blaschke and P. Crochet, Prog. Part. Nucl. Phys. 65 (2010), 209-266 [arXiv:0807.2470 [hep-ph]].
  * [4] R. Larsen, S. Meinel, S. Mukherjee and P. Petreczky, Phys. Lett. B 800 (2020), 135119 [arXiv:1910.07374 [hep-lat]].
  * [5] G. Aarts, S. Kim, M. P. Lombardo, M. B. Oktay, S. M. Ryan, D. K. Sinclair and J. I. Skullerud, Phys. Rev. Lett. 106 (2011), 061602 [arXiv:1010.3725 [hep-lat]].
  * [6] G. Aarts, C. Allton, S. Kim, M. P. Lombardo, M. B. Oktay, S. M. Ryan, D. K. Sinclair and J. I. Skullerud, JHEP 11 (2011), 103 [arXiv:1109.4496 [hep-lat]].
  * [7] G. Aarts, C. Allton, T. Harris, S. Kim, M. P. Lombardo, S. M. Ryan and J. I. Skullerud, JHEP 07 (2014), 097 [arXiv:1402.6210 [hep-lat]].
  * [8] S. Offler, G. Aarts, C. Allton, J. Glesaaen, B. Jäger, S. Kim, M. P. Lombardo, S. M. Ryan and J. I. Skullerud, PoS LATTICE2019 (2019), 076 [arXiv:1912.12900 [hep-lat]].
  * [9] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea and K. Hornbostel, Phys. Rev. D 46 (1992), 4052-4067 [arXiv:hep-lat/9205007 [hep-lat]].
  * [10] G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51 (1995), 1125-1171 [arXiv:hep-ph/9407339 [hep-ph]].
  * [11] Y. Burnier, M. Laine and M. Vepsalainen, JHEP 01 (2008), 043 [arXiv:0711.1743 [hep-ph]].
  * [12] G. Aarts, C. Allton, J. Glesaaen, S. Hands, B. Jäger, S. Kim, M. P. Lombardo, A. A. Nikolaev, S. M. Ryan and J. I. Skullerud, et al. [arXiv:2007.04188 [hep-lat]].
  * [13] C.M. Bishop, Pattern Recognition and Machine Learning, Springer-Verlag New York Inc.
  * [14] $R^{2}$ score, scikit-learn.org/stable/modules/generated/sklearn.metrics.r2_score.html
  * [15] T. Spriggs et al., PoS LATTICE2021 077\.
  * [16] B. Page et al., PoS LATTICE2021 134\.


Oli Danyi Liu, Hao Tang, Sharon Goldwater

# Self-supervised Predictive Coding Models Encode Speaker and Phonetic
Information in Orthogonal Subspaces

###### Abstract

Self-supervised speech representations are known to encode both speaker and
phonetic information, but how they are distributed in the high-dimensional
space remains largely unexplored. We hypothesize that they are encoded in
orthogonal subspaces, a property that lends itself to simple disentanglement.
Applying principal component analysis to representations of two predictive
coding models, we identify two subspaces that capture speaker and phonetic
variances, and confirm that they are nearly orthogonal. Based on this
property, we propose a new speaker normalization method which collapses the
subspace that encodes speaker information, without requiring transcriptions.
Probing experiments show that our method effectively eliminates speaker
information and outperforms a previous baseline in phone discrimination tasks.
Moreover, the approach generalizes and can be used to remove information of
unseen speakers.111This work was supported by the UKRI Centre for Doctoral
Training in Natural Language Processing, funded by the UKRI (grant
EP/S022481/1) and the University of Edinburgh, School of Informatics and
School of Philosophy, Psychology & Language Sciences.

Index Terms: self-supervised learning, unsupervised speech processing, speaker
normalization

##  1 Introduction

Self-supervised learning (SSL) models of speech are trained on large
quantities of unlabeled data to learn useful representations of speech audio.
Empirically, SSL pre-training has been shown to improve performance on
downstream tasks while reducing reliance on annotated data [1, 2]. Researchers
have also started to evaluate SSL models as computational models of phonetic
acquisition and speech perception in humans [3, 4].

It is natural to ask what makes these learned representations so effective. So
far, most work in the area has done so by analyzing what type of information
is encoded, finding that SSL models encode information ranging from acoustic
[5] and contextual [6] cues, to gender [7] and speaker identity [8]. Some
researchers have also analyzed where in these models (at which layer) such
information is encoded [9]. In this paper, we ask a question that has received
little attention so far, namely how different types of information are
distributed across the dimensions of the representation space. Analysis
targeted at this question can potentially improve the interpretability of the
learned representations, inspire new methods for disentangling linguistic and
non-linguistic information, and shed light on properties of the representation
space that can be compared against neural representations or behaviour data of
humans to evaluate the scientific value of the model. We know of no previous
work analyzing this question in detail, although some studies have shown
through two-dimensional visualizations that representations may cluster
according to their corresponding language, gender, or phone class [7]; word-
level context [6]; or speaker identity [8].

In this work, we explicitly investigate how speaker and phonetic information
are distributed in the representation space learned by SSL models. We
hypothesize that a good representation (one that is efficient and works well
for predicting speech) should implicitly disentangle these two sources of
information, since they vary independently in the processes that generate the
speech signal. If so, then the two types of information would be represented
in orthogonal subspaces of the representation space, and it would be possible
to perform speaker normalization by identifying the speaker subspace and
collapsing it.

We test our hypothesis in experiments with two SSL models (APC [10] and CPC
[11]) trained on English LibriSpeech. We use principal component analysis
(PCA) to identify the subspaces that capture most of the variance for speakers
and for phones, and confirm that these subspaces are nearly orthogonal. We
then test whether the speaker subspace generalizes to unseen speakers by
projecting speakers from a test set onto the subspace orthogonal to the one
learned on the training set. We use probing classifiers and an ABX phone
discrimination task to show that with the resulting representations, (a)
little speaker information remains; and (b) within- and across-speaker phone
discrimination improves, outperforming a baseline of utterance-level
standardization [8]. Our results suggest that these two SSL models implicitly
disentangle speaker and phone information into orthogonal subspaces.

##  2 Overview of the approach

We demonstrate the orthogonal speaker and phone subspaces in the SSL
representations using a simple approach, where we first aggregate the
representations by speaker and phone, then perform PCA on the resulting
matrices222Performing PCA on the raw representations is far more compute-
intensive and did not work; the aggregation step seems necessary to overcome
noise. and compare the principal directions. We then show that collapsing the
speaker subspace is an effective speaker normalization technique, for both
seen and unseen speakers. These steps are explained in more detail below.

Aggregating representations by speaker and phone Suppose that the dataset for
analysis contains a set of speakers $S$ and has time-aligned phone
transcriptions (which are used for analysis only, and are not required for the
speaker normalization step). We denote $Z_{s}$ to be the set of frame-level
hidden vectors (each of dimension $D$) from an SSL model of a speaker $s\in
S$. Similarly, we define $Z_{p}$ to be the set of hidden vectors labeled as
phone $p\in P$, where $P$ is a phone set; we define $Z_{s,p}$ to be the set of
hidden vectors labeled as phone $p\in P$ and of a speaker $s\in S$.

For our analysis, we aggregate the representations in three ways: (1) by
speaker (2) by phone, and (3) by each combination of speaker and phone. This
results in three matrices: The first, $M_{\text{spk}}$, is a $|S|\times D$
matrix of speakers, where the $s$-th row $M_{\text{spk}}[s]$ is the average of
all frames of a speaker $s\in S$, i.e., $M_{\text{spk}}[s]=\text{avg}(Z_{s})$.
Analogously, $M_{\text{phn}}$ is a $|P|\times D$ matrix of phones, where the
$p$-th row $M_{\text{phn}}[p]$ is the average of all frames labeled as $p\in
P$, i.e., $M_{\text{phn}}[p]=\text{avg}(Z_{p})$. Finally, $M_{\text{joint}}$
is a $|S||P|\times D$ matrix, where each row corresponds to
$\text{avg}(Z_{s,p})$ for a particular speaker $s$ and a particular phone $p$.

Identifying speaker and phone dimensions Applying PCA to $M_{\text{spk}}$
gives us $|S|$ principal components that describe the directions that account
for the largest variance among the speakers. A speaker subspace is the
subspace spanned by the eigenvectors corresponding to the top eigenvalues.
Similarly, PCA on $M_{\text{phn}}$ and $M_{\text{joint}}$ gives us a phone
subspace and a joint speaker-phone subspace, respectively.

In section 4, we analyze the relationships between these matrices by measuring
the similarity between their principal directions. We also look at the
projection of $M_{\text{joint}}$ on each of its principal directions to
explore what specific information is encoded in each dimension.

Collapsing the speaker subspace If speaker and phonetic information are
represented in orthogonal subspaces, projecting frames to the subspace
orthogonal to the speaker subspace would remove speaker information without
affecting how well phonetic information can be extracted. Specifically, for a
hidden vector $z$ and a principal direction $v$, the vector
$z^{\prime}=z-(z^{\top}v)v$ is orthogonal to $v$. This step can be performed
on multiple principal directions, and we refer to the projection as collapsing
the speaker space when $v$ is a principal component of $M_{\text{spk}}$.
Before doing this, we need to choose the number of principal components to
use, which can be done based on the cumulative variance explained by the
principal components for the training set, or by tuning development set
performance (e.g., on phone discrimination). To test whether a speaker
subspace generalizes across different sets of speakers, we use a speaker
subspace learned from a different set of speakers for collapsing.

##  3 Experimental setup

The two models used in this work are the autoregressive predictive coding
(APC) model [10] and the contrastive predictive coding (CPC) model [11]. Both
CPC and APC are trained to predict one or several future frame(s) given past
context in the same utterance, with the main difference being that CPC is
optimized with a contrastive learning objective, i.e. to distinguish the true
future frame from a set of negative samples, whereas APC is directly optimized
to minimize the distance between the predicted and the target frames. Both
types of models are subsumed by the predictive coding framework [12, 13],
which learns representations by trying to predict unobserved continuous frames
in the future. This stands in contrast to some other SSL models such as HuBERT
[14] and wav2vec 2.0 [15], which perform masked prediction of quantized units,
and which we plan to explore in future work. Apart from their practical
utility, predictive coding models are also interesting from a cognitive
perspective, since the predictive processing mechanism operationalized by
these models also underlies a range of theories on how the brain processes
sensory input [16, 17].

We use two CPC models of different sizes and one APC model. Both CPC-big and
CPC-small are taken from the Zero Speech 2021 baseline333
https://github.com/zerospeech/zerospeech2021_baseline. We use our own
implementation of APC, following [18]. The CPC and APC models differ in
several ways. CPC has a 5-layer convolution that extract 10-ms frames from
wave samples for further processing, while APC operates on 10-ms Log Mel
features. We extract the representations from the LSTM component of the models
(also called context network in CPC). The prediction horizon, as parameterized
by $K$, is 12 for CPC and 3 for APC. The specific hyperparameters are chosen
based on prior work [19, 10] and are listed in Table 1.

Table 1: Summary of model specifications. Models are trained on LibriLight
(LL) or LibriSpeech (LS) train-clean.

Model | LSTM | Extracted | Hidden | $K$ | Training
---|---|---|---|---|---
layers | layer | units | data
CPC-big | 4 | 2 | 512 | 12 | LL 6k hrs
CPC-small | 2 | 2 | 256 | 12 | LS 100 hrs
APC | 3 | 3 | 512 | 3 | LS 360 hrs

For both CPC models, the negative samples are drawn from the same speaker but
not necessarily from the same utterance as the frame to be predicted. It is
possible that a CPC model that draws negative samples without speaker
restrictions would be qualitatively different in terms of speaker and phone
encoding, but we limit the current study to the within-speaker sampling case
since it gives better phone classification results [19] and is thus the more
common setup.

###  3.1 Dataset

We use data from LibriSpeech [20], a corpus of English read speech. We perform
our initial analyses on the dev-clean portion, which contains 8 minutes of
speech for each of 19 male and 21 female speakers.444The documentation claims
that dev-clean contains 20 males and 20 females, but during our analysis, we
found that one of the “male” speakers (ID: 7976, name: JenniferRutters)
appears to be female. For testing generalization, we extract a speaker
subspace from the train-clean-100 portion, which contains 25 minutes of speech
for each of 126 male and 125 female speakers. We then evaluate on dev-clean
and test-clean by collapsing the speaker subspace learned from train-
clean-100. The test-clean portion consists of 8-minute speech per speaker for
20 male and 20 female speakers. The speakers occurred in the three sets are
mutually disjoint.

The phone labels required for our anlayses were obtained by performing forced
alignment with an acoustic model created according to the official Kaldi
recipe for LibriSpeech data555https://github.com/kaldi-
asr/kaldi/blob/master/egs/librispeech/s5/run.sh. We ignore frames aligned to
``silence'' and ``spoken noise'' labels in our analysis, leaving 39 phone
categories (i.e. $|P|=39$).

###  3.2 Evaluation

To evaluate how much speaker information is removed and its effect on the
phonetic information in the representations, we perform two types of tests.
First, we train speaker and phone probing classifiers. We use linear
classifiers, which are trained to predict either speaker or phone labels based
on a single representation frame. We train each classifier on a random half of
each speaker's utterances, using the other half for testing. We also use a
machine ABX phone discrimination test [21]. This test asks whether triphone
$x$ is more similar to triphone $a$ than $b$, where $a$ and $x$ are tokens of
the same type (e.g. `aba') and $b$ is of a different type (`apa'). The final
ABX error rate is computed by aggregating the error rate with ($a$, $b$, $x$)
iterating over all tokens combinations and triphone contrasts. In within-
speaker ABX test, each triplet are tokens produced the same speaker. In the
across-speaker setting, $a$ and $b$ are drawn from the same speaker and $x$
from a different speaker. We follow the Zero Speech challenge splits for
within- and across-speaker ABX tests.

###  3.3 Baselines for speaker normalization experiments

Previous work [8] has shown that utterance-level standardization (centering
plus rescaling) effectively removes speaker information from CPC
representations, improving performance on ABX tests and other tasks
benchmarked in the Zero Speech 2021 challenge. We therefore use utterance-
level standardization as a baseline, but also report centering alone, which we
found to work better. For analyses where data from all speakers is available
in advance, we apply these approaches at the speaker level rather than the
utterance level for a fair comparison to our methods.

##  4 Analysis of subspaces

Figure 1: Similarity (absolute value of dot product) between the principal
directions of $M_{\text{phn}}$ & $M_{\text{joint}}$ (top) and $M_{\text{spk}}$
& $M_{\text{joint}}$ (bottom), computed from dev-clean using representations
from CPC-big. The patterns are similar for APC or CPC-small. Figure 2:
Projection of $M_{\text{joint}}$ onto three of its principal dimensions:
dimension 0 (left), 1 (middle), and 12 (right). Each row in $M_{\text{joint}}$
is plotted as a dot, with those of the same speaker plotted in one column and
their phone identities colour-coded. The first 19 columns represent male and
the rest are female speakers. The legend on top is ordered by each phone's
average projection on dimension 0.

To examine the relationship between speaker and phone subspaces, we compute
the similarity between the top principal directions of
$M_{\text{spk}},M_{\text{phn}},$ and $M_{\text{joint}}$, as extracted from
dev-clean. Since we only care about measuring orthogonality, we use the
absolute value of the dot product as our similarity measure; this ranges from
0 to 1.

When comparing the principal directions of $M_{\text{phn}}$ and
$M_{\text{spk}}$, we found very low similarity: amongst the top 20 speaker
directions, their similarity with the most aligned phone direction is on
average only $0.13$ (variance: $0.002$, maximum: $0.26$), indicating that the
speaker and phone subspaces are nearly orthogonal.

We then compared the principal directions of $M_{\text{phn}}$ and
$M_{\text{spk}}$ to those of $M_{\text{joint}}$, as shown in Figure 1. Looking
at the top 12 directions of $M_{\text{joint}}$, we see that nearly all of them
are highly similar to the top principal directions of $M_{\text{phn}}$, while
only directions 1, 2, and 12 are similar to top directions of
$M_{\text{spk}}$. In other words, except for directions 1-2, all of the top 12
principal directions of $M_{\text{joint}}$ are aligned to either a phone
direction or a speaker direction. Moreover, while directions 1-2 are somewhat
aligned with both speaker direction 0 and phone direction 1, this does not
mean that speaker direction 1 is aligned with phone direction 0—in fact, their
cosine similarity is only $0.169$. Taken together, these observations further
support the orthogonality of speaker and phone directions.

The preceding analysis also informs us of the relative variance of phone and
speaker encoding. While $M_{\text{joint}}$ consists of combinations of 40
speakers and 39 phones, i.e. roughly the same number of speaker and phone
categories, most of the top principal directions are used to encode phones.
That is, the SSL representations encode much less inter-speaker variance and
use more directions to discriminate phones than speakers. This is not a
surprise since the training objective of both models entail learning to
discriminate between phones and not speakers. What remains to be explained is
how the models come to encode speaker information despite only being trained
to predict or discriminate phones within each speaker.

To better understand what is encoded in some of the top principal directions,
we visualize the projection of $M_{\text{joint}}$ onto its dimensions 0, 1,
and 12 (Figure 2). In the projection plot for dimension 0, we can see vowels
lying on the positive end and fricatives and affricates lying on the negative
end, indicating that this dimension discriminates phonetic categories along a
sonority gradient while maintaining no information about speaker differences.
Dimension 1, as well as dimension 2 which is not shown here, differentiate
between genders while the consistent ordering of the coloured dots across all
columns implies that they also capture some phonetic information. In contrast,
dimension 12 contains relatively greater inter-speaker variance, albeit still
much smaller than variance between phones captured in dimension 0 and 1. The
projection on these dimensions are consistent with our observation from the
similarity analysis.

##  5 Application to speaker normalization

Table 2: Probing (speaker & phone) and ABX (within- & across-speaker) error
rates on dev-clean.

Model | Error rate (%) | Original | Speaker-level | Speaker space
---|---|---|---|---
centered | +rescaled | collapsed
CPC-big | Speaker | 0.46 | 76.39 | +12.00 | 96.55
Phone | 24.32 | 24.29 | +2.14 | 24.49
ABX Within | 3.38 | 3.39 | +1.35 | 3.20
ABX Across | 4.11 | 3.98 | +1.34 | 3.75
CPC-small | Speaker | 11.04 | 67.95 | +11.21 | 95.54
Phone | 35.48 | 35.10 | +1.07 | 36.34
ABX Within | 6.22 | 6.14 | +1.97 | 5.11
ABX Across | 8.10 | 7.47 | +1.69 | 6.78
APC | Speaker | 15.87 | 83.42 | +4.95 | 97.23
Phone | 35.93 | 35.78 | +1.20 | 36.44
ABX Within | 6.77 | 6.50 | +0.50 | 6.34
ABX Across | 9.83 | 8.96 | +0.36 | 9.04

Figure 3: Speaker classification (blue) and across-speaker ABX (orange) on
dev-clean after collapsing speaker subspaces learned from either train-
clean-100 (solid lines) or dev-clean (dashed lines). The large dots are where
the number of dimensions covers 95% of the variance.

Collapsing speaker subspace of seen speakers We first explore speaker
normalization for a set of known speakers, by learning the speaker subspace
from dev-clean, and collapsing the speaker directions on the same dataset.
Since the number of speakers is small (40) relative to the number of
dimensions (256 or 512), we collapse all 40 of the speaker directions. (See
below for results with fewer directions collapsed.) Table 2 reports results in
comparison to speaker-level centering and standardization, because all
speakers are known in advance in this setting. Utterance-level centering and
standardization show similar results but is slightly less effective in
removing speaker information.

Table 3: Results of speaker normalization on dev-clean (top) and test-clean
(bottom) by collapsing speaker subspaces learnt from train-clean-100.

| Model | Error (%) | Original | Utt. centered | Collapsed
---|---|---|---|---|---
dev-clean | CPC-big | Speaker | 0.46 | 69.88 | 64.51
ABX Within | 3.38 | 3.43 | 3.20
ABX Across | 4.11 | 4.06 | 3.76
CPC-small | Speaker | 11.04 | 69.20 | 78.12
ABX Within | 6.22 | 5.86 | 5.29
ABX Across | 8.10 | 7.22 | 6.92
APC | Speaker | 15.87 | 77.3 | 79.74
ABX Within | 6.77 | 6.42 | 6.38
ABX Across | 9.83 | 9.02 | 9.14
test-clean | CPC-big | ABX Within | 3.29 | 3.27 | 3.10
ABX Across | 4.22 | 4.11 | 4.01
CPC-small | ABX Within | 5.86 | 5.54 | 4.85
ABX Across | 7.48 | 6.91 | 6.37
APC | ABX Within | 6.62 | 6.13 | 6.19
ABX Across | 9.47 | 8.72 | 8.98

Our method shows the best ABX results in nearly all cases, as well as close to
100% speaker error rates and little change in phone error rates, indicating
that this method removes nearly all the speaker information while improving
the accessibility of phone information for unsupervised tasks. We also see
from the baselines that removing speaker information need not always improve
phone discrimination: in particular, applying rescaling on top of centering
increases both the speaker error rate and the phone classification and ABX
error rates.

Collapsing speaker subspace of unseen speakers We first explore how results
vary depending on the number of dimensions (and % of variance) that are
removed. Figure 3 illustrates the results on CPC-big, comparing the case where
the speaker subspace is learned and applied to the same speakers and the case
where we apply it to new speakers. When we have relatively few speakers as in
dev-clean, and we collapse the speaker subspace of the same speakers,
including all of the speaker dimensions is helpful. However, we see that when
generalizing from the large number of speakers in train-clean-100, many of the
lower principal components appear to be overfitting: after about 50
dimensions, speaker error rate on the dev set increases only slowly while ABX
error also begins to rise.

Finally, after learning the speaker directions from the training set, we
choose the number of directions to collapse based on the best across-speaker
ABX scores on the development set, and present results for both the
development and test sets in Table 3. For CPC-big, CPC-small, and APC, we
collapse speaker subspaces that explain 98%, 95%, and 95% of training variance
(57, 36, and 30 dimensions, respectively). We compare to utterance-level
centering, the baseline that gives the best results for the ABX task (of those
applicable to previously unseen speakers).

Both methods remove a similar amount of speaker information, but our method is
better at the ultimate goal of improving phone discrimination in nearly all
cases (the exception being the APC model on the test set). Moreover, it can be
applied in a fully streaming setting since the speaker subspaces to collapse
are computed in advance.

##  6 Conclusions

In analyses of SSL speech representations based on predictive coding models,
we showed that speaker information and phonetic information are encoded in
orthogonal dimensions of the representation space, indicating that these
models are implicitly disentangling the two sources of information. This
insight led to a simple but effective speaker normalization technique that
requires only speaker-labeled training data, and leads to improved phone
discrimination on the test set. Important avenues for future work include
investigating the extent to which the speaker dimensions generalize to out of
domain data (e.g., other languages or genres of speech), and whether the
orthogonality we found extends to other SSL models based on different
principles (e.g., masked prediction [14, 15]).

## References

  * [1] Y. Zhang, W. Han, J. Qin, Y. Wang, A. Bapna, Z. Chen, N. Chen, B. Li, V. Axelrod, G. Wang, Z. Meng, K. Hu, A. Rosenberg, R. Prabhavalkar, D. S. Park, P. Haghani, J. Riesa, G. Perng, H. Soltau, T. Strohman, B. Ramabhadran, T. N. Sainath, P. J. Moreno, C.-C. Chiu, J. Schalkwyk, F. Beaufays, and Y. Wu, ``Google usm: Scaling automatic speech recognition beyond 100 languages,'' _ArXiv_ , vol. abs/2303.01037, 2023.
  * [2] A. Polyak, Y. Adi, J. Copet, E. Kharitonov, K. Lakhotia, W.-N. Hsu, A. Mohamed, and E. Dupoux, ``Speech Resynthesis from Discrete Disentangled Self-Supervised Representations,'' in _Proc. Interspeech 2021_ , 2021, pp. 3615–3619.
  * [3] M. A. C. Blandón and O. Räsänen, ``Analysis of predictive coding models for phonemic representation learning in small datasets,'' in _ICML 2020 Workshop on Self-supervision in Audio and Speech_ , 2020.
  * [4] M. Lavechin, M. De Seyssel, H. Titeux, H. Bredin, G. Wisniewski, A. Cristia, and E. Dupoux, ``Statistical learning bootstraps early language acquisition,'' 2022.
  * [5] D. Wells, H. Tang, and K. Richmond, ``Phonetic analysis of self-supervised representations of english speech,'' in _23rd Annual Conference of the International Speech Communication Association, INTERSPEECH 2022_. ISCA, 2022, pp. 3583–3587.
  * [6] R. Sanabria, H. Tang, and S. Goldwater, ``Analyzing Acoustic Word Embeddings from Pre-trained Self-supervised Speech Models,'' in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022\.
  * [7] M. de Seyssel, M. Lavechin, Y. Adi, E. Dupoux, and G. Wisniewski, ``Probing phoneme, language and speaker information in unsupervised speech representations,'' in _Proc. Interspeech 2022_ , 2022, pp. 1402–1406.
  * [8] B. van Niekerk, L. Nortje, M. Baas, and H. Kamper, ``Analyzing Speaker Information in Self-Supervised Models to Improve Zero-Resource Speech Processing,'' in _Proc. Interspeech 2021_ , 2021, pp. 1554–1558.
  * [9] A. Pasad, J.-C. Chou, and K. Livescu, ``Layer-wise analysis of a self-supervised speech representation model,'' in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_ , 2021, pp. 914–921.
  * [10] Y.-A. Chung, W.-N. Hsu, H. Tang, and J. Glass, ``An Unsupervised Autoregressive Model for Speech Representation Learning,'' in _Proc. Interspeech 2019_ , 2019, pp. 146–150.
  * [11] A. v. d. Oord, Y. Li, and O. Vinyals, ``Representation Learning with Contrastive Predictive Coding,'' _arXiv:1807.03748 [cs, stat]_ , Jan. 2019, arXiv: 1807.03748.
  * [12] P. Elias, ``Predictive coding–i,'' _IRE transactions on information theory_ , vol. 1, no. 1, pp. 16–24, 1955.
  * [13] B. Atal and M. R. Schroeder, ``Adaptive predictive coding of speech signals,'' _Bell Syst. Tech. J._ , vol. 49, pp. 1973–1986, 1970.
  * [14] W.-N. Hsu, B. Bolte, Y.-H. H. Tsai, K. Lakhotia, R. Salakhutdinov, and A. Mohamed, ``Hubert: Self-supervised speech representation learning by masked prediction of hidden units,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 29, pp. 3451–3460, 2021.
  * [15] A. Baevski, H. Zhou, A. rahman Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _ArXiv_ , vol. abs/2006.11477, 2020.
  * [16] R. P. N. Rao and D. H. Ballard, ``Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.'' _Nature Neuroscience_ , vol. 2, pp. 79–87, 1999.
  * [17] K. J. Friston, ``A theory of cortical responses,'' _Philosophical Transactions of the Royal Society B: Biological Sciences_ , vol. 360, pp. 815 – 836, 2005.
  * [18] G.-P. Yang, S.-L. Yeh, Y.-A. Chung, J. Glass, and H. Tang, ``Autoregressive predictive coding: A comprehensive study,'' _IEEE Journal of Selected Topics of Signal Processing_ , 2022.
  * [19] T. A. Nguyen, M. de Seyssel, P. Rozé, M. Rivière, E. Kharitonov, A. Baevski, E. Dunbar, and E. Dupoux, ``The Zero Resource Speech Benchmark 2021: Metrics and baselines for unsupervised spoken language modeling,'' _arXiv:2011.11588 [cs, eess]_ , Dec. 2020, arXiv: 2011.11588.
  * [20] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, ``Librispeech: An ASR corpus based on public domain audio books,'' in _2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. South Brisbane, Queensland, Australia: IEEE, Apr. 2015, pp. 5206–5210.
  * [21] T. Schatz, V. Peddinti, F. R. Bach, A. Jansen, H. Hermansky, and E. Dupoux, ``Evaluating speech features with the minimal-pair abx task: analysis of the classical mfc/plp pipeline,'' in _Interspeech_ , 2013.