[动态规划]矩阵连乘

矩阵连乘也是经典的动态规划的例子，而且能比较明显的把原来子问题的复杂度从指数级降为O(n^3)，算是效果明显。

【分析】

算法的分析请见参考资料，我就不重复写了。就是写代码的时候需要注意两点：

（1）这个程序是自底向上的动态规划。因为我们在一个段[ i : j ]中进行划分的时候，那么必须下面的每一个小段必须已经生成了。如下图所示：

所以不能采用传统的那种i和j从头到尾的遍历（这种双重循环一般只适合于从前向后生成，比如最长公共子序列等动态规划）。

（2）一开始对书中的递推式m[i][j] = m[i][k] + m[k][j] + pi-1 * pk * pj中的后面三项很不理解。为什么会是三项呢？书中矩阵Ai的行数是pi-1，那么第k项矩阵的行数和列数就分别是pk-1和pk。由于是把从i到j的矩阵分为两段，第一段是从i到k，第二段是从k+1到j。所以第一段计算完之后形成的矩阵应该是pi-1行和pk列；第二段计算完之后形成的矩阵应该是pk行和pj列（开头是第pk+1个矩阵），所以根据矩阵相乘的性质，这两个（分段）矩阵相乘的代价就是pi-1 * pk * pj了。

【代码】

下面是实现的C++代码，注意：代码中的循环方式；以及另一处notice。

#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>

using namespace std;

typedef vector<int> vi;
typedef vector<vi> Matrix;

struct aMatrix
{
     aMatrix( int r, int c ){ row = r; col = c; }
    int row;    // the number of rows
    int col; // the number of columns
};

// add brackets to the matrix multiplication
void addBrackets( Matrix& sep, int start, int end );
// fill in the minimum matrix and the separate point matrix
void MinMatrixMultiply( Matrix& min, vector<aMatrix>& p, Matrix& sep, int size );
//print function
void printMatrix( const Matrix& m );

int main()
{
     vector<aMatrix> p;
     p.push_back( aMatrix( 30, 35 ) );
     p.push_back( aMatrix( 35, 15 ) );
     p.push_back( aMatrix( 15, 5 ) );
     p.push_back( aMatrix( 5, 10 ) );
     p.push_back( aMatrix( 10, 20 ) );
     p.push_back( aMatrix( 20, 25 ) );

    int size = p.size();

     Matrix min( size, vector<int>( size, 0 ) );
     Matrix sep( size, vector<int>( size, 0 ) );

     MinMatrixMultiply( min, p, sep, size );
     printMatrix( min );
     addBrackets( sep, 0, 5 );
    

}

void MinMatrixMultiply( Matrix& min, vector<aMatrix>& p, Matrix& sep, int size )
{
    //because the matrix is initialized to all zero, so here we don't set min[i][i] to zero
    
    /*
     ** Notice!
     ** We can't use the usual way of i,j loop( i from 0 to n, j from i to n )
     ** because we use dynamic programming in a bottom-up way
     ** Now, i is the beginning, j is the end, while r is the distance from i to j
     ** so r is increasing from 1 to n.
    */
    for( int r = 1; r < size ; r++ )
      {
         for( int i = 0 ; i < size - r ; i++ )//i is the beginning
          {
              int j = i + r;// j is the end

              /*
               ** Notice !!!
               ** Now we can try all the separate points from i to j
               ** but we must initialize the min[i][j] by try k = i.
               ** otherwise the t has no value to compare to.
              */
               min[i][j] = min[i][i] + min[i+1][j] + p[i].row * p[i].col * p[j].col;
               sep[i][j] = i;
              for( int k = i + 1; k < j ; k++ )
               {
                   int t = min[i][k] + min[k+1][j] + p[i].row * p[k].col * p[j].col;
                   //save the minimum
                   if( t < min[i][j] )
                    {
                        min[i][j] = t;
                        sep[i][j] = k;
                    }
               }
          }
      }
}

void printMatrix( const Matrix& m )
{
    int size = m.size();
    for( int i = 0; i < size ; i++ )
     {
        for( int j = 0 ; j < size ; j++ )
         {
             cout<<std::setw(5)<<m[i][j]<<" ";
         }
         cout<<endl;
     }
}

void addBrackets( Matrix& sep, int start, int end )
{
    if ( start == end )
        return ;

     addBrackets( sep, start, sep[start][end] );
     addBrackets( sep, sep[start][end]+1, end );
    // print the first half
     cout<< "multiply (A"<<start<<" to A"<<sep[start][end]<<")";
    // print the second half
     cout<<" and (A"<<sep[start][end]+1<<" to A"<<end<<")"<<endl;
}

[参考文献]

王晓东。《计算机算法设计与分析》。P49

以及 http://blog.csdn.net/liguisen/archive/2008/03/08/2158127.aspx