numpy实现高等代数矩阵的求解过程

1 矩阵秩

代数中求矩阵的秩

不为0的维数为2，则秩为2

下面我们看下numpy实现

import numpy as np
A = np.mat([
    [3,2,1,1],
    [1,2,-3,2],
    [4,4,-2,3]],int)

rank = np.linalg.matrix_rank(A)
print(rank)

同样结果是rank=2
利用A.shape[0]>rank 来判断方程组有多少解

求解线性方程组

用矩阵表示方程组

上面一个是系数矩阵，一个是曾广矩阵

答案很明显了，解方程也就是对矩阵的变换
这类方程组分明显是非齐次线性方程组

非齐次线性方程组唯一解

非齐次线性方程组是常数项不全为零的线性方程组
非齐次线性方程组有解的条件是秩相同，也就是rankA=n（阶数）
numpy实现如下

import numpy as np

A = np.mat([
    [1, -2, 1],
    [0, 2, -8],
    [-4, 5, 9]], int)
B = np.array([[[0], [8], [-9]]], int)

rank = np.linalg.matrix_rank(A)
print(rank)  
print(A.shape[0])
print(np.linalg.solve(A,B))

结果输出如下

其中A.shape[0]就是n

无限多解的非齐次线性方程组

有结果我们可以看到rank=n-1，那么矩阵对应的方程组有无穷多解,如果在实战中出现这种情况就说明其中一个因子是无效的，与结果不相关的
如下线性方程组

这个矩阵同样是3阶4列，其中x3可以任意值，就造成了多解的情况。

矩阵分解如下



numpy实现如下

import numpy as np

A = np.mat([
    [1, -1, 1, -1],
    [0, 0, 2, -4]], int)
B = np.array([[[0], [1]]])

rank = np.linalg.matrix_rank(A)
print(rank)
print(A.shape[0])
print(np.linalg.pinv(A).dot(B))

最终输出如下

输出为方程组的最小范数解为 K* [-0.13636364 0.13636364 0.04545455 -0.22727273] ,其中K为任意常数
注意为了满足矩阵满秩的要求，我们只用到了[ [1, -1, 1, -1],[0, 0, 2, -4]] 和[[0], [1]]，其他的舍弃

无解非齐次线性方程组的线性拟合

方程组AX=b无解，说明由矩阵A的列向量进行线性组合，无法得到b，有时候我们退而求其次，欲把A的列向量线性组合到最接近b的值，记作b^{。这时候的解称为变形解X}，AX^=b，b^ <> b但接近b。

此方程组矩阵的秩R(A)=2,R(A,b)=3，因此无解。但可以输入计算机求变形解X^:

numpy实现如下

import numpy as np

A = np.mat([
    [4, 2, -1],
    [3, -1, 2],
    [11, 3, 0]], int)
B = np.array([[2], [10], [8]])

rank = np.linalg.matrix_rank(A)
print(rank)
print(A.shape[0])
print(np.linalg.pinv(A).dot(B))

显示[b-b^]的绝对值总和为4，具体到3个式子，一个式子的误差在1.3左右。

n元齐次线性方程组求解

齐次线性方程组是常数项全为零的线性方程组


求解如下4元齐次线性方程组

numpy实现如下

import numpy as np

A= np.mat([[1,-1,-1,1],
           [1,-1,1,-3],
           [1,-1,-2,3]],int)
print(np.linalg.svd(A))

输出如下

其中第2个矩阵只列出A的奇异值，本例中有3个数，但最后的数很小e-16)应视为0。这样可以确定R(A)=2, (n-r}=2，我们从第3个矩阵截取下面的2个行向量即可。

我们还可以直接求解Ax=0

import numpy as np

A = np.mat([[1, -1, -1, 1],
            [1, -1, 1, -3],
            [1, -1, -2, 3]], int)


def derx(A, eps=1e-15):
    u, s, vt = np.linalg.svd(A)
    n_space = np.compress(s <= eps, vt, axis=0)
    return n_space.T


print(derx(A).T)

输出如下