【啊哈！算法】之七、动态规划---矩阵l链乘

原创

mb61b856e04bb98 2021-12-16 17:39:47 博主文章分类：Algorithm ©著作权

©著作权归作者所有：来自51CTO博客作者mb61b856e04bb98的原创作品，请联系作者获取转载授权，否则将追究法律责任

在学习矩阵链乘之前，先要知道矩阵之间是如何进行乘法运算的，如果你对这个还不是很了解，那么先去看一下线性代数中矩阵的乘法这一节！

OK，我们来说矩阵链乘，这是动态规划算法中的一个基础问题，在算法导论15章中也详细介绍了此算法，本文也是主要参考算法导论。

我们来看一下问题：给你N个矩阵，{A1，A2，A3，A4，......,AN}，这一组矩阵必须是可以相乘的，我们要计算这N个矩阵的乘积，即：A1A2A3..AN。因为我们的矩阵是满足结合律的，所以我们可以在计算矩阵链乘的时候可以有许多不同的计算次序。在这里，以不同的次序做乘法，其乘法的运算次数是不同的。

来举个例子：我们来对{A1，A2，A3}做连乘，假设他们的维数分别是：10*100,100*5,5*50, 那么我们按照（（A1A2）A3）的顺序来做乘法，求A1A2的乘积，也就是10*5的矩阵要做的也能算次数是10*100*5=5000次，而后再乘上A3 5*50，也就是10*5*50=2500次，总的乘法运算次数是 7500次。

我们再按照（A1（A2A3））的次序来计算，先是A2A3 的运算次数：100*5*50=25000次，而后要乘上A1 也就是 10*100*50=50000次，总的运算次数是75000次。 (⊙o⊙)… 通过上面的比较你应该会看出他们是十倍的关系！

简单一点描述矩阵链乘法的问题就是：给你n个矩阵构成的一个链{A1,A2,A3...An}，其中i=1,2,3...n，矩阵Ai的维数是：pi-1 * pi，对于乘积A1A2A3...An用一种最小化标量乘法次数的方式进行加全部括号。

ok，通俗点说，我们就是要寻找最优的解，就是用最少的乘法的次数，也就是流程中乘法次数最小。

我们分三步来解决：

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

动态规划的第一步就是寻找最优子结构，然后利用这个子结构就能根据子问题的最优解构造出原问题的一个最优解。

用记号A[i..j]表示乘积A[i]A[i+1]...A[j]求值的结果，其中i <= j。我们假设存在k，i <= k < j；A[i]A[i+1]...A[j]的任何全部加括号形式乘积在A[k]和A[k+1]这里分开，那么A[i]A[i+1]...A[j]的最优加全部括号的前子链A[i]A[i+1]...A[k]的加全括号必须是A[i]A[i+1]...A[k]的一个最优加全部括号。

这样一来就可以利用多得到的最优子结构来说明能够根据子问题的最优解来构造原问题的一个最优解。

2、一个递归求解

下面我们就要根据子问题的最优解来递归定义一个最优解的代价。

我们先设m[i][j]为计算矩阵A[i..j]所需的标量乘法运算次数的最小值，对于整个问题，计算A[1..n]的最小代价就是m[1][n]。

这里，如果i==j 那么这就是一个普通问题了，矩阵链就有一个矩阵，所以不需要进行任何标量乘法来计算乘积。

但是当i<j的时候，我们必须要来计算m[i...j]，我们就要利用步骤1中得到的最优解的结构了。

假设最优加全部括号将乘积A[i]A[i+1]...A[j]从A[k]和A[k+1]之间分开，i <= k < j。因此m[i...j]就等于计算子乘积Ai...k和Ak+1...j的代价，最后还要加上一个条件是这两个矩阵相乘的代价的最小值。我们计算Ai...kAk+1...j 要做pi-1pkpj 次标量乘法。

最后得出的公式是：

对于第二个公式，实际上是：A[i..j]的最优次序在Ak和Ak+1之间断开，其中i<=k<j。也就是我们要利用最优子结构来计算。所以最后得带m[i,j] = m[i,k] + m[k+1,j] + pi-1 pk pj；

在学习矩阵链乘之前，先要知道矩阵之间是如何进行乘法运算的，如果你对这个还不是很了解，那么先去看一下线性代数中矩阵的乘法这一节！

OK，我们来说矩阵链乘，这是动态规划算法中的一个基础问题，在算法导论15章中也详细介绍了此算法，本文也是主要参考算法导论。

ok，通俗点说，我们就是要寻找最优的解，就是用最少的乘法的次数，也就是流程中乘法次数最小。

我们分三步来解决：

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

动态规划的第一步就是寻找最优子结构，然后利用这个子结构就能根据子问题的最优解构造出原问题的一个最优解。

这样一来就可以利用多得到的最优子结构来说明能够根据子问题的最优解来构造原问题的一个最优解。

2、一个递归求解

下面我们就要根据子问题的最优解来递归定义一个最优解的代价。

我们先设m[i][j]为计算矩阵A[i..j]所需的标量乘法运算次数的最小值，对于整个问题，计算A[1..n]的最小代价就是m[1][n]。

这里，如果i==j 那么这就是一个普通问题了，矩阵链就有一个矩阵，所以不需要进行任何标量乘法来计算乘积。

但是当i<j的时候，我们必须要来计算m[i...j]，我们就要利用步骤1中得到的最优解的结构了。

3、计算最优代价

在这里用的是子问题的自底向上的表格法来计算最优代价。这部分的代码会稍后在文章最后统一给出，其函数名是:MATRIXCHAINORDER。来说一下此函数，函数由3个循环嵌套给出，其时间复杂度是O(n3)。并且还需要O(n2)的空间来保存m和s。

我们可以看到m[][]的值可以根据两个方向的交点看出。其中表项m[i,j]根据乘积pi-1pkpj和所在的m[i,j] 西南方向和东南方向的表项计算得出。

4、构造一个最优解

上个步骤中的函数确定了计算矩阵链乘积所需要的最优标量乘法次数，但是没有直接给出如何对矩阵进行相乘。函数中s[i][j]的数据可以知道计算乘积Ai..Aj 在Ak和Ak+1之间，进行分裂以取得最优加全部括号时的k值。也就是说：最优的方式是（Ai...Ak）（Ak+1...Aj）；

其实在s之中，s[1, s[1, n]]决定计算A1..s[1,n] 中最后的矩阵乘法，而s[s[1,n]+1, n]则决定计算As[1, n]+1...n中最后的矩阵乘法。

在这里有一个函数是输出Ai....Aj 的一个最优加全部括号形式 PRINTOPTIMALPARENS；

下面请看源代码：

#include<iostream>
using namespace std;
/*
 *p[i]:矩阵Ai的列数或Ai-1的行数
 *m[i][j]:纪录Ai - Aj的矩阵连乘的最小代价
 *s[i][j]:纪录Ai - Aj之间得到最小连乘代价时的分割点
 *MAXLENGTH 矩阵链的长度
 */
#define MAXLENGTH 10
void MATRIXCHAINORDER(int p[], int m[][MAXLENGTH], int s[][MAXLENGTH])
{
  int j;
  int temp;
  for(int i = 0; i < MAXLENGTH; i++)   //i=j 的情况
  {
    m[i][i] = 0;
  }
  //i!=j的情况
  for(int len = 2; len <= MAXLENGTH; len++)
  {
    //i 从第i个矩阵开始的len个矩阵如何划分是最优
    for(int i = 0; i < MAXLENGTH - len + 1; i++)
    {
      //j是第i个数len个矩阵就是j，从i到j表示他们之间的len个矩阵如何划分
      j = i + len - 1;
      //m[i][j]是存储从i到j使用最佳划分所得的最优结果
      m[i][j] = INT_MAX; //无穷大
      //k是介于i到j-1之间的数，利用k对i到j之间的矩阵进行试探性划分，
      //从而找到最优解。
      for(int k = i; k < j; k++)
      {
        temp = m[i][k] + m[k+1][j] + p[i]*p[k+1]*p[j+1];
        if(temp < m[i][j])
        {
          m[i][j] = temp;
          s[i][j] = k;
        }
      }
    }
  }
}
void PRINTOPTIMALPARENS(int s[][MAXLENGTH], int i, int j)
{
  if(i == j)
  {
    cout << "A" << i;
  }
  else
  {
    cout << "(";
    PRINTOPTIMALPARENS(s, i, s[i][j]);
    PRINTOPTIMALPARENS(s, s[i][j]+1, j);
    cout << ")";
  }
}
//重叠子结构
int RECURSIVEMATRIXCHAIN(int p[], int m[][MAXLENGTH], int s[][MAXLENGTH], int i, int j)
{
  int temp;
  if(i == j)
    return 0;
  int u = RECURSIVEMATRIXCHAIN(p, m, s, i, i) + RECURSIVEMATRIXCHAIN(p, m, s, i + 1, j) + p[i - 1]*p[i]*p[j];
  s[i][j] = i;
  for(int k = i + 1; k < j; k++)
  {
    temp = RECURSIVEMATRIXCHAIN(p, m, s, i, k) + RECURSIVEMATRIXCHAIN(p, m, s, k + 1, j) + p[i - 1]*p[k]*p[j];
    if(temp < u)
    {
      u = temp;
      s[i - 1][j - 1] = k;
    }
  }
  return u;
}
//备忘录
int LOOKUPCHAIN(int p[], int m[][MAXLENGTH], int s[][MAXLENGTH], int i, int j)
{
  int temp;
  //检查m[i][j]
  if(m[i][j] < INT_MAX)
    return m[i][j];
  if(i == j)
  {
    m[i][j] = 0;
  }
  else
  {
    //递归寻找k，试探划分，最终找到最优解
    for(int k = i; k < j; k++)
    {
      temp = LOOKUPCHAIN(p, m, s, i, k) + LOOKUPCHAIN(p, m, s, k+1, j) + p[i-1]*p[k]*p[j];
      if(temp < m[i][j])
      {
        m[i][j] = temp;
        s[i - 1][j - 1] = k;
      }
    }
  }
  return m[i][j];
}
void MEMOIZEDMATRIXCHAIN(int p[], int m[][MAXLENGTH], int s[][MAXLENGTH], int n)
{
  for(int i = 0; i < n; i++)
  {
    for(int j = i; j < n; j++)
    {
      m[i][j] = INT_MAX;
    }
  }
  //自底向上递归
  m[0][n - 1] = LOOKUPCHAIN(p, m, s, 0, n - 1);
}
int main(void)
{
  int p[MAXLENGTH+1] = {34, 15, 35, 25, 23, 21, 67};
  int m[MAXLENGTH ][MAXLENGTH];
  int s[MAXLENGTH ][MAXLENGTH];
  MATRIXCHAINORDER(p, m, s);
  PRINTOPTIMALPARENS(s, 0, 6);
  cout << endl;
  MEMOIZEDMATRIXCHAIN(p, m, s, 7);
  PRINTOPTIMALPARENS(s, 0, 6);
  cout << endl;
  RECURSIVEMATRIXCHAIN(p, m, s, 1, 7);
  PRINTOPTIMALPARENS(s, 0, 6);
  cout << endl;
  return 0;
}