一种实对称矩阵特征值快速求解方法研究

彭永健 廖晓慧 张鹏飞

摘要：实对称矩阵特征分解在工程项目中经常遇到需要快速求解特征值和特征向量，文章提出了一种基于FPGA的高效并行实对称矩阵特征值分解方法。通过构造旋转矩阵一次消除矩阵多个元素，让矩阵快速收敛为对角阵，从而得到原实对称矩阵的特征值和特征向量。并根据工程误差需求通过Matlab仿真确定了算法循环次数，在FPGA上验证了算法的可行性和快速性。

关键词：实对称矩阵；特征值分解；快速收敛；旋转矩阵

中图分类号：TP3  文献标志码：A

0 引言

在工程应用中，经常需要计算实对称矩阵的特征值和特征向量，特别是信号子空间相关的阵列算法中，特征分解是其中的关键步骤［1］。针对实对称矩阵，通常思路是采用正交变换法计算得到全部特征值和特征向量［2］，即通过多次正交相似变换将实对称矩阵转化为相似矩阵，根据计算过程可以得出原矩阵的特征值和特征向量。

QR分解算法无需将矩阵完全化简，利用Householder变换，将原矩阵化简为对称三对角阵，再通过QR分解迭代，求得特征向量以及特征值［3］。雅可比分解算法通过一系列平面旋转矩阵迭代，最终将原矩阵转化为对角阵［4］。大部分实现算法都是非线性运算，而FPGA的优势主要是并行线性运算带来运算速度的提升。本文介绍一种基于FPGA的高效并行实对称矩阵特征值分解方法。

1 特征值分解计算方法

1.1 单点消除的雅克比分解

两个相似矩阵具有相同的特征值和特征向量［5］。根据这一基本原理，雅可比分解算法通过构造旋转矩阵，多次正交相似变换后，将N阶实对称矩阵A的非对角线元素逐渐变换为零，最终趋近于对角阵，此时对角线元素即为特征值，所有旋转矩阵可计算得到特征向量［6］。

针对N阶实对称矩阵A构造旋转矩阵，对角线上的元素只有pii和pjj为cosθ，其他位置为1；上三角的元素中，只有pij为-sinθ，其他位置为0；下三角的元素中，只有pji为sinθ，其他位置为0；m=-apq，n=0.5×（aqq-app），sinθ=sign（n）×sign（m）×m2/（m2+n2）2［1+n2/（m2+n2）］，cosθ=1-m2/（m2+n2）2［1+n2/（m2+n2）］。

通过矩阵计算A1=PTAP，可知该运算只改变原矩阵A的第i行、第i列、第j行、第j列，并且aij和aji变为0。继续该过程，依次将上三角中的元素消为0，共需要N×（N-1）/2次运算可以遍历一次所有非对角线元素，该过程记为一个循环LOOP。在一个LOOP中，将后面的元素消为0的过程会将前面已经消为0的元素改变，但整个LOOP结束后，对角线元素的平方和增大，非对角线元素的平方和减小。经过若干个LOOP后，矩阵逼近于一个对角阵。此时，对角线上的元素就是原矩阵A的特征值，所有的旋转矩阵P相乘就是特征向量矩阵。

1.2 多点消除的雅克比分解

上述串行分解运算过程中，一个LOOP中每一步仅消去上三角中的一个元素，PTAP涉及两次矩阵乘法运算，共需要N×（N-1）次矩阵乘法，在FPGA中实现该过程耗时较长。

消去上三角中元素aij时，只改变原矩阵A的第i行、第i列、第j行、第j列，而消去上三角中元素akl计算旋转矩阵时，仅需要用到akk，all，akl，因此可以同时消aij和akl，只需要满足i，j，k，l互不相等即可。对于N阶矩阵A，可以构造一个矩阵P，一次将上三角中的N/2个元素消为0，需要（N-1）次运算即可将非对角线元素遍历一次，矩阵乘法的次数降为2（N-1）。

1.3 实对称矩阵特征值和特征向量求解方法

根据上文所述，以16阶实对称矩阵为例，并行雅可比分解法计算特征值和特征向量的实现步骤如下。

（1）根据公式计算旋转矩阵P11。

（2）计算A11=PT11AP11，将8个元素消为0。

（3）重新选定8个元素位置，计算旋转矩阵P12。

（4）计算A12=PT12A11P12。

（5）重复步骤（3）和（4），直到上三角中的120个元素全都完成一次消除，一個循环LOOP内的消除顺序如下表1所示，表中（x，y）代表第x行、第y列。

（6）重复进行M个LOOP，直到对角线元素趋于0。

（7）对角线元素即为特征值，所有的旋转矩阵相乘，最终得到的矩阵每一列都是对角阵相应列特征值的特征向量。

2 Matlab仿真

针对以上算法，编写Matlab程序，仿真每次消单点和消多点在相同的循环迭代次数情况下，收敛速度是否有区别，并以此确定循环迭代次数LOOP的值。收敛速度按照绝对误差为判断标准，将Matlab系统函数Eig的计算结果作为理论值，雅可比分解算法得到的特征值与理论值的差值作为绝对误差。针对不同的LOOP值，随机生成16阶实对称矩阵，仿真次数为10 000。

由仿真结果可知，当LOOP值为6时，消单点的计算结果误差为10^（-6），消多点的计算结果误差为10^（-9），可以看出消多点雅可比算法能够更快的收敛到对角阵，得到与理论值误差更小的特征值。并且，10^（-9）的精度能满足大部分工程应用场景，后续的FPGA实现按照LOOP为6实施。

3 FPGA实现

在FPGA中实现并行雅可比分解算法如图1所示。每次计算旋转矩阵P后，需要计算正交变换A1=PTAP，并且将矩阵P累乘。正交变换和矩阵累乘都需要用到矩阵乘法操作，设计时将矩阵累乘和计算旋转矩阵并行，能够节省整个循环迭代的时间，即计算出旋转矩阵后，先进行正交变换，再同时开展矩阵累乘和计算下一个旋转矩阵的工作。并且，这样操作能够复用矩阵乘法资源，仅例化一套矩阵乘法运算资源，交替进行正交变换运算和变换矩阵累乘运算。

整個实现过程中，资源消耗量大的是矩阵乘法运算。以16阶矩阵乘法B=A×P为例，计算B中的1个点需要16个乘法器、15个加法器，计算所有点共需要4 096个乘法器、3 840个加法器，DSP资源消耗巨大，无法在一片FPGA上实现。因此，需要做适当的串行化，降低资源消耗。经过时间和资源的折中考虑，在本设计中例化16套资源，并行计算B中的16个点，再重复使用16次完成B中所有点的计算。资源消耗为1 248个DSP，计算时间为75个时钟周期。

在Modelsim中对开发的算法程序进行仿真。从输入矩阵A，到计算出特征值和特征向量（LOOP次数为6），共需要29 000个时钟周期，使用300 MHz的计算时钟，所需时间为96.6 us，绝对误差最大仅10^（-6），满足工程应用需求。

4 结语

本文提出了一种高效并行实对称矩阵特征值分解方法，通过Matlab仿真确定了循环次数，并在FPGA上验证了算法的有效性。

参考文献

［1］陈建华.线性代数［M］.4版.北京：机械工业出版   社，2017.

［2］徐士良，马尔妮.常用算法程序集（C/C++描述）［M］.北京：清华大学出版社，2013.

［3］郝英军.MUSIC算法中特征值分解及信源数估计的FPGA实现［D］.哈尔滨：哈尔滨工程大学，2020.

［4］胡乐宇，蔡邢菊.低计算精度下实对称矩阵的特征值分解［J］.高等学校计算数学学报，2021（2）：117-133.

［5］胡乐宇.线性方程组的随机求解算法以及低精度矩阵的特征值分解［D］.南京：南京师范大学，2020.

［6］唐浩伟，罗林，王涛.基于特征值分解的时间反转超声成像技术［J］.信息技术，2019（9）：52-55.

（编辑 沈 强）

Research on a fast method for solving eigenvalues of real symmetric matrix

Peng  Yongjian， Liao  Xiaohui， Zhang  Pengfei

（The 29th Research Institute of China Electronics Technology Group Corporation， Chengdu 610029， China）

Abstract：  Eigendecomposition of real symmetric matrices is often encountered in engineering projects， and it is necessary to quickly solve eigenvalues and eigenvectors. This paper proposes an efficient parallel real symmetric matrix eigenvalue decomposition method based on FPGA. By constructing a rotation matrix to eliminate multiple elements of the matrix at one time， the matrix can quickly converge to a diagonal matrix， and the eigenvalues and eigenvectors of the original real symmetric matrix can be obtained. According to the engineering error requirements， the cycle times of the algorithm are determined by Matlab simulation， and the feasibility and rapidity of the algorithm are verified on FPGA.

Key words： real symmetric matrices; eigenvalue decomposition; quickly converges; rotation matrix