MCM（矩阵链乘法）

Posted 2020-10-24 darkchii

　　这是《算法导论》动态规划中的一个问题。问题简述如下：我们在求解矩阵相乘时通常会有一个最优括号方案来对矩阵进行顺序相乘，这样会减少大量的计算时间。

　　我们知道矩阵A*B相乘，只能是当矩阵A的列数等于矩阵B的行数时才能进行相乘，且假设A·B = C，A为p×q规模的矩阵，B为q×r的矩阵，则乘积C的规模为p×r的矩阵。计算A·B所需的时间由乘法次数决定，即pqr。

　　例如：有三个矩阵的规模分别为：A = 10×100，B = 100×5，C = 5×50。如果按顺序(AB)C计算，则需要10×100×5 + 10×5×50 = 7500次乘法，如果按顺序A(BC)计算，则需要100×5×50 + 10×100×50 = 75000次乘法...所以，按第一种顺序计算矩阵链乘积要比第二种快非常多。

　　矩阵链乘法问题可描述如下：给定n个矩阵的链<A₁,A₂,...,A_n>，矩阵A_i的规模为p_{i - 1}×p_i(1 ≤ i ≤ n)，求完全括号化方案，使得计算乘积A₁A₂...A_n所需标量乘法次数最少。

　　分析问题描述可知道，我们假设在矩阵链乘中找到一个合适的切割点进行括号化，可得到问题的最优括号方案，然后对切割点的两边的子问题用相同的方法独立求解，不断进行下去，则最终将得到原问题的最优解。这就是该问题的最优子结构。

　　构造递归公式：令m[i][j] 表示计算矩阵A_i...j所需标量乘法次数的最小值，那么，原问题的最优解 —— 计算A_i...n所需的最低代价就是m[1][n]。

　　对于i = j 时，矩阵链只包含唯一的矩阵 A_i...i = A_i ，不需要做任何标量乘法运算。若 i < j，利用最优子结构来计算m[i][j]。如前面假设A_i,A_{i + 1},...,A_j的最优括号化方案的分割点在矩阵的A_k和A_{k + 1}之间，其中i ≤ k < j。那么，m[i][j]就等于计算A_i...k和A_{k + 1...j}的代价加上两者相乘的代价的最小值。由于矩阵A_i的大小为p_i×p_{i + 1}，易知A_{k + 1...j}相乘的代价为p_i·p_{k + 1}·p_{j + 1}次标量乘法运算。因此，得到：

m[i][j] = m[i][k] + m[k + 1][j] + p[i]·p[k + 1]·p[j + 1]

　　由于k值取法只有 j - i 种可能的取值，即k = i, i + 1, ..., j。检查所有可能的情况，找到最优，因此递归公式变为：

           (    0,                 i = j;
m[i][j] = {    
           (    min(m[i][k] + m[k + 1][j] + p[i]·p[k + 1]·p[j + 1]),    i < j.

　　O(n^3)迭代解法：

#include <iostream>
#include <vector>

class Solution {
public:
    std::vector<std::vector<int> > MatrixChainOrder(std::vector<int>& p)
    {
        const int n = p.size() - 1;
        std::vector<std::vector<int> > m(n, std::vector<int> (n));

        for(int i = 0; i < n; i++)
        {
            m[i][i] = 0;
        }
        for(int l = 1; l < n; l++)
        {
            for(int i = 0; i < n - l; i++)
            {
                int j = i + l;
                m[i][j] = INT_MAX - 1;
                for(int k = i; k < j; k++)
                {
                    int q = m[i][k] + m[k + 1][j] + p[i]*p[k + 1]*p[j + 1];
                    if(q < m[i][j])
                    {
                        m[i][j] = q;
                    }
                }
            }
        }

        return m;
    }
};

int main()
{
    std::vector<int> p{30, 35, 15, 5, 10, 20, 25};
    Solution solve;
    std::vector<std::vector<int> > res = solve.MatrixChainOrder(p);

    for(int i = res.size() - 1; i >= 0; i--)
    {
        for(int j = 0; j <= i; j++)
            std::cout << res[j][i] << ‘ ‘;
        std::cout << std::endl;
    }

    return 0;
}

　　递归版感觉有点问题...所以再研究一段时间。

　　参考资料：以上皆为《算法导论》 — 动态规划原文