[动态规划]矩阵连乘

矩阵连乘也是经典的动态规划的例子，而且能比较明显的把原来子问题的复杂度从指数级降为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