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

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

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

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

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

问题描述: 给定n个矩阵{A1,A2,...,An},其中Ai与Ai+1是可乘的,i=1,2...,n-1。如何确定计算矩阵连乘积的计算次序,使得依此次序计算矩阵连乘积需要的数乘次数最少。

显而易见,在矩阵个数为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 矩阵 )

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

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

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

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

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

  还拿上面的那个例子来讲: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.然后按照排序顺序(从大到小)进行矩阵乘法。

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

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

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

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

 

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

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

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

  好了,思考暂时便是这么多了,我想也是时候休息一下了(譬如玩玩《KOF》或者《SF4》) :)

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏专知

【论文推荐】最新八篇情感分析相关论文—注意力网络、多模态情感分析、情感分析局限性、跨语言情感分类、多语言情感分析

【导读】专知内容组今天推出最新八篇情感分析(Sentiment Analysis)相关论文,欢迎查看!

79920
来自专栏python读书笔记

《算法图解》NOTE 1-算法的渐近表示法以及二分法1 .渐近表示法2.二分法

15860
来自专栏专知

【论文推荐】最新十篇机器翻译相关论文—记忆增强、词嵌入、形态学信息、非自回归神经序列、用户反馈、语法纠错、单语语料、自注意力机制

26520
来自专栏SeanCheney的专栏

《利用Python进行数据分析·第2版》第13章 Python建模库介绍13.1 pandas与模型代码的接口13.2 用Patsy创建模型描述13.3 statsmodels介绍13.4 sciki

本书中,我已经介绍了Python数据分析的编程基础。因为数据分析师和科学家总是在数据规整和准备上花费大量时间,这本书的重点在于掌握这些功能。 开发模型选用什么库...

1.1K60
来自专栏Java与Android技术栈

常用的像素操作算法:图像加法、像素混合、提取图像中的ROI

图像可以是看成是一个多维的数组。读取一张图片,可以看成是读入了一系列的像素内容。这些像素内容,按照不同的模式具有不同的格式。对于三通道的 RGB 位图来说,每个...

9320
来自专栏老秦求学

从Iris数据集开始---机器学习入门

代码多来自《Introduction to Machine Learning with Python》. 该文集主要是自己的一个阅读笔记以及一些小思考,小总结...

790100
来自专栏AI科技大本营的专栏

如何快速优化机器学习的模型参数

【导读】一般来说机器学习模型的优化没什么捷径可循。用什么架构,选择什么优化算法和参数既取决于我们对数据集的理解,也要不断地试错和修正。所以快速构建和测试模型的能...

12520
来自专栏专知

【论文推荐】最新十篇机器翻译相关论文—自适应机器翻译综述、结构化预测、双向神经机器翻译、图编解码模型、英日机器翻、上下文感知

【导读】专知内容组既昨天推出九篇机器翻译(Machine Translation)相关论文,今天又推出最新十篇机器翻译相关论文,欢迎查看!

21130
来自专栏量化投资与机器学习

从Seq2seq到Attention模型到Self Attention(二)

系列一介绍了Seq2seq和 Attention model。这篇文章将重点摆在Google於2017年发表论文“Attention is all you ne...

67650
来自专栏贾志刚-OpenCV学堂

OpenCV中实现曲线与圆拟合

使用OpenCV做图像处理与分析的时候,经常会遇到需要进行曲线拟合与圆拟合的场景,很多OpenCV开发者对此却是一筹莫展,其实OpenCV中是有现成的函数来实现...

37340

扫码关注云+社区

领取腾讯云代金券