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

本发明属于信号处理领域，尤其涉及种基于jacobi迭代算法的高精度矩阵特征值分解实现方法。

背景技术：

在信号处理中，矩阵的特征值分解evd是一个应用广泛的矩阵运算。如数据压缩、噪声去除、数值分析，包括近几年兴起的机器学习、深度学习其基本核心操作也包括矩阵特征值分解。实现矩阵特征值分解的常用方法有gauss变换、householder变换、jacobi迭代等，其中，jacobi迭代是精度较高的方法，并且很适合在fpga中实现。因此一种基于jacobi迭代算法的高精度矩阵特征值分解实现技术在实际工程中具有很高的应用价值。

经典的jacobi迭代算法计算共轭矩阵a∈c^n×n的特征值分解如图1所示，这种经典的迭代算法虽然有较快收敛速度，但是该算法需要在矩阵a的众多元素中选取aij，使得aij为非对角元素中绝对值最大的一个，再进行后面的计算操作。这样每一步都要寻找绝对值最大的非对角元，比较费时也不适合在fpga实现，因此经典的jacobi迭代算法在实际工程中并不实用。

目前实际工程中多数采用如图2所示的循环jacobi迭代算法，通过逐行扫描遍历法选取aij，这样避免了寻找最大绝对值的非对角元的复杂繁琐步骤。这样选取aij的方式，在aij数值比较大时，fpga中使用cordic算法计算φ、θ误差比较小，可以取得比较好的效果。但当aij比较小甚至接近0时，此时fpga中使用cordic算法计算φ、θ误差比较大，就会导致后面计算a＝q^haq产生误差，其中q∈c^n×n为复数域内的平面旋转矩阵。而计算过程需要多次的迭代运算，如果在迭代过程中出现多次aij比较小甚至接近0的情况，就会产生较大的累计误差，从而使最终计算结果的精度相对较差。

技术实现要素：

发明的目的在于解决在循环jacobi迭代算法进行矩阵特征值分解过程中，因为逐行扫描遍历法选取aij在aij比较小甚至接近于0时，导致fpga中使用cordic算法计算φ、θ误差比较大，进而使迭代过程产生相对较大的累计误差，导致计算结果误差增大。即提供一种基于jacobi迭代算法的高精度矩阵特征值分解实现方法，在不明显增加算法复杂度、fpga实现难度与增加资源消耗的情况下，提高实际工程中基于循环jacobi迭代算法的fpga实现矩阵特征值分解的计算精度。

一种基于jacobi迭代算法的高精度矩阵特征值分解实现方法，包括如下步骤：

s1、设数据矩阵a∈c^n×n为共轭矩阵，并且设定最大遍历次数为t、最小清扫门限a、扩位门限b和算术左移位数m，其中，最小清扫门限a应小于要求计算结果的精度一个数量级，扩位门限b与算术左移位数m与fpga实现里数据位宽size有关，满足b×2^m＜2^size-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×2^m，转入s7；

s7、令a′ij＝aij，进入s8；

s8、计算根据所得计算

s9、计算a＝q^haq和v＝q^hv，其中，q∈c^n×n为复数域内的平面旋转矩阵，

即q的对角元素中除了qii＝e^jφcosθ、qjj＝e^-jφcosθ之外其他均为1，非对角元素中除了qij＝-e^jφ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λ，其中，λ为信号波长，阵列接收噪声是功率为σ²＝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为均值为零、方差σ²＝1的高斯白噪声矢量；远场接收信号s(k)＝ass(k)，其中其幅度as＝10^snr/20；a(γ)＝[1e^-jφ…e^-j(n-1)φ]^t，为空间阵列的n×1维流型矩阵。

a2.由计算数据协方差矩阵r∈c^n×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×2^m＜2^size-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×2^m，进入步骤b26。

b25.不进行位扩展，即a′ij＝aij，进入步骤b26。

b26.计算相角和模值，即和进入步骤b27。

b27.计算平面旋转角，即进入步骤b28。

b28.进行jacobi旋转，即计算a＝q^haq和v＝q^hv，其中q∈c^n×n为复数域内的平面旋转矩阵。

即q的对角元素中除了qii＝e^jφcosθ、qjj＝e^-jφcosθ之外其他均为1，非对角元素中除了qij＝-e^jφ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∈c^n×n，进入步骤c2。

其中，v∈c^(n-1)×(n-1)为前面计算块矩阵r′的特征向量，满足vv^h＝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％。

所以，本发明在基本不增加资源消耗的情况下，不仅可以提高计算精度，还可以提高计算速度，在实际工程中具有重要价值。