趙學(xué)智，葉邦彥

(華南理工大學(xué)機(jī)械與汽車工程學(xué)院 廣州 510640)

近年來奇異值分解(singular value decomposition,SVD)在圖像處理、系統(tǒng)辨識(shí)、控制理論、信號(hào)處理、故障檢測(cè)與診斷等很多領(lǐng)域得到了廣泛的應(yīng)用。將SVD應(yīng)用于工程實(shí)際所需解決的關(guān)鍵是實(shí)現(xiàn)其數(shù)值計(jì)算(一般用QR算法實(shí)現(xiàn))。本文首先系統(tǒng)地總結(jié)了一般矩陣數(shù)值計(jì)算文獻(xiàn)中所介紹的QR算法的特點(diǎn)，給出了一個(gè)該算法在處理由某些小幅度信號(hào)構(gòu)造的大矩陣的奇異值分解時(shí)不收斂的具體實(shí)例，進(jìn)而研究了產(chǎn)生該不收斂現(xiàn)象的本質(zhì)原因。通過對(duì)QR迭代過程中Givens矩陣變化特點(diǎn)的理論分析，對(duì)算法的收斂特性進(jìn)行了深入的解剖，指出該算法在處理大型矩陣的SVD時(shí)不收斂的根本原因在于首行對(duì)角帶元素的衰減，并進(jìn)一步得到了一個(gè)關(guān)于首行元素的衰減對(duì)QR算法收斂速度影響的重要結(jié)論，利用實(shí)例數(shù)據(jù)進(jìn)行證實(shí)，為多次分割、雙向收縮QR算法的提出奠定了基礎(chǔ)。

1 奇異值分解的數(shù)值計(jì)算問題

SVD的數(shù)值計(jì)算并不容易，在一些文獻(xiàn)中，矩陣A的奇異值分解是通過計(jì)算AAT和ATA的特征值和特征向量實(shí)現(xiàn)的[2-5]。由于該二類矩陣是對(duì)稱陣，利用Jacobi方法[2-3]可計(jì)算出其特征值和特征向量。但是，該方法有兩個(gè)缺點(diǎn)：

(1) 計(jì)算AAT或ATA的運(yùn)算量較大；

(2) 在計(jì)算AAT或ATA時(shí)，由于存在舍入誤差，會(huì)導(dǎo)致信息損失，有時(shí)甚至?xí)?dǎo)致計(jì)算所得到的奇異值面目全非。

文獻(xiàn)[1]提出了一種直接計(jì)算SVD的方法，首先利用Householder變換將A化為上雙對(duì)角線矩陣，再利用帶隱式位移的QR算法將上雙對(duì)角線陣化為對(duì)角線陣，可分別得到U、S、V。文獻(xiàn)[6-8]針對(duì)一些特殊的矩陣如正定Toeplitz矩陣、Hessenberg矩陣等，對(duì)該算法提出了一些相應(yīng)的改進(jìn)措施。文獻(xiàn)[9-10]提供了基于文獻(xiàn)[1]提出的QR算法實(shí)現(xiàn)SVD的數(shù)值計(jì)算程序。但是由于該算法本身的原因，這些程序都對(duì)其中的QR迭代總次數(shù)做出了限制，因而其可處理的矩陣階數(shù)絕不可能很大。

2 單向收縮QR算法

該上雙對(duì)角化過程比較簡(jiǎn)單，值得研究的是在完成了上雙對(duì)角化過程后，再對(duì)矩陣B執(zhí)行QR迭代以最終實(shí)現(xiàn)對(duì)角化時(shí)，當(dāng)矩陣階數(shù)很大，在QR迭代時(shí)如果不注意收縮方向的選擇、矩陣的分割、非零元素的不同方向驅(qū)逐等因素，則QR算法就有可能出現(xiàn)不收斂的情況。

QR迭代是通過一系列Givens平面旋轉(zhuǎn)矩陣實(shí)現(xiàn)的，Givens平面旋轉(zhuǎn)矩陣是一個(gè)正交矩陣：

式中c=cosθ；s=sinθ；除了已標(biāo)示出的矩陣元素外，其余位置的元素全部為零。

設(shè)階數(shù)為d的上雙對(duì)角方陣B的通用形式為：

圖1 非零元素的驅(qū)逐過程

式中 0為零矩陣；bd,d為分離出來的一個(gè)奇異值。此時(shí)矩陣可以向上收縮一行，得到矩陣B1，雖然其仍然是一個(gè)上雙對(duì)角陣，但是其階數(shù)已降為d?1。繼續(xù)對(duì)B1執(zhí)行同樣的QR迭代運(yùn)算，矩陣將持續(xù)向上收縮，其階數(shù)不斷降低，新的奇異值不斷地被分離出來，最終可以實(shí)現(xiàn)整個(gè)矩陣的對(duì)角化。

從步驟4可見，執(zhí)行QR迭代時(shí)，矩陣按照從下到上的方向收縮，故可稱此QR迭代為單向收縮QR算法。處理階數(shù)不大的矩陣時(shí)，該算法可實(shí)現(xiàn)較快的SVD計(jì)算。但在基于SVD方法的信號(hào)處理中，經(jīng)常要利用信號(hào)構(gòu)造一個(gè)大階數(shù)的矩陣，其行列數(shù)一般可達(dá)成千上萬[13]。在實(shí)踐中發(fā)現(xiàn)，如果信號(hào)的幅值較小，而構(gòu)造的矩陣階數(shù)又很大，則按照上述算法對(duì)這些大矩陣進(jìn)行SVD運(yùn)算，有時(shí)會(huì)出現(xiàn)不收斂的情況。例如，該算法在處理某些小幅度信號(hào)所構(gòu)造的大型矩陣時(shí)，在矩陣向上收縮了一定的行數(shù)之后，設(shè)此時(shí)的上雙對(duì)角陣為Bi，再對(duì)Bi進(jìn)行QR迭代時(shí)，不管QR迭代過程被執(zhí)行多少次，Bi最右下角的次對(duì)角元卻始終會(huì)停頓在某一個(gè)數(shù)值，再也無法下降，不能滿足式(5)的收斂準(zhǔn)則，矩陣無法收縮下去，這就是以上QR算法在實(shí)際應(yīng)用時(shí)可能出現(xiàn)的問題。實(shí)際上，文獻(xiàn)[9-10]在其提供的計(jì)算SVD的QR算法程序中，也已經(jīng)認(rèn)識(shí)到該算法在處理某些大矩陣時(shí)可能會(huì)出現(xiàn)不收斂的情況，但是該文獻(xiàn)未能解釋出現(xiàn)該問題的原因，也未提出解決措施，僅僅在程序中設(shè)計(jì)了迭代一定次數(shù)后不收斂就強(qiáng)行終止程序的指令。按照最樂觀的估計(jì)，即使一次QR迭代完成后能平均實(shí)現(xiàn)兩行的向上收縮，算法程序所能處理的矩陣階數(shù)最多為QR迭代總次數(shù)的兩倍，實(shí)際處理時(shí)遠(yuǎn)遠(yuǎn)不可能達(dá)到這一目標(biāo)，因此它們都難以用于計(jì)算大型矩陣的奇異值分解。根據(jù)奇異值理論，對(duì)于任何矩陣，不管其階數(shù)有多大、行數(shù)與列數(shù)孰大孰小，其奇異值分解總是存在且總是可以實(shí)現(xiàn)其數(shù)值計(jì)算的，算法出現(xiàn)不收斂的情況只能說明算法存在問題，有待改進(jìn)。

3 單向收縮QR算法的不收斂實(shí)例

本文通過一個(gè)具體實(shí)例證明該QR算法在處理由某些實(shí)際信號(hào)所構(gòu)造的大型矩陣時(shí)會(huì)出現(xiàn)不收斂的情況。

圖2所示為一個(gè)長(zhǎng)度為1 024點(diǎn)的銑削力信號(hào)，以該信號(hào)構(gòu)造一個(gè)行數(shù)512、列數(shù)513的Hankel矩陣。通過Householder變換將其化為上雙對(duì)角陣后，利用QR算法對(duì)其進(jìn)行對(duì)角化。開始時(shí)矩陣向上收縮的速度很快，一般只需執(zhí)行一次QR迭代(最多只執(zhí)行3次)，矩陣就可向上收縮一行。當(dāng)向上收縮了113行即收縮到第399行時(shí)，QR迭代總共被執(zhí)行了168次后，此后再也無法滿足收縮條件。執(zhí)行到第174次QR迭代后，上雙對(duì)角陣Bi最右下角的次對(duì)角元b398,399的值停頓在2.019 732×10?5，再也無法下降，算法陷入無限循環(huán)。

圖2 一個(gè)銑削力信號(hào)

圖3所示為前200次QR迭代時(shí)，矩陣末行向上收縮的變化情況。圖中清晰地反映了剛開始時(shí)矩陣末行快速向上收縮以及之后收縮終止的過程。

圖3 QR迭代時(shí)矩陣末行的收縮過程

4 單向收縮QR算法的收斂特性分析

本文通過對(duì)QR迭代過程中Givens矩陣變化特點(diǎn)的研究和分析，發(fā)現(xiàn)引起上述問題的根源在于第一個(gè)Givens右矩陣。設(shè)矩陣向上收縮了一定行數(shù)之后的上雙對(duì)角陣為Bi，在對(duì)其進(jìn)行QR變換時(shí)，根據(jù)式(2)可得QR變換中的第一個(gè)Givens右矩陣中c和s的值分別為：

從表1可見，s第6個(gè)值的數(shù)量級(jí)是10?10，非常接近零，表明第6次QR迭代開始時(shí)G1(1,2,)θ已接近為單位陣。但是，從圖4可以看出，此時(shí)上雙對(duì)角陣Bi的末行卻依然能以很快的速度向上收縮。可是根據(jù)前面的分析，當(dāng)?shù)谝粋€(gè)Givens右矩陣G1為單位陣時(shí)，后續(xù)的其他一系列左、右Givens變換矩陣肯定是單位陣，此時(shí)式(4)的QR變換將失去意義，Bi再也不可能向上收縮了。

表1 前10次QR迭代時(shí)前兩個(gè)Givens右矩陣的sinθ值

由圖2可見，構(gòu)造矩陣的銑削力信號(hào)中所有數(shù)據(jù)的數(shù)量級(jí)都是一樣的，也即矩陣A中所有元素的數(shù)量級(jí)都相同，而經(jīng)過Householder變換后所得的上雙對(duì)角陣B的兩個(gè)對(duì)角帶元素的數(shù)量級(jí)仍然與矩陣A是一致的，則由式(8)易見，s′和s的數(shù)量級(jí)是相同的。因此，當(dāng)?shù)谝粋€(gè)Givens右矩陣接近單位陣時(shí)，第一個(gè)Givens左矩陣也將接近單位陣，其變換的結(jié)果是在位置(1,3)處產(chǎn)生新元素：

從式(7)可知，第一個(gè)Givens右矩陣的s的數(shù)量級(jí)是由b1,1和b1,2的乘積決定的，當(dāng)s較小時(shí)，表明b1,1或者b1,2的數(shù)量級(jí)也較小。由于b1,1是主對(duì)角帶元素，為未來的奇異值，其衰減的可能性比b1,2小很多，因此更可能是b1,2的數(shù)量級(jí)與s接近，則式(11)中的比值b1,3/b1,2的數(shù)量級(jí)反而會(huì)增大。例如假設(shè)s的數(shù)量級(jí)為10?11，由于構(gòu)造矩陣的信號(hào)數(shù)據(jù)的數(shù)量級(jí)為10?1，可以認(rèn)為此時(shí)b1,1的數(shù)量級(jí)亦為10?1，則由式(7)可得此時(shí)b1,2的數(shù)量級(jí)為10?10；而由前述可知，b1,3的數(shù)量級(jí)將最多比s的數(shù)量級(jí)小10?1，因此可以認(rèn)為b1,3的數(shù)量級(jí)為10?12，此時(shí)b1,3和b1,2都很接近于零；然而b1,3/b1,2的數(shù)量級(jí)卻為10?2，則由式(11)可知第二個(gè)Givens右矩陣已經(jīng)不是單位陣。

本文給出該例中第二個(gè)Givens右矩陣的前10個(gè)sinθ值，將它們列在表1中第3列。可見第10次QR迭代開始時(shí)，第一個(gè)Givens右矩陣的sinθ的數(shù)量級(jí)已下降為10?17，但第二個(gè)Givens右矩陣中sinθ的數(shù)量級(jí)卻還只為10?2，證明以上的分析是正確的。這種結(jié)果將直接導(dǎo)致從第二個(gè)Givens右矩陣開始的后續(xù)一系列左、右Givens變換矩陣都不會(huì)是單位陣，從而使得以后的QR變換可正常進(jìn)行，因此矩陣仍可快速向上收縮。

但是當(dāng)矩陣階數(shù)大到一定程度時(shí)，如果原始矩陣中數(shù)據(jù)的數(shù)量級(jí)較小，有時(shí)在還沒有完成整個(gè)矩陣的對(duì)角化的情況下，隨著QR迭代次數(shù)的增加，第一個(gè)Givens右矩陣的sinθ有可能會(huì)完全下降為零，這會(huì)視原始矩陣數(shù)據(jù)的數(shù)量級(jí)而定。當(dāng)該數(shù)量級(jí)較大時(shí)很難出現(xiàn)這種情況。但是一旦出現(xiàn)這種情況后，所有的左、右Givens變換矩陣都會(huì)變?yōu)閱挝魂?，上雙對(duì)角陣Bi最右下角的次對(duì)角元就會(huì)停頓在某一個(gè)數(shù)值，再也不能下降，算法陷入無限循環(huán)。表2給出了該例中QR迭代到第165～第174次期間，第一、二個(gè)Givens右矩陣中的sinθ值的變化過程。由表可見，在用8個(gè)字節(jié)表示一個(gè)實(shí)數(shù)的double數(shù)據(jù)類型編程的情況下，s的數(shù)量級(jí)可下降為10?323，而此時(shí)第二個(gè)Givens右矩陣的sinθ值的數(shù)量級(jí)卻僅為10?5。迭代到第174次時(shí)，s終于完全降為零，第二個(gè)Givens右矩陣的sinθ值也隨即變?yōu)榱?。此后所有的左、右Givens變換矩陣都變?yōu)閱挝魂?，QR算法失效。

表2 第165～第174次QR迭代時(shí)前兩個(gè)Givens右矩陣的sinθ值

圖4 中的cosθ和sinθ值的變化過程

總結(jié)上述分析，可以得到關(guān)于QR算法的收斂特性的論斷：在對(duì)雙對(duì)角陣執(zhí)行QR迭代時(shí)，除非第一個(gè)Givens右矩陣中的s完全為零，否則即使該s的數(shù)量級(jí)很小。但是，從第二個(gè)Givens右矩陣開始，后續(xù)的一系列左、右Givens變換矩陣卻都不會(huì)相應(yīng)地接近單位陣，從而可以使后面的QR迭代能正常進(jìn)行。這就是在第一個(gè)Givens右矩陣中s的數(shù)量級(jí)相當(dāng)小，即該矩陣非常接近單位陣的情況下，QR算法仍能實(shí)現(xiàn)矩陣快速向上收縮的原因。

5 結(jié) 論

(1) 在處理由某些小幅度信號(hào)構(gòu)造的大型矩陣的奇異值分解時(shí)，單向收縮QR算法會(huì)出現(xiàn)不收斂的情況。

(2) 不收斂的本質(zhì)原因在于矩陣首行對(duì)角帶元素的衰減，當(dāng)構(gòu)造原始矩陣的信號(hào)數(shù)據(jù)的數(shù)量級(jí)較小時(shí)，這種衰減最終會(huì)使QR迭代中的第一個(gè)Givens右矩陣變?yōu)閱挝魂?，從而?dǎo)致后面的一系列Givens變換矩陣全部成為單位陣，使得QR算法失效。

(3) 當(dāng)?shù)谝粋€(gè)Givens右矩陣十分接近單位陣時(shí)，從第二個(gè)Givens右矩陣開始，后續(xù)的一系列左、右變換矩陣卻都不會(huì)相應(yīng)地接近單位陣，使得QR算法仍然有很快的收斂速度。只有當(dāng)?shù)谝粋€(gè)Givens右矩陣完全成為單位陣時(shí)，后續(xù)所有的Givens變換矩陣才會(huì)成為單位陣，QR算法才會(huì)完全失效。

上述工作為多次分割、雙向收縮QR算法的提出指明了方向。該算法還涉及較多的與SVD收斂速度和精度相關(guān)的細(xì)節(jié)和技巧，其具體實(shí)現(xiàn)將另文進(jìn)行探討。

本文的研究工作得到廣州市科技計(jì)劃項(xiàng)目(2008J1-C101)的資助，在此表示感謝。

[1] GOLUB G H, VANLOAN C F. 矩陣計(jì)算[M]. 袁亞湘譯.北京: 科學(xué)出版社, 2001.GOLUB G H, VANLOAN C F. Matrix computations[M].Translated by YUAN Ya-xiang. Beijing: Science Press,2001.

[2] DRMAC Z, VESELIC K. New fast and accurate Jacobi SVD algorithm(I)[J]. SIAM Journal on Matrix Analysis and Application, 2006, 29(4): 1322-1342.

[3] DRMAC Z, VESELIC K. New fast and accurate Jacobi SVD algorithm (II) [J]. SIAM Journal on Matrix Analysis and Application, 2006, 29(4): 1343-1362.

[4] NORDBERG T, GUSTAFSSON I. Using QR factorization and SVD to solve input estimation problems in structural dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(7): 5891-5908.

[5] VANLANDUIT S, CAUBERGHE B, GUILLAUME P.Reduction of large frequency response function data sets using robust singular value decomposition[J]. Computers and Structures, 2006, 84(12): 808-822.

[6] MASTRONARDI N, VANBAREL M, VANDEBRIL R. A fast algorithm for computing the smallest eigenvalue of a symmetric positive-definite Toeplitz matrix[J]. Numerical Linear Algebra with Applications, 2008, 15(4): 327-337.

[7] NIE Y Y, LI Z, HAN J D. Origin-shifted algorithm for matrix eigenvalues[J]. International Journal of Computer Mathematics, 2008, 85(9): 1397-1411.

[8] EIDELMAN Y, GOHBERG I, GEMIGNANIL L. On the fast reduction of a quasiseparable matrix to Hessenberg and tridiagonal forms[J]. Linear Algebra and Its Applications,2007, 420(1): 86-101.

[9] 徐士良. 數(shù)值方法與計(jì)算機(jī)實(shí)現(xiàn)[M]. 北京: 清華大學(xué)出版社, 2006.XU Shi-liang. Numerical methods and computer realization[M]. Beijing: Tsinghua University Press, 2006.

[10] 丁 軍, 楊麗麗. 科學(xué)與工程數(shù)值算法(Java版)[M]. 北京: 清華大學(xué)出版社, 2003.DING Jun, YANG Li-li. Science and engineering numerical algorithms (Java Edition)[M]. Beijing: Tsinghua University Press, 2003.

[11] 奚梅成. 數(shù)值分析方法[M]. 合肥: 中國科學(xué)技術(shù)大學(xué)出版社, 2007.XI Mei-cheng. Numerical analysis methods[M]. Hefei:China Science and Technology University Press, 2007.

[12] 李慶揚(yáng), 王能超, 易大義. 數(shù)值分析[M]. 北京: 清華大學(xué)出版社, 2001.LI Qing-yang, WANG Neng-chao, YI Da-yi. Numerical analysis[M]. Beijing: Tsinghua University Press, 2001.

[13] 趙學(xué)智, 葉邦彥. SVD和小波變換的信號(hào)處理效果相似性及其機(jī)理分析[J]. 電子學(xué)報(bào)，2008, 36(8): 1582-1589.ZHAO Xue-zhi, YE Bang-yan. The similarity of signal processing effect between SVD and wavelet transform and its mechanism analysis[J]. Acta Electronica Sinica, 2008,36(8): 1582-1589.