肖小花

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

0 引言

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

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

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

1 調(diào)和Arnoldi方法

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

Aφi=λiφi(1)

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

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

(2)

(3)

(4)

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

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

Aφi=λiφi(5)

(6)

則有矩陣表達(dá)式：

(7)

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

(8)

(0代表零行向量)

所以有：

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

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

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

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

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

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

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

3.1 數(shù)值算例

圖1

圖2

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

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

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

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

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

3.2 實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

本文對(duì)調(diào)和Arnoldi方法進(jìn)行了改進(jìn).針對(duì)往往調(diào)和Ritz值收斂而相應(yīng)的調(diào)和Ritz向量近似程度很差的情況，保持調(diào)和Ritz值不變，充分利用基向量Vm+1提供的信息求解調(diào)和Ritz向量即改進(jìn)的調(diào)和Ritz向量.理論分析和數(shù)值實(shí)驗(yàn)結(jié)果均表明了該方法的可行性及有效性，并將其用于K-L變換進(jìn)行圖像壓縮.將這種方法引入K-L變換來(lái)進(jìn)行圖像壓縮，解決了K-L變換中變換矩陣過(guò)大而求解困難的問(wèn)題.對(duì)大規(guī)模矩陣特征值問(wèn)題數(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] 趙榮椿，趙忠明，崔蘇生．?dāng)?shù)字圖像處理導(dǎo)論[M]．西安：西北工業(yè)大學(xué)出版社，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].沈陽(yáng)：沈陽(yáng)理工大學(xué)，2008.