scipy.stats.multivariate_normal.pdf is not robust enough to handle singular covariances. In general this package takes care of handling many numerical pitfalls accurately (often in a different way than sklearn's implementation of GMMs by the way) and improved a lot in recent updates. This is a crucial aspect for a robust GMR implementation and it is often not directly supported by scipy and numpy.
pinvh is part of scipy (now?), it has been replaced
... as suggested in https://github.com/openjournals/joss-reviews/issues/3054
Some findings: