生物信息挖掘中ＬＩＳ算法研究
2009-01-01 00:00:00嚴華云李剛張建宏
計算機應用研究 2009年1期
（1.湖州師范學院 信息工程學院， 浙江 湖州 313000； 2.同濟大學 電子與信息工程學院， 上海 201804）
摘 要：探討了生物信息挖掘中ó模式子序列問題的一個特例，即最長遞增子序列(LIS)問題。對于LIS問題，分別用LCS算法、動態規劃、動態規劃結合二分法進行求解，并分析了這三種算法的時間和空間復雜度，對其中兩種算法進行了實現，驗證了時間和空間復雜性理論分析的正確性，最后得出了一種高效的LIS算法。
關鍵詞：最長遞增子序列; 動態規劃; 生物信息挖掘
中圖分類號：TP391;TP301.6 文獻標志碼：A
文章編號：10013695(2009)01006202
LIS algorithm for bioinformatics mining
YAN Huayun1，2， LI Gang1， ZHANG Jianhong1
(1.School of Information Engineering，Huzhou Teachers College， Huzhou Zhejiang 313000， China； 2.College of Electronics Information Engineering，Tongji University， Shanghai 201804， China）
Abstract:This paper introduced a special problem of ó pattern subsequence in bioinformatics mining， that was LIS problem. This paper gave three algorithms of the LIS problem， they were LCSbased， dynamic programming， dynamic programming binary search， and analyzed the advantage and disadvantage of the three algorithms， and tested the performance of the three algorithms’ complexity through experiment. Finally， it educed an effective LIS algorithm.
Key words：LIS(longest increasing subsequence); dynamic programming; bioinformatics mining
最長遞增子序列（LIS）問題是一個組合數學中的問題，它不純粹是理論問題。近年來， LIS等組合學問題已經在生物信息學中廣泛應用，例如LCS（最長公共子序列）問題被應用到對兩個DNA序列進行比較，從而求基因之間的相似性；LIS被應用到DNA序列的BLAST (basic local alignment search tool)數據庫和MUMmer系統中。LIS問題是最長ó模式子序列問題的一個特例[1，2]。
一般情況下， 最長ó模式子序列問題的提法如下:設ó是由字母U和D組成的一個模式， 如ó=UUDUD；設ó∞是由ó組成的無窮序列óóó…，如(UUDUD)∞=UUDUDUUDUD…。 在此例中，UUDUDUU是ó∞的長度為7 的前綴。
設T=t0t1…tm-1是由U和D組成的長度為m的字符串。對于序列V=v0v1…vm，當ti=D時，vi＞vi+1；當ti=U時，vi＜vi+1。此時，稱字符串T=t0t1…tm-1是序列V=v0v1…vm的升降模式。設Wn是1，2，…，n的全體排列組成的集合，對于給定的模式ó和任意a∈Wn，a的以ó∞的前綴為升降模式的子序列稱為a的一個ó模式子序列。當T=Um時，該問題轉換為求ó∞的LIS問題，這正是本文要探討的問題。
動態規劃是求解決策過程的有效最優數學方法之一[3]，它為具有最優子結構性質的實際問題提供了一種重要的解決問題的途徑。該策略最早由Bellman提出，它利用最優性原理和所獲得的遞推關系去解最優決策序列，使問題的計算量急劇下降[4，5]。
LIS問題在國內沒有資料給出其詳細的實現過程。鑒于LIS問題是一個經典的組合數學問題，且在計算機科學中越來越有其實際應用領域。本文準備用多種方法來詳細實現它，并結合實驗進行分析驗證。
最長遞增子序列問題是這樣的：設S=〈a1，a2，…，an〉是n個數的序列，S的一個子序列Ssub定義為將S中的零個或多個元素去掉后剩下的序列。而S的遞增子序列是這樣一個子序列： Sis=〈ak1，ak2，…，akm〉，其中ak1＜ak2＜…＜akm。 LIS問題就是求所有Sis序列中元素最多的那個子序列，從而求得最長遞增子序列Sis。
1 最長單調遞增子序列問題求解
1.1 轉換為LCS問題求解
設序列S′=〈b1，b2，…，bn〉是對序列S=〈a1，a2，…，an〉按遞增排好序的序列，顯然S′與S的最長公共子序列即為S的最長遞增子序列。這樣就把求LIS的問題轉換為求LCS了。轉換為LCS求解的算法描述如下：
算法1
輸入：n個數的數組序列S[1…n]
a)對S[1…n]進行排序，生成一個升序序列S′[1…n]；
b)調用LCS(S，S′)算法，得到LCS序列，從而得到LIS序列;
c)輸出LIS序列。
對于時間來說，算法中的a)步，即對于一般的數組進行排序，可以在O(n log n)時間內解決；算法中的b)步，即LCS問題的時間復雜度為O(n2)。因此總的時間復雜度為O(n log n)+O(n2)，即O(n2)，它表示的是該算法在最壞情況下的運行時間，并不表示一定就是O(n2)時間。后面實驗分析時會討論到這個問題。
對于空間來說，算法中a)步，即對于數組的排序問題來說，可以在原來的數組空間內額外再加一兩個單元就可以解決，即為O(1)；算法中b)步，即對于LCS的動態規劃問題求解中，它需要額外的數組長度n的兩倍空間，即為O(n)。因此總的需要的空間復雜度為O(1)+ O(n)，即O(n)。
1.2 一般動態規劃法求解
對于一個問題可以用動態規劃問題求解的話，一定要滿足最優化原理。下面就來討論LIS問題是否滿足最優化原理。用L(i)表示S中以ai（1≤i＜n）為末元素的最長遞增子序列的長度，則有如下的遞推公式：
L(1)=1L(i)=max{L(k)}+1 if 1≤k＜i， ak≤ai
L(i)=1else(1)
式（1）是在求以ai為末元素的最長遞增子序列時，在S中找到所有序號i前面且小于ai的元素ak，1≤k＜i且ak≤ai。如果這樣的元素存在，那么對所有ak，都有一個以ak為末元素的最長遞增子序列的長度L(k)，把其中最大的L(k)選出來（即在序號i前面的，滿足1≤k＜i并且ak≤ai調節的其中最大的L(k)），那么L(i)=L(k)+1，即得到以ai為末元素的最長遞增子序列，相當于以使L(k)最大的那個ak為末元素的遞增子序列最末再加上ai；如果這樣的元素不存在，那么ai自身構成一個長度為1的以ai為末元素的遞增子序列。
從遞推式（1）可以看出，LIS問題滿足最優化原理。LIS的動態規劃算法實現如下：
算法2
輸入：n個數的數組序列S[1…n]
a)生成數組L[1…n]，初值都設為1，max←1；
b)i←2；
c)k←1，如果max d)如果k=i=n，轉h)； e)如果k=i，則i←i+1，并轉c)； f)如果k g) k←k+1，轉d)； h)j←n；//以下輸出LIS序列 i)如果max=0，轉l)； j)如果L[j]=max并且S[j]小于棧頂元素，則S[j]入棧，max←max-1； k) j←j-1，轉i)； l)取出stack中的元素組成一個LIS序列。 從時間復雜性來說，從a)~g)來看，k每次從1執行到i，而i從2執行到n，可以看出該動態規劃算法是一個二重循環，因此其時間復雜度為O(n2)；而算法從h)~l)只需要掃描數組L一遍就可以得到LIS序列，其時間復雜度為O(n)。總的時間復雜度為O(n)+ O(n2)，即O(n2)。 對于空間來說，由于只是新生成了一個數組，因此總的空間復雜度為O(n)。 算法1和2進行比較，兩者的最壞時間復雜度的階是相同的。但算法1的時間是由排序O(n log n)和LCS（O(n2)）所需的時間組成的；算法2只需要O(n2)的時間，所以算法2比1略優。 經過算法2前七步生成了數組L[1…n]： 以及max=5；根據算法2的h)~l)可以得到進棧順序為8、7、4、2、1，因此出棧后的序列即為LIS的一個序列1、2、4、7、8。 1.3 對算法2的改進 在算法2中，計算每一個L(i)時都要找出最大的L(j)(j＜i)。由于L(j)沒有順序，只能順序查找滿足L(j)＜L(i)最大的L(j)。如果能讓L(j)有序，就可以使用二分查找，這樣算法的時間復雜度就可能降到O(n log n)。受此啟發，可以這樣改進算法2得到。 算法3 輸入：n個數的數組序列S[1…n] a)新生成數組L[1…n]和link[1…n]，max=1；/* link[i] 用于保存長度為i的末元素序列（用鏈表實現），以便高效返回LIS序列*/ b)L[1]←S[1]并將S[1]的值鏈接到link[1]的尾部，i←2； c)如果i＞n，轉l)； d)l←1，h←i；//l和h分別為查詢到該有序數組L的下界和上界 e)如果l＞h，轉i)； f)m←(l+h)/2; g)如果L[M] h)否則，h←m-1，并轉e)； i)L[h]←S[i]，S[i]的值鏈接到link[h]的尾部； j)如果h＞max，則max←max+1； k)i←i+1，轉c)； l)取link[max] 鏈中第一個元素放入stack中 //以下輸出LIS序列 m)如果max=0，轉p)； n)max←max-1； o)在link[max]鏈中從頭開始找第一個滿足小于棧頂元素的元素，并將它壓入棧中，轉m)； p)取出stack中的元素就組成了一個LIS序列。 下面簡單分析該算法的合理性。由算法3可以看出，數組L中存儲LIS序列長度等于下標值末元素的當前最小值，并且LIS長度為i的末元素的值不會比長度為i-1的末元素小，因此數組L是一個遞增序列，從而可以用二分查找進行優化。 算法3從時間上來看，循環次數為n，每次循環二分查找用時步大于log n，所以算法的時間復雜度為O(n log n)。這個算法在算法2的基礎上得到了較好的改進。 算法3從空間上來看，用到了數組L[1…n]和數組鏈表link[1…n]以及堆棧stack，其空間分別都不大于n，因此算法的空間復雜度為O(n)。 例2 以例1中的序列S[1…n]來說，算法3前11步生成了數組鏈表link[1…n]，如圖1所示，以及max=5；根據算法3的l)~p)可以得到進棧順序為8、7、5、2、1，因此LIS的一個序列為1、2、5、7、8。 2 實驗 由于經分析算法1不如算法2，本文只做了算法2和3的實驗。實驗用Java實現，實驗機器的CPU速度為1.6 GHz。因在數組中一個個查找元素的速度很快，本文將每查找一個元素就顯示到控制臺，從而放大時間的對比效果。實驗數據是這樣的，由計算機生成五組數組長度為1萬、10萬、50萬、100萬個的數據（由于內存限制沒辦法得出更大級別數組的實驗）；同時為了效果更好，將所有的生成數據前5千個從1到5千按順序生成，其余數據隨機生成。執行時間取五組的平均值。 從圖2可以看出，算法3所用時間比算法2要少一點，但也不是想象中的那么大的差距。其原因是數組中的絕大部分數據都是隨機生成的(并非最壞情況)。 對這兩個算法的理論分析上的時間復雜度分別是O(n2)和O（n log n），這是最壞情況下的時間，而非壞的情況下則沒有那么大的差距。另一方面，圖示基本上是一條直線，即運行時間是與數組長度成比例的，這可能與生成的數據是隨機數有關系。同時可以看出，隨著數組長度越來越大，算法3更加優越。 3 結束語 動態規劃算法通常用于求解具有某種最優性質的問題。在這類問題中，可能會有許多可行解，從本文所舉兩個例子就可以看出，其解并不相同。 動態規劃算法與分治法類似，其基本思想也是將待求解問題分解成若干個子問題，先求解子問題，然后從這些子問題的解得到原問題的解。與分治法不同的是，適合于用動態規劃求解的問題，經分解得到子問題往往不是互相獨立的。若用分治法來解這類問題，因分解得到的子問題數目太多，有些子問題被重復計算了很多次。如果能夠保存已解決的子問題的答案，而在需要時再找出已求得的答案，這樣就可以避免大量的重復計算，節省時間。可以用一個表來記錄所有已解的子問題的答案，不管該子問題以后是否被用到，只要它被計算過，就將其結果填入表中。這就是動態規劃法的基本思路。 參考文獻： [1]STANLEY R P. Increasing and decreasing subsequences of permutations andtheirvariants[C]//Proc of International Congress of Mathematicians. 2005. [2]FIRRO G， MANSOUR T. Longest alternating subsequences in patternrestricted permutations[J]. Electronic Journal of Combinatorics， 2007，14(1):117. [3]COOPER R. Dynamic programming:an overview[EB/OL]. (20010214). http://econ.bu.edu/faculty /cooper/ dynprog/ introlect.pdf. [4]ALSUWAIYEL M H. Algorithms design techniques and analysis[M]. Beijing:Publishing House of Electronics Industry， 2003. [5]王曉東.計算機算法設計與分析[M].北京：電子工業出版社，2001:4344.
