Strassen矩陣乘法 分治與遞迴
轉自：http://blog.sina.com.cn/s/blog_7e9a88f70100zj2h.html
Strassen矩陣乘法
矩陣乘法是線性代數中最常見的運算之一，它在數值計算中有廣泛的應用。若A和B是2個n×n的矩陣，則它們的乘積C=AB同樣是一個n×n的矩陣。A和B的乘積矩陣C中的元素C[i,j]定義為:
分治與遞迴--strassen矩陣乘法
若依此定義來計算A和B的乘積矩陣C，則每計算C的一個元素C[i,j]，需要做n個乘法和n-1次加法。因此，求出矩陣C的n2個元素所需的計算時間為0(n3)。
60年代末，Strassen採用了類似於在大整數乘法中用過的分治技術，將計算2個n階矩陣乘積所需的計算時間改進到
首先，我們還是需要假設n是2的冪。將矩陣A，B和C中每一矩陣都分塊成為4個大小相等的子矩陣，每個子矩陣都是n/2×n/2的方陣。由此可將方程C=AB重寫為:
分治與遞迴--strassen矩陣乘法(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)應該滿足：
分治與遞迴--strassen矩陣乘法
這個仍然是T(n)=O(n3)。因此，該方法並不比用原始定義直接計算更有效。究其原因，乃是由於式(2)-(5)並沒有減少矩陣的乘法次數。而矩陣乘法耗費的時間要比矩陣加減法耗費的時間多得多。要想改進矩陣乘法的計算時間複雜性，必須減少子矩陣乘法運算的次數。按照上述
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);
分治與遞迴--strassen矩陣乘法;
end;
end;
其中MATRIX-MULTIPLY(A，B，C)是按通常的矩陣乘法計算C=AB的子演算法。
Strassen矩陣乘積分治演算法中，用了7次對於n/2階矩陣乘積的遞迴呼叫和18次n/2階矩陣的加減運算。由此可知，該演算法的所需的計算時間T(n)滿足如下的遞迴方程:
分治與遞迴--strassen矩陣乘法
按照解遞迴方程的套用公式法，其解為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)。因此到目前為止還無法確切知道矩陣乘法的時間複雜性。關於這一研究課題還有許多工作可做。
importjava.util.*;
publicclass Strassen{
public Strassen(){
A = new int[NUMBER][NUMBER];
B = new int[NUMBER][NUMBER];
C = new int[NUMBER][NUMBER];
}
public void input(int a[][]){
Scanner scanner = new Scanner(System.in);
for(int i = 0; i < a.length; i++) {
for(int j = 0; j < a[i].length; j++) {
a[i][j] = scanner.nextInt();
}
}
}
public void output(int[][] resault){
for(int b[] : resault) {
for(int temp : b) {
System.out.print(temp + " ");
}
System.out.println();
}
}
public void Mul(int[][] first, int[][] second, int[][]resault){
for(int i = 0; i < 2; ++i) {
for(int j = 0; j < 2; ++j) {
resault[i][j] = 0;
for(int k = 0; k < 2; ++k) {
resault[i][j] += first[i][k] * second[k][j];
}
}
}
}
public void Add(int[][] first, int[][] second, int[][]resault){
for(int i = 0; i < first.length; i++) {
for(int j = 0; j < first[i].length; j++){
resault[i][j] = first[i][j] + second[i][j];
}
}
}
public void sub(int[][] first, int[][] second, int[][]resault){
for(int i = 0; i < first.length; i++) {
for(int j = 0; j < first[i].length; j++){
resault[i][j] = first[i][j] - second[i][j];
}
}
}
public void strassen(int[][] A, int[][] B, int[][] C){
//定義一些中間變數
int [][] M1=new int [NUMBER][NUMBER];
int [][] M2=new int [NUMBER][NUMBER];
int [][] M3=new int [NUMBER][NUMBER];
int [][] M4=new int [NUMBER][NUMBER];
int [][] M5=new int [NUMBER][NUMBER];
int [][] M6=new int [NUMBER][NUMBER];
int [][] M7=new int [NUMBER][NUMBER];
int [][] C11=new int [NUMBER][NUMBER];
int [][] C12=new int [NUMBER][NUMBER];
int [][] C21=new int [NUMBER][NUMBER];
int [][] C22=new int [NUMBER][NUMBER];
int [][] A11=new int [NUMBER][NUMBER];
int [][] A12=new int [NUMBER][NUMBER];
int [][] A21=new int [NUMBER][NUMBER];
int [][] A22=new int [NUMBER][NUMBER];
int [][] B11=new int [NUMBER][NUMBER];
int [][] B12=new int [NUMBER][NUMBER];
int [][] B21=new int [NUMBER][NUMBER];
int [][] B22=new int [NUMBER][NUMBER];
int [][] temp=new int [NUMBER][NUMBER];
int [][] temp1=new int [NUMBER][NUMBER];
if(A.length==2){
Mul(A, B, C);
}else{
//首先將矩陣A，B 分為4塊
for(int i = 0; i < A.length/2; i++) {
for(int j = 0; j < A.length/2; j++) {
A11[i][j]=A[i][j];
A12[i][j]=A[i][j+A.length/2];
A21[i][j]=A[i+A.length/2][j];
A22[i][j]=A[i+A.length/2][j+A.length/2];
B11[i][j]=B[i][j];
B12[i][j]=B[i][j+A.length/2];
B21[i][j]=B[i+A.length/2][j];
B22[i][j]=B[i+A.length/2][j+A.length/2];
}
}
//計算M1
sub(B12, B22, temp);
Mul(A11, temp, M1);
//計算M2
Add(A11, A12, temp);
Mul(temp, B22, M2);
//計算M3
Add(A21, A22, temp);
Mul(temp, B11, M3);
//M4
sub(B21, B11, temp);
Mul(A22, temp, M4);
//M5
Add(A11, A22, temp1);
Add(B11, B22, temp);
Mul(temp1, temp, M5);
//M6
sub(A12, A22, temp1);
Add(B21, B22, temp);
Mul(temp1, temp, M6);
//M7
sub(A11, A21, temp1);
Add(B11, B12, temp);
Mul(temp1, temp, M7);
//計算C11
Add(M5, M4, temp1);
sub(temp1, M2, temp);
Add(temp, M6, C11);
//計算C12
Add(M1, M2, C12);
//C21
Add(M3, M4, C21);
//C22
Add(M5, M1, temp1);
sub(temp1, M3, temp);
sub(temp, M7, C22);
//結果送回C中
for(int i = 0; i < C.length/2; i++) {
for(int j = 0; j < C.length/2; j++) {
C[i][j]=C11[i][j];
C[i][j+C.length/2]=C12[i][j];
C[i+C.length/2][j]=C21[i][j];
C[i+C.length/2][j+C.length/2]=C22[i][j];
}
}
}
}
public static void main(String[] args){
Strassen demo=new Strassen();
System.out.println("輸入矩陣A");
demo.input(A);
System.out.println("輸入矩陣B");
demo.input(B);
demo.strassen(A, B, C);
demo.output(C);
}
private static int A[][];
private static int B[][];
private static int C[][];
private final static int NUMBER = 4;
}