Floyd-Warshall-Roy, 2

Posted on 2010/02/20 by wdjoyner

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^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^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^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!

Floyd-Warshall-Roy

Posted on 2010/02/08 by wdjoyner

The Floyd-Warshall-Roy algorithm is an algorithm for finding shortest paths in a weighted, directed graph. It allows for negative edge weights and detects a negative weight cycle if one exists. Assuming that there are no negative weight cycles, a single execution of the FWR algorithm will find the shortest paths between all pairs of vertices.

This algorithm is an example of dynamic programming, which allows one to break the computation down to simpler steps using some sort of recursive procedure. If $n=|V|$ then this is a $O(n^3)$ -time algorithm. For comparison, the Bellman-Ford algorithm has complexity $O(|V|\cdot |E|)$ , which is $O(n^3)$ -time for “dense” graphs. However, Bellman-Ford only yields the shortest paths eminating from a single vertex. To achieve comparable output, we would need to iterate Bellman-Ford over all vertices, which would be an $O(n^4)$ -time algorithm for “dense” graphs. Except for “sparse” graphs, Floyd-Warshall-Roy is better than an interated implementation of Bellman-Ford.

Here is a Sage implementation

def floyd_warshall_roy(A):
    """
    Shortest paths

    INPUT: 
        A - weighted adjacency matrix 

    OUTPUT
        dist  - a matrix of distances of shortest paths
        paths - a matrix of shortest paths

    EXAMPLES:
        sage: A = matrix([[0,1,2,3],[0,0,2,1],[20,10,0,3],[11,12,13,0]]); A
        sage: floyd_warshall_roy(A)
        ([[0, 1, 2, 2], [12, 0, 2, 1], [14, 10, 0, 3], [11, 12, 13, 0]], 
          [-1  1  2  1]
          [ 3 -1  2  3]
          [ 3 -1 -1  3]
          [-1 -1 -1 -1])
        sage: A = matrix([[0,1,2,4],[0,0,2,1],[0,0,0,5],[0,0,0,0]])
        sage: floyd_warshall_roy(A)
        ([[0, 1, 2, 2], [+Infinity, 0, 2, 1], [+Infinity, +Infinity, 0, 5],
          [+Infinity, +Infinity, +Infinity, 0]], 
          [-1  1  2  1]
          [-1 -1  2  3]
          [-1 -1 -1  3]
          [-1 -1 -1 -1])
        sage: A = matrix([[0,1,2,3],[0,0,2,1],[-5,0,0,3],[1,0,1,0]]); A
        sage: floyd_warshall_roy(A)
        Traceback (click to the left of this block for traceback)
        ...
        ValueError: A negative edge weight cycle exists.
        sage: A = matrix([[0,1,2,3],[0,0,2,1],[-1/2,0,0,3],[1,0,1,0]]); A
        sage: floyd_warshall_roy(A)
        ([[0, 1, 2, 2], [3/2, 0, 2, 1], [-1/2, 1/2, 0, 3/2], [1/2, 3/2, 1, 0]],
          [-1  1  2  1]
          [ 2 -1  2  3]
          [-1  0 -1  1]
          [ 2  2 -1 -1])
    """
    G = Graph(A, format="weighted_adjacency_matrix")
    V = G.vertices()
    E = [(e[0],e[1]) for e in G.edges()]
    n = len(V)
    dist = [[0]*n for i in range(n)]
    paths = [[-1]*n for i in range(n)]
    # initialization step
    for i in range(n):
        for j in range(n):
            if (i,j) in E:
                paths[i][j] = j
            if i == j:
                dist[i][j] = 0
            elif A[i][j] != 0:
                dist[i][j] = A[i][j]
            else:
                dist[i][j] = infinity
    # iteratively finding the shortest path
    for j in range(n):
        for i in range(n):
            if i != j:
                for k in range(n):
                    if k != j:
                        if dist[i][k]>dist[i][j]+dist[j][k]:
                            paths[i][k] = V[j]
                        dist[i][k] = min(dist[i][k], dist[i][j] +dist[j][k])
    for i in range(n):
        if dist[i][i] < 0:
            raise ValueError, "A negative edge weight cycle exists."
    return dist, matrix(paths)

Sage and the future of mathematics

Posted on 2010/02/04 by wdjoyner

I am not a biologist nor a chemist, and maybe neither are you, but suppose we were. Now, if I described a procedure, in “standard” detail, to produce a result XYZ, you would (based on your reasonably expected expertise in the field), follow the steps you believe were described and either reproduce XYZ or, if I was mistaken, not be able to reproduce XYZ. This is called scientific reproducibility. It is cructial to what I believe is one of the fundamental principles of science, namely Popper’s Falsifiability Criterion.

More and more people are arguing, correctly in my opinion, that in the computational realm, in particular in mathematical research which uses computational experiments, that much higher standards are needed. The Ars Technica article linked above suggests that “it’s time for a major revision of the scientific method.” Also, Victoria Stodden argues one must “facilitate reproducibility. Without this data may be open, but will remain de facto in the ivory tower.” The argument basically is that to reproduce computational mathematical experiments is unreasonably difficult, requiring more that a reasonable expertise. These days, it may in fact (unfortunately) require purchasing very expensive software, or possessing of very sophisticated programming skills, accessibility to special hardware, or (worse) guessing parameters and programming procedures only hinted at by the researcher.

Hopefully, Sage can play the role of a standard bearer for such computational reproducibility. Sage is free, open source and there is a publically available server it runs on (sagenb.org).

What government agencies should require such reproducibility? In my opinion, all scientific funding agencies (NSF, etc) should follow these higher standards of computational accountability.

Bases of representable/vector matroids in Sage

Posted on 2010/02/01 by wdjoyner

A matroid M is specified by a set X of “points” and a collection B of “bases, where $X \subset 2^X$ , satisfying certain conditions (see, for example, Oxley’s book, Matroid theory).

One way to construct a matroid, as it turns out, is to consider a matrix M with coefficients in a field F. The column vectors of M are indexed by the numbers $J = \{1,2,\dots,n\}$ . The collection B of subsets of J associated to the linearly independent column vectors of M of cardinality $r = rank_F(M)$ forms a set of bases of a matroid having points $X=J$ .

How do you compute B? The following Sage code should work.



def independent_sets(mat):
    """
    Finds all independent sets in a vector matroid.

    EXAMPLES:
        sage: M = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])
        sage: M
        [1 0 0 1 1 0]
        [0 1 0 1 0 1]
        [0 0 1 0 1 1]
        sage: independent_sets(M)  
        [[0, 1, 2], [], [0], [1], [0, 1], [2], [0, 2], [1, 2], 
        [0, 1, 4], [4], [0, 4], [1, 4], [0, 1, 5], [5], [0, 5], 
        [1, 5], [0, 2, 3], [3], [0, 3], [2, 3], [0, 2, 5], [2, 5], 
        [0, 3, 4], [3, 4], [0, 3, 5], [3, 5], [0, 4, 5], [4, 5], 
        [1, 2, 3], [1, 3], [1, 2, 4], [2, 4], [1, 3, 4], [1, 3, 5], 
        [1, 4, 5], [2, 3, 4], [2, 3, 5], [2, 4, 5]]
        sage: M = matrix(GF(3), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])
        sage: independent_sets(M)
        [[0, 1, 2], [], [0], [1], [0, 1], [2], [0, 2], [1, 2], 
        [0, 1, 4], [4], [0, 4], [1, 4], [0, 1, 5], [5], [0, 5], 
        [1, 5], [0, 2, 3], [3], [0, 3], [2, 3], [0, 2, 5], [2, 5], 
        [0, 3, 4], [3, 4], [0, 3, 5], [3, 5], [0, 4, 5], [4, 5], 
        [1, 2, 3], [1, 3], [1, 2, 4], [2, 4], [1, 3, 4], [1, 3, 5], 
        [1, 4, 5], [2, 3, 4], [2, 3, 5], [2, 4, 5], [3, 4, 5]]
        sage: M345 = matrix(GF(2), [[1,1,0],[1,0,1],[0,1,1]])
        sage: M345.rank()
        2
        sage: M345 = matrix(GF(3), [[1,1,0],[1,0,1],[0,1,1]])
        sage: M345.rank()
        3
    """
    F = mat.base_ring()
    n = len(mat.columns())
    k = len(mat.rows())
    J = Combinations(n,k)
    indE = []
    for x in J:
        M = matrix([mat.column(x[0]),mat.column(x[1]),mat.column(x[2])])
        if k == M.rank():     # all indep sets of max size
            indE.append(x)
            for y in powerset(x): # all smaller indep sets
                if not(y in indE):
                    indE.append(y)
    return indE

def bases(mat, invariant = False):
    """
    Find all bases in a vector/representable matroid.

    EXAMPLES:
        sage: M = matrix(GF(2), [[1,0,0,1,1,0],[0,1,0,1,0,1],[0,0,1,0,1,1]])
        sage: bases(M)

        [[0, 1, 2],
         [0, 1, 4],
         [0, 1, 5],
         [0, 2, 3],
         [0, 2, 5],
         [0, 3, 4],
         [0, 3, 5],
         [0, 4, 5],
         [1, 2, 3],
         [1, 2, 4],
         [1, 3, 4],
         [1, 3, 5],
         [1, 4, 5],
         [2, 3, 4],
         [2, 3, 5],
         [2, 4, 5]]
        sage: L = [x[1] for x in bases(M, invariant=True)]; L.sort(); L
         [5, 16, 17, 33, 37, 44, 45, 48, 81, 84, 92, 93, 96, 112, 113, 116]
    """
    r = mat.rank()
    ind = independent_sets(mat)
    B = []
    for x in ind:
        if len(x) == r:
            B.append(x)
    B.sort()
    if invariant == True:
        B = [(b,sum([k*2^(k-1) for k in b])) for b in B]
    return B

M	T	W	T	F	S	S
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28

Yet Another Mathblog

Monthly Archives: February 2010

Floyd-Warshall-Roy, 2

Floyd-Warshall-Roy

Sage and the future of mathematics

Bases of representable/vector matroids in Sage