动态规划概论

头像
fogsail
    阅读 18 分钟

    动态规划(dynamic programming)是运筹学的一个分支，是求解决策过程(decision process)最优化的数学方法。20世纪50年代初美国数学家R.E.Bellman等人在研究多阶段决策过程(multistep decision process)的优化问题时，提出了著名的最优化原理(principle of optimality)，把多阶段过程转化为一系列单阶段问题，利用各阶段之间的关系，逐个求解，创立了解决这类过程优化问题的新方法——动态规划。

    钢条切割问题

    算法运行时间分析

    如果采用最平凡的递推算法，由$r_n= underset{1 leq i leq n}{max} (p_i+r_{n-i})$
    令$T(n)$表示算法第二个参数值为n的时候CUT-ROD的调用次数，可以推出：$T(n)=1+ sum_{j=0}^{n-1}T(j)$

    $T(n)=1+T(1)+T(2)+cdots+T(n-1)$
    $T(n-1)=1+T(1)+T(2)+cdots+T(n-2)$

    即$T(n)=2+2T(1)+cdots+2T(n-2)=2T(n-1)$

    $frac{T(n)}{T(n-1)}=2$

    由等比数列的递推公式，$T(n)=2^n$
    所以，平凡的递归算法反复求解相同的子问题，耗费了大量的时间。

    使用动态规划

    图片描述
    求解子问题的两种方法：
    使用memorized_cut_rod方法的时候，memorized_cut_rod_aux(p,n,r)的递归，最多递归到r[]的下标为n-1，因为在循环中执行的是for(int i=0;i<n;i++)

    图片描述
    使用自底向上递归求解。

    图片描述
    由上图所示，在递归过程中：

    r[0]=0;
for(int j=0;j<n;j++)
{
  int result=-INFINITY;
  for(int i=0;i<=j;i++)
  {
    //获取scale=j时候的max值，即最大收益result
    //随着j++，我们获取最大收益的时候要返回r[j+1]的值，这样
    //scale=scale+1时候的最大收益，即r[j+1]，即为上一步求出的
    //result值，即r[j+1]=result
    result=max(result,p[i]+r[j-1]);
  }
  r[j+1]=result;

}
return r[n];
//返回的r[n]这里的n指问题的规模

    算法实现过程：

    #include <iostream>
#define INFINITY 0x7fffffff
using namespace std;

int max(int a,int b)
{
  return a>b?a:b;
}

int memorized_cut_rod_aux(int p[],int n,int r)  //n表示问题的规模
{
  //注意，这里p[i]的意义和r[i]的意义不一样
  //p[i]表示对应的数组元素，i表示下标
  //r[i]表示问题的规模为i所对应的解
  //p[i]对应的问题是第i个，其规模为i+1，剩余的问题规模是n-(i+1)
  //从另一个角度，0<=i<=n-1 0<=r<=n-1，所以递归的时候是(p,n-1-i,r)
  int result;
  if(r[n]>=0)
    return r[n];  //表示问题已经被求解完毕了，返回该值

  if(n==0)
    result=0;  //表示没有切割
  else
  {
    result=-INFINITY;
    for(int i=0;i<n;i++)
    {
      result=max(result,p[i]+memorized_cut_rod_aux(p,n-1-i,r));
    }
  }
  //表示规模为n的问题求解完毕
  r[n]=result;
  return result;

}
int memorized_cut_rod(int p[],int n)
{
  int *r=new int[n];
  for(int i=0;i<n;i++)
  {
    r[i]=-INFINITY;
  }

  return memorized_cut_rod_aux(p,n,r);  //这里(p,n,r)
  //r[]的递归标最多为n-1
}

int BOTTOM_UP_CUT_ROD(int p[],int n)
{
  int *r=new int [n];
  r[0]=0;
  for(int j=0;j<n;j++)
  {
    int result=-INFINITY;
    for(int i=0;i<=j;i++)
    {
      //获取scale=j时候的max值，即最大收益result
      //随着j++，我们获取最大收益的时候要返回r[j+1]的值，这样
      //scale=scale+1时候的最大收益，即r[j+1]，即为上一步求出的
      //result值，即r[j+1]=result
      result=max(result,p[i]+r[j-1]);
    }
    r[j+1]=result;

  }
  return r[n];
  //返回的r[n]这里的n指问题的规模
}

    主函数：

    int main()
{
  const int n=10;
  int p[10]={1,5,8,9,10,17,17,20,24,30};
  cout<<memorized_cut_rod(p,9)<<endl;
  cout<<BOTTOM_UP_CUT_ROD(p,9)<<endl;
}

    子问题图与重构解

    这里特别说明一种表示方式

    array *EXTENDED_BOTTOM_UP_CUT_ROD(int p[],int n)

    这里p[]表示指向数组第一个元素的指针，注意区分以下两种表示方法：当表示第n个元素的时候

    *(p+n)->element
p[n].element

    这两种表示方式不同，一种是指针，一种是结构体表示。

    具体的说明见下图：

    图片描述
    动态规划解决钢条切割问题，问题的输出：注意让规模逐渐减小。

    #include <iostream>
using namespace std;
#define INFINITY 0x7fffffff

struct rod
{
  int result;   //代表最大收益
  int solution;  //代表切割方案
};

rod *EXTENDED_BOTTOM_UP_CUT_ROD(int p[], int n)
{
  rod *res=new rod[n];
  res[0].result=0;  //初始化，钢条规模为0，收益为0
  int q;
  for(int j=0;j<n;j++)
  {
    q=-INFINITY;
    for(int i=0;i<=j;i++)
    {
      if(q<p[i]+res[j-i].result)
      {
        q=p[i]+res[j-i].result;
        res[j+1].solution=i+1;  //保存切割方式
      }        //这里的i表示切割下来的片段，规模有多大？
    }
    //此时已经完成了从1---j的循环
    res[j+1].result=q;
  }
  return res+n;  //注意：这里res表示指向第一个元素的指针
}
//res+n表示指针移动n位，新的指针指向n之后的那个位置
//rod *res=new rod[n] res[i].result  (res+i)->result
//两种用法注意区分

//特别注意：这里res[0].result=0  最后res[j+1].result
//res[j+1].solution，意味这res[]的存储范围如下：
//res[1]，res[2],......,res[n]
//res初始指针位置为res[0]，res末指针的位置为res[n]

void PRINT_CUT_ROD_SOLUTION(int p[], int n)
{
  rod *res=EXTENDED_BOTTOM_UP_CUT_ROD(p,n);
  //注意这里函数的返回值为res+n，即最后一个数组元素的地址
  while(n>0)
  {
    cout<<(*res).solution<<" ";
    n=n-(*res).solution;
    //注意返回的res指的是res[n]的位置
    //返回的是规模为n时候的切割方法,从最大规模往最小规模输出
    //这里从大规模往小规模输出，指针也要调整
    //从规模为n的位置逐渐减小到i，再逐渐减小到1，0
    res=res-(*res).solution;

  }
}

void main()
{
  const int n=10;
  int p[10]={1,5,8,9,10,17,17,20,24,30};
  cout<<(*EXTENDED_BOTTOM_UP_CUT_ROD(p,4)).result<<endl;
  PRINT_CUT_ROD_SOLUTION(p,4);
}

    其他问题

    15.1-2

    反例如下图：

    图片描述
    15.1-3 钢条切割的时候考虑固定成本为c

    特别注意，这里的固定成本为c，在执行循环的时候，存在表达式的不同！

    //当j==i的时候有不同！
for(int j=0;j<n;j++)
{
  //省略的code
  for(int i=0;i<j;i++)
  {
    //i<j时候表示钢条一定会进行切割！
    if(q<=p[i]+res[j-i].result-c)
    {
      q=p[i]+res[j-i].result-c;
      res[j+1].solution=i+1;
      cut_or_not=true;  //表示进行切割了
    }
  }

  //当i==j时候的情况不一样了！
  //此时已经保存的最大收益是q
  //很有可能
  //p[i]+res[j-i].result-c<q<=p[i]+res[j-i].result
  //这个时候j==i可能没有发生切割，没有发生切割和发生切割
  //收益的表达式不同，一个减去成本，一个没有
  if(j==i)
  {
    if(q<=p[i]+res[j-i].result && cut_or_not==false)
    {
      res[j+1].solution=i+1;
    }
  }

  //j求解完毕
  res[j+1].result=q;

}

return res+n;

    具体实现的过程如下：

    #include <iostream>
using namespace std;
#define INFINITY -0x7fffffff
#define cost 2

int max(int a, int b)
{
  return a>b?a:b;
}

struct rod_cost
{
  int result;
  int solution;
  bool cut_or_not;
};

rod_cost *EXTENDED_BOTTOM_UP_CUT_ROD(int p[], int n)
{
  rod_cost *res=new rod_cost[n];
  res[0].result=0;
  int q;

  for(int j=0;j<n;j++)
  {
    q=-INFINITY;
    bool cut_or_not=false;
    for(int i=0;i<j;i++)
    {
      if(q<=p[i]+res[j-i].result-cost)
      {
        q=p[i]+res[j-i].result-cost;
        res[j+1].solution=i+1;
        cut_or_not=true;
      }
    }

    //每一次循环，j=0----j=n-1，表示问题的规模从1---n
    //每一个规模，cut_or_not都初始化为false，最后看是否切割？
    //如果发生了切割，在i的循环中，也就是j的子问题中
    //让cut_or_not=true，最后判断是否切割？

    if(j==i)
    {
      if(q<=p[i]+res[j-i].result && cut_or_not==false)
      {
        res[j+1].solution=i+1;
      }
    }

    //求解完毕
    res[j+1].result=q;
  }

  return res+n;
}

void PRINT_CUT_ROD_SOLUTION(int p[], int n)
{
  rod_cost *res=EXTENDED_BOTTOM_UP_CUT_ROD(p,n);
  while(n)
  {
    cout<<(*res).solution<<" "；
    n=n-(*res).solution;
    res=res-(*res).solution;
  }
}

void main()
{
  const int n=10;
  int p[10]={1,5,8,9,10,17,17,20,24,30};
  cout<<(*EXTENDED_BOTTOM_UP_CUT_ROD(p,10)).result<<endl;
  PRINT_CUT_ROD_SOLUTION(p,10);
}

    15.1-4修改memorized_cut_rod，使之返回最优收益值和切割方案

    实现函数如下：

    #include <iostream>
using namespace std;
#define INFINITY 0x7fffffff

struct rod
{
  int result;
  int solution;
};

int max(int a, int b)
{
  return a>b?a:b;
}

rod *memorized_cut_rod_aux(int p[], int n, rod res[])
{
  int q;
  if(res[n].result>=0)
    return res+n;

  if(n==0)
    q=0;   //规模为0的时候，自然收益为0
  else
  {
    q=-INFINITY;
    for(int i=0;i<n;i++)
    {
      //会犯错误的地方：
      //这里不能够直接写
      //q=max(q,p[i]+(memorized_cut_rod(p,n-1-i,res))->solution)
      //因为还要进行储存切割方案，在子问题规模为i+1的时候，有可能进行切割，也有可能不切割
      int tmp=p[i]+memorized_cut_rod_aux(p,n-1-i,res)->solution;
      //在p[i]处的最优方案是tmp，tmp和q进行比较判断后，推出在i+1处的最优方案

      //这里把n分解成更小的子问题，怎么解决？通过递归解决
      //memorized_cut_rod_aux(p,n-1-i,res)把n-1的子问题转化为n-1-i
      if(q<tmp)
      {
        q=tmp;
        res[n].solution=i+1;
        //我们这里返回的就是问题规模为n时候的切割方案，这里的切割方案是：
        //切割一次or无切割！
        //i<n-1时候，q=tmp表示切割
        //i==n-1的时候，tmp=p[n-1]+memorized_cut_rod(p,0,res)->solution
        //如果此时还满足q<tmp，res[n].solution=n-1+1=n，表示无切割
        //这一步的目的在于保存每一个循环过程中i的最优方案
      }
    }
  }

  res[n].result=q;
  return res+n;
}

rod *memorized_cut_rod(int p[], int n)
{
  rod *res=new rod[n];
  for(int i=0;i<=n;i++)
  {
    res[i].result=-INFINITY;
  }

  return memorized_cut_rod_aux(p,n,res);
}

void PRINT_CUT_ROD_SOLUTION(int p[], int n)
{
  rod *res=memorized_cut_rod(p,n);
  cout<<"max result"<<res->result<<endl;
  cout<<"solution: "<<endl;
  while(n)
  {
    cout<<(*res).solution<<" ";
    n=n-(*res).solution;
    res=res-(*res).solution;
  }
}

//重点：在于区分数组下标和子问题的规模，数组下标为i-1，子问题的规模为i

void main()
{
  const int n=10;
  int p[10]={1,5,8,9,10,17,17,20,24,30};
  PRINT_CUT_ROD_SOLUTION(p,10);
  cout<<endl;
}

    15.1-5斐波那契数列

    子问题图：

    图片描述
    实现方法：

    #include <iostream>
using namespace std;

int fibonacci(int array[],int n)
{
  if(n>=0)
    array[0]=0;
  if(n>=1)
    array[1]=1;

  for(int i=2;i<n;i++)
  {
    array[i]=array[i-1]+array[i-2];
  }
  return array[n];
}

    特别说明：

    for(int j=0;j<n;j++)
{
  //code
  for(int i=0;i<=j;i++)
  {
    q=max(q,p[i]+r[j-i]);
  }
  r[j]=q;     //这样写也是可\ \;\;以的
}

return r[n-1];  //只是返回值会不一样

    矩阵链乘法

    矩阵规模序列

    $A1 qquad ;;|$$A2 \qquad\ \;\;|$$A3 qquad ;;|$
    $underline{10}times100,|$$100\times5\;\;|$$underline{5times50}quad|$

    $A1A2 qquadtimes A3=10times5times50$

    由此可以推出：
    $A_{i ldots k} qquad qquad times qquad qquad A_{k+1 ldots j}$
    $=underline{p_{i-1}} times p_i times p_{i+1} times qquadquad underline{p_k} times p_{k+1}timesunderline{p_j}$
    $=p_{i-1} times p_k times p_j$

    由矩阵规模序列推出矩阵链乘法，如上图所示，相乘的代价是：把有下划线的项相乘的值。
    可以得到递推式如下：

    $$m[i,j]=
    begin{cases}
    0& text{i=j}\
    min(m[i,k]+m[k+1,j]+p_{i-1}p_kp_j)& text{i<j}
    end{cases}$$

    矩阵链乘法的运算过程

    如下图所示：

    图片描述
    可以看出，三角形处的min值取决于两个圆圈处的值，具体的实现过程如下：
    line01:

    $$m[1,2]=
    begin{cases}
    m[1,1]+m[2,2]+p_0p_1p_2\
    =0+0+30times35times15
    =15750
    end{cases}$$

    $$m[2,3]=
    begin{cases}
    m[2,2]+m[3,3]+p_1p_2p_3\
    =0+0+35times15times5
    =2625
    end{cases}$$

    $$m[3,4]=
    begin{cases}
    m[3,3]+m[4,4]+p_3p_4\
    =0+0+15times5times10
    =750
    end{cases}$$

    $$m[4,5]=
    begin{cases}
    m[4,4]+m[5,5]+p_4p_5\
    =0+0+5times10times20
    =1000
    end{cases}$$

    $$m[5,6]=
    begin{cases}
    m[5,5]+m[6,6]+p_5p_6\
    =0+0+10times20times25
    =5000
    end{cases}$$

    line02：括号中表示划分方式

    $$m[1,3]=min
    begin{cases}
    0([1,1])+2625([2,3])+30times35times5=7875\
    15750([1,2])+0+30times15times5=18000
    end{cases}$$

    $$m[2,4]=min
    begin{cases}
    0([2,2])+750([3,4])+35times15times10=6000\
    2625([2,3])+0+35times5times10=4375
    end{cases}$$

    $$m[3,5]=min
    begin{cases}
    0([3,3])+1000([4,5])+15times5times20=2500\
    750([3,4])+0+15times10times20=3750
    end{cases}$$

    $$m[4,6]=min
    begin{cases}
    0([4,4])+5000([5,6])+5times10times25=6250\
    1000([4,5])+0+5times20times25=3500
    end{cases}$$

    line03：

    $$m[1,4]=min
    begin{cases}
    0([1,1])+4375([2,4])+30times35times10=14875\
    15750([1,2])+750([3,4])+30times15times10=4500+750+15750=21000\
    7875([1,3])+0([4,4])+30times5times10=9375
    end{cases}$$

    $$m[3,6]=min
    begin{cases}
    0([3,3])+3500([4,6])+15times5times25=5375\
    750([3,4])+5000([5,6])+15times10times25=9500\
    2500([3,5])+0([6,6])+15times20times25=10000
    end{cases}$$

    line04：

    $$m[1,5]=min
    begin{cases}
    0([1,1])+7125([2,5])+30times35times20=28125\
    15750([1,2])+2500([3,5])+30times15times20=27250\
    7875([1,3])+1000([4,5])+30times5times20=11875\
    9375([1,4])+0([5,5])+30times10times20=15375
    end{cases}$$

    $$m[2,6]=min
    begin{cases}
    0([2,2])+5375([3,6])+35times15times25=18500\
    2625([2,3])+3500([4,6])+35times5times25=10500\
    4375([2,4])+5000([5,6])+35times10times25=18125\
    7125([2,5])+0([6,6])+35times20times25=24625
    end{cases}$$

    line05：

    $$m[1,6]=min
    begin{cases}
    0([1,1])+10500([2,6])+30times35times25=36750\
    15750([1,2])+5375([3,6])+30times15times25=32375\
    7875([1,3])+3500([4,6])+30times5times25=15125\
    9375([1,4])+5000([5,6])+30times10times25=21875\
    11875([1,5])+0([6,6])+30times20times25=26875
    end{cases}$$

    计算最优代价

    图片描述
    具体的实现过程如下：

    #include <iostream>
#include <vector>
#include <stdlib.h>
#include <ctime>
using namespace std;
#define n 7
#define INFINITY 0x7fffffff

struct Matrix
{
    int rows;  //表示行数
    int columns;  //表示列数
    vector<vector<int> > matrix;
    Matrix(int r,int c)
    {
        rows=r;
        columns=c;
        matrix.resize(rows);
        for(int i=0;i<rows;i++)
          matrix[i].resize(columns);
    }
};

struct Matrix_Chain
{
    vector<vector<int> > m;  //表示运算次数
    vector<vector<int> > s;  //划分方式
    Matrix_Chain()
    {
        m.resize(n-1);  //有n-1个括号
        for(int k=0;k<n-1;k++)
        {
            m[k].resize(n-1);    //最后输出的是一个对角阵
                                 //辅助表m[1...n-1, 1...n-1]来保存代价
                                 //s[1...n-1,2...n]记录最优值m[i,j]对应的分割点k
        }
        s.resize(n-1);
        for(int t=0;t<n-1;t++)
        {
            s[t].resize(n-1);
        }
    }
};

Matrix init(Matrix A)
{
    for(int i=0;i<A.rows;i++)
    {
        for(int j=0;j<A.columns;j++)
        {
            A.matrix[i][j]=rand()%10;
        }
    }
    return A;
}

Matrix Matrix_Mutiply(Matrix A,Matrix B)
{
    if(A.columns!=B.rows)
    {
        Matrix D(1,1);
        D.matrix[0][0]=0;
        cerr<<"incompatible dimensions"<<endl;
        return D;
    }
    else
    {
        Matrix C(A.rows, B.columns);
        for(int i=0;i<A.rows;i++)
        {
            for(int j=0;j<B.columns;j++)
            {
                C.matrix[i][j]=0;
                for(int k=0;k<A.columns;k++)
                   C.matrix[i][j]=C.matrix[i][j]+A.matrix[i][k]*B.matrix[k][j];
            }
        }
        return C;
    }
}

Matrix_Chain Matrix_Chain_Order(int p[])
{
    Matrix_Chain T;
    int len=n-1;    //辅助表m[k,k]来保存代价，其中k=n-1，n为矩阵链长度
    for(int i=0;i<len;i++) //s[k,k]用来保存最优值m[i,j]对应的分割点
    {
        T.m[i][i]=0;
    }
    for(int l=2;l<=len;l++)  //l表示矩阵链长度
    {
        for(int i=1;i<=len-l+1;i++)
        {
            int j=i+l-1;
            T.m[i-1][j-1]=INFINITY;   //用来计算[i,i+1..k] [k+1..,j-1,j]这样划分的运算代价
            for(int k=i;k<j;k++)     //循环中的i表示第几个数，T.m[i]中的i表示数组下标
            //千万注意，这里k从i开始！
            {
                int q=T.m[i-1][k-1]+T.m[k][j-1]+p[i-1]*p[k]*p[j];
                if(q<T.m[i-1][j-1])
                {
                    T.m[i-1][j-1]=q;
                    T.s[i-1][j-1]=k-1;
                }
            }
        }
    }
    return T;
}

void print(Matrix A)
{
    for(int i=0;i<A.rows;i++)
    {
         for(int j=0;j<A.columns;j++)
            cout<<A.matrix[i][j]<<" ";
            cout<<endl;
    }
    cout<<endl;
}

int main()
{
    int p[n]={5,10,3,12,5,50,6};
    Matrix_Chain T=Matrix_Chain_Order(p);
    for(int i=0;i<n-1;i++)
    {
        for(int j=0;j<n-1;j++)
        {
            if(T.m[i][j]<0)
              cout<<"-1"<<"\t";
            else cout<<T.m[i][j]<<"\t";
        }
        cout<<endl;
    }
}

    运算结果：

    图片描述

    矩阵链总乘法

    矩阵链总乘法的递归表示：

    图片描述
    递归运算的核心：递归进行到最底层的时候，算法的运行情况？
    这里的运算是矩阵乘法：

    if(j==i+1)
{
  return Matrix_Mutiply(A[i],A[j]);
}

    相当于二路归并排序，递归进行到最末端时候的情况。A[i]*A[i+1]
    在每一次的递归过程中，划分的下标信息，是存储在S[1,...,n]中的，我们在函数中要传递该参数。

    具体实现过程：

    #include <iostream>
#include <vector>
#include <stdlib.h>
#include <ctime>
#include <memory>
using namespace std;
#define n 7
#define INFINITY 0x7fffffff

struct Matrix
{
    int rows;  //表示行数
    int columns;  //表示列数
    vector<vector<int> > matrix;
};

struct Matrix_Chain
{
    vector<vector<int> > m;  //表示运算次数
    vector<vector<int> > s;  //划分方式
    Matrix_Chain()
    {
        m.resize(n-1);  //有n-1个括号
        for(int k=0;k<n-1;k++)
        {
            m[k].resize(n-1);    //最后输出的是一个对角阵
                                 //辅助表m[1...n-1, 1...n-1]来保存代价
                                 //s[1...n-1,2...n]记录最优值m[i,j]对应的分割点k
        }
        s.resize(n-1);
        for(int t=0;t<n-1;t++)
        {
            s[t].resize(n-1);
        }
    }
};

Matrix init(Matrix &A,int r,int c)  //对A的矩阵具体值的初始化
{
    A.rows=r;
    A.columns=c;
    A.matrix.resize(A.rows);
    for(int i=0;i<A.rows;i++)
        A.matrix[i].resize(A.columns);

    srand((unsigned)time(0));
    for(int i=0;i<A.rows;i++)
    {
        for(int j=0;j<A.columns;j++)
        {
            A.matrix[i][j]=rand()%n+1;
        }
    }
    return A;
}

void print(Matrix A)
{
    for(int i=0;i<A.rows;i++)
    {
         for(int j=0;j<A.columns;j++)
            cout<<A.matrix[i][j]<<" ";
            cout<<endl;
    }
    cout<<endl;
}

Matrix Matrix_Mutiply(Matrix A,Matrix B)
{
    if(A.columns!=B.rows)
    {
        Matrix D;
        init(D,1,1);
        D.matrix[0][0]=0;
        cerr<<"incompatible dimensions"<<endl;
        return D;
    }
    else
    {
        Matrix C;
        init(C,A.rows,B.columns);
        for(int i=0;i<A.rows;i++)
        {
            for(int j=0;j<B.columns;j++)
            {
                C.matrix[i][j]=0;
                for(int k=0;k<A.columns;k++)
                   C.matrix[i][j]=C.matrix[i][j]+A.matrix[i][k]*B.matrix[k][j];
            }
        }
        return C;
    }
}

Matrix_Chain Matrix_Chain_Order(int p[])  //这里的i j统一表示是：数组中的第几个数？
{                                         //所以，数组下标是T.m[i-1][j-1]
    Matrix_Chain T;
    int len=n-1;    //辅助表m[k,k]来保存代价，其中k=n-1，n为矩阵链长度
    for(int i=0;i<len;i++) //s[k,k]用来保存最优值m[i,j]对应的分割点
    {
        T.m[i][i]=0;
    }
    for(int l=2;l<=len;l++)  //l表示矩阵链长度
    {
        for(int i=1;i<=len-l+1;i++)
        {
            int j=i+l-1;
            T.m[i-1][j-1]=INFINITY;   //用来计算[i,i+1..k] [k+1..,j-1,j]这样划分的运算代价
            for(int k=i;k<j;k++)     //循环中的i表示第几个数，T.m[i]中的i表示数组下标
            //千万注意，这里k从i开始！
            {
                int q=T.m[i-1][k-1]+T.m[k][j-1]+p[i-1]*p[k]*p[j];
                if(q<T.m[i-1][j-1])
                {
                    T.m[i-1][j-1]=q;
                    T.s[i-1][j-1]=k-1;
                }
            }
        }
    }
    return T;
}


//矩阵链总乘法，相当于二路归并
//在二路归并排序中，递归处理，让区间的长度从l逐渐减小到0
//递归到尾端的时候，是两个数相乘A[i]*A[j] (j==i+1)

//函数参数：矩阵链乘法中，1、传递A[]矩阵的原始数据，2、括号的划分情况，存储在数组T.s[]中
Matrix Matrix_Chain_Multiply(Matrix A[],Matrix_Chain T,int i,int j)
{
    if(j==i)
      return A[i];
    if(j==i+1)
      return Matrix_Mutiply(A[i],A[j]);
    Matrix t1=Matrix_Chain_Multiply(A,T,i,T.s[i][j]);  //矩阵T.s[i][j]表示A[i...j]的划分方式
    Matrix t2=Matrix_Chain_Multiply(A,T,T.s[i][j]+1,j);
    return Matrix_Mutiply(t1,t2);
}

//同理，打印输出矩阵链乘法也是如此：
void Print_Optimal_Parents(Matrix_Chain T,int i,int j)
{
    if(i==j)
      cout<<"A"<<i<<" ";
    else
    {
        cout<<"("<<" ";
        Print_Optimal_Parents(T,i,T.s[i][j]);
        Print_Optimal_Parents(T,T.s[i][j]+1,j);
        cout<<")"<<" ";
    }
}

int main()
{
    Matrix test;
    init(test,4,5);
    print(test);
    cout<<"Program begins: "<<endl;

    int p[n]={5,10,3,12,5,50,6};
    Matrix_Chain T=Matrix_Chain_Order(p);

    for(int i=0;i<n-1;i++)   //矩阵的dimension，等于矩阵链规模n-1
    {
        for(int j=0;j<n-1;j++)  //原矩阵对应的值A[i....j]=A[0...n-2]
        {
            if(T.m[i][j]<0)
              cout<<"-1"<<"\t";
            else
              cout<<T.m[i][j]<<"\t";
        }
        cout<<endl;
    }

    Print_Optimal_Parents(T,0,n-2);
    cout<<endl;

    cout<<"Concret Implement:"<<" "<<endl;

    Matrix A[n];   //这里的A存放的是指向Matrix数组第一个元素的指针
    //A[]的dimension的值,具体一点为p(i-1)*p(i)
    for(int j=1;j<n;j++)
    {
        Matrix t;
        init(t,p[j-1],p[j]);    //A[rows*columns]--->p[j-1]*p[j]
        A[j-1]=t;    //p[]仅仅表示A[]的维度，具体的值用init初始化
        cout<<endl;
        print(A[j-1]);
        cout<<endl;
    }

    Matrix result=Matrix_Chain_Multiply(A,T,0,n-2);
    print(result);

    return 0;
}

    算法运行情况：

    图片描述
    图片描述

    问题解答

    15.2-3 用代入法证明递归公式(15.6)的结果是：$Omega(2^n)$
    $n=2qquadquad P(2)=p^2(1)=1$
    $n=3qquadquad P(3)=2P(1)P(2)$
    $n=4qquadquad P(4)=2^2P^2(1)P(2)+P^2(2)$
    $n=5qquadquad P(5)=2^3P^3(1)P(2)+(2+2^2)P(1)P^2(2)$
    $cdotscdots$

    $P(n)=2^{[n/2]+1}P^{[n/2]+1}(1)P(2)+(2+2^2+cdots+2^{[n/2]})P(1)P^2(2)$
    由等比数列的求和公式，可以推出：
    $P(n)=Omega(2^n)$

    15.2-4 对输入链长度为n的矩阵链乘法问题，包含多少个顶点？包含多少条边？这些边分别连接哪些顶点？
    值得注意的是，矩阵链不等于二叉树，子问题就对应矩阵中的每个点。

    如下图所示：

    图片描述
    子问题图包含$n^2/2=O(n^2)$个顶点。
    边数？第i层的边数为$i(i-1)$

    总边数为$sum_{i=1}^{i=n}i(i-1)=O(n^3)$

    15.2-5 计算其他表项的时候，访问m[i,j]，主要的时间代价是由
    $q=m[i,k]+m[k+1,j]+p_{i-1}p_kp_j$ 产生的，每一次访问都对m[]调用2次。

    总的遍历循环次数，为k可以取到的值的总数，即为矩阵链长度length
    令$l=length$
    $sum_{i=1}^nsum_{j=1}^nR(i,j)=sum_{l=2}^{l=n}sum_{i=1}^{n-l+1}2(l-1)=2sum_{l=2}^n(n-(l-1))(l-1)=2sum_{t=1}^{n-1}(n-t)t$
    $=O(2frac{n(n-1)(2n-1)}{6})=frac{n^3-n}{3}$

    其实本质上还是这张图：

    图片描述
    图中红色线段长度表示$l$，即$sum_{l=1}{l=n}$，每一个红色线段对应两个m[]点，在每一层，有$l-1$个三角形的点
    所以，总代价是：
    $sum_{l=2}^{l=n}2 times l(l-1)=2 times frac{n(n-1)(2n-1)}{6}=O(frac{n^3-n}{3})$

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

    数学归纳法，k个元素要k-1对括号，当元素增加到k+1时，保持k-1对括号不变，再添加新的一个元素，括号数为k

    动态规划算法
    fogsail
    31 声望6 粉丝

    一直认为数学、音乐、诗歌都是世间最美好的语言，程序是用数学这种语言写成的逻辑。于2014年12月被人拉到兰州大学开源社区中去（mirror.lzu.edu.cn），从此与自由软件和编程解下不解之缘，毅然抛弃生物专业从事编...

    « 上一篇
    在c++中选择适当的操作
    下一篇 »
    最长公共子序列问题

    引用和评论

    推荐阅读
    头像
    投资策略规划

    fogsail赞 1阅读 3.1k

    头像
    基础思想总结

    zxl20070701赞 1阅读 932

    头像
    Python用CNN - LSTM、ARIMA、Prophet股票价格预测的研究与分析|附数据代码

    拓端tecdat阅读 1.5k

    头像
    Python深度学习股价预测、量化交易策略：LSTM、GRU深度门控循环神经网络|附代码数据

    拓端tecdat阅读 1.4k

    头像
    CUDA Toolkit常见安装问题一览

    techlead_kris赞 1阅读 537

    头像
    【专题】AI市场趋势洞察（2024年）报告汇总PDF洞察（附原数据表）

    拓端tecdat阅读 1.1k

    头像
    信息熵系列（一）：信息熵的定义及理解

    ZincO阅读 792

    0 条评论
    得票最新