基于OpenCL的自動微分并行實(shí)現(xiàn)及其應(yīng)用

(1.浙江東方職業(yè)技術(shù)學(xué)院 電氣自動化研究室，浙江 溫州 325035；2.溫州大學(xué) 電氣數(shù)字化設(shè)計(jì)技術(shù)國家地方聯(lián)合工程實(shí)驗(yàn)室，浙江 溫州 325035；3.湖南大學(xué) 電氣與信息工程學(xué)院，湖南 長沙 410082)

0 引言

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

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

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

1 自動微分的原理

相對于封閉形式求導(dǎo)結(jié)果或者近似的數(shù)值求導(dǎo)結(jié)果，自動微分可以獲得無截?cái)嗾`差的數(shù)值結(jié)果。這很大程度上得益于鏈?zhǔn)椒▌t和計(jì)算機(jī)編程模式的相結(jié)合。每一個函數(shù)可以被分解為基本運(yùn)算符的組合，其中基本運(yùn)算符包括加、減、乘、除等二元算術(shù)運(yùn)算符，以及像三角函數(shù)和指數(shù)函數(shù)等在內(nèi)的超越函數(shù)。因此函數(shù)可以采用計(jì)算圖的形式來表示。例如考慮如下的函數(shù)：

(1)

式中，可以視為一個向量函數(shù)f:R2|→R2。其分解得到的計(jì)算圖見0。圖中每個運(yùn)算都由一個帶編號的節(jié)點(diǎn)表示?？梢酝ㄟ^前向和反向兩個模式來計(jì)算f關(guān)于x的導(dǎo)數(shù)。

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

1.1 前向計(jì)算模式

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

(2)

(3)

(4)

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

1.2 反向計(jì)算模式

表1 自動微分反向計(jì)算模式的計(jì)算過程

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

2 并行實(shí)現(xiàn)

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

2.1 主機(jī)端編程

為了簡化使用的方式，對待求導(dǎo)函數(shù)的構(gòu)建應(yīng)該盡量的趨近于原生的C/C++風(fēng)格的代碼編寫。式的函數(shù)可以寫成如下的簡潔形式：

DVAR x1,x2;

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

DVAR f2=x1*x2

要實(shí)現(xiàn)如此的函數(shù)構(gòu)建的簡單化，需要完成一些關(guān)鍵的任務(wù)。首先，定義一個用于描述節(jié)點(diǎn)的結(jié)構(gòu)體：

template

struct ADV_Node {

OpType op;

shared_ptr> arg[2];

T val;

T dval;

int id;

}

在該結(jié)構(gòu)體中，OpType是一個用于指示運(yùn)算符的枚舉類型，指示該節(jié)點(diǎn)為其前驅(qū)進(jìn)行某種運(yùn)算的結(jié)果。除了常規(guī)的運(yùn)算之外，還引入“CONST”和“INPUT”來分別表示該節(jié)點(diǎn)是否為常量或者函數(shù)的輸入?yún)?shù)。每個節(jié)點(diǎn)在生成時(shí)還被賦予一個唯一的id值。arg指向該節(jié)點(diǎn)的一個或兩個前驅(qū)節(jié)點(diǎn)，指示該節(jié)點(diǎn)運(yùn)算符的操作數(shù)。dval用于在求導(dǎo)過程中存放函數(shù)對該節(jié)點(diǎn)的導(dǎo)數(shù)值。

然后，實(shí)現(xiàn)了一個Wrapper類ADV來描述最終的數(shù)學(xué)形式上的變量，并用此來完成數(shù)學(xué)表達(dá)式的構(gòu)建。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變量創(chuàng)建節(jié)點(diǎn)可以使得ADV變量能夠不受作用域的限制，像數(shù)學(xué)形式上的函數(shù)那樣完成代碼中的對應(yīng)函數(shù)。例如要實(shí)現(xiàn)一個3元素向量的點(diǎn)積，可以寫為“dot3(ADV *a, ADV*b)”。該函數(shù)能返回一個新建的ADV變量來記錄所有的運(yùn)行結(jié)構(gòu)。ADV類中所有的構(gòu)造函數(shù)和賦值運(yùn)算符重載都用來確保變量的正確生成。所有的關(guān)于ADV變量間的運(yùn)算都用過類定義外的友元的函數(shù)重載來實(shí)現(xiàn)。每個函數(shù)重載創(chuàng)建新ADV變量來記錄這一關(guān)于其前驅(qū)的運(yùn)算。例如加法的函數(shù)重載實(shí)現(xiàn)為：

ADV operator+(const ADV &x,

const ADV &y)

{

ADV adv;

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

adv()->op=ADD;

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

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

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

return adv;

}

程序上的“c=a+b;”表達(dá)式便可以直接地表示數(shù)學(xué)上的c=a+b。通過使用不同運(yùn)算的函數(shù)重載來構(gòu)造函數(shù)，函數(shù)所對應(yīng)的計(jì)算圖也同時(shí)被構(gòu)建出來。另外對于大多數(shù)情況下使用的雙精度浮點(diǎn)型來說，可以將ADV定義為DVAR。

2.2 計(jì)算序列的生成

在構(gòu)建函數(shù)以及對應(yīng)的計(jì)算圖之后，便可以據(jù)此運(yùn)用反向計(jì)算模式來進(jìn)行導(dǎo)數(shù)的計(jì)算。為了正確地計(jì)算各個節(jié)點(diǎn)的導(dǎo)數(shù)值，同時(shí)也是為了便于設(shè)備端的并行計(jì)算。需要將計(jì)算圖轉(zhuǎn)化為計(jì)算序列。

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

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

因此，采用拓?fù)渑判騺硗瓿煞聪蛴?jì)算序列生成。所生成的反向計(jì)算序列中各節(jié)點(diǎn)滿足相互間的依賴關(guān)系。拓?fù)渑判虻膫未a見算法1：

算法1:用于生成反向計(jì)算序列拓?fù)渑判?/p>

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

輸出：反向計(jì)算序列L

初始化計(jì)數(shù)數(shù)值C, 用于記錄各節(jié)點(diǎn)的入度

forvinVdo

令s,t為v的前驅(qū)節(jié)點(diǎn)

增加s,t在C中的入度

end for

whileL非滿 do

找到具有零入度的節(jié)點(diǎn)v

添加到L尾部

令s,t為v的前驅(qū)節(jié)點(diǎn)

減少s,t在C中的入度

end while

2.3 設(shè)備端并行計(jì)算

當(dāng)在主機(jī)端完成計(jì)算序列的生成之后，便可以將需計(jì)算的函數(shù)的參數(shù)傳入設(shè)備端，在設(shè)備端按照計(jì)算序列計(jì)算出函數(shù)值以及函數(shù)的導(dǎo)數(shù)值。

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

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

算法2:計(jì)算導(dǎo)數(shù)的Kernel函數(shù)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 應(yīng)用實(shí)例

四元數(shù)可以表示為q=。在相機(jī)坐標(biāo)系下對3D點(diǎn)的旋轉(zhuǎn)可以表示為：

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

(5)

其中:⊕為Hamilton積，q(X)則用于將X變?yōu)樗脑獢?shù)形式<0,X>。在Rodrigues參數(shù)形式中，旋轉(zhuǎn)表示為繞單位向量k的旋轉(zhuǎn)，旋轉(zhuǎn)角度為θ。因此，X的旋轉(zhuǎn)在Rodrigues參數(shù)形式下表示為：

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

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

(6)

通過將旋轉(zhuǎn)與位移的結(jié)合，可以得到3D點(diǎn)在相機(jī)上的投影。在光束平差中，視3D點(diǎn)坐標(biāo)X和相機(jī)矩陣[R|t]為投影x的參數(shù)。采用Rodrigues參數(shù)形式，編寫如下的代碼完成x相對于X和[R|t]的計(jì)算圖的構(gòu)建：

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];

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

圖3 Rodrigues參數(shù)形式下投影函數(shù)的計(jì)算圖

4 測試與分析

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

表2 各數(shù)據(jù)集的主要數(shù)據(jù)

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

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

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

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

表3 自動微分在各數(shù)據(jù)集上的測試結(jié)果

5 總結(jié)

本文首先展示了自動微分系統(tǒng)的工作原理，包括正向模式和反向模式。對于多參數(shù)的函數(shù)，揭示了反正模式比正向模式具有更高的效率。以O(shè)penCL為并行計(jì)算框架，實(shí)現(xiàn)了針對大型優(yōu)化問題的并行化自動微分。該實(shí)現(xiàn)采用反向計(jì)算模式，以C/C++的風(fēng)格構(gòu)建函數(shù)，并生成計(jì)算圖和計(jì)算序列。以光束平差為應(yīng)用背景，通過測試，該實(shí)現(xiàn)比Ceres Soler快約3.6倍。同時(shí)可以發(fā)現(xiàn)，自動微分在CPU上的實(shí)現(xiàn)要優(yōu)于在GPU上的實(shí)現(xiàn)。