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 \(\vec{r}\) is the vector of ratings and \(\vec{y}\) is the vector of point differences, then we can express this as the matrix equation \[X\vec{r} = \vec{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 so-called normal equations obtained by multiplying both sides by the transpose of \(X\): \[X^T X\vec{r} = X^T \vec{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 \[\vec{p} = X^T y\] with zero. This implies the constraint that the sum of the ratings must be zero. Furthermore, the resulting system \[M \vec{r} = \vec{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 \vec{r} = \vec{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\vec{r} = \vec{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 \(\vec{p}\) are the set of score differences for each team.

For example, if we expand the first term in the product \(M\vec{r}\) and set that equal to the first entry of \(\vec{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.

Real College Football Data

Let’s take a look at a somewhat bigger example. I happen to have a CSV file on my web space that has the scores of major conference college football games from 2017-2024 (excluding 2020), as well as this year’s results for the B1G so far. We can grab that data and examine a sample like so:

cfb_data = pd.read_csv('https://marksmath.org/data/cfb_games_by_conf_and_year.csv')
cfb_data.sample(10, random_state = 1)
year conf WTeamName WScore LTeamName LScore
1204 2021 PAC12 California 24 USC 14
1335 2019 Big12 Texas_Tech 45 Oklahoma_St 35
547 2023 SEC Mississippi 55 LSU 49
1797 2017 ACC Pittsburgh 31 Virginia 14
510 2023 Big12 Oklahoma_St 27 Oklahoma 24
607 2023 PAC12 Stanford 46 Colorado 43
1143 2021 SEC Mississippi 31 Mississippi_St 21
1444 2019 PAC12 Stanford 23 Washington 13
1724 2018 PAC12 Washington 31 UCLA 24
1554 2018 B1G Ohio_St 49 Indiana 26

I’m going to grab the current results for this year from the B1G:

games = cfb_data[
  cfb_data.apply(lambda r: r.year == 2025 and r.conf == "B1G", axis=1)
].copy()
m = len(games)
print(m)
games.head(13)
56
year conf WTeamName WScore LTeamName LScore
0 2025 B1G Oregon 34 Northwestern 14
1 2025 B1G USC 33 Purdue 17
2 2025 B1G Iowa 38 Rutgers 28
3 2025 B1G Indiana 63 Illinois 10
4 2025 B1G Michigan 30 Nebraska 27
5 2025 B1G USC 45 Michigan_St 31
6 2025 B1G Maryland 27 Wisconsin 10
7 2025 B1G Illinois 34 USC 32
8 2025 B1G Indiana 20 Iowa 15
9 2025 B1G Minnesota 31 Rutgers 28
10 2025 B1G Northwestern 17 UCLA 14
11 2025 B1G Oregon 30 Penn_St 24
12 2025 B1G Ohio_St 24 Washington 6

So, there’ve been 56 games played so far and we can see the first 13 games.

In order to apply Massey’s technique, we’re going to need to set up a matrix indexed by the teams. Thus, let’s assign an index to each team.

wteams = games.WTeamName
lteams = games.LTeamName
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)

games['WTeamIdx'] = games.WTeamName.apply(
  lambda t: team_to_idx_dict[t]
)
games['LTeamIdx'] = games.LTeamName.apply(
  lambda t: team_to_idx_dict[t]
)
games = games.reset_index().drop('index', axis=1)
games.tail(6)
year conf WTeamName WScore LTeamName LScore WTeamIdx LTeamIdx
50 2025 B1G Oregon 18 Iowa 16 10 2
51 2025 B1G Indiana 27 Penn_St 24 1 11
52 2025 B1G Ohio_St 34 Purdue 10 9 12
53 2025 B1G Rutgers 35 Maryland 20 13 3
54 2025 B1G Nebraska 28 UCLA 21 7 14
55 2025 B1G Wisconsin 13 Washington 10 17 16

OK, to implement Massey’s algorithm, we’ll initialize a matrix X and a vector y, and iterate through the games setting the values in X and y as we go. We use those to setup and solve the Massey equation.

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

for i,r in games.iterrows():
    X[i, r.WTeamIdx] = 1
    X[i, r.LTeamIdx] = -1
    y[i] = r.WScore - r.LScore

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

r = solve(M,p)
r
array([ -4.01752502,  28.30597823,  15.82568951, -10.32468996,
        -1.54255912, -12.4182134 ,  -8.93963107,  -4.6096638 ,
        -3.31914843,  20.02159984,  16.48575846,   5.55867224,
       -13.62625285, -12.83094573,  -3.96784419,   5.24443956,
        -3.21996918, -12.62569508])

It’s a bit hard to interpret just the vector r, so we’ll match up the values in r with the teams and sort by rating.

rank_teams = pd.DataFrame(teams)
rank_teams.columns = ['TeamName']
rank_teams['Rating'] = rank_teams.TeamName.apply(
  lambda t: r[team_to_idx_dict[t]]
)
rank_teams = rank_teams.sort_values('Rating', ascending=False).reset_index(drop=True)
rank_teams.index = rank_teams.index+1
rank_teams
TeamName Rating
1 Indiana 28.305978
2 Ohio_St 20.021600
3 Oregon 16.485758
4 Iowa 15.825690
5 Penn_St 5.558672
6 USC 5.244440
7 Michigan -1.542559
8 Washington -3.219969
9 Northwestern -3.319148
10 UCLA -3.967844
11 Illinois -4.017525
12 Nebraska -4.609664
13 Minnesota -8.939631
14 Maryland -10.324690
15 Michigan_St -12.418213
16 Wisconsin -12.625695
17 Rutgers -12.830946
18 Purdue -13.626253

Hey, not bad!

Lab problem

OK, now you’re going to give it a try! Specifically, you should rank the teams for a specific year and conference as indicated by choosing your name from the select menu below:

Once you’ve got your year and conference, you can open a Colab notebook and use the code on this page to rank your teams. Once you’ve ranked your teams, you can share your notebook with me.