張 雷

(中國(guó)西南電子技術(shù)研究所，成都610036)

1 引 言

在無(wú)線通信多天線系統(tǒng)中，基于信道矩陣奇異值分解/特征值分解的預(yù)編碼/波束成形技術(shù)［1－2］主要問(wèn)題是計(jì)算半正定Hermitian 矩陣的最大前n 個(gè)特征值對(duì)應(yīng)的特征向量(為方便，將這些特征向量分別稱為第i(i=1，2，…，n)特征向量)。矩陣特征值分解主要有兩大類方法:一類是矩陣變換方法，如Jacobi 方 法［3］;另 一 類 是 向 量 迭 代 法，如 乘 冪法［4－5］、子空間迭代法［6］。Jacobi 方法通常得到矩陣的全部特征值及其對(duì)應(yīng)特征向量，而乘冪法是求解半正定Hermitian 矩陣最大特征值對(duì)應(yīng)特征向量的一種經(jīng)典方法。文獻(xiàn)［5，8］進(jìn)一步指出將乘冪法與壓縮法相結(jié)合可逐個(gè)求出半正定Hermitian 矩陣的所有特征向量。

但是，上述經(jīng)典乘冪法(及其與壓縮法的結(jié)合)存在兩個(gè)主要問(wèn)題:一是迭代次數(shù)固定，使得待分解矩陣的隨機(jī)變化和初始向量的不同選擇導(dǎo)致計(jì)算精度波動(dòng)較大——較小迭代次數(shù)使得部分矩陣計(jì)算精度不高，較大迭代次數(shù)盡管能更好保證計(jì)算精度，但會(huì)顯著增加計(jì)算復(fù)雜度;二是較大特征值對(duì)應(yīng)特征向量的計(jì)算誤差會(huì)傳播至較小特征值對(duì)應(yīng)特征向量的計(jì)算。

著眼于解決以上問(wèn)題，在前述研究基礎(chǔ)上，本文提出了計(jì)算第一特征向量的基于向量距離的改進(jìn)乘冪法，在每次迭代時(shí)計(jì)算本次所得向量與前一次所得向量的距離，當(dāng)該向量距離大于預(yù)先設(shè)定的閾值后即停止迭代;同時(shí)，還證明了所提方法的收斂性。此外，將改進(jìn)乘冪法與壓縮法相結(jié)合，給出了計(jì)算第二特征向量的算法，并證明了存在誤差傳播時(shí)該算法的收斂性。理論計(jì)算結(jié)果表明，對(duì)4 階半正定Hermitian 隨機(jī)矩陣，在相同計(jì)算精度前提條件下，相比經(jīng)典乘冪法(及其與壓縮法的結(jié)合)，所提方法可至少降低一半計(jì)算復(fù)雜度。

符號(hào)說(shuō)明:大小寫(xiě)黑斜體字母分別表示矩陣和(列)向量;(·)T和(·)H分別表示矩陣的轉(zhuǎn)置和Hermitian 轉(zhuǎn)置;‖·‖表示矩陣的Frobenius 范數(shù)。

2 第一特征向量的計(jì)算

2.1 經(jīng)典乘冪法簡(jiǎn)介

對(duì)一m 階半正定Hermitian 矩陣R，設(shè)其m 個(gè)特征值按降序排列為λ1≥λ2≥…≥λm≥0，對(duì)應(yīng)的特征向量分別為u1，u2，…，um。經(jīng)典乘冪法由算法1 給出［4］。

算法1:經(jīng)典乘冪法計(jì)算第一特征向量

Step 1:初始化m 維單位向量x(0);

Step 2:設(shè)k=1，2，…，K1，計(jì)算

則x(k)為所求第一特征向量u1。

考察算法1 的收斂性。因u1，u2，…，um構(gòu)成m維向量空間的一組正交基，x(0)可表示為

因篇幅所限，僅研究λ1為單特征值(對(duì)應(yīng)的特征向量有且僅有一個(gè))的情形。為討論方便，設(shè)

顯然有1 =r1＞r2≥…≥rm≥0 成立。定義x(k)和u1的距離為

它反映了x(k)與u1之間的差異程度，sinθ1，(k)越小，x(k)與u1越接近。一般地，有如下定理［4］:

定理1 若a1≠0，則對(duì)由算法1 求得的x(k)，下式成立(UB 表示上界，下同):

由定理1，只要a1≠0，tanθ1，(0)就為一有限正數(shù)，因此，當(dāng)k→"時(shí)，有sinθ1，(k)，UB→0，即x(k)→ejφu1(φ 為任一實(shí)數(shù))，這表明x(k)將隨著k 的增大而收斂于u1。

2.2 基于向量距離的乘冪法

雖然定理1 保證了算法1 的收斂性，但給定x(0)和K1，不同R 得到的計(jì)算誤差sinθ1，(k)差異較大，主要原因是:x(0)與u1之間距離不同使得|a1|取值不同;R 特征值分布不同使得ri取值不同。為此，可考慮引入一個(gè)合適的“收斂條件”，只要滿足該條件，即停止迭代，并保證對(duì)不同R 得到的計(jì)算誤差sinθ1，(k)均不超過(guò)某一閾值。由此，將算法1 改進(jìn)得算法2。

算法2:基于向量距離的乘冪法計(jì)算第一特征向量

Step 1:初始化m 維單位向量x(0)、d =1、k =1、ε1(ε1＜1);

Step 2:計(jì)算

直到d ＜ε1為止，則x(k)為所求第一特征向量u1。

由定理1 知，只要算法2 中的d 逐漸減小并收斂于某一定值，就能有效控制迭代次數(shù)，并使求出的x(k)逐漸收斂于u1。下面分析d 的收斂性。定義x(k)和x(k－1)之間的距離為

一般地，有如下定理:

定理2 對(duì)由算法2 求得的x(k)和x(k－1)，若a1≠0，則下式成立:

證明:根據(jù)算法2，有

從而

因此，定理2 得證。

由定理2，只要a1≠0，當(dāng)k→"時(shí)，有sinφ1，(k)，UB→0，即x(k)→ejφx(k－1)(φ 為任一實(shí)數(shù))，從而保證了算法2 的收斂性。

3 第二特征向量的計(jì)算

3.1 乘冪法與壓縮法的結(jié)合

由此可見(jiàn)，R 的第二大特征值λ2和第二特征向量u2恰好分別是的第一大特征值和第一特征向量。因此，這種乘冪法與壓縮法的結(jié)合可將“計(jì)算R的第二特征向量”問(wèn)題轉(zhuǎn)化為“計(jì)算的第一特征向量”問(wèn)題，從而再次運(yùn)用算法1 或算法2 求解。與第2 節(jié)類似，本文僅研究λ2為單特征值的情形。此外，考慮到誤差傳播的影響(詳見(jiàn)3.2 節(jié))，在計(jì)算u1時(shí)擬采用算法2，以便對(duì)任意矩陣R 都能將計(jì)算u1的誤差控制在較小水平。

算法3:基于向量距離的乘冪法計(jì)算第二特征向量

Step 1:同算法2;

Step 3:初始化m 維單位向量z(0)、d =1、k =1、ε2ε2()＜1 ;

Step 4:計(jì)算

直到d ＜ε2為止，則z(k)為所求第二特征向量u2。

3.2 誤差傳播的影響

理論上，乘冪法與壓縮法的結(jié)合可依次計(jì)算出R 的所有特征向量，但實(shí)際上總存在誤差，此誤差傳播會(huì)降低后續(xù)特征向量的計(jì)算精度。本小節(jié)分析珘u1≠u1對(duì)計(jì)算u2帶來(lái)的影響。為簡(jiǎn)化推導(dǎo)，假設(shè)珘u1僅含u1和u2兩個(gè)分量(只要ε1值較小，這個(gè)假設(shè)幾乎總是成立)，即

與前文類似，定義z(k)和u2的距離為

設(shè)初始向量

一般地，有如下定理:

定理3 對(duì)由算法3 求得的z(k)(k≥1)，若b1α2≠b2α1，則下式成立(LB 表示下界)

式中，

證明:根據(jù)算法3，有

若b1α2=b2α1，則式(15)退化為

這表明z(k)中不包含u1和u2兩個(gè)分量，有sinθ2，(k)=1?，F(xiàn)考察b1α2≠b2α1情形，由式(15)得

首先，有

其次，有

因此，定理3 得證。

從定理3 可看出:

(1)若λ3＜λ2(即λ2為單特征值)，則當(dāng)k→"時(shí)，有sinθ2，(k)，UB→sinθ2，(k)→sinθ2，(k)，LB;

(2)因sinθ2，(k)，LB為與k 無(wú)關(guān)的定值，故當(dāng)α2≠0 時(shí)，即使k→"，sinθ2，(k)也不會(huì)趨于0;這與2.2 節(jié)求解u1的結(jié)論不同。當(dāng)然，若計(jì)算u1的精度足夠高，使得sinθ2，(k)，LB足夠小，則可使sinθ2，(k)收斂于足夠小的值。

下面再考察算法3 步驟4 中d 的收斂性。定義z(k)和z(k－1)之間的距離為

一般地，有如下定理:

定理4 對(duì)由算法3 求得的z(k)和z(k－1)，下式成立:

式中，

證明:根據(jù)算法3，當(dāng)k≥2 時(shí)，有

式中，p1、c1和c2由式(22)給出，且ci=bi，pi=ri(i=3，4，…，m)。進(jìn)一步，有

由定理4，即使計(jì)算第一特征向量時(shí)存在誤差，sinφ2，(k)，UB也將隨著第二次迭代次數(shù)k 的增加而不斷減小，并最終趨于0，從而保證了算法3 的收斂性。

4 計(jì)算復(fù)雜度分析

在算法3 中，將計(jì)算第一特征向量的Step1 和計(jì)算第二特征向量的Step 4(迭代次數(shù)為K2)都應(yīng)用經(jīng)典乘冪法的算法命名為算法4。設(shè)算法1 和算法2 計(jì)算第一特征向量的迭代次數(shù)分別為K1和L1，算法3和算法4 計(jì)算第二特征向量的迭代次數(shù)分別為L(zhǎng)1(Step 1)和L2(Step 4)、K1(Step 1)和K2(Step 4)，算法1、算法2、算法3 和算法4 的復(fù)數(shù)乘法次數(shù)分別約為K1(m2+3m/4)、L1(m2+7m/4)、(L1+L2+2)m2+7 (L1+L2)m/4 和 (K1+K2+2)m2+3 (K1+K2)m/4，在5.2節(jié)中將結(jié)合隨機(jī)矩陣計(jì)算結(jié)果進(jìn)一步討論。

5 計(jì)算結(jié)果

5.1 確定矩陣計(jì)算結(jié)果

本小節(jié)提供R 為確定矩陣時(shí)的計(jì)算實(shí)例，以驗(yàn)證定理1～定理4 的正確性。R 由式(25)給出:

其特征值從大到小依次為10、6、4、2。

實(shí)例1 用算法1 或算法2 計(jì)算第一特征向量，設(shè)x(0)= [ 1 0 0 0]T。圖1給出了sinθ1，(k)、sinθ1，(k)，UB、sinφ1，(k)和sinφ1，(k)，UB與迭代次數(shù)k 之間的關(guān)系?？梢钥闯?，它們都隨著k 的增大逐漸趨于0，且4 條曲線斜率幾乎一致，這與定理1 和定理2結(jié)論相符。

圖1 sinθ1，(k)、sinθ1，(k)，UB、sinφ1，(k)和sinφ1，(k)，UB與迭代次數(shù)k 之間的關(guān)系Fig.1 sinθ1，(k)，sinθ1，(k)，UB，sinφ1，(k)and sinφ1，(k)，UB vs. iteration number k

實(shí)例2 用算法3 計(jì)算第二特征向量，設(shè)x(0)=z(0)= [ 1 0 0 0]T。圖2給出了閾值ε1取0.1和0. 001 時(shí)sinθ2，(k)、sinθ2，(k)，LB和sinθ2，(k)，UB與Step 4迭代次數(shù)k 之間的關(guān)系?？梢钥闯?對(duì)同一ε1，sinθ2，(k)和sinθ2，(k)，UB都隨著k 的增加而收斂于sinθ2，(k)，LB;對(duì)不同的ε1，sinθ2，(k)，LB隨著ε1的減小而降低，且 大 致 與 ε1的 數(shù) 量 級(jí) 相 當(dāng)，sinθ2，(k)和sinθ2，(k)，UB收斂于sinθ2，(k)，LB所需迭代次數(shù)也逐漸增加，這與定理3 的結(jié)論相符。

圖2 sinθ2，(k)、sinθ2，(k)，LB和sinθ2，(k)，UB與算法3 中Step 4 迭代次數(shù)k 之間的關(guān)系Fig.2 sinθ2，(k)、sinθ2，(k)，LB and sinθ2，(k)，UB vs. iteration number k of Step 4 in Algorithm 3

圖3給出了閾值ε1取0.1 和0.001 時(shí)sinφ2，(k)和sinφ2，(k)，UB與步驟S4 迭代次數(shù)k 之間的關(guān)系?？梢?看 出，在 兩 種 ε1取 值 情 況 下，sinφ2，(k)和sinφ2，(k)，UB都隨著k 的增加而逐漸減小，且斜率幾乎相同，這與定理4 的結(jié)論相符。

圖3 sinφ2，(k)和sinφ2，(k)，UB與算法3 中Step 4 迭代次數(shù)k 之間的關(guān)系Fig.3 sinφ2，(k)and sinφ2，(k)，UB vs. iteration number k of Step 4 in Algorithm 3

5.2 隨機(jī)矩陣計(jì)算結(jié)果

本小節(jié)利用Matlab 隨機(jī)產(chǎn)生半正定Hermitian矩陣，以更好地比較經(jīng)典乘冪法和所提改進(jìn)乘冪法。具體方法為:設(shè)m=4，先用

生成4 階隨機(jī)矩陣A，再用R=AHA 產(chǎn)生R;一共獨(dú)立產(chǎn)生10 000次R。事實(shí)上，當(dāng)應(yīng)用于無(wú)線通信多天線系統(tǒng)時(shí)，A 可視為收發(fā)兩端之間的信道衰落隨機(jī)矩陣，R 的特征向量即為相應(yīng)的預(yù)編碼/波束成形向量。此外，真實(shí)的第一特征向量和第二特征向量由Matlab 函數(shù)eig 求得。

首先考察第一特征向量的計(jì)算。圖4對(duì)比了算法1 和算法2 計(jì)算出的第一特征向量與真實(shí)值距離(用sinθ1表示，由式(3)計(jì)算)的累積分布函數(shù)。對(duì)算法1，由于迭代次數(shù)固定，sinθ1波動(dòng)范圍很大，不易控制計(jì)算精度。K1= 5 和K1= 10 對(duì)應(yīng)的P(sinθ1≥0.1)值分別約為0.18和0.04;當(dāng)K1=20時(shí)，才有P(sinθ1≥0.1)≈0。對(duì)算法2，ε1=0.01和ε1= 0.001 分別對(duì)應(yīng)于P(sinθ1≥0.1)≈0 和P(sinθ1≥0.01)≈0，對(duì)應(yīng)的迭代次數(shù)均值 珔L1分別為7.3和10.9。為比較公平，將P(sinθ1≥0.1)≈0 作為計(jì)算精度要求，由第4 節(jié)分析知，算法1 和算法2的復(fù)數(shù)乘法次數(shù)分別約為380 和168。這表明，相對(duì)于算法1，算法2 不僅能將幾乎所有矩陣的計(jì)算精度控制在較高水平，還降低了約56%的計(jì)算復(fù)雜度。

圖4 sinθ1 的累積分布函數(shù)Fig.4 The cumulative distribution function of sinθ1

接著考察計(jì)算第二特征向量的結(jié)果。圖5給出了算法3 和算法4 求出的第二特征向量與真實(shí)值距離(用sinθ2表示，由式(11)計(jì)算)的累積分布函數(shù)，其特點(diǎn)與圖4類似。對(duì)算法4，sinθ2波動(dòng)范圍很大;K1=K2=5 和K1=K2=10 對(duì)應(yīng)的P(sinθ2≥0.1)值分別升至約0.4 和0.06，當(dāng)K1= K2=20 時(shí)，才有P(sinθ2≥0.1)≈0。對(duì)算法3，ε1=ε2=0.01 和ε1=ε2=0. 001 分別對(duì)應(yīng)于P(sinθ2≥0. 1)≈0 和P(sinθ2≥0.01)≈0，對(duì)應(yīng)的Step 4 的迭代次數(shù)均值珔L2分別為6.1 和8.8。同樣將P(sinθ1≥0.1)≈0 作為計(jì)算精度要求，算法3 和算法4 的復(fù)數(shù)乘法次數(shù)分別約為340 和792，即算法3 比算法4 降低了約57%的計(jì)算復(fù)雜度。

圖5 sinθ2 的累積分布函數(shù)Fig.5 The cumulative distribution function of sinθ2

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

計(jì)算半正定Hermitian 矩陣最大前n 個(gè)特征值對(duì)應(yīng)特征向量是工程中的常見(jiàn)問(wèn)題，在無(wú)線通信多天線系統(tǒng)的預(yù)編碼/波束成形技術(shù)中已有諸多應(yīng)用。相比經(jīng)典乘冪法(及其與壓縮法的結(jié)合)，本文所提基于向量距離的改進(jìn)乘冪法(及其與壓縮法的結(jié)合)能在保證計(jì)算精度的前提下顯著降低計(jì)算復(fù)雜度，更有利于工程實(shí)現(xiàn)。盡管本文僅研究了計(jì)算第一特征向量和第二特征向量情形，但進(jìn)一步通過(guò)3.1節(jié)所述的逐層壓縮，所提方法及理論分析框架也可推廣至計(jì)算所有特征向量情形。

此外，為討論方便，在計(jì)算第i(i =1，2)特征向量時(shí)，本文均假設(shè)第i 特征值為單特征值。文獻(xiàn)［4，8］指出，第i 大特征值為多重特征值的情形并不影響乘冪法及其與壓縮法的結(jié)合。因此，本文所提方法也適用于存在多重特征值的情形。

［1］ 3GPP TS 36.211(V9.1.0)，Evolved Universal Terrestrial Radio Access (E－UTRA);Physical Channels and Modulation (Release 9)［S］.

［2］ 郭建光，王宇欣，孫占委. TD－LTE TM8 傳輸模式分析［J］. 移動(dòng)通信，2014，38(12):33－36.GUO Jianguang，WANG Yuxin，SUN Zhanwei. Research on TD－LTE TM8 Transmission Mode［J］. Mobile Communications，2014，38(12):33－36. (in Chinese)

［3］ Sleupen G，Vander A. A Jacobi－Davidson Iteration Method for Linear Eigenvalue Problem［J］. SIAM Journal on Matrix Analysis and Applications，1996，17(2):401－425.

［4］ Golub G H，Van Loan C F. Matrix Computation［M］.3版.北京:人民郵電出版社，2009.Golub G H，Van Loan C F.Matrix Computation［M］. 3rd ed.Beijing:People's Posts＆Telecom Press，2009.(in English)

［5］ 張賢達(dá). 矩陣分析與應(yīng)用［M］. 北京:清華大學(xué)出版社，2004.ZHANG Xianda. Matrix Analysis and Applications［M］.Beijing:Tsinghua University Press，2004. (in Chinese)

［6］ Bramble J H，Knyazev A V，Pasciak J E. A Subspace Preconditioning Algorithm for Eigenvalue Computation［J］. Advances in Computational Mathematics，1996(6):159－189.

［7］ Paulraj A，Nabar R，Gore D. Introduction to space－time wireless communications (photocopy)［M］. 北京:清華大學(xué)出版社，2005.Paulraj A，Nabar R，Gore D.Introduction to space－time wireless communications (photocopy)［M］. Beijing:Tsinghua University Press，2005.(in English)

［8］ 張青，茍國(guó)楷，呂崇德. 乘冪法的改進(jìn)算法［J］. 應(yīng)用數(shù)學(xué)與計(jì)算數(shù)學(xué)學(xué)報(bào)，1997，11(1):51－55.ZHANG Qing，GOU Guokai，LYU Chongde. An Improvement Scheme for the Power Method［J］. Communication on Applied Mathematics and Computation，1997，11(1):51－55. (in Chinese)