---
author:
- 'Piotr Sankowski[ [^1] ]{}'
- 'Karol Węgrzycki[^2]'
bibliography:
- 'bib.bib'
---

Introduction
============

The *All-Pairs Shortest Paths* (APSP) problem asks to find distances between all pairs of vertices in a graph. For a directed graphs with weights in $\mathbb{R}$, there is a classical $O(n^3)$ time algorithm @floyd [@warshall]. Currently best upper bound for this problem is due to @williams $O(\frac{n^3}{2^{\Omega(\log n)^{0.5}}})$ algorithm. It is asymptotically faster than $O(n^3/\log^c n)$ for any $c>0$ (see survey [@chan07] for earlier algorithms). Showing any algorithm that would work in $O(n^{3-\epsilon})$ time for some $\epsilon > 0$ is a major open problem [@williams].

If we consider unweighted, directed graphs there are subcubic algorithms that exploit fast matrix multiplication. For the undirected graph @seidel presented the optimal ${\widetilde{O}}(n^\omega)$ time algorithm, where $\omega < 2.373$ is the matrix multiplication exponent [@legall]. For the directed case @zwick presented an $O(n^{2.575})$ time algorithm that is based on the fast rectangular matrix multiplication. Moreover, if we are interested in small integer weights from the set $\{-M,\ldots,M\}$ we have $O(M^{0.68}n^{2.575})$ algorithm [@zwick].

Because APSP in undirected graphs can be solved in ${\widetilde{O}}(n^\omega)$, diameter, radius, shortest cycle, etc. can be determined in ${\widetilde{O}}(n^\omega)$ time as well. It is surprising that for a directed case, where merely $O(n^{2.575})$ APSP is known there are also ${\widetilde{O}}(n^\omega)$ algorithms for determining these properties. After a long line of improvements @cygan showed an ${\widetilde{O}}(M n^\omega)$ time algorithms for finding minimum weight perfect matching, shortest cycle, diameter and radius. Also, they showed an application of their techniques that improves upon @cycle-each-vertex ${\widetilde{O}}(M n^\omega t)$ time algorithm for the following problem: *determine the set of vertices that lie on some cycle of length at most $t$*. @cygan managed to solve this problem in ${\widetilde{O}}(M n^\omega)$ time using Baur-Strassen’s theorem.

All of these algorithms are effective only in the case of a dense graphs. For graphs with the small number of edges there are better algorithms (e.g., APSP in $\tilde{O}(|V||E|)$ time [@thorup]). But these algorithms are $\Theta(n^3)$ when $|E|
= \Theta(n^2)$.

Related Work
============

Distance Queries
----------------

@distance-queries considered the weighted, directed graphs with weights in $\{-M,\ldots,M\}$. They showed an algorithm that needs ${\widetilde{O}}(M n^\omega)$ preprocessing time. After preprocessing each distance $\delta (u,v)$ in the graph can be computed exactly in $O(n)$ query time. In the special case $M=1$ they showed ${\widetilde{O}}(n^\omega)$ algorithm that solves *Single Source Shortest Paths* (SSSP).

We will match their bounds (up to the polylogarithmic factors) using Frobenius normal form. Next we will extend their algorithm so it will return more information about a graph in the same query/preprocessing time.

Counting Cycles
---------------

For a given graph $G$ determining whether $G$ contains a simple cycle of length exactly $k$ is NP-hard (in particular determining whether a graph contains a Hamiltonian cycle is NP-complete). However, if we fix $k$ to be a constant this problem can be solved in polynomial time.

@color-coding introduced a color coding technique. For a fixed $k$ if a graph $G(V,E)$ contains a simple cycle of size exactly $k$ then such cycle can be found in ${\widetilde{O}}(|V|^\omega)$ time. Unfortunately, their algorithm depends exponentially $2^{O(k)}$ on the length of the cycle and in consequence is inapplicable for large $k$. In the next years, @finding-cycles showed (using a different technique) that for $k \le 7$ it is possible to count the number of cycles of length exactly $k$ in a graph in ${\widetilde{O}}(|V|^\omega)$ time. In [@even-cycles] it is shown that for any even $k$, cycles of length $k$ can be found in $O(|V|^2)$ time in undirected graphs (if they contain such a cycle). @finding-cycles showed more methods that depend solely on a number of edges in a graph. For example for odd $k$ they showed $O(E^{2-\frac{2}{k+1}})$ algorithm for finding a cycles of length $k$. However, for dense graphs these results are worse than @color-coding.

It appears that to break the exponential dependence on the length of the cycle we can do one of the following:

-   Consider non-simple cycles (the vertices can reoccur) of length exactly $k$,

-   Determine cycles of length at most $k$.

To detect whether a non-simple cycle of length exactly $k$ exists one can use the folklore algorithm. It starts by taking the adjacency matrix $A$ of a graph $G$. Subsequently, in ${\widetilde{O}}(n^\omega)$ time compute $A^k$ by repeated squaring. If $\mathrm{Tr}\left[A^k\right] > 0$ then a non-simple $k$ length cycle exists [^3].

@cycle-each-vertex considered the following problem: *for every vertex in a graph find a shortest cycle that contains it*. He called this problem *All-Nodes Shortest Cycle* (ANSC). He showed a randomized algorithm that solves ANSC for undirected graphs with weights $\{1,\ldots, M\}$ in ${\widetilde{O}}(\sqrt{M} n^{(\omega + 3)/2})$ time. He noted that for simple digraphs (directed graphs with no anti-parallel edges) it is reduced to All-Pairs Shortest Paths problem. The fastest known APSP algorithm for unweighted, directed graphs runs in $O(n^{2.575})$ due to [@zwick]. Here, we will show how to solve ANSC in ${\widetilde{O}}(n^\omega)$ for general, unweighted, directed graphs. Unfortunately, our techniques will allow us only to find the length of such a cycle (not determining it). But we can return the set of points, that lie on some cycle of a given length. Independently to our work @agarwal proved that ANSC can be solved in ${\widetilde{O}}(n^\omega)$ for unweighted, undirected graphs using a completely different technique.

@cycle-each-vertex also considered following problem: *given a graph and an integer $t$. Let $S(k)$ denote the set of all vertices lying in a cycle of length $\le k$. Determine $S(t)$*. He considered directed graphs with weights in $\{ -M, \ldots, M\}$ and showed ${\widetilde{O}}(M n^\omega
t)$ algorithm

Recently, @cygan improved his algorithm. They showed that for a fixed $t
\in [0,n M]$ the set $S(t)$ can be computed in ${\widetilde{O}}(M n^\omega)$ randomized time. We show, that for an unweighted ($M=1$) directed graphs we can compute sets $S(1),S(2),\ldots,S(D)$ in ${\widetilde{O}}(n^\omega)$ time with high probability.

Preliminaries {#prelim}
=============

Let $T(n)$ be the minimal number of algebraic operations needed to compute the product of $n\times n$ matrix by an $n \times n$ matrix. We say that $\omega$ is the exponent of square matrix multiplication. For now the best known upper bound on $\omega$ is due to @legall:

For every $\epsilon>0$, $T(n) < O(n^{\omega +
    \epsilon})$ where $\omega < 2.37287$.

In this paper, we will omit $\epsilon$ in definition and will assume (like in many other papers) that $O(n^\omega)$ operations are needed to multiply two matrices. The best lower bound for the exponent of matrix multiplication is $\omega \ge 2$. For convenience in this paper we will assume that $\omega >
2$. The $\tilde{O}$ notation hides polylogarithmic factors in the complexity. We will use it to emphasize that all our algorithms need the polylogarithmic number of calls to the fast matrix multiplication algorithm.

The Frobenius canonical-form of a matrix can be computed deterministically using ${\widetilde{O}}(n^\omega)$ field operations.

The comprehensive description of Frobenius normal form will be presented in Section \[frobenius\]. The properties of Frobenius normal form used in this paper are well known in literature [@frobenius; @storjohann1; @dummit]. There are also probabilistic algorithms that compute this form in expected ${\widetilde{O}}(n^\omega)$ with high probability over small fields [@eberly]. In this paper, all algorithms are deterministic if we set the upper bound on the number of distinct walks $W$. Then, due to the time of a single field operation we need additional $O(\log{W})$ factor in the complexity. However, if we are only interested if cycle/walk of a given length exists in a graph, we can set sufficiently small field $\mathbb{Z}_p$ ($p$ has $O(\log{n})$ bits). This way when algorithm returns nonzero we are sure that some walk exists there. If algorithm returns zero, then with high probability there is no such walk.

In this paper, randomization occurs only because for some graphs number of walks can be exponential (e.g., $W = O(2^n)$) and the output requires to return them.

For matrices $A \in \mathbb{R}^{n\times k}$ and $B \in \mathbb{R}^{n \times l}$ the $A \oplus B \in \mathbb{R}^{n \times (k+l)}$ denotes the concatenation of their columns. $C_{a,b} \in \mathbb{R}^{n\times (b-a)}$ denotes a matrix constructed by concatenating columns $c_a,c_{a+1},\ldots,c_b$ of the matrix $C \in \mathbb{R}^{n
\times m}$.

Consequences of Frobenius Normal Form {#frobenius}
=====================================

Let $K$ be a commutative field. For any matrix $A \in K^{n\times n}$ there exists an invertible $U$ over $K$ such that:

$$U^{-1} A U = F = 
    \begin{bmatrix}
      C_1 &     &     &        & 0   \\
          & C_2 &     &              \\
          &     & C_3 &        &     \\
          &     &     & \ddots &     \\
        0 &     &     &        & C_k \\
    \end{bmatrix}
    .$$

and $F$ is the Frobenius-canonical-form [^4] of $A$. The diagonal block $C_i$ is called the companion matrix:

$$C_i = 
    \begin{bmatrix}
        0& &\ldots & &0      & -c_0\\
        1&0&       & &0      & -c_1\\
         &1&\ddots & &\vdots & -c_2\\
         & &\ddots &0&       & \vdots \\
         & &       &1&0      & -c_{r-2}\\
        0& &       & &1      & -c_{r-1}\\
    \end{bmatrix}
    \in K^{r\times r}
    .$$

Each companion matrix corresponds to the monic polynomial $C_i(x) = x^r +
c_{r-1}x^{r-1} + \ldots + c_0 \in {K\left[x\right]}$ and is called the *minimal polynomial* of $A$. Each minimal polynomial has a property that $C_i(A) = 0$. To guarantee that matrix has only one decomposition into Frobenius normal form we require that every polynomial must divide the next one, i.e., $C_i(x) |
C_{i+1}(x)$. The final list of polynomials is called the *invariant factors* of matrix $A$ [@frobenius].

Cyclic Subspaces
----------------

Frobenius decomposition can be used to get the desired power of a matrix (analogously to the diagonal decomposition):

$$A^k = (U F U^{-1})^k = U F (U^{-1} U) F \cdots F (U^{-1} U) F U^{-1} = U F^k U^{-1}.$$

Moreover, we will use the property that the power of block diagonal matrix $F$ is block diagonal:

$$F^k = 
    \begin{bmatrix}
      C_1^k &     &     &        & 0   \\
          & C_2^k &     &              \\
          &     & C_3^k &        &     \\
          &     &     & \ddots &     \\
        0 &     &     &        & C_l^k \\
    \end{bmatrix}
    .$$

Now, we need a property of companion matrices that will enable us to power them efficiently.

\[cyclic-property\] Let $v_1, \ldots, v_n$ be the columns of a matrix $C \in K^{n\times
    n}$. Let $v_{n+1},\ldots, v_{2n}$ be the columns of matrix $C^{n+1}$. If, for every $1 \le k \le n$ the columns of matrix $C^k$ are $v_k, v_{k+1}, \ldots,
    v_{k+n}$ then the $C$ has a **cyclic property**.

To further explore an idea of cyclic subspaces and get proper introduction to the theorem below, see Appendix \[appendix-cyclic\].

$$\left[\begin{array}{*9{c}cc}
    {\tikz[overlay,remember picture,baseline=(fuleft.base)] \node (fuleft) {};} & & & {\tikz[overlay,remember picture,baseline=(fkleft.base)] \node (fkleft) {};} & & & {\tikz[overlay,remember picture,baseline=(fnuleft.base)] \node (fnuleft) {};} & & {\tikz[overlay,remember picture,baseline=(fkuright.base)] \node (fkuright) {};}\\
     \\
    v_1 & v_2 & \ldots & v_k & \ldots & v_n & v_{n+1} & \ldots & v_{k+n-1} & \ldots & v_{2n} \\
    \\
    {\tikz[overlay,remember picture,baseline=(fleft.base)] \node (fleft) {};} & & & & & {\tikz[overlay,remember picture,baseline=(fright.base)] \node (fright) {};} & {\tikz[overlay,remember picture,baseline=(fnleft.base)] \node (fnleft) {};} & & {\tikz[overlay,remember picture,baseline=(fkright.base)] \node (fkright) {};} & & {\tikz[overlay,remember picture,baseline=(fnright.base)] \node (fnright) {};} \\
  \end{array}\right]
    \tikz[overlay,remember picture]{
        \draw[decoration={brace,mirror,raise=20pt},decorate] (fleft) -- node[below=23pt] {$C^1$} (fright);
        \draw[decoration={brace,mirror,raise=20pt},decorate] (fnleft) -- node[below=23pt] {$C^{n+1}$} (fnright);
        \draw[decoration={brace,raise=20pt},decorate] (fkleft) -- node[above=23pt] {$C^k$} (fkuright);
    }
    {
    \tikz[overlay,remember picture]{
    \node[highlight,fit=(fuleft.north west) (fright.south east)] () {};}
}
    {
    \tikz[overlay,remember picture]{
    \node[highlight,fit=(fnuleft.north west) (fnright.south east)] () {};}
}
    {
    \tikz[overlay,remember picture]{
    \node[highlight,fill=red!30, inner sep=12pt,fit=(fkleft.north west) (fkright.south east)] () {};}
}$$

\[theorem-cyclic\] Every companion matrix has a cyclic property.

For the numerical example of Theorem \[theorem-cyclic\] see Appendix \[appendix-example\].

Matching Distance Queries on Directed Unweighted Graphs
=======================================================

In this section, we will present a simple algorithm that matches the best known upper bounds of @distance-queries distance queries in directed unweighted graphs.

Answering Distance Queries by Using Frobenius Normal Form
---------------------------------------------------------

We take the adjacency matrix $A$ of a graph $G$ (i.e., $n \times n$ matrix with $a_{u,v} = 1$ when $(u,v) \in G$ and $0$ otherwise). The $k$-th power of the adjacency matrix of the graph $G$ holds the number of walks, i.e., an $a_{u,v}$ element of $A^k$ is *the count of distinct walks from $u$ to $v$ of length $k$ in the graph*. The shortest path between vertices $u,v$ is the smallest $k$ such that $A^k$ has nonzero element $a_{u,v}$. For a brief moment, we will forget about graph theory interpretation and focus only on finding such $k$.

We decompose matrix $A$ into the Frobenius normal form. @frobenius showed an algorithm that returns $U$ and $F$ deterministically in ${\widetilde{O}}(n^\omega)$ time (note that matrix inverse can also be computed in ${\widetilde{O}}(n^\omega)$ time).

Single Invariant Factor {#single-factor}
-----------------------

To better explain the idea, for a start we will consider a simple case when a number of invariant factors of $A$ is exactly $1$. In that situation, the matrix $F$ is a companion matrix $C \in K^{n\times n}$.

First, we compute the $(n+1)$-th power of the companion matrix $F^{n+1}$. This can be done in ${\widetilde{O}}(n^\omega)$ time by repeated squaring (compute $F, F^2,
F^4,\ldots, F^{n+1}$ with $O(\log{n})$ matrix multiplications) [^5].

If the matrix $U F$ has columns $v_1,\ldots, v_n$ and the matrix $U F^{n+1}$ has columns $v_{n+1},\ldots,v_{2n}$, then the columns $v_k,\ldots,v_{k+n-1}$ construct $U F^k$ (see Figure \[mult\]). It is because the matrix $F$ has the cyclic property.

$$\begin{array}{cc}
         & {
  \begin{array}{*9{c}ccccc}

          {\tikz[remember
picture,baseline=(f1a.base)] \node (f1a) {\ensuremath{v_1}};} & \ldots & {\tikz[remember
picture,baseline=(fka.base)] \node (fka) {\ensuremath{v_k}};} & \ldots &
          {\tikz[remember
picture,baseline=(f1b.base)] \node (f1b) {\ensuremath{v_n}};} & {\tikz[remember
picture,baseline=(fna.base)] \node (fna) {\ensuremath{v_{n+1}}};} & \ldots &
          {\tikz[remember
picture,baseline=(fkb.base)] \node (fkb) {\ensuremath{v_{k+n-1}}};} & \ldots & {\tikz[remember
picture,baseline=(fnb.base)] \node (fnb) {\ensuremath{v_{2n}}};} \\
  \end{array}
  {    \tikz[overlay,remember picture]{
    \node[highlight,,fit=(f1a) (f1b)] () {};}
}
  {    \tikz[overlay,remember picture]{
    \node[highlight,,fit=(fna) (fnb)] () {};}
}
  {    \tikz[overlay,remember picture]{
    \node[highlight,fill=red!15, below=1.6cm, minimum height=1.8cm,fit=(fka) (fkb)] () {$U F^k$};}
}
  \vspace{10pt}
}\\
        {
    U = \left[\begin{array}{ccc}
            & \; u_1    & \; \\
            & \; u_2    & \; \\
            & \; \vdots & \; \\
            & \; u_n    & \; 
        \end{array}\right]
}& \\
    \end{array}$$

This step takes just $2$ matrix multiplications, because we need to multiply $U$ times $F$ and $F^{n+1}$. The preprocessing phase takes only ${\widetilde{O}}(n^\omega)$ time.

Now, if a user wants to know the number of distinct walks from vertices $u$ to $v$ of length exactly $k$ he takes:

-   The $u$-th row of matrix $U F^k$ ($n$ numbers),

-   The $v$-th column of matrix $U^{-1}$,

-   Multiplies them in $O(n)$ time (dot product of $n$ dimensional vectors).

This will give us the $u,v$ element of matrix $U F^k U^{-1} = A^k$. To get the length of the shortest path (i.e., the minimal $k$ that $w_{u,v} > 0$), we will modify our matrices slightly to get the number of walks of length $\le k$. At the end, we will fit in $\tilde{O}(n)$ query time (by using binary search) and ${\widetilde{O}}(n^\omega)$ preprocessing time.

Basically, for a given $k$ we need to get the $u,v$ element of matrix $A + A^2 +
\cdots + A^k$. It suffices to add consecutive columns of matrix $U F \oplus
U F^{n+1} = v_1 \oplus v_2 \oplus \ldots \oplus v_{2n}$ in the following manner [^6]:

$$M' = \left[\begin{array}{ccccccc}
            v_1  & v_1 + v_2  & v_1 + v_2 + v_3  & \ldots &  \sum^k_{i=1}
            v_i &  \ldots & \sum^{2n}_{i=1} v_i \\
    \end{array}\right] \in \mathbb{R}^{n \times 2n}
    .$$

Now, to get $A + A^2 + \cdots + A^k$ one needs to multiply $M'_{k,k+n-1}
U^{-1}$ and subtract $M'_{1,n} U^{-1}$ for a balance [^7].

We can transform matrices $U$ and $F$ to matrix $M'$ in $O(n^2)$ time during preprocessing. During query, we need to compute 2 dot products ($u$-th row of $M'_{k,k+n-1}$ times $v$-th column of $U^{-1}$ and $u$-th row of $M'_{1,n}$ times $v$-th column of $U^{-1}$) and subtract them.

We have an algorithm that for a given vertices $u,v \in G$ and integer $k \in \{1,\ldots,n\}$ can answer: *how many walks from $u$ to $v$ of length less or equal $k$ are in the graph $G$* in $\tilde{O}(n)$ query and ${\widetilde{O}}(n^\omega)$ preprocessing.

Because the result of the query is increasing in $k$ we can use binary search. We can determine the first $k$ for which the query will answer nonzero value in $O(\log n)$ tries. Hence, in $\tilde{O}(n)$ we can find the length of the shortest path. This generalized query can easily return the number of walks of length exactly $k$, i.e., $q(u,v,k) - q(u,v,k-1)$.

For now, we only matched the result of @distance-queries for unweighted graphs with a single, invariant factor. In the next section, we will show how to generalize our technique for graphs with any number of invariant factors. Then, we will extend the @distance-queries algorithm. Namely, we will show that in ${\widetilde{O}}(n^\omega)$ preprocessing time and $\tilde{O}(n)$ query time we can get $D$ numbers (where $D$ is the graph diameter): $w_{u,v}^1,w_{u,v}^2,\ldots,w_{u,v}^D$. The number $w_{u,v}^i$ tells how many walks of length exactly $i$ are from vertex $u$ to $v$.

With such an algorithm we can easily implement @distance-queries distance queries by linearly scanning (see Section \[fft-algorithm\]). Thus, for the multiple invariant factors we will skip the description of algorithm that returns the number of walks of length smaller than $k$ (the technique is the same).

Multiple Invariant Factors
--------------------------

Now, we will consider a case when $k \ge 1$, i.e., matrix $F$ has multiple invariant factors. First of all, we need to note that this generalization is not perfect and will allow only the walks of length up to $D$ (the longest distance in a graph, i.e., diameter).

In the real world applications of our framework (detecting cycles, determining distance between nodes, etc.) we do not need walks longer than the longest possible distance in a graph. It is natural that the diameter is considered to be a bound of an output in graph problems [@diameter; @cygan; @vanessa2015; @vanessa2016].

### Relation of the Graph Diameter and Frobenius Normal Form

We begin with relating the graph diameter to the Frobenius normal form. It turns out that the graph diameter is bounded by the degree of a smallest invariant factor.

\[mu\] Given a directed, unweighted graph $G$ with a diameter $D$. Let $\mu$ denote the degree of the smallest invariant factor (i.e., the dimension of the smallest block in the Frobenius matrix $F$) of an adjacency matrix of the graph $G$. Then $D \le \mu$.

The bounds of this inequality are tight. There are graphs with diameter $D =
\mu$ and graphs with $\mu = n$ and arbitrary small diameter [@spectral-graph]. Our algorithms are able to return walks up to the length $\mu$. We use the bound on $D$ solely because it is easier to interpret diameter than the *smallest degree of the invariant factor*.

### Generalization to Multiple Invariant Factors

(5,5) rectangle (50,50); (-5,-5) rectangle (-50,-50) ; (55,5) rectangle (100,50) node\[pos=.5\] [$U^{-1}$]{};

(5,-5) rectangle (50,-50);

(5,50) rectangle (20,35) node\[pos=.5\] [$F_1$]{}; (20,35) rectangle (40,15) node\[pos=.5\] [$F_2$]{}; (40,15) rectangle (50,5) node\[pos=.5\] [$F_3$]{};

(5,-5) rectangle (20,-50) node\[pos=.5\] [$U_1 F_1$]{}; (20,-5) rectangle (40,-50) node\[pos=.5\] [$U_2 F_2$]{}; (40,-5) rectangle (50,-50) node\[pos=.5\] [$U_3 F_3$]{};

(20,5) – (20,50); (40,5) – (40,50);

(-50,-5) rectangle (-35,-50) node\[pos=.5\] [$U_1$]{}; (-35,-5) rectangle (-15,-50) node\[pos=.5\] [$U_2$]{} ; (-15,-5) rectangle (-5,-50) node\[pos=.5\] [$U_3$]{} ;

Let $k$ denote the number of blocks in the Frobenius matrix $F$ and $\mu$ be the number of columns of the smallest block. To multiply matrix $U$ by $F$ we can start by multiplying strips of matrix $U$ by appropriate blocks of $F$ and concatenate them later (see Figure \[fig:block-multi\]).

We start by splitting the matrix $U$ into $k$ strips with rows corresponding to the appropriate blocks of $F$ (strip $U_i$ has as many columns as block $F_i$). Then we multiply $U F$ and have $k$ strips: $U_1 F_1, U_2 F_2, \ldots U_k
F_k$ (each with at least $\mu$ columns). Next, we multiply $U F^\mu$ and we keep $k$ strips: $U_1 F_1^\mu, U_2 F_2^\mu,\ldots,U_k F_k^\mu$. Our goal is to get a data structure such that if we need $U F^k$, we can quickly choose appropriate columns and append them.

The matrix $U_i F_i$ has $l_i$ columns: $v_1,\ldots, v_{l_i}$. Because $F_i$ is a companion matrix, the $U_i F_i^\mu$ has the cyclic property (Definition \[cyclic-property\]). And the matrix $U_i F_i^\mu$ has columns: $v_\mu,\ldots, v_{\mu+l_i}$. There are some duplicate columns in $U_i
F_i$ and $U_i F_i^\mu$, when $\mu < l_i$. Hence, we only need to keep columns $v_1,\ldots,v_{\mu+l_i}$ for this strip. We do this for all strips (see Figure \[combining\]).

(0,0) rectangle (15,10) node\[pos=.5\] [$U_1 F_1$]{}; (10,0) rectangle (25,10) node\[pos=.5, right\] [$U_1 F_1^\mu$]{};

(30,0) rectangle (55,10) node\[pos=.5\] [$U_2 F_2$]{}; (40,0) rectangle (65,10) node\[pos=.5, right, xshift=10pt\] [$U_2 F_2^\mu$]{};

(70,0) rectangle (80,10) node\[pos=.5\] [$U_3 F_3$]{}; (80,0) rectangle (90,10) node\[pos=.5\] [$U_3 F_3^\mu$]{};

We are left with a matrix that has at most $2n$ columns (because $l_1 + \mu + l_2 + \mu + \ldots l_k + \mu = k \mu + \sum^n_{i=1} l_i = n + k\mu
\le 2n$). To generate it we need to power $F$ to $\mu$ and do multiplications $U
\cdot F$ and $U \cdot F^\mu$. This can be computed in ${\widetilde{O}}(n^\omega)$ time.

### Queries with Multiple Invariant Factors {#query-algorithm}

When a query for the number of walks of length $k$ from node $u$ to $v$ comes, we do:

1.  For each strip $i$ take the $u$-th row of $U_i F_i \oplus U_i
            F_i^\mu$ concatenate them (see Figure \[quering-multiple\]) into vector $\bar{u}$,

2.  Take $v$-th column of $U^{-1}$ matrix and denote it $\bar{v}$,

3.  Return the dot product $\bar{u} \cdot \bar{v}$.

Because $l_1 + l_2 + \ldots + l_k = n$ the vector $\bar{u} \in K^n$. Query needs $O(n)$ time.

(0,0) rectangle (25,10) node\[above, xshift=-40pt\] [$U_1 F_1 \oplus U_1 F_1^\mu$]{};

(5,0) – (5,10) node\[pos=-0.25\] [i]{}; (20,0) – (20,10) node\[pos=-0.25\] [$i+l_1$]{}; (5,7) rectangle (20,5);

(30,0) rectangle (65,10) node\[above,xshift=-70pt\] [$U_2 F_2 \oplus U_2 F_2^\mu$]{};

(35,0) – (35,10) node\[pos=-0.25\] [i]{}; (60,0) – (60,10) node\[pos=-0.25\] [$i+l_2$]{}; (35,7) rectangle (60,5);

(70,0) rectangle (90,10) node\[above, xshift=-40pt\] [$U_3 F_3 \oplus U_3 F_3^\mu$]{}; (75,0) – (75,10) node\[pos=-0.25\] [i]{}; (85,0) – (85,10) node\[pos=-0.25\] [$i+l_3$]{}; (75,7) rectangle (85,5);

Finally, this dot product is $a_{u,v}$ element of the matrix $U F^k U^{-1}$, for a fixed $k \le \mu$. Analogously to Section \[single-factor\] one can extend this result to return the number of walks of length less or equal $k$. This matches (up to the polylogarithmic factor) the result of @distance-queries. However in the next section we will show a more general result.

Almost Optimal Query {#fft-algorithm}
====================

In the previous section, we showed how to preprocess a directed, unweighted graph in ${\widetilde{O}}(n^\omega)$ time in such a way that in $O(n)$ query we can return a number of distinct walks of a length $k$ from vertex $u$ to $v$. However, in linear time $O(n)$ we return only a single number. Our goal is to get far richer information in $\tilde{O}(n)$ query time.

\[fft-query\] Let $G = (V,E)$ be a directed, unweighted graph with $n$ vertices and a diameter $D$. There exists an algorithm that after some preprocessing can answer queries for any given $u,v \in V$:

-   Query returns $\{ w_i \; | \; 1\le i \le D \}$, where $w_i$ is the number of distinct walks from $u$ to $v$ of length exactly $i$,

-   Preprocessing takes ${\widetilde{O}}(n^\omega)$ and query takes $\tilde{O}(n)$ field operations.

The algorithm is deterministic (but there are $O(\log W)$ factors in derandomized version, see Section \[prelim\] for explication).

This algorithm is a significant improvement over @distance-queries:

-   One can use Theorem \[fft-query\] to find the distance between $u,v$ by linearly scanning the array and returning the first $k$ such that $w_k > 0$,

-   Theorem \[fft-query\] can count cycles. In contrast the @distance-queries cannot, because the distance from $u$ to $u$ is always $0$ (we will see that in the next section),

-   Theorem \[fft-query\] is almost optimal (when $D=O(n)$ its output and time are $O(n)$).

On the other hand, Theorem \[fft-query\] is only a functional improvement and it does not break the ${\widetilde{O}}(n^\omega)$ of the *Single Source Shortest Paths* (SSSP) for dense, directed graphs.

Hankel Matrix
-------------

Now, we will focus on the proof of Theorem \[fft-query\]. First, we need to introduce the Hankel matrix and its properties.

$$H = \left(\begin{array}{cccccc}
            c_1 & c_2 & \ldots & c_n \\
            c_2 & c_3 & \ldots & c_{n+1} \\
            \vdots & \vdots & & \vdots \\
            c_n & c_{n+1} & \ldots & c_{2n-1} 
    \end{array}\right)$$

Hankel matrix is defined by its first row and last column ($2n-1$ numbers define $n\times n$ Hankel matrix). The numbers from the previous row are left-shifted by one and the new number is added at the end. Hankel matrices have some similarities with Topelitz and Circulant matrices.

The basic property we need is that the product of Hankel matrix and vector can be computed in $O(n\log n)$ time (see [@hankel; @hankel2]) even though explicitly writing the Hankel matrix as $n\times n$ matrix takes $O(n^2)$ time. The algorithm takes $2n-1$ parameters that define the Hankel matrix and $n$ parameters that define the vector. The technique is based on the Fast Fourier Transformation [@hankel; @hankel2].

Using Hankel Matrices to Return Number of Walks
-----------------------------------------------

The algorithm from Section \[query-algorithm\] concatenates the strips $U_i
F_i$ and builds a single vector. Subsequently, that vector is multiplied by a column of matrix $U^{-1}$. But we can also do it in a different order: first we multiply the strip by a section of matrix $U^{-1}$ and sum the results at the end. Thus, we perform $k$ (number of strips) dot products of smaller vectors (see Figure \[fft-strips\]).

(0,0) rectangle (2.5,1.0) node\[below=6pt, pos=.5\] [$U_1 F_1$]{};

(0,0.7) rectangle (2.5,0.5); (0.5,0.7) rectangle (2.0,0.5) node\[above=2pt,pos=.5\] (start\_U1) ;

(3.0,0) rectangle (6.5,1.0) node\[below=6pt, pos=.5\] [$U_2 F_2$]{};

(3.0,0.7) rectangle (6.5,0.5); (3.5,0.7) rectangle (6.0,0.5) node\[above=2pt,pos=.5\] (start\_U2) ;

(7.0,0.0) rectangle (9.0,1.0) node\[below=6pt, pos=.5\] [$U_3 F_3$]{}; (7.0,0.7) rectangle (9.0,0.5); (7.5,0.7) rectangle (8.5,0.5) node\[above=2pt,pos=.5\] (start\_U3) ;

(10.0,1.5) rectangle (11.0,2.5); (10.3,1.5) rectangle (10.5,2.5) node\[left=2pt, pos=.5\] (end\_U3) ;

(10.0,3.0) rectangle (11.0,5.5); (10.3,3.0) rectangle (10.5,5.5) node\[left=2pt, pos=.5\] (end\_U2) ;

(10.0,6.0) rectangle (11.0,7.5); (10.3,6.0) rectangle (10.5,7.5) node\[left=2pt, pos=.5\] (end\_U1) ;

(11.2,1.5) – node\[right=6pt\] [$U^{-1}$]{} (11.2,7.5); (start\_U1) to\[out=90,in=180\] (end\_U1); (start\_U2) to\[out=90,in=180\] (end\_U2); (start\_U3) to\[out=90,in=180\] (end\_U3);

Consider the query for a number of walks of length exactly $k$. The strips in the matrix $U^{-1}$ do not depend on $k$ (vector $(u_0,\ldots,u_l)$). However, the vector taken from $U_i F_i$ (vectors $(x_i,\ldots,x_{i+l})$) will be left shifted if we want to compute the next one.

$$\begin{array}{ccccccccr}
            ( x_0 & x_1 & x_2 & \ldots & x_l )\\
            ( x_1 & x_2 & \ldots & x_l & x_{l+1} )\\
            ( x_2 & \ldots & x_l & x_{l+1} & x_{l+2})\\
            \vdots & & & & \vdots \\
            (x_\mu & &\ldots & & x_{\mu+l}) \\
        \end{array}
        \times
        {
    \global\colveccount5
    \begin{pmatrix}
        \colvecnext
    }{u_0}{}{\vdots}{}{u_l}$$

As you can see, the subsequent rows can be written as the Hankel matrix (we need to add some zeros and discard certain results to get square matrix, but it will not influence asymptotic complexity). By using the *fast Hankel matrix-vector multiplication* we can compute $\mu$ values for every strip $i$ in time $O(l_i \log l_i)$ ($l_i$ was defined as the length of $i$-th strip). At the end, we need to sum all results into a single array. The total complexity is $O(\sum_{i=1}^k l_i \log l_i))$. Because $\sum^k_{i=1} l_i = n$ the algorithm needs $O(n\log n)$ field operations. It proves Theorem \[fft-query\].

We need to address some issues regarding field operations. As mentioned in the Section \[prelim\], the $\tilde{O}$ notation hides the polylogarithmic factors. Here, we silently assume that the number of walks is bounded by $W$. It means that the complexity is multiplied by $O(\log W)$ factor because of the cost of arithmetic operations. If the number of walks is exponential in $n$ then the cost increases by a linear factor. However, at the beginning we could randomly select the prime number $p$ with $O(\log n)$ bits and do arithmetic operations in the field $\mathbb{Z}_p$. In some applications we can use our algorithm to answer whether with high probability there exists a walk of a given length. The problem of large number of distinct walks is more of a theoretical issue than a practical one.

All Pairs All Walks (APAW)
--------------------------

Now we will show the application of Theorem \[fft-query\]. We begin with almost optimal algorithm to compute the number of all walks between all pairs of vertices. We are not aware of any other research concerning this problem.

Given a directed, unweighted graph $G$ with a diameter $D$. The task is to return an array $A$, such that for every pair of vertices $u,v \in G$ and every $k \in \{1,\ldots,D\}$ an element $A[u,v,k]$ is the number of distinct walks from $u$ to $v$ of length $k$.

The only solution to this problem we are aware of needs $O(D
n^{\omega})$ time. The naive approach: *take the adjacency matrix $A$ of graph $G$ and save it in $A[u,v,1]$. Then, square it to get $A^2$ and save it in $A[u,v,2]$. Continue until you fill out complete table*. In the worst case this algorithm needs $D = O(n)$ matrix multiplications, thus its complexity is $O(D n^{\omega})$. At the first glance it is surprising that we can improve it to ${\widetilde{O}}(n^3)$ especially when $\omega > 2$.

The $\tilde{O}(n^3)$ algorithm works as follows. First, preprocess the algorithm from Theorem \[fft-query\] which takes ${\widetilde{O}}(n^\omega)$ time. Then, for every pair of vertices $u,v$ ask a query. A single query takes $\tilde{O}(n)$ time. Then, save it in the table $A[u,v]$ (query gives $D$ numbers $w_1,\ldots,w_D$, such that $w_i$ is the number of walks of length i and save it $A[u,v,i] \coloneqq w_i$).

Because there are $O(n^2)$ pairs and for each pair we need $\tilde{O}(n)$ time, the complexity of our solution is $\tilde{O}(n^3)$. The algorithm is almost optimal because the output in the worst case may be $O(n^3)$.

Counting and Determining the Lengths of Cycles
==============================================

We will use Theorem \[fft-query\] to solve All-Nodes Shortest Cycles (ANSC) problem efficiently.

\[t-cycles\] There exists an algorithm that for a given unweighted digraph $G$ with a diameter $D$:

-   For every vertex $u$ returns $D$ numbers: $c_u^1, c_u^2, \ldots c_u^D$,

-   The $c_u^k$ is the number of non-simple cycles of length exactly $k$, that contain vertex $u$,

-   Algorithm works in ${\widetilde{O}}(n^\omega \log{W})$ time (where $W$ is the maximum $c^k_u$).

We will use Theorem \[fft-query\]. We start by preprocessing the graph $G$ in time ${\widetilde{O}}(n^\omega)$. Theorem \[fft-query\] allows us to ask for a number of walks from $u$ to $v$ and receive $D$ numbers: $w_{u,v}^k$. So, we ask for the number of walks from vertex $u$ to the same vertex $u$. This is exactly the number of non-simple cycles of a given length that contain vertex $u$.

Because we need to ask only $n$ queries (it is the number of vertices in a graph) and each query takes $\tilde{O}(n)$ time we have ${\widetilde{O}}(n^\omega
    + n^2)$ = ${\widetilde{O}}(n^\omega)$ algorithm.

Solving ANSC Faster
-------------------

To solve ANSC problem in ${\widetilde{O}}(n^\omega)$ time and beat @cycle-each-vertex ${\widetilde{O}}(n^{(\omega + 3)/2})$ algorithm we do the linear search on the output. For every vertex we search for the first nonzero element linearly. This with high probability is the length of the shortest cycle that contains it. Because the output contains $O(n^2)$ numbers the complexity is equal to the preprocessing time ${\widetilde{O}}(n^\omega)$.

Similarly, we can scan the output to compute the set $S(c)$ that contains all vertices that lie on some cycle of length $\le c$. Then, by linear scan we can return the sets $S(1),\ldots,S(D)$. This improves upon @cygan.

Conclusion and Open Problems
============================

We introduced the framework based on Frobenius normal form and used it to solve problems on directed, unweighted graphs in matrix multiplication time. The main open question is to use this framework to prove that APSP on such graphs can be solved in ${\widetilde{O}}(n^\omega)$ or at least $O(n^{2.5})$. The promising way is to use the algorithms that determine operators of matrices of polynomials (e.g., determinant, solving linear system [@storjohann2; @storjohann3]). Additionally, algorithms for a black-box polynomial degree determination seem to be a promising path.

Another interesting problem is to use this framework to obtain additive approximation for APSP. Currently, the best additive approximation of APSP is due to [@rodity]. However, none of known additive approximation of APSP works in ${\widetilde{O}}(n^\omega)$ time.

Application in dynamic algorithm also seems to be a promising approach. @dynamic-frobenius showed an algorithm, that dynamically preserves Frobenius normal form in $O(k n^2)$ time. Our algorithms use fast Hankel matrix-vector multiplication. The algorithm behind fast Hankel matrix-vector multiplication is based on Discrete Fourier Transform (DFT). @rief presented an $O(\sqrt{n})$ time per request algorithm for DFT. Can we use those approaches to obtain a faster dynamic algorithm?

Finally, it remains open how to apply the Frobenius normal form in the weighted directed graphs with small, integer weights $\{-M,\ldots,M\}$. @cygan took degree $M$ polynomials and used [@storjohann2] to get radius and diameter in ${\widetilde{O}}(M n^\omega)$ time. We suspect that similar technique can be applied to Frobenius normal form framework.

Acknowledgment
==============

The work of Karol Wgrzycki is part of a project TOTAL that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 677651). Piotr Sankowski is supported by the Polish National Science Center, grant no UMO-2014/13/B/ST6/01811. We would also like to thank Agnieszka Ciepielewska for meaningful discussions and useful comments.

Introduction to the Cyclic Subspaces {#appendix-cyclic}
====================================

We have a constant-recursive sequence

$$a_n = c_0 a_{n-1} + c_1 a_{n-2} + \ldots + c_{r-1} a_{n-r}$$

of order $r$, where all $c_i$ are constants and initial conditions $a_0, \ldots,
a_{r-1}$ are given. Perhaps, the most familiar example of such sequence are the Fibonacci numbers $F_n = F_{n-1} + F_{n-2}$.

When we take the companion matrix:

$$C = 
    \begin{bmatrix}
        0&\ldots&0      & -c_0\\
        1&\ddots&       & -c_1\\
         &\ddots&      0&\vdots \\
        0&      &      1&-c_{r-1}\\
    \end{bmatrix}
    .$$

We can generate the next element of sequence by multiplying it:

$$C^T
    \begin{bmatrix}
        a_{n-r} \\
        \vdots  \\
        a_{n-2} \\
        a_{n-1} \\
    \end{bmatrix}
    =
    \begin{bmatrix}
        a_{n-r+1} \\
        \vdots    \\
        a_{n-1}   \\
        a_{n}     \\
      \end{bmatrix}
      .$$

For example, for Fibonacci sequence we have:

$$\begin{bmatrix}
        0 & 1 \\
        1 & 1 \\
    \end{bmatrix}
    \begin{bmatrix}
        F_k \\
        F_{k+1} \\
    \end{bmatrix}
    =
    \begin{bmatrix}
        F_{k+1} \\
        F_{k+2} \\
    \end{bmatrix}
    .$$

Now, if we want the $a_{n+1}$ element of the series we can square the companion matrix $C^2$ and get:

$$C^T C^T
    \begin{bmatrix}
        a_{n-r} \\
        \vdots  \\
        a_{n-2} \\
        a_{n-1} \\
    \end{bmatrix}
    =
    C^T
    \begin{bmatrix}
        a_{n-r+1} \\
        \vdots  \\
        a_{n-1} \\
        a_{n} \\
    \end{bmatrix}
    =
    \begin{bmatrix}
        a_{n-r+2} \\
        \vdots    \\
        a_{n}   \\
        a_{n+1}     \\
      \end{bmatrix}
      .$$

In linear algebra such transformations are known as the *cyclic subspaces* generated by the vector $a$. For our purposes, we will restrict ourself to a conclusion, that some columns of a companion matrix occur in its powers. These properties are well known in the linear algebra theory (see [@gantmacher; @eberly] for more cyclic properties).

Example of Powering a Companion Matrix {#appendix-example}
======================================

Example:

$$\begin{aligned}
C =
\begin{bmatrix}
   0 & 0 & 0 & 0 & 1 \\
   1 & 0 & 0 & 0 & 2 \\
   0 & 1 & 0 & 0 & 3 \\
   0 & 0 & 1 & 0 & 4 \\
   0 & 0 & 0 & 1 & 5 \\
\end{bmatrix}, \;
C^2 =
\begin{bmatrix}
    0 & 0 & 0 & 1 & 5 \\
    0 & 0 & 0 & 2 & 11 \\
    1 & 0 & 0 & 3 & 17 \\
    0 & 1 & 0 & 4 & 23 \\
    0 & 0 & 1 & 5 & 29 \\
\end{bmatrix}, \;
C^3 = 
\begin{bmatrix}
    0 & 0 & 1 & 5  & 29 \\
    0 & 0 & 2 & 11 & 63 \\
    0 & 0 & 3 & 17 & 98 \\
    1 & 0 & 4 & 23 & 133 \\
    0 & 1 & 5 & 29 & 168 \\
\end{bmatrix}, \end{aligned}$$

$$\begin{aligned}
C^4 = 
\begin{bmatrix}
    0 & 1 & 5 & 29 & 168 \\
    0 & 2 & 11& 63& 365 \\
    0 & 3 & 17& 98& 567 \\
    0 & 4 & 23& 13& 770 \\
    1 & 5 & 29& 16& 973 \\
\end{bmatrix}, \;
C^5 = 
\begin{bmatrix}
    1 & 5  & 29  & 168 & 973  \\
    2 & 11 & 63  & 365 & 2114 \\
    3 & 17 & 98  & 567 & 3284 \\
    4 & 23 & 133 & 770 & 4459 \\
    5 & 29 & 168 & 973 & 5635 \\
\end{bmatrix}
C^6 = 
\begin{bmatrix}
    5  & 29  & 168 & 973  & 5636  \\
    11 & 63  & 365 & 2114 & 12243 \\
    17 & 98  & 567 & 3284 & 19019 \\
    23 & 133 & 770 & 4459 & 25824 \\
    29 & 168 & 973 & 5635 & 32634 \\
\end{bmatrix}\end{aligned}$$

The matrix $C^i$ has $4$ columns identical to matrix $C^{i+1}$. $C$ has coefficients of order dimension (dimension is 5 and maximum coefficient is 5). After powering to the 5th power, the coefficients can be of order $5^5$. Over a finite field $\mathbb{Z}_p$, all those coefficients will have $O(\log p)$ bits.

Proof of Theorem \[mu\] {#proof-mu}
=======================

Lets assume, that $D > \mu$ and w.l.o.g. $\delta(u,v) = D$. Let $A$ be the adjacency matrix of $G$. We know, that there is the minimal polynomial of degree $\mu$:

$$A^\mu = a_0 \mathbb{1} + a_1 A + a_2 A^2 + \ldots + a_{\mu-1} A^{\mu-1}
        .$$

Term $a^k_{i,j}$ denotes the $i,j$ element of the matrix $A^k$. Now, consider the elements $u,v$ of each matrix. The diameter $D > \mu$ and $\delta(u,v) = D$, so for every $k \le \mu$ the elements $a^k_{u,v} = 0$ (because there is no walk of length less than $D$ from $u$ to $v$). Now, if we multiply the minimal polynomial by the matrix $A$ we get:

$$A^{\mu+1} = a_0 A + a_1 A^2 + a_2 A^3 + \ldots + a_{\mu-1} A^{\mu}
        .$$

Hence $a^{\mu+1}_{u,v} = 0$, because every element in the sequence $a^k_{u,v} = 0$ for $k \le \mu$. By repeating this reasoning, we get that for every $k > 0$ the element $a^k_{u,v} = 0$. So, for every achievable pair of vertices, there must be some $k \le \mu$, such that $a^k_{u,v} \ne 0$ and diameter is bounded by $\mu$.

[^1]: Institute of Informatics, University of Warsaw, Poland,

[^2]: Institute of Informatics, University of Warsaw, Poland,

[^3]: $\mathrm{Tr}\left[A\right]$ denotes the trace of a matrix $A$, i.e., the sum of elements on main diagonal.

[^4]: sometimes called the rational-canonical form

[^5]: One can compute $F^{n+1}$ even faster. Namely, in $\tilde{O}(n^2)$ time by applying the cyclic property (see Definition \[cyclic-property\])

[^6]: operation $\oplus$ denotes concatenation

[^7]: $X_{a,b}$ denotes a matrix constructed by concatenating columns $x_a,x_{a+1},\ldots,x_{b}$ of a matrix $X$
