基于滑動窗的Needleman-Wunsch 算法硬件加速方案設計與實現

肖俊杰

（上海交通大學電子信息與電氣工程學院，上海200240）

0 引言

生物序列比對是生物信息學和計算生物學的基本問題之一，它旨在找出一個未知序列與已知序列在結構和功能上是否相關，這對于疾病的早期防治和藥物工程至關重要[1]。Needleman-Wunsch 算法[2]是一種廣泛使用的序列比對算法，用于全局比對兩個序列之間的相似性。這一算法基于動態規劃實現，可以求得比對的最優解，但其時間復雜度與序列長度的乘積成正比，考慮到生物序列數據庫的巨大規模和其指數增長率，需要更快的計算技術來應對這一問題。

高算力的并行結構為解決這一問題提供了機會，但是Needleman-Wunsch 算法內部數據之間存在依賴性，需要各計算節點之間相互通信。現有計算方法大多基于波前進位（wave-front）[3]這一并行方式進行優化刪減[4-6]，但其訪存特征不規則，難以實現高效的訪存調度，有限的存儲資源和過多的片外存儲器訪問導致并行效率下降成為了硬件優化的主要阻礙。因此本文提出了一種基于滑動窗的解決方案，采用了通用化的并行策略，簡化了對Needleman-Wunsch 算法進行并行的難度；并提出一種延時串連的脈動陣列（systolic array）[7]結構，減少了對存儲的訪問。為與不同加速方案下的計算和存儲特征相匹配，本文基于粗粒度可重構陣列（CGRA）[8]對改進前后的方案進行了比較，并與通用計算架構下的解決方案對比，量化本文方案在執行時間和執行效率方面的優勢。

1 Needleman-Wunsch算法

Needleman- Wunsch 算法的目標是在序列X=x1x2…xi…xM（序列長度為M）和序列Y=y1y2…yj…yN（序列長度為N）之間尋找最佳的全局匹配結果。該算法建立了一個(M+1) ×(N+1) 的分數矩陣F，其中矩陣的元素F(i,j)表示X的子序列x1x2…xi和Y的子序列y1y2…yj比對的得分，分值越高表示兩條序列的相似性越大。矩陣的最后一個元素F(M,N)為序列X和序列Y的最優比對得分，最優比對結果通過回溯得出。該算法分為初始化、填充和回溯三個步驟，如圖1 給出的兩條示例序列所示。

分數矩陣的第一行和第一列被初始化，其余元素的值由其左元素、上元素和對角線元素的值確定，其填充規則如下：

上述式子表明了位置(xi,yi)的3 種更新方式：通過（xi-1,yj）的垂直方向，此時序列X需要插入空位，罰分為-d；通過（xi,yj-1）的水平方向，此時序列Y需要插入空位，罰分為-d；通過（xi-1,yi-1）的對角方向，此時序列匹配成功，得分為s（xi,yi）。

圖1 Needleman-Wunsch算法執行步驟

回溯過程從F(M,N)開始，如圖1 所示，以箭頭標記的方向追溯路徑來源，通過箭頭的方向確定是否需要在序列中插入空位，獲得具有最優解特性的匹配結果。

Needleman-Wunsch 的填充階段本質上是二維動態規劃數組的建立階段，其時間和空間復雜度均為O(M×N)，是整個算法計算、訪存最密集的部分。又由于動態規劃的特點，需要各個節點之間相互通信，這一特性增大了Needleman-Wunsch 算法并行化的難度，使其只能以對角線條帶（diagonal strip manner）的形式進行內存訪問，這種不規則的訪問方式破壞了數據的空間局域性。目前的加速平臺通常采用一級或多級高速緩存（cache）來提升訪存的性能，但cache 的容量有限，單純增大并行度并不能顯著提升系統的性能，反而可能降低系統的存儲帶寬和增大緩存缺失率。采用合適的并行計算方法提高計算和訪存的效率成為了加速Needleman-Wunsch 算法的關鍵。

2 加速方案設計

2.1 并行特征提取

依據公式（1），Needleman-Wunsch 算法在執行比對任務時內部的數據存在依賴關系，即元素F(i,j)的計算依賴于相鄰三個元素：F(i,j-1)、F(i-1,j)和F(i-1,j-1)，只有這三個元素計算完成后才開始計算F(i,j)。因此通常的串行方式是按照矩陣的行或者列逐層進行運算，如圖2（a）所示。這種方法在一個時間步內只能處理矩陣的一個元素，計算效率低，而且下一個節點的計算操作數來自前一節點的計算結果，不能簡單地對矩陣不同位置的元素進行并行展開。

圖2 分數矩陣填充步驟的實現和優化

但不難發現當矩陣某一元素F(i,j)的計算完成后，其右方的元素F(i,j+1)和下方的元素F(i+1,j)就可以開始計算，即位于矩陣同一條反對角線上的元素不存在數據相關，從矩陣左上角的第一個元素開始沿著對角線方向通過波前進位的方式進行并行處理，就能夠在一個時間步內同時處理同一條反對角線的元素，如圖2（b）所示。針對波前進位的優化主要是通過不同的分塊策略實現計算和訪存之間的平衡，如圖3 所示，但對角線元素數量不一使得計算和訪存依舊是不規則的，并行效率不高。

圖3 分數矩陣的分塊

圖4 滑動窗數據流圖

2.2 并行策略改進

本文利用串聯延時脈動陣列結構對上述策略進行了改進，提出了一種滑動窗式的并行加速方法，不僅考慮了單個應用的特性，還提取出算法的公共特征，得到的計算結構更易于擴展。如圖2（c）所示，硬件在每一個時間步內只需要處理單個滑動窗內的元素，每個滑動窗遵循相同的計算規則，并且以邊界填充（padding）的形式避免邊界條件的判斷并保證計算規則的統一。計算規則確保了每個滑動窗內部計算的并行性，滑動窗的大小可以依據片上存儲器的容量和計算單元的數量確定。這一方法的好處是滑動窗的計算規則保證了單個滑動窗內的數據不存在相關性，硬件只需要對單個滑動窗的計算進行并行加速。并且邊界填充可以在算法的初始化階段完成，簡化了邊界條件的判斷，提高了資源的利用效率。

以圖2（c）為例，假定滑動窗大小為4×4，滑動窗在單個時間步內可同時計算3 個元素的值，在分數矩陣的邊界處（圖中的t=0 和t=7）額外填充了2 列元素，使得滑動窗可以使用相同的規則進行計算。滑動窗首先進行橫向滑動，在計算得到當前位置的結果之后向右滑，并將之前的結果作為新位置下的輸入，按照此方式處理分數矩陣的0 至3 行后向下滑，以相同的規律處理矩陣的3 至7 行。滑動窗的滑動過程很容易通過脈動陣列結構實現，脈動陣列結構可以在消耗較小的內存帶寬的情況下實現較高的運算吞吐率。而且這種計算方式容易擴展，能夠用于不同規模的輸入。

這一方法得到的數據流圖如圖4 所示，在硬件概念上使用計算處理單元（Processing Element，PE）對滑動窗內的元素實現并行計算。因為某一次計算得到的結果會在接下來兩次計算中使用，所以計算得分的PE會將結果直接或經過延時后轉發給另一個PE，實現局部的數據復用，PE 陣列只有在分數矩陣的上邊界和左邊界處才需要從存儲器中讀取操作數，其余時間可以從計算得分結果的PE 處獲取。單純采用對角線展開的方式處理3 個元素需要進行7 次讀操作，相比之下采用滑動窗實現只需要2 次讀操作，極大地減少了帶寬需求。同時，在計算得到矩陣元素分值的同時也可以知道該矩陣元素是由哪個方向的元素得到的，因此回溯方向的保存可以和分值的計算同時進行。此外由于滑動窗是按照矩陣的行、列方向滑動的，在降低了訪存次數的同時也降低了對角線條帶式訪問帶來的緩存缺失問題。

3 實驗結果與分析

3.1 實驗平臺

本實驗中，我們采用了粗粒度可重構陣列（CGRA）結構來實現滑動窗的結構。粗粒度可重構陣列兼具靈活度和高性能，能夠與不同的訪存和計算結構相適配，從而對不同的算法提供廣泛而靈活的支持。本文基于C++搭建了一個包含典型CGRA 結構的系統級模擬器，用以模擬可重構陣列的執行過程，模擬器包含周期精確的陣列模型和的內存模型[9]，該CGRA 系統結構如圖5 所示。

圖5 粗粒度可重構陣列系統結構

其中，粗粒度可重構陣列負責處理核心計算，主處理器對系統進行控制，對目標應用進行編譯生成帶有配置信息的數據流圖（Data Flow Graph，DFG），任務控制器對數據流圖進行解析，隨后向PE 組成的PE 陣列發送配置信息，PE 陣列依據配置完成任務的執行。本地存儲器和主存儲器構成存儲器層次結構，實現高效的數據存儲和訪問功能，本地存儲器采用16KB 的一級cache。

PE 陣列包含有8×8 規模的64 個PE，每個PE 可依據配置字配置成不同的功能，32 個存儲處理單元（Load Store Element，LS）分布在PE 的外圍負責管理訪存請求。每個PE 可以與同行同列的PE 或LS 互連，這種互連方式支持同行和同列的PE 之間相互轉發數據，從而實現脈動陣列，大量PE 組成的陣列帶來了高并行性，并且PE 計算的結果不需要存回存儲器，可以直接轉發給下一PE 進行處理，提高了計算效率。

3.2 實驗結果分析

為了驗證改進后的方案在CGRA 上實現的效率，本文對上一節的三種填充實現策略進行了映射，分別是串行、波前進位和滑動窗式，堿基序列長度均為256，最終得到的實驗結果的如表1 所示。表中的PE 數量表示執行應用所需要的PE 總數，表明硬件的資源消耗情況；PE 利用率表示PE 執行有效計算的時間占執行時間的比率，表明硬件的執行效率。

其中方案1 不進行并行、分塊等任何優化，作為對比的基準，方案1 中的許多PE 被用于循環控制和邊界判斷，因此利用率不高。方案2 將矩陣分成8×8 大小的塊，采用對角線條帶的并行方式，單個塊最多可同時計算8 個對角線元素，增大了計算的并行度。但是大并行度也導致了更高的緩存缺失率（cache miss rate），多個請求相互爭搶cache 也導致了訪存沖突（cache conflict），不規則的對角線計算也使得部分PE 閑置，PE利用率反而有所降低，方案2 相比方案1 最終獲得了3.5 倍的加速比。方案3 選擇的滑動窗大小為9×9，最多可同時計算8 個對角線元素；由于計算更規則，在相同并行度下減少了用于邊界判斷的PE 數量；方案3 相比方案1 的加速比提升到了7.5 倍，而且由于方案3 采用脈動陣列實現，降低了訪存次數，很大程度上減少了由于緩存訪存缺失和訪存沖突引起的流水線停頓，相比方案2 獲得了2.3 倍的性能提升。

表1 不同映射方式執行時間對比

表2 不同架構執行時間對比

表2 給出了本文方案與CPU 和GPU 對比的結果，以驗證本文的上述方案相對通用加速平臺的性能提升。其中CPU 測試平臺采用Intel Core i7-4770 處理器，16GB DDR3 內存，以執行時間作為性能標準，本文提出的方案獲得了4.5 倍的加速比。GPU 平臺采用NVIDIA GeForce GTX 1080，并使用Rodinia[10]中的加速方案。由于Needleman-Wunsch 算法需要計算節點之間相互通信，GPU 不同線程之間需要頻繁同步，因此整體的并行效率不高，本文方案在使用更少計算資源情況下依然獲得了15.3%的性能提升，執行效率更高。

4 結語

本文提出了一種基于滑動窗的解決方案來加速Needleman-Wunsch 算法，并基于CGRA 對該方案進行了驗證。相比傳統波前進位的并行方式，該方案提取了算法的公共計算特征，簡化了算法在硬件上的實現難度，使得更多的硬件資源可以用于核心計算；采用串連延時脈動陣列結構降低了訪存次數，減少了由于訪存延時帶來的流水線停頓，與改進前相比獲得了的2.3倍加速比，相對CPU 獲得了4.5 倍的加速比，相對GPU，在使用計算資源更少的情況下仍然獲得了15.3%的性能提升。值得一提的是，本文提出的方案核心是對二維動態規劃矩陣的建立過程進行硬件優化，這種優化方法易于拓展到其他具有二維動態規劃特性的算法上。