[發明專利]一種面向硬件實現的快速投影原子選擇正交匹配追蹤重構算法

在審

• 申請號: 202110451299.2
• 申請日: 2021-04-26
• 公開（公告）號: CN113300713A
• 公開（公告）日: 2021-08-24
• 發明（設計）人: 劉素娟;鄭麗麗
• 申請（專利權）人: 北京工業大學
• 主分類號: H03M7/30
• 分類號: H03M7/30
• 代理公司: 北京思海天達知識產權代理有限公司 11203
• 代理人: 劉萍

1.一種面向硬件實現的快速投影原子選擇正交匹配追蹤重構算法，其特征在于：

表1為變量符號的含義說明；

表1變量符號含義說明

具體步驟如下：

A.主函數Fast-POMP

步驟一、輸入輸出數據

輸入:y,A,K,S

輸出:t,I

輸入數據中:y是測量值向量；A是測量矩陣；K是信號稀疏度，表明信號需要稀疏先驗；S是投影原子選擇階段索引并行選擇的個數；

輸出數據中：是重構信號；t是Fast-POMP算法重構信號所用的迭代次數；I是最終支撐集，即所有被選擇的索引的集合；

步驟二、數據初始化

t＝1,L＝K,r⁰＝y,[～,N]＝size(A)

迭代次數t初始化為1；L是索引初步選擇的個數，初始化為K與稀疏度大小相等；r表示殘差，其意義在于表征測量值向量y與信號估計值的貢獻分量之間的差值，r⁰表示殘差初始值，這里初始化為測量值向量y；最終支撐集I初始化為空集size(A)表示取矩陣A的行數和列數，符號“～”表示不輸出，N是測量矩陣A的列數也是原始信號的長度(點數)，即信號長度需要先驗；子函數Incremental_IMGS_QRD分解得到的變量數據包括Qid1，Rid1，Did1，Qid2，Rid2，Did2，均初始化為空集

步驟三、Fast-POMP算法主體，分為以下33Steps：

Step 1-26是Fast-POMP算法中用來確定最終支撐集I的過程；

1.for i＝1:K

2.Qid1＝Qid2,Rid1＝Rid2,Did1＝Did2；

算法在信號完全重構的時候，迭代次數遠小于稀疏度K次，而在不能完全重構的時候，Step 1中主循環體的設置保證迭代次數不超過K次；Step 2用于實現中間數據傳輸，將后面Step 19中得到的Qid2,Rid2,Did2的值分別賦給Qid1,Rid1,Did1；該數據傳輸過程是基于分解基礎的增量矩陣分解結構的重要橋梁；

3.Ip＝argmax_L(|A^T*r^t-1|)；

Step 3實現索引的初步選擇；函數argmax_L()表示取輸入向量的前L個最大的原子所對應的索引；||表示對輸入數據取絕對值，輸入數據是單一變量、向量、或者矩陣；上標T表示對矩陣進行轉置，A^T代表矩陣A的轉置，r^t-1表示上一次迭代后得到的殘差；先計算測量矩陣A和殘差r^t-1的內積，內積向量取絕對值后，取前L個絕對值最大的原子所對應的索引放入集合Ip中；初步選擇的索引集合Ip稱之為潛在支撐集；

4.Iu＝I^t-1∪Ip；

Step 4實現候選支撐集擴充；上標t-1表示上一次迭代過程；I^t-1是最終支撐集I的子集，在I逐漸擴充的過程中，為方便起見，將I^t-1稱之為最終支撐集；Iu是擴充后的集合，將潛在支撐集Ip放入最終支撐集I^t-1得到需要后續處理的集合Iu，稱之為候選支撐集；

5.for j＝1:L

6.[Q,R,D]＝Incremental_IMGS_QRD(A(:,Ip_(j)),Qid1,Rid1,Did1,(i-1)*S+j)；

7.Qid1＝Q,Rid1＝R,Did1＝D；

8.endfor

Step 5-8第一個基于已分解矩陣的增量矩陣分解；Step 5-8是一個小循環，Step 6中調用了Incremental_IMGS_QRD函數完成矩陣分解，Step 7實現矩陣分解后的數據存儲；Incremental_IMGS_QRD函數使用了IMGS-QRD算法并在此基礎上實現增量分解矩陣的功能，B部分將會詳細介紹該函數；輸入數據中，Ip_(j)代表潛在支撐集Ip的第j個元素；A(:,Ip_(j))代表潛在支撐集Ip的第j個元素所對應測量矩陣A中的列向量；Qid1，Rid1，Did1是之前迭代過程分解好的矩陣或者空矩陣；(i-1)*S+j表示本次是第幾次分解；Incremental_IMGS_QRD函數只分解了潛在支撐集Ip所對應矩陣A中的列，而沒有分解候選支撐集Iu所對應測量矩陣A中的列；

9.v＝Q^T*y；

10.m1＝(i-1)*S+L,xu＝zeros(m1,1)；

11.for a＝1:L

12.xu_(m1-a+1)＝(v_(m1-a+1,:)-R_(m1-a+1,：)*xu)；

13.enfor

Step 9-13是部分索引投影的過程；Step 9將分解得到的矩陣Q的轉置矩陣與測量值向量y相乘，相乘結果記為列向量v；Step 10創建一個尺寸和候選支撐集Iu相同的零向量xu，其作用是存放執行后續步驟后得到的投影值；m1代表當前迭代次數下候選支撐集Iu中索引的個數，其具體數值為上一次迭代次數i-1與每次迭代過程中并行原子選擇個數s的乘積再與當前迭代次數下初步原子選擇個數L的和；Step12實現投影計算，v_(m1-a+1,:)代表矩陣v的第m1-a+1行；R_{(m1-a+1,；)}代表矩陣R的第m1-a+1行；使用v_(m1-a+1,:)減去R_(m1-a+1,:)與已得到的估計值xu的乘積得到當前估計值xu_(m1-a+1)；

14.Is＝argmax_S(|xu|)；

Step 14是投影后的索引選擇；函數argmax_S()表示取輸入數據前S個最大元素所對應的索引；所以Step 14意思是將xu絕對值的前S個最大的原子所對應的索引存儲在集合Is中；該過程采用每輪迭代中多個索引并行選擇的方法，極大地減少了總體迭代次數，提高了計算效率；由于迭代次數的大幅度減少，算法整體計算復雜度也會大幅度地降低；

15.I^t＝I^t-1∪Ip_(Is)；

Step 15是支撐集擴充；Ip_(Is)表示集合Is中的索引在Ip中相應位置的元素所組成的集合；取集合Ip_(Is)與上一次迭代得到的最終支撐集I^t-1的并集得到更新后的最終支撐集I^t；

16.for j＝1:S

17.z＝(i-1)*S+j；

18.[Q1,R1,D1]＝Incremental_PMGS_QRD(A(:,Is_(j)),Qid2,Rid2,Did2,z)；

19.Qid2＝Q1,Rid2＝R1,Did1＝D1；

20.r^z＝r^z-1-Q1_(:,z)*(Q1_(:,z))^T*r^z-1*D1_(z)；

21.endfor

Step 16-21實現第二個基于已分解矩陣的增量矩陣分解并完成殘差的計算；Step 18和Step 6的功能相同，通過調用Incremental_PMGS_QRD函數實現矩陣分解，通過Step 2傳遞數據，實現它們在相同的已分解矩陣的基礎上進行繼續增量分解的目的；該過程在分解得到的矩陣Qid2，Rid2，Did2的基礎上，只需要S次增量；Step 19實現數據的存儲；Step20是殘差的計算；通過Q1矩陣的第z列與其轉置及上次的殘差r^z-1，向量D的第z個元素的乘積得到關于殘差的貢獻分量，使用上次迭代產生的殘差r^z-1減去該貢獻分量得到本次迭代的殘差r^z；

22.if||r^t*S||₂1e-6

23.break；

24.endif

25.t＝t+1；

26.endfor

Step 22-26是判斷算法停止的條件；假設x是一個向量，||x||₂表示向量x的l₂范數，具體公式為如果殘差的l₂范數大于10^-6，說明誤差尚未滿足要求，算法繼續運行，迭代次數t＝t+1，新的迭代過程開始；如果殘差的l₂范數小于10^-6，說明誤差已經減小到忽略不計的程度，則跳出大循環，算法停止，迭代次數t不加1；

27.v1＝Q1^T*y；

28.m2＝i*S,xl^t＝zeros(m2,1)；

29.for b＝1:m2

30.xl^t_(m2-b+1)＝(v1_(m2-b+1)-R1_(m2-b+1,:)*xl^t)；

31.endfor

Step 27-31是信號頻譜上非零值的估計過程；Step 27中Q1的轉置矩陣乘以測量值向量y，將結果記為v1；Step 28創建一個和最終支撐集I同樣尺寸的零向量x1，其作用是存放后續計算得到的非零元素估計值，此時最終支撐集I所包含的索引個數為迭代次數與每次迭代過程中并行選取的索引個數S的乘積，將結果記為m2；Step 29-31是計算信號非零元素的過程；Step 30和Step 12計算過程相同，xl^t_(m2-b+1)表示本次(t次)迭代過程中向量xl的第m2-b+1個元素，v1_(m2-b+1)表示v1中的第(m2-b+1)個元素，R1_(m2-b+1,:)表示R1中的第(m2-b+1)行；v1_(m2-b+1)減去R1_(m2-b+1,:)與已估計值的xl^t乘積得到當前估計值xl^t_(m2-b+1)；

32.

33.t,I；

Step 32非零值歸位；Step 32中，首先定義了一個長度為N的零向量然后將估計信號x1的非零值放置到索引集I所對應的的位置上，得到最終的估計信號Step 33列出了要輸出的數據，包括信號估計值迭代次數t和最終索引集I；

B.子函數Incremental_IMGS_QRD

步驟一、輸入輸出數據

輸入:w，Q^i-1，R^i-1，D^i-1，i

輸出:Qⁱ，Rⁱ，Dⁱ

輸入數據中，w代表一個索引對應的測量矩陣A中的一列；Q^i-1、R^i-1和D^i-1表示經過上一次迭代后分解好的矩陣；i表示第i次增量分解；

輸出數據中，Qⁱ、Rⁱ和Dⁱ表示本次分解好的矩陣；

步驟二、第一次分解

1.ifi＝＝1；

2.Dⁱ＝w^T*w；

3.Qⁱ＝w./Dⁱ；

4.Rⁱ＝[1]；

Step 1-4是第一次分解；Dⁱ等于輸入列向量w的內積，該內積通過列向量w的轉置與w相乘得到；Qⁱ通過輸入列向量w的每個元素除以Dⁱ得到，此時Dⁱ只有一個元素；Rⁱ為1；

步驟三、增量分解

5.else

6.Dⁱ＝[D^i-1 0]；

7.

8.for j＝1:i-1

9.

10.

11.endfor

12.

13.

14.end

Step 5-14完成了基于已分解矩陣的增量分解操作；其中Step 6-7是矩陣Dⁱ和矩陣Rⁱ的初始化，目的是為后續增量分解操作預留存儲空間；Step 6中D向量拓展了一個元素長度的空間，并賦值為0；Step 7中，矩陣R的右上角通過(i-1)長度的列零向量來初始化，左下角通過(i-1)長度的行零向量來初始化，右下角的對角線上的單個元素初始化為1；循環體中的Step 9完成上三角矩陣Rⁱ的非對角元素Rⁱ_(j,i)的賦值，表示矩陣Q^i-1的第j列的轉置，非對角元素Rⁱ_(j,i)通過與輸入列向量的w乘積得到；Step 10中，輸入的列w被更新；其中，表示向量D^i-1的第j個元素，輸入列w通過自身減去Rⁱ_(j,i)、與的乘積實現更新；Step 12中，Dⁱ的第i個元素Dⁱ_(i)被計算出，其中表示l₂范數的平方，計算公式為Step 13中，矩陣Qⁱ被計算出，Qⁱ通過Q^i-1并上新增的列產生，每次擴充增加的列通過輸入的列向量w與Dⁱ_(i)相除得到。