# The Nesting Instinct

We have spent some time now[1][2][3] looking at how we might numerically approximate the integrals of functions but have thus far completely ignored the problem of integrating functions with more than one argument, known as multivariate functions.
When solving integrals of multivariate functions mathematically we typically integrate over each argument in turn, treating the others as constants as we do so. At each step we remove one argument from the problem and so must eventually run out of them, at which point we will have solved it.
For example, we can integrate the square root of the sum of $$x$$ and $$y$$ with the steps
\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*}
It is consequently extremely tempting to approximate multivariate integrals by recursively applying univariate numerical integration algorithms to each argument in turn.

### A Bivariate Integrator

We can illustrate the principle with an integrator for bivariate functions, being those that have two arguments. The trick is to implement a helper function that integrates over their second arguments for a fixed value of their first ones and then to integrate it over their first arguments.
In listing 1 that helper function is fx which in turn defines another helper function fy that fixes the value of the first argument of the function f and then integrates it with our ak.rombergIntegral integrator. The biIntegral function then numerically approximates the integral of fx to yield an approximation of the bivariate integral of f.

Listing 1: A Bivariate Integrator
function biIntegral(f, x0, x1, y0, y1) {
function fx(x) {
function fy(y) {return f(x, y);}
return ak.rombergIntegral(fy)(y0, y1);
}
return ak.rombergIntegral(fx)(x0, x1);
}


Program 1 demonstrates that this algorithm is indeed capable of approximating the integral of the square root of $$x$$ plus $$y$$ that we solved above by comparing its result to our exact solution.

Program 1: Using biIntegral

Unfortunately, numerically integrating numerical integrals compounds their approximation errors, which you can see for yourself by setting x0 and y0 to zero in program 1 and noting that the integral fails to converge. At the cost of increased approximation error we can fix this by using a larger convergence threshold, as demonstrated by program 2.

Program 2: Fixing The Convergence Failure

### A Multivariate Integrator

Rather than defining multivariate functions as having long lists of arguments, it is convenient to define them as taking a single vector argument, which we can represent with our ak.vector type. Listing 2 shows how we might do so for the function that we've been using as an example so far.

Listing 2: A Vector Argument Multivariate Function
function f(x) {
return Math.sqrt(x.at(0)+x.at(1));
}


In generalising the bivariate integrator to a multivariate integrator we no longer know in advance how many dimensions we're going to be integrating over and so we must construct the nested integrals at run time. The nested function in listing 3 does so recursively, using an array x as a buffer for the elements of the ak.vector argument of the multivariate function f which is constructed in the function returned in the last step. Note that the functions created at each step populate their elements of that array with their arguments and then integrate over the function created at the next step.
Finally, the multiIntegral function integrates the outermost nested function over the last dimension.

Listing 3: A Multivariate Integrator
function nested(f, x0, x1, x, n) {
if(n===0) {
return function(xn) {
x[0] = xn;
return f(ak.vector(x));
};
}

var g = nested(f, x0, x1, x, n-1);
return function(xn) {
x[n] = xn;
return ak.rombergIntegral(g)(x0.at(n-1), x1.at(n-1));
};
}

function multiIntegral(f, x0, x1) {
var n = x0.dims();
var x = new Array(n);
var g = nested(f, x0, x1, x, n-1);
return ak.rombergIntegral(g)(x0.at(n-1), x1.at(n-1));
}


Program 3 demonstrates that this yields exactly the same result as program 1 for the vector argument version of our test function.

Program 3: Using multiIntegral

Whilst this certainly demonstrates how we can recursively integrate a multivariate function with a univariate integrator, it can hardly be considered a robust implementation given its hard coded use of ak.rombergIntegral and its complete lack of error checking.

### ak.nestedIntegral

Whilst we go about bringing the implementation of the algorithm up to scratch, we might as well put some thought into optimising it. For one thing, we can convert the vector bounds into arrays to slightly reduce the cost of accessing their elements at each step of the recursion. For another, we can move the construction of the integrators outside of the recursively defined functions that use them. The nested function defined in listing 4 exploits both of these optimisations whilst allowing the integrator, its convergence threshold and maximum number of steps to be specified rather than hard coded. Apart from these minor improvements it has exactly the same structure as our prototypical implementation.

Listing 4: An Improved Nested Integral
function nested(integrator, f, threshold, steps, x0, x1, x, n) {
var g, integral;

if(n===0) {
return function(xn) {
x[0] = xn;
return f(ak.vector(x));
};
}

g = nested(integrator, f, threshold, steps, x0, x1, x, n-1);
integral = integrator(g, threshold, steps);
x0 = x0[n-1];
x1 = x1[n-1];

return function(xn) {
x[n] = xn;
return integral(x0, x1);
};
}


Note that if the integrator is ak.trapeziumIntegral, the threshold will be interpreted as the trapezium width and the number of steps will be ignored, so the nested function is compatible with it as well as our other, more sophisticated, numerical integrators.

The multiIntegral function given in listing 5 adds error checking to the arguments. In particular, that the function to be integrated is indeed a function, that the number of dimensions is greater than zero, that the number of steps is a number, that the bounds have the same number of dimensions, and that the integration algorithm produces a function.

Listing 5: An Improved Multivariate Integral
function multiIntegral(integrator, f, threshold, steps, x0, x1) {
var n = x0.dims();
var x, g, integral;

if(ak.nativeType(f)!==ak.FUNCTION_T) {
throw new Error('invalid function in ak.nestedIntegral');
}

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

steps = ak.nativeType(steps)===ak.UNDEFINED_T ? 18 : ak.floor(Math.abs(steps));
if(isNaN(steps)) throw new Error('non-numeric steps passed to ak.nestedIntegral');
steps = ak.ceil(steps/n);

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.nestedIntegral');
}

x0 = x0.toArray();
x1 = x1.toArray();
x = new Array(n);

g = nested(integrator, f, threshold, steps, x0, x1, x, n-1);
integral = integrator(g, threshold, steps);

if(ak.nativeType(integral)!==ak.FUNCTION_T) {
throw new Error('invalid integrator in ak.nestedIntegral');
}

return integral(x0[n-1], x1[n-1]);
}


Here we're providing a default lower bound of all zeros if x1 is undefined and a default maximum number of steps if it isn't specified. Furthermore, if each step of the univariate integrator involves a subdivision into $$k$$ intervals, as they typically will, then in one dimension we have for $$n_s$$ steps
$n_f = k^{n_s}$
function evaluations, whereas in $$d$$ dimensions we have
$n_f = \left(k^{n_s}\right)^d = k^{n_s \times d}$
and so we divide the maximum number of steps by the number of dimensions so that we have the same maximum number of function calls for any number of dimensions.

Finally, the ak.nestedIntegral function given in listing 6 creates a function that uses this to perform the actual integration of multivariate functions.

Listing 6: ak.nestedIntegral
function uniIntegral(integrator, f, threshold, steps, x0, x1) {
var integral = integrator(f, threshold, steps);
if(ak.nativeType(integral)!==ak.FUNCTION_T) {
throw new Error('invalid integrator in ak.nestedIntegral');
}
return integral(x0, x1);
}

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

ak.nestedIntegral = function(integrator) {
if(ak.nativeType(integrator)!==ak.FUNCTION_T) {
throw new Error('invalid integrator in ak.nestedIntegral');
}

return function(f, threshold, steps) {
return nestedIntegral(integrator, f, threshold, steps);
};
};


This checks that the integrator is itself a function and the function that it creates with nestedIntegral is an integrator that switches between the multiIntegral and uniIntegral helper functions depending upon whether the lower bound of integration is a vector or not, with the latter simply integrating the function with the univariate integrator.

This robust implementation of our recursive numerical integration algorithm is defined in NestedIntegral.js and is demonstrated by program 4.

Program 4: Using ak.nestedIntegral

You can see that this can create a multivariate integrator from any of our univariate integrators by changing uni from ak.rombergIntegral to ak.simpsonIntegral, for example, and running the program again.

Unfortunately, recursively applying numerical integrators to numerical integrators gets more and more error prone as the number of dimensions increases and we must consequently increase the number of steps to maintain accuracy. This can get extremely costly for more than say three or four dimensions and so we must turn to alternative approaches for such integrals.

### References

[1] Trapezium Artistry, www.thusspakeak.com, 2015.

[2] The Tip Of The Romberg, www.thusspakeak.com, 2015.

[3] A Jump To The Left And Some Steps To The Right, www.thusspakeak.com, 2016.