肖小花 戴 芳 郭文艷

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

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

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

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

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

對(duì)于任意給定的單位初始向量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過(guò)程生成Qm－1的一組標(biāo)準(zhǔn)正交基，記Qm－1＝ ［q1，q2，…，qm－1］，Qm＝ ［Qm－1，qm］，這一過(guò)程可以用矩陣形式表示成

其 中 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)正交基的過(guò)程稱為位移Arnoldi過(guò)程［7］。

從上面的敘述可以得出：

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

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

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

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

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

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

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

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

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

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

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

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

（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ù)值算例

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

圖1 山水Fig.1 Landscape

圖2 松樹(shù)Fig.2 Pine

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 牽?；ǖ脑瓐D及保留特征值的重建圖像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］.沈陽(yáng)：沈陽(yáng)理工大學(xué)檔案館，2008.