Given a Hermitian operator $H$, I can calculate the variance of the Hamiltonian ${\rm Var}[H]$ as
$$ {\rm Var}[H] = \langle H^2 \rangle -\langle H \rangle^2 $$
Now, $H$ can be decomposed as
$$ H = \sum_i \alpha_i P_i \hspace{1 cm} P_i \in \{ I,X,Y,Z\}^{\otimes n} $$ so by linearity $\langle H \rangle = \sum_i \alpha_i \langle P_i \rangle $. In the Variational Quantum Eigensolver (VQE), we want to minimize the above expression, that is, $\min \langle H \rangle_{|\psi (\theta) \rangle} $.
Note that if $|\psi (\theta)\rangle$ is indeed an eigenvector of $H$ then ${\rm Var}[H] = 0 $. And so one can verify whether they have indeed reached/found an eigenstate by performing such calculation.
The question is how do we calculate $\langle H^2 \rangle$?
Yes, I can decomposed $H^2$ into sum of Pauli terms, that is
$$ H^2 = \sum_k \beta_k P_k \hspace{1 cm} P_i \in \{ I,X,Y,Z\}^{\otimes n} $$
and then using the optimized $|\psi (\theta) \rangle$ I have found through VQE to perform $$\langle H^2 \rangle = \langle \psi(\theta) | \sum_k \beta_k P_k | \psi(\theta) \rangle $$
This would be the same process as I did to calculate $\langle H \rangle$, applying the appropriate rotations to my circuit before measuring in the $Z$ basis... Then once I have this, $\langle H^2 \rangle$, I can put it together with the expectation value I have found through VQE to calculate ${\rm Var}[H]$. This is the naive way of doing it I guess...
What is a better way to do this?
