学步园

    对于"矩阵连乘问题"的一点想法

 在算法设计的学习中，每到“动态规划”一节，一般都会涉及到“矩阵连乘”问题（例如《Algorithms》，中文译名《算法概论》），可想而知该题的经典程度 ：）

  前些天复习动态规划的时候，瞅着这个问题突然有了一点有趣的想法：难道该题只能以动态规划求解吗？能不能用……（暂且卖个关子 ：））？遂而有了以下的一些思考。

  首先，让我不厌其烦的再次回味一遍“矩阵连乘”,算法大牛们可以直接无视：

  显而易见，在矩阵个数为1个或者2个时，该问题易解，在此考虑一种3个矩阵的情况：

A1：10×100，A2：100×5，A3：5×50；并且我们简记两个矩阵相乘为：矩阵1矩阵2，

例如A1乘以A2记为：A1A2，并用 () 表示矩阵之间的相乘顺序，例如A1乘以A2再乘以

A3，我们便记作：((A1A2)A3)，另外，就一个a行b列的矩阵与一个b行c列的矩阵相乘而言（注意，必须满足矩阵可乘条件），其需要的乘法次数为：aXbXc（不相信的朋友可以自己演算：）），对于上面提到的这种情况，如果采用 ((A1A2)A3) 的连乘顺序，那么总的乘法次数为：10X100X5 + 10X5X50 = 7500；而采用 (A1(A2A3))的连乘顺序，那么总的连乘次数为：100X5X50 + 10X100X50 = 75000 ！两种连乘顺序的乘法次数竟然相差十倍之巨！可想而知一个好的矩阵连乘顺序选择是多么重要。

  至于如何解决这个“矩阵连乘”问题，一般都采用动态规划方法，具体思路如下：

对于一连串的矩阵相乘，我们定义问题 P(i,j) ( j >= i ) ：原矩阵链中矩阵Ai至Aj之间的矩阵

连乘最小次数，显而易见，原问题是该问题的一个子问题，P(1,n)即代表原问题的解，并且 

P(i,j)( 1>= j - i >=0 ) 的解都是易解的，或者说平凡的，那么，对于这个自定义的问题，我们很自然的可以总结出以下的递推公式：

                         0 ( i = j )

P(i,j) (j>=i) =      axbxc ( i = j+1 )

                    Min( P(i,k)+p(k+1,j)+axmxc ( k>=i 且 k <=j )) (其他情况)

（ 设Ai为 axb 矩阵，Aj为 bxc 矩阵，Ak为 nxm 矩阵 ）

  显而易见，这是一个天生的动态规划问题（良好的递归问题定义，以及诸多重复的子

问题计算），那么接下来，就让我们继续深入细节，编码来实现这个算法，由于递归公式已

经给出，实际编码其实并无多大问题，需要注意的可能就是子问题的求解顺序了：

/*<br /> Name: matrixs_chain_mutiply_dp.cpp<br /> Copyright: No Copyright<br /> Author: Hugo Yu<br /> Date: 09-01-10 21:47<br /> Description: A Simple Implementation Of matrix chain-multiplication<br /> ( Dynamic Programming )<br /> */</p> <p>#include <iostream><br /> using std::endl;<br /> using std::cout;<br /> using std::cin;<br /> using std::istream;</p> <p>#include <cstddef><br /> using std::size_t;</p> <p>struct SimpleMatrix<br /> {<br /> //friend istream& operator >> ( istream& input, SimpleMatrix& sm );<br /> size_t rows;<br /> size_t cols;<br /> };</p> <p>//the max matrixs count ( just for simple <img src="https://admin.xuebuyuan.com/wp-includes/images/smilies/icon_smile.gif" alt=":)" class="wp-smiley" /> )<br /> const size_t MAX_MATRIXS_COUNT = 512;<br /> //matrix array<br /> SimpleMatrix g_sm_array[MAX_MATRIXS_COUNT];<br /> //buffer<br /> size_t g_buffer[MAX_MATRIXS_COUNT][MAX_MATRIXS_COUNT];<br /> //max size_t<br /> const size_t MAX_SIZE_T = size_t(-1);</p> <p>istream& operator >> ( istream& input, SimpleMatrix& sm )<br /> {<br /> input >> sm.rows >> sm.cols;<br /> return input;<br /> }</p> <p>size_t Calculate( size_t matrix_count )//just use g_sm_array<br /> {<br /> //initiate the g_buffer<br /> for( int i = 0; i < matrix_count; ++i )<br /> {<br /> g_buffer[i][i] = 0;<br /> }<br /> for( int i = 0; i < matrix_count - 1; ++i )<br /> {<br /> g_buffer[i][i+1] = g_sm_array[i].rows * g_sm_array[i].cols * g_sm_array[i+1].cols;<br /> }</p> <p> //do the calculate<br /> for( int k = 2; k < matrix_count; ++k )<br /> {<br /> for( int i = 0; i <= matrix_count - k; ++i )<br /> {<br /> size_t min_count = MAX_SIZE_T;//just for simple<br /> for( int j = i; j < i + k; ++j )<br /> {<br /> size_t result = g_buffer[i][j] + g_buffer[j+1][i+k] + g_sm_array[i].rows * g_sm_array[j].cols * g_sm_array[i+k].cols;<br /> if( result < min_count ) min_count = result;<br /> }<br /> g_buffer[i][i+k] = min_count;<br /> }<br /> }<br /> return g_buffer[0][matrix_count-1];<br /> }</p> <p>int main()<br /> {<br /> size_t matrix_count = 0;<br /> cout << "Please Input The Matrixs Count : " << endl;<br /> cin >> matrix_count;</p> <p> cout << "Please Input The Matrixs : " << endl;<br /> //SimpleMatrix *sm_array = new SimpleMatrix[matrix_count];//should be more safe here <img src="https://admin.xuebuyuan.com/wp-includes/images/smilies/icon_smile.gif" alt=":)" class="wp-smiley" /><br /> for( int i = 0; i < matrix_count; ++i )<br /> {<br /> cin >> g_sm_array[i];<br /> }</p> <p> cout << "The Result Is : "<<<br /> Calculate( matrix_count )<<endl;</p> <p> //delete[] sm_array;//should be more safe here <img src="https://admin.xuebuyuan.com/wp-includes/images/smilies/icon_smile.gif" alt=":)" class="wp-smiley" /><br /> system( "pause" );<br /> return 0;<br /> }

  好了，至此经典的算法已经讲述完毕，接下来该轮到我自个儿瞎扯几句了 ：）

  还拿上面的那个例子来讲：A1：10×100，A2：100×5，A3：5×50，我首先注意到

其中100这个数字，毫无疑问，它是所有矩阵行列数中最大的，并且如果我们采用这种矩阵连乘顺序：(A1(A2A3)) ，那么，100首先作为A2矩阵的行数参加了一次乘法运算(即100×5×50)，然后再次作为A2矩阵的行数(或者说A1矩阵的列数)参加了一次乘法运算(即10×100×50)，也就是说，100作为各矩阵行列数中最大的一个数，在上述的这种矩阵连乘顺序中，总共参加了两次乘法运算，即产生了两次作用！毫无疑问这定会带来很大的乘法次数，那么我的想法是，能不能首先将这个100消除，以使其产生的效果最小，答案很简单：首先计算 A1A2，而这种乘法方案恰恰就是该矩阵情况下矩阵连乘的最优顺序方案！那么很自然的，我便有了以下一个“贪心”的算法：

  1.首先按照各矩阵的共有行列数排序（所谓共有行列数，举例来说，对于axb矩阵A1及bXc矩阵A2的共有行列数即为共有的b，暂且不知道更为准确的叫法，就容我先这么叫吧 ：））

  2.然后按照排序顺序（从大到小）进行矩阵乘法。

  好了，算法就是这么简单，并且显而易见，该算法属于贪心法的范畴，而往往贪心法总是过程简单，证明复杂，我也不出意外的卡在了证明准确性的关卡上，思考了一些方法但感觉仍然不得要领，所以我为此还特意去请教了一下我们学院的陈院长（算法大牛），牛人就是牛人，我刚刚提了“贪心”及“矩阵连乘”这两个关键字，还没有透露任何细节，我们的陈院长便将该算法八九不离十的自己推测出来了，并且还给出了一个解释（类似于之前我的思考），但我认为这般论据仍然不够充分，而陈院长则认为这个算法可能是正确的，所以目前对于更深层的原因，我仍然不得而知，期望某位算法牛人可以不吝赐教，小弟在此先拜谢了：）但不能深层理解并不代表不能加以实现，所以在此我也编写了一段简单的代码（虽说简单，但仍然觉得作为示例还是复杂了些）：

/*<br /> Name: matrixs_chain_mutiply_gd.cpp<br /> Copyright: No Copyright<br /> Author: Hugo Yu<br /> Date: 10-01-10 19:57<br /> Description: A Simple Implementation Of matrix chain-multiplication<br /> ( Greed algorithm )<br /> */</p> <p>#include <iostream><br /> using std::endl;<br /> using std::cout;<br /> using std::cin;<br /> using std::istream;</p> <p>#include <cstddef><br /> using std::size_t;</p> <p>#include <algorithm><br /> using std::sort;</p> <p>struct SimpleMatrix<br /> {<br /> //friend istream& operator >> ( istream& input, SimpleMatrix& sm );<br /> size_t rows;<br /> size_t cols;<br /> };</p> <p>istream& operator >> ( istream& input, SimpleMatrix& sm )<br /> {<br /> input >> sm.rows >> sm.cols;<br /> return input;<br /> }</p> <p>struct Pair<br /> {<br /> bool operator < ( const Pair& pair ) const { return this->element < pair.element; }<br /> size_t element;<br /> size_t index;<br /> };</p> <p>//the max matrixs count ( just for simple <img src="https://admin.xuebuyuan.com/wp-includes/images/smilies/icon_smile.gif" alt=":)" class="wp-smiley" /> )<br /> const size_t MAX_MATRIXS_COUNT = 512;<br /> //matrix array<br /> SimpleMatrix g_sm_array[MAX_MATRIXS_COUNT];<br /> //buffer<br /> Pair g_buffer[MAX_MATRIXS_COUNT];<br /> //max size_t<br /> const size_t MAX_SIZE_T = size_t(-1);</p> <p>size_t Calculate( size_t matrix_count )<br /> {<br /> //initiate the g_buffer<br /> for( int i = 0; i < matrix_count - 1; ++i )<br /> {<br /> g_buffer[i].element = g_sm_array[i].cols;<br /> g_buffer[i].index = i;<br /> }<br /> sort( g_buffer, g_buffer + matrix_count - 1 );<br /> size_t result = 0;<br /> for( int i = matrix_count - 2; i >= 0; --i )<br /> {<br /> size_t index = g_buffer[i].index;<br /> //calculate result<br /> result += ( g_sm_array[index].rows * g_sm_array[index].cols * g_sm_array[index+1].cols );<br /> //update matrixs array<br /> for( int j = index + 2; j < matrix_count - 1; ++j )<br /> {<br /> if( ( g_sm_array[index+1].rows != g_sm_array[j].rows ) ||<br /> ( g_sm_array[index+1].cols != g_sm_array[j].cols ) )<br /> break;<br /> g_sm_array[j].rows = g_sm_array[index].rows;<br /> }<br /> g_sm_array[index+1].rows = g_sm_array[index].rows;</p> <p> for( int j = index - 1; j >= 0; --j )<br /> {<br /> if( ( g_sm_array[index].rows != g_sm_array[j].rows ) ||<br /> ( g_sm_array[index].cols != g_sm_array[j].cols ) )<br /> break;<br /> g_sm_array[j].cols = g_sm_array[index+1].cols;<br /> }<br /> g_sm_array[index].cols = g_sm_array[index+1].cols;<br /> }<br /> return result;<br /> }</p> <p>int main()<br /> {<br /> size_t matrix_count = 0;<br /> cout << "Please Input The Matrixs Count : " << endl;<br /> cin >> matrix_count;</p> <p> cout << "Please Input The Matrixs : " << endl;<br /> //SimpleMatrix *sm_array = new SimpleMatrix[matrix_count];//should be more safe here <img src="https://admin.xuebuyuan.com/wp-includes/images/smilies/icon_smile.gif" alt=":)" class="wp-smiley" /><br /> for( int i = 0; i < matrix_count; ++i )<br /> {<br /> cin >> g_sm_array[i];<br /> }</p> <p> cout << "The Result Is : "<<<br /> Calculate( matrix_count )<<endl;</p> <p> //delete[] sm_array;//should be more safe here <img src="https://admin.xuebuyuan.com/wp-includes/images/smilies/icon_smile.gif" alt=":)" class="wp-smiley" /><br /> system( "pause" );<br /> return 0;<br /> }

  最后我将该程序产生结果与上面程序（动态规划）产生结果加以比较，所使用数据则是由一个简单的随机函数产生：

//max matrix rows or cols<br /> const size_t MAX_ROW_COL = 100;</p> <p>void RandomMatrixs( size_t matrix_count )<br /> {<br /> g_sm_array[0].rows = rand() % ( MAX_ROW_COL - 1 ) + 1;//avoid get 0<br /> g_sm_array[matrix_count-1].cols = rand() % ( MAX_ROW_COL - 1 ) + 1;<br /> for( int i = 0; i < matrix_count - 1; ++i )<br /> {<br /> g_sm_array[i].cols = g_sm_array[i+1].rows = rand() % ( MAX_ROW_COL - 1 ) + 1;<br /> }<br /> }

  很快便得到了一组结果不一致的数据：A1 42X98， A2 98X68，A3 68X63，A4 63X83，A5 83X54。 

  使用动态规划的正规方法，所得结果为：867678（ （（（（A1A2）A3）A4）A5） ），而贪心算法的结果却为：885066（ （（（A1A2）A3）（A4A5）） ），显然，矩阵连乘问题使用贪心法是错误的 ：（ 。

  现在的理解是，贪心的局部最优，在“矩阵连乘”问题中并不会导致全局最优，也就是说我对于本题的看法还是落入了“短视”的窠臼，不过明晰的数学分析抑或缜密的证明推断，现在的我还是无能为力（囧...），再次渴望一下大牛们的谆谆教诲 ：），不过最为“矩阵连乘”问题的近似算法，我想也许这个贪心思路能够带来一点启示 ：）

  好了，思考暂时便是这么多了，我想也是时候休息一下了（譬如玩玩《KOF》或者《SF4》） ：）