基于OpenCL的自動微分并行實現及其應用

(1.浙江東方職業技術學院 電氣自動化研究室，浙江 溫州 325035；2.溫州大學 電氣數字化設計技術國家地方聯合工程實驗室，浙江 溫州 325035；3.湖南大學 電氣與信息工程學院，湖南 長沙 410082)

0 引言

在很多的數值優化[1]方法中都會涉及到對函數導數，雅可比矩陣的計算。在傳統上，這些計算都是通過手工完成的。手工的計算方式對于很多小型的問題來說是比較簡單可行的。但是對于有眾多參數的復雜函數來說，手工計算其導數將是一個龐雜且容易出錯的過程。隨著現代計算機的迅速發展，對導數的計算都轉而采用計算機來實現。當前有三種通過計算機計算導數的方式。第一種為使用有限差分近似的數值微分方式。該方式簡單實用，但存在著較大的舍入和截取誤差。第二種方式為是使用符號計算的符號微分，該方式可以得到精確的封閉形式的求導結果。但存在一個嚴重的問題，那就是隨著函數的復雜度的增加，符號微分產生的符號表達式會呈指數級的增長。這就是所謂的表示式膨脹問題。而第三種方式則是本文所討論的自動微分[5-6]方式。自動微分克服了前兩種方式的缺點，會是數值計算方面主要的導數計算方式。

自動微分可以很有效地應用于計算機視覺中的大規模光束平差[2-3](Bundle Adjustment, BA)優化問題中。光束平差通過不斷地優化2D點的誤差來優化三維重建中相機和3D點的位置。采用的方法通常有Levenberg-Marquardt (LM)、Dog-Leg(DL)等[4]。這些方法都需要計算投影函數x=f(X)的導數來形成雅可比矩陣，其中x為2D點，X為其對應的3D點。由于3D點的數量眾多，可達到上百萬級。因此，對這一過程進行自動微分的并行化計算可以大大提高計算的速度。

自動微分有前向模式和反向模式兩種。對于含有多個參數的函數，從 計算效率上來說更傾向于采用后向計算模式。目前存在這些通用的用于計算自動微分的軟件，如ADOL-C[7]、CPPAD以及Ceres Solver。對于用于光束平差的Ceres Solver，其實現了一個基于OpenMP的自動微分。在本文中，則基于開放的OpenCL構架，提供了一個更高效的并行化自動微分實現。

1 自動微分的原理

相對于封閉形式求導結果或者近似的數值求導結果，自動微分可以獲得無截斷誤差的數值結果。這很大程度上得益于鏈式法則和計算機編程模式的相結合。每一個函數可以被分解為基本運算符的組合，其中基本運算符包括加、減、乘、除等二元算術運算符，以及像三角函數和指數函數等在內的超越函數。因此函數可以采用計算圖的形式來表示。例如考慮如下的函數：

(1)

式中，可以視為一個向量函數f:R2|→R2。其分解得到的計算圖見0。圖中每個運算都由一個帶編號的節點表示。可以通過前向和反向兩個模式來計算f關于x的導數。

圖1 公式(1)的計算流程圖

1.1 前向計算模式

令vi-n=xi,i∈[1,n]為輸入，vi,i∈[1,l]為中間變量，ym-i,i∈[m-1,0]為函數fk,k=m-i的輸出。一個三階段的記法可以用來對導數的前向計算進行形式化，自動微分前向計算模式的計算過程如下所示：

(2)

(3)

(4)

式中，關系ji表示vi與vj是直接相關的。在計算圖也表示節點vi是vj的一個直接后繼。φi是描述基本運算的函數。ui=(vj)ji表示vi的所有前驅的集合。令所有輸入為一個向量x=[x1,…xi,…,xn]。若將輸入的導數設定為則就為希望計算得到的導數?fk/?xi。前向計算模式比較簡潔易懂，并且對于f:R|→Rm這樣的標量函數來說，其計算是非常有效的。因為輸入參數只有一個，那么只需要令便可以計算出所有的導數值。但是對于像f:Rn|→R或f:Rn|→Rm這樣的向量函數，雙述的計算過程必須運行n次或n×m次。而反向的計算模式則可以避免這個問題，提高計算的效率。

1.2 反向計算模式

表1 自動微分反向計算模式的計算過程

該計算過程分為兩部分。第一部分像正向計算模式一樣計算出函數值。第二部分則反向地估計函數對所有輸入參數的導數。對于函數f:Rn|→Rm，所示的計算過程只需要運行m次便可以獲得所有的求導結果，從而得到最終的雅克布矩陣。

2 并行實現

本節中實現的基于OpenCL[8-9]的自動微分采用反向計算模式，可以用于“large-small”問題。所謂的“large-small”問題，就是單一函數的計算圖并不復雜，但是卻存在著大量的重復的計算。光束平差問題便是一個例子。由于大量的3D點和相機的存在，就需要計算大量的投影函數對3D點坐標和相機參數的導數，來形成最終的稀疏雅克布矩陣。基于OpenCL的自動微分的實現分為主機端和設備端兩部分。在主機端運用了C++的函數重載和模板特性[10]來有效地生成計算圖。而設備端則根據計算圖并行地計算出求導結果。

2.1 主機端編程

為了簡化使用的方式，對待求導函數的構建應該盡量的趨近于原生的C/C++風格的代碼編寫。式的函數可以寫成如下的簡潔形式：

DVAR x1,x2;

DVAR f1=ln(x1)+x1*x2-sin(x1)

DVAR f2=x1*x2

要實現如此的函數構建的簡單化，需要完成一些關鍵的任務。首先，定義一個用于描述節點的結構體：

template

struct ADV_Node {

OpType op;

shared_ptr> arg[2];

T val;

T dval;

int id;

}

在該結構體中，OpType是一個用于指示運算符的枚舉類型，指示該節點為其前驅進行某種運算的結果。除了常規的運算之外，還引入“CONST”和“INPUT”來分別表示該節點是否為常量或者函數的輸入參數。每個節點在生成時還被賦予一個唯一的id值。arg指向該節點的一個或兩個前驅節點，指示該節點運算符的操作數。dval用于在求導過程中存放函數對該節點的導數值。

然后，實現了一個Wrapper類ADV來描述最終的數學形式上的變量，并用此來完成數學表達式的構建。ADV類的定義如下：

template

class ADV {

public:

shared_ptr> ADVptr;

ADV();

ADV(shared_ptr> ptr);

ADV(const ADV &adv);

ADV(const T val);

ADV& operator=(const ADV &rhs);

ADV& operator=(const T &val);

}

采用智能指針來為每個ADV變量創建節點可以使得ADV變量能夠不受作用域的限制，像數學形式上的函數那樣完成代碼中的對應函數。例如要實現一個3元素向量的點積，可以寫為“dot3(ADV *a, ADV*b)”。該函數能返回一個新建的ADV變量來記錄所有的運行結構。ADV類中所有的構造函數和賦值運算符重載都用來確保變量的正確生成。所有的關于ADV變量間的運算都用過類定義外的友元的函數重載來實現。每個函數重載創建新ADV變量來記錄這一關于其前驅的運算。例如加法的函數重載實現為：

ADV operator+(const ADV &x,

const ADV &y)

{

ADV adv;

adv()->val=y()->val+x()->val;

adv()->op=ADD;

adv()->var[0]=x.ADVptr;

加大依法治林力度，提升森林資源管護水平，全市將以執行《煙臺市森林防火條例》為總抓手，突出抓好以防火道路網、引水上山水網、預警監測網、森林消防專業隊伍、專職護林隊伍為主體的“三網兩隊”建設，大力構建森林火災應急處置信息化、撲救隊伍專業化、設施裝備現代化體系，開展嚴厲打擊非法野外用火行為專項行動，全面提升森林防火能力和法制化水平，確保全市森林防火形勢持續保持平穩。同時，認真貫徹執行相關法律法規，進一步規范涉林審核審批，持續開展嚴厲打擊亂征濫占林地、亂砍濫伐林木等違法行為專項行動，積極運用飛防、地面噴灑藥物、釋放天敵生物等多種措施進行綜合防治，確保林業有害生物實現持續有效控制。

adv()->var[1]=y.ADVptr;

return adv;

}

程序上的“c=a+b;”表達式便可以直接地表示數學上的c=a+b。通過使用不同運算的函數重載來構造函數，函數所對應的計算圖也同時被構建出來。另外對于大多數情況下使用的雙精度浮點型來說，可以將ADV定義為DVAR。

2.2 計算序列的生成

在構建函數以及對應的計算圖之后，便可以據此運用反向計算模式來進行導數的計算。為了正確地計算各個節點的導數值，同時也是為了便于設備端的并行計算。需要將計算圖轉化為計算序列。

在反向計算模式下，計算序列又分為正向計算序列和反向計算序列，正向計算序列用于計算函數值，而反向計算序列則用于計算函數的導數值。由于每個節點的id在生成時滿足后繼節點的id值一定大于其前驅節點的id值。所以用于計算函數值的正向計算序列可以簡單地取為各節點的升序排列。然而在確定反向計算序列時，必須考慮節點之間的依賴關系。如0所示，如果先從節點7開始處理，在處理節點2時，其所依賴的節點5的導數值還并沒有被計算出來。

圖2 編程模式下的計算圖及正向和反向計算序列

因此，采用拓撲排序來完成反向計算序列生成。所生成的反向計算序列中各節點滿足相互間的依賴關系。拓撲排序的偽代碼見算法1：

算法1:用于生成反向計算序列拓撲排序

輸入：計算圖G=(V,E)

輸出：反向計算序列L

初始化計數數值C, 用于記錄各節點的入度

forvinVdo

令s,t為v的前驅節點

增加s,t在C中的入度

end for

whileL非滿 do

找到具有零入度的節點v

添加到L尾部

令s,t為v的前驅節點

減少s,t在C中的入度

end while

2.3 設備端并行計算

當在主機端完成計算序列的生成之后，便可以將需計算的函數的參數傳入設備端，在設備端按照計算序列計算出函數值以及函數的導數值。

在OpenCL下，多個kernel函數并行地在如CPU、GPU或者FPGA設備上運行，其每個運行的實例稱為一個工作項。其并行的數量取決于設備上計算單元的數量，以及計算單元上局部存儲器大小等因素。主機端需要指定全局工作項的數量，其對應到需要進行求導的函數的數量。同時也可以指定局部工作項的數量，形成工作組。每個工作組中的工作項以SIMT(Single Instruction Multiple Thread)的模式，并且共享一組數量有限的高速局部存儲器。在像光束平差的應用中，所有的導數計算使用相同的計算圖和序列。如果計算圖的尺寸比較小，那么可以將生成的計算序列傳輸到局部存儲器以訪問提高速度。

用于執行每個實際求導過程的kernel函數相對比較簡單。首先按照給定的正向計算序列和函數參數值計算出函數值。然后再按照反向計算序列依次計算每個節點的導數值。最后將作為輸入節點的導數返回給主機端。該計算過程的偽代碼見算法2：

算法2:計算導數的Kernel函數for i=0 to size(forward_seq) do

compute_val(node_op,arg1_val,arg2_val)

end for

for i=0 to size(input_args) do

val_out[base+i]=node_val[i]

end for

for i=0 to size(funcs) do

for j=0 to size(reverse_seq)

compute_diff(node_op,node_val,

arg1_diff,arg2_diff)

end for

for j=0 to size(input_args) do

diff_out[base+i*N+j]=node_diff[j]

end for

end for

3 應用實例

四元數可以表示為q=。在相機坐標系下對3D點的旋轉可以表示為：

Rot(X)=q·q(X)=q⊕q(X)⊕q-1

(5)

其中:⊕為Hamilton積，q(X)則用于將X變為四元數形式<0,X>。在Rodrigues參數形式中，旋轉表示為繞單位向量k的旋轉，旋轉角度為θ。因此，X的旋轉在Rodrigues參數形式下表示為：

Rot(X)=Xcosθ+(k×X)sinθ+

(1-cosθ)(k*v)k

(6)

通過將旋轉與位移的結合，可以得到3D點在相機上的投影。在光束平差中，視3D點坐標X和相機矩陣[R|t]為投影x的參數。采用Rodrigues參數形式，編寫如下的代碼完成x相對于X和[R|t]的計算圖的構建：

DVAR x,y;

DVAR theta2=dot(angle_axis, angle_axis);

DVAR theta=sqrt(theta2);

DVAR costheta=cos(theta);

DVAR sintheta=sin(theta);

DVAR theta_inverse=1.0/theta;

DVAR w[3]={ angle_axis[0]*theta_inverse,

angle_axis[1]*theta_inverse,

angle_axis[2]*theta_inverse

};

DVAR w_cross_pt[3];

cross(w,pt3D,w_cross_pt);

DVAR tmp=dot(w,pt3D)*(1.0-costheta);

DVAR tmp3D[3];

tmp3D[0]=pt3D[0]*costheta+

w_cross_pt[0]*sintheta+w[0]*tmp;

tmp3D[1]=pt3D[1]*costheta +

w_cross_pt[1]*sintheta+w[1]*tmp;

tmp3D[2]=pt3D[2]*costheta +

w_cross_pt[2]*sintheta+w[2]*tmp;

tmp3D[0]=tmp3D[0]+transl[0];

tmp3D[1]=tmp3D[1]+transl[1];

tmp3D[2]=tmp3D[2]+transl[2];

x=focal*tmp3D[0]/tmp3D[2];

y=focal*tmp3D[1]/tmp3D[2];

其對應的計算圖見圖3。該計算圖隨后被轉化為計算序列，并伴隨所有3D點坐標和相機參數送入設備的并行計算。計算出投影點坐標，以及投影點作為函數對3D點坐標和相機參數的導數。最終獲得雅可比矩陣的每個塊Aij和Bij。Aij為投影點坐標(x,y)對3D點坐標(X,Y,Z)的導數。Bij為投影點坐標對向量化的相機參數p的導數。

圖3 Rodrigues參數形式下投影函數的計算圖

4 測試與分析

本節中對實現的并行化自動微分進行測試。將其應用到BAL(Bundle Adjustment in the Large)問題[12]中。所選擇的數據集見表表3。作為對比，同時引入了對Ceres Solver在該問題上的測試。Ceres Solver是Google公司推出的用于解大型數值優化問題的開源軟件。其中的自動微分基于OpenMP多線程框架，在CPU上實現了線程級別的并行計算。

表2 各數據集的主要數據

測試的平臺為一臺兼容PC機。采用的處理器為Intel XEON E5 2643 v2，其共有12個超線程核心，運行頻率為3.2 GHz，配備的內存為1866 MHz的8GB DDR3 RAM。 用于測試OpenCL并行計算的GPU采用AMD的R9 290，其共有40個計算單元，4GB顯存，核心和顯存的工作頻率分別為945 MHz和1 240 MHz。R9 290能夠提供很強大的并行計算能力，其單精度浮點計算性能可達4848 GFLOPS，而雙精度浮點計算能力也可達到606 GFLOPS。對于大多數科學計算來說，例如使用光束平差法的優化問題，通常都使用雙精度浮點數據。

圖4 工作組尺寸對性能的影響

首先測試工作組的大小對計算的性能的影響，采用的數據集為Dubrovnik。工作組的大小從4變化到128。在CPU和GPU上各運行100次取均值，得到的結果見0。從結果可以看出，CPU運行所需時間普遍小于GPU所需時間。隨著工作組尺寸的增長，GPU計算的性能會逐漸的下降。其主要原因在于GPU上工作組內存儲器的訪問沖突。GPU上的存儲器按較大地址塊(如1K字節)的形式訪問，同一個塊上的地址必須訪問。而在CPU上該問題得到很大的緩解，因為CPU上的地址塊大小要小很多。

在對所有的數據進行測試時，結合對工作組尺寸的分析，針對CPU采用的工作組尺寸為32，而針對GPU采用的工作組尺寸為8。所獲得的測試結果見0。使用同樣的CPU，本文基于OpenCL的方法比Ceres Solver基于OpenMP的方法在速度上快了大約3.6倍。而基于GPU的實現比Ceres Solver快了1.6倍。由此可見，CPU實現比GPU實現的性能要好。盡管R9 290的GPU有40個計算單元，但其只運行在1 GHz。而XEON E5的CPU運行在3.2 GHz。同時，自動微分中存在著許多的轉移分支，這并不利于GPU發揮其流水線的特性。相反，CPU是被設計為執行復雜任務的，有著強大的分支預測能力。因此，大規模的自動微分更適合在多核的CPU中執行。

表3 自動微分在各數據集上的測試結果

5 總結

本文首先展示了自動微分系統的工作原理，包括正向模式和反向模式。對于多參數的函數，揭示了反正模式比正向模式具有更高的效率。以OpenCL為并行計算框架，實現了針對大型優化問題的并行化自動微分。該實現采用反向計算模式，以C/C++的風格構建函數，并生成計算圖和計算序列。以光束平差為應用背景，通過測試，該實現比Ceres Soler快約3.6倍。同時可以發現，自動微分在CPU上的實現要優于在GPU上的實現。