程序算法艺术与实践:递归策略之矩阵乘法问题

矩阵预备知识

在数学中，矩阵（Matrix）是指纵横排列的二维数据表格，最早来自于方程组的系数及常数所构成的方阵。这一概念由19世纪英国数学家凯利首先提出。 矩阵是高等代数学中的常见工具，也常见于统计分析等应用数学学科中。并且在ACM竞赛，有很多涉及到矩阵知识的题。许多算法都会结合矩阵来处理，而比较具有代表性的矩阵算法有：矩阵快速幂、高斯消元等等。例如下面的图片就是一个矩阵：

上述矩阵是一个 4 × 3 矩阵：某矩阵 A 的第 i 行第 j 列，或 I , j位，通常记为 A[i,j] 或 A_ij。在上述例子中A[3,3]=37。此外 A_ij = (a_ij)，意为 A[i,j] = a_ij。①  矩阵相乘的规则：矩阵与矩阵相乘 第一个矩阵的列数必须等于第二个矩阵的行数 假如第一个是m*n的矩阵 第二个是n*p的矩阵 则结果就是m*p的矩阵 且得出来的矩阵中元素具有以下特点：第一行第一列元素为第一个矩阵的第一行的每个元素和第二个矩阵的第一列的每个元素乘积的和 以此类推 第i行第j列的元素就是第一个矩阵的第i行的每个元素与第二个矩阵第j列的每个元素的乘积的和。②  单位矩阵： n*n的矩阵  mat ( i , i )=1;  任何一个矩阵乘以单位矩阵就是它本身 n*单位矩阵=n， 可以把单位矩阵等价为整数1。（单位矩阵用在矩阵快速幂中）例如下图就是一个3*3的单位矩阵：

问题描述:

设A1,A2,…,An为矩阵序列，Ai为Pi-1×Pi阶矩阵，i = 1,2,…,n. 确定乘法顺序使得元素相乘的总次数最少.输入：向量P = <P0, P1, … , Pn>

实例： 

P = <10, 100, 5, 50>  A1: 10 × 100, A2: 100 × 5, A3: 5 × 50

乘法次序：

(A1 A2)A3: 10 × 100 × 5 + 10 ×5 × 50 = 7500

        A1(A2 A3): 10 × 100 × 50 + 100 × 5 × 50 = 75000

搜索空间的规模先将矩阵链加括号分为两部分，即P=A1*A2*...*An=（A1*A2...*Ak）*（Ak+1*...An），则有f(n)=f(1)*f(n-1)+f(2)*f(n-2)+...+f(n-1)*f(1)种方法。f(n)为一个Catalan数，所以一般的方法要计算种。

动态规划算法

输入P=< P0, P1, …, Pn>，Ai..j 表示乘积 AiAi+1…Aj 的结果，其最后一次相乘是：

m[i,j] 表示得到Ai..j的最少的相乘次数。

递推方程：

为了确定加括号的次序，设计表s[i,j],记录求得最优时最一位置。

算法递归实现

由上面的递归公式，很容易得到算法的递归实现,用一个N=10, P=<30,35,15,5,40,20,10,8,9,60,80>的简单实例.参考代码如下所示：

int RecurMatrixChain(int *ptrArray,int start,int end){      m[start][end]=100000;      s[start][end]=start;      if(start==end) m[start][end]=0;      else{          for(int k=start;k<end;k++){              int q=RecurMatrixChain(ptrArray,start,k)+RecurMatrixChain(ptrArray,k+1,end)+ptrArray[start]*ptrArray[k+1]*ptrArray[end+1];              if(q<m[start][end]){                  m[start][end]=q;                  s[start][end]=k;              }          }      }      return m[start][end];  }  

测试结果如下所示:

递归实现的复杂性

复杂性满足递推关系：

由数学归纳法可得：

可见递归实现的复杂性虽然较一般算法有改进，但还是较高。分析原因，主要是子问题重复程度高。如下图所示(ps:此图来自网络)：

1..4表示计算Ai..j中i=1,j=4的子问题，其子问题包括A1..1，而A1..2，A1..3中都包括子问题A1..1，所以很多子问题被重复计算了多次。于是，我们想到用自底向上的迭代实现。

算法迭代实现

迭代实现主要思想是子问题由小到大，每个子问题只计算一次，并且把结果保存起来，后来用到这个子问题时，直接代入。

void MatrixChain(int *P,int n){int r,i,j,k,t;for(i=0;i<N;i++)for(j=0;j<N;j++)m[i][j]=0;//r为当前计算的链长（子问题规模）for(r=2;r<=n;r++){  //n-r+1为最后一个r链的前边界for(i=0;i<n-r+1;i++){//计算前边界为r，链长为r的链的后边界j=i+r-1;//将链ij划分为A(i) * ( (A(i+1) ... A(j) )m[i][j]=m[i+1][j]+P[i]*P[i+1]*P[j+1];//记录分割位置s[i][j]=i;for( k=i+1;k<j-1;k++){//将链ij划分为( A(i)...A(k) )* ( (A(k+1) ... A(j) )t=m[i][k]+m[k+1][j]+P[i]*P[i+1]*P[j+1];if(t<m[i][j]){m[i][j]=t;s[i][j]=k;}}}}}

其中: 迭代实现的复杂性.行7,9,16的循环为O(n)，外层循环为O(1)，所以算法复杂度W(n)=O(n^3).迭代过程的一个实例:子问题由小到大的计算过程如下图所示(ps:此图来自网上)：

再写一个打印结果，以及打印优化函数备忘录m和标记函数的s的函数,：

void PrintMatrixChain(int s[][N],int i,int j){if (i==j) { cout<<"A"<<i+1; }  else { cout<<"("; PrintMatrixChain(s, i, s[i][j]); PrintMatrixChain(s, s[i][j]+1, j); cout<<")"; } }void PrintMS(int m[][N],int s[][N],int N){for(int r=0;r<N;r++){for(int i=0;i<N-r;i++){int j=i+r;cout<<"m["<<i+1<<","<<j+1<<"]="<<m[i][j]<<"\t";}cout<<endl;}for(int r=1;r<5;r++){for(int i=0;i<N-r;i++){int j=i+r;cout<<"s["<<i+1<<","<<j+1<<"]="<<s[i][j]+1<<"\t";}cout<<endl;}}

一个简单的测试实例:用一个N=6，P=<{8,30,80,35,5,10,20}>的简单实例，运行上述代码：

HDOJ 4920 Matrix multiplication

Problem Description

Given two matrices A and B of size n×n, find the product of them.bobo hates big integers. So you are only asked to find the result modulo 3.

Input

The input consists of several tests. For each tests:

The first line contains n (1≤n≤800). Each of the following n lines contain n integers — the description of the matrix A. The j-th integer in the i-th line equals Aij. The next n lines describe the matrix B in similar format (0≤Aij,Bij≤109).

Output

For each tests: Print n lines. Each of them contain n integers — the matrix A×B in similar format.

Sample Input

1 0 1 2 0 1 2 3 4 5 6 7

Sample Output

0 0 1 2 1

题目大意：求两个矩阵相乘mod3的结果。这道题就是一个赤裸裸的矩阵相乘的题。但是要注意题目的时间复杂度，一定要进行优化。并且，此题还有二个坑点，如果把握不好就会超时。① 就是Mod 3 时，一定不能在矩阵相乘算法的for循环里，如果放在相乘算法里会超时。②在进行相乘之前把所有元素先Mod 3 这样做也会优化一定的时间。因为题目所给数据并不是很大，所以即使把mod 3 放到结束输出语句里面就可以了，不用担心相乘的过程会超出int型。

源代码参考：

   int X[MAXN][MAXN];    int Y[MAXN][MAXN];    int t[MAXN][MAXN];     int main() {        while(scanf("%d",&N)!=EOF){            for(int i=0;i<N;i++){                 for(int j=0;j<N;j++){                    scanf("%d",&X[i][j]);                     X[i][j]%=3;                 }            }            for(int i=0;i<N;i++){                 for(int j=0;j<N;j++){                    scanf("%d",&Y[i][j]);                     Y[i][j]%=3;                     t[i][j]=0;                 }            }            for(int i=0;i<N;i++){                 for(int k=0;k<N;k++)                     if(X[i][k])                         for(intj=0;j<N;j++)                            t[i][j]=(t[i][j]+X[i][k]*Y[k][j]%3)%3;            }            for(int i=0;i<N;i++){                 printf("%d",t[i][0]);                 for(int j=1;j<N;j++){                     printf("%d",t[i][j]);                 }                 printf("\n");            }        }        return 0;  }

单独考察矩阵相乘的题目并不多见，普遍都是结合矩阵快速幂，高斯消元等矩阵算法。进行这些算法之前我们必须要掌握矩阵相乘算法.我们需掌握简单的相乘规则，才能学习后面的一些矩阵算法。同时，为了以后学习算法以及做题的需要，我们也要记得矩阵相乘时的算法复杂度及优化细节。

关于程序算法艺术与实践更多讨论与交流，敬请关注本博客和新浪微博songzi_tea.