AI 助理
备案控制台
开发者社区 人工智能 文章 正文

【算法导论】矩阵乘法

2015-02-11 1102
版权
版权声明：
本文内容由阿里云实名注册用户自发贡献，版权归原作者所有，阿里云开发者社区不拥有其著作权，亦不承担相应法律责任。具体规则请查看《 阿里云开发者社区用户服务协议》和 《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容，填写 侵权投诉表单进行举报，一经查实，本社区将立刻删除涉嫌侵权内容。
简介： 离过年都不到十天了，还要等到这周五才能回家，想想也一年没回家了。从寒假开始到现在，已经有二十来天，这期间把2014年总结中的寒假计划也大多数完成了：The Element Of Style的阅读，三门数学课《随机过程》、《工程优化》、《数值分析》的算法实现。

离过年都不到十天了，还要等到这周五才能回家，想想也一年没回家了。从寒假开始到现在，已经有二十来天，这期间把2014年总结中的寒假计划也大多数完成了：The Element Of Style的阅读，三门数学课《随机过程》、《工程优化》、《数值分析》的算法实现。回家过年期间肯定不会写博客了，今天一看，这个月只写了三篇，于是乎今天必须再写一篇来完成这个月的基本工作量。言归正传，这篇文章写写选修课《算法设计》作业题中的矩阵乘法的三种方法。

矩阵乘法

  • 传统方法
    • 理论公式
      C=AB
      cij=∑k=1naikbkj
    • 算法实现 
void TraditionalMethod(float A[][N],float B[][N],float C[][N])//传统方法，三重循环
{
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;//之所以每次调用都清零，是因为前面是循环调用，如果只调用一次就不需要
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}
  • 分块相乘法
    • 理论公式
      A=[A11A21A12A22],B=[B11B21B12B22]
      C=AB=[C11C21C12C22]
      C11=A11B11+A12B21C12=A11B12+A12B22C21=A21B11+A22B21C22=A21B12+A22B22
    • 算法实现 
void BlockMatrix()//分块矩阵计算
{     
      for(int i=0;i<N/2;i++)              
         for(int 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];

                C11[i][j]=0;
                C12[i][j]=0;
                C21[i][j]=0;
                C22[i][j]=0;
            }       //将矩阵A和B式分为四块

         MATRIX_Multiply(N/2,A11,B11, AA);
         MATRIX_Multiply(N/2,A12,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A11,B12, AA);
         MATRIX_Multiply(N/2,A12,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B11, AA);
         MATRIX_Multiply(N/2,A22,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B12, AA);
         MATRIX_Multiply(N/2,A22,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z

    for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
        for(int 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];
        }                                            //计算结果送回C[N][N]

    }
  • Strassen法
    • 理论公式
      M1=A11(B12−B22)M2=B22(A11+A12)M3=B11(A21+A22)M4=A22(B21−B11)M5=(A11+A22)(B11+B22)M6=(A12−A22)(B21+B22)M7=(A12−A21)(B11+B12)C11=M4+M5+M6−M2C12=M1+M2C21=M3+M4C22=M1+M5−M3−M7
    • 算法实现
void STRASSEN()  //STRASSEN函数
{
    int i,j;//,x;

    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];
        }       //将矩阵A和B式分为四块




    MATRIX_SUB(N/2,B12,B22,BB);         
    MATRIX_Multiply(N/2,A11,BB,M1);

    MATRIX_ADD(N/2,A11,A12,AA);
    MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22

    MATRIX_ADD(N/2,A21,A22,AA);
    MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11

    MATRIX_SUB(N/2,B21,B11,BB);
    MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)

    MATRIX_ADD(N/2,A11,A22,AA);
    MATRIX_ADD(N/2,B11,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)


    MATRIX_SUB(N/2,A12,A22,AA);
    MATRIX_ADD(N/2,B21,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)

    MATRIX_SUB(N/2,A11,A21,AA);
    MATRIX_ADD(N/2,B11,B12,BB);
    MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
    //计算M1,M2,M3,M4,M5,M6,M7（递归部分）


    MATRIX_ADD(N/2,M5,M4,MM1);                
    MATRIX_SUB(N/2,M2,M6,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6

    MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2

    MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4

    MATRIX_ADD(N/2,M5,M1,MM1);
    MATRIX_ADD(N/2,M3,M7,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7

    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];
        }                                            //计算结果送回C[N][N]



}
  • 完整程序实现
#include <iostream>
#include<ctime>

using namespace std;

const int N=32; //常量N用来定义方阵的大小
void output(int n,float C[][N]); //函数声明部分
void TraditionalMethod(float A[][N],float B[][N],float C[][N]);//传统的矩阵相乘
void BlockMatrix();
void STRASSEN();
void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2]);

float A[N][N];
float B[N][N];
float C[N][N];  //定义三个矩阵A,B,C

float A11[N/2][N/2],A12[N/2][N/2],A21[N/2][N/2],A22[N/2][N/2];
float B11[N/2][N/2],B12[N/2][N/2],B21[N/2][N/2],B22[N/2][N/2];
float C11[N/2][N/2],C12[N/2][N/2],C21[N/2][N/2],C22[N/2][N/2];
float M1[N/2][N/2],M2[N/2][N/2],M3[N/2][N/2],M4[N/2][N/2],M5[N/2][N/2],M6[N/2][N/2],M7[N/2][N/2];
float AA[N/2][N/2],BB[N/2][N/2],MM1[N/2][N/2],MM2[N/2][N/2];
void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]); //矩阵加法函数X+Y—>Z

void main()
{
  //初始化，使相乘的两个矩阵都为全1矩阵
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        {
            A[i][j]=1;
            B[i][j]=1;
        }
    //将结果矩阵C初始化为全0矩阵
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;

    clock_t start, finish; //用于计时
    double   duration; 
    int loop=0;

    cout<<"当矩阵为"<<N<<"×"<<N<<"，循环次数为10000时："<<endl<<endl;

    //--------使用传统方法--------------//
    start = clock();  
    while(loop<10000)//循环10000次,这里可以更改
    {
        loop++;
        TraditionalMethod(A,B,C);   //传统方法计算

    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用传统方法："<<endl;
    cout<<"所需时间为："<<duration<<endl<<endl;

  //  output(N,C);  //输出计算结果

    //---------使用分块矩阵乘法------------//
    start = clock();  
    loop=0;
    while(loop<10000)//循环10000次,这里可以更改
    {
        loop++;
        BlockMatrix();   //分块矩阵计算

    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用分块相乘方法："<<endl;
    cout<<"所需时间为："<<duration<<endl<<endl;
    //  output(N,C);  //输出计算结果

    //-------使用strassen方法-------------------//
    start = clock();  
    loop=0;
    while(loop<10000)//当时间非常小时，需要加大循环次数,这里可以更改
    {
        loop++;
        STRASSEN();   //调用STRASSEN函数计算      
    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用strassen方法："<<endl;
    cout<<"所需时间为："<<duration<<endl;
//  output(N,C);  //输出计算结果
}

void TraditionalMethod(float A[][N],float B[][N],float C[][N])//传统方法，三重循环
{
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;//之所以每次调用都清零，是因为前面是循环调用，如果只调用一次就不需要
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}


void BlockMatrix()//分块矩阵计算
{     
      for(int i=0;i<N/2;i++)              
         for(int 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];

                C11[i][j]=0;
                C12[i][j]=0;
                C21[i][j]=0;
                C22[i][j]=0;
            }       //将矩阵A和B式分为四块

         MATRIX_Multiply(N/2,A11,B11, AA);
         MATRIX_Multiply(N/2,A12,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A11,B12, AA);
         MATRIX_Multiply(N/2,A12,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B11, AA);
         MATRIX_Multiply(N/2,A22,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B12, AA);
         MATRIX_Multiply(N/2,A22,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z

    for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
        for(int 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];
        }                                            //计算结果送回C[N][N]

    }



void output(int n,float C[][N]) //矩阵输出函数
{
    int i,j;
    cout<<"输出矩阵:"<<endl;
    for(i=0;i<n;i++)
    {
        cout<<endl;
        for(j=0;j<n;j++)
            cout<<C[i][j]<<"  ";
    }
    cout<<endl<<endl;

}


void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2])//矩阵加法函数X*Y—>C
{


    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            C[i][j]=0;
            for(int k=0;k<n;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}



void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵加法函数X+Y—>Z
{
    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            Z[i][j]=X[i][j]+Y[i][j];
}

void MATRIX_SUB(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵减法函数X-Y—>Z
{

    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            Z[i][j]=X[i][j]-Y[i][j];

}

void STRASSEN()  //STRASSEN函数
{


    int i,j;//,x;

    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];
        }       //将矩阵A和B式分为四块




    MATRIX_SUB(N/2,B12,B22,BB);         
    MATRIX_Multiply(N/2,A11,BB,M1);

    MATRIX_ADD(N/2,A11,A12,AA);
    MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22

    MATRIX_ADD(N/2,A21,A22,AA);
    MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11

    MATRIX_SUB(N/2,B21,B11,BB);
    MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)

    MATRIX_ADD(N/2,A11,A22,AA);
    MATRIX_ADD(N/2,B11,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)


    MATRIX_SUB(N/2,A12,A22,AA);
    MATRIX_ADD(N/2,B21,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)

    MATRIX_SUB(N/2,A11,A21,AA);
    MATRIX_ADD(N/2,B11,B12,BB);
    MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
    //计算M1,M2,M3,M4,M5,M6,M7（递归部分）


    MATRIX_ADD(N/2,M5,M4,MM1);                
    MATRIX_SUB(N/2,M2,M6,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6

    MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2

    MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4

    MATRIX_ADD(N/2,M5,M1,MM1);
    MATRIX_ADD(N/2,M3,M7,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7

    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];
        }                                            //计算结果送回C[N][N]   
}

运行结果如下图：
三种矩阵乘法的对比

原文：http://blog.csdn.net/tengweitw/article/details/43732253
作者：nineheadedbird

文章标签：
算法
C语言
机器学习/深度学习
人工智能
BI
Serverless
关键词：
算法矩阵
算法乘法
算法矩阵乘法
tengweitw
目录
相关文章
闻缺陷则喜何志丹
|
4月前
|
算法 测试技术 C++
【动态规划】【矩阵快速幂】【滚动向量】C++算法552. 学生出勤记录 II
【动态规划】【矩阵快速幂】【滚动向量】C++算法552. 学生出勤记录 II
闻缺陷则喜何志丹
42 2
红目香薰
|
4月前
|
算法 Java C语言
第十四届蓝桥杯集训——练习解题阶段(无序阶段)-ALGO-48 算法训练 关联矩阵
第十四届蓝桥杯集训——练习解题阶段(无序阶段)-ALGO-48 算法训练 关联矩阵
红目香薰
57 0
yinying293
|
1月前
|
算法 Java
LeetCode经典算法题：矩阵中省份数量经典题目+三角形最大周长java多种解法详解
LeetCode经典算法题：矩阵中省份数量经典题目+三角形最大周长java多种解法详解
yinying293
34 6
算法金
|
3月前
|
机器学习/深度学习 数据采集 人工智能
算法金 | 协方差、方差、标准差、协方差矩阵
**摘要：** 本文介绍了统计学中的基础概念，包括方差、标准差、协方差及其矩阵。方差衡量数据的分散程度，标准差是方差的平方根，提供相同单位下的波动度量。协方差则分析两个变量的关联性，正负值表示正负相关。协方差矩阵扩展到多变量情况，展示多个变量间的关系。这些工具在金融、质量控制、机器学习等领域有广泛应用。文章通过实例和公式清晰解释了每个概念，并强调理解它们之间的关系对于数据分析和统计建模的重要性。
算法金
29 0
算法金 | 协方差、方差、标准差、协方差矩阵
jjczqyjxgdegi
|
3月前
|
存储 算法 Python
技术心得记录：大整数算法【10】Comba乘法（实现）
技术心得记录：大整数算法【10】Comba乘法（实现）
jjczqyjxgdegi
21 0
阿旭算法与机器学习
|
3月前
|
算法
【经典LeetCode算法题目专栏分类】【第6期】二分查找系列：x的平方根、有效完全平方数、搜索二位矩阵、寻找旋转排序数组最小值
【经典LeetCode算法题目专栏分类】【第6期】二分查找系列：x的平方根、有效完全平方数、搜索二位矩阵、寻找旋转排序数组最小值
阿旭算法与机器学习
28 0
阿旭算法与机器学习
|
3月前
|
算法
【经典LeetCode算法题目专栏分类】【第4期】BFS广度优先算法：单词接龙、最小基因变化、二进制矩阵中的最短路径
【经典LeetCode算法题目专栏分类】【第4期】BFS广度优先算法：单词接龙、最小基因变化、二进制矩阵中的最短路径
阿旭算法与机器学习
21 0
电力程序小学童
|
4月前
|
算法
【免费】面向多微网网络结构设计的大规模二进制矩阵优化算法
【免费】面向多微网网络结构设计的大规模二进制矩阵优化算法
电力程序小学童
46 1
闻缺陷则喜何志丹
|
4月前
|
算法 测试技术 Serverless
【二分查找】【C++算法】378. 有序矩阵中第 K 小的元素
【二分查找】【C++算法】378. 有序矩阵中第 K 小的元素
闻缺陷则喜何志丹
41 1
闻缺陷则喜何志丹
|
4月前
|
算法 测试技术 C++
【字符串】【 LCP】【C++算法】2573找出对应 LCP 矩阵的字符串
【字符串】【 LCP】【C++算法】2573找出对应 LCP 矩阵的字符串
闻缺陷则喜何志丹
32 1

热门文章

最新文章

  • 1
    五分钟了解一下什么是「贪心算法 」 | 算法必看系列二十一
  • 2
    算法面试实战(二)- 面试中的复杂度分析
  • 3
    TapTap 算法平台的 Serverless 探索之路
  • 4
    《推荐系统：技术、评估及高效算法》一3.1　简介
  • 5
    某研究院的二叉树深度优先遍历变种的算法面试题以及答案
  • 6
    椭圆曲线密码算法（ECC）：数学之美与安全之钥
  • 7
    MMEditing | 新视频超分算法冠军BasicVSR++来了
  • 8
    数据结构与算法(三) 线性表之双向链表
  • 9
    机器学习算法：补一个k-近邻算法的测试
  • 10
    预约直播 | 流批一体机器学习算法平台Alink介绍及应用
  • 1
    m基于LDPC编译码的matlab误码率仿真,对比SP,MS,NMS以及OMS四种译码算法
    208
  • 2
    【AI 初识】描述遗传算法概念
    135
  • 3
    【AI 初识】人工智能中使用了哪些不同的搜索算法？
    266
  • 4
    机器学习算法原理与应用：深入探索与实战
    87
  • 5
    【C言专栏】递归算法在 C 语言中的应用
    54
  • 6
    基于DCT变换的彩色图像双重水印嵌入和提取算法matlab仿真
    72
  • 7
    排序算法--------计数排序
    32
  • 8
    Matlab决策树、模糊C-均值聚类算法分析高校教师职称学历评分可视化
    48
  • 9
    【Python机器学习专栏】异常检测算法在Python中的实践
    155
  • 10
    圆填充( CIRCLE PACKING)算法圆堆图圆形空间填充算法可视化
    77

    • 相关课程

    更多
  • 智能运维赛（复赛）：利用数据和算法，快速定位系统异常并进行根因分析
  • 智能创作赛（复赛）：相册应用中的视频故事生成算法介绍
  • 智能创作赛（初赛）：相册应用中的故事生成算法介绍
  • 相册服务中的故事生成算法介绍
  • Go语言核心编程 - 数据结构和算法
  • 神经网络概览及算法详解

    • 相关电子书

    更多
  • 数据+算法定义新世界
  • 袋鼠云基于实时计算的反黄牛算法
  • Alink：基于Apache Flink的算法平台

    • 相关实验场景

    更多
  • 推荐系统入门之使用ALS算法实现打分预测
  • 欧拉图的构造性证明与算法实现
  • RSA密码算法设计与实现
  • RSA非对称加密算法
  • 使用Swing算法实现商品推荐
    • 下一篇
    基于LNMP搭建WordPress