首先,我们还是需要假设n是2的幂。将矩阵A,B和C中每一矩阵都分块成为4个大小相等的子矩阵,每个子矩阵都是n/2×n/2的方阵。由此可将方程C=AB重写为:
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
java实现代码如下,仅限于2的n次幂的矩阵
/**
* Strassen矩阵乘法
*
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
16 16 16 16 16 16 16 16
* */
import java.util.*;
public class Strassen {
private static int A[][];
private static int B[][];
private static int C[][];
private final static int NUMBER = 8;
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();
}
}
/**
* 矩阵乘法,此处只是定义了2*2矩阵的乘法
* */
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[][] result) {
for (int i = 0; i < first.length; i++) {
for (int j = 0; j < first[i].length; j++) {
result[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];
}
}
}
/**
* strassen矩阵算法
* */
public void strassen(int[][] A, int[][] B, int[][] C) {
// 定义一些中间变量
int[][] M1 = new int[A.length/2][A.length/2];
int[][] M2 = new int[A.length/2][A.length/2];
int[][] M3 = new int[A.length/2][A.length/2];
int[][] M4 = new int[A.length/2][A.length/2];
int[][] M5 = new int[A.length/2][A.length/2];
int[][] M6 = new int[A.length/2][A.length/2];
int[][] M7 = new int[A.length/2][A.length/2];
int[][] C11 = new int[A.length/2][A.length/2];
int[][] C12 = new int[A.length/2][A.length/2];
int[][] C21 = new int[A.length/2][A.length/2];
int[][] C22 = new int[A.length/2][A.length/2];
int[][] A11 = new int[A.length/2][A.length/2];
int[][] A12 = new int[A.length/2][A.length/2];
int[][] A21 = new int[A.length/2][A.length/2];
int[][] A22 = new int[A.length/2][A.length/2];
int[][] B11 = new int[A.length/2][A.length/2];
int[][] B12 = new int[A.length/2][A.length/2];
int[][] B21 = new int[A.length/2][A.length/2];
int[][] B22 = new int[A.length/2][A.length/2];
int[][] temp = new int[A.length/2][A.length/2];
int[][] temp1 = new int[A.length/2][A.length/2];
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);
strassen(A11, temp, M1);
// 计算M2
Add(A11, A12, temp);
strassen(temp, B22, M2);
// 计算M3
Add(A21, A22, temp);
strassen(temp, B11, M3);
// M4
sub(B21, B11, temp);
strassen(A22, temp, M4);
// M5
Add(A11, A22, temp1);
Add(B11, B22, temp);
strassen(temp1, temp, M5);
// M6
sub(A12, A22, temp1);
Add(B21, B22, temp);
strassen(temp1, temp, M6);
// M7
sub(A11, A21, temp1);
Add(B11, B12, temp);
strassen(temp1, temp, M7);
// 计算C11 = M5+M4-M2+M6
Add(M5, M4, temp1);
sub(temp1, M2, temp);
Add(temp, M6, C11);
// 计算C12 = M1+M2
Add(M1, M2, C12);
// C21 = M3+M4
Add(M3, M4, C21);
// C22 = M5+M1-M3-M7
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);
}
}