一般实矩阵求逆

文章目录

一般实矩阵求逆
- 实矩阵求逆

实矩阵求逆

我们所学习过的矩阵求逆的方法有许多，如高斯-约当消元法、伴随矩阵法、LU分解法等。在此我会一一介绍。

我们先来复习一下逆矩阵：

逆矩阵的定义：

假设 $A$ 是一个方阵，如果存在一个矩阵 $A^{-1}$ ，使得 $A^{-1}A=I$ 并且 $AA^{-1}=I$ ，那么称矩阵 $A$ 是可逆的， $A^{-1}$ 则称为矩阵 $A$ 的逆矩阵。

首先就是逆矩阵的存在性，若逆矩阵都不存在，我们自然就不需要去求这个逆矩阵了。

若矩阵 $A$ 可逆，则有 $\neq 0$

若 $\neq 0$ ，则矩阵 $A$ 可逆，且 $A^{-1} = \frac{1}{|A|} \dot{} A^*$ ，其中 $A^*$ 称为矩阵 $A$ 的伴随矩阵

高斯-约当消元法

这方法听着很高大上，但却是大部分人经常用的消去法，若有一个矩阵 $A_{n\times n}$ ，满足 $AA^{-1}=I$ 。

基于这个方法我们构造一个增广矩阵 $W = [A ∣ I]$ ，则 $W_{n\times 2n}$ 是一个 $n\times 2n$ 阶矩阵，我们对其进行变化，使其变为 $[I ∣ B]$ 类型的矩阵，若变换成功，则称矩阵 $B$ 为矩阵 $A$ 的逆矩阵。

具体操作方法：

（1）将矩阵进行扩展为 $W = [A ∣ I]$

（2）从上到下做行变换，将 $W$ 的左半部分转换为上三角矩阵

（3）将 $W$ 的左半部分转化为对角矩阵

（4）将 $W$ 的左半部分乘以系数使之变为单位矩阵

（5）提取出 $W$ 的右半部分，即为矩阵 $A$ 的逆矩阵

在转化过程中若 $W$ 矩阵前半部分不能单位化，则判定 $A$ 矩阵没有逆矩阵

我们将接口设为：int Guassian(double* a, double* b, const int n)

接口与参数说明：

形参与函数类型	参数意义
`double* a`	存放矩阵 $A$ 的元素
`double* b`	存放 $A$ 的逆矩阵的元素
`const int n`	存放矩阵的阶数
`int Guassian()`	若有逆矩阵，则返回1，若无逆矩阵，则返回0

现在就可以将其代码写出来了：

void Guassian(double* a, double* b, const int n) {
    int i, j, k;
    double temp1, temp2, temp3;
    double w[n][2*n];
    double result[n][n];

    // 对矩阵a进行扩展
    for(i = 0; i != n; ++i) {
        for(j = 0; j != 2*n; ++j) {
            if(j < n){
                w[i][j] = a[n*i+j];
            } else {
                w[i][j] = j-n == i ? 1 : 0;
            }
        }
    }

    // 化为对角线元素为1的上三角矩阵
    for(i = 0; i != n; ++i) {
        // 首元素等于0 则将后面随便一行同一列元素不为0的值加到该行
        if((int)w[i][i] == 0) {
            // 寻找同列首元素不为0的行数
            for(j = i+1; j != n; ++j) {
                if((int)w[j][i] != 0)
                    break;
            }
            if(j == n) {
                std::cout << "无法求逆" << std::endl;
                break;
            }
            // 该行的元素加上找到行的元素
            for(k = 0; k < 2*n; ++k) {
                w[i][k] += w[j][k];
            }
        }

        // 将该行首元素置为1
        temp1 = w[i][i];
        for(j = 0; j != 2*n; ++j) {
            w[i][j] = w[i][j] / temp1;
        }

        // 将后面元素的首元素置为0
        for(j = i+1; j != n; ++j) {
            temp2 = w[j][i];
            for(k = i; k != 2*n;++k) {
                w[j][k] = w[j][k] - temp2 * w[i][k];
            }
        }
    }

    // 将前半部分标准化
    for(i = n-1; i >= 0; --i) {
        for(j = i-1; j >= 0; --j) {
            temp3 = w[j][i];
            for(k = i; k < 2*n; ++k) {
                w[j][k] = w[j][k] - temp3 * w[i][k];
            }
        }
    }

    // 得出逆矩阵
    for(i = 0; i != n; ++i) {
        for(j = n; j != 2*n; ++j) {
            b[i*n+(j-n)] = w[i][j];
        }
    }
}

在这里求一个难一点的矩阵的逆：
$\begin {pmatrix} 0.2368 & 0.2471 & 0.2568 & 1.2671\\ 1.1161 & 0.1254 & 0.1397 & 0.1490\\ 0.1582 & 1.1675 & 0.1768 & 0.1871\\ 0.1968 & 0.2071 & 1.2168 & 0.2271 \end {pmatrix}$
运行程序如下：

#include "Guassian.hpp"

int main() {
    double a[4][4] = {{0.2368, 0.2471, 0.2568, 1.2671},
                      {1.1161, 0.1254, 0.1397, 0.1490},
                      {0.1582, 1.1675, 0.1768, 0.1871},
                      {0.1968, 0.2071, 1.2168, 0.2271}};
    double b[4][4];
    Guassian((double*)a, (double*)b, 4);

    for(int i = 0; i != 4; ++i) {
        for(int j = 0; j != 4; ++j) {
            std::cout << b[i][j] << " ";
        }
        std::cout << std::endl;
    }
    return 0;
}

运行结果如下：

-0.0859208 0.937944 -0.0684372 -0.0796077
-0.10559 -0.0885243 0.905983 -0.0991908
-0.127073 -0.111351 -0.116967 0.878425
0.851606 -0.135456 -0.140183 -0.143807

在这里由于C++的二维数组传参极其严格，因此可以运用模板降低标准，使其方便调用，详见上一篇博客，在此不再赘述。

全选主元高斯-约当消去

我们换一个角度来看高斯-约当消元法：

对于 $k = 1, 2, . . ., n$ 做如下计算：

（1）归一化计算：
$\frac{1}{a_{kk}} \Rightarrow a_{kk}\\ a_{kj}a_{kk} \Rightarrow a_{kj}, j = 1, 2, ..., n;j \neq k$
（2）消元计算：
$KaTeX parse error: Expected 'EOF', got '&' at position 41: …htarrow a_{ij} &̲i=1,2,...,n;i \…$

普通的高斯-约当消去已经可以了，但算法不稳定，容易出现误差与频繁叠加行元素等。所以为了数值计算的稳定性，这整个过程需要全选主元，全选主元的意思就是将绝对值最大的元素交换到主对角线上，在这里即交换到左半部分的主对角线上。全选主元的过程如下所示：

对于矩阵求逆过程中的第 $k$ 步，首先，在 $a_{kk}$ 的右下方（包括 $a_{kk}$ 自己）的 $n - k + 1$ 阶子矩阵中选取绝对值最大的元素，并将该元素所在的行号记录在 $I S (k)$ 中，而列号记录在 $J S (k)$ 中，然后通过行交换与列交换将该绝对值最大的元素交换到主对角线 $a_{kk}$ 的位置上，即做以下两步：
$\Leftrightarrow A(IS(k),l),l=1,2,...,n\\ A(l,k) \Leftrightarrow A(l,JS(k)),l=1,2,...,n\\$
经过全选主元过后，矩阵求逆的过程是稳定的。但最后需要对结果进行恢复。恢复过程如下：
$\Leftrightarrow A(JS(k),l),l=1,2,...,n\\ A(l,k) \Leftrightarrow A(l,IS(k)),l=1,2,...,n\\$

我们将接口设为：int GuassianRinv(double* a, double* b, const int n)

接口与参数说明：

形参与函数类型	参数意义
`double* a`	存放矩阵 $A$ 的元素，返回时存放 $A$ 的逆矩阵
`const int n`	存放矩阵的阶数
`int GuassianRinv()`	若有逆矩阵，则返回1，若无逆矩阵，则返回0

现在就可以将其代码写出来了：

template <int N>
int GuassianRinv(double a[N][N]) {
    int is[N], js[N];
    int i, j, k;
    double d, p;

    for(k = 0; k != N; ++k) {
        d = 0.0;
        // 找出最大元素并存储
        for(i = k; i != N; ++i) {
            for(j = k; j != N; ++j) {
                p = fabs(a[i][j]);
                if(p > d) {
                    d = p;
                    is[k] = i;
                    js[k] = j;
                }
            }
        }

        // 没有逆矩阵
        if(d + 1.0 == 1.0) {
            return 0;
        }

        // 使最大元素在对角线上，行交换
        if(is[k] != k) {
            for(j = 0; j != N; ++j) {
                p = a[k][j];
                a[k][j] = a[is[k]][j];
                a[is[k]][j] = p;
            }
        }

        // 使最大元素在对角线上，列交换
        if(js[k] != k) {
            for(i = 0; i != N; ++i) {
                p = a[k][i];
                a[k][i] = a[js[k]][i];
                a[js[k]][i] = p;
            }
        }

        // 归一计算
        a[k][k] = 1.0/a[k][k];
        for(j = 0; j != N; ++j) {
            if(j != k) {
                a[k][j] *= a[k][k];
            }
        }

        // 消元计算
        for(i = 0; i != N; ++i) {
            if(i != k) {
                for(j = 0; j != N; ++j) {
                    if(j != k) {
                        a[i][j] -= a[i][k] * a[k][j];
                    }
                }
            }
        }

        // 消元计算
        for(i = 0; i != N; ++i) {
            if(i != k) {
                a[i][k] = -a[i][k] * a[k][k];
            }
        }
    }

    // 将结果恢复
    for(k = N-1; k >= 0; --k) {
        if(js[k] != k) {
            for(j = 0; j != N; ++j) {
                p = a[k][j];
                a[k][j] = a[js[k]][j];
                a[js[k]][j] = p;
            }
        }

        if(is[k] != k) {
            for(i = 0; i != N; ++i) {
                p = a[k][i];
                a[k][i] = a[js[k]][i];
                a[js[k]][i] = p;
            }
        }
    }

    return 1;
}

在这里求与上面一样的矩阵的逆：
$\begin {pmatrix} 0.2368 & 0.2471 & 0.2568 & 1.2671\\ 1.1161 & 0.1254 & 0.1397 & 0.1490\\ 0.1582 & 1.1675 & 0.1768 & 0.1871\\ 0.1968 & 0.2071 & 1.2168 & 0.2271 \end {pmatrix}$
运行程序如下：

#include "GuassianRinv.hpp"

int main() {
    double a[4][4] = {{0.2368, 0.2471, 0.2568, 1.2671},
                      {1.1161, 0.1254, 0.1397, 0.1490},
                      {0.1582, 1.1675, 0.1768, 0.1871},
                      {0.1968, 0.2071, 1.2168, 0.2271}};
    GuassianRinv<4>(a);

    for(int i = 0; i != 4; ++i) {
        for(int j = 0; j != 4; ++j) {
            std::cout << a[i][j] << " ";
        }
        std::cout << std::endl;
    }

    return 0;
}

运行结果如下：

-0.0859208 0.937944 -0.0684372 -0.0796077
-0.10559 -0.0885243 0.905983 -0.0991908
-0.127073 -0.111351 -0.116967 0.878425
0.851606 -0.135456 -0.140183 -0.143807

伴随矩阵法

一个方阵 $A$ 如果满足 $\neq 0$ ，则可以认为矩阵 $A$ 可逆，它的逆矩阵可以表示为:
$A^{-1} = \frac{1}{|A|}A^*$
其中求出伴随矩阵是使用该方法求逆的最关键的一步，如何使用矩阵 $A$ 求解它的伴随矩阵 $A *$ 呢？方法如下：

（1）求出该矩阵的余子式 $M_{ij}$ ， $M_{ij}$ 为矩阵 $A$ 去掉第 $i$ 行，第 $j$ 列之后的行列式的值

（2） $A_{ij}$ 称为该矩阵的代数余子式，其中 $A_{ij}=(-1)^{i+j} \dot{} M_{ij}$

（3）运用以下公式求出该矩阵的伴随矩阵：
$A^* = \begin{bmatrix} A_{11} & A_{21} & \dots & A_{n1}\\ A_{12} & A_{22} & \dots & A_{n2}\\ \vdots & \vdots & \ddots & \vdots\\ A_{1n} & A_{2n} & \dots & A_{nn}\\ \end{bmatrix}$
我们对其进行接口设计：int adjoint(double* a, double* b, int n)

接口与参数说明：

形参与函数类型	参数意义
`double* a`	存放矩阵 $A$ 的元素
`double* b`	存放 $A$ 的逆矩阵的元素
`int n`	存放矩阵的阶数
`int adjoint()`	若有逆矩阵，则返回1，若无逆矩阵，则返回0

现在就可以将其代码写出来了：

#include "det.hpp" // 上一篇博客中的代码

template <int N>
int adjoint(double a[N][N], double b[N][N]) {
    int i, j, k, m, n, d;
    double w[N][N], tem[N-1][N-1], result[N][N];
    double detValue, temp;

    // 以下做行列式计算会破坏原行列式，先进行备份
    for(i = 0; i != N; ++i) {
        for(j = 0; j != N; ++j) {
            w[i][j] = a[i][j];
        }
    }

    // 计算矩阵A的行列式|A|
    detValue = det<N>(w);
    if(detValue + 1.0 == 1.0) {
        return 0;
    }

    // 求矩阵的余子式矩阵
    for(i = 0; i != N; ++i) {
        for(j = 0; j != N; ++j) {
            // 计算每个元素
            k = 0;
            for(m = 0; m != N-1; ++m) {
                if(k == i) {
                    ++k;
                }

                d = 0;
                for(n = 0; n != N-1; ++n) {
                    if(d == j) {
                        ++d;
                    }

                    tem[m][n] = a[k][d];
                    ++d;
                }
                ++k;
            }

            // 计算余子式矩阵
            temp = det<N-1>(tem);

            // 进行转置
            result[j][i] = pow(-1, i+j) * temp;
        }
    }

    for(i = 0; i != N; ++i) {
        for(j = 0; j != N; ++j) {
            b[i][j] = result[i][j] / detValue;
        }
    }
    return 1;
}

我们运用同样的数据进行检验：

#include "adjoint.hpp"

int main() {
    double a[4][4] = {{0.2368, 0.2471, 0.2568, 1.2671},
                      {1.1161, 0.1254, 0.1397, 0.1490},
                      {0.1582, 1.1675, 0.1768, 0.1871},
                      {0.1968, 0.2071, 1.2168, 0.2271}};
    double b[4][4];
    adjoint<4>(a, b);

    for(int i = 0; i != 4; ++i) {
        for(int j = 0; j != 4; ++j) {
            std::cout << b[i][j] << " ";
        }
        std::cout << std::endl;
    }

    return 0;
}

输出结果如下：

-0.0859208 0.937944 -0.0684372 -0.0796077
-0.10559 -0.0885243 0.905983 -0.0991908
-0.127073 -0.111351 -0.116967 0.878425
0.851606 -0.135456 -0.140183 -0.143807

LU分解法

相信大名鼎鼎的LU分解法大部分人都是听过的，先来看一下LU分解法求逆的可行性：
$KaTeX parse error: Expected 'EOF', got '&' at position 2: &̲A=(L*U)\\ \Righ…$
众所周知， $L$ 矩阵为下三角矩阵， $U$ 矩阵为上三角矩阵，因此不论是列顺序消去或是行顺序消去都会比普通矩阵简单快捷。具体步骤如下：

1、把矩阵 $A$ 分解为 $L$ 下三角矩阵和 $U$ 上三角矩阵

以 $3\times3$ 的矩阵为例，其LU分解可以写成以下形式：
$\begin{bmatrix} a_{11} & a_{12} &a_{13}\\ a_{21} & a_{22} &a_{23}\\ a_{31} & a_{32} &a_{33}\\ \end{bmatrix}= \begin{bmatrix} l_{11} & 0 & 0\\ l_{21} & l_{22} & 0\\ l_{31} & l_{32} & l_{33}\\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0 & 0 & u_{33}\\ \end{bmatrix}$
LU分解在本质上就是高斯消元法的一种表达形式而已，就是将矩阵 $A$ 通过行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这就是所谓的杜尔里特算法，从上至下地对矩阵 $A$ 做初等行变换，将对角线左下方元素变为0，然后证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积大的逆就是 $L$ 矩，它也是一个单位下三角矩阵。

因此更准确的表达以上LU分解：
$\begin{aligned} A&= \begin{bmatrix} a_{11} & a_{12} &a_{13}\\ a_{21} & a_{22} &a_{23}\\ a_{31} & a_{32} &a_{33}\\ \end{bmatrix}= \begin{bmatrix} 1 & 0 & 0\\ l_{21} & 1 & 0\\ l_{31} & l_{32} & 1\\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0 & 0 & u_{33}\\ \end{bmatrix}\\&= \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ l_{21}u_{11} & l_{21}u_{12}+u_{22} & l_{21}u_{13}+u_{23}\\ l_{31}u_{11} & l_{31}u_{12}+l_{32}u_{22} & l_{31}u_{13}+l_{32}u_{23}+u_{33}\\ \end{bmatrix} \end{aligned}$
按理来说已经可以按照这个公式来求出 $L$ 矩阵和 $U$ 矩阵了：
$\begin{aligned} U_{0, j}&=A_{0, j}&(j = 0, 1, \dots, m-1)\\\\ L_{i, 0}&=\frac{A_{i, 0}}{U_{00}}&(i = 1, 2, \dots, m-1)\\ \end{aligned}$
还可以看出 $U$ 时按行迭代运算， $L$ 是按列迭代运算， $U L$ 交错运算，且 $U$ 先 $L$ 一步，因此有：
$\begin{aligned} & U_{k, j}=\frac{A_{k, j}-\sum_{t=0}^{k-1}(L_{k,t}U_{t,j})}{L_{k,k}}=A_{k, j}-\sum_{t=0}^{k-1}(L_{k,t}U_{t,j})&(k=1,2,\dots m-1,j=k,\dots m-1)\\ & L_{i,k}=\frac{A_{i, k}-\sum_{t=0}^{k-1}(L_{i,t}U_{t,k})}{U_{k, k}}&(k=1,2, \dots m-1,i=k,\dots m-1) \end{aligned}$
此时便可以分别对矩阵 $L$ 和 $U$ 求逆矩阵：

列顺序行顺序： $\dots,m-1,i=j,\dots,m-1)$ :
$L_{inv(i,j)}= \begin{aligned} \begin{cases} &L^{-1}_{i, j}&i=j\\ &0&i<j\\ &-L_{inv(j,j)}\sum_{k=j}^{i-1}(L_{i,k}L_{inv(k,j)})&i>j \end{cases} \end{aligned}$
列顺序行倒序： $\dots,m-1,i=j,\dots,0)$ :
$U_{inv(i,j)}= \begin{aligned} \begin{cases} &L^{-1}_{i, j}&i=j\\ &0&i<j\\ &-\frac{1}{U_{inv(i,j)}}\sum_{k=i+1}^{j}(U_{i, k}U_{inv(k,j)})&i>j \end{cases} \end{aligned}$

此时便可以通过这两个逆矩阵求出原矩阵的逆： $A^{-1}=U^{-1}*L^{-1}$ :

我们对该函数进行接口设计：void luinv(double* a, double* b, int n) :

接口与参数说明：

形参与函数类型	参数意义
`double* a`	存放矩阵 $A$ 的元素
`double* b`	存放 $A$ 的逆矩阵的元素
`int n`	存放矩阵的阶数

现在就可以将其代码写出来了：

template <int N>
void luinv(double a[N][N], double b[N][N]) {
    double w[N][N], l[N][N], ln[N][N], u[N][N], un[N][N];
    int i, j, k, d;
    double s;

    // 初始化
    for(i = 0; i != N; ++i) {
        for(j = 0; j != N; ++j) {
            w[i][j] = a[i][j];
            b[i][j] = 0;
            l[i][j] = 0;
            ln[i][j] = 0;
            u[i][j] = 0;
            un[i][j] = 0;
        }
    }

    // 将l矩阵的对角线元素置为1
    for(i = 0; i != N; ++i) {
        l[i][i] = 1.0;
    }

    // 根据公式计算出u[0][i]
    for(i = 0; i != N; ++i) {
        u[0][i] = w[0][i];
    }

    // 根据公式计算出u[i][0]
    for(i = 1; i != N; ++i) {
        l[i][0] = w[i][0] / u[0][0];
    }

    for(i = 1; i != N; ++i) {
        // 求出u
        for(j = i; j != N; ++j) {
            s = 0;
            for(k = 0; k != i; ++k) {
                s += l[i][k] * u[k][j];
            }
            u[i][j] = w[i][j] - s;
        }

        // 求出l
        for(d = i; d != N; ++d) {
            s = 0;
            for(k = 0; k != i; ++k) {
                s += l[d][k] * u[k][i];
            }
            l[d][i] = (w[d][i] - s) / u[i][i];
        }
    }

    // 求l的逆
    for(j = 0; j != N; ++j) {
        for(i = j; i != N; ++i) {
            if(i == j) {
                ln[i][j] = 1 / l[i][j];
            } else if (i < j) {
                ln[i][j] = 0;
            } else {
                s = 0;
                for(k = j; k != i; ++k) {
                    s += l[i][k] * ln[k][j];
                }
                ln[i][j] = -ln[j][j] * s;
            }
        }
    }

    // 求u的逆
    for(i = 0; i != N; ++i) {
        for(j = i; j >= 0; --j) {
            if(i == j) {
                un[j][i] = 1 / u[j][i];
            } else if (i < j) {
                un[j][i] = 0;
            } else {
                s = 0;
                for(k = j+1; k <= i; ++k) {
                    s += u[j][k] * un[k][i];
                }
                un[j][i] = (-1.0 / u[j][j]) * s;
            }
        }
    }

    // 将l的逆和u的逆相乘保存到b
    for(i = 0; i != N; ++i) {
        for(j = 0; j != N; ++j) {
            for(k = 0; k != N; ++k) {
                b[i][j] += un[i][k] * ln[k][j];
            }
        }
    }
}

以同样的测试案例进行测试：

#include "luinv.hpp"

int main() {
    double a[4][4] = {{0.2368, 0.2471, 0.2568, 1.2671},
                      {1.1161, 0.1254, 0.1397, 0.1490},
                      {0.1582, 1.1675, 0.1768, 0.1871},
                      {0.1968, 0.2071, 1.2168, 0.2271}};
    double b[4][4];
    luinv<4>(a, b);

    for(int i = 0; i != 4; ++i) {
        for(int j = 0; j != 4; ++j) {
            std::cout << b[i][j] << " ";
        }
        std::cout << std::endl;
    }

    return 0;
}

运行结果如下：

-0.0859208 0.937944 -0.0684372 -0.0796077
-0.10559 -0.0885243 0.905983 -0.0991908
-0.127073 -0.111351 -0.116967 0.878425
0.851606 -0.135456 -0.140183 -0.143807

原文链接：https://blog.csdn.net/qq_50952332/article/details/125033082