动态规划之矩阵链乘

问题提出：（代码下载）

对于如下矩阵：

 

 其中各矩阵A[i]下标为

计算其乘积的结果，以及我们需要计算其最小标量乘法次数。

问题分析：

首先我们需要明确的是何为标量：标量即为没有方向的量，而有方向的量即为矢量。（严谨的定义自己百度去）

那么标量乘法就变成了最基本的数字相乘。

其次对于两个矩阵相乘，需满足下示公式所示的形式：（左边矩阵的列数与右边矩阵的行数必须一致）

上述条件可从矩阵相乘的定义中看出：

在计算机，我们可以用一个二维数组来表示矩阵。

一个m行n列的矩阵与一个n行p列的矩阵相乘，会得到一个m行p列的矩阵，具体步骤如下：

其中第i行第j列位置上的数为第一个矩阵的第i行上的n个数与第二个矩阵第j列上的n个数对应相乘后所得的n个乘积之和。

下图展示了一个简单的例子：

 

 

 

从上图我们可以得出两个矩阵相乘城的次数，即为取得结果矩阵每个点所要做的标量乘法次数之和

而每个点的标量乘法次数为左边矩阵的列数n，结果矩阵有m*p个点，那么总数即为m*n*p；

另外我们也应该知道矩阵乘法符合乘法结合律，但不符合交换律；也就是说，我们可以为矩阵链乘设置括号，括号放置的位置不同，产生标量乘法的总量也是不一样的哦，这一点我们必须明确，对这点不清楚的，可以自己举个例子试试看。如果没有这点，那这题就没意义了。

前面科普了一大堆知识，现在进正题了：

对于矩阵链乘，我们可知道不同的括号化方案，会产生不同的计算代价；

算法导论书中指出了求完全括号化方案，所谓矩阵乘积链为完全括号化需满足如下性质：

要么是单一矩阵，要么是两个完全括号化的矩阵乘积链的积，且已加外括号。说实话这个定义让人有点摸不着头脑，待我仔细解释下：

首先：该乘积链中只有一个矩阵，即单一矩阵，该乘积链必是完全括号化，且该矩阵外部不用加括号。

再者：必须是两个括号化的矩阵链的乘积，我们可以把两个括号化的矩阵链看作两个结果矩阵（即两个单一矩阵），该两个矩阵相乘，并加上外括号把他们括住，即这种形式也叫着完全括号化的。

啰嗦了这么多，感觉简单的说法就是：单一矩阵（不加括号）是完全括号化，两个矩阵（单一矩阵）加上外括号也是完全括号化的，完全括号化的乘积链可以当作单一矩阵处理。

下面上一道课后的证明题：

对n个元素的表达式进行完全括号化时，恰好需要n-1对括号。

证明：

注：采用数学归纳法

已知仅有一个元素时，不需要括号，即0对括号。

假设k个元素时，需要k-1对括号；

当有k+1个元素时，我们可以看作k个元素的完全括号化的乘积链外加1个新加的元素。已知我们可以把完全括号化的乘积链看作单一矩阵，那么现有两个单一矩阵，那么需再加1对括号，才能再次完全括号化，而由假设知，k个元素的完全括号化需k-1对括号，那么k+1个元素的完全括号化方案就需要k-1+1对括号；

命题得证！

额，又科普了部分知识，只不过这些科普对后续理解问题都是有用的。

应用动态规划的方法的步骤：（出自书本）

1、刻画一个最优解的结构特征

2、递归地定义最优解的值

3、计算最优解的值，通常采用自底向上的方法

4、利用计算出的信息构造一个最优解的值

第一步：刻画一个最优解的结构特征

对于最开始提出的问题，我们看作为矩阵链乘设计最优的括号化方案。而设置括号的最终目的就是为了打乱计算顺序，通过结合律来找到最小计算代价，就比如下面一个问题33*4*25，是个正常人就会先计算后面两项，因为这样好计算，同样计算矩阵，我们也希望计算的少，计算的快啊，设置括号就变的尤为重要了。

对于完全括号化，每一个括号的作用是将一个链乘转化成两个完全括号化的链乘的积，前面有提到。

因此我最优解的结构特征就如下了：

对于第m个矩阵到第n个矩阵，他的最优解存在于对它的一次划分中，它的划分有n-m种情况，我们对这些情况中代价取最小，不就获得了最优解

而对于每种情况，我们需要划分后的两部分矩阵链都是最优解，乘积后才可能最小啊，（这部分大家可以用“剪切--粘贴”去反证），这样问题就转换成了两个子部分矩阵链的最优括号化方案的问题，以此不断的递归下去。

第二步：递归地定义最优解的值

我们的划分点k可从m到n-1，代价如下所示，每次的代价如下，而最优解只有一个，即它们的最小值：

其一般化的公式如下：

从这里可以看出，这是一个递归的问题，看着很类似上篇所述的钢条切割问题啊，仔细比对比对吧。

若采用递归策略，可能会有很多的重复子问题，无疑提高了计算的复杂度。也可以采取带备忘的递归来降低复杂度。

第三步：计算最优解的值，通常采用自底向上的方法

另外从公式中可以看出，上述计算结果要依赖于m,n（m<n）的所有可能组合。区间[m，n]组合要依赖于它们之间的子区间组合。所以我们采取自下而上的策略，逐渐扩大m,n之间的区间长度，最终算到m=0,n=n，从而获得最优的括号化方案，

具体方法为：matrix_chain_mutiply.cpp中的void matrix_divide_option(int *num,int** cost,int** point,int length)方法。

第四步、利用计算出的信息构造一个最优解的值

具体解析见matrix_chain_mutiply.cpp最后两个方法

代码如下：

matrix.h文件：

#include<iostream>
#ifndef MATRIX_H
#define MATRIX_H
/*
矩阵类
row为其行数
col为该矩阵的列数，
value表示一个二维数组，主要存储该矩阵各点的值
*/
class matrix
{
    public:
        matrix();                                //默认的构造函数
        matrix(int row,int col);                 //给定行数和列数，构造一个矩阵
        virtual ~matrix();                       //析构函数
        matrix(const matrix& other);             //拷贝构造函数
        matrix operator*(const matrix& a);       //矩阵的乘法
        matrix& operator=(const matrix&);        //重写赋值函数
        matrix init();                           //初始化矩阵
        friend std::ostream& operator << (std::ostream&,const matrix &); //重写输出
    protected:
    private:
        double **value;                         //主要用于存储矩阵各点的值
        int col;                                 //列数
        int row;                                 //行数
};
#endif // MATRIX_H

matrix.cpp

#include<iostream>
#include "matrix.h"
#include<cstdlib>
#include<ctime>

matrix::matrix()                                //默认构造函数
{
    col = row = 0;
    value = nullptr;
}

matrix::matrix(int row, int col)                //根据给定行，列数来初始化一个矩阵，并分配相应的内存
{
    this->row = row;
    this->col = col;
    value = new double*[row];
    for(int i=0;i<row;i++)
    {
        value[i] = new double[col];
    }
}

matrix::matrix(const matrix& other)             //复制构造函数，注意此处要用matrix&，重写复制构造函数为了防止值复制造成
{
    row = other.row;
    col = other.col;
    value = new double*[row];                   //此处使用new分配内存
    for(int i=0;i<row;i++)
    {
        value[i] = new double[col];
    }
    for(int i=0;i<row;i++)
    {
        for(int j=0;j<col;j++)
        {
            value[i][j]=(other.value)[i][j];
        }
    }
}

matrix::~matrix()
{
    for(int i=0;i<row;i++)                      //因为构造函数中new分配有内存，所以在析构函数中要对内存进行释放。
    {
        delete[] value[i];
    }
    delete[] value;
}


matrix& matrix::operator=(const matrix& other)  //重写了赋值函数
{
    if(this == &other)                          //防止a=a的情况出现
    {
        return *this;
    }
    for(int i=0;i<row;i++)                      //先释放原有区域的内存，从此处我们可以上面的判断的重要性
    {                                           //若没有上述判断，a=a在赋值的同时，也将内存释放了，造成数据丢失
        delete[] value[i];
    }
    delete[] value;
    row = other.row;                            //将成员变量进行复制
    col = other.col;
    value = new double*[row];                   //重新分配内存
    for(int i=0;i<row;i++)
    {
        value[i] = new double[col];
    }
    for(int i=0;i<row;i++)                      //对二维数组进行复制
    {
        for(int j=0;j<col;j++)
        {
            value[i][j]=(other.value)[i][j];
        }
    }
    return *this;
}

matrix matrix::operator*(const matrix& n)       //矩阵的乘法
{
    matrix result(row,n.col);
    for(int i=0;i<row;i++)
    {
        for(int j=0;j<n.col;j++)
        {
            (result.value)[i][j] = 0;
            for(int k=0;k<n.row;k++)
            {
                (result.value)[i][j] += value[i][k]*(n.value)[k][j];
            }
        }
    }
    return result;
}

matrix matrix::init()                          //初始化矩阵
{
    matrix result(row,col);
    srand((unsigned)time(NULL));
    for(int i=0;i<row;i++)
    {
        for(int j=0;j<col;j++)
        {
            (result.value)[i][j] = rand()%100;
        }
    }
    return result;
}

std::ostream& operator <<(std::ostream& os,const matrix& m)//将矩阵输出
{
    for(int i=0;i<m.row;i++)
    {
        for(int j=0;j<m.col;j++)
        {
            os << (m.value)[i][j] << '	';
        }
        os << std::endl;
    }
    return os;
}

 matrix_chain_mutiply.cpp(核心)

#include<iostream>
#include "matrix.h"
const int SIZE = 5;                                     //SIZE的大小决定了矩阵的数量（SIZE-1），可自行改变
const int max_size = 2147483647;

void matrix_divide_option(int *,int **,int **,int);     //确定最优的括号化方案
void show(int **);                                      //输出二维数组
void print(int**,int,int);                              //输出矩阵链乘的括号化方案
matrix matrix_mutiply(int**,int,int,matrix*);           //获得矩阵链乘的最终结果矩阵
int main()
{
    using namespace std;
    cout << "input " << SIZE << " number(press enter to end your input)：" << endl;
    int *num = new int[SIZE];                           //用于存储矩阵A[i]下标为num[i]*num[i+1];
    int i=0;
    cin >> num[i];                                      //获得矩阵的行列数
    while(cin.get() != '
')
    {
        i++;
        cin >> num[i]; //读取数组维数；
    }
    cout << endl << "what you input is:" << endl;
    for(int j=0;j<SIZE;j++)
    {
        cout << num[j] << "-----";                      //展示输入的是否正确
    }
    cout << endl << endl;
    matrix* test = new matrix[SIZE-1];                  //为矩阵数组分配内存
    for(int j=1;j<SIZE;j++)
    {
        matrix m(num[j-1],num[j]);                      //初始化矩阵大小
        test[j-1] = m.init();                           //初始化该矩阵各点的值
    }
    cout << "矩阵信息如下：" << endl;
    for(int j=0;j<SIZE-1;j++)
    {
        cout << test[j] << endl;                        //输出矩阵信息
    }

    int **cost = new int*[SIZE-1];                      //存储矩阵相乘时，各标量乘法次数，即代价表
    int **point = new int*[SIZE-1];                     //存储矩阵括号化方案时的分割点
    for(int j=0;j<SIZE-1;j++)
    {
        cost[j] = new int[SIZE-1];                      //分配内存
        point[j] = new int[SIZE-1];
    }
    matrix_divide_option(num,cost,point,SIZE-1);        //获取划分方案，获得最小的代价
    cout << "各区间的计算代价：" << endl;
    show(cost);
    cout << endl;
    cout << "区间划分点所在的位置：" << endl;
    show(point);
    cout << endl;
    cout << "输出完全括号化方案：" << endl;
    print(point,0,SIZE-2);
    cout << endl << endl;
    cout << "输出结果矩阵：" << endl;
    matrix result = matrix_mutiply(point,0,SIZE-2,test);
    cout << result;
    for(int j=0;j<SIZE-1;j++)
    {
        delete[] cost[j];
        delete[] point[j];
    }
    delete[] cost;
    delete[] point;
    delete[] test;
    delete[] num;
    return 0;
}


/*
此方法主要获得计算代价，以及括号化方案的划分点
为此问题解决的核心方法。
num存储的矩阵的下标，num[i-1]为第i个矩阵的行数，num[i]表示第i个矩阵的列数
因为数组存储时，最初始的下标为0，所以下述的定义就变的有点不好理解！仔细想想也容易明白。
cost----二维数组，cost[i][j]代表的就是第（i+1）个矩阵到第（j+1）个矩阵的链乘的最小计算代价
point----二维数组，point[i][j]代表的就是第（i+1）个矩阵到第（j+1）个矩阵的链乘的最佳划分处，假如为k，则k就是第k+1个矩阵，
（i+1）到(k+1)矩阵属于左边的完全括号化矩阵，其余为右侧完全括号化矩阵
length代表的是矩阵的个数，即num的长度-1
*/

void matrix_divide_option(int *num,int** cost,int** point,int length)
{
    //首先将代价表中的值重置
    //下标i到下标j之间的矩阵链乘的最小值
    for(int i=0;i<length;i++)
    {
        for(int j=0;j<length;j++)
        {
            if(i>=j)
            {
                cost[i][j] = 0;                 //单个矩阵不存在乘积，代价0，同时也将不合理的下标组合置为0
            }
            else
            {
                cost[i][j] = max_size;          //先将各合理的组合代价设为最大，以便后续比较使用。
            }
            point[i][j] = -1;                   //将所有组合的划分点置为-1，以便后续记录。
        }
    }
/*
    l代表矩阵链中矩阵的个数，我们这里从2个矩阵起步.由公式可知，我们要算遍所有的合理组合，同时大区间的组合
    还要依赖于其子区间的合理组合，因为完全括号化的性质，最终都归根于2个矩阵的乘积，所以我们从2个矩阵起步，
    逐步向上计算，最终就能计算（0，n）的组合的最优括号化的代价。
*/
    for(int l=2;l<=length;l++)
    {
        for(int i=0;i<=length-l;i++)            //此处要注意的i+l-1不能超过数组的个数
        {
            int j = i+l-1;                      //矩阵个数要为l，则j-i需满足j-i+1=l，所以有该等式
            int q = 0;
            for(int k=i;k<j;k++)                //注意k要在区间内移动，利用下述if条件，获得所有可能值中的最小值，
            {
                q = cost[i][k] + cost[k+1][j] + num[i]*num[k+1]*num[j+1];
                if(q < cost[i][j])
                {
                    cost[i][j] = q;
                    point[i][j] = k;
                }
            }
        }
    }

}

void show(int **p)
{
    using namespace std;
    for(int i=0;i<SIZE-1;i++)
    {
        for(int j=0;j<SIZE-1;j++)
        {
            cout << p[i][j] << '	';
        }
        cout << endl;
    }
}

void print(int** point,int i,int j)
{
    if(i == j)
    {
        std::cout << "A" << i+1;
    }
    else
    {
        std::cout << "(";
        print(point,i,point[i][j]);             //输出划分点左侧的完全括号化表达式
        print(point,point[i][j]+1,j);            //输出划分点右侧的完全括号化表达式
        std::cout << ")";
    }
}
/*矩阵链乘的方法与括号化方案表达式输出的方法是同样的道理*/
matrix matrix_mutiply(int** point,int i,int j,matrix* m)
{
    if(i==j)
        return m[i];
    else
    {
        matrix m1 = matrix_mutiply(point,i,point[i][j],m);
        matrix m2 = matrix_mutiply(point,point[i][j]+1,j,m);
        return m1*m2;
    }
}

该问题的解决代码在上述三个文件中。

同时上述代码也存在一定的问题，关于内存是否泄漏以及代码是否符合规范，本人不是非常熟悉。望各位大神，看后觉得有什么问题，直接回复，我也希望提高，毕竟code能力有待提高。

以上问题也属于典型的动态规划问题，也属于递归向迭代的转换。

总结：

应用动态规划的方法的步骤：

1、刻画一个最优解的结构特征

2、递归地定义最优解的值

3、计算最优解的值，通常采用自底向上的方法

4、利用计算出的信息构造一个最优解的值

源代码文件。。。。。。。点击这里