# Examples of graph-theoretic harmonic maps 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 $\Gamma_2$ in Sagemath

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

# Boolean functions from the graph-theoretic perspective

This is a very short introductory survey of graph-theoretic properties of Boolean functions.

I don’t know who first studied Boolean functions for their own sake. However, the study of Boolean functions from the graph-theoretic perspective originated in Anna Bernasconi‘s thesis. More detailed presentation of the material can be found in various places. For example, Bernasconi’s thesis (e.g., see [BC]), the nice paper by P. Stanica (e.g., see [S], or his book with T. Cusick), or even my paper with Celerier, Melles and Phillips (e.g., see [CJMP], from which much of this material is literally copied).

For a given positive integer $n$, we may identify a Boolean function

$f:GF(2)^n\to GF(2),$
with its support

$\Omega_f = \{x\in GF(2)^n\ |\ f(x)=1\}.$

For each $S\subset GF(2)^n$, let $\overline{S}$ denote the set of complements $\overline{x}=x+(1,\dots,1)\in GF(2)^n$, for $x\in S$, and let $\overline{f}=f+1$ denote the complementary Boolean function. Note that

$\Omega_f^c=\Omega_{\overline{f}},$

where $S^c$ denotes the complement of $S$ in $GF(2)^n$. Let

$\omega=\omega_f=|\Omega_f|$

denote the cardinality of the support. We call a Boolean function even (resp., odd) if $\omega_f$ is even (resp., odd). We may identify a vector in $GF(2)^n$ with its support, or, if it is more convenient, with the corresponding integer in $\{0,1, \dots, 2^n-1\}.$ Let

$b:\{0,1, \dots, 2^n-1\} \to GF(2)^n$

be the binary representation ordered with least significant bit last (so that, for example, $b(1)=(0,\dots, 0, 1)\in GF(2)^n$).

Let $H_n$ denote the $2^n\times 2^n$ Hadamard matrix defined by $(H_n)_{i,j} = (-1)^{b(i)\cdot b(j)}$, for each $i,j$ such that $0\leq i,j\leq n-1$. Inductively, these can be defined by

$H_1 = \left( \begin{array}{cc} 1 & 1\\ 1 & -1 \\ \end{array} \right), \ \ \ \ \ \ H_n = \left( \begin{array}{cc} H_{n-1} & H_{n-1}\\ H_{n-1} & -H_{n-1} \\ \end{array} \right), \ \ \ \ \ n>1.$
The Walsh-Hadamard transform of $f$ is defined to be the vector in ${\mathbb{R}}^{2^n}$ whose $k$th component is

$({\mathcal{H}} f)(k) = \sum_{i \in \{0,1,\ldots,2^n-1\}}(-1)^{b(i) \cdot b(k) + f(b(i))} = (H_n (-1)^f)_k,$

where we define $(-1)^f$ as the column vector where the $i$th component is

$(-1)^f_i = (-1)^{f(b(i))},$

for $i = 0,\ldots,2^n-1$.

Example
A Boolean function of three variables cannot be bent. Let $f$ be defined by:

$\begin{array}{c|cccccccc} x_2 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 \\ x_1 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 \\ x_0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 \\ \hline (-1)^f & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 \\ {\mathcal{H}}f & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array}$
This is simply the function $f(x_0,x_1,x_2)=x_0$. It is even because

$\Omega_f = \{ (0,0,1), (0,1,1), (1,0,1), (1,1,1) \},\ \mbox{ so } \ \omega = 4.$

Here is some Sage code verifying this:

sage: from sage.crypto.boolean_function import *
sage: f = BooleanFunction([0,1,0,1,0,1,0,1])
sage: f.algebraic_normal_form()
x0
(0, -8, 0, 0, 0, 0, 0, 0)


(The Sage method walsh_hadamard_transform is off by a sign form the definition we gave.) We will return to this example later.

Let $X=(V,E)$ be the Cayley graph of $f$:

$V = GF(2)^n,\ \ \ \ E = \{(v,w)\in V\times V\ |\ f(v+w)=1\}.$
We shall assume throughout and without further mention that $f(0)\not=1,$ so $X$ has no loops. In this case, $X$ is an $\omega$-regular graph having $r$ connected components, where

$r = |GF(2)^n/{\rm Span}(\Omega_f)|.$

For each vertex $v\in V$, the set of neighbors $N(v)$ of $v$ is given by

$N(v)=v+\Omega_f,$

where $v$ is regarded as a vector and the addition is induced by the usual vector addition in $GF(2)^n$. Let $A = (A_{ij})$ be the $2^n\times 2^n$ adjacency matrix of $X$, so

$A_{ij} = f(b(i)+b(j)), \ \ \ \ \ 0\leq i,j\leq 2^n-1.$

Example:
Returning to the previous example, we construct its Cayley graph.

First, attach afsr.sage from [C] in your Sage session.

     sage: flist = [0,1,0,1,0,1,0,1]
sage: V = GF(2)ˆ3
sage: Vlist = V.list()
sage: f = lambda x: GF(2)(flist[Vlist.index(x)])
sage: X = boolean_cayley_graph(f, 3)
[0 1 0 1 0 1 0 1]
[1 0 1 0 1 0 1 0]
[0 1 0 1 0 1 0 1]
[1 0 1 0 1 0 1 0]
[0 1 0 1 0 1 0 1]
[1 0 1 0 1 0 1 0]
[0 1 0 1 0 1 0 1]
[1 0 1 0 1 0 1 0]
sage: X.spectrum()
[4, 0, 0, 0, 0, 0, 0, -4]
sage: X.show(layout="circular")


In her thesis, Bernasconi found a relationship between the spectrum of the Cayley graph $X$,

${\rm Spectrum}(X) = \{\lambda_k\ |\ 0\leq k\leq 2^n-1\},$

(the eigenvalues $\lambda_k$ of the adjacency matrix $A$) to the Walsh-Hadamard transform $\mathcal H f = H_n (-1)^f$. Note that $f$ and $(-1)^f$ are related by the equation $f=\frac 1 2 (e - (-1)^f),$ where $e=(1,1,...,1)$. She discovered the relationship

$\lambda_k = \frac 1 2 (H_n e - \mathcal H f)_k$

between the spectrum of the Cayley graph $X$ of a Boolean function and the values of the Walsh-Hadamard transform of the function. Therefore, the spectrum of $X$, is explicitly computable as an expression in terms of $f$.

References:

[BC] A. Bernasconi and B. Codenotti, Spectral analysis of Boolean functions as a graph eigenvalue problem, IEEE Trans. Computers
48(1999)345-351.

[C] C. Celerier, github repository, at https://github.com/celerier/oslo/.

[CJMP] Charles Celerier, David Joyner, Caroline Melles, David Phillips, On the Hadamard transform of monotone Boolean functions,
posted to [C].

[S] P. Stanica, Graph eigenvalues and Walsh spectrum of Boolean functions, Integers 7(2007)\# A32, 12 pages.

Here are two excellent videos of Pante Stanica on interesting applications of Boolean functions to cryptography.

This one is 50 minutes:

This one is 30 minutes:

# Remarks on the Hamming [7.4.3] code and Sage

This post simply collects some very well-known facts and observations in one place, since I was having a hard time locating a convenient reference.

Let $C$ be the binary Hamming [7,4,3] code defined by the generator matrix $G = \left(\begin{array}{rrrrrrr} 1 & 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 1 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 1 & 1 \end{array}\right)$ and check matrix $H = \left(\begin{array}{rrrrrrr} 1 & 1 & 1 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 1 \end{array}\right)$. In other words, this code is the row space of G and the kernel of H. We can enter these into Sage as follows:

sage: G = matrix(GF(2), [[1,0,0,0,1,1,1],[0,1,0,0,1,1,0],[0,0,1,0,1,0,1],[0,0,0,1,0,1,1]])
sage: G
[1 0 0 0 1 1 1]
[0 1 0 0 1 1 0]
[0 0 1 0 1 0 1]
[0 0 0 1 0 1 1]
sage: H = matrix(GF(2), [[1,1,1,0,1,0,0],[1,1,0,1,0,1,0],[1,0,1,1,0,0,1]])
sage: H
[1 1 1 0 1 0 0]
[1 1 0 1 0 1 0]
[1 0 1 1 0 0 1]
sage: C = LinearCode(G)
sage: C
Linear code of length 7, dimension 4 over Finite Field of size 2
sage: C = LinearCodeFromCheckMatrix(H)
sage: LinearCode(G) == LinearCodeFromCheckMatrix(H)
True


The generator matrix gives us a one-to-one onto map $G:GF(2)^4\to C$ defined by $m \longmapsto m\cdot G$. Using this map, the codewords are easy to describe and enumerate:

$\begin{tabular}{ccc} decimal & binary & codeword \\ 0 & 0000 & 0000000 \\ 1 & 0001 & 0001011 \\ 2 & 0010 & 0010101 \\ 3 & 0011 & 0011110 \\ 4 & 0100 & 0100110 \\ 5 & 0101 & 0101101 \\ 6 & 0110 & 0110011 \\ 7 & 0111 & 0111000 \\ 8 & 1000 & 1000111 \\ 9 & 1001 & 1001100 \\ 10 & 1010 & 1010010 \\ 11 & 1011 & 1011001 \\ 12 & 1100 & 1100001 \\ 13 & 1101 & 1101010 \\ 14 & 1110 & 1110100 \\ 15 & 1111 & 1111111 \end{tabular}$.

Using this code, we first describe a guessing game you can play with even small children.

Number Guessing game: Pick an integer from 0 to 15. I will ask you 7 yes/no questions. You may lie once.
I will tell you when you lied and what the correct number is.

Question 1: Is n in {0,1,2,3,4,5,6,7}?
(Translated: Is 1st bit of Hamming_code(n) a 0?)
Question 2: Is n in {0,1,2,3,8,9,10,11}?
(Is 2nd bit of Hamming_code(n) a 0?)
Question 3: Is n in {0,1,4,5,8,9,12,13}?
(Is 3rd bit of Hamming_code(n) a 0?)
Question 4: Is n in {0,2,4,6,8,10,12,14} (ie, is n even)?
(Is 4th bit of Hamming_code(n) a 0?)
Question 5: Is n in {0,1,6,7,10,11,12,13}?
(Is 5th bit of Hamming_code(n) a 0?)
Question 6: Is n in {0,2,5,7,9,11,12,14}?
(Is 6th bit of Hamming_code(n) a 0?)
Question 7: Is n in {0,3,4,7,9,10,13,14}?
(Is 7th bit of Hamming_code(n) a 0?)

Record the answers in a vector (0 for yes, 1 for no): $v = (v_1,v_2,...,v_7)$. This must be a codeword (no lies) or differ from a codeword by exactly one bit (1 lie). In either case, you can find n by decoding this vector.

We discuss a few decoding algorithms next.

Venn diagram decoding:

We use a simple Venn diagram to describe a decoding method.

sage: t = var('t')
sage: circle1 = parametric_plot([10*cos(t)-5,10*sin(t)+5], (t,0,2*pi))
sage: circle2 = parametric_plot([10*cos(t)+5,10*sin(t)+5], (t,0,2*pi))
sage: circle3 = parametric_plot([10*cos(t),10*sin(t)-5], (t,0,2*pi))
sage: text1 = text("$1$", (0,0))
sage: text2 = text("$2$", (-6,-2))
sage: text3 = text("$3$", (0,7))
sage: text4 = text("$4$", (6,-2))
sage: text5 = text("$5$", (-9,9))
sage: text6 = text("$6$", (9,9))
sage: text7 = text("$7$", (0,-9))
sage: textA = text("$A$", (-13,13))
sage: textB = text("$B$", (13,13))
sage: textC = text("$C$", (0,-17))
sage: text_all = text1+text2+text3+text4+text5+text6+text7+textA+textB+textC
sage: show(circle1+circle2+circle3+text_all,axes=false)


This gives us the following diagram:

Decoding algorithm:
Suppose you receive $v = ( v_1, v_2, v_3, v_4, v_5, v_6, v_7)$.
Assume at most one error is made.
Decoding process:

1. Place $v_i$ in region i of the Venn diagram.
2. For each of the circles A, B, C, determine if the sum of the bits in four regions add up to 0 or to 1. If they add to 1, say that that circle has a “parity failure”.
3. The error region is determined form the following table.

$\begin{tabular}{cc} parity failure region(s) & error position \\ none & none \\ A, B, and C & 1 \\ B, C & 4 \\ A, C & 2 \\ A, B & 3 \\ A & 5 \\ B & 6 \\ C & 7 \end{tabular}$

For example, suppose v = (1,1,1,1,1,0,1). The filled in diagram looks like

This only fails in circle B, so the table says (correctly) that the error is in the 6th bit. The decoded codeword is $c = v+e_6 = (1,1,1,1,1,1,1).$

Next, we discuss a decoding method based on the Tanner graph.

Tanner graph for hamming 7,4,3 code

The above Venn diagram corresponds to a bipartite graph, where the left “bit vertices” (1,2,3,4,5,6,7) correspond to the coordinates in the codeword and the right “check vertices” (8,9,10) correspond to the parity check equations as defined by the check matrix. This graph corresponds to the above Venn diagram, where the check vertices 8, 9, 10 were represented by circles A, B, C:

sage: Gamma = Graph({8:[1,2,3,5], 9:[1,2,4,6], 10:[1,3,4,7]})
sage: B = BipartiteGraph(Gamma)
sage: B.show()
sage: B.left
set([1, 2, 3, 4, 5, 6, 7])
sage: B.right
set([8, 9, 10])
sage: B.show()


This gives us the graph in the following picture:

Decoding algorithm:
Suppose you receive $v = ( v_1, v_2, v_3, v_4, v_5, v_6, v_7)$.
Assume at most one error is made.
Decoding process:

1. Place $v_i$ at the vertex i on the left side of the bipartite graph.
2. For each of the check vertices 8,9,10 on the right side of the graph, determine of the if the sum of the bits in the four left-hand vertices connected to it add up to 0 or to 1. If they add to 1, we say that that check vertex has a “parity failure”.
3. Those check vertices which do not fail are connected to bit vertices which we assume are correct. The remaining bit vertices
connected to check vertices which fail are to be determined (if possible) by solving the corresponding check equations.

check vertex 8: $v_2+v_3+v_4+v_5 = 0$

check vertex 9: $v_1+v_3+v_4+v_6 = 0$

check vertex 10: $v_1+v_2+v_4+v_7 = 0$

Warning: This method is not guaranteed to succeed in general. However, it does work very efficiently when the check matrix H is “sparse” and the number of 1’s in each row and column is “small.”

For example, suppose v = (1,1,1,1,1,0,1). The check vertex 9 fails its parity check, but vertex 8 and 10 do not. Therefore, only bit vertex 6 is unknown, since vertex 6 is the only one not connected to 8 and 10. This tells us that the decoding codeword is $c = (1,1,1,1,1,v_6,1)$, for some unknown $v_6$. We solve for this unknown using the check vertex equation $v_1+v_3+v_4+v_6 = 0$, giving us $v_6 = 1$. The decoded codeword is $c = (1,1,1,1,1,1,1).$

This last example was pretty simple, so let’s try $v=(0,1,1,1,1,1,1)$. In this case, we know the vertices 9 and 10 fail, so $c = (v_1,1,1,1,1,v_6,v_7)$. We solve using

$v_1+1+1+v_6 = 0$

$v_1+1+1+v_7 = 0$

This simply tells us $v_1=v_6=v_7$. By majority vote, we get $c = (1,1,1,1,1,1,1)$.

# Applications of graphs to Boolean functions

Let f be a Boolean function on $GF(2)^n$. The Cayley graph of f is defined to be the graph

$\Gamma_f = (GF(2)^n, E_f )$,

whose vertex set is $GF(2)^n$ and the set of edges is defined by

$E_f =\{ (u,v) \in GF(2)^n \times GF(2)^n \ |\ f(u+v)=1\}$.

The adjacency matrix $A_f$ is the matrix whose entries are

$A_{i,j} = f(b(i) + b(j))$,

where b(k) is the binary representation of the integer k.
Note $\Gamma_f$ is a regular graph of degree wt(f), where wt denotes the Hamming weight of f when regarded as a vector of values (of length $2^n$).

Recall that, given a graph $\Gamma$ and its adjacency matrix A, the spectrum Spec($\Gamma$) is the multi-set of eigenvalues of A.

The Walsh transform of a Boolean function f is an integer-valued function over $GF(2)^n$ that can be defined as

$W_f(u) = \sum_{x in GF(2)^n} (-1)^{f(x)+ \langle u,x\rangle}.$
A Boolean function f is bent if $|W_f(a)| = 2^{n/2}$ (this only makes sense if n is even). The Hadamard transform of a integer-valued function f is an integer-valued function over $GF(2)^n$ that can be defined as

$H_f(u) = \sum_{x in GF(2)^n} f(x)(-1)^{\langle u,x\rangle}.$
It turns out that the spectrum of $\Gamma_f$ is equal to the Hadamard transform of f when regarded as a vector of (integer) 0,1-values. (This nice fact seems to have first appeared in [2], [3].)

A graph is regular of degree r (or r-regular) if every vertex has degree r (number of edges incident to it).    We say that an r-regular graph $\Gamma$ is a strongly regular graph with parameters (v, r, d, e)  (for nonnegative integers e, d) provided, for all vertices u, v the number of vertices adjacent to both u, v is equal to

e, if u, v are adjacent,
d, if u, v are nonadjacent.

It turns out tht f is bent iff $\Gamma_f$ is strongly regular and e = d (see [3] and [4]).

The following Sage computations illustrate these and other theorems in [1], [2], [3], [4].

Consider the Boolean function $f: GF(2)^4 \to GF(2)$ given by $f(x_0,x_1,x_2) = x_0x_1+x_2x_3$.

sage: V = GF(2)^4
sage: f = lambda x: x[0]*x[1]+x[2]*x[3]
sage: CartesianProduct(range(16), range(16))
Cartesian product of [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
sage: C = CartesianProduct(range(16), range(16))
sage: Vlist = V.list()
sage: E = [(x[0],x[1]) for x in C if f(Vlist[x[0]]+Vlist[x[1]])==1]
sage: len(E)
96
sage: E = Set([Set(s) for s in E])
sage: E = [tuple(s) for s in E]
sage: Gamma = Graph(E)
sage: Gamma
Graph on 16 vertices
sage: VG = Gamma.vertices()
sage: L1 = []
sage: L2 = []
sage: for v1 in VG:
....:         for v2 in VG:
....:             N1 = Gamma.neighbors(v1)
....:         N2 = Gamma.neighbors(v2)
....:         if v1 in N2:
....:                 L1 = L1+[len([x for x in N1 if x in N2])]
....:         if not(v1 in N2) and v1!=v2:
....:                 L2 = L2+[len([x for x in N1 if x in N2])]
....:
....:
sage: L1; L2
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2]


This implies the graph is strongly regular with d=e=2.

sage: Gamma.spectrum()
[6, 2, 2, 2, 2, 2, 2, -2, -2, -2, -2, -2, -2, -2, -2, -2]
sage: [walsh_transform(f, a) for a in V]
[4, 4, 4, -4, 4, 4, 4, -4, 4, 4, 4, -4, -4, -4, -4, 4]
sage: Omega_f = [v for v in V if f(v)==1]
sage: len(Omega_f)
6
sage: Gamma.is_bipartite()
False
sage: Gamma.is_hamiltonian()
True
sage: Gamma.is_planar()
False
sage: Gamma.is_regular()
True
sage: Gamma.is_eulerian()
True
sage: Gamma.is_connected()
True
sage: Gamma.is_triangle_free()
False
sage: Gamma.diameter()
2
sage: Gamma.degree_sequence()
[6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
sage: show(Gamma)
# bent-fcns-cayley-graphs1.png


Here is the picture of the graph:

sage: H = matrix(QQ, 16, 16, [(-1)^(Vlist[x[0]]).dot_product(Vlist[x[1]]) for x in C])
sage: H
[ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1]
[ 1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1]
[ 1  1 -1 -1  1  1 -1 -1  1  1 -1 -1  1  1 -1 -1]
[ 1 -1 -1  1  1 -1 -1  1  1 -1 -1  1  1 -1 -1  1]
[ 1  1  1  1 -1 -1 -1 -1  1  1  1  1 -1 -1 -1 -1]
[ 1 -1  1 -1 -1  1 -1  1  1 -1  1 -1 -1  1 -1  1]
[ 1  1 -1 -1 -1 -1  1  1  1  1 -1 -1 -1 -1  1  1]
[ 1 -1 -1  1 -1  1  1 -1  1 -1 -1  1 -1  1  1 -1]
[ 1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1]
[ 1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1  1]
[ 1  1 -1 -1  1  1 -1 -1 -1 -1  1  1 -1 -1  1  1]
[ 1 -1 -1  1  1 -1 -1  1 -1  1  1 -1 -1  1  1 -1]
[ 1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1]
[ 1 -1  1 -1 -1  1 -1  1 -1  1 -1  1  1 -1  1 -1]
[ 1  1 -1 -1 -1 -1  1  1 -1 -1  1  1  1  1 -1 -1]
[ 1 -1 -1  1 -1  1  1 -1 -1  1  1 -1  1 -1 -1  1]
sage: flist = vector(QQ, [int(f(v)) for v in V])
sage: H*flist
(6, -2, -2, 2, -2, -2, -2, 2, -2, -2, -2, 2, 2, 2, 2, -2)
sage: A = matrix(QQ, 16, 16, [f(Vlist[x[0]]+Vlist[x[1]]) for x in C])
sage: A.eigenvalues()
[6, 2, 2, 2, 2, 2, 2, -2, -2, -2, -2, -2, -2, -2, -2, -2]



Here is another example: $f: GF(2)^3 \to GF(2)$ given by $f(x_0,x_1,x_2) = x_0x_1+x_2$.

sage: V = GF(2)^3
sage: f = lambda x: x[0]*x[1]+x[2]
sage: Omega_f = [v for v in V if f(v)==1]
sage: len(Omega_f)
4
sage: C = CartesianProduct(range(8), range(8))
sage: Vlist = V.list()
sage: E = [(x[0],x[1]) for x in C if f(Vlist[x[0]]+Vlist[x[1]])==1]
sage: E = Set([Set(s) for s in E])
sage: E = [tuple(s) for s in E]
sage: Gamma = Graph(E)
sage: Gamma
Graph on 8 vertices
sage:
sage: VG = Gamma.vertices()
sage: L1 = []
sage: L2 = []
sage: for v1 in VG:
....:         for v2 in VG:
....:             N1 = Gamma.neighbors(v1)
....:         N2 = Gamma.neighbors(v2)
....:         if v1 in N2:
....:                 L1 = L1+[len([x for x in N1 if x in N2])]
....:         if not(v1 in N2) and v1!=v2:
....:                 L2 = L2+[len([x for x in N1 if x in N2])]
....:
sage: L1; L2
[2, 0, 2, 2, 2, 2, 0, 2, 2, 2, 0, 2, 2, 2, 2, 0, 0, 2, 2, 2, 2, 0, 2, 2, 2, 0, 2, 2, 2, 2, 0, 2]
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]



This implies that the graph is not strongly regular, therefore f is not bent.

sage: Gamma.spectrum()
[4, 2, 0, 0, 0, -2, -2, -2]
sage:
sage: Gamma.is_bipartite()
False
sage: Gamma.is_hamiltonian()
True
sage: Gamma.is_planar()
False
sage: Gamma.is_regular()
True
sage: Gamma.is_eulerian()
True
sage: Gamma.is_connected()
True
sage: Gamma.is_triangle_free()
False
sage: Gamma.diameter()
2
sage: Gamma.degree_sequence()
[4, 4, 4, 4, 4, 4, 4, 4]
sage: H = matrix(QQ, 8, 8, [(-1)^(Vlist[x[0]]).dot_product(Vlist[x[1]]) for x in C])
sage: H
[ 1  1  1  1  1  1  1  1]
[ 1 -1  1 -1  1 -1  1 -1]
[ 1  1 -1 -1  1  1 -1 -1]
[ 1 -1 -1  1  1 -1 -1  1]
[ 1  1  1  1 -1 -1 -1 -1]
[ 1 -1  1 -1 -1  1 -1  1]
[ 1  1 -1 -1 -1 -1  1  1]
[ 1 -1 -1  1 -1  1  1 -1]
sage: flist = vector(QQ, [int(f(v)) for v in V])
sage: H*flist
(4, 0, 0, 0, -2, -2, -2, 2)
sage: Gamma.spectrum()
[4, 2, 0, 0, 0, -2, -2, -2]
sage: A = matrix(QQ, 8, 8, [f(Vlist[x[0]]+Vlist[x[1]]) for x in C])
sage: A.eigenvalues()
[4, 2, 0, 0, 0, -2, -2, -2]

sage: show(Gamma)
# bent-fcns-cayley-graphs2.png



Here is the picture:

REFERENCES:
[1] Pantelimon Stanica, Graph eigenvalues and Walsh spectrum of Boolean functions, INTEGERS 7(2) (2007), #A32.

[2] Anna Bernasconi, Mathematical techniques for the analysis of Boolean functions, Ph. D. dissertation TD-2/98, Universit di Pisa-Udine, 1998.

[3] Anna Bernasconi and Bruno Codenotti, Spectral Analysis of Boolean Functions as a Graph Eigenvalue Problem,  IEEE TRANSACTIONS ON COMPUTERS, VOL. 48, NO. 3, MARCH 1999.

[4] A. Bernasconi, B. Codenotti, J.M. VanderKam. A Characterization of Bent Functions in terms of Strongly Regular Graphs, IEEE Transactions on Computers, 50:9 (2001), 984-985.

[5] Michel Mitton, Minimal polynomial of Cayley graph adjacency matrix for Boolean functions, preprint, 2007.

[6] ——, On the Walsh-Fourier analysis of Boolean functions, preprint, 2006.

# Floyd-Warshall-Roy, 3

As in the previous post, let $G=(V,E)$ be a weighted digraph having $n$ vertices and let $A=A_G$ denote its$n\times n$ adjacency matrix. We identify the vertices with the set $\{1,2,\dots, n\}$.

The previous post discussed the following result, due to Sturmfels et al.

Theorem: The entry of the matrix A$\odot$ n-1 in row i and column j equals the length of a shortest path from vertex i to vertex j in G. (Here A$\odot$ n-1 denotes the n-1-st tropical power of A.)

This post discusses an implementation in Python/Sage.

Consider the following class definition.

class TropicalNumbers:
"""
Implements the tropical semiring.

EXAMPLES:
sage: T = TropicalNumbers()
sage: print T
Tropical Semiring

"""
def __init__(self):
self.identity = Infinity

def __repr__(self):
"""
Called to compute the "official" string representation of an object.
If at all possible, this should look like a valid Python expression
that could be used to recreate an object with the same value.

EXAMPLES:
sage: TropicalNumbers()
TropicalNumbers()

"""
return "TropicalNumbers()"

def __str__(self):
"""
Called to compute the "informal" string description of an object.

EXAMPLES:
sage: T = TropicalNumbers()
sage: print T
Tropical Semiring

"""
return "Tropical Semiring"

def __call__(self, a):
"""
Coerces a into the tropical semiring.

EXAMPLES:
sage: T(10)
TropicalNumber(10)
sage: print T(10)
Tropical element 10 in Tropical Semiring
"""
return TropicalNumber(a)

def __contains__(self, a):
"""
Implements "in".

EXAMPLES:
sage: T = TropicalNumbers()
sage: a = T(10)
sage: a in T
True

"""
if a in RR or a == Infinity:
return a==Infinity or (RR(a) in RR)
else:
return a==Infinity or (RR(a.element) in RR)

class TropicalNumber:
def __init__(self, a):
self.element = a
self.base_ring = TropicalNumbers()

def __repr__(self):
"""
Called to compute the "official" string representation of an object.
If at all possible, this should look like a valid Python expression
that could be used to recreate an object with the same value.

EXAMPLES:

"""
return "TropicalNumber(%s)"%self.element

def __str__(self):
"""
Called to compute the "informal" string description of an object.

EXAMPLES:
sage: T = TropicalNumbers()
sage: print T(10)
Tropical element 10 in Tropical Semiring
"""
return "%s"%(self.number())

def number(self):
return self.element

"""
Implements +. Assumes both self and other are instances of
TropicalNumber class.

EXAMPLES:
sage: T = TropicalNumbers()
sage: a = T(10)
sage: a in T
True
sage: b = T(15)
sage: a+b
10

"""
T = TropicalNumbers()
return T(min(self.element,other.element))

def __mul__(self, other):
"""
Implements multiplication *.

EXAMPLES:
sage: T = TropicalNumbers()
sage: a = T(10)
sage: a in T
True
sage: b = T(15)
sage: a*b
25
"""
T = TropicalNumbers()
return T(self.element+other.element)

class TropicalMatrix:
def __init__(self, A):
T = TropicalNumbers()
self.base_ring = T
self.row_dimen = len(A)
self.column_dimen = len(A[0])
# now we coerce the entries into T
A0 = A
m = self.row_dimen
n = self.column_dimen
for i in range(m):
for j in range(n):
A0[i][j] = T(A[i][j])
self.array = A0

def matrix(self):
"""
Returns the entries (as ordinary numbers).

EXAMPLES:
sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
sage: AT = TropicalMatrix(A)
sage: AT.matrix()
[[0, 1, 3, 7], [2, 0, 1, 3], [4, 5, 0, 1], [6, 3, 1, 0]]
"""
m = self.row_dim()
n = self.column_dim()
A0 = [[0 for i in range(n)] for j in range(m)]
for i in range(m):
for j in range(n):
A0[i][j] = (self.array[i][j]).number()
return A0

def row_dim(self):
return self.row_dimen

def column_dim(self):
return self.column_dimen

def __repr__(self):
"""
Called to compute the "official" string representation of an object.
If at all possible, this should look like a valid Python expression
that could be used to recreate an object with the same value.

EXAMPLES:

"""
return "TropicalMatrix(%s)"%self.array

def __str__(self):
"""
Called to compute the "informal" string description of an object.

EXAMPLES:

"""
return "Tropical matrix %s"%(self.matrix())

"""
Implements +. Assumes both self and other are instances of
TropicalMatrix class.

EXAMPLES:
sage: A = [[1,2,Infinity],[3,Infinity,0]]
sage: B = [[2,Infinity,1],[3,-1,1]]
sage: AT = TropicalMatrix(A)
sage: BT = TropicalMatrix(B)
sage: AT
TropicalMatrix([[TropicalNumber(1), TropicalNumber(2), TropicalNumber(+Infinity)],
[TropicalNumber(3), TropicalNumber(+Infinity), TropicalNumber(0)]])
sage: AT+BT
[[TropicalNumber(1), TropicalNumber(2), TropicalNumber(1)],
[TropicalNumber(3), TropicalNumber(-1), TropicalNumber(0)]]
"""
A = self.array
B = other.array
C = []
m = self.row_dim()
n = self.column_dim()
if m != other.row_dim:
raise ValueError, "Row dimensions must be equal."
if n != other.column_dim:
raise ValueError, "Column dimensions must be equal."
for i in range(m):
row = [A[i][j]+B[i][j] for j in range(n)] # + as tropical numbers
C.append(row)
return C

def __mul__(self, other):
"""
Implements multiplication *.

EXAMPLES:
sage: A = [[1,2,Infinity],[3,Infinity,0]]
sage: AT = TropicalMatrix(A)
sage: B = [[2,Infinity],[-1,1],[Infinity,0]]
sage: BT = TropicalMatrix(B)
sage: AT*BT
[[TropicalNumber(1), TropicalNumber(3)],
[TropicalNumber(5), TropicalNumber(0)]]
sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
sage: AT = TropicalMatrix(A)
sage: A = [[0,1,3,7],[2,0,1,3],[4,5,0,1],[6,3,1,0]]
sage: AT = TropicalMatrix(A)
sage: print AT*AT*AT
Tropical matrix [[0, 1, 2, 3], [2, 0, 1, 2], [4, 4, 0, 1], [5, 3, 1, 0]]
"""
T = TropicalNumbers()
A = self.matrix()
B = other.matrix()
C = []
mA = self.row_dim()
nA = self.column_dim()
mB = other.row_dim()
nB = other.column_dim()
if nA != mB:
raise ValueError, "Column dimension of A and row dimension of B must be equal."
for i in range(mA):
row = []
for j in range(nB):
c = T(Infinity)
for k in range(nA):
c = c+T(A[i][k])*T(B[k][j])
row.append(c.number())
C.append(row)
return TropicalMatrix(C)

This shows that the shortest distances of digraph with adjacency matrix $\begin{pmatrix} 0&1&3&7\\ 2&0&1&3\\ 4&5&0&1\\ 6&3&1&0\end{pmatrix},$ is equal to A$\odot$ 3, which is equal to $\begin{pmatrix} 0&1&2&3\\ 2&0&1&2\\ 4&4&0&1\\ 5&3&1&0\end{pmatrix}.$ This verifies an example given in chapter 1 of the book by Maclagan and Sturmfels, Introduction to Tropical Geometry .

# Floyd-Warshall-Roy, 2

Let $G = (V,E)$ be a weighted digraph having $n$ vertices and let $A = A_G$ denote its $n \times n$ adjacency matrix. We identify the vertices with the set $\{1,2,\dots, n\}$. As mentioned in the previous post, the FWR algorithm is an example of dynamic programming. In the book Algebraic Statistics for Computational Biology by L. Pachter and B. Sturmfels, and the book Introduction to Tropical Geometry by D Maclagan and Sturmfels, this algorithm is described using the tropical matrix powers of $A$. In fact, in the latter book, the corresponding section is titled “Dynamic programming”. They formulate it “tropically” as follows.

Theorem: The entry of the matrix A$\odot$ n-1 in row i and column j equals the length of a shortest path from vertex i to vertex j in G. (Here A$\odot$ n-1 denotes the n-1-st tropical power of A.)

How do you compute the tropical power of a matrix? Tropical matrix arithmetic (addition and multiplication) is the obvious extension of the addition and multiplication operations on the tropical semiring. This semiring, $({\mathbb{R}} \cup \{\infty \}, \oplus, \odot)$, is defined with the operations as follows: $x \oplus y = {\rm min}(x,y)$, $x \odot y = x+y$. The ordinary product of two $n\times n$ matrices takes $O(n^3)$ operations (roughly the same number of additions and multiplications). The tropical product of two $n\times n$ matrices still takes $O(n^3)$ operations (but only additions and comparisons, no multiplications). Of course, ordinary powers, as well as tropical powers, can be computed using the repeated squaring algorithm. This implies that the complexity of computing the N-th (tropical) power of an $n \times n$ matrix A is $O(M\log(N))$, where M is the complexity of computing the (tropical) product of two $n \times n$ matrices.

At the current time, the complexity of ordinary matrix multiplication is best estimated by either Strassen’s algorithm (for practical use) or the Coppersmith-Winograd algorithm (for theoretically best known). The CW algorithm can multiply two $n \times n$ matrices in $O(n^{2.376})$ time. However, the implied “big-O” constant is so large that the algorithm is not practical. Strassen’s algorithm can multiply two $n \times n$ matrices in $O(n^{\log_2(7)+o(1)})=O(n^{2.807})$ time.

Question: Do these complexity estimates extend to tropical matrix multiplication?

If so, then the complexity of computing the matrix A$\odot$ n-1 is $O(n^{2.807})$, using Strassen’s algorithm. This is better than the $O(n^3)$ estimate for the FWR algorithm mentioned in the previous post.

In fact, these observations seem so “obvious” that the fact that Sturmfels did not mention them in two different books, makes me wonder if the above question is more difficult to answer than it seems!