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

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.

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 kth 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 ith 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
sage: f.walsh_hadamard_transform()
(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)
     sage: X.adjacency_matrix()
     [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 itsn\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

    def __add__(self, other):
        """
        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())

    def __add__(self, other):
        """
        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!