基于GPU的高光譜遙感主成分分析并行優(yōu)化

柳家福 李歡 賀金平 劉天石 王啟聰 吳澤彬,4

（1 南京理工大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院，南京 210094）（2 北京空間機(jī)電研究所，北京 100094）（3 北京航空航天大學(xué)電子信息工程學(xué)院，北京 100191）（4 南京理工大學(xué)連云港研究院，連云港 222006）

0 引言

高光譜遙感具有較高的光譜分辨率，能夠獲取豐富的地表信息，是繼多光譜遙感技術(shù)之后的一種新型遙感技術(shù)。其光譜分辨率小于10nm，波長(zhǎng)范圍在0.37～2.48μm之間，能從紫外線、可見光、近紅外和短波紅外等光譜區(qū)域獲取大量且光譜連續(xù)的遙感數(shù)據(jù)。一幅高光譜圖像通常包含了同一個(gè)地表區(qū)域內(nèi)上百個(gè)不同波段的光譜信息[1]，如我國(guó)首次自主研制的搭載在HJ-1A星上的干涉型高光譜成像儀（HSI），能夠以 115個(gè)波段同時(shí)對(duì)地表物質(zhì)成像[2]。但豐富的遙感數(shù)據(jù)信息，在實(shí)際應(yīng)用中提高了高光譜遙感數(shù)據(jù)傳輸和處理效率的難度。主成分分析（principal component analysis, PCA）是光譜特征提取常用的線性組合方法[3-4]，能夠有效提取出特征光譜，降低光譜圖像的數(shù)據(jù)量級(jí)，但是在實(shí)際的工程應(yīng)用中存在降維計(jì)算速度慢的問題，尤其在處理海量的高光譜數(shù)據(jù)時(shí)十分耗時(shí)，嚴(yán)重影響高光譜圖像處理的實(shí)時(shí)性，計(jì)算效率已成為其應(yīng)用推廣的瓶頸。如何快速降低高光譜圖像的維數(shù)，提取出特征數(shù)據(jù)，是高光譜遙感圖像處理面臨的重要問題。

隨著圖形處理技術(shù)不斷發(fā)展，尤其是NVIDIA公司使用類C語言開發(fā)的用于通用計(jì)算的CUDA框架問世之后[5]，具備CUDA框架的顯卡可以支持大量的線程并行，而且能夠動(dòng)態(tài)的調(diào)度和執(zhí)行。該框架中的單指令多線程執(zhí)行模型，是對(duì)單指令多數(shù)據(jù)的一種改進(jìn)，這種模型能夠?qū)?zhí)行數(shù)據(jù)的寬度作為硬件細(xì)節(jié)被隱藏起來，使硬件可以自適應(yīng)不同的執(zhí)行寬度，提高編程的靈活性。此外，圖形處理器（graphic processing unit，GPU）的設(shè)計(jì)者將更多的晶體管作為執(zhí)行單元，而不是像中央處理器（central processing unit，CPU）那樣用作復(fù)雜的控制單元和高速緩存，因而在密集型的數(shù)據(jù)計(jì)算中GPU具有明顯的優(yōu)勢(shì)，CPU則更擅長(zhǎng)邏輯控制較多的數(shù)據(jù)運(yùn)算。以上特點(diǎn)，使圖像處理單元具有高度的并行計(jì)算和并行數(shù)據(jù)處理能力，成為能夠輔助CPU計(jì)算的通用計(jì)算單元。有學(xué)者嘗試?yán)肎PU來提高PCA算法的執(zhí)行效率，并取得了一定的成果。文獻(xiàn)[6]比較了基于SSE與GPU兩種方法對(duì)PCA算法的加速效果，結(jié)果證明基于GPU的PCA加速算法的最高加速比達(dá)到20倍，明顯優(yōu)于基于SSE（stream ing SIMD extensions）的加速算法，但效率仍有較大提高；文獻(xiàn)[7]研究了光譜解混的GPU加速流程，其中PCA算法加速比最高可達(dá)122倍，但在特征分解部分采用乘冪法，該方法主要用于求取實(shí)矩陣的主特征值，求矩陣的全部特征值需要對(duì)矩陣進(jìn)行降階，而每降階一次計(jì)算精度就會(huì)損失或降低一些，因此這種PCA算法的計(jì)算精度也較低。由于PCA算法的復(fù)雜性，部分計(jì)算過程并不適合GPU并行執(zhí)行，因而如何設(shè)計(jì)適當(dāng)?shù)牟⑿兴惴鞒?，在保證PCA并行算法的精度基礎(chǔ)上提高計(jì)算效率，已成為快速遙感處理領(lǐng)域中亟待解決的問題。

針對(duì)PCA串行計(jì)算方法效率低、GPU并行算法精度難以保證的特點(diǎn)，本文提出基于GPU+CPU異構(gòu)系統(tǒng)的PCA并行優(yōu)化算法，將適合GPU加速計(jì)算的協(xié)方差矩陣計(jì)算與矩陣投影運(yùn)算部分在GPU上執(zhí)行，而特征分解部分采用改進(jìn)的Jacobi單側(cè)旋轉(zhuǎn)算法在CPU上進(jìn)行計(jì)算[8]，充分利用CPU與GPU各自的計(jì)算特性，在保證了PCA算法的精度的同時(shí)提高了該算法的計(jì)算效率。

1 基于PCA方法的高光譜遙感數(shù)據(jù)特征提取原理

高光譜遙感數(shù)據(jù)特征提取是將光譜數(shù)據(jù)特征按一定準(zhǔn)則由原高維空間變換到較低維數(shù)空間，提取的特征應(yīng)盡可能保留針對(duì)地物光譜差異的有價(jià)值信息[9]。一般通過對(duì)原始n維特征空間進(jìn)行空間變換，然后求得其特征子空間。特征提取需要找到一個(gè)映射關(guān)系Q:X→Y，將原始特征空間的n個(gè)特征X=(X,X,X,…,X)T映射到維數(shù)較低的特征子空間Y中，使其通過這種映射關(guān)系產(chǎn)生m個(gè)新特征

1 2 3n

Y= (Y,Y,Y,…,Y)T，其中m＜n。

1 2 3m

PCA又稱為 K-L變換，是高光譜特征提取常用的線性組合方法[1]。若X=(X1,X2,X3,…,Xn)T是n維隨機(jī)變量，且均值E(X)=μ，協(xié)方差D(X)=Δ，考慮其線性變換:

可得:

式中ai（i= 1,2,3,…,n）為X變換的系數(shù)矩陣；Ci（i= 1,2,3,…,n）為變換矩陣；Var(Ci)為Ci的方差；Cov(Ci,Cj)為Ci、Cj的協(xié)方差。

假設(shè)C1=a1TX為所求的第一主成分，則C1應(yīng)能盡可能多地反映原有n個(gè)變量的信息，其信息可以用C1的方差來表達(dá)。該問題轉(zhuǎn)化為求a1= [a11,a21,…,an1]T，使得在約束條件a1=1下，Var(Ci)達(dá)到最大值。可用拉格朗日乘數(shù)法求解，令

得到:

式中λ為拉格朗日乘子；為對(duì)φ關(guān)于自變量a1求導(dǎo)；為對(duì)φ關(guān)于自變量λ求導(dǎo)；I為單位矩陣。由于a1≠0，因此∑-λI=0，即轉(zhuǎn)化為求解協(xié)方差矩陣Σ特征向量和特征值問題。假設(shè)λ=λ1是∑的最大特征值，則相應(yīng)的單位特征向量a1即為所求值。以此類推，求X的第i個(gè)主成分需要先求出協(xié)方差矩陣∑的第i個(gè)特征值λi對(duì)應(yīng)的單位特征向量ai，為了更有效地表征原來變量的信息，前i-1個(gè)主成分體現(xiàn)的信息不希望在Ci中出現(xiàn)，則需添加約束條件j= 1,2,3,…,i-1；通過i= 1,2,3,…,n進(jìn)行運(yùn)算，即可得到變換后的矩陣。傳統(tǒng)的PCA算法流程如圖1所示。

圖1 傳統(tǒng)PCA算法流程Fig. 1 Traditional PCA algorithm flow

計(jì)算時(shí)通常只需要選擇與最大的m個(gè)特征值相對(duì)應(yīng)的特征向量即可，這樣便可以減少原始矩陣的維數(shù)，還保留了主要的數(shù)據(jù)信息。

2 基于GPU+CPU異構(gòu)系統(tǒng)的主成份分析方法并行優(yōu)化

2.1 PCA并行優(yōu)化設(shè)計(jì)

根據(jù)上述分析，PCA算法的關(guān)鍵步驟在于協(xié)方差矩陣的計(jì)算、特征分解以及光譜數(shù)據(jù)投影。其中協(xié)方差矩陣計(jì)算最為耗時(shí)，是PCA算法并行優(yōu)化的重點(diǎn)[10-11]；特征分解則是該算法的性能瓶頸，對(duì)精度有很大影響；光譜矩陣投影是常規(guī)的矩陣運(yùn)算，十分適合GPU并行加速。

高光譜遙感數(shù)據(jù)的協(xié)方差矩陣計(jì)算，用傳統(tǒng)方法需要對(duì)各波段圖像求均值并對(duì)所有像元逐一去中心化，但由于原始高光譜矩陣的數(shù)據(jù)量大，該方法的執(zhí)行效率較低。本文對(duì)其進(jìn)行改進(jìn)，如式（6）所示，若)是高光譜圖像任意兩個(gè)波段對(duì)應(yīng)的像元向量，∑ij為其協(xié)方差，s為高光譜圖像像元個(gè)數(shù)。使用傳統(tǒng)方法需要對(duì)N i、N j中每個(gè)像元分別求均值，由于s的量級(jí)通常大于 105，因而計(jì)算量很大。改進(jìn)后的方法只需對(duì)結(jié)果去一次均值即可，即先求N i、N j對(duì)應(yīng)像元乘積的累加和，再求出N i、N j對(duì)應(yīng)的均值，二者相減即可求出波段i和j的協(xié)方差∑ij，從而極大地提高運(yùn)行效率。

對(duì)于協(xié)方差矩陣的特征分解，求解算法很多，如 QR算法[12]、乘冪法[7]、MRRR算法[13]、Arnoldi算法[14]、Jacobi-Davidson[15]算法等，但這些算法都有各自的特性，或是通用算法，收斂速度慢；或是精度不高，準(zhǔn)確性低。如文獻(xiàn)[7]中提出的PCA并行優(yōu)化算法（以下稱為SPCA算法），利用乘冪法進(jìn)行特征分解，該方法采用矩陣乘法運(yùn)算，適合并行加速，但是由于乘冪法求矩陣的全部特征值需要對(duì)矩陣進(jìn)行降階運(yùn)算，導(dǎo)致相鄰特征對(duì)數(shù)據(jù)依賴較大，每降階一次計(jì)算精度就會(huì)降低一些，所以降階法實(shí)際上只可使用少數(shù)幾次，僅適用于求矩陣前幾個(gè)按模最大特征值及相應(yīng)特征向量。本文采用了Jacobi單側(cè)旋轉(zhuǎn)算法求解協(xié)方差矩陣的特征分解值，該算法是對(duì)傳統(tǒng)Jacobi雙側(cè)旋轉(zhuǎn)算法的改進(jìn)。Jacobi算法的主要思想是通過正交相似變換將一個(gè)實(shí)對(duì)稱矩陣對(duì)角化，從而求出該矩陣的全部特征值和對(duì)應(yīng)的特征向量。已知協(xié)方差矩陣Σ是對(duì)稱方陣，設(shè)λ為其特征值，a為其特征向量矩陣，則:Σa=λa。由于∑=ΣT，所以ΣTΣa=Σλa=λλa，這說明λ2為ΣTΣ的特征值，a為ΣTΣ的特征向量，即ΣTΣ與Σ具有相同的特征向量a，且特征值具有平方關(guān)系?，F(xiàn)使用Givens旋轉(zhuǎn)變換對(duì)Σ進(jìn)行一系列的列變換，得到方陣Q，使其各列兩兩正交，即ΣV=Q，這里V為正交化過程的變換方陣。由于QTQ為n階對(duì)稱方陣，因此可得即可得到特征值λ，其中V即為所求特征向量。由以上分析可知，Jacobi單側(cè)旋轉(zhuǎn)法僅需對(duì)對(duì)稱方陣實(shí)施列變換，就可使數(shù)據(jù)相關(guān)關(guān)系僅在同列之間，不但計(jì)算速度快，而且精確性較高。

根據(jù)上述方法，將求得的特征向量矩陣a= [a1,a2,a3,…,am]T與光譜矩陣N進(jìn)行矩陣運(yùn)算，即可得到經(jīng)過特征提取的降維矩陣。

本文設(shè)計(jì)的基于GPU+CPU異構(gòu)系統(tǒng)的高光譜遙感主成分分析并行優(yōu)化算法如圖2所示，其中s為高光譜圖像像元個(gè)數(shù)；M為高光譜圖像及其自身轉(zhuǎn)置乘積；M uv為第u、v波段對(duì)應(yīng)像元向量?jī)?nèi)積；Su、Sv分別為第u、v波段對(duì)應(yīng)像元向量累加和。當(dāng)高光譜數(shù)據(jù)加載到內(nèi)存后，CPU端將數(shù)據(jù)由主機(jī)端拷貝到GPU的Global Memory；由設(shè)備端執(zhí)行Kernel函數(shù)，計(jì)算高光譜數(shù)據(jù)的矩陣乘積以及高光譜圖像像元向量的累加和；再由改進(jìn)的協(xié)方差矩陣公式計(jì)算出高光譜圖像的協(xié)方差矩陣；由CPU端輔助計(jì)算出協(xié)方差矩陣的特征向量矩陣，并按特征值的大小進(jìn)行排序；最后再由GPU計(jì)算出特征圖像，并將結(jié)果拷回主機(jī)端內(nèi)存。

圖2 PCA方法的并行優(yōu)化算法流程Fig. 2 Parallel optim ization flow of PCA

2.2 協(xié)方差矩陣計(jì)算并行優(yōu)化

計(jì)算高光譜數(shù)據(jù)的協(xié)方差矩陣Σ，需要對(duì)高光譜數(shù)據(jù)及其轉(zhuǎn)置求矩陣乘積，該過程十分適合 GPU并行加速，可采用CUBLAS庫中的矩陣乘積函數(shù)進(jìn)行計(jì)算，此函數(shù)會(huì)發(fā)布與結(jié)果矩陣相同大小的線程個(gè)數(shù)來執(zhí)行計(jì)算，效率較高。為了計(jì)算矩陣N的協(xié)方差，還需要知道每個(gè)波段圖像的均值，由于各波段求和相互獨(dú)立，可以使用GPU進(jìn)行計(jì)算，而各波段內(nèi)部求和可采用Reduction的方法來進(jìn)行計(jì)算，為此需要發(fā)布一個(gè)名為SumKernel的核來執(zhí)行該操作，其中SumKernel為核函數(shù)名稱，本文實(shí)驗(yàn)中，該核創(chuàng)建與原始光譜圖像波段數(shù)n相同個(gè)數(shù)的線程塊，每個(gè)塊內(nèi)包含256個(gè)線程，計(jì)算一個(gè)波段圖像累加和的線程，在一個(gè)塊內(nèi)，塊內(nèi)線程共享一個(gè)共享存儲(chǔ)器，計(jì)算時(shí)使用合并訪問將數(shù)據(jù)拷貝到該共享內(nèi)存內(nèi)可提高讀取數(shù)據(jù)的速度。如圖3所示，描述了該累加和的縮減過程。計(jì)算結(jié)束后，根據(jù)式（6），將光譜數(shù)據(jù)乘積減去該均值，即可求出光譜矩陣的協(xié)方差，該過程如圖4所示。

圖3 GPU累加和歸約過程Fig. 3 Reduction process of summation on GPU

圖4 協(xié)方差矩陣計(jì)算并行優(yōu)化流程Fig.4 Parallel optim ization flow of covariance matrix calculation

2.3 特征分解計(jì)算及矩陣投影并行優(yōu)化

通過以上方法計(jì)算得到圖像各波段的協(xié)方差矩陣Σ后，需要計(jì)算該矩陣的特征分解值，記特征向量為a= [a,a,a,…,a]T，由以上分析可知，乘冪法適合GPU并行加速，但每降階一次，計(jì)算精度就會(huì)損

1 2 3m失或降低一些，實(shí)際中只可用少數(shù)幾次。本文采用Jacobi單側(cè)旋轉(zhuǎn)算法進(jìn)行求解，該算法是Jacobi雙側(cè)旋轉(zhuǎn)法的改進(jìn)，精確性較高，盡管數(shù)據(jù)之間的依賴性較大，不能在GPU上加速，但這已不是并行加速的瓶頸。得到特征向量矩陣后，按方差大小進(jìn)行排序，取前m個(gè)特征向量a，將其與矩陣N相乘，即得到降維矩陣C。這一步仍然是標(biāo)準(zhǔn)的矩陣相乘運(yùn)算，可以調(diào)用GPU中CUBLUAS庫的內(nèi)置函數(shù)cublasSgemm來執(zhí)行。

與SPCA方法相比，本文提出的優(yōu)化算法有兩方面改進(jìn):1）在波段對(duì)應(yīng)圖像累加和計(jì)算方面，采用線程折半順序累加，減少了bank沖突，提高了訪存效率；2）在計(jì)算協(xié)方差矩陣特征對(duì)時(shí)，使用Jacobi單側(cè)旋轉(zhuǎn)法，不但精度較高，而且速度也接近SPCA方法中的乘冪法。

3 實(shí)驗(yàn)結(jié)果與分析

根據(jù)以上算法，本文采用NASA的AVIRIS高光譜圖像數(shù)據(jù)（美國(guó)內(nèi)華達(dá)州南部沙漠的Cuprite數(shù)據(jù)）進(jìn)行模擬測(cè)試，圖像已經(jīng)經(jīng)過大氣校正并且已轉(zhuǎn)換為光譜反射率，原有224個(gè)波段，去除1～3、105～115、150～170、223～224噪聲段，剩余波段數(shù)為187；然后從中截取3個(gè)不同像元大小的圖像（300×300像元、450×450像元、614×512像元）進(jìn)行試驗(yàn)測(cè)試，其中614×512像元是在NASA網(wǎng)站上下載的原始高光譜圖像。本文實(shí)驗(yàn)環(huán)境為:英特爾Xeon E5603四核CPU；主頻為1.6GHz；內(nèi)存為8 Gbyte。使用的GPU為Nvidia Quadro 600，擁有96個(gè)CUDA cores和1Gbyte的顯存。軟件環(huán)境為:Windows7 專業(yè)版64bit操作系統(tǒng)；Visual Studio 2010集成開發(fā)環(huán)境和CUDA 5.0版本工具包。

3種不同數(shù)據(jù)規(guī)模下，串行算法和并行算法的運(yùn)行時(shí)間以及并行算法的加速比見表1。結(jié)果表明，隨著數(shù)據(jù)規(guī)模的增大，加速效果更加明顯，當(dāng)數(shù)據(jù)大小為614×512像元時(shí)，最高加速比可達(dá)141倍，其并行算法優(yōu)勢(shì)明顯。

表1 不同數(shù)據(jù)的執(zhí)行時(shí)間及加速比Tab. 1 Execution time and speedup of different hyperspectral images

表2比較了本文方法（JPCA算法）、SPCA以及MATLAB計(jì)算特征分解中特征值結(jié)果，由表中的數(shù)據(jù)可知，采用本文的特征值求解算法，其特征值結(jié)果和 MATLAB軟件計(jì)算的結(jié)果完全相同，而 SPCA方法盡管適合GPU并行優(yōu)化，但其計(jì)算精度卻有較大的誤差。

表2 PCA特征值比較Tab. 2 Eigenvalues comparison of PCA

綜合以上實(shí)驗(yàn)數(shù)據(jù)可以驗(yàn)證，本文提出的基于GPU+CPU異構(gòu)系統(tǒng)的高光譜遙感PCA并行優(yōu)化算法，在保證算法精度的同時(shí)，計(jì)算效率也得到了顯著提高。

4 結(jié)束語

本文從高光譜遙感圖像數(shù)據(jù)處理的實(shí)時(shí)性問題出發(fā)，提出了基于GPU+CPU異構(gòu)系統(tǒng)對(duì)高光譜遙感圖像進(jìn)行特征提取的PCA并行優(yōu)化方法，經(jīng)實(shí)驗(yàn)驗(yàn)證，達(dá)到以下結(jié)果:

1）改進(jìn)了協(xié)方差矩陣的計(jì)算步驟，優(yōu)化了高光譜像元去均值的計(jì)算流程，減少了運(yùn)算次數(shù)，避免了多次訪存操作；

2）優(yōu)化了GPU中的累加和計(jì)算，減少了非合并訪問；

3）在特征分解部分，使用精確度更高的Jacobi單側(cè)旋轉(zhuǎn)算法進(jìn)行求解，并將該部分設(shè)計(jì)在CPU端執(zhí)行，不但提高其計(jì)算效率，而且精確度很高；

4）充分利用GPU的存儲(chǔ)層次模型，將部分?jǐn)?shù)據(jù)存儲(chǔ)在共享內(nèi)存，減少對(duì)全局內(nèi)存的訪問，節(jié)省了訪存時(shí)間。

GPU的并行優(yōu)勢(shì)主要體現(xiàn)在密集型數(shù)據(jù)計(jì)算上，將迭代次數(shù)多、需要反復(fù)通信的特征分解放在CPU端執(zhí)行，也體現(xiàn)了基GPU+CPU異構(gòu)計(jì)算的模型思想。本文提出的PCA并行優(yōu)化方法達(dá)到了高光譜數(shù)據(jù)實(shí)時(shí)特征提取的效果。但在實(shí)際應(yīng)用過程中還需考慮GPU顯存的限制，當(dāng)高光譜遙感數(shù)據(jù)超出GPU的顯存能力時(shí)，如何對(duì)數(shù)據(jù)進(jìn)行分塊處理，平衡數(shù)據(jù)塊在存儲(chǔ)設(shè)備間傳輸?shù)臅r(shí)間開銷，需要進(jìn)一步的研究。

（References）

[1] 童慶禧, 張兵, 鄭蘭芳. 高光譜遙感—原理、技術(shù)與應(yīng)用[M]. 北京: 高等教育出版社, 2006.TONG Qingxi, ZHANG Bing, ZHENG Lanfang. Hyperspectral Remote Sensing-Principles, Technology and Applications[M]. Beijing: Higher Education Press, 2006. (in Chinese)

[2] 王愛春, 閔祥軍, 李杏朝, 等.“環(huán)境-1號(hào)”A星高光譜成像儀飛行定標(biāo)[J]. 航天返回與遙感, 2009, 30(3): 34-41.WANG Aichun, M IN Xiangjun, LI Xingchao, et al. In-flight Absolute Calibration of the HJ-1A HSI[J]. Spacecraft Recovery& Remote Sensing, 2009, 30(3): 34-41. (in Chinese)

[3] Kambhatla N, Leen T K. Dimension Reduction by Local Principal Component Analysis[J]. Neural Computation, 1997, 9(7):1493-1516.

[4] Manolakis D G, Marden D B. Dimensionality Reduction of Hyperspectral Imaging Data Using Local Principal Components Transforms[C]//Defense and Security. International Society for Optics and Photonics, 2004: 393-401.

[5] Manavski S A, Valle G. CUDA Compatible GPU Cards as Efficient Hardware Accelerators for Smith-Waterman Sequence Alignment[J]. BMCbioinformatics, 2008, 9(2): 10-19.

[6] Joth R, Antikainen J, Havel J, et al. Real-time PCA Calculation for Spectral Imaging (Using SIMD and GP-GPU)[J]. Journal of Real-time Image Processing, 2012, 7(2): 95-103.

[7] Sánchez S, Ramalho R, Sousa L, et al. Real-time Implementation of Remotely Sensed Hyperspectral Image Unm ixing on GPUs[J]. Journal of Real-Time Image Processing, 2012: 1-15.

[8] 李平. 基于FPGA的矩陣特征值并行計(jì)算研究[D]. 重慶: 重慶大學(xué), 2013: 44-47.LI Ping. Study on FPGA-based Parallel Computing of the Matrix Eigenvalues[D]. Chongqing: Chongqing University, 2013:44-47. (in Chinese)

[9] 何明一, 暢文娟, 梅少輝. 高光譜遙感數(shù)據(jù)特征挖掘技術(shù)研究進(jìn)展[J]. 航天返回與遙感, 2013, 34(1): 1-12.HE M ingyi, CHANG Wenjuan, MEI Shaohui. Advance in Feature M ining from Hyperspectral Remote Sensing Data[J].Spacecraft Recovery & Remote Sensing, 2013, 34(1): 1-12. (in Chinese)

[10] LU J, ZHANG B, HUANG W, et al. IHS Transform Algorithm of Remote Sensing Image Data Fusion Based on GPU[J].Computer Engineering, 2009, 35(7): 261-263.

[11] SONG J, ZHOU S. Fast Image Matching Based on GPU Parallel Computing[J]. Journal of Hubei University for Nationalities:Natural Science Edition, 2011, 29(3): 306-310.

[12] 袁暉坪. 關(guān)于酉對(duì)稱矩陣的QR分解及其算法[J]. 系統(tǒng)科學(xué)與數(shù)學(xué), 2012, 32(2): 172-180.YUAN Huiping. On QR Factorization and Algorithm for Unitary Symmetric Matrix[J]. Journal of Systems Science and Mathematical Sciences, 2012, 32(2): 172-180. (in Chinese)

[13] Petschow M, Bientinesi P. The Algorithm of Multiple Relatively Robust Representations for Multi-core Processors[M].Berlin: Springer Berlin Heidelberg, 2012: 152-161.

[14] Guan L, Gao J L, Wang Z W, et al. A Refined Arnoldi Algorithm Based Krylov Subspace Technique for MEMS Model Order Reduction[J]. Key Engineering Materials, 2012, 503: 260-265.

[15] Andrzejew ski J. On Optimizing Jacobi–Davidson Method for Calculating Eigenvalues in Low Dimensional Structures Using Eight Band k/p Model[J]. Journal of Computational Physics, 2013, 249: 22-35.