Czy Python ma wbudowaną funkcję, która konwertuje matrycę w wiersz Echelon Form (znany również jako górny trójkątny)?

31
user968102 5 październik 2011, 20:00

5 odpowiedzi

Najlepsza odpowiedź

Jeśli możesz użyć sympy, Matrix.rref() może to zrobić:

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])
38
NPE 5 październik 2011, 16:18

Zobacz http://mail.scipy.org/piperMail/numpy- Dyskusja / 2008 - listopad / 038705.html

Zasadniczo: nie rób tego.

Algorytm RREF wytwarza zbyt wiele niedokładności przy wdrożeniu na komputerze. Więc chcesz rozwiązać problem w inny sposób lub sugerować symboliczne symboliczne, podobnie jak @aix.

5
Winston Ewert 5 październik 2011, 17:31

Tak. W scipy.linalg

Istnieją inne czynniki, takie jak qr, rq, svd i więcej, jeśli jesteś zainteresowany.

Dokumentacja.

4
Steve Tjoa 5 październik 2011, 18:43

Zgadzam się z komentarzem @mile do @Winstonewert Odpowiedź nie ma powodu, dla którego komputer nie mógł wykonać RREF z danej precyzji .

Realizacja RREF nie powinna być bardzo skomplikowana, a Matlab jakoś udało się mieć tę funkcję, więc Numpy też powinno mieć.

Zrobiłem bardzo prostą i prostą realizację, która jest bardzo nieefektywna. Ale dla prostych matryc działa całkiem dobrze:

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   

Tutaj kilka prostych testów

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.]]
4
rth 28 wrzesień 2019, 03:44

Możesz go zdefiniować:

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
Peter Francis 21 styczeń 2020, 03:09