Hi, nice and efficient piece of code.
Nevertheless, I believe that the polynomial expression that you compute in a first step is not guaranteed to be one of an ellipsoid. I think that it could be any quadric (hyperboloid, paraboloid, etc), right?
It is probably no problem if you have much data (corresponding to a real ellipsoid) to fit your quadric on.
But in case of sparse data (for instance, data points only on a couple of planes), you could have bad surprises (and in particular negative eigenvalues).