一种基于约化因子上三角矩阵求逆方法与MATLAB仿真

2024-03-23 12:36:04
开发
18

一种基于约化因子上三角矩阵求逆的方法与MATLAB仿真

前言

一、上三角矩阵单位化

前言

矩阵运算广泛应用于实时性要求的各类电路中，其中矩阵求逆运算最难以实现。本文是在阅读文献后，仿真文中采用的一种约化因子求逆的优化算法，将任意一个 n×n 阶上三角矩阵转换成对角线为 1 的上三角矩阵，使得除法运算与乘加运算分离开来，大大简化矩阵求逆运算过程。文献中有些地方表述有误，在撰写本文时已经改正。

提示：以下是本篇文章正文内容，希望能帮助到各位，转载请附上链接。

一、上三角矩阵单位化

假设A为N×N阶上三角矩阵，其对角线元素非0，可以通过变换将其用一个对角矩阵C和对角线为1的上三角矩阵A'的乘积表示。

$\textbf{CA}'=\textbf{A}$

其中

$\textbf{A}_{N\times N}=\begin{pmatrix}c_{11}&c_{12}&c_{13}&\cdots&c_{1n}\\0&c_{22}&c_{23}&\cdots&c_{2n}\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&c_{n-1n-1}&c_{n-1n}\\0&0&0&\cdots&c_{n\times n}\end{pmatrix}$

$\textbf{C}_{N\times N}=\begin{pmatrix}c_{11}&0&0&\cdots&0\\0&c_{22}&0&\cdots&0\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&c_{33}&0\\0&0&0&\cdots&c_{44}\end{pmatrix}$

$\textbf{A}_{N\times N}'=\begin{pmatrix}1&a_{12}&a_{13}&\cdots&a_{1n}\\0&1&a_{23}&\cdots&a_{2n}\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&1&a_{n-1n}\\0&0&0&\cdots&1\end{pmatrix}$

通式为

$a_{mn}=\frac{C_{mn}}{C_{mm}},1\leqslant m\leqslant N ,2\leqslant n\leqslant N$

二、C对角矩阵求逆

对于C这个对角矩阵求逆很简单，只需要对其对角元取倒数即可。

$\textbf{C}_{_{N\times N}}^{-1}=\begin{pmatrix}\frac{1}{c_{11}}&0&0&\cdots&0\\0&\frac{1}{c_{22}}&0&\cdots&0\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&\frac{1}{c_{33}}&0\\0&0&0&\cdots&\frac{1}{c_{44}}\end{pmatrix}$

三、A' 矩阵求逆

对于上三角矩阵A'，每个 2 行 2 列交集所构成的二阶子矩阵，称为操作矩阵，如

$\begin{pmatrix}a_{12}&a_{13}\\1&a_{23}\end{pmatrix}$

构成一个操作矩阵，每次进行行变换的效果体现在操作矩阵的右上角元素上，即令 $b_{13}=a_{13}-a_{12}a_{23}$ ，则 $b_{13}$ 称为约化因子，其通式为

$b_{mn}=a_{mn}-\sum_{k=m+1}^{n-1}(a_{mk}b_{kn}),n\geqslant3$

则可求出A' 的逆矩阵，即为

$\mathbf{A}_{N\times N}^{'-1}=\begin{pmatrix}1&-a_{12}&-b_{13}&\cdots&-b_{1n}\\0&1&-a_{23}&\cdots&-b_{2n}\\\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&\cdots&1&-a_{n-1n}\\0&0&0&\cdots&1\end{pmatrix}$

需要计算(N-1)(N-2)/2个元素。

根据公式，计算 $b_{mn}$ 时，会发现需要先求 $b_{kn}$ ，其他都已知，会发现他俩的下标n一样，区别在于表示行的下标，而且k>m，所以计算时每一列从下至上依次计算，列与列之间互不影响。

谈到这里，是不是会发现如果在FPGA上面用此方法计算上三角矩阵的逆矩阵会大大提高速度，因为关键的计算就在于A' 的逆矩阵的计算，而计算其逆矩阵时列与列又互不影响，直接可以并行计算。

四、A矩阵求逆

由于 $\textbf{CA}'=\textbf{A}$ ，那么

$\textbf{A}^{-1}=(\textbf{CA}')^{-1}=\textbf{A}'^{-1}\textbf{C}^{-1}$

由于 $\textbf{C}^{-1}$ 是个对角阵，所以也可以表示为 $\textbf{C}^{-1}\textbf{A}'^{-1}$ 。

五、计算量分析

限制 $N\geqslant 3$ ，否则根据前面的分析根本不需要求 $b_{mn}$ 。

对N×N阶对角矩阵C求逆需要N次除法。

将一个主对角元素非0且非1的N×N阶矩阵A变为一个主对角元素全1的矩阵A'和对角矩阵C需要运算的乘法次数（将除法看成去乘以C逆矩阵的对角元）为：

$\frac{N(N-1)}{2}$

求 $\mathbf{A}_{N\times N}^{'-1}$ 的第n（n≥3）列需要的乘法次数与加法（在FPGA中实现时，减法可以看成加法）次数相等，均为

$1+2+\cdots +n-2=\frac{(n-1)(n-2)}{2}$

那么计算所有的列需要的乘法与加法次数均为：

$\sum_{n=3}^{N} \frac{(n-1)(n-2)}{2} \\= \frac{1}{2} \sum_{n=1}^{N} (n^2 - 3n + 2)\\=\frac{1}{2} \left( \frac{N(N+1)(2N+1)}{6} - 3 \cdot \frac{N(N+1)}{2} + 2N \right)\\=\frac{N^{3}-3N^{2}+2N}{6}$

其中用到公式

$\sum_{n=1}^{N} n^2=\frac{N(N+1)(2N+1)}{6}$

最后计算 $\textbf{A}'^{-1}\textbf{C}^{-1}$ 需要的乘法次数为

$\frac{N(N-1)}{2}$

因为 $\textbf{A}'^{-1}$ 是对角元为1的上三角矩阵， $\textbf{C}^{-1}$ 是对角矩阵，所以只有上述表达式的乘法次数。

那么总的运算量为：

乘法： $\frac{N(N-1)}{2}+\frac{N^{3}-3N^{2}+2N}{6}+\frac{N(N-1)}{2}=\frac{N^{3}+3N^{2}-4N}{6}$

加法： $\frac{N^{3}-3N^{2}+2N}{6}$

除法： $N$

特殊的，如果待求矩阵A本来就是一个主对角线全为1的上三角矩阵，那么就可以省略上三角矩阵单位化，那么阵个求逆只会用到乘法和加法各 $\frac{N^{3}-3N^{2}+2N}{6}$ 次。

六、MATLAB仿真

以MATLAB自带求逆函数inv为对比，仿真得出以下结果：

误差定义及源代码见参考资料。

七、参考资料

https://download.csdn.net/download/m0_66360845/89011663

总结

以上介绍了一种基于约化因子上三角矩阵求逆的方法与MATLAB仿真，小伙伴们认真看完必定有收获。

原文地址:https://blog.csdn.net/m0_66360845/article/details/136931187 本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：https://www.suanlizi.com/kf/1771395461043851264.html 如若内容造成侵权/违法违规/事实不符，请联系《酸梨子》网邮箱：1419361763@qq.com进行投诉反馈，一经查实，立即删除！

阅读全部