CN106940689A - 基于Jacobi迭代算法的高精度矩阵特征值分解实现方法 - Google Patents

基于Jacobi迭代算法的高精度矩阵特征值分解实现方法 Download PDF

Info

Publication number
CN106940689A
CN106940689A CN201710130519.5A CN201710130519A CN106940689A CN 106940689 A CN106940689 A CN 106940689A CN 201710130519 A CN201710130519 A CN 201710130519A CN 106940689 A CN106940689 A CN 106940689A
Authority
CN
China
Prior art keywords
precision
eigenvalue decomposition
matrix
iterative algorithms
transferred
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201710130519.5A
Other languages
English (en)
Inventor
廖红舒
赵文扬
龙慧敏
甘露
梁先明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201710130519.5A priority Critical patent/CN106940689A/zh
Publication of CN106940689A publication Critical patent/CN106940689A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)

Abstract

本发明属于信号处理领域，尤其涉及种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法。发明提供一种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法，在不明显增加算法复杂度、FPGA实现难度与增加资源消耗的情况下，提高实际工程中基于循环Jacobi迭代算法的FPGA实现矩阵特征值分解的计算精度。

Description

基于Jacobi迭代算法的高精度矩阵特征值分解实现方法
技术领域
本发明属于信号处理领域，尤其涉及种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法。
背景技术
在信号处理中，矩阵的特征值分解EVD是一个应用广泛的矩阵运算。如数据压缩、噪声去除、数值分析，包括近几年兴起的机器学习、深度学习其基本核心操作也包括矩阵特征值分解。实现矩阵特征值分解的常用方法有Gauss变换、Householder变换、Jacobi迭代等，其中，Jacobi迭代是精度较高的方法，并且很适合在FPGA中实现。因此一种基于Jacobi迭代算法的高精度矩阵特征值分解实现技术在实际工程中具有很高的应用价值。
经典的Jacobi迭代算法计算共轭矩阵A∈Cn×n的特征值分解如图1所示，这种经典的迭代算法虽然有较快收敛速度，但是该算法需要在矩阵A的众多元素中选取aij，使得aij为非对角元素中绝对值最大的一个，再进行后面的计算操作。这样每一步都要寻找绝对值最大的非对角元，比较费时也不适合在FPGA实现，因此经典的Jacobi迭代算法在实际工程中并不实用。
目前实际工程中多数采用如图2所示的循环Jacobi迭代算法，通过逐行扫描遍历法选取aij，这样避免了寻找最大绝对值的非对角元的复杂繁琐步骤。这样选取aij的方式，在aij数值比较大时，FPGA中使用Cordic算法计算φ、θ误差比较小，可以取得比较好的效果。但当aij比较小甚至接近0时，此时FPGA中使用Cordic算法计算φ、θ误差比较大，就会导致后面计算A＝QHAQ产生误差，其中Q∈Cn×n为复数域内的平面旋转矩阵。而计算过程需要多次的迭代运算，如果在迭代过程中出现多次aij比较小甚至接近0的情况，就会产生较大的累计误差，从而使最终计算结果的精度相对较差。
发明内容
发明的目的在于解决在循环Jacobi迭代算法进行矩阵特征值分解过程中，因为逐行扫描遍历法选取aij在aij比较小甚至接近于0时，导致FPGA中使用Cordic算法计算φ、θ误差比较大，进而使迭代过程产生相对较大的累计误差，导致计算结果误差增大。即提供一种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法，在不明显增加算法复杂度、FPGA实现难度与增加资源消耗的情况下，提高实际工程中基于循环Jacobi迭代算法的FPGA实现矩阵特征值分解的计算精度。
一种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法，包括如下步骤：
S1、设数据矩阵A∈Cn×n为共轭矩阵，并且设定最大遍历次数为T、最小清扫门限a、扩位门限b和算术左移位数m，其中，最小清扫门限a应小于要求计算结果的精度一个数量级，扩位门限b与算术左移位数m与FPGA实现里数据位宽size有关，满足b×2m＜2size-4保证计算结果不溢出，n为不为零的自然数，共轭矩阵A中的元素为aij，i＝1,2,3,...,n，j＝1,2,3,...,n，1≤t≤T且t为自然数；
S2、初始化遍历次数计数器，令t＝0，
初始化特征向量初始矩阵，令V＝E，其中，E为单位阵；
S3、在S1所述共轭矩阵A中选取aij，初始化清扫元aij行列下标，令i＝1，j＝2；
S4、判断aij是否满足跳过清扫条件|real(aij)|＜a&|imag(aij)|＜a，若满足则转入S10，如不满足则转入S4；
S5、判断aij是否满足扩展条件|real(aij)|＜b&|imag(aij)|＜b，若满足则转入S6，若不满足则转入S7；
S6、进行位扩展，即计算a′ij＝aij×2m，转入S7；
S7、令a′ij＝aij，进入S8；
S8、计算根据所得计算
S9、计算A＝QHAQ和V＝QHV，其中，Q∈Cn×n为复数域内的平面旋转矩阵，
即Q的对角元素中除了qii＝ejφcosθ、qjj＝e-jφcosθ之外其他均为1，非对角元素中除了qij＝-ejφsinθ、qji＝e-jφsinθ其他元素均为0，θ为旋转角度；
S10、判断j＝n是否成立，是则进入S11，否则j＝j+1后跳转到S4；
S11、判断i＝n-1是否成立，是则进入S12，否则i＝i+1，j＝i+1后跳转到S4；
S12、判断t＝T是否成立，是则进入S13，否则t＝t+1后跳转到S3；
S13、输出迭代计算结果A和V，其中，A的对角元数位S1输入数据矩阵A的特征值，V为对应的特征向量矩阵。
进一步地，S1所述遍历次数T越大迭代次数越多，计算越精确，但计算时间越长，为了取得速度与精度的平衡，当n≤8时T＝3，当n＞8时T＝6。
进一步地，S1所述a＝e×10-1。
本发明的有益效果是：
在不明显增加算法复杂度、FPGA实现难度与硬件资源消耗、计算消耗时间的情况下，提高基于循环Jacobi迭代算法的FPGA实现矩阵特征值分解的计算精度和计算速度，在实际工程中具有重要价值。
附图说明
图1为经典Jacobi迭代算法流程。
图2为循环Jacobi迭代算法流程。
图3为本发明算法流程。
具体实施方式
下面将结合实施例和附图，对本发明方法进行进一步说明。
本发明应用于估计信号与噪声对应的盖氏圆盘半径，提高圆盘半径计算精度，和计算速度。
实施例1、
接收阵列为8阵元组成的均匀线阵。
如图3所示，考虑N＝1个载波频率为的BPSK调制的远场信号s(k)，以γ＝0°的方向入射到阵元数n＝8的均匀线阵上，且阵元间距为d＝0.5λ，其中，λ为信号波长，阵列接收噪声是功率为σ2＝1的零均值高斯白噪声，接收信号信噪比SNR＝20dB，快拍数为L＝1024。通过进行特征值分解估计信号与噪声对应的盖氏圆盘半径i＝1,2,…,n-1。
在估计性能包括计算精度、计算速度和资源消耗，具体用下面指标评价：
1.计算精度：
(1).圆盘半径计算精度：i＝1,2,…,n-1其中κi为圆盘半径的理论值。εi越小表示计算精度越高。
(2).圆盘平均计算精度： 越小表示平均计算精度越高。
2.计算速度：
(1).计算消耗的时钟数Nclk，越小表示计算消耗时间越少，计算速度越快。
3.资源消耗：
(1).寄存器消耗数量Nreg，越小对应寄存器资源消耗越少。
(2).逻辑门消耗数量Nlut，越小对应逻辑门资源消耗越少。
应用特征值分解估计信号与噪声对应的盖氏圆盘半径包括，a.仿真接收信号数据建模、b.应用本发明进行特征值分解、c.计算圆盘半径，具体为以下步骤：
a.仿真接收信号数据建模。
a1.由下式产生阵元数n＝8的阵列接收信号向量X(k)＝[x1(k) x2(k) … x8(k)]H进入步骤a2。
X(k)＝a(γ)S(k)+N(k)，k＝1,2,…,L
式中，N(k)为8×1为均值为零、方差σ2＝1的高斯白噪声矢量；远场接收信号S(k)＝Ass(k)，其中其幅度As＝10SNR/20；a(γ)＝[1 e-jφ … e-j(n-1)φ]T，为空间阵列的n×1维流型矩阵。
a2.由计算数据协方差矩阵R∈Cn×n，进入步骤a3。
a3.根据对数据协方差矩阵进行分块，得到块矩阵R′∈C(n-1)×(n-1)，进入步骤b，其中，R′∈C(n-1)×(n-1)，r∈C(n-1)×1，rnn为数据矩阵R第n行第n列的元素。
b.应用本发明对块矩阵R′进行特征值分解，计算R′的特征值矩阵D与特征向量矩阵V。
b1.进行初始化，具体方法为：
b11.设数据矩阵A＝R′为共轭矩阵，并且设定遍历次数T＝5、最小清扫门限a＝10-8、扩位门限b＝10-5和算术左移位数m＝8，进入步骤b12。
其中T越大迭代次数越多计算越精确但计算时间越长，根据矩阵维数n进行选取，当n≤8时T＝4，n＞8时T＝6可以取得速度与精度的平衡；最小清扫门限a与计算精度要求e有关，a＝e×10-1，比如要求计算精度为10-5则a≈10-6。扩位门限b与算术左移位数m与FPGA实现里数据位宽size有关，满足b×2m＜2size-4保证计算结果不溢出。
b12.初始化遍历次数计数器和特征向量初始矩阵，t＝0、V＝E，其中，E为单位阵，进入步骤b13。
b13.初始化清扫元aij的行列下标，i＝1、j＝2，进入步骤b2。
b2.进行Jaocbi旋转，具体方法如下：
b21.在矩阵A中选取aij，进入步骤b22。
b22.判断是否满足跳过清扫条件|real(aij)|＜a&|imag(aij)|＜a，是则跳转到步骤b3，否则进入步骤b23。
b23.判断是否满足扩展条件|real(aij)|＜b&|imag(aij)|＜b，是则进入步骤b24，否则进入步骤b25。
b24.进行位扩展，即a′ij＝aij×2m，进入步骤b26。
b25.不进行位扩展，即a′ij＝aij，进入步骤b26。
b26.计算相角和模值，即和进入步骤b27。
b27.计算平面旋转角，即进入步骤b28。
b28.进行Jacobi旋转，即计算A＝QHAQ和V＝QHV，其中Q∈Cn×n为复数域内的平面旋转矩阵。
即Q的对角元素中除了qii＝ejφcosθ、qjj＝e-jφcosθ之外其他均为1，非对角元素中除了qij＝-ejφsinθ、qji＝e-jφsinθ其他元素均为0，进入步骤b3。
b3.对迭代过程进行判断。
b31.判断j＝n是否成立，是则进入步骤b32，否则j＝j+1后跳转到步骤b21。
b32.判断i＝n-1是否成立，是则进入步骤b33，否则i＝i+1，j＝i+1后跳转到步骤b21。
b33.判断t＝T是否成立，是则进入步骤b4，否则t＝t+1后跳转到步骤b13。
b4.输出迭代计算结果A和V，其中A的对角元就是数据矩阵A的特征值，V为对应的特征向量矩阵，进入步骤c。
c.对数据协方差矩阵R进行酉变换，计算信号与噪声对应的圆盘半径。
c1.由下式构造酉变换矩阵T∈Cn×n，进入步骤c2。
其中，V∈C(n-1)×(n-1)为前面计算块矩阵R′的特征向量，满足VVH＝E，E为单位阵。
c2.进行酉变换得到圆盘半径，即计算下式，进入步骤c3。
式中，λi，i＝1,2,…,n-1为块矩阵R′的特征值。
c3.由ri＝|ρi|，i＝1,2,…,n-1计算圆盘半径ri，进入步骤c4。
c4.计算圆盘半径计算精度：i＝1,2,…,n-1其中κi为圆盘半径的理论值，和圆盘平均计算精度：进入步骤c5。
c5.统计计算消耗的时钟数Nclk、寄存器消耗数量Nreg和逻辑门消耗数量Nlut，算法结束。
仿真结果为：
计算精度：
计算速度:Nclk＝11710
资源消耗:Nreg＝29104、Nlut＝30254
此时，估计信号所对应的圆盘半径计算精度为ε1≈10-9；估计噪声所对应的圆盘半径计算精度为εi≈10-4，i＝2,3,…,7；圆盘平均计算精度
实施例2、
经典方法循环Jacobi算法应用于估计信号与噪声对应的盖氏圆盘半径的估计性能，作为实施例1的对比例。
实施例2的方法如附图2所示，其余仿真条件与实施例1的相同，进行信号与噪声对应的盖氏圆盘半径的估计。
实施例2的评价标准与实施例1的一致。
仿真结果为：
计算精度：
计算速度:N′clk＝17960
资源消耗:N′reg＝29101、N′lut＝29998
本发明此时，估计信号所对应的圆盘半径计算精度为ε1≈10-8；估计噪声所对应的圆盘半径计算精度为εi≈10-1，i＝2,3,…,7；圆盘平均计算精度
综上所述，对比实施例1与实施例2，本发明相对于经典方法在增加(Nreg-N′reg)/N′reg×％≈0.01％寄存器资源消耗，(Nlut-N′lut)/N′lut×％≈0.85％查找表资源消耗的情况下，平均计算精度从10-1提高到了10-4数量级，提高了3个数量级，同时计算速度提高了|Nclk-N′clk|/N′clk×％≈34.8％。
所以，本发明在基本不增加资源消耗的情况下，不仅可以提高计算精度，还可以提高计算速度，在实际工程中具有重要价值。

Claims (3)

1.一种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法，其特征在于，包括如下步骤：
S1、设数据矩阵A∈Cn×n为共轭矩阵，并且设定最大遍历次数为T、最小清扫门限a、扩位门限b和算术左移位数m，其中，最小清扫门限a应小于要求计算结果的精度一个数量级，扩位门限b与算术左移位数m与FPGA实现里数据位宽size有关，满足b×2m＜2size-4保证计算结果不溢出，n为不为零的自然数，共轭矩阵A中的元素为aij，i＝1,2,3,...,n，j＝1,2,3,...,n，1≤t≤T且t为自然数；
S2、初始化遍历次数计数器，令t＝0，
初始化特征向量初始矩阵，令V＝E，其中，E为单位阵；
S3、在S1所述共轭矩阵A中选取aij，初始化清扫元aij行列下标，令i＝1，j＝2；
S4、判断aij是否满足跳过清扫条件|real(aij)|＜a&|imag(aij)|＜a，若满足则转入S10，如不满足则转入S4；
S5、判断aij是否满足扩展条件|real(aij)|＜b&|imag(aij)|＜b，若满足则转入S6，若不满足则转入S7；
S6、进行位扩展，即计算a′ij＝aij×2m，转入S7；
S7、令a′ij＝aij，进入S8；
S8、计算根据所得计算
S9、计算A＝QHAQ和V＝QHV，其中，Q∈Cn×n为复数域内的平面旋转矩阵，
即Q的对角元素中除了qii＝ejφcosθ、qjj＝e-jφcosθ之外其他均为1，非对角元素中除了qij＝-ejφsinθ、qji＝e-jφsinθ其他元素均为0，θ为旋转角度；
S10、判断j＝n是否成立，是则进入S11，否则j＝j+1后跳转到S4；
S11、判断i＝n-1是否成立，是则进入S12，否则i＝i+1，j＝i+1后跳转到S4；
S12、判断t＝T是否成立，是则进入S13，否则t＝t+1后跳转到S3；
S13、输出迭代计算结果A和V，其中，A的对角元数位S1输入数据矩阵A的特征值，V为对应的特征向量矩阵。
2.根据权利要求1所述的一种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法，其特征在于：S1所述遍历次数T越大迭代次数越多，计算越精确，但计算时间越长，为了取得速度与精度的平衡，当n≤8时T＝3，当n＞8时T＝6。
3.根据权利要求1所述的一种基于Jacobi迭代算法的高精度矩阵特征值分解实现方法，其特征在于：S1所述a＝e×10-1。
CN201710130519.5A 2017-03-07 2017-03-07 基于Jacobi迭代算法的高精度矩阵特征值分解实现方法 Pending CN106940689A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710130519.5A CN106940689A (zh) 2017-03-07 2017-03-07 基于Jacobi迭代算法的高精度矩阵特征值分解实现方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710130519.5A CN106940689A (zh) 2017-03-07 2017-03-07 基于Jacobi迭代算法的高精度矩阵特征值分解实现方法

Publications (1)

Publication Number Publication Date
CN106940689A true CN106940689A (zh) 2017-07-11

Family

ID=59468750

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710130519.5A Pending CN106940689A (zh) 2017-03-07 2017-03-07 基于Jacobi迭代算法的高精度矩阵特征值分解实现方法

Country Status (1)

Country Link
CN (1) CN106940689A (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108228536A (zh) * 2018-02-07 2018-06-29 成都航天通信设备有限责任公司 使用FPGA实现Hermitian矩阵分解的方法
CN109245781A (zh) * 2018-09-18 2019-01-18 重庆九洲星熠导航设备有限公司 一种圆环共形阵空时抗干扰方法
CN110110285A (zh) * 2019-04-10 2019-08-09 浙江大学 一种用于FPGA的并行Jacobi计算加速实现方法
CN110222307A (zh) * 2019-06-12 2019-09-10 哈尔滨工程大学 基于fpga的实对称矩阵的特征值分解的并行实现方法
CN111859035A (zh) * 2020-08-12 2020-10-30 华控清交信息科技(北京)有限公司 数据处理方法及装置
CN112528224A (zh) * 2020-12-28 2021-03-19 上海微波技术研究所（中国电子科技集团公司第五十研究所） 一种矩阵特征值分解分组循环迭代流水实现方法及系统

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108228536A (zh) * 2018-02-07 2018-06-29 成都航天通信设备有限责任公司 使用FPGA实现Hermitian矩阵分解的方法
CN108228536B (zh) * 2018-02-07 2021-03-23 成都航天通信设备有限责任公司 使用FPGA实现Hermitian矩阵分解的方法
CN109245781A (zh) * 2018-09-18 2019-01-18 重庆九洲星熠导航设备有限公司 一种圆环共形阵空时抗干扰方法
CN110110285A (zh) * 2019-04-10 2019-08-09 浙江大学 一种用于FPGA的并行Jacobi计算加速实现方法
CN110110285B (zh) * 2019-04-10 2020-05-22 浙江大学 一种用于FPGA的并行Jacobi计算加速实现方法
CN110222307A (zh) * 2019-06-12 2019-09-10 哈尔滨工程大学 基于fpga的实对称矩阵的特征值分解的并行实现方法
CN110222307B (zh) * 2019-06-12 2022-10-28 哈尔滨工程大学 基于fpga的实对称矩阵的特征值分解的并行实现方法
CN111859035A (zh) * 2020-08-12 2020-10-30 华控清交信息科技(北京)有限公司 数据处理方法及装置
CN112528224A (zh) * 2020-12-28 2021-03-19 上海微波技术研究所（中国电子科技集团公司第五十研究所） 一种矩阵特征值分解分组循环迭代流水实现方法及系统

Similar Documents

Publication Publication Date Title
CN106940689A (zh) 基于Jacobi迭代算法的高精度矩阵特征值分解实现方法
CN101763445B (zh) 一种高光谱图像降维芯片
CN110222307B (zh) 基于fpga的实对称矩阵的特征值分解的并行实现方法
Ambikasaran et al. Fast direct methods for Gaussian processes and the analysis of NASA Kepler mission data
CN106919537A (zh) 一种基于FPGA的Jacobi变换的高效实现方法
Bhakthavatchalu et al. A comparison of pipelined parallel and iterative CORDIC design on FPGA
Ling et al. A stochastic extended Rippa’s algorithm for LpOCV
Jiang et al. An efficient FPGA-based direct linear solver
Bärthel et al. Hardware implementation of a latency-reduced sphere decoder with SORN preprocessing
Zhang et al. Efficient structure preserving schemes for the Klein–Gordon–Schrödinger equations
WO2020206716A1 (zh) 一种用于FPGA的并行Jacobi计算加速实现方法
Wahid et al. Hybrid architecture and VLSI implementation of the Cosine–Fourier–Haar transforms
Francescatto et al. A semi‐coarsening strategy for unstructured multigrid based on agglomeration
Chhajed et al. BitMAC: bit-serial computation-based efficient multiply-accumulate unit for DNN accelerator
CN104156953B (zh) 一种基于离散正交矩的图像重构方法
Wu et al. Low-latency low-complexity method and architecture for computing arbitrary Nth root of complex numbers
Ibrahim et al. Assessment of FPGA implementations of one sided Jacobi algorithm for singular value decomposition
Sethumadhavan et al. A case for hybrid discrete-continuous architectures
CN105182299B (zh) 基于复旋转矩阵的非正交联合对角化回波信号处理方法
Chen et al. Hardware efficient massive MIMO detector based on the Monte Carlo tree search method
Bottegal et al. Outlier robust kernel-based system identification using ℓ 1-Laplace techniques
Bangqiang et al. Base-N logarithm implementation on FPGA for the data with random decimal point positions
US8452830B2 (en) Programmable CORDIC processor with stage re-use
Irturk et al. An efficient FPGA implementation of scalable matrix inversion core using QR decomposition
CN107590105B (zh) 面向非线性函数的计算装置及方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WD01 Invention patent application deemed withdrawn after publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20170711