I use the following notation and definition for the (first) variation of some functional $E[\psi]$: \begin{equation} \delta E[\psi][\delta\psi] := \lim_{\varepsilon \rightarrow 0} \frac{E[\psi + \varepsilon \delta\psi] - E[\psi]}{\varepsilon}. \end{equation} For a Hamiltonian that depends on some parameter, $H(\lambda)$, which fulfills the eigenequation $H(\lambda) |\psi(\lambda)\rangle = E(\lambda) |\psi(\lambda)\rangle$, one can quite easily show that \begin{equation} \frac{dE(\lambda)}{d\lambda} = \langle \psi(\lambda) | \frac{dH(\lambda)}{d\lambda} | \psi(\lambda) \rangle \end{equation} which is called the Hellmann-Feynman theorem. Now in quantum chemistry, one often seeks an approximation of the true eigenfunctions by making the functional for the energy expectation value stationary, \begin{equation} E[\psi] = \frac{\langle \psi | H | \psi \rangle}{\langle \psi | \psi \rangle} \end{equation} \begin{equation} \delta E[\psi][\delta \psi] = \frac{1}{\langle \psi | \psi \rangle}\left[ \langle \delta\psi | H - E[\psi] | \psi\rangle + \langle \psi | H - E[\psi] | \delta\psi\rangle \right] = 0 \end{equation} for all variations $\delta\psi$ in some subspace (called the variational space) of the full Hilbert space. In this case, one often reads that the Hellmann-Feynman theorem is still valid for the approximate wave functions and energies.
Unfortunately, I have problems finding a clear proof of this statement in books. Wikipedia tries a proof (section "Alternate proof"), but it assumes that the functional derivative \begin{equation} \frac{\delta E[\psi]}{\delta\psi} \end{equation} exists, which is not the case because of the antilinearity in bra vectors, which leads to \begin{equation} \delta E[\psi][\alpha \cdot \delta \psi] \neq \alpha \cdot \delta E[\psi][ \delta \psi] \end{equation} for complex $\alpha$. Can anyone provide me with a clear justification of why the Hellmann-Feynman theorem works for variational wavefunctions?