矩阵求逆运算 Python实现

原理：应用列主元消去法运算矩阵A的逆矩阵，利用初等矩阵行变换A转化单位矩阵时，同样的行变化可将单位矩阵转化为A的逆矩阵。

步骤：

1. 编制下三角部分消元和上三角部分消元的代码
   a. 从对角线元素往下比较取得这一列的最大值所在的行，与对角线元素所在行进行交换。
   b. 从对角线所在行往下，利用矩阵的行变换将这一列下所有元素消为0，存储消元过程中的系数（三角矩阵）
   c. 类似b步骤从最后一列往前消元消成对角矩阵，同样储存消元过程中的系数
   d. 从第一列开始依次进行a，b，c步骤

2. 对角单位化（用一元数组存储）
   矩阵每一行除对角线元素，将矩阵化为单位矩阵。

3. 生成一个与矩阵A同阶的单位矩阵，对单位矩阵进行上述过程相同的行变化（取上述变换生成的主元交换和消元参数矩阵），即可以得到矩阵A的逆

# -*- coding: utf-8 -*-
"""
@author: machineRandy
"""
import warnings
warnings.filterwarnings("ignore")

matrix1 = [
    [1,2,0,0],
    [3,4,0,0],
    [0,0,4,1],
    [0,0,3,2],
]

matrix2 = [
    [1,0,-1,2,1],
    [3,2,-3,5,-3],
    [2,2,1,4,-2],
    [0,4,3,3,1],
    [1,0,8,-11,4],
]

matrix3 = [
    [1,0,-1,2,1,0,2],
    [1,2,-1,3,1,-1,4],
    [2,2,1,6,2,1,6],
    [-1,4,1,4,0,0,0],
    [4,0,-1,21,9,9,9],
    [2,4,4,12,5,6,11],
    [7,-1,-4,22,7,8,18],
]


def step0(m):
    n = len(m)
    l = []
    for i in range(0,n):
        l.append([])
        for j in range(0,n):
            if i == j:
                l[i].append(1)
            else:
                l[i].append(0)
    return l

def step1(m):
    n = len(m)
    """交换操作记录数组 swap"""
    swap = []
    l = []
    for i in range(0,n):
        swap.append(i)
        l.append([])
        for j in range(0,n):
            l[i].append(0)

    """对每一列进行操作"""
    for i in range(0,n):
        max_row = m[i][i]
        row = i
        for j in range(i,n):
            if m[j][i] >= max_row:
                max_row = m[j][i]
                global row
                row = j
        swap[i] = row

        """交换"""
        if row != i:
            for j in range(0,n):
                m[i][j],m[row][j] = m[row][j],m[i][j]

        """消元"""
        for j in range(i+1,n):
            if m[j][i] != 0:
                l[j][i] = m[j][i] / m[i][i]
                for k in range(0,n):
                    m[j][k] = m[j][k] - (l[j][i] * m[i][k])

    return (swap,m,l)

def step2(m):
    n = len(m)
    long = len(m)-1
    l = []
    for i in range(0,n):
        l.append([])
        for j in range(0,n):
            l[i].append(0)

    for i in range(0,n-1):
        for j in range(0,long-i):
            if m[long-i-j-1][long-i] != 0 and m[long-i][long-i] != 0:
                l[long-i-j-1][long-i] = m[long-i-j-1][long-i] / m[long-i][long-i]
                for k in range(0,n):
                    m[long-i-j-1][k] = m[long-i-j-1][k] - l[long-i-j-1][long-i] * m[long-i][k]


    return (m,l)

def step3(m):
    n = len(m)
    l = []
    for i in range(0,n):
        l.append(m[i][i])
    return l



def gauss(matrix):
    n = len(matrix)
    new = step0(matrix)
    (swap,matrix1,l1) = step1(matrix)
    (matrix2,l2) = step2(matrix1)
    l3 = step3(matrix2)
    for i in range(0,n):
        if swap[i] != i:
            new[i],new[swap[i]] = new[swap[i]],new[i]
        for j in range(i+1,n):
            for k in range(0,n):
                if l1[j][i] != 0:
                    new[j][k] = new[j][k] - l1[j][i] * new[i][k]   
    for i in range(0,n-1):
        for j in range(0,n-i-1):
            if l2[n-1-i-j-1][n-1-i] != 0:
                for k in range(0,n):
                    new[n-1-i-j-1][k] = new[n-1-i-j-1][k] - l2[n-1-i-j-1][n-i-1] * new[n-1-i][k]
    for i in range(0,n):
        for j in range(0,n):
            new[i][j] = new[i][j] / l3[i]
    return new

x1 = gauss(matrix1)
x2 = gauss(matrix2)
x3 = gauss(matrix3)