浮点厄米特矩阵特征值分解的FPGA实现

王传根 吕乐群 林都督

摘要：为了提高厄米特矩阵特征值分解计算的速度和精度，基于QR迭代的方法，在FPGA上实现了一种并行可重构的特征值分解方案。设计在Intel公司型号为10AX066H1F341的FPGA上进行验证，通过调用硬浮点DSP计算资源，实现了8×8阶的浮点厄米特矩阵的特征值分解。验证结果显示：该方案耗时11425个时钟周期，计算误差小于3.74e-7，资源消耗小于总资源的32%，最高运行时钟频率可达311MHz，可支撑MUSIC等复杂算法的高速硬件工程实现。

关键词：厄米特矩阵;特征值分解;可重构;FPGA

中图分类号：TN911.7 文献标识码：A 文章编号：1007-9416（2020）02-0133-03

在阵列信号处理领域中，基于特征分解的超分辨率空间谱估计算法以其高效高精度的特点，一直以来都是研究的焦点，这其中最为经典的就是多重信号分类（MUSIC）算法[1-2]。MUSIC算法具有测向精度高、算法稳定、优异的多信号测向能力等众多优点，然而其运算却尤为复杂，不利于硬件实现。这其中以矩阵特征值分解（EVD）的运算最为复杂，大约占到了MUSIC计算60%的计算量，所以研究EVD在FPGA上的加速运算有着十分重要的意义。当前研究EVD在FPGA上的实现多是针对实对称矩阵，并且多是采用定点的方式实现[3-5]。然而在实际工程应用中，经常需要计算复数自共轭矩阵（即厄米特矩阵）的特征值分解，并且随着矩阵阶数增大，定点数据在计算过程中的截位误差累积，造成最终结果精度大幅下降，难以满足应用需求。

常用的EVD分解的方法主要有基于Jacobi旋转的方法和基于QR迭代的方法，基于Jacobi的方法精度高、收敛速度慢，且不能直接计算复数矩阵的特征值分解。相比Jacobi方法，基于QR迭代的方法精度稍差，但是其收敛速度更快，并行度高，可以直接对复数矩阵进行计算，对于加速计算的效果更为明显，特别适合于FPGA并行实现。

本文以8×8阶的厄米特矩阵特征值分解为目标，采用并行QR迭代的方案实现FPGA设计。通过采取浮点运算机制和算法结构优化，在保证计算精度的同时，加快计算速度，降低资源消耗。

1 算法原理

1.1 基于QR分解的EVD计算

对于非奇异方阵A，进行QR分解，可得A=QR。其中Q为正交矩阵（即QTQ=I），R为上三角矩阵，则：

（1）

也即是说，RQ与A是正交相似的，它们具有相同的特征值。因此可以采用如下的QR迭代格式。

（2）

取A1=A，经过式（2）迭代，可以证明Ak收敛到矩阵A的特征值，Uk=Q1Q2…Qk收敛至矩阵A的特征向量。

1.2 基于QR分解的EVD计算

根据正交理论的思想，矩阵QR分解可以采用Givens旋转变换实现，对矩阵Am×m=（a1，a2，…，am）T，Givens旋转矩阵定义为：

（3）

其中，，，，且Gi，j为酉矩阵。执行Givens旋转变换，可使矩阵按照一定角度旋转，使矩阵i行j列元素置换为0。由此，可以重复执行这个变换过程，生成一系列Givens矩阵Gm，m-1…Gm，2…G3，2Gm，1…G2，1，将对角线以下的元素全部置换为0，得到上三角矩阵R，于是有：

（4）

按照（4）式所示的方法，可以迭代完成矩阵A的QR分解。然而每次使用Givens矩阵置换元素后，新的Givens矩阵必须基于Givens变换后的矩阵生成。由于前后数据的相关性，造成FPGA实现时只能顺序执行，无法发挥FPGA并行执行的优势，难以提升运算速度。由此，对Givens旋转矩阵进行改进，将每列Givens矩阵的计算展开，每次取矩阵A第j列对角线以下元素，令

（5）

（6）

将（5）式和（6）式代入（3）式计算Givens矩阵，一次得到j列的Givens矩阵，计算可将矩阵A的j列对角线以下的元素置换为0。因此可以每次并行计算，一次迭代置换1列元素，经过m-1次迭代即可完成QR分解。这样的计算步骤非常利于FPGA并行实现，从而极大地提高运算速度。

2 FPGA实现方案

为了提升运算效率，设计选用搭载硬浮点DSP计算单元的IntelArria10系列FPGA实现。该DSP单元采用芯片硬件定制浮点设计，遵循IEEE754标准，单个DSP即可实现一个浮点加法和一个浮点乘法，无需额外的逻辑资源。此外通过调整DSP内部流水寄存器级数，可以满足不同的latency和Fmax的需求，无需在timing上做过多努力，就能工作在极高的工作频率。在浮点计算密集的情况下，Arria10比起业界同等级的FPGA在资源、效率、运算速度方面均有明显優势。

在FPGA上实现EVD分解的结构如图1所示，设计包含4个模块：矩阵数据缓存（BUF），Givens矩阵计算单元（GCU）、可重构计算单元（RCU）和计算控制单元（CCU）。

BUF模块实现矩阵A_in输入缓存、迭代结果缓存和EVD输出结果缓存，为加快随机访问速度，BUF模块缓存均使用寄存器实现。CCU模块为整个EVD运算的控制枢纽，根据迭代次数iter_num、启动信号start以及各个模块的状态，控制迭代运算过程，完成EVD的计算，CCU模块主体为控制状态机，其状态转换如图2所示。按照相关性可将EVD计算分为4个计算步骤：

（a）计算Givens旋转矩阵G;

（b）计算旋转矩阵连乘G_prod;

（c）计算G=G_prod*G和A=G_prod*A;

（d）计算Ak+1=Rk*Qk和Uk+1=Uk*Qk。

这4个步骤输入输出上存在关联，无法拆解开来并行实现，但工作时间上是分开的，且计算相似性极高。故在保证步骤内高并行度的情况下，采用分时重构计算资源的设计方法。

GCU和RCU相互协作，根据CCU设置的工作模式，重构完成不同的计算功能。GCU实现计算过程中开平方和除法等个性的计算，同时设置RCU重构参数，调度RCU协处理向量二范数、复矩阵乘法等计算资源密集的共性操作。如此通过循环调度计算资源，即可完成EVD计算的各个步骤。RCU在不同工作模式下，完成如下计算：

在步骤（a）时，取矩阵的列向量按（5）式和（6）式计算各个旋转矩阵的c和s，此时需要求取列向量的二范数，重构RCU并行完成多个的计算;

在步骤（b）时，需要计算一系列Givens矩阵的连乘，对于j列生成的Givens矩阵Gm，j乘以矩阵B，根据（3）式矩阵的特点，可以简化计算，仅需计算乘积的第m行和j行，其余的行等于B，此时RCU完成（7）式的运算过程。

（7）

在步骤（c）和（d）时，需要计算两个8×8矩阵的乘法，此时RCU完成矩阵乘法的计算，即计算。

经过分析剥离，RCU中只保留复数乘和复数加两种操作，并以此为最小重构颗粒，设计RCU模块如图3所示。RCU将重构参数解析为数據选通MUX的选通信号，去选通相应的复数乘法器和加法器，实现以上RCU要求的功能，其中“取反”、“求共轭”操作由格式转换子模块实现。经分析，计算资源消耗最大的步骤为矩阵乘法计算，并行计算一个矩阵乘积元素需要8个乘法器和7个加法器。设计例化64个乘法器和56个加法器，可以一次并行计算8个（一行）的矩阵乘积结果，仅需8+N（N为计算延迟）拍就能完成一次矩阵乘法。其余的运算操作通过重构RCU，选通部分乘法器和加法器即可实现。通过采用RCU重构的方法，在步骤（a）～（d）中资源重构，可以节约60%以上资源。

3 设计实现及结果

本设计选用的器件为Intel公司的10AX066H1F341FPGA，使用Verilog语言进行电路描述。在Quartus上进行编译综合，并使用Modelsim进行仿真。为验证设计正确性，在matlab上生成好激励，同时将激励用于Modelsim仿真和matlab上eig函数计算，然后比对仿真和计算结果。

由于矩阵分解特征向量不唯一，FPGA实现方法和matlab中eig函数实现方法不一致，导致特征向量不同，不能直接数值比较。本文主要进行特征值数值比较，同时根据矩阵特征值分解的特点，验证特征向量是否为酉矩阵、特征值是否为对角线实数矩阵、VDVH是否等于矩阵A。为了评估设计的计算精度，按（8）式对特征值误差进行计算。其中rDi，j为仿真得到的特征值，gDi，j为eig计算得到的特征值。

（8）

通过matlab随机产生21组数据，仿真设置迭代次数为12次，对数据进行了仿真。结果显示，仿真得到的特征向量均为酉矩阵，特征值为对角线实数矩阵，符合矩阵特征值分解的要求。特征值误差如图4所示，计算误差均小于3.74e-7。仿真显示，一次计算需11425个时钟周期，耗时约57us（工作频率200MHz时）。在Quartus上综合报告如表1所示。

4 结语

本文基于Givens变换和QR分解的方法，采用并行可重构的设计策略，在FPGA上实现了厄米特矩阵特征值分解的计算。设计结果表明：在资源节省60%的情况下，仍能保证极高的精度和运算速度。该设计可以满足MUSIC算法实时信号处理的需求，也可以在人工智能、图像处理、盲源分离、频谱分析等有EVD需求的领域作为IP使用。

参考文献

[1] Schmidt R O.Multiple emitter location and signal parameter estimation[J].Antennas and Propagation，1986，34（3）：276-280.

[2] Stoica P，Nehorai A.MUSIC，maximum likelihood，and Cramer-Rao bound[M].IEEE Transactions on Acoustics，Speech，and Signal Processing，1989.

[3] 袁生光，沈海斌.基于Jacobi算法对称矩阵特征值计算的FPGA实现[J].机电工程，2008，25（10）：80-82.

[4] 刘永勤.对称矩阵特征值分解的FPGA实现[J].现代电子技术，2017（12）：15-18.

[5] 王飞，王建业，张安堂，等.实对称矩阵特征值分解高速并行算法的FPGA实现[J].空军工程大学学报（自然科学版），2008，9（6）：71-74.

FPGA Implementation of  Eigenvalue Decomposition for Floating-point

Hermite Matrix

WANG Chuan-gen，LV Le-qun，LIN Du-du

（The 29th Research Institute of  CETC， Chengdu  Sichuan  610000）

Abstract：In order to improve the speed and accuracy of eigenvalue decomposition of Hermite matrix， a parallel reconfigurable FPGA solution based on QR iteration was implemented. For verifying this design， eigenvalue decomposition of 8×8 Hermite matrix was realized by calling hard floating-point DSP on Intels 10AX066H1F341 FPGA. The verification results show that the calculation takes 11425 clock cycles， the error is less than 3.74e-7， the resource consumption is less than 32%， and the maximum running frequency can reach 311MHz， which can support the high-speed hardware engineering implementation of complex algorithms such as MUSIC.

Key words：Hermite matrix;eigenvalue decomposition;reconfigurable;FPGA