Recently I uploaded the paper

- J. Dick, Higher order scrambled digital nets achieve the optimal rate of the root mean square error for smooth integrands. ArXiv, 8 July 2010.

This paper deals with a generalization of Owen’s scrambling algorithm which improves on the convergence rate of the root mean square error for smooth integrands. The bound on the root mean square error is best possible (apart from the power of the factor) and this can also be observed from some simple numerical examples shown in the paper (note that the figures in the paper show the standard deviation (or root mean square error) and not the variance of the estimator). In this post you can also find Matlab programs which generate the quadrature points introduced in this paper and a program to generate the numerical results shown in the paper.

The approach used in this paper is similar to the one used in

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

There are some differences though. The concept of higher order digital -nets, introduced in 1. and 2., cannot be used directly. We shall explain the reason for this below.

**The main result**

Let and . We define a norm

For instance, if then the norm above becomes

where denotes the th derivative of We consider integrands for which

Let be an Owen-scrambled digital -net over (for more information on digital nets see Chapter 4 in the book and for more information on scrambling see Chapter 13 in the book).

Define the digit interlacing function for vectors where is the base expansion of by

For we define Our algorithm then simple computes

Since digital nets and scrambling are already implemented in Matlab, one only needs to implement the digit interlacing function to obtain an implementation in Matlab. A simple Matlab program which generates the quadrature points is included in this post (see below).

The main result is now the following.

Theorem

Let such that Let be defined as above. Then the variance of the estimator satisfieswhere denotes the number of quadrature points.

It is known that this result is best possible (except for the power of the factor).

If the smoothness of the integrand is known we can choose and the above result can be slightly improved to

See also the comment in the MathSciNet report of Erich Novak on Art Owen’s paper on local antithetic sampling with scrambled nets, which can be downloaded here.

**Higher order scrambling**

The scrambling method proposed here uses Owen’s scrambling together with the interlacing of digits.

Let be the generating matrices of a digital -net over Then the interlacing of the digits can also be described using the generating matrices. The following picture illustrates how the new first generating matrix is obtained from

The algorithm in the paper first uses a scrambled digital -net and then interlaces the digits of the points. In principle one can of course also use a higher order digital net, use the inverse of the interlacing function as defined in the manuscript, scramble this point set and then apply the interlacing function

The inverse of the interlacing principle can also be applied to the generating matrices of the higher order digital nets. In this case one has generating matrices and then obtains matrices The following picture illustrates how one obtains matrices from

The matrices generate a digital -net over However, even if are the generating matrices of a digital -net, not much can be said about the quality parameter of the digital net generated by Hence, if one would first apply the inverse interlacing function to a higher order digital net, scramble this point set and then apply the interlacing function then one looses the information about the quality parameter of the higher order digital net. For this reason, the scrambling method of the manuscript uses a scrambled digital -net and applies the interlacing function to these sample points directly. (Whether there is a scrambling which can be applied directly to digital higher order nets remains open.)

**Matlab function which generates higher order scrambled digital nets based on Sobol sequence**

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

% Creates a higher order scrambled net based on 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;

P = scramble(P,’MatousekAffineOwen’); % Scramble Sobol points;

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

% Create binary representation of digits;

u=52; % Number of digits used in higher order representation;

depth=floor(u/d); % Number of digits used of original point set;

W = sobolpoints * pow2(depth);

Z = floor(transpose(W));

Y = zeros(s,N);

for j = 1:s,

for i = 1:depth,

for k = 1:d

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

end;

end;

end;

Y = Y * pow2(-depth*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

**Numerical Example**

We consider the function

The following Matlab program estimates the variances using higher order scrambling.

s=2; % Dimension;

mmin=1; % Minimum number of points;

mmax=17; % Maximum number of points;

dmin=1; % Minimum of interlacing factor;

dmax=3; % Maximum of interlacing factor;

rmax=100; % Number of repeats;

%Notice that we must have mmax * dmax < 53.

for d=dmin:dmax

var=zeros(mmax,1);

for m=mmin:mmax

z=zeros(rmax,1);

for r=1:rmax

x=HOscrambleSobol(m,s,d);

z(r)= sum(x(:,2).*exp(x(:,1).*x(:,2))/(exp(1)-2))/(pow2(m)*(exp(1)-2)); % Evaluation of function;

end;

var(m) = sum((z – sum(z)/rmax).^2)/(rmax*(rmax-1)); % Estimation

% of the variance;

d

m

var(m)

end;

loglog(2.^(mmin:mmax),sqrt(var));

hold on;

end;

for d=dmin:dmax

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

end;

hold off;

**Further research**

The construction of higher order polynomial lattice rules is not based on the construction principle of 1. and 2. Hence the question arises how suitable higher order scrambled polynomial lattice rules can be constructed for which the root mean square error converges at the optimal rate. An approach similar to the one used in the recent paper on scrambled polynomial lattice rules, see here, might also work in this case.