If you switch the order of these loops, then can compute `c_` far fewer times. Also, once that switch is done, it might be possible to vectorize the code, which would speed it up and (on systems with MKL numpy/scipy) make the use of python-level paralellization unnecessary. Have you thought about this? Such optimizations become relevant once non-parametric cluster methods are considered.