[math-fun] Another multi-D numerical integration idea -- orthogonal arrays and product rules
It turns out the same idea I recently thought of and will describe here, was invented earlier by Greg Kuperberg, and also he found out a version of his idea had been invented even earlier by Nicolas Victoir: Asymmetric cubature formulae with few points in high dimension for symmetric measures, SIAM J. Numer. Anal. 42,1 (2004) 209-227 http://www.math.hu-berlin.de/~finance/papers/victoir2.pdf Greg Kuperberg: Numerical cubature using error-correcting codes, http://arxiv.org/abs/math/0402047 My idea was the following. A "product rule" for integrating a function over an D-dimensional hypercube, is simply to integrate in each dimension one at a time, causing a K-node rule to become a K^D-node rule, with the same degree of exactness t as the 1-dimensonal rule had to start with (e.g. t=2*K-1 for a Gauss rule). View each of the K^D points as indexed by a D-letter word in a K-letter alphabet. Now discard points. If the words corresponding to the points you keep form, combinatorially, an "orthogonal array of strength>=t" (in radix K), then the new slimmed down rule will still have degree of exactness t. (Of course it needs to have all weights rescaled by an overall constant factor to compensate for the discardings.) There are various further games one can play, but that is the basic construction. As a result: THEOREM: the D-cube enjoys (for any fixed t>0) an exactness-degree-t quadrature formula with all weights EQUAL and O(D^floor(t/2)) nodes, all interior to the cube. The constant (which depends upon t) hidden in the O may be large, but for any particular fixed t this result is best possible. The same result is also true, but with a better constant inside the O, if all weights are positive but not necessarily equal. PROOF SKETCH: I'd actually originally shown this for O(D^(c*t)) for various constants c with 0<c. That comes by combining known results on 1-dimensional equal-weight quadrature rules (specifically Gauss rules for unequal weights, and Chebyshev and related rules for equal weights) and known results from the theory of orthogonal arrays as duals of linear error correcting codes within prime fields (specifically BCH codes), and finally known results about gaps between consecutive primes. The fact one can get c down to 1/2, which is minimum possible, was shown by Kuperberg via his own extensions of Chebyshev rules (which I had not known about but which are elegantly suited for this purpose) which enable always regarding the radix as 2 or a power of 2 and thus enable the use of BINARY BCH codes (or any binary codes better than BCH). Some results I have that were not mentioned by Kuperberg include: * You can product up, not 1-dimensional rules, but rather 2-diml or 3-diml basic rules. For example there is a 7-point degree-5 rule in the 2-dimensional square by Albrecht & Collatz (beating the naive 9 point count), yielding after D-fold producting a 2*D-dimensional rule of degree 5 with 7^D points (before slimming). Using an wordsize=D orthogonal array of radix=7 and strength=5 we slim this to a number of points bounded by a QuadraticPolynomial(D). J.Albecht & L.Collatz: Z.Angew.Math.Mech. 38 (1958) 1-15 (German). There also is a 32-point degree-5 rule for the 3-dimensional cube with all point-weights equal, by L.N. Dobrodeev: Cubature rules with equal coefficients for integrating functions with respect to symmetric domains, U.S.S.R. Comput. Math. and Math. Phys. 18,4 (1978) 27-34 (in English) I'm unsure how much this idea can pay off. * You can use the product+slimming method to create a rule "almost" of degree of exactness t, but not quite. That is, it will be exact for all degree-t polynomials in D variables as integrand, except for (say) the univariate monomials x^t or bivariate monomials x^A * y^B for A+B=t. Then you can handle these few exceptions by adjoining a comparatively small number of additional integration nodes (with careful design of their locations and weights). Warning: the weights of these extra nodes can sometimes be negative. For example, Genz & Malik devised a rule in D dimensions with 2^D + 2*D^2 + 2*D + 1 points (centrosymmetric) and degree of exactness 7, where the 2^D points in the leftmost term form the vertices of a hypercube. The 2^D may safely be slimmed down to only the number of words in a radix=2 orthogonal array with strength=6. (Septics are integrated exactly to 0 automatically just by centrosymmetry, so strength=7 not needed.) For example a "Kerdock code" accomplishes that. The 2^D points alone will then integrate every degree<=7 monomial exactly except for x^6 and x^4 * y^2, and Genz & Malik's extra points deal with correcting things in those cases. These extra points can have negative weights but the total sum of |weights| is bounded by 0.082*D*D + O(D). A.C. Genz and A.A. Malik. Remarks on algorithm 006: An adaptive algorithm for numerical integration over an n-dimensional rectangular region, J. of Computational and Applied Mathematics 6,4 (1980) 295-302 A.C. Genz and A.A. Malik: An imbedded family of fully symmetric numerical integration rules, SIAM Journal on Numerical Analysis 20,3 (June 1983) 580-588. The Genz and Malik extra points and weights are: 1 point at (0,0,0...,0) with weight= V*(12824 - 9129*D + 400*D^2 ) / 19683 each 2*D points (+-A, 0, 0, ..., 0) with weight = 980*V/6561 each 2*D points (+-B, 0, 0, ..., 0) with weight = V*(1820-400*D)/19683 each [negative if D>=5] D*(D-1) points +-(C,C,0,0,...,0) with weight = V*200/19683 each A=sqrt(9/70) B=sqrt(9/10) C=sqrt(9/19) V=2^D=volume of hypercube [-1,1]^D. The Genz-Malik "extra points" actually form a degree-5 integration rule by themselves (albeit with different weights needed) and the single point (0,0,0...0) is a degree-1 integration rule by itself, which is very convenient for software error-estimation purposes. I believe, however, that the "slimmed" Genz-Malik rule will be far superior for the purposes they used it for, efficiency-wise, than their original "fat" rule. For example, in D=256 dimensions, there is a binary orthogonal array of strength 6 with 2^23 (about 8 million) words arising from a linear code [256,233,6] over GF(2), which I looked up at http://codetables.de . Hence 2^23 + 2*256^2 + 2*256 + 1 = 8520193 points suffice to get an integration rule in a 256-dimensional hypercube exact for polynomials of degree<=7. (Could we do a bit better with a nonlinear code? Have not checked.) This blows away all previous practical 256-dimensional cubature formulas, far as I can tell (at least if you are willing to live with some negative weights). Apparently the numerical integration software community has not used these ideas (Genz & Malik's unslimmed rule is used in some widely used software, though, e.g. CERN library), which is a pity because they've thereby lost a lot of efficiency. -- Warren D. Smith http://RangeVoting.org <-- add your endorsement (by clicking "endorse" as 1st step)
participants (1)
-
Warren D Smith