I'm looking into the circuit for the VQE, but am stumped at how we can identify the expectation value of the Pauli series. Essentially, how do we find:
$$ \langle \psi | H_i | \psi \rangle $$
Given $H_i$ can be decomposed into Paulis?
The first strategy that comes to mind is simply applying the adjoint of the state prep routine and finding the average number of 0s measured, but this doubles the number of state prep calls necessary. Is there an alternative circuit I'm missing?
Let's denote $P$ as a Pauli tensor product term (Pauli term: e.g. $Z \otimes Z$, $X \otimes X$ or $Y \otimes X \otimes X$), then note that $P$ is a Hermitian matrix whose expectation value is real, the eigenvalues are $+1$ and $-1$ with corresponding eigenspaces. The expectation value of $P$ will be equal to:
$$ \langle \psi | P |\psi\rangle = p_{+} - p_{-}$$
where $p_{+}$ is the probability of $|\psi\rangle$ being in the $+1$ eigenspace and $p_{-}$ is the probability of $|\psi\rangle$ being in the $-1$ eigenspace of $P$. We will prove this. Note that $p_{+}$ and $p_{-}$ can be calculated by measuring in $P$ eigenbasis: $p_{+} = \frac{N_{+}}{N}$ and $p_{-} = \frac{N_{-}}{N}$, where $N_+$ is the number of measured eigenvectors that have $+1$ eigenvalue, $N_-$ is the number of measured eigenvectors that have $-1$ eigenvalue and $N$ is the total number of measurements.
Now the prove. Arbitrary $|\psi \rangle$ can be expressed as superposition of eigenbasis vectors of arbitrary Hermitian matrix (but we are interested in $P$) that acts in the same Hilbert space, so:
$$|\psi\rangle = \sum_{j_+} a_{j_+}|j_+\rangle + \sum_{j_-} b_{j_-}|j_-\rangle$$
where $a_{j_+}$ and $b_{j_-}$ are amplitudes of eigenvectors $|j_+\rangle$ and $|j_-\rangle$ of the $P$. We have separated the sum for our convenience: $|j_+\rangle$ are eigenvectors that have $+1$ eigenvalue (that means $P|j_+\rangle = |j_+\rangle$) and $|j_-\rangle$ are eigenvectors that have $-1$ eigenvalue (that means $P|j_-\rangle = -|j_-\rangle$). So:
$$P |\psi\rangle = \sum_{j_+} a_{j_+}|j_+\rangle - \sum_{j_-} b_{j_-}|j_-\rangle$$
By taking into account that all eigenbasis vectors $|j_+\rangle$, $|j_-\rangle$ are orthogonal to each other (e.g. $\langle j_+| k_+\rangle = 0$ when $j \ne k$ or $\langle j_+| k_-\rangle = 0$ ):
$$\langle \psi |P |\psi\rangle = \sum_{j_+} |a_{j_+}|^2 - \sum_{j_-} |b_{j_-}|^2= p_{+} - p_{-}$$
In this answer, one can find Qiskit implementation for $Z\otimes Z \otimes ... \otimes Z$ Pauli term. The answer also includes a discussion on how one can use the same procedure for calculating expectation value for an arbitrary Pauli term. Also, note that if the Hamiltonian $H = \sum_i c_i P_i$, where $P_i$ are Pauli terms and $c_i$ are real numberes, then $\langle \psi | H| \psi \rangle = \sum_i c_i \langle \psi | P_i | \psi \rangle$.
Replying to the comments of the answer
In Microsoft's docs it was proposed (not only) an approach of measuring $\langle ZZ \rangle$ and the question was how one can generalize their approach. Actually the generalization of their approach is not obvious to me and I suggest asking a separate question focused on generalizing their approach. Nevertheless, I want to mention that their approach is not the only way of finding $\langle ZZ \rangle$ and here are three alternative ways for doing it:
If the measured bitstring has odd (even) $1$s then we have measured the state in $-1$ ($+1$) eigenspace. The generalized Qiskit implementation of this approach can be found in this answer. Also, this approach is used (if I have understood the paper right) in this paper from IBM (Fig. 1 shows not only measurements but also simple one qubit unitary transformations ($I$ or $X_{-\frac{\pi}{2}}$ or $Y_{-\frac{\pi}{2}}$) that are specific for each Pauli term).
If we measure the ancillary qubit in $| 0\rangle$ ($|1\rangle$) state then we have measured it in $+1$ ($-1$) eigenspace. This circuit is also used for the Hadamard test / scattering circuit. Note that, Pauli terms have real expectation values, so the described circuit for calculating imaginary part of the expectation value in the Hadamard test is not relevant here. The generalization of this approach is fairly simple: for example, this is the circuit for $\langle XZZY \rangle$.
If measure the target qubit (targeted by CNOT) $|0\rangle$ ($|1\rangle$) then we have measured the state in the $+1$ ($-1$) eigenspace. Although its working, but the generalization of this approach is not obvious for me and I think this question about generalizing the approach shown in Microsoft's docs deserves a separate question.
P.S. There exist also other approaches for calculating $\langle ZZ \rangle$.
