如何求有向圖的鄰接矩陣的冪斂指數和周期？

01-03

一個有向圖的鄰接矩陣A是一個n階布爾矩陣，它的冪組成的序列會具有周期性。
（注意，計算冪的時候要用布爾代數，即乘法是與，加法是或）
設k,d是滿足 $A^k = A^{k+d}$ 的最小正整數，則k稱為A的冪斂指數（index），d稱為A的周期（period）。
舉個例子：

上面這個圖的鄰接矩陣為
$A = left[ egin{array}{ccccc} 0 1 0 0 0 \ 0 0 1 0 0 \ 0 0 0 1 0 \ 0 0 0 0 1 \ 0 0 1 0 0 end{array} ight]$
它的各次冪為：
$A^2 = left[ egin{array}{ccccc} 0 0 1 0 0 \ 0 0 0 1 0 \ 0 0 0 0 1 \ 0 0 1 0 0 \ 0 0 0 1 0 end{array} ight]$
$A^3 = left[ egin{array}{ccccc} 0 0 0 1 0 \ 0 0 0 0 1 \ 0 0 1 0 0 \ 0 0 0 1 0 \ 0 0 0 0 1 end{array} ight]$
$A^4 = left[ egin{array}{ccccc} 0 0 0 0 1 \ 0 0 1 0 0 \ 0 0 0 1 0 \ 0 0 0 0 1 \ 0 0 1 0 0 end{array} ight]$
$A^5 = left[ egin{array}{ccccc} 0 0 1 0 0 \ 0 0 0 1 0 \ 0 0 0 0 1 \ 0 0 1 0 0 \ 0 0 0 1 0 end{array} ight]$
此時發現 $A^5 = A^2$ ，故k=2, d=3。

給定A，有什麼演算法可以求k和d？
我搜到了這麼幾篇文獻，但都並不是直接解決上述問題：
[1] COMPUTING INDEXES AND PERIODS OF ALL
BOOLEAN MATRICES UP TO DIMENSION N = 8
http://www.cai.sk/ojs/index.php/cai/article/viewFile/1310/485
[2] On the semigroup of binary relations on a finite set
http://dml.cz/bitstream/handle/10338.dmlcz/100989/CzechMathJ_20-1970-4_9.pdf
[3] 布爾矩陣的冪斂指數集
http://advmath.pku.edu.cn/CN/article/downloadArticleFile.do?attachType=PDFid=9421

經過跟 @七月的討論，現在得出了這麼一種求周期d的方法：

用Tarjan演算法求出圖中所有強連通分量；
求每個強連通分量的周期，這等於強連通分量中所有環的長度的最大公約數；
整個圖的周期等於每個強連通分量的周期的最小公倍數。

這裡2、3兩步的原理出現在下面文獻[4]中定理2.5下面。

上面的第2步，評論中的 @吼吼吼給出了如下的演算法：

任取強連通分量中的一個點x作為起點進行DFS，並記錄每個點的層號；
在DFS過程中，若發現第i層的某個點到第j層的某個點有邊，則說明圖中有一個長度為 $d=|i+1-j|$ 的環，或者有兩個長度相差d的環；
對上一步中得到的所有d求最大公約數，就是這個強連通分量的周期。

對於一個有n個頂點的稠密圖，Tarjan演算法和求每個強連通分量的周期的複雜度為O(n^2)；對V個點、E條邊的稀疏圖則是O(V + E)。（求最大公約數和最小公倍數的時間看作常數）

關於冪斂指數k， @劉城提出，可以用二分法來求：

首先，用矩陣快速冪演算法算出 $A^d$ ；
然後，從 $p=1$ 開始，每次把p增大一倍，直到 $A^p = A^{p+d}$ ，這樣就可以知道 $p/2 < k le p$ ；
在這個區間內，用二分法尋找使 $A^k = A^{k+d}$ 成立的最小k值。

每次矩陣乘法的複雜度是O(n^3)。

第1步需要O(log d)次矩陣乘法，第2步需要O(log k)次矩陣乘法。

如果在第2步中存儲了 $A^1, A^2, A^4, A^8, ...$ 的值，則第3步每次區間折半隻需要1次矩陣乘法，驗證 $A^k = A^{k+d}$ 也只需要1次矩陣乘法，總共需要O(log k)次矩陣乘法。

題目描述中的文獻[1]以及下面的文獻[4]指出，k的上界是(n-1)^2 + 1，即O(n^2)。

而d則可能比較大，比如一個圖中若含有4個不相交的環，環長分別為2,3,5,7，則d = 2*3*5*7 = 210。我估算O(log d) = $O(sqrt{n} log n)$ （可能不夠緊）。

這樣，求k的總複雜度就是 $O(n^{3.5} log n)$ 。

[4] On the sequence of consecutive powers of a
matrix in a Boolean algebra

http://www.dcsc.tudelft.nl/~bdeschutter/pub/rep/97_67.pdf

糾正一下吧……至少有這麼個結論：k為不經過重複點的最長路長度，d為圖中所有環的長度的最小公倍數……

好吧有反例……

————————————————————

其實這跟鏈表找環沒啥區別……大概最簡單的方法是所謂的步差法（追趕法、快慢指針法）。空間O(N^2)，時間O(N^2M)，其中M=k+d。多個$N^2是因為要驗證矩陣

說白了就是逐個檢查是否有滿足 $A^{k+1}=A^{2k+1}$ 的k。具體做的時候我們可以考慮這樣：

首先來個『慢矩陣』S0=A，以及一個『快矩陣』F0=A，然後循環做Fi=Fi-1 * A^2，Si=Si-1 * A。如果有環，必定會遇到Fi=Si的情況（相遇）。可以證明在Si走完一圈時就必然相遇。記錄相遇時的下標為p。

接下來求d：繼續迭代往下算Si，如果出現Sq=Fp，說明Si從p點開始又走了一圈回到了p，此時可以得到我們要的d=q-p。

那求k怎麼求呢？

來個S"0=A^d, F"0=A，然後再找重複最初的迭代。初次相遇的位置（滿足S"k=F"k）正好是圓環行程差，此時的k就是環開始的位置了。

手機就偽代碼了變數方便起見都int

struct info {k,d}

array of info infos[n]

for i in 1…n

array mink[n] = {0}

pos = i

j =0

loop

pos = nextpos(pos)

j++

until mink[pos]!=0 or pos==i

otherwise mink[pos]= j continue loop

infos[i] = {mink[pos], j - mink[pos]}

continue forloop

k = max k in infos

d = least common multiple of d in infos

兩層loop,如果靜態化一個nextpos的查表就是O(n^2),否則nextpos也是一層循環則共計O(n^3)

=======

實際一個因數n(inner loop)應為max d。不過max d ≦ n

給一個 naive 地利用 Floyd判圈演算法（Floyd"s cycle detection algorithm）的 C++11 實現。

#include & #include &


template &

struct BoolMatrix {

    static BoolMatrix identity() {

        BoolMatrix r;

        for (size_t i = 0; i &< n; i++)
            for (size_t j = 0; j &< n; j++)
                r.m[i][j] = (i == j);
        return r;
    }

    BoolMatrix operator*(const BoolMatrix rhs) const {
        BoolMatrix r;
        for (size_t i = 0; i &< n; i++)
            for (size_t j = 0; j &< n; j++) {
                bool e = false;
                for (size_t k = 0; k &< n; k++)
                    if (m[i][k]  rhs.m[k][j]) {
                        e = true;
                        break;
                    }
                r.m[i][j] = e;
            }
        return r;
    }

    bool operator==(const BoolMatrix rhs) const {
        for (size_t i = 0; i &< n; i++)
            for (size_t j = 0; j &< n; j++)
                if (m[i][j] != rhs.m[i][j])
                    return false;
        return true;
    }

    bool m[n][n];
};

using namespace std::rel_ops; // defines operator!=() based on operator==()

// https://en.wikipedia.org/wiki/Cycle_detection#Tortoise_and_hare
template &

void floyd(F f, const T x0, int mu, int lam) {

    T tortoise = f(x0);

    T hare = f(tortoise);

    while (tortoise != hare) {

        tortoise = f(tortoise);

        hare = f(f(hare));

    }
    mu = 0;

    tortoise = x0;

    while (tortoise != hare) {

        tortoise = f(tortoise);

        hare = f(hare);

        mu++;

    }
    lam = 1;

    hare = f(tortoise);

    while (tortoise != hare) {

        hare = f(hare);

        lam++;

    }

}
int main() {

    const BoolMatrix&<5&> a = {{

        { 0, 1, 0, 0, 0 },

        { 0, 0, 1, 0, 0 },

        { 0, 0, 0, 1, 0 },

        { 0, 0, 0, 0, 1 },

        { 0, 0, 1, 0, 0 }

    }};

int mu, lam; floyd( [a](decltype(a) x) { return x * a; }, decltype(a)::identity(), mu, lam); std::cout &<&< mu &<&< " " &<&< lam &<&< std::endl; }

Floyd演算法複雜度是O(k + d)，矩陣乘法複雜度 O(n^3)，因此整體時間複雜度為 O((k + d)n^3)。由於只需存儲兩個臨時矩陣，空間複雜度為 O(n^2)。

這個其實百度搜索解釋結構模型就知道了。

演算法複雜度就是 O(V + E) + O(V )

用trajan 演算法解釋結構模型方法運用--一步到位的計算方法，無須可達矩陣，速度飛快圖形表示見這個的鏈接