标签: 矩阵乘法

由于笔者马上就要毕业了，想到公司里面锻炼一下，最近的一个多月都在公司实习。而公司也没有跟我客气，上来就叫我做算法，当然是一些已经成熟的技术，可查找的资料比较多。比如：矩阵求逆。由于公司的保密性质，我以下写的内容都是可以公开，可以查到的，而我的工作是将这些资料做一些整合，加上我的理解，让大家更容易读懂。 刚接到老大说要做矩阵求逆，内心是没什么波澜的，因为没接触过，不懂而无畏（O(∩_∩)O哈哈~） 。接下来上网查一些资料，弄清楚矩阵求逆的整体架构，包括哪些部分（包括矩阵分解、三角矩阵求逆、矩阵乘法三个部分），心中有个整体的把握之后，开始着手矩阵乘法（选它是因为矩阵乘法的资料最多，觉得是个比较好的突破口）。看了一天资料，第二天开始设计结构（结构设计是最难的，这个是论文里面设计好的，好在矩阵运算的大部分的结构设计都可以找的到），编程加仿真，ok。元旦三天假期回来开始矩阵分解，并不像矩阵乘法那样顺利，两三天过去了，仍然没有进展，内心是有点着急，开始求助于各个FPGA群，希望能从经验丰富的大牛们那里得到些灵感，而只得到了“你还太年轻了”的回答（原因是对方觉得一个月的时间绝对做不出，想当年他工作的时候还做了一个多月，而基于保密性，不能提供给我帮助）。不知道大家是不是这种，你越是说不可能，我就越要证明给你看。周末跑去学校借了几本书，从头开始研究脉动阵列（脉动阵列是矩阵运算的关键，以后会讲到），论文看了十几篇，终于把矩阵分解做了出来。由于这段时间的研究打下的基础，没多久，整个矩阵求逆也顺利的得到了。。内心还是有些小兴奋的。 总结一下，我想说的是1.“基础一个要打牢，做事要踏踏实实的准没错”2.“不逼自己一把，你不会知道你有多大的潜力”。上面那些讲的是我的学习经历，大家可以看看，也可以当废话略过。 接下来进入正题，浮点转定点是个比较基础的知识点吧，所以作为开篇，简单的举几个小例子，通过例子，相信大家都能掌握它。 简单说明一下，浮点包括 符号位| 指数位|小数位 。浮点的类型包括 单精度浮点数|双精度浮点数 。这里用到的是单精度浮点数。单精度浮点数： 1 位符号位， 8 位指数位， 23 位尾数位。也有说是24位尾数位，这里笔者认为这样划分，更便于说明（我的地盘听我的，嘿嘿）。 浮点转定点的步骤如下：a）将浮点数划分 符号位| 指数位|小数位; b)计算指数位与偏差位的值;单精度浮点数的偏差值固定为127. c)计算并得到定点数。 看例子 ： 例1： A =01000000.010000000000000000000000 划分如下： 01000000010000000000000000000000 符号位 =0 ，即为正 指数位 =128 ，偏差值 =127 ，指数位 - 偏差位 =1 尾数位为 1 意味着实际尾数为 1.1 ，包括小数点前面隐藏的 1 所以实际的实数为 +1.1*2^1=(2^1+2^-1*2^1)=3 再举个例子，例子 2 ： B =01000000.101000000000000000000000 划分如下： 01000000101000000000000000000000 符号位 =0 ，即为正 指数位 =129 ，偏差值 =127 ，指数位 - 偏差位 =2 尾数位为 01 意味着实际尾数为 1.01 ，包括小数点前面隐藏的 1 所以实际的实数为 +1.01*2^2=(2^2+2^-2*2^2)=5 到了这里，相信大家基本掌握了单精度浮点转定点的方法，接下来的两篇计划介绍矩阵乘法以及LU矩阵分解。由于保密性，里面的内容都是可以公开的。 连接到下一篇《 矩阵系列之矩阵乘法 》 下载视频

上一篇说到一个基本的小知识点浮点到定点的转换，这一篇来说说矩阵乘法。矩阵乘法和下一篇要说的矩阵LU分解是矩阵求逆的重要组成部分，所以就算大家不需要做矩阵求逆，对其先有个整体的认识也是好的。（矩阵求逆的整体框图还是很好理解的 ，甚至你只要瞟一眼图就好）。 1 矩阵求逆的整体框图 矩阵求逆的步骤如下： 1. 原始矩阵 A 通过 LU 分解为上三角矩阵 L 与单位下三角矩阵 U 。 2. 分别通过三角矩阵的求逆运算得到 L 逆和 U 逆。 3. 最后通过矩阵乘法得到 A 逆。 显而易见，矩阵求逆由如下三个部分组成： 1. LU 矩阵分解 2. 三角矩阵求逆 3. 矩阵乘法 2 脉动阵列介绍 做矩阵运算，脉动阵列绝对是个神器，笔者大部分的时间是花在脉动阵列上面的，脉动阵列主要完成算法到结构的映射，而这种结构的设计往往是最难的， 值得一提的是，大部分矩阵运算的脉动阵列都已经被设计出，有的甚至不止一种结构。下面简单介绍脉动阵列，如果感兴趣的可以深入研究。 脉动阵列：多个相同的处理单元（ PE ），按照一定的规则组成的网络，成为脉动阵列。脉动阵列可以是一维线性、二维三角形、二维矩形、二维六边形等等。 特点： 1. 每一个节点，即 PE ，都是相同的，（个别也可以不同）。 2. 每个 PE 只与其邻近的 PE 进行通信，也就是说 PE 之间的通信具有局部性，而且通信是规则的。这点非常重要！！！ 3. 每个 PE 都有其局部的储存器，也就是 PE 的某些边带有延时，延时在硬件上对应于寄存器。 由于以上特点，造成 PE 之间的高度流水化、规则化，因此系统吞吐率非常大且易于 VLSI 的实现。 3 矩阵乘法模块 3.1 矩阵乘法的数学表达式 从公式上可以看出，计算后的矩阵的每个单元（ i ， j ）是矩阵 A 的第 i 行矩阵与矩阵 B 的第 j 列对应元素的乘积和。 可以看出，矩阵乘法的计算规则化，数据间的通讯具有局部化的特点，因此很适合用脉动阵列来实现。 3.2 乘法脉动阵列的设计 根据算法的特点，设计二维方形的脉动阵列，如下图所示：从上到下，从左到右分别流经矩阵 A 、 B 的数据。处理单元 PE 则保存结果矩阵 C ，处理单元以矩阵 C(i,j) 命名，初始值均为 0 。设计好的脉动这列图与以及每个时刻详细的步骤如下： 乘法脉动阵列图 注： * 号代表延时一拍 下面给出计算数据流向与结果的示意图 1. 时刻 1 数据 a11 、 b11 进入 PE 单元，计算并更新 PE ，将数据保存在 PE 中。 2. 时刻 2 数据 a11 继续向右流入 C12 ，与向下流入的 b22 做运算将结果更新保存在 C12 。数据 a21 流入处理单元 C21 ， b11 继续向下流动，与 a21 做同样的运算并更新保存 C21 。而此时 C11 保存的是 a11*b11+a12*b21 的结果。 依次向下向右流动，经过 3*n-2 个时钟周期得到结果矩阵 C ，相比串行执行的时间复杂度 (n^3) ，阶数越多，时间缩短的越多。值得说明的是，该阵列中所有处理单元所完成的功能相同，均为乘加运算。 3.3模块的 测试 测试模块包括：矩阵生成模块，定点转浮点模块，矩阵乘积模块（核心模块），浮点转定点。通过matlab仿真与modelsim的结果对比，验证其正确性。 上一篇《 矩阵系列之开篇浮点转定点 》

脉动阵列（ Systolic Array ）计算矩阵乘法（ Array Multiplication ） 下一个目标是实现流水线输出，提升硬件资源的利用率。 脉动阵列 (Systolic Array) ：数据流同步流过相邻的二维阵列单元的处理器结构，一般不同方向流过不同数据。 结构： 矩阵计算： C 语言描述每个输出矩阵中的值： For I = 1 to N      For J = 1 to N           For K = 1 to N                C = C + A * B ; 运用 N x N processing units ，输入数据呈批次输入： 二维不同数据在同一时钟下依次输入每个处理单元，而后完成乘法并存在其寄存器中。     其中每个 PE （处理单元）结构如下： 是一个乘加单元   c=c+(a*b)   例子：计算两个３×３矩阵的乘积 结构：     在 CLK 驱动下的每一个步骤如下 : Clk1:   Clk2: Clk3: Clk4:          Clk5:     Clk6: Clk7: Clk8: 输出 功能仿真图： 在 start 上升沿到来后的第一个 CLK 上升沿开始计数 Count_start 高电平期间 Cout =1 时，准备 a11 和 b11; Cout =2 时，将数据打入寄存器，并计数出 a11*b11; Cout =3 时，计数 a11*b11+a12*b21 Cout =4 时，计数 a11*b11+a12*b21+a13*b31 Cout =5 时，用寄存器打一拍输出 Y11 。 其他类似。 时序仿真图： 连续运算，中间忘了将乘加单元寄存器清零的情况，功能仿真： 每次计算出结果后清零寄存器，修改后的功能仿真图： 数据在送入运算单元之前，采用寄存器打一拍，功能仿真图：       状态机便于实现控制。 状态机控制：功能仿真 时序仿真图：  

二维流水线结构矩阵乘法（ Array Multiplication ） 上一篇文中建立了矩阵乘法运算的数据路径，从仿真结构中可以看出整个计算方案的可行性，但是存在一个问题，就是硬件运算单元的资源利用率不高。这是什么意思呢？就是说，在每次计算两个 3*3 矩阵的乘法之前，需要将整个运算单元中的每个寄存器都清零，但是 9 个输出结果不是在同一个时刻输出，有先有后。当最先出结果的 P11 计算出第一个结果之后，它就不再输出新值了。其实这时硬件电路是存在的，而且是在不断计算出新值，只是这些值不是我们需要的有效值。那么如果将这些硬件资源完全被利用起来，让它们在最短的时间间隔里都有有效正确数据输出呢。 其实要说的就是流水线设计思想，我们只需要当 P11 单元计算出新的值，下个 CLK 将其计算结果输出（只要有另一个机制接收这些值），然后将其清零（如果不清零，那个会累加了上次的计算结果），再然后就可以将下个要计算的两个 3*3 矩阵对应位置上面的新值输入给 P11 单元，让其计算两个数的乘积。 其他单元类似，思想就是在输出结果后立即做清零和赋新值，不必等整个计算单元都出结果之后才开始进行下个矩阵乘法的计算。 那么从大局来看，就出现个这样的情况：在上个矩阵计算还没有完全计算出 9 个值的时候（右下角部分还在计算），左上角的单元就已经开始下个矩阵乘法的计算。关于输出就是，每个 CLK 都有有效的数据输出，产生了流水输出。而这些输出是从这整个二维矩阵计算单元的某些地方同时输出的，就像水流即向下流也向右流，只是在这里水流对应了数据输出。 回顾一下这样做的目的是什么：为了增加整个运算单元的利用率。做了改进，我们可以看出，每个 CLK 到来时，每个计算单元中的乘法和加法器的运算都是有效的计算。我们一定要记住 FPGA 做运算和 CPU 做计算的不同点。 FPGA 做运算，那些设计好的运算单元，不管你对这些运算单元操作或者不操作，它们都已经以硬件电路的形式存在了，也就是说，每个 CLK 到来时，它们都进行了一次计算。那么如何有效的利用这些运算单元，就是不要让它们做无效的计算，让它们每一次的计算都对我们来说是有效的，是一个对整个单元输出有影响的计算。而想实现这样的功能流水线设计方法是必须要采用的。只是有时候单向流水（一维）更容易被我们所理解。其实这二维运算单元也是如此，也可以产生流水。 相对于不采用流水线操作，数据路径不必做太多修改，主要是控制单元更为复杂，要考虑何时对某个单元清零，何时再对其赋值，何时把数据读走等。这些需要比较精准的调度。有人会想，那么这样是不是需要更多状态机状态？答案是肯定的。那么是不是需要更多的硬件资源的开销呢？答案也是肯定的。那么这种流水线的优势在哪呢？其实是速度上面的优势，有时候在不能提高整个系统的工作频率的情况下，再消耗一点资源来提高其它部分（运算单元）的利用率也是一个很好的提高运算速度的办法。（注：这并不等同于面积与速度的法则）其实这也是一个比较好的，硬件算法优化方案。我们可以将其归纳为以优化算法路径的算法优化方法。 下面是仿真结果：   矩阵还是之前的矩阵。可以看出每个 CLK 都有 2-3 个有效输出。 增加了 Valid 信号，当其为高电平时，说明对应运算单元输出是有效值。其实在算法成熟之后，这些 Valid 信号是可以撤销的，因为我们完全知道了，输出的规律，只需要在特定的时间读走数据即可。而且我们可以将这些数据合并到几根连续有效的总线上面。让每个总线上面每个 CLK 都是有效值。（笔者将按照此方法，继续优化输出总线） 下面是对输入最大值时的输出仿真，主要是看数据会不会益处： 测试 8bits 输入最大值 255 对应的输出值： 时序仿真： 255*255+255*255+255*255=195075 255*254+255*254+255*254=194310