矩阵链乘法

ResysChina 2018-02-17

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

我们知道矩阵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个矩阵的链<A1,A2,...,An>,矩阵Ai的规模为pi - 1×pi(1 ≤ i ≤ n),求完全括号化方案,使得计算乘积A1A2...An所需标量乘法次数最少。

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

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

对于i = j 时,矩阵链只包含唯一的矩阵 Ai...i = Ai,不需要做任何标量乘法运算。若 i < j,利用最优子结构来计算m[i][j]。如前面假设Ai,Ai + 1,...,Aj的最优括号化方案的分割点在矩阵的Ak和Ak + 1之间,其中i ≤ k < j。那么,m[i][j]就等于计算Ai...k和Ak + 1...j的代价加上两者相乘的代价的最小值。由于矩阵Ai的大小为pi×pi + 1,易知Ak + 1...j相乘的代价为pi·pk + 1·pj + 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;
}

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

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

相关推荐