矩阵连乘求解优化

前言

从旭东的博客 看到一篇博文：矩阵连乘最优结合 动态规划求解，挺有意思的，这里做个转载【略改动】。

问题

矩阵乘法满足结合律，但不满足交换律。例如矩阵$A_{ab}, B_{bc}, C_{cd}$ 连乘得到矩阵$S_{ad}$

[ S_{ad}=A_{ab} B_{bc} C_{cd} ]

实际运算可以先将矩阵A和B相乘，S=(AB)C，也可以先将矩阵B和C相乘，S=A(BC) 。这两种算法的计算量是否一样呢？

先看A×B的计算量。A×B是a行c列矩阵，共a×c个元素。b是A、B两矩阵的乘法接口，积矩阵的每个元素都需经过b次乘法运算，那么A×B需要a×b×c次乘法运算。这里注意喽，两个b维向量点乘需要的计算次数不是$b^2$次，而是b次。

(AB)C乘法运算次数为：abc+acd=ac(b+d)

A(BC)乘法运算次数为：bcd+abd=bd(a+c)

由于a、b、c、d的取值是任意的， 所以ac(b+d)和bd(a+c)不可能保持相同，也即结合律可以用于矩阵乘法优化。那么，对于连乘A_iA_i+1…A_j，如何找出最优的运算顺序？

分析

这个问题长得像递归，但是用递归似乎不好做，因为子问题里面有很多分支，那子问题的返回值是什么，你返回哪个分支的返回值？而且就算能实现，肯定会对一些相同的子问题重复计算。

由于计算n个矩阵的连乘，总是通过计算更少个数的矩阵连乘实现的，最终化归到计算两个矩阵的相乘，所以考虑从两个矩阵相乘着手。这样自底向上，就把问题解决了。当我们把所有的两个矩阵相乘的乘法运算次数求出后，便可以去计算三个矩阵相乘的乘法运算次数。这时，三个矩阵的连乘相当于两个矩阵相乘。这样，所有的矩阵连乘都是两个矩阵相乘，其乘法运算次数为两个矩阵各自的乘法运算次数之和再加上两个矩阵相乘的乘法运算次数。对于两个矩阵A、B相乘的乘法运算次数$eta$，上面A×B的示例已经讲过了。进一步可以抽象为

[eta=行*接口*列]

我们将矩阵连乘表达式视为一系列子式，为确定最佳计算方案，我们需要知道每个子式的情况。每一个子式均由起点和跨度唯一决定，起点和跨度刚好构成两个维度，于是可以用n阶方阵存储子式的一些情况。比如索引[1][2]表示起点为1跨度为2的子式，里面存的值专属于这个子式。起点至多有n个，最小为0，最大为n-1。跨度至多n种，最小为0，最大为n-1。用两个n阶方阵分别存储每个子式的最少计算次数及分割点。

C++提供了vector模板类，作为数组的一个替代品，我们用vector实现n阶方阵。首先确定跨度，再确定起点，构成一个二级嵌套循环，通过起点的平移确定子式。子式确定后，再在内部嵌套一级循环，遍历子式的可能的分割点，并保存子式的最少计算次数及对应分割点。最后我们得到最大的子式，也就是连乘本身的最少运算次数及分割方案。

代码如下：

#include <iostream>
#include <vector>
#include <algorithm>
#include <limits.h>
#include <string>

using namespace std;

//计算矩阵的分割方案（基于向量）
int calculate_M(vector<vector<int> >&num, vector<pair<int, int> > &data, vector<vector<int> > &points) {
    int matrixNum = data.size();
    int span;//起止点间隔距离，表示跨度
    int start;//起点
    int end;//终点
    int spiltPoint;//分割点
    int multiplyTimes;//临时存储子式的乘法运算次数
    for (span = 1; span < matrixNum; span++) {  ///间隔距离 相邻矩阵的间隔距离为1
        for (start = 0; start < matrixNum - span; start++) {  ///起点平移
            //子式确立，下面开始计算这个子式的最少乘法运算次数，初值先设为最大整数
            end = start + span;
            num[start][end] = INT_MAX;
            for (spiltPoint = start; spiltPoint < end; spiltPoint++) {
                
                multiplyTimes = num[start][spiltPoint] + num[spiltPoint + 1][end] + data[start].first*data[spiltPoint].second*data[end].second;
                if (multiplyTimes < num[start][end]) {
                    points[start][end] = spiltPoint;  ///记录分割点
                    num[start][end] = multiplyTimes;   ///记录最少乘法次数
                }
            }
        }
    }
    return 0;
}

//根据记录的分割点，生成最后的矩阵相乘表达式
string make_result(vector<vector<int> > &points, int t1, int t2) {
    if (t1 == t2)
        return string(1, 'A' + t1);
    int split = points[t1][t2];
    return "(" + make_result(points, t1, split) + "*" + make_result(points, split + 1, t2) + ")";
}


int main()
{
    int matrixNum=0;///连乘的矩阵个数
    vector<pair<int, int>> data; ///存储矩阵的尺寸

    //固定输入
    data.push_back(make_pair(10, 100));  //A
    data.push_back(make_pair(100, 5));   //B
    data.push_back(make_pair(5, 25));    //C
    data.push_back(make_pair(25, 15));   //D
    data.push_back(make_pair(15, 20));   //E
    matrixNum = 5;

    //为记录乘法次数和分割点的n阶矩阵分配空间并初始化
    vector<vector<int> > num(matrixNum, vector<int>(matrixNum));
    vector<vector<int> > points(matrixNum, vector<int>(matrixNum));
    for (int i = 0; i < matrixNum; i++) {
        for (int j = 0; j < matrixNum; j++) {
            points[i][j] = -1;
            num[i][j] = 0;
        }
    }

    calculate_M(num, data, points);
    cout << make_result(points, 0, matrixNum - 1) << "	最少乘法次数为:" << num[0][matrixNum - 1] << endl;
    
    cin >> matrixNum;
    return 0;
}

 代码输出为：

((A*B)*((C*D)*E))       最少乘法次数为:9375

如果从控制台输入矩阵尺寸的话，可以输入第一个矩阵的行数和各矩阵的列数。直接使用这个序列存储矩阵尺寸。

#include <iostream>
#include <vector>
#include <algorithm>
#include <limits.h>
#include <string>

using namespace std;

int* uncertainInput(int& inputNum);

//计算矩阵的分割方案（基于数组）
int calculate_M1(int*data, int matrixNum, vector<vector<int> >&num, vector<vector<int> > &points) {
    
    int span;//起止点间隔距离，表示跨度 最大跨度比总的矩阵数小1
    int start;//起点
    int end;//终点
    int spiltPoint;//分割点 分割线在分割点与下一点之间
    int multiplyTimes;//临时存储子式的乘法运算次数
    int rows, columns, interfaces;//行 列 接口

    for (span = 1; span < matrixNum-1; span++) {                ///间隔距离 相邻矩阵的间隔距离为1
        for (start = 1; start + span < matrixNum; start++) {    ///起点从第一个矩阵的列开始
            //子式确立，下面开始计算这个子式的最少乘法运算次数，初值先设为最大整数
            end = start + span;
            num[start][end] = INT_MAX;
            for (spiltPoint = start; spiltPoint < end; spiltPoint++) {
                rows = *(data + start-1);
                columns = *(data + end);
                interfaces = *(data + spiltPoint);

                multiplyTimes = num[start][spiltPoint] + num[spiltPoint + 1][end] + rows * interfaces * columns;
                if (multiplyTimes < num[start][end]) {
                    points[start][end] = spiltPoint;  ///记录分割点
                    num[start][end] = multiplyTimes;   ///记录最少乘法次数
                }
            }
        }
    }
    return 0;
}

//根据记录的分割点，生成最后的矩阵相乘表达式
string make_result(vector<vector<int> > &points, int t1, int t2) {
    if (t1 == t2)
        return string(1, 'A' + t1 -1);
    int split = points[t1][t2];
    return "(" + make_result(points, t1, split) + "*" + make_result(points, split + 1, t2) + ")";
}


int main()
{
    int matrixNum = 0;///连乘的矩阵个数
    vector<pair<int, int>> data; ///存储矩阵的尺寸

    //手动输入
    //输入的第一个值是第一个矩阵的行，其余是各矩阵的列
    int* matrixSizeSeq;
    matrixSizeSeq = uncertainInput(matrixNum);
    if (matrixSizeSeq == NULL)
        return false;
    
    //重构为向量
    /*for (int ptrOffset = 0; ptrOffset < matrixNum - 1; ptrOffset++)    ///遍历到倒数第二个输入
        data.push_back(make_pair(*(matrixSizeSeq + ptrOffset), *(matrixSizeSeq + ptrOffset + 1)));*/

    //为记录乘法次数和分割点的n阶矩阵分配空间并初始化
    vector<vector<int> > num(matrixNum, vector<int>(matrixNum));
    vector<vector<int> > points(matrixNum, vector<int>(matrixNum));
    for (int i = 0; i < matrixNum; i++) {
        for (int j = 0; j < matrixNum; j++) {
            points[i][j] = -1;
            num[i][j] = 0;
        }
    }

    calculate_M1(matrixSizeSeq, matrixNum, num, points);
    cout << make_result(points, 1, matrixNum - 1) << "	最少乘法次数为:" << num[1][matrixNum - 1] << endl;
    cin >> matrixNum;
    return 0;
}

int* uncertainInput(int& inputNum)
{
    //输入0表示输入完成 输入负数认为手误
    using namespace std;

    int allocateSpace = 5;
    int* inputSequence = (int *)malloc(allocateSpace * sizeof(int));    ///存储矩阵的尺寸输入


    int ptrOffset = 0;        ///指针偏移
    int tempSize = 0;        ///临时存储输入值

    do {
        cin >> tempSize;
        if (tempSize <= 0)
            continue;
        *(inputSequence + ptrOffset++) = tempSize;
        if (ptrOffset >= allocateSpace)
        {
            allocateSpace += 5;
            inputSequence = (int*)realloc(inputSequence, allocateSpace * sizeof(int));
        }
    } while (tempSize != 0);
    inputNum = ptrOffset;    ///偏移量即输入个数
    return inputSequence;
}

比如求$S_{ad}$，则在控制台输入a b c d 0，然后回车即可。当然 a b c d 必须替换为数字。

代码中用到的一些知识

C++提供模版类string，其中一个构造方法可将字符转化为字符串。如 string(1, 'A'+1)，第一个参数是源字符延拓次数，这个构造函数将‘B’转化为"B"。