Computes the jacobian of a real-valued function, using the complex step method. The result is accurate nearly down to the machine epsilon with a very simple code, at the expense of uncessary computations, and assuming that the function to differentiate is defined in the complex domain

complex-step-jacobian.py

def get_jacobian(func, X, h = 1e-7):
	J = []	
	for i in range(X.shape[0]):
		Xp = X.astype(complex)
		Xp[i] += h * 1j	
		J.append([numpy.imag(func(Xp))])

	J = numpy.concatenate(J, axis = 0)
	J = numpy.transpose(J)
	J /= h
	return J