LCIS 最長公共上升子序列問題DP算法及優化

---

一. 知識簡介

1. 學習 LCIS 的預備知識： 動態規划基本思想, LCS, LIS

2. 經典問題：給出有 n 個元素的數組 a[] , m 個元素的數組 b[] ,求出它們的最長上升公共子序列的長度.

3. 例如:

a[] data:
5
1 4 2 5 -12

b[] data:
4
-12 1 2 4

LCIS is 2
LCIS 所含元素為 1 4

---

二.LCIS問題分析

1. 確定狀態
     可以定義 dp[i][j] 表示以 a 數組的前 i 個整數與 b 數組的前 j 個整數且以 b[j] 為結尾構成的公共子序列的長度。
     對於解決DP問題,第一步定義狀態是很重要的!
     需要注意，以 b[j] 結尾構成的公共子序列的長度不一定是最長公共子序列的長度！

2. 確定狀態轉移方程

   • 當 a[i] == b[j] 時,我們只需要在前面找到一個能將 b[j] 接到后面的最長的公共子序列.
     • 之前最長的公共子序列在哪呢？首先我們要去找的 dp[][] 的第一維必然是 i - 1 ,因為 i 已經拿去和 b[j] 配對去了，不能用了。並且也不能是 i - 2 ，因為 i - 1 必然比 i - 2 更優。
     • 第二維呢？那就需要枚舉 b[1]...b[j-1] 了，因為你不知道這里面哪個最長且哪個小於 b[j] 。
     • 這里還有一個問題，可不可能不配對呢？也就是在 a[i] == b[j] 的情況下，需不需要考慮 dp[i][j] == dp[i-1][j] 的決策呢？答案是不需要。因為如果 b[j] 不和 a[i] 配對，那就是和之前的 a[1]...a[j-1] 配對（假設 dp[i-1][j]>0 ，等於0不考慮），這樣必然沒有和 a[i] 配對優越。（為什么必然呢？因為 b[j] 和 a[i] 配對之后的轉移是 max(dp[i-1][k])+1 ，而和之前的 i1 配對則是 max(dp[i1-1][k])+1 。
   • 當 a[i] != b[j] 時, 因為 dp[i][j] 是以 b[j] 為結尾的LCIS，如果 dp[i][j] > 0 那么就說明 a[1]..a[i] 中必然有一個整數 a[k] 等於 b[j] ,因為 a[k] != a[i] ，那么 a[i] 對 dp[i][j] 沒有貢獻，於是我們不考慮它照樣能得出 dp[i][j] 的最優值。所以在 a[i] != b[j] 的情況下必然有 dp[i][j] == dp[i-1][j] 。

• 綜上所述,我們可以得到狀態轉移方程:
  ① a[i] != b[j], dp[i][j] = dp[i-1][j]
② a[i] == b[j], dp[i][j] = max(dp[i-1][k]+1) (1 <= k <= j-1 && b[j] > b[k])

---

三. LCIS算法實現及遞進的優化方法

• O(N * M^2) 算法
  分析： 根據狀態轉移方程可以得到非常容易的寫出 O(N * M^2) 的算法

// 復雜度 O(N * M^2)
int lengthOfLCIS(int *a, int n, int *b, int m) {
    int ans = 0;
    int dp[505][505] = {0};
    int max_dp = 0;
    for (int i = 0 ; i < n ; ++i) {
        for (int j = 0 ; j < m ; ++j) {
            dp[i + 1][j + 1] = dp[i][j + 1];
            if (a[i] == b[j]) {
                max_dp = 0;
                for (int k = 0 ; k < j ; ++k) {
                    if (b[j] > b[k] && max_dp <= dp[i][k + 1]) {
                        max_dp = dp[i][k + 1];
                    }
                }
                dp[i + 1][j + 1] = max_dp + 1;
            }
            ans = ans > dp[i + 1][j + 1] ? ans : dp[i + 1][j + 1];
        }
    }
    return ans;
}

• **O(N * M * log(M))**算法
  

  分析： 觀察O(N * M^2)的算法，我們可以發現當 a[i] == b[j] 時在前面找到一個能將 b[j] 接到后面的最大長度，這一個查找過程是可以優化的。
  設置新的輔助數組 len[]，len[i] = value 代表長度為 i 的所有公共子序列中最后一個數字為的最小值為value。
  • 為什么要這樣設計輔助數組 len[] ？
      首先我們優化的過程是一個查找過程( dp[i][j] = max(dp[i-1][k]+1) (1 <= k <= j-1 && b[j] > b[k]) )，實際上這個查找過程是在一個一維數組中查找滿足條件 b[j] > b[k] 的一個最大值。
      然后，既然是優化查找過程，可以采用 hash 或者是二分來將查找問題從O(N)降到一個可以接受的復雜度，大體判斷一下二分似乎更加適合該問題。
      由於二分查找的對象必須具有單調性，而在最長公共上升子序列中，所以的子序列是具有單調性的，所以我們可以將子序列的結尾數字作為數組中的值。
      為什么用數組 len[i] 來維護長度為 i 的所有公共子序列中最后一個數字為的最小值？因為維護一個最小值能夠判斷 b[j] 到底能夠接在哪個位置，試想一下，如果維護一個最大值，本來 b[j] 能接在長度為 i 的某一個子序列后面，因為 b[j] < 最大值導致查找失敗，這是不是很 **: ) **...
  • len[] 更新策略 ( 請參考下面代碼 )
      由於本人能力有限，len[] 數組的更新策略不能夠足夠嚴謹、清晰的表述出來，其中有一些細節着實不知道該怎么說，如果有好的見解請一定提出！
      首先用二分查找找到 b[j] 能“嫁接”到之前子序列中的最長長度( ind是二分出來的len[]數組下標，根據上面的定義，ind就是以b[j]結尾的最長公共上升子序列長度 )。
      當a[i] == b[j]，這時候 b[j] 將“嫁接”到之前子序列中的合適位置，同時更新 len[ind] = b[j] 。
      當a[i] != b[j]，我們也需要更新 len[] 的信息，因為后面的更新需要前面的數據。當 b[j] 與 a[i] 適配，那么我們需要查看 b[j] 與前 i - 1個數字所構成的以 b[j] 為結尾數字的最長的公共子序列。如果有多個跟 dp[i][j + 1] ( 因為所有 i , j 都 + 1 所以是dp[i][j + 1])相同長度的子序列，查看 len[dp[i][j + 1]] 的值是否比 d[j] 小，如果比 d[j] 小則更新 len[dp[i][j + 1]]的值。

// 復雜度 O(N * M * log(M))
int binary_search(int *len, int length, int value) {
    int l = 0, r = length - 1, mid;
    while (l < r) {
        mid = (l + r) >> 1;
        if (len[mid] >= value) {
            r = mid;
        } else {
            l = mid + 1;
        }
    }
    return l;
}
int lengthOfLCIS(int *a, int n, int *b, int m) {
    int ans = 0;
    int dp[505][505] = {0};
    int len[505];
    for (int i = 0 ; i < n ; ++i) {
        len[0] = INT_MIN;  
        for (int j = 0 ; j < m ; ++j) {
            len[j + 1] = INT_MAX;
            int ind = binary_search(len, j + 2, b[j]);
            if (a[i] != b[j]) {
                if (b[j] < len[dp[i][j + 1]]) { // 需要注意這里的更新策略,只有當b[j]的值小於當前len[dp[i][j + 1]]的最小值時才能更新
                    len[dp[i][j + 1]] = b[j];
                }
                dp[i + 1][j + 1] = dp[i][j + 1];
            } else {
                dp[i + 1][j + 1] = ind;
                len[ind] = b[j];
            }
            ans = ans > dp[i + 1][j + 1] ? ans : dp[i + 1][j + 1];
        }
    }
    return ans;
}

• O(N * M)算法
  分析：到這里好像沒有更好的方法來優化查找了，那讓我們再仔細的分析一下問題。
      當a[i] != b[j], dp[i][j] = dp[i-1][j],這對問題並沒有什么影響。
      當a[i] == b[j]，的時候 dp[i][j] = max(dp[i-1][k]+1) (1 <= k <= j-1 && b[j] > b[k])，我們發現一個特點，其實 a[i] ( b[j] )與 b[k] 的關系早就在之前就可以確定了！( i 是最外層循環， j` 是內層循環，當 j` 遍歷到 k 時，就足以判斷 b[j] ? b[j`] 的大小關系了)，因此我們只需要在內層循環與外層循環直接維護一個 max_dp 的值，等到 a[i] == b[j] 的時候，直接令 dp[i][j] = max_dp + 1即可，時間復雜度降到了 O(N * M)。
  注意：這里的優化由問題的特性決定的，即上一段黑色加粗部分文字。對於這種情況下，普通優化技巧已經很無力了，那我們要做的就是對問題不斷的思考，不斷的尋找探究問題的本質以及其獨有的特性，學而不思則惘！

// 復雜度 O(N * M)
int lengthOfLCIS(int *a, int n, int *b, int m) {
    int ans = 0;
    int dp[505][505] = {0};
    for (int i = 0 ; i < n ; ++i) {
        int max_dp = 0;
        for (int j = 0 ; j < m ; ++j) {
            dp[i + 1][j + 1] = dp[i][j + 1];
            if (a[i] > b[j] && max_dp < dp[i + 1][j + 1]) max_dp = dp[i + 1][j + 1];
            if (a[i] == b[j]) {
                dp[i + 1][j + 1] = max_dp + 1;
            }
            ans = max(ans, dp[i + 1][j + 1]);
        }
    }
    return ans;
}

• 將空間復雜度優化至O(M)
  分析：現在空間復雜度是O(N * M)，經過觀察可以發現，對於當 dp[i + 1][j + 1] 來說，我只需要上一層的數據和當前這一層左側的數據，因此 dp[][] 數組第一維開為 2 就足夠使用了，這種不斷利用兩層或者幾層數組滾動求解的技巧叫滾動數組。
  以下樣例以O(N * M)算法為基礎進行了空間上的優化

int lengthOfLCIS(int *a, int n, int *b, int m) {
    int ans = 0;
    int dp[2][505] = {0};
    for (int i = 0 ; i < n ; ++i) {
        int max_dp = 0;
        for (int j = 0 ; j < m ; ++j) {
            dp[(i + 1) % 2][j + 1] = dp[i % 2][j + 1];
            if (a[i] > b[j] && max_dp < dp[(i + 1) % 2][j + 1]) max_dp = dp[(i + 1) % 2][j + 1];
            if (a[i] == b[j]) {
                dp[(i + 1) % 2][j + 1] = max_dp + 1;
            }
            ans = max(ans, dp[(i + 1) % 2][j + 1]);
        }
    }
    return ans;
}

---

四.LCIS練習題

傳送門：HDU 1423 Greatest Common Increasing Subsequence (http://acm.hdu.edu.cn/showproblem.php?pid=1423)

AC 代碼 ：

/*************************************************************************
    > File Name: LCIS-1.cpp
    > Author: 
    > Mail: 
    > Created Time: 2017年09月03日 星期日 14時26分47秒
 ************************************************************************/

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

#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))

#define N2

#ifdef N2
int lengthOfLCIS(int *a, int n, int *b, int m) {
    int ans = 0;
    int dp[2][505] = {0};
    for (int i = 0 ; i < n ; ++i) {
        int max_dp = 0;
        for (int j = 0 ; j < m ; ++j) {
            dp[(i + 1) % 2][j + 1] = dp[i % 2][j + 1];
            if (a[i] > b[j] && max_dp < dp[(i + 1) % 2][j + 1]) max_dp = dp[(i + 1) % 2][j + 1];
            if (a[i] == b[j]) {
                dp[(i + 1) % 2][j + 1] = max_dp + 1;
            }
            ans = max(ans, dp[(i + 1) % 2][j + 1]);
        }
    }
    return ans;
}
#endif

int main() {
    int n, m, T;
    int a[505], b[505];
    scanf("%d", &T);
    while (T--) {
        scanf("%d", &n);
        for (int i = 0 ; i < n ; ++i) {
            scanf("%d", &a[i]);
        }
        scanf("%d", &m);
        for (int j = 0 ; j < m ; ++j) {
            scanf("%d", &b[j]);
        }
        printf("%d\n", lengthOfLCIS(a, n, b, m));
        if (T) printf("\n");
    }
    return 0;
}