I see. But since the grad calculation requires the eigenvectors -- and thus a call to `scipy.linalg.eigh(a, b, lower=self.lower)` -- I would also need to make an Op that wraps `scipy.linalg.eigh` that I would call inside `Eigvalsh.grad`, I think?

There is already an Op for numpy.linalg.eigh, but not one for scipy.linalg.eigh, which lets you pass the `b` matrix.

Pending a name for the scipy.linalg.eigh Op, I think I can do this, but I'm not sure if it's worth it.