Gaussian Elimination with Pivots
For some inputs where there is no solution, gaussian_reduce returns an array of NaNs:

&gt;&gt;&gt; gaussian_reduce(np.zeros((2, 2)), np.ones((2, 1)))
__main__:63: RuntimeWarning: invalid value encountered in double_scalars
array([[nan, nan]])


A better way to handle an exceptional situation is to raise an exception, like np.linalg.solve:

&gt;&gt;&gt; np.linalg.solve(np.zeros((2, 2)), np.ones((2, 1)))
Traceback (most recent call last):
  File "&lt;stdin&gt;", line 1, in &lt;module&gt;
  File "numpy/linalg/linalg.py", line 390, in solve
    r = gufunc(a, b, signature=signature, extobj=extobj)
  File "numpy/linalg/linalg.py", line 89, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.linalg.LinAlgError: Singular matrix


It looks as though you attempted to do this:

matrix = np.append(matrix, b, axis=1)
matrix = to_row_echelon(matrix)
if all(matrix[-1,:] == 0):
    raise Exception("System of equations is inconsistent")


but it doesn't work, because in this case (a zero matrix) the last row gets filled with NaNs rather than zeros:

&gt;&gt;&gt; m = np.append(np.zeros((2, 2)), np.ones((2, 1)), axis=1)
&gt;&gt;&gt; to_row_echelon(m)
array([[ 0.,  0.,  1.],
       [nan, nan, nan]])


The NaNs result from the division in subtract_rows:

return row - (row[0]/top_row[0])*top_row


This division needs to be skipped if top_row[0] is zero. For example, in pivot you would have:

if matrix[0, 0]:


before the call to np.apply_along_axis. (But see below for further improvements here.)
In mathematical code, you should be on the lookout for division by zero. Every division operation is suspicious: what if the divisor is zero? There are some more division operations in backsolve which lead to this kind of incorrect result:

&gt;&gt;&gt; gaussian_reduce(np.ones((2, 2)), np.array([[2], [3]]))
__main__:70: RuntimeWarning: divide by zero encountered in double_scalars
array([[-inf, inf]])


In this case you need to raise an appropriate exception, again like  np.linalg.solve:

&gt;&gt;&gt; np.linalg.solve(np.ones((2, 2)), np.array([[2], [3]]))
Traceback (most recent call last):
  File "&lt;stdin&gt;", line 1, in &lt;module&gt;
  File "numpy/linalg/linalg.py", line 390, in solve
    r = gufunc(a, b, signature=signature, extobj=extobj)
  File "numpy/linalg/linalg.py", line 89, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.linalg.LinAlgError: Singular matrix

gaussian_reduce returns the result as a row vector (an array with one row):

&gt;&gt;&gt; A = np.diag((1., 2.))
&gt;&gt;&gt; b = np.array([[4.], [5.]])
&gt;&gt;&gt; gaussian_reduce(A, b)
array([[4. , 2.5]])


But in the equation \$Ax = b\$ the solution \$x\$ is a column vector (an array with one column). And this is what np.linalg.solve returns:

&gt;&gt;&gt; np.linalg.solve(A, b)
array([[4. ],
       [2.5]])


Returning a column vector means that we can verify the solution by multiplying:

&gt;&gt;&gt; A @ np.linalg.solve(A, b)
array([[4.],
       [5.]])


which is equal to b as required. Whereas in the gaussian_reduce case, it doesn't work:

&gt;&gt;&gt; A @ gaussian_reduce(A, b)
Traceback (most recent call last):
  File "&lt;stdin&gt;", line 1, in &lt;module&gt;
ValueError: shapes (2,2) and (1,2) not aligned: 2 (dim 1) != 1 (dim 0)

The name gaussian_reduce could be improved. It's better to name a function based on what it computes, not after the algorithm it uses to compute it. For example, if you want to find the strongly connected components of a graph, you'd rather have a function named strongly_connected_components than one named tarjans_algorithm (Robert Tarjan invented many algorithms). Naming a function after its result rather than its method also allows you to change the method if you find a better one.

Since gaussian_reduce is a function that returns a solution to a system of linear equations, a name like linear_system_solution (or some suitable abbreviation) would be clearer.
When writing NumPy code, it's a good idea to be clear about when you are modifying arrays in-place and when you are returning new arrays. If you don't keep this distinction clear, then you are at risk of either (i) modifying an array that you didn't mean to, for example one passed by the caller; or (ii) making an unnecessary copy of an array when you could have just modified the original.

There are a couple of things you can do to keep this distinction clear. First, distinguish between functions that construct new arrays and functions that modify their input arrays, by having the former return the new array, and the latter return nothing (consider the difference between Python's built-in sorted, which returns a new list, and list.sort, which returns nothing). In this case, since to_row_echelon modifies its input, it should return nothing.

Second, avoid reusing a variable name to refer to a newly constructed array, because this will lead to confusion about whether it is safe to modify the array. For example, in gaussian_reduce you have:

matrix = np.append(matrix, b, axis=1)


The function np.append always constructs a new array, and so the original matrix and b are unchanged. This means that the copy and the deletion are unnecessary in the following lines:

b = matrix[:,-1].copy()
matrix = np.delete(matrix, -1, axis =1)
b = backsolve(matrix, b)


So to improve this code I would give the full matrix a new variable name to avoid confusion:

augmented_matrix = np.append(matrix, b, axis=1)


Then since to_row_echelon modifies its input, it should not return a result:

to_row_echelon(augmented_matrix)


Then the copy and deletion are unnecessary, because I can just index into augmented_matrix:

solution = backsolve(augmented_matrix[:,:-1], augmented_matrix[:, -1])


Note the use of the new variable name solution to avoid confusion with b.

Finally, instead of constructing yet another array:

return np.array([b])


we could reshape the array we already have:

return solution.reshape(1, -1)


but in fact as discussed above we should return a column vector here, not a row vector, so we actually need:

return solution.reshape(-1, 1)

Similarly, pivot modifies the input and so should not also return it.
The docstring for pivot says:

... to set all other elements of the first row to be zero


but this should say "column", not "row".
pivot needs to find the row with the largest absolute value in the first column, which it does like this:

idx_of_max = np.argmax(list(map(abs, matrix.T[0])))


It is clearer, I think, to get the first column of the matrix using matrix[:, 0] instead of transposing. Then you can just call abs(matrix[:, 0]), which returns a new array containing absolute values element-wise. This avoids the need for map and list.
Instead of calling subtract_rows for each row:

matrix[1:,] = np.apply_along_axis(subtract_rows, 1, matrix[1:], matrix[0])


do all the subtractions in one go, avoiding the need for subtract_rows:

matrix[1:] -= (matrix[0] / matrix[0, 0]) * matrix[1:, 0]


But to avoid division by zero you will need to guard this, as discussed above:

pivot_value = matrix[0, 0]
if pivot_value:
    matrix[1:] -= (matrix[0] / pivot_value) * matrix[1:, 0]

In swap_rows, it is misleading to name the arguments a and b, since elsewhere in the code you are using b for a column vector, but here b needs to be an index. It helps the reader to follow the code if you are consistent about the use of variable names. In this case I would use the names i and j.
swap_rows can be implemented more simply like this, using NumPy's "fancy indexing":

def swap_rows(a, i, j):
    "Swap rows i and j in the array a."
    a[[i, j]] = a[[j, i]]


This is now so simple that you would probably want to inline it at its single point of use.