一種改進的調(diào)和Arnoldi算法及其應(yīng)用

肖小花

(西安理工大學 理學院， 陜西 西安 710054)

0 引言

20世紀50年代以來，求解大規(guī)模矩陣特征值問題主要采用的方法為投影類方法.目前，投影類方法中出現(xiàn)了很多比較有代表性的算法，其中Arnoldi-type算法[1]被公認為是投影類方法中最成功的方法之一.傳統(tǒng)Arnoldi方法[2]雖然對求解大規(guī)模矩陣的端部特征對有著很好的效果，但是對矩陣的內(nèi)部特征值問題常常失效.后來paige等人提出的調(diào)和Arnoldi方法已經(jīng)成為求解大規(guī)模矩陣特征對，尤其是內(nèi)部部分特征對的最重要的方法之一.

本文根據(jù)投影子空間中所含的特征信息越豐富，該子空間上Ritz對的近似程度越好的指導思想，結(jié)合調(diào)和Ritz值收斂而調(diào)和Ritz向量卻難以收斂的情況，充分利用Arnoldi過程產(chǎn)生的第m+1個基向量提供的有用信息，在m+1維的Krylov子空間中來求解近似向量，使Ritz對的殘量范數(shù)達到極小.并對這種方法進行了理論分析，給出了數(shù)值實驗.理論分析顯示這種方法的可行性；數(shù)值實驗顯示這種方法的有效性.最后將本文的算法應(yīng)用于K-L變換的變換矩陣求解中，K-L變換的核心過程是計算特征值和特征向量.由于待處理矩陣維數(shù)高的情況，一般的方法很難求出其特征值及特征向量，甚至無法求出.而K-L變換[5]的一些優(yōu)化處理過程復(fù)雜，很難滿足實時性的要求，使得用K-L變換進行圖像壓縮難以廣泛應(yīng)用.而本文的算法正好就是一種較快求解大規(guī)模矩陣特征值及特征向量的方法，實驗表明將其用于K-L變換來進行圖像壓縮是可行的.

文中Km(A,v1)=span{q0,Aq0,…,Am-1q0}表示m維的Krylov子空間，上標*表示矩陣或向量的共軛轉(zhuǎn)置，I表示N階單位矩陣，em表示m階單位矩陣IM的第M列，σmin(G)表示G的最小奇異值.文中的范數(shù)‖*‖如無特殊說明均為2范數(shù).

1 調(diào)和Arnoldi方法

大規(guī)模矩陣特征值問題：

Aφi=λiφi(1)

的數(shù)值求解，其中A為n階實或復(fù)矩陣，(λi,φi)i=1,2,…,n為A的特征對(‖φi‖=1).

給定m維Krylov子空間Km(A,v1)由Arnoldi過程產(chǎn)生它的一組標準正交基，這一過程的矩陣表達式為：

(2)

(3)

(4)

2 改進的調(diào)和ARNOLDI算法

考慮大規(guī)模矩陣特征值問題：

Aφi=λiφi(5)

(6)

則有矩陣表達式：

(7)

給(7)式兩邊同時左乘(A-τI)得到：

(8)

(0代表零行向量)

所以有：

從上面的討論，可得以下結(jié)論:

根據(jù)以上過程，下面給出整個算法如下:

(1)給定子空間維數(shù)m，需要計算特征向量的個數(shù)k以及要求達到的精度tol，選擇一個單位初始向量q0以及位移點τ；

(5)重新啟動：構(gòu)造新的初始向量q0，轉(zhuǎn)向第(2)步.

注:算法第(5)步采用顯式重新啟動策略，即重新啟動后的初始向量q0取作：

其中α為規(guī)范化因子.

3 數(shù)值算例與實驗結(jié)果

3.1 數(shù)值算例

圖1

圖2

數(shù)值算例1．以圖1中的圖像矩陣作為待計算特征值和特征向量的矩陣，矩陣的規(guī)模為1000×1 000.按本文方法計算得到矩陣按實部最大的3個特征值依次近似為：λ1≈12.3285，λ2≈9.5228,λ3≈8.2896.表1給出了計算結(jié)果，其中m表示子空間的維數(shù)，iter表示重新啟動的次，time表示運算時間(秒).mv表示矩陣與向量乘積的個數(shù).tol=10-7,位移τ=1，初始向量q0是按均勻分布生成的隨機向量.

表1 數(shù)字圖像轉(zhuǎn)換生成的1 200×1 200矩陣

數(shù)值算例2．此例中是以圖2中的圖像矩陣作為待計算特征值和特征向量的矩陣,矩陣的規(guī)模為2 400×2 400.按本文方法算得矩陣按實部最大的三個特征值依次近似為：λ1≈13.358 1，λ2≈10.807 9，λ3≈7.394 2.表2給出了計算的結(jié)果，其中m、iter、 time、mv的表示同數(shù)值算例1.位移τ=1，精度tol=10-6，初始向量q0是按均勻分布生成的隨機向量.

表2 對2 400階方陣用不同方法計算的結(jié)果比較

由以上兩例的計算結(jié)果可以得出：當子空間維數(shù)越小時，本文算法達到收斂迭代的步數(shù)越少，優(yōu)越性比較明顯.隨著子空間維數(shù)的增加，子空間中含有的特征信息比較豐富，則兩種方法的收斂速度都比較快且比較接近.

3.2 實驗

在用K-L變換對圖像進行壓縮的實驗中，需考慮計算機處理器的性能，如果計算機處理器的性能比較好，則用本文的算法就可以直接對整幅圖像進行K-L變換來壓縮圖像；如果計算機處理器的性能不是太好，則可以把圖像矩陣分成大小相同的幾個矩陣，分別對每個矩陣進行K-L變換，再把得到的每個壓縮圖像矩陣相應(yīng)的合成為壓縮后的整幅圖像矩陣.對于將圖像分成若干小塊而對每個小塊分別進行K-L變換的方法[9]，一般選擇大小為8×8的塊.

實驗1對大小為200×200×8 bit的灰度圖像shanshui.jpg進行壓縮.比較了以下兩種方法：第一種是本文方法即將本文的算法直接用于圖像矩陣進行K-L變換；第二種是將圖像分成若干小塊，而對每個小塊分別進行變換[10]的方法，每個小塊選為8×8.在實驗中，本文的算法選取Krylov子空間 的維數(shù)m=30 ,V表示算法的執(zhí)行速度(以分為單位).對特征值的保留個數(shù)分別取為1、2、4、8、16，將特征向量量化為8位二進制數(shù)，用兩種方法分別進行測試的情況如表3所示.

表3 200×200×8 bit圖像的壓縮比和峰值信嗓比及算法執(zhí)行速度

用本文算法對shanshui.jpg進行壓縮的原圖和保留特征值個數(shù)分別為1，2，4，8，16的壓縮后的重建圖像如圖3所示.

將圖像分塊的方法對shanshui.jpg進行壓縮的原圖和保留特征值個數(shù)分別為1，2，4，8，16的壓縮后的重建圖像如圖4所示.

圖3 shanshui.jpg原圖及保留特征值個數(shù)為1、2、4、8、16的重建圖像

圖4 shanshui.jpg原圖及保留特征值個數(shù)為1、2、4、8、16的重建圖像

圖5 fengye.jpg原圖及保留特征值個數(shù)為1、2、4、8、16的重建圖像

圖6 fengye.jpg原圖及保留特征值個數(shù)為1、2、4、8、16的重建圖像

實驗2對大小為1 000×1 000×8 bit的灰度圖像fengye.jpg進行壓縮,比較了以下兩種方法：第一種是本文方法即將本文的算法直接用于圖像矩陣進行K-L變換,將圖像矩陣分成25個200×200的方陣；第二種是將圖像分成若干小塊而對每個小塊分別進行K-L變換,每個數(shù)據(jù)塊選為8×8.在實驗中，本文的算法選取Krylov子空間的維數(shù)m=30 ,V表示算法的執(zhí)行速度(以分為單位).對特征值保留個數(shù)分別取為1，2，4，8，16時，將特征向量量化為8位二進制數(shù)，用兩種方法分別進行測試的情況如表4所示.

表3 200×200×8 bit圖像的壓縮比和峰值信嗓及算法執(zhí)行速度

用本文算法對fengye.jpg進行壓縮的原圖和保留特征值個數(shù)分別為1，2，4，8，16的壓縮后的重建圖像如圖5所示.

將圖像分塊的方法對fengye.jpg進行壓縮的原圖和保留特征值個數(shù)分別為1，2，4，8，16的壓縮后的重建圖像如圖6所示.

從以上兩個實驗的結(jié)果可以看出：從壓縮比、峰值信嗓比以及算法的執(zhí)行時間上，用本文的方法來使用K-L變換進行圖像壓縮要優(yōu)于將圖像數(shù)據(jù)矩陣分成若干小塊而在每個小塊上使用K-L變換的方法.

4 結(jié)束語

本文對調(diào)和Arnoldi方法進行了改進.針對往往調(diào)和Ritz值收斂而相應(yīng)的調(diào)和Ritz向量近似程度很差的情況，保持調(diào)和Ritz值不變，充分利用基向量Vm+1提供的信息求解調(diào)和Ritz向量即改進的調(diào)和Ritz向量.理論分析和數(shù)值實驗結(jié)果均表明了該方法的可行性及有效性，并將其用于K-L變換進行圖像壓縮.將這種方法引入K-L變換來進行圖像壓縮，解決了K-L變換中變換矩陣過大而求解困難的問題.對大規(guī)模矩陣特征值問題數(shù)值求解方法的討論和研究，將有助于K-L變換在圖像壓縮中的廣泛應(yīng)用.

[1] Jia Z．Arnoldi type algorithms for large unsymmetric multipleeig envalue problems[J]．J Comput Math，1999，17：257-274．

[2] Saad Y．Variations on Arnoldi′s method for computing eigenelements of large unsymmetric matrices[J]．Linear Algebra Appl，1980，34：269-295．

[3] Sorensen D C．Implicit application of polynomial filters in a k-step Arnoldi method[J]．Siamj Matrix Anal Appl,1992,13：357-385．

[4] Jia Z．The convergence of harmonic Ritz values，harmonic Ritz vectors and refined harmonic Ritz vectors[J]．Math Comput,2005，74：1 441-1 456．

[5] 趙榮椿，趙忠明，崔蘇生．數(shù)字圖像處理導論[M]．西安：西北工業(yè)大學出版社，1999.

[6] Morgan R B．Zeng M．Harmonic projection methods for large non-symmetric eigenvalue problems[J]．Numer Linear Algebra Appl, 1998，5：33-55．

[7] Jia Z．The refined harmonic Arnoldi method and an implicitly restarted refined algorithm for computing interior eigenpairs of large matrices[J]．Appl Numer Math，2002，42：489-512．

[8] Morgan R B. Implicitiy restarted GMRES and Arnoldi methods for nonsymmetrics-tems of equations[J].Siamj Matrix Anal Appl,2000,21:1 112-1 135.

[9] 王文峰.K-L變換的研究及其在圖像壓縮編碼中的應(yīng)用[D].沈陽：沈陽理工大學，2008.