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.
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
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.
% Digital nets in R^s numerical experiment.
mmin = 13;
mmax = 23;
s = 3;
tSobol = 1;
Pointcounter=zeros(1,mmax-mmin+1); % Counts the number of points which fall into intervals for which
for m = mmin:mmax
N = pow2(m); % Number of points;
P = sobolset(s); % Get Sobol sequence;
sobolpoints = net(P,N); % Get net from Sobol sequence with N points;
% 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
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;
R = R + 2^(m-1-k);
Z(N) = -Y(m+1);
W(N) = Y(m+1)-Y(m);
%Points in R^s;
% m_i values
Msum = sum(Mi,2);
Mpoint = m*(1-s)+Msum;
if Mpoint(k) < tSobol
Mpoint(k)=tSobol; % Gives a small weight to points which lie outside hyperbolic cross
Pointcounter(m-mmin+1) = Pointcounter(m-mmin+1)+1;
% Example function f(x,y,z) = exp(2 \sqrt(Pi) * (x+y+z));
correct_result = exp(3);
QRs = 0;
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]);
error_Rs(m-mmin+1) = QRs – correct_result
% Other method
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));
Result = Q/N;
error_other_method(m-mmin+1) = Result-correct_result
© 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.