矩阵乘法的算法实现
2016-11-30 18:25
330 查看
一般矩阵乘法算法:
原理:矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的栏数(column)和第二个矩阵的列数(row)相同时才有定义。一般单指矩阵乘积时,指的便是一般矩阵乘积。若A为m×n矩阵,B为n×p矩阵,则他们的乘积AB会是一个m×p矩阵。其乘积矩阵的元素如下面式子得出:代码如下:
struct mat{ int n, m; double data[MAXN][MAXN]; }; int mul(mat& c, const mat& a, const mat& b){ int i, j, k; if (a.m != b.n) return 0; c.n = a.n; c.m = b.m; for (i = 0; i < c.n; i++) for (j = 0; j < c.m; j++) for (c.data[i][j] = k = 0; k < a.m; k++) c.data[i][j] += a.data[i][k] * b.data[k][j]; return 1; }
时间复杂度:O(n3)
传统分治方法:
C = AB将A, B, C分成相等大小的方块矩阵:
C11=A11B11+A12B21
(2)
C12=A11B12+A12B22
(3)
C21=A21B11+A22B21
(4)
C22=A21B12+A22B22
(5)
如果n=2,则2个2阶方阵的乘积可以直接用(2)-(5)式计算出来,共需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(2) = b
T(n) = 8T(n / 2) + cn2 (n > 2)
由上式可知:分治法的运用,方阵的乘法的算法效率并没有提高!
分治法的设计思想:将一个难以直接解决的大问题,分割成一些规模较小的相同问题,以便各个击破,分而治之。
任何一个可以用计算机求解的问题所需的计算时间都与其规模有关。问题的规模越小,越容易直接求解,解题所需的计算时间也越少。例如,对于n个元素的排。n=2时,只要作一次比较即可排好序。n=3时序问题,当n=1时,不需任何计算只要作3次比较即可,…。而当n较大时,问题就不那么容易处理了。要想直接解决一个规模较大的问题,有时是相当困难的。
分治法所能解决的问题一般具有以下几个特征:
1.可缩性。问题的规模缩小到一定的程度就可以容易地解决;
2.最有子结构性。问题可以分解为若干个规模较小的相同问题;
3.可合性。利用该问题分解出的子问题的解可以合并为该问题的解;
4.独立性。该问题所分解出的各个子问题是相互独立的,即子问 题之间不包含公共的子子问题。
分治思想与递归就像一对孪生兄弟,经常同时应用在算法设计中,并由此产生高效的算法!
Strassen算法
传统方法2个2阶方阵相乘需要8次乘法。按照上述分治法的思想可以看出,要想减少乘法运算次数,关键在于计算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
算法效率分析:
Strassen矩阵乘积分治算法中,用了7次对于n/2阶矩阵乘积的递归调用和18次n/2阶矩阵的加减运算。由此可知,该算法的所需的计算时间T(n)满足如下的递归方程:
T(2) = b
T(n) = 7T(n / 2) + cn2 (n > 2)
按照解递归方程的套用公式法,其解为T(n)=O(nlog7)≈O(n2.81)。由此可见,Strassen矩阵乘法的计算时间复杂性比普通矩阵乘法有阶的改进。
数据的存储结构:
二维数组(弊端:必须宏定义一常量N)
二重指针 (摆脱宏定义的限制,但程序变得繁琐)
原始矩阵:
A,B,C(C=AB)
分块矩阵:
A11,A12,A21,A22(对方阵A分成四块)
B11,B12,B21,B22(对方阵B分成四块)
七块关键阵:
M1,M2,M3,M4,M5,M6,M7
中间方阵:
AA,BB(计算加法时作一个桥梁)
二级指针定义说明:
int **c; c = (int**)malloc(n * sizeof(int)); for(int i = 0; i < n; i++) c[i] = (int *)malloc(n * sizeof(int));
两矩阵相加问题:
//矩阵加法: void add(int n, int A[] , int B[] , int R[] ) { int i, j; for(i = 0; i < n; i++) for(j = 0; j < n; j++) R[i][j] = A[i][j] + B[i][j]; }
两矩阵相减问题:
//矩阵减法: void sub(int n,int A[] ,int B[] ,int R[] ) { int i,j; for(i = 0; i < n; i++) for(j = 0; j < n; j++) R[i][j] = A[i][j] - B[i][j]; }
两方阵相乘问题:
//阶为2的方阵相乘: void multiply(int A[] ,int B[] ,int R[] ) { int M1 = A[0][1] * (B[0][1] - B[1][1]); int M2 = (A[0][0] + A[0][1]) * B[1][1]; int M3 = (A[1][0] + A[1][1]) * B[0][0]; int M4 = A[1][1] * (B[1][0] - B[0][0]); int M5 = (A[0][0] + A[1][1]) * (B[0][0] + B[1][1]); int M6 = (A[0][1] - A[1][1]) * (B[1][0] + B[1][1]); int M7 = (A[0][0] - A[1][0]) * (B[0][0] + B[0][1]); R[0][0] = M5 + M4 - M2 + M6; R[0][1] = M1 + M2; R[1][0] = M3 + M4; R[1][1] = M5 + M1 - M3 - M7; }
将大的方阵划分的问题
//分块函数 void divide(int A[] , int n, int X11[] , int X12[] , int X21[] , int X22[] ) { int i, j; for(i = 0; i < (n / 2); i++) for(j = 0; j < (n / 2); j++) { X11[i][j] = A[i][j]; X12[i][j] = A[i][j + (n / 2)]; X21[i][j] = A[i + (n / 2)][j]; X22[i][j] = A[i + (n / 2)][j + (n / 2)]; } }
完整代码如下:
#include<stdio.h>
#include<math.h>
#define N 4
//二阶矩阵相乘
void Matrix_Multiply(int A[]
, int B[]
, int C[]
) {
for(int i = 0; i < 2; i++) {
for(int j = 0; j < 2; j++) {
C[i][j] = 0;
for(int t = 0; t < 2; t++) {
C[i][j] = C[i][j] + A[i][t] * B[t][j];
}
}
}
}
//矩阵加法: void add(int n, int A[] , int B[] , int R[] ) { int i, j; for(i = 0; i < n; i++) for(j = 0; j < n; j++) R[i][j] = A[i][j] + B[i][j]; }
//矩阵减法:
void sub(int n, int A[]
, int B[]
, int R[]
)
{
int i,j;
for(i = 0; i < n; i++)
for(j = 0; j < n; j++)
R[i][j] = A[i][j] - B[i][j];
}
void strassen(int n, int A[]
, int B[]
, int C[]
)
{
int i, j;
int A11
, A12
, A21
, A22
;
int B11
, B12
, B21
, B22
;
int C11
, C12
, C21
, C22
;
int AA
, BB
;
int M1
, M2
, M3
, M4
, M5
, M6
, M7
;
if(n == 2) {
Matrix_Multiply(A, B, C);
} else {
for(i = 0; i < n / 2; i++) {
for(j = 0; j < n / 2; j++) {
A11[i][j] = A[i][j];
A12[i][j] = A[i][j + n / 2];
A21[i][j] = A[i + n / 2][j];
A22[i][j] = A[i + n / 2][j + n / 2];
B11[i][j] = B[i][j];
B12[i][j] = B[i][j + n / 2];
B21[i][j] = B[i + n /2][j];
B22[i][j] = B[i + n /2][j + n / 2];
}
}
sub(n / 2, B12, B22, BB);
strassen(n / 2, A11, BB, M1);
add(n / 2, A11, A12, AA);
strassen(n / 2, AA, B22, M2);
add(n / 2, A21, A22, AA);
strassen(n / 2, AA, B11, M3);
sub(n / 2, B21, B11, BB);
strassen(n / 2, A22, BB, M4);
add(n / 2, A11, A22, AA);
add(n / 2, B11, B22, BB);
strassen(n / 2, AA, BB, M5);
sub(n / 2, A12, A22, AA);
add(n / 2, B21, B22, BB);
strassen(n / 2, AA, BB, M6);
sub(n / 2, A11, A21, AA);
add(n / 2, B11, B12, BB);
strassen(n / 2, AA, BB, M7);
//C11 = M5 + M4 - M2 + M6
add(n / 2, M5, M4, AA);
sub(n / 2, M6, M2, BB);
add(n / 2, AA, BB, C11);
//C12 = M1 + M2
add(n / 2, M1, M2, C12);
//C21 = M3 + M4
add(n / 2, M3, M4, C21);
//C22 = M5 + M1 - M3 - M7
sub(n / 2, M5, M3, AA);
sub(n / 2, M1, M7, BB);
add(n / 2, AA, BB, C22);
for(i = 0; i < n / 2; i++) {
for(j = 0; j < n / 2; j++) {
C[i][j] = C11[i][j];
C[i][j + n / 2] = C12[i][j];
C[i + n / 2][j] = C21[i][j];
C[i + n / 2][j + n / 2] = C22[i][j];
}
}
}
}
int main(void)
{
int A
, B
, C
;
printf("input A: \n");
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
scanf("%d", &A[i][j]);
printf("input B: \n");
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++)
scanf("%d", &B[i][j]);
strassen(N, A, B, C);
printf("C:\n");
for(int i = 0; i < N; i++)
for(int j = 0; j < N; j++) {
printf("%d ", C[i][j]);
if(j > 0 && j % 3 == 0)
printf("\n");
}
return 0;
}
Strassen算法虽然把时间代价从O(n3)降到O(n2.81),但是进一步分析,后者的系数比1大很多,当n比较小(如n<45)且矩阵中非零元素较少时,不宜采用此方法。Strassen算法仅当n很大的时候才有价值,因此,可以说这方面的研究更偏重于理论价值。
相关文章推荐
- 一些算法的MapReduce实现——矩阵-向量乘法实现
- 【数据结构与算法】数组应用1:矩阵乘法(Java实现)
- 【甘道夫】Mapreduce实现矩阵乘法的算法思路
- MPI群通信与矩阵乘法的Fox算法实现
- 转置矩阵的分块并行乘法(C语言实现),计算矩阵C[rawn][rawn]=A[rawm][rawn]'*B[rawm][rawn],子块大小为S*T,其算法实现原理参加本代码的附件。
- 【甘道夫】Mapreduce实现矩阵乘法的算法思路
- MPI群通信与矩阵乘法的Fox算法实现
- 基于矩阵分解推荐算法之交替最小二乘法(ALS)--附实现代码
- 实现两个N*N矩阵的乘法,矩阵由一维数组表示
- 快速立即除法的乘法实现(通用算法)
- 求助:矩阵求逆、矩阵乘法的Delphi实现(20个积分的奖励)
- 递归+优化 一维数组实现矩阵乘法
- Strassen矩阵乘法 + 快速计算乘方的算法 + 矩阵的次幂
- 矩阵乘法 C++实现
- hadoop/mapreduce-矩阵乘法的实现代码下载
- MapReduce框架中矩阵相乘的算法思路及其实现
- 实现两个N*N矩阵的乘法,矩阵由一维数组表示
- 矩阵链乘法(最小乘法数) (动态规划算法实现)算法导论p201
- 实现N*N矩阵的乘法,矩阵有一维数组表示
- 矩阵乘法算法