Fast random field generation with $H$-matrices

TL;DR

This paper introduces a linear-cost method using $H$-matrices to efficiently generate normal and log-normal random fields on arbitrary point sets, with proven error bounds and broad applicability.

Contribution

It presents a novel $H$-matrix based approach for approximate square root computation of covariance matrices, enabling fast random field generation on general point sets.

Findings

Linear computational cost for random field generation

Rigorous error estimates demonstrating convergence

Applicable to non-stationary covariance functions and irregular point sets

Abstract

We use the $H$-matrix technology to compute the approximate square root of a covariance matrix in linear cost. This allows us to generate normal and log-normal random fields on general point sets with optimal cost. We derive rigorous error estimates which show convergence of the method. Our approach requires only mild assumptions on the covariance function and on the point set. Therefore, it might be also a nice alternative to the circulant embedding approach which applies only to regular grids and stationary covariance functions.

Tables

Table 1. Table 1 : Condition numbers of 𝑪 𝑪 \boldsymbol{C} for the covariance function from ( 2 ) with 𝒩 𝒩 {\mathcal{N}} being a Sobol point set with 2 m superscript 2 𝑚 2^{m} points.

$m=$	5	6	7	8	9	10
$\lambda =1$	2.0e+09	6.1e+16	8.6e+17	2.6e+19	1.8e+20	1.4e+20
$\lambda ={10}^{-1}$	3.9e+07	5.5e+14	1.8e+17	8.4e+18	4.8e+20	4.6e+20
$\lambda ={10}^{-2}$	6.5e+06	2.6e+12	2.7e+17	1.2e+19	3.3e+19	2.8e+20
$\lambda ={10}^{-3}$	4.2e+06	9.4e+11	6.1e+17	4.2e+18	2.6e+19	1.1e+20

Equations

$$\displaystyle\mathcal{Z}({\boldsymbol{x}}_{i},\omega)=\big{(}\boldsymbol{C}^{1/2}{\boldsymbol{z}}(\omega)\big{)}_{i}\quad\text{for all }i\in\{1,\ldots,N\},$$

$$\mathcal{Z}({\mathbf{x}}_i,\omega )\approx {\mathcal{Z}}_{k,p}(\mathbf{z}(\omega ){)}_i\text{for all }i\in \{1,\dots ,N\}.$$

$${\lambda }_{\max }(\mathbf{M}):=\underset{\mathbf{z}\in {\mathbb{R}}^N\setminus \{0\}}{\sup }\frac{|\mathbf{M}\mathbf{z}|}{|\mathbf{z}|}\text{and}{\lambda }_{\min }(\mathbf{M}):=\underset{\mathbf{z}\in {\mathbb{R}}^N\setminus \{0\}}{\inf }\frac{(\mathbf{M}\mathbf{z}{)}^T\mathbf{z}}{|\mathbf{z}{|}^2}.$$

$$\mathcal{Z}(\mathbf{x},\omega )\text{or}\exp (\mathcal{Z}(\mathbf{x},\omega ))\text{for all }\omega \in \Omega ,\mathbf{x}\in D$$

$$|{\partial }_{\mathbf{x}}^{\alpha }{\partial }_{\mathbf{y}}^{\beta }\varrho (\mathbf{x},\mathbf{y})|\le c_1(c_2|\mathbf{x}-\mathbf{y}|{)}^{-|\alpha {|}_1-|\beta {|}_1}|\alpha +\beta {|}_1!\text{for all }\mathbf{x}\ne \mathbf{y}\in D,$$

$$\displaystyle\Big{(}\mathcal{Z}({\boldsymbol{x}},\omega)\Big{)}_{{\boldsymbol{x}}\in{\mathcal{N}}}\in{\mathbb{R}}^{N}\quad\text{or}\quad\Big{(}\exp(\mathcal{Z}({\boldsymbol{x}},\omega))\Big{)}_{{\boldsymbol{x}}\in{\mathcal{N}}}\in{\mathbb{R}}^{N}$$

$$\displaystyle\varrho({\boldsymbol{x}},{\boldsymbol{y}})=\sigma^{2}\frac{2^{1-\mu}}{\Gamma(\mu)}\Big{(}\sqrt{2\mu}\frac{|{\boldsymbol{x}}-{\boldsymbol{y}}|_{p}}{\lambda}\Big{)}^{\mu}K_{\mu}\Big{(}\sqrt{2\mu}\frac{|{\boldsymbol{x}}-{\boldsymbol{y}}|_{p}}{\lambda}\Big{)},$$

$$\displaystyle\varrho({\boldsymbol{x}},{\boldsymbol{y}})=\sigma^{2}\exp\Big{(}-\frac{|{\boldsymbol{x}}-{\boldsymbol{y}}|_{p}}{\lambda}\Big{)}$$

$$\displaystyle\varrho({\boldsymbol{x}},{\boldsymbol{y}})=\sigma^{2}\exp\Big{(}-\frac{|{\boldsymbol{x}}-{\boldsymbol{y}}|_{p}^{2}}{2\lambda^{2}}\Big{)}.$$

$$\displaystyle\varrho({\boldsymbol{x}},{\boldsymbol{y}}):=\sigma^{2}\frac{{\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}})^{1/4}{\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{y}})^{1/4}}{\sqrt{2}{\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}}+\boldsymbol{\Sigma}_{\boldsymbol{y}})^{1/2}}\exp\Big{(}-({\boldsymbol{x}}-{\boldsymbol{y}})^{T}\frac{(\boldsymbol{\Sigma}_{\boldsymbol{x}}+\boldsymbol{\Sigma}_{\boldsymbol{y}})^{-1}}{2}({\boldsymbol{x}}-{\boldsymbol{y}})\Big{)}.$$

$$\underset{\alpha \in {\mathbb{N}}^d}{\sup }\underset{\mathbf{x}\in D}{\sup }\Vert {\partial }_{\mathbf{x}}^{\alpha }{\mathbf{\Sigma }}_{\mathbf{x}}\Vert <\infty .$$

$$\max \{\mathrm{diam}(B),\mathrm{diam}(B^{\prime })\}\le \eta \mathrm{dist}(B,B^{\prime }),$$

$$\mathrm{dist}(B,B^{\prime }):=\underset{\mathbf{x}\in B,\mathbf{y}\in B^{\prime }}{\inf }|\mathbf{x}-\mathbf{y}|.$$

$$\displaystyle\left\{\;\qquad\parbox{346.89731pt}{\begin{itemize} \par\itemize@item@[\textbf{if}] $B_{X}$ and $B_{Y}$ satisfy\leavevmode\nobreak\ \eqref{eq:adm} or if ${\rm sons}(X)=\emptyset={\rm sons}(Y)$ {set} ${\rm sons}(X,Y)=\emptyset$ \par\itemize@item@[\textbf{else if}]${\rm sons}(Y)\neq\emptyset$ and ${\rm sons}(X)=\emptyset$ {set} ${\rm sons}(X,Y)=\{X\}\times{\rm sons}(Y)$ \par\itemize@item@[\textbf{else if}] ${\rm sons}(X)\neq\emptyset$ and ${\rm sons}(Y)=\emptyset$ {set} ${\rm sons}(X,Y)={\rm sons}(X)\times\{Y\}$ \par\itemize@item@[\textbf{else}]${\rm sons}(X)\neq\emptyset$ and ${\rm sons}(Y)\neq\emptyset$ {set} ${\rm sons}(X,Y)={\rm sons}(X)\times{\rm sons}(Y)$ \par\end{itemize}}\right.$$

$$\displaystyle{\mathbb{T}}_{\rm far}:=\big{\{}(X,Y)\in{\mathbb{T}}\,:\,{\rm sons}(X,Y)=\emptyset\text{ and }B_{X},B_{Y}\text{ satisfy }\leavevmode\nobreak\ \eqref{eq:adm}\big{\}}$$

$$\displaystyle{\mathbb{T}}_{\rm near}:=\big{\{}(X,Y)\in{\mathbb{T}}\,:\,{\rm sons}(X,Y)=\emptyset\text{ and }B_{X},B_{Y}\text{ do not satisfy }\leavevmode\nobreak\ \eqref{eq:adm}\big{\}}.$$

$$X\times Y=\bigcup _{(X^{\prime },Y^{\prime })\in \mathrm{sons}(X,Y)}{X}^{\prime }\times Y^{\prime }$$

$$\varrho (\mathbf{x},\mathbf{y})\approx c_p^{XY}(\mathbf{x},\mathbf{y}):=\sum _{n,m=1}^{p^d}\varrho (q_n^X,q_m^Y)L_n^X(\mathbf{x})L_m^Y(\mathbf{y})\text{for all }\mathbf{x}\in X,\mathbf{y}\in Y.$$

$${\mathbf{V}}_{in}^X:=L_n^X({\mathbf{x}}_i),{\mathbf{V}}_{jm}^Y:=L_m^Y({\mathbf{x}}_j),\text{ and}{\mathbf{M}}_{nm}^{XY}:=\varrho (q_n^X,q_m^Y)$$

$$\mathbf{C}{|}_{I(X)\times I(Y)}\approx {\mathbf{V}}^X{\mathbf{M}}^{XY}({\mathbf{V}}^Y{)}^T.$$

$${\mathbf{V}}_{in}^X:=L_n^X({\mathbf{x}}_i)=\sum _{m=1}^{p^d}{L}_n^X(q_m^{X^{\prime }})L_m^{X^{\prime }}({\mathbf{x}}_i)=\sum _{m=1}^{p^d}{L}_n^X(q_m^{X^{\prime }}){\mathbf{V}}_{im}^{X^{\prime }}=({\mathbf{V}}^{X^{\prime }}{\mathbf{T}}^{X^{\prime }X}{)}_{in}.$$

$$C_{\mathrm{u}}^{-1}{N}^{-1/d}\le \underset{\mathbf{x},{\mathbf{x}}^{\prime }\in \mathcal{N}}{\min }|\mathbf{x}-{\mathbf{x}}^{\prime }|\le \underset{\mathbf{x}\in D}{\sup }\underset{{\mathbf{x}}^{\prime }\in \mathcal{N}}{\min }|\mathbf{x}-{\mathbf{x}}^{\prime }|\le C_{\mathrm{u}}{N}^{-1/d}.$$

$$\displaystyle\|\boldsymbol{C}-\boldsymbol{C}_{p}\|_{2}\leq\|\boldsymbol{C}-\boldsymbol{C}_{p}\|_{F}:=\Big{(}\sum_{i,j=1}^{N}|\boldsymbol{C}-\boldsymbol{C}_{p}|_{ij}^{2}\Big{)}^{1/2}\leq C_{H}N(\log(p)+1)^{2d-1}\Big{(}\frac{\eta}{4c_{2}}\Big{)}^{p}.$$

$$\displaystyle C_{H}N(\log(p)+1)^{2d-1}\Big{(}\frac{\eta}{4c_{2}}\Big{)}^{p}<\lambda_{\rm min}(\boldsymbol{C}).$$

$${\mathcal{K}}_k:=\mathrm{span}\{\mathbf{z},\mathbf{M}\mathbf{z},{\mathbf{M}}^2\mathbf{z},\dots ,{\mathbf{M}}^{k-1}\mathbf{z}).$$

$$\mathbf{U}:={\mathbf{Q}}^T\mathbf{M}\mathbf{Q}.$$

$${\mathbf{M}}^{1/2}\mathbf{z}=\mathbf{Q}{\mathbf{U}}^{1/2}{\mathbf{Q}}^T\mathbf{z}.$$

$$\displaystyle\frac{|\boldsymbol{C}^{1/2}{\boldsymbol{z}}-\mathcal{Z}_{k,p}({\boldsymbol{z}})|}{|{\boldsymbol{z}}|}\leq\delta_{k_{0},k}{\sqrt{2\|\boldsymbol{M}\|_{2}}}\frac{4r^{{2}}}{r-1}r^{-k}+\frac{2C_{\rm H}N(\log(p)+1)^{2d-1}\Big{(}\frac{\eta}{4c_{2}}\Big{)}^{p}}{\max\{\lambda_{\rm min}(\boldsymbol{C}),\lambda_{\rm min}(\boldsymbol{C}_{p})\}^{1/2}},$$

$$r:=\frac{{\lambda }_{\max }({\mathbf{C}}_p)+{\lambda }_{\min }({\mathbf{C}}_p)}{{\lambda }_{\max }({\mathbf{C}}_p)-{\lambda }_{\min }({\mathbf{C}}_p)}>1.$$

$$|{\lambda }_i-{\lambda }_j|\le {\lambda }_{\max }({\mathbf{C}}_p)C_{\kappa }{\kappa }^{\min \{i,j\}}\text{for all }1\le i,j\le M$$

Full text

1 Introduction

Generating samples of random fields is a common bottleneck in simulation and modeling of real life phenomena as, e.g., structural vibrations ap2 , groundwater flow flow , and composite material behavior ap1 . A standard approach is to truncate the Karhunen-Loève expansion of the random field. This can, particularly for rough fields with short correlation length, be very expensive, as many summands of the expansion have to be evaluated to compute a decent approximation. Often, it suffices to evaluate the random field only on some particular (quadrature) nodes. If the random field $\mathcal{Z}({\boldsymbol{x}},\omega)$ is Gaussian with given covariance function $\varrho(\cdot,\cdot)$, it is well-known that the evaluation at the quadrature nodes ${\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{N}$ can be done by computing the square-root of the corresponding covariance matrix $\boldsymbol{C}=(\varrho({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j}))_{i,j\in\{1,\ldots,N\}}\in{\mathbb{R}}^{N\times N}$, i.e.,

[TABLE]

where ${\boldsymbol{z}}(\omega)\in{\mathbb{R}}^{N}$ is a vector of i.i.d. standard normal random numbers. Since each evaluation requires a matrix-vector multiplication with $\boldsymbol{C}^{1/2}$, a direct approach requires $\mathcal{O}(N^{2})$ operations for the multiplication plus $\mathcal{O}(N^{3})$ operations for computing the square-root itself and thus is prohibitively expensive. An efficient method first proposed in circ1 ; circ2 is circulant embedding, which employs fast FFT techniques to realize the factorization and the matrix-vector multiplication in $\mathcal{O}(N\log(N))$ operations. This approach, however, works solely for stationary covariance functions $\varrho({\boldsymbol{x}},{\boldsymbol{y}})=\rho(|{\boldsymbol{x}}-{\boldsymbol{y}}|)$ and regular grids of quadrature nodes. Since non-stationary covariance functions are of great interest for the modeling of natural structures (e.g., porous rock, wood,…), and since finite element methods often use irregular grids, we propose a new method which removes both restrictions.

The idea is to approximate the covariance matrix $\boldsymbol{C}$ by an $H^{2}$-matrix, as described in, e.g, h2mat , and to use an iterative method to compute an approximation $\mathcal{Z}_{k,p}({\boldsymbol{z}})$ ($k$ and $p$ are parameters of the methods, see below) to $\boldsymbol{C}^{1/2}{\boldsymbol{z}}$ for any ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$. We therefore obtain the approximation to the random field by feeding the algorithms with i.i.d. standard normal random vectors ${\boldsymbol{z}}(\omega)\in{\mathbb{R}}^{N}$, i.e.,

[TABLE]

This is feasible since matrix-vector multiplication with $H^{2}$-matrices can be done in $\mathcal{O}(N)$ operations. The only assumption on the covariance function of the random field is that it is asymptotically smooth. We propose two iterative algorithms, each with individual advantages for smooth or rough random fields. This algorithms might also be of interest for the approximation of random fields with covariance kernels of random solutions of certain stochastic operator equations, as considered in schwab .

The idea to use $H$-matrices for random field approximation has already been used indirectly in h2KL ; harbrecht , where the authors efficiently compute eigenfunctions of the covariance operator by use of $H$-matrix techniques.

1.1 Notation

Throughout the text, $\alpha\lesssim\beta$ denotes $\alpha\leq C\beta$ for some generic constant $C>0$ and $\alpha\simeq\beta$ means $\alpha\lesssim\beta$ and $\beta\lesssim\alpha$. The notation $|\cdot|$ has several unambiguous meanings: for vectors, it denotes the euclidean norm, while for sets, $|\cdot|$ is the natural measure, which is the Lebesgue measure (volume, area) for continuous sets and the counting measure (cardinality) for finite sets. The notation $\|\cdot\|_{2}$ is used for the spectral matrix norm and $|{\boldsymbol{z}}|_{p}:=(\sum_{j=1}^{N}|{\boldsymbol{z}}_{i}|^{p})^{1/p}$ for all ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$ denotes the $\ell_{p}$-norm. By ${\mathcal{P}}^{k}$ we denote the set of polynomials of maximal degree $k$. For brevity, we write $|\cdot|:=|\cdot|_{2}$. We denote the maximal and minimal eigenvalues of a positive definite and symmetric matrix $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ by

[TABLE]

We denote the $k$-th component of a vector ${\boldsymbol{v}}\in{\mathbb{R}}^{N}$ by ${\boldsymbol{v}}_{k}$, whereas sequences of vectors are denoted by ${\boldsymbol{v}}^{1},{\boldsymbol{v}}^{2},\ldots$.

2 Model Problem

Let $(\Omega,\Sigma,{\mathbb{P}})$ be a probability space and let $D\subseteq{\mathbb{R}}^{d}$, $d\in{\mathbb{N}}$ be a Lipschitz domain. We consider a random field which is normal or log-normal,

[TABLE]

for some zero-mean Gaussian random field $\mathcal{Z}(\cdot,\cdot)$ (note that the assumption on the mean is purely for brevity of presentation). The covariance function $\varrho\colon D\times D\to{\mathbb{R}}$ of $\mathcal{Z}(\cdot,\cdot)$ is assumed asymptotically smooth: that is, \varrho\in C^{\infty}\big{(}\big{\{}({\boldsymbol{x}},{\boldsymbol{y}})\in D\times D\,:\,{\boldsymbol{x}}\neq{\boldsymbol{y}}\big{\}}\big{)} and there exist constants $c_{1},c_{2}>0$ such that

[TABLE]

for all multi-indices $\alpha,\beta\in{\mathbb{N}}_{0}^{d}$ with $|\alpha|_{1}+|\beta|_{1}\geq 1$. (The expert reader will notice that the original definition of asymptotically smooth includes a singularity order. As our covariance functions are always finite in value, we do not consider this.) The goal of this work is to derive an efficient method which evaluates the random field at certain (quadrature) points ${\mathcal{N}}\subseteq D$, where ${\mathcal{N}}=\{{\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{N}\}$ is a finite set, i.e., we aim to approximate

[TABLE]

for given $\omega\in\Omega$.

2.1 Examples of valid covariance functions

The condition above includes the important class of isotropic stationary covariance functions of Matérn form, e.g.,

[TABLE]

where $\Gamma(\cdot)$ is the gamma function, $K_{\mu}$ is the modified Bessel function of second kind, and $\lambda,\sigma>0$, $\mu\in(0,\infty]$, $p\in{\mathbb{N}}$ are parameters. For $\mu=1/2$, the above function takes the form

[TABLE]

and the limit case $\mu=\infty$ satisfies

[TABLE]

Also much more general non-stationary, non-isotropic covariance functions, e.g.,

[TABLE]

satisfy the assumptions. Here, $\boldsymbol{\Sigma}_{(\cdot)}\colon D\to{\mathbb{R}}^{d\times d}$ is a smooth mapping into the symmetric positive definite matrices and $\sigma>0$ is a parameter. This covariance function was first suggested in cov to model spatially dependent anisotropies in a material.

Lemma 1

The covariance functions from (2) satisfy (1). Assume the mapping ${\boldsymbol{x}}\mapsto\boldsymbol{\Sigma}_{\boldsymbol{x}}$ satisfies (for any matrix norm $\|\cdot\|$)

[TABLE]

Then, the covariance function from (3) is asymptotically smooth (1).

We postpone the proof of the lemma to Appendix A.

3 Sampling the random field

By definition, $\mathcal{Z}({\boldsymbol{x}},\cdot)$, ${\boldsymbol{x}}\in{\mathcal{N}}$ is a Gaussian random field with covariance matrix $\boldsymbol{C}\in{\mathbb{R}}^{N\times N}$, $N=|{\mathcal{N}}|$, and $\boldsymbol{C}_{ij}=\varrho({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j})$, where we write ${\mathcal{N}}:=\{{\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{N}\}$. The main goal of this section is to establish a new way to efficiently approximate $\boldsymbol{C}^{1/2}{\boldsymbol{z}}$ for given ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$. Roughly, the strategy is to approximate $\boldsymbol{C}$ by an $H^{2}$-matrix and to benefit from the fast matrix-vector multiplication provided by it. This allows us to efficiently approximate $\boldsymbol{A}{\boldsymbol{z}}$ (without actually factorizing the matrix $\boldsymbol{C}$).

3.1 $H^{2}$-matrix approximation of the covariance matrix

Given the finite set of evaluation points ${\mathcal{N}}:=\{{\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{N}\}\subset D$, we approximate the covariance matrix $\boldsymbol{C}\in{\mathbb{R}}^{N\times N}$, $\boldsymbol{C}_{ij}:=\varrho({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j})$ by an $H^{2}$-matrix $\boldsymbol{C}_{p}$ via interpolation of order $p\in{\mathbb{N}}$.

In the following, we recall the definition of $H^{2}$-matrices and the approximation process as laid out in, e.g., h2mat . The rough idea is to partition the index set of the covariance matrix into far-field blocks, which can be approximated efficiently by interpolation of the covariance function, and near-field blocks, which are stored exactly.

3.1.1 Block partitioning

For each subset $X\subseteq{\mathcal{N}}$, we denote by $B_{X}\subseteq{\mathbb{R}}^{d}$, the smallest axis-parallel box such that $X\subseteq B_{X}$. We build a binary tree of clusters in the following way. Let $X_{\rm root}:={\mathcal{N}}=\{{\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{N}\}$ denote the root of the tree which has level zero ${\rm level}(X_{\rm root})=0$ by definition. For each node of the tree $X$ with $|X|>C_{\rm leaf}$ for some cut-off constant $C_{\rm leaf}{\geq 2}$ (usually $C_{\rm leaf}\approx 20$), we define two sons of $X$ as follows: Split $B_{X}$ in half along its longest edge into $B_{0}\cup B_{1}=B_{X}$. Define ${\rm sons}(X):=\{X_{0},X_{1}\}$ with $X_{0}:=X\cap B_{0}$ and $X_{1}:=X\setminus X_{0}$ and set ${\rm level}(X_{i})={\rm level}(X)+1$ for $i=0,1$. For a node $X$ with $|X|\leq C_{\rm leaf}$, we define ${\rm sons}(X):=\emptyset$. This procedure generates a binary tree denoted by ${\mathbb{T}}_{\rm cl}$ (where ${\rm cl}$ stands for cluster) and guarantees that its leaves satisfy $|X|\leq C_{\rm leaf}$.

For a parameter $\eta>0$, we consider the admissibility condition for axis parallel boxes $B,B^{\prime}\subseteq{\mathbb{R}}^{d}$

[TABLE]

where the euclidean distance between the bounding boxes is defined by

[TABLE]

The condition (5) will be used to build the block-cluster tree ${\mathbb{T}}\subseteq{\mathbb{T}}_{\rm cl}\times{\mathbb{T}}_{\rm cl}$ as follows. The root of ${\mathbb{T}}$ is $(X_{\rm root},X_{\rm root})$. For each node $(X,Y)\in{\mathbb{T}}$ of the tree, define ${\rm sons}(X,Y)$, the set of sons, as:

[TABLE]

We also define the level as ${\rm level}(X_{\rm root},X_{\rm root})=0$ and ${\rm level}(X,Y)={\rm level}(X^{\prime},Y^{\prime})+1$ for $(X,Y)\in{\rm sons}(X^{\prime},Y^{\prime})$. Further, we define

[TABLE]

as well as

[TABLE]

Note that by definition of the block-cluster tree ${\mathbb{T}}$, the set ${\mathbb{T}}_{\rm near}\cup{\mathbb{T}}_{\rm far}$ contains all the leaves of ${\mathbb{T}}$. Moreover, we see that for each $(X,Y)\in{\mathbb{T}}\setminus({\mathbb{T}}_{\rm near}\cup{\mathbb{T}}_{\rm far})$, there holds

[TABLE]

Therefore, ${\mathbb{T}}_{\rm near}\cup{\mathbb{T}}_{\rm far}$ is a partition of ${\mathcal{N}}\times{\mathcal{N}}$ in the sense that each pair of points $({\boldsymbol{x}}_{i},{\boldsymbol{x}}_{j})\in{\mathcal{N}}\times{\mathcal{N}}$ for $1\leq i,j\leq N$ is contained in exactly one $(X,Y)\in{\mathbb{T}}_{\rm near}\cup{\mathbb{T}}_{\rm far}$.

3.1.2 Interpolation

The blocks $(X,Y)\in{\mathbb{T}}_{\rm far}$ satisfy (5) and hence interpolation of the kernel function is highly accurate. This allows us to store the matrix very efficiently. Let I(X):=\big{\{}i\in{\mathbb{N}}\,:\,{\boldsymbol{x}}_{i}\in X\big{\}} denote the index set of $X$. The basic idea now is to replace $\boldsymbol{C}|_{I(X)\times I(Y)}$ by a low-rank approximation $\boldsymbol{V}^{X}\boldsymbol{M}^{XY}(\boldsymbol{V}^{Y})^{T}$ with $\boldsymbol{V}^{X}\in{\mathbb{R}}^{|X|\times p^{d}}$, $\boldsymbol{M}^{XY}\in{\mathbb{R}}^{p^{d}\times p^{d}}$, and $\boldsymbol{V}^{Y}\in{\mathbb{R}}^{|Y|\times p^{d}}$, where $p$ is the interpolation order. The three matrices are defined by Chebychev interpolation of the covariance function. To that end, let $\{q_{1}^{X},\ldots,q_{p^{d}}^{X}\}$ denote transformed, tensorial Chebychev nodes in $B_{X}$ with the corresponding Lagrange basis functions $L_{1}^{X},\ldots,L_{p^{d}}^{X}\colon B_{X}\to{\mathbb{R}}$. Given $(X,Y)\in{\mathbb{T}}_{\rm far}$, we may approximate

[TABLE]

For $i,j\in\{1,\ldots,N\}$ and $n,m\in\{1,\ldots,p^{d}\}$, this leads to

[TABLE]

and hence

[TABLE]

The admissibility condition (5) guarantees that the approximation error converges to zero exponentially in $p$, as we prove in Proposition 1 below. Further note that the Chebychev interpolation described above is exact on polynomials of degree $p$. Thus, for $X\in{\mathbb{T}}_{\rm cl}$ and ${\boldsymbol{x}}_{i}\in X^{\prime}\in{\rm sons}(X)$, there holds with the transfer matrices $\boldsymbol{T}^{X^{\prime}X}:=(L^{X}_{n}(q_{m}^{X^{\prime}}))_{mn}\in{\mathbb{R}}^{p^{d}\times p^{d}}$

[TABLE]

Thus, it suffices to store $\boldsymbol{V}^{X}$ only for the leaves of ${\mathbb{T}}_{\rm cl}$ together with the transfer matrices $\boldsymbol{T}^{X^{\prime}X}$. This enables very efficient storage and arithmetics for $H^{2}$ matrices.

The capabilities of $H^{2}$-matrices which we employ in this work are summarized below in Proposition 1. To that end, we assume that the points ${\mathcal{N}}$ are approximately uniformly distributed, in the following sense.

Assumption 1** (quasi-uniform distribution)**

We say that ${\mathcal{N}}$ is quasi-uniformly distributed if there exists a constant $C_{\rm u}>0$ such that

[TABLE]

Proposition 1

Suppose we have a covariance matrix $\boldsymbol{C}\in{\mathbb{R}}^{N\times N}$ and an asymptotically smooth kernel $\varrho(\cdot,\cdot)$ and recall Assumption 1 on approximate uniform distribution of ${\mathcal{N}}$. Then, there exists a constant $C_{H}>0$ such that, for all $p\in{\mathbb{N}}_{0}$, the $H^{2}$-matrix $\boldsymbol{C}_{p}\in{\mathbb{R}}^{N\times N}$ constructed as above satisfies

[TABLE]

(The constant $c_{2}$ is defined in (1).) The $H^{2}$-matrix $\boldsymbol{C}_{p}$ is symmetric and can be stored using less than $C_{H}p^{2d}N$ memory units. Moreover, given any vector ${\boldsymbol{x}}\in{\mathbb{R}}^{N}$, it is possible to compute $\boldsymbol{C}_{p}x\in{\mathbb{R}}^{N}$ in less than $C_{H}p^{2d}N$ arithmetic operations. The constant $C_{H}$ depends only on $C_{\rm leaf}$ and $d$. The matrix $\boldsymbol{C}_{p}$ is positive definite if $p$ is sufficiently large such that

[TABLE]

We postpone the proof of the lemma to Appendix B.

3.2 Computing the square-root (Method 1)

Since $\boldsymbol{C}$ is positive definite in our case, a standard method is to compute the Cholesky factorization $\boldsymbol{L}\boldsymbol{L}^{T}=\boldsymbol{C}$. This can be done using $H^{2}$-matrices in almost linear cost (analyzed in hmatrices for $H$-matrices, but the method transfers to $H^{2}$-matrices). However, to the authors’ best knowledge, there is no complete error analysis available, and due to the complicated structure of the algorithm, the worst-case error estimate may be overly pessimistic. Therefore, we propose an iterative algorithm based on a variant of the Lanczos iteration. Note that polynomial or rational approximations of the square root (as pursued in, e.g., rational ) are doomed to fail since smooth random fields result in very badly conditioned covariance matrices $\boldsymbol{C}$ (see also the numerical experiments below). This implies that a polynomial approximation of the square root over the spectrum of $\boldsymbol{C}$ is very costly, whereas a rational approximation requires the inverse of $\boldsymbol{C}$ which is hard to compute due to the bad condition number.

The idea behind the algorithm below is as follows. Given a positive definite symmetric matrix $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ and a vector ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$, the aim is to compute efficiently an approximation to $\boldsymbol{M}^{1/2}{\boldsymbol{z}}$. For arbitrary $k\leq N$ define the order-$k$ Krylov subspace of $\boldsymbol{M}$ and ${\boldsymbol{z}}$ as

[TABLE]

Assuming ${\mathcal{K}}_{k}$ is $k$-dimensional, consider the orthogonal matrix $\boldsymbol{Q}\in{\mathbb{R}}^{N\times k}$ whose columns are the orthonormal basis vectors of the Krylov subspace, i.e., $\boldsymbol{Q}^{T}\boldsymbol{Q}=\boldsymbol{I}_{k}$ and ${\rm range}(\boldsymbol{Q})={\mathcal{K}}_{k}$. Now define $\boldsymbol{U}\in{\mathbb{R}}^{k\times k}$ by

[TABLE]

If $k=N$ then $\boldsymbol{Q}\boldsymbol{Q}^{T}=\boldsymbol{I}_{N}$ and $\boldsymbol{Q}\boldsymbol{U}\boldsymbol{Q}^{T}=\boldsymbol{M}$, from which it follows that

[TABLE]

The algorithm relies on explicit matrix multiplication to construct $\boldsymbol{U}$ and then a direct factorization of $\boldsymbol{U}$, thus for large $N$ it is feasible only when $k\ll N$, in which case (9) does not hold exactly. However, as we show later it may hold to a good enough approximation. The following Lanczos type algorithm builds up progressively the columns of $\boldsymbol{Q}$ without fully computing ${\mathcal{K}}_{k}$ first.

Remark 1

In the following, we make frequent use of the $QR$-factorisation of matrices and therefore recall the most important facts: For a matrix $\boldsymbol{A}\in{\mathbb{R}}^{n\times k}$ with $k\leq n\in{\mathbb{N}}$, there exists a $QR$-factorization $\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{R}$ such that $\boldsymbol{Q}\in{\mathbb{R}}^{n\times k}$ and $\boldsymbol{R}\in{\mathbb{R}}^{k\times k}$. The columns of $\boldsymbol{Q}$ are orthonormal and for $1\leq j\leq{\rm rank}(\boldsymbol{A})$, the first $j$ columns of $\boldsymbol{Q}$ span the same linear space as the first $j$ columns of $\boldsymbol{A}$. Moreover, $\boldsymbol{R}$ is upper triangular. If we restrict to positive diagonal entries of $\boldsymbol{R}$, the factorization is unique if ${\rm rank}(\boldsymbol{A})=k$.

Algorithm 1

Input: positive definite symmetric matrix $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$, vector ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$, and maximal number of iterations $k\in{\mathbb{N}}$.

Compute Krylov subspace: Set $\boldsymbol{Q}_{1}:={\boldsymbol{z}}/|{\boldsymbol{z}}|\in{\mathbb{R}}^{N\times 1}$ and $k_{0}=k$. For $j=2,\ldots,k$ do:

1. (a)

Compute $\widetilde{{\boldsymbol{q}}}:=\boldsymbol{M}{\boldsymbol{q}}^{j-1}\in{\mathbb{R}}^{N}$, where ${\boldsymbol{q}}^{j-1}$ is the $(j-1)$-th column of $\boldsymbol{Q}_{j-1}\in{\mathbb{R}}^{N\times(j-1)}$. 2. (b)

Compute $QR$-factorization $\boldsymbol{Q}_{j}\in{\mathbb{R}}^{N\times j}$ (with orthonormal columns), $\boldsymbol{R}_{j}\in{\mathbb{R}}^{j\times j}$ (upper triangular) such that $\boldsymbol{Q}_{j}\boldsymbol{R}_{j}=(\boldsymbol{Q}_{j-1},\widetilde{{\boldsymbol{q}}})\in{\mathbb{R}}^{N\times j}$. 3. (c)

If $(\boldsymbol{R}_{j})_{jj}=0$, set $k_{0}=j-1$ and goto Step 2. 2. 2.

Compute $\boldsymbol{U}_{k_{0}}:=\boldsymbol{Q}_{k_{0}}^{T}\boldsymbol{M}\boldsymbol{Q}_{k_{0}}\in{\mathbb{R}}^{{k_{0}}\times{k_{0}}}$. 3. 3.

Compute $\boldsymbol{U}_{k_{0}}^{1/2}$ directly. 4. 4.

Return ${\boldsymbol{y}}=\boldsymbol{Q}_{k_{0}}\boldsymbol{U}_{k_{0}}^{1/2}\boldsymbol{Q}_{k_{0}}^{T}{\boldsymbol{z}}$.

Output: Approximation ${\boldsymbol{y}}\approx\boldsymbol{M}^{1/2}{\boldsymbol{z}}$ and number of steps $k_{0}$.

Remark 2

Obviously, the orthogonal basis ${\boldsymbol{q}}^{1},\ldots,{\boldsymbol{q}}^{k}$ could also be generated by Gram-Schmidt orthogonalization. However, numerical experiments show that this is not stable with respect to roundoff errors. Moreover, also the classical Lanczos algorithm seems to be prone to rounding errors, especially for ill-conditioned matrices. Therefore, we propose to use the $QR$-factorization as above.

Remark 3

As proved in Lemma 4 below (and as is easily verified), a generic $QR$-algorithm produces $\boldsymbol{Q}_{j}$ which coincides with the first $j$ columns of $\boldsymbol{Q}$ up to signs. For simplicity, we assume in the following that the $QR$-algorithm ensures that the diagonal entries of $\boldsymbol{R}_{j}$ are always non-negative. This guarantees that the first $j$ columns of $\boldsymbol{Q}_{j+1}$ coincide with $\boldsymbol{Q}_{j}$. Thus, it suffices to store only the new column ${\boldsymbol{q}}^{j}$.

Theorem 3.1

Let $0<\eta<4c_{2}$ and let $p$ be sufficiently large such that $\boldsymbol{C}_{p}$ constructed from $\boldsymbol{C}$ as in Section 3.1 is positive definite (condition (7) is sufficient), and suppose Assumption 1 holds. Given ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$, call Algorithm 1 with $\boldsymbol{M}=\boldsymbol{C}_{p}$, ${\boldsymbol{z}}$, and a maximal number of iterations $k\in{\mathbb{N}}$. The output of Algorithm 1 contains the approximation $\mathcal{Z}_{k,p}({\boldsymbol{z}}):={\boldsymbol{y}}\in{\mathbb{R}}^{N}$ to $\boldsymbol{C}^{1/2}{\boldsymbol{z}}$ and the step number $k_{0}\leq k$.

• (i)

There holds with Kronecker’s delta $\delta_{i,j}$

[TABLE]

where $C_{\rm H}$, $\eta$, $c_{2}$, and $p$ are as in Proposition 1, and

[TABLE]

• (ii)

Let $\lambda_{\rm max}(\boldsymbol{C}_{p})=\lambda_{1}>\lambda_{2}>\ldots>\lambda_{M}>0$ denote the distinct eigenvalues of $\boldsymbol{C}_{p}$ for some $M\leq N$ and assume

[TABLE]

for some $C_{\kappa}>0$ and $0<\kappa<1$, then

[TABLE]

The algorithm completes in $\mathcal{O}(k^{3}p^{2d}N)$ arithmetic operations and uses less than $\mathcal{O}(kN)$ storage.

Remark 4

The theorem covers two regimes of covariance matrices. Whereas case (i) is the classical Lanczos convergence analysis for well-conditioned matrices, case (ii) considers ill-conditioned matrices with rapidly decaying eigenvalues. The numerical examples in Section 4 suggest that the error estimates might be more or less sharp, since Algorithm 1 performs remarkably well for smooth random fields (with rapidly decaying eigenvalues) and very rough random fields (with well-conditioned covariance matrices). Note that $k_{0}<k$ (hence $\delta_{k_{0},k}=0$) implies that the condition in the if-clause 1(c) is true. This however is an exotic case, meaning that ${\boldsymbol{z}}$ lies some non-trivial invariant subspace of $\boldsymbol{C}_{p}$ with fewer than $k$ dimensions. In this situation the algorithm computes $\boldsymbol{C}_{p}^{1/2}{\boldsymbol{z}}$ exactly and only the $H$-matrix approximation error remains. We note that by use of (16) instead of (17) in the proof below, it is possible to replace $(k+1)/4$ by $(k+1)/2$ and $p/2$ by $p$ in the exponents in (ii) at the price of including the square-root of the minimal eigenvalue in the denominator as in (i).

Proof (Proof of Theorem 3.1)

The cost estimate is proved as follows. The Krylov subspace loop of Algorithm 1 completes at most $k$ iterations. In each iteration, we have one $H^{2}$-matrix-vector multiplication which needs $\mathcal{O}(p^{2d}N)$ operations. Moreover, the $QR$-factorization needs $\mathcal{O}(Nk^{2})$ arithmetic operations. After the matrix $\boldsymbol{Q}_{k}$ is set up, we have $k$ $H^{2}$-matrix-vector multiplications to compute $\boldsymbol{M}\boldsymbol{Q}_{k}$ and $k^{2}$ scalar products to compute $\boldsymbol{U}_{k_{0}}$. In total, this needs $\mathcal{O}(N(k+k^{2}))$ arithmetic operations. The computation of $\boldsymbol{U}_{k_{0}}^{1/2}$ can be done in $\mathcal{O}(k^{3})$ operations (see, e.g., sqrtm for the algorithm and the corresponding analysis). Finally, to compute ${\boldsymbol{y}}$, we have $k$ scalar products, a matrix vector multiplication with a $(k\times k)$ matrix and a matrix-matrix multiplication of $(N\times k)$ and $(k\times k)$ matrices, all of which can be done in $\mathcal{O}(Nk^{2})$ arithmetic operations.

To see (i), we employ the triangle inequality

[TABLE]

For the first term on the right-hand side, Lemma 6 below proves

[TABLE]

As shown in (16) of Lemma 2 below, the second term on the right-hand side of (10) is bounded by

[TABLE]

Hence, (i) follows from Proposition 1. For (ii), we note that the combination of both estimates in Proposition 2 below shows for $\boldsymbol{U}_{j}:=\boldsymbol{Q}_{j}^{T}\boldsymbol{M}\boldsymbol{Q}_{j}$

[TABLE]

We may eliminate the minimum in the error estimate since Algorithm 1 is essentially (up to roundoff errors) of Lanczos type, and for this algorithm, (monotone, , Example 5.1) shows that the approximation error $|\boldsymbol{C}_{p}^{1/2}{\boldsymbol{z}}-\boldsymbol{Q}_{j}(\boldsymbol{U}_{j}^{1/2})\boldsymbol{Q}_{j}^{T}{\boldsymbol{z}}|$ decreases monotonically in $j$. Since $\boldsymbol{Q}_{k_{0}}(\boldsymbol{U}_{k_{0}}^{1/2})\boldsymbol{Q}_{k_{0}}^{T}{\boldsymbol{z}}=\mathcal{Z}_{k,p}({\boldsymbol{z}})$, the remainder of the proof then follows as for (i) but we use (17) instead of (16) of Lemma 2 below.

3.3 Computing the square-root (Method 2)

The main drawback of Algorithm 1 is the additional storage requirements due to the necessity to store the matrix $\boldsymbol{Q}_{k}$. For this reason, we here follow a different approach, proposing a second algorithm that improves this situation.

The matrix sign function is defined for all square matrices $\widetilde{\boldsymbol{M}}$ with no pure imaginary eigenvalues as

[TABLE]

The sign function ${\rm sgn}(\widetilde{\boldsymbol{M}})$ can be computed using the Schultz iteration via

[TABLE]

The iterates $\boldsymbol{M}_{k}$ converge quadratically towards ${\rm sgn}(\widetilde{\boldsymbol{M}})$ if $\|\boldsymbol{I}-\widetilde{\boldsymbol{M}}^{2}\|_{2}<1$ in any matrix norm (see (sqit2, , Theorem 5.2)). It is observed in sqit1 , that all matrices $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ with only positive real eigenvalues satisfy

[TABLE]

where $\boldsymbol{I}\in{\mathbb{R}}^{N\times N}$ denotes the identity matrix, which opens the possibility to compute $\boldsymbol{M}^{1/2}$ via the sign function of the matrix By inserting

[TABLE]

By inserting this choice of $\widetilde{\boldsymbol{M}}$ into (12), we see that all iterates have the form

[TABLE]

As already observed in sqit1 , this leads to the iteration

[TABLE]

starting with $\boldsymbol{A}_{0}=\boldsymbol{M}$ and $\boldsymbol{B}_{0}=\boldsymbol{I}\in{\mathbb{R}}^{N\times N}$. The iterates $\boldsymbol{A}_{k}$ converge towards $\boldsymbol{M}^{1/2}$, which is what we aim to compute. The considerations above lead us to the following recursive form of the Schulz algorithm above, which uses only matrix vector multiplication. The subroutines PartA and PartB compute $\boldsymbol{A}_{k}{\boldsymbol{z}}$ and $\boldsymbol{B}_{k}{\boldsymbol{z}}$ respectively.

Algorithm 2

**Input: positive definite symmetric matrix $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$, vector ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$, maximal number of iterations $k\in{\mathbb{N}}$, temporary storage vectors ${\boldsymbol{z}}^{j}\in{\mathbb{R}}^{N}$, $j\in\{1,\ldots,k\}$, and scaling factor $0<s<2\|\boldsymbol{C}_{p}\|_{2}^{-1}$ (the scaling factor ensures convergence of the algorithm).

Main:**

Compute ${\boldsymbol{y}}=\texttt{PartA}(s\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j})_{j=1}^{k},k)$. 2. 2.

Return ${\boldsymbol{y}}/\sqrt{s}$.

**Output: the approximation ${\boldsymbol{y}}\approx\boldsymbol{M}^{1/2}{\boldsymbol{z}}$.

Subroutines:

$\texttt{PartA}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j}),k)$:**

• (i)

If $k=0$, return $\boldsymbol{M}{\boldsymbol{z}}$.

• (ii)

Compute ${\boldsymbol{z}}^{k}:=\texttt{PartA}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)$ and ${\boldsymbol{z}}^{k}:=\texttt{PartB}(\boldsymbol{M},{\boldsymbol{z}}^{k},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)$.

• (iii)

Compute ${\boldsymbol{z}}:=3{\boldsymbol{z}}-{\boldsymbol{z}}^{k}$.

• (iv)

Return $\frac{1}{2}\texttt{PartA}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)$.

$\texttt{PartB}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j}),k)$:

• (i)

If $k=0$, return ${\boldsymbol{z}}$.

• (ii)

Compute ${\boldsymbol{z}}^{k}:=\texttt{PartB}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)$ and ${\boldsymbol{z}}^{k}:=\texttt{PartA}(\boldsymbol{M},{\boldsymbol{z}}^{k},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)$.

• (iii)

Compute ${\boldsymbol{z}}:=3{\boldsymbol{z}}-{\boldsymbol{z}}^{k}$.

• (iv)

Return $\frac{1}{2}\texttt{PartB}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)$.

Remark 5

The extra storage vectors $({\boldsymbol{z}}^{j})_{j=1}^{k}$ are needed to avoid allocation of a new temporary storage vector in each call of either PartA are PartB. This would result in $\mathcal{O}(3^{k})$ additional allocations. By supplying the additional storage vectors, we can exploit the fact that each level of recursion can share a single storage vector.

Theorem 3.2

Suppose Assumption 1 holds and and let ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$. If $0<\eta<4c_{2}$ and $p$ is sufficiently large such that $\boldsymbol{C}_{p}$ constructed from $\boldsymbol{C}$ as in Section 3.1 is positive definite (condition (7) is sufficient), Algorithm 1 called with $\boldsymbol{M}=\boldsymbol{C}_{p}$ and $0<s<2\|\boldsymbol{C}_{p}\|_{2}^{-1}$ computes the approximation $\mathcal{Z}_{k,p}({\boldsymbol{z}}):={\boldsymbol{y}}\in{\mathbb{R}}^{N}$ such that

[TABLE]

where $\kappa:=\max\{|1-s\lambda_{\rm max}(\boldsymbol{C}_{p})|,|1-s\lambda_{\rm min}(\boldsymbol{C}_{p})|\}<1$. The algorithm completes in $\mathcal{O}(3^{k}p^{2d}N)$ arithmetic operations and uses less than $kN$ extra storage. The constant $C_{\rm H}$ is defined in Proposition 1.

Remark 6

In contrast to Algorithm 1 which needs $\mathcal{O}(|\log_{\kappa}(\varepsilon)|N)$ extra storage (at least in case (ii)), we see that Algorithm 2 requires only $\mathcal{O}(\log|\log(\varepsilon)|N)$ additional storage for an error request of $\varepsilon>0$.

Proof (Proof of Theorem 3.2)

First, we prove that PartA and PartB from Algorithm 2 correctly compute $\boldsymbol{A}_{k}{\boldsymbol{z}}$ and $\boldsymbol{B}_{k}{\boldsymbol{z}}$ from (13). This is done by induction on $k$. First, for $k=0$, the output of PartA is obviously $\boldsymbol{M}{\boldsymbol{z}}=\boldsymbol{A}_{0}{\boldsymbol{z}}$ and the output of PartB is ${\boldsymbol{z}}=\boldsymbol{B}_{0}{\boldsymbol{z}}$. This confirms the case $k=0$. Assume that PartA and PartB work correctly for $k\in{\mathbb{N}}$. By substitution of $\texttt{PartA}(\boldsymbol{M},{\boldsymbol{z}},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)=\boldsymbol{A}_{k-1}{\boldsymbol{z}}$ and $\texttt{PartB}(\boldsymbol{M},{\boldsymbol{z}}^{k},({\boldsymbol{z}}^{j})_{j=1}^{k-1},k-1)=\boldsymbol{B}_{k-1}{\boldsymbol{z}}^{k}$ in PartA, the variable ${\boldsymbol{z}}^{k}$ before step (iii) is given by ${\boldsymbol{z}}^{k}=\boldsymbol{B}_{k-1}\boldsymbol{A}_{k-1}{\boldsymbol{z}}$. Thus, step (iii)–(iv) correctly compute $\tfrac{1}{2}\boldsymbol{B}_{k-1}(3{\boldsymbol{z}}-\boldsymbol{B}_{k-1}\boldsymbol{A}_{k-1}{\boldsymbol{z}})=\boldsymbol{A}_{k}{\boldsymbol{z}}$. During the execution of $\texttt{PartA}(\cdot,\cdot,\cdot,k)$, extra storage vector ${\boldsymbol{z}}^{k}$ is not accessed by other instances of the subroutines (the function calls to $\texttt{PartA}(\cdot,\cdot,\cdot,k-1)$ and $\texttt{PartB}(\cdot,\cdot,\cdot,k-1)$ access only $({\boldsymbol{z}}^{k})_{j=1}^{k-1}$). This ensures that the correct value of ${\boldsymbol{z}}^{k}$ is used at each point of the execution. Analogously, we argue that PartB works correctly and thus conclude the induction.

For the computational cost estimate, we prove by induction that each subroutine $\texttt{PartA}(\cdot,\cdot,\cdot,k)$ and $\texttt{PartB}(\cdot,\cdot,\cdot,k)$ requires less than

[TABLE]

operations for some universal constant $C\geq 1$ and all $k\in{\mathbb{N}}$. For $k=0$, subroutine PartA performs an $H^{2}$-matrix-vector multiplication which, according to Proposition 1, costs less than $\mathcal{O}(p^{2d}N)$. Subroutine PartB just returns the vector ${\boldsymbol{z}}$. This shows (14) for $k=0$ for both subroutines. Assume that (14) is correct for both subroutines for some $k>0$. The fact that each subroutine $\texttt{PartA}(\cdot,\cdot,\cdot,k+1)$ and $\texttt{PartB}(\cdot,\cdot,\cdot,k+1)$ performs one scalar-vector multiplication and one vector addition as well as three calls to $\texttt{PartA}(\cdot,\cdot,\cdot,k)$ or $\texttt{PartB}(\cdot,\cdot,\cdot,k)$ shows that the cost of each subroutine $\texttt{PartA}(\cdot,\cdot,\cdot,k+1)$ and $\texttt{PartB}(\cdot,\cdot,\cdot,k+1)$ is bounded by

[TABLE]

This concludes the proof of (14), which proves the cost estimate since

[TABLE]

To see the error estimate, we use (10) and note that Algorithm 2 is nothing else than a recursive version of the iteration (13). The scaling $s<2\|\boldsymbol{C}_{p}\|_{2}^{-1}$ ensures $\kappa<1$, since $\lambda\in\{\lambda_{\rm min}(\boldsymbol{C}_{p}),\lambda_{\rm max}(\boldsymbol{C}_{p})\}$ satisfies $1-s\lambda<1$ (since $s,\lambda>0$) as well as $s\lambda-1\leq s\|\boldsymbol{C}_{p}\|_{2}-1<2-1=1$. Thus, Lemma 8 shows

[TABLE]

We conclude the proof with the aid of (11) and Proposition 1.

4 Numerical experiments

All numerical experiments where computed in Matlab, by use of a Matlab-$H^{2}$-matrix library which can be downloaded under software.michaelfeischl.net. The authors are well aware that the Matlab implementation prohibits high-end performance. However, we wanted to demonstrate the feasibility of our algorithms and show the correct convergence rates, for which purpose the Matlab implementation is sufficient.

For the first example, we consider a covariance function of the form (3) with

[TABLE]

We use Algorithm 1 to generate six samples on the unit square $D=[0,1]^{2}$ of the corresponding normal random field $\mathcal{Z}$ shown in Figure 1. Figure 2–3 show samples of the covariance functions from (2) with different parameters.

To illustrate the challenging nature of handling these covariance matrices, Table 1 shows condition numbers of $\boldsymbol{C}$ for different problem sizes and the Matérn covariance function (2).

For a performance comparison of Algorithm 1 and Algorithm 12, we consider the covariance function of the form (2) with $p=2$, $\sigma=1$, and varying $\mu\in\{1/2,\infty\}$, $\lambda\in\{1,10^{-1},10^{-2},10^{-3}\}$. We compute samples of $\mathcal{Z}({\boldsymbol{x}},\omega)$ on a Sobol pointset with $2^{10}$ points. The results are plotted in Figure 4 where we see the relative approximation error versus the computation time in seconds. We observe that with respect to computational time, Algorithm 1 is superior in almost all cases (particularly for smooth fields). However, keep in mind that according to Theorem 3.1, Algorithm 1 needs up to $\mathcal{O}(\log_{\kappa}(\varepsilon)N)$ extra storage, while Algorithm 2 uses only $\mathcal{O}(\log(\log(\varepsilon))N)$ extra storage units. (See Theorem 3.2, where the quadratic convergence shows that $k\simeq\log(\log(\varepsilon))$ is sufficient to reach a given accuracy $\varepsilon>0$. However, we have to mention that $k$ iterations of Algorithm 2 require $\mathcal{O}(3^{k})$ arithmetic operations.)

Figure 5 compares the two algorithms with the direct matrix square root provided by Matlab. We evaluate $\mathcal{Z}({\boldsymbol{x}},\omega)$ on a Sobol pointset with size $2^{m}$ for $m\in\{1,\ldots,14\}$. The number of iterations in both algorithm is set such that the relative error is smaller than $10^{-10}$ for the example from above with $p=2$, and varying $\mu\in\{1/2,\infty\}$, $\lambda\in\{1,10^{-1},10^{-2},10^{-3}\}$. We see that both, Algorithm 1–2, perform in linear time, whereas the direct approach comes closer to $\mathcal{O}(N^{3})$. Even though our $H^{2}$-matrix library is programmed entirely in Matlab (and thus nowhere near optimal performance), the breakthrough point at around $N=10^{3}$ shows that also small problems benefit from the speed up.

5 Lemmas for the proof of Theorem 3.1

First, we state a slight generalization of a well-known result.

Lemma 2

Let $\boldsymbol{A},\boldsymbol{B}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite. Then, there holds

[TABLE]

as well as

[TABLE]

Proof

The estimate (16) is proved in (matrixsq, , Lemma 2.2). To obtain (17), let $\boldsymbol{U}\in{\mathbb{R}}^{N\times N}$ denote the orthonormal matrix that diagonalizes $\boldsymbol{A}$, i.e., $\boldsymbol{U}^{T}\boldsymbol{A}\boldsymbol{U}=\boldsymbol{D}$ for a positive diagonal matrix $\boldsymbol{D}\in{\mathbb{R}}^{N\times N}$. With $\boldsymbol{U}\boldsymbol{D}^{1/2}\boldsymbol{U}^{T}=\boldsymbol{A}^{1/2}$ and $\boldsymbol{U}\boldsymbol{U}^{T}=\boldsymbol{I}$, there holds for arbitrary $\alpha\geq 0$

[TABLE]

where we used $x+y\leq(\sqrt{x}+\sqrt{y})^{2}$ and hence $\sqrt{x+y}\leq\sqrt{x}+\sqrt{y}$ for $x,y\geq 0$ in the last estimate. With $\alpha:=\|\boldsymbol{A}-\boldsymbol{B}\|_{2}$, (16) shows

[TABLE]

The combination of the last two estimates concludes the proof of (17).

Lemma 3

Let $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite and assume that $0<\kappa<1$ and $C_{\kappa}>0$ are such that the sequence of all distinct eigenvalues $\lambda_{1}>\ldots>\lambda_{M}>0\in{\mathbb{R}}$ (for some $M\leq N$) of $\boldsymbol{M}$ satisfies $|\lambda_{i}-\lambda_{j}|\leq\lambda_{1}C_{\kappa}\kappa^{\min\{i,j\}}$ for all $1\leq i,j\leq M$. Given $1\leq k\leq M$ and ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$, define $\boldsymbol{Z}\in{\mathbb{R}}^{N\times k}$ by

[TABLE]

Consider the $QR$-factorization $\boldsymbol{Z}=\boldsymbol{Q}\boldsymbol{R}$, with $\boldsymbol{Q}\in{\mathbb{R}}^{N\times k}$ satisfying $\boldsymbol{Q}^{T}\boldsymbol{Q}=\boldsymbol{I}_{k}$ and $\boldsymbol{R}\in{\mathbb{R}}^{k\times k}$ upper triangular with non-negative diagonal entries (note that if $\boldsymbol{Z}$ has full rank, this ensures uniqueness of $\boldsymbol{Q}$ and $\boldsymbol{R}$). Then the diagonal entries of $\boldsymbol{R}$ satisfy

[TABLE]

Proof

Let ${\boldsymbol{q}}^{i}$, $i=1,\ldots,k$ denote the orthonormal columns of $\boldsymbol{Q}$. By definition of the $QR$-factorization, there holds for $1\leq n\leq k$

[TABLE]

Since the ${\boldsymbol{q}}^{i}$ are orthogonal, the best approximation (with respect to $|\cdot|$) of $\lambda_{1}^{-(n-1)}\boldsymbol{M}^{n-1}{\boldsymbol{z}}$ in ${\rm span}\{{\boldsymbol{q}}^{1},\ldots,{\boldsymbol{q}}^{n-1}\}$ is given by $\sum_{i=1}^{n-1}\boldsymbol{R}_{in}{\boldsymbol{q}}^{i}$ for all $n\geq 2$. Therefore, we obtain

[TABLE]

where we used ${\rm span}\{{\boldsymbol{z}},\ldots,\boldsymbol{M}^{n-2}{\boldsymbol{z}}\}\subseteq{\rm span}\{{\boldsymbol{q}}^{1},\ldots,{\boldsymbol{q}}^{n-1}\}$ by definition of the $QR$-factorization (see also Remark 1). We may choose ${\boldsymbol{v}}=p(\boldsymbol{M}){\boldsymbol{z}}$, where $p(x)$ is the polynomial of degree $n-2$ interpolating $f(x):=(x/\lambda_{1})^{n-1}$ at the points $x=\lambda_{1},\ldots,\lambda_{n-1}$. Since $\boldsymbol{M}$ is symmetric and positive definite, we may diagonalize it with an orthogonal matrix $\boldsymbol{U}\in{\mathbb{R}}^{N\times N}$, i.e., $\boldsymbol{U}^{T}\boldsymbol{M}\boldsymbol{U}=\boldsymbol{D}$ with a diagonal matrix $\boldsymbol{D}\in{\mathbb{R}}^{N\times N}$ containing the eigenvalues of $\boldsymbol{M}$. This allows us to conclude

[TABLE]

The function $f(x)-p(x)$ is a polynomial of degree $n-1$ with known zeros $\lambda_{1},\ldots,\lambda_{n-1}$ and thus reads

[TABLE]

for some leading coefficient $\alpha\in{\mathbb{R}}$. Differentiation reveals $\alpha(n-1)!=f^{(n-1)}(x)=(n-1)!\lambda_{1}^{-(n-1)}$ and hence $\alpha=\lambda_{1}^{-(n-1)}$. This shows

[TABLE]

By the decay assumption on the $\lambda_{i}$ it follows that

[TABLE]

This concludes the proof.

The next lemma shows that the matrices $\boldsymbol{Q}_{j}$ from Algorithm 1 are strongly tied to the matrices $\boldsymbol{Z}=\boldsymbol{Q}\boldsymbol{R}$ defined in Lemma 3.

Lemma 4

Given ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$ and let $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite. Call Algorithm 1 with $\boldsymbol{M}$, ${\boldsymbol{z}}$, and $k\in{\mathbb{N}}$ to compute $k_{0}\leq k$ and $\boldsymbol{R}_{j}$, $\boldsymbol{Q}_{j}$ for all $1\leq j\leq k_{0}$. Define $\boldsymbol{Z},\boldsymbol{Q},\boldsymbol{R}$ satisfying $\boldsymbol{Z}=\boldsymbol{Q}\boldsymbol{R}$ as in Lemma 3. Then, $\boldsymbol{Q}_{j}$ (as defined in Algorithm 1) for $1\leq j\leq k_{0}$ satisfies $\boldsymbol{Q}_{j}=\boldsymbol{Q}|_{\{1,\ldots,N\}\times\{1,\ldots,j\}}$, i.e., the first $j$ columns coincide and

[TABLE]

for all $1\leq j\leq k_{0}$. Moreover, $\boldsymbol{Z}$ has full rank if and only if $k_{0}=k$.

Proof

Let ${\boldsymbol{q}}^{j}$ denote the $j$-th column of $\boldsymbol{Q}_{j}$ and note that by definition of Algorithm 1 we have

[TABLE]

In order to prove (21), we first show

[TABLE]

for all $1\leq j\leq k_{0}$ by induction. To that end, note that $\boldsymbol{Q}_{1}={\boldsymbol{q}}^{1}={\boldsymbol{z}}/|{\boldsymbol{z}}|$ and consequently (23) holds for $j=1$. Assume (23) holds for all $1\leq j<j_{0}\leq k_{0}$. By construction of the matrices in Algorithm 1, we have

[TABLE]

By the induction assumption, ${\boldsymbol{q}}^{j_{0}-1}\in{\rm span}\{{\boldsymbol{z}},\ldots,\boldsymbol{M}^{j_{0}-2}{\boldsymbol{z}}\}$. Thus, (24) and the fact that $\boldsymbol{R}_{j_{0}}$ is regular (by (22)) imply

[TABLE]

The fact that $\boldsymbol{Q}_{j_{0}}$ is orthogonal (and hence its range is $j_{0}$ dimensional) shows even equality, that is

[TABLE]

This concludes the induction, and proves (23) for all $1\leq j\leq k_{0}$. The second equation in (21) follows by definition of $\boldsymbol{Z}$.

To see the remainder of the statement, we first assume $k_{0}=k$ and proceed to prove that $\boldsymbol{Z}$ has full rank. To that end, we apply (21) with $j=k$ to see that ${\rm range}(\boldsymbol{Z})={\rm range}(\boldsymbol{Q}_{k})$ is $k$-dimensional and therefore $\boldsymbol{Z}$ has full rank.

For the converse implication, assume that $\boldsymbol{Z}$ has full rank. We prove $k_{0}=k$ by induction. By construction, we have $(\boldsymbol{R}_{1})_{11}=1$ and thus $k_{0}\geq 1$. Assume $k_{0}\geq j_{0}$ for some $j_{0}<k$. Then, since $(\boldsymbol{R}_{j})_{jj}\neq 0$ for all $1\leq j<j_{0}$, the identity (21) shows ${\rm range}(\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,j\}})={\rm range}(\boldsymbol{Q}_{j})$ for all $j<j_{0}$. From this, we argue that

[TABLE]

which, by definition of $\boldsymbol{Z}=({\boldsymbol{z}},\lambda_{1}^{-1}\boldsymbol{M}{\boldsymbol{z}},\ldots,\lambda_{1}^{-(k-1)}\boldsymbol{M}^{k-1}{\boldsymbol{z}})$, shows that ${\boldsymbol{q}}^{j_{0}-1}=\sum_{i=0}^{j_{0}-2}\alpha_{i}\boldsymbol{M}^{i}{\boldsymbol{z}}$ for some $\alpha_{i}\in{\mathbb{R}}$ with $\alpha_{j_{0}-2}\neq 0$. Consequently, we obtain $\boldsymbol{M}{\boldsymbol{q}}^{j_{0}-1}=\sum_{i=0}^{j_{0}-2}\alpha_{i}\boldsymbol{M}^{i+1}{\boldsymbol{z}}\in{\rm range}(\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,j_{0}\}})\setminus{\rm range}(\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,j_{0}-1\}})$. Since ${\rm range}(\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,j_{0}-1\}})={\rm range}(\boldsymbol{Q}_{j_{0}-1})$, this implies the identity ${\rm range}((\boldsymbol{Q}_{j_{0}-1},\boldsymbol{M}{\boldsymbol{q}}^{j_{0}-1}))={\rm range}(\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,j_{0}\}})$ and therefore the matrix $(\boldsymbol{Q}_{j_{0}-1},\boldsymbol{M}{\boldsymbol{q}}^{j_{0}-1})$ has full rank. Hence, (24) implies that $\boldsymbol{R}_{j_{0}}$ has full rank, which in particular implies $(\boldsymbol{R}_{j_{0}})_{j_{0}j_{0}}\neq 0$ and thus $k_{0}\geq j_{0}+1$. This concludes the induction and shows $k_{0}=k$.

The following result proves that if Algorithm 1 terminates in less than $k$ steps (due to the criterion in step 1(c)), the quantity $\boldsymbol{M}^{1/2}{\boldsymbol{z}}$ is computed exactly.

Lemma 5

Let ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$ and let $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite. Call Algorithm 1 with $\boldsymbol{M}$, ${\boldsymbol{z}}$, and $k\in{\mathbb{N}}$ to compute $k_{0}\leq k$ as well as $\boldsymbol{Q}_{j}$ for all $1\leq j\leq k_{0}$. Define $\boldsymbol{U}_{k_{0}}=\boldsymbol{Q}_{k_{0}}^{T}\boldsymbol{M}\boldsymbol{Q}_{k_{0}}$ as in Algorithm 1. If $k_{0}<k$, there holds

[TABLE]

Proof

If $k_{0}<k$ then Lemma 4 shows that $\boldsymbol{Z}$ as defined in Lemma 3 does not have full rank. Moreover, the identity  (21) shows that $\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,k_{0}\}}$ has full rank. By definition of $\boldsymbol{Z}$, this implies ${\rm range}(\boldsymbol{Z}|_{\{1,\ldots,N\}\times\{1,\ldots,k_{0}\}})={\rm range}(\boldsymbol{Z})$. Therefore, (21) shows

[TABLE]

Let $\overline{\boldsymbol{Q}}\in{\mathbb{R}}^{N\times N}$ be an orthonormal matrix such that its first $k_{0}$ columns coincide with $\boldsymbol{Q}_{k_{0}}$, i.e., $\overline{\boldsymbol{Q}}=(\boldsymbol{Q}_{k_{0}},\boldsymbol{Q}_{\perp})$ for some orthonormal $\boldsymbol{Q}_{\perp}\in{\mathbb{R}}^{N\times(N-k_{0})}$. We obtain

[TABLE]

There holds

[TABLE]

The invariance property (26) shows $\boldsymbol{Q}_{\perp}^{T}\boldsymbol{M}\boldsymbol{Q}_{k_{0}}=\boldsymbol{0}$, and by symmetry also $\boldsymbol{Q}_{k_{0}}^{T}\boldsymbol{M}\boldsymbol{Q}_{\perp}=\boldsymbol{0}$. Therefore, we have

[TABLE]

This and (27), together with ${\boldsymbol{z}}\in{\rm range}(\boldsymbol{Q}_{k_{0}})$, show $\boldsymbol{M}^{1/2}{\boldsymbol{z}}=\overline{\boldsymbol{Q}}\,(\overline{\boldsymbol{Q}}^{T}\boldsymbol{M}\overline{\boldsymbol{Q}})^{1/2}\,\overline{\boldsymbol{Q}}^{T}{\boldsymbol{z}}=\boldsymbol{Q}_{k_{0}}\boldsymbol{U}_{k_{0}}^{1/2}\boldsymbol{Q}_{k_{0}}^{T}{\boldsymbol{z}}$ and conclude the proof.

The following result is the main tool to prove Theorem 3.1 (i).

Lemma 6

Let ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$ and let $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite. Call Algorithm 1 with $\boldsymbol{M}$, ${\boldsymbol{z}}$, and $k\in{\mathbb{N}}$ to compute $k_{0}\leq k$ as well as $\boldsymbol{Q}_{j}$ for all $1\leq j\leq k_{0}$. Let $\boldsymbol{U}_{k_{0}}=\boldsymbol{Q}_{k_{0}}^{T}\boldsymbol{M}\boldsymbol{Q}_{k_{0}}$ be defined as in Algorithm 1. Then, there holds

[TABLE]

where

[TABLE]

Proof

The case $k_{0}<k$ is covered in Lemma 5. Assume $k_{0}=k$. Note that $\boldsymbol{Q}_{k}\boldsymbol{Q}_{k}^{T}$ is the identity on ${\rm range}(\boldsymbol{Q}_{k})$. Lemma 4 shows that $\boldsymbol{M}^{j}{\boldsymbol{z}}\in{\rm range}(\boldsymbol{Q}_{k})$ for all $0\leq j\leq k-1$. Moreover, ${\boldsymbol{z}}\in{\rm range}(\boldsymbol{Q}_{k})$ by construction. Hence, we have

[TABLE]

Thus, any polynomial $p\in{\mathcal{P}}^{k-1}$ of degree $k-1$ satisfies

[TABLE]

This implies for all $p\in{\mathcal{P}}^{k-1}$

[TABLE]

With $f(x):=\sqrt{(x+1)(\lambda_{\rm max}(\boldsymbol{M})-\lambda_{\rm min}(\boldsymbol{M}))/2+\lambda_{\rm min}(\boldsymbol{M})}$, the result (h2mat, , Lemma 4.14) proves

[TABLE]

with $r>1$ from (28) and

[TABLE]

Since $x\in{\mathbb{C}}_{r}$ implies $|x|\leq r$, straightforward calculations show

[TABLE]

which implies the estimate $\min_{p\in{\mathcal{P}}^{k-1}}\|f-p\|_{L^{\infty}([-1,1])}\leq{\sqrt{2\|\boldsymbol{M}\|_{2}}}\frac{2r^{{2}}}{r-1}r^{-k}$. Hence, we obtain for $g(x):=\sqrt{x}$ (note that $x\mapsto(x+1)(\lambda_{\rm max}(\boldsymbol{M})-\lambda_{\rm min}(\boldsymbol{M}))/2+\lambda_{\rm min}(\boldsymbol{M})$ maps $[-1,1]$ onto $[\lambda_{\rm min}(\boldsymbol{M}),\lambda_{\rm max}(\boldsymbol{M})]$) also

[TABLE]

Let $\boldsymbol{U}\in{\mathbb{R}}^{N\times N}$ denote the orthonormal matrix ($\boldsymbol{U}\boldsymbol{U}^{T}=\boldsymbol{I}$) that diagonalizes $\boldsymbol{M}$, i.e., $\boldsymbol{M}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{T}$ with $\boldsymbol{D}\in{\mathbb{R}}^{N\times N}$ being the diagonal matrix containing the eigenvalues of $\boldsymbol{M}$. There holds $\boldsymbol{M}^{1/2}=\boldsymbol{U}\boldsymbol{D}^{1/2}\boldsymbol{U}^{T}$ as well as $p(\boldsymbol{M})=\boldsymbol{U}p(\boldsymbol{D})\boldsymbol{U}^{T}$. This, (30), and invariance of the spectral norm $\|\cdot\|_{2}=\|\boldsymbol{U}(\cdot)\boldsymbol{U}^{T}\|_{2}$ show

[TABLE]

Since $\boldsymbol{U}_{k}$ is an orthogonal projection of $\boldsymbol{M}$, we have $\lambda_{\rm min}(\boldsymbol{M})\leq\lambda_{\rm min}(\boldsymbol{U}_{k})\leq\lambda_{\rm max}(\boldsymbol{U}_{k})\leq\lambda_{\rm max}(\boldsymbol{M})$. Thus, repeating the above argument for $\boldsymbol{U}_{k}$ instead of $\boldsymbol{M}$ yields

[TABLE]

This in combination with (29) and Lemma 4 conclude the proof.

The next result quantifies the distance of ${\rm range}(\boldsymbol{Q}_{j})$ to ${\rm range}(\boldsymbol{M}\boldsymbol{Q}_{j})$ in terms of the projection $\boldsymbol{Q}_{j}\boldsymbol{Q}_{j}^{T}$ onto ${\rm range}(\boldsymbol{Q}_{j})$.

Lemma 7

Assume the requirements of Lemma 3. Call Algorithm 1 with $\boldsymbol{M}$, ${\boldsymbol{z}}$, and $k\in{\mathbb{N}}$ to compute $k_{0}\leq k$ as well as $\boldsymbol{Q}_{j}$ for all $1\leq j\leq k_{0}$. Let ${\boldsymbol{q}}^{j}$ be the last column of $\boldsymbol{Q}_{j}$ for all $1\leq j\leq k_{0}$. There holds for all $1\leq j<k_{0}$

[TABLE]

as well as

[TABLE]

Proof

Recall $\boldsymbol{Z},\boldsymbol{Q},\boldsymbol{R}$ satisfying $\boldsymbol{Z}=\boldsymbol{Q}\boldsymbol{R}$ from Lemma 3 with $k$ replaced by $k_{0}$ in the call to Algorithm 1. By Lemma 4, ${\boldsymbol{q}}^{j}$ coincides with the $j$-th column of $\boldsymbol{Q}$ for all $1\leq j\leq k_{0}$. Moreover, let ${\boldsymbol{r}}^{j}$ be the $j$-th column of $\boldsymbol{R}$ and let $\widetilde{{\boldsymbol{r}}}^{j}$ be the $j$-th column of $\boldsymbol{R}^{-1}$ from Lemma 3. All quantities are well-defined since $\boldsymbol{Z}$ has maximal rank $k_{0}$ by Lemma 4. For the first statement (31), note that ${\rm range}(\boldsymbol{M}\boldsymbol{Q}_{j})\subseteq{\rm range}(\boldsymbol{Q}_{j+1})$ implies $\boldsymbol{M}\boldsymbol{Q}_{j}=\boldsymbol{Q}_{j+1}\boldsymbol{Q}_{j+1}^{T}\boldsymbol{M}\boldsymbol{Q}_{j}$. Moreover, due to Lemma 4, we have $\boldsymbol{Q}_{j+1}=(\boldsymbol{Q}_{j},{\boldsymbol{q}}^{j+1})$, and hence $\boldsymbol{Q}_{j+1}\boldsymbol{Q}_{j+1}^{T}={\boldsymbol{q}}^{j+1}({\boldsymbol{q}}^{j+1})^{T}+\boldsymbol{Q}_{j}\boldsymbol{Q}_{j}^{T}$. Altogether, this shows

[TABLE]

where the last step follows because ${\boldsymbol{q}}^{j+1}$ is orthogonal to $\boldsymbol{M}{\boldsymbol{q}}^{i}$, $i=1,\ldots,j-1$. This proves (31).

To see the remaining statement, note that the definition of $\boldsymbol{Z}$ in (18) implies

[TABLE]

as well as

[TABLE]

The last two identities, and the fact that $({\boldsymbol{q}}^{j+1})^{T}(\boldsymbol{M}\boldsymbol{Z})|_{\{1,\ldots,N\}\times\{i\}}=0$ for all $1\leq i\leq j-1$, imply

[TABLE]

The triangular structure of $\boldsymbol{R}$ implies $(\boldsymbol{R}^{-1})_{jj}=1/\boldsymbol{R}_{jj}$ and hence $(\widetilde{\boldsymbol{r}}^{j})_{j}=1/\boldsymbol{R}_{jj}$ (where $\boldsymbol{R}_{jj}\neq 0$ by assumption). This shows

[TABLE]

With Lemma 3, we have

[TABLE]

Moreover, we know $\boldsymbol{R}_{11}=|{\boldsymbol{q}}^{1}\boldsymbol{R}_{11}|=|{\boldsymbol{z}}|$. This implies that at least one of the fractions on the left-hand side of (33) must be smaller than the $j$-th root of the right hand side of (33) divided by $|{\boldsymbol{z}}|$ and hence

[TABLE]

With this, (32), and (31), we obtain

[TABLE]

This concludes the proof.

The following proposition is the main tool to prove Theorem 3.1 (ii).

Proposition 2

Let ${\boldsymbol{z}}\in{\mathbb{R}}^{N}$ and let $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite. Call Algorithm 1 with $\boldsymbol{M}$, ${\boldsymbol{z}}$, and $k\in{\mathbb{N}}$ to compute $k_{0}\leq k$ as well as $\boldsymbol{Q}_{j}$ for all $1\leq j\leq k_{0}$. Then, $\boldsymbol{U}_{j}:=\boldsymbol{Q}_{j}^{T}\boldsymbol{M}\boldsymbol{Q}_{j}$ satisfies the error bound

[TABLE]

and we have the a priori estimate

[TABLE]

for all $1\leq j<k_{0}$.

Proof

The case $k_{0}<k$ and $j=k_{0}$ is trivially covered in Lemma 5. For the other cases, let $\overline{\boldsymbol{Q}}\in{\mathbb{R}}^{N\times N}$ be orthonormal such that the first $j$ columns coincide with $\boldsymbol{Q}_{j}$, i.e., $\overline{\boldsymbol{Q}}=(\boldsymbol{Q}_{j},\boldsymbol{Q}_{\perp})$ for some orthonormal $\boldsymbol{Q}_{\perp}\in{\mathbb{R}}^{N\times(N-j)}$. Then, we write

[TABLE]

for matrices $\boldsymbol{S}=\boldsymbol{Q}_{\perp}^{T}\boldsymbol{M}\boldsymbol{Q}_{j}\in{\mathbb{R}}^{(N-j)\times j}$, $\boldsymbol{T}\in{\mathbb{R}}^{(N-j)\times(N-j)}$. This means that

[TABLE]

Lemma 2 then implies

[TABLE]

Since $\boldsymbol{I}-\boldsymbol{Q}_{j}\boldsymbol{Q}_{j}^{T}=\boldsymbol{Q}_{\perp}\boldsymbol{Q}_{\perp}^{T}$, we have

[TABLE]

With $(\overline{\boldsymbol{Q}}^{T}\boldsymbol{M}\overline{\boldsymbol{Q}})^{1/2}=\overline{\boldsymbol{Q}}^{T}\boldsymbol{M}^{1/2}\overline{\boldsymbol{Q}}$ and since the ranges of $\boldsymbol{Q}_{j}$ and $\boldsymbol{Q}_{\perp}$ are orthogonal, we have $\boldsymbol{Q}_{\perp}^{T}\boldsymbol{Q}_{j}\boldsymbol{Q}_{j}^{T}=0$ and

[TABLE]

The combination of (34) and (35) shows

[TABLE]

We conclude the proof with ${\boldsymbol{z}}=\boldsymbol{Q}_{j}\boldsymbol{Q}_{j}^{T}{\boldsymbol{z}}$ due to ${\boldsymbol{z}}\in{\rm range}(\boldsymbol{Q}_{j})$ and Lemma 7.

6 Lemma for the proof of Theorem 3.2

The following lemma is the main tool for the proof of Theorem 3.2.

Lemma 8

Let $\boldsymbol{M}\in{\mathbb{R}}^{N\times N}$ be symmetric positive definite. Then, the iteration (13) with initial values $\boldsymbol{A}_{0}=s\boldsymbol{M}$ and $\boldsymbol{B}_{0}=\boldsymbol{I}$ satisfies

[TABLE]

for all $k\in{\mathbb{N}}$ and all $s>0$. The minimum bound is attained at $s=2/(\lambda_{\rm min}(\boldsymbol{M})+\lambda_{\rm max}(\boldsymbol{M}))$ such that $\max\{|1-s\lambda_{\rm max}(\boldsymbol{M})|,|1-s\lambda_{\rm min}(\boldsymbol{M})|\}=1-2\lambda_{\rm min}(\boldsymbol{M})/(\lambda_{\rm min}(\boldsymbol{M})+\lambda_{\rm max}(\boldsymbol{M}))$.

Proof

Straightforward calculations show

[TABLE]

The result (sqit2, , Theorem 5.2) shows $\|\boldsymbol{I}-\boldsymbol{X}_{n}^{2}\|_{2}<\|\boldsymbol{I}-\boldsymbol{X}_{0}^{2}\|_{2}^{(e_{1}+e_{2}+1)^{n}}$ for all $n\in{\mathbb{N}}$, where $\boldsymbol{X}_{n+1}=-\boldsymbol{X}_{n}P_{e_{1}e_{2}}(\boldsymbol{I}-\boldsymbol{X}_{n}^{2})Q_{e_{1}e_{2}}^{-1}(\boldsymbol{I}-\boldsymbol{X}_{n}^{2})$ and $\boldsymbol{X}_{0}$ has no purely imaginary eigenvalues. Here $P_{e_{1}e_{2}}/Q_{e_{1}e_{2}}$ is the $(e_{1}/e_{2})$-Padé approximant to $(1-x)^{-1/2}$. We obtain from (sqit2, , Table 1) that for $e_{1}=1$ and $e_{2}=0$, $\boldsymbol{X}_{n}$ satisfies the Schultz iteration (12) and thus we may use the result with

[TABLE]

to show

[TABLE]

for all $k\in{\mathbb{N}}$. By scaling of $\boldsymbol{M}$, we may minimize the right-hand side. To that end, we observe that the spectrum satisfies $\sigma(\boldsymbol{I}-s\boldsymbol{M})\subset[1-s\lambda_{\rm max}(\boldsymbol{M}),1-s\lambda_{\rm min}(\boldsymbol{M})]$. The fact $\|\boldsymbol{I}-s\boldsymbol{M}\|_{2}\leq\max\{|1-s\lambda_{\rm max}(\boldsymbol{M})|,|1-s\lambda_{\rm min}(\boldsymbol{M})|\}$ proves (36). A straightforward optimization of $s>0$ concludes the proof.

Appendix A Proof of Lemma 1

The following lemma is an elementary statement on holomorphic functions

Lemma 9

Let $f\colon O\to{\mathbb{C}}$ be a continuous function on the domain $O\subset{\mathbb{C}}^{n}$ which is holomorphic in $O$ in all variables ${\boldsymbol{x}}_{i}$, $i\in\{1,\ldots,n\}$, i.e.,

[TABLE]

is holomorphic in \big{\{}{\boldsymbol{x}}_{i}\in{\mathbb{C}}\,:\,({\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{i},\ldots,{\boldsymbol{x}}_{n})\in O\big{\}} for all ${\boldsymbol{x}}_{1},\ldots,{\boldsymbol{x}}_{i-1},{\boldsymbol{x}}_{i+1},\ldots,{\boldsymbol{x}}_{n}\in{\mathbb{C}}$. Then, for all multi-indices $\alpha\in{\mathbb{N}}_{0}^{n}$, the function $\partial_{\boldsymbol{x}}^{\alpha}f$ is holomorphic in $O$ in all variables ${\boldsymbol{x}}_{i}$, $i\in\{1,\ldots,n\}$ as defined above.

Proof

The result is proved by induction on $|\alpha|_{1}$. Obviously, for $|\alpha|_{1}=0$, $\partial_{\boldsymbol{x}}^{\alpha}f=f$ and the statement is true. Assume the statement holds for all $|\alpha|_{1}\leq k$ and choose some $\alpha\in{\mathbb{N}}_{0}^{n}$ with $|\alpha|_{1}=k+1$. Then, we have for some $i\in\{1,\ldots,n\}$ and some $\alpha_{0}\in{\mathbb{N}}_{0}^{n}$ with $|\alpha_{0}|_{1}=k$ that

[TABLE]

Since, $\partial_{\boldsymbol{x}}^{\alpha_{0}}f$ is holomorphic in $O$ in all variables by the induction hypothesis, obviously $\partial_{\boldsymbol{x}}^{\alpha}f$ is holomorphic in $O$ at least in ${\boldsymbol{x}}_{i}$ (derivatives of holomorphic functions are holomorphic). To prove the statement for all other variables, we may employ Cauchy’s integral formula to obtain

[TABLE]

for some $\varepsilon>0$ with $B_{\varepsilon}({\boldsymbol{x}}_{i})\subset{\mathbb{C}}$ being the ball with radius $\varepsilon$. The integrand is holomorphic in all variables ${\boldsymbol{x}}_{j}$, $j\neq i$. Hence, we conclude that $\partial_{\boldsymbol{x}}^{\alpha}f({\boldsymbol{x}})$ is holomorphic in all variables and prove the assertion.

The following result is elementary but technical.

Lemma 10

For $n,p\in{\mathbb{N}}$, define the set M:=\big{\{}{\boldsymbol{x}}\in{\mathbb{C}}^{n}\,:\,{\rm real}(\sum_{i=1}^{n}{\boldsymbol{x}}_{i}^{p})\leq 0\big{\}}. Then, there holds ({\mathbb{R}}^{n})_{+}:=\big{\{}{\boldsymbol{x}}\in{\mathbb{R}}^{n}\setminus\{0\}\,:\,{\boldsymbol{x}}_{i}\geq 0\big{\}}\cap M=\emptyset and

[TABLE]

Proof

Let ${\boldsymbol{x}}\in({\mathbb{R}}^{n})_{+}$, then we have $\sum_{i=1}^{n}{\boldsymbol{x}}_{i}^{p}>0$ and hence ${\boldsymbol{x}}\notin M$. It is easy to see that the cone C_{p}:=\big{\{}r\exp(i\phi)\,:\,r>0,\,\phi\in(-\frac{\pi}{2p},\frac{\pi}{2p})\big{\}}\subset{\mathbb{C}} satisfies ${\rm real}(x^{p})>0$ for all $x\in C_{p}$. Thus, we have that

[TABLE]

satisfies $C_{p}^{n}\cap M=\emptyset$.

Moreover, a simple geometric argument (see Figure 6) shows that all $x>0$ satisfy

[TABLE]

Since $({\mathbb{R}}^{n})_{+}\subseteq C_{p}^{n}$, this implies

[TABLE]

This concludes the proof.

Products of asymptotically smooth functions are again asymptotically smooth. This is shown in the next lemma.

Lemma 11

Given two functions $f,g\colon D\times D\to{\mathbb{R}}$ which are asymptotically smooth (1). Then, also their product $fg$ satisfies (1).

Proof

To simplify the notation, we consider $f,g$ as functions of one variable ${\boldsymbol{z}}=({\boldsymbol{x}},{\boldsymbol{y}})\in D\times D\subset{\mathbb{R}}^{2d}$. For multi-indices $\alpha,\beta\in{\mathbb{N}}^{2d}$, define

[TABLE]

Note that there holds $\binom{\alpha}{\beta}\leq\binom{|\alpha|_{1}}{|\beta|_{1}}$. This follows from the basic combinatorial fact that the number of possible choices of $\beta_{i}$ elements out of a set of $\alpha_{i}$ elements for all $i=1,\ldots,2d$ is smaller than the number of choices of $|\beta|_{1}$ elements out of a set of $|\alpha|_{1}$ elements.

The Leibniz formula together with the definition of asymptotically smooth function (1) show for $\alpha\in{\mathbb{N}}^{2d}$

[TABLE]

where we used $(|\alpha|_{1}+1)^{2d}\leq(2d\exp(2d))^{|\alpha|_{1}}$ and $\widetilde{c}_{2}=c_{2}/(2d\exp(2d))$. This concludes the proof.

The final lemma of this section proves the concatenations of certain asymptotically smooth functions are asymptotically smooth.

Lemma 12

Let $g\colon D\times D\to{\mathbb{R}}$ be asymptotically smooth (1) with constants $c_{1},c_{2}>0$.

• (i)

If $c_{g}:=\sup_{{\boldsymbol{x}}\in D\times D}g({\boldsymbol{x}})<\infty$. Then, $\exp\circ g$ satisfies (1) with constants $\widetilde{c}_{1}:=\exp(c_{g})$ and $\widetilde{c}_{2}:=c_{2}/(2\max\{1,c_{1}\})$.

• (ii)

If $g$ satisfies $\partial_{\boldsymbol{x}}^{\alpha}\partial_{\boldsymbol{y}}^{\alpha}g({\boldsymbol{x}},{\boldsymbol{y}})\leq C_{g}$ for all $\alpha,\beta\in{\mathbb{N}}_{0}^{d}$ and some $C_{g}<\infty$ as well as $g({\boldsymbol{x}},{\boldsymbol{y}})\geq C_{g}^{-1}|{\boldsymbol{x}}-{\boldsymbol{y}}|$, then, $g^{1/q}$ satisfies (1) with $\widetilde{\varrho}_{1}=1/2$ and $\widetilde{\varrho}_{2}=C_{g}^{-1}$ for all $q\in{\mathbb{N}}$.

• (iii)

If $g$ satisfies the assumptions from (ii) and additionally $g({\boldsymbol{x}},{\boldsymbol{y}})\geq c_{0}>0$ for all ${\boldsymbol{x}},{\boldsymbol{y}}\in D$, then $g^{-1/q}$ satisfies (1) for all $q\in{\mathbb{N}}$.

Proof

To simplify the notation, we consider $g$ as a function of one variable ${\boldsymbol{z}}=({\boldsymbol{x}},{\boldsymbol{y}})\in D\times D\subset{\mathbb{R}}^{2d}$. Define the set of all partitions of $\{1,\ldots,n\}$ as

[TABLE]

For a multi-index $\alpha\in{\mathbb{N}}^{2d}$, we define $\widetilde{\alpha}\in\{1,\ldots,2d\}^{n}$ by $\widetilde{\alpha}_{i}=j$ for all $1+\sum_{k=1}^{j-1}\alpha_{k}\leq i\leq\sum_{k=1}^{j}\alpha_{k}$ and all $1\leq j\leq 2d$ (e.g., $\alpha=(2,3,1,1)$ yields $\widetilde{\alpha}=(1,1,2,2,2,3,4)$). With $n=|\alpha|_{1}$ and some $S\in P\in\Pi(n)$, we define

[TABLE]

(the definition implies $\partial_{{\boldsymbol{z}}}^{\{1,\ldots,n\}}g({\boldsymbol{z}})=\partial_{\boldsymbol{z}}^{\alpha}g({\boldsymbol{z}})$.) With those definitions and given a function $f\colon{\mathbb{R}}\to{\mathbb{R}}$, Faà di Bruno’s formula reads for a multi-index $\alpha\in{\mathbb{N}}^{2d}$

[TABLE]

For (i), Faà di Bruno’s formula (37) and $\partial_{x}^{|P|}\exp=\exp$ show for all multi indices $\alpha\in{\mathbb{N}}^{2d}$ with $n=|\alpha|_{1}$ that

[TABLE]

The definition of asymptotically smooth (1) and $\|g\|_{L^{\infty}(D\times D)}=c_{g}$ imply

[TABLE]

With $f(x):=(1-x)^{-1}$, $x\in{\mathbb{R}}\setminus\{1\}$, we have $\partial_{x}^{k}f(x)=k!(1-x)^{-1-k}$. Hence, the last factor can be written, using Faà di Bruno’s formula again, as

[TABLE]

As the function $h(x):=\exp((1-x)^{-1})$, $x\in{\mathbb{C}}$ is holomorphic at least for $|x|<1$, Cauchy’s integral formula shows

[TABLE]

Altogether, we conclude the proof of (i) by

[TABLE]

For (ii), Faà di Bruno’s formula (37) shows again for $q>1$

[TABLE]

where we used $f(x):=x^{1/q}$ and $|\partial_{x}^{|P|}f(x)|=|(1/q)(1/q-1)(1/q-2)\cdots(1/q-|P|+1)||x|^{1/q-|P|}\leq|P|!|x|^{1/q-|P|}$ as well as the boundedness assumption on the derivatives of $g$ from (ii). With $r(x):=\exp(x)-1$ and $f(x):=(1-x)^{-1}$, $x\in{\mathbb{R}}$, the last factor satisfies

[TABLE]

The function $h(x):=f\circ r(x)=(2-\exp(x))^{-1}$, $x\in{\mathbb{C}}$ is holomorphic at least for $|x|\leq 1/2$. As above, this implies

[TABLE]

and thus concludes the proof of (ii).

For (iii), we conclude the proof as for (ii) by use of the estimate $g(z)^{{-}1/q-|P|}\leq c_{0}^{-1-n}$.

At last, we are ready to prove Lemma 1 which states that the covariance functions from (2) and (3) are asymptotically smooth (1).

Proof (Proof of Lemma 1)

To see (1), consider $\varrho(\cdot,\cdot)$ from (2). We define for complex variables ${\boldsymbol{x}}_{i},{\boldsymbol{y}}_{i}\in{\mathbb{C}}$

[TABLE]

whenever $(\cdot)^{1/p}$ is defined in ${\mathbb{C}}$. and consider $\widetilde{\varrho}({\boldsymbol{x}},{\boldsymbol{y}})$ which is $\varrho({\boldsymbol{x}},{\boldsymbol{y}})$ from (2) but with $d({\boldsymbol{x}}-{\boldsymbol{y}})$ instead of $|{\boldsymbol{x}}-{\boldsymbol{y}}|_{p}$. With the notation of Lemma 10, the above sum has positive real part in O:=\big{\{}({\boldsymbol{x}},{\boldsymbol{y}})\in{\mathbb{C}}^{2d}\,:\,{\boldsymbol{x}}-{\boldsymbol{y}}\notin M\big{\}}. Thus, the function $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto d({\boldsymbol{x}}-{\boldsymbol{y}})$ is holomorphic in each variable in $O$. Since for $a>0$, ${\boldsymbol{x}}\mapsto{\boldsymbol{x}}^{\mu}K_{\mu}(ax)$ is a holomorphic function on ${\mathbb{C}}\setminus({\mathbb{R}}_{-}\cup\{0\})$, and $d({\boldsymbol{x}}-{\boldsymbol{y}})$ has positive real part, we deduce that $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto\widetilde{\varrho}({\boldsymbol{x}},{\boldsymbol{y}})$ is holomorphic in each variable in $O$. Thus, Lemma 9 proves that $\partial_{\boldsymbol{x}}^{\alpha}\partial_{\boldsymbol{y}}^{\beta}\widetilde{\varrho}({\boldsymbol{x}},{\boldsymbol{y}})$ is holomorphic in $O$ in all variables ${\boldsymbol{x}}_{i}$ and ${\boldsymbol{y}}_{i}$. Therefore, Cauchy’s integral formula applied in all variables shows

[TABLE]

The balls $B_{{\boldsymbol{x}},i}$ and $B_{{\boldsymbol{y}},i}$ have to be chosen such that $\prod_{i=1}^{d}B_{{\boldsymbol{x}},i}\times\prod_{i=1}^{d}B_{{\boldsymbol{y}},i}\subset O$. With Lemma 10, and for $({\boldsymbol{x}},{\boldsymbol{y}})\in{\mathbb{R}}^{2d}$ such that ${\boldsymbol{x}}-{\boldsymbol{y}}\in({\mathbb{R}}^{n})_{+}$ (note that Lemma 10 implies $({\boldsymbol{x}},{\boldsymbol{y}})\in O$), this can be achieved by setting $B_{{\boldsymbol{x}},i}:=B_{\varepsilon}({\boldsymbol{x}}_{i})$ and $B_{{\boldsymbol{y}},i}:=B_{\varepsilon}({\boldsymbol{y}}_{i})$ with $\varepsilon:=\sin(\pi/(2p))|{\boldsymbol{x}}-{\boldsymbol{y}}|/(2d+1)$. From this, we obtain the estimate

[TABLE]

for all $({\boldsymbol{x}},{\boldsymbol{y}})\in{\mathbb{R}}^{2d}$ such that ${\boldsymbol{x}}-{\boldsymbol{y}}\in({\mathbb{R}}^{n})_{+}$, where the first equality follows from $d({\boldsymbol{x}}-{\boldsymbol{y}})=|{\boldsymbol{x}}-{\boldsymbol{y}}|_{p}$ for all ${\boldsymbol{x}}-{\boldsymbol{y}}\in({\mathbb{R}}^{n})_{+}$. To remove the restriction ${\boldsymbol{x}}-{\boldsymbol{y}}\in({\mathbb{R}}^{n})_{+}$, consider $b\in\{0,1\}^{d}$ and define the function

[TABLE]

Since we consider $\varrho(\cdot,\cdot)$ from (2), there holds $\varrho\circ F_{b}=\varrho$. Since for all ${\boldsymbol{x}},{\boldsymbol{y}}\in{\mathbb{R}}^{d}$ with ${\boldsymbol{x}}\neq{\boldsymbol{y}}$, there exists some $b\in\{0,1\}^{d}$ such that $({\boldsymbol{x}}_{b},{\boldsymbol{y}}_{b}):=F_{b}({\boldsymbol{x}},{\boldsymbol{y}})$ satisfies ${\boldsymbol{x}}_{b}-{\boldsymbol{y}}_{b}\in({\mathbb{R}}^{n})_{+}$, we prove (38) for all ${\boldsymbol{x}},{\boldsymbol{y}}\in{\mathbb{R}}^{d}$ with ${\boldsymbol{x}}\neq{\boldsymbol{y}}$. Finally, the fact $\alpha!\beta!\leq|\alpha+\beta|_{1}!$, proves that $\varrho(\cdot,\cdot)$ from (2) is asymptotically smooth (1).

Next, consider the covariance function $\varrho(\cdot,\cdot)$ from (3). By definition $\boldsymbol{\Sigma}_{\boldsymbol{x}}$ is continuous on $\overline{D}$. Hence, ${\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}})\geq c_{0}>0$ for all ${\boldsymbol{x}}\in D$. The assumption (4) implies that also ${\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}})$ has bounded derivatives in the sense of (4) (since ${\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}})$ is a polynomial in the matrix entries of $\boldsymbol{\Sigma}_{\boldsymbol{x}}$). Thus, Lemma 12 shows that the functions $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto{\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}})^{1/4}$, $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto{\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{y}})^{1/4}$, and $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto{\rm det}(\boldsymbol{\Sigma}_{\boldsymbol{x}}+\boldsymbol{\Sigma}_{\boldsymbol{y}})^{-q}$, $q\in\{1/2,1\}$ satisfy (1). With $\boldsymbol{\Sigma}_{\boldsymbol{x}}$, also all functions $\widetilde{\boldsymbol{\Sigma}}_{\boldsymbol{x}}$ defined by considering only sub-matrices of $\boldsymbol{\Sigma}_{\boldsymbol{x}}$ satisfy (4). Thus, Cramer’s rule and Lemma 11 show that the map $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto((\boldsymbol{\Sigma}_{\boldsymbol{x}}+\boldsymbol{\Sigma}_{\boldsymbol{y}})^{-1})_{i,j}$ for all $i,j\in\{1,\ldots,d\}$ satisfies (1). From this, we conclude (again with Lemma 11), that $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto({\boldsymbol{x}}-{\boldsymbol{y}})^{T}(\boldsymbol{\Sigma}_{\boldsymbol{x}}+\boldsymbol{\Sigma}_{\boldsymbol{y}})^{-1}({\boldsymbol{x}}-{\boldsymbol{y}})$ as sum and product of asymptotically smooth functions is asymptotically smooth (1). Finally, Lemma 12 shows that $\varrho({\boldsymbol{x}},{\boldsymbol{y}})$ satisfies (1). This concludes the proof.

Appendix B Proof of Proposition 1

The following lemmas state facts about the $H^{2}$-matrix block partitioning, which are well-known but cannot be found explicitly in the literature.

Lemma 13

Under Assumption 1, there exists a constant $C_{B}>0$ which depends only on $d$, $C_{\rm u}$, $D$, and $B_{X_{\rm root}}$ such that all $X\in{\mathbb{T}}_{\rm cl}$ satisfy

[TABLE]

Moreover, all $(X,Y)\in{\mathbb{T}}$ satisfy

[TABLE]

where $C_{BB}>0$ depends only on $C_{B}$, $C_{\rm leaf}$, and $D$.

Proof

The first estimate (39a) follows from the fact that always the longest edge of a bounding box is halved. This means that the ratio $L_{\rm max}/L_{\rm min}$ of the maximal and the minimal side length of a bounding box $B_{X}$ stays bounded in terms of the corresponding ratio for $B_{X_{\rm root}}$. Therefore, we have

[TABLE]

To see the second estimate (39b), consider a given bounding box $B$ with side lengths $L_{1},\ldots,L_{d}$. Due to Assumption 1 the balls $Q_{\boldsymbol{x}}$ with centre ${\boldsymbol{x}}$ and radius $C_{\rm u}^{-1}N^{-1/d}/2$ for all ${\boldsymbol{x}}\in{\mathcal{N}}$ do not overlap. All balls $Q_{\boldsymbol{x}}$ with ${\boldsymbol{x}}\in B$ are containted in a box with sidelengths $L_{\rm max}+C_{\rm u}^{-1}N^{-1/d}$. Thus, the number $m_{B}$ of ${\boldsymbol{x}}\in{\mathcal{N}}$ contained in $B$ can be bounded by

[TABLE]

Since $m_{B}\leq 1$ if $L_{\rm max}<C_{\rm u}^{-1}N^{-1/d}/2$ and since $L_{\rm max}^{d}\simeq|B|$, we may improve the estimate to

[TABLE]

where $C_{B}$ depends only on $d$ and $C_{\rm u}$. On the other hand, Assumption 1 implies that any ball with radius $C_{\rm u}N^{-1/d}$ contains at least one point ${\boldsymbol{x}}\in{\mathcal{N}}$. Since each such ball fits inside a box with sidelength $2C_{\rm u}N^{-1/d}$, we obtain

[TABLE]

points of ${\mathcal{N}}$. This allows us to estimate $m_{B}\geq C_{B}^{-1}|B|N-1$ and conclude (39b). The estimate (39c) follows from the fact $|B_{X}|=|B_{X^{\prime}}|/2$ for all $X\in{\rm sons}(X^{\prime})$. For (40), we observe with (39b) that

[TABLE]

for all $X\in{\mathbb{T}}_{\rm cl}$ with hidden constants depending only on $C_{B}$. Thus, with (39c), we have for all $X\in{\mathbb{T}}_{\rm cl}$ with $X\in{\rm sons}(X^{\prime})$ that

[TABLE]

Moreover, if additionally ${\rm sons}(X)=\emptyset$, we have even $2^{-{\rm level}(X)}\simeq|B_{x}|\lesssim C_{\rm leaf}/N$. By definition of the block-tree ${\mathbb{T}}$, a level difference between $X$ and $Y$ for $(X,Y)\in{\mathbb{T}}$ can only happen, if ${\rm sons}(X)=\emptyset$ or ${\rm sons}(Y)=\emptyset$. Assume ${\rm sons}(X)=\emptyset$. In this case, we have ${\rm level}(Y)\geq{\rm level}(X)$. Then, we have

[TABLE]

with hidden constants depending only on $C_{B}$ and $D$. This implies ${\rm level}(Y)\leq{\rm level}(X)+C$ for some constant $C>0$ which depends only on $C_{\rm leaf}$, $D$, and $C_{B}$ from (39). From this we derive (40) by use of (39).

Lemma 14

Given the definition of ${\mathbb{T}}_{\rm far}$ in Section 3.1, there exists a constant $C>0$ such that all $(X,Y)\in{\mathbb{T}}_{\rm far}$ satisfy

[TABLE]

Proof

By Lemma 13, we have

[TABLE]

For $(X,Y)\in{\rm sons}(X^{\prime},Y^{\prime})$, we obtain additionally

[TABLE]

By definition of the block-partitioning, for $(X,Y)\in{\mathbb{T}}_{\rm far}$ there holds that $B_{X},B_{Y}$ satisfy (5) and $B_{X^{\prime}},B_{Y^{\prime}}$ do not satisfy (5). Altogether, this implies

[TABLE]

where we used $|{\rm level}(X^{\prime})-{\rm level}(X)|\leq 1$. This concludes the proof.

The following lemma gives some basic facts about tensorial Chebychev-interpolation (see, e.g., (h2mat, , Section 4.4))

Lemma 15

Let $f\colon B\to{\mathbb{R}}$ for an axis parallel box $B\subseteq{\mathbb{R}}^{2d}$ such that $\partial_{j}^{k}f\in L^{\infty}(B)$ for all $j=1,\ldots,d$ and all $0\leq k\leq p$. Then, the tensorial Chebychev-interpolation operator of order $p$, $I_{p}\colon C(B)\to{\mathcal{P}}^{p}(B)$ satisfies

[TABLE]

where

[TABLE]

is the operator norm of the one dimensional Chebychev interpolation operator

Proof

It is well-known that the one dimensional Chebychev interpolation operator $I_{p}^{{\boldsymbol{x}}}$ satisfies the error estimate for any $f\in C([-1,1])$

[TABLE]

with an operator norm given in (43). Consider $B:=[-1,1]^{2d}$. Then, there holds with $I_{p}^{{\boldsymbol{x}}_{i}}$ denoting interpolation in the ${\boldsymbol{x}}_{i}$-variable $i\in\{1,\ldots,2d\}$

[TABLE]

Since, for any affine transformation $A\colon{\mathbb{R}}^{2d}\to{\mathbb{R}}^{2d}$, we have $I_{p}(f\circ A)=I_{p}(f)\circ A$, a standard scaling argument concludes the proof.

Proof (Proof of Proposition 1)

We start by proving that $\lambda_{\rm min}(\boldsymbol{C}_{p})>0$ if $p$ satisfies (7). To that end, note

[TABLE]

since the Frobenius norm is an upper bound for the spectral norm. By use of (6) (which is proved below) and (7), we conclude $\lambda_{\rm min}(\boldsymbol{C}_{p})>0$.

To see (6), we first estimate the maximal depth of the tree ${\mathbb{T}}_{\rm cl}$. With (39b)–(39c), we obtain $C_{\rm leaf}\leq|X|\lesssim 2^{-{\rm level}(X)}$ for all $X\in{\mathbb{T}}_{\rm cl}$ with ${\rm sons}(X)\neq\emptyset$. Thus, there holds

[TABLE]

Second, we bound the so-called sparsity constant

[TABLE]

The $H$-matrix case can be found in (hmatorig, , Lemma 4.5). For the $H^{2}$-matrix case, the combination of (40) and (41) (from Lemma 14) shows that $(X,Y)\in{\mathbb{T}}_{\rm far}$ only if $B_{Y}$ touches the (hyper-) annulus with center $B_{X}$ and radii $C^{-1}{\rm diam}(B_{X})$ and $C{\rm diam}(B_{X})$. By comparing the volumes of this annulus and of $B_{Y}$ and using the fact that all the bounding boxes are disjoint, we see that the number of $Y$ such that $(X,Y)\in{\mathbb{T}}_{\rm far}$ is bounded in terms of $C$ and the constants in (39).

For $Y\in{\mathbb{T}}_{\rm cl}$ such that $(X,Y)\in{\mathbb{T}}_{\rm near}$, we have with (39)–(40)

[TABLE]

Again, comparing the volumes of the ball with radius ${\rm diam}(B_{X})$ and of $B_{Y}$, we see that the number of $Y$ such that $(X,Y)\in{\mathbb{T}}_{\rm near}$ is bounded in terms of the constants in (39). Altogether, we bound $C_{\rm sparse}$ uniformly in terms of the constants of Lemma 13. Now, (h2mat, , Lemma 3.38) proves the estimate for storage requirements and (h2mat, , Theorem 3.42) proves the estimate for matrix-vector multiplication.

It remains to prove the error estimate (see also (h2mat, , Section 4.6) for the integral operator case). To that end, note that since the near field ${\mathbb{T}}_{\rm near}$ is stored exactly, there holds

[TABLE]

Given, $(i,j)\in I(X)\times I(Y)$, we have with the interpolation operator $I_{p}$ from Lemma 15 and (1)

[TABLE]

With the admissibility condition (5), we get

[TABLE]

and hence

[TABLE]

The combination of the above estimates concludes the proof.