<figure><img src="image_1.png" /></figure>

HAL
open science

# Error Estimates for Sparse Tensor Products of B-spline Approximation Spaces

Clément Guillet

► To cite this version:

Clément Guillet. Error Estimates for Sparse Tensor Products of B-spline Approximation Spaces. RR-9598, Centre Inria de l'Université de Bordeaux. 2025. hal-05326739

HAL Id: hal-05326739

https://hal.science/hal-05326739v1

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inria

Error Estimates for
Sparse Tensor Products
of B-spline
Approximation Spaces

Clément Guillet

RESEARCH
REPORT

Nº 9598

October 2025

Project-Teams Concace
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

# Error Estimates for Sparse Tensor Products of B-spline Approximation Spaces

Clément Guillet*

Project-Teams Concace

Research Report n° 9598 — October 2025 — 27 pages

**Abstract:** This work introduces and analyzes B-spline approximation spaces defined on general geometric domains obtained through a mapping from a parameter domain. These spaces are constructed as sparse-grid tensor products of univariate spaces in the parameter domain and are mapped to the physical domain via a geometric parametrization. Both the univariate approximation spaces and the geometric mapping are built using maximally smooth B-splines. We construct two such spaces, employing either the sparse-grid combination technique or the hierarchical subspace decomposition of sparse-grid tensor products, and we prove their mathematical equivalence. Furthermore, we derive approximation error estimates and inverse inequalities that highlight the advantages of sparse-grid tensor products. Specifically, under suitable regularity assumptions on the solution, these spaces achieve the same approximation order as standard tensor product spaces while using significantly fewer degrees of freedom. Additionally, our estimates indicate that, in the case of non-tensor-product domains, stronger regularity assumptions on the solution—particularly concerning isotropic (non-mixed) derivatives—are required to achieve optimal convergence rates compared to sparse-grid methods defined on tensor-product domains.

**Key-words:** sparse grids, combination technique, approximation spaces, error estimates, interpolation, inverse inequality, B-splines, maximal smoothness, isogeometric analysis

<footnote>* Inria, Concace joint project Airbus CR & T, Cerfacs and Inria, France</footnote>
# Estimations d'Erreur pour les Espaces d'Approximation *Sparse-grid* par B-splines

**Résumé** : Ce travail introduit et analyse des espaces d'approximation par B-splines définis sur des domaines géométriques généraux obtenus par transformation géométrique depuis un domaine paramétrique. Ces espaces sont construits comme produits tensoriels *sparse grids* d'espaces unidimensionnels définis dans le domaine paramétrique, puis transférés vers le domaine physique au moyen d'une paramétrisation géométrique. Tant les espaces d'approximation unidimensionnels que la transformation géométrique sont construits à partir de B-splines de régularité maximale. Nous construisons deux de ces espaces, en utilisant soit la technique de combinaison *sparse grid*, soit la décomposition hiérarchique en sous-espaces des produits tensoriels *sparse grid*, et nous en démontrons l'équivalence mathématique. Nous établissons des estimations d'erreur d'approximation et des inégalités inverses qui mettent en évidence les avantages des produits tensoriels *sparse grid*. Plus précisément, sous des hypothèses de régularité appropriées sur la solution, ces espaces atteignent le même ordre d'approximation que les espaces à produit tensoriel complet, tout en nécessitant un nombre de degrés de liberté bien inférieur. Enfin, nos estimations montrent que, dans le cas de domaines non tensoriels, des hypothèses de régularité plus fortes sur la solution — notamment en ce qui concerne les dérivées isotropes (non mixtes) — sont nécessaires pour obtenir des taux de convergence optimaux, comparativement aux méthodes *sparse grids* définies sur des domaines tensoriels.

**Mots-clés** : sparse grids, technique de combinaison, espaces d'approximation, estimations d'erreur, interpolation, inégalité inverse, B-splines, régularité maximale, analyse isogéométrique
# Contents

<table><tr><td>1 Introduction</td><td>3</td></tr><tr><td>2 Preliminaries</td><td>5</td></tr><tr><td>2.1 Sobolev Spaces and Norms</td><td>5</td></tr><tr><td>2.2 Parameter Space Meshes</td><td>6</td></tr><tr><td>3 Standard Approximation Spaces</td><td>7</td></tr><tr><td>3.1 Univariate B-splines Spaces</td><td>7</td></tr><tr><td>3.2 Anisotropic d-dimensional B-spline Spaces</td><td>8</td></tr><tr><td>4 Sparse-Grid Approximation Spaces</td><td>12</td></tr><tr><td>4.1 Approximation Spaces in Tensor-Product Domains</td><td>12</td></tr><tr><td>4.1.1 Hierarchical Subspace Decomposition</td><td>13</td></tr><tr><td>4.1.2 Combination Technique</td><td>14</td></tr><tr><td>4.1.3 Equivalence Between Spaces</td><td>14</td></tr><tr><td>4.1.4 Approximation Error Estimate</td><td>16</td></tr><tr><td>4.1.5 Inverse Inequality</td><td>19</td></tr><tr><td>4.2 Extension to Non-Tensor-Product Domains</td><td>21</td></tr><tr><td>4.2.1 Approximation Error</td><td>22</td></tr><tr><td>4.2.2 Inverse inequality</td><td>24</td></tr><tr><td>4.2.3 Complexity of the Spaces</td><td>24</td></tr></table>

# 1 Introduction

Many problems in computational science involve the approximation of high-dimensional functions. However, representing such functions on isotropic, uniformly refined meshes quickly becomes computationally intractable as the dimensionality increases beyond three or four. This phenomenon, known as the curse of dimensionality [6], refers to the exponential growth of computational cost with respect to the number of dimensions. Sparse tensor-product constructions, commonly referred to as sparse-grid (SG) methods, were introduced primarily to alleviate this curse by reducing the computational complexity while maintaining high approximation accuracy.

Sparse-grid methods [9, 12] were originally developed for the interpolation of high-dimensional functions and have since been extended to the approximation of partial differential equations [13, 14, 15]. Under the assumption that the solution has bounded mixed derivatives, an accurate approximation, achieving nearly the same order of accuracy as standard tensor-product methods— using isotropic uniform meshes— can be computed at a reduced cost. Convergence rates of sparse-grid approximations using piecewise linear functions are well established [9, 22]. Specifically, the $L^2$ error of sparse-grid approximations scales as $O(h^2|\log h|^{d-1})$, where $h$ denotes the mesh size and $d$ the dimension. Constructing the sparse-grid approximation requires only $O(h^{-1}|\log h|^{d-1})$ basis functions whereas standard tensor-product approximations require $O(h^{-d})$ basis functions to achieve an $L^2$ error of $O(h^2)$. These results were extended to polynomials [25] and piecewise Lagrange polynomials [1, 7, 8, 9] with an $L^2$ error scaling as $O(h^{p+1}|\log h|^{d-1})$, for polynomials of degree $p$. Similar results have also been proven for tensor-product operators [18, 16] and for non-tensor-product operators [17].

Sparse-grid methods were initially restricted to generic tensor-product domains. Some early works proposed to use of transfinite interpolation to generate curvilinear domains [8, 11]. However, the regularity assumptions on the solution may no longer hold if the geometric mapping is insufficiently smooth, as is the case with transfinite interpolation. More recently, to overcome this limitation, sparse grids were
introduced in the context of isogeometric analysis (IGA) [3], which can be viewed as a natural extension of finite element analysis using high-order, highly regular B-spline or non-uniform rational B-spline (NURBS) basis functions. The main idea of IGA, initially introduced in [20] is to employ the same type of functions—B-splines or NURBS—to define both the computational domain and the unknown solution of the differential equations. This approach allows for more accurate representation of complex geometries and enables exact descriptions of common geometries such as circles or ellipses. In [3], the authors introduce a sparse-grid isogeometric solver based on the so-called sparse-grid combination technique[19], which consists of forming a linear combination of several inexpensive tensor-product approximations to recover an approximation of nearly the same order as that of full tensor-product discretizations. Numerical experiments demonstrated that the expected convergence rate of $O(h^{p+1}|\log h|^{d-1})$ is achieved for sparse-grid approximations using NURBS basis functions of degree $p$.

However, rigorous approximation error estimates for sparse-grid spaces on general, non-tensor-product domains are still lacking. Understanding such approximation properties is of fundamental importance in the context of Galerkin methods applied to the variational formulations of differential equations. Indeed, error estimates in appropriate Sobolev norms guarantee the convergence of the methods through the application of Céa's or Strang's lemma. The first error estimates for IGA approximation spaces were derived in the case of h-refinement with fixed polynomial degree [2]. The case of anisotropic meshes, involving different mesh sizes in each coordinate direction, was addressed in [5]. This analysis was later extended to include the explicit dependence on the polynomial degree $p$ [4], although it did not cover the most relevant case of maximally smooth B-splines. This gap was filled in [24], where sharp error estimates and inverse inequalities with constants independent of the spline degree and the mesh size were established for maximally smooth B-splines.

In this work, we introduce and analyze approximation spaces based on sparse-grid tensor products of B-splines. This extends the work of [24], where such estimates were established for univariate B-splines of maximal smoothness and for $d$-dimensional spaces constructed from full tensor products. Here, we generalize these results to sparse-grid tensor-product spaces in the parameter domain and their images in the physical domain via a geometric B-spline mapping. For simplicity, we restrict our analysis to B-spline functions; however, the extension to NURBS functions, as in [2], is straightforward. Following [24], we restrict our analysis to single-patch geometric domains; however, an extension to multi-patch domains, along the lines of [23], could also be considered.

We propose two distinct approaches to construct such spaces. The first approach is based on the sparse-grid combination technique, while the second relies on the hierarchical subspace decomposition of sparse-grid tensor products. Here, we focus exclusively on the approximation properties of these spaces, leaving implementation and performance considerations to future work.

We further show that, despite their different constructions, the two approximation spaces are identical and can thus be considered a single, unified space.

## The Main Results

We summarize in the following the two main results of this paper, Theorem 4.2 and Theorem 4.3. These results provide interpolation error estimates and inverse inequalities in Sobolev norms for the sparse-grid approximation spaces, with constants independent of the mesh size, the spline degree, and the B-spline regularity.

We consider the physical domain constructed by the application of a geometric mapping $\mathbf{F}$ to a tensor-product parameter domain $\hat{\Omega} := (0, 1)^d$, i.e.,

$$ \Omega := \mathbf{F}(\hat{\Omega}). \qquad (1) $$

Let $n \in \mathbb{N}$ be the maximum mesh level in all directions, associated with the dyadic mesh size $h_n = 2^{-n}$. We denote by $S_{p,h_n}^{\mathcal{H}}(\Omega)$ the sparse-grid approximation space.
**Theorem.** Let $u \in H^{dq}(\Omega)$ with $q \in \mathbb{N}_0$. For any mesh level $n \in \mathbb{N}$, spline degree $p \in \mathbb{N}_0$, and any $r \in \mathbb{N}_0$ satisfying $0 \le r \le q \le p + 1$, there exists a spline approximation $u_{p,h_n}^{\mathcal{H}} \in S_{p,h_n}^{\mathcal{H}}(\Omega)$ such that:

$$|u - u_{p,h_n}^{\mathcal{H}}|_{H^r(\Omega)} \lesssim h_n^{q-r} |\log h_n|^{d-1} \|u\|_{H^{dq}(\Omega)}.$$

**Theorem.** For any mesh level $n \in \mathbb{N}$, spline degree $p \in \mathbb{N}_0$ and any $q \in \mathbb{N}_0$ satisfying $0 \le q \le p$, the inverse inequality:

$$|u_{p,h_n}^{\mathcal{H}}|_{H^q(\Omega)} \lesssim h_n^{-q} |\log h_n|^{d/2} \|u_{p,h_n}^{\mathcal{H}}\|_{L^2(\Omega)}$$

holds true for all spline approximations $u_{p,h_n}^{\mathcal{H}} \in \tilde{S}_{p,h_n}^{\mathcal{H}}(\Omega)$.

The symbol $\lesssim$ refers to an inequality up to a constant that does not depend on the mesh size $h_n$ or the spline degree $p$, but may depend on the shape and the size of the domain $\Omega$.

This paper is organized as follows. In Section 2, we introduce the main notation related to meshes, Sobolev spaces, and the associated norms. In Section 3, we describe the construction of standard approximation spaces based on univariate B-splines and their extension to anisotropic $d$-dimensional splines via tensor products. In Section 4, we present the sparse-grid approximation spaces, first in tensor-product domains and then extended to general geometric domains. The main theoretical results are stated and proved therein.

# 2 Preliminaries

In this preliminary section, we introduce the main notations for spaces, norms, and meshes that will be used in our analysis.

## 2.1 Sobolev Spaces and Norms

Let $d \in \mathbb{N}$ denote the dimension of the physical domain of interest, which we take as a generic open set $\Omega \subset \mathbb{R}^d$. We denote by $W^{q,p}(\Omega)$ the standard $L^p$ Sobolev space of order $q \in \mathbb{N}_0$ with $p \in [1, +\infty]$, and by $H^q(\Omega)$ when $p = 2$. Specifically, these spaces are defined as

$$W^{q,p}(\Omega) := \{u : \Omega \to \mathbb{R} \mid D^{\alpha} u \in L^p(\Omega), \ |\alpha|_1 \le q\},$$

$$H^q(\Omega) := W^{q,2}(\Omega).$$

In the above definitions, the derivative operator, defined for functions $u \in W^{q,p}(\Omega)$ and multi-indices $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{N}^d$, is given by

$$D^{\alpha} u := \frac{\partial^{\alpha_1}}{\partial x_1^{\alpha_1}} \cdots \frac{\partial^{\alpha_d}}{\partial x_d^{\alpha_d}} u, \quad \text{and} \quad D^q u := D^{\alpha} u \quad \text{for } |\alpha|_1 = q.$$

The $\ell^1$- and $\ell^\infty$-norms are defined by

$$|\alpha|_1 := \sum_{i=1}^{d} |\alpha_i|, \quad |\alpha|_{\infty} := \max_{i=1, \dots, d} |\alpha_i|.$$

For a smooth vector function $F = (F_1, \dots, F_d) \in (\mathcal{C}^\infty(\Omega))^d$, we denote by $\nabla^q F$ the $q$-th order multi-derivative tensor of $F$, defined componentwise as

$$(\nabla^q F)_i := \nabla^\alpha F_i, \quad i = 1, \dots, d, \quad |\alpha|_1 = q.$$
In this work, we focus on the L²-based Sobolev space $H^q(\Omega)$. We define the inner product on $H^q(\Omega)$ for $u, v \in H^q(\Omega)$ by

$$ (u, v)_{H^q(\Omega)} := \sum_{|\boldsymbol{\alpha}|_1 \le q} \int_{\Omega} D^{\boldsymbol{\alpha}} u(\mathbf{x}) D^{\boldsymbol{\alpha}} v(\mathbf{x}) d\mathbf{x}. $$

The associated seminorm and norm are defined, respectively, by

$$ |u|_{H^q(\Omega)} := \left( \sum_{|\boldsymbol{\alpha}|_1=q} \|D^{\boldsymbol{\alpha}} u\|_{L^2(\Omega)}^2 \right)^{1/2}, $$

$$ \|u\|_{H^q(\Omega)} := (u, u)_{H^q(\Omega)}^{1/2} = \left( \sum_{m=0}^{q} |u|_{H^m(\Omega)}^2 \right)^{1/2}. $$

Sparse-grid methods require stronger assumptions on the mixed regularity of the solutions than standard tensor-product methods. To this end, we consider function spaces suitable for their analysis. We define the L²-based Sobolev space with dominating mixed smoothness of order $q \in \mathbb{N}_0$, also called the mixed Sobolev space, by

$$ H_{\text{mix}}^{q}(\Omega) := \{u : \Omega \to \mathbb{R} \mid D^{\boldsymbol{\alpha}} u \in L^2(\Omega), \ |\boldsymbol{\alpha}|_{\infty} \le q\}. $$

The associated seminorm and norm are defined by

$$ |u|_{H_{\text{mix}}^q(\Omega)} := \left( \sum_{|\boldsymbol{\alpha}|_{\infty}=q} \|D^{\boldsymbol{\alpha}} u\|_{L^2(\Omega)}^2 \right)^{1/2}, $$

$$ \|u\|_{H_{\text{mix}}^q(\Omega)} := \left( \sum_{m=0}^{q} |u|_{H_{\text{mix}}^m(\Omega)}^2 \right)^{1/2}. $$

The standard Sobolev space is a subspace of the mixed Sobolev space, and the following relation holds between the full norms:

$$ \|u\|_{H^q(\Omega)} \le \|u\|_{H_{\text{mix}}^q(\Omega)} \le \|u\|_{H^{dq}(\Omega)}. \qquad (2) $$

## 2.2 Parameter Space Meshes

We consider a family of meshes $\{\hat{\Omega}_h\}_h$ on a tensor-product parameter domain $\hat{\Omega} := (0, 1)^d$, where $h$ denotes the family index and represents the global mesh size. In this work, we consider a dyadic subdivision of the domain. Let $\ell = (\ell_1, \dots, \ell_d) \in \mathbb{N}^d$ be a multi-index specifying the mesh level in each coordinate direction. The corresponding mesh sizes are given by the $d$-dimensional vector

$$ h_{\ell} := (2^{-\ell_1}, \dots, 2^{-\ell_d}) \in \mathbb{R}^d. $$

These mesh sizes partition the $d$-dimensional parameter domain into elements constructed by tensor products of unidimensional dyadic meshes:

$$ \hat{\Omega}_{h_{\ell}} := \bigotimes_{i=1}^{d} \hat{\Omega}_{h_{\ell_i}}, \quad \text{where} \quad \hat{\Omega}_{h_{\ell_i}} := \bigcup_{j=0}^{h_{\ell_i}^{-1}-1} (h_{\ell_i} j, h_{\ell_i}(j+1)). $$

The family $\{\hat{\Omega}_{h_{\ell}}\}_{\ell}$ consists of anisotropic meshes indexed by the mesh levels $\ell$, where we allow the mesh size to vary along different directions. This flexibility is essential for defining sparse-grid approximations.
# 3 Standard Approximation Spaces

In this section, we introduce finite-dimensional approximation spaces

$$S_{p,h_l}(\hat{\Omega}) \subset H^q(\hat{\Omega})$$

defined on the tensor-product domain $\hat{\Omega} = (0, 1)^d$. These spaces provide increasingly accurate approximations in $H^q(\hat{\Omega})$ as $h_l \to 0$. These spaces will later serve as building blocks for the construction of sparse-grid approximation spaces. We begin by defining the univariate spaces and then extend the construction to the $d$-dimensional setting using tensor-product formulations.

## 3.1 Univariate B-splines Spaces

In this section, we consider the unidimensional case corresponding to $d = 1$, i.e., $\hat{\Omega} = (0, 1)$. The mesh level is denoted by $l \in \mathbb{N}$.

**Definition 3.1.** The space of univariate B-splines is defined by

$$S_{p,k,h_l}(\hat{\Omega}) := \{ u \in H^k(\hat{\Omega}) \mid u|_{(h_l,j,h_l(j+1))} \in \mathbb{P}^p, \forall j = 0, \dots, h_l^{-1} - 1 \}.$$

Here, $\mathbb{P}^p$ is the space of polynomials of degree $p \in \mathbb{N}_0$, and the parameter $k$, with $0 \le k \le p$, specifies the spline regularity.

We restrict ourselves to the case of B-splines with maximal smoothness, i.e., $k = p$, and denote

$$S_{p,h_l}(\hat{\Omega}) := S_{p,p,h_l}(\hat{\Omega}).$$

For $r \in \mathbb{N}_0$ verifying $0 \le r \le p + 1$, we introduce the seminorm $H^r$-orthogonal projection operator onto the spline space and its complementary operator:

$$\pi_{p,h_l}: H^r(\hat{\Omega}) \to S_{p,h_l}(\hat{\Omega}), \\ (I - \pi_{p,h_l}): H^r(\hat{\Omega}) \to H^r(\hat{\Omega})/S_{p,h_l}(\hat{\Omega}),$$

defined such that

$$|(I - \pi_{p,h_l}) u|_{H^r(\hat{\Omega})} = \min_{v \in S_{p,h_l}(\hat{\Omega})} |u - v|_{H^r(\hat{\Omega})}.$$

The main result of this section, originally presented in [24], Theorem 7, provides an estimate of the approximation error in the $H^r$ seminorm for B-spline spaces with regularity up to $k = p$, with a constant independent of both the spline degree and the mesh size. We recall this result here, expressed with our notations.

**Lemma 3.1.** Let $u \in H^q(\hat{\Omega})$ with $q \in \mathbb{N}_0$. For all mesh level $l \in \mathbb{N}$, each spline degree $p \in \mathbb{N}_0$ such that $h_l p < 1$, and any $r \in \mathbb{N}_0$ satisfying $0 \le r \le q \le p + 1$, the following estimate holds true:

$$|(I - \pi_{p,h_l}) u|_{H^r(\hat{\Omega})} \le C_1(q,r) h_l^{q-r} |u|_{H^q(\hat{\Omega})},$$

with $C_1(q,r) = (\sqrt{2})^{q-r}$.

**Proof of Lemma 3.1.** According to [24], Theorem 7, under the stated assumptions, there exists a spline approximation $u_{p,h_l} \in S_{p,h_l}(\hat{\Omega})$ such that:

$$|u - u_{p,h_l}|_{H^r(\hat{\Omega})} \le (\sqrt{2}h_l)^{q-r}|u|_{H^q(\hat{\Omega})}.$$

Applying the definition of the seminorm $H^r$-orthogonal projection, we obtain the result. $\square$

<figure><img src="image_5.png" /></figure>
Note that, since

$$S_{p,h_l}(\hat{\Omega}) \subseteq S_{p,k,h_l}(\hat{\Omega}),$$

for all $0 \le k < p$, Lemma 3.1 is also valid in that case.

We also introduce, as in [24], a suitable subspace of $S_{p,h_l}(\hat{\Omega})$ to derive inverse inequalities.

**Definition 3.2.** The space of $q$-vanishing derivative univariate B-splines is defined by

$$\tilde{S}_{p,h_l}^{(q)}(\hat{\Omega}) := \left\{ u_{p,h_l} \in S_{p,h_l}(\hat{\Omega}) \left| \frac{\partial^{2l+q}}{\partial x^{2l+q}} u_{p,h_l}(0) = \frac{\partial^{2l+q}}{\partial x^{2l+q}} u_{p,h_l}(1) = 0 \quad \text{for all } l \in \mathbb{N}_0 \text{ with } 2l+q < p \right. \right\}$$

Note that according to [24], Theorem 5, Lemma 3.1 also holds for the space $\tilde{S}_{p,h_l}^{(q)}(\hat{\Omega})$. An inverse inequality, which is a result from [24], Theorem 6, can be derived within this space.

**Lemma 3.2.** For all mesh levels $l \in \mathbb{N}$, each $p \in \mathbb{N}_0$ such that $h_l p < 1$ and each $q \in \mathbb{N}_0$ with $0 \le q \le p$,

$$|u_{p,h_l}|_{H^q(0,1)} \le C_2(q) h_l^{-q} \|u_{p,h_l}\|_{L^2(0,1)}$$

is satisfied for all spline approximation $u_{p,h_l} \in \tilde{S}_{p,h_l}^{(q)}(\hat{\Omega})$ with $C_2(q) = (2\sqrt{3})^q$.

*Proof of Lemma 3.2.* From [24], Theorem 6, by induction, we obtain the result. $\square$

<figure><img src="image_6.png" /></figure>

## 3.2 Anisotropic $d$-dimensional B-spline Spaces

We extend the results of the previous section to anisotropic meshes in multiple dimensions, where $d \ge 2$. The tensor-product parameter domain is $\hat{\Omega} = (0, 1)^d$.

**Definition 3.3.** Let $l = (l_1, \dots, l_d) \in \mathbb{N}^d$ be the mesh level in each direction. The anisotropic $d$-dimensional B-spline space is defined by the tensor products of the univariate spaces:

$$S_{p,h_l}(\hat{\Omega}) := \bigotimes_{i=1}^{d} S_{p,h_{l_i}}(0,1).$$

Here, we consider B-splines of the same degree $p$ in each direction, while allowing the mesh size to vary along different directions.

To extend the notion of projection introduced in the previous section to this $d$-dimensional space, we begin by defining the projection along a single direction.

**Definition 3.4.** The $H^r$-orthogonal projection onto the spline space $S_{p,h_{l_i}}(0,1)$ in the single direction $i \in \{1, \dots, d\}$ is defined by

$$\Pi_{p,h_{l_i}} := I \otimes \dots \otimes \Pi_{p,h_{l_i}} \otimes \dots \otimes I$$

Note that the unidimensional projections are commutative, i.e., for $i_1, i_2 \in \{1, \dots, d\}$,

$$\Pi_{p,h_{l_1}} \circ \Pi_{p,h_{l_2}} = \Pi_{p,h_{l_2}} \circ \Pi_{p,h_{l_1}}.$$

We introduce the projection onto the B-spline space in a subset of dimension indices by successively applying the unidimensional projection.

**Definition 3.5.** Let $k \in \mathbb{N}_0$ with $0 \le k \le d$ and $J \subset \{1, \dots, d\}$ be a subset of cardinality $|J| = k$. The $k$-dimensional projection onto the B-spline space in the directions $J$ is defined by

$$\Pi_{p,h_l}^J u := \left( \prod_{i \in J} \Pi_{p,h_{l_i}} \right) u.$$
We also define the complementary projection by:

$$ (I - \Pi_{p,h_l})^J u := \left( \prod_{i \in J}^{\circ} (I - \Pi_{p,h_{l_i}}) \right) u. $$

Here the notation $\prod_{i \in J}^{\circ}$ indicates composition of the operators with respect to the indices in $J$, the order being irrelevant since the operators commute.

Based on these definitions, the projection of a function $u$ onto the $d$-dimensional spline space $S_{p,h_l}(\hat{\Omega})$ can be defined as the composition of the univariate $H^r$-orthogonal projections applied successively along each coordinate direction, i.e.,

$$ \Pi_{p,h_l} u := \prod_{p,h_l}^{\{1,...,d\}} u. \quad (3) $$

We introduce an approximation error estimate for the B-spline spaces on anisotropic meshes.

**Proposition 3.1.** Let $u \in H_{\text{mix}}^q(\hat{\Omega})$ with $q \in \mathbb{N}_0$. For all mesh levels $\ell \in \mathbb{N}^d$, each $p \in \mathbb{N}_0$, and any $r \in \mathbb{N}$ such that $0 \le r \le q \le p+1$,

$$ \|(I - \Pi_{p,h_l}) u\|_{H^r(\hat{\Omega})} \le C_3(d,q,r) \left( \sum_{i=1}^{d} h_{\ell_i}^{q-r} \right) \|u\|_{H^q(\hat{\Omega})} $$

is satisfied with $C_3(d,q,r) = (r+1)^{d/2} \sqrt{2}^{q-r}$

*Proof of Proposition 3.1.* The proof is strictly similar to [24], Theorem 8, but considering different mesh sizes $h_{\ell_i}$, $i = 1, \dots, d$, in each direction. $\square$

The result of Proposition 3.1 is not suitable for proving sharp approximation properties of the sparse-grid spaces, as it does not provide sufficient information in the mixed directions. In the remainder of this section, we introduce three lemmas that will be used later to study the approximation properties of these spaces.

First, in Lemma 3.4, we decompose the $d$-dimensional projection defined in Eq. (3) into a telescopic sum of partial projections over subsets of coordinate directions. This decomposition will allow us to exploit cancellations of coarse errors in the sparse-grid combination technique. To prove Lemma 3.4, we introduce the combinatorial result of Lemma 3.3.

**Lemma 3.3.** For all $i \in \mathbb{N}_0$ such that $0 \le i \le d-2$, we have:

$$ \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} l^i = 0. \quad (4) $$

*Proof of Lemma 3.3.* We demonstrate the result by induction. For $i=0$, the relation is immediate using the binomial theorem. We assume that Eq. (4) holds for $0 \le i \le d-3$, then we have:

$$ \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} l^{i+1} = \sum_{l=1}^{d-1} (-1)^l \binom{d-1}{l} l^{i+1} \quad (5) $$

$$ = (d-1) \sum_{l=1}^{d-1} (-1)^l \binom{d-2}{l-1} l^i \quad (6) $$

$$ = -(d-1) \sum_{l=0}^{d-2} (-1)^l \binom{d-2}{l} (l+1)^i \quad (7) $$

$$ = -(d-1) \sum_{j=0}^{i} \binom{i}{j} \sum_{l=0}^{d-2} (-1)^l \binom{d-2}{l} l^j \quad (7) $$

$$ = 0. \quad (8) $$
In the above, we have used the relation

$$ \binom{d-1}{l} = \frac{d-1}{l} \binom{d-2}{l-1} $$

in Eq. (6); the binomial theorem in Eq. (7); and the induction hypothesis in Eq. (8).

**Lemma 3.4.** For all mesh levels $l \in \mathbb{N}^d$ and all $u \in L^2(\hat{\Omega})$ the following decomposition holds true:

$$ (I - \Pi_{p,h_l})u = \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} (-1)^{k-1} (I - \Pi_{p,h_l})^J u. \quad (9) $$

*Proof of Lemma 3.4.* By expanding each of the following products as

$$ (I - \Pi_{p,h_l})^J u = I + \sum_{l=1}^{k} (-1)^l \sum_{\substack{\tilde{J} \subset J \\ |\tilde{J}|=l}} \Pi_{p,h_l}^{\tilde{J}}, $$

the right hand side of Eq. (9) can be rewritten as

$$ \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} (-1)^{k-1} (I - \Pi_{p,h_l})^J u = I \left[ \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} (-1)^{k-1} \right] \\ + \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} (-1)^{k-1} \sum_{l=1}^{k} (-1)^l \sum_{\substack{\tilde{J} \subset J \\ |\tilde{J}|=l}} \Pi_{p,h_l}^{\tilde{J}}. \quad (10) $$

The first term of the right hand side of Eq. (10) is equal to I because:

$$ \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} (-1)^{k-1} = \sum_{k=1}^{d} (-1)^{k-1} \binom{d}{k} \quad (11) $$

$$ = - \sum_{k=0}^{d} (-1)^k \binom{d}{k} - (-1) \quad (12) $$

$$ = 1. \quad (13) $$

Eq. (11) is obtained by counting the number of elements in the sum on $J \subset \{1, \dots, d\}$, $|J|=k$; in Eq. (12), we have added and subtracted the term $k=0$ to the sum; in Eq. (13), we have used Eq. (4) for $i=0$, which is a direct consequence of the binomial theorem. We recast the second term of the right hand side of Eq. (10) by rearranging the sums as

$$ \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} (-1)^{k-1} \sum_{l=1}^{k} (-1)^l \sum_{\substack{\tilde{J} \subset J \\ |\tilde{J}|=l}} \Pi_{p,h_l}^{\tilde{J}} = \sum_{l=1}^{d} \sum_{k=l}^{d} (-1)^{l+k-1} \underbrace{\sum_{J \subset \{1, \dots, d\} \atop |J|=k} \sum_{\tilde{J} \subset J \atop |\tilde{J}|=l} \Pi_{p,h_l}^{\tilde{J}}}_{P(l)} $$

and we look at each term $P(l)$, for $1 \le l \le d$, of the outermost sum. First for $l=d$, we directly obtain:

$$ P(d) = -(\Pi_{p,h_l})^{\{1, \dots, d\}}. $$
For $1 \le l \le d-1$, we rearrange the sums appearing in $P(l)$ as

$$P(l) = \sum_{\substack{\tilde{l} \subset \{1, \dots, d\} \\ |\tilde{l}| = l}} \prod_{p,h_l}^{\tilde{l}} \sum_{k=l}^{d} (-1)^{l+k-1} \binom{d-l}{k-l} \\ = 0.$$

We obtain the last equation because the innermost sum vanishes. Indeed, by substituting $\kappa = k-l$, we have:

$$\sum_{k=l}^{d} (-1)^{l+k-1} \binom{d-l}{k-l} = - \sum_{\kappa=0}^{d-l} (-1)^{\kappa} \binom{d-l}{\kappa} = 0.$$

□

Second, in Lemma 3.5, we give an estimate of the approximation error of each B-spline space defined by the partial projections appearing in the decomposition Eq. (9).

**Lemma 3.5.** Let $u \in H_{\text{mix}}^q(\hat{\Omega})$ with $q \in \mathbb{N}_0$. For all mesh levels $\ell \in \mathbb{N}^d$, each $p \in \mathbb{N}_0$, each $k \in \mathbb{N}$ with $1 \le k \le d$, each index set $J \subset \{1, \dots, d\}$ with $|J| = k$ and each $r \in \mathbb{N}$ such that $0 \le r \le q \le p+1$,

$$\|(I - \Pi_{p,h_\ell})^J u\|_{H_{\text{mix}}^r(\hat{\Omega})} \le C_3(d, k, q, r) \left(\prod_{i \in J} h_i^{q-r}\right) \|u\|_{H_{\text{mix}}^q(\hat{\Omega})}$$

is satisfied with $C_3(d, k, q, r) = (r+1)^{d/2}(\sqrt{2})^{k(q-r)}$.

*Proof of Lemma 3.5.* We demonstrate the result in two dimensions, for $d=2$. The generalization to higher dimensions is straightforward and uses the same arguments.

We start by estimating the term corresponding to the projection in both directions, i.e., $k=2$ and $J = \{1,2\}$. Let $s, m \in \mathbb{N}$ be such that $0 \le s, m \le r$. Since $I - \Pi_{p,h_{l_2}}$ is the $H^r$-orthogonal complementary projection, it is a linear operator and it verifies:

$$\| I - \Pi_{p,h_{l_2}} \|_{\mathcal{L}(H^q(\hat{\Omega}))} = 1.$$

Then, for all $0 \le s \le r$,

$$\left\| \frac{\partial^s}{\partial y^s} \left[ (I - \Pi_{p,h_{l_2}}) u(x,y) \right] \right\|_{H^r(\hat{\Omega})} \le \left\| \frac{\partial^s}{\partial y^s} u(x,y) \right\|_{H^r(\hat{\Omega})}, \quad \forall x, y \in (0,1).$$

As a result,

$$\frac{\partial^s}{\partial y^s} \left[ (I - \Pi_{p,h_{l_2}}) u(\cdot, y) \right] \in H^r(0,1) \subset H^m(0,1), \quad \forall y \in (0,1),$$

because $u \in H^q(\hat{\Omega})$ and $0 \le m \le r \le q$. Because of linearity we have:

$$\frac{\partial^s}{\partial y^s} \left[ (I - \Pi_{p,h_l})^{\{1,2\}} u(\cdot, y) \right] = (I - \Pi_{p,h_{l_1}}) \left[ \frac{\partial^s}{\partial y^s} (I - \Pi_{p,h_{l_2}}) u(\cdot, y) \right].$$

According to Lemma 3.1,

$$\left| \frac{\partial^s}{\partial y^s} (I - \Pi_{p,h_l})^{\{1,2\}} u(\cdot, y) \right|_{H^m(0,1)} \le (\sqrt{2}h_{l_1})^{q-m} \left| \frac{\partial^s}{\partial y^s} (I - \Pi_{p,h_{l_2}}) u(\cdot, y) \right|_{H^q(0,1)} . \quad (14)$$
By squaring and taking the integral in y on both sides, we obtain:

$$ \left\| \frac{\partial^m}{\partial x^m} \frac{\partial^s}{\partial y^s} (I - \Pi_{p,h_l})^{\{1,2\}} u \right\|_{L^2(\hat{\Omega})}^2 \le (\sqrt{2}h_{l_1})^{2(q-m)} \int_0^1 |(I - \Pi_{p,h_l}) \frac{\partial^q}{\partial x^q} u(x, \cdot)|_{H^s(0,1)}^2 dx. \quad (15) $$

Since $\frac{\partial^q}{\partial x^q} u(x, \cdot) \in H^s(0, 1)$ for all $x \in (0, 1)$, we can apply again Lemma 3.1 to obtain, for all $0 \le s, m \le r$,

$$ \begin{aligned} \left\| \frac{\partial^m}{\partial x^m} \frac{\partial^s}{\partial y^s} (I - \Pi_{p,h_l})^{\{1,2\}} u \right\|_{L^2(\hat{\Omega})}^2 & \le (\sqrt{2}h_{l_1})^{2(q-m)} (\sqrt{2}h_{l_2})^{2(q-s)} \int_0^1 \left| \frac{\partial^q}{\partial x^q} u(x, \cdot) \right|_{H^q(0,1)}^2 dx \\ & \le (2h_{l_1}h_{l_2})^{2(q-r)} \|u\|_{H_{\text{mix}}^q(\hat{\Omega})}^2. \end{aligned} \quad (16) $$

The second inequality is obtained by using the relations:

$$ \sqrt{2}h_{l_i} < 1 \quad \forall i = 1, \dots, d, \quad \text{and} \quad \left\| \frac{\partial^q}{\partial x^q} \frac{\partial^q}{\partial y^q} u \right\|_{L^2(\hat{\Omega})} \le \|u\|_{H_{\text{mix}}^q(\hat{\Omega})} $$

Each term in the $H_{\text{mix}}^r(\hat{\Omega})$ norm can be bounded using Eq. (16), with $0 \le s, m \le r$. Since there are $(r+1)^d$ terms in the seminorm, we obtain the result for $k=2$.

Concerning the terms corresponding to projections along a single direction, i.e., $k=1$, since:

$$ \frac{\partial^s}{\partial y^s} u(\cdot, y) \in H^m(0, 1), \forall y \in (0, 1) \quad \text{and} \quad \frac{\partial^s}{\partial x^s} u(x, \cdot) \in H^m(0, 1), \forall x \in (0, 1), $$

we can apply Lemma 3.1 to these two functions and use the same steps as in Eq. (14)-Eq. (16) to obtain:

$$ \left\| \frac{\partial^m}{\partial x^m} \frac{\partial^s}{\partial y^s} (I - \Pi_{p,h_j}) u \right\|_{L^2(\hat{\Omega})} \le (\sqrt{2}h_{j})^{q-r} \|u\|_{H_{\text{mix}}^q(\hat{\Omega})}, \quad \text{for } j=1,2. \quad (17) $$

We conclude by bounding all of the $(r+1)^d$ terms in the norm using Eq. (17). $\square$

<figure><img src="image_7.png" /></figure>

# 4 Sparse-Grid Approximation Spaces

In this section we introduce finite-dimensional approximation spaces defined by the sparse-grid tensor products of the unidimensional spaces $S_{p,h_l}(0, 1)$. We first introduce in Section 4.1 the sparse-grid spaces in the tensor-product parameter domain $\hat{\Omega} = (0, 1)^d$ before extending in Section 4.2 the results to more general non-tensor-product domains defined by $\Omega = F(\hat{\Omega})$, where $F$ is a geometric B-spline-based mapping.

## 4.1 Approximation Spaces in Tensor-Product Domains

Let $n \in \mathbb{N}$ be the maximum mesh level, i.e., $h_{l_i} \le h_n = 2^{-n}$, for all $i=1, \dots, d$.

We first introduce the so-called full tensor-product approximation space, constructed by taking uniform mesh size in all directions.

**Definition 4.1.** The full tensor-product approximation space is defined by

$$ S_{p,h_n}(\hat{\Omega}) := S_{p,h_l}(\hat{\Omega}) \quad \text{where } l_i = n, \forall i = 1, \dots, d. $$
In Section 4.1.1, we introduce a d-dimensional sparse-grid approximation space, denoted by $S_{p,h_n}^{\mathcal{H}}(\hat{\Omega})$, which is a subspace of the full tensor-product space:

$$S_{p,h_n}^{\mathcal{H}}(\hat{\Omega}) \subset S_{p,h_n}(\hat{\Omega}).$$

The superscript $\mathcal{H}$ indicates that the space is constructed using a hierarchical formulation. The construction proceeds by expressing the full tensor product space as a direct sum of hierarchical subspace increments, and by neglecting the less significant contributions under the assumption that the target functions possess sufficiently smooth mixed derivatives. In this setting, the approximation acts on the space $H_{\text{mix}}^q(\hat{\Omega})$ rather than on $H^q(\hat{\Omega})$.

In Section 4.1.2, we introduce a second sparse-grid approximation space, denoted by $S_{p,h_n}^{\mathcal{L}}(\hat{\Omega})$, based on the sparse-grid combination technique and first proposed in [3]. The combination technique can be viewed as an approximation of hierarchical-based sparse-grid methods and is known, under suitable regularity assumptions, to achieve the same approximation order as the latter [13, 18, 21].

### 4.1.1 Hierarchical Subspace Decomposition

A key element in the construction of sparse-grid tensor-product spaces is the definition of hierarchical subspaces, also called hierarchical increments, yielding a direct sum decomposition. Specifically, we want to express the anisotropic approximation spaces as the direct sum decomposition

$$S_{p,h_{\ell}}(\hat{\Omega}) = \bigoplus_{k \le \ell} W_{p,h_k}(\hat{\Omega}), \qquad (18)$$

where the multi-index notation $k \le \ell$ stands for component-wise inequalities, i.e., $k_j \le \ell_j$ for all $j = 1, \dots, d$.

A crucial property to obtain Eq. (18) is that the B-spline spaces are nested under refinement via knot insertion:

$$S_{p,h_{\ell-e_j}}(\hat{\Omega}) \subset S_{p,h_{\ell}}(\hat{\Omega}), \quad \forall j = 1, \dots, d,$$

where $e_j$ denotes the unit vector in the j-th direction. Since knot insertion preserves nestedness, every new function introduced during refinement is linearly independent of the previous ones. The hierarchical subspaces are then defined as

$$W_{p,h_k}(\hat{\Omega}) := S_{p,h_k}(\hat{\Omega}) / \bigoplus_{j=1}^{d} S_{p,h_{k-e_j}}(\hat{\Omega}), \qquad (19)$$

where

$$S_{p,h_k}(\hat{\Omega}) := 0 \quad \text{if } \exists j \in \{1, \dots, d\} \text{ s.t. } k_j = -1.$$

By defining the hierarchical increments as in Eq. (19), Eq. (18) is satisfied. Intuitively, the hierarchical increments $W_{p,h_k}(\hat{\Omega})$ spans the new basis functions introduced at the refinement step.

Using the definition of these hierarchical subspaces, we introduce the hierarchical-based sparse-grid approximation space.

**Definition 4.2.** The hierarchical sparse-grid approximation space is defined by

$$S_{p,h_n}^{\mathcal{H}}(\hat{\Omega}) := \{ u \in H_{\text{mix}}^p(\hat{\Omega}) \mid u = \sum_{\ell \in \mathcal{H}} w_{p,h_{\ell}}, \ w_{p,h_{\ell}} \in W_{p,h_{\ell}}(\hat{\Omega}) \},$$

where

$$\mathcal{H} = \{ \ell \in \mathbb{N}^d \mid |\ell|_1 \le n + (d-1)\lambda_p, \ \ell_j \ge \lambda_p \},$$
and

$$ \lambda_p := \begin{cases} 1, & \text{if } p=1, \\ \frac{\log p}{\log 2}, & \text{otherwise.} \end{cases} $$

The constraint $\ell_j \ge \lambda_p$ ensures that $h_{\ell_j} p < 1$ for all $j = 1, \dots, d$, a condition required to apply Theorem 7 from [24], which underlies our analysis.

### 4.1.2 Combination Technique

The sparse-grid combination technique is an alternative formulation for representing functions from hierarchical sparse-grid spaces, such as $S_{p,h_n}^{\mathcal{H}}(\hat{\Omega})$. In general, it yields an approximation of $H_{\text{mix}}^q(\hat{\Omega})$ of the same order of accuracy [13, 21]. The method relies on forming a linear combination of contributions that live in standard tensor-product approximation spaces defined with admissible levels $\ell$ satisfying $|\ell|_1 \approx n$.

Specifically, the admissible levels are such that $\ell \in \mathcal{L}$. We introduce the index set of admissible levels as

$$ \mathcal{L} := \bigcup_{0 \le l \le d-1} \mathcal{L}_l, \quad \text{where } \mathcal{L}_l := \left\{ \ell \in \mathbb{N}^d \mid |\ell|_1 = n + (d-1)\lambda_p - l, \ell_j \ge \lambda_p \right\}. \quad (20) $$

We define the sparse-grid combination technique approximation space by combining the spaces $S_{p,h_\ell}(\hat{\Omega})$ associated with the admissible levels.

**Definition 4.3.** The sparse-grid combination technique approximation space is defined by

$$ S_{p,h_n}^{\mathcal{L}}(\hat{\Omega}) := \left\{ u \in H_{\text{mix}}^p(\hat{\Omega}) \mid u = \sum_{\ell \in \mathcal{L}} c_{\ell} u_{p,h_{\ell}}, \quad u_{p,h_{\ell}} \in S_{p,h_{\ell}}(\hat{\Omega}) \right\}, $$

with the combination coefficients given by

$$ c_{\ell} := (-1)^l \binom{d-1}{l}, \quad \text{for } \ell \in \mathcal{L}_l. $$

We extend the definition of the $d$-dimensional projection introduced in Eq. (3) to this space by defining the projection as the combination of projections associated with the admissible levels:

$$ \Pi_{p,h_n}^{\mathcal{L}} := \left[ \sum_{\ell \in \mathcal{L}} c_{\ell} \Pi_{p,h_{\ell}} \right]. $$

### 4.1.3 Equivalence Between Spaces

In the previous sections, we have introduced two approximation spaces based on sparse-grid tensor products of univariate B-spline spaces. We now establish that the two approximation spaces are identical.

**Theorem 4.1.** The hierarchical subspace decomposition and the combination technique formulations yield identical spaces, i.e.,

$$ S_{p,h_n}^{\mathcal{L}}(\hat{\Omega}) = S_{p,h_n}^{\mathcal{H}}(\hat{\Omega}). $$

To give the proof of this theorem, we have to introduce two auxiliary lemmas. Lemma 4.2 establishes that the sum of the combination coefficients over the admissible levels equals one, while Lemma 4.1 is a combinatorial result used in the proof of Lemma 4.2.
**Lemma 4.1.** Let $\lambda_p$ be defined as in Eq. (20), $\ell, k \in \mathbb{N}$ be such that $0 \le \ell \le n$ and $0 \le k \le d-1$, then:

$$ \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \binom{n+d-1-\lambda_p-l-\ell}{d-1-k} = \begin{cases} 1 & \text{if } k=0, \\ 0 & \text{else.} \end{cases} \qquad (21) $$

*Proof of Lemma 4.1.* Eq. (21) can be recast as:

$$ \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \binom{n+d-1-\lambda_p-l-\ell}{d-1-k} = \sum_{l=0}^{d-1} \binom{d-1}{l} \frac{(-1)^l}{(d-1-k)!} \prod_{j=1}^{d-1-k} (n-l-\ell-\lambda_p+k+j). $$

By expanding the products in the right side we obtain a sum of terms with powers $(-l)^i$ like the left side of Eq. (4) in Lemma 3.3, with $0 \le i \le d-1$. As a result, the expression vanishes for $k > 0$. If $k = 0$, only one term remains, the one with $(-l)^{d-1}$. Then by applying $d-2$ times the same calculations than in Eq. (5)-Eq. (7), we obtain:

$$ \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} (-l)^{d-1} = (-1)^{d-2} (-1)^{d-1} \frac{(d-1)!}{(d-1)!} \sum_{l=0}^{1} (-1)^l \binom{1}{l} l = 1. \qquad (22) $$

□

**Lemma 4.2.** Let $\mathcal{L}$ be the set of admissible levels in the combination, and $c_\ell$ the combination coefficients, then:

$$ \sum_{\ell \in \mathcal{L}} c_{\ell} = 1. \qquad (23) $$

*Proof of Lemma 4.2.* The sum can be recast into:

$$ \begin{aligned} \sum_{\ell \in \mathcal{L}} c_{\ell} &= \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \sum_{\ell \in \mathcal{L}_l} 1 \\ &= \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \binom{n+d-1-\lambda_p-l}{d-1} \end{aligned} \qquad (24) $$

$$ = 1. \qquad (25) $$

We obtain Eq. (24) by counting the number of elements in $\mathcal{L}_l$; and Eq. (25) by using Lemma 4.1.

□

*Proof of Theorem 4.1.* Let $u_{p,h_n} \in S_{p,h_n}^{\mathcal{L}}(\hat{\Omega})$, then by definition:

$$ u_{p,h_n} = \sum_{\ell \in \mathcal{L}} c_{\ell} u_{p,h_{\ell}}, \quad \text{where } u_{p,h_{\ell}} \in S_{p,h_{\ell}}(\hat{\Omega}). $$

Using the definition of the combination coefficients and the direct sum decomposition in Eq. (18), we have:

$$ u_{p,h_n} = \sum_{l=0} (-1)^l \binom{d-1}{l} \sum_{\ell \in \mathcal{L}_l} \sum_{k \le \ell} w_{p,h_k} \qquad (26) $$

$$ = \sum_{\mathbf{k}} w_{p,h_k} N(\mathbf{k}, n, d, p). \qquad (27) $$
We have rearranged the sums in Eq. (26) and introduced in Eq. (27) the following quantity:

$$N(\mathbf{k}, n, d, p) = \begin{cases} \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \sum_{\ell \in \mathcal{L}_l} 1 & \text{if } |\mathbf{k}|_1 \le n + (d-1)\lambda_p, \\ 0 & \text{if } |\mathbf{k}|_1 > n + (d-1)\lambda_p. \end{cases}$$

By counting the number of elements in $\mathcal{L}_l$, and using Lemma 4.1, we have:

$$N(\mathbf{k}, n, d, p) = \begin{cases} 1 & \text{if } |\mathbf{k}|_1 \le n + (d-1)\lambda_p, \\ 0 & \text{if } |\mathbf{k}|_1 > n + (d-1)\lambda_p. \end{cases}$$

As a result, $u_{p,h_n} \in S_{p,h_n}^{\mathcal{H}}(\hat{\Omega})$ and $S_{p,h_n}^{\mathcal{L}}(\hat{\Omega}) \subset S_{p,h_n}^{\mathcal{H}}(\hat{\Omega})$. By performing the same substitutions in the sums in reverse, we can also show that $S_{p,h_n}^{\mathcal{H}}(\hat{\Omega}) \subset S_{p,h_n}^{\mathcal{L}}(\hat{\Omega})$, which conclude the proof. $\square$

As a result of Theorem 4.1, the two spaces —and thus their approximation properties— are the same. We shall then derive properties of this space using either the combination technique or the hierarchical formulation, choosing the approach that makes a particular proof more straightforward.

However, although the two spaces are mathematically identical, their computational properties may differ due to the distinct constructions. In particular when applied to Galerkin methods, the size, condition number, and sparsity patterns of the resulting discretization matrices, such as the mass and stiffness matrices, are not the same. A detailed investigation of these computational aspects is deferred to future work.

### 4.1.4 Approximation Error Estimate

In this section, we provide an approximation error estimate in mixed Sobolev norm for the sparse-grid space. To prove the result, we consider the combination technique formulation of the sparse-grid space. Before stating the main result of this section, we introduce Lemma 4.3 explicating the cancellations of the coarse contributions in the combination technique.

**Lemma 4.3.** Let $k \in \mathbb{N}$ with $1 \le k \le d$, and $J \subset \{1, \dots, d\}$ be an index set of cardinality $|J| = k$. Suppose $\{a_{J,\ell}\}_{J \subset \{1, \dots, d\}}$ is a family of $d$-dimensional functions depending on the directions indexed by $J$. Then the following identity holds:

$$\sum_{\ell \in \mathcal{L}} c_{\ell} \sum_{k=1}^{d} \sum_{\substack{J \subset \{1, \dots, d\} \\ |J|=k}} a_{J, \ell} = \sum_{k=1}^{d-1} \sum_{l=0}^{d-2} C_3(d, k, l) \sum_{\substack{J \subset \{1, \dots, d\} \\ |J|=k}} \sum_{\ell \in \mathcal{L}_{l,k}^J} a_{J, \ell} + \sum_{l=0}^{d-1} C_4(d, l) \sum_{\ell \in \mathcal{L}_l} a_{\{1, \dots, d\}, \ell}, \quad (28)$$

where the notation

$$\mathcal{L}_{l,k}^J := \left\{ \ell \in \mathbb{N}^d \mid \sum_{i \in J} |\ell_i| = n + (k-1)\lambda_p - l, \ \ell_{J_j} \ge \lambda_p, \text{ for } j = 1, \dots, k \right\}$$

is introduced. The constants are given by

$$C_3(d, k, l) = \sum_{\kappa=0}^{l} (-1)^{\kappa} \binom{d-1}{\kappa} \binom{d-k-1+l-\kappa-k\lambda_p}{d-1-k}, \quad C_4(d, l) = (-1)^l \binom{d-1}{l}.$$
Proof of Lemma 4.3. We recast the expression of the left side of Eq. (28) by introducing the combination coefficient expression:

$$ \sum_{\ell \in \mathcal{L}} c_{\ell} \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} a_{J, \ell} = \sum_{k=1}^{d} \sum_{J \subset \{1, \dots, d\} \atop |J|=k} \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \sum_{\ell \in \mathcal{L}_l} a_{J, \ell}. \quad (29) $$

Note that, in the above, the admissible levels verify a condition on the $\ell^1$ norm of all components and not only the $J$ directions.

If $k=d$, we directly obtain the second term of Eq. (28), for which there is no cancellation. If $k \le d-1$, cancellations between the different admissible levels occur in the combination. The first step is to the rearrange the innermost sums of the right side of Eq. (29), so that we can extract all the terms that will be cancelled. Specifically, we have:

$$ \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \sum_{\ell \in \mathcal{L}_l} a_{J, \ell} \\ = \sum_{l=0}^{d-1} (-1)^l \binom{d-1}{l} \sum_{\ell \in \mathcal{L}_{d-1,k}^J} \binom{n+d-1-\lambda_p-l-(\ell_{J_1}+\dots+\ell_{J_k})}{d-1-k} a_{J, \ell} \quad (30) $$

$$ = \sum_{\ell \in \mathcal{L}_{d-1,k}^J} a_{J, \ell} \sum_{\kappa=0}^{d-1} (-1)^\kappa \binom{d-1}{\kappa} \binom{n+d-1-\kappa-\lambda_p-(\ell_{J_1}+\dots+\ell_{J_k})}{d-1-k}, \\ + \sum_{l=0}^{d-2} \sum_{\ell \in \mathcal{L}_{l,k}^J} a_{J, \ell} \sum_{\kappa=0}^{l} (-1)^\kappa \binom{d-1}{\kappa} \binom{n+d-1-\kappa-\lambda_p-(\ell_{J_1}+\dots+\ell_{J_k})}{d-1-k}, \quad (31) $$

We have introduced the set of levels

$$ \tilde{\mathcal{L}}_{l,k}^J := \{ \ell \in \mathbb{N}^d \mid \sum_{i \in J} |\ell_i| \le n + (k-1)\lambda_p - l, \ell_{J_i} \ge \lambda_p, \text{ for } j=1, \dots, k \}. $$

In Eq. (30), we have rearranged the sum on the levels, by summing in the $J$ directions and letting the remaining directions constraint by the $l^1$-norm equality. We have rearranged again the sum on the levels into two terms: whether if the level appears in each term $0 \le l \le d-1$ of the sum in Eq. (30), which correspond to the first term in Eq. (31); or not, which correspond to the second one.

According to Lemma 4.1, the first term in Eq. (31) is equal to zero. This corresponds to coarse contributions that are cancelled in the combination technique. Substituting $(\ell_{J_1} + \dots + \ell_{J_k})$ by its expression in Eq. (31), we obtain the result for $1 \le k \le d-1$. $\square$

With this lemma in hand, we can now state Lemma 4.4, which is the main result of this section.

**Lemma 4.4.** Let $u \in H_{\text{mix}}^q(\hat{\Omega})$ with $q \in \mathbb{N}_0$. For all mesh size $n \in \mathbb{N}$, each $p \in \mathbb{N}_0$ and each $r \in \mathbb{N}$ with $0 \le r \le q \le p+1$,

$$ \| (I - \Pi_{p,h_n}^{\mathcal{L}}) u \|_{H_{\text{mix}}^r(\hat{\Omega})} \le C_{10}(d,q,r) h_n^{q-r} |\log h_n|^{d-1} \|u\|_{H_{\text{mix}}^q(\hat{\Omega})} \quad (32) $$

is satisfied with

$$ C_{10}(d, q, r) = \frac{(d-1)^{d-1}}{(d-1)!(\log 2)^{d-1}} \sum_{l=0}^{d-2} 2^{-(q-r)(d-1-l)} (-1)^l \binom{d-1}{l} (r+1)^{d/2} (\sqrt{2})^{d(q-r)}, $$

a constant that does not depend on the mesh size $h_n$ or the degree $p$.
*Proof of Lemma 4.4.* The local interpolation error can be written as

$$ (I - \Pi_{p,h_n}^{\mathcal{L}})u = (I - \sum_{\ell \in \mathcal{L}} c_{\ell} \Pi_{p,h_{\ell}})u \\ = \sum_{\ell \in \mathcal{L}} c_{\ell} (I - \Pi_{p,h_{\ell}})u \qquad (33) $$

$$ = \sum_{\ell \in \mathcal{L}} c_{\ell} \sum_{k=1}^{d} \sum_{\substack{J \subset \{1, \dots, d\} \\ |J|=k}} (-1)^k (I - \Pi_{p,h_{\ell}})^J u \qquad (34) $$

$$ = \sum_{k=1}^{d-1} \sum_{l=1}^{d-2} C_4(d,k,l) (-1)^{k-1} \sum_{\substack{J \subset \{1, \dots, d\} \\ |J|=k}} \sum_{\ell \in \mathcal{L}_{l,k}^J} (I - \Pi_{p,h_{\ell}})^J u \\ + \sum_{l=0}^{d-1} C_5(d,l) \sum_{\ell \in \mathcal{L}_l} (I - \Pi_{p,h_{\ell}})^{\{1, \dots, d\}} u. \qquad (35) $$

In the above, Lemma 4.2 has been used in Eq. (33); Lemma 3.4 has been used in Eq. (34); and Lemma 4.3 has been used in Eq. (35). Then, using the triangular inequality and applying Lemma 3.5 to each of the terms, we obtain:

$$ \| (I - \Pi_{p,h_n}^{\mathcal{L}}) u \|_{H_{\text{mix}}^r(\hat{\Omega})} \le \sum_{k=1}^{d-1} \sum_{l=0}^{d-2} C_6(d,r,k,l) \sum_{\substack{J \subset \{1, \dots, d\} \\ |J|=k}} \sum_{\ell \in \mathcal{L}_{l,k}^J} \left( \prod_{i \in J} h_{\ell_i}^{q-r} \right) \|u\|_{H_{\text{mix}}^q(\hat{\Omega})} \\ + \sum_{l=0}^{d-1} C_7(d,r,l) \sum_{\substack{\ell \in \mathcal{L}_l \\ l \in \{1, \dots, d\}}} \left( \prod_{i \in \{1, \dots, d\}} h_{\ell_i}^{q-r} \right) \|u\|_{H_{\text{mix}}^q(\hat{\Omega})}, \qquad (36) $$

where the constants are

$$ C_6(d,r,k,l) = C_3(d,k,l)C_2(d,k,q,r), \quad C_7(d,r,l) = C_4(d,l)C_2(d,d,q,r). $$

The number of elements in each sum is given by

$$ \sum_{\substack{J \subset \{1, \dots, d\} \\ |J|=k}} 1 = \binom{d}{k}, \quad \sum_{\substack{\ell \in \mathcal{L}_{l,k}^J \\ l \in \{1, \dots, d\}}} 1 = \binom{n+k-1-\lambda_p-l}{k-1}, \quad \sum_{\substack{\ell \in \mathcal{L}_l \\ l \in \{1, \dots, d\}}} 1 = \binom{n+d-1-\lambda_p-l}{d-1}. $$

Since the last two sums are maximal for $l=0$, Eq. (36) can be bounded by

$$ \| (I - \Pi_{p,h_n}^{\mathcal{L}}) u \|_{H_{\text{mix}}^r(\hat{\Omega})} \le \sum_{k=1}^{d-1} C_8(d,q,r,k) \binom{n+k-1-\lambda_p}{k-1} h_n^{q-r} \|u\|_{H_{\text{mix}}^q} \\ + C_9(d,q,r) \binom{n+d-1-\lambda_p}{d-1} h_n^{q-r} \|u\|_{H_{\text{mix}}^q}, $$

where

$$ C_8(d,q,r,k) = \binom{d}{k} \sum_{l=0}^{d-2} 2^{-(q-r)(k-1-l)} C_5(d,r,k,l), \quad C_9(d,q,r) = \sum_{l=0}^{d-2} 2^{-(q-r)(d-1-l)} C_6(d,r,l). $$

Eventually, we conclude with the bound

$$ \binom{n+d-1-\lambda_p}{d-1} \le \binom{n+d-2}{d-1} \le \frac{(d-1)^{d-1}}{(d-1)!(\log 2)^{d-1}} |\log h_n|^{d-1}, $$

so that the constant does not depend on $p$.

<figure><img src="image_8.png" /></figure>
On tensor-product domains such as $\hat{\Omega}$, it is well established that sparse-grid approximations using piecewise linear functions achieve an $L^2$ approximation error of order $\mathcal{O}(h_n^2|\log h_n|^{d-1})$, with a dependence on the $H_{\text{mix}}^2$ Sobolev norm of the solution [9, 19]. This result has been generalized to piecewise Lagrange polynomials of degree $p$, for which the $L^2$ error scales as $\mathcal{O}(h_n^{p+1}|\log h_n|^{d-1})$ and depends on the $H_{\text{mix}}^{p+1}$ Sobolev norm [9]. Lemma 4.4 further extends these results to B-splines of maximal smoothness, i.e., B-splines of degree $p$ with regularity $H^p(0,1)$ in each direction. Considering Lemma 4.4 with $r=0$ and $q=p+1$, the $L^2$ approximation error for such B-splines exhibits the same scaling as for Lagrange polynomials. Notably, as in the case of isotropic full tensor-product spaces treated in [24], the approximation error is independent of the spline degree $p$.

### 4.1.5 Inverse Inequality

In this section, we derive an inverse inequality in mixed Sobolev norms for the sparse-grid approximation space. To this end, we extend the definition of the unidimensional space of $q$-vanishing derivative univariate B-splines Definition 3.2 to the sparse-grid spaces. Specifically, we define the spaces $\tilde{S}_{p,h_n}^{\mathcal{H},(q)}(\hat{\Omega})$ and $\tilde{S}_{p,h_n}^{\mathcal{L},(q)}(\hat{\Omega})$ by replacing the univariate spaces $S_{p,h_l}(0,1)$ with their counterparts $\tilde{S}_{p,h_l}^{(q)}(0,1)$ in Definition 4.2 and Definition 4.3, respectively.

Using the same arguments as in the proof of Theorem 4.1, we can demonstrate that these two spaces are identical. Moreover, according to [24], Corollary 1, the estimate Eq. (32) of Lemma 4.4 also holds for spaces constructed with $q$-vanishing-derivative B-splines.

**Lemma 4.5.** For each $0 \le q \le p$, the inverse inequality

$$ \|u_{p,h_n}^{\mathcal{H}}\|_{H_{\text{mix}}^q(\hat{\Omega})} \le C_{11}(d,q)h_n^{-q}|\log h_n|^{d/2}\|u_{p,h_n}^{\mathcal{H}}\|_{L^2(\hat{\Omega})} \quad (37) $$

is satisfied for all $u_{p,h_n}^{\mathcal{H}} \in \tilde{S}_{p,h_n}^{\mathcal{H},(q)}(\hat{\Omega})$, with

$$ C_{11}(d,q) = (q+1)^{d/2}(2\sqrt{3})^{dq} \frac{2^{(d-1)2^{d/2}}}{d!(\log 2)^{d/2}}, \quad (38) $$

a constant that does not depend on the mesh size $h_n$ or the degree $p$.

To prove this result, we consider the hierarchical sparse-grid formulation. We introduce Lemma 4.6, an auxiliary result required to prove Lemma 4.5.

**Lemma 4.6.** Let $s \in \mathbb{N}^d$ with $0 \le s_1, \dots, s_d \le q \le p$, and for all mesh level $l \in \mathcal{H}$, each $p \in \mathbb{N}_0$, $w_{p,h_l} \in W_{p,h_l}(\hat{\Omega})$, then

$$ \left\| \frac{\partial^{s_1}}{\partial x_1^{s_1}} \cdots \frac{\partial^{s_d}}{\partial x_d^{s_d}} w_{p,h_l} \right\|_{L^2(\hat{\Omega})} \le (2\sqrt{3})^{|s|_1} \left( \prod_{i=1}^d h_{l_i}^{-s_i} \right) \|w_{p,h_l}\|_{L^2(\hat{\Omega})}. \quad (39) $$

*Proof of Lemma 4.6.* Like in the proof of Lemma 3.5, we demonstrate the result for $d=2$, as the extension to higher dimensions uses the same arguments. Let $0 \le s, m \le q$ be two integers. The hierarchical increment verifies:

$$ w_{p,h_l} \in W_{p,h_l}(\hat{\Omega}) \subset H^p(\hat{\Omega}) \subset H^q(\hat{\Omega}), $$

so that we have:

$$ \begin{gathered} \frac{\partial^s}{\partial y^s} w_{p,h_l}(\cdot, y) \in H^m(0,1), \quad \forall y \in (0,1), \\ w_{p,h_l}(x, \cdot) \in H^s(0,1), \quad \forall x \in (0,1). \end{gathered} $$
Therefore, we can demonstrate the bounds

$$ \begin{aligned} \left\| \frac{\partial^m}{\partial x^m} \frac{\partial^s}{\partial y^s} w_{p,h_l} \right\|_{L^2(\hat{\Omega})}^2 &= \int_0^1 \left| \frac{\partial^s}{\partial y^s} w_{p,h_l}(\cdot, y) \right|_{H^{m(0,1)}}^2 dy \\ &\le (2\sqrt{3})^{2m} h_{\ell_1}^{-2m} \int_0^1 \int_0^1 \left( \frac{\partial^s}{\partial y^s} w_{p,h_l}(x,y) \right)^2 dxdy \end{aligned} \quad (40) $$

$$ \le (2\sqrt{3})^{2m} h_{\ell_1}^{-2m} \int_0^1 |w_{p,h_l}(x, \cdot)|_{H^{s(0,1)}}^2 dx $$

$$ \le (2\sqrt{3})^{2(m+s)} h_{\ell_1}^{-2m} h_{\ell_2}^{-2s} \|w_{p,h_l}\|_{L^2(\hat{\Omega})}^2 \quad (41) $$

We have used Lemma 3.2 twice: in Eq. (40) and in Eq. (41). □

Proof of Lemma 4.5. The following inequalities are true:

$$ \left\| \frac{\partial^{s_1}}{\partial x_1^{s_1}} \cdots \frac{\partial^{s_d}}{\partial x_d^{s_d}} u_{p,h_n}^{\mathcal{H}} \right\|_{L^2(\hat{\Omega})} \le \sum_{\ell \in \mathcal{H}} \left\| \frac{\partial^{s_1}}{\partial x_1^{s_1}} \cdots \frac{\partial^{s_d}}{\partial x_d^{s_d}} w_{p,h_\ell} \right\|_{L^2(\hat{\Omega})} \quad (42) $$

$$ \le \sum_{\ell \in \mathcal{H}} (2\sqrt{3})^{|\mathcal{S}|_1} \left( \prod_{i=1}^d h_{\ell_i}^{-s_i} \right) \|w_{p,h_\ell}\|_{L^2(\hat{\Omega})} \quad (43) $$

$$ \le (2\sqrt{3})^{dq} \left( \sum_{\ell \in \mathcal{H}} \prod_{i=1}^d h_{\ell_i}^{-2s_i} \right)^{1/2} \|u_{p,h_n}^{\mathcal{H}}\|_{L^2(\hat{\Omega})} \quad (44) $$

We have used the triangular inequality in Eq. (42); applied Lemma 4.6 in Eq. (43); and used Cauchy-Schwarz inequality in Eq. (44). We then have the relation

$$ \sum_{\ell \in \mathcal{H}} \prod_{i=1}^{d} h_{\ell_i}^{-2s_i} \le h_n^{-2q} 2^{2(d-1)} \lambda_p \left( \sum_{\ell \in \mathcal{H}} 1 \right), \quad (45) $$

and we bound the last term of the product by counting the number of elements in $\mathcal{H}$:

$$ \sum_{\substack{\ell \in \mathcal{H} \\ |\mathcal{K}|_1 \le n - \lambda_p}} 1 \le \sum_{\substack{k \in \mathbb{N}^d \\ |k|_1 \le n - \lambda_p}} 1 \le \binom{n - \lambda_p + d}{d}. \quad (46) $$

Gathering Eq. (45) and Eq. (46) together, and since $\lambda_p \ge 1$, we obtain:

$$ \sum_{\ell \in \mathcal{H}} \prod_{i=1}^{d} h_{\ell_i}^{-2s_i} \le h_n^{-2q} 2^{2(d-1)} \binom{n+d-1}{d} \le h_n^{-2q} \log |h_n|^d \frac{2^{2(d-1)} 2^d}{d!(\log 2)^d}. \quad (47) $$

Eventually, we substitute the side of Eq. (47) relation into Eq. (44) and bound all the $(q+1)^d$ in the $H_{\text{mix}}^q$ full norm to get:

$$ \| u_{p,h_n}^{\mathcal{H}} \|_{H_{\text{mix}}^q(\hat{\Omega})} \le (q+1)^{d/2} (2\sqrt{3})^{dq} \frac{2^{(d-1)2^d/2}}{d!(\log 2)^{d/2}} h_n^{-q} \log |h_n|^{d/2}. \quad \square $$
## 4.2 Extension to Non-Tensor-Product Domains

We now extend the previous results to more general geometric domains. Let the physical domain be defined by a smooth mapping

$$F : \hat{\Omega} \to \Omega,$$

which is a diffeomorphism between the tensor-product parameter domain $\hat{\Omega}$ and the physical domain $\Omega$, i.e.,

$$\Omega := F(\hat{\Omega}).$$

The mapping $F$ provides a geometric parametrization of the physical domain.

Given meshes of the parameter and physical domains, denoted respectively by $\hat{\Omega}_{h\ell}$ and $\Omega_{h\ell}$, and an approximation space $S_{h\ell}(\hat{\Omega})$ defined on the parameter domain, we define the corresponding approximation space in the physical domain as

$$S_{h\ell}(\Omega) := \{ u_{h\ell} : \Omega \to \mathbb{R} \mid u_{h\ell} = \hat{u}_{h\ell} \circ F_{h\ell}^{-1}, \ \hat{u}_{h\ell} \in S_{h\ell}(\hat{\Omega}) \}.$$

Here, $F_{h\ell}$ denotes a $\mathscr{C}^1$-diffeomorphism that maps the parameter mesh $\hat{\Omega}_{h\ell}$ to the corresponding physical mesh $\Omega_{h\ell}$:

$$\Omega_{h\ell} := F_{h\ell}(\hat{\Omega}_{h\ell}).$$

Several strategies exist in the context of sparse grids to define $F_{h\ell}$, and consequently to construct the spaces $S_{h\ell}(\Omega)$. In [1], $F_{h\ell}$ is defined as the interpolant of $F$ evaluated at the nodes of $\Omega_{h\ell}$. In [8, 11], curvilinear meshes are employed and $F_{h\ell}$ is defined by transfinite interpolation. In [3], $F_{h\ell}$ is represented using B-splines or NURBS (rational B-splines with associated weights), where the geometry is defined as a linear combination of B-spline basis functions weighted by control points.

In this work, we adopt the isogeometric analysis representation, that is, we define $F_h$ as a B-spline mapping. We restrict ourselves to B-spline functions for simplicity; however, the extension to NURBS functions, as in [2], could also be considered. In practical applications, the geometry of the physical domain is typically described on a coarse mesh with relatively few elements, while the solution is approximated on successively finer meshes to achieve the desired accuracy. This is possible because B-spline or NURBS mappings can represent the geometry exactly, even on the coarsest mesh.

We assume that the geometry is described on the coarsest parameter mesh $\hat{\Omega}_{h0}$, corresponding to unit mesh size ($h_0 = 1$) in each coordinate direction. For any mesh size $h_\ell$, the physical mesh is then defined by

$$\Omega_h := F_{h0}(\hat{\Omega}_{h0}),$$

where the geometric mapping is given by

$$F_{h0}(\xi_1, \dots, \xi_d) = \sum_{i_1=1}^{n_1} \dots \sum_{i_d=1}^{n_d} P_{i_1, \dots, i_d} s_{p,h_0,i_1}(\xi_1) \dots s_{p,h_0,i_d}(\xi_d).$$

Here, $s_{p,h_0,i_j} \in S_{p,h_0}(0,1)$ denote the univariate B-spline basis functions defined on the coarsest mesh with $n_j = h_0^{-1}$ for $j = 1, \dots, d$, and $P_{i_1, \dots, i_d} \in \mathbb{R}^d$ are the control points defining the geometric mapping.

Following [2, 20], when the approximation spaces are refined, the control points are adjusted so that the mapping $F_{h\ell}$ remains unchanged. Thus, the geometry and its parametrization are kept fixed as the mesh is refined, i.e., when considering $\hat{\Omega}_{h\ell}$ for increasing $\ell$. Specifically,

$$F_{h\ell} = F_{h0},$$

where $F_{h\ell}$ is defined with control points $P_{i_1, \dots, i_d}^{(h\ell)}$ that are adjusted according to $h_\ell$.
We define the isogeometric analysis approximation spaces by composing the parameter spaces with the geometric mapping $\mathbf{F}_{h_0}$. This composition is also called the push-forward operation of the parameter space.

**Definition 4.4.** The physical domain sparse-grid approximation spaces are defined by

$$S_{p,h_n}^{\mathcal{H}}(\Omega) := \{u \in H_{\text{mix}}^p(\Omega) \mid u = u_{p,h_\ell}^{\mathcal{H}} \circ \mathbf{F}_{h_0}^{-1}, \quad u_{p,h_\ell}^{\mathcal{H}} \in S_{p,h_\ell}^{\mathcal{H}}(\hat{\Omega})\},$$

and

$$S_{p,h_n}^{\mathcal{L}}(\Omega) := \{u \in H_{\text{mix}}^p(\Omega) \mid u = u_{p,h_\ell}^{\mathcal{L}} \circ \mathbf{F}_{h_0}^{-1}, \quad u_{p,h_\ell}^{\mathcal{L}} \in S_{p,h_\ell}^{\mathcal{L}}(\hat{\Omega})\}.$$

Note that according to Theorem 4.1, the two spaces are identical.

### 4.2.1 Approximation Error

In this section we extend the results from the parameter approximation space to the physical domain approximation space. Before stating the main result of this paper, we introduce a lemma providing bounds for the composition of a function with the geometric mapping $\mathbf{F}_{h_0}$, or its inverse $\mathbf{F}_{h_0}^{-1}$. This result is an extension of [2], Lemma 3.5, to functions of mixed Sobolev spaces; it uses an argument similar to the result from [10], Lemma 3.

**Lemma 4.7.** For all $u \in H^{dq}(\Omega)$, with $0 \le q \le p + 1$, it holds that

$$|u|_{H^q(\Omega)} \le \sum_{m=0}^{q} C_{12} |u \circ \mathbf{F}_{h_0}|_{H^m(\hat{\Omega})}, \quad (48)$$

$$\|u \circ \mathbf{F}_{h_0}\|_{H_{\text{mix}}^q(\hat{\Omega})} \le \sum_{m=0}^{dq} C_{13} |u|_{H^m(\Omega)}, \quad (49)$$

where

$$C_{12} = C_{12}(m, q, \|\nabla \mathbf{F}_{h_0}^{-1}\|_{W^{q,\infty}(\Omega)}, \|\det \nabla \mathbf{F}_{h_0}\|_{L^{\infty}(\hat{\Omega})}),$$

$$C_{13} = C_{13}(m, q, \|\nabla \mathbf{F}_{h_0}\|_{W^{q,\infty}(\hat{\Omega})}, \|\det \nabla \mathbf{F}_{h_0}^{-1}\|_{L^{\infty}(\hat{\Omega})}),$$

are constants that depend in particular on the derivatives up to order $q$ of the geometric mapping $\mathbf{F}_{h_0}$ and its inverse $\mathbf{F}_{h_0}^{-1}$.

**Proof of Lemma 4.7.** The proof uses similar arguments than [10], Lemma 3, and is based on the multivariate Faá di Bruno decomposition formula. Let $\boldsymbol{\alpha}$ be such that $|\boldsymbol{\alpha}|_{\infty} = q$. Using the Faá di Bruno formula, the triangular inequality, the relation $\|f \cdot (\prod_j g_j)\|_{L^2} \le \|f\|_{L^2} \prod_j \|g_j\|_{L^{\infty}}$, and putting together the factors which correspond to identical values of indices of derivatives of $\mathbf{F}_{h_0}$, we obtain:

$$\begin{aligned} \|D^{\boldsymbol{\alpha}}(u \circ \mathbf{F}_{h_0})\|_{L^2(\hat{\Omega})} & \\ & \le C(q) \sum_{m=0}^{dq} \|D^m u(\mathbf{F}_{h_0})\|_{L^2(\hat{\Omega})} \sum_{i \in I(m,q,\boldsymbol{\alpha})} \prod_{j=1}^{q} \|\nabla^j \mathbf{F}_{h_0}\|_{L^{\infty}(\hat{\Omega})}^{i_j} \\ & \le C(q) \|\det \nabla \mathbf{F}_{h_0}^{-1}\|_{L^{\infty}(\hat{\Omega})}^{1/2} \sum_{m=0}^{dq} \|\nabla^m u\|_{L^2(\hat{\Omega})} \sum_{i \in I(m,q,\boldsymbol{\alpha})} \prod_{j=1}^{q} \|\nabla^j \mathbf{F}_{h_0}\|_{L^{\infty}(\hat{\Omega})}^{i_j}. \end{aligned} \quad (50)$$

In the above, we have introduced the following index set:

$$I(m, q, \boldsymbol{\alpha}) = \{i \in \mathbb{N}^q \mid |i|_1 = m, 1 \cdot i_1 + 2 \cdot i_2 + \dots + q \cdot i_q = |\boldsymbol{\alpha}|_1\}.$$
Eq. (50) is obtained with a change of variable. Then, using the Hölder inequality, and gathering all the terms in a constant $C_{12}$, we get:

$$ \|u \circ \mathbf{F}_{h_0}\|_{H_{\text{mix}}^q(\hat{\Omega})} \le \sum_{m=0}^{dq} C_{12}(m, q, \|\nabla \mathbf{F}_{h_0}\|_{W^{q,\infty}(\hat{\Omega})}, \|\det \nabla \mathbf{F}_{h_0}^{-1}\|_{L^\infty(\Omega)}) |u|_{H^m(\Omega)}. $$

We obtain Eq. (48) by applying the same arguments to $\|D^\alpha u\|_{L^2(\hat{\Omega})} = \|D^\alpha(u \circ \mathbf{F}_{h_0} \circ \mathbf{F}_{h_0}^{-1})\|_{L^2(\hat{\Omega})}$, with $|\alpha|_1 = q$. □

Note that, according to Eq. (49), in order to bound the $H_{\text{mix}}^q(\Omega)$ norm of the composition of a function $u$ with the geometric mapping $\mathbf{F}_{h_0}$, it is necessary to control the derivatives of $u$ up to order $dq$ in each direction. This arises because the mapping $\mathbf{F}_{h_0}$ is not a separable function. Although $\mathbf{F}_{h_0}$ is expressed in a tensor-product form, its control points couple the contributions of the individual basis functions, each of which is separable on its own. As a result, the mixed derivatives of the composed function $D^\alpha(u \circ \mathbf{F}_{h_0})$ with $|\alpha|_\infty = q$, appearing in the $H^q(\Omega)$ seminorm, do not vanish, thereby destroying the coordinatewise independence.

The main result of the paper is the following.

**Theorem 4.2.** Let $u \in H^{dq}(\Omega)$ with $q \in \mathbb{N}_0$. For any mesh level $n \in \mathbb{N}$, spline degree $p \in \mathbb{N}_0$, and any $r \in \mathbb{N}$ satisfying $0 \le r \le q \le p + 1$, the following estimate holds:

$$ |(I - \Pi_{p,h_n}^L) u|_{H^r(\Omega)} \le C h_n^{q-r} |\log h_n|^{d-1} \|u\|_{H^{dq}(\Omega)}, $$

where

$$ C = C_{14}(d, q, r, \|\nabla \mathbf{F}_{h_0}^{-1}\|_{W^{q,\infty}(\Omega)}, \|\nabla \mathbf{F}_{h_0}\|_{W^{q,\infty}(\hat{\Omega})}, \|\det \nabla \mathbf{F}_{h_0}\|_{L^\infty(\hat{\Omega})}, \|\det \nabla \mathbf{F}_{h_0}^{-1}\|_{L^\infty(\Omega)}) $$

is a constant independent of the mesh size $h_n$ and the spline degree $p$.

**Proof of Theorem 4.2.** For simplicity of notations, we omit the parameters dependence of the constants. The interpolation error in the $H^r(\Omega)$ seminorm is bounded by:

$$ |(I - \Pi_{p,h_n}^L) u|_{H^r(\Omega)} = \left| \left[ (I - \Pi_{p,h_n}^L) (u \circ \mathbf{F}) \right] \circ \mathbf{F}^{-1} \right|_{H^r(\Omega)} \\ \le \left( \max_{0 \le m \le dq} C_{12} \right) \sum_{m=0}^{r} \left| (I - \Pi_{p,h_n}^L) (u \circ \mathbf{F}_{h_0}) \right|_{H^m(\Omega)} \quad (51) $$

$$ \le \left( \max_{0 \le m \le dq} C_{12} \right) \sqrt{r} \|(I - \Pi_{p,h_n}^L)(u \circ \mathbf{F}_{h_0})\|_{H^r(\Omega)} \quad (52) $$

$$ \le \left( \max_{0 \le m \le dq} C_{12} \right) \sqrt{r} \|(I - \Pi_{p,h_n}^L)(u \circ \mathbf{F}_{h_0})\|_{H_{\text{mix}}^r(\Omega)} \quad (53) $$

$$ \le \left( \max_{0 \le m \le dq} C_{12} \right) C_{10} \sqrt{r} h_n^{q-r} |\log h_n|^{d-1} \|u \circ \mathbf{F}_{h_0}\|_{H_{\text{mix}}^q(\Omega)} \quad (54) $$

$$ \le \left( \max_{0 \le m \le dq} C_{12} \right) C_{10} \sqrt{r} h_n^{q-r} |\log h_n|^{d-1} \sum_{m=0}^{dq} C_{13} |u|_{H^m(\Omega)} \quad (55) $$

$$ \le \left( \max_{0 \le m \le dq} C_{12} \right) C_{10} \left( \max_{0 \le m \le dq} C_{13} \right) \sqrt{dqr} h_n^{q-r} |\log h_n|^{d-1} \|u\|_{H^{dq}(\Omega)}. \quad (56) $$

In the above, Eq. (48) of Lemma 4.7 has been used to obtain Eq. (51); Cauchy-Schwarz inequality has been used in Eq. (52); and the norm inequalities Eq. (2) has been used in Eq. (53). Then, we used Lemma 4.4 to obtain Eq. (54) and we apply Eq. (49) of Lemma 4.7 to get Eq. (55). Eventually applying again the Cauchy-Schwarz inequality we obtain the result of Eq. (56). □
Note that the result from Theorem 4.2 differs from the convergence results obtained for sparse-grid approximations defined on tensor-product domains, such as $\hat{\Omega}$. In those cases, Lemma 4.4—which can equivalently be expressed in terms of the seminorm instead of the full norm—shows that the $H_{\text{mix}}^r(\hat{\Omega})$ seminorm of the interpolation error is controlled by the $H_{\text{mix}}^q(\hat{\Omega})$ seminorm of the solution, for $0 \le r \le q$. However, for non-tensor-product domains constructed through a geometric mapping $\mathbf{F}$, controlling the $H^r(\Omega)$ seminorm of the interpolation error requires controlling the $H^{dq}(\Omega)$ norm of the solution, with $0 \le r \le q$. This implies that higher-order derivatives of the solution must be bounded—up to order $dq$ in a single direction. Consequently, employing sparse-grid tensor products to construct approximation spaces for non-tensor-product domains necessitates stronger regularity assumptions, particularly on isotropic (non-mixed) derivatives, compared to standard methods.

### 4.2.2 Inverse inequality

We extend the inverse inequality of Section 4.1.5 to the physical domain approximation space. To this end, we introduce the spaces constructed from $q$-vanishing derivative B-splines. Specifically, we define the space $\tilde{S}_{p,h_n}^{\mathcal{H},(q)}(\Omega)$ by taking the push-forward operation from the parameter space $\tilde{S}_{p,h_n}^{\mathcal{H},(q)}(\hat{\Omega})$. Thanks to [24], Corollary 1, we can prove that Theorem 4.2 is also valid for this space.

**Theorem 4.3.** Let $0 \le q \le p$, the following inverse inequality holds true:

$$|u_{p,h_n}^{\mathcal{H}}|_{H^q(\Omega)} \le C h_n^{-q} |\log h_n|^{d/2} \|u_{p,h_n}^{\mathcal{H}}\|_{L^2(\Omega)}$$

for all $u_{p,h_n}^{\mathcal{H}} \in \tilde{S}_{p,h_n}^{\mathcal{H},(q)}(\Omega)$ and where:

$$C = C_{15}(q, \|\nabla \mathbf{F}_{h_0}^{-1}\|_{W^{q,\infty}(\Omega)}, \|\nabla \mathbf{F}_{h_0}\|_{W^{q,\infty}(\hat{\Omega})}, \|\det \nabla \mathbf{F}_{h_0}\|_{L^{\infty}(\hat{\Omega})}, \|\det \nabla \mathbf{F}_{h_0}^{-1}\|_{L^{\infty}(\Omega)})$$

is a constant that does not depend on the mesh size $h_n$ or the degree $p$.

**Proof of Theorem 4.3.** For simplicity of notations, we omit the parameters dependence of the constants. Let $u_{p,h_n}^{\mathcal{H}} \in \tilde{S}_{p,h_n}^{\mathcal{H},(q)}(\Omega)$, then its $H^q(\Omega)$ seminorm error is bounded by:

$$|u_{p,h_n}^{\mathcal{H}} \circ|_{H^q(\Omega)} = |(u_{p,h_n}^{\mathcal{H}} \circ \mathbf{F}_{h_0}) \circ \mathbf{F}^{-1}|_{H^q(\Omega)} \\ \le \sum_{m=0}^{q} C_{12} |u_{p,h_n}^{\mathcal{H}} \circ \mathbf{F}_{h_0}|_{H^m(\hat{\Omega})} \qquad (57)$$

$$\le \sqrt{q} C_{12} \|u_{p,h_n}^{\mathcal{H}} \circ \mathbf{F}_{h_0}\|_{H^q(\hat{\Omega})} \qquad (58)$$

$$\le \sqrt{q} C_{12} \|u_{p,h_n}^{\mathcal{H}} \circ \mathbf{F}_{h_0}\|_{H_{\text{mix}}^q(\hat{\Omega})} \qquad (59)$$

$$\le \sqrt{q} C_{12} C_{11} h_n^{-q} |\log h_n|^{d/2} \|u_{p,h_n}^{\mathcal{H}} \circ \mathbf{F}\|_{L^2(\hat{\Omega})} \qquad (60)$$

$$\le \sqrt{q} C_{12} C_{11} C_{13} h_n^{-q} |\log h_n|^{d/2} \|u_{p,h_n}^{\mathcal{H}}\|_{L^2(\Omega)} \qquad (61)$$

In the above, Eq. (48) of Lemma 4.7 has been used to obtain Eq. (57); Cauchy-Schwarz inequality has been used in Eq. (58); and the norm inequalities Eq. (2) has been used in Eq. (59). Then, we used Lemma 4.5 to obtain Eq. (60) and eventually we apply Eq. (49) of Lemma 4.7 with $q=0$ to get Eq. (61).

□

### 4.2.3 Complexity of the Spaces

The main benefit of sparse-grid spaces is that they retain nearly the same approximation properties as standard tensor-product spaces, as demonstrated in Theorem 4.2, while significantly reducing the number of degrees of freedom. This can be formalized by the following result concerning the dimension of the sparse-grid approximation space.
**Proposition 4.1.** The dimension of the sparse-grid B-spline approximation space satisfies

$$ \dim S_{p,h_n}^{\mathcal{L}}(\hat{\Omega}) \lesssim (h_n^{-1}) |\log h_n|^{d-1}. $$

*Proof of Proposition 4.1.* The dimension of the physical domain approximation space coincides with that of the parameter domain space. To estimate this dimension, we count the number of basis functions involved in the combination:

$$ \dim S_{p,h_n}^{\mathcal{L}}(\hat{\Omega}) \le \sum_{l=0}^{d-1} |(-1)^l| \left| \begin{pmatrix} d-1 \\ l \end{pmatrix} \right| \sum_{\substack{|\ell|=n+(d-1)\lambda_p-l, \\ \ell_j \ge \lambda_p}} \dim S_{p,h_\ell}(\hat{\Omega}). \quad (62) $$

Since all of the basis functions of $S_{p,h_\ell}(\hat{\Omega})$, for a given level $\ell \in \mathbb{N}^d$, are linearly independent we have:

$$ \dim S_{p,h_\ell}(\hat{\Omega}) = \prod_{l=1}^{d} (h_{\ell l}^{-1} + p). \quad (63) $$

Substituting Eq. (63) into Eq. (62), we obtain:

$$ \dim S_{p,h_n}^{\mathcal{L}}(\hat{\Omega}) \le h_n^{-1} \sum_{l=0}^{d-1} \begin{pmatrix} d-1 \\ l \end{pmatrix} \begin{pmatrix} n+d-1-\lambda_p-l \\ d-1 \end{pmatrix} (1+p)^d 2^{-[(d-1)\lambda_p-l]}. $$

We conclude with the estimation

$$ \begin{pmatrix} n+d-1-\lambda_p-l \\ d-1 \end{pmatrix} \le \begin{pmatrix} n+d-2-l \\ d-1 \end{pmatrix} \lesssim |\log h_n|^{d-1}. $$

<figure><img src="image_9.png" /></figure>

For comparison, consider the full tensor-product space obtained by refining the univariate spline spaces uniformly in all coordinate directions. From Proposition 3.1 and Lemma 4.7, for all $u \in H^q(\Omega)$, we have

$$ \|(I - \Pi_{p,h_n}) u\|_{H^r(\Omega)} \lesssim h_n^{q-r} \|u\|_{H^q(\Omega)}, \quad (64) $$

where $\Pi_{p,h_n}$ denotes the $H^r$-orthogonal projection onto $S_{p,h_n}(\Omega)$.

When sparse-grid spaces are employed instead of full tensor-product spaces, the approximation error deteriorates only mildly—by a factor of $|\log h_n|^{d-1}$—which becomes negligible as the mesh is refined ($n \to \infty$). This mild deterioration is more than compensated by the significant reduction in the number of degrees of freedom. Indeed, the dimension of the full tensor-product space is given by

$$ \dim S_{p,h_n}(\Omega) = (h_n^{-1} + p)^d, \quad (65) $$

so the relative computational advantage of sparse grids increases rapidly with the spatial dimension $d$.

## References

[1] S. Achatz. Higher order sparse grid methods for elliptic partial differential equations with variable coefficients. *Computing*, 71(1):1–15, 2003.

[2] Y. Bazilevs, L. Beirão Da Veiga, J. A. Cottrell, T. J. R. Hughes, and G. Sangalli. Isogeometric analysis: Approximation, stability and error estimates for h-refined meshes. *Mathematical Models and Methods in Applied Sciences*, 16(07):1031–1090, 2006.
[3] J. Beck, G. Sangalli, and L. Tamellini. A sparse-grid isogeometric solver. *Computer Methods in Applied Mechanics and Engineering*, 335:128-151, 2018.

[4] L. Beirão da Veiga, A. Buffa, J. Rivas, and G. Sangalli. Some estimates for h-p-k-refinement in isogeometric analysis. *Numerische Mathematik*, 118(2):271-305, 2011.

[5] L. Beirão da Veiga, D. Cho, and G. Sangalli. Anisotropic NURBS approximation in isogeometric analysis. *Computer Methods in Applied Mechanics and Engineering*, 209-212, 2012.

[6] R. E. Bellman. *Adaptive Control Processes: A Guided Tour*. Princeton University Press, Princeton, NJ, 1961.

[7] H.-J. Bungartz. *Finite elements of higher order on sparse grids*. Berichte aus der Informatik. Shaker, 1998.

[8] H.-J. Bungartz and T. Dornseifer. Sparse grids: Recent developments for elliptic partial differential equations. In *Multigrid Methods V*, pages 45-70. Springer Berlin Heidelberg, 1998.

[9] H.-J. Bungartz and M. Griebel. Sparse grids. *Acta Numerica*, 13:147-269, 2004.

[10] P. Ciarlet and P. Raviart. Interpolation theory over curved elements, with applications to finite element methods. *Computer Methods in Applied Mechanics and Engineering*, 1(2):217-249, 1972.

[11] T. Dornseifer and C. Pflaum. Discretization of elliptic differential equations on curvilinear bounded domains with sparse grids. *Computing*, 56(3):197-213, 1996.

[12] J. Garcke. Sparse grids in a nutshell. *Lecture Notes in Computational Science and Engineering*, 88:57-80, 2013.

[13] M. Griebel. The combination technique for the sparse grid solution of PDE's on multiprocessor machines. *Parallel Processing Letters*, 02(01):61-70, 1992.

[14] M. Griebel. Adaptive sparse grid multilevel methods for elliptic PDEs based on finite differences. *Computing*, 61(2):151-179, 1998.

[15] M. Griebel and J. Hamaekers. Sparse grids for the Schrödinger equation. *ESAIM: M2AN*, 41(2):215-247, 2007.

[16] M. Griebel and H. Harbrecht. A note on the construction of L-fold sparse tensor product spaces. *Constructive Approximation*, 38(2):235-251, 2013.

[17] M. Griebel and H. Harbrecht. On the construction of sparse tensor product spaces. *Mathematics of Computation*, 82:975-994, 2013.

[18] M. Griebel and H. Harbrecht. On the convergence of the combination technique. In Jochen Garcke, editor, *Sparse Grids and Applications - Munich 2012*, pages 55-74, Cham, 2014. Springer International Publishing.

[19] M. Griebel, M. Schneider, and C. Zenger. A combination technique for the solution of sparse grid problems. *Forschungsberichte, TU Munich*, TUM I 9038:1-24, 1990.

[20] T. Hughes, J. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. *Computer Methods in Applied Mechanics and Engineering*, 194(39):4135-4195, 2005.
[21] C. Pflaum. Convergence of the combination technique for second-order elliptic differential equations. *SIAM Journal on Numerical Analysis*, 34(6):2431–2455, 1997.

[22] C. Pflaum and A. Zhou. Error analysis of the combination technique. *Numerische Mathematik*, 84(2):327–350, 1999.

[23] S. Takacs. Robust approximation error estimates and multigrid solvers for isogeometric multi-patch discretizations. *Mathematical Models and Methods in Applied Sciences*, 28(10):1899–1928, 2018.

[24] S. Takacs and Takacs T. Approximation error estimates and inverse inequalities for B-splines of maximum smoothness. *Mathematical Models and Methods in Applied Sciences*, 26:1411–1445, 2016.

[25] B. Volker, N. Erich, and R. Klaus. High dimensional polynomial interpolation on sparse grids. *Advances in Computational Mathematics*, 12(4):273–288, 2000.
Inria

RESEARCH CENTRE
Centre Inria de l'Université de Bordeaux

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399