Strassen矩陣乘法 + 快速計算乘方的演算法 + 矩陣的次冪

1.Strassen矩陣乘法

矩陣乘法是線性代數中最常見的運算之一，它在數值計算中有廣泛的應用。若A和B是2個n×n的矩陣，則它們的乘積C=AB同樣是一個n×n的矩陣。A和B的乘積矩陣C中的元素C[i,j]定義為:　

若依此定義來計算A和B的乘積矩陣C，則每計算C的一個元素C[i,j]，需要做n個乘法和n-1次加法。因此，求出矩陣C的n2個元素所需的計算時間為0(n3)。

60年代末，Strassen採用了類似於在大整數乘法中用過的分治技術，將計算2個n階矩陣乘積所需的計算時間改進到O(nlog7)=O(n2.18)。

首先，我們還是需要假設n是2的冪。將矩陣A，B和C中每一矩陣都分塊成為4個大小相等的子矩陣，每個子矩陣都是n/2×n/2的方陣。由此可將方程C=AB重寫為:

　(1)

由此可得:　

C11=A11B11+A12B21                           (2)

C12=A11B12+A12B22                           (3)

C21=A21B11+A22B21                           (4)

C22=A21B12+A22B22                           (5)

如果n=2，則2個2階方陣的乘積可以直接用(2)-(3)式計算出來，共需8次乘法和4次加法。當子矩陣的階大於2時，為求2個子矩陣的積，可以繼續將子矩陣分塊，直到子矩陣的階降為2。這樣，就產生了一個分治降階的遞迴演算法。依此演算法，計算2個n階方陣的乘積轉化為計算8個n/2階方陣的乘積和4個n/2階方陣的加法。2個n/2×n/2矩陣的加法顯然可以在c*n2/4時間內完成，這裡c是一個常數。因此，上述分治法的計算時間耗費T(n)應該滿足：

這個遞迴方程的解仍然是T(n)=O(n3)。因此，該方法並不比用原始定義直接計算更有效。究其原因，乃是由於式(2)-(5)並沒有減少矩陣的乘法次數。而矩陣乘法耗費的時間要比矩陣加減法耗費的時間多得多。要想改進矩陣乘法的計算時間複雜性，必須減少子矩陣乘法運算的次數。按照上述分治法的思想可以看出，要想減少乘法運算次數，關鍵在於計算2個2階方陣的乘積時，能否用少於8次的乘法運算。Strassen提出了一種新的演算法來計算2個2階方陣的乘積。他的演算法只用了7次乘法運算，但增加了加、減法的運算次數。這7次乘法是:　

M1=A11(B12-B22)

M2=(A11+A12)B22

M3=(A21+A22)B11

M4=A22(B21-B11)

M5=(A11+A22)(B11+B22)

M6=(A12-A22)(B21+B22)

M7=(A11-A21)(B11+B12)

做了這7次乘法後，再做若干次加、減法就可以得到:　

C11=M5+M4-M2+M6

C12=M1+M2

C21=M3+M4

C22=M5+M1-M3-M7

以上計算的正確性很容易驗證。例如:　

C22=M5+M1-M3-M7

   =(A11+A22)(B11+B22)+A11(B12-B22)-(A21+A22)B11-(A11-A21)(B11+B12)

   =A11B11+A11B22+A22B11+A22B22+A11B12

     -A11B22-A21B11-A22B11-A11B11-A11B12+A21B11+A21B12

   =A21B12+A22B22　

由(2)式便知其正確性。

至此，我們可以得到完整的Strassen演算法如下：

procedure STRASSEN(n,A,B,C);
begin
  if n=2 then MATRIX-MULTIPLY(A，B，C)
         else begin
                將矩陣A和B依(1)式分塊;
                STRASSEN(n/2,A11,B12-B22,M1);
                STRASSEN(n/2,A11+A12,B22,M2);
                STRASSEN(n/2,A21+A22,B11,M3);
                STRASSEN(n/2,A22,B21-B11,M4);
                STRASSEN(n/2,A11+A22,B11+B22,M5);
                STRASSEN(n/2,A12-A22,B21+B22,M6);
                STRASSEN(n/2,A11-A21,B11+B12,M7);                    ;               end;end;

其中MATRIX-MULTIPLY(A，B，C)是按通常的矩陣乘法計算C=AB的子演算法。

Strassen矩陣乘積分治演算法中，用了7次對於n/2階矩陣乘積的遞迴調用和18次n/2階矩陣的加減運算。由此可知，該演算法的所需的計算時間T(n)滿足如下的遞迴方程:

按照解遞迴方程的套用公式法，其解為T(n)=O(nlog7)≈O(n2.81)。由此可見，Strassen矩陣乘法的計算時間複雜性比普通矩陣乘法有階的改進。

有人曾列舉了計算2個2階矩陣乘法的36種不同方法。但所有的方法都要做7次乘法。除非能找到一種計算2階方陣乘積的演算法，使乘法的計算次數少於7次，按上述思路才有可能進一步改進矩陣乘積的計算時間的上界。但是Hopcroft和Kerr(197l)已經證明，計算2個2×2矩陣的乘積，7次乘法是必要的。因此，要想進一步改進矩陣乘法的時間複雜性，就不能再寄希望於計算2×2矩陣的乘法次數的減少。或許應當研究3×3或5×5矩陣的更好演算法。在Strassen之後又有許多演算法改進了矩陣乘法的計算時間複雜性。目前最好的計算時間上界是O(n2.367)。而目前所知道的矩陣乘法的最好下界仍是它的平凡下界Ω(n2)。因此到目前為止還無法確切知道矩陣乘法的時間複雜性。關於這一研究課題還有許多工作可做。

 

 

 

 

 

 

 

2.快速冪:

unsigned int power(unsigned int n, unsigned int p) </p><p>{ </p><p> // 計算n的p次方 </p><p> unsigned int tmp = 1; </p><p> while (p > 1) </p><p> { </p><p> // 判斷p是否奇數，偶數的最低位必為0 </p><p> if (( p & 1 )!=0) </p><p> { </p><p> tmp *= n; // 若p為奇數，則把“剩下的”乘起來,如果不明白請想想2進位的末尾1是不是應該乘在結果裡面 </p><p> } </p><p> n *= n; </p><p> p >>= 1; </p><p> }</p><p> return n * tmp; // 最後把主體和“剩下的”乘起來作為結果 </p><p>}

圖示:

 

 

 

3.矩陣次冪

 只需注意單位矩陣和零矩陣在數學系統中的對應