Quadrature on sparse grids:

Code to generate and readily evaluated nodes and weights

• Matlab code
• Stata code
• Readily evaluated nodes and weights

Updates

General remarks

To approximate a D-variate integral of the form
∫_Ωg(x)f(x)dx
where Ω⊆R^D and f(x) is a weight function, a quadrature rule prescribes a RxD matrix of nodes X and a Rx1 vector of weights w depending on Ω and f. The approximation is calculated as
g(X)'w.

In all three versions provided here, four quadrature rules for D-dimensional numerical integration are implemented:

• Rules for the unweighted integration: Ω=[0,1]^D, f(x)=1
• GQU: Univariate Gaussian quadrature rules as basis
• KPU: Univariate nested quadrature rules as basis - delayed Kronrod-Patterson rules, see Knut Petras (2003): "Smolyak cubature of given polynomial degree with few nodes for increasing dimension." Numerische Mathematik 93, 729-753.
• Rules for the integration over Ω=R^D with Gaussian weight f(x)=exp(-x'x/2)/sqrt(2π)
• GQN: Univariate Gaussian quadrature rules as basis
• KPN: Univariate nested quadrature rules as basis, see A. Genz and B. D. Keister (1996): "Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight." Journal of Computational and Applied Mathematics 71, 299-309.

As a general rule, the versions with the nested univariate rules are more efficient in the sense that they require fewer nodes to achieve the same level of polynomial exactness. The Matlab code and the Stata code allow the user to use own one-dimensional quadrature rules.

The accuracy level in each of the implementations is denoted as k in the paper. Each rule with accuracy level k integrates complete polynomials of total order 2k-1 exactly.

Matlab code

The Matlab code can be used to extend any univariate quadrature rule to multiple dimensions. The four rules discussed above are included in the code.

To install the package, perform the following steps:

• Download nwspgr.m or download and unpack matlab.zip.
• Make sure Matlab will find the file: Put it into the current working directory, or somewhere into the Matlab search path. You can view and modify it in the menu, see File→Set Path...

The syntax is the following:
[n,w]=nwspgr(type,dim,acc,sym)

• type is a string choosing the underlying univariate quadrature rule. It is either 'GQU', 'KPU', 'GQN', or 'KPN' for the built-in rules, or the name of a function for user-supplied univariate rules. In the latter case, the user must make sure that this function is known to Matlab and that it
  • accepts one scalar argument as the accuracy level
  • returns a column vector of nodes and a column vector of weights.
  for example this function could be defined as function [n w] = myquadrule(acc). In this case, the argument type would be the string 'myquadrule'.
• dim is the number of dimensions.
• acc is the accuracy level.
• sym is optional and only interpreted if a user-supplied univariate quadrature rule is used. In this case, if sym is used and a value other than 0 is supplied, this has two effects: (1) the code runs much faster with high dim and/or acc and (2) the result will only be valid if the underlying univariate rule is "symmetric". This is the case if
  • it uses one node for accuracy level 1, call it m
  • in higher accuracy levels, the nodes cosist of m and/or a set of other nodes. For each of these other nodes m+d, there is also a node m-d with the same weight.
• The results will be returned as a (Rxdim) matrix of nodes n and a (Rx1) vector of corresponding weights w, where R depends on dim and acc, see the examples below.

An example of the program in action is in nwSpGr_demo.m.

Stata code

The Stata code will run in Stata versions 9 and above. There are three ways the nodes and weights can be obtained:

• As a Mata matrix
• As a Stata matrix
• As variables

To install the package, perform the following steps:

• Download stata.zip and unpack its contents to a place where Stata will find the files. A natural place is the personal ado directory. If you don't know where it is, Stata's command sysdir will tell you.
• In order to re-index the installed functions, either close and reopen Stata or use the command mata: mata mlib index.

The syntax is the following:

• For generating nodes and weights as Stata matrices and/or variables, use the command
  nwspgr, type(string) dimensions(integer) accuracy(integer) [ matnodes(string) matweights(string) varnodes(string) varweights(string) ]
  Example: nwspgr, type(KPN) dim(2) acc(4) varn(x) varw(w)
  The first three options determine what is calculated:
  • type(string) chooses the underlying univariate quadrature rule, where string is either GQU, KPU, GQN, or KPN.
  • dim(integer) determines the number of dimensions.
  • acc(integer) determines the accuracy level.
  The other options determine how the results are returned:
  • matn(string) requests that the matrix of nodes is returned as a matrix named string.
  • matw(string) requests that the vector of weights is returned as a matrix named string.
  • varn(string) requests that the matrix of nodes is returned as variables named string1,...,stringD .
  • varw(string) requests that the vector of weights is returned as a variable named string.
  Not specifying any of those options is equivalent to specifying matn(nodes) matw(weights).
• For generating a matrix of nodes and weights from within Stata's matrix language Mata, use the Mata function
  nwspgr(type,dim,acc,sym)
  Example: mata: nw = nwspgr("KPN",2,4)
  • type chooses the underlying univariate quadrature rule. It can be one of the strings "GQU", "KPU", "GQN", or "KPN" if the built-in rules are used. If the user wants to use an own underlying univariate quadrature rule, type has to be a pointer to a function. This function
    • accepts one scalar argument as the accuracy level
    • returns a matrix with two columns, where the first column represents nodes and the second column represents weights.
    for example this function could be declared as real matrix myquadrule(real scalar acc). In this case, the argument type would be the pointer &myquadrule().
  • dim is an integer determining the number of dimensions.
  • acc is an integer determining the accuracy level.
  • sym is optional and only interpreted if a user-supplied univariate quadrature rule is used. In this case, if sym is used and a value other than 0 is supplied, this has two effects: (1) the code runs much faster with high dim and/or acc and (2) the result will only be valid if the underlying univariate rule is "symmetric". This is the case if
    • it uses one node for accuracy level 1, call it m
    • in higher accuracy levels, the nodes cosist of m and/or a set of other nodes. For each of these other nodes m+d, there is also a node m-d with the same weight.
  The result is a Rxdim+1 matrix where the last column represents the weights.

One example for each of the types of use is here:

• Using Mata: demo_mata.do
• Using Stata variables: demo_var.do
• Using Stata Matrices: demo_matrix.do

The source code of the compiled Mata library lnwspgr.mlib contained in the installation ZIP file is not necessary to use the program. For the interest users, here it is: build_nwspgr.do

Readily evaluated nodes and weights

For a variety of combinations of the number of dimensions and the accuracy level, readily generated nodes and weights can be downloaded as comma-seperated ASCII files. The file names are
TYPE_dDIM_lLEV.asc
where TYPE denotes the quadrature rule and is GQU, KPU, GQN, or KPN, DIM is the number of dimensions with 1≤DIM≤20, and LEV is the accuracy level with 1≤LEV≤25. The ASCII files for DIM dimensions have DIM+1 columns, with the last column representing the integration weight.

The following tables show the number of nodes required by the combinations of TYPE, DIM, and LEV. Note that since the number of nodes increases with the number of dimensions, rules for higher dimensions are only implemented for lower accuracy levels. The files are provides as ZIPped archives - one for each TYPE and for each combination of TYPE and DIM.

Rule = KPU

Rule = KPN

Rule = GQU

Rule = GQN

Contact: