OPENCL设计基础 基于Arria10硬浮点的矩阵分解
【时间：2017 年 11 月 01 日, 来源：本站原创, 作者：Jerry Zheng】
一、概述
QR矩阵分解( qrdQRD)，也称为正交矩阵三角化，基本原理是将矩阵( A )的分解为正交矩阵(Q )和上三角矩阵( R)。QRD是求解最小二乘问题和联立方程的有效方法。QRD在无线应用有很多的有效应用场景，如在OFDM的MIMO系统通常需要对对个低阶矩阵（如4*4矩阵）求逆，这种类型的系统通常需要采用典型的非递归技术如QRD。
二、QRD算法介绍
参考如下方程：
AY = Z；
其中：
A，Y，Z均为矩阵，A为N*N的方阵，Y，Z为N阶向量。A与Z是已知的，Y是未知的，最终目的就是确定Y 矩阵的值。
将矩阵A用QR替换得到如下结果：
(QR)Y= Z；
将Q移到方程右边得到如下结果：
RY= Q-1 Z；
Q是一个正交(酉)矩阵，因此Q -1等于复数Q的共轭转置，因此便可得到如下方程：
RY= Z'；
其中
Z' = Q-1 Z；
实现QRD算法包括gram - Schmidt正交，Householder变换和Givens转换，由于Givens转换具有并行化的优势，下面主要介绍基于Givens转换的QRD实现。
我们可以采用每个单元均实现Givens转换的结构实现QRD，图1显示了不同输入如何反馈到不同收缩阵列的不同列。
当所有输入数据输入到收缩阵列并通过数组传播，数组所有列中的值(不包括最后一列)
提供上三角R矩阵。最后一列提供Z’的值，及Q-1Z的值。确定R跟Z’后便可以求解方程RY = Z' 得到Y.
三、A10 DSP模块特性
Arria10 可变精度 DSP Block同时支持浮点运算及定点运算 ，浮点 运算时特性：
1、 完整 硬件结构支持乘法、 加法、 减法、乘加、乘减等浮点运算；
2、 带累加及动态累加器复位控制浮点乘法；
3、 带加法 或减法 级联的浮点乘法；
4、 浮点复数乘法运算；
5、 直接向量点积运算；
6、 脉动 FIR滤波器 ；
四、QRD算法实现
下面是4阶方阵QRD Matlab的实现代码，并给出了直接通过Matlab QR分解函数和该算法的结果对比：
m=4;
A = randn(m,m)+i*randn(m,m);
[q,r] =qr(A) %%Matlab QR分解函数
Q = zeros(m,m);
R = zeros(m,m);
%% k=1
R(1,1) = sqrt(sum(abs(A(:,1)).^2));
Q(:,1) = A(:,1)/R(1,1);
%% k= 2
R(1,2) = Q(:,1)'*A(:,2);
A(:,2) = A(:,2)-Q(:,1)*R(1,2);
R(2,2) = sqrt(sum(abs(A(:,2)).^2));
Q(:,2) = A(:,2)/R(2,2);
%% k= 3
R(1,3) = Q(:,1)'*A(:,3);
A(:,3) = A(:,3)-Q(:,1)*R(1,3);
R(2,3) = Q(:,2)'*A(:,3);
A(:,3) = A(:,3)-Q(:,2)*R(2,3);
R(3,3) = sqrt(sum(abs(A(:,3)).^2));
Q(:,3) = A(:,3)/R(3,3);
%% k= 4
R(1,4) = Q(:,1)'*A(:,4);
A(:,4) = A(:,4)-Q(:,1)*R(1,4);
R(2,4) = Q(:,2)'*A(:,4);
A(:,4) = A(:,4)-Q(:,2)*R(2,4);
R(3,4) = Q(:,3)'*A(:,4);
A(:,4) = A(:,4)-Q(:,3)*R(3,4);
R(4,4) = sqrt(sum(abs(A(:,4)).^2));
Q(:,4) = A(:,4)/R(4,4);
Q
R
q =
-0.2696 + 0.2758i 0.0926 -0.1033i 0.3199 - 0.6630i -0.0668 + 0.5344i
-0.2042 + 0.2798i -0.4199 -0.0996i 0.4549 + 0.3194i 0.6203 + 0.0077i
-0.7428 - 0.2179i -0.3284 +0.3758i -0.3217 - 0.1656i 0.0147 - 0.1438i
0.3629 + 0.0195i -0.4886 +0.5539i -0.0138 + 0.1345i -0.1966 + 0.5153i
r =
3.9640 + 0.0000i -2.1874 +1.5695i -0.0841 + 0.8639i 0.1994 + 0.4748i
0.0000 + 0.0000i 1.8159 +0.0000i 0.4109 - 1.0269i 0.6065 - 1.2118i
0.0000 + 0.0000i 0.0000 +0.0000i 0.5670 + 0.0000i 0.5431- 1.2020i
0.0000 + 0.0000i 0.0000 +0.0000i 0.0000 + 0.0000i -0.1807 + 0.0000i
Q =
-0.2696 + 0.2758i 0.0926 -0.1033i 0.3199 - 0.6630i 0.0668 - 0.5344i
-0.2042 + 0.2798i -0.4199 -0.0996i 0.4549 + 0.3194i -0.6203 - 0.0077i
-0.7428 - 0.2179i -0.3284 +0.3758i -0.3217 - 0.1656i -0.0147 + 0.1438i
0.3629 + 0.0195i -0.4886 +0.5539i -0.0138 + 0.1345i 0.1966 - 0.5153i
R =
3.9640 + 0.0000i -2.1874 +1.5695i -0.0841 + 0.8639i 0.1994 + 0.4748i
0.0000 + 0.0000i 1.8159 + 0.0000i 0.4109 - 1.0269i 0.6065 - 1.2118i
0.0000 + 0.0000i 0.0000 +0.0000i 0.5670 + 0.0000i 0.5431 - 1.2020i
0.0000 + 0.0000i 0.0000 +0.0000i 0.0000 + 0.0000i 0.1807 + 0.0000i
通过以上算法，实现QRD包含了浮点复数模值，浮点复数矩阵乘法，浮点复数矩阵加乘等预算，通过合理调用硬符点DSPblock实现相关复数浮点运算。下面是基于A10硬浮点结构实现4*4矩阵QR分解的部分代码，基于该参考代码用户实现了8阶的矩阵的QR分解，基于A10硬浮点的优势，实现8阶的QR分解仅需要用到66个DSPblock，并且延时仅为500个时钟周期，相比软浮点设计大幅度提升。
如欲了解更多关于该方案的信息，请联系骏龙科技各地的办事处，或发邮件至stonemao@macnica.com。
骏龙科技微信公众号：
骏龙Altera FPGA微信公众号：
