python built-in function to do matrix reduction

Yes. In scipy.linalg, lu does LU decomposition which will essentially get you row-echelon form.

There are other factorizations such as qr, rq, svd, and more, if you're interested.

Documentation.

If you can use sympy, Matrix.rref() can do it:

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])

see http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038705.html

Basically: don't do it.

The rref algorithm produces too much inaccuracy when implemented on a computer. So you either want to solve the problem in another way, or use symbolics like @aix suggested.