Date: 4/01/2016 23:25:01
From: mollwollfumble
ID: 825209
Subject: Integration question?

Am I right in thinking that there are only two basic numerical methods for multidimensional integration of a function?

One, used in finite element and similar, is to divide space up into small pieces and use polynomials on those pieces to integrate and sum.

The second, used in quantum mechanics, is to find the maximum of the function and iterate randomly (Monte-Carlo) around that maximum.

I’m looking into the possibility of two other methods, a fast approximate method (exact when the function is a Gaussian pdf in arbitrarily many dimensions) and a deterministic variant of the Monte-Carlo.

The approximate method relies on the shape of isosurfaces near the function maximum being n-D ellipsoids. The principal axes of these ellipsoids pop out of a Conjugate Gradient method for finding the minimum. So far so good. Now calculate the integral and the n-1th moment along each of the principle axes. That’s only 1-D integration so is extremely fast. This actually gives more than enough information for estimating the integral. The integral can be estimated using the product of the 1-D integrals, with a correction based on the surface area on an n-sphere an a factor based on the average n-1th moment. Or the integral can be estimated in a simpler way – the n-1th moment times the surface area of a unit n-sphere gives the exact integral for every spherically symmetric function – so work from that.

The exact method starts with the Monte-Carlo method (eg. a variant of simulated annealing, or adaptive refinement) and replaces random deviations by deterministic deviations along vectors found from sets of orthogonal polynomials. A deterministic search to given accuracy converges as the square root of the time required for a random search. (A quasi-random search such as a Sobol search takes an intermediate time). So if a random search requires 100 million function evaluations then a deterministic algorithm operates roughly 10 thousand times as fast.

Reply Quote

Date: 5/01/2016 21:21:53
From: The Rev Dodgson
ID: 826016
Subject: re: Integration question?

Did you see:
https://newtonexcelbach.wordpress.com/2014/12/29/numerical-integration-tanh-sinh-quadrature-v-4-3/

The Tanh-Sinh method is an alternative approximate integration method to the Gauss. It’s very well documented in the download file. (I’m allowed to say that because it was written by someone else).

For your second alternative, what advantages do you see to this method over the Gaussian (or similar) method?

Reply Quote

Date: 6/01/2016 18:41:20
From: mollwollfumble
ID: 826690
Subject: re: Integration question?

The Rev Dodgson said:


Did you see:
https://newtonexcelbach.wordpress.com/2014/12/29/numerical-integration-tanh-sinh-quadrature-v-4-3/

The Tanh-Sinh method is an alternative approximate integration method to the Gauss. It’s very well documented in the download file. (I’m allowed to say that because it was written by someone else).

From your link:

“To the best of the author’s knowledge, this Tanh-Sinh and DE Quadrature workbook is the fastest, most powerful, most accurate, most robust and most comprehensive general-purpose quadrature package available today at no cost. It includes full open source code and extensive documentation. The speed and the accuracy of the Tanh-Sinh method in particular and the Double-Exponential (DE) method in general are simply astounding! The DE technique (Ref. 1) combines the fastest speed and highest accuracy of all numerical integration techniques developed thus far.

The quadrature programs are provided as custom functions (UDFs) for performing numerical integration of single-dimension analytic functions.
The programs are grouped below according to the interval of integration, as follows:

1. Finite Interval (a, b)

a. Tanh-Sinh Quadrature. The UDF is named: QUAD_TANH_SINH, and is demonstrated on the TANH-SINH worksheet. Refer to the T-S README, T-S DATA and T-S EXAMPLES worksheets for information on using the program. Use this program as the general-purpose integrator for the finite interval (a, b). b. Gauss-Kronrod Quadrature. The UDF is named: QUAD_GAUSS_KRONROD, and is demonstrated on the G-K worksheet. Refer to the T-S README, T-S DATA and T-S EXAMPLES worksheets for information on using the program. Try using this program as the first alternative should the Tanh-Sinh program fail to integrate a function. c. Recursive Monotone Stable (RMS) Quadrature. The UDF is named QUAD_RMS and is demonstrated on the RMS worksheet. Refer to the T-S README, T-S DATA and T-S EXAMPLES worksheets for information on using the program. Try using this program should the above two programs fail to integrate a function. d. Romberg Quadrature. The UDF is named: QUAD_ROMBERG, and is demonstrated on the ROMBERG worksheet. Refer to the T-S README, T-S DATA and T-S EXAMPLES worksheets for information on using the program. This program is included for historical reasons and for comparisons with the above programs. It is the slowest and least accurate of the four finite-interval programs.

4. Infinite Interval (-inf, inf)

a. Double Exponential (DE) Quadrature for Non-Periodic functions. The UDF name is: QUAD_DE_3, and is demonstrated on the DE-3 worksheet. Refer to the DE README and DE EXAMPLES worksheets for information on using the program. NOTE: This program does not support periodic functions (ie, those with sine, cosine etc terms). Use this program as the general-purpose integrator for the infinite interval (-inf, inf) for non-periodic functions. b. Double Exponential (DE) Quadrature for Periodoc functions. The UDF name is: QUAD_DE_OSC_3, and is demonstrated on the DE-OSC-3 worksheet. Refer to the DE README and DE EXAMPLES worksheets for information on using the program. NOTE: This program does support periodic functions (ie, those with sine, cosine etc terms). Use this program as the general-purpose integrator for the infinite interval (-inf, inf) for periodic functions. c. QAGI Quadrature. The UDF is named: QUAD_QAGI_3 and is demonstrated on the QAGI-3 worksheet. Try this program should either of the above two DE programs fail to integrate a function. ————————————————————

Thanks a million Rev D. I’ll look into these in detail. Particularly those on an infinite domain. Which of these are multi-dimensional and how many dimensions?

The Rev Dodgson said:


For your second alternative, what advantages do you see to this method over the Gaussian (or similar) method?

It’s used in problems that have a very large number of dimensions, it’s widely used in computing path integrals in lattice quantum dynamics. Let’s suppose that you need to integrate something along a path that is free to fluctuate in position. Further suppose that between start and end there are 100 waypoints (in 1-D or 50 in 2-D, 33 in 3-D) that can vary. Then that becomes an integration in 100 dimensions. If you were to set up a grid to do the integration then to get anywhere near sufficient accuracy would require at least 5 grid points in each dimension. The total number of points that need to be sampled is then 5^100, which at 10^70 points is way beyond the capacity of even the most powerful supercomputers, and 100/50/33 isn’t that many waypoints.

The location of the path in lattice quantum dynamics is constrained by a probability function that varies sort of as e^(-x^2), so their strategy is to randomly vary all dimensions at once in such a way that changes have a 50% chance of satisfying e^(-x^2)>0.5. This reduces the integration domain by a huge amount without sacrificing significant accuracy. “x” here is a 100 dimensional vector, and the function is more complicated than -x^2, but you get the idea. To put it another way, the function is very rarely >>0.

For a civil engineering equivalent, suppose you wanted to calculate the average terrain height for a very fast train between Sydney and Canberra when you don’t know the exact route, you only know that the route should have less than a certain amount of cut and fill and a limit on the maximum gradient and minimum bend radius.

Reply Quote

Date: 6/01/2016 19:33:47
From: mollwollfumble
ID: 826766
Subject: re: Integration question?

> The Tanh-Sinh method is an alternative approximate integration method to the Gauss.

How different is this to the use of the type of orthogonal polynomials known as “Hermite Polynomials”? http://mathworld.wolfram.com/OrthogonalPolynomials.html
This also uses a double exponential tail for integrating over the domain from -infinity to infinity.

Or, to put it another way, would the grid points of Tanh-Sinh quadrature be the same as the zeros of Hermite polynomials with the same number of unknowns?

Reply Quote

Date: 7/01/2016 16:25:15
From: mollwollfumble
ID: 827199
Subject: re: Integration question?

mollwollfumble said:


I’m looking into the possibility of two other methods, a fast approximate method (exact when the function is a Gaussian pdf in arbitrarily many dimensions)

The approximate method relies on the shape of isosurfaces near the function maximum being n-D ellipsoids. The principal axes of these ellipsoids pop out of a Conjugate Gradient method for finding the minimum. So far so good. Now calculate the integral and the n-1th moment along each of the principle axes. That’s only 1-D integration so is extremely fast. This actually gives more than enough information for estimating the integral. The integral can be estimated using the product of the 1-D integrals, with a correction based on the surface area on an n-sphere an a factor based on the average n-1th moment. Or the integral can be estimated in a simpler way – the n-1th moment times the surface area of a unit n-sphere gives the exact integral for every spherically symmetric function – so work from that.

Aha!

I’ve found a name for this method. (I think, not confirmed yet). It’s called Sag–Szekeres quadrature. My version may be more or less advanced than the original Sag-Szekeres, but various modified versions already exist. I hadn’t heard of this, its missing from Wikipedia for example.

T.W. Sag, G. Szekeres
Numerical evaluation of high-dimensional integrals
Math. Comp., 18 (1964), pp. 245–253
http://www.ams.org/journals/mcom/1964-18-086/S0025-5718-1964-0165689-X/S0025-5718-1964-0165689-X.pdf

Reply Quote

Date: 7/01/2016 18:42:12
From: mollwollfumble
ID: 827319
Subject: re: Integration question?

mollwollfumble said:


mollwollfumble said:

I’m looking into the possibility of two other methods, a fast approximate method (exact when the function is a Gaussian pdf in arbitrarily many dimensions)

The approximate method relies on the shape of isosurfaces near the function maximum being n-D ellipsoids. The principal axes of these ellipsoids pop out of a Conjugate Gradient method for finding the minimum. So far so good. Now calculate the integral and the n-1th moment along each of the principle axes. That’s only 1-D integration so is extremely fast. This actually gives more than enough information for estimating the integral. The integral can be estimated using the product of the 1-D integrals, with a correction based on the surface area on an n-sphere an a factor based on the average n-1th moment. Or the integral can be estimated in a simpler way – the n-1th moment times the surface area of a unit n-sphere gives the exact integral for every spherically symmetric function – so work from that.

Aha!

I’ve found a name for this method. (I think, not confirmed yet). It’s called Sag–Szekeres quadrature. My version may be more or less advanced than the original Sag-Szekeres, but various modified versions already exist. I hadn’t heard of this, its missing from Wikipedia for example.

T.W. Sag, G. Szekeres
Numerical evaluation of high-dimensional integrals
Math. Comp., 18 (1964), pp. 245–253
http://www.ams.org/journals/mcom/1964-18-086/S0025-5718-1964-0165689-X/S0025-5718-1964-0165689-X.pdf

No.

The paper by Sag and Szekeres seems to only specify a way to generate a uniform grid in an n-sphere up to 20 dimensions, which was the maximum possible by computers in 1964. Then a transformation from the n-sphere onto a hypercube and a simplex. It’s valuable, but more closely resembles my second approach (grid in n-D using either Fourier series or orthogonal polynomials) than my first method (multiplication of 1-D integrals along principle axes).

Reply Quote