矩陣乘法的順序安排問題 Python簡單實現

矩陣乘法的順序安排問題

問題背景

設矩陣 A、B 大小分別 \(p\times q\) , \(q \times r\) ，則矩陣乘積 AB 需要做的標量乘法次數為 \(p\times q \times r\) 。我們知道矩陣的乘法運算是不可交換的，但它是可結合的。因此對於多個矩陣的連乘，我們可以以任意順序添加括號改變其中相鄰矩陣乘法的優先級。不同計算順序下總的標量乘法運算次數是不同的，我們的目標是找到一個最優的矩陣乘法計算順序。

給定矩陣乘法序列 \(A_1, A_2, ..., A_n\)，將乘法序列以第 \(i\) 個矩陣分為前后兩部分，則方案數為前后兩部分方案數之積。因此乘法計算的順序個數為

\[T(n) = \Sigma_{i=1}^{n-1} T(i) \cdot T(n-i) \]

T(n) 的解為 Catalan數，這里不加證明給出結果為

\[T(n) = \text{C}_{2n}^n - \text{C}_{2n}^{n+1} = \frac {\text{C}_{2n}^n} {n+1} \]

由此可見的矩陣乘法順序個數為問題規模 \(n\) 的指數級，顯然通過枚舉找到最優的乘法順序是不合適的。

 

暴力算法

首先還是試探一下如何用最朴素的方式解決。

設 \(M_{i, j}\) 表示 第 \(i\) 個矩陣到第 \(j\) 個矩陣的最少乘法運算次數，用數學化的語言表達我們的目標，即

\[M_{1, n} = \min_i \{ M_{1, i} + M_{i+1, n} + p \times q \times r \} \]

其中 p、q、r為最后兩個矩陣的大小。

代碼很容易實現：

def minMatrixMultiplication(Mats):
	"""
	:param Mats: Mat類型的list
	:return:	 矩陣乘法的最小乘法次數，及對應的括號位置
	"""
    
	
	if len(Mats)==1:
		return 0, '[%d,%d]' % (Mats[0].n, Mats[0].m)

	import math
	minCost = math.inf
	bestSeq = ''        # 記錄添加的括號位置
	for i in range(0, len(Mats)-1):
		leftCost, leftSeq = minMatrixMultiplication(Mats[:i+1])
		rightCost, rightSeq = minMatrixMultiplication(Mats[i+1:])
		tmpCost = leftCost + rightCost +  Mats[0].n * Mats[i].m * Mats[-1].m

		if tmpCost < minCost:
			minCost = tmpCost
			bestSeq = '(' + leftSeq + '*' + rightSeq + ')'

	return minCost, bestSeq

測試用的矩陣類型Mat定義如下：

class Mat:
	def __init__(self, mat=None):
		if mat and isinstance(mat[0], list):
			self.mat = mat
			self.n = len(mat)
			self.m = len(mat[0])
		else:
			self.mat = [[]]
			self.n = 0
			self.m = 0

	def __init__(self, n, m):
		self.mat = [[]]
		self.n = n
		self.m = m

以上算法的函數調用次數 \(f(n) = 1 + f(1)+f(n-1) + f(2)+f(n-2) + ... + f(n-1)+f(1)\)，

容易驗證得到\(f(n) = 3^n\), 即該算法的復雜度為 O(\(3^n\))，這是不可接受的。

 

記憶化

分析一番可以發現，對於矩陣序列 i~j 之間乘法的最優結果 \(M_{i, j}\) 只有 \(\text{C}_n^2\) 種，那么上述代碼的中間很多段都進行了重復計算。如果把中間得到的答案記錄下來，可以大大減少計算量。

在不改變上述算法的框架下，將 i~j 之間的結果 \(M_{i, j}\) 定義Python嵌套的內部函數。新增了變量 invokeCnt 統計遞歸函數需要重新計算 \(M_{i, j}\) 的次數。

def minMatrixMultiplication2(Mats):
	siz = len(Mats) + 1
	# 血的教訓：不要使用下面的方法定義二維數組
	# minCostMem = [[-1]*siz]*siz
	# bestSeqMem = [['']*siz]*siz
	minCostMem = [[-1]*siz for i in range(siz)]
	bestSeqMem = [['']*siz for i in range(siz)]

	invokeCnt = 0  # 統計遞歸函數重新執行次數
	def helper(s, t):
		if s==t:
			return 0, '[%d,%d]' % (Mats[s].n, Mats[s].m)

		if minCostMem[s][t]!=-1:
			return minCostMem[s][t], bestSeqMem[s][t]

		nonlocal invokeCnt
		invokeCnt += 1

		import math
		minCost = math.inf
		bestSeq = ''
		for i in range(s, t):
			leftCost, leftSeq = helper(s, i)
			rightCost, rightSeq = helper(i+1, t)
			tmpCost = leftCost + rightCost +  Mats[s].n * Mats[i].m * Mats[t].m
			if tmpCost < minCost:
				minCost = tmpCost
				bestSeq = '(' + leftSeq + '*' + rightSeq + ')'

		minCostMem[s][t] = minCost
		bestSeqMem[s][t] = bestSeq

		return minCost, bestSeq

	return helper(0, len(Mats)-1), invokeCnt

 

動態規划

(待補充。。。）

 

運行對比

Mats = [Mat(2,3), Mat(3,5), Mat(5,8), Mat(8,2), Mat(2,3), Mat(3,2), Mat(2,5), Mat(5, 3)]
print(minMatrixMultiplication(Mats))
# (184, '((([2,3]*([3,5]*([5,8]*[8,2])))*([2,3]*[3,2]))*([2,5]*[5,3]))')
# 調用次數 3^8 = 2187
print(minMatrixMultiplication2(Mats))
# ((184, '((([2,3]*([3,5]*([5,8]*[8,2])))*([2,3]*[3,2]))*([2,5]*[5,3]))'), 28)

 

注意事項

Python 定義二維矩陣，千萬不要使用注釋寫法。調試了很久才發現問題。 T^T
正確的寫法為

• matrix = [[0]*m for i in range(n)]
• 或使用numpy庫
  import numpy
  matrix = numpy.zeros((n, m))

原因可以簡單理解為

n = 5
m = 3
matrix = [[0]*m]*n

# 相當於
"""
array = [0 0 0]
matrix = [array]*5
# matrix內的5個元素都是同一個列表引用
# 當使用 matrix[3][2] = 1 賦值
# 則 array[2] = 1
# 所以 matrix[0~4][2]都為 1
"""