As pointed out by@pmli in#27, thecurrent implementation involves unnecessarily computing the characteristic polynomial, and then finding its roots. This is inefficient and numerically unstable. Much better to usenp.linalg.eigvals.