Higher order digital nets and sequences are quasi-Monte Carlo point sets where for digital nets and for digital sequences (for more background information on this topic see Chapters 1-4 and Chapter 8 in the book; several sample chapters of this book can be downloaded here), which satisfy the following property:

where is the smoothness of the integrand and is a constant which only depends on and (but not on For instance, if has square integrable partial mixed derivatives up to order in each variable, then we get a convergence rate of where is the number of quadrature points used in the approximation and can be chosen arbitrarily small. If has derivatives then we get a convergence of and so on.

This method has been introduced in the papers:

- J. Dick, Explicit constructions of quasi-Monte Carlo rules for the numerical integration of high dimensional periodic functions. SIAM J. Numer. Anal., 45, 2141–2176, 2007.
- J. Dick, Walsh spaces containing smooth functions and quasi-Monte Carlo rules of arbitrary high order. SIAM J. Numer. Anal., 46, 1519–1553, 2008.

The first paper only deals with periodic functions, whereas the second paper also includes nonperiodic functions.

For some heuristic explanation how higher order digital nets and sequences work see the paper:

- J. Dick, On quasi-Monte Carlo rules achieving higher order convergence. In: Proceedings of the MCQMC’08 conference, Montreal, Canada, P. L’Ecuyer and A. Owen (eds.), pp. 73–96, 2009. doi: 10.1007/978-3-642-04107-5_5 An earlier version can be found here.

In this entry we provide a Matlab program which generates point sets which satisfy the property (1).

Before we do so, we give some further explanations. The construction of the point set is based on Sobol sequences (which is a particular low-discrepancy sequence) and requires several parameters. First notice that for a Sobol sequence, the base is always The first parameter which needs to be specified in the program is the natural number the number of points will be The second parameter is the dimension (which is a natural number). The third parameter is the interlacing factor Again, is a natural number and should be chosen according to the smoothness of the integrand. So if the integrand has absolutely integrable partial mixed derivatives up to order then choose Of course, might not be known. In this case the following happens:

- If then the method yields a convergence rate of order
- If then the method yields a convergence rate of order

Notice that choosing too large might be disadvantageous (high numerical precision is required and the constants increase rapidly, which can have a stronger effect than the improved rate of convergence for reasonable sample sizes).

Added 29 June 2010: today, the following preprint was made available:

- D. Nuyens and R. Cools, Higher order quasi-Monte Carlo methods: A comparison. Preprint, UNSW Applied Mathematics Reports 2010.

**Matlab function for generating higher order Sobol points**

function [X] = HOSobol(m,s,d)

% Higher order Sobol sequence

% Create a higher order Sobol sequence.

% 2^m number of points

% s dimension of final point set

% d interlacing factor

% X Output Sobol sequence

N = pow2(m); % Number of points;

P = sobolset(d*s); % Get Sobol sequence;

sobolpoints = net(P,N); % Get net from Sobol sequence with N points;

% Create binary representation of digits;

W = sobolpoints* N;

Z = transpose(W);

Y = zeros(s,N);

for j = 1:s,

for i = 1:m,

for k = 1:d

Y(j,:) = bitset( Y(j,:),(m*d+1) – k – (i-1)*d,bitget( Z((j-1)*d+k,:),(m+1) – i));

end;

end;

end;

Y = Y * pow2(-m*d);

X=transpose(Y); % X is matrix of higher order Sobol points,

% where the number of columns equals the dimension

% and the number of rows equals the number of points;

end

* Note* Choosing in the program HOSobol simply yields the classical Sobol sequence.

* Note* The precision required for the optimal performance of this point set is bits. Hence if we have bits available, then we require that

**How to improve higher order digital nets and sequences **

Simple examples show that the method achieves the theoretically predicted convergence rate (see below).

The method runs into problems when the dimension increases (see also the curse of dimensionality). One way to improve the implementation above is by replacing the Sobol points with Niederreiter-Xing sequences (Niederreiter-Xing points can be downloaded from this page). For these points the theoretical dependence on the dimension is better.

Improving higher order digital nets and sequences is an ongoing area of research and quite a number of papers have been published recently on this topic, which can be found here.

In the following we present some very simple numerical examples.

**Example 1: Dimension , Function **

- Matlab program
mmax=13;

dmax=4;for d=1:dmax

z=zeros(mmax);

for m=1:mmax

x=HOSobol(m,1,d);

z(m)=abs((2^(-m)*sum(x)) -1/2);

end;

loglog(2.^(1:mmax),z);

hold on;

end;for d=1:dmax

loglog(2.^(1:mmax),2.^(-d*(1:mmax)),’green’);

end; - This program generates the following figure:
Here, the blue line shows the integration error and the green lines are for One can see that in each case we obtain a convergence rate of as predicted by the theory.

**Example 2: Dimension , Function **

- Matlab program
mmin=1;

mmax=13;

d=2;

dmin=d;

dmax=d;for d=dmin:dmax

z=zeros(mmax);

for m=mmin:mmax

x=HOSobol(m,1,d);

z(m)=abs(2^(-m)*sum(sin(2*pi*x))-0);

end;

loglog(2.^(mmin:mmax),z);

hold on;

end;for d=dmin:dmax

loglog(2.^(mmin:mmax),2.^(-d*(mmin:mmax)),’green’);

end;

hold off; - Numerical results for ; Comparison of error with line (The point Sobol points have some symmetry which yields the small integration error in this case.)
- Numerical results for ; Comparison of error with line
- Numerical results for ; Comparison of error with line
- Numerical results for ; Comparison of error with line

In each case one can approximately see the convergence of the error of order

**Example: Dimension , Function **

This function is interesting as it has absolutely integrable derivatives up to order approximately and hence we can expect a convergence rate of at most This can also be observed in the numerical experiment. Increasing beyond does not have any advantage in this case.

- Matlab program
s=1; % Dimension

mmin=1; % Minimum number of points is 2^(mmin);

mmax=17; % Maximum number of points is 2^(mmax);dmin=1; % Minimum d;

dmax=3; % Maximum d;for d=dmin:dmax

z=zeros(mmax);

for m=mmin:mmax

x=HOSobol(m,s,d);

z(m)=abs(2^(-m)*5/2*sum(x.^(3/2))-1);

end;

loglog(2.^(mmin:mmax),z);

hold on;

end;for d=dmin:dmax

loglog(2.^(mmin:mmax),2.^(-d*(mmin:mmax)),’green’);

end;

hold off; - Numerical results for ;
For we obtain a convergence rate of approximately

- Numerical results for ;
The convergence rate does not increase by increasing beyond We still only get a convergence rate of approximately in this case.

**Example: Dimension , Function **

- Matlab program
s=2; % Dimension

mmin=1; % Minimum number of points is 2^(mmin);

mmax=17; % Maximum number of points is 2^(mmax);dmin=1; % Minimum d;

dmax=3; % Maximum d;for d=dmin:dmax

z=zeros(mmax);

for m=mmin:mmax

x=HOSobol(m,s,d);

z(m)=abs(2^(-m)*5*sum(x(:,1).^(3/2).*x(:,2))-1);

end;

loglog(2.^(mmin:mmax),z);

hold on;

end;for d=dmin:dmax

loglog(2.^(mmin:mmax),2.^(-d*(mmin:mmax)),’green’);

end;

hold off; - Numerical results for ;
Again we obtain a convergence of at most as predicted by the theory.

**Example: Dimension , Function **

- Matlab program
s=5; % Dimension

mmin=1; % Minimum number of points is 2^(mmin);

mmax=23; % Maximum number of points is 2^(mmax);dmin=2; % Minimum d;

dmax=2; % Maximum d;for d=dmin:dmax

z=zeros(mmax);

for m=mmin:mmax

x=HOSobol(m,s,d);

z(m)=abs(2^(-m)*sum(32*x(:,1).*x(:,2).*x(:,3).*x(:,4).*x(:,5))-1);

end;

loglog(2.^(mmin:mmax),z);

hold on;

end;for d=dmin:dmax

loglog(2.^(mmin:mmax),2.^(-d*(mmin:mmax)),’green’);

end;

hold off; - Numerical results for ;
We obtain the expected rate of convergence of (notice that even though the function is infinitely smooth, we have chosen and hence the rate of convergence can at most only be ).

**Plots of higher order Sobol points in dimension **

The following figures show some plots of higher order Sobol sequences in dimension

In Matlab, these plots can be generated for instance using the commands:

m=10; % Change the parameter m here;

s=2; % Change the parameter s here;

d=3; % Change the parameter d here;

X= HOSobol(m,s,d);

plot(X(:,1),X(:,2),’.’);

That’s great explanation. I wonder if you can show us how to plot the sobol points on the surface of a sphere. Thanks a lot!