Exploration of a Cosine Expansion Lattice Scheme

TL;DR

This paper introduces a novel approximation scheme combining Quasi-Monte Carlo lattice sequences with Fourier-cosine methods for expectation computation, analyzing its error and comparing it with wavelet-based approaches through numerical experiments.

Contribution

It presents a new lattice-based scheme that integrates Fourier-cosine techniques, advancing numerical expectation computation methods.

Findings

The scheme effectively approximates expectations with quantifiable error.

Numerical experiments demonstrate the scheme's competitiveness against wavelet methods.

Abstract

In this article, we combine a lattice sequence from Quasi-Monte Carlo rules with the philosophy of the Fourier-cosine method to design an approximation scheme for expectation computation. We study the error of this scheme and compare this scheme with our previous work on wavelets. Also, some numerical experiments are performed.

Tables

Table 1. Table 1: Reference solution for Uniform Distribution

Dimension	1	2	3	4	5	6	7	8
Ref. Values	1	2.167	2.167	4.424	4.424	8.893	8.893	17.81

Equations

$$\mathbb{E}[\mathfrak{f}(\mathbf{Y})]={\int }_{{\mathbb{R}}^s}\mathfrak{f}(\mathbf{y})\mu (d\mathbf{y}).$$

$${\int }_{[0,1{]}^s}\mathfrak{f}(\mathbf{y})d\mathbf{y}$$

$$\frac{1}{N}\sum _{n=0}^{N-1}\mathfrak{f}({\mathbf{y}}_n).$$

$${\left(\left\{\frac{n\mathbf{g}}{N}\right\}\right)}_{0\le n<N-1},$$

$${\bf 1}_{\mathcal{D}}({\bf y})=\left\{\begin{array}[]{ll}1,&\mbox{ if }{\bf y}\in\mathcal{D};\\ 0,&\mbox{ otherwise;}\end{array}\right.$$

$$\bar{\mathfrak{f}}(y):=\sum _{k=1}^{N-1}{}^{\prime }{\tilde{\mathfrak{f}}}_{\cos }(k)\cos \left(k\pi \frac{y-a}{b-a}\right),$$

$${\tilde{\mathfrak{f}}}_{\cos }(k):=\frac{2}{b-a}{\int }_a^b\mathfrak{f}(y)\cos \left(k\pi \frac{y-a}{b-a}\right)dx.$$

$$\mathbb{E}[\mathfrak{f}(Y)]$$

$$=\sum _{k=1}^{N-1}{}^{\prime }{\tilde{\mathfrak{f}}}_{\cos }(k){\int }_{\mathbb{R}}\cos \left(k\pi \frac{y-a}{b-a}\right)\mu (dy)$$

$$=\sum _{k=1}^{N-1}{}^{\prime }{\tilde{\mathfrak{f}}}_{\cos }(k)\Re \left\{\mathcal{F}\left(\frac{k\pi }{b-a}\right)\exp \left(-ık\pi \frac{a}{b-a}\right)\right\}.$$

$$<\mathfrak{f},\mathfrak{g}{>}_{{\mathcal{K}}_{\alpha ,\mathbf{\gamma },s}^{\cos }(\mathcal{D})}=\sum _{\mathbf{k}\in {\mathbb{N}}^s}{\tilde{\mathfrak{f}}}_{\cos }(\mathbf{k}){\tilde{\mathfrak{g}}}_{\cos }(\mathbf{k})r_{\alpha ,\mathbf{\gamma },s}^{-1}(\mathbf{k}),$$

$${\int }_{[0,1{]}^s}\mathfrak{f}(\mathbf{y}\times (\mathbf{b}-\mathbf{a})+\mathbf{a})2^{|\mathbf{k}{|}_0/2}\prod _{j=1}^s\cos (\pi k_j{y}_j)d\mathbf{y},\text{mbox}for\mathbf{k}=(k_1,\cdots ,k_s)\in {\mathbb{N}}^s.$$

$$r_{\alpha,\gamma}(k):=\left\{\begin{array}[]{ll}1&\mbox{if }k=0;\\ \gamma|k|^{-2\alpha}&\mbox{if }k\neq 0,\end{array}\right.$$

$$r_{\alpha ,\mathbf{\gamma },s}(\mathbf{k}):=\prod _{j=1}^s{r}_{\alpha ,{\gamma }_j}(k_j),$$

$$||\mathfrak{f}|{|}_{{\mathcal{K}}_{\alpha ,\mathbf{\gamma },s}^{\cos }(\mathcal{D})}^2=\sum _{\mathbf{k}\in {\mathbb{N}}^s}\frac{|{\tilde{\mathfrak{f}}}_{\cos }(\mathbf{k}){|}^2}{r_{\alpha ,\gamma ,s}(\mathbf{k})}=\sum _{\mathbf{h}\in {\mathbb{Z}}^s}{2}^{-|\mathbf{h}{|}_0}\frac{|{\tilde{\mathfrak{f}}}_{\cos }(\mathbf{h}){|}^2}{r_{\alpha ,\gamma ,s}(\mathbf{h})}.$$

$${\mathcal{K}}_{\alpha ,\mathbf{\gamma },s}^{\cos }(\mathcal{D},\mathbf{x},\mathbf{y}):=$$

$$=$$

$${\mathcal{K}}_{\alpha ,\gamma }^{\cos }(\mathcal{D},x,y):=$$

$$=$$

$$\mathfrak{f}(\mathbf{y})=<\mathfrak{f},{\mathcal{K}}_{\mathcal{D},\alpha ,\mathbf{\gamma },s}^{\cos }(\cdot ,y){>}_{{\mathcal{K}}_{\alpha ,\mathbf{\gamma },s}^{\cos }(\mathcal{D})}.$$

$${\mathcal{K}}_s^{\cos }(\mathbf{x},\mathbf{y}):=$$

$$\mathfrak{f}(\mathbf{y})=\sum _{\mathbf{k}\in {\mathbb{N}}^s}{\tilde{\mathfrak{f}}}_{\cos }(\mathbf{k})(\sqrt{2}{)}^{|\mathbf{k}{|}_0}\prod _{j=1}^s\cos \left(\pi k_i\frac{y_i-a_i}{b_i-a_i}\right),$$

$$\bar{\mathfrak{f}}(\mathbf{y}):=\sum _{\mathbf{k}\in {\mathbb{N}}^s}{\tilde{\mathfrak{f}}}_{\cos }(\mathbf{k})(\sqrt{2}{)}^{|\mathbf{k}{|}_0}\prod _{j=1}^s\cos \left(\pi k_i\frac{y_i-a_i}{b_i-a_i}\right),$$

$${\tilde{\mathfrak{f}}}_{\cos }(\mathbf{k}):={\int }_{[0,1{]}^s}\mathfrak{f}(\mathbf{y}\times (\mathbf{b}-\mathbf{a})+\mathbf{a})2^{|\mathbf{k}{|}_0/2}\prod _{j=1}^s\cos (\pi k_j{y}_j)d\mathbf{y}.$$

$$\mathbb{E}[\mathfrak{f}(Y)]\approx \frac{1}{N}\sum _{n=1-N}^N\mathfrak{f}\left(\frac{r}{2^{\mathtt{m}}}\right)\mathbb{E}[{\phi }_{N,n}(Y)],$$

$${\int }_{{\mathbb{R}}^s}\bar{\mathfrak{f}}(\mathbf{y})\mu (d\mathbf{y}).$$

$${\int }_{{\mathbb{R}}^s}\frac{1}{N}\sum _{n=0}^{N-1}\mathfrak{f}\left({\mathbf{p}}_n\right){\mathcal{K}}_s^{\cos }\left({\mathbf{p}}_n,\mathbf{y}\right)\mu (d\mathbf{y}),$$

$$P(\mathbf{g},N):=\left\{\left\{\frac{n\mathbf{g}}{N}\right\}:0\le n<N\right\}.$$

$$\varphi (x):=1-|2x-1|,$$

$$P_{\mathbf{\varphi }}(\mathbf{g},N):=\left\{\mathbf{\varphi }\left(\left\{\frac{n\mathbf{g}}{N}\right\}\right):0\le n<N\right\}.$$

$$P_{\mathbf{\varphi }}(\mathbf{g},N,\mathcal{D})=$$

$${\left({\int }_{{\mathbb{R}}^s}\left(\prod _{j=1}^s\cos \left(\pi k_j\frac{y_j-a_j}{b_j-a_j}\right)\right)\mu (d\mathbf{y})\right)}^2\le r_{\beta ,\mathbf{\rho },s}(\mathbf{k}),$$

Full text

1 Introduction

In this work, we derive a numerical integration formula for a function $\mathfrak{f}:\mathbb{R}^{s}\rightarrow\mathbb{R}$ with respect to a probability measure $\mu$. Namely, we aim to approximate the following quantity:

[TABLE]

The inspiration of this work comes from two promising methods in numerical integration, the COS method and lattice rules.

The COS method is a numerical integration method developed in the context of option pricing, see [10]. Within this framework, the fair value of a financial option can be expressed as the expectation of the discounted payoff function. Using the fact that the payoff function can be expressed as a cosine series, an approximation of the integration in terms of the Fourier transform of the risk neutral measure and the cosine coefficient of the payoff function were presented. This method is efficient in terms of calculation with only simple addition operations. It also makes use of the characteristic function of a probability measure, which is available more often than the density function. Ever since, there has been numerous extensions of the COS method, including pricing Bermudan options [11], finding ruin probability [4], solving backward stochastic differential equations [22] and computation of valuation adjustment [2].

There are also several applications for which it is necessary to extend the COS method beyond the one-dimensional situation, for example, in [21]. However, the tensor extension used in previous studies suffers from the curse of dimensionality, with the number of summation terms increasing exponentially when the number of dimensions increases. Therefore, further input is required for such extensions to be feasible in practice.

Thus, we took the insight from Quasi-Monte Carlo (QMC) rules, which approximate integrals of the form

[TABLE]

with equal weight quadrature rules

[TABLE]

Readers are referred to [8], [18] and references therein for further details.

In particular, we are interested in the rank 1 lattice rule construction method for quadrature points, where for a given positive number $N$ under dimension $s$, one chooses a vector ${\bf g}\in\{1,\ldots,N-1\}^{s}$ called the generating vector and generates a points set:

[TABLE]

where the notation $\{\}$ denotes the fractional part of the real number in each dimension. The quadrature points are the result of applying some fixed function on this points set.

There have been numerous research papers supporting the implementation of lattice rule QMC, in identifying a suitable generating vector [6, 16, 24, 23], in efficient algorithms to generate such vectors [19, 20], and in extensible lattice rules [5, 9, 13, 14]. In [7], the authors derived an error bound for QMC with tent transformed lattice rules for functions within the half-period cosine space. Noting the connection between the half-period cosine space and the tensor extension of the cosine series, we are inspired to combine the two approaches and transfer the rich results from the lattice literature to Fourier expansion schemes.

This work is organized as follows. In Section 2, we present the two components of the cosine expansion lattice scheme, the half-period cosine expansion and the tent transformed lattice. We also present the full scheme in this section. In Section 3, we give further details regarding our current results by providing an alternative formation of its error bound and connecting the scheme to periodic wavelets. Numerical experiments are presented in Section 4 and we conclude our findings in Section 5.

Before we begin, we mention some conventions used in this work. We assume all the integrals in the computation to be finite, therefore Fubini’s theorem can be applied and we exchange the order of integration without notice. The operations $\times$ and $/$ act component-wise when used on vectors.

Finally, some notation we use in this article:

• •

The natural number set $\mathbb{N}:=\{0,1,2,\ldots\}$;

• •

The positive integer set $\mathbb{Z}^{+}:=\{i\in\mathbb{Z}|i>0\}$, similarly for $\mathbb{R}^{+}$;

• •

The truncated integer set, for $s\in\mathbb{Z}^{+}$ we write $[s]:=\{1,\ldots,s\}$;

• •

The indicator function ${\bf 1}_{\mathcal{D}}:\mathbb{R}^{s}\rightarrow\{0,1\}$

[TABLE]

• •

The ceil function $\lceil\cdot\rceil:\mathbb{R}\rightarrow\mathbb{Z}$, $\lceil y\rceil=\min\{i\in\mathbb{Z}|i\geq y\}$.

2 Cosine Expansion Lattice Scheme

In this section, we introduce the cosine expansion lattice scheme, whose construction consists of two parts: a projection of the original integrand on a reproducing kernel space and a numerical integration technique based on a tent-transformed lattice rule. We will briefly describe the intuition behind our derivation.

2.1 The half-period cosine space

In a similar framework as the COS method from [10], we would like to define a cosine based periodic expansion of function $\mathfrak{f}$ in $\mathbb{R}^{s}$. This allows us to connect the expectation problem (1.1) to a Fourier transform.

The common practice to transform Equation (1.1) into a finite problem for computational purposes is restricting the domain of integration to a predefined box $\mathcal{D}:=[a_{1},b_{1}]\times\cdots\times[a_{s},b_{s}]$. This step is justified as long as $\mathcal{D}$ contains the majority of the mass of the measure. In fact, this is equivalent to replacing the original integrand $\mathfrak{f}({\bf y})$ by $\mathfrak{f}({\bf y}){\bf 1}_{\mathcal{D}}({\bf y})$, or setting the function values outside $\mathcal{D}$ to zero.

The COS method uses this concept to its advantage by replacing an originally non-periodic one-dimensional function $\mathfrak{f}$ by a cosine series projection that coincides with $\mathfrak{f}$ on the domain $[a,b]$, but is periodic throughout the whole real line.

To be precise, the integrand $\mathfrak{f}$ is replaced by111The notation $\sum{}^{\prime}$ denotes the first term of the summation is weighted by one half.

[TABLE]

where

[TABLE]

In this case, the original 1-D expectation is approximated by

[TABLE]

Note that within $[a,b]$, $\bar{\mathfrak{f}}$ is a projection of $\mathfrak{f}$ onto a finite cosine series space and $\bar{\mathfrak{f}}$ converges to $\mathfrak{f}$ in the $L^{2}([a,b])$ norm, when $N$ tends to infinity. However, when considering the function $\bar{\mathfrak{f}}$ on the whole domain $\mathbb{R}$, it is a periodic function and it thus deviates from $\mathfrak{f}$ itself. Once again, the error of this transformation is kept under control as the probability mass outside $[a,b]$ is small and we make use of the Fourier transform $\mathcal{F}$, which is typically available, in the approximation.

The main goal of this work is to apply the same idea in a higher-dimensional setting. Thus, we will consider a functional space built on cosine functions and use the series projection of $\mathfrak{f}$ onto such space as our replacement integrand. The space we pick is a modification to the half-period cosine space introduced in Section 2.3 of [7] and their article serves as an inspiration of the current work.

We define an inner product $<\cdot,\cdot>_{\mathcal{K}^{\cos}_{\alpha,{\boldsymbol{\gamma}},s}(\mathcal{D})}$ as

[TABLE]

for some real number $\alpha>1/2$ and vector $\boldsymbol{\gamma}\in(\mathbb{R}^{+})^{s}$, where the multi-dimensional cosine coefficients are given by

[TABLE]

Here we define $|{\bf k}|_{0}:=\#\{j\in[s]:k_{j}\neq 0\}$ to be the number of non-zero components in ${\bf k}$. Note that only the portion of $\mathfrak{f}$ within the predefined domain $\mathcal{D}$ is used here.

For $\alpha>1/2,k\in\mathbb{Z}$ and $\gamma>0$, the one-dimensional $r$ function is defined as:

[TABLE]

and the multidimensional $r$ function is set as,

[TABLE]

for ${\bf k}=(k_{1},\ldots,k_{s})\in\mathbb{Z}^{s}$ and ${\boldsymbol{\gamma}}=(\gamma_{1},\ldots,\gamma_{s})\in(\mathbb{R}^{+})^{s}$. The $r$ function is introduced to the norm here to assess the decay rate of the cosine coefficients, as it is closely related to the approximation error.

We define the corresponding norm as $\sqrt{<\mathfrak{f},\mathfrak{f}>_{\mathcal{K}^{\cos}_{\alpha,{\boldsymbol{\gamma}},s}(\mathcal{D})}}$ and denote it by $||\mathfrak{f}||_{\mathcal{K}^{\cos}_{\alpha,{\boldsymbol{\gamma}},s}(\mathcal{D})}$. In particular, we have

[TABLE]

The dummy variable is changed from ${\bf k}$ to ${\bf h}$ here in preparation for future computations. We denote any function $\mathfrak{f}$ such that $||\mathfrak{f}||_{\mathcal{K}^{\cos}_{\alpha,{\boldsymbol{\gamma}},s}(\mathcal{D})}<\infty$ as $\mathfrak{f}\in\mathcal{H}_{\mathcal{K}^{\cos}_{\alpha,{\boldsymbol{\gamma}},s}(\mathcal{D})}$.

The half-period cosine space is an example of a reproducing kernel Hilbert space with the corresponding reproducing kernel,

[TABLE]

and for the one-dimensional version,

[TABLE]

For any ${\bf y}\in\mathcal{D}$ and $\mathfrak{f}\in\mathcal{H}_{\mathcal{K}^{\cos}_{\alpha,\boldsymbol{\gamma},s}(\mathcal{D})}$, we have the reproducing property,

[TABLE]

Readers are referred to [1] and [8] for further information on reproducing kernel Hilbert spaces.

In this work, we use an alternative kernel which drops the $r$ function. It is defined as

[TABLE]

We suppress the $\mathcal{D}$ part here to simplify our notation.

Using the reproducing property, we have the following equation:

[TABLE]

for any $\mathfrak{f}\in\mathcal{H}_{\mathcal{K}^{\cos}_{\alpha,{\boldsymbol{\gamma}},s}(\mathcal{D})}$ and ${\bf y}\in\mathcal{D}$. We use the expansion at the right-hand side of (2.1) as the replacement integrand in Equation (1.1), denoted by $\bar{\mathfrak{f}}$.

Definition 2.1** (Half-period Cosine Expansion).**

For any given function $\mathfrak{f}:\mathbb{R}^{s}\rightarrow\mathbb{R}$ and given domain $\mathcal{D}\subset\mathbb{R}^{s}$, the half-period cosine expansion, $\bar{\mathfrak{f}}:\mathbb{R}^{s}\rightarrow\mathbb{R}$, is defined as

[TABLE]

where the cosine coefficients are defined as

[TABLE]

2.2 Lattice rule approximations

There are essentially two drawbacks when extending the COS method to higher dimensions. First, the cosine coefficient $\widetilde{\mathfrak{f}}_{\cos}$ may be difficult to calculate, especially in a recurring situation. In our previous work [3], we aimed to remedy this shortcoming by adopting a wavelet basis such that we can approximate Equation (1.1) as a weighted sum of local values, in the form of

[TABLE]

with $\varphi$ being the wavelet basis functions. We will compare our current scheme to that work in Section 3.2.

The second point of concern is that if we simply apply the tensor product of cosine spaces in a higher-dimensional case, as we have done in the last section, we will face the curse of dimensionality.

Here, we aim to address the above two issues by constructing quadrature rules for the integration of

[TABLE]

In particular, the approximant takes the specific form,

[TABLE]

where $\mathit{P}:=\{{\bf p}_{0},\ldots,{\bf p}_{N-1}\}$ is some predetermined quadrature points set.

For the half-period cosine space that we adopted in the previous section, Dick et. al. showed in [7] that a combination of rank-1 lattice rules and tent transformations converges well in the context of quasi-Monte Carlo methods. We will now introduce the quadrature points proposed in [7] for the half-period cosine space, which we should adapt for the approximation of Equation (2.3).

A lattice point set with $N\geq 2$ points and generating vector ${\bf g}\in[N-1]^{s}$ is given by

[TABLE]

The tent transformation, $\phi:[0,1]\rightarrow[0,1]$ is defined as

[TABLE]

and the higher-dimensional version $\boldsymbol{\phi}:[0,1]^{s}\rightarrow[0,1]^{s}$ is obtained by applying the function component-wise. The tent-transformed lattice point set in [7] is given by

[TABLE]

However, since we have to apply the quadrature points on a general box $\mathcal{D}$, instead of just on $[0,1]^{s}$, we have to transform the lattice rules:

[TABLE]

We wish to quantify the integration error of applying Equation (2.3) to approximate the expectation in the form of equation (2.2) for functions in the half-period cosine space $\mathcal{H}_{\mathcal{K}^{\cos}_{\alpha,\boldsymbol{\gamma},s}(\mathcal{D})}$. In addition to the above condition, we also enforce some convergence requirements on the cosine transform of the measure $\mu$.

Theorem 2.2**.**

Consider any function $\mathfrak{f}\in\mathcal{H}_{\mathcal{K}^{\cos}_{\alpha,\boldsymbol{\gamma},s}(\mathcal{D})}$ and any probability measure $\mu$ such that

[TABLE]

for some real number $\beta>s$ and vector $\boldsymbol{\rho}\in(\mathbb{R}^{+})^{s}$. Furthermore, we assume that $\beta-\alpha>1/2$, namely, the cosine transform of the measure $\mu$ decays at least algebraically and quicker than the cosine coefficients of $\mathfrak{f}$. The error for the numerical integration using the tent-transformed lattice rule on $\mathcal{D}$ is then bounded by

[TABLE]

for some constant $C_{i}$, where $L^{\bot}:=\{{\bf h}\in\mathbb{Z}^{s}:{\bf h}\cdot{\bf g}\equiv 0\;(\bmod\;N)\}$ is the dual lattice.

Proof.

The general flow of this proof follows the proof of Theorem 2 in [7]. Let $\mathfrak{f}\in\mathcal{H}_{\mathcal{K}^{\cos}_{\alpha,\boldsymbol{\gamma},s}(\mathcal{D})}$ and $\bar{\mathfrak{f}}$ be its half-period cosine expansion.

For any $k\in\mathbb{N}$, we have

[TABLE]

and hence

[TABLE]

The reproducing kernel can be rewritten as

[TABLE]

Therefore, we obtain

[TABLE]

The sum $\frac{1}{N}\sum^{N-1}_{n=0}e^{2\pi\imath n[({\bf h}+{\bf k})\cdot{\bf g}]/N}$ is a character sum over the group $\mathbb{Z}/N\mathbb{Z}$, which is equal to one if ${\bf(h+k)}\cdot{\bf g}$ is a multiple of $N$ and zero otherwise. From this, we get

[TABLE]

Based on this formula and an application of the Cauchy-Schwarz inequality, we obtain

[TABLE]

Making use of the smoothness assumption on the probability measure, we find

[TABLE]

In order to control the error, we need to control the sum $I_{j}$ in Equation (2.6) for all $h_{j}\in\mathbb{Z}$. We have three different cases to consider.

Case 1: $h_{j}=0.$

In this case,

[TABLE]

in which $\zeta$ denotes the Riemann zeta function.

Case 2: $h_{j}>0.$

In this case,

[TABLE]

The above inequality simply follows from the definition and the orderings $|h_{j}+k_{j}|^{-1}<|h_{j}|^{-1}$, $|h_{j}|^{-\beta}<|h_{j}|^{-\alpha}$ and $(h_{j}-1)\times 1<(h_{j}-2)\times 2<\cdots\left(\frac{h_{j}}{2}\right)$. Note that we use the convention $\sum^{0}_{k_{j}=1}=0$ here, so the inequality also holds for $h_{j}=1$.

Case 3: $h_{j}<0.$

Finally, in this case,

[TABLE]

The derivation of Case 3 is similar to Case 2. The only difference is that we have the following orderings: $|-h_{j}+k_{j}|>|h_{j}|$ for $k_{j}>0$ and $(-h_{j}-1)\times 1<(-h_{j}-2)\times 2<\cdots<\left(\frac{h_{j}}{2}\right)^{2}$.

Finally, let $C_{j}=\max\{1+2\gamma_{j}\rho_{j}\zeta(2\alpha+2\beta),1+\rho_{j}\left(\zeta(2\alpha)+\zeta(2\beta)+\frac{1}{\gamma_{j}}+2^{2\alpha}\zeta(2(\beta-\alpha))\right)\}$, which is independent of $h_{j}$. Then,

[TABLE]

In this proof, we break down the sum in Equation (2.5) and derive a bound for each dual lattice point ${\bf h}$ along each direction, namely, the $I_{j}(h_{j})$ terms in (2.6). By considering the three possible cases for $h_{d}$, positive, negative and zero, we provide a bound in each case. It is necessary to separate the three cases as we have to identify the section where both $|h_{j}-k_{j}|$ and $|k_{j}|$ are smaller than $|h_{j}|$ when $k_{j}$ is moving along the real line and apply a separated bound on these sections.

Taking the maximum out of the three cases and putting $I_{j}(h_{j})$ back into the sum in Equation (2.3) finishes the proof. ∎

Remark 2.3*.*

The integrals of the form, $\int_{[0,1]^{s}}\mathfrak{f}({\bf y})d{\bf y},$ can be seen as special cases of Equation (2.2) where the probability measure is an uniform distribution over $\mathcal{D}=[0,1]^{s}$. In this case,

[TABLE]

Therefore, we may simplify the square of the term in Equation (2.5), as follows

[TABLE]

which is the same as the right-hand side of Equation (6) in [7]. Therefore, our numerical integration can be seen as an extension of the Quasi-Monte Carlo (QMC) rules in [7].

The above proof is rather rough, note the $2^{2\alpha}$ term in the last two cases. However, it suggests we can still relate the integration error to the smoothness of the target function $\mathfrak{f}$ in terms of the cosine coefficients. The error can be controlled by the sum $\left(\sum_{{\bf h}\in L^{\bot}\backslash\{{\bf 0}\}}r_{\alpha,\boldsymbol{\gamma},s}({\bf h})\right)^{1/2}$. Thus all the results for the worst-case error for the QMC integration in the half-period cosine space using tent-transformed lattice rules can be used here. This opens the door for applying the generating vector for extensible lattice rules from [13] and other results from the QMC literature to the half-period cosine expansion scheme. We believe this result is a promising starting point for the study of lattice expansions for higher-dimension numerical integration with general finite measure.

2.3 Full Approximation Schemes and Errors

Now we can introduce the full approximation scheme. The final obstacle lies in the reproducing kernel, which is defined as an infinite sum. Instead of calculating the full sum, a truncated version is used in the actual algorithm. We define

[TABLE]

For any expectation that satisfies the conditions in the previous section, we have our full approximant:

[TABLE]

where $\mathcal{F}_{\mu}$ is the Fourier transform of the measure $\mu$. Recall from Equation (2.1) that $\mathfrak{f}$ and $\bar{\mathfrak{f}}$ coincide in $\mathcal{D}$. We can simply replace $\bar{\mathfrak{f}}$ by $\mathfrak{f}$ in the above scheme.

In practice, one would first calculate the expected reproducing kernel value $\int_{\mathbb{R}^{s}}\mathcal{K}^{\cos}_{s,K}\left({\bf p}_{\boldsymbol{\phi},n},{\bf y}\right)\mu(d{\bf y})$ for each lattice point, possibly in an offline setting and then calculate the main sum.

Define the absolute approximation error under this setting as $\epsilon$ and we see that

[TABLE]

For the third term, a bit more detail is required. Again we have the convergence assumption on the cosine transform of measure $\mu$ with the additional condition that $\beta>s$. With the inequality $(K-1)\times 1<(K-2)\times 2<\ldots<\left(\frac{K}{2}\right)^{2}$, for all $K$ with $K\in\mathbb{N}$ and $K>2$ and mathematical induction, one can show that if $\sum|k_{i}|=K$,

[TABLE]

There are in total $\frac{1}{(s-1)!}\prod^{s-1}_{j=1}(K+j)$ s-tuples of natural numbers satisfying $\sum^{s}_{i=1}k_{i}=K$ for $s\geq 2$. (see [17] for further explanation.) Taking into account all the possible combinations of positive and negative signs among all dimensions when we consider all integers, we have

[TABLE]

The final inequality also holds when $s=1$.

Therefore, we find the following error bound for the lattice expansion scheme:

[TABLE]

To conclude, we can apply the telescoping technique and describe our error by three separate pieces: the projection error from a non-periodic to a periodic function, the lattice integration error and the kernel truncation error. With this proof, we successfully extended the lattice integration in [7] to a more general setting and made use of both the smoothness of the integrand in terms of cosine coefficients and the convergence of the cosine transform with respect to the probability measure. However, with our “number theory based” derivation, the model dimension $s$ still heavily impacts the complete error bound, which may be an obstacle when applying the method to higher-dimensional problems. Moreover, the assumption of the cosine transform’s decay is rough. Results from Fourier analysis may be incorporated to further improve the error bound. Further study for the error control is therefore recommended.

3 Discussion

In Section 2, we have provided a justification for the cosine expansion lattice scheme. We extended the tent transformed lattice rule in the half-cosine space to general probability measures, with its error controlled by the decay rate of the cosine transform of the probability measures as well as the cosine coefficients of the integrands.

However, there are remaining questions. Is there a better way to control the approximation error than simple algebraic manipulation? How does our current scheme connect to previous work on wavelets? Here we aim to present some discussion on the above issues.

3.1 Alternative Error Formulation

The first possibility we would like to consider is whether the error of the approximation in Section 2.2 can be written in an alternative form. By doing so, we aim to avoid bounding the term $I({\bf h})$ along each dimension, as this estimation is rough.

We revisit the derivation of $\epsilon_{2}$ in Section 2.2, specifically the right-hand side of Equation (2.4). Note that the innermost sum in the expression can be rewritten as follows:

[TABLE]

by using the definition of cosine coefficients. We separate the sum here in three parts according to whenever $k_{1}$ is positive, negative or zero. Next, we combine the terms whose $|k_{1}|$ value is the same.

[TABLE]

Then we repeat the similar steps of separating terms and combining those with the same absolute value for $k_{2},\cdots,k_{s}$, which leads to:

[TABLE]

This can be seen as the half-period cosine expansion of the function

[TABLE]

Therefore, the right-hand side of Equation (2.4) can be rewritten as

[TABLE]

To summarize, the error of this approximation is just the sum of functions $\bar{\mathfrak{f}}_{{\bf h}}$ integrated over the probability measure $\mu$ over the dual lattice. The approximation error is controlled by a series with individual term $\int_{\mathbb{R}^{s}}\bar{\mathfrak{f}}_{{\bf h}}({\bf y})\mu(d{\bf y})$, a cosine transform with respect to the cosine series of $\mathfrak{f}$ under the probability measure $\mu$. The decay rate of such integral, depending on the combined smoothness of the measure $\mu$ and the function $\mathfrak{f}$, is the key to limit the error of our scheme.

3.2 Cosine Wavelets

We place our work in the context of the construction of wavelets in [3] to see the similarities and differences between the two approaches. Instead of using the full Fourier series for the wavelet construction, as in [3], we choose a “cosine function only” basis set as our starting point, corresponding to the half-period cosine space we used in Section 2.

Under this setting, we start our derivation from the set

[TABLE]

with $[a,b]\in\mathbb{R}$ a finite range. We define an inner product

[TABLE]

and the corresponding norm $||\cdot||_{L^{2}([a,b])}$. We denote any function $\mathfrak{f}$ such that $||\mathfrak{f}||_{L^{2}([a,b])}<\infty$, as $\mathfrak{f}\in L^{2}([a,b])$.

Equipped with the above definitions, we construct an approximation space together with a localized basis. Consider the following function $\mathcal{K}^{wl}_{N^{\prime}}:\mathbb{R}\times[N^{\prime}]\rightarrow\mathbb{R}$,

[TABLE]

where $r=1,2,\ldots,N^{\prime}$. This definition is a special case of the scaling functions given in Equation (2.13) of [12], in which the authors presented an uniform approach for the construction of wavelets based on orthogonal polynomials. The properties of $\mathcal{K}^{wl}_{N^{\prime}}$, that are relevant to our numerical method are listed in the next proposition.

Proposition 3.1**.**

The function $\mathcal{K}^{wl}_{N^{\prime}}$, which is defined in Equation (3.1), satisfies the following properties:

1. (a)

The inner product of two scaling functions is given by the following equation:

[TABLE]

Thus, $\{\mathcal{K}^{wl}_{N^{\prime}}(x,r)|r=1,\ldots,N^{\prime}\}$ is an orthogonal set. 2. (b)

The scaling function $\mathcal{K}^{wl}_{N^{\prime}}(\cdot,r)$ is localized around $a+\frac{2r-1}{2N^{\prime}}(b-a)$. By this we mean that for the subspace

[TABLE]

we have

[TABLE] 3. (c)

$\{\mathcal{K}^{wl}_{N^{\prime}}(\cdot,r)|r=1,2,\ldots,N^{\prime}\}$* is a basis for $V_{N^{\prime}}$.* 4. (d)

The scaling function $\mathcal{K}^{wl}_{N^{\prime}}$ is also a kernel polynomial in the sense that, for any function $v$ in $V_{N^{\prime}}$, we have

[TABLE]

Readers are referred to [12] for further properties of such functions and for a condensed proof of the above properties one may follow the derivation of Theorem 3.1 in [3].

Applying the results obtained above, we may define the wavelet expansion $\mathfrak{f}_{wl}:\mathbb{R}\rightarrow\mathbb{R}$ for any function $\mathfrak{f}\in L^{2}([a,b])$ by

[TABLE]

and approximate the expectation $\mathbb{E}[\mathfrak{f}(Y)]$ by,

[TABLE]

Comparing this expression with Equation (2.8) when $s=1$, it is clear that these two formulas are of the same form. The main difference is that the wavelet formula can be seen as an expansion scheme using the lattice points

[TABLE]

with the kernel also bounded at the corresponding value $N^{\prime}-1$.

Considering a two-dimensional function $\mathfrak{f}:\mathbb{R}^{2}\rightarrow\mathbb{R}$, we can extend the cosine wavelet to two dimensions by applying the expansion to each dimension separately.

[TABLE]

Again this is in a similar form as in Equation (2.8). However, the total number of terms ${N^{\prime}}^{2}$ increases with the dimension and there is no clear way to select the more significant ones.

One of the key advantages of the cosine expansion lattice scheme is that we have removed the link between the quadrature points and the number of summation terms $N$. This allows us to apply an online/offline construction in our algorithm. We can first compute the expectations of reproducing kernels $\mathbb{E}\left[\mathcal{K}^{\cos}_{s,K}\left({\bf p}_{\boldsymbol{\phi},n},{\bf Y}\right)\right]$ for all points in the lattice sequence and use the stored values in the approximation scheme. This construction is not feasible in the cosine wavelets setting and it is one of the advancements from the cosine expansion lattice scheme.

4 Numerical Experiments

In this section we perform numerical experiments to test the cosine expansion lattice scheme.

In order to demonstrate that our scheme can inherit the results from the previous literature, we use the lattice sequence from [13] and perform the tent-transformation on them in all our tests. This sequence has been used in the numerical results section of [7] and readers are referred to the references therein for further information.

All figures displaying numerical results in this Section present the log absolute error $\log_{10}(|\epsilon|)$ against the log number of lattice points $\log_{10}(N)$.

For the first two experiments, we compute the expectation with the following test function:

[TABLE]

4.1 Uniform Distribution

Whenever the reproducing kernel’s expectation $\mathbb{E}[\mathcal{K}^{\cos}_{s}\left({\bf p}_{\boldsymbol{\phi},n},{\bf Y}\right)]$ is known explicitly, we can simplify the algorithm to its original form in Equation (2.3) instead of the full scheme in Equation (2.8). This results in an algorithm that has no difference to the original QMC rule in terms of computational complexity. As stated in Remark 2.3, this is indeed the case for the uniform distribution.

In this subsection, we approximate

[TABLE]

where ${\bf Y}$ follows an uniform distribution on the domain $[0,1]\times[-1,1]\times[0,1]\times[-1,1]\times\cdots$ for dimension $s\in[8]$. The reference value is $\prod_{j=1}^{s}\left(2+w^{j}\frac{2}{3}\right)^{(j\mod 2)}$.

Note that we are not establishing a new numerical result or proof of convergence here as this problem can be solved by the original QMC algorithm through a change of variables. These tests serve two purposes. They establish the base line result as a comparison to the results in the next section and allow us to provide comments on the actual implementation of the QMC algorithm.

The result for $w=0.5$ can be seen in Figure 1. All tests use $2^{20}$ evaluations for their final result. The reader should note that instead of having 1 as the common reference value throughout all dimensions, as in [7], the reference values of our setting follow those in Table 1. This is by design to demonstrate some key properties of lattice expansions. Under this setting, results from neighboring dimensions form a pair, for example, dimensions 4 and 5 share the same reference value. As demonstrated in Figure 1, the absolute errors converge to similar limits for these pairs, while the results worsen when the reference increases. This demonstrates that the QMC rules can be generalized to higher dimensions, but their error depends on the value of the target integration, which fits the result from Theorem 2.2.

Finally, it should be noted that the result for dimension 1 is particular strong. This is because the function evaluation for $y_{j}\in[0,1)$, the monomial $y^{k}_{j}$ would always be around zero which helps stabilizing the function evaluation. When we use a projection domain that is larger than the unit box when approximating a polynomial, as we do in the upcoming examples, the results may seem to be worse than standard QMC result and this is one of the contributing factors.

4.2 Normal Distribution

In this subsection, we approximate

[TABLE]

where ${\bf Y}$ follows a multivariate normal distribution with mean [math] for all dimensions and covariance matrix

[TABLE]

Further, we restrict our lattice to the box $[-4.5,4.5]^{s}$ and we study the cases of dimension $s\in[3]$. We applied the full scheme in Equation (2.8) here with $K=2^{7}$. The reference values listed in vector form are $(1.2324,1.4901,1.7705)$.

The convergence result can be seen in Figure 2. All tests use $2^{18}$ evaluations for their final result. The approximation value converges to the references value in all tests, as we would expect. We can observe the structure of the error plateauing for a while with respect to the number of lattice points and suddenly improving when a critical number is reached in Figure 2. This seems to be related to the lattice sequence.

In the case of the normal distribution, there is a final level of error for each dimension, where increasing the number of lattice points further would not improve the approximation result. This is related to the error of projecting a non-periodic function to the half-period cosine space and the truncation of reproducing kernel. These errors dominate when the lattice rule error has converged and limit the final result.

4.2.1 Further studies on parameters

With the successful implementation of our scheme, we conducted further tests for the selection of the kernel truncation parameter and projection domain. We did not consider the generation and construction of different lattice sequences in this article considering the large amount of previous research and work on this topic. It is difficult to draw a satisfying picture on this topic from a simple test. However, we would further consider the choices of projection domain and kernel truncation limit $K$ and their impact on the approximation error. In the following tests, we focus on the 2 dimensional case.

In Figure 3, we see the result of the same test we used for the normal distribution but instead of changing the test dimension $s$, we fixed the test dimension to 2 and varied the projection domain through a variable $L$. The projection domain is defined as $[-0.5L,0.5L]^{2}$ and the kernel truncation $K$ is fixed at $2^{7}$.

From the result, we can see that the convergence behavior is similar among all tests for $L\geq 5$. Our result implies that the lattice sequence has the most significant influence on the convergence behavior with respect to the number of lattice points of our scheme. When most of the probability mass is included in the projection range ($L\geq 5$ here), the projection range only affects the final error. As one can see, the error is the lowest when $L=9$, which can possibly be explained by the expression in Equation (2.9). The first term at the right hand side decreases when the probability mass outside of the projection range increases, while the rest of the terms relate to the averages of the integrand and increases with the projection range. We have to find the balance for these two competing effects when we choose the projection range. Detailed results can be seen in the appendix.

Next, we can see the result of changing the kernel truncation factor $K$ in Figure 4. The projection domain is fixed as $[-4.5,4.5]^{2}$. Once again, the approximation converges to the final limit when the number of lattice points reaches $2^{12}$. This again shows that the kernel truncation is only related to the final limit of approximation but not to the convergence behavior.

Another point of interest for us is the exponential decay of the Fourier transforms (the building blocks of the reproducing kernel). When $K\geq 64$, including more terms in the reproducing kernel does not improve the approximation. This was one of the main advantages of the COS method, making use of the smoothness of the probability measure to achieve better convergence results. It seems that we can arrive at similar results for this combined scheme.

However, the construction for the truncated kernel, defined in Equation (2.7), still suffers from the curse of dimensionality. Under our current setting, this issue can be addressed through software engineering by making use of an online/offline construction, in which we calculate and store the reproducing kernel in an offline stage and then use the stored results in actual calculation. Parallel computing or GPU computing are also possible means to speed up the computation.

We believe that this scheme is valuable for specifically making use of the smoothness of the probability measure and the clear separation of the measure from the integrand, which implies we may switch one part without affecting the other. Possible remedies to be studied for this scheme may include hyperbolic cross or non-linear basis constructions for the truncated kernel, or alternatively, different means for approximating the reproducing kernel.

To conclude, our additional parameter of kernel truncation and domain projection does not change the convergence behavior but only affects the final approximation error when it is sufficiently large. However, the kernel construction still gives rise to the curse of dimensionality.

4.3 Asymmetric Multivariate Laplace Distribution

Finally, we implement our scheme to a typical asymmetric 2-dimensional Laplace distribution ${\bf Y}$, whose character function is given by

[TABLE]

in this subsection. In our setting, $\bar{\boldsymbol{\mu}}$ is a 2-dimensional vector given by $(0.3,-0.1)^{\top}$, while $\boldsymbol{\Sigma}$ is a $2$ matrix defined as

[TABLE]

The mean and covariance of this distribution are given by $\bar{\boldsymbol{\mu}}$ and $\boldsymbol{\Sigma}+\bar{\boldsymbol{\mu}}\bar{\boldsymbol{\mu}}^{\top}$, respectively. For further information, readers are referred to [15]. In this test, we compare the algorithm in Section 2 and the extended COS wavelet scheme in Section 3.2.

In this example, the integrand is simply $Y_{1}Y_{2}$ and the analytic solution for the expectation is $\bar{\mu}_{1}\bar{\mu_{2}}+(\boldsymbol{\Sigma}+\bar{\boldsymbol{\mu}}\bar{\boldsymbol{\mu}}^{\top})_{1,2}=-0.21$.

We pick the projection domain $[\bar{\mu}_{1}-20\Sigma_{1,1},\bar{\mu}_{1}+20\Sigma_{1,1}]\times[\bar{\mu}_{2}-20\Sigma_{2,2},\bar{\mu}_{2}+20\Sigma_{2,2}]$ and set the kernel truncation parameter $K$ as $2^{6}$ for the cosine expansion lattice scheme. For the COS wavelet scheme, the order of the kernel and the number of lattice points are linked. Therefore, when we use $N$ lattice points for the lattice scheme, we pick $N^{\prime}=\left\lceil\sqrt{N}\right\rceil$ for the respective wavelet test. We report in Figure 5.

It is clear that both schemes converge to the same limit. The approximation error is dominated by the measure truncation error when the number of lattice points is high enough and therefore they have similar end results. Comparing both schemes, it seems that the COS wavelet performs better with a similar number of lattice points at the early stage. This may suggest that its quadrature point distribution is more efficient. The difference in quadrature points for the two schemes can be seen in Figure 6.

However, there are still some shortcomings in the construction of the COS wavelet that makes it a less attractive option comparing to the lattice scheme. For example, the COS wavelet is not extendable in the sense that we cannot simply add an extra evaluation point to improve the approximation. A new evaluation from the beginning is required in such case.

Nevertheless, this test suggests that it would be of great interest to implement other quadrature rules with the half-cosine reproducing kernel.

Another point of interest is that in fact the cosine transform decay rate $\beta=2$ for this example, so the rough error bound for the last term in Equation 2.9 fails to converge. However, the success in applying our algorithm clearly shows that we underestimate the decay rate of the cosine transform and a tighter error bound may be derived from advanced results in Fourier analysis. This motivates us to derive better error bound for our scheme.

5 Conclusion

Equipped with the Fourier-cosine projection technique, we extended the lattice sequence from [13] to the broader space of probability measures (instead of just the uniform distribution in [7]). The resulting approximation is reproducible and, in theory, this scheme can be generalized to higher-dimensional space by the extension of the results from the previous lattice literature.

The fact that our scheme concentrates all the information from the probability measure into the reproducing kernel and simply evaluates the integrand at some preset quadrature points gives us flexibility. We may use the same set of quadrature points for different probability measures as our derivation is not measure specific. Moreover, our scheme can be adjusted and possibly improved by simply replacing the lattice sequence by another one.

The main remaining issue is the approximation of the reproducing kernel which still suffers from the curse of dimensionality. However, with the rapid decay of the Fourier transform and the specific form of the reproducing kernel, the authors believe that further research on this topic would be fruitful.

Appendix A Detail results for some numerical tests

A.1 Results for projection domain test

[TABLE]

A.2 Results for kernel truncation test

[TABLE]