Massey ratings

This web page outlines a technique based in linear algebra to rate teams playing games in a common league based on the score differences that we see in the games. The basic idea is to perform a linear regression where the input consists of games and the output is the set of score differences.

This technique dates back to Ken Massey’s 1997 senior honors thesis at Bluefield College. My understanding of the technique is based mostly on Chapter 2 of Langville and Meyer’s text on sports ranking.

Code

We’ll describe the mathematics of Massey’s algorithm and also shows how to implement it in Python. Each code block has a clipboard icon in the upper right that you can click on to copy the code. You should be able to run the code in a scientific Python environment like Colab.

Imports

Let’s go ahead and import the libraries that we’ll need straight away:

# Data manipulation
import pandas as pd

# Numerics
import numpy as np

# SciPy's linear system solver
from scipy.linalg import solve

Description of the algorithm

Let’s describe the algorithm and take a look at a relatively simple example.

The general problem

In general, we’ve got \(m\) games played between \(n\) teams. Let’s suppose that in the \(k^{\text{th}}\) game, team \(i\) beat team \(j\) by \(y_k\) points. We’d like the teams’ ratings \(r_i\) and \(r_j\) to reflect this fact, i.e. we’d like \[r_i - r_j = y_k.\] To put it another way, we’d like our ratings of teams \(i\) and \(j\) to indicate the expected score difference between the teams, if they play again.

If \(\vect{r}\) is the vector of ratings in \(\mathbb{R}^n\) and \(\vect{y}\) is the vector of point differences in \(\mathbb{R}^m\), then we can express this as the matrix equation \[X\vect{r} = \vect{y},\] where \(X_{ki} = 1\), \(X_{kj} = -1\), and all other entries are zero.

Now, there are generally many more games than teams so that the system we’ve written down is likely overdetermined. Nonetheless, we can attempt to find an approximate solution in the least squares sense and an efficient way to do so is to solve the normal equations obtained by multiplying both sides by the transpose of \(X\): \[X^T X\vect{r} = X^T \vect{y},\]

The matrix \(M = X^T X\) formed in this process is called the Massey matrix and has some interesting properties. Note that each entry in column \(i\) of the matrix \(X\) is either \(1\) (if team \(i\) won that game), \(-1\) (if team \(i\) lost that game), or zero (if team \(i\) didn’t play in that game). As a result, \(M_{ii}\) (which is obtained by computing the dot product of the \(i^{\text{th}}\) column of \(X\) with itself) denotes the number of games played by team \(i\). Similarly, \(M_{ij}\) is the opposite of the number of games played by teams \(i\) and \(j\). Thus, the sum of each row is zero. Put another way, the columns sum to the zero vector and are therefore linearly dependent. This implies that \(M\) is singular so that the normal equations don’t have a unique solution.

A way to get around this problem is to set all the elements in some row of \(M\) equal to 1. This is conventionally done in the last row. We similarly replace the last element of \[\vect{p} = X^T y\] with zero. This implies the constraint that the sum of the ratings must be zero. Furthermore, the resulting system \[M \vect{r} = \vect{p}\] can be shown to have a unique solution.

A simple working example

As a simple example, let’s suppose that we have three teams labeled \(A\), \(B\), and \(C\) and that each team plays the other teams twice each. We might visualize the collection of games via a labeled, directed graph like so:

Note that a directed edge from team \(i\) to team \(j\) indicates that team \(i\) defeated team \(j\) in that game; the edge labels indicate the score differences.

Given the directed graph shown above, how would you rank the teams?

To rank the teams the Massey way, we start with the general equation \[X \vect{r} = \vect{y}\] which takes the specific form

\[ \begin{pmatrix}1&-1&0\\1&-1&0\\0&1&-1\\0&-1&1\\1&0&-1\\-1&0&1\end{pmatrix} \begin{pmatrix}r_0\\r_1\\r_2\end{pmatrix} = \begin{pmatrix}3\\2\\1\\2\\3\\1\end{pmatrix}. \]

If we multiply both sides by \(X^T\) and modify the last row according Massey’s strategy, the equation becomes

\[ \begin{pmatrix}4&-2&-2\\-2&4&-2\\1&1&1\end{pmatrix} \begin{pmatrix}r_0\\r_1\\r_2\end{pmatrix} = \begin{pmatrix}7\\-6\\0\end{pmatrix}. \]

The solution of this system is

\[ \begin{pmatrix}r_0\\r_1\\r_2\end{pmatrix} = \begin{pmatrix}7/6\\-1\\-1/6\end{pmatrix}. \]

In particular, note that \(r_1 < r_2 < r_0\), which is probably in accordance with the order you expected.

The Massey equation

It’s worth noting that the Massey equation \[M\vect{r} = \vect{p}\] can itself be written down and interpreted without reference to \(X\) or \(y\). As we’ve already noted, \(M_{ii}\) denotes the number of games played by team \(i\) and \(M_{ij}\) is the opposite of the number of games played by teams \(i\) and \(j\). Similarly, the entries of \(\vect{p}\) are the set of score differences for each team.

For example, if we expand the first term in the product \(M\vect{r}\) and set that equal to the first entry of \(\vect{p}\), we get \[4r_0 - 2r_1 - 2r_2 = 7.\] This indicates that team \(0\) played four games, two against each other team, and that it outscored its opponents by a combined score of seven. If we solve for \(r_0\), we get \[r_0 = \frac{1}{2}r_1 + \frac{1}{2}r_2 + \frac{7}{2}.\] In general, the rating of team \(i\) is a combination of the ratings of the other teams together with the number of points that team \(i\) outscored those other teams by.

Implementation

Now, we’ll implement the technique in code for both the simple example and a more realistic example.

The simple example

Here’s an implementation of the technique for the simple example described above in Python:

X = np.array([
    [1, -1, 0],
    [1, -1, 0],
    [0, 1, -1],
    [0, -1, 1],
    [1, 0, -1],
    [-1, 0, 1]
  ])
M = X.transpose().dot(X)
M[-1] = np.ones(len(M))

y = np.array([[3], [2], [1], [2], [3], [1]])
y = np.array([3, 2, 1, 2, 3, 1])
y = X.transpose().dot(y)
y[-1]=0

solve(M,y)

array([ 1.16666667, -1.        , -0.16666667])

This is exactly the [7/6,-1,-1/6] ratings that we mentioned above.

The Big South

Let’s take a look at a somewhat bigger example, namely the current Big South. I happen to have a CSV file on my web space that has the scores of all Big South games played this season up through last night (Thur, Feb 12). We can download it right off my web space and take a look:

games = pd.read_csv('https://marksmath.org/data/BigSouthBasketball2026.csv')  
m = len(games)
games.tail(8)

	name1	score1	name2	score2	date	idx1	idx2
43	High_Point	86	Radford	77	2026-02-07	2	5
44	Winthrop	79	Longwood	74	2026-02-07	8	3
45	Presbyterian	68	Gardner_Webb	62	2026-02-07	4	1
46	SC_Upstate	100	Charleston_So	94	2026-02-07	6	0
47	Winthrop	103	Gardner_Webb	85	2026-02-12	8	1
48	UNC_Asheville	79	Longwood	74	2026-02-12	7	3
49	Charleston_So	84	Presbyterian	67	2026-02-12	0	4
50	High_Point	95	SC_Upstate	70	2026-02-12	2	6

It looks like there have been 51 games played so far. Let’s see how many teams there are:

wteams = games.name1
lteams = games.name2
teams = pd.concat([wteams, lteams]).drop_duplicates().sort_values()
team_to_idx_dict = dict([(v,i) for i,v in enumerate(teams)])
n = len(teams)
n

We could plot the directed graph but it’s a bit busy. The same information is displayed by a chord diagram, though:

So, we’ve got 9 teams who played 51 games. Our design matrix \(X\) will have 51 rows and 9 columns. It will map 9 dimensional space into 51 dimensional space and the solution to the regression problem will be the nine-dimensional vector of ratings. Here’s how to setup and solve the regression problem using the scipy.linalg.solve function:

X = np.zeros([m,n])
y = np.zeros([m])
X.transpose().dot(y)

for i,r in games.iterrows():
    X[i, r.idx1] = 1
    X[i, r.idx2] = -1
    y[i] = r.score1 - r.score2

M = X.transpose().dot(X)
M[-1,:] = 1
p =  X.transpose().dot(y)
p[-1] = 0

r = solve(M,p)
r

array([ -1.85118093, -17.62671491,  13.2013467 ,   2.18874022,
        -2.07113002,   2.62433563,  -6.34438783,   3.51718027,
         6.36181086])

Finally, we’ll sort those teams by their ratings in descending order to get a ranking:

teams = pd.DataFrame(teams)
teams.columns = ['TeamName']
teams['Rating'] = teams.TeamName.apply(
  lambda t: r[team_to_idx_dict[t]]
)
teams = teams.sort_values('Rating', ascending=False)
teams.reset_index(inplace=True)
teams.drop('index', inplace=True, axis=1)
teams

	TeamName	Rating
0	High_Point	13.201347
1	Winthrop	6.361811
2	UNC_Asheville	3.517180
3	Radford	2.624336
4	Longwood	2.188740
5	Charleston_So	-1.851181
6	Presbyterian	-2.071130
7	SC_Upstate	-6.344388
8	Gardner_Webb	-17.626715

Not bad!