施育鑫,魯信金,孫藝夫,雷 菁,李玉生
(1.國防科技大學(xué)第六十三研究所,江蘇 南京 210000;2.國防科技大學(xué) 電子科學(xué)學(xué)院,湖南 長沙 410000)
矩陣常常被用于表示無線通信、圖像視頻處理、計算機視覺、相控陣?yán)走_(dá)的數(shù)據(jù)表示。隨著用戶需求的不斷增加,數(shù)據(jù)規(guī)模、通信陣列的持續(xù)擴大,矩陣的大小和維度也隨之快速增加,這給矩陣的數(shù)據(jù)存儲、算法計算帶來了很大的困難。矩陣的關(guān)聯(lián)性是指矩陣數(shù)據(jù)在某些維度上的相關(guān)特性,例如視頻中時間相鄰的幀具有很強的矩陣關(guān)聯(lián)性。因此,充分挖掘矩陣間關(guān)聯(lián)性,以實現(xiàn)低復(fù)雜度的計算具有十分重要的價值和意義。
對于所給定復(fù)數(shù)矩陣H={Hj,k},Hj,k∈M×N,j=1,2,…,J;k=1,2,…,K。其中,矩陣之間以及同一矩陣的元素之間有一定的相關(guān)性,包括:相同j下標(biāo)、不同k下標(biāo)的矩陣間存在一定的關(guān)聯(lián),即{Hj,1,Hj,2,Hj,3,…,Hj,K}間存在關(guān)聯(lián);且矩陣的各個元素間也存在關(guān)聯(lián),矩陣Hj,k∈M×N,j=1,2,…,J;k=1,2,…,K可表示為:
(1)
此外定義矩陣組H={Hj,k}的一組數(shù)學(xué)運算,其中間結(jié)果V={Vj,k}由式(2)給出:
(2)
式中,j=1,2,…,J;k=1,2,…,K,svd(·)為矩陣奇異值分解(Singular Value Decomposition,SVD)中求解右奇異向量的過程;Vj,k是由Hj,k的前L個右奇異向量構(gòu)成的矩陣,維度為N×L。
為得到最終輸出結(jié)果W={Wj,k},先將不同j下標(biāo)、相同k下標(biāo)的Vj,k進(jìn)行橫向的拼接,得到維度N×LJ的Vk=[V1,k,…Vj,k,…VJ,k],然后根據(jù)式(3)獲取Wk:
(3)
式中,σ2為固定常數(shù);Wk維度同Vk;I為單位矩陣,維度為LJ×LJ。
最后將各Wk按如式(4)進(jìn)行拆解:
Wk=[W1,k,…Wj,k,…WJ,k],
(4)
式中,Wj,k為Wk中順序排列的子矩陣,維度為N×L。
為了降低計算和儲存的復(fù)雜度,分析相關(guān)矩陣組的關(guān)聯(lián)性,通過建模對輸出結(jié)果進(jìn)行估計,建模過程可用式(5)表示:
(5)
(6)
定義W的建模估計精度為:
(7)
為描述方便,額外定義W的最低建模精度為ρmin(W):
(8)
計算復(fù)雜度定義為由矩陣組H計算得到結(jié)果矩陣組W所需要的總計算復(fù)雜度。復(fù)數(shù)矩陣運算可拆解為基本的復(fù)數(shù)運算,而基本的復(fù)數(shù)運算又可進(jìn)一步拆解為基本的實數(shù)運算。例如,復(fù)數(shù)乘法(a+bj)(c+dj)=(ac-bd)+(ad+bc)j的復(fù)雜度為4次實數(shù)乘法和2次實數(shù)加(減)法。實數(shù)基本運算的復(fù)雜度按照表1計算。
表1 實數(shù)基本運算的計算復(fù)雜度
2021研究生數(shù)學(xué)建模A題提供的數(shù)據(jù)集附件(Data1~Data6)給出的詳細(xì)數(shù)據(jù),包括輸入矩陣組、標(biāo)準(zhǔn)中間矩陣組和標(biāo)準(zhǔn)輸出矩陣組的數(shù)據(jù)及其維度,其中M=4,N=64,L=2,J=4,K=384,σ2=0.01,數(shù)據(jù)為十進(jìn)制格式。根據(jù)所給數(shù)據(jù)Data1~Data6中的M=4,N=64,J=4,可對應(yīng)通信模型中共有J=4個用戶,每個用戶的發(fā)射天線數(shù)為M=4,基站的接收天線數(shù)為N=64。L=2表示取信道衰落程度最低的2個信道向量,即對矩陣進(jìn)行壓縮。K表示信道探測時隙數(shù)。σ2表示信道中高斯白噪聲的噪聲方差。
建立基于矩陣的多輸入多輸出(Multiple-Input Multiple-Output,MIMO)通信模型[1]如圖1所示,J個用戶發(fā)送信息,信號經(jīng)過信道H到達(dá)基站,基站有N根接收天線。
圖1 矩陣關(guān)系的通信模型建立示意圖Fig.1 Schematic diagram of establishing communication model of matrix relationship
圖2 單個用戶和基站的通信示意圖Fig.2 Schematic diagram of communication between a single user and a base station
此時,可使用壓縮矩陣Vj,k來表示H矩陣。進(jìn)一步,引入最小均方誤差(Minimum Mean Square Error,MMSE)的概念來解釋題設(shè)條件。如圖3所示的信道模型中,信號x經(jīng)過信道V=Vj,k,由于白噪聲n的影響,接收信號y可表示為:
y=Vx+n。
(9)
圖3 信號傳輸模型(求解MMSE流程)Fig.3 Signal transmission model (process of solving MMSE)
(10)
此時的MMSE為:
MMSE=E{eHe}。
(11)
假設(shè)接收到的數(shù)據(jù)y和誤差e是不相關(guān)的,即
E{e·yH}=0。
(12)
將式(10)代入式(12)可得:
E{(Wy-x)·yH}=0。
(13)
將式(13)左邊進(jìn)一步展開可得:
E{(Wy-x)·yH}=E{WyyH}-E{xyH}=
WE{yyH}-E{xyH}。
(14)
由式(13)和式(14)可得:
W=E{xyH}E{yyH}-1。
(15)
接下來對E{yyH}和E{xyH}進(jìn)行處理,首先對于E{yyH},將其進(jìn)一步展開:
E{yyH}=E{(Vx+n)(Vx+n)H}=
E{VxxHVH+VxnH+nxHVH+nnH}。
(16)
此處假設(shè)輸入信號x和噪聲n不相關(guān),則nxH與nxH值為0,可得:
E{yyH}=E{VxxHVH+nnH}=
VE{xxH}VH+E{nnH}=
V(P·I)VH+σ2·I,
(17)
式中,P為發(fā)送信號x的能量,σ2為噪聲n的方差。
其次對于E{xyH},展開如下:
E{xyH}=E{x(Vx+n)H}=
E{xxHVH+xnH}=
E{xxHVH}=
E{xxH}VH=
(P·I)VH。
(18)
得到E{yyH}和E{xyH}后,將其代入式(15)可得到W的表達(dá)式:
W=E{xyH}E{yyH}-1=
(P·I)VH(V(P·I)VH+σ2·I)-1=
P·VH(PVVH+σ2·I)-1=
VH(VVH+·I)-1。
(19)
當(dāng)發(fā)送信號x的能量P為1時,則可得:
W=VH(VVH+σ2·I)-1。
(20)
綜合上述分析,MIMO模型中利用信道數(shù)據(jù)計算信道的MMSE均衡器接收矩陣的復(fù)雜度主要來源于SVD分解與式(20)中的矩陣求逆。
由前文可知,H、V和W之間的關(guān)系可以由圖4表示。
圖4 H、V和W的關(guān)系示意圖Fig.4 Relationship of H,V and W
利用相關(guān)矩陣組的關(guān)聯(lián)性以降低計算復(fù)雜度,其具體分析及操作如下。
3.1.1 矩陣數(shù)據(jù)的關(guān)聯(lián)性
對于信道系數(shù)復(fù)數(shù)矩陣H={Hj,k},其中Hj,k∈M×N,j=1,2,…,J;k=1,2,…,K。因此,該矩陣是一個M×N×J×K維的信道矩陣,其中M表示接收天線的數(shù)量,N表示發(fā)射天線的數(shù)量,J表示用戶數(shù)量,K表示時隙個數(shù)。
由前文建立的通信模型,對H矩陣內(nèi)在的關(guān)聯(lián)性進(jìn)行分析。首先,Hj,k內(nèi)部的關(guān)系可表示為不同天線構(gòu)建出的不同信道之間的相關(guān)性。在一般高斯白噪聲信道條件下,天線陣列之間固定的距離和入射角關(guān)系將帶來一定的規(guī)律,但數(shù)據(jù)集中未能發(fā)現(xiàn)Hj,k內(nèi)部可靠的相關(guān)特性。這可能是由于天線之間的方向性、距離之間的差異較大,使得信道在空間上的相關(guān)特性不再明顯。
考慮不同k時隙,同一用戶j的信道系數(shù)情況,即MN個信道的時間相關(guān)性。圖5給出了MN個信道在時隙k=1,2,3情況下的幅度響應(yīng)和相位響應(yīng)。可以看出,在不同k下的幅度和角度的變化程度不大,這可以理解為在相關(guān)時間內(nèi),信道的變化很小,這進(jìn)一步驗證通信建模的可行性。
(a) 不同信道系數(shù)的幅度響應(yīng)
(b) 不同信道系數(shù)的相位響應(yīng)圖5 同j不同k的H塊之間的幅度和相位響應(yīng)關(guān)系Fig.5 Amplitude and phase response relationship between blocks H with the same j and different k
基于上述兩層分析,得出Hj,k塊在時間維上的相關(guān)性。而在時間維上利用SVD與MMSE求W矩陣是獨立的,無法直接用于算法簡化,為此,本文利用時間相關(guān)性,并運用插值算法直接估計W矩陣。
圖6 W矩陣的插值示意圖Fig.6 Schematic diagram of interpolation of W matrix
3.1.2 矩陣數(shù)據(jù)W的插值算法
對于矩陣數(shù)據(jù)W的插值算法,采用linear插值法、spline插值法與Pchip插值法進(jìn)行建模插值[2-4]。此處,引入“最大插值比”作為評估方法來評價插值性能,參數(shù)尋優(yōu)的過程可以表示為:
(21)
式中,R表示每隔R個點進(jìn)行一次插值運算。因此,式(21)表示插值結(jié)果滿足最小建模精度的約束下,使得插值數(shù)量越多的尋優(yōu)目標(biāo)。
因此,為滿足題設(shè)對于擬合后W矩陣對最小建模精度的要求,需選擇合適的插值方法并計算最大插值比,為此進(jìn)一步提出了最大插值比搜索方法,以評估在不同信道條件數(shù)據(jù)集下可用的插值參數(shù),最大插值比搜索算法的詳細(xì)過程在算法1中給出,其基本思想為通過給定的初始插值比Rate,和選定的插值類型,不斷迭代和逼近給定插值類型下的滿足要求的最大插值比。當(dāng)取得的插值比越大時,意味著W矩陣的更多部分可以通過在K維度上的相關(guān)性插值得到,不需要通過SVD與MMSE求逆過程,這將大大減少計算的復(fù)雜度。此外,Data集最大插值比的計算過程可以理解為適應(yīng)信道的訓(xùn)練過程。在后續(xù)過程中,在外部信道條件未劇烈改變的情況下,不需要再次執(zhí)行,因此最大插值比搜索可以在線下執(zhí)行,不會影響算法的復(fù)雜度。
算法1 所提出的最大插值比搜索算法輸入:訓(xùn)練數(shù)據(jù)集Data,包含通過MMSE計算的標(biāo)準(zhǔn)W矩陣;初始化:初始插值比 Rate=1/R,插值間隔R的初始值可以設(shè)定為R=2。R的中止值可以設(shè)定為Rmax=200選定的插值類型:Linear,Spline,Pchip。執(zhí)行過程:ForR
表2給出了3種常見插值方法在6個Data集的最大插值比。
由于Data1~Data6來自不同的信道條件,因此同一插值方法在插值過程中計算出的最大插值比有較大區(qū)別。例如Data 3數(shù)據(jù)集的最大插值比顯著較小,這可以理解為信道的時變性強或受到干擾噪聲較大,插值算法在此時難以滿足要求,需要降低最大插值比。
進(jìn)一步,橫向?qū)Ρ?種插值方法,可見在多數(shù)的數(shù)據(jù)集下,Linear插值的最大插值比最小,性能最差,這是由于簡單的Linear插值精度較低。Spline插值與Pchip插值的最大插值比性能相近,Spline插值在Data 1與Data 3上表現(xiàn)比Pchip插值較好。從原理上分析,可以理解為當(dāng)基礎(chǔ)函數(shù)振蕩時,Spline比Pchip更好地捕獲點之間的移動,后者會在局部極值附近急劇扁平化,這在該場景下帶來了更好的插值性能[5-6]。
表2 3種插值方法在6個Data集的最大插值比
復(fù)數(shù)矩陣運算可拆解為基本的復(fù)數(shù)運算,而基本的復(fù)數(shù)運算又可進(jìn)一步拆解為基本的實數(shù)運算,實數(shù)基本運算的復(fù)雜度按照表1計算。
對于Linear插值法,根據(jù)上述描述,可得需要加減法6次、乘法2次、倒數(shù)1次,且對于復(fù)數(shù),幅度和相位要分別插值,總復(fù)雜度是實數(shù)插值的兩倍。然而,特殊的是這里插值點的橫坐標(biāo)是均勻的,因此計算的復(fù)雜度可以大大簡化,僅需要3次加法、2次乘法和1次倒數(shù),因此總復(fù)雜度為68。每個插值時刻k,共需要NJL次插值,因此需要復(fù)雜度68NJL。本文的參數(shù)取值為N=64,J=4,L=2,因此計算復(fù)雜度為17 408。
對于Spline插值法,其計算步驟為:求三次函數(shù)的系數(shù),然后將插值點橫坐標(biāo)代入三次函數(shù),計算又需要6次加法、6次乘法。且對于復(fù)數(shù),幅度和相位要分別插值,總復(fù)雜度是實數(shù)插值的2倍,因此總復(fù)雜度為124,每個插值時刻k,計算復(fù)雜度124NJL。取本文參數(shù),計算復(fù)雜度為31 744。
對于Pchip插值,由于其性能不如Spline且計算復(fù)雜度相似,因此不在此處進(jìn)行考慮。
首先,對于矩陣分塊直接求逆,假設(shè)一個N階(這里的N被重新定義)的方陣A,分塊如下:
(22)
設(shè)A的逆矩陣分塊如下:
(23)
則根據(jù)矩陣分塊求逆的原理有:
(24)
式中,對于N×N階矩陣,需要的矩陣相乘的運算量級為N3。
對于n>1,m>1,直接分塊求逆算法在具體實現(xiàn)中需要利用遞歸實現(xiàn),具體的運算量按照復(fù)數(shù)乘加次數(shù)統(tǒng)計,對于上述的直接求逆算法,以乘加次數(shù)統(tǒng)計理論運算量,式(24)各部分運算量:c11的運算量為4m2n+4mn2-mn,c12和c21的運算量相等,均為4m2n+4mn2-6mn,同理,c22的運算量為4m2n+4mn2-m2,為此,矩陣A利用一次分塊求逆的總的運算量為:
T(1)(N)=T(m)+T(n)+16m2n+
12mn2-13mn-m2+4m3,
(25)
式中,T(1)(N)表示利用一次矩陣分塊求逆算法計算矩陣求逆的總計算量,T(m)和T(n)分別表示對m和n階復(fù)矩陣求逆所需的運算量。
采用改進(jìn)的Strassen矩陣求逆算法,結(jié)合式(22)~(23),Strassen算法應(yīng)用到求逆運算有如下公式:
(26)
按照前面相同的運算量統(tǒng)計方法,根據(jù)式(25),矩陣?yán)靡淮畏謮K求逆的總運算量為:
T(1)=T(m)+T(N)+13m2n+
11mn2-4mn+3m2。
(27)
對于n>1,m>1,Strassen矩陣求逆算法也是利用遞歸實現(xiàn)的,但因為Strassen算法減少了矩陣復(fù)乘次數(shù),所以相比直接分塊的常規(guī)算法運算量有明顯降低。
(28)
式中,對于N×N階矩陣,需要的矩陣相乘的運算量級為N2。
為了便于比較新求逆算法的性能改善,按照前面相同的運算量統(tǒng)計方法,根據(jù)式(28),矩陣A利用一次分塊求逆的總的運算量為:
T(1)(N)=T(m)+T(N)+8m2n+8mn2+mn,
(29)
式中,T(1)(N)表示利用改進(jìn)算法計算一次矩陣求逆的運算量,T(m)和T(n)部分表示對m和n階復(fù)矩陣求逆所需的運算量。
和矩陣直接分塊求逆算法相比,新的求逆算法雖然增加了一些加減運算,但復(fù)乘次數(shù)降低。對于維數(shù)較高的矩陣,其中有大量的復(fù)矩陣運算,復(fù)乘消耗的運算量將遠(yuǎn)大于加減法,而這個運算量隨著矩陣維數(shù)增加將有顯著增大,因此新算法對于復(fù)乘次數(shù)的降低將顯著改善求逆運算的實時性能。和常規(guī)Strassen矩陣求逆算法相比,改進(jìn)的算法由于利用了求逆矩陣的特點,即對Hermite矩陣進(jìn)行求逆運算,所以在運算量和算法復(fù)雜度上都有明顯的降低。
常規(guī)算法計算一次矩陣求逆的計算復(fù)雜度如表3所示,改進(jìn)算法計算一次矩陣求逆的計算復(fù)雜度如表4所示。
表3 常規(guī)算法計算一次矩陣求逆的計算復(fù)雜度
表4 改進(jìn)算法計算一次矩陣求逆的計算復(fù)雜度
在本文中,當(dāng)矩陣維度為8(J=4)時,改進(jìn)算法總共復(fù)雜度為4 776,常規(guī)算法總共5 255,復(fù)雜度降低10.03%。
進(jìn)一步,圖7給出了不同用戶數(shù)量J時,改進(jìn)的Strassen算法與常規(guī)Strassen算法的復(fù)雜度比較??梢钥闯?,隨著用戶數(shù)量增加,矩陣求逆時的維度增加,采用改進(jìn)的Strassen算法對復(fù)雜度的降低更加明顯。
圖7 在不同用戶個數(shù)J下,改進(jìn)的Strassen算法與 常規(guī)Strassen算法的復(fù)雜度比較Fig.7 Under different number of users J,the complexity comparison of the improved Strassen algorithm and the conventional Strassen algorithm
圖8 采用改進(jìn)的矩陣求逆算法后對各數(shù)據(jù)集 最小建模精度的影響Fig.8 Influence of the improved matrix inversion algorithm on the minimum modeling accuracy of each data set
對于不采用插值算法的情況,每個插值時刻k,由H矩陣到W矩陣,需要進(jìn)行J(J=4)次SVD和1次MMSE均衡(主要復(fù)雜度在于求逆)的計算。其中4次SVD需要3 574 400次計算。MMSE均衡需要521 112次計算,因此共需要計算復(fù)雜度C1=4 095 512。采用Spline插值時,每個插值時刻的復(fù)雜度為31 744??梢杂嬎愠雒看尾逯档膹?fù)雜度收益為C2=3574400+521112-31744 = 4 063 768。因此采用插值的最終復(fù)雜度收益為:
ΔC=(C1-C2)×K×Rate。
(30)
假設(shè)Rate=1/3時,ΔC=524 225 536。可見,插值對于計算復(fù)雜度的降低比較明顯。同樣,計算復(fù)雜度可以降低為:
(31)
當(dāng)Rate=1/3時,計算復(fù)雜度降低為原來的66.93%。當(dāng)Rate=1/10時,計算復(fù)雜度降低為原來的90.08%。
綜上,當(dāng)采用改進(jìn)的Strassen矩陣求逆算法時,復(fù)雜度降低了10.03%。進(jìn)一步采用插值算法后,計算復(fù)雜度能夠在上述的基礎(chǔ)上再降低9.92%(插值比為1/10)、33.07%(插值比為1/3)。
本論文主要解決MIMO場景下的相關(guān)矩陣組的低復(fù)雜度計算問題,首先利用H矩陣在時間相關(guān)性推導(dǎo)了W矩陣的相關(guān)性,通過對已有W矩陣的相關(guān)性直接插值出部分缺失的W。這使得在接收H矩陣時,在求取部分W矩陣后通過相關(guān)性重建完整的W矩陣;也避免了一部分H矩陣的存儲以及這部分H矩陣計算SVD與求逆獲得W矩陣的過程。相比SVD與求逆的復(fù)雜度,插值的復(fù)雜度要低得多。此外,采用了一種改進(jìn)的Strassen矩陣求逆算法來降低求逆過程的復(fù)雜度。該算法結(jié)合了Strassen矩陣求逆的高效性以及Hermite正定陣的共軛對稱性特點,結(jié)構(gòu)更簡化。