演算法題系列：矩陣鏈乘問題

矩陣鏈乘問題是最典型的動態規劃問題，要理解下面的內容請先閱讀這篇動態規劃的總結。

1.問題描述

矩陣鏈乘問題的描述如下，就是說要確定一個完全加括號的形式使得矩陣鏈乘需要進行的標量計算數目最少，矩陣Ai的維數為pi−1×pi，如果窮舉所有可能形式的話，時間複雜度是指數級的！因為該問題滿足最優子結構，並且子問題存在重疊，所以我們可以藉助動態規劃來求解。

2.問題分析

我們需要確定一個遞迴式來將我們要求解的問題表示出來，下面摘自演算法導論，介紹地非常詳細

最後給出的遞迴式如下，就是說我們要如何確定從第i個矩陣到第j個矩陣組成的矩陣鏈的最優解。如果i和j相等，那麼就是一個矩陣，不需要運算；如果i小於j，那麼肯定要從它們中間的某個位置分開來，那從哪裡分開來呢? 這個我們可以嘗試下所有可能的選擇，也就是嘗試不同的位置k，k滿足條件(i <= k < j)，在位置k將矩陣鏈進行分開，看看它需要的計算次數，然後我們從這些可能的k中選擇使得計算次數最小的那個k進行分開，分開了之後我們的問題就變成了2個小問題，確定矩陣鏈從i到k 和另一個矩陣鏈從k+1到j的最優解。如果我們一開始設定i=1(第一個矩陣)，j=n(最後一個矩陣)，那麼，經過上面的遞迴即可得到我們需要的解。這就是遞迴的思想！

3.程式碼實現

根據上面的思想我們很快就可以寫出一個遞迴版本的矩陣鏈承法的實現程式碼，輸出的結果也沒有錯，給出的加括號的方式是( ( A1 ( A2 A3 ) ) ( ( A4 A5 ) A6 ) )。[問題的資料是演算法導論中的問題的資料，值是30,35,15,5,10,20,25]。

def matrixchain_rec(p,i,j):
    if i==j:
        return 0
    for k in range(i,j):
        q=matrixchain_rec(p,i,k)+matrixchain_rec(p,k+1,j)+p[i-1]*p[k]*p[j]
        if q<m[i][j]:
            m[i][j]=q
            s[i][j]=k
    return m[i][j]

def showmatrixchain(s,i,j):
    if i==j:
        print 'A%d'%(i),
    else:
        print '(',
        showmatrixchain(s,i,s[i][j])
        showmatrixchain(s,s[i][j]+1,j)
        print ')',

n=6
p=[30,35,15,5,10,20,25]
m=[[sys.maxint for i in range(n+1)] for j in range(n+1)]
s=[[0 for i in range(n+1)] for j in range(n+1)]
# pprint.pprint(m)
result=matrixchain_rec(p,1,6)
print(result) #15125
showmatrixchain(s,1,6) #( ( A1 ( A2 A3 ) ) ( ( A4 A5 ) A6 ) )

上面的程式碼執行沒有問題，但是，它不夠完美！為什麼呢? 很明顯，矩陣鏈乘問題子問題存在重疊，下面這張圖很形象地顯示了哪些子問題被重複計算了，所以我們需要改進，改進的方法就是使用帶備忘錄的遞迴形式！

要改成帶備忘錄的很簡單，但是，這次我們不能直接使用原來的裝飾器，因為Python中的dict不能對list物件進行hash，所以我們要簡單地修改下我們key值的構建，也很簡單，看下程式碼就明白了：

from functools import wraps

def memo(func):
    cache={}
    @wraps(func)
    def wrap(*args):
        #build new key!!!
        key=str(args[1])+str(args[2])
        if key not in cache:
            cache[key]=func(*args)
        return cache[key]
    return wrap

@memo
def matrixchain_rec(p,i,j):
    if i==j:
        return 0
    for k in range(i,j):
        q=matrixchain_rec(p,i,k)+matrixchain_rec(p,k+1,j)+p[i-1]*p[k]*p[j]
        if q<m[i][j]:
            m[i][j]=q
            s[i][j]=k
    return m[i][j]

def showmatrixchain(s,i,j):
    if i==j:
        print 'A%d'%(i),
    else:
        print '(',
        showmatrixchain(s,i,s[i][j])
        showmatrixchain(s,s[i][j]+1,j)
        print ')',

n=6
p=[30,35,15,5,10,20,25]
m=[[sys.maxint for i in range(n+1)] for j in range(n+1)]
s=[[0 for i in range(n+1)] for j in range(n+1)]
# pprint.pprint(m)
result=matrixchain_rec(p,1,6)
print(result) #15125
showmatrixchain(s,1,6) #( ( A1 ( A2 A3 ) ) ( ( A4 A5 ) A6 ) )

接下來的一個問題是，我們怎麼實現迭代版本呢? 迭代版本關鍵在於順序！我們怎麼保證我們在計算$A{i…j}的最優解時，所有可能的k的選擇需要求解的子問題A{i…k}以及A_{(k+1)…j}$是已經求解出來了的呢? 一個簡單但是有效的想法就是看矩陣鏈的長度，我們先計算矩陣鏈短的最優解，然後再計算矩陣鏈長的最優解，後者計算時所需要求解的子問題肯定已經求解完了，對不對? 於是就有了迭代版本的實現，需要注意的就是其中的i,j,k的取值範圍。

import sys
def matrixchain_iter(p):
    n=len(p)-1 #total n matrices 6
    #to solve the problem below, so initialize to n+1!!!
    m=[[0 for i in range(n+1)] for j in range(n+1)]
    s=[[0 for i in range(n+1)] for j in range(n+1)]
    # for i in range(n): #for matrix with len=1
        # m[i][i]=0
    # pprint.pprint(m)
    for l in range(2,n+1): #iterate the length, max is n
        for i in range(1,n-l+2): #i max is n-l+1
            j=i+l-1 #j is always l away from i
            m[i][j]=sys.maxint #initial to infinity
            for k in range(i,j):
                #attention to python array when index < 0!!!
                #solution is using more space with useless values
                q=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j]
                if q<m[i][j]:
                    m[i][j]=q
                    s[i][j]=k
        # print('when len is %d ' % (l))
        # pprint.pprint(m)
    return m,s

print('')
m,s=matrixchain_iter(p)
print(m[1][6]) #15125
showmatrixchain(s,1,6) #( ( A1 ( A2 A3 ) ) ( ( A4 A5 ) A6 ) )

實現的時候需要注意一點，在Python中取list中的值時，如果索引是負值的話會從後面往前數返回對應的元素，而以前我們用其他語言的時候肯定是提示越界了，所以程式碼中用來儲存結果的數陣列是(n+1)x(n+1)，而不是nxn的，這樣的話就能夠保證返回的是0，而不是從後往前數得到的結果。

得到的陣列m如下，m[1,6]就是我們需要的解。

[[0, 0, 0, 0, 0, 0, 0],
 [0, 0, 15750, 7875, 9375, 11875, 15125],
 [0, 0, 0, 2625, 4375, 7125, 10500],
 [0, 0, 0, 0, 750, 2500, 5375],
 [0, 0, 0, 0, 0, 1000, 3500],
 [0, 0, 0, 0, 0, 0, 5000],
 [0, 0, 0, 0, 0, 0, 0]]

陣列s如下：

[[0, 0, 0, 0, 0, 0, 0],
 [0, 0, 1, 1, 3, 3, 3],
 [0, 0, 0, 2, 3, 3, 3],
 [0, 0, 0, 0, 3, 3, 3],
 [0, 0, 0, 0, 0, 4, 5],
 [0, 0, 0, 0, 0, 0, 5],
 [0, 0, 0, 0, 0, 0, 0]]

將這個兩個陣列旋轉下，並且只看上三角部分的數字，就可以得到演算法導論中給出的那張三角圖形了，非常類似楊輝三角