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

arXiv:1705.08743v2 [cs.DM] 7 Nov 2017

Michiel de Bondt
November 8, 2017
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.

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.
I ∨ A ∨ (A ∧ A) ∨ (A ∧ A ∧ A) ∨ · · ·

(1)

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.
A ∨ (A ∧ A) ∨ (A ∧ A ∧ A) ∨ · · ·

(2)

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
1

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 FloydWarshall-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 ∗ 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 B = {0, 1} denote the Boolean set in numerical fashion. We define an
addition on B, namely
∨:B×B →B

by

a ∨ b = max{a, b}

Furthermore, we define a multiplication on B, namely
∧:B×B →B

by

a ∧ b = min{a, b}

Using this addition and multiplication, we can define matrix addition and multiplication over B, namely in the same way as usual. So if we take A, B ∈
MatR,C (B), then
(A ∨ B)ij = Aij ∨ Bij

(1 ≤ i ≤ R, 1 ≤ j ≤ C)

If we take A ∈ MatR,C (B) and B ∈ MatC,C ′ (B), then
(A ∗ B)ij = (Ai1 ∧ B1j ) ∨ · · · ∨ (AiC ∧ BCj )

(1 ≤ i ≤ R, 1 ≤ j ≤ C ′ )

Let Alk denote the k-th column of A and Bk↔ denote the k-th row of B.
Write 0R×C for the zero matrix with R rows and C columns. We can compute
P = A ∗ B by way of the following code:
P := 0R×C
for k := 1 to C do
2

P := P ∨ (Alk ∗ Bk↔ )
end for
return P
We can interpret a square Boolean matrix, say A ∈ MatD,D (B), as the
adjacency matrix of a directed graph. Let B ∈ MatD,D (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 ∨ B is adjacency matrix of
the union of A and B. The zero matrix 0D×D is the neutral element of ∨, and
the interpretation of 0D×D is the empty graph.
The interpretation of A ∗ 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 ID of size
D is the neutral element of ∗, and the interpretation of ID 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
A ∨ (A ∗ A) ∨ (A ∗ A ∗ A) ∨ · · ·
(3)
which is actually just
′

′

A ∨ (A ∗ A) ∨ (A ∗ A ∗ A) ∨ · · · ∨ A∗D = A ∗ (ID′ ∨ A)∗(D −1)
′

for every D′ ≥ D, where A∗D is the product of D′ copies of A.
′
(ID′ ∨ A)∗(D −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 ∨ (Tlk ∗ Tk↔ )
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 ∨ B, 0D×D , A ∗ B,
and ID as much as possible. For that reason, a 1 indicates distance 0, and a
0 indicates distance ∞. Any number in between 0 and 1 indicates a distance
between 0 and ∞, i.e. a proper distance.
So let us replace the Boolean set B by the interval C = [0, 1]. Let δ : C →
[0, ∞] denote the distance associated to an element of C. Then δ(0) = ∞ and
δ(1) = 0, so it is natural to impose that δ is decreasing.
The addition operator on C should take the minimum distance, which is
accomplished by taking the maximum value, because δ is decreasing. So we can
take ∨ as the addition operator for C, just as for B.
3

The multiplication operator on C should add distances, but ∧ only does this
for the distances 0 and ∞, which are in fact the only distances in the Boolean
case. So we cannot preserve ∧. No, we use another multiplication on C, namely
we define
∗:C×C →C
by
a ∗ b = δ −1 (δ(a) + δ(b))
If we extend ∧ to [0, ∞], then

(a ∨ b) ∗ c = δ −1 δ(a ∨ b) + δ(c)


= δ −1 δ(a) ∧ δ(b) + δ(c)


= δ −1 δ(a) + δ(c) ∧ δ(b) + δ(c)


= δ −1 δ(a) + δ(c) ∨ δ −1 δ(b) + δ(c)
= (a ∗ c) ∨ (b ∗ c)

and similarly c ∗ (a ∨ b) = (c ∗ a) ∨ (c ∗ b), so the distributive law is fulfilled.
From this, the distributive law for matrices follows.
We can take δ(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, ∞) for some positive number S, or we are just
not interested in distinguishing distances in the interval [S, ∞], we can take
δ(a) = S · (1 − a).
We can go even further. Since continuous intervals are not so convenient
in computer practice, we can replace C = [0, 1] by D = {0, S1 , S2 , . . . , 1}. Then
δ(D) = {0, 1, 2, . . . , S}. We call S is the saturation distance, and define
∗:D×D →D

by

a ∗ b = max{a + b − 1, 0}

which is just as above with δ(a) = S · (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 δ(C) = [0, ∞] instead of C, with ∧ 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


0 ∞ ··· ∞

. 
..
 ∞ 0
. .. 



 . .
.. ... ∞ 
 ..
∞ ··· ∞ 0
The neutral matrix for addition by way of ∧ is the matrix of which all entries
are ∞ 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 ∗ B, both for Boolean matrices and antidistance matrices, where A and B are square matrices of the same size.

4

Note that sub-cubic multiplication techniques such as Strassen’s algorithm
cannot be applied directly here, because the addition operator (∨ or ∧) 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 mr := (∗this) ∗ 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 D =
{0, S1 , S2 , . . . , 1}, but we used S D = {0, 1, 2, . . . , S}. This is to be compatible
with standard unsigned integer types. We made implementations for three values of S, namely 28 −1, 216 −1, and 232 −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 mr := (∗this)∗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; ) {

5

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 ∗ 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 D = {0, 1, 2, . . . , S}, with as many coordinates as the register size permits.
S

integer type

vector size

28 − 1

unsigned char

16

unsigned short

8

unsigned int

4

2

16

−1

232 − 1

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 mr := (∗this) ∗ 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

6

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)
);

7

}

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 ∨ 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 ∧ ¬b) + b, if ¬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 ∧ instead of ∨, 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.
¬(A ∨ B) = (¬A) ∧ (¬B)

¬(A ∧ B) = (¬A) ∨ (¬B)

for anti-distance matrix addition ∨ and distance matrix addition ∧, and
¬(A ∗ B) = (¬A) + (¬B)

¬(A + B) = (¬A) ∗ (¬B)

for anti-distance matrix multiplication ∗ and distance matrix multiplication +.
8

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 ∗) 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.
mathematical
operator

c++
operator

entrywise
computation

∧

&

min{a, b}

∨

|

max{a, b}

⊻

^

|a − b|

¬

~

S−a

References
[1] Alfred V. Aho, John E. Hopcroft, and Jeffrey D. Ullman. The design and
analysis of computer algorithms. Addison-Wesley Publishing Co., Reading,
Mass.-London-Amsterdam, 1975. Second printing, Addison-Wesley Series in
Computer Science and Information Processing.
[2] Noga Alon, Zvi Galil, and Oded Margalit. On the exponent of the all pairs
shortest path problem. J. Comput. System Sci., 54(2, part 1):255–262, 1997.
32nd Annual Symposium on Foundations of Computer Science (San Juan,
PR, 1991).
9

[3] V. L. Arlazarov, E. A. Dinic, M. A. Kronrod, and I. A. Faradžev. On
economical construction of the transitive closure of a directed graph. Soviet
Mathematics—Doklady, 11(5):1209–1210, 1970.
[4] Sung-Chul Han, Franz Franchetti, and Markus Püschel. Program generation
for the all-pairs shortest path problem. In Proceedings of the 15th International Conference on Parallel Architectures and Compilation Techniques,
PACT ’06, pages 222–232, New York, NY, USA, 2006. ACM.
[5] Pawan Harish and P. J. Narayanan. Accelerating large graph algorithms
on the GPU using CUDA. In High Performance Computing – HiPC 2007,
pages 197–208. Springer Berlin Heidelberg.
[6] Gary J. Katz and Joseph T. Kider, Jr. All-pairs shortest-paths for
large graphs on the gpu.
In Proceedings of the 23rd ACM SIGGRAPH/EUROGRAPHICS Symposium on Graphics Hardware, GH ’08,
pages 47–55, Aire-la-Ville, Switzerland, Switzerland, 2008. Eurographics Association.

10

