第3章动态规划——矩阵连乘最优计算方式查找

问题：如何得到n个矩阵连乘 $A_{p_1p_2}A_{p_2p_3}A_{p_3p_4}...A_{p_{n-1}p_n}A_{p_{n}p_{n+1}}$ 的最少计算次数的计算顺序？先计算 $A_{p_1p_2}A_{p_2p_3}$ ，还是先计算 $A_{p_2p_3}A_{p_3p_4}$ ？其中， $p_k$ 为矩阵的维度。

1、两个矩阵相乘 $A_{p_1p_2}A_{p_2p_3}$

两个矩阵相乘 $A_{p_1p_2}A_{p_2p_3}$ 需要多少次运算呢？需要做 $p_1\times p_2\times p_3$ 次乘法。具体来说，它将得到矩阵 $(p_1\times p_3)$ 大小的矩阵，该矩阵中每个元素由 $p_2$ 次乘法得到。

2、多个矩阵连乘

多个矩阵连乘会有多少种计算方式呢？

穷举法
$A_{p_1p_2}A_{p_2p_3}A_{p_3p_4}...A_{p_{n-1}p_n}A_{p_{n}p_{n+1}}$ ，假设从第k个矩阵断开，分为
$(A_{p_1p_2}A_{p_2p_3}...A_{p_kp_{k+1}})(A_{p_{k+1}p_{k+2}}...A_{p_{n-1}p_n}A_{p_{n}p_{n+1}})$ ，先计算前半部分，再计算后半部分，最后前后两个矩阵相乘。则前半部分有P(k)种计算方式，后半部分P(n-k)种，总共P(k)P(n-k)种计算方式。
k从1到n-1，则总共有P(n)种：
$P(n)=\sum_{k=1}^{n-1}{P(k)P(n-k)}, P(1)=1$
据说P(n)=C(n-1)= $\frac{1}{n} \left( \begin{matrix} 2(n-1) \\ (n-1) \end{matrix} \right) = \Omega \left(4^{n-1}/(n-1)^{3/2}\right)$
上述方法中，会涉及到递归，同时有很多重复运算。
动态规划就是把重复计算的结果保存下来，下次遇到时只需要查找到已经计算好的结果。即上述穷举过程中，从i到j个矩阵相乘 $A_{p_{i}p_{i+1}}A_{p_{i+1}p_{i+2}}...A_{p_{j}p_{j+1}}$ ，记为 $A[i:j]$ ，会出现很多次，记录下来后，就不用重复计算。去除重复计算后，整个连乘只需要计算 $n^2$ 次。
程序
以 $A_{30\times 35}A_{35\times 15}A_{15\times 5}A_{5\times 10}A_{10\times 20}A_{20\times 25}$ 矩阵连乘为例，优化后递归实现方式的Python程序如下：

# 输入 一组矩阵的shape，A[[p0, p1], [p1, p2], [p2, p3],  ..., [p_n-1, p_n]]
# 中间存储 A[i,:]到A[j,:]到最小计算次数m[i,j]、最优断开点s
# 输出  A[0,:]到A[n-1,:]到最小计算次数、最优断开点s

import numpy as np 
A = [[30,35], [35,15], [15,5], [5,10], [10,20], [20,25]]
n = np.shape(A)[0] #获得连乘矩阵的个数n
m = np.zeros((n, n, 2), dtype=int) #存储[i,j]从Ai到Aj的[最小计算次数, 最优断开点s]

# A里面相邻矩阵：前一个矩阵列数必须等于后一个矩阵的行数，否则返回错误
def isAillegal(A, i, j):
    # 略
    return False

# 主要函数matrixChain
# 输入  矩阵链A
#      i,j表示从Ai到Aj的 
# 返回  从Ai到Aj的[最小计算次数, 最优断开点s],并记录在m中
def matrixChain(A, i, j):
    # 输入合法性判断
    if i==j:
        return [0,0] #i==j时，m[i,j]的最小计算次数为0
    if i>j or j>=n or i<0:
        return [-1,-1] 
    if isAillegal(A, i, j):
        return [-1,-1]

    # 查找m是否已经有值
    if m[i,j,0]!=0 :
        return m[i,j]

    #相邻矩阵相乘，例如A(P1,P2)*A(P2,P3)需要P1*P2*P3次乘法操作
    if j==(i+1):
        m[i][j][0] = A[i][0]*A[i][1]*A[j][1]
        m[i][j][1] = i
        return m[i,j]

    # 矩阵序列从s处分成前后两部分，循环判断s的位置使得矩阵连乘计算次数最少
    minCount = 0
    minS = i
    for s in range(i,j):
        first = matrixChain(A,i,s)
        last = matrixChain(A,s+1,j)
        count = first[0] + last[0] + A[i][0]*A[s][1]*A[j][1]
        if (s==i) or (count<minCount):
            minCount = count
            minS = s
    m[i][j] = [minCount, minS]
    return m[i,j]

# 输出Ai到Aj的最优计算方式，类似(A0(A1A2))((A3A4)A5))
def traceS(mm, i, j):
    if i==j:
        return 'A'+str(i)  #只剩余一个矩阵，则返回‘Ai’
    if (i+1)==j:
        return 'A'+str(i)+'A'+str(j)+'' #剩余2个矩阵，则返回‘AiAj’
    
    # 否则，分成左右两个部分
    leftStr = '('+traceS(mm, i, mm[i,j,1])+')'
    RightStr = '('+traceS(mm, mm[i,j,1]+1, j)+')'
    return leftStr+RightStr
    
matrixChain(A, 0, n-1)
print(m[:,:,1])
print(traceS(m, 0, n-1))

结果如下：

最优计算次数：
[[    0 15750  7875  9375 11875 15125]
 [    0     0  2625  4375  7125 10500]
 [    0     0     0   750  2500  5375]
 [    0     0     0     0  1000  3500]
 [    0     0     0     0     0  5000]
 [    0     0     0     0     0     0]]
最优断开处：
[[0 0 0 2 2 2]
 [0 0 1 2 2 2]
 [0 0 0 2 2 2]
 [0 0 0 0 3 4]
 [0 0 0 0 0 4]
 [0 0 0 0 0 0]]
最优计算方式：
((A0)(A1A2))((A3A4)(A5))

直接动态优化算法的Python程序如下：

import numpy as np 
A = [[30,35], [35,15], [15,5], [5,10], [10,20], [20,25]]
n = np.shape(A)[0] #获得连乘矩阵的个数n
m = np.zeros((n, n, 2), dtype=int) #存储[i,j]从Ai到Aj的[最小计算次数, 最优断开点s]

def matrixChainDirect(A):
    for k in range(1, n):  #从2个矩阵连乘开始，计算到n个矩阵连乘。
        for i in range(0,n-k): #从第i个矩阵开始连乘k个矩阵
            minCount = 0
            minS = i
            for s in range(i,i+k): # A_i * A_i+1 *... *A_i+k-1 矩阵连乘从s处断开
                count = m[i][s][0] + m[s+1][i+k][0] + A[i][0]*A[s][1]*A[i+k][1]
                if (s==i) or (count<minCount):
                    minCount = count
                    minS = s
            m[i][i+k] = [minCount, minS]
    return m

matrixChainDirect(A)
print('最优计算次数：')
print(m[:,:,0])
print('最优断开处：')
print(m[:,:,1])
print('最优计算方式：')
print(traceS(m, 0, n-1))

结果同递归调用方式。

第3章 动态规划——矩阵连乘最优计算方式查找

1、两个矩阵相乘

2、多个矩阵连乘

推荐阅读更多精彩内容

第3章动态规划——矩阵连乘最优计算方式查找

1、两个矩阵相乘 $A_{p_1p_2}A_{p_2p_3}$