2012-09-07 29 views
16

単純な除算および征服アルゴリズムよりもM^n(Mは行列、nは整数)を計算するための行列累乗の方法はありません。高速行列累乗

+1

ねえ、私はstackoverflowの内の1つのリンクを見つけただけで、それをチェックアウト http://stackoverflow.com/questions/12268516/matrix-exponentiation -using-fermats-theorem –

+0

Expokitは、行列累乗を実行するためのよく知られたパッケージです。 http://fortranwiki.org/fortran/show/Expokit – Sayan

答えて

22

行列を固有値と固有ベクトルに分解できます。そして、あなたは得る。

M = V^-1 * D * V 

ここで、Vは固有ベクトル行列であり、Dは対角行列である。すべてのVおよびV^-1項はキャンセルので

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

:N乗にこれを上げるために、あなたのような何かを得ます。

Dは対角線であるため、完全な行列ではなく、n乗の(実際の)数を増やすだけです。あなたはnの対数時間でそれを行うことができます。

固有値と固有ベクトルの計算はr^3(rはMの行数/列数)です。 rとnの相対的な大きさに応じて、これは速くてもいいかもしれません。

+0

私の知る限り、この方法はSquaringによる累乗と同じ複雑さを持っています。だから、もっと速い方法がありますか? –

+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; 
}