Multivariable Gradient Descent in Numpy
Without sample inputs I can't run your whole code.  And I prefer not to guess.

The use of np.matrix suggests it was translated from MATLAB/Octave code.  That array subclass, in numpy, is always 2d, which makes it behave more like MATLAB matrices, especially old versions.  Transpose always has effect; row and column indexing returns 2d matrices; and * is matrix multiplication (as opposed to element wise, the .* of MATLAB).

I'll focus on the scaling function.  I don't see it being used, but it's simple and typical of the other functions.

import numpy as np
&lt;your code&gt;

X = np.arange(12).reshape(3,4)
print(X)
Xm = np.matrix(X)
print(scaling(Xm))    


produces:

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
[[-0.5 -0.5 -0.5 -0.5]
 [ 0.   0.   0.   0. ]
 [ 0.5  0.5  0.5  0.5]]


scaling doesn't work with X an array, because x = X[:, k] would be 1d, which conflicts with the x[:, 0] use.

But I can perform this same operation without row iteration.

def ascaling(X, n=0):
    # array scaling
    ctr = X.mean(axis=n, keepdims=True)
    rge = X.max(axis=n, keepdims=True)-X.min(axis=n, keepdims=True)
    return (X - ctr)/rge
print(ascaling(X,0))
print(ascaling(X,1))


producing

[[-0.5 -0.5 -0.5 -0.5]
 [ 0.   0.   0.   0. ]
 [ 0.5  0.5  0.5  0.5]]
[[-0.5        -0.16666667  0.16666667  0.5       ]
 [-0.5        -0.16666667  0.16666667  0.5       ]
 [-0.5        -0.16666667  0.16666667  0.5       ]]


In this case, the equivalent code, assumng X is np.matrix is simpler

def mscaling(X, n=0):
    # matrix scaling
    ctr = X.mean(axis=n)
    rge = X.max(axis=n)-X.min(axis=n)
    return (X - ctr)/rge

print(mscaling(Xm,0))
print(mscaling(Xm,1))


Maybe this example will make these operations clear:

Make a small 2d array:

In [861]: X=np.arange(6).reshape(2,3)
In [862]: X
Out[862]: 
array([[0, 1, 2],
       [3, 4, 5]])


Sum rows, resulting in a 1d array of length 3 (the columns)

In [863]: X.sum(axis=0)
Out[863]: array([3, 5, 7])


Sum columns - again a 1d result

In [864]: X.sum(axis=1)
Out[864]: array([ 3, 12])


But if I add keepdims, the result is 2d, (2,1) shape:

In [866]: X.sum(axis=1,keepdims=True)
Out[866]: 
array([[ 3],
       [12]])


sum applied to matrix does the same thing.

In [867]: np.matrix(X).sum(axis=1)
Out[867]: 
matrix([[ 3],
        [12]])


The nice thing about keeping dims is that I can do math like

In [868]: X-X.sum(axis=1,keepdims=True)
Out[868]: 
array([[-3, -2, -1],
       [-9, -8, -7]])


Without keepdims I'd have to do X-X.sum(axis=1)[:,None].

Your frequent use of transpose suggests that the mix of dimensions haven't been fully thought out.  MATLAB code does use x' or x.' quite a bit.  On the other hand beginner numpy coders try to apply transpose to 1d arrays, and wonder why nothing happens.