ACM基础之动态规划DP:矩阵链乘法Matrix-chain multiplication


一、思路

1. 数学原理

矩阵相乘时不同的结合顺序,虽然结果都一样,但会导致不同的cost(标量乘法代价)。

在这里插入图片描述

可以从其得出规律:

  • 不用考虑A1,A2,…,An的Ai的前后顺序。意思是说不用考虑这种情况A1 x A3 x A2(违背矩阵乘法的基本原则)。
  • 只用考虑先结合谁,再结合谁的结合顺序。而且,大的结合顺序可以划分成两个更小的结合顺序,如A1 x A3划分成了A1A2 x A3这两个子结合顺序。所以,整体的代价就是两个子结合顺序内部的代价,和将两个子部分再结合的代价。如,A1 x (A2 x A3)分别是0,1000,200

那么如何知道该怎么结合,答案很简单,用for循环暴力试。

2. 核心公式

在这里插入图片描述
符号:

  • A i A_i Ai:第i个矩阵,从 A 1 A_1 A1开始。

  • m[i,j]:表示子结合顺序内部的即 A i × A i + 1 . . . × A j A_i\times A_{i+1} ... \times A_j Ai×Ai+1...×Aj的代价,而且是 最优 代价。其中,m[i,i]表示 A i × A i A_i \times A_i Ai×Ai,代价为0。
    如:A1 x (A2 x A3)的最优总代价就是m[1,3]。

  • k:表示划分,左边的子结合顺序是到Ak,右边的子结合顺序从Ak+1开始。
    i ≤ k < j i \leq k < j ik<j这个等于不等于是根据 m i n . . . min{...} min...那个公式来划分的。

    • k ≥ i k \geq i ki是因为当k=i时,因为可以有m[i,i]。
    • k < j k < j k<j是因为m[k+1,j],k+1必须小于等于j。
  • p:表示矩阵的维度描述,用n+1个数描述n个矩阵,用来直接计算代价。如:A1 , A2 , A3,就是p=<10,20,50,1>
    从0开始A1 x A2p0 x p1 x p2

  • p i − 1 p k p j p_{i-1}p_{k}p_{j} pi1pkpj:表示将两个子部分 ( A i × . . . × A k ) × ( A k + 1 × . . . × A j ) (A_i \times ... \times A_k)\times(A_{k+1} \times ... \times A_j) (Ai×...×Ak)×(Ak+1×...×Aj)再结合的代价。
    p 0 p 1 p 3 p_0p_1p_3 p0p1p3就是 ( A 1 ) ( A 2 A 3 ) (A_1)(A_2A_3) (A1)(A2A3) p 0 p 2 p 3 p_0p_2p_3 p0p2p3就是 ( A 1 A 2 ) ( A 3 ) (A_1A_2)(A_3) (A1A2)(A3)

3. 图示

【两个右上的三角表格】:

  • m表:表示代价。n行n列。
  • s表:表示k。因为m[i,i]不需要也不能划分,所以去掉最外面这条对角线,n-1行n-1列。

在这里插入图片描述
在这里插入图片描述

使用动态程序填写表m[i,j]

  • 首先将m[i,i]=0,其中i=1,…,n。(对应单个的矩阵代价,如例A1 x (A2 x A3)划分的A1代价)
  • 然后计算m[1,2]m[2,3],…,m[n-1,n]。(对应最小结合代价,两个矩阵的结合代价)
  • 然后m[1,3]m[2,4],…,m[n-2,n],…(对应三个矩阵的结合总代价)
  • …一直到我们可以计算m[1,n]。(这就是从下到上的体现)

PS1:计算m[i,j]的k展开形状:
在这里插入图片描述

PS:如何从S推导出矩阵的结合顺序?
要划分Ai到Aj,就去查s[i,j],根据其值k划分为(Ai...Ak)(Ak+1...Aj)
在这里插入图片描述

例题:
在这里插入图片描述
在这里插入图片描述

p=<30,35,15,5,10,20,25>

  • m [ 1 , 2 ] = m [ 1 , 1 ] + m [ 2 , 2 ] + p 0 p 1 p 2 = 0 + 0 + ( A 1 ) ( A 2 ) = 30 × 35 × 15 = 15750 m[1,2]=m[1,1]+m[2,2]+p_0p_1p_2=0+0+(A_1)(A_2)=30\times35\times15=15750 m[1,2]=m[1,1]+m[2,2]+p0p1p2=0+0+(A1)(A2)=30×35×15=15750
    k=1,s[1,2]=1。
  • m [ 1 , 3 ] m[1,3] m[1,3]
    • k=1时,即 ( A 1 ) ( A 2 A 3 ) (A_1)(A_2A_3) (A1)(A2A3) m [ 1 , 3 ] = m [ 1 , 1 ] + m [ 2 , 3 ] + p 0 p 1 p 3 = 0 + 2625 + 30 × 35 × 5 = 7875 m[1,3]=m[1,1]+m[2,3]+p_0p_1p_3=0+2625+30\times35\times5=7875 m[1,3]=m[1,1]+m[2,3]+p0p1p3=0+2625+30×35×5=7875
    • k=2时,即 ( A 1 A 2 ) ( A 3 ) (A_1A_2)(A_3) (A1A2)(A3) m [ 1 , 3 ] = m [ 1 , 2 ] + m [ 3 , 3 ] + p 0 p 2 p 3 = 15750 + 0 + 30 × 15 × 5 = 18000 m[1,3]=m[1,2]+m[3,3]+p_0p_2p_3=15750+0+30\times15\times5=18000 m[1,3]=m[1,2]+m[3,3]+p0p2p3=15750+0+30×15×5=18000
      所以,选择k=1,s[1.3]=1。

4. 伪代码

在这里插入图片描述

// p比如<30,35,15>,表示A1(30*35)和A2(35*15)
MATRIX-CHAIN-ORDER(p):
	// p用n+1个数表示n个矩阵相乘
	n ← length[p] - 1
	// 对表中第一条对角线初始化当i=j的情况
	for i ← 1 to n:
		do	m[i, i]0
	// 控制每条红色对角线的行号区间递减,每条红色对角线列号在随行号递增的同时随对角线向上加一个数递增
	for l ← 2 to n:
			// 每条红色对角线的行号区间,[1,n-1],[1,n-2],...,[1,1]
		do	for i ← 1 to n - l + 1:
					// 每条红色对角线的每行对应的列号
				do	j ← i + l - 1
					// 先初始化正无穷
					m[i, j] ← ∞
					// 遍历k,i≤k<j,即[i,j-1]:找出最小的m[i,j]
					for k ← i to j - 1:
							// 公式
						do	q ← m[i, k] + m[k + 1, j] + (pi-1)(pk)(pj)
							// 比较新试的k和原来哪个更小
							if	q < m[i, j]
							then	m[i, j] ← q
									s[i, j] ← k 	

二、cpp实现

#include <iostream>
#include <vector>
using namespace std;

#define INF 0x3f3f3f3f

void matrix_chain_order(vector<int> &p, vector<vector<int>> &m, vector<vector<int>> &s)
{
    // p用n+1个数表示n个矩阵相乘
    int n = p.size() - 1;
    // 对表中第一条对角线初始化当i=j的情况
    // 表行列总数都是n个矩阵个,[1,n]
    for (int i = 1; i <= n; i++)
    {
        m[i][i] = 0;
    }
    // 控制每条红色对角线的行号区间递减,每条红色对角线列号在随行号递增的同时随对角线向上加一个数递增
    for (int l = 2; l <= n; l++)
    {
        // 每条红色对角线的行号区间,[1,n-1],[1,n-2],...,[1,1]
        for (int i = 1; i <= n - l + 1; i++)
        {
            // 每条红色对角线的每行对应的列号
            int j = i + l - 1;
            // 先初始化正无穷
            m[i][j] = INF;
            // 遍历k,i≤k<j,即[i,j-1]:找出最小的m[i,j]
            for (int k = i; k <= j - 1; k++)
            {
                // 公式,表示新试k的代价
                int q = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
                // 比较新试的k和原来哪个更小
                if (q < m[i][j])
                {
                    m[i][j] = q;
                    s[i][j] = k;
                }
            }
        }
    }
}

// void print_optimal_parens(vector<vector<int>> &s, int i, int j)
// {
//     if (i == j)
//     {
//         printf("A%d\n", i);
//     }
//     else
//     {
//         printf("(");
//         print_optimal_parens(s, i, s[i, j]);
//         print_optimal_parens(s, s[i, j] + 1, j);
//         printf(")\n");
//     }
// }

int main(void)
{
    // 矩阵链,7个数表示6个矩阵
    vector<int> p{30, 35, 15, 5, 10, 20, 25};

    // 这里我们不从0开始,从1开始使用,[1,6]范围
    // m,表示记录结果
    vector<vector<int>> m(p.size(), vector<int>(p.size()));
    // s,表示选择的k
    vector<vector<int>> s(p.size(), vector<int>(p.size()));

    // 调用
    matrix_chain_order(p, m, s);

    // 输出m
    for (int i = 1; i < p.size(); i++)
    {
        for (int j = 1; j < p.size(); j++)
        {
            printf("%5d ", m[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    // 输出s
    for (int i = 1; i < p.size(); i++)
    {
        for (int j = 1; j < p.size(); j++)
        {
            printf("%5d ", s[i][j]);
        }
        printf("\n");
    }
    return 0;
}

Reference

(最大矩阵链乘)Matrix-chain product

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值