Jul 20, 2017 by Bennet Carstensen •

GCA Part 1 - Green's second identity and quadrature

To receive a data-sparse approximation of a discretized matrix, which is both fast and reliable, applying cross approximation to the matrix directly would lead either to high runtime complexity of full pivoting or to potentially unreliable results of for example heurisitc pivoting strategies and error estimators.

For this reason, in a first step, we use Green’s second identity to gain a degenerate approximation of the kernel function and thus a $\mathcal{H}$-matrix representation of the Galerkin discretization [1].

Let $d \in \{2, 3\}$, $\omega \subseteq \mathbb{R}^d$ be a Lipschitz domain and $u \colon \bar{\omega} \to \mathbb{R}$ be harmonic in $\omega$. Green’s second identity (cf., e.g. [11, Theorem 2.2.2]) states

$$u(x) = \int\limits_{\partial \omega} g(x, z) \frac{\partial u_{z}}{\partial n}(z) dz - \int\limits_{\partial \omega} \frac{\partial g_{z}}{\partial n}(x, z) u(z) dz \qquad \text{for all} \qquad x \in \omega,$$

with $\frac{\partial u_{z}}{\partial n} and \frac{\partial g_{z}}{\partial n}, n \in \mathbb{R}^d,$ denoting the directional derivative in $z$ with regard to $n$, and $g$ still denoting a fundamental solution of the negative Laplace operator $- \Delta$.

Since the function $u(x) = g(x, y)$ is harmonic for any $y \notin \bar{\omega}$, replacing $u$ with $g$ results in

$$g(x, y) = \int\limits_{\partial \omega} g(x, z) \frac{\partial g_{z}}{\partial n}(z, y) dz - \int\limits_{\partial \omega} \frac{\partial g_{z}}{\partial n}(x, z) g(z, y) dz \qquad \text{for all} \qquad x \in \omega, y \notin \omega.$$

We pay particular attention to the right-hand side, where the variables $x$ and $y$ no longer appear together as arguments of $g$ and $\frac{\partial g}{\partial n}$. If $x$ and $y$ are sufficiently far from the boundary $\partial \omega$, the integrals are smooth [1]. This way, we can approximate the integrals by a quadrature rule. Denoting the quadrature weights by $(w_{i})_{i = 1}^{n}$ and the corresponding quadrature points by $(z_{i})_{i = 1}^{n}$, $n \in \mathbb{N}$, $w_{i} \in \mathbb{R}$ and $z_{i} \in \mathbb{R}^{d}$, $1 \leq i \leq n$, we can approximate Green’s second identity with

$$g(x, y) \approx \sum\limits_{i = 1}^n w_{i} g(x, z_{i}) \frac{\partial g_{z_{i}}}{\partial n}(z_{i}, y) - w_{i} \frac{\partial g_{z_{i}}}{\partial n}(x, z_{i}) g(z_{i}, y) \qquad \text{for all } x \in \omega, \ y \notin \bar{\omega}.$$

Since this is a degenerate approximation of the kernel function, approximating the Galerkin discretization leads to a hierarchical matrix.

In the interest of ensuring uniform convergence of the approximation, we have to choose a suitable admissiblity condition which ensures that $x$ and $y$ are sufficiently far from the boundary $\partial \omega$. One approach relies on bounding boxes: for a given cluster $t \in \mathcal{T}_{\mathcal{I}}$, we construct an axis-parallel box

$$\mathcal{B}_{t} = [a_{t}, b_{t}] = [a_{t, 1}, b_{t, 1}] \times \cdots \times [a_{t, d}, b_{t, d}], \qquad a_{t} = (a_{t, 1} \cdots a_{t, d})^{T}, b_{t} = (b_{t, 1} \cdots b_{t, d})^{T} \in \mathbb{R}^d,$$

containing the supports of all basis functions corresponding to the indices in $\hat{t}$ such that

$$supp \ \varphi_{i} \subseteq \mathcal{B}_{t} \qquad \text{for all } i \in \hat{t}.$$

Given a cluster $t \in \mathcal{T}_{\mathcal{I}}$, we define the diameter $\text{diam}_{\infty}(\mathcal{B}_{t})$ and distance $\text{dist}_{\infty}(\mathcal{B}_t, y)$ of a bounding box $\mathcal{B}_t$ to a vector $y \in \mathbb{R}^d$ with regard to the maximum norm

$$\begin{align*} \text{diam}_{\infty}(\mathcal{B}_{t} &:= \max{\{ \left\lVert x - y \right\rVert_{\infty} : x \in \mathcal{B}_{t}\}}\\ \text{dist}_{\infty}(\mathcal{B}_t, y) &:= \max{\{ \left\lVert x - y \right\rVert_{\infty} : x, y \in \mathcal{B}_{t} \}} \end{align*}.$$

This way, the farfield

$$\mathcal{F}_{t} := \{ y \in \mathbb{R}^d : \text{diam}_{\infty}(\mathcal{B}_{t}) \leq \text{dist}_{\infty}(\mathcal{B}_t, y)\},$$

describes the space in $\mathbb{R}^d$ where each $y \in \mathbb{F}_{t}$ is sufficiently far away from each $x \in \mathcal{B}_{t}$ for the approximation to converge, thus describing our admissiblity condition: a block $b = (t, s)$ is admissible if $\mathcal{B}_{s}$ is in the farfield of $t$, i.e., if $\mathcal{B}_{s} \subseteq \mathcal{F}_{t}$ holds true.

Now we apply Green’s second identity to the domain $\omega_{t}$ given by

$$\delta_{t} := \text{diam}_{\infty}(\mathcal{B}_{t}) / 2, \qquad \omega_{t} := [a_{t, 1} - \delta_{t}, b_{t, 1} + \delta_{t}] \times \cdots \times [a_{t, d} - \delta_{t}, b_{t, d} + \delta_{t}]$$

We represent $\omega_{t}$ with regard to the reference cube $[-1, 1]^{d}$, using the affine mapping

$$\Phi_{t} : [-1, 1]^d \to \omega_t, \hat{x} \mapsto \frac{b + a}{2} + \frac{1}{2} \begin{pmatrix} b_{t, 1} - a_{t, 1} + 2 \delta_{t} & & \\ & \ddots & \\ & & b_{t, d} - a_{t, d} + 2 \delta_{t} \\ \end{pmatrix} \hat{x},$$

which leads to the affine parametrizations

$$\gamma_{\iota}(\hat{z}) := \Phi_{t}(\hat{z}_{1}, \cdots, \hat{z}_{\kappa - 1}, \hat{z}_{\kappa}, \hat{z}_{\kappa + 1}, \cdots, \hat{z}_{d - 1}), \qquad \hat{z}_{\kappa} := \begin{cases} -1, & \kappa = 2 \iota - 1 \\ 1, & \kappa = 2 \iota \end{cases} \\ \text{for all } \iota \in \mathbb \{1, \cdots, 2d\}, \kappa \in \{1, \cdots, d\} \text{ and } \hat{z} \in \mathcal{Q},$$

of the boundary $\partial \omega$ with

$$\mathcal{Q} := [-1, 1]^{d - 1}$$

being the parameter domain for one side of the boundary such that

$$\partial \omega = \bigcup\limits_{\iota = 1}^{2d} \gamma_{\iota}(\mathcal{Q}), \qquad \int\limits_{\partial \omega} f(z) dz = \sum\limits_{\iota = 1}^{2d} \int\limits_{\mathcal{Q}} \sqrt{\det(D \gamma_{\iota}^{*} D \gamma_{\iota})} f(\gamma_{\iota}(\hat{z})) d \hat{z}.$$

We approximate the integral on the right side by a tensor quadrature formula: let $m \in \mathbb{N}$, $\xi_{1}, \cdots, \xi_{m} \in [-1, 1]$ the points and $w_{1}, \cdots, w_{m} \in \mathbb{R}$ the weights of the one-dimensional $m$-point Gauss quadrature formula for the reference interval [-1, 1]. Defining

$$\hat{z}_{\mu} := (\xi_{\mu 1}, \cdots, \xi_{\mu d - 1})^{T} \in \mathbb{R}^{d - 1}, \qquad \hat{w}_{\mu} := \prod\limits_{j = 1}^{d - 1} w_{\mu j} \in \mathbb{R} \qquad \text{for all } \mu \in M := \{1, \cdots, m\}^{d - 1},$$

introduces us to the tensor quadrature formula

$$\int\limits_{\mathcal{Q}} \hat{f}(\hat{z}) d \hat{z} \approx \sum\limits_{\mu \in M} \hat{w}_{\mu} \hat{f}(\hat{z}).$$

Applying this formula to all surfaces of $\partial \omega$ yields

$$\int\limits_{\partial \omega} f(z) dz \approx \sum\limits_{\iota = 1}^{2d} \sum\limits_{\mu \in M} \hat{w}_{\mu} \sqrt{\det(D \gamma_{\iota}^{*} D \gamma_{\iota})} f(\gamma_{\iota}(\hat{z}_{\mu})) = \sum\limits_{(\iota, \mu) \in K} w_{\iota \mu} f(z_{\iota \mu}),$$

with

$$K := {1, \cdots, 2d} \times M, \qquad w_{\iota \mu} := \hat{w}_{\mu} \sqrt{\det(D \gamma_{\iota}^{*} D \gamma_{\iota})}, \qquad z_{\iota \mu} := \gamma_{\iota}(\hat{z}_{\mu}).$$

Using this quadrature formula on the approximation of Green’s second identity yields

$$\tilde{g}_{t}(x, y) := \sum\limits_{(\iota, \mu) \in K} w_{\iota \mu} g(x, z_{\iota \mu}) \frac{\partial g_{z_{\iota \mu}}}{n_{\iota}}(z_{\iota \mu}, y) - w_{\iota \mu} \frac{\partial g_{z_{\iota \mu}}}{n_{\iota}}(x, z_{\iota \mu}) g(z_{\iota \mu}, y) \\ \text{for all } x \in \mathcal{B}_{t}, \ y \in \mathcal{F}_{t},$$

where $n_{\iota}$ denotes the outer normal vector of the face $\gamma_{\iota}(\mathcal{Q})$ of $\omega_{t}$.

Be an admissible block $b = (t, s)$ given, replacing the kernel function $g$ with $\tilde{g}$ results in a low-rank approximation

$$G|_{\hat{t} \times \hat{s}} \approx A_{t} B_{ts}^{*} \qquad A_t = (A_{t+} \ A_{t-}), \qquad, B_{ts} = (B_{ts+} \ B_{ts-}), \\ A_{t+}, A_{t-} \in \mathbb{R}^{\hat{t} \times K} \text{ and } B_{ts+}, B_{ts-} \in \mathbb{R}^{\hat{s} \times K},$$

with

$$\begin{align} a_{t+, i \nu} &:= \int\limits_{\Omega} \Phi_{i}(x) g(x, z_{\nu}) dx, b_{ts+, j \nu} &&:= w_{\nu} \int\limits_{\Omega} \Psi_{j}(x) \frac{\partial g_{z_{\nu}}}{n_{\iota}}(z_{\nu}, y) dy, \\ a_{t-, i \nu} &:= w_{\nu} \int\limits_{\Omega} \Phi_{i}(x) \frac{\partial g_{z_{\nu}}}{n_{\iota}}(x, ) dx b_{ts-, j \nu} &&:= - \int\limits_{\Omega} \Psi_{j}(x) g(z_{\nu}, y) dy \end{align} \\ \text{for all } \nu = (\iota, \mu) \in K, i \in \hat{t}, j \in \hat{s}.$$

Important to note is the fact that $A_{t}$ solely depends on $t$ and not on $s$. This will allow us to extend this construction to a $\mathcal{H}^2$-matrix later on.

Remark (Axis-parallel bouding boxes and directional derivatives) In our case, $\omega_{t}$ is still an axis-parallel box. As a result, we can describe the outer normal vectors $n_{\iota}$ with

$$n_{\iota} := (n_1 \cdots n_d)^{T}, \qquad n_i = \begin{cases} z_{\kappa} & \kappa = 2 \iota - 1 \text{ or } \kappa = 2 \iota \\ 0 & \text{otherwise} \end{cases}, \qquad \kappa \in \{1, \cdots, d\}.$$

This way, we can reduce the directional derivative $\frac{\partial g_{z_{\nu}}}{n_{\iota}}$ to a scaled partial derivative in the $\kappa$-th component of $z_{\nu}$: $\hat{z}_{\kappa} \frac{\partial g}{\partial z_{\nu, \kappa}}$. This will change the entries of the low-rank approximation as follwoing:

$$\begin{align} a_{t+, i \nu} &:= \int\limits_{\Omega} \Phi_{i}(x) g(x, z_{\nu}) dx, &&b_{ts+, j \nu} := w_{\nu} \int\limits_{\Omega} \Psi_{j}(x) \hat{z}_{\kappa} \frac{\partial g}{\partial z_{\nu, \kappa}}(z_{\nu}, y) dy, \\ a_{t-, i \nu} &:= w_{\nu} \int\limits_{\Omega} \Phi_{i}(x) \hat{z}_{\kappa}\frac{\partial g}{\partial z_{\nu, \kappa}}(x, z_{\nu}) dx &&b_{ts-, j \nu} := - \int\limits_{\Omega} \Psi_{j}(x) g(z_{\nu}, y) dy \end{align} \\ \text{for all } \nu = (\iota, \mu) \in K, i \in \hat{t}, j \in \hat{s}.$$

An applications with a more complex kernel function $g$ might consider this method to save evaluations of the derivatives. Yet, one restricts themself to a specific kind of bouding boxes. Since the fundamental solution of the negative Laplace operator is rather simple to evaluate and since we are using the implementation of the BEM-Modules of the H2Lib$^1$, we continue with the method using the directional derivative.

---

$^1$ The BEM-Modules are available for 2D- and 3D-Problems