算法导论第15章 动态规划-矩阵链乘法

前言：今天接着学习动态规划算法，学习如何用动态规划来分析解决矩阵链乘问题。首先回顾一下矩阵乘法运算法，并给出C++语言实现过程。然后采用动态规划算法分析矩阵链乘问题并给出C语言实现过程。

1、矩阵乘法

　

　　

 

 

 

　　

　　从定义可以看出：只有当矩阵A的列数与矩阵B的行数相等时A×B才有意义。一个m×r的矩阵A左乘一个r×n的矩阵B，会得到一个m×n的矩阵C。在计算机中，一个矩阵说穿了就是一个二维数组。一个m行r列的矩阵可以乘以一个r行n列的矩阵，得到的结果是一个m行n列的矩阵，其中的第i行第j列位置上的数等于前一个矩阵第i行上的r个数与后一个矩阵第j列上的r个数对应相乘后所有r个乘积的和。采用C++语言实现完整的两个矩阵乘法，程序如下所示：

 1 #include <iostream> 2 using namespace std; 3 #define A_ROWS        3 4 #define A_COLUMNS     2 5 #define B_ROWS        2 6 #define B_COLUMNS     3 7 void matrix_multiply(int A[A_ROWS][A_COLUMNS],int B[B_ROWS][B_COLUMNS],int C[A_ROWS][B_COLUMNS]); 8 int main() 9 {10     int A[A_ROWS][A_COLUMNS] = {1,0,11                                 1,2,12                                 1,1};13     int B[B_ROWS][B_COLUMNS] = {1,1,2,14                                 2,1,2};15     int C[A_ROWS][B_COLUMNS] = {0};16     matrix_multiply(A,B,C);17     for(int i=0;i<A_ROWS;i++)18     {19         for(int j=0;j<B_COLUMNS;j++)20             cout<<C[i][j]<<" ";21         cout<<endl;22     }23     return 0;24 }25 void matrix_multiply(int A[A_ROWS][A_COLUMNS],int B[B_ROWS][B_COLUMNS],int C[A_ROWS][B_COLUMNS])26 {27     if(A_COLUMNS != B_ROWS)28         cout<<"error: incompatible dimensions."<<endl;29     else30     {31         int i,j,k;32         for(i=0;i<A_ROWS;i++)33             for(j=0;j<B_COLUMNS;j++)34             {35                 C[i][j] = 0;36                 for(k=0;k<A_COLUMNS;k++)37                     C[i][j] += A[i][k] * B[k][j]; //将A的每一行的每一列与B的每一列的每一行的乘积求和38             }39     }40 }

程序测试结果如下所示：

2、矩阵链乘问题描述

　　给定n个矩阵构成的一个链<A₁,A₂,A₃,.......A_n>，其中i=1,2,...n，矩阵A的维数为p_i-1p_i，对乘积 A₁A₂...A_n 以一种最小化标量乘法次数的方式进行加全部括号。

　　注意：在矩阵链乘问题中，实际上并没有把矩阵相乘，目的是确定一个具有最小代价的矩阵相乘顺序。找出这样一个结合顺序使得相乘的代价最低。

3、动态规划分析过程

1）最优加全部括号的结构

　　动态规划第一步是寻找一个最优的子结构。假设现在要计算A_iA_i+1....A_j的值，计算A_i...j过程当中肯定会存在某个k值（i<=k<j）将A_i...j分成两部分，使得A_i...j的计算量最小。分成两个子问题A_i...k和A_k+1...j，需要继续递归寻找这两个子问题的最优解。

　　有分析可以到最优子结构为：假设A_iA_i+1....A_j的一个最优加全括号把乘积在A_k和A_k+1之间分开，则A_i..k和A_k+1..j也都是最优加全括号的。

2）一个递归解

　　设m[i,j]为计算机矩阵A_i...j所需的标量乘法运算次数的最小值，对此计算A_1..n的最小代价就是m[1,n]。现在需要来递归定义m[i,j]，分两种情况进行讨论如下：

当i==j时：m[i,j] = 0，（此时只包含一个矩阵）

当i<j 时：从步骤1中需要寻找一个k（i≤k＜j）值，使得m[i,j] ＝min{m[i,k]+m[k+1,j]+p_i-1p_kp_j} （i≤k＜j）。

3）计算最优代价

　　虽然给出了递归解的过程，但是在实现的时候不采用递归实现，而是借助辅助空间，使用自底向上的表格进行实现。设矩阵Ai的维数为pi-1pi，i=1,2.....n。输入序列为：p=<p0，p1,...pn>，length[p] = n+1。使用m[n][n]保存m[i,j]的代价，s[n][n]保存计算m[i,j]时取得最优代价处k的值，最后可以用s中的记录构造一个最优解。书中给出了计算过程的伪代码，摘录如下：

 1 MAXTRIX_CHAIN_ORDER(p) 2   n = length[p]-1; 3   for i=1 to n 4       do m[i][i] = 0; 5   for t = 2 to n  //t is the chain length 6        do for i=1 to n-t+1 7                      j=i+t-1; 8                      m[i][j] = MAXLIMIT; 9                      for k=i to j-110                             q = m[i][k] + m[k+1][i] + qi-1qkqj;11                             if q < m[i][j]12                                then m[i][j] = q;13                                     s[i][j] = k;14   return m and s;

MATRIX_CHAIN_ORDER具有循环嵌套，深度为3层，运行时间为O(n³)。如果采用递归进行实现，则需要指数级时间Ω(2ⁿ)，因为中间有些重复计算。递归是完全按照第二步得到的递归公式进行计算，递归实现如下所示：

 1 int recursive_matrix_chain(int *p,int i,int j,int m[N+1][N+1],int s[N+1][N+1]) 2 { 3     if(i==j) 4        m[i][j] = 0; 5     else 6     { 7         int k; 8         m[i][j] = MAXVALUE; 9         for(k=i;k<j;k++)10         {11             int temp = recursive_matrix_chain(p,i,k,m,s) +recursive_matrix_chain(p,k+1,j,m,s) + p[i-1]*p[k]*p[j];12             if(temp < m[i][j])13             {14                 m[i][j] = temp;15                 s[i][j] = k;16             }17         }18     }19     return m[i][j];20 }

 对递归算计的改进，可以引入备忘录，采用自顶向下的策略，维护一个记录了子问题的表，控制结构像递归算法。完整程序如下所示：

 1 int memoized_matrix_chain(int *p,int m[N+1][N+1],int s[N+1][N+1]) 2 { 3     int i,j; 4     for(i=1;i<=N;++i) 5         for(j=1;j<=N;++j) 6         { 7            m[i][j] = MAXVALUE; 8         } 9     return lookup_chain(p,1,N,m,s);10 }11 12 int lookup_chain(int *p,int i,int j,int m[N+1][N+1],int s[N+1][N+1])13 {14     if(m[i][j] < MAXVALUE)15         return m[i][j]; //直接返回，相当于查表16     if(i == j)17         m[i][j] = 0;18     else19     {20         int k;21         for(k=i;k<j;++k)22         {23             int temp = lookup_chain(p,i,k,m,s)+lookup_chain(p,k+1,j,m,s) + p[i-1]*p[k]*p[j];  //通过递归的形式计算，只计算一次，第二次查表得到24             if(temp < m[i][j])25             {26                 m[i][j] = temp;27                 s[i][j] = k;28             }29         }30     }31     return m[i][j];32 }

4)构造一个最优解

第三步中已经计算出来最小代价，并保存了相关的记录信息。因此只需对s表格进行递归调用展开既可以得到一个最优解。书中给出了伪代码，摘录如下：

1 PRINT_OPTIMAL_PARENS(s,i,j)2   if i== j 3      then print "Ai"4   else5      print "(";6      PRINT_OPTIMAL_PARENS(s,i,s[i][j]);7      PRINT_OPTIMAL_PARENS(s,s[i][j]+1,j);8      print")";

4、编程实现

　　采用C++语言实现这个过程，现有矩阵A₁(30×35)、A₂(35×15)_、A3(15×5)、A4(5×10)、A5(10×20)、A6(20×25)，得到p=<30,35,15,5,10,20,25>。实现过程定义两个二维数组m和s，为了方便计算其第一行和第一列都忽略，行标和列标都是1开始。完整的程序如下所示：

 1 #include <iostream> 2 using namespace std; 3  4 #define N 6 5 #define MAXVALUE 1000000 6  7 void matrix_chain_order(int *p,int len,int m[N+1][N+1],int s[N+1][N+1]); 8 void print_optimal_parents(int s[N+1][N+1],int i,int j); 9 10 int main()11 {12     int p[N+1] = {30,35,15,5,10,20,25};13     int m[N+1][N+1]={0};14     int s[N+1][N+1]={0};15     int i,j;16     matrix_chain_order(p,N+1,m,s);17     cout<<"m value is: "<<endl;18     for(i=1;i<=N;++i)19     {20         for(j=1;j<=N;++j)21             cout<<m[i][j]<<" ";22         cout<<endl;23     }24     cout<<"s value is: "<<endl;25     for(i=1;i<=N;++i)26     {27         for(j=1;j<=N;++j)28             cout<<s[i][j]<<" ";29         cout<<endl;30     }31     cout<<"The result is:"<<endl;32     print_optimal_parents(s,1,N);33     return 0;34 }35 36 void matrix_chain_order(int *p,int len,int m[N+1][N+1],int s[N+1][N+1])37 {38     int i,j,k,t;39     for(i=0;i<=N;++i)40         m[i][i] = 0;41     for(t=2;t<=N;t++)  //当前链乘矩阵的长度42     {43         for(i=1;i<=N-t+1;i++)  //从第一矩阵开始算起，计算长度为t的最少代价44         {45             j=i+t-1;//长度为t时候的最后一个元素46             m[i][j] = MAXVALUE;  //初始化为最大代价47             for(k=i;k<=j-1;k++)   //寻找最优的k值，使得分成两部分k在i与j-1之间48             {49                 int temp = m[i][k]+m[k+1][j] + p[i-1]*p[k]*p[j];50                 if(temp < m[i][j])51                 {52                     m[i][j] = temp;   //记录下当前的最小代价53                     s[i][j] = k;      //记录当前的括号位置，即矩阵的编号54                 }55             }56         }57     }58 }59 60 //s中存放着括号当前的位置61 void print_optimal_parents(int s[N+1][N+1],int i,int j)62 {63     if( i == j)64         cout<<"A"<<i;65     else66     {67         cout<<"(";68         print_optimal_parents(s,i,s[i][j]);69         print_optimal_parents(s,s[i][j]+1,j);70         cout<<")";71     }72 73 }

程序测试结果如下所示：

5、总结

　　动态规划解决问题关键是分析过程，难度在于如何发现其子问题的结构及子问题的递归解。这个需要多多思考，不是短时间内能明白。在实现过程中遇到问题就是数组，数组的下标问题是个比较麻烦的事情，如何能够过合理的去处理，需要一定的技巧。