Python解非满秩的线性方程组

之前也是需要用的时候发现Numpy只能解满秩的，解不了非满秩有无穷解的，也就是特解+核空间线性组合的方法什么的。后面百度到了结果，下面给出知乎大佬的代码，这个函数复制下就能直接用啦。

Python 如何计算线性方程组 Ax=b 的通解？ - 用我双手乾坤转过的回答 - 知乎 https://www.zhihu.com/question/36865787/answer/2669793744

这个也就是线性代数中增广矩阵高斯消元解线性方程组的规范方法。

import numpy as np
from numpy.linalg import matrix_rank, inv, multi_dot

def mylinalgsolve(A,b):

    m, n = A.shape
    rank_a = matrix_rank(A)
    rank_a_b = matrix_rank(np.concatenate((A, b), axis=1))
    print("-----------START-----------")
    print("The coefficient matrix A is %s, rank = %d." % (A.shape, rank_a))
    if rank_a == m:
        print("The column space of A encompasses R^%d." % m)
    else:
        print("The column space of A does not encompass R^%d." % m)
    if rank_a < rank_a_b:
        print("rank(A) = %d < %d = rank(A|b), the equation has no solution." % (rank_a, rank_a_b))
    else:
        if rank_a == n:
            print("rank(A) = rank(A|b) = %d, the equation has only one solution." % n)
            Ar = A
            br = b
            if n < m:
                Ar = Ar[:n]
                br = br[:n]
            Ar_inv = inv(Ar)
            if n < m:
                print("One of the largest sub square of A is T =\n%s" % Ar)
                print("T^-1 =\n%s" % Ar_inv)
                print("x = T^-1 * b = ", end="")
            else:
                print("A^-1 =\n%s" % Ar_inv)
                print("x = A^-1 * b = ", end="")
            x = np.round(multi_dot((Ar_inv, br)), 4).flatten()
            print("%s^T" % x)
        else:
            print("rank(A) = rank(A|b) = %d < %d, the equation has infinitely many solutions." % (rank_a, n))
            # find the solution of T*y = b, where T is compounded of linearly independant column vectors of A
            Ar = A[:, :rank_a]
            br = b
            if rank_a < m:
                Ar = Ar[:rank_a]
                br = br[:rank_a]
            Ar_inv = inv(Ar)
            # find a specific solution
            y = np.concatenate((multi_dot((inv(Ar), br)), np.zeros((n - rank_a, 1))), axis=0)
            # find every linear combination of T representating redundant column vectors of A
            co_independant = np.zeros((n - rank_a, rank_a))
            co_redundant = -np.eye(n - rank_a)
            for i in range(rank_a, n):
                c = A[:rank_a, i]
                co_independant[i - rank_a] = multi_dot((Ar_inv, c)).T
            # find the solutions of homogeneous equation Ax=0
            z = np.concatenate((co_independant, co_redundant), axis=1).T
            y = y
            z = z
            result = "%s^T" % y.flatten()
            for i in range(n - rank_a):
                result += " + k_%d * %s^T" % (i + 1, z[:, i])
            print("x = %s (k_n ∈ R)" % result)
    print("-----------END-----------")
    return result

下面给出一个算例解AX=b，A不满秩的情况。注意下numpy的结构和转换就好。

x=[2,-3,-4,-6,5,3,-2,4]
y=[3,4,-5,2,3,2,-6,8]
x=np.array(x)
y=np.array(y)

b=np.zeros_like(y)
A=np.array([x**3,x**2*y,x*y**2,y**3,x**2,y**2,x*y,x,y,np.ones_like(y)]).T
mylinalgsolve(A,b.reshape(8,1))


'''
-----------START-----------
The coefficient matrix A is (8, 10), rank = 8.
The column space of A encompasses R^8.
rank(A) = rank(A|b) = 8 < 10, the equation has infinitely many solutions.
x = [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]^T + k_1 * [-0.11287802  0.46301482 -0.3657235   0.08385066 -1.08818233 -0.46186757
  1.56357684  0.89698889 -1.         -0.        ]^T + k_2 * [-0.0198804   0.05238849 -0.05231127  0.00908104 -0.10905413 -0.01730029
  0.17617317  0.38083805 -0.         -1.        ]^T (k_n ∈ R)
-----------END-----------
'''