最长递增子序列 (Longest Increasing Subsequence)

我转载两篇文章，同时给出一道poj题目的链接 http://poj.org/problem?id=2533

转自：http://blog.csdn.net/linulysses/article/details/5559262

问题描述： 给定一个序列 A_n = a₁ ,a₂ ,  ... , a_n ，找出最长的子序列使得对所有 i < j ，a_i < a_j 。

显然，暴力算法的时间复杂度是 O(2ⁿ ) ，因为搜索空间呈指数级增长。对于这种问题，如果要找复杂度为多项式时间的算法，自然而然地会想到动态规划。首先，要找出一种方法把该问题分解成只有多项式个子问题。考虑 a _n。如果最长递增子序列不 包含 a_n ，则问题变成要在 A_n-1 中找最长递增子序列；否则，在 A_n-1 之中找最长子序列，且该子序列的最大值不能超 a_n ，然后再加上 a_n 。按照这个思路，我们可以用 OPT(i,x _i ) 来表示A_i 之中最大值不超 x 的最长递增子序列的长度，则

      OPT(i,x) = max { OPT(i-1,x), a_i <= x ? OPT(i-1, a_i )+1 : OPT(i-1, x) }

而最终的解则是 OPT(n, max{a₁ , a₂ , ..., a_n }) 。现在的问题是，到底有多少个子问题？这取决于 x 有多少个取值。观察上述递归式可知， x 只可能取 A_n 中的值，因此最多有 n 可能的值。所以，子问题的数量为 O(n² ) 个。而每个子问题都可以通过 O(1) 的时间获解，从而总的时间复杂度是 O(n² ) 。实现这个算法的时候，为了能够使用一个二维的数组 S[n][n] 来存储状态，我们可以用数组 B 来保存排过序之后的 序列 b₁ , b₂ , ..., b_n ，从而对所有i < j ，b_i < b_j 。而 S[i][j] 表示 A_i 中最大值不超过 b_j 的最长递增子序列的长度。

[c-sharp]  view plain copy

1. copy A to B;  
2. sort B by increasing order;  
3. // initialize S[1][*]  
4. for j from 1 to n  
5.     if A[1] > B[j] then S[1][j] := 0  
6.     else S[1][j] := 1  
7. end  
8. for i from 2 to n  
9.     for j from 1 to n  
10.         if A[i] > B[j] then S[i][j] := S[i-1][j]  
11.         else   
12.             find the index k of A[i] in B[j]; // k <= j  
13.             S[i][j] := max( S[i-1][j], S[i-1][k]+1 )  
14.         end  
15.     end  
16. end  
17. return S[n][n]  

上述的实现中，第12行代码，我们可以使用一个哈希表来实现。我们可以看到，这个实现并不优美，既要一个额外的数组 B 和一次排序，还要动用一个哈希表。可以说，我们对这个算法并不满意。

那么，有没有更好一些的算法呢？注意到排序之后，我们得到新的序列 B，那么，显然，A 中最长的递增子序列也是序列 B 的子序列。因此，我们可以应用求两个序列的最长公共子序列的经典算法来求解，时间复杂也是 O(n²)。这样，我们就省去了使用哈希表的“不雅之举”。然而，最长公共子序列是一个更一般的算法，它不要求在序列的元素之间有序关系 。那么，我们能不能利用本问题中元素之间的序关系来设计一个优美的算法呢？答案是肯定的。注意中在上面的算法中（姑且称为算法1吧），一共有 O(n² ) 个子问题。事实上，我们从另外的角度看这个问题，从而获得只有 O(n) 中子问题的算法，只不过，计算每个子问题需要 O(n) 的时间。

首先，假设 A_n 中最长递增子序列 L 的最后一个元素是 a_t 。考虑 L 中在 a_t 之前的元素 a_k ，则 a_k <= a_t ，并且 L 由在 A_k 中包含 a_k 的最长递增子序列和 a _t 组成。因此，用 OPT(i) 来表示 A_i 中包含 a_i 的最长递增子序列。因为L 的最后一个元素必定在 A_n 中，因此，L 的长度为 OPT(1), OPT(2), ..., OPT(n) 中的最大值。这也证明了此算法的正确性。接下的问题就是如何计算OPT(i)。事实上，

      OPT(i) = max { OPT(j) |  j < i 且 a_j < a_i }

也就是说，通过遍历一次 已经计算好的 OPT(1), OPT(2), ..., OPT(i-1) 就可以计算出 OPT(i)，其时间复杂度为 O(i)。总的时间复杂度为 O(1) + O(2) + ... + O(n-1) = O(n² )。

[cpp]  view plain copy

1. // Let S[i] be OPT(i)  
2. S[1] := 1  
3. L := S[1]  
4. for i from 2 to n  
5.     S[i] := 1 // at least contain A[i]  
6.     for j from 1 to i-1  
7.         if A[j] < A[i] then S[i] := max( S[i], S[j]+1 )  
8.     end  
9.     L := max( L, S[i] )  
10. end  
11. return L  

相比之一，这个算法的实现就很干净和优美，且不容易出错。

O(n² ) 的时间复杂度似乎已经是很不错了。那么，有没有更快的算法？事实上，存在 O(n logk ) 的算法，其中 k为最长递增子序列的长度。为了达到这个时间复杂度，就需要费点脑筋了。考虑 A_i = a₁ , a₂ , ..., a_i 。记 Tail(X)为递增序列 X 的最后一个元素（尾元素），令 R_i,j 表示 A_i 中所有长度为 j 的递增子序列的集合。在所有属于 R _i,j的序列的尾元素中，必有一个最小值 ，记为用 m_i,j 。则

观察一 ：   对任何 i，m_i,1 <= m_i,2 <= ... <= m_i,j 。 

因此，如果我们想找以 a_i+1 结尾的最长递增子序列，则只要找到 k ，使得 m_i,k < a_i+1 <= m_i,k+1 ，并且该最长递增子序列的长度为 k+1 。对于这个搜索过程，利用上述观察一，可以使用二分法搜索 (binary search)。同时，我们注意到

观察二 ：   m_i+1,k+1 = a_i+1 , 且对于所有 t 不等于 k+1 ， m_i+1,t = m_i,t 。

同时，注意到计算 m_i+1,* 只需要用到 m_i,* ，因此，我们可以用 K[j] 来表示 m_i,j 在 A_i 中的下标。当计算在 A_i+1 上进行时，我们只需在 K 中找到一个元素 k ，使得 m_i,k < a_i+1 <= m_i,k+1 ，然后更新 K[k+1] ，这时，K[j] 就可以表示 m_i+1,j 在 A_i+1 中的下标了。

[cpp]  view plain copy

1. L := 0  
2. for i from 1 to n  
3.     find j in K such that A[K[j]] < A[i] <= A[K[j+1]] by binary search; if no such j, then set j := 0  
4.     P[i] := K[j]  
5.     if j == L or A[i] < A[K[j+1]] then  
6.         K[j+1] := i  
7.         L := max( L, j+1 )  
8.     end  
9. end  
10. return L  

其中，P[k] 记录了以 a_k 结尾的最长递增子序列的前一个元素在 A_n 中的下标。通过数组 P，我们就可以求出最长递增子序列，为

        ..., A[P[P[L]]], A[P[L]], A[L]。从上述伪代码可以看出，总的时间复杂度为 O(n logL )。

 

注：这是一道经典的题目。本人解这道题时，想到了第一种算法，在脑子里写出了大概的递归方程。然而，当真正在纸上写下来并建立一个二维的状态表时，才意识到问题没那简单。一开始，我的递归方程是 OPT(i,x) = max { OPT(i-1,x),  OPT(i-1, min(x, a_i ))+1  } 。这是不对的，但似乎不太容易发现，直到建立状态时，我才发现。而且，这个算法对于有重复元素的序列也存在问题，需要进一步修正。而第二个算法则是比较普遍的动态规划解，该算法优美，且对任何序列都不会有问题。最后一个算法则是参考了Algorithmist 上的一篇文章 。同时，该文章也给出 C/C++ 实现。以前解类似的动态规划问题时，总只是在脑子里把递归方程过一遍就觉得OK了，看来，以后要尽量把方程精确写出来并画画状态表。同时，也可以进一步思考更优美的动太规划算法，例如第二个算法。前两个算法都是典型的动态规划思路，但显然，第二个算法的状态要少，也更好。自勉之。

转自：http://www.wutianqi.com/?p=1850
引出：
问题描述:给出一个序列a1,a2,a3,a4,a5,a6,a7….an,求它的一个子序列（设为s1,s2,…sn），使得这个子序列满足这样的性质，s1<s2<s3<…<sn并且这个子序列的长度最长。输出这个最长的长度。（为了简化该类问题，我们将诸如最长下降子序列及最长不上升子序列等问题都看成同一个问题，其实仔细思考就会发现，这其实只是<符号定义上的问题，并不影响问题的实质）
例如有一个序列:1  7  3  5  9  4  8,它的最长上升子序列就是 1 3 4 8 长度为4.
分析：
这题目是经典的DP题目，也可叫作LIS（Longest Increasing Subsequence）最长上升子序列 或者 最长不下降子序列。很基础的题目，有两种算法，复杂度分别为O(n*logn)和O(n^2) 。
算法1:
时间复杂度:O(n^2):
我们依次遍历整个序列，每一次求出从第一个数到当前这个数的最长上升子序列，直至遍历到最后一个数字为止，然后再取dp数组里最大的那个即为整个序列的最长上升子序列。我们用dp[i]来存放序列1-i的最长上升子序列的长度，那么dp[i]=max(dp[j])+1,(j∈[1, i-1]); 显然dp[1]=1,我们从i=2开始遍历后面的元素即可。


// Author: Tanky Woo
// Blog: www.WuTianQi.com
int dp[1000];
int LIS(int arr[1000], int n)
{
for(int i=1; i<=n; ++i)
dp[i] = 0;
int ans;
dp[1] = 1;
for(int i=2; i<=n; ++i)
{
ans = dp[i];
for(int j=1; j<i; ++j)
{
if(arr[i]>arr[j] && dp[j]>ans)
ans = dp[j];
}
dp[i] = ans+1;
}
ans = 0;
for(int i=1; i<=n; ++i)
{
if(dp[i] > ans)
ans = dp[i];
}
return ans;
}


算法2:
时间复杂度:(NlogN):
除了算法一的定义之外，增加一个数组b，b[i]用以表示长度为i最长子序列的最后一个数最小可以是多少。易证：i<j时，b[i]<b[j]。
在二分查找时，一直更新b[]内容，设此时b的总长度为k，
若1. arr[i] >= b[k], 则b[k+1] = arr[i];
若2. arr[i] <  b[k], 则在b[1..k]中用二分搜索大于arr[i]的最小值，返回其位置pos，然后更新b[pos]=arr[i]。


// Author: Tanky Woo
// Blog: www.WuTianQi.com
// num为要查找的数,k是范围上限
// 二分查找大于num的最小值，并返回其位置
int bSearch(int num, int k)  
{  
    int low=1, high=k;  
    while(low<=high)  
    {  
        int mid=(low+high)/2;  
        if(num>=b[mid])  
            low=mid+1;  
        else   
            high=mid-1;  
    }  
    return low;  
}  
 
int LIS()
{
int low = 1, high = n;
int k = 1;
b[1] = p[1];
for(int i=2; i<=n; ++i)
{
if(p[i]>=b[k])
b[++k] = p[i];
else
{
int pos = bSearch(p[i], k);
b[pos] = p[i];
}
}
return k;
}




以下是证明b[]的单调递增性:
b序列是严格递增的，即b[1] < b[2] < … < b[t]。
证明：
若b[i] >= b[i + 1],b[i + 1] 是长度为i+1的递增子序列的尾项的最小值，设此序列为x[1]..x[i+1]，x[1]..x[i]即构成长度为i的递增子序列，x[i] < x[i+1] = b[i+1] <= b[i]，与b[i]定义不符。
最后，给出两个有代表性的题目：
1.HDOJ 1257 最少拦截系统
题目传送门：http://acm.hdu.edu.cn/showproblem.php?pid=1257
解题报告传送门：http://www.wutianqi.com/?p=1841
此题用O(n^2)解法做即可。
2.HDOJ 1025  Constructing Roads In JGShining’s Kingdom
题目传送门：http://acm.hdu.edu.cn/showproblem.php?pid=1025
解题报告传送门：http://www.wutianqi.com/?p=1848
此题数据量较大，所以要用O(NlogN)的解法做。