QMC rules over R^s: Matlab code and numerical example

In this post you can find a Matlab code for constructing digital nets on \mathbb{R}^s which was recently proposed in J. Dick, Quasi-Monte Carlo numerical integration on \mathbb{R}^s: digital nets and worst-case error. Submitted, 2010. See the previous post where an explanation of the method and a link to the paper can be found. In the numerical example we consider a simple three-dimensional integral. In this example the computation time with the new method is reduced by a factor of ten and additionally the integration error is also reduced. The numerical result in the paper shows that, for the example considered there, that the computation time can be reduced from two and a half minutes to less than two seconds for a certain given error level.

The algorithm

In the following we describe the algorithm. A Matlab implementation is given in the next section. The algorithm is of the form

\displaystyle Q(f) = \sum_{n=0}^{N-1} \lambda_n f(\boldsymbol{x}_n),

where the integrand f:\mathbb{R}^s \to \mathbb{R}, \lambda_n \in [0,\infty), 0 \le n < N are weights and \boldsymbol{x}_n \in \mathbb{R}^s, 0 \le n < N are the quadrature points. We now describe how to obtain the quadrature points and the corresponding weights.

To obtain the quadrature points, we use a digital (t,m,s)-net \boldsymbol{y}_0,\ldots, \boldsymbol{y}_{b^m-1} over a finite field \mathbb{Z}_b of prime order {}b. Let \boldsymbol{y}_n = (x_{y,1},\ldots, x_{y,s}). In the following we consider only the case b=2. Choose the number of quadrature points, i.e. choose m.

The algorithm is illustrated by the following picture.

In step a), we choose a partitioning of the real line \mathbb{R}. How to choose this partition depends on to the problem at hand. For instance, in the numerical example below we chose points a_l = \mathrm{erfinv}(1-2^{-l}) X, where \mathrm{erfinv} denotes the inverse error function. Then the points -a_{m-1}, \ldots, -a_2, -a_1,a_0=0,a_1, a_2,\ldots, a_{m-1} partition the real line into 2m intervals. We now put 2^{m-2} points in the interval [a_0, a_1) and [-a_1,a_0). For the interval [a_0,a_1) these are given by z_k = a_0 + k 2^{-m+2} (a_1-a_0), where 0 \le k < 2^{m-2}. Analoguously for the other interval, which we label by z_{2^{m-2}}, \ldots, z_{2^{m-1}-1} from left to right. We continue putting points in the intervals [a_1,a_2) and [-a_2,-a_1) using the same method. Continuing doing so we obtain points z_0,\ldots, z_{2^m-1} \in \mathbb{R}. This is step b) and c) above.

The last step d) now describes how to obtain the point \boldsymbol{x}_n = (x_{n,1},\ldots, x_{n,s}) \in \mathbb{R}^s. For 0 \le n < 2^m=N and 1 \le i \le s we define

\displaystyle x_{n,i} = z_{2^m y_{n,i}}.

Example
Consider the picture above. Here we have z_0 = 0, z_1 = 1/2, z_2 = -1, z_3 = -1/2, z_4 = 1, z_5 = -2, z_6 = 2, z_7 = -4. Assume we are given the point \boldsymbol{y}_n =  (3/4,5/8) of a digital (0,3,2)-net over \mathbb{Z}_2. Then y_{n,1} = 3/4 and y_{n,2} = 5/8 and therefore 2^m y_{n,1} = 8 \cdot 3/4 = 6 and 2^m y_{n,2} = 8 \cdot 5/8 = 5. Thus x_{n,1} = z_6 = 2 and x_{n,2} = z_5 = -2. Therefore \boldsymbol{x}_n = (2,-2). \Box

Now we consider the construction of the weights \lambda_n. With the construction of the points, each point \boldsymbol{x}_n falls into an interval of the form J = \prod_{i=1}^s [\pm a_{l_i}, \pm a_{l_i \pm 1}). Let N_J denote the number of points which are in the interval {}J, that is N_J = |J \cap \{\boldsymbol{x}_0,\ldots, \boldsymbol{x}_{2^m-1} \}|. Then the weight \lambda_n is simply given by

\displaystyle \lambda_n = \frac{\mathrm{Volume}(J)}{N_J}.

With this definition, constant functions in the interval {}J are integrated exactly. However, it is maybe not so obvious how to actually compute the weights. The volume of J can be computed since it is clear in which interval a point \boldsymbol{x}_n lies. The number of points N_J in {}J is a bit more tricky to compute. To do so, we need Theorem 15 in the paper. For each interval [\pm a_{l_i}, a_{l_i \pm 1}) we know that there are 2^{m_{l_i}} points z_k inside. The number of points in the interval J is then given by N = 2^{m_{\boldsymbol{d}}}, where

\displaystyle m_{\boldsymbol{d}} = m - \sum_{i=1}^s (m-m_{l_i}),

provided that m_{\boldsymbol{d}} \ge t, where {}t is the quality parameter of underlying digital net. Using this formulae we can calculate the number of points in {}J, provided that the interval contains at least 2^t points.

Now we discuss the remaining case of intervals where there are less than 2^t points. The analysis of the integration error shows that points which fall into intervals with less than 2^t points are not significant. Thus, conceivably, giving them weight {}0 would not alter the outcome. More generally, any weight in the interval [0,\mathrm{Volume}(J) N^{-1}_J] is possible. Below we choose those weights to be \mathrm{Volume}(J) 2^{-t}.

In the following section we present a Matlab implementation of the algorithm and we use it to estimate the integral

\displaystyle I = \int_{-\infty}^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty \mathrm{e}^{2\sqrt{\pi} [x+y+z]} \mathrm{e}^{-\pi [x^2+y^2+z^2]} \,\mathrm{d} x \,\mathrm{d} y \,\mathrm{d} z.

This integral can also be computed analytically

\displaystyle I = \left(\int_{-\infty}^\infty \mathrm{e}^{1-(1-\sqrt{\pi} x)^2} \,\mathrm{d} x\right)^3 = \mathrm{e}^3.

This allows us to compute the integration error. Using a substitution, the above integral can also be written as

\displaystyle I = \int_0^1 \int_0^1 \int_0^1 \mathrm{e}^{2 \sqrt{\pi} [\mathrm{erfinv}(2x-1) + \mathrm{erfinv}(2y-1) + \mathrm{erfinv}(2z-1)]} \,\mathrm{d} x \,\mathrm{d} y \,\mathrm{d} z.

We compare our algorithm with the algorithm

\displaystyle \frac{1}{2^m} \sum_{n=0}^{2^m-1}  \mathrm{e}^{2 \sqrt{\pi} [\mathrm{erfinv}(2x_n-1) + \mathrm{erfinv}(2y_n-1) + \mathrm{erfinv}(2z_n-1)]},

where (x_n,y_n,z_n) for 0 \le n < 2^m are the first 2^m points of a Sobol sequence.

Matlab code

% Digital nets in R^s numerical experiment.

format long;

mmin = 13;
mmax = 23;
X=6;

s = 3;
tSobol = 1;

trs=zeros(1,mmax-mmin+1);
tother=zeros(1,mmax-mmin+1);
error_Rs=zeros(1,mmax-mmin+1);
error_other_method=zeros(1,mmax-mmin+1);
Pointcounter=zeros(1,mmax-mmin+1); % Counts the number of points which fall into intervals J for which N_J \ge 2^t.

for m = mmin:mmax

m

N = pow2(m); % Number of points;
P = sobolset(s); % Get Sobol sequence;
sobolpoints = net(P,N); % Get net from Sobol sequence with N points;

tic

% Elementary intervals on R;
Y = erfinv(1-2.^(-(0:1:m)))*X;

% One dimensional projections;
Z=zeros(1,N); % Points in R
W=zeros(1,N); % Weights of interval length in R
M=zeros(1,N); % Stores the values of m_i

R=1;
for k=0:m-2
for l =0:2^(m-2-k)-1
Z(R+l) = Y(k+1) + l*2^(-(m-2-k))*[Y(k+2)-Y(k+1)];
Z(R+2^(m-2-k)+l) = -Y(k+2)+l*2^(-(m-2-k))*[Y(k+2)-Y(k+1)];
W(R+l) = Y(k+2)-Y(k+1);
W(R+2^(m-2-k)+l) = Y(k+2)-Y(k+1);
M(R+l) = m-2-k;
M(R+2^(m-2-k)+l) = m-2-k;
end;
R = R + 2^(m-1-k);
end;
Z(N-1)= Y(m);
Z(N) = -Y(m+1);
W(N-1)= Y(m+1)-Y(m);
W(N) = Y(m+1)-Y(m);
M(N-1)= 0;
M(N)= 0;

%Points in R^s;
sobolindex=sobolpoints*N;
Rspoints=zeros(N,s);
for k=1:N
for l=1:s
Rspoints(k,l)=Z(sobolindex(k,l)+1);
end;
end;

% Volume;
Volume=ones(1,N);
for k=1:N
for l=1:s
Volume(k)=Volume(k)*W(sobolindex(k,l)+1);
end;
end;

% m_i values
Mi=zeros(N,s);
for k=1:N
for l=1:s
Mi(k,l)=M(sobolindex(k,l)+1);
end;
end;

Msum = sum(Mi,2);
Mpoint = m*(1-s)+Msum;

for k=1:N
if Mpoint(k) < tSobol
Mpoint(k)=tSobol; % Gives a small weight to points which lie outside hyperbolic cross
else
Pointcounter(m-mmin+1) = Pointcounter(m-mmin+1)+1;
end;
end;
weight=transpose(Volume).*pow2(-Mpoint);

% Example function f(x,y,z) = exp(2 \sqrt(Pi) * (x+y+z));
correct_result = exp(3);

QRs = 0;
for k=1:N
QRs = QRs + weight(k)*exp(2*sqrt(pi)*[Rspoints(k,1)+Rspoints(k,2)+Rspoints(k,3)]) * exp(-pi*[Rspoints(k,1)^2 + Rspoints(k,2)^2 + Rspoints(k,3)^2]);
end;

error_Rs(m-mmin+1) = QRs – correct_result

trs(m-mmin+1) =toc

% Other method
tic
Q=0;
for k=1:N
Q=Q+ exp(2*sqrt(pi)*[erfinv(2*sobolpoints(k,1)-1) + erfinv(2*sobolpoints(k,2)-1)+ erfinv(2*sobolpoints(k,3)-1)]/sqrt(pi));
end;
Result = Q/N;
error_other_method(m-mmin+1) = Result-correct_result
tother(m-mmin+1) =toc

end

error_Rs
trs
error_other_method
tother

© Josef Dick

This algorithm and program is freely available, under the condition that it cannot be patented by any other party. Please cite the paper ‘J. Dick,  Quasi-Monte Carlo numerical integration on \mathbb{R}^s: digital nets and worst-case error. Submitted, 2010.’ if you make use of it. It is well understood that this algorithm is useful in applications and it is not permitted to obtain a patent for applications of this algorithm or program or modifications thereof.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s