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!
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
We have already implemented a quasi random number generator with our
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.
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.
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
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
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.
The
Here,
The
Finally,
Phew!
Listing 2 shows how we can use the
The
For multivariate integration we again need a helper function to calculate the volume of integration, provided by
Note that we are again using our
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.
The
Note that we first convert the
With these functions in place we're ready to use the
The first thing that we do here is check that the lower bound
The
Here we're checking that the function to be integrated is, in fact, a function, giving a default value to the convergence
Finally, we return a function that switches between our multivariate and univariate integrators depending upon whether or not the lower bound
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.
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
Again, the quasi random algorithm typically yields much better results than the Monte Carlo algorithm and consequently
[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.
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}}
\]
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 ourak.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*s2s1*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+\lefts_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 + (x1x0)*((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), x1x0); }
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('nonnumeric 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 nonzero 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 ourak.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