一种厄米特矩阵特征值并行求解方法

严新文 彭永健 张续莹

摘 要：厄米特矩阵特征值特征向量快速求解是在实际工程中常常会碰到的问题，但大部分基于计算机实现的算法都是利用串行思想进行编程求解，不能满足实时求解需求。本文提出基于FPGA的并行算法来求解本类问题。首先，将厄米特矩阵实数化。其次，通过并行雅克比分解和特征值快速排序来恢复厄米特矩阵的特征值特征向量。最后，通过MATLAB仿真验证了算法的有效性，并通过蒙特卡洛仿真确定了雅克比sweep的迭代次数，同时通过FPGA实现验证了算法的实时性。

关键词：厄米特矩阵 特征值 并行雅克比分解 蒙特卡洛仿真

中图分类号：O129 文献标识码：A 文章编号：1672-3791（2018）04（b）-0108-02

矩阵特征值特征向量求解一直是线性代数中最重要的操作之一，在科学计算中有着非常广泛的应用[1]。工程技术和物理等领域中的很多问题在数学上都归结为求矩阵的特征值特征向量问题。例如振动问题和物理学中某些临界值的确定，均归结为求矩阵的特征值问题。

矩阵特征值和特征向量的求解方法有很多，幂法可以求解一般矩阵绝对值最大的特征值，雅可比方法可以求解对称矩阵的全部特征值，QR方法可以求解一般实矩阵的全部特征值[2]。在一般的工程当中，常用的矩阵特征分解方法大部分都是针对实矩阵的，而在实际工程中，我们常需要处理复矩阵，而且大部分情况下为一个厄米特矩阵，由于厄米特矩阵具有对称性，能简化其特征值特征向量求解过程。

但是上述算法的计算机实现大部分都是利用串行思想进行编程求解，并不适合用于并行计算，本文详细介绍一种可用于FPGA并行高效实现的厄米特矩阵特征值特征向量求解方法[3]。

1 厄米特矩阵特征值求解方法

1.1 实数化

如果直接对厄米特矩阵进行特征分解显然是不合适的，也是没有效率的，因此我们在对其进行特征分解前先进行一项预处理过程，即将它实数化。

1.2 并行雅克比分解

对于任意的一个实对称矩阵A，只要能够求得一个正交矩阵U，使得 成为一个对角阵D，就得到了A的所有特征值和对应的特征向量。雅可比法就是基于这一思想，用一系列的初等正交变换逐步消去A的非对角线元素，从而使矩阵A对角化。

基本雅可比法需要逐次寻找非主对角元中绝对值最大元素，这种方法并不适合于并行运算，因此我们可以按照一定的顺序将A中的所有上三角元素全部消去一遍，这样的过程称为一次扫描（sweep）。由于对称性，在一次扫描中，A的下三角元素同时被消去了一遍，被消为0的元素在消去后面元素的旋转过程中可能又变为非0。但经过若干轮之后，可使S收敛于一个对角阵。

并行雅克比分解，通过构造旋转矩阵，可以一次将矩阵N阶矩阵A上三角中的N/2个点消为0，经过（N-1）次旋转变换可将上三角中的每个点都消一次，使矩阵A接近对角阵，该过程作为一个迭代过程。为满足一定的精度要求，可循环迭代多次，迭代次数通过蒙特卡洛仿真确定。

1.3 特征值排序

排序逻辑是结合冒泡排序和归并排序的分组处理思想而设计的。冒泡排序是每次选出当前序列中的最大或最小元素，重复这个过程，直到输入序列中的元素全部有序输出时结束。归并排序算法首先将输入序列分为若干个组，每组内排序，对这些已排序的有序子序列再次分组，同组内的有序序列归并为一个有序序列，重复这个过程直到所有元素有序排列。

1.4 厄米特矩阵特征值求解方法

根据上文所述，按照并行雅可比法计算厄米特矩阵特征值和特征向量的整体实现流程步骤如下：

（1）将输入的N阶复数矩阵B升维为2*N阶实数对称矩阵A。

（2）根据雅可比基本原理的公式计算旋转矩阵R_temp。

（3）计算特征值，即计算R_temp'*A*R_temp，该步骤可将矩阵A的N个点消为0，经过多次循环迭代后，矩阵A接近对角阵，对角阵上的元素即是特征值。

（4）计算特征向量，将每次迭代的旋转矩阵R_temp累乘，最终得到特征向量。

（5）将特征值按照从小到大顺序排序，选择奇数位置的特征值，同时选择对应的特征向量，恢复出N阶复数特征向量。

2 仿真

本文使用matlab蒙特卡洛仿真来验证并行雅可比算法的有效性。仿真中构造随机8*8复数赫尔米特矩阵R，通过设置不同的迭代次数，分别使用上述并行雅可比算法进行特征值分解，以及matlab系统自带函数eig特征值分解，计算两者特征值的差值，并取差值的最大值作为绝对误差。设置仿真循环次数10万次，统计每次仿真的绝对误差，绘制绝对误差直方图。

设置雅可比sweep迭代次数为4，在此情况下，绝对误差最大可达30，大部分误差都集中在2以内。设置雅可比sweep迭代次数为5，在此情况下，绝对误差最大为5×10^（-6），大部分误差都集中在10^（-6）以内。设置雅可比sweep迭代次数为6，在此情况下，绝对误差最大为7×10^（-10），大部分误差都集中在10^（-9）以内。为使绝对误差达到10^（-6），在后续实现过程中迭代次数设置为5。

在Vivado 2016.2的环境下，使用VHDL按照单精浮点模型完成代码开发，然后编写测试激励信号，在Modelsim环境下仿真如图5。从输入矩阵A，到完成特征值和特征向量计算，共需要11576个时钟节拍。在系统时钟为150MHz情況下，所需时间为11576×6.66ns=77.1us。

将Modelsim仿真得到的结果与matlab计算结果对比，可以看出结果基本一致，误差为10^（-6）。

3 结语

本文提出了一种通过将厄米特矩阵实数化，然后再通过并行雅克比分解，特征值快速排序来恢复厄米特矩阵的特征值特征向量的方法。最后通过MATLAB仿真验证了算法的有效性，并通过蒙特卡洛仿真确定了雅克比sweep的迭代次数，同时通过FPGA实现验证了算法的实时性。

参考文献

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

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

[3] 王涛.求矩阵特征值的GPU并行算法的研究[D].黑龙江：黑龙江大学，2012.