You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'We give the first mathematically rigorous analysis of an emerging approach to finite element analysis (see, e.g., Bauer et al. \[Appl. Numer. Math., 2017\]), which we hereby refer to as *the surrogate matrix methodology*. This methodology is based on the piece-wise smooth approximation of the matrices involved in a standard finite element discretization. In particular, it relies on the projection of smooth so-called stencil functions onto high-order polynomial subspaces. The performance advantage of the surrogate matrix methodology is seen in constructions where each stencil function uniquely determines the values of a significant collection of matrix entries. Such constructions are shown to be widely achievable through the use of locally-structured meshes. Therefore, this methodology can be applied to a wide variety of physically meaningful problems, including nonlinear problems and problems with curvilinear geometries. Rigorous *a priori* error analysis certifies the convergence of a novel surrogate method for the variable coefficient Poisson equation. The flexibility of the methodology is also demonstrated through the construction of novel methods for linear elasticity and nonlinear diffusion problems. In numerous numerical experiments, we demonstrate the efficacy of these new methods in a matrix-free environment with geometric multigrid solvers. In our experiments, up to a twenty-fold decrease in computation time is witnessed over the classical method with an otherwise identical implementation.'
author:
- 'Daniel Drzisga[^1]'
- Brendan Keith
- Barbara Wohlmuth
bibliography:
- 'main.bib'
title: 'The surrogate matrix methodology: [a priori]{} error estimation'
---

Surrogate numerical methods, finite element methods, matrix-free, high performance computing, a priori analysis, low order, geometric multigrid.

65D05, 65M60, 65N30, 65Y05, 65Y20.

Introduction {#sec:introduction}
============

In the field of computational science, major funding initiatives in North America, Europe, and Asia have thrust high performance computing (HPC) to ascendancy. In anticipation of future exascale computers, much work in this discipline involves the deep and careful reconstruction of long-established computing practices. An important characteristic of numerical algorithms aimed for these computers, is the floating-point operation (FLOP) per byte ratio. In order to achieve optimal performance and power efficiency on future machines, the time spent on FLOPs relative to memory transfer needs to be substantial.

Most traditional finite element softwares assemble global stiffness matrices by looping over elements and adding the corresponding local contributions to the global matrix. Storing the resulting sparse matrices requires significantly more memory than just storing the degrees of freedom. However, memory consumption is certainly not the only obstacle at the computational frontier. Indeed, at such scales, the memory traffic and latency involved in loading indices and entries for matrix vector products (MVPs) also presents critical challenges.

Since iterative solvers only require MVPs, it is not necessary to store all of the nonzeros of the global matrix in memory. Instead, it is sufficient to compute the nonzero entries on-the-fly, i.e., matrix-free. Different tactics exist to implement matrix-free methods, but the predominant candidate for low-order finite elements is the element-by-element approach [@Arbenz:2008:IJMME; @Bielak:2005:CMES; @Carey:1986:CANM; @Flaig:2012:LSSC; @Rietbergen:1996:IJMME], wherein local stiffness matrices are multiplied by local vectors and later added to the global solution vector. These local stiffness matrices may either be stored in memorywhich actually requires [more]{} memory than storing the global matrixor computed on-the-fly. When using high-order finite elements, the weak forms can be integrated on-the-fly using standard or reduced quadrature formulas [@Brown:2010:JSC; @Kronbichler:2012:CAF; @Ljungkvist:2017:HPC17; @Ljungkvist:2017:TechRep; @May:2015:CMAME]. This is a well-suited tactic for future machines because of its large arithmetic intensity [@loffeld2017arithmetic].

Significant performance gains are often attributed to a problem, scale, and architecture-specific balance between FLOPs and memory traffic. As a matter of course, exploiting symmetries in a problem or discretization can significantly improve the time to solution. This is often the cause of enormous speed-ups in computations on structured meshes. Likewise, high performance of matrix-free methods can be most easily achieved in homogeneous problems with simple geometries. Nevertheless, most geometries coming from significant real-world problems cannot be adequately approximated by fully-structured meshes. A possible trade-off is to use locally-structured meshes like hierarchical hybrid grids (HHG) where initially unstructured coarse grids are locally refined in a uniform way. This local structure allows the application of stencil-based finite element procedures which operate similar to finite difference methods. By using these grids, efficient stencil-based methods have been successfully applied to a wide range of problems [@bergen2005hierarchical; @bergen2004hierarchical; @bergen2007hierarchical; @Engwer:2017:CCPE; @gmeiner2015towards]. A related approach, suitable for low-order finite element discretizations of scalar elliptic partial differential equations (PDEs) with variable coefficients, based on scaling of reference stencils is discussed in [@bauer2017stencil].

In this paper, we revisit the classical lowest-order Bubnov–Galerkin finite element method and analyze a modification of it which is strongly amenable to stencil-based matrix-free computation. In our approach, a macro-mesh, which is not required to have any global structure, is used to triangulate the model geometry. This macro-mesh is then uniformly refined a large number of times, resulting in a fine-scale locally-structured mesh. For each macro-element, a *local approximation* of the fine-scale global matrix delivers a fine-scale *surrogate matrix* which maintains the convergence properties of the fine-scale discrete solution, up to the original order of the approximation. A related investigation [@bauer2017two] illustrated the promise of this methodology and provided numerical evidence for the convergence rates which are proven here rigorously. This work considered only Poisson’s equation. Later on, Stokes flow (with variable viscosity) was considered in two follow-up articles [@bauer2018large; @bauer2018new]. Each of these initial studies focused on the massively parallel high performance computing aspects of their respective methods. These studies used the HHG software framework [@bergen2005hierarchical; @bergen2004hierarchical; @bergen2007hierarchical] in their experiments. Here, the finite-element software framework HyTeG [@Kohl2018HyTeGfinite] is used.

In this paper, we recast the central features of the original work as a methodology complete with a mathematical framework suitable for rigorous analysis. The principal novelty is the mathematical foundation developed here, which can be used to analyze further incarnations of the methodology. In total, we consider three specific mathematical models; namely, the variable coefficient Poisson equation, linear elastostatics, and $p$-Laplacian diffusion. Although our presentation demonstrates that the surrogate matrix methodology applies to each of these models equally well, we only employ a complete *a priori* analysis of the simplest model, the variable coefficient Poisson equation. In our numerical experiments, we carefully verify the proven *a priori* convergence rates with the variable coefficient Poisson equation. We also include proof-of-concept demonstrations from numerical experiments with the linear elastostatics and $p$-Laplacian diffusion problems.

Notation and outline {#ssub:overview}
====================

Let $V$ be a reflexive Banach space over ${\mathbb{R}}$, the field of real numbers, and let $V_h\subsetneq V$ be a finite-dimensional subspace. Consider a continuous and weakly coercive bilinear form $a:V\times V \to {\mathbb{R}}$ and a bounded linear functional $F\in V^\ast$, the topological dual of $V$.

In this paper, we are concerned with the solutions $u$, $u_h$, and $\tilde{u}_h$ of the following three abstract variational problems.

\[eq:VariationalFormulations\] $$\begin{aligned}
{3}
    &
    \text{Find } u \in V \text{ satisfying }
    &
    \quad
    a(u,v)
    &=
    F(v)
        \quad
    &&\text{for all } v\in V
    \,.
\label{eq:ContinousVF}
    \\
    &
    \text{Find } u_h \in V_h \text{ satisfying }
    &
    \quad
    a(u_h,v_h)
    &=
    F(v_h)
        \quad
    &&\text{for all } v_h\in V_h
    \,.
\label{eq:DiscreteVF}
    \\
    &
    \text{Find } \tilde{u}_h \in V_h \text{ satisfying }
    &
    \quad
    \tilde{a}(\tilde{u}_h,v_h)
    &=
    F(v_h)
        \quad
    &&\text{for all } v_h\in V_h
    \,.
\label{eq:SurrogateVF}\end{aligned}$$

In \[eq:SurrogateVF\], a *surrogate* bilinear form $\tilde{a}:V_h\times V_h \to {\mathbb{R}}$ has been introduced. In order to properly define $\tilde{a}(\cdot,\cdot)$, some additional assumptions on $a(\cdot,\cdot)$ are still required; see .

The discrete variational problems \[eq:DiscreteVF,eq:SurrogateVF\] induce matrix equations for coefficients ${\mathsf{u}},\tilde{{\mathsf{u}}}$ in some ${\mathbb{R}}^N$, $${\mathsf{A}}{\mathsf{u}}= {\mathsf{f}}\qquad \text{ and } \qquad\tilde{{\mathsf{A}}}\tilde{{\mathsf{u}}} = {\mathsf{f}},
\label{eq:DiscreteVariationalProblems}$$ respectively. In the first case, fix a basis for $V_h$, say $\{\phi_i\}_{i=1}^N$. For this basis, each $(i,j)$-component of the stiffness matrix ${\mathsf{A}}$ is simply ${\mathsf{A}}_{ij} = a(\phi_j,\phi_i)$. In the following section, we present a methodology to construct a *surrogate* stiffness matrix $\tilde{{\mathsf{A}}}\approx {\mathsf{A}}$ which can be used in place of the *true* stiffness matrix ${\mathsf{A}}$. This methodology stands apart from technical details, such as differences in quadrature formulas. provides a short (non-comprehensive) list of examples fitting into our framework. In , we discuss the incorporation of non-homogeneous boundary conditions and what we hereon refer to as *the zero row sum property*. In , sufficient conditions for the discrete stability of surrogate bilinear forms $\tilde{a}(\cdot,\cdot)$ are briefly discussed. Next, in , we perform a rigorous *a priori* error analysis of our approach applied to the variable coefficient Poisson equation. A brief description of our implementation is given in . Then, in , we document several numerical experiments. Here, a thorough verification of each error estimate in is given. This is complemented by performance measurements for the additional examples.

Throughout this article, we assume that $\Omega\operatorname{\subseteq}{\mathbb{R}}^n$ is a bounded Lipschitz domain. For matrices ${\mathsf{M}}\in {\mathbb{R}}^{l\times m}$, define the $\ell^\infty$-, and ${\operatorname*{\vphantom{p}max}}$-norms, $\|{\mathsf{M}}\|_\infty = \operatorname*{\vphantom{p}max}_{i}\sum_{j}|{\mathsf{M}}_{ij}|$ and $\|{\mathsf{M}}\|_{\operatorname*{\vphantom{p}max}} = \operatorname*{\vphantom{p}max}_{i,j}|{\mathsf{M}}_{ij}|$. Likewise, for any function $v:\Omega\to {\mathbb{R}}$, we will use the similar notation, $\|v\|_0$, $\|v\|_1$, and $\|v\|_2$, for the canonical $L^2(\Omega)$-, $H^1(\Omega)$-, and $H^2(\Omega)$-norms, respectively. When dealing with a subset $T\operatorname{\subseteq}\Omega$, denote the related $L^2(T)$-, $H^1(T)$-, and $H^2(T)$-norms by $\|v\|_{0,T}$, $\|v\|_{1,T}$, and $\|v\|_{2,T}$, respectively. For any simplex $T$ and integer $0\leq q <\infty$, we denote the space of polynomials of degree at most $q$ as ${\mathcal{P}}_q(T)$. All remaining notation will be defined as it arises.

Surrogate stiffness matrices {#sec:surrogate_stiffness_matrix}
============================

In this section, we present the constitutive elements of the surrogate matrix methodology. Our approach here is to gradually introduce the necessary concepts, all the while maintaining a clear sense of generality. In order to arrive at a tractable framework for our problems of interest, we gradually refine the presentation from general $n$-dimensional spaces to only $n=1,2$ or $3$ and from general Banach spaces to only $W^{1,p}(\Omega)$ (or products thereof), where $1<p<\infty$. The intention of proceeding in this way is to indicate that the methodology can be applied to an extremely broad set of problems and, specifically, to most problems where finite element methods are traditionally applied.

Preliminary assumptions {#sub:assumptions_and_preliminaries}
-----------------------

Given a bounded domain $\Omega\operatorname{\subseteq}{\mathbb{R}}^n$, assume that the true bilinear form can be expressed as

$$a(u,v) = \int_\Omega G(x,u(x),v(x)) {\,\mathrm{d}}x
    \qquad \text{for all } u, v\in V.
\label{eq:BilinearFormStructure}$$

Additionally, upon defining $\operatorname*{supp}(u) = \{ y \in \Omega : u(y) \neq 0 \}$ for smooth functions, make the following sparsity assumption: $$G(x,u(y),v(y)) = 0
    \qquad
    \text{whenever }
            y\notin\operatorname*{supp}(u)\cap\operatorname*{supp}(v)
        .
\label{eq:SupportAssumption}$$

These assumptions permit us to consider the discretization of most classical differential operators. Indeed, in the assumptions above, the integrand $G(x,u,v)$ may induce distributional derivatives on its second and third arguments. Meanwhile, the first argument can be identified with the spatial argument of any associated variable coefficients. For example, in the weak form of a Poisson-type equation, $-\operatorname{div}(K\nabla u) = f$, with a variable, symmetric positive-definite diffusion tensor $K(x)$ (cf. ), we have the bilinear form $$a_1(u,v)
    =
    \int_\Omega \nabla u(x)^\top K(x){\hspace{0.25pt}}\nabla v(x) {\,\mathrm{d}}x
        \qquad \text{for all } u, v\in V=H^1(\Omega)
    \,.
\label{eq:BilinearFormPoisson}$$ Here, taking any point $x\in\Omega$, the integrand in \[eq:BilinearFormStructure\] reduces to $G(x,u,v) = {\nabla u^\top K(x){\hspace{0.25pt}}\nabla v}$. Evidently, this $G$ satisfies the sparsity assumption \[eq:SupportAssumption\].

Stencil functions {#sub:stencil_functions}
-----------------

Let $\phi\in V$ be a test function with compact support in $\Omega$ and, for any fixed $y\in {\mathbb{R}}^n$, define $\phi_y(x) = \phi(x-y)$. Now, consider any fixed set of ordered points ${\mathbb{X}}= \{x_i\}$ in $\Omega$ and recall \[eq:BilinearFormStructure\]. Assuming that both $\phi_{x_i},\phi_{x_j}\in V$, observe (via a simple change of variables) that $$\begin{aligned}
                        a(\phi_{x_j},\phi_{x_i})
                =
        \int_\Omega G(y,\phi_{x_j}(y),\phi_{x_i}(y)) {\,\mathrm{d}}y
                                                &=
        \int_{\Omega_{\delta}} \!\!G(x_i+y,\phi_{\delta}(y),\phi(y)) {\,\mathrm{d}}y
        \,,
    \end{aligned}
\label{eq:ChangeOfVars}$$ where $\delta = x_j-x_i$ and $\Omega_{\delta} = \operatorname*{supp}(\phi)\cap\operatorname*{supp}(\phi_{\delta})$. In the second equality, passing from an integral over $\Omega$ to an integral over the subset $\Omega_{\delta}\operatorname{\subseteq}\Omega$ follows immediately from the sparsity assumption \[eq:SupportAssumption\]. For each fixed $x_i$, the affine structure of the identity above may be illuminated by collecting each contributing translation $\delta$ into the set ${\mathscr{D}}(x_i) = \{x_j-x_i \,:\, x_j\in{\mathbb{X}},\, a(\phi_{x_j},\phi_{x_i}) \neq 0\}$ and defining a *stencil function* $$\Phi_i^\delta(x)
        =
    \int_{\Omega_\delta} G(x+y,\phi_{\delta}(y),\phi(y)) {\,\mathrm{d}}y
        \qquad
    \text{for each } \delta \in {\mathscr{D}}(x_i).
    \label{eq:StencilFunction}$$

We have just reduced the computation of any $a(\phi_{x_j},\phi_{x_i})$ to the evaluation of scalar-valued functions enumerated by affine coordinates $(x_i,x_j-x_i)$. Indeed, $$a(\phi_{x_j},\phi_{x_i})
    =
    \begin{cases}
    \Phi_i^\delta(x_i)
        \,,\quad & \text{if } \delta = x_j-x_i \in{\mathscr{D}}(x_i),\\
    0\,, & \text{otherwise.}
    \end{cases}
    \label{eq:StencilFunction2StiffnessMatrixGENERAL}$$ In the present scenario, there may be a different set of translations ${\mathscr{D}}(x_i)$ for every point $x_i$. However, if each point is drawn from a point lattice, most of the sets ${\mathscr{D}}(x_i)$ are identical. This observation is the subject of the following subsection and a foundational principle in our approach.

\[rem:Symmetry\] In the scenario that the bilinear form is symmetric, $a(u,v) = a(v,u)$, it is natural to assume that the stencil functions \[eq:StencilFunction\] will inherit a similar symmetry. Indeed, under the equivalent symmetry condition $G(x,u,v) = G(x,v,u)$ almost everywhere, one may easily verify that if $\delta = x_j - x_i$, then $$\Phi^\delta_i(x_i)
        =
        \int_{\Omega_{\delta}} \!\!G(x_i+y,\phi(y),\phi_{\delta}(y)) {\,\mathrm{d}}y
        =
                \int_\Omega G(x_j+y,\phi_{-\delta}(y),\phi(y)) {\,\mathrm{d}}y
        =
        \Phi^{-\delta}_j(x_j)
                                    \label{eq:SymmetricStencilFunction2StiffnessMatrixGENERAL}$$ or, equivalently, $\Phi^\delta_i(x_i) = \Phi^{-\delta}_j(x_i + \delta)$.

Local stencil functions and locally-structured meshes {#sub:stencil_functions_on_uniformly_structured_meshes}
-----------------------------------------------------

An affine point lattice ${\mathbb{L}}$, from here on referred to only as a *lattice*, is a regularly spaced array of points in ${\mathbb{R}}^n$ where every point $x_i\in{\mathbb{L}}$ belongs to a neighborhood containing no other points in ${\mathbb{L}}$. In this paper, each (possibly finite) lattice is determined by a finite linearly independent set of translations in ${\mathbb{R}}^n$; i.e., ${\mathbb{L}}\operatorname{\subseteq}\{ \delta_0 + a_1 \delta_1 + \cdots + a_l \delta_l \,:\, a_1,\ldots,a_l\in{\mathbb{Z}}\}$.

Assuming that the test function $\phi\in V$ is sufficiently localized and each point $x_i$ is drawn from a lattice ${\mathbb{L}}\operatorname{\subseteq}\Omega$, then each ${\mathscr{D}}(x_i)$ is a subset of a small number of admissible translations ${\mathscr{D}}({\mathbb{L}}) = \bigcup\{{\mathscr{D}}(x_i): x_i\in{\mathbb{L}}\}$, determined solely by the lattice structure. In such a scenario, every stencil function is closely related; i.e., $\Phi_i^\delta(x) = \Phi_j^\delta(x)$, whenever both are defined. Therefore, it is prudent to drop the subscript and define only one common stencil function $\Phi^\delta(x)$ for each $\delta \in {\mathscr{D}}({\mathbb{L}})$. Clearly, $$a(\phi_{x_j},\phi_{x_i})
    =
    \begin{cases}
    \Phi^\delta(x_i)\,,\quad & \text{if } \delta = x_j-x_i \in{\mathscr{D}}({\mathbb{L}}),\\
    0\,, & \text{otherwise.}
    \end{cases}
    \label{eq:StencilFunction2StiffnessMatrixSTRUCTURED}$$

We are interested in exploiting \[eq:StencilFunction2StiffnessMatrixSTRUCTURED\] for solving a wide variety of PDEs with curvilinear geometries. Toward this end, the following examples help motivate our construction further.

\[rem:ReflexiveSobolevSpaces\] From now on, it is useful to let $V$ be a closed subset of $W^{1,p}(\Omega)$, for some $1<p<\infty$. This will allow us to use a basis for $V_h \operatorname{\subseteq}V$ consisting of finite element vertex functions [@ern2013theory]. Other problems where $V\operatorname{\subseteq}\big[W^{1,p}(\Omega)\big]^k$, $k\in{\mathbb{N}}$, can be handled similarly (see, e.g., ), by employing a basis consisting of the same vertex functions in each component.

### The one-dimensional setting {#ssub:the_one_dimensional_setting}

Let $V=H^1_0(\Omega)$, where $\Omega=(0,1)\operatorname{\subseteq}{\mathbb{R}}$, and fix a small translation $dx = 1/(N+1)$. Consider the scenario where each point, ${x_i = x_{i-1} + dx}$, evenly divides $\Omega$ and $\phi$ is the piecewise-linear hat function defined $\phi(x) = \operatorname*{\vphantom{p}max}(1-\frac{|x|}{dx},0)$. Let $V_h = \{v \in H^1_0(\Omega) : v|_t \in {\mathcal{P}}_1(t),\,t = (x_i, x_{i+1}), \text{ for each } 1\leq i\leq N \}$ and identify each shifted hat function with the standard basis, $\phi_{x_i} = \phi_i\in V_h$. Here, we may define ${\mathbb{L}}= \{x_i\}$. In this case, for each $i\geq 1$, the value $a(\phi_i,\phi_j)$ can either be computed directly from , in the standard way, or evaluated using \[eq:StencilFunction2StiffnessMatrixSTRUCTURED\], assuming that each $\Phi^\delta$ is available at the onset of computation. Note that ${\mathscr{D}}(x_1) = \{0,dx\}$ and ${\mathscr{D}}(x_N) = \{-dx,0\}$, but ${\mathscr{D}}(x_i) = \{-dx,0,dx\}$ for each $2\leq i \leq N-1$. Therefore, ${\mathscr{D}}({\mathbb{L}}) = \{-dx,0,dx\}$. Ultimately, defining each structured stencil function $\Phi^\delta(x) = \Phi((x_i,\delta);x)$ from an arbitrary candidate point $x_i$, one may verify that $$a(\phi_j,\phi_i)
    =
    \begin{cases}
    \Phi^\delta(x_i)\,,\quad & \text{if } \delta = x_j-x_i \in \{dx,0,dx\},\\
    0\,, & \text{otherwise,}
    \end{cases}$$ which is clearly the same format as \[eq:StencilFunction2StiffnessMatrixSTRUCTURED\].

Recall \[rem:Symmetry\]. If $a(\cdot,\cdot)$ is symmetric, then, by \[eq:SymmetricStencilFunction2StiffnessMatrixGENERAL\], $\Phi^{dx}(x_i) = \Phi^{-dx}(x_i+dx)$ and the number of required stencil functions can be reduced to two. In some situationse.g., when the zero row sum property can be employed (see )only a single stencil function is actually required.

### Locally-structured meshes with triangles {#ssub:hierarchical_hybrid_grids_with_triangles}

Let $m\in{\mathbb{N}}$. Beginning with a scaled Cartesian lattice ${\mathbb{L}}_m = 2^{-m}{\mathbb{Z}}^n$, it is useful to define its intersection with the closure of the right-angled reference simplex $\hat{T} = \{\hat{x}\in{\mathbb{R}}^n : \|\hat{x}\|_1 < 1,\, \hat{x}\cdot e_i > 0,\,\forall i=1,\ldots,n\}$. This simplicial lattice, $\hat{T}_m = {\mathbb{L}}_m \cap \overline{\hat{T}}$, can easily be transformed into a similar simplicial lattice $T_m$ for any arbitrary simplex $T\operatorname{\subseteq}\Omega$ via an affine transformation (see, e.g., ). Indeed, fixing the unique $A\in{\mathbb{R}}^{n\times n}$ and $b\in{\mathbb{R}}^{n}$ such that $T = \{A\hat{x}+b : \hat{x}\in\hat{T}\}$, the corresponding lattice is clearly $T_m = \{ A\hat{x}_i+b : \hat{x}_i\in \hat{T}_m\}$. Note that there are no interior points, $\mathring{T}_m = T_m\cap T = \emptyset$, if $m<2$. We also define the set of boundary points ${\partial}T_m = T_m \setminus \mathring{T}_m$, which is always non-empty.

Let $\Omega\operatorname{\subseteq}{\mathbb{R}}^n$, where $n=2,3$. The utility of this transformation is evident upon considering a *macro*-triangulation of $\Omega$, say ${{\mathcal{T}}}_H$, where each *macro*-element $T\in{{\mathcal{T}}}_H$ is endowed with a simplicial lattice $T_m$, as defined above. Here, $H = \operatorname*{\vphantom{p}max}_{T\in{{\mathcal{T}}}_H} H_T$, where each $H_T = \mathrm{diam}(T)$ denotes the diameter of $T$. Notice that, for any fixed level $m\geq0$, all interface points $x_i\in {\partial}T_m\cap\Omega$ are coincident with an interface lattice point on some neighbouring simplex. We may now define the set of all vertices on level $m$: $${\mathbb{X}}_m = \bigcup \{T_m:T\in{{\mathcal{T}}}_H\}.$$

For every ${\mathbb{X}}_m$, there is a corresponding finite element mesh such that each vertex function within a fixed macro-element $T\in{{\mathcal{T}}}_H$ is self-similar. In the cases $n=2$ or $3$, we are left to define each triangle or tetrahedron whose vertices coincide with points in ${\mathbb{X}}_m$. Let $[y_0,y_1,\ldots,y_k]\operatorname{\subseteq}{\mathbb{R}}^n$ denote the convex combination of the points $y_0,y_1,\ldots,y_k\in\Omega$. When $n=2$, the natural construction begins by considering the following uniform subdivision of a triangle $T = [y_0,y_1,y_2] \in{{\mathcal{T}}}_H$ into a set of four equal-volume triangles: $${\mathcal{S}}(T)
    =
    \big\{[y_0,\tfrac{y_0+y_1}{2},\tfrac{y_0+y_2}{2}],[\tfrac{y_0+y_1}{2},y_1,\tfrac{y_1+y_2}{2}],[\tfrac{y_0+y_2}{2},\tfrac{y_1+y_2}{2},y_2],[\tfrac{y_0+y_1}{2},\tfrac{y_1+y_2}{2},\tfrac{y_0+y_2}{2}]\big\}
    .$$ For an illustration of this $n=2$ case, see \[fig:Subdivision\]. For $n=3$, see the construction in [@bey1995tetrahedral]. Further subdivisions can then be defined recursively, *viz.*, $${\mathcal{S}}^{m+1}(T)
    = 
    \bigcup\{{\mathcal{S}}(t) : t \in {\mathcal{S}}^{m}(T)\}
        \qquad \text{ for all } m\geq 1\,.$$

![\[fig:Subdivision\]Illustration of three refinement steps of a single macro-element $T$ for $n=2$ with the corresponding vertex lattices $T_0$, $T_1$, $T_2$, and $T_3$. The interior lattice points $\mathring{T}_m$ are colored blue and the boundary lattice points $\partial T_m = T_m \setminus \mathring{T}_m$ are colored in orange. Additionally, the support of an exemplary vertex function $\phi_i$ is shaded in green.](triangle_refinement_level_1.pdf){width="\textwidth"}

The set of all vertices in a given subdivision ${\mathcal{S}}^{m}(T)$ forms an evenly spaced set of points inside $T$. The set of vertices in ${\mathcal{S}}^{1}(T)$ clearly coincides with $T_1$ and one can easily verify from the recursive definition that the set of vertices in ${\mathcal{S}}^{m}(T)$ also coincides with $T_m$. We may finally define the sequence of *locally-structured meshes*: $${\mathcal{S}}^m({{\mathcal{T}}}_H)
    =
    \bigcup \{{\mathcal{S}}^m(T) : T\in{{\mathcal{T}}}_H\},
    \qquad \text{ for all } m\geq 1\,.
\label{eq:FineScaleMesh}$$ Notice that for each $m$, ${\mathbb{X}}_m$ coincides with the set of all vertices in ${\mathcal{S}}^m({{\mathcal{T}}}_H)$. Therefore, each point $x_i$ in ${\mathbb{X}}_m$ can be identified with a vertex function $\phi_i$ supported by only the neighboring elements appearing in ${\mathcal{S}}^m({{\mathcal{T}}}_H)$; see \[fig:Subdivision\].

Assume that the fine mesh level, $m$, is chosen large enough that one may find several points in ${\mathbb{X}}_m$ which lie in the interior of some macro-element $T$. If $x_i\in \mathring{T}_m$ is any such point, then the translation set ${\mathscr{D}}(x_i)={\mathscr{D}}(\mathring{T}_m)$ will only contain translations aligned with edges appearing in the original subdivision ${\mathcal{S}}^1(T)$. Identifying $\phi_{x_i}=\phi_i$ we see that for every $x_i\in \mathring{T}_m$ and $x_j\in T_m$, $$a(\phi_{j},\phi_{i})
    =
    \begin{cases}
    \Phi^\delta_{T}(x_i)\,,\quad & \text{if } \delta = (x_j-x_i) \in{\mathscr{D}}(\mathring{T}_m),\\
        0\,, & \text{otherwise,}
    \end{cases}
    \label{eq:StencilFunction2StiffnessMatrixHHG}$$ where $\Phi^\delta_{T}:\mathrm{conv}(\mathring{T}_m)\to{\mathbb{R}}$ is a local stencil function for the current level $m$ and macro-element $T$ and $\mathrm{conv}(\mathring{T}_m)$ is the convex hull of $\mathring{T}_m$. In general, notice that $\Phi^\delta_{T} \neq \Phi^{2\delta}_{T}$ are two different stencil functions, corresponding to the same direction but different mesh levels.

Consider the bilinear form \[eq:BilinearFormPoisson\] with a variable diffusion coefficient. A visualization of several corresponding local stencil functions, coming from locally-structured meshes used in our numerical experiments, is given in . It is clear from this figure that each $\Phi^\delta_{T}$ has the potential to be a smooth function. We now come to the final essential component of our surrogate methodology; the approximation of $\Phi^\delta_{T}$.

![\[fig:stencilweight\_surfaceplot\]Left: Surface plots of the local stencil functions $\Phi^\delta_{T}(x)$, for the degenerate direction $\delta=0$ and level $m = 5$, from the numerical experiments recounted in . Here, the function is plotted over each subset $\mathrm{conv}(\mathring{T}_m)\operatorname{\subseteq}T\in{{\mathcal{T}}}_H$. In this case, it clearly appears that each stencil function can be related to the restriction of a globally continuous function $\Phi^\delta_{T}(x)$. This is a result of the structure of the macro-mesh only. Right: Surface plots of the stencil functions $\Phi^\delta_{T}(x)$ after the first time step from the experiment recounted in , for the eastern direction $\delta$ relative to each macro-element and level $m = 5$. Moreover, although they are clearly related, it is evident that the corresponding stencil functions lack any global smoothness property.](stencilweight_surfaceplot.png){width="\linewidth"}

![\[fig:stencilweight\_surfaceplot\]Left: Surface plots of the local stencil functions $\Phi^\delta_{T}(x)$, for the degenerate direction $\delta=0$ and level $m = 5$, from the numerical experiments recounted in . Here, the function is plotted over each subset $\mathrm{conv}(\mathring{T}_m)\operatorname{\subseteq}T\in{{\mathcal{T}}}_H$. In this case, it clearly appears that each stencil function can be related to the restriction of a globally continuous function $\Phi^\delta_{T}(x)$. This is a result of the structure of the macro-mesh only. Right: Surface plots of the stencil functions $\Phi^\delta_{T}(x)$ after the first time step from the experiment recounted in , for the eastern direction $\delta$ relative to each macro-element and level $m = 5$. Moreover, although they are clearly related, it is evident that the corresponding stencil functions lack any global smoothness property.](stencilweight_surfaceplot_unstructured.png){width="\linewidth"}

\[rem:Geometry\] The lattice structure of locally-structured meshes is destroyed under smooth, non-affine transformations $\hat{T} \to \tilde{T}\operatorname{\subseteq}\Omega$. This offers a possible impediment to our construction in the case of non-polygonal domains $\overline{\Omega} \neq \overline{\Omega}_H = \bigcup_{T\in{{\mathcal{T}}}_H} \overline{T}$. In fact, if a globally continuous transformation $\varphi:\overline{\Omega}_H \to \overline{\Omega}$ is available, with $\varphi|_T$ a smooth bijection for every $T\in{{\mathcal{T}}}_H$, then an equivalent method can be found using local pull-backs of $\varphi$. For example, the bilinear form in \[eq:BilinearFormPoisson\], $a_1:H^1(\Omega)\times H^1(\Omega) \to {\mathbb{R}}$, simply transforms to $a_{1,H}:H^1(\Omega_H)\times H^1(\Omega_H) \to {\mathbb{R}}$, where $$a_{1,H}(u,v)
        =
        \int_{\Omega_H} \nabla u(x)^\top K_{H}(x) \nabla v(x) {\,\mathrm{d}}x
        \,,
                        \qquad
        K_{H} = \frac{D\varphi^{-1}{\hspace{0.25pt}}(K\circ\varphi){\hspace{0.25pt}}D\varphi^{-\top}}{|\det{\left(D\varphi^{-1}\right)}|}
        \,.
    \label{eq:BilinearFormPoissonTensor}$$ Therefore, from now on, we operate under the assumption that the macro-mesh ${{\mathcal{T}}}_H$ is geometrically conforming, $\overline{\Omega} = \overline{\Omega}_H$.

The inherited regularity of stencil functions {#sub:the_inherited_regularity_of_stencil_functions}
---------------------------------------------

Our approach is to locally project each stencil function $\Phi^\delta_{T}$ in \[eq:StencilFunction2StiffnessMatrixHHG\] onto a high-dimensional space of polynomials and later use this projection to compute approximate values of the stiffness matrix ${\mathsf{A}}$. Let ${{\mathcal{T}}}_H$ be a shape-regular simplicial mesh and consider a locally-structured mesh of level $m$ subordinate to this mesh, ${\mathcal{S}}^m({{\mathcal{T}}}_H)$. Before moving on, observe that definition \[eq:StencilFunction2StiffnessMatrixHHG\] in fact holds for any $x_i,x_j\in T_m$ if $x_i$ or $x_j \in \mathring{T}_m$. Therefore, due to the structure of the vertex functions in a locally-structured mesh, the domain of each $\Phi^\delta_{T}$ can actually be extended to a set $\overline{T}_\delta$ lying between $\mathrm{conv}(\mathring{T}_m)$ and $\overline{T}$. Indeed, identifying the test function in \[eq:StencilFunction\] with the $i$-th vertex function, $\phi(x) = \phi_i(x+x_i)$, for an arbitrary vertex $x_i\in\mathring{T}_m$, define $$T_\delta
    =
        \{x\in T : x + y \in T, \text{ for all } y\in\Omega_\delta = \operatorname*{supp}(\phi)\cap\operatorname*{supp}(\phi_\delta)\}
    .$$ From now on, we assume $\Phi^\delta_{T}:\overline{T}_\delta\to{\mathbb{R}}$. See \[fig:DomainsofTheStencilFunctions\] for a depiction of the sets in a triangular mesh and note that $T_{-\delta} = T_\delta + \delta$, for every $\delta \in {\mathscr{D}}(T_m)$.

![\[fig:DomainsofTheStencilFunctions\]Illustrations of the domains $T_\delta$ in gray and $\Omega_{\delta}$ for six exemplary directions $\delta$. Left: Northern and north-western direction. Middle: Eastern and south-eastern direction. Right: Western and southern direction.](triangle_direction_domains.pdf){width="80.00000%"}

For any $T\in{{\mathcal{T}}}_H$, let ${\mathcal{P}}_q(T_\delta)$ denote the space of polynomials of degree at most $q$ on the simplex $\overline{T}_\delta$ and let $\Pi^\delta_T:C^0(\overline{T}_\delta)\to {\mathcal{P}}_q(T_\delta)$ be an $L^\infty$-continuous projection operator, $\Pi^\delta_T\circ \Pi^\delta_T = \Pi^\delta_T$. For each macro-element $T\in{{\mathcal{T}}}_H$ and level $m\in{\mathbb{N}}$, define the *surrogate stencil function* ${\tilde{\Phi}^\delta_{T}\!:{T}_\delta\to{\mathbb{R}}}$ to be the corresponding polynomial projection of $\Phi^\delta_{T}$. Namely,

\^\_[T]{} = \^\_T \^\_[T]{} . \[eq:StencilInterpolation\]

In order to correctly argue that a polynomial approximation of $\Phi^\delta_{T}$ is feasible, it is necessary to classify its regularity depending on the problem at hand. In the following proposition, we show, under the modest assumptions above, that if $G(\cdot|_T,\cdot,\cdot)$ is a polynomial in its first argument, then $\Phi^\delta_{T}$ is also a polynomial of the same degree.

\[prop:Interpolation\] Fix a simplex $T\in{{\mathcal{T}}}_H$. Assume that the bilinear form $a(\cdot,\cdot)$ in \[eq:ContinousVF\] satisfies assumptions \[eq:BilinearFormStructure,eq:SupportAssumption\], where the integrand $G(\cdot|_T,\cdot,\cdot)$ is a polynomial of at most degree $q$ in its first argument. Then, for any locally-structured mesh ${\mathcal{S}}^m({{\mathcal{T}}}_H)$, as defined above, every local stencil function $\Phi^\delta_{T} \in {\mathcal{P}}_q(T_\delta)$ is a polynomial of the same degree.

Recall definition \[eq:StencilFunction\]. For every $x_i,x_j\in{\mathbb{X}}_m$, every stencil function $\Phi((x_i,x_j-x_i),x)$ is defined with a test function $\phi\in V$. Fixing any arbitrary vertex $x_i\in\mathring{T}_m$, identify this test function with the $i$-th vertex function, $\phi(x) = \phi_i(x+x_i)$. Let $\alpha = (\alpha_1,\ldots,\alpha_n) \in {\mathbb{N}}_0$ denote a standard multi-index, $|\alpha| = \sum_{i=1}^n |\alpha_i|$ and $x^\alpha = x_{\scriptscriptstyle 1}^{\alpha_1}\cdots x_{\scriptscriptstyle n}^{\alpha_n}$. By assumption, we may express $G(x, \phi_\delta(y),\phi(y)) = \sum_{|\alpha|\leq l} c_\alpha(y) x^\alpha$, where each coefficient function $c(y)$ has support only in $\Omega_\delta = \operatorname*{supp}(\phi)\cap\operatorname*{supp}(\phi_\delta)$. Moreover, if $x+y\in T$, then $$G(x+y, \phi_\delta(y),\phi(y))
        =
        \sum_{|\alpha|\leq l} c_\alpha(y) (x+y)^\alpha
        =
        \sum_{|\alpha|\leq l}\sum_{|\nu|\leq \alpha}\binom{\alpha}{\nu} c_\alpha(y)y^\nu x^{\alpha-\nu}
            \label{eq:PolyExpansion}$$ is clearly an equal degree polynomial in the variable $x$. The proof is completed by noting that the integral in \[eq:StencilFunction\] is performed only over the variable $y\in\Omega_\delta$ and so the stencil function acts like a convolution. Indeed, under the assumption $x\in T$, only the subset of points $x\in T_\delta = \{x\in T : x+y \in T \text{ for all } y\in \Omega_\delta\}$ guarantee that \[eq:PolyExpansion\] holds at every point of integration $y$. In this case, by linearity of integration, $\Phi^\delta_{T}$ is a member of ${\mathcal{P}}_q(T_\delta)$, with its coefficients defined by the associated integrals of the $y$-dependent functions in the right-hand side of \[eq:PolyExpansion\].

\[cor:Interpolation\] In the setting of , $
        \tilde{\Phi}^\delta_{T}
        =
        \Phi^\delta_{T}
    $.

Since $\Phi^\delta_{T} \in {\mathcal{P}}_q(T_\delta)$, we immediately see that $\Pi^\delta_T\, \Phi^\delta_{T} = \Phi^\delta_{T}$.

Surrogate stiffness matrices {#sub:the_surrogate_stiffness_matrix}
----------------------------

The main goal of the entire effort above is to guide us in reducing the vast majority of the finite element assembly process to the evaluation of a small set of functions which can, in fact, be locally approximated by polynomials. As in , let ${{\mathcal{T}}}_H$ be a shape-regular triangulation of a domain $\Omega$ into disjoint *macro*-simplices $T$.

Our construction of the surrogate matrix $\tilde{{\mathsf{A}}}$ is obviously built to exploit the local lattice structure of locally-structured meshes. As argued previously, a surrogate stencil function $\tilde{\Phi}^\delta_{T}$ can be used to approximate any matrix entry ${\mathsf{A}}_{ij}$ coming from a locally-structured mesh if at least one of the corresponding vertices $x_i$ or $x_j$ belongs to $\mathring{T}_m$. This leaves us to define only the nonzero matrix entries coming from the mutual interaction of vertex functions at the boundaries of the macro-elements. Although these entries could also be approximated by surrogate stencil functions in this case, these additional functions would be defined on each subsimplex of the macro-mesh ${{\mathcal{T}}}_H$ let us assume that they are computed directly. Because the growth of the macro-mesh boundary and interface interactions grow at an order of magnitude less than the interior interactions, computing these matrix entries directly does not affect to the asymptotic performance of the methodology. Finally, letting ${\partial}{\mathbb{X}}_m= \bigcup_{T\in{{\mathcal{T}}}_H} {\partial}T_m$ denote the union of all macro-mesh boundary vertices, we define the general surrogate stiffness matrix $$\tilde{{\mathsf{A}}}_{ij}
    =
    \begin{cases}
    \int_\Omega G(y,\phi_j(y),\phi_i(y)) {\,\mathrm{d}}y\,,\quad & \text{if both } x_i \text{ and } x_j \in {\partial}{\mathbb{X}}_m \\
    \tilde{\Phi}^\delta_{T}(x_i)\,,\quad & \text{if } \delta = (x_j-x_i) \in{\mathscr{D}}(\mathring{T}_m) \text{ and } x_i \text{ or } x_j \in \mathring{T}_m,\\
    0\,, & \text{otherwise.}
    \end{cases}
    \label{eq:SurrogateStiffnessMatrix}$$

\[rem:SurrogateSymmetry\] Due to the presence of the surrogate stencil functions $\tilde{\Phi}^\delta_T$, even if ${\mathsf{A}}$ is symmetric, $\tilde{{\mathsf{A}}}$ will generally not be. However, recalling \[eq:SymmetricStencilFunction2StiffnessMatrixGENERAL\], observe that if $a(\cdot,\cdot)$ is symmetric, then $\Phi^\delta_{T}(x) = \Phi^{-\delta}_{T}(x+\delta)$. Therefore, if we use related projection operators $$\Big[\Pi^\delta_T \Phi^\delta_{T}\Big](x) = \Big[\Pi^{-\delta}_T \Phi^{-\delta}_{T}\Big](x+\delta)\,,
\label{eq:SymmetricInterpolants}$$ for each opposing direction $\delta$ and $-\delta$, then $\tilde{{\mathsf{A}}}$ will be symmetric. Indeed, if $\delta = x_j-x_i$, then $$\tilde{{\mathsf{A}}}_{ij}
    =
    \Pi^\delta_T \Phi^\delta_{T}(x_i)
    =
    \Pi^{-\delta}_T \Phi^{-\delta}_{T}(x_i+\delta)
    =
    \Pi^{-\delta}_T \Phi^{-\delta}_{T}(x_j)
    =
    \tilde{{\mathsf{A}}}_{ji}
    \,.$$

Examples {#sec:examples}
========

In this section, we present three example problems which easily fit into the framework above.

The variable coefficient Poisson equation {#sub:poisson_s_equation}
-----------------------------------------

Consider the Poisson-type equation ${-\mathrm{div}(K\nabla u) = f}$ in $\Omega$, $u = 0$ on ${\partial}\Omega$, with a load $f\in L^2(\Omega)$ and a variable, symmetric positive-definite tensor $K$. Furthermore, assume that for each index $a,b$, $K_{ab}\in {\mathcal{P}}_q(\Omega)$. Recall \[eq:BilinearFormPoisson\] and note that we have already shown that the weak form of this problem can be cast into the framework above.

Assume that, within some set $T\operatorname{\subseteq}\Omega$, each vertex function $\phi_i$ is a translation of a fixed test function $\phi(x) = \phi_i(x-x_i)$. Then for each $\phi_i,\phi_j$, the stiffness matrix entry $${\mathsf{A}}_{ij}
        =
    \int_{\Omega} \nabla \phi_i(x)^\top K(x){\hspace{0.25pt}}\nabla \phi_j(x)\hspace{0.5pt} {\,\mathrm{d}}x
        \label{eq:TrueStiffnessMatrixPoisson}$$ can equally well be expressed as the evaluation (at the point $x_i$) of a stencil function, which, by , is simply a polynomial of the same degree as the diffusion tensor $K$. In the case of locally-structured meshes, there is a locally defined stencil function $\Phi^{\delta}_{T}:\Omega\to{\mathbb{R}}$ for each macro-element $T$ and level $m$. In this case, each $\Phi^{\delta}_{T}:\overline{T}_\delta\to{\mathbb{R}}$ is a polynomial (of degree at most $q$) on $T_\delta$.

Linearized elasticity {#sub:linearized_elasticity}
---------------------

Let $\vec{\nabla}$ and $\mathrm{Div}$ denote the row-wise distributional gradient and divergence, respectively. Now define $\epsilon(u) = \frac{1}{2}[\vec{\nabla} u + (\vec{\nabla}u)^\top]$ to be the symmetric gradient operator $\epsilon:H^1(\Omega)^n \to L^2(\Omega)^n$, where $n\geq 2$. Consider the following standard PDE model for the displacement $u\in H^1_0(\Omega)^n$ of a linearly elastic isotropic material: $-\mathrm{Div}\, \sigma = \vec{f}$, where the stress $\sigma = 2\mu \epsilon(u) + \lambda I \mathrm{div}\,u$ and the load $f\in L^2(\Omega)^n$.

The weak form of this equation is well known in the literature [@ciarlet1988three] and the associated bilinear form is simply $$a_2(u,v) 
    =
    \int_{\Omega} 2\mu \epsilon(u):\epsilon(v) + \lambda\hspace{0.75pt} \mathrm{div}(u)\, \mathrm{div}(v)\hspace{0.5pt} {\,\mathrm{d}}x
    \qquad \text{for all } u, v\in \big[H^1_0(\Omega)\big]^n
    .
\label{eq:BilinearFormLinearElasticity}$$ This bilinear form obviously satisfies assumptions \[eq:BilinearFormStructure,eq:SupportAssumption\]. If we assume that the Lamé parameters $\mu,\lambda:\Omega\to{\mathbb{R}}$ are piecewise polynomials on a collection of disjoint subdomains $T\in{\mathcal{T}}_H$, then each associated stencil function is also a piecewise polynomial.

$p$-Laplacian diffusion {#sub:p-Laplacian_diffusion}
-----------------------

For any $1<p<\infty$, let $\Delta_p\, u = \operatorname{div}(|\nabla u|^{p-2}\nabla u)$ be the $p$-Laplacian operator. Fix a valid parameter $p$ and consider the nonlinear diffusion equation $\frac{\partial u}{\partial t} - \Delta_p u = f$, where $f\in L^p(\Omega)$. A simple Euler time-stepping scheme replaces the time derivative $\frac{\partial u}{\partial t}$ by the quotient $\frac{u_{k+1}-u_k}{d t}$, where $dt>0$ is a fixed time-step parameter. Choosing backward Euler time-stepping and defining $f_k = f(k\cdot dt)$, we arrive at a semi-discrete nonlinear elliptic PDE for the solution variable $u_k$, which must be solved at each step $k\in{\mathbb{N}}$: $u_{k} - dt\hspace{0.5pt}  \Delta_p u_{k} = dt\hspace{0.5pt} f_k + u_{k-1}$. Upon fixed point linearization of the weak form of this equation, we uncover the following bilinear form: $$b(u,v)
    =
    \int_{\Omega} dt\, |\nabla \tilde{u}|^{p-2}\, \nabla u\cdot\nabla v + u \, v\hspace{0.5pt} {\,\mathrm{d}}x
    ,
    \qquad
    \text{for all } u,v\in W^{1,p}(\Omega)
    \,.
\label{eq:BilinearFormpLaplace}$$ Here, the variable coefficient $\tilde{u}\in W^{1,p}(\Omega)$ is usually identified with the previous solution iteration in the associated fixed point algorithm (cf. ).

The bilinear form $b(\cdot,\cdot)$ can easily be placed into the form of \[eq:BilinearFormStructure\] and each matrix entry can therefore be superceded by stencil function evaluations, $\Phi^\delta_T(x)$, as in \[eq:StencilFunction\]. Alternatively, one may split $b(\cdot,\cdot)$ into a mass term $m(\cdot,\cdot)$ and a $dt$-weighted stiffness term $a_3(\cdot,\cdot)$. Specifically, $b(u,v) = m(u,v) + dt\cdot a_3(u,v)$, where $$m(u,v)
    =
    \int_{\Omega} u \, v\hspace{0.5pt} {\,\mathrm{d}}x
    \qquad
    \text{and}
    \qquad
    a_3(u,v)
    =
    \int_{\Omega} |\nabla \tilde{u}|^{p-2}\, \nabla u\cdot\nabla v\hspace{0.5pt} {\,\mathrm{d}}x
    \,.
\label{eq:SplitBilinearFormspLaplace}$$ With this observation in hand, we see that $b(u,v)$ may be discretized by a linear combination of independent surrogates; one for $m(\cdot,\cdot)$ and one for $a_3(\cdot,\cdot)$ (cf. ). In either approach, the variable coefficient $|\nabla \tilde{u}(x)|^{p-2}$ will generally not remain a polynomial in a subdomain of $\Omega$ and the accuracy of a surrogate stencil function $\tilde{\Phi}^\delta_T$ will reflect the local regularity of the solution from the previous iteration, $\tilde{u}$.

Boundary conditions and the zero row sum property {#sec:the_zero_row_sum_property}
=================================================

It is generally appropriate to define the surrogate stiffness matrix component-wise by the rule given in \[eq:SurrogateStiffnessMatrix\]. Nevertheless, in some problems the operator to be discretized has a kernel which is not guaranteed to be respected by the surrogate. In such scenarios, it is possible that better performance and accuracy can be achieved if elements of this kernel are incorporated into the construction of the surrogate. This occurrence is most easily illustrated with the Poisson example from .

Consider the bilinear form $a_1:H^1(\Omega)\times H^1(\Omega)\to {\mathbb{R}}$, defined in \[eq:BilinearFormPoisson\]. Define $V_h^{{\mathrm{ext}}}= \{ v \in H^1(\Omega) : v|_t \in {\mathcal{P}}_1(t) \text{ for each } t\in{\mathcal{S}}^m({{\mathcal{T}}}_H)\}$ and $V_{h} = \{ v \in H^1_0(\Omega) : v|_t \in {\mathcal{P}}_1(t) \text{ for each } t\in{\mathcal{S}}^m({{\mathcal{T}}}_H)\} \operatorname{\subseteq}V_h^{{\mathrm{ext}}}$. Let the corresponding vertex function bases be $\{\phi_i\} \operatorname{\subseteq}\{\phi_i^{{\mathrm{ext}}}\}$, with $\phi_i = \phi^{{\mathrm{ext}}}_i$ for $1\leq i\leq N$. In most finite element software, a space like $V_h^{{\mathrm{ext}}}$ is used to impose Dirichlet boundary conditions. Indeed, a “lift” of the Dirichlet data, say $u_h^{{\mathrm{ext}}}= \sum_i {\mathsf{u}}^{{\mathrm{ext}}}_i\phi_i^{{\mathrm{ext}}}$, is generally constructed from a linear combination of the set $\{\phi_i^{{\mathrm{ext}}}\}\setminus\{\phi_i\}$. Then, taking ${\mathsf{A}}^{{\mathrm{ext}}}_{ij} = a(\phi_j^{{\mathrm{ext}}},\phi_i)$ for each valid $i,j$, a modified load vector ${\mathsf{f}}^{{\mathrm{ext}}}= {\mathsf{f}}- {\mathsf{A}}^{{\mathrm{ext}}}{\mathsf{u}}^{{\mathrm{ext}}}$ is used in computation.

It is obvious that $\nabla\, 1 = 0$ and so $a_1(1,v) = 0$ for any $v\in H^1\Omega)$. Therefore, by the partition of unity property $\sum_i \phi_i^{{\mathrm{ext}}}= 1$, the zero row sum of the matrix ${\mathsf{A}}^{{\mathrm{ext}}}$ also vanishes. Namely, $
    \sum_j {\mathsf{A}}_{ij}^{{\mathrm{ext}}}=
    0
    $. This property may be induced in the corresponding surrogate matrix if we simply define $
            \tilde{{\mathsf{A}}}_{ii}^{{\mathrm{ext}}}=
    - \sum_{j\neq i} \tilde{{\mathsf{A}}}_{ij}^{{\mathrm{ext}}}$, for every $x_i \in {\mathbb{X}}_m$, where $\tilde{{\mathsf{A}}}_{ij}^{{\mathrm{ext}}}= {\mathsf{A}}_{ij}^{{\mathrm{ext}}}$ for every $j$ where ${\mathsf{A}}_{ij}$ is not defined, and $\tilde{{\mathsf{A}}}_{ij}^{{\mathrm{ext}}}= \tilde{{\mathsf{A}}}_{ij}$ otherwise. With this extra condition, the surrogate matrix \[eq:SurrogateStiffnessMatrix\] actually requires one fewer independent stencil function; i.e., $\Phi^0_T = -\sum_{\delta \in {\mathscr{D}}(T_m)\setminus\{0\}} \Phi^\delta_T$, for every $T\in{{\mathcal{T}}}_H$. By this definition, although $\tilde{{\mathsf{A}}}$ does not satisfy the zero row sum property, the matrix ${\mathsf{A}}-\tilde{{\mathsf{A}}}$ does. Indeed, $$\label{eq:ZeroRowSumPropertySurrogateDifference}
    {\mathsf{A}}_{ii} - \tilde{{\mathsf{A}}}_{ii}
            =
    -
    \sum_{j\neq i} \big({\mathsf{A}}_{ij}^{{\mathrm{ext}}}- \tilde{{\mathsf{A}}}_{ij}^{{\mathrm{ext}}}\big)
    =
    -
    \sum_{j\neq i} \big({\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}\big)
    \,.$$ In this way, the stiffness matrix coming from linearized elasticity \[eq:BilinearFormLinearElasticity\] is similar to the stiffness matrix coming from the Laplacian. Indeed, the zero row sum property can be incorporated into its surrogate via a straight-forward generalization.

Analyzing the surrogate discretization {#sec:stability_of_the_surrogate_operator}
======================================

In this section, we define and motivate what we see as some the most essential features in the analysis of our surrogate methods. We begin with a review of discrete stability in the context of \[eq:SurrogateVF\]. We then touch on the concept of spectral convergence of the surrogate matrix $\tilde{{\mathsf{A}}}\to{\mathsf{A}}$, which helps us motivate the need to control $\|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}$. This specific quantity will repeatedly appear in the *a priori* error analysis in .

Discrete stability {#sub:stability}
------------------

Let $S = \{v\in V \,:\, \|v\|_V = 1\}$ be the surface of the unit ball in $V$. Recall \[eq:DiscreteVariationalProblems\] and assume that the discretization ${\mathsf{A}}{\mathsf{u}}= {\mathsf{f}}$ is stable. In the present context, this is equivalent to the existence of a constant $\alpha>0$ such that $$\sup_{v_h\in V_h\cap S} a(w_h,v_h) \geq \alpha \|w_h\|_V
    \quad
    \text{for all } w_h\in V_h
    \,.$$ Likewise, in order for the surrogate discretization $\tilde{{\mathsf{A}}}\tilde{{\mathsf{u}}} = {\mathsf{f}}$ to be stable, we must show that there exists a constant $\tilde{\alpha}>0$ such that $$\sup_{v_h\in V_h\cap S} \tilde{a}(w_h,v_h) \geq \tilde{\alpha}\|w_h\|_V
    \quad
    \text{for all } w_h\in V_h
    \,.
            \label{eq:DiscreteStability}$$ Inequality \[eq:DiscreteStability\] guarantees that $\tilde{{\mathsf{A}}}\tilde{{\mathsf{u}}} = {\mathsf{f}}$ has a unique solution and that $\|\tilde{u}_h\|_V \leq \tilde{\alpha}^{-1}\|F\|_{V^\ast}$. Equally important, however, it is a necessary precursor to Strang’s First Lemma, which in some cases can be used, in part, to show that $\tilde{u}_h$ converges to the exact solution $u$ (see, e.g., ).

Spectral convergence {#sub:general_principles}
--------------------

By analyzing the singular values of ${\mathsf{A}}$ directly, \[eq:DiscreteStability\] can sometimes be proven by showing that the spectrum of $\tilde{{\mathsf{A}}}$ converges to the spectrum of ${\mathsf{A}}$ at a fast enough rate. The main takeaway from this section is that spectral convergence can be guaranteed by showing that $\tilde{{\mathsf{A}}}\to {\mathsf{A}}$ in the matrix maximum norm, $\|\cdot\|_{\operatorname*{\vphantom{p}max}}$. Before moving on, denote the $k$-smallest eigenvalue of a matrix ${\mathsf{M}}\in{\mathbb{R}}^{N\times N}$ as $\lambda_k({\mathsf{M}})$ and let $\ell({\mathsf{M}}) = \operatorname*{\vphantom{p}max}_{1\leq i \leq N}\,\#\{{\mathsf{M}}_{ij} \neq 0 \text{ where } 1\leq j \leq N\}$ be the maximum number of nonzero components in ${\mathsf{M}}$, taken across all individual rows.

\[prop:SpectralConvergence\] Let ${\mathsf{M}},{\mathsf{N}}\in{\mathbb{R}}^{N\times N}$ be symmetric matrices. Then, for each $k=1,\ldots,N$, it holds that $$|\lambda_k({\mathsf{M}}) - \lambda_k({\mathsf{N}})|
        \leq
        \|{\mathsf{M}}-{\mathsf{N}}\|_{\infty}
        \,.$$

The proof is placed in . Because ${\mathsf{A}}$ and $\tilde{{\mathsf{A}}}$ have the same sparsity pattern, the next result follows trivially. Note that when $n=2$, $\ell({\mathsf{A}}-\tilde{{\mathsf{A}}}) \leq 7$, and when $n=3$, $\ell({\mathsf{A}}-\tilde{{\mathsf{A}}}) \leq 15$, whereas $\ell({\mathsf{A}})$ can be larger depending on the structure of the macro-mesh.

\[cor:BandwidthBound\] Let ${\mathsf{A}}$ and $\tilde{{\mathsf{A}}}$ be the true and surrogate stiffness matrices in \[eq:DiscreteVariationalProblems\], respectively. If both ${\mathsf{A}}$ and $\tilde{{\mathsf{A}}}$ are real symmetric matrices, then $$|\lambda_k({\mathsf{A}}) - \lambda_k(\tilde{{\mathsf{A}}})|
        \leq
        \ell({\mathsf{A}}-\tilde{{\mathsf{A}}})\cdot
        \|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}
        \,,
        \qquad
                k=1,\ldots,N
        .
    \label{eq:BandwidthBound}$$

As stated above, if $\|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}} \to 0$ fast enough, then can be used in proving the stability condition \[eq:DiscreteStability\]. However, \[eq:BandwidthBound\] is generally a very pessimistic bound and, when available, we recommend using more direct means to prove discrete stability (see, e.g., ). Nevertheless, this result illustrates the importance of controlling $\|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}$, which is a central feature in all of the coming analysis.

Controlling $\|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}$ with the variable coefficient Poisson equation
------------------------------------------------------------------------------------------------------------------------------------

Before we begin, some new notation is required. For any tensor $K:\Omega \to {\mathbb{R}}^{n\times n}$, define $\|K\|_{L^\infty(\Omega)} = \operatorname*{\vphantom{p}max}_{a,b} \|K_{ab}\|_{L^\infty(\Omega)}$, likewise, for any $r\geq 0$, define $|K|_{W^{r+1,\infty}(T)} = \operatorname*{\vphantom{p}max}_{a,b} |K_{ab}|_{W^{r+1,\infty}(T)}$. From now on, the notation $A \lesssim B$ will be used when two mesh-dependent quantities $A$ and $B$ satisfy an inequality $A \leq C B$, where $C$ some positive $H$-independent constant. Likewise, when $A \lesssim B$ and $B \lesssim A$, we write $A \eqsim B$. Recall that for a macro-mesh ${{\mathcal{T}}}_H$, the diameter of a single element $T\in{{\mathcal{T}}}$ is denoted $H_T$ and the mesh size is denoted $H = \operatorname*{\vphantom{p}max}_{T\in{{\mathcal{T}}}_H} H_T$. We also denote the fine-scale element diameter $h_T = 2^{-m}H_T$, for each $T \in {{\mathcal{T}}}_H$ and $h = 2^{-m}H$.

\[lem:InfinityNormConvergence\] Let ${\mathsf{A}}$ and $\tilde{{\mathsf{A}}}$, respectively, be the true and surrogate stiffness matrices corresponding to the bilinear form \[eq:BilinearFormPoisson\]. Namely, let each component of ${\mathsf{A}}$ be given by \[eq:TrueStiffnessMatrixPoisson\] and each component of $\tilde{{\mathsf{A}}}$ be defined by \[eq:SurrogateStiffnessMatrix\] with $G(x,u(y),v(y)) := {\nabla u(y)^\top K(x) \nabla v(y)}$. Fix $T\in {{\mathcal{T}}}_H$ and $0\leq r\leq q$. If $K_{ab}|_T\in W^{r+1,\infty}(T)$ for each index $a,b$, then

$$\big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max},T_m}
        \lesssim
        h_T^{n-2} H_T^{r+1}|K|_{W^{r+1,\infty}(T)}
        \,,
    \label{eq:InfinityNormBoundLocal}$$

where $\|{\mathsf{C}}\|_{{\operatorname*{\vphantom{p}max}},T_m} = \operatorname*{\vphantom{p}max}\big\{ |{\mathsf{C}}_{ij}| \,:\, x_i, x_j \in T_m \big\}$, for any matrix ${\mathsf{C}}$. Moreover, if each component $K_{ab}\in W^{r+1,\infty}({{\mathcal{T}}}_H) = \prod_{T\in{{\mathcal{T}}}_H} W^{r+1,\infty}(T)$, then $$\big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max}}
        \lesssim
        h^{n-2} H^{r+1}|K|_{W^{r+1,\infty}({{\mathcal{T}}}_H)}
        \,.
    \label{eq:InfinityNormBoundGlobal}$$

We prove only \[eq:InfinityNormBoundLocal\], \[eq:InfinityNormBoundGlobal\] then follows immediately. Recall \[eq:SurrogateStiffnessMatrix\] and fix $T\in{{\mathcal{T}}}_H$. Let $i$ and $j$ be the indices of the maximal value $|{\mathsf{A}}_{ij}-\tilde{{\mathsf{A}}}_{ij}| = \big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max},T_m}$. Next, because the theorem trivially holds in the degenerate case $\big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max},T_m} = 0$, we proceed under the assumption that ${\mathsf{A}}_{ij} \neq \tilde{{\mathsf{A}}}_{ij}$. Notably, it follows from \[cor:Interpolation\] that if each $K_{ab}|_T \in {\mathcal{P}}_q(T)$, then $\tilde{{\mathsf{A}}}_{ij} = {\mathsf{A}}_{ij}$ and, therefore, we find ourselves in the scenario where the diffusion tensor $K|_T$ is not a polynomial (of degree at most $q$). Here, we may also freely assume that $i\neq j$ because of \[eq:ZeroRowSumPropertySurrogateDifference\]. Indeed, for each $i$, $
        |{\mathsf{A}}_{ii}-\tilde{{\mathsf{A}}}_{ii}|
        \leq
                                \sum_{j\neq i}
        |{\mathsf{A}}_{ij}-\tilde{{\mathsf{A}}}_{ij}|
        \leq
        \ell({\mathsf{A}}-\tilde{{\mathsf{A}}})\cdot \operatorname*{\vphantom{p}max}_{j\neq i}{|{\mathsf{A}}_{ij}-\tilde{{\mathsf{A}}}_{ij}|}
            $.

To fix notation in the remainder of the proof, we take $\phi = \phi_i$, and express $$\tilde{{\mathsf{A}}}_{ij}
                        =
        \bigg[
        \Pi^\delta_T
        \int_{\Omega_{\delta}} \nabla \phi_{\delta}(y)^\top K(\cdot+y)\, \nabla \phi(y)\hspace{0.5pt} {\,\mathrm{d}}y
        \bigg](x_i)
                        ,$$ for some nonzero $\delta\in{\mathscr{D}}(T_m)$. Let ${\mathcal{I}}_T:C^0(\overline{T})\to {\mathcal{P}}_r(T)$ be the local Lagrange interpolant and define $\big[{\mathcal{I}}_T^{n\times n} K\big]_{ab} = {\mathcal{I}}_T K_{ab}$, for each index $a,b$. Splitting the two matrix entries ${\mathsf{A}}_{ij}$ and $\tilde{{\mathsf{A}}}_{ij}$ into polynomial and non-polynomial parts and rewriting $$\begin{aligned}
        \int_\Omega \nabla \phi_j(x)^\top \big[K-{\mathcal{I}}_T^{n\times n} K\big]{\hspace{-0.25pt}}(x){\hspace{0.25pt}}\nabla \phi_i(x) {\,\mathrm{d}}x
        &=
        \int_{\Omega_{\delta}} \nabla \phi_{\delta}(y)^\top \big[K-{\mathcal{I}}_T^{n\times n} K\big]{\hspace{-0.25pt}}(x_i+y){\hspace{0.25pt}}\nabla \phi(y) {\,\mathrm{d}}y
        \,,
        \\
        \int_\Omega \nabla \phi_j(x)^\top \big[{\mathcal{I}}_T^{n\times n} K\big](x){\hspace{0.25pt}}\nabla \phi_i(x) {\,\mathrm{d}}x
        &=
        \int_{\Omega_{\delta}} \nabla \phi_{\delta}(y)^\top \big[{\mathcal{I}}_T^{n\times n} K\big] (x_i+y){\hspace{0.25pt}}\nabla \phi(y)\hspace{0.5pt} {\,\mathrm{d}}y
        \,,
    \end{aligned}$$ we find that $$\begin{aligned}
        {\mathsf{A}}_{ij} &= \int_\Omega \nabla \phi_j^\top \,{\mathcal{I}}_T^{n\times n} K\; \nabla \phi_i {\,\mathrm{d}}x + \int_{\Omega_{\delta}} \nabla \phi_{\delta}(y)^\top \big[K-{\mathcal{I}}_T^{n\times n} K\big](x_i+y){\hspace{0.25pt}}\nabla \phi(y)\hspace{0.5pt} {\,\mathrm{d}}y,\\
        \tilde{{\mathsf{A}}}_{ij} &= \int_\Omega \nabla \phi_j^\top \,{\mathcal{I}}_T^{n\times n} K\; \nabla \phi_i {\,\mathrm{d}}x + \bigg[\Pi^\delta_T\int_{\Omega_{\delta}} \nabla \phi_{\delta}(y)^\top \big[K-{\mathcal{I}}_T^{n\times n} K\big](\cdot+y){\hspace{0.25pt}}\nabla \phi(y)\hspace{0.5pt} {\,\mathrm{d}}y\bigg](x_i)
                .
    \end{aligned}
    \label{eq:InterpolationStep}$$ Upon canceling the first two terms in the expressions above, we arrive at the inequality $$\big|{\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}\big| \leq |\beta_{ij}(x_i)| + |(\Pi^\delta_T\beta_{ij})(x_i)|
        \,,$$ where $\beta_{ij}(x) = \int_{\Omega_{\delta}} \nabla \phi_{\delta}(y)^\top \big[K-{\mathcal{I}}_T^{n\times n} K\big](x+y){\hspace{0.25pt}}\nabla \phi(y)\hspace{0.5pt} {\,\mathrm{d}}y$. Recall that the projection $\Pi^\delta_T:C^0(\overline{T}_\delta)\to{\mathcal{P}}(T_\delta)$ is continuous in the $L^\infty(\Omega)$ norm. Therefore, $$\begin{gathered}
                |(\Pi^\delta_T\beta_{ij})(x_i)| \leq \|\Pi^\delta_T\beta_{ij}\|_{L^\infty(T_\delta)} \lesssim  \|\beta_{ij}\|_{L^\infty(T_\delta)}\lesssim
        \|K-{\mathcal{I}}_T^{n\times n} K\|_{L^\infty(T)} \|\nabla\phi_{\delta} \cdot \nabla\phi\|_{L^1(\Omega_\delta)}
                \,.
    \end{gathered}$$ A standard scaling argument shows that $\|\nabla\phi_{\delta} \cdot \nabla\phi\|_{L^1(\Omega_\delta)} \lesssim h_T^{n-2}$. This, together with the well-known property $\|K_{ab}-{\mathcal{I}}_T K_{ab}\|_{L^\infty(T)} \lesssim H^{r+1}_T|K_{ab}|_{W^{r+1,\infty}(T)}$, yields the sufficient result $$\begin{aligned}
    |{\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}|
    \lesssim
    h_T^{n-2}H^{r+1}_T |K|_{W^{r+1,\infty}(T)}
    \,.
    \end{aligned}$$

The proof of can be read as a blueprint which extends to the settings of the bilinear forms $a_2(\cdot,\cdot)$, $a_3(\cdot,\cdot)$, defined in \[eq:BilinearFormLinearElasticity,eq:SplitBilinearFormspLaplace\]. Indeed, when $a_2(\cdot,\cdot)$ is considered, the only significant modification to the proof above is that an interpolation operator ${\mathcal{I}}_T:C^0(\overline{T})\to {\mathcal{P}}_r(T)$ must be introduced for each Lamé parameter $\mu$ and $\lambda$. The adaption to the setting $a(\cdot,\cdot) = a_3(\cdot,\cdot)$, is obvious. Ultimately, $$\big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max}}
    \lesssim
    h^{n-2} H^{r+1}
    \cdot
    \left\{
    \begin{alignedat}{3}
                                                                        \,&         \big|\lambda\big|_{W^{r+1,\infty}({{\mathcal{T}}}_H)} + \big|\mu\big|_{W^{r+1,\infty}({{\mathcal{T}}}_H)}
                        \qquad
        &&
        \text{if } 
        a(\cdot,\cdot) = a_2(\cdot,\cdot)
        \,,
        \\
        & \big|{\hspace{0.25pt}}|\nabla\tilde{u}|^{p-2}\big|_{W^{r+1,\infty}({{\mathcal{T}}}_H)}
                \qquad
        &&
        \text{if } 
        a(\cdot,\cdot) = a_3(\cdot,\cdot)
        \,.
    \end{alignedat}
    \right.
\label{eq:OtherMaxNormBounds}$$ Moreover, the proof may be easily modified to permit surrogates $\tilde{{\mathsf{A}}}$ without the zero row sum property. Likewise, scenarios involving fewer derivatives (which generally do not possess the zero row sum property), e.g., $a(\cdot,\cdot) = m(\cdot,\cdot)$, have similar bounds but invoke a different scaling in $h$.

*A priori* error estimation for the variable coefficient Poisson equation {#sec:a_priori_error_estimation}
=========================================================================

In this section, we present a thorough analysis of a surrogate discretization of the variable coefficient Poisson equation. Given a load $f\in L^2(\Omega)$ and symmetric positive-definite tensor $K:\Omega \to {\mathbb{R}}^{n\times n}$, the corresponding weak form may be written as follows. $$\text{Find } u\in H^1_0(\Omega)
    \text{ satisfying }
    \quad
                    a(u,v) = F(v)
    \quad
    \text{for all }
    v\in H^1_0(\Omega)
    \,,
\label{eq:PoissonWeakForm}$$ where $a(u,v) = \int_\Omega \nabla u^\top K{\hspace{0.25pt}}\nabla v {\,\mathrm{d}}x$ and $F(v) = \int_\Omega f\hspace{0.25pt} v {\,\mathrm{d}}x$. As done in , define $V_{h} = \{ v \in H^1_0(\Omega) : v|_t \in {\mathcal{P}}_1(t) \text{ for each } t\in{\mathcal{S}}^m({{\mathcal{T}}}_H)\}$ and let the corresponding vertex function basis be $\{\phi_i\}$.

Coercivity {#sub:coercivity}
----------

In the present setting, observe that each $v_h\in V_h$ can be expressed as $v_h(x) = \sum_i v_h(x_i)\phi_i(x)$. Therefore, due to the zero row/column sum property \[eq:ZeroRowSumPropertySurrogateDifference\], we find $$\label{eq:DifferenceIdentity}
    \begin{aligned}
        {\tilde{a}}({{v}_{h}}, {{w}_{h}}) - {a}({{v}_{h}}, {{w}_{h}})
        &=
        \sum_{i,j} 
        \big(\tilde{{\mathsf{A}}}_{ij} - {\mathsf{A}}_{ij}\big){\hspace{0.25pt}}v_h(x_i){\hspace{0.25pt}}w_h(x_j)
                                                \\
        &=
        \frac{1}{2}\sum_{i\neq j} \big({\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}\big) (v_h(x_i) - v_h(x_j)){\hspace{0.25pt}}(w_h(x_i) - w_h(x_j))
        \,.
    \end{aligned}$$ Due to the mutual sparsity of the matrices $\tilde{{\mathsf{A}}}$ and ${\mathsf{A}}$, every nonzero term in the sum above can be rewritten as $
    \big({\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}\big) (v_h(x_i) - v_h(x_i+\delta)){\hspace{0.25pt}}(w_h(x_i) - w_h(x_i+\delta))
    ,
$ for some nonzero $\delta$. Because $|\delta| \eqsim h$ by construction, one easily arrives at the following upper bound: $$\begin{aligned}
    {a}({{v}_{h}}, {{w}_{h}}) - {\tilde{a}}({{v}_{h}}, {{w}_{h}})
    \lesssim
    h^{2-n}
    \big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{{\operatorname*{\vphantom{p}max}}}{\hspace{0.25pt}}\|\nabla v_h\|_{0}\hspace{1pt} \|\nabla w_h\|_{0}
                \,.
                        \label{eq:FiniteDifferenceBound}\end{aligned}$$ We now arrive at the main result of this subsection.

\[cor:Coercivity\] Let $0\leq r\leq q$. Assume that $a(\cdot,\cdot)$ is coercive and $K\in \big[W^{r+1,\infty}({{\mathcal{T}}}_H)\big]^{n\times n}$. Then, for any fine enough macro-mesh ${{\mathcal{T}}}_H$, the surrogate bilinear form $\tilde{a}:V_h\times V_h \to {\mathbb{R}}$ is also coercive.

Let $S = \{v\in H^1 \,:\, \|v\|_1 = 1\}$ be the surface of the unit ball in $H^1$. Recall that since $a(\cdot,\cdot)$ is coercive, there exists a coercivity constant $\alpha>0$ such that $a(v,v) \geq \alpha$ for all $v\in S$. Notice that $\alpha \leq a(v_h,v_h) \leq \tilde{a}(v_h,v_h) + |a(v_h,v_h) - \tilde{a}(v_h,v_h)|$ for all $v_h \in V_h\cap S$ and, therefore, $$\alpha
        -
                |a(v_h,v_h) - \tilde{a}(v_h,v_h)|
        \leq
                \tilde{a}(v_h,v_h)
        \quad
        \text{for all } v_h \in V_h\cap S
        \,.
    \label{eq:Coercivity1}$$ Here, the second term on the left may be bounded from above using \[eq:FiniteDifferenceBound,lem:InfinityNormConvergence\] as follows, $$|a(v_h,v_h) - \tilde{a}(v_h,v_h)|
                                                \lesssim
        h^{2-n}\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}} \|\nabla v_h\|_0^2
        \lesssim
        H^{r+1}|K|_{W^{r+1,\infty}(\Omega)}
        \,.
        \label{eq:Coercivity2}$$ Thus, for any small enough $H$, we see that $0 < \alpha - |a(v_h,v_h) - \tilde{a}(v_h,v_h)| \leq \tilde{a}(v_h,v_h)$, as necessary.

Convergence of the surrogate solution in the $H^1$ norm {#sub:convergence_rates_of_the_solution_in_the_h_1_norm}
-------------------------------------------------------

The purpose of this subsection is to derive a mesh-dependent upper bound on the error in the surrogate solution $\tilde{u}_h$ of the form $\|u - {{\tilde{u}}_{h}}\|_1 \leq C(K,\Omega,u){\hspace{0.25pt}}h + \tilde{C}(K,\Omega,u)H^{r+1}$. In doing so, we choose to emphasize the primary difference from the classical $\|u - u_h\|_1 \leq C(K,\Omega,u){\hspace{0.25pt}}h$ error estimate by absorbing the coercivity and continuity constants (which depend on both $K$ and $\Omega$) into the $\lesssim$ symbol. We begin with a particular version of the First Strang Lemma [@strang1972variational].

\[lem:Strang1\] Let $S = \{v\in H^1 \,:\, \|v\|_1 = 1\}$ be the surface of the unit ball in $H^1$. Assume that $\tilde{a}:V_{h}\times V_{h}\to {\mathbb{R}}$ is coercive. The following error estimate holds for the surrogate solution $\tilde{u}_h$ of the variable coefficient model problem \[eq:PoissonWeakForm\]: $${{\left\|{u} - {{\tilde{u}}_{h}}\right\|_{1}}}
        \lesssim
        \operatorname*{\vphantom{p}inf}_{w_h\in V_h}
        \Big[
                        \|{u} - w_h\|_1
            +
            \sup_{v_h\in V_h\cap S}
            |\tilde{a}(w_h,v_h) - a(w_h,v_h)|
        \Big]
        \,.
    \label{eq:Strang1}$$

\[thm:APrioriBoundH1\] Let $0\leq r\leq q$ and assume that $K \in \big[W^{r+1,\infty}(\Omega)\big]^{n\times n}$ is symmetric and positive definite with $\lambda_1(K)$ bounded away from zero almost everywhere. Let $u\in H^1(\Omega)$ and $\tilde{u}_h\in V_h$ be the unique solutions to \[eq:ContinousVF,eq:SurrogateVF\], respectively, where $a(u,v) = \int_\Omega \nabla u^\top K{\hspace{0.25pt}}\nabla v {\,\mathrm{d}}x$ and $F(v) = \int_\Omega f\hspace{0.25pt} v {\,\mathrm{d}}x$. Then, for any sufficiently fine macro-mesh ${{\mathcal{T}}}_H$, the following upper bound holds: $$\|u - {{\tilde{u}}_{h}}\|_{1} \lesssim h{\hspace{0.25pt}}|u|_{2} + H^{r+1}|K|_{W^{r+1,\infty}(\Omega)} |u|_{1}
        \,.
    \label{APrioriBoundH1Local}$$

With the assumptions above, $a(\cdot,\cdot)$ is coercive. Therefore, by , if the macro-mesh ${{\mathcal{T}}}_H$ is taken fine enough, then $\tilde{a}:V_{h}\times V_{h}\to {\mathbb{R}}$ is coercive. We now bound the right-hand side of \[eq:Strang1\]. Invoking \[eq:FiniteDifferenceBound\], we find $${{\left\|{u} - {{\tilde{u}}_{h}}\right\|_{1}}} \lesssim \|{u} - w_h\|_1 + h^{2-n} \big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max}} \|\nabla w_h\|_{0}
        \,,$$ for every $w_h\in V_h$. Setting $w_h = {\mathcal{S}}\hspace{-0.4pt}{\mathcal{Z}}_h u$, the Scott–Zhang interpolant of $u$ [@scott1990finite], we see that $${{\left\|{u} - {{\tilde{u}}_{h}}\right\|_{1}}}
        \lesssim
        \|{u} - {\mathcal{S}}\hspace{-0.4pt}{\mathcal{Z}}_h u\|_1 + h^{2-n} \big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max}} |{\mathcal{S}}\hspace{-0.4pt}{\mathcal{Z}}_h u|_{1}
        \lesssim
        h{\hspace{0.25pt}}|u|_{2} + h^{2-n}\big\|{\mathsf{A}}- \tilde{{\mathsf{A}}}\big\|_{\operatorname*{\vphantom{p}max}} |u|_{1}
                \,.$$ In order to finish the proof, recall that $h^{2-n} \|{\mathsf{A}}- \tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}} \lesssim H^{r+1}|K|_{W^{r+1,\infty}(\Omega)}$, by . 

Convergence of the surrogate solution in the $L^2$ norm {#sub:convergence_of_the_solution_in_the_}
-------------------------------------------------------

In this subsection, we prove an $L^2$ error estimate of the form $\|u - {{\tilde{u}}_{h}}\|_0 \leq C(K,\Omega,u){\hspace{0.25pt}}h^2 + \tilde{C}(K,\Omega,u)H^{r+1}$. A second result, which elicits accelerated $H$-convergence, is also proved under the additional assumption $\sum_{x_i,\in T^\delta_m} \big[\Phi^\delta_T-\Pi^\delta_T\Phi^\delta_T\big](x_i) = 0$, where $T^\delta_m = T_m \cap \overline{T}_\delta$. This is a property which naturally arises for the specific class of least-squares projections introduced in . Again, we emphasize the primary differences from the corresponding classical error estimate by absorbing the standard constants into the $\lesssim$ symbol.

\[thm:APrioriBoundL2\] Under the conditions of , if $\Omega\operatorname{\subseteq}{\mathbb{R}}^n$ is a convex domain, then the following additional upper bound on the error in the surrogate solution holds:

\[eq:APrioriBoundL2\] $$\|u - {{\tilde{u}}_{h}}\|_0
        \lesssim
        h^2{\hspace{0.25pt}}|u|_{2} + H^{r+1}|K|_{W^{r+1,\infty}(\Omega)} |u|_{1}
                \,.
    \label{eq:APrioriBoundL2_1}$$ Moreover, if $r>0$ and $\sum_{x_i\in T^\delta_m} \big[\Phi^\delta_T-\Pi^\delta_T\Phi^\delta_T\big](x_i) = 0$, for each $T\in{{\mathcal{T}}}_H$ and $\delta\in{\mathscr{D}}(\mathring{T}_m)$, then $$\|u - {{\tilde{u}}_{h}}\|_0
        \lesssim
        h^2{\hspace{0.25pt}}|u|_{2} + H^{r+2}|K|_{W^{r+1,\infty}(\Omega)} \|\nabla u\|_{1}
        \,.
    \label{eq:APrioriBoundL2_2}$$

By the triangle inequality, $\|u-\tilde{u}_h\|_0 \leq \|u-u_h\|_0 + \|u_h-\tilde{u}_h\|_0$, where $u_h\in V_h$ is the discrete solution coming from \[eq:DiscreteVF\]. It can be shown that if $\Omega$ is convex, then $u\in H^2(\Omega)\cap H^1_0(\Omega)$; see, e.g., [@grisvard2011elliptic]. It then follows from standard arguments that $\|u-u_h\|_0 \lesssim h^2{\hspace{0.25pt}}|u|_2$; see, e.g., [@brenner2007mathematical Theorem 5.7.6]. Therefore, we only need to analyze the term $\|u_h-\tilde{u}_h\|_0$. Since, $\|u_h-\tilde{u}_h\|_0 \leq \|u_h-\tilde{u}_h\|_1$, proceeding as in the proof of , we quickly arrive at \[eq:APrioriBoundL2\_1\].

In order to prove \[eq:APrioriBoundL2\_2\], first define $w_h\in V_h$ satisfying $a(w_h,v_h) = (u_h-\tilde{u}_h,v_h)_\Omega$, for all $v_h\in V_h$. Observe that the exact solution of the problem $a(w,v) = (u_h-\tilde{u}_h,v)_\Omega$, for all $v\in H^1_0(\Omega)$, belongs to the space $H^2(\Omega)\cap H^1_0(\Omega)$, $\|w\|_2\lesssim\|u_h-\tilde{u}_h\|_{0}$. Moreover, $$\begin{aligned}
        \|u_h-\tilde{u}_h\|_0^2
        &=
        a(w_h,u_h-\tilde{u}_h)
        =
        \tilde{a}(w_h,\tilde{u}_h) - a(w_h,\tilde{u}_h)
        \\
        &=
        \frac{1}{2}
        \sum_{i\neq j} ({\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}) (\tilde{u}_h(x_i) - \tilde{u}_h(x_j)) (w_h(x_i) - w_h(x_j))
        \,,
        \end{aligned}
    \label{eq:SecondTermManipulationAubinNitsche}$$ where the final line follows from \[eq:DifferenceIdentity\]. As remarked previously, each nonzero term in this sum can be written as $\big({\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij}\big) (\tilde{u}_h(x_i) - \tilde{u}_h(x_i+\delta)){\hspace{0.25pt}}(w_h(x_i) - w_h(x_i+\delta))$, for some $T\in{{\mathcal{T}}}_H$ and nonzero $\delta \in {\mathscr{D}}(T_m)$. We can make better use of this expression with the identity $v_h(x_i) - v_h(x_i+\delta) = \nabla v_h(y_{i,\delta})$, wherein each $y_{i,\delta}$ is chosen from the edge connecting $x_i$ and $x_i+\delta$, and with the relationship ${\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij} = \big[\Phi^\delta_T-\Pi^\delta_T\Phi^\delta_T\big](x_i)$, since ${\mathsf{A}}_{ij} - \tilde{{\mathsf{A}}}_{ij} \neq 0$ and $i\neq j$. With these observations in hand, we have $$\begin{aligned}
        2\|u_h&-\tilde{u}_h\|_0^2
        =
        \sum_{T\in{{\mathcal{T}}}_H}
        \sum_{\delta\in {\mathscr{D}}(\mathring{T}_m)}
        \sum_{x_i\in T^\delta_m}
        \big[\Phi^\delta_T-\Pi^\delta_T\Phi^\delta_T\big](x_i)\,
        (\nabla \tilde{u}_h(y_{i,\delta}) \cdot \delta)(\nabla w_h(y_{i,\delta})\cdot \delta)
        \\
        &\leq
        \sum_{T\in{{\mathcal{T}}}_H}
        \sum_{\delta\in {\mathscr{D}}(\mathring{T}_m)}
        \sum_{x_i\in T^\delta_m}
        \big[\Phi^\delta_T-\Pi^\delta_T\Phi^\delta_T\big](x_i)
        \Big(
        (\nabla \tilde{u}_h(y_{i,\delta}) \cdot \delta)(\nabla w_h(y_{i,\delta})\cdot \delta) - C
        \Big)
        \\
        &\lesssim
        h^{-n} 
        \|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}
        \sum_{T\in{{\mathcal{T}}}_H}
        \sum_{\delta\in {\mathscr{D}}(\mathring{T}_m)}
        \|(\nabla \tilde{u}_h \cdot \delta)(\nabla w_h\cdot \delta)-C\|_{L^1(T)}
        \,.
    \end{aligned}$$ If we set the constant above to the following average value, $C := \frac{1}{\mathrm{vol}(T)} \int_T(\nabla u \cdot \delta)(\nabla w\cdot \delta) {\,\mathrm{d}}x$, then $\|(\nabla u \cdot \delta)(\nabla w\cdot \delta) - C\|_{L^1(T)} \lesssim H_T{\hspace{0.25pt}}|(\nabla u \cdot \delta)(\nabla w\cdot \delta)|_{W^{1,1}(T)}$. Therefore, $$\begin{aligned}
        \|(\nabla \tilde{u}_h \cdot \delta)&(\nabla w_h\cdot \delta)-C\|_{L^1(T)}
        \leq
        \|(\nabla \tilde{u}_h \cdot \delta)(\nabla w_h\cdot \delta) - (\nabla u \cdot \delta)(\nabla w_h\cdot \delta)\|_{L^1(T)}
        \\
        &+
        \| (\nabla u \cdot \delta)(\nabla w_h\cdot \delta) - (\nabla u \cdot \delta)(\nabla w\cdot \delta)\|_{L^1(T)}
        +
        \|(\nabla u \cdot \delta)(\nabla w\cdot \delta) - C\|_{L^1(T)}
        \\
        &\lesssim
        \|\nabla (u - \tilde{u}_h) \cdot \delta\|_{0,T}\,\|\nabla w_h\cdot \delta\|_{0,T}
        +
        \|\nabla (w - w_h)\cdot \delta\|_{0,T}\, \|\nabla u \cdot \delta\|_{0,T}
        \\
        &+
        H_T{\hspace{0.25pt}}\|\nabla u \cdot \delta\|_{1,T}{\hspace{0.25pt}}\|\nabla w\cdot \delta\|_{1,T}
        \,.
    \end{aligned}$$ After summing over all $T\in {{\mathcal{T}}}_H$ and $\delta\in{\mathscr{D}}(T_m)$ and taking into account $|\delta| \eqsim h$, we arrive at the bound $$\begin{aligned}
        \|u_h-\tilde{u}_h\|_0^2
        &\lesssim
        h^{2-n} 
        \|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}
        \big(
        \|u - \tilde{u}_h\|_{1}\, | w_h|_{1}
        +
        \|w - w_h\|_{1}\, |u|_{1}
        +
        H\|\nabla u\|_{1} \|\nabla w\|_{1}
        \big)
        \\
        &\lesssim
        h^{2-n} 
        \|{\mathsf{A}}-\tilde{{\mathsf{A}}}\|_{\operatorname*{\vphantom{p}max}}
        \big(
        h\,|u|_2 + H^{r+1}|K|_{W^{r+1,\infty}(\Omega)}|u|_1 + H\, \|\nabla u\|_{1}
        \big)
        {\hspace{0.25pt}}\|w\|_{2}
        \\
        &\lesssim
        H^{r+1}|K|_{W^{r+1,\infty}(\Omega)}
        \big(
        H\, \|\nabla u\|_{1} + H^{r+1}|K|_{W^{r+1,\infty}(\Omega)}|u|_1
        \big)
        {\hspace{0.25pt}}\|w\|_{2}
        \,.
    \end{aligned}$$ The proof is completed by recalling that $\|w\|_2\lesssim\|u_h-\tilde{u}_h\|_{0}$.

As stated previously, the error estimates in not do track the constants appearing in the corresponding classical estimates. When $r>0$, none of the constants in either classical estimate depend on the higher order norms $|K|_{W^{r+1,\infty}(\Omega)}$.

The main ingredients in the error analysis presented in this section are , Strang’s First Lemma, and \[eq:DifferenceIdentity\], which simply follows from the zero row sum property and symmetry. Simple generalizations of , for the bilinear forms $a_2(\cdot,\cdot)$ and $a_3(\cdot,\cdot)$ have been stated in \[eq:OtherMaxNormBounds\]. Meanwhile, \[eq:DifferenceIdentity\] simply follows from the zero row sum property of $\tilde{{\mathsf{A}}}$ and symmetry. Therefore, we see no reason to doubt that our analysis here can be generalized to the other problems of interest in this paper. Hence, we proceed with numerical verification and demonstration.

Implementation {#sec:implementation}
==============

The performance of a numerical method largely depends on its implementation. Therefore, in this section, we highlight the important features of ours. We use the HyTeG finite-element software framework [@Kohl2018HyTeGfinite] as the core framework for all the numerical experiments in . It offers efficient distributed data structures for simplicial meshes in 2D and 3D, which serve as a basis for the implementation of massively parallel fast iterative solvers. Its main concept is based on the idea that a coarse input mesh is split into its geometrical primitives, i.e., vertices, edges, and faces, and each of these primitives is uniformly refined. Because the primitives of the same dimension are decoupled from the others, all primitives of the same dimension may be processed in parallel. This partitioning and the hierarchy of locally structures meshes allows for efficient parallel implementations of geometric multigrid methods. More importantly, these data structures fit perfectly to the concept of macro-elements introduced in . The problems in this paper are mainly solved by employing a geometric multigrid solver using V-cycles with a hybrid Gauss–Seidel smoother. On the coarsest grid, either a preconditioned conjugate gradient method or the direct solver MUMPS [@MUMPS01; @MUMPS02], as provided by the PETSc interface [@petsc-user-ref; @petsc-efficient], is used. For improved parallel scalability of the coarse grid solver, agglomeration techniques as provided by `PCTELESCOPE` [@MaySananRuppKnepleySmith2016] are used in runs with many processes.

Polynomial least squares regression {#sub:PLSR}
-----------------------------------

An important factor in the performance of the surrogate approach is the approximation of the stencil functions $\Phi^\delta_T$ by polynomials $\tilde{\Phi}^\delta_T$. This step in the solver process must be very fast and is usually done in a pre-process step before the actual solve. After various preparatory experiments, we have seen satisfactory performance and accuracy from simply computing $\tilde{\Phi}^\delta_T = \Pi^\delta_T \Phi^\delta_T$ via solving a simple least-squares problem, which we now describe.

Let $T\in{{\mathcal{T}}}_H$ be a macro-element and recall that $T_m$ is the associated lattice on level $m$. Suppose that $\Phi_T^\delta$ is the stencil function in direction $\delta \in {\mathscr{D}}(T_m)$ which we want to approximate. For the least-squares regression, we fix a level $m_{\mathrm{LS}}$ with $m \geq m_{\mathrm{LS}} \geq 2$ and define the set of least-squares points ${T}{}^\delta_{\mathrm{LS}} \coloneqq {T}_{m_{\mathrm{LS}}}\cap \overline{T}_\delta$. Furthermore, let $\{p_k\}_{k=1}^{M}$ be a basis of $\mathcal{P}_q(T)$, the space of polynomials with maximal degree $q$. Assume that $m_{\mathrm{LS}}$ is chosen large enough such that $\big|{T}{}^\delta_{\mathrm{LS}}\big| \geq M$ and introduce the following norm on $\mathcal{P}_q(T)$: $\|p\|^2_{{T}{}^\delta_{\mathrm{LS}}} \coloneqq \sum_{x_i \in {T}{}^\delta_{\mathrm{LS}}} p(x_i)^2$. The least-squares regression problem, which in turn defines $\Pi^\delta_T$, is formalized as follows: $$\begin{aligned}
\label{eqn:leastsquaresproblem}
\text{Find } {\mathsf{c}}\in \mathbb{R}^M \text{ satisfying }
\quad
{\mathsf{c}}= \operatorname*{arg\,min}_{{\mathsf{d}}\in \mathbb{R}^M} \left\|\Phi^\delta_{T} - \sum_{k=1}^{M} d_k p_k\right\|^2_{{T}{}^\delta_{\mathrm{LS}}}
\,.\end{aligned}$$ The approximated stencil function is then defined as $\tilde{\Phi}^\delta_{T} \coloneqq \sum_{k=1}^{M} c_k p_k$. This problem is equivalent to solving the possibly overdetermined linear system of equations ${\mathsf{B}}{\mathsf{c}}= {\mathsf{f}}$ in a least-squares sense, where ${\mathsf{B}}_{ij} = p_j(x_i)$ and ${\mathsf{f}}_i = \Phi^\delta_{T}(x_i)$ for $1 \leq i \leq |{T}{}^\delta_{\mathrm{LS}}|$ and $1 \leq j \leq N$. The choice of the polynomial basis is arbitrary. However, for an easier implementation, we employ the monomial basis, even knowing that the resulting linear system is ill conditioned. Since it is crucial to get numerically precise results, a stable solver for this problem has to be chosen. For this purpose, we apply the `colPivHouseholderQr` method from the Eigen 3.3.5 library [@eigenweb], which offers a good balance between speed and accuracy. Obviously, each of these linear systems is independent of others, therefore they may be solved in parallel.

\[rem:OptimalityCondition\] Taking into account $\tilde{\Phi}^\delta_{T} = \sum_{k=1}^{M} c_k p_k$, the first order optimality condition for \[eqn:leastsquaresproblem\] can be stated as $\sum_{x_i\in{T}{}^\delta_{\mathrm{LS}}} \big[\Phi^\delta_T-\Pi^\delta_T\Phi^\delta_T\big](x_i) = 0$. If $m=m_{\mathrm{LS}}$, then the secondary assumption in is satisfied and we see higher order convergence in $H$, as stated in \[eq:APrioriBoundL2\_2\]. Usually, when $m_{\mathrm{LS}}$ is close but not equal to $m$, we see preasymptotic $H$-convergence in between the two estimates given in \[eq:APrioriBoundL2\]; see .

\[remark:symmetricregression\] In the case where the bilinear form $a(\cdot,\cdot)$ is symmetric, we need only approximate a single stencil function $\Phi^\delta_T$ for both directions $\delta$ and $-\delta$. Indeed, as observed in , the corresponding stencil functions are identical, up to a shift by $\delta$. Furthermore, the symmetry requirement \[eq:SymmetricInterpolants\], from , is satisfied with the projection operator defined above. Indeed, one may verify that for every $\delta$, $T_\delta = T_{-\delta} - \delta$. Therefore, $$\Big\|\Phi^\delta_{T} - \tilde{\Phi}^\delta_{T}\Big\|^2_{{T}{}^\delta_{\mathrm{LS}}}
    =
    \Big\|\Phi^{-\delta}_{T} - \tilde{\Phi}^{-\delta}_{T}\Big\|^2_{{T}{}^{-\delta}_{\mathrm{LS}}}
    \qquad
    \text{and}
    \qquad
    \tilde{\Phi}^\delta_{T}(x) = \tilde{\Phi}^{-\delta}_{T}(x+\delta)
    \,.$$ Thus, on simplicial meshes in 2D, only four instead of seven polynomials per macro-element have to be determined and stored in memory. In some cases, where the zero row sum property holds, the number of required polynomials may be even reduced to three.

Fast polynomial evaluation
--------------------------

An even more important factor with respect to the performance of our implementation is the fast evaluation of the surrogate stencil functions $\tilde{\Phi}^\delta_{T}$. Contrary to the computation of each $\tilde{\Phi}^\delta_{T}$, which will happen only once per solve, these evaluations will be made during every matrix-vector multiplication. Therefore, the costs of evaluating the stiffness matrix entries associated to a degree of freedom, may not exceed the costs of evaluating the bilinear forms with the respective ansatz functions. In this case not only the reduction of floating point operations per degree of freedom is of importance, but also the required memory traffic has to be taken into account.

When performing a matrix-vector multiplication in HyTeG, the degrees of freedom in a macro-element are processed in a row-wise fashion as illustrated in the left of \[fig:poly\_eval\_1d\]. In each row, the stencil function may be interpreted as a 1D function. We assume without loss of generality, that the 1D stencil functions are aligned with the $x$-axis. This property is also inherited by the approximated stencil function.

![\[fig:poly\_eval\_1d\] Illustration of a loop through the degrees of freedom in a macro-element. In each row of the loop, the 2D stencil function may be interpreted as a 1D function (left). Seven stencil functions have to be evaluated in order to obtain the whole stencil for a degree of freedom (right).](poly_eval.pdf "fig:"){width="0.33\linewidth"} ![\[fig:poly\_eval\_1d\] Illustration of a loop through the degrees of freedom in a macro-element. In each row of the loop, the 2D stencil function may be interpreted as a 1D function (left). Seven stencil functions have to be evaluated in order to obtain the whole stencil for a degree of freedom (right).](stencileval.pdf "fig:"){width="0.375\linewidth"}

To further optimize the evaluation of the 1D polynomial, we exploit that the stencil functions have to be evaluated on a line subdivided into uniformly sized intervals of length $h$. Let $(x_i, y_j)$ be a vertex node in the lattice and let $p_{y_j}(\cdot) \coloneqq \tilde{\Phi}_\delta(\cdot, y_j)$ be the approximated 1D stencil function associated to row $y_j$.

Assuming that we already have evaluated the stencil function $p_{y_j}$ at a point $x_i$, we want to evaluate it at the next point $x_i + h$ as efficiently as possible. Since the grid points are equidistantly distributed, we can use a special case of the divided differences, called forward differences [@RichardL.Burden2015Numerical pg. 126].

First, we need $q+1$ helper variables $\left\{\Delta_{x_0}^{(k)}\right\}$ for $k \in \{0,1,\ldots,q\}$ which are defined in a preprocessing step as follows: $$\begin{aligned}
\Delta^{(0)}_{x_0} &\coloneqq p_{y_j}(x_0)\,,\\
\Delta^{(k)}_{x_0} &\coloneqq \Delta_{x_0+h}^{(k-1)} - \Delta_{x_0}^{(k-1)}, \quad k \in \{1,\ldots,q\}\,.\end{aligned}$$

The value at position $p_{y_j}(x_0)$ is then given by $\Delta_{x_0}^{(0)}$. In order to obtain the value at $p_{y_j}(x_0+h)$, one has to update all the helper variables in the following way: $$\begin{aligned}
\Delta_{x_0}^{(k)} &\coloneqq \Delta_{x_0}^{(k)} + \Delta_{x_0}^{(k+1)}, \quad k \in \{0,1,\ldots,q\}\,.\end{aligned}$$ After that, the value of $p_{y_j}(x_0+h)$ is given by $\Delta_{x_0}^{(0)}$. Doing this recursively yields the approximated stencil function values at all mesh points on a single row using only $q+1$ helper variables and $q+1$ floating point additions. When iterating through a row, in general seven polynomial evaluations, one for each direction, are required, cf. right of \[fig:poly\_eval\_1d\]. In the symmetric case this may be reduced to six polynomial evaluations, since the western stencil weight may be obtained from the previous eastern evaluation. Keep also in mind that in the symmetric case the polynomials of approximated stencil functions in opposite directions are the same but only evaluated at different positions, cf. \[remark:symmetricregression\]. Therefore, $6\cdot(q+1)$ helper variables are required for a single row. For $q=8$, these results in $54\cdot8$ bytes of memory which fits easily into a modern L1 CPU cache. When moving from one lattice point to another, $6\cdot(q+1)$ vectorizable floating point additions have to be performed, to obtain the updated polynomial evaluations. Furthermore, in our implementation, the polynomial degree is realized as a C++ template parameter, therefore, all loops concerning the evaluation of a polynomial of a certain degree may be optimized at compile time. Since our focus lies in the theoretical analysis of the surrogate approach, thorough performance studies employing performance models should be considered beyond the scope of this paper. Similar performance studies have been carefully completed in [@bauer2018new; @bauer2017two]. Thus, in the next section, we only report on relative run-times of the surrogate approach compared to the standard method, also implemented on HyTeG, using on-the-fly quadrature of the integrals stemming from the bilinear forms.

Numerical experiments {#sec:numerical_experiments}
=====================

In this section, we numerically verify , both related to the variable coefficient Poisson equation. Additionally, we present proof-of-concept results for a linearized elasticity application and a simple $p$-Laplacian diffusion problem. While not covered by the theory, we include these latter examples to demonstrate the breadth of generality of the methodology.

All run-time measurements in this sections were obtained on a machine equipped with two Intel^^ Xeon^^ Gold 6136 processors with a nominal base frequency of 3.0 GHz. Each processor has 12 physical cores which results in a total of 24 physical cores. The total available memory of is split into two NUMA domains, one for each socket. We use the GCC 7.3.0 compiler and specify the following compiler arguments: `-O3 -march=native`. All the examples in this section were executed in parallel using all available 24 physical cores.

When comparing run times from the standard and the surrogate approaches, many factors are responsible for the relative speed-up of the surrogate approach. Increasing the polynomial order $q$ of the surrogate stencil functions not only increases the run time of a multigrid iteration but also the time spent in the setup phase (i.e., computing each $\Phi_T^\delta$). The cost of the setup phase, however, is mostly dominated by the sampling level $m_{\mathrm{LS}}$. Therefore, when the ratio of time spent in the iterative solver to the time spent in the setup phase(s) for solving a particular problem is large, the setup cost is almost negligible and we see the best performance. Since the problems in the following subsections differ in complexity and have different ratios of solver to setup time, the observed relative speed-ups are not directly comparable. Nonetheless, the reported speed-ups for all tested examples range between a factor of 14 and 20. Such significant speed-ups are in particular important in case of dynamic or stochastic applications. Most stochastic applications demand an enormous number of deterministic solves resulting quite often in extreme long run times. Having such a surrogate approach at hand can help to make stochastic approaches such as, e.g., multilevel Monte Carlo and its variants, more accessible for complex applications.

Quantitative benchmark problem {#sub:quantitativebenchmark}
------------------------------

In this subsection, we examine the surrogate method for the variable coefficient Poisson equation which has been described and analyzed above. The strong form of the problem is $$\begin{alignedat}{2}
-\mathrm{div}\left(K \nabla u\right) &= f &&\quad \text{in } \Omega,\\
u &= g &&\quad \text{on } \partial \Omega.
\end{alignedat}
\label{eqn:strongform}$$ We consider both the bilinear form coming from the scalar coefficient scenario (i.e., $K = k\cdot \mathrm{Id}$), introduced in \[eq:BilinearFormPoisson\], and the tensorial coefficient scenario, introduced in \[eq:BilinearFormPoissonTensor\]. In the scalar coefficient experiments, we use the unit-square domain $\Omega = (0,1)^2$. In the tensorial coefficient experiments, the domain $\Omega$ has a curvilinear boundary.

![\[fig:unitsquare\_meshes\] Coarse macro-meshes of the unit-square with mesh sizes $H = H_0$ (left), $H = \nicefrac{H_{0}}{2}$ (middle), and $H = \nicefrac{H_{0}}{4}$ (right). Meshes with a smaller $H$ follow the same uniform refinement pattern.](unitsquare_macro_elements_level0.pdf){width="\linewidth"}

![\[fig:unitsquare\_meshes\] Coarse macro-meshes of the unit-square with mesh sizes $H = H_0$ (left), $H = \nicefrac{H_{0}}{2}$ (middle), and $H = \nicefrac{H_{0}}{4}$ (right). Meshes with a smaller $H$ follow the same uniform refinement pattern.](unitsquare_macro_elements_level1.pdf){width="\linewidth"}

![\[fig:unitsquare\_meshes\] Coarse macro-meshes of the unit-square with mesh sizes $H = H_0$ (left), $H = \nicefrac{H_{0}}{2}$ (middle), and $H = \nicefrac{H_{0}}{4}$ (right). Meshes with a smaller $H$ follow the same uniform refinement pattern.](unitsquare_macro_elements_level2.pdf){width="\linewidth"}

### Scalar coefficient on unit square {#sec:scalar_coefficient_benchmark}

In the first benchmark problem ($K = k\cdot \mathrm{Id}$), we take $\Omega = (0,1)^2$ and employ the scalar coefficient function $$\begin{aligned}
k(x,y) = \exp{\left(xy\right)} + \sin{\left (3\pi x y \right )} + \cos{\left (\pi x^{2} y \right )} + 1
\label{eq:ScalarCoefficient}\end{aligned}$$ in problem \[eqn:strongform\]. The manufactured solution $u$ is chosen as $u(x,y) = \sin(x)\sinh(y)$. The restriction of $u$ to the boundary is chosen as Dirichlet datum $g$. The right-hand-side $f$ is directly computed by inserting $u$ into the equation. In this benchmark, we fix the finest mesh size $h$ and report on the errors depending on $H$ and $q$ to show the proven $\mathcal{O}\left(H^{q+1}\right)$ estimate in the $H^1$ norm and $\mathcal{O}\left(H^{q+2}\right)$ in the $L^2$ norm. For this purpose, $h$ is chosen to be very small in order for the error to be mostly dominated by the surrogate part.

The reference macro-mesh size is given by $H_0$, as illustrated in the left of \[fig:unitsquare\_meshes\]. All finer macro-meshes, with associated mesh sizes $H<H_0$, stem from uniformly refining this reference mesh; see middle and right of \[fig:unitsquare\_meshes\]. The fine mesh, with associated mesh size $h\ll H$, is the 13 times uniformly refined reference macro-mesh, i.e., $h = 2^{-13} H_0$. This fine mesh, has about degrees of freedom. The approximation of the stencil functions through least-squares regression, is done on the mesh associated to mesh size $H_{\mathrm{LS}} = 2^{-8} H$. Note that this keeps the number of sampling points in each macro-element constant to . Each linear system is solved by applying geometric multigrid V(2,2) iterations until a relative residual of is obtained.

In \[tab:capitalhconvergence\_scalar\_h1,tab:capitalhconvergence\_scalar\_l2\], the relative $H^1$ and $L^2$ errors for decreasing mesh sizes $H$ are shown. Both tables show the expected convergence rates. In the case of the $L^2$ norm for $q=3$ and $q=4$, the convergence rate deteriorates for small macro-mesh sizes $H$ because the discretization error is dominating the total error.

In order to show the dependence of the least-squares approach on the sampling level, we provide . Here, we show two plots of the relative $L^2$ errors for fixed $q \in \{1,3\}$, $h=2^{-13}H_0$, and varying $H_\mathrm{LS}$. From this figure, one can see that it is crucial to tune the sampling level fine enough in order to achieve optimal $\mathcal{O}\left(H^{q+2}\right)$ convergence in the $L^2$-norm.

The choice of sampling level $m_{\mathrm{LS}}$ or, equivalently, $H_{\mathrm{LS}}$ is very important, since the cost of the polynomial regression grows exponentially with $m_{\mathrm{LS}}$. However, choosing a too large $H_{\mathrm{LS}}$ may violate the discrete $L^2$ projection property required in \[thm:APrioriBoundL2\] in order to obtain an increased order of convergence. Therefore, it is crucial to choose a suitable $H_{\mathrm{LS}}$ for an optimal ratio between the accuracy of the solution and the run time of the stencil function approximation. In each of our experiments, setting $H_{\mathrm{LS}}$ two to four times larger than the fine mesh size $h$ yielded satisfactory results with respect to accuracy and run time.

Furthermore, in \[fig:stencilplot\_var\_q\_H\] we want to illustrate the dependence of the polynomial degree $q$ and the macro-mesh size $H$ within the surrogate approach. For this purpose, we plot the central true and surrogate stencil functions over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ for different pairings of $q$ and $H$. It can be observed, that there is no visible difference of both functions when either the pairing $H = H_0$ and $q = 8$, or the pairing $H = \nicefrac{H_0}{8}$ and $q = 2$ is chosen. Obviously, the quality of $\tilde{{\mathsf{A}}}$ can be improved by either increasing $q$ or decreasing $H$. For smooth coefficients $K$, increasing $q$ is the more efficient option, like in the $hp$-FEM context.

table \[x index = [0]{}, y index=[1]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ; table \[x index = [0]{}, y index=[2]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ; table \[x index = [0]{}, y index=[3]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ; table \[x index = [0]{}, y index=[4]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ; table \[x index = [0]{}, y index=[5]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ; table \[x index = [0]{}, y index=[8]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ; table \[x index = [0]{}, y index=[9]{}, col sep=comma\] [./results/benchmarks/l2errors.csv]{}; ;

plot coordinates [ (1.7, 4.1e-9) (6.3, 4.1e-9)]{}; ; table \[x index = [0]{}, y index=[2]{}, col sep=comma,restrict x to domain=1:6\] [./results/benchmarks/l2errors\_q3.csv]{}; ; table \[x index = [0]{}, y index=[3]{}, col sep=comma,restrict x to domain=1:6\] [./results/benchmarks/l2errors\_q3.csv]{}; ; table \[x index = [0]{}, y index=[4]{}, col sep=comma,restrict x to domain=1:6\] [./results/benchmarks/l2errors\_q3.csv]{}; ; table \[x index = [0]{}, y index=[5]{}, col sep=comma,restrict x to domain=1:6\] [./results/benchmarks/l2errors\_q3.csv]{}; ; table \[x index = [0]{}, y index=[8]{}, col sep=comma\] [./results/benchmarks/l2errors\_q3.csv]{}; ; table \[x index = [0]{}, y index=[9]{}, col sep=comma\] [./results/benchmarks/l2errors\_q3.csv]{}; ;

![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](q2.png "fig:"){width="25.00000%"}![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](q4.png "fig:"){width="25.00000%"}![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](q6.png "fig:"){width="25.00000%"}![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](q8.png "fig:"){width="25.00000%"}\
![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](level_0.png "fig:"){width="25.00000%"}![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](level_1.png "fig:"){width="25.00000%"}![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](level_2.png "fig:"){width="25.00000%"}![\[fig:stencilplot\_var\_q\_H\]Plots of true stencil functions in orange and surrogate stencil functions in blue for $\delta = 0$ over the subdomain $\{(x,y)^\top \in \Omega : x + y \geq 1\}$ in the case of the variable coefficient Poisson equation. Top row: Fixed $H = H_0$ and varying $q = 2$, $4$, $6$, and $8$ from left to right. Bottom row: Fixed $q=2$ and varying $H = H_0$, $\nicefrac{H_0}{2}$, $\nicefrac{H_0}{4}$, and $\nicefrac{H_0}{8}$ from left to right.](level_3.png "fig:"){width="25.00000%"}

### Tensor coefficient on domain with curved boundaries {#sec:tensor_coefficient_benchmark}

In the second benchmark problem, we study problem \[eqn:strongform\] with the symmetric and positive definite tensor coefficient $$\begin{aligned}
\label{eqn:benchmark_tensor_coefficient}
K(x,y) = \left[\begin{matrix}3 x^{2} + 2 y^{2} + 1 & - x^{2} - y^{2}\\- x^{2} - y^{2} & 4 x^{2} + 5 y^{2} + 1\end{matrix}\right]
\,.\end{aligned}$$ Moreover, we consider the domain $\Omega$ with the curved boundary illustrated in \[fig:benchmark\_curved\_domain\]. In the following scenarios, $a = 0.1$ is used as the amplitude of the boundary perturbation. The mapping from the reference unit-square to the perturbed domain is defined by $\varphi$ in \[eq:benchmark\_mapping\_perturbed\]. To map the coefficient onto the perturbed domain, we replace the coefficient $K$ in \[eqn:strongform\] by a new coefficient, $K_0$, induced by the domain transformation, *viz.*, $$\begin{aligned}
\label{eq:benchmark_mapping_perturbed}
K_0 = \frac{D\varphi^{-1} (K\circ \varphi) {\hspace{0.25pt}}D\varphi^{-\top}}{|\det{\left(D\varphi^{-1}\right)}|}
\,,
\qquad\text{where}\qquad
\varphi(x,y) &= \left[\begin{matrix}
x\\
\left(2ay - a\right) \sin^2{\left( 2\pi x\right)} + y
\end{matrix}\right]
\,.\end{aligned}$$ The manufactured solution $u$ is chosen to be $u(x,y) = \sin(\varphi_1(x,y))\sinh(\varphi_2(x,y))$. The restriction of $u$ to the boundary is chosen as Dirichlet datum $g$ and the right-hand-side $f$ is directly computed by inserting $u$ into the strong form of the equation \[eqn:strongform\].

Here, we perform the same verification as in the previous subsection. That is, fixing $h$ and varying $q$ and $H$ with the same meshes and solver settings. In \[tab:capitalhconvergence\_tensor\_h1,tab:capitalhconvergence\_tensor\_l2\], the relative $H^1$ and $L^2$ errors for decreasing mesh sizes $H$ are shown. Both tables show the expected convergence rates. In the case of $q=4$, the convergence rate deteriorates for small macro-mesh sizes $H$, because the discretization error is dominating.

Additionally, we present results for fixed $H = 2^{-3}H_0$ and varying $h$ and $q$. shows the relative $L^2$ errors and convergence rates of the standard approach and the surrogate approach with $q \in \{3,5,7\}$. Only for $q = 7$, the $L^2$ error coincides with the errors from the standard approach for all $h$. The relative time-to-solution (rtts) shown for the surrogate approaches is defined as the time-to-solution (tts) including the setup-phase of the surrogate approach divided by the time-to-solution of the standard approach. In the case with the smallest $h$, the surrogate approach took at most only $7\%$ of the time of the standard approach. That is, a speed-up by more than a factor of $14$.

![\[fig:benchmark\_curved\_domain\]Illustration of the mapping $\varphi$ from the unit-square to the perturbed unit-square. The top boundary is parametrized by $y = a\cdot\sin\left(2\pi x\right)^2 + 1$ and the bottom boundary by $y = -a\cdot\sin\left(2\pi x\right)^2$.](benchmark_curved_boundary_domain.pdf){width="0.6\linewidth"}

Linearized elasticity example
-----------------------------

In this subsection, we compare a standard method with a surrogate method applied to the linearized elasticity problem presented in . In our surrogate method, we employ the zero row sum property described in . We choose an annular domain composed of two distinct and concentric materials under uniform pressure loading. The problem is inspired by a similar 3D experiment documented in [@fuentes2017coupled]. Let $\mathcal{B}_r \operatorname{\subseteq}\mathbb{R}^{2}$ be the two-dimensional open ball of radius $r$ with the midpoint at the origin. The computational domain is then defined as $\Omega = \mathcal{B}_{R_{\mathrm{out}}} \setminus \mathcal{B}_{R_{\mathrm{in}}}$. We split this domain into two disjoint sets $\Omega_I = \mathcal{B}_{R_{\mathrm{mid}}} \setminus \mathcal{B}_{R_{\mathrm{in}}}$ and $\Omega_O = \mathcal{B}_{R_{\mathrm{out}}} \setminus \mathcal{B}_{R_{\mathrm{mid}}}$ corresponding to each material. In our experiments, we fix $R_{\mathrm{in}} = \SI{1}{\centi\meter}$, $R_{\mathrm{mid}} = \SI{1.75}{\centi\meter}$, and $R_{\mathrm{out}} = \SI{2}{\centi\meter}$. Refer to the leftmost diagram in \[fig:linear\_elasticity\_domain\] for an illustration of the setup. Here, the macro-elements adjacent to the boundary and the material interface are mapped to the physical geometry by using the transformation described in [@zlamal1973curved].

![\[fig:linear\_elasticity\_domain\]Linear elasticity problem setup (left) and initial macro-mesh ${{\mathcal{T}}}_H$ (right).](linear_elasticity_annulus.pdf){width="\linewidth"}

![\[fig:linear\_elasticity\_domain\]Linear elasticity problem setup (left) and initial macro-mesh ${{\mathcal{T}}}_H$ (right).](annulus_twomaterial_coarse_mesh.pdf){width="\linewidth"}

Let $r(x) = |x|$. The strong form of the problem is $$\begin{alignedat}{2}
-\mathrm{Div}\left(\sigma\right) &= \vec{f} &&\quad \text{in } \Omega,\\
u_{\theta} &= 0 &&\quad \text{on } \{(-R_{\mathrm{out}},0)^\top, (0,-R_{\mathrm{out}})^\top, (R_{\mathrm{out}},0)^\top\},\\
\sigma \cdot \hat{n} &= p_{\mathrm{in}} \hat{e}_r &&\quad \text{on }\{x \in \partial\Omega : r = R_{\mathrm{in}} \},\\
\sigma \cdot \hat{n} &= -p_{\mathrm{out}} \hat{e}_r &&\quad \text{on } \{x \in \partial\Omega : r = R_{\mathrm{out}} \}.
\end{alignedat}
\label{eqn:linearelasticity}$$ Here, the stress tensor $\sigma$ is given by Hooke’s law for isotropic materials as defined in . The unit vector in radial direction is denoted $\hat{e}_r$ and the outward pointing unit normal vector is denoted $\hat{n}$. We neglect body forces and therefore set $\vec{f} = \vec{0}$. The displacement is described in polar coordinates where $u_r$ is the radial displacement and $u_\theta$ is the tangential displacement. In order to make the system uniquely solvable, we enforce the tangential displacement $u_{\theta}$ to be zero at three points; see \[fig:linear\_elasticity\_domain\]. The materials are chosen to be cork in the inner domain $\Omega_I$ with an A36 steel layer in the outer domain $\Omega_O$. Poisson’s ratio and Young’s modulus for this scenario are $E_I = \SI{0.02}{\giga\pascal}$, $E_O = \SI{200.0}{\giga\pascal}$, $\nu_I = 0$, and $\nu_O = 0.26$. The Lamé parameters are obtained by the expressions $\mu = \frac{E}{2(1+\nu)}$ and $\lambda = \frac{E\nu}{(1-2\nu)\cdot(1+\nu)}$. Note that while these expressions induce piecewise-constant Lamé parameters, $\lambda = \lambda(r)$ and $\mu = \mu(r)$, once the smooth domain is mapped to the computational domain depicted on the right of \[fig:linear\_elasticity\_domain\], these parameters will not be piecewise-constant anymore due to the transformation. The pressure on the outer boundary is set to $p_{\mathrm{out}} = \SI{0}{\mega\pascal}$ and on the inner boundary $p_{\mathrm{in}} = \SI{1}{\mega\pascal}$ is prescribed.

In this particular scenario, there is an analytic solution available for the radial displacement $u_r$ which has the form $$\label{eqn:linear_elasticity_analytic_solution}
\begin{aligned}
u_r(r) = \begin{cases} 
A \cdot r + B \cdot r^{-1} &\quad \text{if }r \in [R_{\mathrm{in}}, R_{\mathrm{mid}}]\,, \\
C \cdot r + D \cdot r^{-1} &\quad \text{if }r \in (R_{\mathrm{mid}}, R_{\mathrm{out}}]\,.
\end{cases}
\end{aligned}$$ The tangential displacement $u_\theta$ is zero everywhere due to the symmetry of the problem. In \[eqn:linear\_elasticity\_analytic\_solution\], the coefficients $A$, $B$, $C$, and $D$ are uniquely determined by the following system of linear equations: $$\label{eqn:linear_elasticity_analytic_system}
\begin{aligned}
\begin{bmatrix}
E_{I}R_{\mathrm{in}}^{2} & E_{I} \left(2\nu_{I} - 1\right) & 0 & 0\\
0 & 0 & E_{O}R_{\mathrm{out}}^{2} & E_{O} \left(2\nu_{O}-1\right)\\
R_{\mathrm{mid}}^2 & 1 & -R_{\mathrm{mid}}^2 & -1\\
- E_{I} R_{\mathrm{mid}}^{2} d_{O} & -d_{O}  E_{I} \left(2\nu_{I} - 1\right) & E_{O} R_{\mathrm{mid}}^{2} d_{I} & d_{I} E_{O} \left(2\nu_{O} - 1\right)
\end{bmatrix} \begin{bmatrix}
A\\
B\\
C\\
D
\end{bmatrix} &= \begin{bmatrix}
p_{in} R_{\mathrm{in}}^{2} d_{I}\\
p_{out} R_{\mathrm{out}}^{2} d_{O}\\
0\\
0
\end{bmatrix},
\end{aligned}$$ where $d_{I} \coloneqq 2 \nu_{I}^{2} + \nu_{I} - 1$ and $d_{O} \coloneqq 2 \nu_{O}^{2} + \nu_{O} - 1$. This system is derived after deducing the continuity of the displacement and surface traction at the material interface, then by incorporating the prescribed external forces at the boundaries.

![\[fig:linear\_elasticity\_results\]Plots of radial displacement $u_r$ (left) and the tangential stress $\sigma_{\theta\theta}$ (right) computed on the fine mesh ${\mathcal{S}}^6({{\mathcal{T}}}_H)$, corresponding to $h = 2^{-6}H$, with $q = 4$.](displacement_cropped.png){width="\linewidth"}

![\[fig:linear\_elasticity\_results\]Plots of radial displacement $u_r$ (left) and the tangential stress $\sigma_{\theta\theta}$ (right) computed on the fine mesh ${\mathcal{S}}^6({{\mathcal{T}}}_H)$, corresponding to $h = 2^{-6}H$, with $q = 4$.](sigma_cropped.png){width="\linewidth"}

table \[x index = [0]{}, y index=[1]{}, col sep=comma\] [./results/linearelasticity/plot.csv]{}; ; table \[x index = [0]{}, y index=[2]{}, col sep=comma\] [./results/linearelasticity/plot.csv]{}; ;

table \[x index = [0]{}, y index=[3]{}, col sep=comma\] [./results/linearelasticity/plot.csv]{}; table \[x index = [0]{}, y index=[3]{}, col sep=comma\] [./results/linearelasticity/plot.csv]{}; ; table \[x index = [0]{}, y index=[4]{}, col sep=comma, restrict x to domain=1:1.75\] [./results/linearelasticity/plot.csv]{}; table \[x index = [0]{}, y index=[4]{}, col sep=comma, restrict x to domain=1.75:2, forget plot\] [./results/linearelasticity/plot.csv]{}; ;

In order to verify the accuracy of the surrogate method, we select the polynomial degree $q=7$ and the macro-mesh ${{\mathcal{T}}}_H$, illustrated on the right of \[fig:linear\_elasticity\_domain\]. Note that the discontinuity in the material parameters lies along the macro-element interfaces and so $\lambda,\mu \in \prod_{T\in{{\mathcal{T}}}_H}C^\infty(T) \subsetneq W^{r+1,\infty}({{\mathcal{T}}}_H)$, for any $r>0$.

Each linear system is solved by applying geometric multigrid iterations with V(3,3) cycles until the relative residual is reduced by the factor . On the coarsest level used in the multigrid hierarchy, we employ MUMPS as a direct solver. In \[tab:littlehconvergence\_linearelasticity\], we report on the results for varying $h$ and present the relative $L^2$ errors for the standard and surrogate approach, respectively. On the finest mesh involving about degrees of freedom, the surrogate approach required only 5% of the time required by the standard approach while having the same accuracy. That is a speed up by a factor of $20$. \[fig:linear\_elasticity\_results\] shows the radial displacement $u_r$ and the tangential stress $\sigma_{\theta\theta}$ computed with the surrogate approach on the fine mesh ${\mathcal{S}}^6({{\mathcal{T}}}_H)$, corresponding to $h = 2^{-6}H$, with $q = 4$. This is illustrated further by the plots in \[fig:linear\_elasticity\_lineplots\] which allow a visual comparison between $u_r$ and $\sigma_{\theta\theta}$ in the surrogate and analytical solutions.

[r|r|r|c|c|c|c|c]{} & & & &\
& & & rel. $L^2$ err. & eoc & rel. $L^2$ err. & eoc & rtts\

$p$-Laplacian diffusion example {#sec:plaplacian_diffusion_example}
-------------------------------

In this subsection, we consider the time-dependent example introduced in . Here, we solve the non-linear $p$-Laplacian diffusion problem \[eqn:nonlinearheatstrong\], given in strong form as $$\begin{alignedat}{2}
\frac{\partial u}{\partial t} - \mathrm{div}\left(|\nabla u|^{p-2} \cdot u\right) &= f &&\quad \text{in } \Omega \times (0,T]\,,\\
u &= 0 &&\quad \text{on } \partial \Omega \times (0,T]\,,\\
u &= u_0 &&\quad \text{in } \Omega \times \{0\}\,.
\end{alignedat}
\label{eqn:nonlinearheatstrong}$$ The computational domain is set to the unit disk, i.e., $\Omega \coloneqq \mathcal{B}_1$ and the right-hand-side is set to a specific constant, $f(x) = 2\hat{q}^{\nicefrac{p}{\hat{q}}}$, where $\hat{q} = \frac{p}{p-1}$. The initial solution is set to $u_0(x) = 0.1 \cdot \left(1 - |x|^2\right)$. For this particular problem, the stationary limit $u_\infty$ has an analytic solution exists with unit magnitude, namely $u_\infty(x) = 1 - |x|^{\hat{q}}$ [@barrett1993finite Example 3.1].

Our discretization follows a standard approach where a mass matrix ${\mathsf{M}}_{ij} = \int_\Omega \phi_i \phi_j{\hspace{0.25pt}}{\,\mathrm{d}}x$ and a stiffness matrix ${\mathsf{A}}_{ij}(\tilde{{\mathsf{u}}}) = \int_\Omega |\nabla \tilde{u}_h|^{p-2}{\hspace{0.25pt}}\nabla\phi_j\cdot\nabla \phi_i{\hspace{0.25pt}}{\,\mathrm{d}}x$ are introduced. At this point, $\tilde{{\mathsf{u}}}$ is the coefficient vector, in the $\{\phi_i\}$ basis, of an arbitrary discrete function $\tilde{u}_h$. The time derivative is discretized by a backward Euler scheme, and the non-linearity in each time step is resolved by Picard fixed-point iterations. Let ${\mathsf{u}}^{l}_k$ be the coefficient vector of the discrete solution at the $k$-th time step and $l$-th fixed point iteration. Employing the bilinear form \[eq:BilinearFormpLaplace\] and fixing a time step size $dt>0$, the discrete problem in each time step $k>0$ and fixed point iteration $l>0$ reads as follows: $$\Big({\mathsf{M}}+ dt{\hspace{0.25pt}}{\mathsf{A}}{\hspace{-0.25pt}}\big({\mathsf{u}}^{l-1}_k\big)\Big){\hspace{0.25pt}}{\mathsf{u}}^{l}_k = {\mathsf{M}}{\mathsf{u}}_{k-1} + dt{\hspace{0.25pt}}{\mathsf{M}}{\mathsf{f}}\,,
\label{eqn:nonlinearheatdiscrete}$$ where ${\mathsf{u}}_{k-1}$ is the final coefficient vector from the previous time step. In each time step, this system is solved multiple times (once for each fixed-point iteration) by the application of five V(2,2) multigrid cycles. The fixed-point iterations continue until the relative increment $\frac{\|{\mathsf{u}}^l_k - {\mathsf{u}}^{l-1}_k\|_2}{\|{\mathsf{u}}^l_k\|_2}$ is smaller than the fixed tolerance . Then $k$ is incremented and a new ${\mathsf{u}}_{k-1} = {\mathsf{u}}_{k-1}^l$ is defined.

In our surrogate method, the stencil functions of the stiffness matrix ${\mathsf{A}}{\hspace{-0.25pt}}\big({\mathsf{u}}^{k-1}\big)$ are approximated by solving the least-squares problems after every fixed-point iteration, all the while enforcing the zero row sum property (cf. ). Meanwhile, the stencil function of the mass matrix ${\mathsf{M}}$ is only approximated once in a pre-processing step because it does not depend on any free variables in the computation. The time step surrogate polynomials of both operators are then simply summed together to obtain the time step matrix ${\mathsf{M}}+ dt{\hspace{0.25pt}}{\mathsf{A}}{\hspace{-0.25pt}}\big({\mathsf{u}}^{k-1}\big)$. This particular splitting of the surrogate matrices, which reproduces the zero row sum property in the stiffness matrix, allows for faster re-approximation of the time step matrix stencil function and appears to improve the stability of the method. In this example, we did not enforce the symmetry condition featured in . Instead, whenever a vertex $x_i\in{\mathbb{X}}_m$ was on the boundary of a macro-element ${\partial}T_m$, we set the surrogate stiffness matrix to the exact value stiffness matrix $\tilde{{\mathsf{A}}}_{ij}= {\mathsf{A}}_{ij}$. This minor asymmetry, is more amenable to computation because there is less data transfer and it did not affect the behavior of our multigrid solver. In fact, this choice improved our results with this problem, which we believe is due to better accuracy in the surrogate near the singularity in the coefficient; i.e., at the origin $x=(0,0)$. The success of this approach suggests that the definition given in \[eq:SurrogateStiffnessMatrix\] may be relaxed in other applications as well.[^2] In the proximity of this singularity and for $p > 2$, the coefficient depending on the solution of the previous fixed-point iteration is getting very close to zero which serves as a challenge for the approximated off-diagonal stencil functions. Depending on the polynomial degree $q$, they might erroneously take on positive values due to overshoots which possibly results in a loss of positive definiteness of the surrogate matrix. However, this drawback could not be observed in the scenario considered in the following example.

table \[x index = [4]{}, y index=[1]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_surrogate.0.csv]{}; table \[x index = [4]{}, y index=[0]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_surrogate.0.csv]{}; table \[x index = [4]{}, y index=[0]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_standard.0.csv]{};

table \[x index = [4]{}, y index=[1]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_surrogate.10.csv]{}; table \[x index = [4]{}, y index=[0]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_surrogate.10.csv]{}; table \[x index = [4]{}, y index=[0]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_standard.10.csv]{};

table \[x index = [4]{}, y index=[1]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_surrogate.100.csv]{}; ; table \[x index = [4]{}, y index=[0]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_surrogate.100.csv]{}; ; table \[x index = [4]{}, y index=[0]{}, col sep=comma\] [./results/nonlinearheatequation/pthree\_standard.100.csv]{}; ;

![\[fig:nonlinear\_surfaceplots\]Surface plot of non-stationary $p$-Laplacian surrogate solution for $t = T_{\mathrm{end}}$ with $p=3$ and $q=6$ (left). Absolute difference between the discrete standard and surrogate solution at the time $t = T_{\mathrm{end}}$ (right).](surfaceplot_p3_u_cropped.png){width="\textwidth"}

![\[fig:nonlinear\_surfaceplots\]Surface plot of non-stationary $p$-Laplacian surrogate solution for $t = T_{\mathrm{end}}$ with $p=3$ and $q=6$ (left). Absolute difference between the discrete standard and surrogate solution at the time $t = T_{\mathrm{end}}$ (right).](surfaceplot_p3_diff_cropped.png){width="\textwidth"}

The unit-disk is discretized by the macro-mesh ${{\mathcal{T}}}_H$ featured on the right of \[fig:stencilweight\_surfaceplot\]. Note that the vertices of the central macro-elements meet at the origin $x=(0,0)$; i.e., exactly where the singularity occurs in the stationary limit $u_\infty$. The simulations are conducted on the mesh ${\mathcal{S}}^9({{\mathcal{T}}}_H)$, which involves about degrees of freedom. The macro-elements adjacent to the boundary are mapped to the physical geometry by using the mapping described in [@zlamal1973curved]. Furthermore, the least-squares regressions are carried out on the mesh corresponding to $H_{\mathrm{LS}} = 4h$ and the polynomial degree of the approximated stencil functions is fixed to $q = 6$. In this scenario, we consider the $p$-Laplacian operator with $p = 3$, fix the time step size $\Delta t = \num{1e-2}$, and solve until time $T_{\mathrm{end}} = 1$. illustrates the standard and surrogate solutions plotted over the line $[0,1] \times \{0\}$ for different times $t$. In the left of \[fig:nonlinear\_surfaceplots\], the surface plot of the surrogate solution is depicted. Since the difference of the solutions is very small, we added in the right of \[fig:nonlinear\_surfaceplots\] a surface plot of the absolute difference of the surrogate and standard solution at the final time $t = T_{\mathrm{end}}$. The surrogate approach required only $5.4$% of the time required by the standard approach. That is, a speed-up by more than a factor of $18$.

Proofs {#app:appendix}
======

By the min-max theorem [@courant1953methods], the $k$-th eigenvalue of ${\mathsf{M}}$ is $$\begin{aligned}
    \lambda_k({\mathsf{M}})
    &=
    \operatorname*{\vphantom{p}min}_{W\operatorname{\subseteq}{\mathbb{R}}^N}
    \Bigg\{
        \operatorname*{\vphantom{p}max}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{M}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        :
        \mathrm{dim}\,W=k
    \Bigg\}
    \\
    &=
    \operatorname*{\vphantom{p}max}_{W\operatorname{\subseteq}{\mathbb{R}}^N}
    \Bigg\{
        \operatorname*{\vphantom{p}min}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{M}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        :
        \mathrm{dim}\,W=N-k+1
    \Bigg\}
    \,.\end{aligned}$$

Define ${\mathsf{D}}= {\mathsf{M}}- {\mathsf{N}}$. We first show that $\lambda_1({\mathsf{D}}) \leq  \lambda_k({\mathsf{M}}) - \lambda_k({\mathsf{N}}) \leq \lambda_N({\mathsf{D}})$. Indeed, $$\begin{aligned}
    \lambda_k({\mathsf{M}})
                                                &\leq
    \operatorname*{\vphantom{p}min}_{W\operatorname{\subseteq}{\mathbb{R}}^N}
    \Bigg\{
        \operatorname*{\vphantom{p}max}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{N}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        +
        \operatorname*{\vphantom{p}max}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{D}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        :
        \mathrm{dim}\,W=k
    \Bigg\}
    \\
    &\leq
    \operatorname*{\vphantom{p}min}_{W\operatorname{\subseteq}{\mathbb{R}}^N}
    \Bigg\{
        \operatorname*{\vphantom{p}max}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{N}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        :
        \mathrm{dim}\,W=k
    \Bigg\}
    +
    \operatorname*{\vphantom{p}max}_{\|{\mathsf{x}}\|_2 = 1}
    \Big\{
        {\mathsf{x}}^\top {\mathsf{D}}{\mathsf{x}}: {\mathsf{x}}\in {\mathbb{R}}^N
    \Big\}
    \\
    &=
    \lambda_k({\mathsf{N}}) + \lambda_N({\mathsf{D}})\end{aligned}$$ and, likewise, $$\begin{aligned}
    \lambda_k({\mathsf{M}})
                                                &\geq
    \operatorname*{\vphantom{p}max}_{W\operatorname{\subseteq}{\mathbb{R}}^N}
    \Bigg\{
        \operatorname*{\vphantom{p}min}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{N}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
                                                +
        \operatorname*{\vphantom{p}min}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{D}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        :
        \mathrm{dim}\,W=N-k+1
    \Bigg\}
    \\
    &\geq
    \operatorname*{\vphantom{p}max}_{W\operatorname{\subseteq}{\mathbb{R}}^N}
    \Bigg\{
        \operatorname*{\vphantom{p}min}_{\|{\mathsf{x}}\|_2 = 1}
        \Big\{
            {\mathsf{x}}^\top {\mathsf{N}}{\mathsf{x}}: {\mathsf{x}}\in W
        \Big\}
        :
        \mathrm{dim}\,W=N-k+1
    \Bigg\}
                        +
    \operatorname*{\vphantom{p}min}_{\|{\mathsf{x}}\|_2 = 1}
    \Big\{
        {\mathsf{x}}^\top {\mathsf{D}}{\mathsf{x}}: {\mathsf{x}}\in {\mathbb{R}}^N
    \Big\}
    \\
    &=
    \lambda_k({\mathsf{N}}) + \lambda_1({\mathsf{D}})
        \,.\end{aligned}$$ This immediately leads us to the inequality $|\lambda_k({\mathsf{M}}) - \lambda_k({\mathsf{N}})| \leq \operatorname*{\vphantom{p}max}\{|\lambda_1({\mathsf{D}})|,|\lambda_N({\mathsf{D}})|\}$. Now, for at least one $i$, $|\lambda_N({\mathsf{D}})| - |{\mathsf{D}}_{ii}| \leq |\lambda_N({\mathsf{D}}) - {\mathsf{D}}_{ii}| \leq \sum_{j\neq i} |{\mathsf{D}}_{ij}|$, by the Gershgorin circle theorem. Therefore, $|\lambda_N({\mathsf{D}})| \leq \sum_{j} |{\mathsf{D}}_{ij}| \leq \|{\mathsf{D}}\|_{\infty}$. Similarly, $|\lambda_1({\mathsf{D}})| \leq \|{\mathsf{D}}\|_{\infty}$.

Acknowledgments {#acknowledgments .unnumbered}
===============

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 800898. This work was also partly supported by the German Research Foundation through the Priority Programme 1648 “Software for Exascale Computing” (SPPEXA) and by grant WO671/11-1. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (GCS, [www.gauss-centre.eu](www.gauss-centre.eu)) for funding this project by providing computing time on the GCS supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, [www.lrz.de](www.lrz.de)).

[^1]: Lehrstuhl für Numerische Mathematik, Fakultät für Mathematik (M2), Technische Universität München, Garching bei München (, , )

[^2]: See [@bauer2017two] for further evidence.