动态规划 - 矩阵链乘法

前言：今天接着学习动态规划算法，学习如何用动态规划来分析解决矩阵链乘问题。首先回顾一下矩阵乘法运算法，并给出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++语言实现完整的两个矩阵乘法，程序如下所示：

#include<iostream>
#include<cstdlib>
using namespace std;

#define A_ROWS        3
#define A_COLUMNS     2
#define B_ROWS        2
#define B_COLUMNS     3

void matrix_multiply(int A[A_ROWS][A_COLUMNS],int B[B_ROWS][B_COLUMNS],int C[A_ROWS][B_COLUMNS])
{
    if(A_COLUMNS!=B_ROWS)
    {
        cout<<"incompatible dimensions"<<endl;
        exit(1);
    }
    int i,j,k;
    for(i=0;i<A_ROWS;i++)
        for(j=0;j<B_COLUMNS;j++)
    {
        C[i][j]=0;
        for(k=0;k<A_COLUMNS;k++)
            C[i][j]+=A[i][k]*B[k][j];
    }
}

int main()
{
    int C[A_ROWS][B_COLUMNS];
    int A[A_ROWS][A_COLUMNS]={{1,2},{3,4},{5,6}};
    int B[B_ROWS][B_COLUMNS]={1,2,3,4,5,6};
    matrix_multiply(A,B,C);
    int i,j;
    for(i=0;i<A_ROWS;i++)
    {
        for(j=0;j<B_COLUMNS;j++)
            cout<<C[i][j]<<" ";
        cout<<endl;
    }
}

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中的记录构造一个最优解。书中给出了计算过程的伪代码，摘录如下：

MAXTRIX_CHAIN_ORDER(p)
  n = length[p]-1;
  for i=1 to n
      do m[i][i] = 0;
  for t = 2 to n  //t is the chain length
       do for i=1 to n-t+1
                     j=i+t-1;
                     m[i][j] = MAXLIMIT;
                     for k=i to j-1
                            q = m[i][k] + m[k+1][i] + qi-1qkqj;
                            if q < m[i][j]
                               then m[i][j] = q;
                                    s[i][j] = k;
  return m and s;

C++代码：

#include<iostream>
using namespace std;
#define N 6
#define MAXVALUE 100000000
void matrix_chain_order(int *p,int m[N+1][N+1],int s[N+1][N+1])
{
    int i,j,l,k;
    for(i=1;i<=N;i++)
        m[i][i]=0;
    for(l=2;l<=N;l++)
    {
        for(i=1;i<=N-l+1;i++)
        {
            j=i+l-1;
            m[i][j]=MAXVALUE;
            for(k=i;k<=j-1;k++)
            {
                int temp=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
                if(temp<m[i][j])
                {
                    m[i][j]=temp;
                    s[i][j]=k;
                }
            }
        }
    }
}

void print_optimal_parens(int s[N+1][N+1],int i,int j)
{
    if(i==j)
        cout<<"A"<<i;
    else
    {
        cout<<"(";
        print_optimal_parens(s,i,s[i][j]);
        print_optimal_parens(s,s[i][j]+1,j);
        cout<<")";
    }
}

int main()
{
    int p[N+1] = {30,35,15,5,10,20,25};
    int m[N+1][N+1]={0};
    int s[N+1][N+1]={0};
    int i,j;
    matrix_chain_order(p,m,s);
    cout<<"m value is: "<<endl;
    for(i=1;i<=N;++i)
    {
        for(j=1;j<=N;++j)
            cout<<m[i][j]<<" ";
        cout<<endl;
    }
    cout<<"s value is: "<<endl;
    for(i=1;i<=N;++i)
    {
        for(j=1;j<=N;++j)
            cout<<s[i][j]<<" ";
        cout<<endl;
    }
    cout<<"The result is:"<<endl;
    print_optimal_parens(s,1,N);
    return 0;
}