基于浮點數(shù)的Cholesky 分解FPGA 實現(xiàn)*

朱 鵬 葉樹霞 楊曉飛

（江蘇科技大學(xué)電子信息學(xué)院 鎮(zhèn)江 212003）

1 引言

求解逆矩陣的方法有許多，如高斯消元，矩陣分解等方法。通過矩陣分解實現(xiàn)矩陣求逆的運算量小且易于實現(xiàn)，在實際工程中被廣泛運用。常用的矩陣分解方法有Cholesky分解、LU 分解和QR 分解。其中，Cholesky 分解適用于共軛對稱正定矩陣；LU 分解適用于順序主子式不為0 的任意矩陣；QR分解適用于任意矩陣。針對協(xié)方差矩陣的共軛對稱厄米特（Hermitian）正定的特性，使用Cholesky分解可以節(jié)約FPGA資源使用，提高運算速度［1～3］。

本設(shè)計采用了浮點數(shù)來實現(xiàn)Cholesky 分解，相比定點數(shù)浮點數(shù)在運算過程中無需考慮截位效應(yīng)帶來的精度損失，大大提高了結(jié)果的精度。在計算浮點數(shù)開方運算上，設(shè)計了基于平方根倒數(shù)方法的硬件電路，相比傳統(tǒng)的牛頓迭代法，大大減少了迭代次數(shù)。

2 Cholesky分解

Cholesky 分解矩陣方法是利用協(xié)方差矩陣R厄米特（Hermitian）正定的特性，將協(xié)方差矩陣R分解為上/下三角矩陣G及其共軛矩陣的乘積。通過求取上/下三角矩陣的逆矩陣的共軛矩陣GH及矩陣G的乘積，即可得到原協(xié)方差矩陣R的逆矩陣［4～5］。

R是共軛對稱厄米特（Hermitian）正定，R=G*GH是矩陣的Choleksy 分解，其中G是一個具有正的對角線元素的下三角矩陣［6～7］。

由式（1）可得：

整理式（2）可得遞推式：

3 平方根倒數(shù)

平方根倒數(shù)方法（Fast inverse square root），經(jīng)常和一個十六進制的常量0x5f3759df 聯(lián)系起來。該算法被用來快速運算平方根倒數(shù)，在20 世紀90年代末由硅圖公司開發(fā)出來，后見于John Carmark的Quake III Arena的源碼中。

首先在IEEE 754 標準下的單精度浮點數(shù)表示如圖1。

圖1 單精度浮點數(shù)表示格式

由已知：

其中，ex表示浮點數(shù)的指數(shù)大小，mx表示浮點數(shù)的余數(shù)大小。

定義：Sx為符號位，當x 為正數(shù)是0，負數(shù)為1。Ex=ex+B為偏移指數(shù)，B=127 。Mx=mx+L，L=223。

在實數(shù)域中，被開方數(shù)為正數(shù)，則Sx=0。對式（4）作以2為底的對數(shù)運算得到：

因為mx∈[0，1)，所以有

其中，α為調(diào)節(jié)精度近似值。當α=0.0430357 時，可得最優(yōu)結(jié)果。

根據(jù)式（5）和式（6）：

這個數(shù)可以表示為

由式（8）可得：

平方根倒數(shù)方法最重要的一步就是找出了一個迭代初值，推導(dǎo)出的初值越接近真實結(jié)果，算法也就收斂越快。

將式（8）的結(jié)論帶入得

將數(shù)據(jù)帶入式（11），并轉(zhuǎn)為十六進制的：

在得到迭代初值后，只需要進行牛頓迭代法就能得到平方根倒數(shù)。

迭代公式如下：

綜合式（13），式（14）得到：

在實際工程中，由于加入了平方根倒數(shù)的方法求初值，一次迭代就可以基本所需要的平方根倒數(shù)。

4 Cholesky分解模塊設(shè)計

實現(xiàn)Cholesky 分解模塊，需要一個數(shù)據(jù)控制模塊和一個Cholesky分解計算模塊，如圖2。

圖2 Cholesky分解模塊總體結(jié)構(gòu)

由于協(xié)方差矩陣為共軛對稱正定，所以FPGA的ram 只需要存儲一半元素和主對角線元素。數(shù)據(jù)控制模塊儲存這些數(shù)據(jù)，并按要求發(fā)送給Cholesky分解計算塊，同時接收和儲存Cholesky分解出來的G 矩陣。數(shù)據(jù)控制模塊數(shù)據(jù)的發(fā)送接收流程，如圖3。

圖3 數(shù)據(jù)控制模塊工作流程

4階Cholesky分解計算模塊，如圖4。

圖4 4階Cholesky分解計算模塊結(jié)構(gòu)示意圖

因為Rsqrt 模塊為求平方根倒數(shù)模塊，所以主對角線元素g11，g22，g33，g44為是實際值的倒數(shù)，但是在Cholesky分解之后所作的高斯消元法、下三角矩陣求逆都需要的主對角線的倒數(shù)，可以節(jié)約后續(xù)開發(fā)資源使用［8～11］。由于每一列的Cholesky分解計算，每列只有首位元素需要求解平方根倒數(shù)，為了節(jié)約FPGA的資源，只需要一個Rsqrt模塊。

如果更高階的Cholesky 分解，只需要例化對應(yīng)的PE 模塊，并設(shè)置好對應(yīng)位置的參數(shù)，將其于Rsqrt模塊相連即可。

其中PE1 只需要計算平方根倒數(shù)為一個特殊節(jié)點，其余PE結(jié)構(gòu)如圖5。每一個PE節(jié)點，包含一個復(fù)數(shù)的浮點數(shù)乘法器、一個復(fù)數(shù)的浮點數(shù)加法器和存儲數(shù)據(jù)的Ram。為了節(jié)約資源，當需要作減法時，只需要將減數(shù)的符號位取反，送入浮點數(shù)加法器。存儲數(shù)據(jù)的Ram也是基于參數(shù)化配置，根據(jù)式（3），靠后的PE 節(jié)點需要儲存前面PE 節(jié)點計算的結(jié)果，所以當前PE 節(jié)點需要前一個PE 節(jié)點加一的Ram 資源。復(fù)數(shù)的浮點數(shù)乘法器和復(fù)數(shù)的浮點數(shù)加法器都為流水線型，提高運算效率。

圖5 PE的內(nèi)部結(jié)構(gòu)

Rsqrt 模塊結(jié)構(gòu)如圖6。綜合式（11），式（15）U1為一個預(yù)處理單元，其中，U1_y=0x5f3759df-(Fx?1)。利用浮點數(shù)的性質(zhì)，求F 的一半，只需要對其階數(shù)加1。U2為一個定制單元用來計算U2_x= .5*U1_x。U3也是一個定制單元用來計算U3_x=U1_x*U1_y*U1_y。最后U4為一個浮點數(shù)減法器F_rsqrt=U2_x-U3_x。Rsqrt 模塊是一次迭代的平方根倒數(shù)方法實現(xiàn)，流水線型實現(xiàn)，4 個時鐘周期就可以得到輸入值的平方根倒數(shù)。

圖6 Rsqrt模塊的結(jié)構(gòu)

5 結(jié)果分析

完成算法設(shè)計之后，對算法的可行性進行驗證。本文所使用的開發(fā)環(huán)境為Vivado 2018.3，使用的硬件編程語言為Verilog HDL，實驗平臺是Xilinx的Xc7k325t。

實驗對4 階、5 階、6 階、7 階和8 階的Cholesky分解模塊進行了仿真測試。每次實驗通過Matlab［12］隨機生成20 組信號協(xié)方差矩陣，將生成的數(shù)據(jù)通過coe 文件保存到FPGA 的rom 里，再將FPGA 運算好的數(shù)據(jù)通過串口讀出。FPGA運算的數(shù)據(jù)與Matlab 運算數(shù)據(jù)進行比較，結(jié)果如表1 所示。運行周期為1組數(shù)據(jù)需要的運行時長，結(jié)果精度為20組數(shù)據(jù)的平均精度［13］。

表1 實驗結(jié)果

在Vivado 中，通過時序分析PE 模塊在4 階的情況下主頻率可達186MHz，隨著矩陣階數(shù)的提高會逐漸降低。實驗均在150MHz的主頻下進行。

6 結(jié)語

本文提出了一種協(xié)方差矩陣的Cholesky 分解的并行算法和相應(yīng)的硬件結(jié)構(gòu)，并使用FPGA 進行了實現(xiàn)與驗證。相比傳統(tǒng)Cholesky 分解實現(xiàn)方法使用定點數(shù)相比，本文采用的浮點數(shù)精度更高，不會丟失小數(shù)位便于后續(xù)處理。采用模塊化，通過PE 節(jié)點配置參數(shù)就能快速實現(xiàn)更高階的Cholesky分解算法，相比全并行結(jié)構(gòu)在高階配置和資源消耗更有優(yōu)勢［14～15］。