【算法导论】矩阵乘法

离过年都不到十天了，还要等到这周五才能回家，想想也一年没回家了。从寒假开始到现在，已经有二十来天，这期间把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

时间： 2024-08-28 01:30:43

【算法导论】矩阵乘法的相关文章

算法导论-矩阵乘法-strassen算法

目录 1.矩阵相乘的朴素算法 2.矩阵相乘的strassen算法 3.完整测试代码c++ 4.性能分析 5.参考资料 内容 1.矩阵相乘的朴素算法 T(n) = Θ(n3) 朴素矩阵相乘算法,思想明了,编程实现简单.时间复杂度是Θ(n^3).伪码如下 1 for i ← 1 to n 2 do for j ← 1 to n 3 do c[i][j] ← 0 4 for k ← 1 to n 5 do c[i][j] ← c[i][j] + a[i][k]⋅ b[k][j] 2.矩阵相乘的stra

蓝桥杯- 算法训练 矩阵乘法

算法训练 矩阵乘法 时间限制:1.0s   内存限制:512.0MB 问题描述 输入两个矩阵,分别是m*s,s*n大小.输出两个矩阵相乘的结果. 输入格式 第一行,空格隔开的三个正整数m,s,n(均不超过200). 接下来m行,每行s个空格隔开的整数,表示矩阵A(i,j). 接下来s行,每行n个空格隔开的整数,表示矩阵B(i,j). 输出格式 m行,每行n个空格隔开的整数,输出相乘後的矩阵C(i,j)的值. 样例输入 2 3 21 0 -11 1 -30 31 23 1 样例输出 -3 2-8

模板C++ 02数论算法 4矩阵乘法

矩阵乘法:用来求某种 递推关系. 矩阵相乘只有在第一个矩阵的列数和第二个矩阵的行数相同时才有意义. 定义 设A为A*M的矩阵,B为M*B的矩阵,那么矩阵C为矩阵A与B的乘积,其中矩阵C中的第i行第j列元素可以表示为: 如下所示: 开一个2*2的矩阵:主要是为了快速幂的方便,一个可以和自己乘上许多次(>=2)的矩阵只有可能是正方形的,所以要开这样一个矩阵. [题目描述] a[1]=a[2]=a[3]=1 a[x]=a[x-3]+a[x-1]  (x>3) 求a数列的第n项对1000000007(

算法提高 矩阵乘法 区间dp

问题描述 有n个矩阵,大小分别为a0*a1, a1*a2, a2*a3, ..., a[n-1]*a[n],现要将它们依次相乘,只能使用结合率,求最少需要多少次运算. 两个大小分别为p*q和q*r的矩阵相乘时的运算次数计为p*q*r. 输入格式 输入的第一行包含一个整数n,表示矩阵的个数. 第二行包含n+1个数,表示给定的矩阵. 输出格式 输出一个整数,表示最少的运算次数. 样例输入 3 1 10 5 20 样例输出 150 数据规模和约定 1<=n<=1000, 1<=ai<=1

算法训练 矩阵乘法

问题描述 输入两个矩阵,分别是m*s,s*n大小.输出两个矩阵相乘的结果. 输入格式 第一行,空格隔开的三个正整数m,s,n(均不超过200). 接下来m行,每行s个空格隔开的整数,表示矩阵A(i,j). 接下来s行,每行n个空格隔开的整数,表示矩阵B(i,j). 输出格式 m行,每行n个空格隔开的整数,输出相乘後的矩阵C(i,j)的值. 样例输入 2 3 2 1 0 -1 1 1 -3 0 3 1 2 3 1 样例输出 -3 2 -8 2 提示 矩阵C应该是m行n列,其中C(i,j)等于矩阵A

算法导论--矩阵链相乘

#include<iostream> using namespace std; /* 计算括号化方案数:标量乘法作为代价衡量,应该使标量乘法尽可能少. m[i,j]表示Ai.....Aj所需标量乘法的最小值. i=j 时只有一个矩阵,无需分割 m[i,i]=0; 采用自底向上的方式: */ int m[100][100]; int p[]={30,35,15,5,10,20,25}; int bottomcut(int n){ int t; for(int l=2;l<=n;l++){

第四章 分治策略 4.2 矩阵乘法的Strassen算法

package chap04_Divide_And_Conquer; import static org.junit.Assert.*; import java.util.Arrays; import org.junit.Test; /** * 矩阵相乘的算法 * * @author xiaojintao * */ public class MatrixOperation { /** * 普通的矩阵相乘算法,c=a*b.其中,a.b都是n*n的方阵 * * @param a * @param b

Spark中的矩阵乘法分析

前言: 矩阵乘法在数据挖掘/机器学习中是常用的计算步骤,并且在大数据计算中,shuffle过程是不可避免的,矩阵乘法的不同计算方式shuffle的数据量都不相同.通过对矩阵乘法不同计算方式的深入学习,希望能够对大数据算法实现的shuffle过程优化有所启发.网上有很多分布式矩阵乘法相关的文章和论文,但是鲜有对Spark中分布式矩阵乘法的分析.本文针对Spark中分布式矩阵乘法的实现进行必要的说明讨论. 分布式矩阵乘法原理: 矩阵乘法计算可以分为内积法和外积法.根据实现颗粒度的不同,也可以分为普通

算法导论--动态规划（矩阵链乘法）

矩阵链乘法问题 给定一个n个矩阵的序列?A1,A2,A3...An?,我们要计算他们的乘积:A1A2A3...An.因为矩阵乘法满足结合律,加括号不会影响结果.可是不同的加括号方法.算法复杂度有非常大的区别: 考虑矩阵链:?A1,A2,A3?.三个矩阵规模分别为10×100.100×5.5×50 假设按((A1A2)A3)方式,须要做10?100?5=5000次,再与A3相乘,又须要10?5?50=2500,共须要7500次运算: 假设按(A1(A2A3))方式计算.共须要100?5?50+10

《算法导论》读书笔记之动态规划—矩阵链乘法

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