矩阵的乘法算法

一般矩阵乘法算法:

原理:矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的栏数(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) + cn  (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[][N], int B[][N], int 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 sub(int n,int A[][N],int B[][N],int 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];
}

两方阵相乘问题:

//阶为2的方阵相乘:
void multiply(int A[][N],int B[][N],int R[][N])
{
    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[][N], int n, int X11[][N], int X12[][N], int X21[][N], int X22[][N] )
{
    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[][N], int B[][N], int C[][N]) {
     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[][N], int B[][N], int 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 sub(int n, int A[][N], int B[][N], int 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 strassen(int n, int A[][N], int B[][N], int C[][N])
{
    int i, j;
    int A11[N][N], A12[N][N], A21[N][N], A22[N][N];
    int B11[N][N], B12[N][N], B21[N][N], B22[N][N];
    int C11[N][N], C12[N][N], C21[N][N], C22[N][N];
    int AA[N][N], BB[N][N];
    int M1[N][N], M2[N][N], M3[N][N], M4[N][N], M5[N][N], M6[N][N], M7[N][N];
    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[N][N], B[N][N], C[N][N];
    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很大的时候才有价值,因此,可以说这方面的研究更偏重于理论价值。

时间: 2024-12-31 02:51:21

矩阵的乘法算法的相关文章

【算法导论】动态规划之矩阵链乘法

       所谓矩阵链乘法是指当一些矩阵相乘时,如何加括号来改变乘法顺序从而来降低乘法次数.例如有三个矩阵连乘:A1*A2*A3,其维数分别为:10*100,100*5,5*50.如果按照((A1*A2)*A3)来计算的话,求(A1*A2)要10*100*5=5000次乘法,再乘以A3需要10*5*50=2500次乘法,因此总共需要7500次乘法.如果按照(A1*(A2*A3))来计算的话,求(A2*A3)要100*5*50=25000次乘法,再乘以A1需要10*100*50=50000次乘法

第十五章 动态规划——矩阵链乘法

前言:今天接着学习动态规划算法,学习如何用动态规划来分析解决矩阵链乘问题.首先回顾一下矩阵乘法运算法,并给出C++语言实现过程.然后采用动态规划算法分析矩阵链乘问题并给出C语言实现过程. 1.矩阵乘法       从定义可以看出:只有当矩阵A的列数与矩阵B的行数相等时A×B才有意义.一个m×r的矩阵A左乘一个r×n的矩阵B,会得到一个m×n的矩阵C.在计算机中,一个矩阵说穿了就是一个二维数组.一个m行r列的矩阵可以乘以一个r行n列的矩阵,得到的结果是一个m行n列的矩阵,其中的第i行第j列位置上的

Ruby实现的矩阵连乘算法

  这篇文章主要介绍了Ruby实现的矩阵连乘算法,本文直接给出实现代码,需要的朋友可以参考下 动态规划解决矩阵连乘问题,随机产生矩阵序列,输出形如((A1(A2A3))(A4A5))的结果. 代码: ? 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 5

在小波变换中,两矩阵求差算法中,+128是什么意思?

问题描述 在小波变换中,两矩阵求差算法中,+128是什么意思? /********************************************************************** 函数名称: MatrixSub() 参数: unsigned char *matrixA 待求差矩阵A数组指针 unsigned char *matrixB 待求差矩阵B数组指针 int height 高度 int width 宽度 unsigned char *result 差矩阵数组指针

矩阵乘法-c++求两个矩阵的乘法,但运行出来的结果是一些不可预料的值

问题描述 c++求两个矩阵的乘法,但运行出来的结果是一些不可预料的值 #include #include using namespace std; class CreateMatrix { private: int m,n,**p; public: void SetMatrix(int m,int n) { p=new int *[m]; //build m pointers and save m rows for(int i=0;i if((p[i]=new int [n])==NULL) e

c语言-关于矩阵快速转置算法问题

问题描述 关于矩阵快速转置算法问题 请问圈出的部分是什么意思啊 解决方案 图片看不清,应该和这个类似,而且有详细注释:http://programdesign.baike.com/article-379765.html 解决方案二: 7. 矩阵的快速转置算法矩阵的快速转置矩阵的快速转置

用Spark学习矩阵分解推荐算法

在矩阵分解在协同过滤推荐算法中的应用中,我们对矩阵分解在推荐算法中的应用原理做了总结,这里我们就从实践的角度来用Spark学习矩阵分解推荐算法. 1. Spark推荐算法概述 在Spark MLlib中,推荐算法这块只实现了基于矩阵分解的协同过滤推荐算法.而基于的算法是FunkSVD算法,即将m个用户和n个物品对应的评分矩阵M分解为两个低维的矩阵: 其中k为分解成低维的维数,一般远比m和n小.如果大家对FunkSVD算法不熟悉,可以复习对应的原理篇. 2. Spark推荐算法类库介绍 在Spar

C++实现大数乘法算法代码_C 语言

C++实现大数乘法算法代码 复制代码 代码如下: //大数乘法算法 #include<iostream> #include<string> #include<cstring> using namespace std; int main() {     string num1,num2;     cin >> num1 >> num2;     //cout << num1.size() << " " &

《数学与泛型编程:高效编程的奥秘》一2.1 埃及乘法算法

2.1 埃及乘法算法 与所有的古文明一样,埃及的计数系统也没有按位置计数这一概念,而且无法表示0.于是,乘法计算起来就特别困难,只有少数受过训练的专家才会做.(你可以想象一下:如果自己只能像使用罗马数字那样来做运算,而且要计算的数字又很大,那么算起来是相当困难的.)怎样来定义乘法呢?宽泛地说,我们可以认为乘法就是"把某物多次加到它自己上面",如果说得严谨一些,那么可以分两种情况来定义:一种情况是乘以1,另一种情况是乘以一个大于1的数.我们将乘以1的乘法运算,定义如下:1a = a(2.