动态规划之矩阵链乘 C++实现

原理

在上一次的文章当中，作者讲解了什么是动态规划，以及动态规划的一个举例应用，这次，我们来看看如何运用动态规划来解决矩阵链乘问题。

关于矩阵的乘法，运用如下公式：

C=A×B

其中

cij=∑k=1naikbkj

并不是任意两个矩阵都能相乘，必须满足A的列数na等于B的行数mb，即有Anama,Bnbmb且na=mb。相乘之后，矩阵C的行列分别为ma、nb。

如果不只是两个矩阵相乘，而是n个矩阵相乘呢？同样可以用上面的公式算出结果，但是存在着一个乘积次数的问题，若有Apq,Bqr且na=mb，则其相乘次数为pqr。假设有三个矩阵⟨A1,A2,A3⟩，其规模分别为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×pi(1≤i≤n)，求完全括号化方案，使得计算乘积A1A2…An所需标量乘法次数最少。

令m[i,j]表示计算矩阵Aij所需标量乘法次数的最小值，那么A1,A2…An相乘，最小代价括号化方案的递归求解公式变为：

m[i,j]=⎧⎩⎨0maxi≤k<j{m[i,k]+m[k+1,j]+pi−1pkpj}(i=j)(i<j)

以下附上自底向上方法O(n3)，递归版Ω(2n)与带备忘的自顶向下法O(n3)的实现。

源代码：

#include <iostream>#include <utility>#include <vector>using namespace std;vector<int> temp_Vec{ 30,35,15,5,10,20,25 }; //A1(30X35), A2(35X15.).....A6(20X25)//Bottom-up.pair<vector<vector<int>>, vector<vector<int>>> Matrix_Chain_Order(const vector<int> &temp_Vec) {    auto temp_n = temp_Vec.size() - 1;    vector<vector<int>> temp_VecM, temp_VecS;    temp_VecM.resize(temp_n + 1);    temp_VecS.resize(temp_n + 1);    for(auto &i : temp_VecM) {        i.resize(temp_n + 1);    }    for (auto &i : temp_VecS) {        i.resize(temp_n + 1);    }    for(auto i = 1; i <= temp_n; ++i) {        temp_VecM[i][i] = 0;    }    for(auto l = 2; l <= temp_n; ++l) {        for (auto i = 1; i <= temp_n - l + 1; ++i) {            auto j = i + l - 1;            temp_VecM[i][j] = INT_MAX;            for(auto k = i; k <= j - 1; ++k) {                auto temp_q = temp_VecM[i][k] + temp_VecM[k + 1][j] + temp_Vec[i - 1] * temp_Vec[k] * temp_Vec[j];                if(temp_q < temp_VecM[i][j]) {                    temp_VecM[i][j] = temp_q;                    temp_VecS[i][j] = k;                }            }        }    }    return make_pair(temp_VecM, temp_VecS);}//Memoized of part1.int Lookup_Chain(vector<vector<int>> &temp_VecM, const vector<int> &temp_Vec, const int &i, const int &j) {    if(temp_VecM[i][j] < INT_MAX) {        return temp_VecM[i][j];    }    if(i == j) {        temp_VecM[i][j] = 0;    }    else {        for(auto k = i; k <= j - 1; ++k) {            auto temp_q = Lookup_Chain(temp_VecM, temp_Vec, i, k) + Lookup_Chain(temp_VecM, temp_Vec, k + 1, j) + temp_Vec[i - 1] * temp_Vec[k] * temp_Vec[j];            if(temp_q < temp_VecM[i][j]) {                temp_VecM[i][j] = temp_q;            }        }    }    return temp_VecM[i][j];}//Memoized of part2.int Memoized_Matrix_Chain(const vector<int> &temp_Vec) {    auto temp_n = temp_Vec.size() - 1;    vector<vector<int>> temp_VecM;    temp_VecM.resize(temp_n + 1);    for (auto &i : temp_VecM) {        i.resize(temp_n + 1);    }    for(auto &i : temp_VecM) {        for(auto &j : i) {            j = INT_MAX;        }    }    return Lookup_Chain(temp_VecM, temp_Vec, 1, temp_n);}//Recusive.int Recursive_Matrix_Chain(const vector<int> &temp_Vec, const int &i, const int &j) {    if (i == j) {        return 0;    }    auto temp_n = temp_Vec.size() - 1;    vector<vector<int>> temp_VecM;    temp_VecM.resize(temp_n + 1);    for (auto &x : temp_VecM) {        x.resize(temp_n + 1);    }    for (auto &x : temp_VecM) {        for (auto &y : x) {            y = INT_MAX;        }    }    for (auto k = i; k <= j - 1; ++k) {        auto temp_q = Recursive_Matrix_Chain(temp_Vec, i, k) + Recursive_Matrix_Chain(temp_Vec, k + 1, j) + temp_Vec[i - 1] * temp_Vec[k] * temp_Vec[j];        if (temp_q < temp_VecM[i][j]) {            temp_VecM[i][j] = temp_q;        }    }    return temp_VecM[i][j];}//Print optimal bracket scheme.void Print_Optimal_Parens(const vector<vector<int>> &temp_VecS, const int &i, const int &j) {    if(i == j) {        cout << "A" << i << " ";    }    else {        cout << "(" << " ";        Print_Optimal_Parens(temp_VecS, i, temp_VecS[i][j]);        Print_Optimal_Parens(temp_VecS, temp_VecS[i][j] + 1, j);        cout << ")" << " ";    }}int main() {    auto temp_Pair = Matrix_Chain_Order(temp_Vec);    auto temp_VecM = temp_Pair.first, temp_VecS = temp_Pair.second;    cout << temp_VecM[2][5] << endl;  //A2 to A5 by bottom-up.    Print_Optimal_Parens(temp_VecS, 2, 5);  //A2 to A5 by bottom-up.    cout << endl << endl;    for(auto &i : temp_VecM) {  //All of the value.        for(auto &j : i) {            cout << j << " ";        }        cout << endl;    }    cout << endl << Recursive_Matrix_Chain(temp_Vec, 1, 6) << endl; // A1 to A6 by recursive.    cout << endl << Memoized_Matrix_Chain(temp_Vec) << endl;  //A1 to A6 by memoized.    return 0;}