矩阵连乘 最优计算次序 动态规划 图文详解

矩阵连乘最优次序求解

没有公式，放心食用。

问题概要

给定矩阵 A1 A2 A3 A4 A5 A6 及其行列维度
求解这些矩阵的最优计算次序，使得乘法运算次数最少
其中 A1 A2 A3 A4 A5 A6 从左至右每个矩阵两两满足“左边的列维度等于右边的行数”

前导知识

一、矩阵运算规则
	设矩阵 A 和 B 的行列维数分别为 (a, b)，(c, d)；
	A和B能够进行矩阵乘法运算的前提为 b==c；（※）
	且运算过后得到的矩阵 C 行列维度为 (a, d)。
	【举例】
	A=[[1, 2, 3],         B=[[1, 2],
	   [1, 2, 2]]            [1, 3],
	                         [1, 4]]  
	则 A·B=[[1*1+2*1+3*1, 1*2+2*3+3*4],
	        [1*1+2*1+2*1, 1*2+2*3+2*4]]
	        
二、矩阵乘法中，乘法的运算次数
	假设A的维度为 Sa=(2, 3)，B的维度为 Sb=(3, 2)；
	则其总的乘法运算次数为 Sa[0]*Sb[0]*Sb[1]=2*3*2=12 次；
	即 “A的行数” 乘 “B的行数” 乘 “B的列数”。
	
	推广到多个矩阵时，
	对于A(i)·A(i+1)，其运算次数为 A(i)的行数 乘 A(i+1)的行数 乘 A(i+1)的列数。

	在这里再强调一下，问题所给出的若干个矩阵，
	相邻两个矩阵“均满足能够进行矩阵运算的前提”（※）。

	所以，对于 A1 A2 A3 A4 A5 A6 ，
	我们想取某个矩阵的行数，只需要去寻找与它相邻的前一个矩阵的列数不就可以了吗？
	因为(a, b)，(c, d)，有b==c，
	那么假设每个矩阵的形状为(a, b) (c, d) (e, f) (g, h) (i, j) (k, l)，
	可以得到b==c，d==e，f==g，h==i，j==k，
	现在我们将每个矩阵的行 用 与其相邻的前一个矩阵的列代替，有
	(a, b) (b, d) (d, f) (f, h) (h, j) (j, l)
	既然除了第一个矩阵，任意一个矩阵的行数都能通过取与其相邻的前一个矩阵的列来获得，
	那么我们不就首先简化了在计算机存储矩阵形状的方式啦！

	定义一个int型数组 p[7] 存取每个矩阵的列，
	但是第一个矩阵无法从它“前一个矩阵”获取行数，
	所以，需要将第一个矩阵的行存储在该数组的第一个位置，即 p[0]位置。
	依次存储得到数组 p[7]={30, 35, 15, 5, 10, 20, 25}
	                      ↑
	                第一个矩阵的行

动态规划

依然引用上面的例子，
对于这六个矩阵 A1 A2 A3 A4 A5 A6 ，
想要获得计算的最优次序使得乘法运算次数最少，怎么办？

容易知道，如果枚举，将会有很高的时间复杂度，在限制时间的前提下确实不是一个好的选择。
所以我们得考虑一下动态规划了。

但是在讲到动态规划实现方法之前，我们先来了解一下如何枚举获得最优计算次序。

一、枚举
	矩阵乘法无非就是以两个矩阵为单位进行运算，
	那我们首先要找到能够让计算次数最少的两相邻矩阵。
	是先算 A1A2？还是 A2A3？还是A3A4？A4A5？A5A6？
	
	（1）第一层：找出第一对运算次数最少的两个矩阵
		 (A1 A2) A3 A4 A5 A6，
		 A1 (A2 A3) A4 A5 A6，
		 A1 A2 (A3 A4) A5 A6，
		 A1 A2 A3 (A4 A5) A6，
		 A1 A2 A3 A4 (A5 A6)，
		记第i个矩阵和第j个矩阵相乘为A(i,j)，对上面的循环得出的序列进行替代，得到
		 A(1,2) A3 A4 A5 A6，
		 A1 A(2,3) A4 A5 A6，
		 A1 A2 A(3,4) A5 A6，
		 A1 A2 A3 A(4,5) A6，
		 A1 A2 A3 A4 A(5,6)，
	（2）第二层：在前一层得到的第一个运算次数最少的两个矩阵的基础上，找第二对
		 A((1,2),3) A4 A5 A6，
		 A(1,2) A(3,4) A5 A6，
		 A(1,2) A3 A(4,5) A6，
		 A(1,2) A3 A4 A(5,6)，
		 		...
		 A1 A((2,3),4) A5 A6，
		 		...
		 A1 A2 A((3,4),5) A6，
		 		...
		 A1 A2 A3 A((4,5),6)，
		 		...
	（3）第三层：找第三对
		 A(((1,2),3),4) A5 A6，
		 		...
		 A1 A(((2,3),4),5) A6，
		 		...
		 A1 A2 A(((3,4),5),6)，
		 		...
	（4）第四层：找第四对
	     A(((1,2),3),4),5) A6，
	     		...
		 A1 A((((2,3),4),5),6)，
		 		...
	（5）第五层：找第五对
		 A(((((1,2),3),4),5),6)
	以此类推。

二、动态规划-自底向上方法
	【变量解释】
	（1）m[i][j]：第i个矩阵到第j个矩阵的最少运算次数
	（2）col[i]：第i个矩阵的列数，其中col[0]存的是第一个矩阵的行数
	（3）k：第i个矩阵与第j个矩阵的分界下标
	（4）b：第i个矩阵与第j个矩阵之间（不包括i和j）有（b-1）个矩阵
	（5）s[i][j]：存第i个矩阵和第j个矩阵之间的最优划分点
    【基本步骤】
    （1）分别求出1个矩阵相乘需要的运算次数、2个矩阵相乘需要的运算次数
         ...
         一直到num个矩阵相乘需要的运算次数，
         用数组m[num][num]来存储第i个矩阵到第j个矩阵的最少运算总次数
         注意：下标从m[1][1]开始.
         
    （2）其中，每次仅仅求出n个矩阵相乘所需要的运算次数后，得到的矩阵m还并不是最少计算次数结果，
         所以，单次求n个矩阵相乘次数后，还需要再以k为第i个矩阵到第j个矩阵的分界点，
         分别比较m[i][j]与m[i][k]+m[k+1][j]的大小，这里的k范围是[i+1, j-1]闭区间。
         意思是第i个矩阵和第k个矩阵相乘后，再与第k+1个矩阵到第j个矩阵相乘得到的矩阵相乘。
         即 （Ai*Ai+1*...Ak）*（Ak+1*...Aj），由于是自底向上算法，（Ai*Ai+1*...Ak）、（Ak+1*...Aj）分别已经是最优计算次序。

图文实例

注意第一层的每轮循环，求n个矩阵相乘得到的计算次数时，是默认计算从左到右矩阵相乘得到的次数，就是 A1 · A2
后再乘A3，再乘A4…乘到Anum。第二层循环才开始比较不同划分点计算的次数，得到最少的计算次数存入数组m，最优计算次序存入数组s。

①初始化数组m[num][num]

②第一次循环，求相邻两个矩阵相乘得到的 最少计算次数 和 最优计算次序
可以发现，第二层循环（如下）没有进行，那是因为两个矩阵之间没有矩阵，无法进行划分，已经是最优计算次序了。

for (int k = i + 1; k < j; k++) {
    int t = m[i][k] + m[k + 1][j] + col[i - 1] * col[k] * col[j];
    if (t < m[i][j]) {
        m[i][j] = t;
        s[i][j] = k; //矩阵最优划分为(i,k)和(k+1,j)两个部分相乘次数最少
    }
}

②第二次循环，计算相邻三个矩阵的 最少计算次数 和 最优计算次序
注意：这次要进行比较，即三个矩阵当中，是选择先算前两个矩阵还是选择先算后两个矩阵，这两种计算次序会影响三个矩阵的总计算次数。

先分别计算三个相邻矩阵的计算次数（还没进行第二层循环，此时数组m中的计算次数还不是最优），得到下表

为了方便图示，这部分计算并没有算上中间嵌套的第三层for循环，而是直接循环完第二层for循环得到的结果。那如果第二层for的每轮循环同时算上第三层for呢？第三层for实际上是用来检查是否还存在比当前计算次数更小的计算次序的（拿三个矩阵的最优连乘次序举例，比如我们现在在未进入第三层for之前算了m[1][3]的计算次数=m[2][3]+p[0]*p[1]*p[2]，这个计算公式的意思是先计算矩阵1和矩阵2，最后计算前两个矩阵计算得到的矩阵与矩阵3，那在第三层for的时候就得检查先计算矩阵2和矩阵3，再计算矩阵1和矩阵2、3相乘得到的矩阵的连乘结果）。第一层for的第二轮循环完毕后可得到下表：

其余循环以此类推。

完整代码

#include <bits/stdc++.h>
using namespace std;
//相乘的矩阵个数为(num-1)个
#define num 7
//col[0]存的是第一个矩阵的行数，其余分别为第1、第2...第num个矩阵的列数
//只存列数是因为两个矩阵可以运算的前提是：前一个矩阵的 列数 等于后一个矩阵的 行数。所以如果想取当前矩阵的 行数，只需要查找矩阵的 列数 即可。
//例如A1.shape=(30,35)，A2.shape=(35,15)，此时A1列数等于A2行数，所以才能进行矩阵乘法运算，并且可以通过取A1的列数得知A2的行数。
int col[num] = {30, 35, 15, 5, 10, 20, 25};
int m[num][num], s[num][num];

void matrixChain() {
    //b表示当前计算的矩阵连乘个数
    for (int b = 2; b < num; b++) {
        //下标从1开始
        //(num-b+1)表示按照从左到右的顺序，相邻b个矩阵连乘的组合方式有多少个
        //例如对于A1 A2 A3矩阵连乘，假设当前b=2，那么b个相邻的矩阵连乘组合方法有A1A2和A2A3
        //对于A1 A2 A3 A4矩阵连乘，设b=2，那么b个相邻的矩阵连乘组合方法有A1A2、A2A3、A3A4
        for (int i = 1; i < num - b + 1; i++) {
            //j与第i个矩阵相隔(b-1)个矩阵，例如A1A2A3，设b=2，则A1与A3相隔(1+2-1)个矩阵
            int j = i + b - 1;
            //获取所有可能的(num-b+1)个相邻b个矩阵连乘组合的计算次数
            //例如，设b=2，两个两个矩阵连乘，那么对于A1A2A3A4的矩阵，两两相邻矩阵连乘的方式有A1A2 A2A3 A3A4，
            //i为b个矩阵中排在第一个的矩阵的下标，j为与i相隔了(b-1)个矩阵的矩阵，设i=1，j=1+2-1=2，计算的就是A1A2两个矩阵连乘的计算次数
            //再例如设i=2，j=3+2-1=3，计算的就是A2A3两个矩阵连乘的次数；如果b=3，那么j=2+3-1=4，计算的就是A2A3A4三个矩阵连乘的计算次数
            m[i][j] = m[i + 1][j] + col[i - 1] * col[i] * col[j];
            s[i][j] = i;//表示此时第i到j个矩阵以第i个矩阵为划分点划分成i~i与(i+1)~j两个部分
            for (int k = i + 1; k < j; k++) {
                int t = m[i][k] + m[k + 1][j] + col[i - 1] * col[k] * col[j];
                if (t < m[i][j]) {
                    m[i][j] = t;
                    s[i][j] = k; //矩阵最优划分为(i,k)和(k+1,j)两个部分相乘次数最少
                }
            }
        }
    }
}

int main() {
    matrixChain();
    for (int i = 1; i < num; i++) {
        for (int j = 1; j < num - 1; j++)
            cout << m[i][j] << " ";
        cout << m[i][num - 1] << endl;
    }
    cout << endl;
    for (int i = 1; i < num; i++) {
        for (int j = 1; j < num - 1; j++)
            cout << s[i][j] << " ";
        cout << s[i][num - 1] << endl;
    }
    return 0;
}