[Performance] Optimise calculation of eigenvalues/vectors in finite shot sampling

[atiyo]

Instead of calculating eigenvalues and eigenvectors of Hamiltonians acting on the entire state vector, we can instead calculate them on the matrix representation of the observable on a limited support and then kron to get the correct representation.

---

This issue it to investigate the performance implementations of this and implement the optimised version if necessary. It's plausible that the performance benefit might be small since XLA compiler might cache the result of the eigh call.

...isn't it just the eigenvector matrix of the small hamiltonian kronned with the identity matrix or something? At least it seems to me like you can get a valid eigendecomposition that way:

$$ H_\text{big} = H_\text{small} \otimes I = \left(\sum_\lambda \lambda|\lambda\rangle\langle\lambda|\right)\otimes I $$

where each eigenvalue of $H_\text{small}$ will show up repeated by the size of the identity. Check Theorem 13.12 here.

Originally posted by @jpmoutinho in #23 (comment)