A Greedy Method for Solving Classes of PDE Problems

TL;DR

This paper introduces a greedy algorithm for efficiently solving classes of PDE problems by selecting optimal data in a Hilbert space, with convergence properties linked to Kolmogoroff N-widths and demonstrated through numerical experiments.

Contribution

It proposes a novel greedy method for PDEs that directly operates in the dual Hilbert space and connects its convergence to Kolmogoroff N-widths, extending greedy algorithms to PDE applications.

Findings

The greedy method performs comparably to P-greedy interpolation for elliptic PDEs.

Numerical results suggest the method achieves near-optimal approximation rates.

The approach is theoretically grounded in Hilbert space properties and N-widths.

Abstract

Motivated by the successful use of greedy algorithms for Reduced Basis Methods, a greedy method is proposed that selects N input data in an asymptotically optimal way to solve well-posed operator equations using these N data. The operator equations are defined as infinitely many equations given via a compact set of functionals in the dual of an underlying Hilbert space, and then the greedy algorithm, defined directly in the dual Hilbert space, selects N functionals step by step. When N functionals are selected, the operator equation is numerically solved by projection onto the span of the Riesz representers of the functionals. Orthonormalizing these yields useful Reduced Basis functions. By recent results on greedy methods in Hilbert spaces, the convergence rate is asymptotically given by Kolmogoroff N-widths and therefore optimal in that sense. However, these N-widths seem to be unknown in PDE applications. Numerical experiments show that for solving elliptic second-order Dirichlet problems, the greedy method of this paper behaves like the known P-greedy method for interpolation, applied to second derivatives. Since the latter technique is known to realize Kolmogoroff N-widths for interpolation, it is hypothesized that the Kolmogoroff N-widths for solving second-order PDEs behave like the Kolmogoroff N-widths for second derivatives, but this is an open theoretical problem.

Equations

$$\begin{array}[]{rcl}Lu&=&f\;\hbox{ on }\overline{\Omega}\subset\mathbb{R}^{d}\\ u&=&g\;\hbox{ on }\Gamma:=\partial\Omega\end{array}$$

$$P_{{\Lambda }_N}(\lambda ):=\mathrm{dist}(\lambda ,span({\Lambda }_N))$$

$${\sigma }_{{\Lambda }_N}:=\underset{\lambda \in \Lambda }{\max }{P}_{{\Lambda }_N}(\lambda )$$

$${\lambda }_{N+1}:=arg\; max{P}_{{\Lambda }_N}(\lambda )$$

$${\rho }_{{\Lambda }_N}:=\underset{x\in \overline{\Omega }}{\max }{P}_{{\Lambda }_N}({\delta }_x)$$

$$\begin{array}[]{rcl}\Lambda_{1}&:=&\{\delta_{x}\circ L\;:\;x\in\overline{\Omega}\}\\ \Lambda_{2}&:=&\{\delta_{x}\;:\;x\in\Gamma:=\partial{\Omega}\}\end{array}$$

$$\Vert u{\Vert }_{\infty ,\overline{\Omega }}\le \Vert u{\Vert }_{\infty ,\partial \Omega }+C\Vert Lu{\Vert }_{\infty ,\overline{\Omega }}\le (C+1)\underset{\lambda \in \Lambda }{\sup }|\lambda (u)|\; for\; allu\in C^2(\overline{\Omega })\cup C(\Gamma ).$$

$$\lambda (u)=0\; for\; all\lambda \in \Lambda \; impliesu=0.$$

$$D_{\Lambda }(u):=\{\lambda (u){\}}_{\lambda \in \Lambda }for\; allu\in \mathcal{H},$$

$$\Vert u{\Vert }_{WP}\le C_{WP}\Vert D_{\Lambda }(u){\Vert }_{\infty }=C_{WP}\underset{\lambda \in \Lambda }{\sup }|\lambda (u)|\; for\; allu\in \mathcal{H}$$

$$L_N:=span\{{\Lambda }_N\}\subseteq {\mathcal{H}}^{\ast }$$

$${\lambda }_k(u)=\sum _{j=1}^N{\alpha }_j{\lambda }_k(v_{{\lambda }_j})=\sum _{j=1}^N{\alpha }_j(v_{{\lambda }_k},v_{{\lambda }_j}{)}_{\mathcal{H}}=\sum _{j=1}^N{\alpha }_j({\lambda }_k,{\lambda }_j{)}_{{\mathcal{H}}^{\ast }},1\le k\le N$$

$${\tilde{u}}_N:=\sum _{j=1}^N{\alpha }_j{v}_{{\lambda }_j}.$$

$$\Vert u{\Vert }_{\mathcal{H}}^2=\Vert u-{\tilde{u}}_N{\Vert }_{\mathcal{H}}^2+\Vert {\tilde{u}}_N{\Vert }_{\mathcal{H}}^2$$

$$P_{{\Lambda }_N}(\mu ):=\underset{\lambda \in L_N}{\min }\Vert \mu -\lambda {\Vert }_{{\mathcal{H}}^{\ast }}=:\mathrm{dist}(\mu ,L_N{)}_{{\mathcal{H}}^{\ast }}for\; all\mu \in {\mathcal{H}}^{\ast }.$$

$${\sigma }_{{\Lambda }_N}:=\underset{\lambda \in \Lambda }{\sup }{P}_{{\Lambda }_N}(\lambda )=\underset{\lambda \in \Lambda }{\sup }\mathrm{dist}(\lambda ,L_N{)}_{{\mathcal{H}}^{\ast }}.$$

$$\begin{array}[]{rcl}\|u-\tilde{u}_{N}\|_{WP}&\leq&C_{WP}\sigma_{\Lambda_{N}}\|u\|_{\cal H}.\end{array}$$

$$\begin{array}[]{rcl}\|u-\tilde{u}_{N}\|_{WP}&\leq&C_{WP}\sup_{\lambda\in\Lambda}|\lambda(u)-\lambda(\tilde{u}_{N})|\\ &=&C_{WP}\sup_{\lambda\in\Lambda}|(\lambda-\mu)(u-\tilde{u}_{N})|\hbox{ for all }\mu\in L_{N}\end{array}$$

$$|\mu (u)-\mu ({\tilde{u}}_N)|\le P_{{\Lambda }_N}(\mu )\Vert u{\Vert }_{\mathcal{H}},$$

$$\Vert D_{\Lambda }(u-{\tilde{u}}_N){\Vert }_{\infty }=\underset{\lambda \in \Lambda }{\sup }|\lambda (u-{\tilde{u}}_N)|\le {\sigma }_{{\Lambda }_N}\Vert u{\Vert }_{\mathcal{H}},$$

$$\begin{array}[]{rcl}|u(x)-\tilde{u}_{N}(x)|&\leq&P_{\Lambda_{N}}(\delta_{x})\|u\|_{\cal H}\end{array}$$

$${\rho }_{{\Lambda }_N}:=\underset{x\in \overline{\Omega }}{\max }{P}_{{\Lambda }_N}({\delta }_x)$$

$$\Vert u-{\tilde{u}}_N{\Vert }_{\infty ,\overline{\Omega }}\le {\rho }_{{\Lambda }_N}\Vert u{\Vert }_{\mathcal{H}}$$

$$\Vert u{\Vert }_{WP}=\underset{\mu \in M}{\sup }|\mu (u)|,u\in \mathcal{H}$$

$${\rho }_{M,{\Lambda }_N}:=\underset{\mu \in M}{\max }{P}_{{\Lambda }_N}(\mu ).$$

$${\rho }_{M,{\Lambda }_N}\le C_{WP}{\sigma }_{{\Lambda }_N}.$$

$$\begin{array}[]{rcl}\rho_{M,\Lambda_{N}}&=&P_{\Lambda_{N}}(\tilde{\mu}_{N})\\ &=&\mathrm{dist}(\tilde{\mu}_{N},\hfil\lx@intercol\hbox{ span }\Lambda_{N})\\ &=&\displaystyle{\left\|\tilde{\mu}_{N}-\sum_{j=1}^{N}(\tilde{\mu}_{N},\mu_{j})\mu_{j}\right\|_{{\cal H}^{*}}}.\end{array}$$

$$\begin{array}[]{rcl}\left|\tilde{\mu}_{N}(u)-\sum_{j=1}^{N}(\tilde{\mu}_{N},\mu_{j})\mu_{j}(u)\right|&=&\left|\tilde{\mu}_{N}(u)-\tilde{\mu}_{N}(\tilde{u}_{N})\right|\\ &\leq&\sup_{\mu\in M}\left|\mu(u)-\mu(\tilde{u}_{N})\right|\\ &=&\|u-\tilde{u}_{N}\|_{WP}\\ &\leq&C_{WP}\sup_{\lambda\in\Lambda}|\lambda(u)-\lambda(\tilde{u}_{N})|\\ &\leq&C_{WP}\sigma_{\Lambda_{N}}\|u\|_{\cal H}\end{array}$$

$$P_{{\Lambda }_N}({\lambda }_{N+1})=\underset{\lambda \in \Lambda }{\sup }{P}_{{\Lambda }_N}(\lambda ).$$

$$d_N(\Lambda ):=\underset{all{H}_N}{\inf }\underset{\mu \in \Lambda }{\sup }\underset{\lambda \in H_N}{\inf }\Vert \lambda -\mu {\Vert }_{{\mathcal{H}}^{\ast }}=\underset{all{H}_N}{\inf }\underset{\mu \in \Lambda }{\sup }\mathrm{dist}(\mu ,H_N{)}_{{\mathcal{H}}^{\ast }}$$

Full text

A Greedy Method for Solving Classes of PDE Problems

Robert Schaback

Draft of

Motivated by the successful use of greedy algorithms for Reduced Basis Methods, a greedy method is proposed that selects $N$ input data in an asymptotically optimal way to solve well-posed operator equations using these $N$ data. The operator equations are defined as infinitely many equations given via a compact set of functionals in the dual of an underlying Hilbert space, and then the greedy algorithm, defined directly in the dual Hilbert space, selects $N$ functionals step by step. When $N$ functionals are selected, the operator equation is numerically solved by projection onto the span of the Riesz representers of the functionals. Orthonormalizing these yields useful Reduced Basis functions. By recent results on greedy methods in Hilbert spaces, the convergence rate is asymptotically given by Kolmogoroff $N$-widths and therefore optimal in that sense. However, these $N$-widths seem to be unknown in PDE applications. Numerical experiments show that for solving elliptic second-order Dirichlet problems, the greedy method of this paper behaves like the known $P$-greedy method for interpolation, applied to second derivatives. Since the latter technique is known to realize Kolmogoroff $N$-widths for interpolation, it is hypothesized that the Kolmogoroff $N$-widths for solving second-order PDEs behave like the Kolmogoroff $N$-widths for second derivatives, but this is an open theoretical problem.

Keywords: Operator equations, Greedy methods, Reduced Basis Methods, Kolmogoroff $N$-widths, Partial Differential Equations, meshless methods, collocation, discretization, error bounds, well-posedness, stability, convergence

Mathematics Subject Classification (2000): 65M12, 65M70, 65N12, 65N35, 65M15, 65M22, 65J10, 65J20, 35D30, 35D35, 35B65, 41A25, 41A63

1 Introduction

For illustration of the application scope of this paper, consider the class of all second-order elliptic boundary value problems

[TABLE]

with arbitrary Dirichlet data and a fixed second-order strongly elliptic operator $L$ on a fixed bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$. We keep this problem class in strong form and pose it in Sobolev space $W_{2}^{m}(\mathbb{R}^{d})$ with $m>2+d/2$ and spaces of data functions $f$ and $g$ of corresponding smoothness. Similar to Reduced Basis (e.g. [14, 1, 3, 6]) or Proper Orthogonal Decomposition methods (e.g. [15, 13, 10]), we focus on a class of PDE problems, not on single problems. The output of this paper will be connected to both areas, since a “reduced” orthonormal basis is produced that is adapted to the given class of PDE problems.

The next section will generalize such problems to operator equations defined by sets $\Lambda$ of infinitely many functionals on an underlying Hilbert space of functions, e.g. a Sobolev space. In this context, well-posedness can be formulated, and for $N$ selected functionals in a set $\Lambda_{N}\subset\Lambda$, numerical solutions can be obtained by Hilbert space projection on the Riesz representers of these functionals. This is the well-known Rayleigh-Ritz idea. For kernel-based spaces, it coincides with Symmetric Collocation and yields the optimal recovery technique in Hilbert space for the given data functionals.

Section 3 analyzes the error in terms of the Generalized Power Function

[TABLE]

on the dual of the Hilbert space and introduces

[TABLE]

that controls the error of the numerical solution in Hilbert space.

The greedy method of Section 4 now selects

[TABLE]

and follows the results of the literature on Reduced Basis Methods, linking the decay of $\sigma_{\Lambda_{N}}$ to Kolmogoroff $N$-widths.

For operator equations defined on Hilbert spaces of functions on a bounded domain $\Omega$, the additional quantity

[TABLE]

directly controls the pointwise and uniform error in the domain, but is not useful for greedy methods. For well-posed problems, Section 3.2 will show that $\rho_{\Lambda_{N}}$ is bounded form above by $\sigma_{\Lambda_{N}}$ up to a factor, and numerical results in the final Section 6 suggest that this bound is asymptotically sharp.

Before that, Section 5 gives a partial analysis of expectable Kolmogoroff $N$-widths for second-order Dirichlet problems in Sobolev spaces $W_{2}^{m}$. It is hypothesized that both $\sigma_{\Lambda_{N}}$ and the $N$-width behave like ${\cal O}(N^{-\frac{m-2-d/2}{d}})$ for $N\to\infty$, which is the Kolmogoroff $N$-width for interpolation problems in $W_{2}^{m-2}$ with respect to the supremum norm. Section 6 shows supporting examples and concludes the paper.

2 Hilbert Space Theory

Following [11, 21] the problem class is written in terms of infinitely many constraints, each defined by a linear functional. The functional sets then are

[TABLE]

combined into $\Lambda:=\Lambda_{1}\cup\Lambda_{2}$. Since all single functionals are continuous on $W_{2}^{m}(\mathbb{R}^{d})$, the above sets are images of compact sets by continuous maps, thus compact. This brings us into line with the literature on reduced basis methods.

The above problems are well-posed in the sense that there is a standard well-posedness inequality of the form [2, 1.5, p. 30]

[TABLE]

2.1 Abstract Problem

Generalizing this case, and following [11, 21], we assume a Hilbert space ${\cal H}$ and a subset $\Lambda$ of its dual ${\cal H}^{*}$ that is total in the sense that

[TABLE]

If we formally introduce the linear data map $D_{\Lambda}\;:{\cal H}\to\mathbb{R}^{\Lambda}$ with

[TABLE]

this means that elements $u\in{\cal H}$ are uniquely identifiable from their data $D_{\Lambda}(u)$. The central background problem in this paper is to recover elements $u$ from their data $D_{\Lambda}(u)$ in practice, i.e. the approximate numerical inversion of the data map. In view of the preceding example we assume that the set $\Lambda$ is compact.

The invertibility of the data map is quantified by assuming a well-posedness inequality

[TABLE]

in some well-posedness norm $\|.\|_{WP}$ on ${\cal H}$ that usually is weaker than the norm on ${\cal H}$. Recall that [21] allows also to handle weakly formulated problems as well this way. Furthermore, the framework applies to general operator equations, including the case of interpolation if the operator is the identity.

3 Error Analysis of Projection Methods

For a finite subset $\Lambda_{N}:=\{\lambda_{1},\ldots,\lambda_{N}\}$ of $\Lambda$ we can define the subspace

[TABLE]

and use the Riesz representers $v_{\lambda_{1}},\ldots,v_{\lambda_{N}}$ of $\lambda_{1},\ldots,\lambda_{N}$ as trial functions. They span a space $V_{\Lambda_{N}}$ that will lead later to reduced bases.

The standard optimal recovery of an element $u\in{\cal H}$ from finite data $\lambda_{1}(u),\ldots,\lambda_{N}(u)$ then proceeds by Hilbert space projection, i.e. by solving the linear system

[TABLE]

to get the numerical approximation

[TABLE]

This satisfies the orthogonality relation

[TABLE]

that implies uniform stability in Hilbert space. In case of the example in the beginning, the method is known as Symmetric Collocation. This is a numerical technique [7] based on [23] with certain optimality properties [20] and a convergence theory [9, 8]. Within the Hilbert space framework, it produces pointwise optimal approximations to the true solution under all possible methods that use the same data [20].

3.1 Power Function

The standard error analysis in Hilbert Spaces uses the generalized Power Function defined as

[TABLE]

This is continuous and attains its maximum on the compact set $\Lambda$. In particular, we are interested in

[TABLE]

3.2 Error Analysis

This quantity leads to error bounds for the discretized recovery problem that we described in the beginning of this section.

Lemma 1**.**

Let $u\in{\cal H}$ supply the finite data $\lambda_{1}(u),\ldots,\lambda_{N}(u)$ that is used to construct $\tilde{u}_{N}$ by projection, and assume well-posedness in the sense of (4). Then

[TABLE]

Proof.

Since $u$ and $\tilde{u}_{N}$ share the same data, the assertion follows from

[TABLE]

using orthogonality (5) and well-posedness (4). ∎

Therefore we are interested to find sets $\Lambda_{N}$ that minimize $\sigma_{\Lambda_{N}}$ under all sets of $N$ functionals from $\Lambda$.

By the same argument, for all continuous test functionals $\mu$ the inequality

[TABLE]

holds and we can check derivative errors once we evaluate the Generalized Power Function on derivative functionals. Taking all functionals of $\Lambda$ here, we get

[TABLE]

i.e. $\sigma_{\Lambda_{N}}$ bounds the error in the data norm.

In case of a Hilbert space ${\cal H}$ containing continuous functions on some compact set $\overline{\Omega}$, there is a special case of (9) in the pointwise form

[TABLE]

that may be checked by evaluation of the Generalized Power Function over all delta functionals for all points of the domain. If we define

[TABLE]

we get a uniform bound

[TABLE]

that is numerically available, and similar to the well-posedness inequality (3).

By definition, the Generalized Power Function decreases at all functionals when we extend the set $\Lambda_{N}$. All of these error bounds will then improve. For fixed $\Lambda_{N}$ and ${\cal H}$, they are optimal, but here we are interested in finding good finite subsets $\Lambda_{N}$ of $\Lambda$.

Note that both Lemma 1 and (12) furnish similar error bounds with associated convergence rates, and they might have the same behavior if the well-posedness norm $\|.\|_{WP}$ coincides with $\|.\|_{L_{\infty}(\overline{\Omega})}$. We shall have a closer look at this now, delaying experiments to Section 6.

To find a connection between $\rho_{\Lambda_{N}}$ and $\sigma_{\Lambda_{N}}$ in well-posed situations, we assume that the well-posedness norm $\|.\|_{WP}$ of (4) can be written as

[TABLE]

with a compact set $M\subset{\cal H}^{*}$. This yields (3) when taking $M$ as the set of delta functionals on $\overline{\Omega}$. The generalization of (11) then is

[TABLE]

Lemma 2**.**

Assuming (4) and (13), we have

[TABLE]

Proof.

The maximum in (14) is attained at some $\tilde{\mu}_{N}\in M$ and then

[TABLE]

Applied to an arbitrary $u\in{\cal H}$ this yields

[TABLE]

where the last line follows like (9) using a generalization of the argument within (10). ∎

4 Greedy Method

Given a set $\Lambda_{N}:=\{\lambda_{1},\ldots,\lambda_{N}\}$, the $P$-greedy algorithm [5, 18] in its abstract form defines $\lambda_{N+1}$ recursively via

[TABLE]

In the context of Reduced Basis Methods in Hilbert spaces, this algorithm coincides with the one in [14, 1, 3, 6]. The main difference is that we work in the dual here, focusing on a single class of PDEs in applications instead of a parametrized family. A first application of Reduced Basis Methods within a Reproducing Kernel Hilbert Space setting is in [4], implementing a greedy method based on discrete least-squares, not on the dual of the basic Hilbert space of functions.

The cited papers provide a useful error analysis of the method that we state here for completeness. The main ingredient is the Kolmogorov $N$-width

[TABLE]

where the first infimum is taken over all $N$-dimensional subspaces $H_{N}\subseteq{\cal H}^{*}$, not only those spanned by $N$ functionals from $\Lambda$. Then [6] proves

[TABLE]

while [1] has

[TABLE]

with suitable constants. This links the behaviour of the $P$-greedy algorithm to Kolmogoroff $N$-widths. The paper [18] exploits this connection for function recovery by interpolation, while this paper extends [18] to classes of operator equations, including PDE solving.

Consequently, the greedy method converges roughly like the Kolmogoroff $N$-widths. Asymptotically, there are no better choices for selecting $N$ functionals out of $\Lambda$. Note that this works for all well-posed problems stated in abstract form in Hilbert space via infinitely many constraints. The differential and boundary operators can be easily generalized.

On the downside, the literature does not provide much information about the Kolmogoroff $N$-widths in such situations. This is why on has to look at special cases. We postpone this to Section 5.

If we run the greedy algorithm numerically, we should get very good candidates for reduced bases. Choosing them with additional orthogonality properties will then bring us close to Proper Orthogonal Decomposition methods. The next section will explain how to do that.

4.1 Implementation

Throughout, we assume that the Hilbert space ${\cal H}$ has a kernel $K$ such that inner products in ${\cal H}^{*}$ can be numerically calculated via

[TABLE]

where the upper index stands for the variable the functional acts on. In case of Sobolev spaces $W_{2}^{m}(\mathbb{R}^{d})$ with $m>d/2$, we use the standard Whittle-Matérn kernel

[TABLE]

with the modified Bessel function $K_{\nu}$ of second kind. In what follows, we treat $\Lambda$ as being very large and finite, but extensions to infinite compact $\Lambda$ will be possible if functions on $\Lambda$ are discretized somehow.

A direct way to assess the generalized Power Function at $\lambda$ for given $\Lambda_{N}:=\{\lambda_{1},\ldots,\lambda_{N}\}$ is to use its definition via the approximation problem

[TABLE]

We assume that we have turned $\lambda_{1},\ldots,\lambda_{N}$ into an ${\cal H}^{*}$-orthonormal basis $\mu_{1},\ldots,\mu_{N}$ already, and then the solution is

[TABLE]

Furthermore, we store the orthonormalization system for getting the $\mu_{j}$ from the $\lambda_{j}$ as

[TABLE]

Besides the $N\times N$ triangular matrix $C$, we store the values $\{(\lambda,\lambda)_{{\cal H}^{*}}\}_{\lambda\in\Lambda}$ and $\{(\lambda,\mu_{k})_{{\cal H}^{*}}\}_{\lambda\in\Lambda}$ for $k=1,\ldots,N$. These make up the method’s bulk storage of order $(N+1)|\Lambda|$. Note that $\{(\lambda,\lambda)_{{\cal H}^{*}}\}_{\lambda\in\Lambda}$ simplifies considerably for translation-invariant kernels.

After maximizing (16) over all $\lambda\in\Lambda$, we assume to have some $\lambda_{N+1}$ at which a nonzero maximum is attained. This can then not be one of the old $\lambda_{j}$, and we retrieve the values $(\lambda_{N+1},\mu_{k})_{{\cal H}^{*}}$ from what we have. We now orthonormalize

[TABLE]

via $c_{N+1,j}=-(\lambda_{N+1},\mu_{k})_{{\cal H}^{*}}c_{N+1,N+1},\;1\leq k\leq N$ and

[TABLE]

to update the $C$ matrix. Finally

[TABLE]

can be calculated from what we have, if we first calculate all $\{(\lambda,\lambda_{N+1})_{{\cal H}^{*}}\}_{\lambda\in\Lambda}$ and overwrite them with $\{(\lambda,\mu_{N+1})_{{\cal H}^{*}}\}_{\lambda\in\Lambda}$ after use. This extends the Newton basis technique in [16] to general functionals.

4.2 Bases and Postprocessing

Useful reduced bases for PDE solving are the Riesz representers $v_{\mu_{k}}$ of the $\mu_{k}$, being orthonormal in Hilbert space. From Riesz representers $v_{\lambda_{j}}(\cdot)=\lambda_{j}^{x}K(x,\cdot)$ of the $\lambda_{j}$, we can calculate them recursively via

[TABLE]

on whatever point sets we like, using our triangular matrix $C$. Via (16) in the form

[TABLE]

this allows to calculate the pointwise error bounds (10) described in Section 3.2 explicitly, up to the term $\|u\|_{\cal H}$.

In view of Proper Orthogonal Decomposition, one can also apply a Singular Value Decomposition to the partial Gramian matrix with entries $(\lambda_{j},\lambda_{k})_{{\cal H}^{*}},\;1\leq j,k\leq N$ and construct a different ${\cal H}$-orthonormal basis. This possibility is not pursued here, because we want to keep the recursive structure of the algorithm. In view of [17], this basis may be closer to what happens for Kolmogoroff $N$-widths, because the spaces for the latter will not necessarily have a recursive structure.

Like in reduced basis techniques, the orthonormal bases are a simple tool to solve a variety of similar problems, namely all problems of the form (1). If data $\lambda_{j}(u)$ are known for an unknown function $u$, we go over to $\mu_{j}(u)$ via the triangular system (17) and then form the optimal projection

[TABLE]

using the ${\cal H}$-orthonormal basis we have constructed.

Note that this algorithm neither stores nor solves a large $|\Lambda|\times|\Lambda|$ system. It works “on–the–fly”. For $N$ steps,  storage is of order ${\cal O}(N\cdot|\Lambda|+N^{2})$ and calculations are of order ${\cal O}(N^{2}\cdot|\Lambda|)$. For a given accuracy requirement, the number $N$ of steps will depend on the Kolmogoroff $N$-width for the set $\Lambda$ and the Hilbert space ${\cal H}^{*}$.

There are no square systems to be solved. Instead, there are $N$ orthonormalization steps on vectors of length $|\Lambda|$ that may require standard stabilization precautions for the basic Gram-Schmidt technique. The resulting triangular $N\times N$ matrix $C$ is not explicitly inverted, but used via (17) to transform input data in terms of the $\lambda_{j}$ functionals into data in terms of the orthonormalized functionals $\mu_{j}$. The error behavior of this is comparable to backsubstitution after an $LR$ or $QR$ factorization, but it will pay the price when some of the $\lambda_{j}$ are strongly correlated. This will be unavoidable for large $N$, but (6) and (15) show that $\lambda_{N+1}$ will be kept away from the zeros $\lambda_{1},\ldots,\lambda_{N}$ of $P_{\Lambda_{N}}$ and the space they span, by construction.

4.3 Extended Greedy Method

One can get somewhat closer to the error analysis in Section 3.2 and in particular to (9) by a modification of the selection strategy of functionals. Given a set $\Lambda_{N}:=\{\lambda_{1},\ldots,\lambda_{N}\}$, and a set $M=\{\mu_{1},\ldots,\mu_{M}\}$ with nonempty intersection, we calculate two maxima

[TABLE]

and set

[TABLE]

This can be called an extended $P$-greedy algorithm. It tries to keep some additional control of $\rho_{M,\Lambda_{N}}$ from (14) by selecting functionals from $\Lambda\cap M$ whenever $\rho_{M,\Lambda_{N}}$ is attained on them.

For solving Dirichlet problems, the set $M$ will consist of delta functionals in $\overline{\Omega}$, the intersection of $M$ and $\Lambda$ being $\Lambda_{2}$, the set of functionals for Dirichlet boundary values. The new technique will make sure that if $P_{\Lambda_{N}}(M)$ attains its maximum on the boundary, the corresponding functional is preferred over the functional where $P_{\lambda_{N}}(\Lambda)$ attains its maximum.

Now for some inplementation details for the Dirichlet case. Define the set $Z$ of boundary points via

[TABLE]

and add some other point set $Y\subset\overline{\Omega}$ to get

[TABLE]

We need the additional value $\|P_{\Lambda_{N}}\|_{\infty,Y}$, while the usual Greedy method provides $\|P_{\Lambda_{N}}\|_{\infty,Z}$ as part of the calculation of $\rho_{\Lambda_{N}}$. If we use (19) on $Y$, we have the necessary data, but for that we have to evaluate (18) on $Y$ as well, and have to store the values of the orthonormal basis on $Y$. This requires an additional storage of size $|Y|\cdot N$ for $N$ steps, and additional calculations of order $|Y|\cdot N^{2}$. On can choose $|Y|$ smaller than $|\Lambda_{1}|$ to keep the complexity at bay.

5 Sobolev case

We now go back to the example at the beginning. We need information on the Kolmogoroff $N$-width

[TABLE]

where the infimum is taken over all $N$-dimensional subspaces $H_{N}$ of ${\cal H}^{*}={W_{2}^{m}(\Omega)}^{*}$, and the set $\Lambda=\Lambda_{1}\cup\Lambda_{2}\subset{\cal H}^{*}$ is formed by (2). There should be an $m$- and $d$-dependent decay rate $\kappa(m,d)$ in the sense

[TABLE]

but no explicit results on this were found yet.

If we restrict attention to spaces $H_{N}$ being generated by $N$ functionals from $\Lambda=\Lambda_{1}\cup\Lambda_{2}\subset{\cal H}^{*}$, we only get an upper bound for $d_{N}(\Lambda)$, and we do not know the splitting $N=N_{1}+N_{2}$ if an optimal choice of $N$ functionals from $\Lambda$ takes $N_{1}$ functionals out of $\Lambda_{1}$ and $N_{2}$ functionals out of $\Lambda_{2}$. The rest of the chapter will give some arguments supporting the hypothesis

[TABLE]

that will be observed in the numerical behavior of the $P$-greedy method in Section 6.

If $\Lambda$ is the union of two disjoint compact sets $\Lambda_{1}$ and $\Lambda_{2}$, then

[TABLE]

where $H_{N_{i}}$ is used for approximation of $\Lambda_{i}$, and the sum of spaces is direct. This splitting argument is similar to the technique in [8]. However, we do not know a priori how the greedy algorithm selects functionals from either set, and how the dimension splits into $N_{1}+N_{2}\leq N$.

We first aim at $d_{N}(\Lambda_{1})$ in the space ${W_{2}^{m}(\Omega)}^{*}$. We have

[TABLE]

and can majorise it by choosing $N$ asymptotically uniformly placed points $x_{1},\ldots x_{N}$ in $\overline{\Omega}$ at fill distance $h_{\Omega}$ and taking $H_{N}$ to be the span of the corresponding functionals $\delta_{x_{j}}\circ L$. Then

[TABLE]

holds and we can expect

[TABLE]

due to standard results on error bounds for interpolation [22], and with generic constants depending on $m$, $d$ and the domain.

On the boundary $\Gamma$, we can argue similarly to expect

[TABLE]

for a fill distance $h_{\Gamma}$ on the boundary, either by working in ${W_{2}^{m}(\Omega)}^{*}$ directly, or via trace theorems, which would give the same rate due to $m-1/2-(d-1)/2=m-d/2$.

If we now consider splittings $N=N_{1}+N_{2}$, we roughly have the upper bound

[TABLE]

and the crude split $N_{1}\approx N/2\approx N_{2}$ will already support (21).

For purposes of asymptotics, we can minimize the sum of the above quantities instead of the maximum. Using standard optimization arguments under the constraint $N_{1}+N_{2}\leq N$, the result after some calculations is that one should expect

[TABLE]

with constants depending on the domain and the space dimension, but not on $N$ and $m$. For the typical case $m=4,\;d=2$ this implies $N_{2}\approx N_{1}^{3/8}$. For large $m$ and $d$ we get $N_{1}\approx N_{2}$, i.e. the domain and the boundary require roughly the same degrees of freedom. This is no miracle, because the volumes of balls of high dimension are increasingly accounted for by the boundary layer. In general, the above argument implies $N_{2}\leq c\,N_{1}$ for $N$ large enough, and then (21) will hold. Section 6 will reveal that the greedy algorithm selects unexpectedly small values of $N_{2}$, and that the hypothesis (21) is supported.

6 Numerical Results

No matter which PDE examples are selected, the most interesting question is the behavior of the greedy method as a function of the $N$ steps it takes. If $\Lambda$ is chosen large enough but still finite, the crucial quantities are $\sigma_{\Lambda_{N}}$ and $\rho_{\Lambda_{N}}$ as defined in (7) and (11). By Section 4 we can expect that $\sigma_{\Lambda_{N}}$ decays at an optimal rate comparable to the Kolmogoroff $N$-width with respect to $\Lambda$, but since the latter is still unknown, we can not yet assess how close we come to it. Since $\rho_{\Lambda_{N}}$ is controlling the error in the sup norm, without being usable for greedy refinement, and in view of Lemma 2, we would like to confirm that $\rho_{\Lambda_{N}}$ decays as fast as $\sigma_{\Lambda_{N}}$. And, the decay rates should improve with smoothness of the functions in the basic Hilbert space, i.e. with $m$ if we work in $W_{2}^{m}(\mathbb{R}^{d})$, our hypothesis being (21). Another interesting question is how the greedy method chooses between boundary functionals from $\Lambda_{2}$ and domain functionals from $\Lambda_{1}$, and whether the corresponding points are roughly uniformly distributed in both cases. Finally, the shape and the behaviour of the basis functions $v_{\mu_{j}}$ should be demonstrated.

6.1 Observations for varying $N$

We start on the 2D unit disk, with $L$ being the Laplace operator, carrying the greedy method out for up to 500 steps, offering 17570 functionals for $\Lambda_{1}$ and 150 functionals for $\Lambda_{2}$. The Hilbert space will be $W_{2}^{m}(\mathbb{R}^{2})$, but we fix $m>2+d/2$ first, to study the behaviour of the greedy method for varying the number $N$ of steps. Classical results on kernel-based interpolation lets us expect rates for $\sigma_{\Lambda_{N}}$ and $\rho_{\Lambda_{N}}$ that are determined by fill distances. If $h_{\Gamma}$ and $h_{\Omega}$ are fill distances for points on the boundary $\Gamma$ and in the domain $\Omega$, we can compare $\sigma_{\Lambda_{N}}$ with plain interpolation of $\Delta u=f$ with smoothness $m-2$ and behavior

[TABLE]

on the domain. The error on the boundary in $L_{\infty}$ should, if it were a plain interpolation, behave like

[TABLE]

if $N_{\Gamma}$ points are asymptotically equally spaced on the boundary. But we do not know a priori how the greedy algorithm chooses between boundary and domain functionals.

The first experiments are for $m=4>2+d/2=3$, and we ignore the split $N=N_{\Omega}+N_{\Gamma}$ in the beginning. The scale of the Whittle-Matérn kernel is chosen to be 1, and then both types of functionals happen to have the same norm. Scaling changes the relation between function value evaluation and Laplace operator evaluation and must be used with care. We wanted to eliminate scaling effects for what follows, in order to let the choice between boundary and domain functionals be unbiased by possibly different norms.

Figure 1 shows $\sigma_{\Lambda_{N}}$ (left) and $\rho_{\Lambda_{N}}$ (right) as functions of $N$, with the observed rates $-0.45$ and $-0.54$, respectively. We see that $\rho_{\Lambda_{N}}$ decays as fast as $\sigma_{\Lambda_{N}}$ as functions of $N$.

The strange drops of $\rho_{\Lambda_{N}}$ at certain $N$ occur exactly when the fill density $h_{\Gamma}$ on the boundary drops, namely after the greedy method has chosen another boundary functional. This can be read off the two plots of Figure 2 that show both $\sigma_{\Lambda_{N}}$ and $\rho_{\Lambda_{N}}$ (left) and both $h_{\Gamma}(N)$ and $h_{\Omega}(N)$ as functions of the $N$ functionals that were greedily selected. Figure 4 shows the strange fact that the greedy method selects only rather few boundary functionals compared to domain functionals (6 versus 494). It allows two large peaks on the boundary, because it still fights for getting small on the domain functionals. These effects were observed in many other cases that we suppress here for brevity. Choosing different weights for domain and boundary functionals makes some sense in view of well-posedness inequalities like (3), and then the effect will be less apparent.

Figure 3 shows decays like

[TABLE]

in comparison to the expectations in (23) and (22) that suggest rates $-0.5$ and $-3$, respectively, for $m=4$. This shows again that the greedy method focuses on domain points, and is able to maintain a small error on the boundary by adding “domain” functionals there and close to the boundary.

If one looks at the Power Function on $\delta$ functionals all over the domain, like for the $\rho_{\Lambda_{N}}$ calculation, the results are in Figure 5. It shows that for a better overall $L_{\infty}$ bound it would be useful to pick boundary functionals at the boundary peaks of the left plot. This calls for a variation of the greedy method that monitors $P_{\Lambda_{N}}$ on all $\delta$ functionals as well, and picks a boundary functional as soon as the maximum is on the boundary. This will work for Dirichlet problems, but not in general circumstances in Hilbert spaces. See Sections 4.3 for theory and 6.4 for numerical results, respectively.

The basis functions $v_{\mu_{N}}$ are orthonormal in Hilbert space but not in $L_{2}$, decay with $N$, and show a sharp bell-shape for large $N$. Figure 6 shows a case like Figure 5, but for $N=150$. The new point, marked with an asterisk, is preferred over any boundary point by the greedy method, though it is close to the boundary. The greedy method, as is, does not need many $\delta$ functionals on the boundary, if it has plenty of $\delta\circ L$ functionals on or near the boundary.

The norm and the condition estimate of the transformation matrix that takes the $\lambda_{j}(u)$ values into the $\mu_{j}(u)$ values, based on (17), are given in Figure 7, behaving roughly like $N^{0.69}$ and $N^{1.7}$, respectively. This is like $h_{\Omega}^{-1/3}$ and $h_{\Omega}^{-2/3}$, respectively, in terms of the fill distance $h_{\Omega}$ in the domain.

To get an idea of the decay of the basis functions $v_{\mu_{j}}$ in $L_{2}$ context, the matrix of the values in a fine point set on the domain is calculated and a singular value decomposition is done on that matrix. The result is shown in Figure 8, the decay behavior of singular values being roughly like $N^{-2.4}$.

Note that each basis function is a worst case for the preceding steps, because it has zero data for them and is approximated by the zero function. Thus $L_{\infty}$ or $L_{2}$ norms of the basis functions are closely related to the worst-case $L_{\infty}$ or $L_{2}$ norms of solutions with Hilbert space norm one. Figure 9 shows the RMSQ and $L_{\infty}$ norms of the basis functions $v_{\mu_{j}}$ as functions of $j$, the estimated decay being like $N^{-1.65}$ and $N^{-1.44}$, respectively, but with serious roundoff pollution for large $j$. The peaks are exactly where boundary points are chosen by the greedy method, and this is explained via (19) by the identity

[TABLE]

that calls for a large $v_{\mu_{j}}$ when there is a sharp drop from $\rho_{\Lambda_{j-1}}$ to $\rho_{\Lambda_{j}}$.

6.2 Observations for varying smoothness

We now check the behavior of the greedy method when the smoothness parameter $m$ of the Hilbert spaces $W_{2}^{m}(\mathbb{R}^{d})$ changes. In view of (22) and (23), and since we saw before that the greedy method focuses on the differential operator and the domain, not on boundary values, Figure 10 shows rates $-0.67$ and $-0.71$ as functions of $m$ that confirm the $-(m-3)/2$ rate of (22) for both $\sigma_{\Lambda_{500}}$ and $\rho_{\Lambda_{500}}$ as functions of $m$ after 500 steps of the greedy method.

6.3 Use of Basis Functions

To get examples of solving specific Poisson problems via (20) using the basis obtained by the Greedy Method, we selected a run for smoothness $m=6$ to generate the basis first, the other parameters being as in the previous examples. Then the functions $f$ and $g$ of (1) were defined to let the true solution be of the form $K(\cdot,z)$ for a different kernel $K$ and a point $z=(-\pi/10,0)$. This allows to check cases with solutions of different smoothness.

The first case is the infinitely smooth situation where $K$ is a Gaussian. Figure 11 shows a very fast decay of the error, a fast increase of the cumulative sum of the coefficients $\mu_{j}^{2}(u)$ from (20), and the final maximal error $8\cdot 10^{-6}$, roughly. A less smooth solution is of the form $r^{2.5}$ with $r=\|\cdot-z\|_{2}$ with a derivative singularity at $z$, and the corresponding results are in Figure 12. The $\mu_{j}^{2}$ decay much more slowly, and the solver has to fight with the derivative singularity at $z$. The two cases were aligned by factors to start with an error of roughly one using the same precalculated basis. Recall that classical Harmonic Analysis shows that decay rates of coefficients of orthonormal expansions depend on smoothness and determine convergence rates.

6.4 The Extended Greedy Method

The Extended Greedy Method from Section 4.3 was run in the same situations as in Sections 6.1 and 6.2. The results are given from Figure 13 on.

The rates in the left parts of Figures 1 and 13 are similar, but not in the right parts. Since the extended method picks more boundary points than the original method, the convergence rate on the delta functionals is now much better. Similarly, Figures 2 and 14 differ considerably. The error in the interior decays much better, and the fill distances in the domain and on the boundary show a better alignment.

Figure 15 adds a plot of the Power Function on the domain to the three plots of Figure 4, and it should be compared to Figure 5. One can see that the improved selection of boundary functionals now avoids large values of the Power Function on the boundary. The other results are very similar, and plots are omitted, except for the condition. The extended method is less stable, if Figures 7 and 16 are compared.

From Lemma 2 one might hope that the extended Greedy Method performs better as a function of smoothness, but this cannot be supported by experiments. The rates for the situation corresponding to Figure 10 come out to be roughly -0.69 and -0.71, respectively, but the right-hand plot is much smoother. The actual plot is suppressed.

7 Summary and Open Problems

Roughly, the $P$-greedy method for solving Dirichlet problems for second-order elliptic operators in $W_{2}^{m}(\mathbb{R}^{d})$ seems to behave like the comparable $P$-greedy method for interpolation of functions in $W_{2}^{m-2}(\mathbb{R}^{d})$, with all its pros and cons. It focuses on the domain, not on the boundary, and it tends to produce an asymptotically uniform distribution of evaluation points there, with an unexpectedly small number of points for sampling the boundary values.

On the theoretical side, this opens the quest for a thorough analysis of Kolmogoroff $N$-widths for such PDE problems. A reasonable hypothesis is that these behave like those without differential operators, but for spaces of functions with lower-order smoothness.

Another observation is that the maximum of the generalized Power Function taken on all delta functionals, being a central quantity for pointwise error bounds, shows the same asymptotics as the maximum of the generalized Power Function taken on all chosen PDE data functionals. This is no surprise for well-posed problems, but it opens a way for explicitly computable factors for error bounds in terms of the Hilbert space norm of the true solution.

We supplied an Extended Greedy Method that cares for the delta functionals in a better way, but it falls out of the Hilbert space foundation, so far. Its analysis is another open problem.

Using the standard Hilbert space background [19, 12], we also can formulate a $P$-greedy method for solving Dirichlet problems for harmonic functions in 2D or 3D. It will follow the Kolmogoroff $N$-width theory for such cases, but the latter seems to be open.

There are other greedy techniques on the market (“$f$–greedy” and “$f/P$-greedy”) that apply to specific problems of the form (1), not uniformly to the whole class. But their convergence analysis is less far developed, see [18].

Acknowledgment

Special thanks go to Gabriele Santin for several helpful remarks.