Fast algorithms for anti-distance matrices as a generalization of Boolean matrices

TL;DR

This paper generalizes Boolean matrices to matrices over a semiring, extending the connection between Boolean matrix multiplication and graph algorithms like Warshall's and Floyd-Warshall, with efficient vectorized implementations.

Contribution

It introduces a generalized matrix framework over a semiring that extends Boolean matrices and maintains algorithmic equivalence with Floyd-Warshall, along with optimized implementation techniques.

Findings

Generalized matrices preserve Floyd-Warshall algorithm properties.

Efficient vectorized operations improve computational performance.

Framework unifies Boolean and weighted graph matrix computations.

Abstract

We show that Boolean matrix multiplication, computed as a sum of products of column vectors with row vectors, is essentially the same as Warshall's algorithm for computing the transitive closure matrix of a graph from its adjacency matrix. Warshall's algorithm can be generalized to Floyd's algorithm for computing the distance matrix of a graph with weighted edges. We will generalize Boolean matrices in the same way, keeping matrix multiplication essentially equivalent to the Floyd-Warshall algorithm. This way, we get matrices over a semiring, which are similar to the so-called "funny matrices". We discuss our implementation of operations on Boolean matrices and on their generalization, which make use of vector instructions.

Equations

$$I\vee A\vee (A\wedge A)\vee (A\wedge A\wedge A)\vee \cdots$$

$$A\vee (A\wedge A)\vee (A\wedge A\wedge A)\vee \cdots$$

$$\vee :\mathcal{B}\times \mathcal{B}\to \mathcal{B}\text{mbox}bya\vee b=\max \{a,b\}$$

$$\wedge :\mathcal{B}\times \mathcal{B}\to \mathcal{B}\text{mbox}bya\wedge b=\min \{a,b\}$$

$$(A\vee B{)}_{ij}=A_{ij}\vee B_{ij}(1\le i\le R,1\le j\le C)$$

$$(A\ast B{)}_{ij}=(A_{i1}\wedge B_{1j})\vee \cdots \vee (A_{iC}\wedge B_{Cj})(1\le i\le R,1\le j\le C^{\prime })$$

$$A\vee (A\ast A)\vee (A\ast A\ast A)\vee \cdots$$

$$A\vee (A\ast A)\vee (A\ast A\ast A)\vee \cdots \vee A^{\ast D^{\prime }}=A\ast (I_{D^{\prime }}\vee A{)}^{\ast (D^{\prime }-1)}$$

$$\ast :\mathcal{C}\times \mathcal{C}\to \mathcal{C}\text{mbox}bya\ast b={\delta }^{-1}(\delta (a)+\delta (b))$$

$$(a\vee b)\ast c$$

$$\displaystyle=\delta^{-1}\big{(}\big{(}\delta(a)\wedge\delta(b)\big{)}+\delta(c)\big{)}$$

$$\displaystyle=\delta^{-1}\big{(}\big{(}\delta(a)+\delta(c)\big{)}\wedge\big{(}\delta(b)+\delta(c)\big{)}\big{)}$$

$$\displaystyle=\delta^{-1}\big{(}\delta(a)+\delta(c)\big{)}\vee\delta^{-1}\big{(}\delta(b)+\delta(c)\big{)}$$

$$=(a\ast c)\vee (b\ast c)$$

$$\ast :\mathcal{D}\times \mathcal{D}\to \mathcal{D}\text{mbox}bya\ast b=\max \{a+b-1,0\}$$

$$\left(\begin{array}[]{cccc}0&\infty&\cdots&\infty\\ \infty&0&\ddots&\vdots\\ \vdots&\ddots&\ddots&\infty\\ \infty&\cdots&\infty&0\end{array}\right)$$

$$\neg (A\vee B)=(\neg A)\wedge (\neg B)\neg (A\wedge B)=(\neg A)\vee (\neg B)$$

$$\neg (A\ast B)=(\neg A)+(\neg B)\neg (A+B)=(\neg A)\ast (\neg B)$$

Full text

Warshall’s algorithm can be generalized to Floyd’s algorithm for computing the distance matrix of a graph with weighted edges. We will generalize Boolean matrices in the same way, keeping matrix multiplication essentially equivalent to the Floyd-Warshall algorithm. This way, we get matrices over a semiring, which are similar to the so-called “funny matrices”.

We discuss our implementation of operations on Boolean matrices and on their generalization, which make use of vector instructions.

1 Introduction

In [1, pp. 200-206], it has been proved that Boolean matrix multiplication has the same time complexity as computing the reflexive transitive closure matrix which corresponds to a Boolean matrix $A$ as adjacency matrix of a graph, i.e.

[TABLE]

where $I$ is the identity matrix. More generally, it has been proved that matrix multiplication over a semiring has the same time complexity as computing (1)

The computation of the transitive closure matrix which corresponds $A$, i.e.

[TABLE]

has the same complexity as well. Namely, we can get (1) from (2) by addition of $I$, and we can obtain (2) from (1) by multiplication with $A$. More generally, matrix multiplication over a semiring has the same time complexity as computing (2).

We connect Boolean matrix multiplication and taking transitive closure in another way. Namely, we show that Boolean matrix multiplication, computed as a sum of products of column vectors with row vectors, is essentially the same as Warshall’s algorithm for computing the transitive closure matrix.

Warshall’s algorithm can be generalized to Floyd’s algorithm for computing the distance matrix of a graph with weighted edges. We will generalize Boolean matrices in the same way, keeping matrix multiplication essentially similar to the Floyd-Warshall algorithm. This way, we get matrices over a semiring, which we call “anti-distance matrices”.

Our anti-distance matrices are similar to the distance-matrices in [1, pp. 200-206] and the “funny matrices” in [2]. This is not very surprising, because these matrices are used to compute the shortest paths for all pairs of vertices, briefly APSP (All Pairs shortest Path).

Furthermore, we discuss our implementations of matrix operations on Boolean matrices and anti-distance matrices. This includes the Floyd-Warshall algorithm. Alternatives to the Floyd-Warshall algorithm may involve some sort of graph traversal. This leads to overhead, and is harder to parallellize. But parallellization is indeed possible, namely with GPU shaders, see e.g. [5].

But we will implement our matrix operations in a simpler way, namely by way of vector instructions of the CPU. Using such instructions, the Floyd-Warshall-algorithm and our matrix multiplication have been implemented more or less in [4], where our matrix multiplication is called Floyd-Warshall-abc, with the three letters referring to the three distinct parameters in $C=A\ast B$.

The authors of [4] wrote a program which writes the instructions for the matrix operation. We choose a simpler approach, namely we wrote g++ template header files with inline assembly instructions. Furthermore, the authors of [4] use signed edge weights without saturation, where we use unsigned edge weights with saturation.

2 Boolean matrix multiplication and Warshall’s algorithm

Let $\mathcal{B}=\{0,1\}$ denote the Boolean set in numerical fashion. We define an addition on $\mathcal{B}$, namely

[TABLE]

Furthermore, we define a multiplication on $B$, namely

[TABLE]

Using this addition and multiplication, we can define matrix addition and multiplication over $\mathcal{B}$, namely in the same way as usual. So if we take $A,B\in\operatorname{Mat}_{R,C}(\mathcal{B})$, then

[TABLE]

If we take $A\in\operatorname{Mat}_{R,C}(\mathcal{B})$ and $B\in\operatorname{Mat}_{C,C^{\prime}}(\mathcal{B})$, then

[TABLE]

Let $A_{{\updownarrow}k}$ denote the $k$-th column of $A$ and $B_{k{\leftrightarrow}}$ denote the $k$-th row of $B$. Write $0_{R\times C}$ for the zero matrix with $R$ rows and $C$ columns. We can compute $P=A\ast B$ by way of the following code:

  $P:=0_{R\times C}$

  for $k:=1$ to $C$ do

     $P:=P\vee(A_{{\updownarrow}k}\ast B_{k{\leftrightarrow}})$

  end for

  return  $P$

We can interpret a square Boolean matrix, say $A\in\operatorname{Mat}_{D,D}(\mathcal{B})$, as the adjacency matrix of a directed graph. Let $B\in\operatorname{Mat}_{D,D}(\mathcal{B})$ be the adjacency matrix for a directed graph with the same vertices, but different edges as the graph of which $A$ is the adjacency matrix. Then $A\vee B$ is adjacency matrix of the union of $A$ and $B$. The zero matrix $0_{D\times D}$ is the neutral element of $\vee$, and the interpretation of $0_{D\times D}$ is the empty graph.

The interpretation of $A\ast B$ is a little more difficult. In the associated graph, there is an edge from $i$ to $j$, if and only if there is a vertex $k$ such that there is a $B$-edge from $i$ to $k$ and an $A$-edge from $k$ to $j$. The identity matrix $I_{D}$ of size $D$ is the neutral element of $\ast$, and the interpretation of $I_{D}$ is a graph with $D$ edges, namely an edge of every vertex to itself.

The matrix which corresponds to the transitive closure of the graph is the infinite sum

[TABLE]

which is actually just

[TABLE]

for every $D^{\prime}\geq D$, where $A^{\ast D^{\prime}}$ is the product of $D^{\prime}$ copies of $A$.

$(I_{D^{\prime}}\vee A)^{\ast(D^{\prime}-1)}$ can be computed by way of repeated squaring, but a better method to compute (3) is Warshall’s algorithm:

  $T:=A$

  for $k:=1$ to $D$ do

     $T:=T\vee(T_{{\updownarrow}k}\ast T_{k{\leftrightarrow}})$

  end for

  return  $T$

Notice that Warshall’s algorithm is very similar to the multiplication algorithm.

3 Anti-distance matrices and Floyd’s algorithm

If we compare the Warshall algorithm with the multiplication algorithm, we see that for the two input matrices $A$ and $B$ and the output matrix $P$ in the multiplication algorithm, the same matrix $T$ is taken in the Warshall algorithm. We will generalize the Booleans matrices in such a way, that Warshall’s algorithm becomes Floyd’s algorithm for distances.

We do this in a way to preserve the interpretations of $A\vee B$, $0_{D\times D}$, $A\ast B$, and $I_{D}$ as much as possible. For that reason, a $1$ indicates distance [math], and a [math] indicates distance $\infty$. Any number in between [math] and $1$ indicates a distance between [math] and $\infty$, i.e. a proper distance.

So let us replace the Boolean set $\mathcal{B}$ by the interval $\mathcal{C}=[0,1]$. Let $\delta:\mathcal{C}\rightarrow[0,\infty]$ denote the distance associated to an element of $\mathcal{C}$. Then $\delta(0)=\infty$ and $\delta(1)=0$, so it is natural to impose that $\delta$ is decreasing.

The addition operator on $\mathcal{C}$ should take the minimum distance, which is accomplished by taking the maximum value, because $\delta$ is decreasing. So we can take $\vee$ as the addition operator for $\mathcal{C}$, just as for $\mathcal{B}$.

The multiplication operator on $\mathcal{C}$ should add distances, but $\wedge$ only does this for the distances [math] and $\infty$, which are in fact the only distances in the Boolean case. So we cannot preserve $\wedge$. No, we use another multiplication on $\mathcal{C}$, namely we define

[TABLE]

If we extend $\wedge$ to $[0,\infty]$, then

[TABLE]

and similarly $c\ast(a\vee b)=(c\ast a)\vee(c\ast b)$, so the distributive law is fulfilled. From this, the distributive law for matrices follows.

We can take $\delta(a)=-\log(a)$, where $\log$ is the logarithm with respect to some base. But infinite intervals and functions like $\log$ and $\exp$ are not the most convenient things in computer practice. If we know in advance that distances will not be in the interval $[S,\infty)$ for some positive number $S$, or we are just not interested in distinguishing distances in the interval $[S,\infty]$, we can take $\delta(a)=S\cdot(1-a)$.

We can go even further. Since continuous intervals are not so convenient in computer practice, we can replace $\mathcal{C}=[0,1]$ by $\mathcal{D}=\{0,\tfrac{1}{S},\tfrac{2}{S},\ldots,1\}$. Then $\delta(\mathcal{D})=\{0,1,2,\ldots,S\}$. We call $S$ is the saturation distance, and define

[TABLE]

which is just as above with $\delta(a)=S\cdot(1-a)$.

In [1, pp. 200-206] and [2], similars things have been done. But the deduction from adjacency matrices is missing. Furthermore, operations are done directly on $\delta(\mathcal{C})=[0,\infty]$ instead of $\mathcal{C}$, with $\wedge$ as the addition operator and $+$ as the multiplication operator. That $+$ plays the role of multiplication leads to “funny matrix multiplication”, where the associated identity matrix is as funny as

[TABLE]

The neutral matrix for addition by way of $\wedge$ is the matrix of which all entries are $\infty$ instead of the zero matrix.

4 Fast algorithms for anti-distance matrices

In [1, pp. 200-206], it is proved that the computation of (3) has the same time complexity as the computation of $A\ast B$, both for Boolean matrices and anti-distance matrices, where $A$ and $B$ are square matrices of the same size.

Note that sub-cubic multiplication techniques such as Strassen’s algorithm cannot be applied directly here, because the addition operator ($\vee$ or $\wedge$) has no inverse. This can be counteracted, leading to algorithms which may have better asymptotic behavior, but which may not be very fast in practice, especially if the matrices are not extemely large.

For arithmetic with Boolean matrices, the first optimization is due to the fact that $64$ bits fit into a register of a modern CPU. Namely, we can do arithmetic on $64$ matrix entries in parallel. In matBool.h, matrix operations are implemented, where each matrix row is divided into blocks of $64$ bits each. Below is our code for the matrix multiplication, as ${\tt mr}:=({\tt*this})\ast{\tt m}$:

template <int R, int C>template <int CC>matBool<R,CC> matBool<R,C>::operator * (const matBool<C,CC> &m) const{    matBool<R,CC> mr = zeroMatrix;    for (int c=C; --c>=0; ) {        blockType *m_c = m.block + c * m.Cblocks; // row c of m        for (int r=R; --r>=0; ) {            bool e = ((unsigned char *)(block + r * Cblocks))[c>>3] &                        (1 << (c & 7)); // entry r,c of *this            if (e) {                blockType *mr_r = mr.block + r * mr.Cblocks; // row r of mr                for (int ccb=m.Cblocks; --ccb>=0; ) {                    mr_r[ccb] |= m_c[ccb]; // block ccb of row c of m                }            }        }    }    return mr;}The code for the transitive closure is similar.

A subsequent optimization for multiplication is to use the “four Russians”-algorithm. See [3] or [1, pp. 243-247]. We did not implement this optimization. The “four Russians”-algorithm cannot be applied on the transitive closure algorithm.

In our implementation of the anti-distance matrices, we did not use $\mathcal{D}=\{0,\tfrac{1}{S},\tfrac{2}{S},\ldots,1\}$, but we used $S\,\mathcal{D}=\{0,1,2,\ldots,S\}$. This is to be compatible with standard unsigned integer types. We made implementations for three values of $S$, namely $2^{8}-1$, $2^{16}-1$, and $2^{32}-1$, corresponding to the standard integer types unsigned char, unsigned short, and unsigned int respectively.

We use template parameter T, being one of unsigned char, unsigned short, and unsigned int, to select the corresponding value of $S$. In matAntidist.h, matrix operations are implemented without optimization by way of vector instructions. Below is our code for the matrix multiplication, again as ${\tt mr}:=({\tt*this})\ast{\tt m}$:

template <int R, int C, class T>template <int CC>matAntidist<R,CC,T> matAntidist<R,C,T>::operator * (const matAntidist<C,CC,T> &m) const{    matAntidist<R,CC,T> mr = zeroMatrix;    for (int c=C; --c>=0; ) {        T *m_c = m.entry + c * CC; // row c of m        for (int r=R; --r>=0; ) {            T e = (entry + r * C)[c]; // entry r,c of *this            if (e) {                T *mr_r = mr.entry + r * CC; // row r of mr                e = ~e;                for (int cc=CC; --cc>=0; ) {                    T ee = m_c[cc]; // entry c,cc of m                    if (ee > e) {                        ee -= e;                        if (mr_r[cc] < ee) mr_r[cc] = ee;                    }                }            }        }    }    return mr;}Notice that $a+b-1$ in the definition of $\ast$ is scaled to $a+b-S$, which is computed as $b-(S-a)$ to reduce the number of operations. The code for the transitive closure is again similar.

In the optimization with vector instructions, we use Intel’s SSE instructions. The SSE instructions on Intel compatible hardware are the successor of MMX instructions. The instructions act on registers of 128 bits. The unsigned variants of integer type SSE instructions interpret such a register as a vector of elements of $S\,\mathcal{D}=\{0,1,2,\ldots,S\}$, with as many coordinates as the register size permits.

[TABLE]

In matAntidist_sse.h, matrix operations are implemented with optimization by way of vector instructions. Below is our code for the matrix multiplication, again as ${\tt mr}:=({\tt*this})\ast{\tt m}$:

template <int R, int C, class T>template <int CC>matAntidist<R,CC,T> matAntidist<R,C,T>::operator * (const matAntidist<C,CC,T> &m) const{    matAntidist<R,CC,T> mr = zeroMatrix;    for (int c=C; --c>=0; ) {        block_sse<T> *m_c = m.block + c * m.Cblocks; // row c of m        for (int r=R; --r>=0; ) {            T e = ((T *) (block + r * Cblocks))[c]; // entry r,c of *this            if (e) {                block_sse<T> *mr_r = mr.block + r * mr.Cblocks; // row r of mr                block_sse<T> b;                b.clonenot (e);                for (int ccb=m.Cblocks; --ccb>=0; ) {                    block_sse<T> bb = m_c[ccb]; // block ccb of row c of m                    bb.subsat (b);                    mr_r[ccb].max (bb);                }            }        }    }    return mr;}The above code uses three functions with SSE instructions, namely clonenot, subsat, and max. The code for the transitive closure is again similar.

In clonenot, the argument is negated logically and copied to all coordinates of the vector, which corresponds to the SSE register. In subsat, two SSE registers are interpreted as vectors, and the second vector is subtracted from the first with saturation. In max, again two SSE registers are interpreted as vectors, and the first vector is replaced by the second vector on spots where the second vector is larger.

Below follows the code for clonenot, subsat, and max with template argument T = unsigned char.

template <>inline void block_sse<unsigned char>::clonenot (unsigned c){    // clone from byte c to doubleword c    c *= 0x01010101;    // negate doubleword c logically    c = ~c;    // clone from doubleword c to block this    asm (        "movd %1, %0\n\t"        "pshufd $0, %0, %0"        : "=x" (val)        : "r" (c)    );}template <>inline void block_sse<unsigned char>::subsat (const block_sse<unsigned char> &b){    asm (        "psubusb %1, %0"        : "+x" (val)        : "xm" (b.val)    );}template <>inline void block_sse<unsigned char>::max (const block_sse<unsigned char> &b){    asm (        "pmaxub %1, %0"        : "+x" (val)        : "xm" (b.val)    );}In [4], the authors use distance matrices instead of anti-distance matrices, so their cloning function does not need negation. Furthermore, in their code for unsigned short matrix entries, the authors use punpck instructions to duplicate the number of coordinates with the clone value, until the whole SSE register is filled, which is after three such instructions.

For template argument T = unsigned int, there is no instruction for subtraction with saturation. To overcome this, we use the formula $a:=(a\vee b)-b$ for subtraction with saturation.

template <>inline void block_sse<unsigned int>::subsat (const block_sse<unsigned int> &b){    // psubusd does not exist, so one extra instruction is needed    // using max and regular subtraction    asm (        "pmaxud %1, %0\n\t"        "psubd %1, %0"        : "+x" (val)        : "xm" (b.val)    );}For addition of unsigned integers with saturation, it is recommended to use the formula $a:=(a\wedge\neg b)+b$, if $\neg b=S-b$ can be precomputed. But a more funny addition of four unsigned integers with saturation is given below.

template <>inline void block_sse<unsigned int>::addsat (const block_sse<unsigned int> &b){    // paddusd does not exist, so extra instructions are needed    // using regular addition and fix afterwards    asm (        "movdqa %1, %%xmm0\n\t"        "paddd %1, %0\n\t"        "pcmpgtd %0, %%xmm0\n\t"        "por %%xmm0, %0"        : "+x" (val)        : "xm" (b.val)    );}

Using distance matrices instead of anti-distance matrices is possible in combination with saturation as well, namely by way of $\wedge$ instead of $\vee$, and addition with saturation instead of subtraction with saturation. Any anti-distance matrix operation and its corresponding distance matrix operation satisfy the rules of De Morgan, e.g.

[TABLE]

for anti-distance matrix addition $\vee$ and distance matrix addition $\wedge$, and

[TABLE]

for anti-distance matrix multiplication $\ast$ and distance matrix multiplication $+$.

5 Conclusion

Most of our ideas were already discovered earlier, namely in [4], [1, pp. 200-206], and [2]. The most important new idea is to use saturation instead of signed edge values.

The saturation idea came in a natural way from the Boolean matrices as adjacency matrices. Due to saturation, we can deal with graphs in which the distance between vertices may be infinite, without the need of using floating point arithmetic.

Inspired by the Boolean matrices as adjacency matrices, we used anti-distance matrices instead of distance matrices. Since others used distance matrices, we added operations for distance matrices as well, which are $+$ (instead of $\ast$) for “funny matrix multiplication”, and dransclose and dransclosure (instead of transclose and transclosure) for computing the distance matrix with minimum loop lengths on the diagonal.

With signed matrix entries, saturation is not possible. This is not a problem if the matrix entry type is float, since saturation does not make any sense in that case. Furthermore, the correspondence with the Boolean case cannot be maintained in a natural way, so distance matrices are recommended with signed matrix entries.

Our implementation of the matrix multiplication and the Floyd-Warshall algorithm is not very suitable for very large matrices, since one has to take caching into account for such matrices. The authors of [4] do take caching into account. For even larger matrices, it is faster to make use of GPU shaders, see [6].

We used g++ template headers with inline assembly instructions, which can easily be added to projects which make use of the matrix operations. In addition to matrix multiplication and the Floyd-Warshall algorithm, we implemented the following entrywise matrix operations.

[TABLE]