矩阵链相乘顺序

• 矩阵链相乘顺序
  • 问题描述
    • 举例说明
  • 问题分析
  • 程序

问题描述

A1,A2,..,An 表示n个矩阵的序列，其中Ai为\(P_{i−1}×P_i\)阶矩阵，i=1,2,...,n。
向量P=<P0,P1,P2..Pi>表示矩阵链的输入，其中P0是A1的行数，P1是A1的列数，P1是A2的行数，以此类推。
计算这个矩阵需要做n−1次两个矩阵的相乘运算，可以用n−1对括号表示运算次序。
因为矩阵乘法满足结合律，所以无论采用那种顺序，最后结果都一样，但是采用不同的顺序计算的工作量不同。如何定义两个矩阵相乘的工作量呢？
所以假设A1有i行k列，A2有k行j列。所以A1A2相乘后的矩阵有i行j列，含ij个元素。
以元素相乘作为基本运算，乘积中每个元素的计算都需要做j次乘法，于是计算A1A2总共需要ijk次乘法。

举例说明

假设输入的是P=<10,100,5,50>,说明有3个矩阵相乘。其中，
A1:10×100
A2:100×50
A3:5×50

有两种乘法次序：
(A1A2)A3
A1(A2A3)

执行第一种运算的基本运算次序：10×100×5+10×5×50=7500
执行第二种运算的基本运算次序：10×100×50+100×5×50=75000
工作量相差达10倍！
所以我们的问题是：给定向量P，确定一种乘法次序，使得基本运算的总次数最少。

问题分析

1、寻找最优子结构：
假设我们已经找到父矩阵链最优解，当我们划分到最后一步时都是两个子矩阵链（分别被括号包围）相乘，如（A1A2A3A4）（A5A6A7），此时显然外括号为父矩阵的最优括号划分。继续往下划分，（（A1A2A3）A4）（A5（A6A7）），则两个子矩阵链的划分括号必须也为他们本身的最优括号划分。因为如果子矩阵链的划分括号不是他们本身的最优括号划分时，两个子矩阵链就有另一种最优括号划分，如（A1（A2A3A4））（A5（A6A7））。那么就与我们的假设相悖。（（A1A2A3）A4）（A5（A6A7））不是我们父矩阵最优解。所以拥有最优划分的父矩阵的子矩阵链显然也拥有最优划分。

2、子问题重叠：
在分解子结构的过程中，会多次访问同一个子问题，比如 (A1(A2A3A4))(A5(A6A7)) 和 (A1(A2A3A4)A5)(A6A7)。
这正是动态规划解决的问题。

程序

#include <iostream>
#include <vector>
#include <new>
#include <climits>

using namespace std;

int solve(const vector<pair<int, int>>& mat_chain);
void print_mat(int** mat, int len1, int len2);


int main (int argc, char** argv) {
    vector<pair<int, int>> mat_chain{
        {10,100},
        {100,5},
        {5,50},
        {50,9},
    };

    int ret = solve(mat_chain);
    cout << "ret: " << ret << endl;

    return 0;
}

int solve(const vector<pair<int, int>>& mat_chain) {
    int len = mat_chain.size();
    int i, j, k, diff, sum;
    int** order = new int*[len];
    // mount 右上存放计算量：[i,j]矩阵链的最小计算量为mount[i,j]
    // mount 左下存放计算顺序：mount[i][j]=k 表示[i,j]矩阵链分割为[i,k]和[k+1,j]两部分
    //   先分别计算这两部分，再对其相乘
    int** mount = new int*[len];
    for (i = 0; i < len; ++i) {
        mount[i] = new int[len];
        for (j = 0; j < len; ++j) {
            mount[i][j] = INT_MAX;
        }
        // 对角线表示自己和自己相乘，计算量为0
        mount[i][i] = 0;
    }

    for (diff = 1; diff < len; ++diff) {
        for (i = 0; i < len - diff; ++i) {
            j = i + diff;
            for (k = i; k < j; ++k) {
                // 计算量=分割的两部分各自的计算量加上两部分相乘的计算量
                sum = mount[i][k] + mount[k+1][j]
                    + mat_chain[i].first * mat_chain[k].second * mat_chain[j].second;
                if (sum < mount[i][j]) {
                    mount[i][j] = sum;
                    mount[j][i] = k;
                }
            }
        }
    }
    cout << "mount:\n";
    print_mat(mount, len, len);
    sum = mount[0][len-1];

    for (i = 0; i < len; ++i) {
        delete[] mount[i];
    }
    delete[] mount;

    return sum;
}

void print_mat(int** mat, int len1, int len2) {
    for (int i = 0; i < len1; ++i) {
        for (int j = 0; j < len2; ++j) {
            cout << mat[i][j] << "\t";
        }
        cout << endl;
    }
}

测试：

mount:
0	5000	7500	7700
0	0	25000	6750
1	1	0	2250
1	1	2	0
ret: 7700