CODEWITHSUNDEEP
pythonmatrixnumpy
flonk
flonk

Reputation: 3876

Is there a standard solution for Gauss elimination in Python?

Is there somewhere in the cosmos of scipy/numpy/... a standard method for Gauss-elimination of a matrix?

One finds many snippets via google, but I would prefer to use "trusted" modules if possible.

Upvotes: 35

Views: 55794

Answers (3)

rsant
rsant

Reputation: 73

You can use the symbolic mathematics python library sympy

import sympy as sp

m = sp.Matrix([[1,2,1],
           [-2,-3,1],
           [3,5,0]])

m_rref, pivots = m.rref() # Compute reduced row echelon form (rref).

print(m_rref, pivots)

This will output the matrix in reduced echelon form, as well as a list of the pivot columns

Matrix([[1, 0, -5],
        [0, 1,  3],
        [0, 0,  0]])

(0, 1)

Upvotes: 5

Bernardo Costa
Bernardo Costa

Reputation: 519

One function that can be worth checking is _remove_redundancy, if you wish to remove repeated or redundant equations:

import numpy as np
import scipy.optimize

a = np.array([[1.,1.,1.,1.],
              [0.,0.,0.,1.],
              [0.,0.,0.,2.],
              [0.,0.,0.,3.]])
print(scipy.optimize._remove_redundancy._remove_redundancy(a, np.zeros_like(a[:, 0]))[0])

which gives:

[[1. 1. 1. 1.]
 [0. 0. 0. 3.]]

As a note to @flonk answer, using a LU decomposition might not always give the desired reduced row matrix. Example:

import numpy as np
import scipy.linalg

a = np.array([[1.,1.,1.,1.],
              [0.,0.,0.,1.],
              [0.,0.,0.,2.],
              [0.,0.,0.,3.]])

_,_, u = scipy.linalg.lu(a)
print(u)

gives the same matrix:

[[1. 1. 1. 1.]
 [0. 0. 0. 1.]
 [0. 0. 0. 2.]
 [0. 0. 0. 3.]]

even though the last 3 rows are linearly dependent.

Upvotes: 8

flonk
flonk

Reputation: 3876

I finally found, that it can be done using LU decomposition. Here the U matrix represents the reduced form of the linear system.

from numpy import array
from scipy.linalg import lu

a = array([[2.,4.,4.,4.],[1.,2.,3.,3.],[1.,2.,2.,2.],[1.,4.,3.,4.]])

pl, u = lu(a, permute_l=True)

Then u reads

array([[ 2.,  4.,  4.,  4.],
       [ 0.,  2.,  1.,  2.],
       [ 0.,  0.,  1.,  1.],
       [ 0.,  0.,  0.,  0.]])

Depending on the solvability of the system this matrix has an upper triangular or trapezoidal structure. In the above case a line of zeros arises, as the matrix has only rank 3.

Upvotes: 42

Related Questions