动态规划---例题1.矩阵连乘问题

2024-02-09 11:29:04

动态规划---例题1.矩阵连乘问题

一.问题描述

矩阵A和B可乘的条件是矩阵A的列数等于矩阵B的行数.若A是一个p×q的矩阵，B是一个q×r的矩阵，则其乘积C=AB是一个p×r的矩阵.其标准计算公式为：

计算C=AB总共需要pqr次的数乘.
给定n个矩阵{A1,A2,…,An}.其中Ai与Ai+1是可乘的，i=1,2,…,n-1.要求计算出这n个矩阵的连乘积A1A2…An.

二.解题思路

矩阵乘法满足结合律，故连乘积的计算可以有许多不同的计算次序.这种计算次序可以用加括号的方式来确定.若一个矩阵连乘积的计算次序已完全确定，也就是说该连乘积已完全加括号，则我们可以通过反复调用两个矩阵相乘的标准算法计算出矩阵连乘积.完全加括号的矩阵连乘积可递归地定义为：

1.单个矩阵是完全加括号的；
2.若矩阵连乘积A是完全加括号的，则A可表示为两个完全加括号的矩阵连乘积B和C的乘积并加括号，即A=(BC).

例如，矩阵连乘积A1A2A3 A4可以有以下5种不同的完全加括号方式：

(A1(A2(A3A4)))
(A1((A2A3)A4))
((A1A2)(A3A4))
((A1(A2A3))A4)
(((A1A2)A3)A4)

每一种完全加括号方式对应于一种矩阵连乘积的计算次序，而这种计算次序与计算矩阵连乘积的计算量有着密切的关系.
为了说明在计算矩阵连乘积时加括号方式对整个计算量的影响，我们来看一个计算3个矩阵{A1，A2，A3}的连乘积的例子.
设这3个矩阵的维数分别为10×100，100×5和5×50.若按第一种加括号方式((A1A2)A3)来计算，总共需要
10×100×5+10×5×50=7500次的数乘.

若按第二种加括号方式(A1(A2A3))来计算，则需要的数乘次数为
100×5×50+10×100×50=75000

后者的计算量是前者的10倍.可见，在计算矩阵连乘积时，计算次序对计算量影响很大.§矩阵连乘积的最优计算次序问题，即对于给定的相继n个矩阵{A1,A2,…,An}(其中Ai的维数为pi-1×p****i ，i=1,2,…,n)，如何确定计算矩阵连乘积A1A2…An的一个计算次序(完全加括号方式)，使得依此次序计算矩阵连乘积需要的数乘次数最少.

解这个问题的最容易想到的方法是穷举搜索法.也就是列出所有可能的计算次序，并计算出每一种计算次序相应需要的计算量，然后找出最小者.
然而，这样做计算量太大.

事实上，对于n个矩阵的连乘积，设有P(n)个不同的计算次序.由于可以首先在第k个和第k+1个矩阵之间将原矩阵序列分为两个矩阵子序列，k=1,2,…,n-1；然后分别对这两个矩阵子序列完全加括号；最后对所得的结果加括号，得到原矩阵序列的一种完全加括号方式.所以关于P(n)，我们有递推式如下：

解此递归方程可得,P(n)实际上是Catalan数，
即P(n)=C(n-1)，其中，

也就是说，P(n)随着n的增长是指数增长的.因此，穷举搜索法不是一个有效算法.
下面我们来考虑用动态规划法解矩阵连乘积的最优计算次序问题.此问题是动态规划的典型应用之一.

1.分析最优解的结构

将矩阵连乘积AiAi+1…Aj简记为Aij.先看计算A1 n的一个最优次序.设这个计算次序在矩阵Ak和Ak+1之间将矩阵链断开，1<=k<n，则完全加括号方式为((A1…Ak)(Ak+1…An)).照此，我们要先计算A1…k和Ak+1…n，然后，将所得的结果相乘才得到A1 n.显然其总计算量为计算A1…k的计算量加上计算Ak+1…n的计算量，再加上A1…k与Ak+1…n相乘的计算量.

问题的关键特征：计算A1 n的一个最优次序所包含的计算A1…k的次序也是最优的.事实上，若有一个计算A1…k的次序需要的计算量更少，则用此次序替换原来计算A1…k的次序，得到的计算A1…n的次序需要的计算量将比最优次序所需计算量更少，这是一个矛盾.同理，计算A1…n的一个最优次序所包含的计算矩阵子链Ak+1…n的次序也是最优的.

最优解包含其子问题的最优解，这种性质称为最优子结构性质.另外，该问题显然满足无后向性，因为前面的矩阵链的计算方法和后面的矩阵链的计算方法无关.

2.建立递归关系

对于矩阵连乘积的最优计算次序问题，设计算Ai…j ,1≤i≤j≤n，所需的最少数乘次数为m[i,j]，原问题的最优值为m[1,n].

当i=j时，Ai…j=Ai为单一矩阵，无需计算，因此m[i,i]=0，i=1,2,…,n ；
当i<j时，可利用最优子结构性质来计算m[i,j].事实上，若计算Ai…j的最优次序在Ak和Ak+1之间断开，i≤k<j，则：m[i,j] = m[i,k] + m[k+1,j] + pi-1pkpj

断开点A的位置还未确定.不过k的位置只有j-i个可能，即k∈{i,i+1,…,j-1}.因此k是这j-i个位置中计算量达到最小的那个位置.从而m[i,j]可递归地定义为：

m[i,j]给出了最优值，即计算Ai…j所需的最少数乘次数.同时还确定了计算Ai…j的最优次序中的断开位置k，也就是说，对于这个k有m[i,j] = m[i,k] + m[k+1,j] + pi-1pkpj .若将对应于m[i,j]的断开位置k记录在s[i,j]中，则相应的最优解便可递归地构造出来.

三.计算递归值

根据m[i,j]的递归定义(2.1)，容易写一个递归程序来计算m[1,n].然而简单地递归计算将耗费指数计算时间.注意：在递归计算过程中，不同的子问题个数只有θ(n2)个.事实上，对于1≤i≤j≤n不同的有序对(i,j)对应于不同的子问题.因此，不同子问题的个数最多只有

个.可见，许多子问题被重复计算多次.

用动态规划算法解此问题，可依据递归式(2.1)以自底向上的方式进行计算，在计算过程中，保存已解决的子问题答案，每个子问题只计算一次，而在后面需要时只要简单查一下，从而避免大量的重复计算，最终得到多项式时间的算法.
下面所给出的计算m[i,j]动态规划算法中，输入是序列P={p0,p1,…,pn}，输出除了最优值m[i,j]外，还有使
m[i,j] = m[i,k] + m[k+1,j] + pi-1pkpj
达到最优的断开位置k=s[i,j]，1≤i≤j≤n;

该算法按照

m[1,1]
m[2,2] m[1,2]
m[3,3] m[2,3] m[1,3]
... ... ...
m[n,n] m[n-1,n] ... .... ... m[1,n]

的顺序根据公式(2.1)计算m[i, j].
该算法的计算时间上界为O(n3).算法所占用的空间显然为O(n2).由此可见，动态规划算法比穷举搜索法要有效得多.

4.构建最优解

算法MatrixChain只是计算出了最优值，并未给出最优解.也就是说，通过MatrixChain的计算，我们只知道计算给定的矩阵连乘积所需的最少数乘次数，还不知道具体应按什么次序来做矩阵乘法才能达到数乘次数最少.
然而，它己记录了构造一个最优解所需要的全部信息.

因为s[i, j]中的数k表明计算矩阵链Ai…j的最佳方式应在矩阵Ak和Ak+1之间断开，即最优的加括号方式应为(A1...k)(Ak+1…n).因此，从s[i,j]记录的信息可知计算A1…n的最优加括号方式为 (A1…s[1,n])(As[1,n]+1…n).
而计算A1…s[1,n]的最优加括号方式为(A1…s[1,s[1,n]])(As[1,s[1,n]]+1…s[1,n]).同理可以确定计算As[1,n]+1…n的最优加括号方式在s[s[1,n]+1,n]处断开.…照此递推下去，最终可以确定As[1,n]+1…n的最优完全加括号方式，即构造出问题的一个最优解.

代码如下:

// 矩阵乘法链动态规划算法
// 矩阵连乘,动态规划迭代实现(自底向上)
// A1 30*35 A2 35*15 A3 15*5 A4 5*10 A5 10*20 A6 20*25
// p[0-6] = {30,35,15,5,10,20,25}
#include<bits/stdc++.h>
using namespace std;
const int L = 7;
int MatrixChain(int *p, int n, int ** m, int **s);
void Traceback(int i, int j, int **s);  //构造最优解

int main()
{
    int p[L] = {30,35,15,5,10,20,25};
    int **s = new int *[L];
    int **m = new int *[L];
    for(int i=0; i<L; i++)
    {
        s[i] = new int[L];
        m[i] = new int[L];
    }

    cout<<"矩阵的最少计算次数为"<<MatrixChain(p, 6, m, s)<<endl;
    cout<<"矩阵最优计算次序为:"<<endl;
    Traceback(1, 6, s);
    system("pause");
    return 0;
}
int MatrixChain(int *p, int n, int ** m, int ** s)  //T(n) = O(n^3)
{
    for(int i=1; i<=n; i++) m[i][i] = 0;
    for(int r=2; r<=n; r++) //r个矩阵连乘.r为当前计算的链长(子问题规模)
        for(int i=1; i<=n-r+1; i++) //n-r+1为最后一个r链的前边界
        {
            int j=i+r-1;  //计算前边界为i,链长为r的后边界
            m[i][j] = m[i+1][j] + p[i-1]*p[i]*p[j];  //将链ij划分为A(i)*(A[i+1:j]). ]m[i][i]=0省略.
            s[i][j] = i;
            for(int k=i+1; k<j; k++) //分割位置从i+1不断向后移动,将链ij划分为(A[i:k]*A[k+1:j])
            {
                int t = m[i][k] +m[k+1][j] + p[i-1]*p[k]*p[j];
                if(t < m[i][j])
                {
                    m[i][j] = t;
                    s[i][j] = k;
                }
            }
        }
    return m[1][L-1];
}
void Traceback(int i, int j, int **s)
{
    if(i==j) return ;
    Traceback(i, s[i][j], s);
    Traceback(s[i][j]+1, j, s);
    cout<<"Mutiply A"<<i<<","<<s[i][j];
    cout<<" and A"<<(s[i][j]+1)<<","<<j<<endl;
}

参考我的老师毕方明《算法设计与分析》课件.

欢迎大家访问个人博客网站---乔治的编程小屋,和我一起为大厂offer努力!

码农公寓