一種調(diào)和Ritz向量的精化算法及應(yīng)用

肖小花 戴 芳 郭文艷

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

在大量的應(yīng)用科學(xué)和工程計算中，如計算流體力學(xué)、量子物理、化學(xué)工程、經(jīng)濟(jì)模擬、材料科學(xué)、結(jié)構(gòu)工程和網(wǎng)絡(luò)排隊的Markov鏈模擬等領(lǐng)域，許多問題最后都?xì)w結(jié)為大規(guī)模矩陣特征值問題。由于大規(guī)模矩陣的階數(shù)很高，而且往往是稀疏矩陣，直接求解其全部特征值不太現(xiàn)實(shí)，在實(shí)際中需要的往往只是矩陣的部分特征對。近幾年來，大規(guī)模矩陣的內(nèi)部特征值問題在許多實(shí)際問題中的應(yīng)用越來越多，成為一個新的研究熱點(diǎn)。

Scott提出的位移求逆的Arnoldi方法對于大規(guī)模矩陣內(nèi)部特征值問題的求解是非常有效的［1］。該方法能夠達(dá)到很快的收斂速度，但由于計算過程中需要對大規(guī)模矩陣進(jìn)行分解，這一方面會破壞原矩陣的稀疏結(jié)構(gòu)，增加存儲量；另一方面，由于Arnoldi算法在迭代過程中對矩陣的擾動是非常敏感的［2］，而求逆過程又很難保證達(dá)到可靠的精度，因此這種方法被認(rèn)為不實(shí)用。由Morgan［3，4］提出的調(diào)和 Arnoldi方法目前被認(rèn)為是解大規(guī)模矩陣內(nèi)部特征問題最有效的方法之一［5］，已成為一種更為實(shí)用的求解方法。

本文利用調(diào)和Arnoldi算法的一種等價形式來求解調(diào)和Ritz對，并應(yīng)用精化Arnoldi算法的思想給出了一種精化變形，進(jìn)一步對調(diào)和Ritz向量進(jìn)行精化求解，尋求使殘量范數(shù)達(dá)到極小的近似特征向量，并對這種方法進(jìn)行了理論分析，給出了數(shù)值實(shí)驗(yàn)。理論分析顯示了這種方法的可行性，數(shù)值實(shí)驗(yàn)顯示這種方法的優(yōu)越性。最后將本文的算法應(yīng)用于K－L變換的變換矩陣求解中。K－L［6］變換的核心過程是計算特征值和特征向量。由于待處理矩陣維數(shù)高，一般的方法很難求出其特征值及特征向量，甚至無法求出。K－L變換的一些優(yōu)化處理過程復(fù)雜，很難滿足實(shí)時性的要求，使得用K－L變換進(jìn)行圖像壓縮難以廣泛應(yīng)用。而本文的算法正好就是一種快速求解大規(guī)模矩陣特征值及特征向量的方法，將其用于K－L變換來進(jìn)行圖像壓縮是可行的。

1 調(diào)和Arnoldi方法的一種實(shí)現(xiàn)形式

對于任意給定的單位初始向量q0∈Rn和實(shí)定點(diǎn)τ，構(gòu)造m維的Krylov子空間

用Q表示V的位移子空間，即Q＝（A－τI）V ，所以有

設(shè)dimV＝dimQ＝m，令r11＝‖（A－τI）q0‖，q1＝ （A－τI）／r11，則有

用Arnoldi過程生成Qm－1的一組標(biāo)準(zhǔn)正交基，記Qm－1＝ ［q1，q2，…，qm－1］，Qm＝ ［Qm－1，qm］，這一過程可以用矩陣形式表示成

其 中 Gm＝ （gij）m×m和 Rm＝ （rij）m×m均為上Hessenberg矩陣，則q1，…，qm構(gòu)成了Q 的一組標(biāo)準(zhǔn)正交基。我們將上述產(chǎn)生Q的一組標(biāo)準(zhǔn)正交基的過程稱為位移Arnoldi過程［7］。

從上面的敘述可以得出：

q0，q1，…，qm－1構(gòu)成了V 的一組基，令Um＝ ［q0，q1，…，qm－1］，則由位移 Arnoldi過程可得

其中Qm＝ ［q1，…，qm］。關(guān)系式（6）揭示了子空間Q的一組標(biāo)準(zhǔn)正交基Qm和V的一組基Um之間的關(guān)系。

利用上述的位移Arnoldi過程和調(diào)和投影得特征值問題［3，8］

2 調(diào)和Ritz向量的精化變形算法

用調(diào)和投影方法來解決大規(guī)模矩陣的部分內(nèi)部特征值問題時，理論分析表明Ritz向量的收斂條件更嚴(yán)格，往往調(diào)和Ritz值作為特征值的近似非常適合，而相應(yīng)的調(diào)和Ritz向量近似程度很差。Jia提出了對于確定的特征值，在近似特征向量的計算上使用滿足某種最優(yōu)條件的精化近似向量來代替?zhèn)鹘y(tǒng)的Ritz向量或者調(diào)和Ritz向量［9，10］，分析表明精化近似特征向量能夠從理論上克服調(diào)和Ritz值收斂到特征值時，調(diào)和Ritz向量經(jīng)常不收斂的問題；同時，他還證明了當(dāng)使用精化調(diào)和Ritz向量來作為特征向量的近似時，能夠得到與相應(yīng)的Ritz值相同的收斂性。本文利用位移Arnoldi算法生成的位移Krylov子空間，結(jié)合精化的思想給出一種精化變形，在位移Krylov子空間中進(jìn)行正交投影來對調(diào)和Ritz向量進(jìn)行精化求解，使改進(jìn)后調(diào)和Ritz向量對應(yīng)的殘量范數(shù)在求解空間上達(dá)到極小。

由前文的介紹可知，位移Krylov子空間Q為m維，則在Q上進(jìn)行正交投影

根據(jù)精化思想，并由關(guān)系式（6）和（10），有以下關(guān)系式成立

結(jié)論：若zi是矩陣Gm，i的最小奇異值σmin（Gm，i）所對應(yīng)的右奇異向量，且ui＝Qmzi，那么，其中－τ）Qm。

由式（10）可知這種精化變形使調(diào)和Ritz向量對應(yīng)的殘量范數(shù)在求解空間上達(dá)到極小，若設(shè)，則，因此精化變形方法比原來的方法能達(dá)到更好的收斂性質(zhì)和更快的收斂速度。下面給出調(diào)和Ritz向量的精化變形算法。

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

（2）生成子空間Q的一組標(biāo)準(zhǔn)正交基向量，得到Um，Qm，Rm。

（3）計算近似特征對。通過式（8）計算近似特征對作為A的特征值的近似，利用ui＝Qmzi計算精化后的調(diào)和Ritz向量，按要求從中選擇k個作為準(zhǔn)確特征對 （λi，φi）的近似（i＝1，2，…，k），其 中zi為 矩 陣Gm，i的 最 小 奇 異 值σmin（Gm，i）所對應(yīng)的右奇異向量；

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

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

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

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

3.1 數(shù)值算例

在INTEL PENTIUM微機(jī)上使用MATLAB7.0軟件包對本文的整個算法進(jìn)行數(shù)值實(shí)驗(yàn)。算法中停機(jī)準(zhǔn)則設(shè)為，實(shí)驗(yàn)中設(shè)精度為tol，當(dāng)stopcrit≤tol時停機(jī)。數(shù)值實(shí)驗(yàn)中所用的大規(guī)模矩陣為圖像矩陣，圖像為圖1和圖2所示的以jpg文件格式存儲的山水和松樹，山水的大小為1 000×1 000，松樹的大小為2 000×2 000。

圖1 山水Fig.1 Landscape

圖2 松樹Fig.2 Pine

數(shù)值算例1以圖1的圖像矩陣作為待計算特征值和特征向量的矩陣，矩陣的規(guī)模為1 000×1 000。m表示Krylov子空間的維數(shù)，m分別取10、20和30；n表示重新啟動的次數(shù)；t表示運(yùn)算時間，以秒為單位；M表示矩陣與向量乘積的個數(shù)。“－”表示算法200步迭代未收斂。tol＝10－7，位移τ＝1，初始向量q0是按均勻分布生成的隨機(jī)向量。按本文算法計算得到矩陣按實(shí)部最大的3個特征值依次近似為λ1≈14.423 5，λ2≈10.211 8，λ3≈8.367 3。按調(diào)和 Arnoldi算法的等價變形計算得到矩陣按實(shí)部最大的3個特征值依次近似為λ1≈14.393 2，λ2≈10.192 6，λ3≈8.287 5。表1給出了計算的結(jié)果。

表1 對1 000階方陣用不同方法計算的結(jié)果比較Table 1 Comparison of the calculated results of 1 000 order matrixes by different methods

數(shù)值算例2此例是以圖2的圖像矩陣作為待計算特征值和特征向量的矩陣，矩陣的規(guī)模為2 000×2 000。m，τ，q0的取值同數(shù)值算例1，精度tol＝10－6。按本文算法計算得到矩陣按實(shí)部最大的3個特征值依次近似為λ1≈15.358 1，λ2≈12.807 9，λ3≈10.394 2。按調(diào)和Arnoldi算法的等價變形計算得到矩陣按實(shí)部最大的3個特征值依次近似為λ1≈15.317 2，λ2≈12.799 1，λ3≈10.295 5。表2給出了計算的結(jié)果。

表2 對2 000階方陣用不同方法計算的結(jié)果比較Table 2 Comparison of the calculated results of 2 000order matrixes by different methods

由以上2例的計算結(jié)果可以得出：對于不同的子空間維數(shù)，本文算法更有效，無論是在運(yùn)算時間還是在達(dá)到收斂時的迭代步數(shù)上均優(yōu)于精化前的調(diào)和Arnoldi算法；特別是當(dāng)子空間維數(shù)越小時，達(dá)到收斂迭代的步數(shù)更少，計算所需的時間更短，優(yōu)越性更明顯。

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

在用K－L變換對圖像進(jìn)行壓縮的實(shí)驗(yàn)中，考慮到計算機(jī)處理器的性能，本文采用了如下方法：首先把圖像矩陣分成了大小相同的幾個矩陣，再用本文算法算出每個矩陣的協(xié)方差矩陣的特征值和特征向量，進(jìn)而對每個矩陣進(jìn)行K－L變換，并把變換后得到的每個壓縮圖像矩陣相應(yīng)地合成為壓縮后的整幅圖像矩陣，最后顯示出壓縮后的整幅圖像。如果計算機(jī)處理器的性能比較好，則用本文的算法就可以直接計算整幅圖像矩陣的協(xié)方差矩陣的特征值和特征向量，對整幅圖像矩陣進(jìn)行K－L變換來壓縮圖像。對于將圖像分成若干小塊而對每個小塊分別進(jìn)行K－L變換的方法，一般選擇大小為8×8的塊。

實(shí)驗(yàn)1對大小為256×256×8bit的以png文件格式存儲的米粒灰度圖像進(jìn)行壓縮。我們比較了以下2種方法：第一種是本文的圖像壓縮方法，將圖像矩陣分成16個64×64的方陣。第二種是將圖像分成若干小塊而對每個小塊分別進(jìn)行變換的方法［11］，記為JK－L，每個小塊選為8×8。t表示算法的執(zhí)行的時間（以分鐘為單位）。在實(shí)驗(yàn)中，本文的算法選取Krylov子空間的維數(shù)m＝30。對特征值的保留個數(shù)分別取為1、2、4、8、16，將特征向量量化為8位二進(jìn)制數(shù)，用2種方法分別進(jìn)行測試的情況如表3所示。

表3 256×256×8bit圖像的壓縮比和峰值信噪比及算法執(zhí)行時間Table 3 256×256×8bit image compression ratio and peak signal／noise ratio and the algorithm time

用本文方法對米粒圖進(jìn)行壓縮的原圖和保留特征值個數(shù)分別為1、2、4、8、16的壓縮后的重建圖像如圖3所示。

將圖像分塊的JK－L方法對米粒圖進(jìn)行壓縮的原圖和保留特征值個數(shù)分別為1、2、4、8、16的壓縮后的重建圖像如圖4所示。

實(shí)驗(yàn)2對大小為512×512×8bit的以jpg文件格式存儲的牽牛花灰度圖像進(jìn)行壓縮。比較了以下2種方法：第一種是本文的圖像壓縮方法，將圖像矩陣分成16個128×128的方陣。第二種是將圖像分成若干小塊而對每個小塊分別進(jìn)行變換的方法，記為JK－L，每個小塊選為8×8。t表示算法的執(zhí)行時間（以分鐘為單位）。在實(shí)驗(yàn)中，本文的算法選取Krylov子空間的維數(shù)m＝30。對特征值保留個數(shù)分別取為1、2、4、8、16，將特征向量量化為8位二進(jìn)制數(shù)，用2種方法分別進(jìn)行測試的情況如表4所示。

圖3 米粒的原圖及保留特征值的重建圖像Fig.3 Original image of rice and the reconstructed image of the eigenvalue

圖4 米粒原圖及保留特征值的重建圖像Fig.4 Original of rice image and the reconstructed image of the eigenvalue

表4 512×512×8bit圖像的壓縮比和峰值信噪比及算法執(zhí)行時間Table 4 512×512×8bit image compression ratio and peak signal／noise ratio and the algorithm time

用本文方法對牽牛花圖進(jìn)行壓縮的原圖和保留特征值個數(shù)分別為1、2、4、8、16的壓縮后的重建圖像如圖5所示。

將圖像分塊的JK－L方法對牽牛花圖進(jìn)行壓縮的原圖和保留特征值個數(shù)分別為1、2、4、8、16的壓縮后的重建圖像如圖6所示。

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

4 結(jié)束語

本文研究了大規(guī)模矩陣特征值問題的一種數(shù)值求解方法，對調(diào)和Arnoldi方法進(jìn)行了改進(jìn)，在精化思想的基礎(chǔ)上給出了一種精化變形方法。保持調(diào)和Ritz值不變，對調(diào)和Ritz向量進(jìn)行了精化求解，理論分析和數(shù)值實(shí)驗(yàn)結(jié)果均表明了該方法的可行性及有效性，并將其用于K－L變換進(jìn)行圖像壓縮。這種將大規(guī)模矩陣特征值問題的數(shù)值求解方法引入K－L變換來進(jìn)行圖像壓縮，解決了KL變換中變換矩陣過大而求解困難的問題。進(jìn)一步對大規(guī)模矩陣特征值問題數(shù)值求解方法的研究，將有助于用K－L變換對圖像進(jìn)行實(shí)時壓縮和傳輸。

圖5 牽牛花的原圖及保留特征值的重建圖像Fig.5 The original image of morning glory and the reconstructed image of the eigenvalue

圖6 牽牛花原圖及保留特征值的重建圖像Fig.6 The original image of morning glory and the reconstructed image of the eigenvalue

［1］Scott D C.The advantages of inverted operators in Rayleigh－Ritz approximations SIAMJ.Sci.Statist［J］.Comput，1982，3：68－75.

［2］Bouras A，F(xiàn)raysse V.A relaxation strategy for the Arnoldi method in eigenproblems，technical report TR／PA／00／16［R］.France：CERFACS，Toul－ouse，2000.

［3］Morgan R B.Computing interior eigenvalues of large matrices［J］.Linear Algebra Appl，1991，154：289－309.

［4］Morgan R B，Zeng M.Harmonic projection methods for large non－symmetric eigenvalue problems［J］.Numer Linear Algebra Appl，1998，5：33－55.

［5］Bai J Z，Barret R，Day D，et al.Test matrix collection for non－Hermitian eigenvalue problems，Technical Report，Department of Mathematic，University of Kentucky，US，1995 ［EB／OL］.（2003－08－15），［2009－11－07］.http／／math.nist.gov／MatrixMarket／.

［6］趙榮椿，趙忠明，崔蘇生.數(shù)字圖像處理導(dǎo)論［M］.西安：西北工業(yè)大學(xué)出版社，1999：66－70.

［7］Walker H F，Zhou L.A simpler GMRES［J］.Numer Linear Algebra APPl，1994，1：57l－581.

［8］Morgan R B.Implicitiy restarted GMRES and Arnoldi methods for nonsymmetrics－tems of equations［J］.SIAMJ Matrix Anal Appl，2000，21：1112－1135.

［9］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.

［10］Jia Z.The convergence of harmonic Ritz values，harmonic Ritz vectors and refined harmonic Ritz vectors［J］.Math Comput，2005，74：1441－1456.

［11］王文峰.K－L變換的研究及其在圖像壓縮編碼中的應(yīng)用［D］.沈陽：沈陽理工大學(xué)檔案館，2008.