Integration, Quasi Mode

| 0 Comments

Last time[1] we saw how it was possible to use uniformly distributed random variables to approximate the integrals of univariate and multivariate functions, being those that take numbers and vectors as arguments respectively. Specifically, since the integral of a univariate function is equal to the net area under its graph within the interval of integration it must equal its average height multiplied by the length of that interval and, by definition, the expected value of that function for a uniformly distributed random variable on that interval is the average height and can be approximated by the average of a large number of samples of it. This is trivially generalised to multivariate functions with multidimensional volumes instead of areas and lengths.
We have also seen[2] how quasi random sequences fill areas more evenly than pseudo random[3] sequences and so you might be asking yourself whether we could do better by using the former rather than the latter to approximate integrals.

Clever you!

Quasi Versus Pseudo Random Means And Standard Deviations

Derivation 1: The Mean And Standard Deviation
By definition, the mean \(\mu\) and variance \(\sigma^2\) are given by the expectations
\[ \begin{align*} \mu &= \mathrm{E}[x]\\ \sigma^2 &= \mathrm{E}\left[(x-\mu)^2\right] = \mathrm{E}\left[x^2 - 2 \mu x + \mu^2\right] \end{align*} \]
The probability density function of the standard uniform distribution is one between zero and one and zero everywhere else and so, by the definition of expected values, we have
\[ \begin{align*} \mu &= \int_0^1 x \; \mathrm{d}x = \left[\frac12x^2\right]_0^1 = \frac12\\ \sigma^2 &= \int_0^1 x^2 - 2 \mu x +\mu^2 \; \mathrm{d}x = \left[\frac13x^3 - \mu x^2 + \mu^2 x\right]_0^1\\ &= \frac13 - \frac12 + \frac14 = \frac4{12} - \frac6{12} + \frac3{12} = \frac1{12} \end{align*} \]
Finally, the standard deviation is the square root of the variance
\[ \sigma = \sqrt{\sigma^2} = \frac1{\sqrt{12}} \]
To get some sense of how effective quasi random sequences might be in the approximation of integrals we can compare how closely the mean and standard deviation of samples of them match those of truly uniformly distributed random variables against how closely those of pseudo random numbers do.
The mean \(\mu\) and standard deviation \(\sigma\) of a standard uniform random variable, having a range of zero to one, are given by
\[ \begin{align*} \mu &= \frac12\\ \sigma &= \frac1{\sqrt{12}} \end{align*} \]
as proven in derivation 1.

We have already implemented a quasi random number generator with our ak.haltonRnd function and program 1 uses samples of numbers generated by both it and Math.random to perform this comparison using the formulae
\[ \begin{align*} \mu &= \frac1n \sum_{i=1}^n u_i\\ \sigma &= \sqrt{\frac1n \sum_{i=1}^n u_i^2 - \mu^2} \end{align*} \]
that we derived last time for the sample mean and standard deviation, in which \(\sum\) is the summation sign and the \(u_i\) are the samples.

Program 1: Means And Standard Deviations

Evidently the quasi random number generator does a much better job of approximating the mean and standard deviation of the standard uniform distribution than the pseudo random number generator does, which should give us some confidence that it will also do better at approximating integrals.

Applying The Central Limit Theorem

For our ak.monteCarloIntegral integrator we estimated the approximation error with the central limit theorem which states that the mean of \(n\) observations of a random variable with a mean of \(\mu\) and a standard deviation of \(\sigma\) will be normally distributed with a mean of \(\mu\) and a standard deviation of \(\sigma\) divided by the square root of \(n\).
Unfortunately, since quasi random sequences don't approximate truly random sequences, we can't use it to estimate the error in the means of samples of them as we would for pseudo random sequences, as demonstrated by program 2 which compares the errors in the means of standard uniform pseudo and quasi random sequences to the standard deviation predicted by the central limit theorem.

Program 2: Pseudo And Quasi Random Errors

If we were to use the central limit theorem to measure how far the mean of a sample of a function of the values of a quasi random sequence might deviate from that function's actual expected value we'd be in danger of greatly overestimating it. As a consequence we'd risk making very many more calls than necessary to any function whose integral we might wish to approximate to some given level of error.

Fortunately Cranley and Patterson[4] came up with a rather neat trick to estimate the errors of a class of integration algorithms that Joe[5] later adapted to algorithms using deterministic sequences like quasi random sequences. The idea is to use a set of pseudo random numbers \(o_j\) to transform a single quasi random sequence \(s_i\) into several with
\[ s^\prime_{j,i} = \left(o_j + s_i\right) \, \% \, 1 \]
where \(x \, \% \, 1\) represents the remainder of \(x\) upon division by one, or in other words its fractional part.
We can now compute expectations using each of these sequences in turn and, since they will be pseudo randomly distributed, we can use the central limit theorem to approximate the error of their mean. For example, we can compute \(m\) approximations of the mean of the standard uniform distribution with
\[ \mu_j = \frac1n \sum_{i=1}^n s^\prime_{j,i} \]
and then use their sample mean and standard deviation of
\[ \begin{align*} \mu &= \frac1m \sum_{j=1}^m \mu_j\\ \sigma &= \sqrt{\frac1m \sum_{j=1}^m \mu_j^2 - \mu^2} \end{align*} \]
to estimate the true mean and the error of that estimation, as done by program 3.

Program 3: Deviation Of Randomised Sequences

Despite using only ten aggregate samples this is clearly significantly more accurate than naively trying to use the central limit theorem directly to estimate the error.

A Quasi Random Integrator

Just as we did for Monte Carlo integration we shall exploit the relationship between integrals and expectations given by
\[ \int_S f(x) \; \mathrm{d}x = \mathrm{E}[x] \times |S| = \mathrm{E}[x \times |S|] \]
where \(S\) represents the space over which the function \(f\) is being integrated, \(\mathrm{E}\) stands for expectation and the vertical bars stand for the size, or volume, of the space enclosed between them.
The integral function, given in listing 1, uses this relationship together with Cranley and Patterson's trick to approximate the integral of the function f.

Listing 1: The integral Function
var N = 16;

function integral(f, threshold, steps, qrnd, prnd, map, v) {
 var i = 0;
 var o = new Array(N);
 var s = new Array(N);
 var j, k, u, s1, s2, fx, s2n2, e1n1;

 if(!isFinite(v)) throw new Error('invalid bounds in ak.quasiRandomIntegral');

 steps = ak.ceil(steps/N);

 for(k=0;k<N;++k) {
  o[k] = prnd();
  s[k] = 0;
 }

 do {
  for(j=0;j<32;++j) {
   u = qrnd();
   for(k=0;k<N;++k) s[k] += f(map(o[k], u))*v;
  }
  i += j;

  s1 = 0;
  s2 = 0;
  for(k=0;k<N;++k) {
   fx = s[k]/i;
   s1 += fx;
   s2 += fx*fx;
  }
  s2n2 = N*s2-s1*s1;
  e1n1 = threshold*(N+Math.abs(s1));
 }
 while(i<steps && s2n2>N*e1n1*e1n1);

 if(!(s2n2<=N*e1n1*e1n1)) {
  throw new Error('failure to converge in ak.quasiRandomIntegral');
 }
 return s1/N;
}

Here, N is the number of pseudo randomly offset quasi random sequences that we use to compute the aggregate samples. Note that we divide the maximum number of steps by N so that we don't treat it as a limit for each of the aggregate samples but rather as a limit for all of them.
The prnd function is a standard uniform pseudo random number generator that is used to initialise their offsets in the array o in the first loop of the function whilst their sums are initialised to zero in the array s. The map function will combine these offsets with values drawn from the qrnd standard uniform quasi random number generator by taking the fractional part of their sum and mapping it to the region of integration. The sums of the results of the function for the arguments that it returns for each offset multiplied by the volume of the region of integration v are accumulated in s in batches of 32 at the start of the main loop of the function. The sample means fx are derived from these sums in the loop that follows it and their sum and sum of their squares are stored in s1 and s2 respectively.
Finally, s1, s2, N and the convergence threshold are used to compute the terms s2n2 and e1n1 that are used to decide whether or not the quasi random integration algorithm has converged by exactly the same argument that we used for our Monte Carlo integration algorithm, given again in derivation 2.

Derivation 2: The Convergence Criterion (Redux)
Given the \(N\) aggregate sample means \(\mu_j\), if we define
\[ \begin{align*} s_1 &= \sum_{j=1}^N \mu_j\\ s_2 &= \sum_{j=1}^N \mu_j^2\\ \end{align*} \]
then their mean and variance are equal to
\[ \begin{align*} \mu &= \frac{s_1}N\\ \sigma^2 &= \frac{s_2}N - \mu^2 \end{align*} \]
and we can use
\[ \sqrt{\max\left(\frac{\sigma^2}N,0\right)} \]
as a measure of the approximation error of \(\mu\). We can again use the normalised error to test for convergence with
\[ \frac{\sqrt{\max\left(\frac{\sigma^2}N,0\right)}}{1+|\mu|} \leqslant \epsilon \]
where the vertical bars stand for the absolute value of the term that they enclose. Rearranging yields
\[ \begin{align*} \sqrt{\max\left(\frac{\sigma^2}N,0\right)} &\leqslant \epsilon \times (1+|\mu|)\\ \max\left(\frac{\sigma^2}N,0\right) &\leqslant \epsilon^2 \times (1+|\mu|)^2 \end{align*} \]
We don't need to set the minimum value of the left hand side to zero since the right hand side is a square and so cannot be negative, and so
\[ \begin{align*} \frac{\sigma^2}N &\leqslant \epsilon^2 \times (1+|\mu|)^2\\ \sigma^2 &\leqslant N \times \epsilon^2 \times (1+|\mu|)^2 \end{align*} \]
Substituting the above definitions of \(\mu\) and \(\sigma^2\) yields
\[ \begin{align*} \frac{s_2}N - \frac{s_1^2}{N^2} &\leqslant N \times \epsilon^2 \times \left(1+\left|\frac{s_1}N\right|\right)^2\\ N \times s_2 - s_1^2 &\leqslant N \times \epsilon^2 \times \left(N+\left|s_1\right|\right)^2 \end{align*} \]
with the final step following from the fact that \(N\) is positive and so is trivially equal to \(|N|\).

Phew!

Listing 2 shows how we can use the integral function for univariate integration.

Listing 2: Univariate Integration
function uniMap(x0, x1) {
 return function(o, u) {
  return x0 + (x1-x0)*((o+u)%1);
 }
}

function uniIntegral(f, x0, x1, threshold, steps, qrnd, prnd) {
 if(ak.nativeType(x1)===ak.UNDEFINED_T) {
  x1 = x0;
  x0 = 0;
 }

 if(ak.nativeType(x0)!==ak.NUMBER_T) {
  throw new Error('invalid lower bound in ak.quasiRandomIntegral');
 }
 if(ak.nativeType(x1)!==ak.NUMBER_T) {
  throw new Error('invalid upper bound in ak.quasiRandomIntegral');
 }

 if(ak.nativeType(qrnd)===ak.UNDEFINED_T) qrnd = ak.haltonRnd(2);

 return integral(f, threshold, steps, qrnd, prnd, uniMap(x0, x1), x1-x0);
}

The uniMap function implements the pseudo random offsetting trick and maps its resulting values from the interval zero to one onto the interval x0 to x1. The uniIntegral follows our usual scheme of changing the bounds of integration to zero to x0 if x1 is undefined. It then checks that both x0 and x1 are numbers and replaces qrnd with a default ak.haltonRnd quasi random generator if it hasn't been defined before finally forwarding to the integral function to create the integrator.

For multivariate integration we again need a helper function to calculate the volume of integration, provided by volume in listing 3.

Listing 3: The volume Helper Function
function volume(x0, x1) {
 var n = x0.dims();
 var v = 1;
 var i;

 for(i=0;i<n;++i) v *= x1.at(i)-x0.at(i);
 return v;
}

Note that we are again using our ak.vector type to represent the bounds of integration x0 and x1 and consequently use its dims and at methods to retrieve the number of dimensions an elements respectively.
Now our multivariate quasi random generator returns arrays and so it's convenient to create an array valued pseudo random number generator with which to offset its results, given in listing 4.

Listing 4: The arrayRnd Generator
function arrayRnd(n, rnd) {
 return function() {
  var x = new Array(n);
  var i;

  for(i=0;i<n;++i) x[i] = rnd();
  return x;
 };
}

The multiMap function given in listing 5 performs the offsetting and mapping for multivariate integration on an element by element basis of the arrays of pseudo and quasi random numbers o and u.

Listing 5: The multiMap Function
function multiMap(x0, x1) {
 var n = x0.dims();

 x0 = x0.toArray();
 x1 = x1.toArray();

 return function(o, u) {
  var x = new Array(n);
  var i;

  for(i=0;i<n;++i) x[i] = x0[i] + (x1[i]-x0[i])*((o[i]+u[i])%1);
  return ak.vector(x);
 }
}

Note that we first convert the ak.vector bounds x0 and x1 to arrays and finally convert the offset and mapped array to an ak.vector to be passed to the function being integrated.

With these functions in place we're ready to use the integral function for multivariate integration, as done in listing 6.

Listing 6: Multivariate Integration
function multiIntegral(f, x0, x1, threshold, steps, qrnd, prnd) {
 var n = x0.dims();
 var base, i;

 if(n===0) {
  throw new Error('invalid lower bound in ak.quasiRandomIntegral');
 }

 if(ak.nativeType(x1)===ak.UNDEFINED_T) {
  x1 = x0;
  x0 = ak.vector(n, 0);
 }
 else if(ak.type(x1)!==ak.VECTOR_T || x1.dims()!==n) {
  throw new Error('invalid upper bound in ak.quasiRandomIntegral');
 }

 if(ak.nativeType(qrnd)===ak.UNDEFINED_T) {
  base = new Array(n);
  for(i=0;i<n;++i) base[i] = ak.primeSequence(i);
  qrnd = ak.haltonRnd(base);
 }

 return integral(f, threshold, steps, qrnd, arrayRnd(n, prnd),
                 multiMap(x0, x1), volume(x0, x1));
}

The first thing that we do here is check that the lower bound x0 has at least one dimension. Next, if the upper bound x1 is undefined then we change the bounds of integration to be from the zero vector to x0, otherwise we check that it is a vector with the correct number of dimensions. If the quasi random generator hasn't been defined, then we set it to a default ak.haltonRnd seeded with consecutive prime numbers generated by our ak.primeSequence function[6]. Finally, we again use the integral function to actually create the integrator.

The ak.quasiRandomIntegral function given in listing 7, and defined in QuasiRandomIntegral.js together with all of the helper functions, performs the final task of selecting between the univariate and multivariate integrals.

Listing 7: ak.quasiRandomIntegral
var DEFAULT_STEPS = Math.pow(2, 22);
var DEFAULT_THRESHOLD = Math.pow(2, -15);

ak.quasiRandomIntegral = function(f, threshold, steps, qrnd, prnd) {
 if(ak.nativeType(f)!==ak.FUNCTION_T) {
  throw new Error('invalid function in ak.quasiRandomIntegral');
 }

 threshold = ak.nativeType(threshold)===ak.UNDEFINED_T ? DEFAULT_THRESHOLD
                                                       : Math.abs(threshold);
 if(isNaN(threshold)) {
  throw new Error('non-numeric threshold passed to ak.quasiRandomIntegral');
 }

 steps = ak.nativeType(steps)===ak.UNDEFINED_T ? DEFAULT_STEPS
                                               : ak.floor(Math.abs(steps));
 if(!(steps>0)) {
  throw new Error('invalid steps passed to ak.quasiRandomIntegral');
 }

 if(ak.nativeType(qrnd)!==ak.UNDEFINED_T && ak.nativeType(qrnd)!==ak.FUNCTION_T) {
  throw new Error('invalid quasi random generator in ak.quasiRandomIntegral');
 }

 if(ak.nativeType(prnd)===ak.UNDEFINED_T) prnd = Math.random;
 else if(ak.nativeType(prnd)!==ak.FUNCTION_T) {
  throw new Error('invalid pseudo random generator in ak.quasiRandomIntegral');
 }

 return function(x0, x1) {
  return ak.type(x0)===ak.VECTOR_T ? multiIntegral(f,x0,x1,threshold,steps,qrnd,prnd)
                                   : uniIntegral(f,x0,x1,threshold,steps,qrnd,prnd);
 };
};

Here we're checking that the function to be integrated is, in fact, a function, giving a default value to the convergence threshold if necessary and checking that it's a number, providing a default number of steps if necessary and checking that it's a non-zero number, checking that the quasi random number generator qrnd is either undefined or a function and setting the pseudo random number generator prnd to Math.random if it isn't defined and otherwise checking that it's a function.
Finally, we return a function that switches between our multivariate and univariate integrators depending upon whether or not the lower bound x0 is an ak.vector.
Note that, in the expectation that quasi random integration will be more accurate than Monte Carlo integration, the default convergence threshold is significantly smaller than the \(2^{-9}\) that we used last time.

Using ak.quasiRandomIntegral

To see our ak.quasiRandomIntegral in action we can compare its results to those of our ak.monteCarloIntegral. Program 4 does so for the exponential function. which has an integral of
\[ \int_{x_0}^{x_1} e^x \; \mathrm{d}x = \Big[e^x\Big]_{x_0}^{x_1} = e^{x_1} - e^{x_0} \]
Program 4: Quasi Versus Pseudo Univariate

Well, the quasi random integration algorithm certainly seems to be living up to our expectations!

For a multivariate example we can again turn to the square root of the sum of the elements a two dimensional vector, which has an exact integral of
\[ \begin{align*} \int_{x_0}^{x_1} \int_{y_0}^{y_1} \sqrt{x + y} \; \mathrm{d}y \; \mathrm{d}x &=\int_{x_0}^{x_1} \left(\int_{y_0}^{y_1} \left(x + y\right)^\tfrac12 \; \mathrm{d}y\right) \; \mathrm{d}x\\ &=\int_{x_0}^{x_1} \bigg[\tfrac23\left(x+y\right)^\tfrac32\bigg]_{y_0}^{y_1} \; \mathrm{d}x\\ &=\int_{x_0}^{x_1} \tfrac23\left(x+y_1\right)^\tfrac32 - \tfrac23\left(x+y_0\right)^\tfrac32 \; \mathrm{d}x\\ &=\bigg[\tfrac4{15}\left(x+y_1\right)^\tfrac52 - \tfrac4{15}\left(x+y_0\right)^\tfrac52\bigg]_{x_0}^{x_1}\\ &=\left(\tfrac4{15}\left(x_1+y_1\right)^\tfrac52 - \tfrac4{15}\left(x_1+y_0\right)^\tfrac52\right) - \left(\tfrac4{15}\left(x_0+y_1\right)^\tfrac52 - \tfrac4{15}\left(x_0+y_0\right)^\tfrac52\right) \end{align*} \]
Program 5 compares the results of our quasi random and Monte Carlo integrators to this exact result.

Program 5: Quasi Versus Pseudo Multivariate

Again, the quasi random algorithm typically yields much better results than the Monte Carlo algorithm and consequently ak.quasiRandomIntegral will be our integrator of choice for high dimensional integrals.

References

[1] Monte Carlo Or Bust, www.thusspakeak.com, 2017.

[2] Halton, Here's A Pseud, www.thusspakeak.com, 2016.

[3] Pseudo Algorithmical, www.thusspakeak.com, 2016.

[4] Cranley, R. and Patterson, T., Randomization of Number Theoretic Methods for Multiple Integration, SIAM Journal on Numerical Analysis, Vol. 13, No. 6, 1976.

[5] Joe, S., Randomization of lattice rules for numerical multiple integration, Journal of Computational and Applied Mathematics, Vol. 31, No. 2, 1990.

[6] A Few Sequences More, www.thusspakeak.com, 2016.

Leave a comment