李俊，周薇娜

　　（上海海事大學 信息工程學院,上海 201306）

　　摘要：為了實現快速運動目標檢測，利用低秩矩陣恢復原理進行視頻前景檢測，主要針對低秩矩陣恢復算法存在的耗費大部分運算時間且運算較為復雜的奇異值分解問題，應用統一計算結構裝置（CUDA）第三方庫實現加速計算奇異值分解的低秩矩陣恢復算法優化，得到快速且高效的前景檢測方法。基于開源視頻序列實驗，與原有的低秩矩陣恢復算法進行各項參數的比較，其中加速倍數達一倍以上。實驗結果證明，經過優化的算法運算時間變短，具有更高效率。

　　關鍵詞：低秩矩陣恢復；前景檢測；統一計算結構裝置

0引言

　　目標檢測是把視頻序列中感興趣的目標分割出來，作為許多視覺應用中最基礎的部分，深入探討其相關算法研究以提高其檢測性能與精度，顯得相當重要。傳統目標檢測算法計算較為復雜，檢測精度相對較低，性能也較低下。隨著計算機技術的高速發展，所拍攝的視頻序列圖像分辨率也越來越高，這使得視頻序列中圖像數據量過大和數據維數過高，增加目標檢測計算難度。對高維數據降維，用降維得到的低維數據表征原有的高維數據從而降低數據計算復雜度成了視頻圖像數據處理的重要方法。主成分分析(Principal Component Analysis,PCA)是一種很有效的分析高維數據的工具，把高維數據投影到低維空間，奇異值分解(Singular Value Decomposition,SVD)使得PCA有了穩定高效的計算能力，但當矩陣元素受到嚴重破壞時，PCA則對矩陣估計不準確。為了改善PCA對這種情況的處理，CANDES E J提出了魯棒主成分分析［1］(Robust Principle Component Analysis,RPCA)，也稱低秩矩陣恢復［2］，指元素受到嚴重破壞的矩陣也得以恢復。

　　一般說來由視頻數據構造成的觀測矩陣可以分解為稀疏和低秩兩個矩陣，矩陣恢復時可以用凸優化的方法，視頻前景檢測擬合這種分解特性。前景檢測是從拍攝的視頻中分離出背景和前景，由于視頻的背景基本是不變的，因此如果把每幀當做矩陣的各列，那么此矩陣將是低秩的，低秩矩陣恢復可以同時完成背景和前景分離，得到的前景即為需要檢測的目標。

　　低秩矩陣恢復算法運行在MATLAB平臺上，MATLAB雖易用性強，但是有著計算慢、效率低等缺點。經過對算法和CUDA［3］深入分析，使用MATLAB和C++ 的MEX混合編程，進行了在MATLAB內調用CUDA第三方庫加速原有算法的運算，經加速優化后，低秩矩陣恢復前景檢測方法效率更高。最后基于真實數據實驗，進行經過優化后的算法與原算法各項參數的比較。

1低秩矩陣恢復

　　觀測得到的數據組成的數據矩陣D秩一般很低，為了恢復它的低秩結構，把矩陣D分解成一個低秩矩陣和一個稀疏矩陣的和，即D=A+E，A未知是低秩的，E也未知是稀疏的。當矩陣E的元素服從獨立同分布的高斯分布時，PCA可解SVD得到最優的矩陣A，即求解下面的最優化問題：

　　min∫‖E‖Fs.trank(A)≤c,D=A+E （1）

　　其中，‖.‖F為F范數，c為子空間維數。當E不理想時，經典PCA對A的估計就相當不準，WRIGHT J 等人提出的RPCA很好地解決了上述問題。

　　同PCA一樣，RPCA的問題可描述為：觀察數據矩陣D為已知，D=A+E，A、E未知，A為低秩的，但是E稀疏且非零元素E可能無限大，要從中恢復出A，需得到觀察數據D的最小的秩，且干擾誤差是稀疏的。通過以下優化問題來描述：

　　minA,Erank(A)+λE0s.tA+E=D（2）

　　其中，rank(A)為A的秩；.0為求解零范數，表示矩陣中不為零元素數目；λ為正常數。因為沒有有效方法解決這個NPHard問題，所以決定對式（2）做松弛變化轉換成容易求解的凸優化問題，將l0范數由l1范數替代，用核范數A=∑iσi(A)代替A的秩，即：

　　minA,EA*+λE1s.tA+E=D（3）

　　該問題也被稱之為主成分追蹤(Principal Component Pursuit,PCP)，在一定條件下，求解式(3)可以得到理想且唯一的A和E。將以上的凸優化問題稱為RPCA，該方程式在實際應用中性能良好，能從高達1/3嚴重腐蝕的觀測值中恢復真實的低秩矩陣，并且具有優秀的理論保證。例如，對一般的矩陣A，它可以顯示成功糾正錯誤恒比，甚至當 A 的秩是近似地正比于環境度:rank(A)=Cm/log(m)。參考文獻［1］也給出了結論定理：如果A0為n×n矩陣，其秩rank(A0)≤ρrnμ－1·(logn)－2，E0是隨機稀疏的n×n矩陣，基數m≤ρsn2ρrρs，其中ρr為低秩率，ρs為稀疏率，當λ=1/n時，PCP能夠準確恢復出L∧=L0,S∧=S0的概率為1－O(n－10)。對于任意矩形矩陣，對λ=1/max(m,n)有與上相同的結論。

2算法

　　前景檢測是把目標由視頻里提取出來。因為視頻的背景基本是不變的，假設視頻有n幀圖像，把每幀當做矩陣D每一列，由每幀圖像的像素構成的m維列向量作為矩陣D每一行，得到大小為m×n低秩矩陣D，因此可用低秩矩陣恢復來恢復出背景。低秩矩陣A可作為視頻中有很大相關性的背景部分，稀疏矩陣E則可作為分布稀疏的前景部分也即感興趣的目標。

　　有許多方法可求解RPCA問題［47］，本文對非精確增廣拉格朗日乘子法(IALM)中的奇異值分解計算進行加速優化，使得能夠加速計算奇異值分解，算法運算效率提升，運行時間縮減。

　　構造增廣拉格朗日函數：

　　L(A,E,Y,μ)=‖A‖+λ‖E‖1,1+

　　<Y,D－A－E>+μ‖D－A－E‖2F/2（4）

　　其中，Y∈Rm×n為約束乘子，μ>0為懲罰參數，·表示計算內積，‖·‖F表示F范數。每一步最小化式(4)時使用輪流更新的方法，先保持Y與E不變，求得使L最小的A,然后保持Y與A不變，求使L最小的E，這樣迭代次數就可收斂到這個子問題的最優解。更新A時：

　　arg minAλ‖A‖*  + μ2‖D-A-E + μ-1Y‖2F

　　=Θμ-1 （D-E + μ-1Y）（5）

　　更新E時：

　　arg minEλ‖E‖1  + μ2‖D-A-E + μ-1Y‖2F

　　=  Λμ-1 （D-A + μ-1Y）（6）

　　輪流更新，直到對子問題求解收斂時為止，該算法被稱為精確增廣拉格朗日乘子法(EALM)，見算法1［4,8］。

　　算法1

　　（1）初始化Y0，E0=0，μ0，k=0

　　（2）while not converged do

　　（3）E0k+1=E*k，j=0

　　（4）while not converged do

　　（5）(U,S,V)=svd(D－Ejk+1+μ－1kY*k)

　　（6）Aj+1k+1=UΛμ－1k［S］VT

　　（7）Ej+1k+1=Λλμ－1k［D－Aj+1k+1+μ－1kYk］

　　（8）j=j+1

　　（9）end while

　　(10)Y*k+1=Yk+μk(D－A*k+1－E*k+1),μk+1=ρμk

　　(11)k=k+1

　　(12)end while

　　稱算法1為精確增廣拉格朗日乘子法，原因在于最外層的循環并不需要求出子問題的精確解，事實上，只要將A和E都更新一次得子問題的一個近似解，這樣得到的最優解就可使算法收斂到原問題，因此可以得到簡潔且收斂更快的算法——IALM，見算法2［4,8］。

　　算法2

　　(1)初始化Y0，E0，μ0，k=0

　　(2)while not converged do

　　(3)(U,S,V)=svd(D－Ek+μ－1kYk)

　　(4)Ak+1=UΛμ－1k［S］VT

　　(5)Ek+1=Λλμ－1k［D－Ak+1+μ－1kYk］

　　(6)Yk+1=Yk+μk(D－Ak+1－Ek+1)

　　(7)更新μk

　　(8)k=k+1

　　(9)end while

　　以上算法會耗費大量時間計算奇異值分解，所以找到快速處理奇異值分解計算的方法，對于整個算法的執行效率將會大有改善。本文基于CUDA 架構第三方CULA［9 10］庫在MATLAB中實現加速優化算法中奇異值分解的計算。MATLAB中CUDA加速的原理為：CPU執行流程控制的操作代碼，一些需要并行運算的操作代碼放在GPU中處理，這樣算法運算時間將大幅度減少。CUDA計算的流程如圖1所示。

　　經過仔細閱讀CULA參考手冊和CULA編程指南，分析低秩矩陣恢復算法，編寫可以在MATLAB里面通過一些接口使用CULA庫加速的culaSvd mex源文件。首先，需要在MATLAB 里面處理多種CULA的數據類型，MATLAB和CULA 支持的4個主要的數據類型為: 單精度實數、雙精度實數、單精度復數和雙精度復數。由于MATLAB MEX編譯器與CULA 支持C++語言，可以通過代碼來實現數據類型的處理，用C++模板來處理所有的4個MATLAB數據類型，使用MATLAB函數mxGetClassID() 和 mxIsComplex()確定輸入矩陣的精度和復雜度。上述代碼主要包括3個部分，SVD函數的CULA頭文件和MATLAB數據類型和接口的頭文件、支持4種MATLAB數據類型的模板函數、在MATLAB中調用的C++模板所需要的網關函數。其次，因為CULA和 MATLAB以不同的類型存儲復雜的數據，創建helper函數，實現兩個類型之間的轉換。 創建好helper函數后，通過查閱CULA的xGESVD函數的文檔可以發現，MATLAB中對角陣返回奇異值，而CULA返回一個向量，MATLAB返回V，而CULA返回V的轉置，為了對接 MATLAB 的接口，必須從奇異值向量轉換成一個對角矩陣并轉置。至此完成了culaSvd程序，在MATLAB中用MEX命令編譯culaSvd.cpp源文件生成二進制mex文件，完成MATLAB中實現在底層調用經CULA 加速優化后的奇異值分解。

3實驗結果及分析

　　本節將原有算法與基于CUDA優化后的算法進行運行效果比較，測試是基于拍攝的真實數據進行的，試驗運行環境配置為：Intel Core i5 M 520 2.40 GHz CPU，4 GB 內存，Nvidia 公司的GT218 GPU，編程環境為MATLAB R2010a+Visual Studio 2010+CUDA Toolkit Version 4.0+CULA R12，選取的視頻是拍攝的學校路口某個時刻的場景。視頻序列總共為200幀，每幀大小為576 768，由此組成的觀測矩陣為M（27 648×200）。

　　圖2為選取的第150幀3種算法運行效果對比圖。表1為取150幀進行前景檢測，兩種算法及其加速改進后算法的運算時間、迭代次數、SVD時間及SVD次數的比較。圖2第180幀3種算法仿真對比圖

　　表2為分別取180幀、160幀、140幀、120幀、100幀和80幀進行前景檢測，為了得到精確運行時間，分別進行100次運算，然后取平均運行時間，計算算法的運行時間、平均加速倍數。對視頻進行處理，將其裁剪為不同分辨率的視頻，對其進行測試，表3 為分別采用不同分辨率進行視頻測試時兩種算法的運行時間及其加速倍數。

　　從仿真效果可以看出，經CULA優化后的算法分離出的背景A3比其他兩種原有的算法得到的背景A1和背景A2相比，前景目標更加分明，前景目標檢測效果更好。

　　表1表明，優化改進后的算法運行效率有所提升。表2取不同的視頻幀數進行前景目標檢測，表明實現了至少1倍以上加速比。最后表3給出了不同測試視頻的運行結果，表明隨著圖像尺寸的加大，加速比增大。綜合得知，本文的優化改進算法效率更高。

4結論

　　本文應用CULA加速優化后的低秩矩陣恢復算法實現前景目標檢測，加速優化后的算法與原有算法相比，奇異值分解次數減少，表明算法里的奇異值分解計算得到加速改進，因此算法的運行時間得以較大幅度縮短。仿真結果表明，本文加速優化后的算法有更好的前景目標檢測效率。本文的前景目標檢測算法在計算效率上得到改善，發揮了軟硬件結合的效果，這種學習成本低廉、使用代價較小、運行效率較高的優點彌補了MATLAB計算速度上的缺點，也使得視覺應用的研究有了良好的開端。

參考文獻

　　［1］ CANDES E J, Li Xiaodong, Ma Yi, et al. Robust principalcomponent analysis［J］. Journal of the ACM,2009,58(3):233279.

　　［2］  WRIGHT J，GANESH A，RAO S，et al．Robust principal component analysis:exact recovery of corrupted low rank matrices via convex optimization［C］.Proceedings of Neural Information Processing Systems (NIPS), 2009,87(4):20:320:56.

　　［3］  NVIDIA Corporation. NVIDIA CUDA programming guide［Z］. 2009.

　　［4］  Lin Zhouchen,Chen Minming,Ma Yi．The augmented Lagrange multiplier method for exact recovery of a corrupted lowrank matrix［R］. Eprint Arxiv,2010.

　　［5］ GANESH A, Lin Zhouchen, WRIGHT J, et al. Fast algorithms for recovering a corrupted low rank matrix［C］.International Workshop on Computational Advances in MultiSensor Adaptive Processing,2009:213216.

　　［6］ Lin Zhouchen, GANESH A,WRIGHT J,et al.Fast convex optimization algorithms for exact recovery of a corrupted low rank matrix［C］.Proceedings of Computational Advances in MultiSensor Adaptive Processing (CAMSAP),2009.

　　［7］ YUAN X, YANG J. Sparse and low rank matrix decomposition via alternating direction methods［R］．Hong Kong: Hong Kong Baptist University，2009．

　　［8］ Chen Minming．Algorithms and implementations of matrix reconstruction ［D］ . Beijing: Graduate School Chinese Academy of Science,2010.

　　［9］ EM Photonics Corporation. CULA reference guide［Z］.  2014.

　　［10］ EM Photonics Corporation. CULA programmers guide［Z］. 2014.