矩阵链乘法

Posted 2021-03-01 zhhhb

整理自博客：https://blog.csdn.net/c18219227162/article/details/50412333

什么是矩阵链乘法？

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

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

 

解析

两个矩阵A和B只有相容(compatible)，即A的列数等于B的行数时，才能相乘。如果A是p×q的矩阵，B是q×r的矩阵，那么乘积C是p×r的矩阵。计算C所需要时间由第8行的标量乘法的次数决定的，即pqr。

以矩阵链<A1,A2,A3>为例，来说明不同的加括号方式会导致不同的计算代价。假设三个矩阵的规模分别为10×100、100×5和5×50。

如果按照((A1A2)A3)的顺序计算，为计算A1A2(规模10×5)，需要做10*100*5=5000次标量乘法，再与A3相乘又需要做10*5*50=2500次标量乘法，共需7500次标量乘法。

如果按照(A1(A2A3))的顺序计算，为计算A2A3(规模100×50)，需100*5*50=25000次标量乘法，再与A1相乘又需10*100*50=50000次标量乘法，共需75000次标量乘法。因此第一种顺序计算要比第二种顺序计算快10倍。

矩阵链乘法问题(matrix-chain multiplication problem)可描述如下：给定n个矩阵的链<A1,A2,...,An>，矩阵Ai的规模为p(i-1)×p(i) (1<=i<=n)，求完全括号化方案，使得计算乘积A1A2...An所需标量乘法次数最少。

因为括号方案的数量与n呈指数关系，所以通过暴力搜索穷尽所有可能的括号化方案来寻找最优方案是一个糟糕策略。

 

分析

1.刻画一个最优解的结构特征

2.递归地定义最优解的值

3.计算最优解的值，通常采用自底向上的方法

4.利用计算出的信息构造一个最优解

 

1.刻画一个最优解的结构特征

动态规划的第一步是寻找最优子结构，然后就可以利用这种子结构从子问题的最优解构造出原问题的最优解。在矩阵链乘法问题中，我们假设A(i)A(i+1)...A(j)的最优括号方案的分割点在A(k)和A(k+1)之间。那么，继续对“前缀”子链A(i)A(i+1)..A(k)进行括号化时，我们应该直接采用独立求解它时所得的最优方案。

    我们已经看到，一个非平凡(i≠j)的矩阵链乘法问题实例的任何解都需要划分链，而任何最优解都是由子问题实例的最优解构成的。为了构造一个矩阵链乘法问题实例的最优解，我们可以将问题划分为两个子问题(A(i)A(i+1)...A(k)和A(k+1)A(k+2)..A(j)的最优括号化问题)，求出子问题实例的最优解，然后将子问题的最优解组合起来。我们必须保证在确定分割点时，已经考察了所有可能的划分点，这样就可以保证不会遗漏最优解。

2.递归地定义最优解的值

下面用子问题的最优解来递归地定义原问题最优解的代价。对于矩阵链乘法问题，我们可以将对所有1<=i<=j<=n确定A(i)A(i+1)...A(j)的最小代价括号化方案作为子问题。令m[i,j]表示计算矩阵A(i..j)所需标量乘法次数的最小值，那么，原问题的最优解—计算A(1..n)所需的最低代价就是m[1,n]。

   我们可以递归定义m[i,j]如下。对于i=j时的平凡问题，矩阵链只包含唯一的矩阵A(i..j)=A(i)，因此不需要做任何标量乘法运算。所以，对所有i=1,2,...,n，m[i,i]=0。若i<j，我们利用步骤1中得到的最优子结构来计算m[i,j]。我们假设A(i)A(i+1)...A(j)的最优括号化方案的分割点在矩阵A(k)和A(k+1)之间，其中i<=k<j。那么，m[i,j]就等于计算A(i..k)和A(k+1..j)的代价加上两者相乘的代价的最小值。由于矩阵Ai的大小为p(i-1)*pi，易知A(i..k)和A(k+1..j)相乘的代价为p(i-1)p(k)p(j)次标量乘法运算。因此，我们得到

              m[i,j]=m[i,k]+m[k+1,j]+ p(i-1)p(k)p(j)

    此递归公式假定最优分割点k是已知的，但实际上我们是不知道。不过，k只有j-i种可能的取值，即k=i,i+1,...,j-1。由于最优分割点必在其中，我们只需检查所有可能情况，找到最优者即可。

    因此，A(i)A(i+1)...A(j)的最小代价括号化方案的递归求解公式变为：

    ①如果i=j，m[i,j]=0

    ②如果i<j，m[i,j]=min{m[i,k]+m[k+1,j]+p(i-1)p(k)p(j)}  i<=k<j

 

m[i,j]的值给出了子问题最优解的代价，但它并未提供足够的信息来构造最优解。为此，我们用s[i,j]保存最优括号化方案的分割点位置k，即使得m[i,j]=m[i,k]+[k+1,j]+p(i-1)p(k)p(j)成立的k值。

3.计算最优解的值，通常采用自底向上的方法

我们采用自底向上表格法代替递归算法来计算最优代价。此过程假定矩阵Ai的规模为p(i-1)×pi(i=1,2,...,n)。它的输入是一个序列p=<p0,p1,...,pn>，其长度为p.length=n+1。过程用一个辅助表m[1..n,1..n]来保存代价m[i,j]，用另一个辅助表s[1..n-1,2..n](s[1,2]..s[n-1,n]这里i<j)记录最优值m[i,j]对应的分割点k。我们就可以利用表s构造最优解。

   对于矩阵A(i)A(i+1)...A(j)最优括号化的子问题，我们认为其规模为链的长度j-i+1。因为j-i+1个矩阵链相乘的最优计算代价m[i,j]只依赖于那么少于j-i+1个矩阵链相乘的最优计算代价。因此，算法应该按长度递增的顺序求解矩阵链括号化问题，并按对应的顺序填写表m。

4.利用计算出的信息构造一个最优解

因此，我们A(1..n)的最优计算方案中最后一次矩阵乘法运算应该是以s[1,n]为分界的A(1..s[1,n])*A(s[1,n]+1..n)。我们可以用相同的方法递归地求出更早的矩阵乘法的具体计算过程，因为s[1,s[1,n]]指出了计算A(1..s[1,n])时应进行的最后一次矩阵乘法运行；s[s[1,n]+1,n]指出了计算A(s[1,n]+1..n)时应进行的最后一次矩阵乘法运算。下面给出的递归过程可以输出<A(i),A(i+1),...,A(j)>的最优括号化方案。

 

代码

#include<iostream>
using namespace std;
#define INT_MAX 2147483647;
int const M=7;
void MATRIX_CHAIN_ORDER(int *p,int Length,int m[][M],int s[][M])
{
    int q,n=Length-1;
    for(int i=1;i<=n;i++) m[i][i]=0;
    for(int l=2;l<=n;l++)     /* 矩阵链的长度 */
    {
        for(int i=1;i<=n-l+1;i++) 
        {
            int j=i+l-1;         /* 等价于 l=j-i+1 */
            m[i][j]=INT_MAX;
            for(int k=i;k<=j-1;k++)
            {
                q=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
                if(q<m[i][j])
                {
                    m[i][j]=q;
                    s[i][j]=k;
                }
            }
        }
    }
}
void PRINT_OPTIMAL_PARENS(int s[][M],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[M]={30,35,15,5,10,20,25};
   int m[M][M],s[M][M];
   MATRIX_CHAIN_ORDER(p,M,m,s);
   cout<<"当n=6时最优解为: 
"<<m[1][6];
   cout<<"
括号化方案为：
"; 
   PRINT_OPTIMAL_PARENS(s,1,6);
   return 0;
}