Defining the restriction operator in nested multigrid FEM

I'm new in the multigrid approach and I'm trying to implement an algorithm from this paper: https://www.researchgate.net/publication/242913687_A_Multigrid_Algorithm_for_the_p-Laplacian

I'm stuck with the construction of the restriction operator $$R$$. I see people recommending taking $$R = I^T$$, where $$I$$ is my interpolation operator (I have no trouble constructing it). I'm not sure, however, if it works for this algorithm, since the authors define $$R$$ via $$$$\langle Rw, v \rangle_{k} \equiv \langle w, Iv \rangle_{k+1},$$$$ noting: I is the transpose of R with respect to inner products $$\langle \cdot \rangle_{k}$$ and $$\langle \cdot \rangle_{k+1}$$. Here $$\langle u,v \rangle_{k}$$ is defined as $$\langle u, v\rangle_k = \sum_{T \in {\cal T}^k} \int_{T} uv \: dT,$$ the sum being taken over the elements of the current finite element mesh $${\cal T}^k$$.

It is not obvious to me if $$I^T$$ satisfies this identity for $$R$$. It's known that in the fixed bases $$\{\phi^k\}$$ and $$\{\psi^{k+1}\}$$ the transpose of $$I$$ is represented by the matrix $$R = \Gamma_k^{-1} I^T \Gamma_{k+1},$$ where $$\Gamma_k$$ and $$\Gamma_{k+1}$$ are the Gram matrices of $$\langle \cdot \rangle_{k}$$ in $$\{\phi^k\}$$ and $$\langle \cdot \rangle_{k+1}$$ in $$\{\psi^{k+1}\}$$ respectively. I tried using this formula to get $$R$$, but the resulting matrix is in no way sparse which seems suspicious to me. Could someone tell me which is the proper way to choose $$R$$? Thank you

The operator $$R$$ is not the operator you would use to "restrict" a function from one level to another (which you would typically do via interpolation to the nodes of the coarser mesh). Rather, it is defined as the transpose of the interpolation operator, $$R := I^T$$ where you define $$I$$ in the obvious way as the interpolation onto a finer mesh.
You choose $$R$$ in this way because you want that for a symmetric matrix $$A$$, the correction matrix $$IA_\text{coarse}R$$ is also symmetric, because that's what you need for certain theoretical guarantees (and also so that you can use this matrix, or its inverse, as a preconditioner for CG).
• In the linked paper they do not construct the subproblems by using the interpolation and restriction operators, i.e. $A_{k-1} \ne I^{k-1}_k A_k I^k_{k-1}$. It's a geometric multigtid formulation with hierarchical meshes where they discretise the problem at each level, and the restriction and interpolation serve to just transfer the residual and correction between levels. The adjointness wrt the inner products on the meshes is to account for a potentially non-uniform mesh, otherwise it would be as if you set the mass matrix to the identity. Oct 10, 2023 at 7:55
$$\newcommand{\bm}[1]{\boldsymbol{#1}}$$ In the paper that you mentioned they define the coarse-to-fine/interpolation operator $$I^{k}_{k-1}:V^0_{k}\to V^0_{k-1}$$ as the natural injection: $$I^k_{k-1}v = v, \, v\in V^{0}_{k-1}.$$ Note that $$V^0_{k} \leq V^0_{k-1}$$ so there is no issue with that. For Lagrange elements this would just consist of setting the new fine grid node values to the interpolated values from the coarse grid mesh. That is, suppose you are given the coefficients vector $$\bm{v}^{k-1}\in\mathbb{R}^{n_{k-1}}$$ such that $$v = \sum_{j=1}^{n_{k-1}} v^{k-1}_j \phi^{k-1}_j$$ then you can find the coefficients vector $$\bm{v}^k\in\mathbb{R}^{n_k}$$ for $$v = \sum_{j=1}^{n_k} v^k_j \phi^k_j$$ through a matrix $$\bm{I}^{k}_{k-1}\in\mathbb{R}^{n_k \times n_{k-1}}$$: $$\bm{v}^{k} = \bm{I}^k_{k-1}\bm{v}^{k-1}$$, where the new node values are just the interpolated values along the edges. It's clear that $$\bm{I}^k_{k-1}$$ is sparse.
The fine-to-coarse/restriction operator $$I^{k-1}_k : V^0_{k-1} \to V^0_{k}$$ is chosen so that: $$\langle I^{k-1}_k w, v \rangle_{k-1} = \langle w, v\rangle_k = \langle w, I^k_{k-1}v \rangle_k, \, w\in V^0_k, \, v\in V^0_{k-1}.$$ If we were to expand both using the basis functions we get: \begin{align} \langle I^{k-1}_k w,v\rangle_{k-1} &= \left\langle \sum_{i=1}^{n_{k-1}}w^{k-1}_i\phi^{k-1}_i, \sum_{j=1}^{n_{k-1}}v^{k-1}_j \phi^{k-1}_j\right\rangle = (\bm{w}^{k-1})^T\bm{M}_{k-1} \bm{v}^{k-1} = (\bm{I}^{k-1}_k \bm{w}^k)^T\bm{M}_{k-1}\bm{v}^{k-1}, \\ \langle w,I^k_{k-1}v\rangle_{k} &= \left\langle \sum_{i=1}^{n_{k}}w^{k}_i\phi^{k}_i, \sum_{j=1}^{n_{k}}v^{k}_j \phi^{k}_j\right\rangle = (\bm{w}^{k})^T\bm{M}_k \bm{v}^{k} = (\bm{w}^k)^T\bm{M}_{k}(\bm{I}^k_{k-1}\bm{v}^{k-1}). \end{align} Equating the two yields: $$(\bm{w}^k)^T(\bm{I}^{k-1}_k)^T\bm{M}_{k-1}\bm{v}^{k-1} = (\bm{w}^k)^T\bm{M}_k\bm{I}_{k-1}^k\bm{v}^{k-1}.$$ So I believe that your derivation is correct and $$\bm{I}_{k}^{k-1} = \bm{M}_{k-1}^{-T}(\bm{I}_{k-1}^k)^T\bm{M}_{k}^T$$, and if $$\bm{M}_{k-1}$$ and $$\bm{M}_k$$ are symmetric we can drop the transposes for those. If you use mass matrices that are lumped (i.e. diagonal), then this is as sparse as $$(\bm{I}_{k-1}^k)^T$$ and only the rows and columns are rescaled. If you use non-lumped matrices then the matrix $$\bm{M}^{-1}_{k-1}$$ is in general non-sparse. However, there you can apply the matrix $$\bm{I}_{k}^{k-1}$$ as three (or two) consecutive matrices, in which case you can still exploit the sparsity at the expense of potentially some more memory reads and computations (you can solve the linear system $$\bm{M}_{k-1}\bm{y} = \bm{x}$$ in order to compute $$\bm{y} = \bm{M}^{-1}_{k-1}\bm{x}$$, which would still keep things sparse).