（通信与信息系统专业论文）基于fpga的矩阵运算实现.pdf

硕士论文基于f p g a 的矩阵运算实现 摘要 密集型的矩阵运算在信号处理和图像处理中被广泛应用，而且往往需要系统进行 实时运算，这就需要系统具有很高的吞吐率。因此寻找矩阵运算的高速实现方法是很 有意义的。f p g a 的运算速度快并且可以并行运算，和其它矩阵运算的实现方式相比， f p g a 有其独特的优势。本文主要设计并实现了基于f p g a 的各种矩阵运算模块。 本文首先介绍了矩阵运算的特点和原理，接着讨论了f p g a 浮点运算单元的 v h d l 设计方法，在此基础上，设计了矩阵相乘累加、三角矩阵求逆和一般矩阵分解 求逆的运算模块，给出矩阵阶数扩大时各种矩阵运算的分块实现方法。然后在 m o d e l s i m 环境下仿真了一般矩阵的求逆模块，与m a n a b 仿真结果比较，分析了运算精 度、时间复杂度和资源占用情况，在v i r t e x 4 系列f p g a 硬件平台上进行了调试和测试， 并通过u s b 接口将矩阵运算结果送入p c 机，验证了基于f p g a 矩阵运算的正确性和可 行性。最后对矩阵求逆模块在雷达信号中的应用作了简单介绍。 关键词：矩阵运算，f p g a ，并行结构，浮点运算 硕士论文 基于f p g a 的矩阵运算实现 i n t e n s i v em a l r i xo p e r a t i o n sa r ew i d e l yu s e di ns i g n a lp r o c e 醛i n ga n di 1 1 硷g pp r o c e s s i n g 1 惋so p e r a t i o n sa l w a y sr e q u k et h es y s t e mf o rr e a l - t i m e a n dt h ea r c h i t e c m r e 、v i t l lh i 班 t h r o u g h p u tr a t ei st h e r e f o r eh i g h l ys o u g h tf o r s oi t su s e f u lt of i n dw a y st or e a 比bm a t r i x o p e r a t i o n sf o rh i g hs p e 虻d c o m p a r e dw i t ht h eo t h e rm e t h o d , u s i n gf p g a h a si t so w n u n i q u ea d v a n t a g e sb e c a u s eo f i t sh i 9 1 ls p e e d , a n di t se a s yt or e a l i z ep a r a l l e lc o m p u t a t i o n s h l 】c t l 搬b yf p g a 1 1 1 i sd i s s e r t a t i o nm a i n l yd i s c u s s e st h er e a l i z a t i o no fm a t r i xo p e r a t i o n b a s e do nf p g a f i r s t , t h ec h a r a c t e r i s t i ca n dt h et h e o r e t i co fm a t r i xo p e r a t i o n sa 托i n t r o d u c e d a n dt h e d e s i g no ff l o a t i n g - p o i n to p e r a t i o nu n i tu s i n gv h d l i sd e s c r i b e d o nt h i sb a s i s ，t h em a t r i x m u l t i p l i c a t i o nc u m u l a t i v em o d u l e ，t h et r i a n g u l a rm a t r i xi n v e r s i o nm o d u l ea n dt h eg e n e r a l m a t r i xi n v e r s i o nm o d u l ea r ed e s i g n e d i no r d e rt os o l v et h ep r o b l e m sw h i c hb r o u g h tb yt h e e x p a n s i o no fm a t r i xo r d e r , t h ep a r t i t i o n e dm a t r i xa l g o r i t h mt or e a l i z eh i g ho r d e rm a t r i x o p e r a t i o ni sd i s c u s s e d t h e nt h eg e n e r a lm a t r i xi n v e r s i o nm o d u l ei ss i m u l a t e di nm o d e l s i m 1 1 l cp r e c i s i o no fr e s u l tc o m p a r e dw i t hm a t l a b t i m ea n dr e s o u r c e so c c u p i e db yf p g aa 坞 a n a l y z e d t h e nt h ep r o g r a m m e ri sd o w n l o a d e do nv i r t e x - 4f p g ah a r d w a r ep l a t f o r mt o r e a l i z eh i i xi n v e r s i o n a n dt h er e s u l ti st r a n s f e r e dt op cb yu s b a c c o r d i n g l y , t h e c o r r e c t n e s sa n d 龟a s i b i l i t yo fm a t r i xo p e r a t i o nb 髂e do nf p g aa r ev a l i d a t e d 。a tl a s t , t h e u s eo f m a t r i xi n v e r s i o ni nr a d a rs i g n a lp r o c e s s i n gi ss i m p l yi n t r o d u c e d k e y w o r d s ：m a t r i xo p e r a t i o n , f p g a ，p a r a l l e ls t r u c t u r e ，f l o a t i n g - p o i n t o p e r a t i o n i i 声明 本学位论文是我在导师的指导下取得的研究成果，尽我所知，在本 学位论文中，除了加以标注和致谢的部分外，不包含其他人已经发表或 公布过的研究成果，也不包含我为获得任何教育机构的学位或学历而使 用过的材料。与我一同工作的同事对本学位论文做出的贡献均已在论文 中作了明确的说明。 研究生签名： 臃7 a 占日 学位论文使用授权声明 南京理工大学有权保存本学位论文的电子和纸质文档，可以借阅或 上网公布本学位论文的部分或全部内容，可以向有关部门或机构送交并 授权其保存、借阅或上网公布本学位论文的部分或全部内容。对于保密 论文，按保密的有关规定和程序处理。 硬士论文 基于f p g a 的矩阵运算实现 1 绪论 1 1 研究背景 近年来由于计算机的发展和普及，使矩阵运算的重要性愈加显著，应用日益广泛。 这是因为用矩阵理论和方法来解决现代工程技术中的各种问题，不仅表述简洁，便于 进行研究，而且具有适合计算机处理的特点。矩阵运算如矩阵分解、矩阵求逆等，在 简化和解决很多问题上是关键环节。如求解代数线性方程，分析综合线性电路系统， 以及参数优化等，都涉及到了矩阵运算删。 当代许多信号处理、图像处理，以及通信中的操作和运算，都要求有很高的吞吐 量，大部分的操作都要求实时完成，这样对于算法的实现速度就有很高的要求。在数 字信号处理中，有许多信号和图像处理的算法具有运算局部化、计算密集以及大多数 是矩阵运算等特点。因此寻找高速矩阵运算的方法对于很多复杂的算法操作是很有意 义的。 例如在雷达成像和阵列信号处理算法中，很多要求计算权系数的地方都要用到矩 阵求逆。但由于矩阵求逆算法的复杂性等原因，大多数情况下使用一些迭代的算法， 如最小均方算法、最小二乘法来求近似的权系数，这样在运算精度上不够高，收敛速 度不够快。为此，鉴于f p g a 高速和并行的特点，我们开展了基于f p g a 矩阵运算 的设计研究，以有效实现快速的矩阵求逆等运算。 1 2 研究概况 一般来说实现矩阵运算的方法可以有以下几种：在通用处理器( g p p ) 上实现； 用微控制器( m c u ) 实现；在通用数字信号处理器( d s p ) 上实现：用专用芯片( a s i c 或f p g a ) 实现1 9 1 。 在通用的计算机上选择一种语言来编制算法的软件包实现，这种方法由于通用计 算机在体系结构上采用冯诺依曼结构，缺乏专用的硬件乘法器、循环条件判断等硬 件支持，执行速度太慢，不适用于实时系统，而多用于事后处理和仿真研究工作。 微控制器m c u 又称单片机，将计算机的c p u 、r a m 、r o m 、定时器和多种i o 接口集成在一片芯片上，形成芯片级的计算机。其接口性能好，易于实现人机接口， 各种应用程序也很成熟。但由于其结构不是专门为复杂运算设计的，尤其是乘法速度 慢，因此只能用于一些不太复杂的处理，如数字控制等。 采用d s p 实现算法的方法可移植性和灵活性较好，并且d s p 有定点处理器和浮 点处理器，在应用中浮点d s p 能满足大动态范围的要求。目前市场上主流的浮点d s p 芯片包括1 r i 公司的t m s 3 2 0 c 系列和a d 公司的a d s p t s 系列等都具有良好的性能， 1 硕士论文基于f p g a 的矩阵运算实现 目前基于d s p 矩阵运算的系统已有成功的应用。 然而d s p 实现的方法在运算效率、执行速度等方面有时还不能完全满足工程的 需要，特别是对于大量复杂运算的实时应用方面。而采用f p g a 实现方法，有执行效 率好、速度快、集成度高等很多优点。近年来随着f p g a 工艺不断的进步，国内外的 数字信号处理算法研究者和开发者也更多地采用基于f p g a 芯片来实现自己专用的 算法。随着动态可配置技术的发展，用f p g a 来实现众多复杂算法的灵活性已经有了 长足进展。 f p g a ( 现场可编程门阵列) 是一种可由用户编程来实现所需逻辑功能的数字集 成电路【1 。随着e d a ( 电子设计自动化) 技术和微电子技术的进步，f p g a 的时钟 延迟可达到n s 级，结合其并行工作方式，在超高速、实时测控方面有非常广阔的应 用前景；并且f p g a 具有高集成度、高可靠性，几乎可将整个设计系统下载于同一芯 片中，实现所谓片上系统，从而大大缩小其体积，因此得到了广泛的应用。 用f p g a 来实现信号处理算法主要有以下几点优势【l l 】：( 1 ) f p g a 非常适合于各 种算术运算，包括含有大量乘加运算的数字信号处理算法。用f p g a 设计数字信号处 理系统时，可以利用并行结构和算术运算的优点，充分利用硬件资源，并在性能上超 过d s p 器件。在现有的数字信号处理方案中，利用f p g a 中阵列乘法的分步式运算， 可以将数据带宽和通过量提高几个数量级。( 2 ) 相对通用处理器和d s p 处理器而言， f p g a 可以由设计者根据算法的内在结构设计合适的处理阵列，避免前者串行执行指 令的低效。( 3 ) 相对a s i c 而言，采用f p g a 可避免初期巨大的开发投资，同时拥有 微处理器的通用性和灵活性。在算法修改时，短时间内就将新的算法付诸实际。由于 f p g a 仍然属于通用器件，效率要比a s i c 低，但其灵活性的优势仍然可以在很大程 度上弥补其缺点。因此利用f p g a 完成信号处理算法是一种方便、快捷、最具优势的 优化设计方案。 f p g a 元器件最具吸引力的实现数字信号处理算法的功能就是快速进位逻辑的 能力，能够以超过5 0 m h z 的速度实现3 2 位( 非流水线) 的加法。x i l i n x 元器件具有 f p g a 中典型的宽泛的路由选择级，例如x i l i n xx c 4 0 0 0 系列的基本逻辑单元可配置 逻辑模块( c l b ) 具有两个独立的4 输入l 输出的u j t 和快速进位，另外一个3 输 入l 输出的l u t 将两个独立的u j t 连接起来，还有两个触发器。x i l i n x 元器件具有 5 层路由，从c l b 到c l b ，再跨过整个芯片的长线，每一个c l b 都可以用作1 6 x 2 或3 2 x1 位的r a m 或r o m 。由此可见，x i l i n x 的方法拥有更多的局部路由资源而全 局资源则较少，这对信号处理算法的使用是有促进作用的，因为绝大部分数字信号处 理算法都是处理局部数据的l “。 1 3 论文结构安捧 2 硕士论文基于f p g a 的矩阵运算实现 本文主要对矩阵运算中的矩阵相乘累加、三角矩阵求逆以及一般矩阵求逆等运算 的f p g a 实现做了较为深入的探讨，包括算法的选择、设计、仿真和验证，最后将程 序下载到带u s b 接口的f p g a 板卡上运算并将结果送入主机查看验证。 论文结构安排如下： 第一章：介绍了本课题的研究背景和概况，阐述了矩阵运算的意义，实现矩阵运 算的方法比较以及用f p g a 实现矩阵运算的优势。 第二章：有关矩阵运算的一些基本概念的介绍，包括矩阵乘法，矩阵求逆和矩阵 分解方法的比较，讨论了适合于f p g a 实现的矩阵分解求逆方法。 第三章：对矩阵运算中基本运算单元及逻辑单元的v i - i i ) l 设计与仿真进行了介绍， 分析了定点数与浮点数相互转换模块，以及浮点乘加器和除法器模块的设计过程，给 出仿真结果，最后说明了流水线结构中基本逻辑模块的设计方法。 第四章：详细介绍了各种矩阵运算模块的设计实现方法，对矩阵相乘累加、三角 矩阵求逆模块进行了仿真分析，讨论了一般矩阵求逆的实现方法，给出l u 分解和三 角矩阵相乘的硬件实现结构，最后对高阶矩阵运算的分块实现方法作了说明。 第五章：讨论了基于x i l i n x 公司的v l r t e x - 4s x 系列f p g a 硬件平台上的矩阵求 逆运算过程的实现，介绍了设计中所使用的硬件平台，说明了u s b 接口部分程序设 计方法，并对一般的矩阵的求逆过程进行仿真分析，将程序下载到硬件平台上调试， 并将结果通过u s b 送入p c 机，证实了该算法的正确性和可行性。 硕士论文基于f p g a 的矩阵运算实现 2 矩阵运算的基本概念 2 1 引言 矩阵运算是描述许多工程问题中的数学关系所不可缺少的工具，矩阵乘积、逆矩 阵等都是信号处理和系统理论中最基本的数学工具。本章重点介绍这些常用矩阵运算 的基本概念，并讨论适合于f p g a 实现的矩阵运算方法。 2 2 矩阵乘法1 2 l 矩阵的乘法运算包括标量与矩阵的乘法、h a d a m a r d 积、k r o n e c h e r 积和矩阵与矩 阵的乘法，下面分别介绍这几种乘法。 标量与矩阵的乘法简称数乘，是标量与矩阵各元素相乘得到与原矩阵相同维数的 新矩阵，属于矩阵的线性运算。h a d a m a r d 积是按矩阵的对应分量相乘的矩阵乘法。 k r o n e c k e r 积又称为直积，它是信号处理与系统理论中的随机静态分析、随机向量和 随机向量过程分析等的一种基本分析工具。p x q 矩阵a 和m x r l 矩阵b 的k r o n e c k 盯 积记作a o b ，它是一个p m x q n 矩阵，定义为 | a i l b a i q bi a 圆b = i ji l ( 2 2 1 ) l a p l b a p q b j 由此可见，a b 是以a ；b 为子块的分块矩阵，虽然a o b 与b o a 的阶数相同， 但它们并不相等，不满足交换律。 以上三种矩阵的乘法运算都可以看成是标量与矩阵的运算，不涉及乘累加操作， 重点在结果矩阵的排列上，因此用硬件实现较容易，本文中不作重点介绍。 矩阵与矩阵的乘法定义为，设a = 【a i i 】是一个m s 矩阵，b = b i 3 】是一个s x n 矩 阵，那么规定矩阵a 与矩阵b 的乘积记作c = a b ，c 是一个m n 矩阵 厂1 c = 【c i j 】= i a m b ( i - - 1 ，妒，m ；j = 1 ，2 ，n ) ( 2 2 2 ) l k ij 只有当a 的列数等于b 的行数时，a b 才有意义，并且a b 是m n 矩阵，而乘 积a b 的第i 行第j 列元素是矩阵a 的第i 行各元素分别与矩阵b 的第j 列各对应元 素乘积之和。一般矩阵的乘法不满足交换律，因此，矩阵乘法运算一定要注意先后次 序。 2 3 矩阵求逆 4 硕士论文基于f f g a 的矩阵运算实现 2 3 1 逆矩阵 2 1 矩阵求逆算法是矩阵乘法的逆运算，在线性代数中，逆矩阵的定义是：对于i i 阶方阵a ，如果存在n 阶方阵b ，使a b = b a = e ，则称方阵a 是可逆的，并把b 称为 a 的逆矩阵，其中e 为n 阶单位阵。如果方阵a 可逆，则其逆矩阵是唯一的。 判断矩阵是否可逆的方法可以看它是否为奇异矩阵。一线性变换或正方矩阵a ， 若它只对零输入产生零输出，或者说a 的行列式值不为零，贝u 称为非奇异的，否则 是奇异的。如果一个矩阵非奇异，则必存在逆矩阵，奇异矩阵不存在逆矩阵。矩阵奇 异性的另外一种常用表示是行列式，一个n x n 矩阵的行列式用( n - 1 ) x ( n - 1 ) 矩阵的行 列式定义，即 d c t ( a ) = ( 一1 ) a l jd e t ( a i j ) ( 2 3 1 ) j - l 其中a l ；是除去4 的第1 行和第j 列后得到的( n 1 ) ( n 一1 ) 矩阵。 一个n x n 矩阵彳是奇异的，当且仅当d e t ( a ) = 0 。因此我们可以得到如式( 2 3 2 ) 所示的判断关系。 d e t ( a ) o a 是非奇异的a _ 可逆 ( 2 3 2 ) 2 3 2 一般矩阵求逆的方法 求逆矩阵的方法有很多，如线性代数中介绍的伴随矩阵法、初等变换法、分块矩 阵法等，还有一些工程中常用的求逆方法，如g a u s s - j o r d a n 消去法求逆阵、矩阵分解 求逆等。 l 、伴随矩阵法( 公式法) 由线性代数中n 阶方阵a 的逆矩阵公式，如式( 2 3 3 ) ，可以直接求出矩阵的逆 矩阵。 1 a 一= 击a ( 2 3 3 ) j a i 其中a = a “a 1 2 a 2 i a 2 2 a n la n 2 a l 。 a 2 n a 。 为矩阵a 的伴随矩阵，a 。是除去a 的第i 行和 第j 列后得到的( n 1 ) ( n - 1 ) 矩阵的行列式。 2 、初等变换法 由n 阶方阵和n 阶单位阵构成n x 2 n 增广矩阵( a ；e ) ，然后对此矩阵作仅限于行 5 硕士论文基于f p g a 的矩阵运算实现 的初等变换，使子块a 化为e ，同时子块e 则转化为a - 1 3 、分块矩阵法( 降阶法) 1 2 令n x n 矩阵y 被子矩阵a ，b ，c ，d 分块为 y = ( a ：) 式中a 是m x m 矩阵，b 是( n m ) x ( n m ) 矩阵，c 是( n m ) x m 矩阵，而d 是 m x ( n 一神矩阵。矩阵y 的逆矩阵由式( 2 4 5 ) 给出 v - lf a - + a - l d q c a - 1 - a - 1 d a - 11 l 一q c a 1- 1 j i a - l- a - i d b - l i g 3 5 ) l b - 1 c a - 1b - 1 + b - 1 c a - 1 d b - 1j 假定a 和a = b c a 。d 的逆矩阵存在，或b 和a = a d b - 1 c 的逆矩阵存在。式 中，矩阵和人称为矩阵a 的s c h u r 补。 4 、g a u s s - j o r d a n 消去法脚】 对于n x n 矩阵，g a u s s j o r d a n 法求逆的方法是：对k 从1 到n 作如下操作，从第 k 行、第k 列开始的右下角子阵中选取绝对值最大的元素，并记住此元素所在的行号 和列号，再通过行交换和列交换将它交换到主元素位置上，这一步称为全选主元。因 此g a u s s - j o r d a n 消去法又称为全选主元法。最后根据在全选主元过程中所记录的行、 列交换的信息进行恢复，恢复时应注意在全选主元过程中先交换的行( 列) 后进行恢 复，原来的行( 列) 交换用列( 行) 交换来恢复。 5 、矩阵分解法【2 1 】 把矩阵分解为较简单的一些矩阵的乘积，如三角矩阵或酉矩阵这种具有某种特 性，较容易得到逆矩阵的特殊矩阵，对分解得到的矩阵求逆后相乘，得到原矩阵的逆 矩阵。若将矩阵a 作三角分解成两个三角矩阵，即a = l u ，则a 一= u 。1 l - 1 。 上述各种矩阵求逆的方法中，伴随矩阵法中需要求大量的行列式，每个行列式都 几乎要计算到所有的矩阵元素，计算量大，对存储空间的需求也较大；初等变换法需 要增加矩阵的阶数、且涉及复杂的矩阵变换操作，不利于硬件实现；分块矩阵法主要 用于阶数很大的矩阵变换为阶数较低的矩阵，以避免阶数大时计算量和资源需求都呈 平方增长的缺陷；g a u s s j o r d a n 消去法由于运用了大量的排序操作，计算量较大， 实现过程不宜采用并行计算，运算时间较长，较为适合于软件实现；矩阵分解实现求 逆则客服了上述方法的缺点，三角矩阵求逆硬件实现简单，分解和求逆都可以采用并 行结构实现，运算速度快，且采用分块的方法可以实现高阶的矩阵求逆。因此本文主 要讨论三角分解矩阵求逆方法的f p g a 实现。 6 硕士论文基于f p g a 的矩阵运算实现 2 3 3 三角矩阵求逆的方法 采用初等变换法可以得到三角矩阵逆矩阵的一般公式。对于n 阶上三角矩阵u ， 求逆过程中与n 阶单位阵e 构成增广矩阵如式( 2 3 6 ) ( u l e ) - - u l l b 1 2 0 u n oo u l i l u 2 i ；l u 。i ( 2 3 6 ) 对第n 列的操作步骤如下：从第n 行开始，先除以u 。；然后乘以一u 。加到第i 行( i = n l ，n 一2 ，1 ) 消去u 矩阵中第1 1 列除了主对角线外的所有元素 u 。u ：n ，“，u 。，依此类推，直至处理完第一列，此时，u 矩阵变换成了单位矩阵， 而单位矩阵就变换为u 的逆阵u 一。 设v = u ，综合上面的步骤可给出u 的逆矩阵v 的计算公式如式( 2 3 7 ) 所示。 v “= 亩i ( i = 1 ，2 ，。，n ) v q u m v 日2 一午( i = n - l , n - 2 , - 1 ；j = 1 + 1 ，n ) 对于下三角阵，可以对其作如式( 2 3 8 ) 的操作。 儿( ( l t ) 7 卜( ( l t ) 1 ) t 先计算下三角阵l 的转置、上三角阵l t 的逆，然后再转置求出p 。 2 4 矩阵分解 ( 2 3 7 ) ( 2 t 3 。8 ) 所谓矩阵的分解就是通过线性变换，将某个给定或已知的矩阵分解为二或三个矩 阵标准型的乘积( 个别情况下分解为两个矩阵标准型之和) 。 根据矩阵分解后得到的矩阵标准型以及是对单个矩阵还是两个矩阵组成的矩阵 束或矩阵对进行分解来划分矩阵的分解类别。单个矩阵的分解包括对角化分解、三角 化分解、三角一对角化分解和三对角化分解等。矩阵束的分解涉及两个矩阵的同时分 解，主要用于求解矩阵束广义特征值分解( g e v d ) 问题的q z 方法。主要有s c h u r 分解，实现广义s c h u r 分解需要先使用h e s s e n b c r g 三对角化分解【2 j 。 本节重点讨论单个矩阵的三角化分解，即将矩阵分解为正交矩阵与三角矩阵之 积，或分解为两个三角矩阵之积，主要介绍的是c h o l e s l c y 分解、q r 分解和l u 分解。 2 4 1c h o l e s k y 分解1 2 1 7 硕士论文 基于f p g a 的矩阵运算实现 设a = ( a 日) e r “。是对称正定矩阵，a f g g 7 称为c h o l e s k y 分解，其中g r 一是 一个具有正的对角线元素的下三角矩阵： g = g i i 9 2 19 1 2 g n ig 量2 0 g 。 ( 2 4 1 ) 所求的下三角矩阵g 称为c h o l e s k y 三角，可视为矩阵a 的“平方根”。一个非 奇异矩阵a 的逆矩阵a 。可以通过c h o l e s k y 分解求得，即a - 1 = g _ f g ，其中 g - t = f g 7 r 1 。 比较a = g g 7 两边的下三角部分元素，易得 f a l l = l g , , 1 2 ， l a i l2 9 i l g l ：( i 。2 ，n ) - ( 2 4 2 ) i a 址- - i g 。1 2 + l g 。：1 2 + + i g 。1 2 ( k = 2 ，3 ，n ) ， i a 址= g i l 瓦l + g i 2 虱2 + + g 址吾娃( i = k + l ，n ；k = 2 , 3 ，n ) 从而得到求对称正定矩阵a 的c h o l e s k y 分解的紧凑计算格式： g “= a i _ _ l ( i = 2 ，3 ，n ) ， g n 酽、鹰k - 面i 2 ( k - - 2 , 3 , - - , n ) ， q a 3 y t f f i l 驴亡( 善嘱 ( i = 呲”川k 喇，神 2 4 2l u 分解1 u 设a c ”，如果存在下三角矩阵l c 一和上三角矩阵r c “，使得a = l r ， 则称a 可以作三角分解。a 可以作三角分解的充分必要条件是a k 0 ( k = l ，2 ，n 1 ) ， 其中a 。= d e t a 。为a 的k 阶顺序主子式，a 。为a 的k 阶顺序主子阵。因此不是每个 可逆矩阵都可以作三角分解。另外，当矩阵a 的秩为r 的时候，当a 的前r 个顺序主 子式不为零时，a 可以作三角分解。 由于三解分解不是唯一的，为了规范化三角分解，定义当a c ”时，如果a 可 分解成a = l r ，其中l 是对角元素为1 的下三角矩阵( 单位下三角矩阵) ，r 是上三 角矩阵，则称之为a 的d o o l i t t l e 分解，也就是上述的l u 分解；如果l 是下三角矩 阵，r 是对角元素为1 的上三角矩阵( 单位上三角矩阵) ，则称之为a 的c r o u t 分解； 8 硕士论文基于f p g a 的矩阵运算宴现 如果a 可分解为a = l d r ，其中l 是单位下三角矩阵，d 是对角矩阵，r 是单位上三 角矩阵，则称之为a 的l d r 分解。 综上，l u 分解就是将矩阵a 分解成一个单位下三角矩阵l 和一个上三角矩阵u 的乘积，即a = l u ，如式( 2 4 4 ) 示。由此可以得到l u 分解的计算方法。 于是 a l ia 1 2 a l nl a ：- a n a ：| - ；j【 a 。a 。：a j a l j - - - - u l j 0 = 1 ，2 ，n ) ； a i l = l i l u l l ( i = 2 , 3 ，n ) ； a q = 1 h u 日+ u q ( j = k ，k + l ，n ；k = 2 ，3 ，n ) ； ( 2 4 5 ) t = l k - i a m = k t u d c + k u h ( i = k + 1 ，k + 2 ，n ；k = 2 ，n ) 由上式可得到求a 的l u 分解的紧凑计算格式为 u l j2 a l j ( j5 l ，2 ，n 】； 1 i 1 ：旦生( i ：2 州3 ，n ) ； u l l u ：a q 一k - t l h u 日0 ：k ，k + l ，n ；k ：2 ，3 ，n ) 2 4 6 k = 亡( 尊u 。) 雌扎吼，n 嗡一m 具体计算时可按图2 1 所示一框一框地逐步进行。对同一框中的元素，u 。必须 在计算k 之前先算出，其余元素的计算先后没有影响，可以并行计算，并且在计算出 u i j 或k 后，就不再使用了 毡iq 2 1 3 巩 第l 框 第2 框 第n 框 图2 1l u 三角分解的紧凑计算格式示意图 2 4 3q r 分解1 j 设a e c ，如果存在n 阶酉矩阵q 和n 阶上三角矩阵r ，使得a = q r ，则称 9 442 、l，； 伦 盟 u u u ，j。l 、1liii， b l k ；k 硕士论文基于f p g a 的矩阵运算实现 之为a 的q r 分解或酉- 三角分解。任意a c 都可以作q r 分解。当a 是可逆矩 阵时，a 的q r 分解是惟一的，其中r c ：”是具有正对角元的上三角矩阵。 下面简单介绍三种主要的q r 分解算法，分别是采用h o u s e h o l d e r 变换的q r 分 解、采用g i v e n s 旋转的q r 分解和s c h m i d t 正交化的方法。 l 、采用h o u s e h o l d e r 变换 将矩阵a 按歹q 分块为a = ( a l ，a 2 ，a 。) ，存在h o u s e h o l d e r 矩阵h = i - 2 u u “( 其 中u c “是单位向量，即u “u = 1 ) ，使得h a ，= a l e ，于是 h 1 a = ( h l a l ：h 1 a 2 ，h l a 。) = a t + 1 7 艮。i 包4 力 0 其中b o 1 c ”1 m “，再将e l 作上述h o u s e h o l d e r 变换，继续到第n - 1 步得 f a 一事1 h 。- 1 h 2 h 1 a = l | - r ( 2 4 8 ) 【a 。j 其中h k ( k = l ，2 ，行一1 ) 都是n 阶h o u s e h o l d e r 矩阵。由于h k 均是自逆矩阵，则 有a = h h 2 h 。r = q r ，这里q = h l h 2 h 。一l 是酉矩阵，r 是上三角矩阵，于是得 到a 的q r 分解。 2 、采用g i v e n s 变换 将矩阵a 按列分块为a = “，a ：，a 。) ，存在n 阶c f i v c n s 矩阵t 1 2 ，王。，使得 玉。玉：a 。- - b 。i l ：e 1 ，于是 小私= ( 1 i 四2o ：i ( b k 蚪，呦旺4 m 对于其第二列再作上述g i v e n s 变换，最后得 t 。t 2 。k i 。玉2 a = r ( 2 4 1 0 ) 于是a = 碟碟t 昙碟堪，。r = q r ，其中q = 碟磺磕堪，是酉 矩阵，r 是上三角矩阵，得到a 的q r 分解。 3 、采用s c h m i d t 正交化方法求可逆矩阵的q r 分解 将矩阵a 按列分块为a = ( a i ，a 2 ，a 。) ，由于a 可逆，所以a l ，a 2 ，a 。线性无关， 用s c h m i d t 正交化法将其正交化： i o 硕士论文基于f l e a 的矩阵运算实现 仨一k ， 其中气= 石( a 孟i , p 万j ) ，再将p k ( k = l ，2 ，n ) 单位化得 q t 。奇( k _ l 2 ，n ) 故 ( 2 4 1 1 ) ( 2 4 1 2 ) a ，= p l - - l i p 。i i ：q 。， a ：= 五p 。+ p 2 = 五。n p 。0 ：q 。+ 9 p ：0 ：q 2 ， ( 2 4 1 3 ) a 。= 五l p l + + 九正一i p - i + p 。 = 九1 i p 。l i ：q 。+ + 五，。i i p 。一。l ：q 。一。+ l i p 。l i ：q 。， a = ( a l , a 2 ，a 。) = ( q 1 ，q 2 ，q 。) i i p l | i ：五l l | p 1 i i ：钆| l p l l i ： | i p 2 l l ：五：i l p 2 | l i i l p n i i ： = q g ( 2 4 1 4 ) 其中q 是酉矩阵，r 是具有正对角元的上三角矩阵，由此得到可逆矩阵4 的q p 分解。 2 4 4 三种分解方式的比较 上述三种矩阵分解都是将单个矩阵分解成至少有一个三角矩阵的三角化分解方 式，下面分别从其分解后矩阵的特点、适用矩阵范围和运算的复杂程度对这三种分解 方式进行比较。 1 、分解后的矩阵特点：c h o l e s k y 分解后的矩阵是两个三角矩阵，其中上三角矩 阵为下三角矩阵的转置，即分解为一个下三角阵和它的转置的乘积；i 耵分解后是一 个单位下三角矩阵和一个上三角矩阵；q r 分解后是一个酉矩阵和一个具有正对角元 的上三角矩阵。 2 、分解适用矩阵范围：c h o l e s k y 分解是对h e r m i t e 正定矩阵的三角分解；l u 分解是对所有顺序主子式不为零的矩阵的三角分解；q r 分解对任意a c “都可以 进行。因此q r 分解在信号处理中有着最广泛的应用，特别是自适应信号处理中的参 堡主丝苎茎主! 堡垒箜堑堕墨蔓壅墨 l 数估计问题可以通过q r 分解进行分析。另两种分解则是对特殊矩阵进行的分解，在 应用上有局限性。 3 、运算的复杂程度：q r 分解的多种分解计算方式中，g i v e n s 变换里g i v e n s 矩 阵中c 和s 的选择必须为实数且c 2 + s 2 = l ，因此构造g i v e n s 矩阵过程较复杂。 h o u s e h o l d r 变换法和s c h m i d t 法q r 分解在数学上是等价的，都有向量的单位化，运 算量较大。c h o l e s k y 分解和l u 分解的计算格式较简单，程序化实现容易，计算量较 q r 分解小。 综上的比较，c h o l e s k y 分解虽然计算量不大，但适用矩阵范围太小；q r 分解适 用于任意矩阵，但由于计算过程复杂，不适合于硬件实现；l u 分解虽然不适用于任 意矩阵，但由于大多数可逆矩阵都满足所有顺序主子式不为零的条件，因此l u 分解 在矩阵求逆中仍然得到广泛的运用。本文中基于f p g a 的矩阵求逆中矩阵分解采用的 是l u 分解的算法。 2 5 本章小节 本章主要介绍了矩阵运算的基本概念，首先给出矩阵的几种乘法形式，然后比较 矩阵求逆的各种计算方法，讨论适合于f p g a 实现的矩阵分解求逆方法，在此基础上 说明了三角矩阵求逆的计算格式。最后分别介绍了几种矩阵三角化分解形式，包括 c h o l e s k y 分解、l u 分解和q r 分解，在分析了这三种分解的差异后，选择l u 分解作 为矩阵求逆中矩阵分解的方法。 硕士论文基于f p g a 的矩阵运算实现 3 矩阵运算中基本运算单元和逻辑模块的设计 3 1 引言 在矩阵运算中，涉及大量乘累加操作，由此带来的“位增长”率使每一级运算的 最大值可能会逐级加倍( 但是这样的增长率并不是每一级都会持续发生) i 阍。因此如 果不精心地规划设计，这些值就会溢出，结果会因为精度不够而无法使用。浮点运算 则克服了这些问题，因为它具有极大的动态范围而不会发生溢出，还因为它可以自动 地缩放，无需担心数值会溢出或者当算法需要缩放时会丢失尾数位的问题。因此考虑 采用自定义的浮点数进行基本的运算单元设计。 用v h b l 实现矩阵运算的模块，其中最基本的运算单元包括乘加器和除法器，逻 辑模块包括数据选择器、数据分配器以及用于流水线操作的锁存器等。本章主要讨论 浮点运算中定点与浮点之间转换模块的设计，以及浮点乘加器和除法器的设计过程， 并介绍了流水线中数据选择分配器的设计方法。 3 2 浮点运算数制 浮点数制可以在更大的动态范围内提供更高的分辨率，由于当前f p g a 具备大量 的门数量，使浮点算法的设计成为可行选择。本节主要介绍浮点数格式，以及浮点数 与定点数相互转换的模块设计。 3 2 1 浮点数格式l a l 标准浮点数字长由一个符号位s 、指数e 和无符号( 小数) 的规格化尾数m 构成。 其格式如图3 1 所示。 广1 t l 符号s l 指数el 尾数n li ijj1j 图3 1 标准浮点数格式 最高位为符号位，0 表示正数，l 表示负数。其后是用原码表示的指数e ，若指 数为e 位，则表示范围为：0 2 e 一1 ，指数的偏移量b i a s 由式( 3 2 1 ) 表示。 b i a s = e “一l( 3 2 1 ) 因此指数的表示范围为一2 “1 1 ) 2 “。最后为原码表示的尾数m 部分，为了使 尾数部分能表示更多一位的有效值，采用隐含尾数最高位l 的方法，假设尾数用m 位表示，则实际上尾数一共有m + i 位。由此可见，浮点数字长的代数表达式如式 ( 3 2 2 ) 所示。 x = ( 一1 ) 5 1 r e x 2 。妇 ( 3 2 2 ) 本文考虑的1 8 位浮点数采用( 1 , 7 ，1 0 ) 浮点格式，即是由1 位符号位，e = 7 位指 1 3 硕士论文基于f l e a 的矩阵运算实现 数宽度和m = 1 0 位的尾数( 不包括隐藏的1 ) 组成的浮点数表达式。由式( 3 2 1 ) 得 到此格式浮点数的偏移量是b i a s = e “一1 = 6 3 。它所能表示的实数范围为： _ 2 “( 2 2 1 0 ) ：2 “( 2 2 - 1 0 ) ，能表示的最小实数是- 2 。，小于- 2 。的实数按0 处理，动 态范围为8 2 5 d b ，而1 8 位的有符号定点数的动态范围仅为1 0 2 d b 。 假设一个1 8 位( 1 , 7 。1 0 ) 格式的浮点数1 1 0 0 0 1 0 1 1 1 0 0 1 0 1 1 0 1 要表示成十进制数格 式，则进行如下转换：首先符号位s = 1 ，则表示该数为负；指数位为( 1 0 0 0 1 0 1 ) 2 = ( 6 9 ) i o ， b i a s = ( 6 巩o ，因此式( 3 2 2 ) 中指数为( 6 9 - 6 3 ) 1 0 = 6 l o ；尾数位为m - - ( 1 1 0 0 1 0 1 1 0 1 ) ，。 根据式( 3 2 2 ) 可以将该浮点数表示为式( 3 2 3 ) 所示。 x = ( 一1 ) 8 x 1 m x 2 。- b i - = ( 一1 ) 1x 1 1 1 0 0 1 0 1 1 0 1 x 2 6 ( 3 2 3 ) = ( 一111 0 0 1 0 11 0 1 ) 2 = ( - 11 4 8 1 2 5 ) i o 可见该浮点数所表示的十进制数为一11 4 8 1 2 5 。下面将讨论的定点与浮点之间转 换就是基于这种方式的转换。 另外，该浮点数算法还定义了一些其它有用的特殊数。如指数e = e 一= 1 1 ：与 0 尾数m = 0 组合是为。保留的。这样在运算中如果出现溢出，或是除数为0 的- i t 算， 结果都可以用c o 表示。0 是用0 指数e = e 曲= o 和0 尾数m = o 编码。 下面讨论定点数与浮点数转换的实现方法，当采用不同的定点或浮点格式时，将 此方法作相应的改动即可。 3 2 2 定点一浮点转换模块设计 由于设计中运算过程采用浮点数制，而输入输出数据采用的是定点数制，因此要 考虑定点数与浮点数的转换模块实现。这里所使用的定点数具有l 位符号位，1 5 位 整数位和1 6 位小数位的3 2 位定点数制。考虑1 8 位浮点数可表示数的范围大于3 2 位定点数的表示范围，定点数中的最大值和最小值在浮点数中都可以表示，因此转换 中不考虑溢出的情况。 假设一个3 2 位( 1 , 1 5 ，1 6 ) 格式定点数1 0 0 1ll 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 11 0 0 1 0 0 1 0 0 ，该 定点数表示的十进制数由式( 3 2 4 ) 计算，由于是负数，首先要将该定点数转换成补 码形式表示。 1 0 0 1l l1 0 0 0 1 0 0 l 0 0 0 1 0 0 0 1 11 0 0 1 0 0 1 0 0 = - 1 1 0 0 0 0 1 1 1 0 1 1 0 1 1 1 0 1 1 1 0 0 0 1 1 0 1 1 1 0 0(