strassen算法

❶ C语言实现Strassen算法计算两个矩阵积

        如果你学过线性代数,说到计算矩阵乘法,那么我们一般常规操作一般公式是 ;即A的每一行乘B的每一列，但是需要注意的是A的列数要与B的行数要相等，否则A与B不可做计算操作，对应的项乘积进行相加处理，就可得到 储存在在C中对应的位置，直到每一项被计算完成。这便是矩阵乘法的定义。

        上面是对矩阵成积计算方法定义的陈述，但是如果要把它转为程序，计算机读懂的语言，我们首先会想到使用三层for循环，一层用于 赋值操作，两层用于计算 ，那么程序基本会定义为如下程序：

for(i=0;i<n;i++)

    for(j=0;j<n;j++)

        for(m=0;m<n;m++)

            +=

        由于有三个for循环，那么这个计算方法时间复杂度 ，这个计算方法效率相对有些低了，为此出现了Strassen计算方法，但是这个算法的要求是两个矩阵的size ，n要满足： 。Strassen 算法类似于分治的思想，先将矩阵分解为更小的矩阵，直到n=1，即矩阵的size ，再计算两个数乘积返回数值即可。

        根据Strassen 计算法则可以知道先将矩阵A，B进行拆分，使得子矩阵size满足 ，如下所示：

根据Strassen计算法则有以下公式要计算.

        根据上诉计算规则,我们只要对分解到符合条件的矩阵,都可使用上诉公式进行计算.分解示例图如下：

        说到这里，想必你应该对Strassen计算矩阵规则有所了解了，其实就是对矩阵分解后作加，减，乘三种运算，那么接下去就是用程序实现了。源程序如下所示：

#include <stdio.h>

#include<stdlib.h>

/******************************************

*

* function add_matrix()

*

* args sub_ma1,sub_ma2 inttype **pointers

* n the size of sub_ma1 sub_ma2 n*n

*

* return matrix **inttype

*

* *****************************************/

int** add_matrix(int** sub_ma1, int** sub_ma2, int n) {

    int** temp = init_matrix(n);

    for(int i=0; i<n; i++)

        for(int j=0; j<n; j++)

            temp[i][j] = sub_ma1[i][j] + sub_ma2[i][j];

    return temp;

}

/******************************************

*

* function subtract_matrix()

*

* args sub_ma1,sub_ma2 inttype **pointers

* n the size of sub_ma1 sub_ma2 n*n

*

* return matrix **inttype

*

* *****************************************/

int** subtract_matrix(int** sub_ma1, int** sub_ma2, int n) {

    int** temp = init_matrix(n);

    for(int i=0; i<n; i++)

        for(int j=0; j<n; j++)

            temp[i][j] = sub_ma1[i][j] - sub_ma2[i][j];

    return temp;

}

/**********************************************

*

* function square_matrix_mutiply_recursive()

*

* args

* A,B, inttype ** pointer

* n inttype matrix size n*n

*

* return C inttype array

*

************************************************/

int** square_matrix_strassen_mutiply(int **A,int **B,int n){

    int i,j;

    //only one element

    if (n == 1) {

        int** C = init_matrix(1);

        C[0][0] = A[0][0] * B[0][0];

        return C;

    }

    else{

        //init C,A,B

        int** C = init_matrix(n);

        int k = n/2;

        int** A11 = init_matrix(k);

        int** A12 = init_matrix(k);

        int** A21 = init_matrix(k);

        int** A22 = init_matrix(k);

        int** B11 = init_matrix(k);

        int** B12 = init_matrix(k);

        int** B21 = init_matrix(k);

        int** B22 = init_matrix(k);

        //resolve A,B matrixs to A11...A22,B11...B22

        for(i=0; i<k; i++)

            for(j=0; j<k; j++) {

                A11[i][j] = A[i][j];

                A12[i][j] = A[i][k+j];

                A21[i][j] = A[k+i][j];

                A22[i][j] = A[k+i][k+j];

                B11[i][j] = B[i][j];

                B12[i][j] = B[i][k+j];

                B21[i][j] = B[k+i][j];

                B22[i][j] = B[k+i][k+j];

            }

        //calculate P[1-7]

        int** P1 = square_matrix_strassen_mutiply(A11, subtract_matrix(B12, B22, k), k);

        int** P2 = square_matrix_strassen_mutiply(add_matrix(A11, A12, k), B22, k);

        int** P3 = square_matrix_strassen_mutiply(add_matrix(A21, A22, k), B11, k);

        int** P4 = square_matrix_strassen_mutiply(A22, subtract_matrix(B21, B11, k), k);

        int** P5 = square_matrix_strassen_mutiply(add_matrix(A11, A22, k), add_matrix(B11, B22, k), k);

        int** P6 = square_matrix_strassen_mutiply(subtract_matrix(A12, A22, k), add_matrix(B21, B22, k), k);

        int** P7 = square_matrix_strassen_mutiply(subtract_matrix(A11, A21, k), add_matrix(B11, B12, k), k);

        //calculate C11.....C22

        int** C11 = subtract_matrix(add_matrix(add_matrix(P5, P4, k), P6, k), P2, k);

        int** C12 = add_matrix(P1, P2, k);

        int** C21 = add_matrix(P3, P4, k);

        int** C22 = subtract_matrix(subtract_matrix(add_matrix(P5, P1, k), P3, k), P7, k);

        // C11,C12,C13,C14 to C

        for(i=0; i<k; i++)

            for(j=0; j<k; j++) {

                C[i][j] = C11[i][j];

                C[i][j+k] = C12[i][j];

                C[k+i][j] = C21[i][j];

                C[k+i][k+j] = C22[i][j];

            }

        //free memory

        for(i=0;i<k;i++){

            // free subarrays of A,B

            free(A11[i]);

            free(A12[i]);

            free(A21[i]);

            free(A22[i]);

            free(B11[i]);

            free(B12[i]);

            free(B21[i]);

            free(B22[i]);

            //free subarray of P

            free(P1[i]);

            free(P2[i]);

            free(P3[i]);

            free(P4[i]);

            free(P5[i]);

            free(P6[i]);

            free(P7[i]);

            //free subarray of C

            free(C11[i]);

            free(C12[i]);

            free(C21[i]);

            free(C22[i]);

        }

        //free rows of A,B

        free(A11);

        free(A12);

        free(A21);

        free(A22);

        free(B11);

        free(B12);

        free(B21);

        free(B22);

        //free rows of P

        free(P1);

        free(P2);

        free(P3);

        free(P4);

        free(P5);

        free(P6);

        free(P7);

        //free rows of C

        free(C11);

        free(C12);

        free(C21);

        free(C22);

        //make NULL

        A11=NULL;

        A12=NULL;

        A21=NULL;

        A22=NULL;

        B11=NULL;

        B12=NULL;

        B21=NULL;

        B22=NULL;

        P1=NULL;

        P2=NULL;

        P3=NULL;

        P4=NULL;

        P5=NULL;

        P6=NULL;

        P7=NULL;

        C11=NULL;

        C12=NULL;

        C21=NULL;

        C22=NULL;

        return C;

    }

}

int main()

{

    int n=2;

    int i,j;

    // init A,B

    int**A=init_matrix(2);

    int**B=init_matrix(2);

    A[0][0]=1;

    A[0][1]=3;

    A[1][0]=7;

    A[1][1]=5;

    B[0][0]=6;

    B[0][1]=8;

    B[1][0]=4;

    B[1][1]=2;

    //use square_matrix_strassen_mutiply()

    int **C=square_matrix_strassen_mutiply(A,B,n);

    //output C

    printf("the result of C is:\n");

    for(i=0;i<n;i++){

        for(j=0;j<n;j++)

            printf("%d ",*(*(C+i)+j));

        printf("\n"); //next line

    }

    return 0;

}

/****************************************

*

* input A={{1,3},{7,5}} B={{6,8},{4,2}}

*

* output C={{18,14},{62,66}};

*

* **************************************/

以上就是对Strassen算法实现的C源程序。

运行结果截图如下：

有兴趣的朋友，可以一起交流，共同进步。

❷ 矩阵乘法的strassen算法的直观理解

如果想要直观的理解,那么好的方式是用图形代替代数公式,下面的内容就是基于这个想法
1.代数公式和图形的联系
对于a b,我们可以把它看作是一个线段,a和b是它的两个端点
那么a b+a*c,即a(b+c)就是由两个线段组成的,它们共同的点是a
进一步,(a+b)(c+d)可以看作是一个四边形,四点是a,b,c,d,四条边就是它们的项.如下图:

当然a,b,c,d的正负都行
2.strassen算法的直观效果
根据以上的想法,我们可以把strassen算法过程画成以下的效果

其中实线代表我们想要得到的矩阵乘法的项,如下
C11​=A11​⋅B11​+A12​⋅B21​
C12​=A11​⋅B12​+A12​⋅B21
​C21​=A21​⋅B11​+A22​⋅B21
​C22​=A21​⋅B12​+A22​⋅B22​

A22B22是两个四边形的重复线段,可以使它在两个四边形的正负不同相抵
A11B22和B22A12是由两条线段组成,可以看作是B22(A11+A12)
同样B11A22,A22B21可以看作是A22(B11+B21)
所以最终是由(A22+A12)(B21+B22)+(A11+A22)(B11+B22)+B22(A11+A12)+A22(B11+B21)通过调整正负号可以得到A11​⋅B11​+A12​⋅B21

通过类比,有时能够更方便的理解和记忆

❸ 大数据最常用的算法有哪些

奥地利符号计算研究所(Research Institute for Symbolic Computation，简称RISC)的Christoph Koutschan博士在自己的页面上发布了一篇文章，提到他做了一个调查，参与者大多数是计算机科学家，他请这些科学家投票选出最重要的算法，以下是这次调查的结果，按照英文名称字母顺序排序。

大数据等最核心的关键技术：32个算法

1、A* 搜索算法——图形搜索算法，从给定起点到给定终点计算出路径。其中使用了一种启发式的估算，为每个节点估算通过该节点的最佳路径，并以之为各个地点排定次序。算法以得到的次序访问这些节点。因此，A*搜索算法是最佳优先搜索的范例。

2、集束搜索(又名定向搜索，Beam Search)——最佳优先搜索算法的优化。使用启发式函数评估它检查的每个节点的能力。不过，集束搜索只能在每个深度中发现最前面的m个最符合条件的节点，m是固定数字——集束的宽度。

3、二分查找(Binary Search)——在线性数组中找特定值的算法，每个步骤去掉一半不符合要求的数据。

4、分支界定算法(Branch and Bound)——在多种最优化问题中寻找特定最优化解决方案的算法，特别是针对离散、组合的最优化。

5、Buchberger算法——一种数学算法，可将其视为针对单变量最大公约数求解的欧几里得算法和线性系统中高斯消元法的泛化。

6、数据压缩——采取特定编码方案，使用更少的字节数(或是其他信息承载单元)对信息编码的过程，又叫来源编码。

7、Diffie-Hellman密钥交换算法——一种加密协议，允许双方在事先不了解对方的情况下，在不安全的通信信道中，共同建立共享密钥。该密钥以后可与一个对称密码一起，加密后续通讯。

8、Dijkstra算法——针对没有负值权重边的有向图，计算其中的单一起点最短算法。

9、离散微分算法(Discrete differentiation)。

10、动态规划算法(Dynamic Programming)——展示互相覆盖的子问题和最优子架构算法

11、欧几里得算法(Euclidean algorithm)——计算两个整数的最大公约数。最古老的算法之一，出现在公元前300前欧几里得的《几何原本》。

12、期望-最大算法(Expectation-maximization algorithm，又名EM-Training)——在统计计算中，期望-最大算法在概率模型中寻找可能性最大的参数估算值，其中模型依赖于未发现的潜在变量。EM在两个步骤中交替计算，第一步是计算期望，利用对隐藏变量的现有估计值，计算其最大可能估计值;第二步是最大化，最大化在第一步上求得的最大可能值来计算参数的值。

13、快速傅里叶变换(Fast Fourier transform，FFT)——计算离散的傅里叶变换(DFT)及其反转。该算法应用范围很广，从数字信号处理到解决偏微分方程，到快速计算大整数乘积。

14、梯度下降(Gradient descent)——一种数学上的最优化算法。

15、哈希算法(Hashing)。

16、堆排序(Heaps)。

17、Karatsuba乘法——需要完成上千位整数的乘法的系统中使用，比如计算机代数系统和大数程序库，如果使用长乘法，速度太慢。该算法发现于1962年。

18、LLL算法(Lenstra-Lenstra-Lovasz lattice rection)——以格规约(lattice)基数为输入，输出短正交向量基数。LLL算法在以下公共密钥加密方法中有大量使用：背包加密系统(knapsack)、有特定设置的RSA加密等等。

19、最大流量算法(Maximum flow)——该算法试图从一个流量网络中找到最大的流。它优势被定义为找到这样一个流的值。最大流问题可以看作更复杂的网络流问题的特定情况。最大流与网络中的界面有关，这就是最大流-最小截定理(Max-flow min-cut theorem)。Ford-Fulkerson 能找到一个流网络中的最大流。

20、合并排序(Merge Sort)。

21、牛顿法(Newton’s method)——求非线性方程(组)零点的一种重要的迭代法。

22、Q-learning学习算法——这是一种通过学习动作值函数(action-value function)完成的强化学习算法，函数采取在给定状态的给定动作，并计算出期望的效用价值，在此后遵循固定的策略。Q-leanring的优势是，在不需要环境模型的情况下，可以对比可采纳行动的期望效用。

23、两次筛法(Quadratic Sieve)——现代整数因子分解算法，在实践中，是目前已知第二快的此类算法(仅次于数域筛法Number Field Sieve)。对于110位以下的十位整数，它仍是最快的，而且都认为它比数域筛法更简单。

24、RANSAC——是“RANdom SAmple Consensus”的缩写。该算法根据一系列观察得到的数据，数据中包含异常值，估算一个数学模型的参数值。其基本假设是：数据包含非异化值，也就是能够通过某些模型参数解释的值，异化值就是那些不符合模型的数据点。

25、RSA——公钥加密算法。首个适用于以签名作为加密的算法。RSA在电商行业中仍大规模使用，大家也相信它有足够安全长度的公钥。

26、Sch?nhage-Strassen算法——在数学中，Sch?nhage-Strassen算法是用来完成大整数的乘法的快速渐近算法。其算法复杂度为：O(N log(N) log(log(N)))，该算法使用了傅里叶变换。

27、单纯型算法(Simplex Algorithm)——在数学的优化理论中，单纯型算法是常用的技术，用来找到线性规划问题的数值解。线性规划问题包括在一组实变量上的一系列线性不等式组，以及一个等待最大化(或最小化)的固定线性函数。

28、奇异值分解(Singular value decomposition，简称SVD)——在线性代数中，SVD是重要的实数或复数矩阵的分解方法，在信号处理和统计中有多种应用，比如计算矩阵的伪逆矩阵(以求解最小二乘法问题)、解决超定线性系统(overdetermined linear systems)、矩阵逼近、数值天气预报等等。

29、求解线性方程组(Solving a system of linear equations)——线性方程组是数学中最古老的问题，它们有很多应用，比如在数字信号处理、线性规划中的估算和预测、数值分析中的非线性问题逼近等等。求解线性方程组，可以使用高斯—约当消去法(Gauss-Jordan elimination)，或是柯列斯基分解( Cholesky decomposition)。

30、Strukturtensor算法——应用于模式识别领域，为所有像素找出一种计算方法，看看该像素是否处于同质区域( homogenous region)，看看它是否属于边缘，还是是一个顶点。

31、合并查找算法(Union-find)——给定一组元素，该算法常常用来把这些元素分为多个分离的、彼此不重合的组。不相交集(disjoint-set)的数据结构可以跟踪这样的切分方法。合并查找算法可以在此种数据结构上完成两个有用的操作：

查找：判断某特定元素属于哪个组。

合并：联合或合并两个组为一个组。

32、维特比算法(Viterbi algorithm)——寻找隐藏状态最有可能序列的动态规划算法，这种序列被称为维特比路径，其结果是一系列可以观察到的事件，特别是在隐藏的Markov模型中。

以上就是Christoph博士对于最重要的算法的调查结果。你们熟悉哪些算法?又有哪些算法是你们经常使用的?

❹ 如何修改strassen算法

#define N 8
//用到的函数：
//矩阵加法：
void MatrixAdd(int n,float A[][N],float B[][N],float R[][N])
{ 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 MatrixSub(int n,float A[][N],float B[][N],float R[][N])
{ int j,i;
for(i=0;i<n;i++)
for(j=0;i<n;j++)
{
R[i][j]=A[i][j]+B[i][j];
}
}
//矩阵相乘：
//传递数值应该可以用FOR循环
void MatrixMultiply(float A[][N],float B[][N],float R[][N])
{
int i,j,t;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
R[i][j]=0;
for(t=0;t<2;t++)
R[i][j]=R[i][j]+A[i][t]*B[t][j];
}

}
//strassen算法！！！
void strassen(int n,float A[][N],float B[][N],float R[][N])
{
int i,j;
//for控制变量
//分块矩阵：
float A11[N][N];
float A12[N][N];
float A21[N][N];
float A22[N][N];
float B11[N][N];
float B12[N][N];
float B21[N][N];
float B22[N][N];
//中间矩阵：
float Mid[N][N];
float Midd[N][N];
//M块（7个）
float M1[N][N];
float M2[N][N];
float M3[N][N];
float M4[N][N];
float M5[N][N];
float M6[N][N];
float M7[N][N];
if(n==2)
{
MatrixMultiply(A ,B ,R );
}
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]=A[i][j];
B12[i][j]=A[i][j+(n/2)];
B21[i][j]=A[i+(n/2)][j];
B22[i][j]=A[i+(n/2)][j+(n/2)];
}
}
//分别求出M1---M7；
MatrixSub(n/2,B12,B22,Mid);
strassen(n/2,A11,Mid,M1);

MatrixAdd(n/2,A11,A12,Mid);
strassen(n/2,Mid,B22,M2);

MatrixAdd(n/2,A21,A22,Mid);
strassen(n/2,Mid,B11,M3);

MatrixSub(n/2, B21,B11,Mid);
strassen(n/2,A22,Mid,M4);

MatrixAdd(n/2,A11,A22,Mid);
MatrixAdd(n/2,B11,B22,Midd);
strassen(n/2,Mid,Midd,M5);

MatrixSub(n/2,A12,A22,Mid);
MatrixAdd(n/2,B21,B22,Midd);
strassen(n/2,Mid,Midd,M6);

MatrixSub(n/2,A11,A22,Mid);
MatrixAdd(n/2,B11,B12,Midd);
strassen(n/2,Mid,Midd,M7);
//组合矩阵：中间矩阵
float C11[N][N];
float C12[N][N];
float C21[N][N];
float C22[N][N];
MatrixAdd(N/2,M5,M4,Mid);
MatrixAdd(N/2,Mid,M2,Mid);
MatrixAdd(N/2,Mid,M6,C11);

MatrixAdd(N/2,M1,M2,C12);

MatrixAdd(N/2,M3,M4,C21);

MatrixAdd(N/2,M5,M1,Mid);
MatrixAdd(N/2,Mid,M3,Mid);
MatrixSub(N/2,Mid,M7,C11);

for(i=0;i<n/2;i++)
for(j=0;j<n/2;j++)
{
R[i][j]=C11[i][j];
R[i][j+n/2]=C12[i][j];
R[i+n/2][j]=C21[i][j];
R[i+n/2][j+n/2]=C22[i][j];
} //计算结果送回C[N][N]

}

int _tmain(int argc, _TCHAR* argv[])
{

float A[N][N],B[N][N];

for(int i=1;i<=N;i++)
{
for(int j=1;j<=N;j++)
{

A[i-1][j-1]=float(1)/float(i+j-1);
printf("%f ",A[i-1][j-1]);
}
}

strassen(8,A,A,B);
getchar( );
return 0;
}

❺ 怎样修改矩阵乘法的strassen算法，它也可以用于大小不必为2的幂的矩阵

strassen 之所以要2分的原因在于是矩阵2分后，两个子矩阵加减法要满足相同规模。
首先最简单的：把矩阵补全成为2的幂次规模即可。由于矩阵乘法性质，就算扩大矩阵（补0），也会保留原有的结果，而其他部分为0，也就是说算完之后再从结果矩阵将需要部分扣下来即可。
其次稍微动脑子的：规定一个最小乘法规模m*m，当子问题小于该规模的时候，用普通矩阵乘法，于是问题转化为，把当前矩阵规模扩大为 p*p，满足 n*2^r = p, 且n < m，我们的目标就成了找到一组比较小的 n , r
但是，毕竟要追求卓越，实际问题中补成2次幂的代价往往很大，比如总共1.5T的空间,矩阵们 总共占了1T，最差情况我们要补成2T，显然不合理，因此扩充矩阵的条件只考虑判断奇偶性即可，在每一次递归计算的时候如果发现矩阵规模是奇数，我们把它+1补成偶数即可。