【算法】动态规划

最近看了一下动态规划，发现其实还是很容易的，不知道为什么当时上课的时候觉得很难。脑子瓦特了。网上找了一些资料：

幽默的解释，感性的理解动态规划：http://blog.youkuaiyun.com/woshioosm/article/details/7438834

一、基本概念

    动态规划过程是：每次决策依赖于当前状态，又随即引起状态的转移。一个决策序列就是在变化的状态中产生出来的，所以，这种多阶段最优化决策解决问题的过程就称为动态规划。

二、基本思想与策略

    基本思想与分治法类似，也是将待求解的问题分解为若干个子问题（阶段），按顺序求解子阶段，前一子问题的解，为后一子问题的求解提供了有用的信息。在求解任一子问题时，列出各种可能的局部解，通过决策保留那些有可能达到最优的局部解，丢弃其他局部解。依次解决各子问题，最后一个子问题就是初始问题的解。

    由于动态规划解决的问题多数有重叠子问题这个特点，为减少重复计算，对每一个子问题只解一次，将其不同阶段的不同状态保存在一个二维数组中。

    与分治法最大的差别是：适合于用动态规划法求解的问题，经分解后得到的子问题往往不是互相独立的（即下一个子阶段的求解是建立在上一个子阶段的解的基础上，进行进一步的求解）。

 

---

三、适用的情况

能采用动态规划求解的问题的一般要具有3个性质：

    (1) 最优化原理：如果问题的最优解所包含的子问题的解也是最优的，就称该问题具有最优子结构，即满足最优化原理。

    (2) 无后效性：即某阶段状态一旦确定，就不受这个状态以后决策的影响。也就是说，某状态以后的过程不会影响以前的状态，只与当前状态有关。

   （3）有重叠子问题：即子问题之间是不独立的，一个子问题在下一阶段决策中可能被多次使用到。（该性质并不是动态规划适用的必要条件，但是如果没有这条性质，动态规划算法同其他算法相比就不具备优势）

 

---

四、求解的基本步骤

     动态规划所处理的问题是一个多阶段决策问题，一般由初始状态开始，通过对中间阶段决策的选择，达到结束状态。这些决策形成了一个决策序列，同时确定了完成整个过程的一条活动路线(通常是求最优的活动路线)。如图所示。动态规划的设计都有着一定的模式，一般要经历以下几个步骤。

    初始状态→│决策１│→│决策２│→…→│决策ｎ│→结束状态

                      图1 动态规划决策过程示意图

    (1)划分阶段：按照问题的时间或空间特征，把问题分为若干个阶段。在划分阶段时，注意划分后的阶段一定要是有序的或者是可排序的，否则问题就无法求解。

    (2)确定状态和状态变量：将问题发展到各个阶段时所处于的各种客观情况用不同的状态表示出来。当然，状态的选择要满足无后效性。

    (3)确定决策并写出状态转移方程：因为决策和状态转移有着天然的联系，状态转移就是根据上一阶段的状态和决策来导出本阶段的状态。所以如果确定了决策，状态转移方程也就可写出。但事实上常常是反过来做，根据相邻两个阶段的状态之间的关系来确定决策方法和状态转移方程。

    (4)寻找边界条件：给出的状态转移方程是一个递推式，需要一个递推的终止条件或边界条件。

    一般，只要解决问题的阶段、状态和状态转移决策确定了，就可以写出状态转移方程（包括边界条件）。

实际应用中可以按以下几个简化的步骤进行设计：

    （1）分析最优解的性质，并刻画其结构特征。

    （2）递归的定义最优解。

    （3）以自底向上或自顶向下的记忆化方式（备忘录法）计算出最优值

    （4）根据计算最优值时得到的信息，构造问题的最优解

 

---

五、算法实现的说明

    动态规划的主要难点在于理论上的设计，也就是上面4个步骤的确定，一旦设计完成，实现部分就会非常简单。

     使用动态规划求解问题，最重要的就是确定动态规划三要素：

    （1）问题的阶段 （2）每个阶段的状态

    （3）从前一个阶段转化到后一个阶段之间的递推关系。

     递推关系必须是从次小的问题开始到较大的问题之间的转化，从这个角度来说，动态规划往往可以用递归程序来实现，不过因为递推可以充分利用前面保存的子问题的解来减少重复计算，所以对于大规模问题来说，有递归不可比拟的优势，这也是动态规划算法的核心之处。

    确定了动态规划的这三要素，整个求解过程就可以用一个最优决策表来描述，最优决策表是一个二维表，其中行表示决策的阶段，列表示问题状态，表格需要填写的数据一般对应此问题的在某个阶段某个状态下的最优值（如最短路径，最长公共子序列，最大价值等），填表的过程就是根据递推关系，从1行1列开始，以行或者列优先的顺序，依次填写表格，最后根据整个表格的数据通过简单的取舍或者运算求得问题的最优解。

          f(n,m)=max{f(n-1,m), f(n-1,m-w[n])+P(n,m)}

f(n,m)通常表示为一个nxm的数组，保存f(n,m)子问题的最优值，然后自底向上的求解，减少递归带来的重复计算问题。需要注意的是计算顺序的问题。

---

六、动态规划算法基本框架

代码

   

     1
     for(j=1; j<=m; j=j+1) // 第一个阶段
 2    xn[j] = 初始值;
 3 
 4  for(i=n-1; i>=1; i=i-1)// 其他n-1个阶段
 5    for(j=1; j>=f(i); j=j+1)//f(i)与i有关的表达式
 6      xi[j]=j=max（或min）{g(xi-1[j1:j2]), ......, g(xi-1[jk:jk+1])};
 8 
 9 t = g(x1[j1:j2]); // 由子问题的最优解求解整个问题的最优解的方案
10 
11 print(x1[j1]);
12 
13 for(i=2; i<=n-1; i=i+1）
15 {  
17      t = t-xi-1[ji];
18 
19      for(j=1; j>=f(i); j=j+1)
21         if(t=xi[ji])
23              break;
25 }
   

七、练习

7.1 矩阵连乘问题

　一个给定的矩阵序列A1A2...An计算连乘乘积，有不同的结合方法，并且在结合时，矩阵的相对位置不能改变，只能相邻结合。根据矩阵乘法的公式，10*100和100*5的矩阵相乘需要做10*100*5次标量乘法。那么对于维数分别为10*100、100*5、5*50的矩阵A、B、C，用(A*B)*C来计算需要10*100*5 + 10*5*500 =7500次标量乘法；而A*(B*C)则需要100*5*50+10*100*50=75000次标量乘法。

　　那么对于由n个矩阵构成的链<A1,A2,...,An>，对i-1,2,...n，矩阵Ai的维数为pi-1*pi，对乘积A1A2...An求出最小化标量乘法的加括号方案。

数学基础：

1. [n][m]矩阵乘[m][k]矩阵，乘积次数为 n*m*k,乘的结果是一个n x k的矩阵

2. 这个矩阵列，前一个矩阵的列数等于后一个矩阵的行数。Ai x... x Ak 的结果是一个 p[i - 1]p[k]的矩阵。

输入矩阵的数目和一系列可以连乘的矩阵，输出最佳计算次数

样例输入：6 30 35 15 5 10 20 25

样例输出：

先分析转移方程：

其中k是一个在[i, j)之间游走的值，表示从第i个乘到第k个，加上括号，第k+1个乘到第j个，加上括号，这两部分左边乘的结果是一个 p(i-1) x p(k)矩阵，右边的结果是一个p(k) x p(j) 的矩阵，这两个矩阵再乘，又需要 p(i - 1) x p (k) x p(j) 次运算。

代码实现时需要注意的问题：计算顺序！！！

因为你要保证在计算m[i][j]查找m[i][k]和m[k+1][j]的时候，m[i][k]和m[k+1][j]已经计算出来了。

#include<iostream>
using namespace std;
const int MAX = 100;
//p用来记录矩阵的行列，main函数中有说明
//m[i][j]用来记录第i个矩阵至第j个矩阵的最优解
//s[][]用来记录从哪里断开的才可得到该最优解
int p[MAX + 1], m[MAX][MAX], s[MAX][MAX];
int n;//矩阵个数

void matrixChain() {
    // 对角线初始化为0
    for (int i = 1; i <= n; i++)m[i][i] = 0;

    for (int r = 2; r <= n; r++)//对角线循环,r是连乘矩阵的个数
        for (int i = 1; i <= n - r + 1; i++) {//行循环
            int j = r + i - 1;//列的控制
            //找m[i][j]的最小值，先初始化一下，令k=i
            m[i][j] = m[i][i] + m[i + 1][j] + p[i - 1] * p[i] * p[j];
            s[i][j] = i;
            //k从i+1到j-1循环找m[i][j]的最小值
            for (int k = i + 1; k<j; k++) {
                int temp = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
                if (temp<m[i][j]) {
                    m[i][j] = temp;
                    //s[][]用来记录在子序列i-j段中，在k位置处
                    //断开能得到最优解
                    s[i][j] = k;
                }
            }
        }
}

//根据s[][]记录的各个子段的最优解，将其输出
void traceback(int i, int j) {
    if (i == j)return;

    traceback(i, s[i][j]);
    traceback(s[i][j] + 1, j);
    cout << "Multiply A" << i << "," << s[i][j] << "and A" << s[i][j] + 1 << "," << j << endl;
}

int main() {
    cin >> n;
    for (int i = 0; i <= n; i++)cin >> p[i];
    //测试数据可以设为六个矩阵分别为
    //A1[30*35],A2[35*15],A3[15*5],A4[5*10],A5[10*20],A6[20*25]
    //则p[0-6]={30,35,15,5,10,20,25}
    //输入：6 30 35 15 5 10 20 25
    matrixChain();

    traceback(1, n);
    //最终解值为m[1][n];
    cout << m[1][n] << endl;
    return 0;
}

7.2 发工资问题

–问题描述

作为老师，最盼望的日子就是每月的10号了，因为这一天是发工资的日子，养家糊口就靠它了，但是对于学校财务处的工作人员来说，这一天则是很忙碌的一天，财务处的小胡老师最近就在考虑一个问题：如果每个老师的工资额都知道，最少需要准备多少张人民币，才能在给每位老师发工资的时候都不用老师找零呢？这里假设老师的工资都是正整数，单位元，人民币一共有100元、50元、10元、5元、2元和1元六种。 
–输入

输入数据包含多个测试实例，每个测试实例的第一行是一个整数n（n<100），表示老师的人数，然后是n个老师的工资。n=0表示输入的结束，不做处理。 
–输入样例

3
1 2 3
0

–输出样例

4

这个问题的关键在于计算每个老师的工资最少需要多少张人民币，加起来就行了。实际上问题是求给定的工资数，有六种人民币100, 50, 10, 5, 2，最少需要几张加起来为工资额。转移方程也很简单，用f(m) = min(f(m - v[i] +１)其中f(m)表示m钱数时需要的最小值。计算顺序是从低向上，肯定是先计算小值。

#include <iostream>
using namespace std;
#define INF 100000
int main()
{
    int n;
    int t[100];
    int m[6] = { 100, 50, 10, 5, 2, 1 }; 
    int f[1000] = {0}; // 1000块以内的每种钱的最优解

    for (int i = 1; i <= 1000; i++) // 注意计算顺序
    {
        int min = INF;
        for (int j = 0; j < 6; j++) 
        {
            if (i - m[j] >= 0) // 边界
            {
                int count = f[i - m[j]] + 1;
                if (count < min)
                    min = count;
            }
        }
        f[i] = min; 
    }

    while (cin >> n)
    {
        if (n == 0) break;
        int count = 0; 
        for (int i = 0; i < n; i++)
        {
            cin >> t[i];
            count += f[t[i]];
        }
        cout << count;
    }
    getchar();
}

当然，用常规想法也是可以做的，比如对于一定钱数，肯定是先用100的，再用50的，依次下去，反正有1块的来填坑。有贪心的思想。

#include<iostream>
using namespace std;
int f(int n)
{
    int a,b,c,d,e,f;
    a=n/100;
    b=(n-100*a)/50;
    c=(n-100*a-50*b)/10;
    d=(n-100*a-50*b-10*c)/5;
    e=(n-100*a-50*b-10*c-5*d)/2;
    f=n-100*a-50*b-10*c-5*d-2*e;
    return (a+b+c+d+e+f);
}
int main()
{
    int n;
    int a[100];
    while(cin>>n)
    {
        if(n==0) break;
        int t=0;
        for(int i=0;i<n;i++)
        cin>>a[i];
        for(int i=0;i<n;i++)
            t=t+f(a[i]);
            cout<<t<<endl;
    }

    return 0;
}

7.3 买书

有一书店引进了一套书，共有3卷，每卷书定价是60元，书店为了搞促销，推出一个活动，活动如下：

如果单独购买其中一卷，那么可以打9.5折。
如果同时购买两卷不同的，那么可以打9折。
如果同时购买三卷不同的，那么可以打8.5折。

如果小明希望购买第1卷x本，第2卷y本，第3卷z本，那么至少需要多少钱呢？（x、y、z为三个已知整数）。

输入：x, y, z 

输出需要的最少金钱

这个问题最好的买法肯定是一次一次买，一次可以买一本、两本或者三本。

 转移方程 

f(x, y, z) = min ( 

f( x - 1, y, z)  + 60*0.95 , 

f(x, y - 1, z) + 60*0.95, 

f(x, y, z) + 60*0.95, 

f(x - 1, y - 1, z) + 120 * 0.9, 

f(x - 1, y, z- 1) + 120 * 0.9, 

f(x, y - 1, z - 1) + 120 * 0.9, 

f( x - 1,y -1 ,z -1 ) + 180 * 0.85

)

f(x, y, z)表示购买x本1卷，y本2卷，z本3卷的最优解。共有x * y * z 个问题，每个问题又有7中选择，所以时间复杂度是O(x ×ｙ×ｚ×７) => 忽视7 复杂度为O（x×y×z）。

#include <iostream>
using namespace std;
#define MAX 20
#define INF 10000
int f[MAX][MAX][MAX] = { { INF },{ INF },{ INF }}; // 购买x本1卷，y本2卷，z本3卷的最优解

int main()
{
    // 先初始化
    for (int x = 0; x < MAX; x++)
        for (int y = 0; y < MAX; y++)
            for (int z = 0; z < MAX; z++)
                f[x][y][z] = INF; 
    f[0][0][0] = 0;
    f[1][0][0] = 60 * 0.95;
    f[0][1][0] = 60 * 0.95;
    f[0][0][1] = 60 * 0.95;
    for(int x = 0; x < MAX; x++)
        for(int y = 0; y < MAX; y++)
            for (int z = 0; z < MAX; z++)
            {
#define T(X, V)\
do{\
if(X)\
{\
if(V < f[x][y][z])\
    f[x][y][z]= V;\
}\
}while(0)
                T(x - 1 >= 0, f[x - 1][y][z] + 60* 0.95);
                T(y - 1 >= 0, f[x][y -1 ][z] + 60 * 0.95);
                T(z - 1 >= 0, f[x][y][z-1] + 60 * 0.95);
                T(x - 1 >= 0 && y - 1 >=0, f[x - 1][y - 1][z] + 120*0.9);
                T(x - 1 >= 0 && z - 1 >= 0, f[x - 1][y ][z - 1] + 120 * 0.9);
                T(y - 1 >= 0 && z - 1 >= 0, f[x][y - 1][z - 1] + 120 * 0.9);
                T(y - 1 >= 0 && z - 1 >= 0 && x - 1 >= 0, f[x - 1][y - 1][z - 1] + 180 * 0.85);
            }
    for (;;)
    {
        int x, y, z;
        cin >> x >> y >> z; 
        cout << f[x][y][z] << endl;
    }
}

7.4 最长公共子序列Longest Common Subsequence,lcs

难度：1

对于序列S和T，求它们的最长公共子序列。例如X={A,B,C,B,D,A,B}，Y={B,D,C,A,B,A}则它们的lcs是{B,C,B,A}和{B,D,A,B}。求出一个即可。

对于两个序列X[0] .. X[i], Y[0] ... Y[j],  i,j分别表示两个序列的最后一个元素的索引。我们可以用一个二维数组记录序列Xi，Yi的LCS的长度。当i == 0 || j == 0时，公共子序列长度为0，这是边界。有以下转移方程：

c[i][j] => (

if(i ==0 || j == 0)  return 0; 

if(X[i] == Y[j]) return c[i - 1][j - 1] + 1; 

if(X[i] != Y[j]) max( c[i -1][j] , c[i][j-1] );

)

然后用一个数组b[i][j] 是由哪个子问题得到的最优解，用来构造最长公共子序列。

以下用string来表达序列。
#include <string>
#include <iostream>
using namespace std;

#define SIZE 20
// 注意i和j是从1开始的
void LCS(const string& str, int b[][SIZE] , int i, int j)
{
// 递归终点
if (i == 0 || j == 0) return;
if (b[i][j] == 1)
{
LCS(str, b, i - 1, j - 1);
cout << str[i - 1];
}
else if (b[i][j] == 2)
{
LCS(str, b, i -1, j);
}
else if (b[i][j] == 3)
{
LCS(str, b, i, j - 1);
}
}

int main()
{
for (;;)
{
int c[SIZE][SIZE] = {}; // 最优解
int b[SIZE][SIZE] = {}; // 构造最长公共子序列
string x;
string y;
cin >> x >> y;

// 边界，两者中一方是空的，则公共子序列是空的
// 注意c和d并不是从0开始的，0留给了边界
for (int i = 0; i < SIZE; i++) c[0][i] = 0; 
for (int i = 0; i < SIZE; i++) c[i][0] = 0;

for (int i = 1; i <= x.size(); i++) // 从最短的开始
{
for (int j = 1; j <= y.size(); j++)
{
#define max(x,y) (x >y ?x:y)
if (x[i - 1] != y[j - 1]) // 注意这里
{
if (c[i - 1][j] > c[i][j - 1])
{
c[i][j] = c[i - 1][j];
b[i][j] = 2;// 标记为2
}
else if (c[i - 1][j] < c[i][j - 1])
{
c[i][j] = c[i][j - 1]; 
b[i][j] = 3;// 标记为3
}
}
else if (x[i - 1] == y[j - 1])
{
c[i][j] = c[i - 1][j - 1] + 1;
b[i][j] = 1; // 标记为1
}
}
}
LCS(x, b, x.size(), y.size());
}
}

7.5 最长递增子序列

可以转化为最长公共子序列问题，将序列排序形成新序列， 然后找两个序列的最长公共子序列就是结果。或者用动态规划：

对于一个序列如1，-1，2，-3，4，-5，6，-7，其最长第增子序列为1,2,4,6。

设长度为N的数组为{a0，a1, a2, ...an-1)，则假定以aj结尾的数组序列的最长递增子序列长度为L(j)，则L(j)={ max(L(i))+1, i<j且a[i]<a[j] }。也就是说，我们需要遍历在j之前的所有位置i(从0到j-1)，找出满足条件a[i]<a[j]的L(i)，求出max(L(i))+1即为L(j)的值。最后，我们遍历所有的L(j)（从0到N-1），找出最大值即为最大递增子序列。时间复杂度为O(N^2)。
例如给定的数组为{5，6，7，1，2，8}，则L(0)=1, L(1)=2, L(2)=3, L(3)=1, L(4)=2, L(5)=4。所以该数组最长递增子序列长度为4，序列为{5，6，7，8}。算法代码如下：

#include <iostream>  
using namespace std;  
#define len(a) (sizeof(a) / sizeof(a[0])) //数组长度  
int lis(int arr[], int len)  
{  
    int longest[len];  
    for (int i=0; i<len; i++)  
        longest[i] = 1;  
  
    for (int j=1; j<len; j++) {  
        for (int i=0; i<j; i++) {  
            if (arr[j]>arr[i] && longest[j]<longest[i]+1){ //注意longest[j]<longest[i]+1这个条件，不能省略。  
                longest[j] = longest[i] + 1; //计算以arr[j]结尾的序列的最长递增子序列长度  
            }  
        }  
    }  
  
    int max = 0;  
    for (int j=0; j<len; j++) {  
        cout << "longest[" << j << "]=" << longest[j] << endl;  
        if (longest[j] > max) max = longest[j];  //从longest[j]中找出最大值  
    }  
    return max;  
}  
  
int main()  
{  
    int arr[] = {1, 4, 5, 6, 2, 3, 8}; //测试数组  
    int ret = lis(arr, len(arr));  
    cout << "max increment substring len=" << ret << endl;  
    return 0;  
}  

7.6 0-1背包

难度评级：★★

　　一个贼在偷窃一家商店时发现了n件物品，其中第i件值vi元，重wi磅。他希望偷走的东西总和越值钱越好，但是他的背包只能放下W磅。请求解如何放能偷走最大价值的物品，这里vi、wi、W都是整数。

解法：

　　如果每个物品都允许切割并只带走其一部分，则演变为部分背包问题，可以用贪心法求解。0-1背包问题经常作为贪心法不可解决的实例（可通过举反例来理解），但可以通过动态规划求解。

　　为了找出子结构的形式，粗略地分析发现，对前k件物品形成最优解时，需要决策第k+1件是否要装入背包。但是此时剩余容量未知，不能做出决策。因此把剩余容量也考虑进来，形成的状态由已决策的物品数目和剩余容量两者构成。这样，所有状态可以放入一个n*(W+1)的矩阵c中，其值为当前包中物品总价值，这时有

#include <stdio.h>
#include <stdlib.h>

int package_dp(int *v,int *w,int n,int total) {
    int i,j,tmp1,tmp2;
    int **c = (int **)malloc((n+1)*sizeof(int *));
    for(i=0;i<n+1;i++)
        c[i]=(int *)malloc((total+1)*sizeof(int));
    for(i=0,j=0;j<total;j++)
        c[i][j] = 0;
    for(i=1;i<=n;i++) {
        c[i][0] = 0;
        for(j=1;j<=total;j++) {
            if(w[i]>j)
                c[i][j] = c[i-1][j];
            else {
                tmp1 = v[i]+c[i-1][j-w[i]];
                tmp2 = c[i-1][j];
                c[i][j]=(tmp1>tmp2?tmp1:tmp2);
            }
        }
    }
    printf("c[%d][%d]:%d\n",n,total,c[n][total]);
    return 0;
}

int main() {
    int v[] = {0,10,25,40,20,10};
    int w[] = {0,40,50,70,40,20};
    int total = 120;
    package_dp(v,w,sizeof(v)/sizeof(int)-1,total);
    return 0;
}