I am considering a two-grid finite element scheme for eigenvalue problems arising from PDEs. The toy problem I am addressing first is the 2D Laplace eigenvalue equation
$$\Delta u(\mathbf{x}) = \lambda u(\mathbf{x}), \quad \mathbf{x}\in\Omega$$ $$ u(\mathbf{x}) = 0, \quad \mathbf{x}\in\Gamma,$$
where $\Omega = [0,1]\times[0,1]$, with boundary denoted by $\Gamma$. I am using two uniform triangular meshes $\mathcal{T}_H$ and $\mathcal{T}_h$ (coarse and fine, respectively), such that $\mathcal{T}_h$ is obtained by one uniform refinement of $\mathcal{T}_H$. Denote the corresponding piecewise linear finite element spaces by $V_H$ and $V_h$, respectively. The associated matrix eigenvalue problems are then
$$A_H u_H = \lambda_H M_H u_H \quad \text{and} \quad A_h u_h = \lambda_h M_h u_h,$$
respectively. Here $A_H$ (resp. $A_h$) is the stiffness matrix given by:
$$ a_{ij} = a(\phi_i,\phi_j) = \int_{\Omega} \nabla\phi_i \nabla\phi_j dx dy,$$
and $M_H$ (resp. $M_h$) is the mass matrix given by:
$$ m_{ij} = (\phi_i,\phi_j) = \int_{\Omega} \phi_i \phi_j dx dy,$$
where $\{\phi_1,\dots,\phi_{N_H}\}$ (resp. $\{\phi_1,\dots,\phi_{N_h}\}$) are the standard linear basis functions.
For this particular problem $\{A_H,M_H\}$ and $\{A_h,M_h\}$ are symmetric positive-definite pencils, and consequently the eigenvectors of the coarse-grid problem are both $A_H$- and $M_H$-orthogonal (by which I mean, e.g., $u_{i,H}^{T} A_H u_{j,H} = 0$, whenever $i \neq j$), and similarly for the fine-grid problem.
Suppose I have two "exact" (machine precision) eigenpairs $(\lambda_{1,H},u_{1,H}) \in \mathbb{R}\times V_H$ and $(\lambda_{2,H},u_{2,H}) \in \mathbb{R}\times V_H$ for the coarse problem, so that $u_{1,H}$ and $u_{2,H}$ are $A_H$-orthogonal, and $\lambda_{i,H}$ taken to be the generalised Rayleigh quotient of $u_{i,H}$ with respect to $A_H$ and $B_H$. I want to construct two approximate eigenpairs, say $(\mu_{1,h},v_{1,h})\in \mathbb{R}\times V_h$ and $(\mu_{2,h},v_{2,h})\in \mathbb{R}\times V_h$, for the fine-grid problem by taking
$$ v_{i,h} = I(u_{i,H}) \quad i=1,2,$$
where $I(u_{i,H})$ denotes some interpolation of the coarse-grid solution onto the fine grid. $\lambda_{i,h}$ would then be taken to be the generalised Rayleigh quotient of $v_{i,h}$ with respect to $A_h$ and $B_h$.
Question: What assumptions are required on the interpolation method, in order to ensure that $v_{1,h}$ and $v_{2,h}$ are $A_h$-orthogonal?
My thoughts:
Having looked at the simple 1D scenario, it seems intuitive to me that linear interpolation will yield the required orthogonality, due to the fact that the coarse-grid solution is piecewise linear (so, $I(u_{i,h})$ is essentially $u_{i,H}$ sitting in $V_h$). In a similar manner, I can deduce that the Rayleigh quotient approximations $\mu_{i,h}$ are exactly equal to $\lambda_{i,H}$.
It does not seem obvious to me that the same will hold for higher degree interpolation. In particular, I am interested in using a cubic interpolant, as this gives far superior approximations, but my a priori error bounds are dependent on this orthogonality.
Having implemented all of the above in MATLAB, both linear and cubic interpolation seem to preserve the orthogonality. This suggests that I am overlooking something with the cubic interpolation.
Actually, if you consider $y_{h,i} = P y_{H,i}$ and take the scalar product $$ (y_{h,i}, y_{h,j})_{A_h} = y_{H,i}^T P^T A_h P y_{H,j} = 0 $$ then you see that fine-grid vectors will be always orthogonal if you take coarse-grid matrix $A_H$ to be related to fine-grid matrix by $$ A_H = P^T A_h P $$ which is quite a standard thing to be assumed in multigrid methods. So, it is more about "how to define properly coarse-grid operator after choosing the interpolation" than "how to define interpolation for given $A_h$ and $A_H$".