Python内置函数进行矩阵化简

Question

Python有没有内置的功能可以把一个矩阵转换成行阶梯形（也叫上三角形）？

Answer 1

可以查看这个链接：http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038705.html

简单来说：别这么做。

在电脑上实现的rref算法会产生很多不准确的结果。所以你要么换个方法来解决问题，要么像@aix建议的那样使用符号计算。

Answer 2

我同意@Mile对@WinstonEwert回答的评论，计算机是可以以给定的精度来执行RREF的。

实现RREF其实并不复杂，matlab已经有这个功能了，所以numpy也应该有。

我做了一个非常简单直接的实现，虽然效率不高，但对于简单的矩阵来说效果还不错：

更新：

看起来@JustMe发现了这个rref实现中的一个bug，这个问题出现在输入矩阵的第一列全是零的情况下。所以这里有一个更新版本。

from numpy import *

def rref(mat,precision=0,GJ=False):
    m,n = mat.shape
    p,t = precision, 1e-1**precision
    A = around(mat.astype(float).copy(),decimals=p )
    if GJ:
        A = hstack((A,identity(n)))
    pcol = -1 #pivot colum
    for i in range(m):
        pcol += 1
        if pcol >= n : break
        #pivot index
        pid = argmax( abs(A[i:,pcol]) )
        #Row exchange
        A[i,:],A[pid+i,:] = A[pid+i,:].copy(),A[i,:].copy()
        #pivot with given precision
        while pcol < n and abs(A[i,pcol]) < t:
            pcol += 1
            if pcol >= n : break
            #pivot index
            pid = argmax( abs(A[i:,pcol]) )
            #Row exchange
            A[i,:],A[pid+i,:] = A[pid+i,:].copy(),A[i,:].copy()
        if pcol >= n : break
        pivot = float(A[i,pcol])
        for j in range(m):
            if j == i: continue
            mul = float(A[j,pcol])/pivot
            A[j,:] = around(A[j,:] - A[i,:]*mul,decimals=p)
        A[i,:] /= pivot
        A[i,:] = around(A[i,:],decimals=p)
        
    if GJ:
        return A[:,:n].copy(),A[:,n:].copy()
    else:
        return A

这里有一些简单的测试

print("/*--------------------------------------/")
print("/             Simple TEST               /")
print("/--------------------------------------*/")
A = array([[1,2,3],[4,5,6],[-7,8,9]])
R = rref(A,precision=6)
print("A:\n",A)
print("R:\n",R)
print()
print("With GJ ")
R,E =   rref(A,precision=6,GJ=True)
print("R:\n",R)
print("E:\n",E)
print("AdotE:\n",around( dot(A,E), decimals=0))
print()
A = array([[0, 0, 1], [0, 1, 0]]) 
R = rref(A, precision=1)
print("A:\n",A)
print("R:\n",R)
print()
A = array([[1,2,2,2],[2,4,6,8],[3,6,8,10]])
R = rref(A,precision=6)
print("A:\n",A)
print("R:\n",around(R, decimals=0))
print()

print("/*--------------------------------------/")
print( "/        Not Invertable TEST            /")
print( "/--------------------------------------*/")
A = array([
    [2,2,4, 4],
    [3,1,6, 2],
    [5,3,10,6]])
R = rref(A,precision=2)
print("A:\n",A)
print("R:\n",R)
print()
print("A^{T}:\n",A.T)
R = rref(A.T,precision=10)
print("R:\n",R)

/*--------------------------------------/
/             Simple TEST               /
/--------------------------------------*/
A:
 [[ 1  2  3]
  [ 4  5  6]
  [-7  8  9]]
R:
 [[1. 0. 0.]
  [0. 1. 0.]
  [0. 0. 1.]]

With GJ 
R:
 [[1. 0. 0.]
  [0. 1. 0.]
  [0. 0. 1.]]
E:
 [[-0.071428  0.142857 -0.071429]
  [-1.857142  0.714285  0.142857]
  [ 1.595237 -0.523809 -0.071428]]
AdotE:
 [[ 1.  0.  0.]
  [ 0.  1.  0.]
  [-0.  0.  1.]]

A:
 [[0 0 1]
  [0 1 0]]
R:
 [[0. 1. 0.]
  [0. 0. 1.]]

A:
 [[ 1  2  2  2]
  [ 2  4  6  8]
  [ 3  6  8 10]]
R:
 [[ 1.  2.  0. -2.]
  [ 0.  0.  1.  2.]
  [ 0.  0.  0.  0.]]

/*--------------------------------------/
/        Not Invertable TEST            /
/--------------------------------------*/
A:
 [[ 2  2  4  4]
  [ 3  1  6  2]
  [ 5  3 10  6]]
R:
 [[ 1.  0.  2.  0.]
  [-0.  1. -0.  2.]
  [ 0.  0.  0.  0.]]

A^{T}:
 [[ 2  3  5]
  [ 2  1  3]
  [ 4  6 10]
  [ 4  2  6]]
R:
 [[ 1.  0.  1.]
  [-0.  1.  1.]
  [ 0.  0.  0.]
  [ 0.  0.  0.]]

Answer 3

如果你可以使用 sympy 这个库的话，Matrix.rref() 这个功能可以帮你完成这个任务：

In [8]: sympy.Matrix(np.random.random((4,4))).rref()
Out[8]: 
([1, 1.42711055402454e-17, 0, -1.38777878078145e-17]
[0,                  1.0, 0,  2.22044604925031e-16]
[0, -2.3388341405089e-16, 1, -2.22044604925031e-16]
[0, 3.65674099486992e-17, 0,                   1.0],
 [0, 1, 2, 3])

Python内置函数进行矩阵化简

6 个回答

撰写回答