  We investigate contributions of excited states to nucleon matrix elements
computed in lattice QCD by employing, in addition to the standard nucleon
interpolating operator, pion-nucleon ($\pi$-$N$) operators. We solve a
generalized eigenvalue problem (GEVP) to obtain an optimal interpolating
operator that minimizes overlap with the $\pi$-$N$ states. We derive a variant
of the standard application of the GEVP method, which allows for constructing
3-point correlation functions using the optimized interpolating operator
without requiring the computationally demanding combination that includes
$\pi$-$N$ operators in both sink and source. We extract nucleon matrix elements
using two twisted mass fermion ensembles, one ensemble generated using pion
mass of 346 MeV and one ensemble tuned to reproduce the physical value of the
pion mass. Especially, we determine the isoscalar and isovector scalar,
pseudoscalar, vector, axial, and tensor matrix elements. We include results
obtained using a range of kinematic setups, including momentum in the sink. Our
results using this variational approach are compared with previous results
obtained using the same ensembles and multi-state fits without GEVP
improvement. We find that for the physical mass point ensemble, the
improvement, in terms of suppression of excited states using this method, is
most significant for the case of the matrix elements of the isovector axial and
pseudoscalar currents.
