动态规划之矩阵链乘

问题提出:(代码下载

对于如下矩阵:

clip_image001

 

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

clip_image001[18]

 

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

问题分析:

首先我们需要明确的是何为标量:标量即为没有方向的量,而有方向的量即为矢量。(严谨的定义自己百度去)

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

其次对于两个矩阵相乘,需满足下示公式所示的形式:(左边矩阵的列数与右边矩阵的行数必须一致)

clip_image001

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

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

一个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] << '\t';
        }
        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() != '\n')
    {
        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] << '\t';
        }
        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、利用计算出的信息构造一个最优解的值

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

posted @ 2014-11-06 14:07  天地由我欣  阅读(3234)  评论(0编辑  收藏  举报