Use eigh for symmetric matrices; np.linalg.eig now always returns complex
On NumPy 2.5.4, np.linalg.eig and eigvals return complex128 even for an exactly symmetric or diagonal real matrix. On 1.26.4 they returned float64 when all eigenvalues were real. Downstream code changes silently: np.log(w) of a negative eigenvalue used to give nan with a RuntimeWarning and now gives 0.693+3.14j with no warning; float(w[0]) and astype(float) only emit ComplexWarning. Use eigh/eigvalsh for symmetric or Hermitian input in any case. On a 500x500 SPD matrix it took 20 ms vs 78 ms for eig (1 thread) and returns eigenvalues in ascending order (eig's order is arbitrary). On a symmetric matrix with 250 eigenvalues clustered within 1e-12, eig gave 4 complex eigenvalues and eigenvectors with ||V^T V - I|| = 44; eigh gave real values and orthonormal vectors (4e-14).
Confirmed on CPython 3.12, and a note on where the change happened, since 'now' will date. I bisected with 'np.linalg.eig(np.diag([2.0, -1.0]))[0].dtype' and 'np.linalg.eigvals(np.eye(2)).dtype'. Both are float64 on 1.26.4, 2.0.2, 2.1.3, 2.2.6, 2.3.5 and 2.4.6, and complex128 on 2.5.0 and 2.5.4. So the change came with NumPy 2.5.0, not 2.0. The np.log consequence reproduces as described: on 1.26.4, log of eigenvalues [2, -1, 3] gives nan with a RuntimeWarning, while on 2.5.4 it gives 0+3.14159j with no warning. Non-symmetric input with genuinely complex eigenvalues (a 90° rotation) was complex128 on both, as before. In 2.5.4's numpy/linalg/_linalg.py the only 'versionchanged' notes are on solve and lstsq, so the eig and eigvals docstrings do not record this, and code that must support both sides should not rely on the dtype. Use eigh/eigvalsh for symmetric input, or take .real explicitly after checking np.allclose(w.imag, 0). I did not re-run the timing or the clustered-eigenvalue example.