Comaprsion between Natural Gradient Descent and Stochastic Gradient Descent Natural Gradient Descent Algorithm NGD was introdued by Amari. It has been used in Deep Networks and in ICA.
Although in this Comparison of NGD with GD they discussed the comparison between the three algorithms. 
Can somebody explain the difference between the two algorithm using a small numerical example on some data in Python and discuss how they perform in detail?
 A: The example task: Given the $(x,y)$, compute the equivalent polar coordinates $(r, \theta)$ using (natural) gradient descent. This makes a good candidate for showing off NGD, because we are concerned about non-euclidean parameter spaces.
First some boilerplate code. Our target is $(x,y) = (1,1)$.
import numpy as np
from scipy.optimize import approx_fprime as gradient
import matplotlib.pyplot as plt

xy = np.array([1.0, 1.0])

def xy_for_rtheta(rtheta):
    r = rtheta[0]
    theta = rtheta[1]
    return np.array([r*np.cos(theta), r*np.sin(theta)])

We're going to attempt to minimize the squared error     
def err(rtheta):
    x, y = xy_for_rtheta(rtheta)
    dx = x-xy[0]
    dy = y-xy[1]
    return dx*dx + dy*dy

Here's the function which computes the gradients. For a derivation of the riemannian metric tensor $G$ see this earlier paper. Note that I normalized the gradients.
def compute_grad(rtheta, natural = True):
    grad = gradient(rtheta, err, epsilon = 1E-6)
    if natural:
        G = np.array([[1.0, 0.0],[0.0, rtheta[0]**2]])
        grad = np.matmul(np.linalg.inv(G), grad)
    return grad / np.linalg.norm(grad)

Note that the metric can also be computed using the jacobian, so alternatively we could write:
    J = np.stack([gradient(rtheta, lambda rt: xy_for_rtheta(rt)[0], epsilon = 1E-6),
                 gradient(rtheta, lambda rt: xy_for_rtheta(rt)[1], epsilon = 1E-6)], axis = 0)
    G = np.matmul(J.T, J)

Then a straightforward implementation of gradient descent, starting with $(r,\theta) = (0.5, 0.75\ \pi)$:
def descend(color, natural):
    rtheta = np.array([0.5, np.pi*3/4.])
    rthetas = [rtheta]
    stepsize = 0.001
    tol = 0.0001

    while err(rtheta) > tol:
        rtheta = rtheta - stepsize * compute_grad(rtheta, natural)
        rthetas.append(rtheta)

    xs, ys = zip(*map(xy_for_rtheta, rthetas))
    plt.plot(xs, ys, c = color)

Finally we run the code: In the red is ordinary gradient descent and in the blue is natural gradient descent. While ordinary GD takes a winding path, due to the particular parameterization of problem (polar), natural gradient descent will make a beeline for the target no matter which coordinate system is used.

if __name__ == '__main__':
    descend('r', False)
    descend('b', True)
    plt.show()

What happens if the mapping is changed to $x = 2r \cos(\theta), y = r\sin(\theta)$?

Even in a very simple setting such as $x' = 2x, y' = y$, we can observe behaviour like this:

The reason for this behaviour is that normal GD aims for the direction of steepest descent in parameter space, whereas NGD aims to find the best direction on the underlying manifold. 
When you're optimizing a neural network with NDG, each point in our search space is a statistical model, and NDG tries to take a straight line to the best statistical model regardless of what "coordinate system" (parameterization) that model uses.
