In this post you can find a Matlab code for constructing digital nets on which was recently proposed in J. Dick, Quasi-Monte Carlo numerical integration on : 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

where the integrand are weights and 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 -net over a finite field of prime order Let In the following we consider only the case Choose the number of quadrature points, i.e. choose

The algorithm is illustrated by the following picture.

In step a), we choose a partitioning of the real line How to choose this partition depends on to the problem at hand. For instance, in the numerical example below we chose points where denotes the inverse error function. Then the points partition the real line into intervals. We now put points in the interval and For the interval these are given by where Analoguously for the other interval, which we label by from left to right. We continue putting points in the intervals and using the same method. Continuing doing so we obtain points This is step b) and c) above.

The last step d) now describes how to obtain the point For and we define

*Example*

Consider the picture above. Here we have Assume we are given the point of a digital -net over . Then and and therefore and Thus and Therefore

Now we consider the construction of the weights With the construction of the points, each point falls into an interval of the form Let denote the number of points which are in the interval that is Then the weight is simply given by

With this definition, constant functions in the interval are integrated exactly. However, it is maybe not so obvious how to actually compute the weights. The volume of can be computed since it is clear in which interval a point lies. The number of points in is a bit more tricky to compute. To do so, we need Theorem 15 in the paper. For each interval we know that there are points inside. The number of points in the interval is then given by where

provided that where is the quality parameter of underlying digital net. Using this formulae we can calculate the number of points in provided that the interval contains at least points.

Now we discuss the remaining case of intervals where there are less than points. The analysis of the integration error shows that points which fall into intervals with less than points are not significant. Thus, conceivably, giving them weight would not alter the outcome. More generally, any weight in the interval is possible. Below we choose those weights to be

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

This integral can also be computed analytically

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

We compare our algorithm with the algorithm

where for are the first 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 for which

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 : 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.