1.矩阵相乘的朴素算法
时间复杂度T(n)=Θ(n3),朴素矩阵相乘算法,思想明了,编程实现简单。时间复杂度是Θ(n^3)。伪码如下
for i ← 1 to n
do for j ← 1 to n
do c[i][j] ← 0
for k ← 1 to n
do c[i][j] ← c[i][j] + a[i][k]⋅ b[k][j]
2.strassen算法
矩阵乘法是典型采用了分治思想的算法,将一个 n × n的矩阵乘法分割为 2 × 2个(n/2) × (n/2)大小的子矩阵相乘。
P1=a ⋅(f−h)
P2= (a+ b)⋅h
P3= (c+ d)⋅e
P4=d ⋅(g−e)
P5= (a+ d)⋅(e+h)
P6= (b−d)⋅(g+h)
P7= (a−c)⋅(e+f )
r =P5+ P4−P2+P6
s =P1+ P2
t =P3+ P4
u =P5+ P1−P3−P7
在这样的加乘中,进行了7次乘法操作,而常规的朴素算法是8次乘法操作,在矩阵很大的时候8次乘法的复杂度远高于7次,关于strassen算法的正确性在算法导论的分治这节给出了证明。
伪代码如下:
Strassen (N,MatrixA,MatrixB,MatrixResult)
//splitting input Matrixes, into 4 submatrices each.
for i <- 0 to N/2
for j <- 0 to N/2
A11[i][j] <- MatrixA[i][j]; //a矩阵块
A12[i][j] <- MatrixA[i][j + N / 2]; //b矩阵块
A21[i][j] <- MatrixA[i + N / 2][j]; //c矩阵块
A22[i][j] <- MatrixA[i + N / 2][j + N / 2];//d矩阵块
B11[i][j] <- MatrixB[i][j]; //e 矩阵块
B12[i][j] <- MatrixB[i][j + N / 2]; //f 矩阵块
B21[i][j] <- MatrixB[i + N / 2][j]; //g 矩阵块
B22[i][j] <- MatrixB[i + N / 2][j + N / 2]; //h矩阵块
//here we calculate M1..M7 matrices .
//递归求M1
HalfSize <- N/2
AResult <- A11+A22
BResult <- B11+B22
Strassen( HalfSize, AResult, BResult, M1 ); //M1=(A11+A22)*(B11+B22)p5=(a+d)*(e+h)
//递归求M2
AResult <- A21+A22
Strassen(HalfSize, AResult, B11, M2); //M2=(A21+A22)B11 p3=(c+d)*e
//递归求M3
BResult <- B12 - B22
Strassen(HalfSize, A11, BResult, M3); //M3=A11(B12-B22) p1=a*(f-h)
//递归求M4
BResult <- B21 - B11
Strassen(HalfSize, A22, BResult, M4); //M4=A22(B21-B11) p4=d*(g-e)
//递归求M5
AResult <- A11+A12
Strassen(HalfSize, AResult, B22, M5); //M5=(A11+A12)B22 p2=(a+b)*h
//递归求M6
AResult <- A21-A11
BResult <- B11+B12
Strassen( HalfSize, AResult, BResult, M6); //M6=(A21-A11)(B11+B12) p7=(c-a)(e+f)
//递归求
M7 AResult <- A12-A22
BResult <- B21+B22
Strassen(HalfSize, AResult, BResult, M7); //M7=(A12-A22)(B21+B22) p6=(b-d)*(g+h)
//计算结果子矩阵
C11 <- M1 + M4 - M5 + M7;
C12 <- M3 + M5;
C21 <- M2 + M4;
C22 <- M1 + M3 - M2 + M6;
//at this point , we have calculated the c11..c22 matrices, and now we are going to //put them together and make a unit matrix which would describe our resulting Matrix.
for i <- 0 to N/2
for j <- 0 to N/2
MatrixResult[i][j] <- C11[i][j];
MatrixResult[i][j + N / 2] <- C12[i][j];
MatrixResult[i + N / 2][j] <- C21[i][j];
MatrixResult[i + N / 2][j + N / 2] <- C22[i][j];
性能分析:
矩阵大小 | 朴素矩阵算法(秒) | Strassen算法(秒) |
32 | 0.003 | 0.003 |
64 | 0.004 | 0.004 |
128 | 0.021 | 0.071 |
256 | 0.09 | 0.854 |
512 | 0.782 | 6.408 |
1024 | 8.908 | 52.391 |
可以看到使用Strassen算法时,耗时不但没有减少,反而剧烈增多,在n=512时计算时间就无法忍受,效果没有朴素矩阵算法好
1)采用Strassen算法作递归运算,需要创建大量的动态二维数组,其中分配堆内存空间将占用大量计算时间,从而掩盖了Strassen算法的优势
2)于是对Strassen算法做出改进,设定一个界限。当n<界限时,使用普通法计算矩阵,而不继续分治递归。需要合理设置界限,不同环境(硬件配置)下界限不同