前言 從旭東的博客 看到一篇博文：矩陣連乘最優結合 動態規劃求解，挺有意思的，這裡做個轉載【略改動】。 問題 矩陣乘法滿足結合律，但不滿足交換律。例如矩陣$A_{ab}, B_{bc}, C_{cd}$ 連乘得到矩陣$S_{ad}$ \[ S_{ad}=A_{ab} B_{bc} C_{cd} \] ...
前言
從旭東的博客 看到一篇博文：矩陣連乘最優結合 動態規劃求解，挺有意思的，這裡做個轉載【略改動】。
問題
矩陣乘法滿足結合律，但不滿足交換律。例如矩陣$A_{ab}, B_{bc}, C_{cd}$ 連乘得到矩陣$S_{ad}$
\[ S_{ad}=A_{ab} B_{bc} C_{cd} \]
實際運算可以先將矩陣A和B相乘，S=(AB)C，也可以先將矩陣B和C相乘，S=A(BC) 。這兩種演算法的計算量是否一樣呢？
先看A×B的計算量。A×B是a行c列矩陣，共a×c個元素。b是A、B兩矩陣的乘法介面，積矩陣的每個元素都需經過b次乘法運算。那麼A×B需要a×b×c次乘法運算。
(AB)C乘法運算次數為：abc+acd=ac(b+d)
A(BC)乘法運算次數為：bcd+abd=bd(a+c)
由於a、b、c、d的取值是任意的， 所以ac(b+d)和bd(a+c)不可能保持相同，也即結合律可以用於矩陣乘法優化。那麼，對於連乘AiAi+1…Aj，如何找出最優的運算順序？
分析
這個問題長得像遞歸，但是用遞歸似乎不好做，因為子問題裡面有很多分支，那子問題的返回值是什麼，你返回哪個分支的返回值？而且就算能實現，肯定會對一些相同的子問題重覆計算。
由於計算n個矩陣的連乘，總是通過計算更少個數的矩陣連乘實現的，最終化歸到計算兩個矩陣的相乘，所以考慮從兩個矩陣相乘著手。這樣自底向上，就把問題解決了。當我們把所有的兩個矩陣相乘的乘法運算次數求出後，便可以去計算三個矩陣相乘的乘法運算次數。這時，三個矩陣的連乘相當於兩個矩陣相乘。這樣，所有的矩陣連乘都是兩個矩陣相乘，其乘法運算次數為兩個矩陣各自的乘法運算次數之和再加上兩個矩陣相乘的乘法運算次數。對於兩個矩陣A、B相乘的乘法運算次數$\beta$，上面A×B的示例已經講過了。進一步可以抽象為
\[\beta=行*介面*列\]
按上述需要，我們將矩陣連乘表達式視為一系列子式。每一個子式，由起點和跨度唯一決定。起點和跨度剛好構成兩個維度，於是可以用n階方陣索引子式，比如索引[1][2]表示起點為1跨度為2的子式。起點至多有n個，最小為0，最大為n-1。跨度至多n種，最小為0，最大為n-1。用兩個n階方陣分別存儲每個子式的最少計算次數及分割點。C++提供了vector模板類，作為數組的一個替代品，我們用vector實現n階方陣。首先確定跨度，再確定起點，構成一個二級嵌套迴圈，通過起點的平移確定子式。再嵌套一級遍歷子式的分割點求出子式的最少計算次數。最後我們得到最大的子式，也就是連乘本身的最少運算次數及分割方案。
代碼如下：
1 #include <iostream> 2 #include <vector> 3 #include <algorithm> 4 #include <limits.h> 5 #include <string> 6 7 using namespace std; 8 9 //計算矩陣的分割方案（基於向量） 10 int calculate_M(vector<vector<int> >&num, vector<pair<int, int> > &data, vector<vector<int> > &points) { 11 int matrixNum = data.size(); 12 int span;//起止點間隔距離，表示跨度 13 int start;//起點 14 int end;//終點 15 int spiltPoint;//分割點 16 int multiplyTimes;//臨時存儲子式的乘法運算次數 17 for (span = 1; span < matrixNum; span++) { ///間隔距離 相鄰矩陣的間隔距離為1 18 for (start = 0; start < matrixNum - span; start++) { ///起點平移 19 //子式確立，下麵開始計算這個子式的最少乘法運算次數，初值先設為最大整數 20 end = start + span; 21 num[start][end] = INT_MAX; 22 for (spiltPoint = start; spiltPoint < end; spiltPoint++) { 23 24 multiplyTimes = num[start][spiltPoint] + num[spiltPoint + 1][end] + data[start].first*data[spiltPoint].second*data[end].second; 25 if (multiplyTimes < num[start][end]) { 26 points[start][end] = spiltPoint; ///記錄分割點 27 num[start][end] = multiplyTimes; ///記錄最少乘法次數 28 } 29 } 30 } 31 } 32 return 0; 33 } 34 35 //根據記錄的分割點，生成最後的矩陣相乘表達式 36 string make_result(vector<vector<int> > &points, int t1, int t2) { 37 if (t1 == t2) 38 return string(1, 'A' + t1); 39 int split = points[t1][t2]; 40 return "(" + make_result(points, t1, split) + "*" + make_result(points, split + 1, t2) + ")"; 41 } 42 43 44 int main() 45 { 46 int matrixNum=0;///連乘的矩陣個數 47 vector<pair<int, int>> data; ///存儲矩陣的尺寸 48 49 //固定輸入 50 data.push_back(make_pair(10, 100)); //A 51 data.push_back(make_pair(100, 5)); //B 52 data.push_back(make_pair(5, 25)); //C 53 data.push_back(make_pair(25, 15)); //D 54 data.push_back(make_pair(15, 20)); //E 55 matrixNum = 5; 56 57 //為記錄乘法次數和分割點的n階矩陣分配空間並初始化 58 vector<vector<int> > num(matrixNum, vector<int>(matrixNum)); 59 vector<vector<int> > points(matrixNum, vector<int>(matrixNum)); 60 for (int i = 0; i < matrixNum; i++) { 61 for (int j = 0; j < matrixNum; j++) { 62 points[i][j] = -1; 63 num[i][j] = 0; 64 } 65 } 66 67 calculate_M(num, data, points); 68 cout << make_result(points, 0, matrixNum - 1) << "\t最少乘法次數為:" << num[0][matrixNum - 1] << endl; 69 70 cin >> matrixNum; 71 return 0; 72 }View Code
代碼輸出為：
| ((A*B)*((C*D)*E)) 最少乘法次數為:9375 |
如果從控制台輸入矩陣尺寸的話，可以輸入第一個矩陣的行數和各矩陣的列數。直接使用這個序列存儲矩陣尺寸。
1 #include <iostream> 2 #include <vector> 3 #include <algorithm> 4 #include <limits.h> 5 #include <string> 6 7 using namespace std; 8 9 int* uncertainInput(int& inputNum); 10 11 //計算矩陣的分割方案（基於數組） 12 int calculate_M1(int*data, int matrixNum, vector<vector<int> >&num, vector<vector<int> > &points) { 13 14 int span;//起止點間隔距離，表示跨度 最大跨度比總的矩陣數小1 15 int start;//起點 16 int end;//終點 17 int spiltPoint;//分割點 分割線在分割點與下一點之間 18 int multiplyTimes;//臨時存儲子式的乘法運算次數 19 int rows, columns, interfaces;//行 列 介面 20 21 for (span = 1; span < matrixNum-1; span++) { ///間隔距離 相鄰矩陣的間隔距離為1 22 for (start = 1; start + span < matrixNum; start++) { ///起點從第一個矩陣的列開始 23 //子式確立，下麵開始計算這個子式的最少乘法運算次數，初值先設為最大整數 24 end = start + span; 25 num[start][end] = INT_MAX; 26 for (spiltPoint = start; spiltPoint < end; spiltPoint++) { 27 rows = *(data + start-1); 28 columns = *(data + end); 29 interfaces = *(data + spiltPoint); 30 31 multiplyTimes = num[start][spiltPoint] + num[spiltPoint + 1][end] + rows * interfaces * columns; 32 if (multiplyTimes < num[start][end]) { 33 points[start][end] = spiltPoint; ///記錄分割點 34 num[start][end] = multiplyTimes; ///記錄最少乘法次數 35 } 36 } 37 } 38 } 39 return 0; 40 } 41 42 //根據記錄的分割點，生成最後的矩陣相乘表達式 43 string make_result(vector<vector<int> > &points, int t1, int t2) { 44 if (t1 == t2) 45 return string(1, 'A' + t1 -1); 46 int split = points[t1][t2]; 47 return "(" + make_result(points, t1, split) + "*" + make_result(points, split + 1, t2) + ")"; 48 } 49 50 51 int main() 52 { 53 int matrixNum = 0;///連乘的矩陣個數 54 vector<pair<int, int>> data; ///存儲矩陣的尺寸 55 56 //手動輸入 57 //輸入的第一個值是第一個矩陣的行，其餘是各矩陣的列 58 int* matrixSizeSeq; 59 matrixSizeSeq = uncertainInput(matrixNum); 60 if (matrixSizeSeq == NULL) 61 return false; 62 63 //重構為向量 64 /*for (int ptrOffset = 0; ptrOffset < matrixNum - 1; ptrOffset++) ///遍歷到倒數第二個輸入 65 data.push_back(make_pair(*(matrixSizeSeq + ptrOffset), *(matrixSizeSeq + ptrOffset + 1)));*/ 66 67 //為記錄乘法次數和分割點的n階矩陣分配空間並初始化 68 vector<vector<int> > num(matrixNum, vector<int>(matrixNum)); 69 vector<vector<int> > points(matrixNum, vector<int>(matrixNum)); 70 for (int i = 0; i < matrixNum; i++) { 71 for (int j = 0; j < matrixNum; j++) { 72 points[i][j] = -1; 73 num[i][j] = 0; 74 } 75 } 76 77 calculate_M1(matrixSizeSeq, matrixNum, num, points); 78 cout << make_result(points, 1, matrixNum - 1) << "\t最少乘法次數為:" << num[1][matrixNum - 1] << endl; 79 cin >> matrixNum; 80 return 0; 81 } 82 83 int* uncertainInput(int& inputNum) 84 { 85 //輸入0表示輸入完成 輸入負數認為手誤 86 using namespace std; 87 88 int allocateSpace = 5; 89 int* inputSequence = (int *)malloc(allocateSpace * sizeof(int)); ///存儲矩陣的尺寸輸入 90 91 92 int ptrOffset = 0; ///指針偏移 93 int tempSize = 0; ///臨時存儲輸入值 94 95 do { 96 cin >> tempSize; 97 if (tempSize <= 0) 98 continue; 99 *(inputSequence + ptrOffset++) = tempSize; 100 if (ptrOffset >= allocateSpace) 101 { 102 allocateSpace += 5; 103 inputSequence = (int*)realloc(inputSequence, allocateSpace * sizeof(int)); 104 } 105 } while (tempSize != 0); 106 inputNum = ptrOffset; ///偏移量即輸入個數 107 return inputSequence; 108 }
比如求$S_{ad}$，則在控制台輸入a b c d 0，然後回車即可。當然 a b c d 必須替換為數字。
代碼中用到的一些知識
C++提供模版類string，其中一個構造方法可將字元轉化為字元串。如 string(1, 'A'+1)，第一個參數是源字元延拓次數，這個構造函數將‘B’轉化為"B"。
