inner = x_j.T.dot(theta)
            gradient = np.exp(inner)/np.sum(np.exp(inner)) - y_j
            hessian = (np.sum(np.exp(inner))*np.exp(inner) - (np.exp(inner))*(np.exp(inner)))/(np.sum(np.exp(inner)))*(np.sum(np.exp(inner)))
            hessianInverse = 1/hessian
            result = hessianInverse*gradient
            grad_step += np.outer(result, x_j).T / len(X)