Sports ranking methods, 1

This is the first of a series of expository posts on matrix-theoretic sports ranking methods. This post, which owes much to discussions with TS Michael, discusses Massey’s method.

Massey’s method, currently in use by the NCAA (for football, where teams typically play each other once), was developed by Kenneth P. Massey
while an undergraduate math major in the late 1990s. We present a possible variation of Massey’s method adapted to baseball, where teams typically play each other multiple times.

There are exactly 15 pairing between these teams. These pairs are sorted lexicographically, as follows:

(1,2),(1,3),(1,4), …, (5,6).

In other words, sorted as

Army vs Bucknell, Army vs Holy Cross, Army vs Lafayette, …, Lehigh vs Navy.

The cumulative results of the 2016 regular season are given in the table below. We count only the games played in the Patriot league, but not including the Patriot league post-season tournament (see eg, the Patriot League site for details). In the table, the total score (since the teams play multiple games against each other) of the team in the vertical column on the left is listed first. In other words, ”a – b” in row $i$ and column $j$ means the total runs scored by team i against team j is a, and the total runs allowed by team i against team j is b. Here, we order the six teams as above (team 1 is Army (USMI at Westpoint), team 2 is Bucknell, and so on). For instance if X played Y and the scores were 10-0, 0-1, 0-1, 0-1, 0-1, 0-1, then the table would read 10-5 in the position of row X and column Y.

X\Y Army Bucknell Holy Cross Lafayette Lehigh Navy
Army x 14-16 14-13 14-24 10-12 8-19
Bucknell 16-14 x 27-30 18-16 23-20 10-22
Holy Cross 13-14 30-27 x 19-15 17-13 9-16
Lafayette 24-14 16-18 15-19 x 12-23 17-39
Lehigh 12-10 20-23 13-17 23-12 x 12-18
Navy 19-8 22-10 16-9 39-17 18-12 x
sm261_baseball-ranking-application_teams-digraph

Win-loss digraph of the Patriot league mens baseball from 2015

In this ordering, we record their (sum total) win-loss record (a 1 for a win, -1 for a loss) in a 15\times 6 matrix:

M = \left(\begin{array}{cccccc} -1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 0 & 0 \\ -1 & 0 & 0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 & 1 & 0 \\ -1 & 0 & 0 & 0 & 0 & 1 \\ 0 & -1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 \\ 0 & 1 & 0 & 0 & -1 & 0 \\ 0 & -1 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & -1 & 0 & 0 & 1 \\ 0 & 0 & 0 & -1 & 1 & 0 \\ 0 & 0 & 0 & -1 & 0 & 1 \\ 0 & 0 & 0 & 0 & -1 & 1 \end{array}\right).

We also record their total losses in a column vector:

{\bf b}= \left(\begin{array}{c} 2 \\ 1 \\ 10 \\ 2 \\ 11 \\ 3 \\ 2 \\ 3 \\ 14 \\ 4 \\ 14 \\ 10 \\ 11 \\ 22 \\ 6 \\ \end{array}\right).

The Massey ranking of these teams is a vector {\bf r} which best fits the equation

M{\bf r}={\bf b}.

While the corresponding linear system is over-determined, we can look for a best (in the least squares sense) approximate solution using the orthogonal projection formula

P_V = B(B^tB)^{-1}B^t,

valid for matrices B with linearly independent columns. Unfortunately, in this case B=M does not have linearly independent columns, so the formula doesn’t apply.

Massey’s clever idea is to solve

M^tM{\bf r}=M^t{\bf b}

by row-reduction and determine the rankings from the parameterized form of the solution. To this end, we compute

M^tM= \left(\begin{array}{rrrrrr} 5 & -1 & -1 & -1 & -1 & -1 \\ -1 & 5 & -1 & -1 & -1 & -1 \\ -1 & -1 & 5 & -1 & -1 & -1 \\ -1 & -1 & -1 & 5 & -1 & -1 \\ -1 & -1 & -1 & -1 & 5 & -1 \\ -1 & -1 & -1 & -1 & -1 & 5 \end{array}\right)

and

M^t{\bf b}= \left(\begin{array}{r} -24 \\ -10 \\ 10 \\ -29 \\ -10 \\ 63 \\ \end{array}\right).

Then we compute the rref of

A= (M^tM,M^t{\bf b}) = \left(\begin{array}{rrrrrr|r} 5 & -1 & -1 & -1 & -1 & -1 & -24 \\ -1 & 5 & -1 & -1 & -1 & -1 & -10 \\ -1 & -1 & 5 & -1 & -1 & -1 & 10 \\ -1 & -1 & -1 & 5 & -1 & -1 & -29 \\ -1 & -1 & -1 & -1 & 5 & -1 & -10 \\ -1 & -1 & -1 & -1 & -1 & 5 & 63 \end{array}\right),

which is

rref(M^tM,M^t{\bf b})= \left(\begin{array}{rrrrrr|r} 1 & 0 & 0 & 0 & 0 & -1 & -\frac{87}{6} \\ 0 & 1 & 0 & 0 & 0 & -1 & -\frac{73}{6} \\ 0 & 0 & 1 & 0 & 0 & -1 & -\frac{53}{6} \\ 0 & 0 & 0 & 1 & 0 & -1 & -\frac{92}{3} \\ 0 & 0 & 0 & 0 & 1 & -1 & -\frac{73}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right).

If {\bf r}=(r_1,r_2,r_3,r_4,r_5,r_6) denotes the rankings of Army, Bucknell, Holy Cross, Lafayette, Lehigh, Navy, in that order, then

r_1=r_6-\frac{87}{6},\ \ r_2=r_6-\frac{73}{6},\ \ r_3=r_6-\frac{53}{6},\ \ r_4=r_6-\frac{92}{6},\ \ r_5=r_6-\frac{73}{6}.

Therefore

Lafayette < Army = Bucknell = Lehigh < Holy Cross < Navy.

If we use this ranking to predict win/losses over the season, it would fail to correctly predict Army vs Holy Cross (Army won), Bucknell vs Lehigh, and Lafayette vs Army. This gives a prediction failure rate of 20\%.

The minimum upset ranking problem

Suppose n teams play each other, and let Team r_1 < Team r_2 < \dots < Team r_n denote some fixed ranking (where r_1,\dots,r_n is some permutation of 1,\dots,n). An upset occurs when a lower ranked team beats an upper ranked team. For each ranking, {\bf r}, let U({\bf r}) denote the total number of upsets. The minimum upset problem is to find an “efficient” construction of a ranking for which U({\bf r}) is as small as possible.

In general, let A_{ij} denote the number of times Team i beat team $j$ minus the number of times Team j beat Team i. We regard this matrix as the signed adjacency matrix of a digraph \Gamma. Our goal is to find a Hamiltonian (undirected) path through the vertices of \Gamma which goes the “wrong way” on as few edges as possible.

  1. Construct the list of spanning trees of \Gamma (regarded as an undirected graph).
  2. Construct the sublist of Hamiltonian paths (from the spanning trees of maximum degree 2).
  3. For each Hamiltonian path, compute the associated upset number: the total number of edges transversal in \Gamma going the “right way” minus the total number going the “wrong way.”
  4. Locate a Hamiltonian for which this upset number is as large as possible.

Use this sagemath/python code to compute such a Hamiltonian path.

def hamiltonian_paths(Gamma, signed_adjacency_matrix = []):
    """
    Returns a list of hamiltonian paths (spanning trees of 
    max degree <=2).

    EXAMPLES:
        sage: Gamma = graphs.GridGraph([3,3])
        sage: HP = hamiltonian_paths(Gamma)
        sage: len(HP)
        20
        sage: A = matrix(QQ,[
        [0 , -1 , 1  , -1 , -1 , -1 ],
        [1,   0 ,  -1,  1,  1,   -1  ],
        [-1 , 1 ,  0 ,  1 , 1  , -1  ],
        [1 , -1 , -1,  0 ,  -1 , -1  ],
        [1 , - 1 , - 1 , 1 , 0 , - 1  ],
        [1 ,  1  ,  1  , 1  , 1  , 0 ]
        ])
        sage: Gamma = Graph(A, format='weighted_adjacency_matrix')
        sage: HP = hamiltonian_paths(Gamma, signed_adjacency_matrix = A)
        sage: L = [sum(x[2]) for x in HP]; max(L)
        5
        sage: L.index(5)
        21
        sage: HP[21]                                 
        [Graph on 6 vertices,
         [0, 5, 2, 1, 3, 4],
         [-1, 1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1]]
        sage: L.count(5)
        1

    """
    ST = Gamma.spanning_trees()
    if signed_adjacency_matrix == []:
        HP = []
        for X in ST:
            L = X.degree_sequence()
            if max(L)<=2:
                #print L,ST.index(X), max(L)
                HP.append(X)
        return HP
    if signed_adjacency_matrix != []:
        A = signed_adjacency_matrix
        HP = []
        for X in ST:
            L = X.degree_sequence()
            if max(L)<=2:
                #VX = X.vertices()
                EX = X.edges()
		if EX[0][1] != EX[-1][1]:
                    ranking = X.shortest_path(EX[0][0],EX[-1][1])
		else:
		    ranking = X.shortest_path(EX[0][0],EX[-1][0])
		signature = [A[ranking[i]][ranking[j]] for i in range(len(ranking)-1) for j in range(i+1,len(ranking))]
                HP.append([X,ranking,signature])
        return HP

Wessell describes this method in a different way.

  1. Construct a matrix, M=(M_{ij}), with rows and columns indexed by the teams in some fixed order. The entry in the i-th row and the j-th column is defined bym_{ij}= \left\{ \begin{array}{rr} 0,& {\rm if\ team\ } i {\rm \ lost\ to\ team\ } j,\\ 1,& {\rm if\ team\ } i {\rm\ beat\ team\ } j,\\ 0, & {\rm if}\ i=j. \end{array} \right.
  2. Reorder the rows (and corresponding columns) to in a basic win-loss order: the teams that won the most games go at the
    top of M, and those that lost the most at the bottom.
  3. Randomly swap rows and their associated columns, each time checking if the
    number of upsets has gone down or not from the previous time. If it has gone down, we keep
    the swap that just happened, if not we switch the two rows and columns back and try again.

An implementaiton of this in Sagemath/python code is:

def minimum_upset_random(M,N=10):
    """
    EXAMPLES:
        sage: M = matrix(QQ,[
        [0 , 0 , 1  , 0 , 0 , 0 ],
        [1,   0 ,  0,  1,  1,   0  ],
        [0 , 1 ,  0 ,  1 , 1  , 0  ],
        [1 , 0 , 0,  0 ,  0 , 0  ],
        [1 , 0 , 0 , 1 , 0 , 0  ],
        [1 ,  1  ,  1  , 1  , 1  , 0 ]
        ])
        sage: minimum_upset_random(M)
        (
        [0 0 1 1 0 1]                    
        [1 0 0 1 0 1]                    
        [0 1 0 0 0 0]                    
        [0 0 1 0 0 0]                    
        [1 1 1 1 0 1]                    
        [0 0 1 1 0 0], [1, 2, 0, 3, 5, 4]
        )

    """
    n = len(M.rows())
    Sn = SymmetricGroup(n)
    M1 = M
    wins = sum([sum([M1[j][i] for i in range(j,6)]) for j in range(6)])
    g0 = Sn(1)
    for k in range(N):
        g = Sn.random_element()
        P = g.matrix()
        M0 = P*M1*P^(-1)
        if sum([sum([M0[j][i] for i in range(j,6)]) for j in range(6)])>wins:
            M1 = M0
            g0 = g*g0
    return M1,g0(range(n))

Mathematics of zombies

What do you do if there is a Zombie attack? Can mathematics help? This page is (humorously) dedicated to collecting links to papers or blog posted related to the mathematical models of Zombies.

 

George Romero’s 1968 Night of the Living Dead, now in the public domain, introduced reanimated ghouls, otherwise known as zombies, which craved live human flesh. Romero’s script was insired on Richard Matheson’s I Am Legend. In Romero’s version, the zombies could be killed by destroying the zombie’s brain. A dead human could, in some cases be “reanimated,” turning into a zombie. These conditions are modeled mathematically in several papers, given below.

  1. When Zombies Attack! Mathematical Modelling of a Zombie Outbreak!, paper by Mathematicians at the University of Ottawa, Philip Munz, Ioan Hudea, Joe Imad and Robert J. Smith? (yes, his last name is spelled “Smith?”).
  2. youtube video 28 Minutes Later – The Maths of Zombies , made by Dr James Grime (aka, “siningbanana”), which references the above paper.
  3. Epidemics in the presence of social attraction and repulsion, Oct 2010 Zombie paper by Evelyn Sander and Chad M. Topaz.
  4. Statistical Inference in a Zombie Outbreak Model, slides for a talk given by Des Higman, May 2010.
  5. Mathematics kills zombies dead!, 08/17/2009 blog post by “Zombie Research Society Staff”.
  6. The Mathematics of Zombies, August 18, 2009 blog post by Mike Elliot.
  7. Love, War and Zombies – Systems of Differential Equations using Sage, April 2011 slides by David Joyner. Sage commands for Love, War and Zombies talk. This was given as a Project Mosaic/M-cast broadcast.
  8. Public domain 1968 film Night of the Living Dead by George Romero.

Splitting fields of representations of generalized symmetric groups, 7

In this post, we discover which representations of the generalized symmetric group G = S_n\ wr\ C_\ell = C_\ell^n\, >\!\!\lhd \, S_n can be realized over a given abelian extension of {\mathbb{Q}}.

Let \theta_{\mu,\rho}\in G^* be the representation defined previously, where \rho\in ((S_n)_\mu)^*.

Let K\subset {\mathbb{Q}}(\zeta_\ell) be a subfield, where \zeta_\ell is a primitive \ell^{th} root of unity. Assume K contains the field generated by the values of the character of \theta_{\mu,\rho}. Assume K/{\mathbb{Q}} is Galois and let \Gamma_K=Gal({\mathbb{Q}}(\zeta_\ell)/K). Note if we regard C_\ell as a subset of {\mathbb{Q}}(\zeta_\ell) then there is an induced action of \Gamma_K on C_\ell,

\sigma:\mu \longmapsto \mu^\sigma, \ \ \ \ \ \ \ \ \ \mu\in (C_\ell)^*,\ \ \sigma\in \Gamma_K,

where \mu^\sigma(z)=\mu(\sigma^{-1}(z)), z\in C_\ell. This action extends to an action on (C_\ell^n)^*=(C_\ell^*)^n.

Key Lemma:
In the notation above, \theta_{\mu,\rho}\cong\theta_{\mu,\rho}^\sigma if and only if \mu is equivalent to \mu^\sigma under the action of S_n on (C_\ell^n)^*.

Let

n_\mu(\chi)=|\{i\ |\ 1\leq i\leq n,\ \mu_i=\chi\}|,

where \mu=(\mu_1,...,\mu_n)\in (C_\ell^n)^* and \chi\in C_\ell^*.

Theorem: The character of \theta_{\mu,\rho}\in G^* has values in K if and only if n_\mu(\chi)=n_\mu(\chi^\sigma),
for all \sigma\in \Gamma_K and all \chi\in C_\ell^*.

This theorem is proven in this paper.

We now determine the splitting field of any irreducible character of a generalized symmetric group.

Theorem: Let \chi=tr(\theta_{\rho,\mu}) be an irreducible character of G=S_n\ wr\ C_\ell. We have

Gal({\mathbb{Q}}(\zeta_\ell)/{\mathbb{Q}}(\chi))= Stab_\Gamma(\chi).

This theorem is also proven in this paper.

In the next post we shall give an example.

Linear systems of graphs in Sage

Let \Gamma be a graph. A divisor on \Gamma is an element of the free group generated by the vertices V, {\mathbb{Z}}[V].

We say that divisors D and D^\prime are linearly equivalent and write D \sim D^\prime if D^\prime-D is a principal divisor, i.e., if D^\prime = D + \text{div}(f) for some function f : V \rightarrow {\mathbb{Z}}. Note that if D and D^\prime are linearly equivalent, they must have the same degree, since the degree of every principal divisor is 0. Divisors of degree 0 are linearly equivalent if and only if they determine the same element of the Jacobian. If D is a divisor of degree 0, we denote by [D] the element of the Jacobian determined by D. A divisor D is said to be effective if D(v) \geq 0 for all vertices v. We write D \geq 0 to mean that D is effective. The linear system associated to a divisor D is the set

|D| = \{ D^\prime \in \text{Div}(\Gamma ) : D^\prime \geq 0 \text{ and } D^\prime \sim D\},

i.e., |D| is the set of all effective divisors linearly equivalent to D. Note that if D_1 \sim D_2, then |D_1| = |D_2|. We note also that if \text{deg}(D)<0, then |D| must be empty.

Sage can be used to compute the linear system of any divisor on a graph.

def linear_system(D, Gamma):
    """
    Returns linear system attached to the divisor D.

    EXAMPLES:
        sage: Gamma2 = graphs.CubeGraph(2)
        sage: Gamma1 = Gamma2.subgraph(vertices = ['00', '01'], edges = [('00', '01')])
        sage: f = [['00', '01', '10', '11'], ['00', '01', '00', '01']]
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
        True
        sage: PhiV = matrix_of_graph_morphism_vertices(Gamma1, Gamma2, f); PhiV
        [1 0 1 0]
        [0 1 0 1]
        sage: D = vector([1,0,0,1])
        sage: PhiV*D
        (1, 1)
        sage: linear_system(PhiV*D, Gamma1)
        [(2, 0), (1, 1), (0, 2)]
        sage: linear_system(D, Gamma2)
        [(0, 2, 0, 0), (0, 0, 2, 0), (1, 0, 0, 1)]
        sage: [PhiV*x for x in linear_system(D, Gamma2)]
        [(0, 2), (2, 0), (1, 1)]

    """
    Q = Gamma.laplacian_matrix()
    CS = Q.column_space()
    N = len(D.list())
    d = sum(D.list())
    #print d
    lin_sys = []
    if d < 0:
        return lin_sys
    if (d == 0) and (D in CS):
        lin_sys = [CS(0)]
        return lin_sys
    elif (d == 0):
        return lin_sys
    S = IntegerModRing(d+1)^N
    V = QQ^N
    for v in S:
        v = V(v)
        #print D-v,v,D
        if D-v in CS:
            lin_sys.append(v)
    return lin_sys

 

Differential calculus using Sagemath

Granville’s classic text book Elements of the Differential and Integral Calculus fell into the public domain and then much of it (but not all, at the time of this writing) was scanned into wikisource primarily by R. J. Hall. Granville’s entire book contains material on differential, integral, and even multivariable calculus. The material in the subset here is restricted to differential calculus topics, though contains some material which might properly belong to an elementary differential geometry course. The above-mentioned wikisource document uses mathml and latex and some Greek letter fonts.

In particular, the existence of this document owes itself primarily to three great open source projects: TeX/LaTeX, Wikipedia, and Sagemath (http://www.sagemath.org). Some material from Sean Mauch’s public domain text on Applied Mathematics, was also included.

The current latex document is due to David Joyner, who is responsible for re-formatting, editing for readability, the correction (or introduction) of typos from the scanned version, and any extra material added (for example, the hyperlinked cross references, and the Sagemath material). Please email corrections to wdjoyner@gmail.com.

Though the original text of Granville is public domain, the extra material added in this version is licensed under the GNU Free Documentation License (please see the FDL) as is most of Wikipedia.

Acknowledgements: I thank the following readers for reporting typos: Mario Pernici, Jacob Hicks.

Now available from amazon.com for $20 (not including shipping).

Discrete Fourier transforms using Sagemath

Here are some Sagemath examples for DFTs, DCTs, and DST’s. You can try copying and pasting them into the Sagemath cloud, for example.

The Sagemath dft command applies to a sequence S indexed by a set J computes the un-normalized DFT: (in Python)

[sum([S[i]*chi(zeta**(i*j)) for i in J]) for j in J]Here are some examples which explain the syntax:

sage: J = range(6)
sage: A = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: s.dft(lambda x:x^2)
   Indexed sequence: [6, 0, 0, 6, 0, 0]
   indexed by [0, 1, 2, 3, 4, 5]
sage: s.dft()
   Indexed sequence: [6, 0, 0, 0, 0, 0]
   indexed by [0, 1, 2, 3, 4, 5]
sage: G = SymmetricGroup(3)
sage: J = G.conjugacy_classes_representatives()
sage: s = IndexedSequence([1,2,3],J) # 1,2,3 are the values of a class fcn on G
sage: s.dft()   # the "scalar-valued Fourier transform" of this class fcn
    Indexed sequence: [8, 2, 2]
    indexed by [(), (1,2), (1,2,3)]
sage: J = AbelianGroup(2,[2,3],names='ab')
sage: s = IndexedSequence([1,2,3,4,5,6],J)
sage: s.dft()   # the precision of output is somewhat random and arch. dependent.
    Indexed sequence: [21.0000000000000, -2.99999999999997 - 1.73205080756885*I, -2.99999999999999 + 1.73205080756888*I, -9.00000000000000 + 0.0000000000000485744257349999*I, -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I]                
     indexed by Multiplicative Abelian Group isomorphic to C2 x C3
sage: J = CyclicPermutationGroup(6)
sage: s = IndexedSequence([1,2,3,4,5,6],J)
sage: s.dft()   # the precision of output is somewhat random and arch. dependent.
    Indexed sequence: [21.0000000000000, -2.99999999999997 - 1.73205080756885*I, -2.99999999999999 + 1.73205080756888*I, -9.00000000000000 + 0.0000000000000485744257349999*I, -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I]
     indexed by Cyclic group of order 6 as a permutation group
sage: p = 7; J = range(p); A = [kronecker_symbol(j,p) for j in J]
age: s = IndexedSequence(A,J)
sage: Fs = s.dft()
sage: c = Fs.list()[1]; [x/c for x in Fs.list()]; s.list()
     [0, 1, 1, -1, 1, -1, -1]
     [0, 1, 1, -1, 1, -1, -1]

 

The DFT of the values of the quadratic residue symbol is itself, up to a constant factor (denoted c on the last line above).

Here is a 2nd example:

sage: J = range(5)
sage: A = [ZZ(1) for i in J]
sage: s = IndexedSequence(A,J)
sage: fs = s.dft(); fs
  Indexed sequence: [5, 0, 0, 0, 0]
   indexed by [0, 1, 2, 3, 4]
sage: it = fs.idft(); it
  Indexed sequence: [1, 1, 1, 1, 1]
   indexed by [0, 1, 2, 3, 4]
age: it == s
True
sage: t = IndexedSequence(B,J)
sage: s.convolution(t)
 [1, 2, 3, 4, 5, 4, 3, 2, 1]

Here is a 3rd example:

sage: J = range(5)
sage: A = [exp(-2*pi*i*I/5) for i in J]
sage: s = IndexedSequence(A,J)
sage: s.dct()    # discrete cosine
   Indexed sequence: [2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I, 2.50000000000011 + 0.00000000000000582867087928207*I]
    indexed by [0, 1, 2, 3, 4]
sage: s.dst()        # discrete sine
  Indexed sequence: [0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I, 0.0000000000000171529457304586 - 2.49999999999915*I]
   indexed by [0, 1, 2, 3, 4]

Here is a 4th example:

sage: I = range(3)
sage: A = [ZZ(i^2)+1 for i in I]
sage: s = IndexedSequence(A,I)
sage: P1 = s.plot()
sage: P2 = s.plot_histogram()

P1 and P2 are displayed below:

The plots of P1

The plot of P1

The plot of P2

The plot of P2

Examples of graph-theoretic harmonic morphisms using Sage

In the case of simple graphs (without multiple edges or loops), a map f between graphs \Gamma_2 = (V_2,E_2) and \Gamma_1 = (V_1, E_1) can be uniquely defined by specifying where the vertices of \Gamma_2 go. If n_2 = |V_2| and n_1 = |V_1| then this is a list of length n_2 consisting of elements taken from the n_1 vertices in V_1.

Let’s look at an example.

Example: Let \Gamma_2 denote the cube graph in {\mathbb{R}}^3 and let \Gamma_1 denote the “cube graph” (actually the unit square) in {\mathbb{R}}^2.

This is the 3-diml cube graph

This is the 3-diml cube graph \Gamma_2 in Sagemath

The cycle graph on 4 vertices

The cycle graph \Gamma_1 on 4 vertices (also called the cube graph in 2-dims, created using Sagemath.

We define a map f:\Gamma_2\to \Gamma_1 by

f = [[‘000’, ‘001’, ‘010’, ‘011’, ‘100’, ‘101’, ‘110’, ‘111’], [“00”, “00”, “01”, “01”, “10”, “10”, “11”, “11”]].

Definition: For any vertex v of a graph \Gamma, we define the star St_\Gamma(v) to be a subgraph of \Gamma induced by the edges incident to v. A map f : \Gamma_2 \to \Gamma_1 is called harmonic if for all vertices v' \in V(\Gamma_2), the quantity

|\phi^{-1}(e) \cap St_{\Gamma_2}(v')|

is independent of the choice of edge e in St_{\Gamma_1}(\phi(v')).

 
Here is Python code in Sagemath which tests if a function is harmonic:

def is_harmonic_graph_morphism(Gamma1, Gamma2, f, verbose = False):
    """
    Returns True if f defines a graph-theoretic mapping
    from Gamma2 to Gamma1 that is harmonic, and False otherwise. 

    Suppose Gamma2 has n vertices. A morphism 
              f: Gamma2 -> Gamma1
    is represented by a pair of lists [L2, L1],
    where L2 is the list of all n vertices of Gamma2,
    and L1 is the list of length n of the vertices
    in Gamma1 that form the corresponding image under
    the map f.

    EXAMPLES:
        sage: Gamma2 = graphs.CubeGraph(2)
        sage: Gamma1 = Gamma2.subgraph(vertices = ['00', '01'], edges = [('00', '01')])
        sage: f = [['00', '01', '10', '11'], ['00', '01', '00', '01']]
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
        True
        sage: Gamma2 = graphs.CubeGraph(3)
        sage: Gamma1 = graphs.TetrahedralGraph()
        sage: f = [['000', '001', '010', '011', '100', '101', '110', '111'], [0, 1, 2, 3, 3, 2, 1, 0]]
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
        True
        sage: Gamma2 = graphs.CubeGraph(3)
        sage: Gamma1 = graphs.CubeGraph(2)
        sage: f = [['000', '001', '010', '011', '100', '101', '110', '111'], ["00", "00", "01", "01", "10", "10", "11", "11"]]
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
        True
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f, verbose=True)
        This [, ]] passes the check: ['000', [1, 1]]
        This [, ]] passes the check: ['001', [1, 1]]
        This [, ]] passes the check: ['010', [1, 1]]
        This [, ]] passes the check: ['011', [1, 1]]
        This [, ]] passes the check: ['100', [1, 1]]
        This [, ]] passes the check: ['101', [1, 1]]
        This [, ]] passes the check: ['110', [1, 1]]
        This [, ]] passes the check: ['111', [1, 1]]
        True
        sage: Gamma2 = graphs.TetrahedralGraph()
        sage: Gamma1 = graphs.CycleGraph(3)
        sage: f = [[0,1,2,3],[0,1,2,0]]
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f)
        False
        sage: is_harmonic_graph_morphism(Gamma1, Gamma2, f, verbose=True)
        This [, ]] passes the check: [0, [1, 1]]
        This [, ]] fails the check: [1, [2, 1]]
        This [, ]] fails the check: [2, [2, 1]]
        False

    """
    V1 = Gamma1.vertices()
    n1 = len(V1)
    V2 = Gamma2.vertices()
    n2 = len(V2)
    E1 = Gamma1.edges()
    m1 = len(E1)
    E2 = Gamma2.edges()
    m2 = len(E2)
    edges_in_common = []
    for v2 in V2:
        w = image_of_vertex_under_graph_morphism(Gamma1, Gamma2, f, v2)
        str1 = star_subgraph(Gamma1, w)
        Ew = str1.edges()
        str2 = star_subgraph(Gamma2, v2)
        Ev2 = str2.edges()
        sizes = []
        for e in Ew:
            finv_e = preimage_of_edge_under_graph_morphism(Gamma1, Gamma2, f, e)
            L = [x for x in finv_e if x in Ev2]
            sizes.append(len(L))
            #print v2,e,L
        edges_in_common.append([v2, sizes])
    ans = True
    for x in edges_in_common:
        sizes = x[1]
        S = Set(sizes)
        if S.cardinality()>1:
            ans = False
            if verbose and ans==False:
                print "This [, ]] fails the check:", x
        if verbose and ans==True:
            print "This [, ]] passes the check:", x
    return ans
            

For further details (e.g., code to

star_subgraph

, etc), just ask in the comments.

A curious conjecture about self-reciprocal polynomials

Let p be a polynomial
p(z) = a_0 + a_1z+\dots + a_nz^n\ \ \ \ \ \ a_i\in {\bf C},
and let p^* denote the reciprocal polynomial
p^*(z) = a_n + a_{n-1}z+\dots + a_0z^n=z^np(1/z).
We say p is self-reciprocal if p=p^*. For example,
1+2z+3z^2+2z^3+z^4
and
1+2z+3z^2+3z^3+2z^4+z^5
are self-reciprocal. Suppose
p(z) = a_0 + a_1z+\dots + a_dz^d+a_{d+1}z^{d+1}+a_dz^{d}+\dots +a_{1}z^{2d+1}+a_0z^{2d+2},\ \ \ \ \ \ a_i\in {\bf C},
or
p(z) = a_0 + a_1z+\dots + a_dz^d+a_{d+1}z^{d+1}+a_{d+1}z^{d+2}+a_dz^{d+3}+\dots +a_{1}z^{2d+2}+a_0z^{2d+3},\ \ \ \ \ \ a_i\in {\bf C},

The question is this: for which increasing sequences a_0<a_1<\dots a_{d+1} do the polynomial roots p(z)=0 lie on the unit circle |z|=1?

Some examples:

symmetric-increasing-coeff-plot1
This represents, when d=1 and a_0=1 and a_1=1.1 the largest a_{d+1} which has this roots property is as in the plot.

symmetric-increasing-coeff-plot2
This represents, when d=1 and a_0=1 and a_1=1.2 the largest a_{d+1} which has this roots property is as in the plot.

symmetric-increasing-coeff-plot3
This represents, when d=1 and a_0=1 and a_1=1.3 the largest a_{d+1} which has this roots property is as in the plot.

symmetric-increasing-coeff-plot4
This represents, when d=1 and a_0=1 and a_1=1.4 the largest a_{d+1} which has this roots property is as in the plot.

symmetric-increasing-coeff-plot5
This represents, when d=1 and a_0=1 and a_1=1.5 the largest a_{d+1} which has this roots property is as in the plot.

Conjecture 1: Let s:{\bf Z}_{>0}\to {\bf R}_{>0} be a ”slowly increasing” function.

  1. Odd degree case.
    If g(z)= a_0 + a_1z + \dots + a_dz^d, where a_i=s(i), then the roots of p(z)=g(z)+z^{d+1}g^*(z) all lie on the unit circle.
  2. Even degree case.
    The roots of
    p(z)=a_0 + a_1z + \dots + a_{d-1}z^{d-1} + a_{d}z^{d} + a_{d-1}z^{d+1} + \dots + a_{1}z^{2d-1}+a_0z^{2d}
    all lie on the unit circle.

The next conjecture gives an idea as to how fast a “slowly increasing” function can grow.

Conjecture 2:* Consider the even degree case only. The polynomial
p(z)=a_0 + a_1z + \dots + a_{d-1}z^{d-1} + a_{d}z^{d} + a_{d-1}z^{d+1} + \dots + a_{1}z^{2d-1}+a_0z^{2d}  ,
with a_0=1 and all a_i>0, is symmetric increasing if and only if it can be written as a product of quadratics of the form x^2-2\cos(\theta)x+1, where all but one of the factors satisfy 2\pi/4\leq \theta\leq 4\pi/3. One of the factors can be of the form x^2+\alpha x+1, for some \alpha \geq 0.

* It was pointed out to me by Els Withers that this conjecture is false.

More odd examples of p-ary bent functions

In an earlier post I discussed bent functions f:GF(3)^2\to GF(3). In this post, I’d like to give some more examples, based on a recent paper with Caroline Melles, Charles Celerier, David Phillips, and Steven Walsh, based on computations using Sage/pythpn programs I wrote with Charles Celerier.

We start with any function f:GF(p)^n\to GF(p). The Cayley graph of f is defined to be the edge-weighted digraph

\Gamma_f = (GF(p)^n, E_f ),
whose vertex set is V=V(\Gamma_f)=GF(p)^n and the set of edges is defined by

E_f =\{(u,v) \in GF(p)^n\ |\ f(u-v)\not= 0\},
where the edge (u,v)\in E_f has weight f(u-v). However, if f is even then we can (and do) regard \Gamma_f as an edge-weighted (undirected) graph.

We assume, unless stated otherwise, that f is even.

For each u\in V, define

  • N(u)=N_{\Gamma_f}(u) to be the set of all neighbors of u in \Gamma_f,
  • N(u,a)=N_{\Gamma_f}(u,a) to be the set of all neighbors v of u in \Gamma_f for which the edge (u,v)\in E_f has weight a (for each a\in GF(p)^\times = GF(p)-\{0\}),
  • N(u,0)=N_{\Gamma_f}(u,0) to be the set of all non-neighbors v of u in \Gamma_f (i.e., we have (u,v)\notin E_f),
  • the support of f is
    {\rm supp}(f)=\{v\in V\ |\ f(v)\not=0\}

Let \Gamma be a connected edge-weighted graph which is regular as a simple (unweighted) graph. The graph \Gamma is called strongly regular with parameters v, k=(k_a)_{a\in W}, \lambda=(\lambda_a)_{a\in W^3}, \mu=(\mu_a)_{a\in W^2}, denoted SRG_{W}(v,k,\lambda,\mu), if it consists of v vertices such that, for each a=(a_1,a_2)\in W^2

|N(u_1,a_1) \cap N(u_2,a_2)| =  \left\{  \begin{array}{ll}  k_{a}, & u_1=u_2,\\  \lambda_{a_1,a_2,a_3}, & u_1\in N(u_2,a_3),\ u_1\not= u_2,\\  \mu_{a}, &u_1\notin N(u_2),\ u_1\not= u_2,\\  \end{array}  \right.
where k, \lambda, \mu are as above. Here, W= GF(p) is the set of weights, including 0 (recall an “edge” has weight 0 if the vertices are not neighbors).

This example is intended to illustrate the bent function b_8 listed in the table below
\begin{array}{c|ccccccccc}  GF(3)^2 & (0, 0) & (1, 0) & (2, 0) & (0, 1) & (1, 1) & (2, 1) & (0,  2) & (1, 2) & (2, 2) \\ \hline  b_8 & 0  & 2  & 2  & 0  & 0  & 1  & 0  & 1  & 0 \\  \end{array}

Consider the finite field
GF(9) = GF(3)[x]/(x^2+1) = \{0,1,2,x,x+1,x+2,2x,2x+1,2x+2\}.
The set of non-zero quadratic residues is given by
D = \{1,2,x,2x\}.
Let \Gamma be the graph whose vertices are GF(9) and whose edges e=(a,b) are those pairs for which a-b\in D.

The graph looks like the Cayley graph for b_8 in the Figure below

Bent function b_8 on GF(3)^2

Bent function b_8 on GF(3)^2

except

8\to 0, 0\to 2x+2, 1\to 2x+1, 2\to 2x,
3\to x+2, 4\to x+1, 5\to x, 6\to 2,  7\to 1, 8\to 0.
This is a strongly regular graph with parameters (9,4,1,2).

\begin{array}{c|ccccccccc}  v       & 0       & 1       & 2       & 3       & 4       & 5       & 6       & 7 & 8 \\ \hline  N(v,0)  & 3,4,6,8 & 4,5,6,7 & 3,5,7,8 & 0,2,6,7 & 0,1,7,8 & 1,2,6,8 & 0,1,3,5 & 1,2,3,4 & 0,2,4,5 \\     N(v,1) & 5,7 & 3,8 & 4,6 & 1,8 & 2,6 & 0,7 & 2,4 & 0,5 & 1,3 \\  N(v,2) & 1,2 & 0,2 & 0,1 & 4,5 & 3,5 & 3,4 & 7,8 & 6,8 & 6,7 \\   \end{array}

The axioms of an edge-weighted strongly regular graph can be directly verified using this table.

Let S be a finite set and let R_0, R_1, \dots, R_s denote binary relations on S (subsets of S\times S). The dual of a relation R is the set

R^* = \{(x,y)\in S\times S\ |\ (y,x)\in R\}.
Assume R_0=\Delta_S= \{ (x,x)\in S\times S\ |\ x\in S\}. We say (S,R_0,R_1,\dots,R_s) is an s-class association scheme on S if the following properties hold.

  • We have a disjoint union

    S\times S = R_0\cup R_1\cup \dots \cup R_s,
    with R_i\cap R_j=\emptyset for all $i\not= j$.

  • For each i there is a j such that R_i^*=R_j (and if R_i^*=R_i for all i then we say the association scheme is symmetric).
  • For all i,j and all (x,y)\in S\times S, define

    p_{ij}(x,y) = |\{z\in S\ |\ (x,z)\in R_i, (z,y)\in R_j\}|.
    For each k, and for all x,y\in R_k, the integer p_{ij}(x,y) is a constant, denoted p_{ij}^k.

These constants p_{ij}^k are called the intersection numbers of the association scheme.

For this example of b_8, we compute the adjacency matrix associated to the members R_1 and R_2 of the association scheme (G,R_0,R_1,R_2,R_3), where G = GF(3)^2,

R_i = \{(g,h)\in  G\times G\ |\ gh^{-1} \in D_i\},\ \ \ \ \ i=1,2,
and D_i = f^{-1}(i).

Consider the following Sage computation:

sage: attach "/home/wdj/sagefiles/hadamard_transform2b.sage"
sage: FF = GF(3)
sage: V = FF^2
sage: Vlist = V.list()
sage: flist = [0,2,2,0,0,1,0,1,0]
sage: f = lambda x: GF(3)(flist[Vlist.index(x)])
sage: F = matrix(ZZ, [[f(x-y) for x in V] for y in V])
sage: F  ## weighted adjacency matrix
[0 2 2 0 0 1 0 1 0]
[2 0 2 1 0 0 0 0 1]
[2 2 0 0 1 0 1 0 0]
[0 1 0 0 2 2 0 0 1]
[0 0 1 2 0 2 1 0 0]
[1 0 0 2 2 0 0 1 0]
[0 0 1 0 1 0 0 2 2]
[1 0 0 0 0 1 2 0 2]
[0 1 0 1 0 0 2 2 0]
sage: eval1 = lambda x: int((x==1))
sage: eval2 = lambda x: int((x==2))
sage: F1 = matrix(ZZ, [[eval1(f(x-y)) for x in V] for y in V])
sage: F1
[0 0 0 0 0 1 0 1 0]
[0 0 0 1 0 0 0 0 1]
[0 0 0 0 1 0 1 0 0]
[0 1 0 0 0 0 0 0 1]
[0 0 1 0 0 0 1 0 0]
[1 0 0 0 0 0 0 1 0]
[0 0 1 0 1 0 0 0 0]
[1 0 0 0 0 1 0 0 0]
[0 1 0 1 0 0 0 0 0]
sage: F2 = matrix(ZZ, [[eval2(f(x-y)) for x in V] for y in V])
sage: F2
[0 1 1 0 0 0 0 0 0]
[1 0 1 0 0 0 0 0 0]
[1 1 0 0 0 0 0 0 0]
[0 0 0 0 1 1 0 0 0]
[0 0 0 1 0 1 0 0 0]
[0 0 0 1 1 0 0 0 0]
[0 0 0 0 0 0 0 1 1]
[0 0 0 0 0 0 1 0 1]
[0 0 0 0 0 0 1 1 0]
sage: F1*F2-F2*F1 == 0
True
sage: delta = lambda x: int((x[0]==x[1]))
sage: F3 = matrix(ZZ, [[(eval0(f(x-y))+delta([x,y]))%2 for x in V] for y in V])
sage: F3
[0 0 0 1 1 0 1 0 1]
[0 0 0 0 1 1 1 1 0]
[0 0 0 1 0 1 0 1 1]
[1 0 1 0 0 0 1 1 0]
[1 1 0 0 0 0 0 1 1]
[0 1 1 0 0 0 1 0 1]
[1 1 0 1 0 1 0 0 0]
[0 1 1 1 1 0 0 0 0]
[1 0 1 0 1 1 0 0 0]
sage: F3*F2-F2*F3==0
True
sage: F3*F1-F1*F3==0
True
sage: F0 = matrix(ZZ, [[delta([x,y]) for x in V] for y in V])
sage: F0
[1 0 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 1 0 0]
[0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 1]
sage: F1*F3 == 2*F2 + F3
True

The Sage computation above tells us that the adjacency matrix of R_1 is

A_1 =   \left(\begin{array}{rrrrrrrrr}  0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\  0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 \\  0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 \\  0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\  0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \\  1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\  0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\  1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\  0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0  \end{array}\right),
the adjacency matrix of R_2 is

A_2 =   \left(\begin{array}{rrrrrrrrr}  0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\  1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\  1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 \\  0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\  0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0  \end{array}\right),
and the adjacency matrix of R_3 is

A_3 =   \left(\begin{array}{rrrrrrrrr}  0 & 0 & 0 & 1 & 1 & 0 & 1 & 0 & 1 \\  0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 0 \\  0 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 \\  1 & 0 & 1 & 0 & 0 & 0 & 1 & 1 & 0 \\  1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\  0 & 1 & 1 & 0 & 0 & 0 & 1 & 0 & 1 \\  1 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\  0 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\  1 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 0  \end{array}\right)
Of course, the adjacency matrix of R_0 is the identity matrix. In the above computation, Sage has also verified that they commute and satisfy

A_1A_3 = 2A_2+A_3
in the adjacency ring of the association scheme.

Conjecture:
Let f:GF(p)^n\to GF(p) be a bent function, with p>2. If the level curves of f give rise to a weighted partial difference set then f is weakly regular, and the corresponding (unweighted) partial difference set is of (positive or negative) Latin square type.

For more details, see the paper [CJMPW] with Caroline Melles, Charles Celerier, David Phillips, and Steven Walsh.