This work is concerned with numerical methods for matrix eigenvalue problems that are nonlinear in the eigenvalue parameter. In particular, we focus on eigenvalue problems for which the evaluation of the matrix-valued function is computationally expensive. Such problems arise, e.g., from boundary integral formulations of elliptic PDE eigenvalue problems and typically exclude the use of established nonlinear eigenvalue solvers. Instead, we propose the use of polynomial approximation combined with non-monomial linearizations. Our approach is intended for situations where the eigenvalues of interest are located on the real line or, more generally, on a pre-specified curve in the complex plane. A first-order perturbation analysis for nonlinear eigenvalue problems is performed. Combined with an approximation result for Chebyshev interpolation, this shows exponential convergence of the obtained eigenvalue approximations with respect to the degree of the approximating polynomial. Preliminary numerical experiments demonstrate the viability of the approach in the context of boundary element methods.