趙學智，葉邦彥

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

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

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

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

(1) 計算AAT或ATA的運算量較大；

(2) 在計算AAT或ATA時，由于存在舍入誤差，會導致信息損失，有時甚至會導致計算所得到的奇異值面目全非。

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

2 單向收縮QR算法

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

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

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

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

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

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

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

3 單向收縮QR算法的不收斂實例

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

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

圖2 一個銑削力信號

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

圖3 QR迭代時矩陣末行的收縮過程

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

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

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

表1 前10次QR迭代時前兩個Givens右矩陣的sinθ值

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

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

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

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

表2 第165～第174次QR迭代時前兩個Givens右矩陣的sinθ值

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

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

5 結(jié) 論

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

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

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

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

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

[1] GOLUB G H, VANLOAN C F. 矩陣計算[M]. 袁亞湘譯.北京: 科學出版社, 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ù)值方法與計算機實現(xiàn)[M]. 北京: 清華大學出版社, 2006.XU Shi-liang. Numerical methods and computer realization[M]. Beijing: Tsinghua University Press, 2006.

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

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

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

[13] 趙學智, 葉邦彥. SVD和小波變換的信號處理效果相似性及其機理分析[J]. 電子學報，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.