标签: lu分解

由于笔者马上就要毕业了，想到公司里面锻炼一下，最近的一个多月都在公司实习。而公司也没有跟我客气，上来就叫我做算法，当然是一些已经成熟的技术，可查找的资料比较多。比如：矩阵求逆。由于公司的保密性质，我以下写的内容都是可以公开，可以查到的，而我的工作是将这些资料做一些整合，加上我的理解，让大家更容易读懂。 刚接到老大说要做矩阵求逆，内心是没什么波澜的，因为没接触过，不懂而无畏（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矩阵分解。由于保密性，里面的内容都是可以公开的。 连接到下一篇《 矩阵系列之矩阵乘法 》 下载视频

1 矩阵 LU 分解模块 1.1 LU 分解数学表达 首先要明确的是，矩阵的 LU 分解是有局限性的，即 LU 分解只针对非奇异矩阵。那么什么是非奇异矩阵呢？即各阶顺序主子式不为零。 （1） 高斯消去法 LU 分解的思想来源于高斯消去法，拿方阵为例（因为本项目中要处理的就是方阵）。将一个 n*n 的方阵 A ，通过左乘一系列消去矩阵（笔者自己起的名字，便于理解）。使得 (((( n 个 L1*(((L2*((....*(Ln*A)=UU 为上三角矩阵。 之所以这样化简是因为上三角矩阵便于求解方程组（当然这不是本文的主题，只是为了说明 LU 分解的来历）。将所有的消去矩阵按顺序相乘，得到 L ’ *A=U 这样就把 Ax=b 的方程组，变换为求解 Ux=b 三角方程的形式。之所以说 LU 与高斯消去法有紧密的联系，是因为 LU 分解在分解的过程中同样用到了高斯消去法。 （2） LU 分解 LU 分解是将一个 n 阶的非奇异矩阵 A 分解成一个单位下三角矩阵 L 和上三角矩阵 U 的乘积，称作矩阵 A 的 LU 分解。如下图所示： L 和 U 的计算公式如下： 1.2 LU 分解脉动阵列的设计 为了提高 PE 的利用率，让 A 、 L 、 U 三个矩阵的元素中选取 L 矩阵的元素存于 PE 之中，使另外两个矩阵的元素沿不同的方向流动，即 A 、 U 中的元素分别沿水平向右与垂直向下的方向流动，构成如下图所示的三角状阵列。阵列中的正方形 PE 为内积步 PE 。当 A 的主对角线元素到达 PE （ i,j ）（ i 不等于 j ）时，产生 L 矩阵的元素 L ij ，且从此保存与 PE 中。阵列中的圆形 PE 仅起到改变方向的作用。 （1） 数据流向说明 A 的元素沿水平方向流经各正方形 PE 并被处理后，其下三角元素在 PE 中被除以 a jj 后变为 L ij 而存于 PE （ i,j ）中； A 的上三角元素到达圆形 PE 时，变为 u ij 而改为垂直向下流动，参加后面 L 、 U 值得计算。最后 U 的元素由阵列底部输出。 LU 分解脉动设计图 注： 1 上述阵列需要 n （ n+1 ） /2 个 PE ，所需时间为 3N-2 个时钟周期。 1.3 结果验证 （1） matlab 分解结果 矩阵A 矩阵 U 矩阵 L （2）Modelsim 分解结果： 矩阵 U 矩阵 L 大家应该注意到，矩阵L的结果是使用单精度的浮点数表示的，大家可以参照《 矩阵开篇之浮点转定点 》，练习练习，然后跟matlab的结果做对比，看看自己是否掌握了该转换方法。

上一篇说到一个基本的小知识点浮点到定点的转换，这一篇来说说矩阵乘法。矩阵乘法和下一篇要说的矩阵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的结果对比，验证其正确性。 上一篇《 矩阵系列之开篇浮点转定点 》