2012-09-07 170 views
16

是否有任何更快的矩陣求冪方法來計算M^n(其中M是一個矩陣,n是一個整數)比簡單的除法和征服算法。快速矩陣求冪

+1

嘿,我發現了一個計算器只鏈接檢查出來 http://stackoverflow.com/questions/12268516/matrix-exponentiation -using-fermats-theorem – 2012-09-07 04:57:52

+0

Expokit是一個衆所周知的用於執行矩陣求冪的包。 http://fortranwiki.org/fortran/show/Expokit – Sayan

回答

22

您可以將矩陣分解爲特徵值和特徵向量。那麼你得到

M = V^-1 * D * V 

其中V是特徵向量矩陣和D是一個對角矩陣。爲了提高到N次方,您得到如下結果:

M^n = (V^-1 * D * V) * (V^-1 * D * V) * ... * (V^-1 * D * V) 
    = V^-1 * D^n * V 

因爲所有的V和V^-1項都被取消了。由於D是對角線,你只需要將一堆(實數)數字提升到第n次方,而不是滿矩陣。你可以在n的對數時間內做到這一點。

計算特徵值和特徵向量是r^3(其中r是M的行數/列數)。取決於r和n的相對大小,這可能會更快或者不更快。

+0

據我所知,這種方法與Squaring的Exponentiation具有相同的複雜度。那麼有沒有更快的方法? –

+0

@AkashdeepSaluja:這比通過平方冪取冪更快。這是O(r^3)時間,平方的指數是O(r^3 logn)時間。 –

+0

爲更好地解釋上述方法http://www.google.co.in/url?sa=t&rct=j&q=pdf%20nth%20power%20of%20matrix&source=web&cd=1&cad=rja&ved=0CCAQFjAA&url=http% 3A%2F%2Fwww.qc.edu.hk%2Fmath%2FTeaching_Learning%2FNth%2520power%2520of%2520a%2520square%2520matrix.pdf&ei = Jf9JULrwFsi8rAejh4C4DQ&usg = AFQjCNE7yqQce5jdtyyVLFpSZmYUnoWyVA –

4

Exponentiation by squaring經常用來獲得矩陣的高次冪。

+0

我知道這種方法,但需要進一步加速。 –

+0

您最好在問題中添加此算法名稱以避免類似的答案:) – MBo

+0

更快的算法要複雜得多。 – Ari

0

我會推薦用於計算matrix form中Fibbonacci序列的方法。 AFAIK,其效率是O(log(n))。

+1

您必須乘以矩陣相乘的成本。整體運行時間爲O(n^3 log n)。 – saadtaame

6

使用歐拉快速冪算法很簡單。使用下一個算法。

#define SIZE 10 

//It's simple E matrix 
// 1 0 ... 0 
// 0 1 ... 0 
// .... 
// 0 0 ... 1 
void one(long a[SIZE][SIZE]) 
{ 
    for (int i = 0; i < SIZE; i++) 
     for (int j = 0; j < SIZE; j++) 
      a[i][j] = (i == j); 
} 

//Multiply matrix a to matrix b and print result into a 
void mul(long a[SIZE][SIZE], long b[SIZE][SIZE]) 
{ 
    long res[SIZE][SIZE] = {{0}}; 

    for (int i = 0; i < SIZE; i++) 
     for (int j = 0; j < SIZE; j++) 
      for (int k = 0; k < SIZE; k++) 
      { 
       res[i][j] += a[i][k] * b[k][j]; 
      } 

    for (int i = 0; i < SIZE; i++) 
     for (int j = 0; j < SIZE; j++) 
      a[i][j] = res[i][j]; 
} 

//Caluclate a^n and print result into matrix res 
void pow(long a[SIZE][SIZE], long n, long res[SIZE][SIZE]) 
{ 
    one(res); 

    while (n > 0) { 
     if (n % 2 == 0) 
     { 
      mul(a, a); 
      n /= 2; 
     } 
     else { 
      mul(res, a); 
      n--; 
     } 
    } 
} 

下面請找出等值數字:

long power(long num, long pow) 
{ 
    if (pow == 0) return 1; 
    if (pow % 2 == 0) 
     return power(num*num, pow/2); 
    else 
     return power(num, pow - 1) * num; 
}