2
$\begingroup$

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 \begin{equation} \langle Rw, v \rangle_{k} \equiv \langle w, Iv \rangle_{k+1}, \end{equation} 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

$\endgroup$

2 Answers 2

2
$\begingroup$

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).

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – lightxbulb
    Oct 10, 2023 at 7:55
  • 1
    $\begingroup$ @lightxbulb I did not claim that one needs to construct the coarse matrix by restriction of the fine matrix. In geometric multigrid, it is typically assembled on the coarse mesh, like you say. But you still need that the multigrid correction operator is a symmetric operator. $\endgroup$ Oct 10, 2023 at 14:39
  • $\begingroup$ I misunderstood your comment about the correction matrix, since you had sandwiched the matrix between the restriction and interpolation. In the cited paper the requirement that the restriction should be the transpose of the interpolation is generalised to them being adjoints wrt the inner product with the mass matrix playing the role of metric tensor. Also note that here this is not used as a preconditioner for CG, instead a nonlinear version of CG is used as a "smoother" (irrespective of the fact that CG is really a rougher) in their multigrid approach. $\endgroup$
    – lightxbulb
    Oct 10, 2023 at 16:21
1
$\begingroup$

$$\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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.