大規(guī)模非對稱線性方程組Lanczos 算法和精化Lanczos算法的對比

張亞蕾，楊少靜

1.仰恩大學(xué)數(shù)學(xué)系，福建 泉州 362014；2.河北科技學(xué)院公共課部，河北 保定 071000

0 引言

現(xiàn)代社會(huì)當(dāng)中許多科學(xué)和工程技術(shù)等實(shí)際問題都可以轉(zhuǎn)換歸納成為非對稱線性方程組的特征值和特征向量求解問題，也被稱為非對稱性方程組求解問題，例如醫(yī)學(xué)上的數(shù)據(jù)分析問題、圖像以及信號等問題的處理、化學(xué)反應(yīng)問題等，最后都?xì)w結(jié)為解決矩陣的特征問題，通過求解大規(guī)模的非對稱線性方程組

Rm×n(Cm×n)αi=λiαi

(1)

來解決實(shí)際問題.其中，Rn(Cn) 表示n維實(shí)(復(fù))向量空間，Rm×n(Cm×n)表示m×n實(shí)(復(fù))矩陣的全體，(λi,αi)i=1,2,...,n為Rm×n(Cm×n)的特征對( ‖αi‖=1).因此研究非對稱線性方程組的特征求解問題的數(shù)值方法具有非常重要的作用，關(guān)于探討算法設(shè)計(jì)的好壞，開發(fā)出更有效的算法也成為一個(gè)非常重要的問題.

為了求出更優(yōu)，精度更高的特征對，需要對算法的精度和實(shí)用性進(jìn)行反復(fù)地研究.一個(gè)好的算法不僅需要具備數(shù)值收斂、占用的存儲空間小、運(yùn)算量少，還需要易于計(jì)算機(jī)實(shí)現(xiàn).傳統(tǒng)方式解決此類問題應(yīng)用Lanczos 算法，當(dāng)原始矩陣規(guī)模龐大時(shí)，利用這種方法就存在著計(jì)算復(fù)雜，存儲量過大，實(shí)際問題中我們產(chǎn)生的矩陣都是上萬階，甚至于幾十萬階，因此它的速度和內(nèi)存問題就會(huì)格外突出.Lanczos 算法在求解時(shí)采用的是正交投影方法[1]，雖然只需要存儲4個(gè)向量，但是其特征值的精度不夠高.除此之外，隨著算法生成的Krylov子空間維度的增大，向量會(huì)逐漸失去正交性，直接導(dǎo)致算法的收斂速度和收斂性降低.為了解決這一系列的問題，對原始的Lanczos 算法進(jìn)行改進(jìn)和精化，使得求解出的特征對精度更高，速度更快.改變算法的解題方式，從正交投影的方式轉(zhuǎn)變成為雙正交的方式，由于矩陣空間維度擴(kuò)大會(huì)影響向量的雙正交性，所以借助壓縮技術(shù)，在得到高精度特征值的同時(shí)維持其雙正交性，且能減少硬件的內(nèi)存消耗[2].

1 精細(xì)化對比

1.1 構(gòu)造三角矩陣

首先對非線性方程組進(jìn)行初始化處理，將方程組轉(zhuǎn)換成等價(jià)的三角矩陣.非線性方程組的矩陣形式為

(2)

其中，大寫字母V,W,Y表示矩陣，小寫字母v,w表示向量.在三角化中，首先設(shè)定u為啟動(dòng)向量，對于i=1,2,...,n,計(jì)算式(3)：

Wk+1(β1e1)=bAWk=Wk+1Yk+1(∶,k+1)

(3)

其中，β1是非負(fù)的且滿足Wk+1∈Rk+1,k+1，Wk+1為非對稱三角矩陣，Yk+1(∶k+1)表示的是將Wk+1的最后一列去除掉.通常情況下，要是想找到矩陣A的奇異三元組，則必須對Wk進(jìn)行奇異值分解，接著用Wk的奇異值近似A的奇異值，將Wk的奇異值向量與Lanczos向量相結(jié)合并對A的奇異值向量進(jìn)行逼近，則可進(jìn)行三角化處理，設(shè)A為非對稱線性方程的等價(jià)矩陣表現(xiàn)形式，且A∈Rm×n(Cm×n)，b表示正向量且b∈Rm×n，k為整數(shù).執(zhí)行循環(huán)算法：

y=norm(b);v=b/y;g=A*v;x=v*g;g=g-x*v;T(1,1)=x;
fori=2:k;Y(i,i-1)=Y(i-1,i)=y;end

1.2 精化向量子空間的Ritz值

(4)

由多項(xiàng)式的連續(xù)性可以得出，位移點(diǎn)對于不需要奇異值的近似程度更好.用每一次計(jì)算出的那些不需要的Ritz值作為位移，即為準(zhǔn)確位移.取那些不需要的近似奇異值(Ritz值)作為位移，由此得到重新啟動(dòng)方法.

1.3 精細(xì)化投影

(5)

(6)

通過算法計(jì)算出近似值與真實(shí)值進(jìn)行對比，對比兩種算法計(jì)算出的近似值與真實(shí)值之間的誤差，從而得出精化度對比結(jié)果.誤差估計(jì)以子空間迭代法為標(biāo)準(zhǔn)進(jìn)行對比，計(jì)算是在IBM-PC機(jī)上進(jìn)行的.兩種算法都是采用的多個(gè)初始向量進(jìn)行迭代計(jì)算，精化Lanczos算法對所有初始向量采用了平行迭代格式，且借助了逐個(gè)加入初始向量的迭代格式.從而在取相同截?cái)嘀档那闆r下，精化Lanczos算法的初始向量實(shí)際參加迭代的次數(shù)更少.迭代次數(shù)越多則收斂越快，計(jì)算結(jié)果精度也就越高.針對誤差值的對比分析進(jìn)行數(shù)據(jù)統(tǒng)計(jì)并繪制對比曲線如圖1所示.

圖1 精細(xì)度對比曲線Fig.1 The curve of precision contrastion

結(jié)合圖中曲線分析，受收斂作用的影響，兩種算法的精度隨著求解數(shù)量的增加而呈上升趨勢.數(shù)量較少時(shí)精化Lanczos算法的精度較低，比Lanczos算法的精度低102，待數(shù)量增多后精化Lanczos算法的誤差急劇上升，最后維持在102以上.

2 收斂效果對比

表1 Lanczos算法與精化Lanczos算法的收斂效果Tab.1 The convergence effect of Lanczos algorithm and refined Lanczos algorithm

從表1中可知，兩種算法在運(yùn)行過程中得到的收斂結(jié)果相似，但由于在算法運(yùn)行當(dāng)中的收斂效果主要體現(xiàn)在算法的運(yùn)行速度上，由于精化算法自身具備收斂檢驗(yàn)和重啟操作，所以精化Lanczos算法計(jì)算相對更加簡單快速.

3 時(shí)間消耗對比

時(shí)間消耗是檢測一個(gè)算法質(zhì)量和效率的重要參考依據(jù)，在保證準(zhǔn)確率的情況下，算法計(jì)算的運(yùn)行時(shí)間越短，該算法的質(zhì)量也就越高，性能也就越強(qiáng).由于精化Lanczos算法有收斂檢驗(yàn)的功能，所以可以自身控制收斂速度，進(jìn)而影響了計(jì)算的時(shí)間.針對Lanczos算法和精化Lanczos算法都可以順利地進(jìn)行大規(guī)模非對稱線性方程組的求解，并得到一定誤差范圍內(nèi)的準(zhǔn)確結(jié)果，可以確保此兩種算法為可運(yùn)行的算法，即可進(jìn)行計(jì)算時(shí)間的探究與分析，可以統(tǒng)計(jì)出兩種算法計(jì)算方程組消耗的時(shí)間并整理成對比圖，如圖2所示.

從圖2中結(jié)果來看，對于大規(guī)模非對稱線性方程組的求解算例，兩種Lanczos算法計(jì)算得到的參數(shù)估值均值出入不大，而從觀測數(shù)據(jù)上可以看出，精化Lanczos算法的計(jì)算時(shí)間在每次計(jì)算時(shí)均小于Lanczos算法30秒的時(shí)間.從圖中可以看出隨著數(shù)據(jù)量的增加，利用Lanczos算法計(jì)算時(shí)所涉及系數(shù)矩陣的協(xié)因數(shù)陣的維度擴(kuò)大，分解出的向量構(gòu)造也逐漸復(fù)雜，而精化Lanczos算法提取了系數(shù)矩陣中重復(fù)出現(xiàn)的隨機(jī)元素，進(jìn)行收斂操作和檢驗(yàn)，此時(shí)兩種Lanczos算法的計(jì)算時(shí)間差會(huì)迅速增大，精化的Lanczos算法大大節(jié)省了時(shí)間消耗.因此，實(shí)驗(yàn)同樣證明了對大量數(shù)據(jù)進(jìn)行求解處理時(shí)，采用精化Lanczos算法可以明顯地提高計(jì)算效率，減少計(jì)算時(shí)間.

圖2 兩種算法時(shí)間消耗對比圖Fig.2 The comparison of time consumption between two algorithms

4 占用內(nèi)存空間對比

在精化Lanczos算法當(dāng)中引用了隱式壓縮重啟技術(shù)[10]，其操作分為兩個(gè)步驟，一個(gè)是隱式壓縮，該技術(shù)的特點(diǎn)是選擇若干個(gè)位移，對Lanczos過程產(chǎn)生的低階上矩陣進(jìn)行分解，使得計(jì)算量和存儲量大大減少了，從而加快了Lanczos算法的收斂速度.隱式壓縮重啟技術(shù)已廣泛應(yīng)用于特征值問題，例如隱式重啟Lanczos算法、Lanczos壓縮重啟算法、循環(huán)收縮Lanczos算法、隱式重啟塊Lanczos算法、隱式重啟全局Lanczos算法[10]等.該技術(shù)在全局的應(yīng)用中首先令k為要求特征值的個(gè)數(shù)，i=k+p，考慮i=k+p步全局Lanczos算法過程，有

(7)

其中,e為單位向量，Is為最小奇異值的奇異向量.x為位移，對Wk+p做帶位移的QR分解，有

Wk+p-xI=QR

(8)

其中，Q是正交矩陣，R是三角矩陣.經(jīng)過p步帶位移隱式重啟之后，可以得出結(jié)果：

(9)

5 計(jì)算結(jié)果對比

計(jì)算結(jié)果分為計(jì)算結(jié)果的準(zhǔn)確性和殘量值[5].設(shè)定CPU表示計(jì)算時(shí)間單位為秒，Iter表示重啟次數(shù)，m為投影了空間的維數(shù)，允許誤差界tol=10-5，全局算法的殘量范數(shù)取res，初始矩陣或向量均為隨機(jī)選取.調(diào)用Matlab R2010a中的eig命令求得矩陣A的4個(gè)最大特征值為

?1=7.992 413 314 948 18
?2=7.981 047 676 817 97
?3=7.981 047 676 817 96
?4=7.699 682 038 687 75

用精化Lanczos算法計(jì)算初始矩陣的4個(gè)最大特征值，使用兩種解題算法進(jìn)行計(jì)算得出計(jì)算結(jié)果如表2所示.

表2 兩種算法計(jì)算結(jié)果Tab.2 The results of two methods

從計(jì)算結(jié)果中可以看出 Lanczos算法得出的結(jié)果為小數(shù)點(diǎn)后兩位，而精化 Lanczos算法特征值為小數(shù)點(diǎn)后五位，與設(shè)定的最大特征值對比兩種算法的計(jì)算結(jié)果均在準(zhǔn)確范圍之內(nèi)，但精化的算法在結(jié)果的體現(xiàn)上更加準(zhǔn)確，即精化 Lanczos算法的準(zhǔn)確性更高.至于殘量值根據(jù)殘量定理：當(dāng)精化向量收斂到原矩陣的特征值當(dāng)中，則精化向量的殘量收斂值為零.利用Ipsen引理可以得出精化向量和特征向量的殘量值.

(10)

圖3 殘量對比結(jié)果圖Fig.3 The figure of residual comparison

6 結(jié)束語

相對非對稱線性方程組求解算法及應(yīng)用而言，非對稱線性方程組求解初始化策略的研究相對較少，通過提出精化Lanczos算法并與Lanczos算法的各項(xiàng)性能進(jìn)行對比，可以充分地驗(yàn)證出精化算法的優(yōu)勢，從理論上說明精化Lanczos方法比Lanczos方法更優(yōu)越，可以在實(shí)際問題的解決中得到廣泛應(yīng)用.