基于多重序列所有公共子序列的啟發式算法度量多圖的相似度

王 旭,歐陽繼紅,陳桂芬

(1.吉林大學 計算機科學與技術學院,長春130012;2.吉林大學 符號計算與知識工程教育部重點實驗室,長春130012;3.吉林農業大學 信息技術學院,長春130118)

0 引 言

圖是最重要和值得深入研究的數據結構之一,例如在生物信息學和化學領域中,DNAs和分子都是以圖的結構出現的[1,2]。在分析這種圖數據類型時,度量圖的相似度承擔了一個重要角色[3,4],圖的相似度度量方法能夠處理這些領域的圖數據,對DNAs和分子進行提取和分類,有效地解決分類問題。隨著圖數據日益增加及其復雜性,需要更加有效的方法處理任意多個圖的相似度方法[5,6],而現有的大部分度量圖的相似度方法都是度量兩個圖的相似度[7,8]。圖相似度度量方法可將圖表示為序列[9-11],用度量序列的方法度量圖的相似度。度量序列相似度方法可以采用所有公共子序列方法(All common subsequences,ACS)[12,13]。ACS方法通過計算所有公共子序列數度量序列的相似度,與最長公共子序列方法相比,ACS不僅包含了最長的公共子序列,而且包含了第二長、第三長、……的公共子序列,最大化地獲取公共子序列的信息,更能反映出序列間的相似程度。所有公共子序列數越大,序列相似度越大;反之,序列相似度越小。度量序列的相似度也可以采用啟發式方法,在查找多重序列的匹配時,最大化啟發估計值,以便最少地擴展節點,進而找到最優解。但現有的啟發式算法提出的公共子序列只包含一個字符,不適應實際的應用[14],查找最長公共子序列的啟發式算法[15],不如ACS方法更能反映出序列之間的相似程度。為度量多圖的相似度,如何將圖表示為包含更多原圖信息的對應序列,并應用啟發函數減少擴展節點的個數,成為度量多圖相似度的重要問題。

本文提出了多重序列所有公共子序列的啟發式算法度量多圖的相似度,將多圖表示多重序列,當多重序列的一個匹配出現時,該算法遞歸地計算在匹配點上的所有公共子序列數,通過下標比較查找多重序列的公共子序列、剔除重復的匹配,通過后綴序列的啟發函數減小計算多重序列的所有公共子序列數的節點個數,有效地解決任意多個圖的相似度問題。

1 圖表示為序列

1.1 圖表示序列方法

將圖表示為序列,序列要保留原圖的信息,這樣通過度量序列的相似度度量圖的相似度時,度量結果才能保證圖的相似度的準確性。本文采用圖的深度優先搜索方法,以任一頂點作為序列的初始節點,按照該頂點與其它頂點間的路徑搜索,將圖表示為頂點序列。頂點序列體現了圖中頂點訪問的層次順序,反映了頂點間的路徑信息以及連通性,較好地保留了圖的信息。在搜索頂點時,若某一頂點已訪問過,就不再從該頂點出發進行搜索,搜索圖的過程實質上是對每個頂點查找其鄰接頂點的過程。為保證度量的準確性,選擇圖中相同或相似位置的頂點作為序列的初始節點。應用圖深度優先搜索方法可將d個圖G1,G2,…,G d表示為d個序列,這些d個序列組成集合S＝{s1,s2,…,s d},S為字母表Σ上的多重序列(multiple sequences),d≥2。本文研究的圖是頂點帶標號的有向圖。

以兩個圖G1和G2(如圖1所示)為例,選取具有相同結構的頂點(e和a)為序列的初始節點,應用圖深度優先搜索方法將G1和G2表示為對應頂點序列s1和s2,s1＝expect,s2＝accept。

圖1 有向圖G1和G2Fig.1 The directed graphs G1 and G2

1.2 所有公共子序列方法

所有公共子序列方法(ACS)可以度量序列的相似度,所有公共子序列數越大,序列的相似度越大;反之,序列的相似度越小。下面給出d個序列的所有公共子序列的定義,d≥2。

定義1 從序列s中刪除0個或更多的元素所獲得的序列稱為s的子序列:

(1)若序列x為s i的子序列,1≤i≤d,則稱x為多重序列S的公共子序列;

(2)滿足(1)的情況下所有的x稱為S的所有公共子序列。用|MACS|表示S的所有公共子序列數。

為度量多重序列S的相似度,應用所有公共子序列方法計算|M ACS|,|M ACS|越大,多重序列的相似度越大;反之,多重序列的相似度越小。以兩個序列s1,s2為例,長度分別為n1,n2,ACS(x,y)表示s1和s2所有公共子序列數,有下式成立。

式中:當x＝0或y＝0時,s1和s2的公共子序列僅有空集?,所以ACS(0,y)＝1,ACS(x,0)＝1和ACS(0,0)＝1;當s1[x]＝s2[y]時,s1[x]或s2[y]成為s1和s2新的公共子序列,s1[x]或s2[y]與ACS(x－1,y－1)已得到的公共子序列組成新的公共子序列,新的公共子序列數為ACS(x－1,y－1),在s1[x]或s2[y]點上的所有公共子序列數是ACS(x－1,y－1)的2倍,所以ACS(x,y)＝ACS(x－1,y－1)×2,0≤x≤n1,0≤y≤n2。公式(1)迭代地建立n1×n2矩陣M,計算s1和s2的所有公共子序列數,ACS(x,y)值越大,s1和s2越相似。

對于圖1中G1和G2的對應序列s1＝expect和s2＝accept,由公式(1)計算得到的所有公共子序列數如表1所示,s1和s2的所有公共子序列為:{?,e,p,c,t,ct,ep,et,pt,ept},s1和s2所有公共子序列數ACS(s1,s2)＝10。

表1 s1和s2的所有公共子序列數Table 1 The number of all common subsequences between s1 and s2

2 ACS的啟發式方法

多重序列表示為S＝{s1,s2,…,s d},d≥2,對應的長度分別為n1,n2,…,n d。S中長度分別為n i,n j任意兩個序列s i和s j的后綴序列為:s i[x＋1,…,n i]和s j[y＋1,…,n j],0≤x≤n i,0≤y≤n j,0≤i,j≤d,ACSsuf(x,y)ij表示兩個后綴序列的所有公共子序列數,有下式成立:

式中:當x＝n i或y＝n j,s i和s j沒有后綴序列,所以ACSsuf(n i,y)ij＝0,ACSsuf(x,n j)ij＝0和ACSsuf(n i,n j)ij＝0;當s i[x＋1]＝s j[y＋1]時,s i[x＋1]或s j[y＋1]與ACSsuf(x＋1,y＋1)ij得到的公共子序列組成新的公共子序列,所以ACSsuf(x,y)ij＝ACSsuf(x＋1,y＋1)ij×2。公式(2)迭代地建立(n i－x)×(n j－y)矩陣Msuf,計算s i和s j后綴序列的所有公共子序列數,ACSsuf(x,y)ij值越大,s i和s j后綴序列越相似。

p表示多重序列S上的點,p＝(p1,p2,…,p i,…,p d),1≤i≤d,1≤p i≤n i,p i表示序列n i上的點,n i表示n i的長度。對于點p的d個后綴序列s i[p i＋1,…,n i],可以應用公式(2)計算這些后綴序列的所有公共子序列數,h?(p)表示任意兩個后綴序列s i[p i＋1,…,n i]和s j[p j＋1,…,n j]的所有公共子序列數,1≤i,j≤d,有下式成立。

定理1 若h(p)表示d個后綴序列s1[p1＋1,…,n1],s2[p2＋1,…,n2],…,s i[p i＋1,…,n i],…,s d[p d＋1,…,n d]的所有公共子序列數,1≤i≤d,則h(p)≤h?(p)。

證明:d個后綴序列的所有公共子序列是d中任意兩個序列最少的所有公共子序列的子集,d個后綴序列的所有公共子序列數小于等于d中任意兩個序列的所有公共子序列數的最小值,由公式(3)可知,h?(p)表示任意兩個序列的所有公共子序列數的最小值,h?(p)肯定大于等于d個后綴序列的所有公共子序列數h(p),即h(p)≤h?(p),得證。

若s1[p1]＝s2[p2]＝…＝s i[p i]＝…＝s d[p d],稱點p＝(p1,p2,…,p i,…,p d)為多重序列S上的一個匹配。對于圖1中G1和G2的對應序列s1＝expect和s2＝accept,點(4,4)為s1和s2在字符e的匹配,(5,3)為s1和s2在c的匹配,s1和s2的所有匹配出現的位置如表2所示。

表2 s1和s2匹配出現的位置Table 2 The locations of the matches between s1 and s2

定義2 對于多重序列S上的兩個點p＝(p1,p2,…,p i,…,p d),q＝(q1,q2,…,q i,…,q d),如果滿足p i＜q i,1≤i≤d,則稱p的所有下標對應地小于q的所有下標,記Subscript(p)＜Subscript(q)。

若p和q為S上的兩個匹配,Combine(p,q)表示p和q的合并后的序列:

對于圖1中G1和G2的對應序列s1和s2,(1,4)和(3,5)分別為字符e和p的匹配,且Subscript(e)＜Subscript(p),合并(1,4)和(3,5)得到的序列為:Combine((1,4),(3,5))＝Combine(e,p)＝(ep)。

定理2 若點p和q為S上的兩個匹配,則Combine(p,q)為S的公共子序列。

證明:根據定義2和公式(4),當q的所有下標對應地大于p的所有下標,即p i＜q i,或q的所有下標對應地小于p的所有下標,即q i＜p i,才可以合并p和q。點p和q為S上的匹配,則p和q分別為S的公共子序列,對于s i上的點p i和q i,p i和q i分別為S的公共子序列,當p i＜q i或q i＜p i,則序列(p iq i)或(q ip i)也是S的公共子序列,1≤i≤d。所以合并p和q的序列Combine(p,q)為S的公共子序列,得證。

合并兩個匹配可以得到公共子序列,當新的匹配出現,如果新的匹配的下標大于或小于合并后序列的下標,就會產生新的公共子序列,所以需要標記合并后序列的下標。合并后的序列Com bine(p,q)的下標為p和q中最大的下標,有下式成立。

對于圖1中G1和G2的對應序列s1和s2,(1,4)為字符e的匹配,(3,5)為字符p的匹配,合并后的序列(ep)的下標Subscript(ep)＝Subscript(p)＝(3,5)。

由定理2可知,當Combine(p,q)為S的公共子序列,Combine(p,q)可以看作S上新的匹配,Combine(p,q)下標為p和q中最大的下標,所以有:

推論1 點p,q,r為S上的匹配,若Subscript(r)＞Subscript(q)＞Subscript(p),則(pqr)為S的公共子序列。

對于圖1中的s1和s2,(1,4)為字符e的匹配,(3,5)為字符p的匹配,(6,6)為字符t的匹配,合并后的序列(ept)為S的公共子序列。

定義3 設f(n)＝g(n)＋h(n),g(n)表示從初始節點到節點n付出的實際代價,h(n)表示從節點n到目標節點的最優路徑的估計代價,稱使用f(n)作為估價函數的GRAPHSEARCH算法為A。若算法A中使用的啟發函數h(n)對任何節點都有h(n)≤h?(n),則稱其為A?算法。

|MACS|表示多重序列(s1,s2,…,s d)的所有公共子序列數,d≥2,f(p)表示d個序列(s1,s2,…,s d)的|M ACS|的估計函數,g(p)表示d個前綴序列s1[1,…,p1],s2[1,…,p2],…,s i[1,…,p i],…,s d[1,…,p d]的|MACS|,g(p)可由公式(1)計算,h(p)表示d個后綴序列s1[p1＋1,…,n1],s2[p2＋1,…,n2],…,s i[p i＋1,…,n i],…,s d[p d＋1,…,n d]的|M ACS|,1≤i≤d,h(p)可由公式(2)計算,所以|M ACS|的估計函數為:

h?(p)表示d中任意兩個后綴序列s i[p i＋1,…,n i]和s j[p j＋1,…,n j]的|M ACS|的估計函數,1≤i,j≤d,h?(p)可由公式(3)計算。用h?(p)代替h(p)計算后綴序列的|MACS|,|MACS|的估計函數變成:

定理3 計算|MACS|的估計函數f?(p)是A?算法。

證明:由公式(7)可知,f?(p)＝g?(p)＋h?(p),f?(p)表示從初始節點出發經由節點p到達目標節點的所有公共子序列數的估計函數,g?(p)表示從初始節點到節點p的所有公共子序列數,h?(p)表示從節點p到目標節點的所有公共子序列數的估計函數。由定理1可知h(p)≤h?(p),根據定義2,f?(p)是A?算法,得證。

用MACS-A?表示計算|MACS|的估計函數f?(p),與A?算法查找最短搜索路徑不同,MACS-A?是在d維矩陣里計算d個序列的所有公共子序列數。M ACS-A?是A?算法的變形,在每次計算點p位置的所有公共子序列數時,最大化啟發函數h?(p)。h?(p)越大,包含的啟發信息越多,所需計算的節點越少,更快地找到目標節點,計算所有公共子序列數。若h?(p)＝0,點p到達目標節點,f?(p)＝g?(p)＋h?(p)＝g?(p)＝|M ACS|,所以有:

推論2 若h?(p)＝0,則g?(p)＝|M ACS|。

3 度量多圖相似度算法MACS-A?

通過計算多重序列的所有公共子序列數可以度量多重序列的相似度,根據將多圖表示為多重序列的圖深度優先搜索方法,度量多重序列的相似度可以度量多圖的相似度,所以,計算多重序列的所有公共子序列數可以度量多圖的相似度,度量圖的相似度問題簡化為計算所有公共子序列數的問題。多重序列的所有公共子序列數越大,多圖的相似度越大;反之,相似度越小。因此,本文提出了計算多重序列所有公共子序列數的啟發式算法MACS-A?,通過計算d個序列的前綴序列和后綴序列的所有公共子序列數度量d個圖的相似度,d≥2。

MACS-A?算法從初始點p0＝(0,0,…,0)開始,計算d個序列的所有單個字符的匹配,放在集合Q中。從表1可以看出,當一個匹配出現時,該匹配為多重序列的公共子序列,且該匹配與之前的公共子序列組成新的公共子序列。但對于表2圓圈中標記的匹配位置,該匹配已經出現,與之前的公共子序列組成新的公共子序列也已經出現,而且公式(2)在該匹配計算的所有公共子序列數并沒有變化。因為圓圈中匹配的下標不大于之前出現過的匹配的下標,所以當圓圈中的匹配出現時,不能組成新的公共子序列。對于圓圈中的匹配,在查找所有公共子序列和計算所有公共子序列數時,通過下標比較,不對該匹配進行合并序列、標記下標操作。

MACS-A?算法首先應用圖深度優先搜索方法將d個圖表示為d個序列;接著從Q中提取一個匹配p,T為d個序列的所有公共子序列的集合,比較p與Q、T中的序列的下標,合并序列后組成新的公共子序列,并對新的公共子序列標記下標,將不重復新的公共子序列放在T中;然后用公式(7)計算d個序列的所有公共子序列數,其中f?(p)表示d個序列的所有公共子序列數估計函數,g?(p)表示d個前綴序列的所有公共子序列數,h?(p)表示d個后綴序列的所有公共子序列數的啟發函數;最后返回d個序列的所有公共子序列和所有公共子序列數。

算法1:MACS-A?(s1,s2,…,s d)算法

輸入:圖G1,G2,…,G d

輸出:d個圖的所有公共子序列和所有公共子序列數

步驟1:應用圖深度優先搜索方法將d個圖表示為d個序列(s1,s2,…,s d),d≥2;

步驟2:初始p0＝(0,0,…,0),g?(p0)＝0,f?(p0)＝h?(p0),集合T表示d個序列的所有公共子序列,初始T＝{?},空集的下標Subscript(?)＝(0,0,…,0);

步驟3:計算d個序列的所有匹配,并放入Q中;

步驟4:|Q|表示d個序列的所有匹配數,從Q中取出一個匹配p,q表示Q中的任一匹配,t表示T中任一公共子序列,若|Q|＞0,轉到步驟5,否則,轉到步驟10;

步驟5:比較p和q的下標,若Subscript(p)＞Subscript(q),否則,轉到步驟6,合并p和q,Combine(q,p)＝(qp),并用p的下標標記(qp)的下標,Subscript(qp)＝Subscript(p),若(qp)不在T中,Label(qp)?Label(t),將(qp)放入T中,T＝T＋{qp},轉到步驟7;

步驟6:若Subscript(p)＜Subscript(q),合并p和q,Combine(q,p)＝(pq),并用q的下標標記(pq)的下標,Subscript(pq)＝Subscript(q),若(pq)不在T中,Label(pq)?Label(t),將(pq)放入T中,T＝T＋{pq},轉到步驟7;

步驟7:若p不在T中,label(p)?label(t),將p并入T中,T＝T＋{p},計算d個序列的所有公共子序列數,f?(p)＝g?(p)＋h?(p),轉到步驟8;

步驟8:比較p和t的下標,若Subscript(p)＞Subscript(t),將合并p和t的序列(tp)并入T中,T＝T＋{tp},標記(tp)的下標,Subscript(tp)＝Subscript(p),轉到步驟9;

步驟9:從Q中刪除p,Q＝Q–{p},轉到步驟4;

步驟10:輸出T和f?(p),T為d個序列的所有公共子序列集合,包括空集,f?(p)為d個序列的所有公共子序列數。

算法1包括兩部分:將d個圖G1,G2,…,G d表示為d個序列(s1,s2,…,s d),d≥2;計算d個序列的所有公共子序列數,并輸出所有公共子序列。為了方便分析算法1的時間復雜度,設d個圖的頂點數均為n,則d個序列(s1,s2,…,s d)的長度均為n。

算法1的步驟1應用圖深度優先搜索方法將頂點帶標號的有向圖G表示為序列s,時間復雜度為O(n＋|E|),其中n和|E|分別表示G的頂點數和邊數。圖深度優先搜索方法將d個圖表示為d個序列的時間復雜度為O(dn＋d|E|)。若|E|＝n×(n－1)＝(n2－n),則G為完全有向圖。最壞情況下,將d個圖表示為d個序列的時間復雜度為O(dn＋d(n2－n))＝O(dn2)。

算法1的步驟3為計算d個序列的所有匹配,從表1可以看出,建立2維矩陣M ij需要O(n2)的時間復雜度,1≤i,j≤d。從矩陣M ij可看出,計算2個序列中的所有匹配,需要O(n2)的時間復雜度,計算d個序列的所有匹配,需要O(dn2)的時間復雜度。

在算法1的步驟4中,Q為d個序列的所有單個字符匹配的集合,|Q|最大為n,n表示序列s i的長度,1≤i≤d。步驟5、步驟6和步驟8為比較p和Q中匹配的下標、p和T中的公共子序列的下標,時間復雜度為O(d)。步驟7為避免重復序列出現在T中,比較新組成的公共子序列與T中序列的字符,判斷新的序列是否存在T中,新的序列最長為n,最多為|Σ|個字符,Σ表示字母表,比較序列字符的時間復雜度為O(n|Σ|)。當新組成的公共子序列出現時,對p與Q中的每個匹配、T中的公共子序列進行合并和標記下標,其操作的時間復雜度為O(1)。所以,處理點p的時間復雜度總共為O(dn|Σ|)。步驟9從Q中刪除點p,當執行完步驟4,Q為空,所需的時間復雜度為O(dn2|Σ|)。

點p0從(0,0,…,0)到(n,n,…,n),當處理d個序列的每個匹配時,用公式(1)和公式(3)遞歸地計算d個序列在點p的所有公共子序列數f?(p)＝g?(p)＋h?(p),其時間復雜度也為O(dn2|Σ|)。結合步驟3的時間復雜度,用MACS-A?算法查找d個序列的所有公共子序列的時間復雜度為O(dn2＋dn2|Σ|),遞歸地計算d個序列的所有公共子序列數f?(p)的時間復雜度也為O(dn2＋dn2|Σ|)。

所以,結合步驟1的時間復雜度,用M ACS-A?算法度量d個圖的相似度的時間復雜度為O(dn2＋dn2＋dn2|Σ|)。

4 結束語

本文提出的啟發式算法MACS-A?通過計算多重序列的所有公共子序列數度量多圖的相似度,由于所有公共子序列數的變化都是在多重序列的匹配出現之后,所以MACS-A?算法遞歸地計算在匹配點上的所有公共子序列數,不必計算所有點的公共子序列數,避免了在非匹配點上冗余計算。該算法在處理匹配的過程中最大化后綴序列的啟發函數值h?(p),將訪問的匹配點p限制在矩陣M(公式(1))的匹配的子集,通過下標比較剔除不能組成新的公共子序列的匹配,進一步減少了計算節點的個數,能夠快速地度量多圖的相似度。

[1]Hattori M,Okuno Y,Goto S,et al.Development of a chemical structure comparison method for integrated analysis of chemical and genomic information in the metabolic pathways[J].Journal of the American Chemical Society,2003,125(39):11853-11865.

[2]Shanavas N,Wang H,Lin Z,et al.Supervised graph-based term weighting scheme for effective Text Classification[C]∥Proceedings of the 22nd European Conference on Artificial Intelligence,Hague,Netherlands,2016:1710-1711.

[3]Elzinga C,Wang H.Kernels for acyclic digraphs[J].Pattern Recognition Letters,2012,33(16):2239-2244.

[4]Sugiyama M,Llinares F,Kasenburg N,et al.Significant subgraph mining with multiple testing correction[C]∥Proceedings of the SIAM International Conference on Data Mining,Vancouver,Canada,2015:37-45.

[5]Li Tao,Dong Han,Shi Yong-tang,et al.A comparative analysis of new graph distance measures and graph edit distance[J].Information Sciences,2017,403:15-21.

[6]Zhu Gang-gao,Iglesias C.Computing semantic similarity of concepts in knowledge graphs[J].IEEE Transactions on Knowledge and Data Engineering,2017,29(1):72-85.

[7]Sugiyama M,Borgwardt K.Halting in random walk Kernels[C]∥28th International Conference on Neural Information Processing Systems,Montreal,Canada,2015:1639-1647.

[8]Borgwardt K,Kriegel H.Shortest-path Kernels on graphs[C]∥Proceedings of IEEE International Conference on Data Mining,Houston,USA,2005:74-81.

[9]Szilagyi S,Szilagyi L.A fast hierarchical clustering algorithm for large-scale protein sequence data sets[J].Computers in Biology and Medicine,2014,48:94-101.

[10]Forster D,Bittner L,Karkar S,et al.Testing ecological theories with sequence similarity networks:marine ciliates exhibit similar geographic dispersal patterns as multicellular organisms[J].BMC Biology,2015,13(1):1-16.

[11]Yanardag P,Vishwanathan S.A structural smoothing framework for robust graph comparison[C]∥Proceedings of Neural Information Processing Systems,Montreal,Canada,2015:2134-2142.

[12]Wang H.All common subsequences[C]∥Proceeding 20th International Joint Conference on Artificial Intelligence,Hyderabad,India,2007:635-640.

[13]Lin Z,Wang H,McClean S.A multidimensional sequence approach to measuring tree similarity[J].IEEE Transactions on Knowledge and Data Engineering,2012,24(2):197-208.

[14]Chin F,Poon C.Performance analysis of some simple heuristics for computing longest common subsequences[J].Algorithmica,1994,12:293–311.

[15]Wang Q,Korkin D,Shang Y.A fast multiple longest common subsequence(MLCS)algorithm[J].IEEE Transactions on Knowledge and Data Engineering,2011,23(3):321-334.