【发布时间】：2011-12-01 15:07:42
【问题描述】：
python有没有将矩阵转换为行梯形形式（也称为上三角）的内置函数？
【问题讨论】：
【发布时间】：2011-12-01 15:07:42
【问题描述】：
如果你可以使用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])
【讨论】：
-
如何解释结果？
-
4x4 矩阵是您在 RRE 中的矩阵（注意浮点精度），1x4 矩阵列出了您的枢轴变量的索引。
-
手动，约 5-10 分钟。同情，1.13 毫秒：-3
-
感谢您的有用回答。不过，转换回 numpy 数组并不是那么简单。我用
np.array(sympy.lambdify((),result[0])())。也许有更好的方法？ -
没关系，
np.array(result[0].tolist(), dtype=float)更简单。
见http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038705.html
基本上：不要这样做。
rref 算法在计算机上实施时会产生过多的不准确性。因此，您要么想以其他方式解决问题，要么使用 @aix 建议的符号。
【讨论】：
-
只需将您的数字包装在一个有理数据类型中，例如 python 的
fraction.Fraction，它已经为算术运算重载。你甚至可以用 numpy 向量化构造函数，即VecFraction = np.vectorize(fraction.Fraction)和rational_array = VecFraction(numeric_array)。计算机没有理由不能执行精确的 RREF。 -
答案中的链接对我不起作用，numpy-discussion.10968.n7.nabble.com/…
-
“或使用@aix 建议的符号”，这是在哪里写的？
-
如果 rref 在计算机上不准确，那么为什么我们要在 scipy 上实现类似的算法 LU 以供生产使用？
【讨论】：
-
LU 分解与行梯形不同。见math.stackexchange.com/a/1614952
-
是的，可以从这些函数计算减少的行梯队，但为什么要让用户跳过箍。大多数明智的应用程序，如 scilab、Matlab、Mathematica 等，都内置了这些。相反，Python 的使用必须通过符号代数库。
我同意@Mile 的评论to @WinstonEwert answer 没有理由计算机不能以给定的精度 执行 RREF。
RREF 的实现应该不会很复杂，matlab 莫名其妙地有这个功能，所以 numpy 应该也有。
我做了一个很简单直接的实现，效率很低。但对于简单的矩阵，它工作得很好：
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 xrange(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:
#pivot index
pid = argmax( abs(A[i:,pcol]) )
#Row exchange
A[i,:],A[pid+i,:] = A[pid+i,:].copy(),A[i,:].copy()
pcol += 1
if pcol >= n : break
pivot = float(A[i,pcol])
for j in xrange(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]])
print "A:\n",R
R = rref(A,precision=6)
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)
/*--------------------------------------/
/ Simple TEST /
/--------------------------------------*/
A:
[[ 1. 0. 1.]
[-0. 1. 1.]
[ 0. 0. 0.]
[ 0. 0. 0.]]
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.]]
print "/*--------------------------------------/"
print "/ Not Invertable TEST /"
print "/--------------------------------------*/"
A = array([
[2,2,4, 4],
[3,1,6, 2],
[5,3,10,6]])
print "A:\n",A
R = rref(A,precision=2)
print "R:\n",R
print
print "A^{T}:\n",A.T
R = rref(A.T,precision=10)
print "R:\n",R
/*--------------------------------------/
/ 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.]]
【讨论】：
你可以自己定义：
def rref(matrix):
A = np.array(matrix, dtype=np.float64)
i = 0 # row
j = 0 # column
while True:
# find next nonzero column
while all(A.T[j] == 0.0):
j += 1
# if reached the end, break
if j == len(A[0]) - 1 : break
# if a_ij == 0 find first row i_>=i with a
# nonzero entry in column j and swap rows i and i_
if A[i][j] == 0:
i_ = i
while A[i_][j] == 0:
i_ += 1
# if reached the end, break
if i_ == len(A) - 1 : break
A[[i, i_]] = A[[i_, i]]
# divide ith row a_ij to make it a_ij == 1
A[i] = A[i] / A[i][j]
# eliminate all other entries in the jth column by subtracting
# multiples of of the ith row from the others
for i_ in range(len(A)):
if i_ != i:
A[i_] = A[i_] - A[i] * A[i_][j] / A[i][j]
# if reached the end, break
if (i == len(A) - 1) or (j == len(A[0]) - 1): break
# otherwise, we continue
i += 1
j += 1
return A
【讨论】：
-
这里有几个错误，但我仍然觉得这个过程很优雅： 1. if j == len(A[0]) - 1 应该突破外部而 2. 你应该在规范化 3 之前检查 A[i][j] 是否不为零（理想情况下不为零）。如果 A[i][j] 不为零，您应该只乘其他行
这是一个工作版本，它几乎只是 MATLAB 的 rref 函数的 numpy 版本：
def rref(A, tol=1.0e-12):
m, n = A.shape
i, j = 0, 0
jb = []
while i < m and j < n:
# Find value and index of largest element in the remainder of column j
k = np.argmax(np.abs(A[i:m, j])) + i
p = np.abs(A[k, j])
if p <= tol:
# The column is negligible, zero it out
A[i:m, j] = 0.0
j += 1
else:
# Remember the column index
jb.append(j)
if i != k:
# Swap the i-th and k-th rows
A[[i, k], j:n] = A[[k, i], j:n]
# Divide the pivot row i by the pivot element A[i, j]
A[i, j:n] = A[i, j:n] / A[i, j]
# Subtract multiples of the pivot row from all the other rows
for k in range(m):
if k != i:
A[k, j:n] -= A[k, j] * A[i, j:n]
i += 1
j += 1
# Finished
return A, jb
例子：
A = np.array([[16.0, 2, 3, 13], [5, 11, 10, 8],
[9, 7, 6, 12], [4, 14, 15, 1]])
Areduced, jb = rref(A)
print(f"The matrix as rank {len(jb)}")
print(Areduced)
【讨论】：