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

news2024/11/15 17:33:08

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

前言

一、上三角矩阵单位化

二、C对角矩阵求逆

三、A' 矩阵求逆

四、A矩阵求逆

五、计算量分析

六、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}$