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補成偶數即可。