Using SageMath to compute Massey products

Published: Published December 18, 2025.

Background

In my habilitation thesis, I included an appendix that used Mathematica to compute Massey products in the cohomology of configuration spaces of an oriented surface. I used Mathematica because it was what I knew, and its rewriting engine is well-suited for commutative algebra computations. However, I had to reimplemented myself a lot of things.

After some prompting from a Strasbourgeois colleague, I gave a go at reimplementing the appendix with SageMath. I have previously written on why I prefer Mathematica to SageMath for teaching. When reimplementing the appendix, I stumbled upon the same roadblocks as I had before:

the mere installation of SageMath is a chore,
documentation is a mess,
the notebook interface is archaic compared to Mathematica or VS Code (and there is some weird drama around “classic” notebook vs JupyterLab): autocomplete is limited and simplistic, no hover documentation, etc.,
debugging is inexistent; or at least I couldn’t make it work and had to rely on old-school print statements inside my code.

However, after a lot of fiddling and trial and error, I managed to do it. I have to admit, for all its warts, SageMath allows for much cleaner code in the end.

It “knows” what a vector space is and what operations I can do on subspaces.
It also “knows” about rings, ideals, and can easily make computations, although I’m concerned about optimization in some cases. (To be fair, Mathematica is very good about Gröbner basis computations too.)
Commutative differential graded algebras are already implemented in the standard library; Mathematica only knows about commutative algebras and nothing about differentials.
There are some useful helpers such as @cached_function.
NumPy is easily accessible from SageMath and has a lot of useful matrix tools, such as einsum.

Without further ado, here is the notebook in HTML format; you can also download the notebook directly.

A trick

Note that compared to my habilitation thesis, I use a fancy trick to find the nontrivial Massey product. Indeed, with Mathematica, I could just call FindInstance with a system of equations and inequations. So, since I’m looking for vectors $u,v,w \in H^1$ such that $m_2(u,v) = m_2(v,w) = 0$ and $m_3(u,v,w) \neq 0$ , I can just ask Mathematica to find a solution to that. SageMath can compute the kernel of a matrix, but I couldn’t get it to solve tensor equations; $m_2$ is bilinear, $m_3$ is trilinear, so that won’t do. Instead, I treat $v_{\mathrm{alg}} = (v_1, \dots, v_{b_1})$ as a vector with coefficients in $R = \mathbb{Q}[v_1, \dots, v_{b_1}]$ and try to make it so that $L : u \mapsto m_2(u, v_{\mathrm{alg}})$ is highly singular. I want it to have a kernel of dimension at least $2$ , to find my independent vectors $u,v$ .

To do that, I compute the zeroth and first Fitting ideals of $L$ ; this gives me an ideal $J \subset R$ of dimension $2$ . I then successively refine this ideal, guessing that $v_1 = 1$ or $v_1 = 0$ , then $v_2 = 1$ or $v_2 = 0$ will eventually lead to an ideal of dimension zero, i.e., a discrete set of points. I can then ask SageMath to solve the system of equations.

In the end, it finds the following Massey product in the cohomology of $G_A(3)$ , for $A = H^*(\Sigma_1)$ :

u = \alpha_1 + \alpha_2 - 2 \alpha_3, \quad v = \alpha_1 + \beta_1 + \alpha_2 + \beta_2 - 2 \alpha_3 - 2 \beta_3, \quad w = \beta_1 + \beta_2 - 2 \beta_3,

m_2(u,v) = m_2(v,w) = 0,

m_3(u,v,w) = -2 \alpha_2 \omega_{12} - 2 \beta_2 \omega_{12} + 2 \alpha_3 \omega_{12} + 2 \beta_3 \omega_{12} + 2 \alpha_2 \omega_{13} + 2 \beta_2 \omega_{13} - 2 \alpha_3 \omega_{13} - 2 \beta_3 \omega_{13} + 2 \alpha_1 \omega_{23} + 2 \beta_1 \omega_{23} - 2 \alpha_3 \omega_{23} - 2 \beta_3 \omega_{23}.

Massey products in $\mathrm{Conf}_n(\Sigma_g)$ [the notebook]

import numpy as np
import scipy as sp
import itertools

Definition of the model

g = 1
r = 3

abGen = [c + "_" + str(u) + "_" + str(i) for i in range(1, r+1) for u in range(1, g+1) for c in ("a", "b")]
wGen = ["w_" + str(i) + "_" + str(j) for i in range(1, r+1) for j in range(i+1, r+1)]
allGen = abGen + wGen
A = GradedCommutativeAlgebra(QQ, names = allGen)

abGen = [f"{c}_{u}_{i}" for i in range(1, r+1) for u in range(1, g+1) for c in ("a", "b")]
wGen = [f"w_{i}_{j}" for i in range(1, r+1) for j in range(i+1, r+1)]
allGen = abGen + wGen
A = GradedCommutativeAlgebra(QQ, names = allGen)

a = {(u,i) : A.gens()[allGen.index(f"a_{u}_{i}")] for i in range(1, r+1) for u in range(1, g+1)}
b = {(u,i) : A.gens()[allGen.index(f"b_{u}_{i}")] for i in range(1, r+1) for u in range(1, g+1)}
w = {(i,j) : A.gens()[allGen.index(f"w_{i}_{j}")] for i in range(1, r+1) for j in range(i+1, r+1)}

B = A.quotient(A.ideal(
    [a[u,i] * b[u,i] - a[1, i] * b[1, i] for u in range(2, g+1) for i in range(1, r+1)] +
    [(a[u,i]-a[u,j])*w[i,j] for u in range(1, g+1) for i in range(1, r+1) for j in range(i+1, r+1)] +
    [(b[u,i]-b[u,j])*w[i,j] for u in range(1, g+1) for i in range(1, r+1) for j in range(i+1, r+1)] +
    [w[i,j] * w[j, k] + w[j, k] * w[i, k] + w[i, k] * w[i, j] for i in range(1, r+1) for j in range(i+1, r+1) for k in range (j+1, r+1)] +
    [a[u,i] * b[v, i] for u in range(1, g+1) for v in range(1, g+1) for i in range(1, r+1) if u != v] +
    [a[u,i] * a[v, i] for u in range(1, g+1) for v in range(1, g+1) for i in range(1, r+1) if u != v] +
    [b[u,i] * b[v, i] for u in range(1, g+1) for v in range(1, g+1) for i in range(1, r+1) if u != v]
))

G = B.cdg_algebra({
    w[i,j]: a[1,i]*b[1,i] + a[1,j]*b[1,j] - sum((a[u,i] * b[u,j] + a[u,j] * b[u,i] for u in range(1, g+1))) for i in range(1, r+1) for j in range(i+1, r+1)
})

[len(G.basis(n)) for n in range(8)]

[1, 9, 29, 42, 29, 9, 1, 0]

[G.homology(n).dimension() for n in range(8)]

[1, 6, 14, 14, 5, 0, 0, 0]

Deformation retraction

@cached_function
def Q(n):
    return G.cocycles(n).complement().basis_matrix()

[Q(n).rank() for n in range(8)]

[0, 3, 12, 16, 8, 1, 0, 0]

dG = G.differential()

@cached_function
def Qt(n):
    return Q(n) * dG.differential_matrix(n)

[Qt(n).rank() for n in range(8)]

[0, 3, 12, 16, 8, 1, 0, 0]

@cached_function
def inj(n):
    """Basis for the homology, chosen as an orthogonal complement of Q^n + Qt^{n-1}"""
##    return G.cocycles(n).quotient(G.coboundaries(n)).lift_map().image().basis_matrix()
    return Q(n).stack(Qt(n-1)).row_space().complement().basis_matrix()

@cached_function
def rtr(n):
    basis = inj(n)
    proj = basis.transpose() * (basis * basis.transpose()).inverse() * basis
    return basis.solve_left(proj)

@cached_function
def diff(n):
    return dG.differential_matrix(n)

@cached_function
def htp(n):
    qt = Qt(n-1)
    q = Q(n-1)
    #proj_qt = qt.transpose() * (qt * qt.transpose()).inverse() * qt
    #proj_q = q.transpose() * (q * q.transpose()).inverse() * q
    #s = qt.solve_right(q)
    return qt.transpose() * (qt * qt.transpose()).inverse() * q

all(diff(n) * htp(n+1) + htp(n) * diff(n-1) == identity_matrix(len(G.basis(n))) - rtr(n)*inj(n) for n in range(8))

True

Computing the transferred $A_\infty$ -structure

@parallel
def mulcoeff(i, j, x, y, dim_prod):
    z = x * y
    if z != 0:
        return z.basis_coefficients()
    else:
        return [0]*dim_prod

@cached_function
def mult(m,n):
    bas1 = G.basis(m)
    bas2 = G.basis(n)
    len_b1 = len(bas1)
    len_b2 = len(bas2)
    dim_prod = len(G.basis(m + n))

    inputs = [
        (i, j, x, y, dim_prod)
        for i, x in enumerate(bas1)
        for j, y in enumerate(bas2)
    ]
    output = np.empty((len_b1, len_b2, dim_prod), dtype="object")

    for in_args, result in mulcoeff(inputs):
        args = in_args[0]
        i, j = args[0], args[1]
        output[i, j] = result
    return output

@cached_function
def m2(d1, d2):
    r = rtr(d1+d2).numpy(dtype="object")
    i1 = inj(d1).numpy(dtype="object")
    i2 = inj(d2).numpy(dtype="object")
    mu = mult(d1, d2)
    return np.einsum("ji,klj,mk,nl->mni", r, mu, i1, i2, optimize=True)

@cached_function
def m3(d1, d2, d3):
    i1 = inj(d1).numpy(dtype="object")
    i2 =  inj(d2).numpy(dtype="object")
    i3 = inj(d3).numpy(dtype="object")
    r = rtr(d1+d2+d3-1).numpy(dtype="object")

    # right part
    m_r = mult(d2, d3)
    h = htp(d2+d3).numpy(dtype="object")
    m_l = mult(d1, d2+d3-1)
    m3_plus = np.einsum('ji,klj,nl,opn,mk,qo,rp->mqri', r, m_l, h, m_l, i1, i2, i3, optimize=True)

    # left part
    m_l = mult(d1, d2)
    h = htp(d1+d2).numpy(dtype="object")
    m_r = mult(d1+d2-1, d3)
    m3_minus = np.einsum('ji,klj,mk,nom,qn,ro,pl->qrpi', r, m_r, h, m_l, i1, i2, i3, optimize=True)

    return m3_plus - m3_minus

m3(1,1,1).shape

(6, 6, 6, 14)

Finding the Massey product

We are looking for vectors $u,v,w \in H^1$ such that $m_2(u,v) = 0$ , $m_2(v,w) = 0$ , but $m_3(u,v,w) \neq 0$ .

Solving tensor equations is hard. We could simplify the problem by picking a random $v_0$ ; the first two equations become linear as we just need to compute the nullspaces of two matrices. Then we can pick random $u$ and $v$ from these nullspaces and check whether $m_3(u,v_0,w) \neq 0$ . If this doesn’t work, we try again.

Well actually this doesn’t work as the space of singular $v$ ‘s is of measure zero. Let’s do some algebraic geometry instead.

Ambient polynomial ring for the vector $v$

a,b,c,z = m3(1,1,1).shape

R = PolynomialRing(QQ, names=[f'v{i}' for i in range(b)])

v_vars = R.gens()

Construct Symbolic Matrix $L_A(v)$

The condition $m_2(u, v) = 0$ is equivalent to $u * L(v) = 0$ . $L_A$ has entries that are linear polynomials in $v$ .

I am working in a $C_\infty$ -algebra, so $m_2(u,v) = 0 \iff m_2(v,u) = 0$ .

L = matrix(np.einsum('j,ijk->ik', v_vars, m2(1,1)))

The ideal of solutions

We compute the Fitting ideals of $L$ , generated by minors of size $a$ and $a-1$ .

I0 = L.transpose().fitting_ideal(0)
I1 = L.transpose().fitting_ideal(1)
J = I0 + I1

len(J.gens())

4321

The problem is homogeneous so let us set $v = (1, \dots)$ . If this does not work, we can then try $v = (0, 1, \dots)$ until we find a solution.

J = J + R.ideal(v_vars[0]-1)

GB = J.groebner_basis()
JGB = R.ideal(GB)
len(JGB.gens())

4

if 1 in GB:
    print("No solution found! Try again.")

if JGB.dimension() > 0:
    print(f"Variety of dimension {JGB.dimension()}, must add conditions to find a specific instance.")

Variety of dimension 2, must add conditions to find a specific instance.

current_J = JGB
fixed_vars = {}
for i in range(1,b):
    print(f"Trying to set v_{i} = 1")
    test_ideal = current_J + R.ideal(v_vars[i] - 1)
    if test_ideal.dimension() >= 0: # It's still valid
        print(f"-> still valid, adding v_{i} = 1")
        current_J = test_ideal
        fixed_vars[i] = 1
        if current_J.dimension() == 0:
            break
    else:
        print(f"Try setting v_{i} = 0")
        test_ideal = current_J + R.ideal(v_vars[i])
        if test_ideal.dimension() >= 0:
            print(f"-> still valid, adding v_{i} = 0")
            current_J = test_ideal
            fixed_vars[i] = 0
            if current_J.dimension() == 0:
                break
J_specialized = current_J

Trying to set v_1 = 1
-> still valid, adding v_1 = 1
Trying to set v_2 = 1
-> still valid, adding v_2 = 1

JGB_specialized = R.ideal(J_specialized.groebner_basis())

pts = JGB_specialized.variety()
print(f"Found {len(pts)} solutions")

Found 1 solutions

pt = pts[0]
v_vec = vector(QQ(f.subs(pt)) for f in v_vars)

num_L = L.apply_map(lambda x: QQ(x.subs(pt)))

ker = num_L.left_kernel()
ker

Vector space of degree 6 and dimension 2 over Rational Field
Basis matrix:
[ 1 0 1 0 -2 0]
[ 0 1 0 1 0 -2]

u_vec = ker.basis()[0]
w_vec = ker.basis()[1]

massey = vector(np.einsum('i,j,k,ijkl->l', u_vec, v_vec, w_vec, m3(1,1,1)))

if massey != 0:
    print("Solution found!")

Solution found!

np.einsum('i,ij,j', u_vec, inj(1), G.basis(1)) # u

a_1_1 + a_1_2 - 2*a_1_3

np.einsum('i,ij,j', v_vec, inj(1), G.basis(1)) # v

a_1_1 + b_1_1 + a_1_2 + b_1_2 - 2*a_1_3 - 2*b_1_3

np.einsum('i,ij,j', w_vec, inj(1), G.basis(1)) # w

b_1_1 + b_1_2 - 2*b_1_3

np.einsum('i,ij,j', massey, inj(2), G.basis(2))

-2*a_1_2*w_1_2 - 2*b_1_2*w_1_2 + 2*a_1_3*w_1_2 + 2*b_1_3*w_1_2 + 2*a_1_2*w_1_3 + 2*b_1_2*w_1_3 - 2*a_1_3*w_1_3 - 2*b_1_3*w_1_3 + 2*a_1_1*w_2_3 + 2*b_1_1*w_2_3 - 2*a_1_3*w_2_3 - 2*b_1_3*w_2_3

Checking that $m_3(u,v,w)$ is not in $(u,w)$

ideal_u = matrix(np.einsum('i,ijk', u_vec, m2(1,1))).row_space()
ideal_w = matrix(np.einsum('j,ijk', w_vec, m2(1,1))).row_space()
ideal_tot = ideal_u + ideal_w
ideal_tot

Vector space of degree 14 and dimension 7 over Rational Field
Basis matrix:
[ 1 0 0 2 0 0 0 0 2 -1 0 0 0 0]
[ 0 1 0 0 0 0 0 2 0 0 0 0 0 0]
[ 0 0 1 1 0 0 0 0 2 0 0 0 0 0]
[ 0 0 0 0 1 0 0 0 0 0 0 2 0 0]
[ 0 0 0 0 0 1 0 1 0 0 0 0 0 0]
[ 0 0 0 0 0 0 1 0 1 1 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 1 1 0 0]

massey in ideal_tot

False