This can be implemented through broadcasting, to avoid the allocation of the whole diagonal matrices. This is faster:

    diag1 = 1. / np.sum(X, axis=1)
    diag2 = 1. / np.sum(X, axis=0)
    an = diag1[:, np.newaxis] * X * diag2

Or, with less operations but maybe more obscure, you can even broadcast directly with division:

    an = X / diag1[:, np.newaxis] / diag2

I'm not sure if our scalers should be reused for this as here the scaling is along both axes, but either way, I think this tip is cool :)
Also I don't know if this is explicitly written down but I am used to seeing np.dot(A, B) rather than A.dot(B).