This paper presents a computationally efficient approach to evaluate multidimensional expectation integrals. Specifically, certain nonproduct cubature points are constructed that exploit the symmetric structure of the Gaussian and uniform density functions. The proposed cubature points can be used as an efficient alternative to the Gauss–Hermite (GH) and Gauss–Legendre quadrature rules, but with significantly fewer number of points while maintaining the same order of accuracy when integrating polynomial functions in a multidimensional space. The advantage of the newly developed points is made evident through few benchmark problems in uncertainty propagation, nonlinear filtering, and control applications.
The problem of uncertainty quantification and uncertainty propagation through nonlinear system for various engineering applications requires the evaluation of multidimension expected value integrals with respect to an appropriate probability density function (PDF). The computation of multidimension expected value integrals constitute an important component in defining the computational complexity associated with many uncertainty quantification and propagation methods like stochastic collocation [1,2], polynomial chaos [3,4], Gaussian mixture model [5–7], etc. Analytical expressions for these multidimension integrals in general exist for linear systems. However, in most other cases, one often does not have analytical solution for these integrals and has to approximate the integral values numerically.
Several computation techniques exist in the literature to approximate the expectation integral with respect to a Gaussian PDF and uniform PDF, the most popular being Monte Carlo (MC) methods , Gauss–Hermite (GH) and Gauss–Legendre (GLgn) quadrature rules , sparse grid quadratures such as Gauss–Hermite Smolyak (GHS) Gauss–Legendre Smolyak (GLgnSM) quadratures , unscented transformation (UT) , and many other Cubature methods . All these methods involve approximating the expectation integral as a weighted sum of integrand values at specified points within the domain of integration. These methods basically differ from each other in the generation of these specific points. For example, MC methods involve random samples from the specified PDF while the Gaussian quadrature scheme involves deterministic points carefully chosen to reproduce exactly the integrals for polynomials.
Quadrature rules in higher dimensions are usually referred to as cubature rules. A cubature rule is said to be exact to degree d, if it can exactly integrate all polynomials with degree ≤ d . The Gaussian quadrature methods provide the minimal number of quadrature points to integrate the polynomial function in a one-dimensional () space. One needs N quadrature points according to the Gaussian quadrature scheme to exactly reproduce the expectation integrals of polynomials with degree 2N − 1 or less in a space . However, one needs to take the tensor product of quadrature points to evaluate multidimensional integrals involving polynomial functions in general space, resulting in a total of Nn quadrature points. Even for a moderate-dimension system involving six random variables, the number of points required to evaluate the expectation integral with only five points along each direction is 56 = 15,625. This is a nontrivial number of points that might make the calculation of the integral computationally expensive, especially when evaluation of the function at each cubature point itself is an expensive procedure, e.g., one may need to solve a system of partial differential equations at each cubature point to compute the function of interest. The sparse grid quadratures or Smolyak quadrature schemes in particular take the sparse product of quadrature rules and thus have fewer points than the equivalent Gaussian quadrature rules at the cost of introducing negative weights. The negative weights are not desired as the error in computing the integral of a nonpolynomial function by making use of any cubature scheme is proportional to the sum of absolute value of weights . In Refs.  and , the authors especially emphasize the necessity of cubature points with positive weights for stable and accurate computations of expectation integrals. Fortunately, the Gaussian quadrature rule is not minimal for n ≥ 2 and there exists cubature rules with reduced number of points . For example, the UT is a degree 3 cubature method with only 2n + 1 points . This forms the motivation for the work presented in this paper.
An extensive amount of work has been performed to develop nonproduct cubature rules, which can exactly reproduce integrals involving polynomials of degree less than or equal to d = 2N − 1 with fewer points, i.e., less than Nn. A good description of nonproduct Gaussian cubature rules, particularly for symmetric density functions, that are second, third, and fifth degree exact in any dimension is provided in Ref. . In general, the development of a minimal cubature rule with positive weights of particular degree d that is applicable to any dimension is still an open problem and often the rules are derived individually for certain degrees and dimensions. In Ref. , the authors provide a method to construct fully symmetric integration formulas of arbitrary degree in space, specifically the formulas for degrees 3, 5, 7, 9 were developed in detail, but all these rules have large negative weights. A cubature rule of degree 2 with n + 1 points and a cubature method of degree 3 with 2n points for a general centrally symmetric weight function (such as the Gaussian and Uniform PDF) in is developed in Ref. . Fully symmetric integration rules with minimal points for that are exact for degree 9–15 were developed in Ref. . A 19-point cubature rule was developed, for symmetric PDFs in that are exact for degree 9 or less in Ref. . For a integral, with symmetric weight function, a 12-point cubature rule for degree 7 was developed in Ref. . It is even claimed that this is the minimal number of points required and there exist many such 12-point cubature rules for degree 7 in systems. A complete account of the cubature rules to date is out of the scope of this paper but a collection of minimal cubature rules are documented in Ref.  and an updated summary of cubature rules can be found in Refs.  and .
With regard to nonlinear filtering, where the integrals with Gaussian probability density functions frequently arise, cubature methods such as the unscented Kalman filter (UKF) with reduced number of points have tremendous potential especially in a real-time scenario. The highly acclaimed UT with only 2n + 1 points is a degree 3 cubature method, where these cubature points are called sigma points . The 2n + 1 UT points can also capture all higher-order moments with odd exponents due to the symmetrical structure of the points. Julier and Uhlmann  derive the simplex sigma points that require only n + 1 points in to capture the first two moments. Unlike the conventional UT, these points do not capture higher-order (≥3) odd moments. Similar to the UT method, a more recent development is the cubature Kalman filter (CKF), which is again exact to degree 3 but uses only 2n points . It can be noticed that the sigma/cubature points of 2n + 1 UT and 2n CKF can be observed as special cases of a more general structure of cubature points previously presented in Ref.  for symmetrical probability density functions. In Ref. , a summary of application sigma points to nonlinear filtering is described.
It is interesting to observe that most of these nonproduct cubature rule possess certain similarities that can be broadly summarized as: (1) they exploit the symmetry of the PDF such as in the Gaussian and Uniform PDFs, (2) a special structure for the cubature points is assumed, and finally (3) a set of nonlinear equations, which are polynomial in nature, are solved. Correspondingly, they also possess some similar drawbacks such as: (1) inconsistency of the set of nonlinear equations, (2) numerical solutions often possess undesirable negative weights and complex, and (3) some cubature points tend to lie outside the support of the PDF. These drawbacks make it very challenging to have a general cubature rules that can be easily constructed for any dimension and order; hence, most cubature rules are constructed for the desired dimension and order. The valuable insight gained by observing the drawbacks of these rules has influenced the development of some new cubature rules with positive weights in this paper.
Building on our prior work [1,22–24], the objective of this paper is to present an unified approach, which we label as conjugate unscented transformation (CUT), to find nonproduct sigma/cubature points with positive weights that can exactly integrate polynomial functions of desired order with respect to the Gaussian and uniform probability density functions. It is desired that the new sigma/cubature rule possesses fewer points than the equivalent Gaussian quadrature product rules. The CUT approach can be considered as the extension of the unscented transformation to provide higher-order sigma points with positive weights for multivariate Gaussian and uniform probability density functions. It is explicitly shown in this work that one cannot capture higher-order moments by considering additional points along the principal direction used in the unscented transformation. This work presents a systematic way to choose conjugate directions to complement principal directions to synthesize symmetric higher-order cubature rules. According to the CUT methodology, the cubature points are chosen to lie symmetrically along specially defined directions. For each cubature point, two unknown variables, a weight and a scaling parameter, are assigned. The moment constraints equations (MCEs) for the desired order are derived in terms of the unknown scaling and weight variables. The primary contribution of this paper is the derivation and documentation of few novel cubature rules for degrees 4, 6, and 8. An analytical solution of degree 4 quadrature points is provided for general space while numerical solutions are presented for degree 6 and 8 quadrature points in a specific dimension for both Gaussian and uniform PDFs. The points have the desired characteristics of all weights being positive and all points lying within the support of the PDF.
The organization of the paper can be summarized as follows: Section 2 introduces the problem of computing multidimensional expectation integrals with quadrature methods. A brief introduction to few popular quadrature rules, particularly the unscented transform, is presented in Sec. 3. An intuitive understanding of how these methods work provides invaluable insight in the development of the CUT methodology. Next, in Sec. 4, the general framework for the CUT approach is described and the newly developed higher-order cubature points are presented. Finally, Sec. 5 presents numerical results illustrating the effectiveness of the CUT approach in evaluating multidimensional expectation integrals relative to conventional approaches in terms of accuracy and computational cost. The application of the CUT methodology to nonlinear filtering and stochastic control is also discussed. Section 6 concludes the paper.
Expectation Integral and Cubature Methods
where represents the order/degree of the moment of the PDF. Alternatively, the MCE in Eq. (4) can be determined by stating that the sigma point set (X, w) can integrate all monomials (hence polynomials) of degree up to and including d. This is known as the moment matching method.
Notice that these constraint equations simply convey that the sigma point set (X, w) should ideally capture all the infinite moments of the PDF. Often, this is not required and one seeks to find a sigma point set (X, w) that can reproduce the first d moments of the PDF exactly to guarantee the exact evaluation of the expectation integral of Eq. (1) when f(x) is a polynomial function of total degree d or less. On the other hand, when integrating nonpolynomial functions, one usually does not know the number of moments to be satisfied but one can expect that higher-order cubature points would have lower error than the lower-order cubature points. These moment constraint equations are high degree multivariate polynomial equations whose solution is not trivial and in some cases even intractable. The Gauss quadrature product rule in indeed satisfies these moment constraint equations theoretically to any prescribed order d by taking tensor product of quadrature points. Hence, Gauss quadrature methods suffer from the curse of dimensionality with the number of points increasing exponentially with dimension. Nevertheless, it is interesting to observe that Gaussian quadratures are symmetric about the mean and thereby inherently satisfy all moments involving any odd exponent, i.e., Ni is odd.
As an alternative to the Gauss quadrature scheme, the nonproduct cubature methods solve these moment constraint equations by judiciously selecting quadrature points also known as sigma points in and offer similar numerical accuracy with fewer points. However, the development of a cubature rule of any given degree d that is applicable to any dimension is still an open problem. The Smolyak quadrature scheme and UT are among the most widely used nonproduct quadrature schemes. The Smolyak quadrature scheme is a general cubature rule of any given degree d and involves sparse tensor product of Gauss quadrature rule. However, it can result in quadrature points with negative weights . The UT is only a third-degree quadrature rule and involves the judicious selection of sigma set by constraining the points to lie on prescribed directions in . The main focus of this work is to extend the UT rules to higher-order quadrature rules by selecting additional axes or directions in . Due to the symmetry of the Gaussian and Uniform probability density functions, the odd moments are always zero and the remaining even order moments are tabulated in Tables 1 and 2. These even order moments can be simply evaluated by taking the product of moments.
These points are selected in such a way that one of the point, x0, lies at the mean of the Gaussian density function and the remaining 2n points, xi, are constrained to lie symmetrically at a distance about the mean on n-symmetrical axes defined by columns of . When P = I (identity matrix), these axes are simply the n-orthogonal axes that are generally known as “principal axes.” Figure 1(a) shows the sigma points and the corresponding length of the stems represent the weights for each point while Fig. 1(b) shows the role of κ that adjusts the distance of the points from the mean and their corresponding weights. Note that the size of the circle represents the weight associated with the sigma point, i.e., as κ increases, the sigma points move farther away from the mean and the associated weights decrease.
Furthermore, the tuning parameter κ is generally chosen such that one of the fourth moment constraint equation error can be minimized, which leads to the constraint of n + κ = 3. This method works well for evaluation of integrals with dimensions 3 or less, after which κ < 0 and the weight corresponding to the central point becomes negative. It is known that cubature rules with negative weight have large errors . Hence, it is always desirable to have a cubature rule with all positive weights.
Tenne and Singh  have developed higher-order unscented transformation methods for a general PDF in , which can satisfy the moment equations up to eighth-order. Lévesque generalizes this to develop unscented transformation for generic integral . However, this method fails to capture cross moments. Next, it is shown that no cross moment of any order can be satisfied by the selection of symmetric points constrained to lie only on the principle axes.
Unscented Transform and Higher-Order Moment Constraint Equations.
It can be seen that the choice of in Eq. (10) leads to the expression for weights in Eqs. (6) and (7). It should be noticed that moments of Eq. (8) containing odd exponents are already satisfied due to the symmetry of 2n + 1 sigma points. Furthermore, the fourth-order moment in Eq. (11) leads to n + κ = 3. The cross moment in Eq. (12) cannot be satisfied by this particular selection of symmetric sigma points that are constrained to lie only on the principal axes.
In fact, no cross moment of any order can be satisfied by the selection of symmetric points constrained to lie only on the orthogonal Cartesian axes. Hence, by merely adding more points or scaling points, they can only minimize the error in the higher-order moments but cannot capture even order cross moments. For a symmetric PDF such as a Gaussian or uniform PDF, moments such as are nonzero when N1, N2,…, Nn are all even numbers. But points on the orthogonal Cartesian axes, which are of the form (±r, 0, 0,…, 0), always produce a zero when substituted into the left-hand side of the MCE for all cross moments.
Also, notice that the CKF is a special case of UT when κ = 0. Hence, the UT and CKF have discrepancy in the fourth-order moments. In summary, the set of nonlinear moment constraint equations is tedious to solve for a general system. Either one needs to break the symmetry of the sigma points or one needs to look for alternative axes to define new sigma points. Breaking the symmetry of sigma points is not desirable since this will require one to include moment equations containing odd exponents in the formulation. In this work, the second option is pursued by including additional axes, which are labeled as conjugate axes such that the resultant points are more “spread out” symmetrically in space.
Conjugate Unscented Transformation
In this section, the conjugate unscented transform is introduced as an extension of the conventional UT to define new sigma points, which can capture higher-order moments. The CUT points have been developed for the Gaussian PDF and Uniform PDF. The points are labeled as CUT-G for the Gaussian PDF and CUT-U for the Uniform PDF. As discussed in Sec. 3, the UT provides the insight that “the points can be constrained to lie on some carefully selected axes” and one only needs to solve for the distances of these points from the origin and their corresponding weights. The set of moment constraint equations along with these additional constraints can make the system tractable for higher-orders and dimensions. However, choosing appropriate axes is still not a trivial problem. Hence, as a first step, the following important axes are defined to select higher-order sigma points.
where and (I)k represents the kth row or column of the n × n dimensional identity matrix I. Each point on the principle axes is at a unit distance from the origin.
The various axes for the case are shown in Fig. 2(a) while Figs. 2(b) and 2(c) show two different perspective views of the first octant for a case. It should be mentioned that all the eight octants in the case are symmetrical.
The next step involves constraining the sigma points to lie on these special axes so that they satisfy higher-order MCE:
As the Gaussian and uniform PDF are symmetric, choosing points on the symmetric axes inherently satisfies all the moment constraint equations involving odd exponents. If the set is a fully symmetric set , then it is closed under the operations of coordinate position and sign permutations. For example, in if , then . The various structures of symmetric sigma point sets used to develop the CUT method are summarized in Table 3. Only one sample point is shown for each type of axes. Given a single point, the fully symmetric set for each type can be easily completed by taking all possible permutations of coordinate position and sign.An example for a system is provided for which the symmetric points in Table 3 are explicitly illustrated in the below equation:Points on the same set of symmetric axes are equidistant from the mean and have equal weights. For the ith set of axes, the distance scaling variables are labeled as ri and weight variables are labeled as wi.(14)
The moment constraint equations up to the desired order are derived in terms of the unknown variables ri and wi by making use of the moments of a Gaussian and Uniform PDF. Each sigma point set in Table 3 can be scaled by a variable ri, assigned a weight wi and substituted into the left hand side of Eq. (4) resulting in a polynomial equation. Hence, each moment constraint equation would contain monomials with variables ri and wi corresponding to each sigma point set. The coefficients of these monomials are tabulated in Table 4 and they can be used directly in generating the moment constraint equations. The triplet [a, b, c] represents the monomial for the ith set of axes, where ri and wi are the unknown variables. Here, n is the dimension of the integral and is the binomial coefficient.
This nonlinear set of equations is solved for ri and wi, yielding the required sigma point set. Usually, for lower-order moment constraint equations, the polynomial equations can be solved analytically or by the help of symbolic computation software. For higher-order moment constraint equations, efficient polynomial solvers such as bertini  can be used.
where yi,j and xi,j are the jth coordinates of the ith point yi and xi, respectively. In Sec. 4.1, certain minimal cubature points are developed for specific order and dimension. Starting with the lowest-order, the axes are judiciously and progressively added to capture the higher-order moments. The method is illustrated primarily for the Gaussian PDF, i.e., for the CUT-G points. As the essential procedure remains the same, only a brief extension is provided for the development of the CUT-U points. The numerical solutions to both CUT-G and CUT-U points are provided in the Appendix.
Fourth-Order Conjugate Unscented Transformation (CUT4).
The unscented transformation provides the minimal number of second degree equivalent quadrature points for the Gaussian PDF. In Ref. , a general second-order rule for any symmetric PDF is provided, which includes UT as the special case. In this section, we derive conjugate unscented transformation sigma points that are fourth moment equivalent, represented by CUT4 and can exactly integrate all polynomials with total degree 5 or less. In this section, we first develop the CUT4 points for the Gaussian PDF and then for the uniform PDF.
Gaussian Probability Density Function.
As the scaling variable r and weight w appear in pairs, the equations in Eq. (17) can be solved by introducing two set of axes or equivalently five variables r1, r2, w1, w2 and the central weight w0. Note that the sum of the weights should be equal to 1, so the central weight, which does not show up in any MCE, is used to sum all the weights to 1. This results in an under-determined problem and an optimization problem can be framed or one can arbitrarily set w0 to 0 and solve for the remaining parameters. Since the conjugate axes provide a nonunique selection of points to solve the moment constraint equations (17), potential benefits of selection of specific axes need to be analyzed.
Table 5 illustrates the growth in the number of points as a function of dimension for various conjugate axes. The number with an underscore refers to the principal and second conjugate axes, which were used by Julier and Uhlmann . Note that for n = 3, selection of the principal axis in conjunction with the third conjugate axis results in the smallest number of points necessary to satisfy the required constraints. Julier and Uhlmann have described an approach in the Appendix of Ref.  to calculate the sigma points that can capture all the fourth-order moments. But this method suffers from the presence of a negative/zero weight for dimensions higher than or equal to 4. It is to be noted that this cubature rule with the similar drawback has been listed in Ref. . To gain more insight, these equations are rederived and solved. Let us consider that one sigma point of weight w0 lies on the origin, 2n points of weight w1 lie symmetrically on each principal axes at a distance of r1, and 2n(n − 1) points of weight w2 lie symmetrically on the second-conjugate axes (c2) at a distance scaled by r2. As the points are symmetric about the origin, the moments involving odd exponents need not be included. The moment constraint equations up to order 4 for this particular selection of points are given as
As a consequence of this simplification, there is only one constraint equation (23) to solve in terms of the variables r1 and r2. Even though it can be solved by minimizing the sixth moment constraint violation, the weight w1 in Eq. (22) would always result in a negative weight for dimension greater than 4. One can address this problem of negative weight by choosing a different set of conjugate axes.
From Table 5, it can be seen that cn axes are the next set of axes with fewer points, which will be used to satisfy the moment constraint equations. The new sigma points are selected such that 2n points of weight w1 lie symmetrically on each principal axis at a distance of r1 and 2n points of weight w2 lie on the nth conjugate axes at a distance of . Note that the boldface number correspond to the fewest quadrature points required to satisfy all the constraints. Substitution of this particular selection of sigma points as summarized in Table 6 leads to the following fourth-order moment constraint equations
Notice that r2 is undefined for n ≤ 2 and thereby this solution is valid for n ≥ 3. Alternatively, one can also find the central weight w0 by minimizing the error in one of the sixth moment constraint equation namely . The solution for dimension n = 2 is shown in Table 7 where the sixth moment constraint equation error has been minimized. In addition, for n = 1, there is no cross order moment (26) and this equation can be replaced by the sixth moment constraint equation leading to the standard Gauss–Hermite quadrature with four points. Note that, according to the CUT4-G, one would need 14 and 1044 sigma points to exactly evaluate the expectation of a polynomial function of degree 4 in three- and ten-dimensional space where as the GH Product rule would require a minimal of 27 and 59,049 points, respectively.
Uniform Probability Density Function.
Following a similar procedure, the CUT4-U points for the uniform PDF can also be developed. But generating the points for the uniform PDF involves another challenge: all the points have to be within the support of the uniform PDF, i.e., within . Hence, rather than developing a general rule for all dimensions, the points and weights have been developed for each specific dimension. As there are three moment constraint equations, two sets of axes have to be chosen. With reference to Table 8, one can begin by choosing the axes with the lowest number of points, until the MCE are satisfied along with the desired characteristics of positive weights. For example, points along the principal axes and nth conjugate axes, cn, constitute minimal number of points for n ≤ 5 while points along principal axes and c4 constitute minimal number of points for n = 6. The appropriate choice for axes is highlighted in Table 8 and is also listed in Table 9. The moment constraint equations up to fourth-order with selection of points underlined in Table 9 for dimensions 2–5 are given as follows:
Unfortunately, after dimension 5, the variable r1 becomes greater than 1 and hence moves out of the support of the uniform PDF. For higher dimensions, a different set of axes are chosen as shown in Table 9. In all these cases, the axes with fewer points were sequentially chosen to satisfy the moment constraint equations up to order 4 and such that the points are within and all have positive weights. The MCE can be easily derived and their analytical solutions are given in Ref. . The final solution of the scaling variables and weights is given in Table 29 of the Appendix.
Sixth-Order Conjugate Unscented Transformation: CUT6.
This section describes the CUT procedure to solve all the moment constraint equations up to sixth-order for both Gaussian and uniform density functions.
Gaussian Probability Density Function.
These equations can be solved by introducing at least six variables r1, r2, r3, w1, w2, w3 or three set of axes to make the system of equation consistent. Table 10 lists the points along different axes with underlined points used by CUT6-G. The axes were chosen sequentially in order to keep the total number of points to a minimum. The sigma points on the principal axes have been assigned a weight of w1 and are constrained to lie symmetrically at a distance r1 from the origin. Points on the nth conjugate axes have been chosen with a weight w2 and are constrained to lie symmetrically at a distance scaled by r2. Finally, the third set of points are selected with weight w3 and are constrained to lie symmetrically at a distance scaled by r3 along the second-conjugate axes. These points are enumerated in Table 11. The set of moment constraint equations using points in Table 11 are given as
This reduced system of equations in Eq. (45) is simpler to solve than the original system of equations given by Eqs. (36)–(42). From Eq. (43), it is evident that w1 < 0 for n > 8. Furthermore, the central weight becomes negative for n > 6. Hence, the solution to Eq. (45) providing quadrature points with positive weights is valid for n ≤ 6. These points are listed in Table 27 of the Appendix.
Furthermore, the problem of negative weight can be avoided by choosing third-conjugate axes instead of second-conjugate axes, which makes the proposed approach valid for 7 ≤ n ≤ 9. The set of sixth-order moment constraint equations in terms of these new sigma points, as enumerated in Table 12, are given as
The solution for r1, r2, and r3 can be computed from the remaining three equations (46)–(48). Notice that this aforementioned set of Eqs. (46)–(51) is preferred till dimension 9 as the central weight becomes negative for n > 9. But if one allows the presence of negative weight, the equations can be solved up to dimension n = 13 after which real roots cannot be guaranteed. The results of evaluating the expectation integral by the CUT6-G method show a contrasting difference in the number of points required. For example, the CUT6-G would need only 49 and 1203 function evaluations while the GH product rule would require a minimal of 256 and 262,144 function evaluations for the exact computation of the expectation integral of a polynomial function of degree 6 in four- and nine-dimensional space, respectively. The solution to these MCEs (36)–(41) and (46)–(51) are provided in Table 27 of the Appendix. Thereby, along with Tables 11 and 12, one can directly generate the sigma point set of order 6. In Table 10, the numbers indicted by an underscore correspond to the number of CUT6-G points required to exactly integrate all polynomials with total degree 7 or less and correspondingly the specific axes, which are required to ensure positive weights for the CUT points.
Uniform Probability Density Function.
Table 13 shows the choice of axes, for the uniform PDF, that would lead to the least number of points within which have positive weights. The axes are selected sequentially from Table 14 to satisfy the moment constraint equations with minimal number of points. The underlined number in Table 14 corresponds to points along the specific axes for which solution to MCEs is found with all positive weights. For n = 2, there are only five moment constraint equations up to order 6 because one of the sixth-order cross moment does not exist in . An example of the MCE, using the choice of axes from Table 13, in the case is
As there are six equations and seven variables, an optimized solution is sought such that the square of the error in the eighth moment is minimized subject to the constraints in Eqs. (53)–(58). The solution to this optimization problem is given in Table 30 of the Appendix. For n ≥ 3, there are six moment constraint equations and one equation for the weights to add up to one. One additional set of axes was used to satisfy the positive weight constraint and guarantee that all the points lie within . When the points are selected as in Table 13, it can be seen that the system is underdetermined with more variables than equations. One can again pose an optimization problem with the objective to minimize the error in higher moment constraint equations. But due to the polynomial nature of the set of MCE, the higher the degree of the equations the higher is the computational complexity. One cannot always guarantee solutions with real or the positive weight or even a solution with all points within the support of the uniform PDF. Hence, appropriate values for some variables are assumed to reduce the computation complexity. The resultant system of equations can easily be solved by numerical polynomial solvers such as bertini . The final numerical solutions for the particular selection of points for dimensions n = 3, 4, 5, 6, 7, 8, 9 in Table 13 are provided in Table 30 of the Appendix.
Eighth-Order Conjugate Unscented Transform: CUT8.
In this section, sigma points are selected such that they satisfy all the moment constraint equations up to order 8 for both Gaussian and uniform PDFs.
Gaussian Probability Density Function.
To solve the system of 11 moment constraint equations, one would require at least 11 variables. As the distance and weight variables appear in pairs, six sets of points including the scaled conjugate axes parameter h would lead to 13 variables. Table 15 shows the axes, which are used for the CUT8-G points as a function of the dimension of the random variables. The specific structure of points chosen to capture these 11 MCEs is shown in Table 16. For instance, this particular selection of points in space leads to the following 11 moment constraint equations:
It should be noticed that there are only eight moment constraint equations in and hence the sigma points corresponding to r5 and w5 are unnecessary; in effect, there are only eight variables r1, r2, r3, r4 and w1, w2, w3, w4. Similarly, there are ten moment constraint equations for system and hence the sigma points corresponding to r5 and w5 are dropped. Using Tables 16 and 28 of the Appendix, the CUT8-G method would need only 355 and 745 function evaluations to compute the expectation integral for a polynomial function of degree 8 in the and space, respectively, whereas the GH product rule would need 3125 and 15,625 function evaluations for the same and space, respectively. It is clear that the number of function evaluations required by the CUT methods is significantly less than that required by the GH quadrature rule while having the same order of accuracy when integrating multivariate polynomial functions.
Uniform Probability Density Function.
The eighth-order conjugate unscented transform points for the uniform PDF are shown in Table 17. The axes are selected sequentially in a judicious way to satisfy the moment constraint equations, have positive weights, and have all the points within the hypercube . There are 11 MCEs corresponding to the 11 even order moments as shown in Table 2. There are only eight MCE and10 MCE for and cases, respectively. For example, eighth-order moment constraint equations for the uniform PDF in space are given as follows:
Comparison to Other Cubature Methods.
For the dimensions and order described in this paper, the CUT-G is compared to equivalent Gauss–Hermite quadrature, Gauss Hermite Smolyak quadrature, and Kronrod Patterson Smolyak quadrature (KPS) [29,30]. The recently developed high-order cubature Kalman filter  for fifth (denoted as HCKF-5) and seventh (denoted as HCKF-7) degrees is also considered for a thorough comparison.
Tables 18–20 show the number of points, N, used by each of the methods considered for specific dimensions. In addition, the stability factor for each method is also shown. The stability factor, as described in Refs.  and , is defined as , i.e., the sum of absolute values of the weights. If some weights are negative, then the stability factor would be greater than 1 leading to a cubature rule or formula that is inaccurate. It is clear that GHS has all negative weights for the dimensions and orders considered. The KPS method has negative weights for n ≥ 5 for fifth-order quadrature rules and negative weights from n ≥ 3 for seventh-order quadrature rules. From Tables 18–20 and Figs. 3(a)–3(c), it is evident that CUT4-G, CUT6-G, and CUT8-G provide minimal number of quadrature points with all positive weights.
In the case of the uniform PDF, Figs. 4(a)–4(d) compare the number of points required to exactly evaluate the expectation integrals of polynomials up to total degree 5, 7, and 9. Figure 4(c), in particular, compares only the CUT6 points and the equivalent Sparse Gauss–Legendre quadrature with four points in each dimension. It can be seen that the CUT-U points, similar to the CUT-G points, are much lower than the equivalent Gauss–Legendre quadrature. Though in some cases the sparse grid (Smolyak scheme) Gauss–Legendre quadrature points achieve fewer points than the CUT-U points, they do so by introducing negative weights. On the other hand, the CUT-U points have all positive weights.
In this section, a few benchmark nonpolynomial functions are considered to illustrate the effectiveness of the proposed CUT method in numerical evaluation of expectation integrals involving Gaussian as well as uniform PDFs. Furthermore, a nonlinear filtering problem and an optimal control problem with uncertain parameters are also considered to show the diversity of applications which can benefit from CUT. The notations such as imply Gauss–Hermite product rule with quadrature points in a single dimension and hence the tensor product would result in a total of points, respectively. Similar notation is adopted for the GHS, GLgn, and GLgnSM points.
Expectation Integral Evaluations
Example 1: Polynomial Function: Gaussian Probability Density Function.
The true value of this integral is 63, and Table 21 lists the numerical value of the integral computed with the help of different methods. Table 21 also lists the number of quadrature points (N), summation of weights, and minimum value of quadrature weight. As the integrand is always non-negative, the integral value should also be always non-negative. As expected, the evaluation of the integral becomes more and more accurate as the order of the cubature rule is increased and cubature rule of order 9 (i.e., GH5, CUT8 and GHS5) exactly reproduces the true integral value. Notice that the fifth-order cubature rule, HCKF–5, in Ref.  estimates a negative value −28.8 due to negative quadrature weights. The presence of negative weights may reduce the accuracy, but can even pose a more significant disadvantage when they are used to integrate positive functions (such as estimating the covariance matrix). This example clearly illustrates the importance of quadrature points with positive weights.
Example 2: Nonpolynomial Function: Gaussian Probability Density Function.
The integral is of dimension n = 6. The true solution as given in Ref.  is −0.543583844. Table 22 lists the integral value obtained by the various methods. The fifth degree cubature rule HCKF–5 with 73 points has negative weights and has 14.72% relative error. In Table 1 of Ref. , the seventh-order higher cubature rule with 584 points achieves 1% relative error while CUT4-G achieves the 1% error with only 76 points as all the weights are positive. Further, CUT6-G achieves 0.3% error with only 137 points, again because all the weights are positive. Sparse grid quadrature that is accurate up to ninth degree achieves the 1% relative error with 1433 points due to the presence of negative weights. The superior performance of the CUT method can again be attributed to all positive weights.
Polar to Cartesian Coordinates (x=r cos θ,y=r sin θ)
This problem of uncertainty transformation under a coordinate change  was introduced to illustrate the effectiveness of the conventional UT. The example is resimulated in Fig. 5(a) while making use of same parameters as described in Ref. . The mean and 1-σ contours are plotted for each method in Fig. 5(a). The truth is taken to be the result of Monte Carlo integration using 3 × 106 random samples in two-dimensional space. All the methods, except CKF, seem to do equally well compared to the MC simulations. CKF predicts the mean correctly but covariance is slightly under estimated due to error in all the fourth moments but UT can capture one fourth-order moment for dimensions ≤ 3. For the same example when the parameters are changed as , the inclusion of higher moment constraints tends to do better than the conventional UT as seen in Fig. 5(b). The CKF slightly under estimates the covariance while the UT and GH3 slightly over estimate the covariance. CUT4-G in 2 can additionally capture one of the sixth-order moments and hence does better than GH3 and UT. The mean and covariance can be evaluated analytically  by the following expressions where and
In Table 23, the various methods are compared to the analytical result and the %-absolute relative error is tabulated for the means and variances clearly illustrating the improved accuracy of the CUT method relative to the other techniques.
Nonpolynomial Function: Uniform Probability Density Function.
This problem was introduced in Ref.  where it had been discussed that computing the expectation of f(x) for negative values of α is a challenging task since α < 0 leads to delta-sequence functions. The convergence of the integral using the GLgn rules for various dimensions is shown in Fig. 6(a) and the converged value for each dimension is used to compute the relative error for other methods. Figure 6(b) shows the relative error using the sparse-GLgn method (GLgnSM). The convergence in the integral is very slow and in higher dimensions such as 8 and 9, the integral does not seem to converge. The error in the integral value first rapidly increases to huge values and then gradually decreases for a very large number of points. This could be due to the presence of large negative weights for many points. Figure 6(c) shows the number of CUT-U points required to achieve at least 10% relative error in contrast to GLgn. For GLgnSM method, the relative error is below 10% only for dimensions 2 and 3 and hence the results are not shown in this figure.
Nonlinear Filtering Examples.
where is the discrete model for the dynamics of the system and can be thought of as the flow of nonlinear continuous time dynamics of the system. is the measurement model, nx and ny being the dimensions of the process model and measurement model, respectively. and are assumed to be zero-mean Gaussian white noise sequences with covariances Qk and Rk+1, respectively. The filtering algorithm is initiated with uncertain initial condition , a Gaussian PDF with mean and covariance matrix and corresponds to computing posterior mean, , and covariance, Pk/k given the measurements up to time k. The notation used in this section is adopted from Ref. .
It should be emphasized that nonlinear filters like the extended Kalman filter (EKF), UKF, CKF, and quadrature Kalman filter (QKF) all make use of the minimum variance formulation. However, each of these methods makes different approximations to evaluate the expectations of the quantities involved. The EKF involves the linearization of both, the system dynamics and measurement model, about the mean and assumes the state PDF to be Gaussian. The EKF has been used successfully in many applications; however, it may result in poor performance due to large errors in state estimates or sparse measurement data . Furthermore, the linearization involves the calculation of Jacobians that could be computationally expensive procedure. On the other hand, nonlinear filters such as UKF, CKF, and QKF do not involve linearization but numerically evaluate the expectation integrals by making use of the specific quadrature scheme. When UT-derived sigma points are used for numerical integration, the resultant filter is called UKF  and when κ = 0, the filter is equivalent to CKF . When Gauss–Hermite quadrature is used, the filter corresponds to QKF  and when Sparse grid quadrature such as Smolyak quadrature is used, the resultant filter is equivalent to the filter developed in Ref. . As discussed earlier, UKF and CKF are only third-order cubature rules and can evaluate the state mean exactly if the process model and measurement model h(xk+1) are polynomials of degree 3 or less. As the degree of the and h(xk+1) increases above 3, errors will be introduced at each stage. Ideally, one would desire to use the Gauss–Hermite product rule as in QKF, but the number of quadrature points grows exponentially with dimension. The minimal cubature rules and CUT sigma points achieve higher accuracy than UKF or CKF with only a fraction of the points used by QKF. Unlike the Smolyak quadratures, the CUT points have all positive weights making it as ideal choice for nonlinear filtering. After each measurement update, the state density function is approximated as Gaussian PDF with mean and covariance Pk+1/k. A new set of quadrature points are generated from this new Gaussian density function, which are propagated through the nonlinear system equations until the next measurement instant. The filter algorithm that directly incorporates these cubature/sigma points by evaluating the expectation integral as a weighted average is described in Fig. 8. In case there are no measurements available at time k + 1, the measurement step is avoided altogether and the algorithm proceeds to the time evolution step. It is interesting to see that the sigma points can be generated from any quadrature method and the filter algorithm would remain unchanged. In Ref. , a polynomial chaos-based filtering algorithm is presented, which makes use of quadrature schemes to update higher-order moments without any assumption on initial density function.
It should be noticed that all the cubature/sigma points can be generated offline and stored at time k = 0 for zero mean and identity covariance Gaussian PDF of appropriate dimension. At each given time k, the required cubature/sigma points can be generated by a simple affine transformation as in Eq. (15).
Air Traffic Tracking Scenario.
A typical air traffic control scenario  is simulated using various cubature/sigma points. The simulation consists of an aircraft that executes the following maneuvers in order: Horizontal motion for 125 s at a speed of 120 m/s toward west, then a coordinated turn (CT) of 90 deg toward south for 90 s with a turn rate of 1 deg/s. The aircraft continues to head south at a speed of 120 m/s for 125 s, finally takes another coordinated turn toward west with turn rate of −3 deg/s for 30 s, thus completing a 90 deg turn and then heading west at a speed of 120 m/s for 125 s. The trajectory of the flight is shown in Fig. 7(a). The kinematics of the turning motion is modeled by the set of nonlinear equations called coordinated turn as described in Ref. . The CT model is characterized by constant speed and constant turn rate. The turn rate Ω is usually unknown and is hence appended to the state vector making the model nonlinear. The system equations for CT model, where the state vector is and radar model that measures the range and bearing are:
In Table 24, the 2-norm of the RMSE in position, velocity, and angular rate is shown for the case when T = 5 s. With significant time between consecutive measurements, the lower-order sigma point filters CKF and UKF seem to diverge with large velocity errors. Further, CUT4-G and CUT6-G also seem to have large velocity errors rendering the estimates unusable. But, CUT8-G has performed significantly better with similar order of errors compared to the particle filter. Further there is a gradual improvement as the order of the filters is increased. CKF and UKF being third-order rules have the maximum error.
A sample trajectory for the Lorenz system with these parameters is shown in Fig. 10. The uncertainty lies in σ, ρ, and the initial conditions. The parameters σ and ρ are considered to be uncertain and hence appended to the state vector [x, y, z]. The uncertainty in the initial condition x = [x, y, z, σ, ρ] is modeled by a Gaussian PDF where
The measurement model consists of the three states h(x) = [x, y, z]T, where the measurements are made every T seconds. The measurement noise covariance is R = Diag([3, 3, 3]). The process noise covariance for the continuous dynamics is Q = 2I5×5. As measurement times are discrete, the process noise covariance for discrete implementation of the dynamics is approximated as Q = 2TI5×5. Between the measurements, the continuous dynamics are used to propagate the particle/sigma/cubature points. The filters PF, CKF, UKF, CUT4-G, CUT6-G, and CUT8-G are simultaneously initiated from the same mean and initial Gaussian PDF. A value of κ = 1 is chosen for the UKF. The bootstrap particle filter  is initiated with 5000 samples and has resampling with a threshold of 0.6 Np where Np is the number of particles. A total of NR = 100 runs are performed for each filter with each run having a different set of random measurements. The RMSE error for each time step is evaluated as where is the estimate from the ith run and is the corresponding truth. The 2-norm and max of RMSE are evaluated as in Eq. (73).
In Fig. 11, the 2-norm of the RMSE and the max norm of the RMSE are plotted for each filter. It can be observed that the particle filter performs the best. The low-order filters UKF and CKF have large errors compared to the higher-order filters. As the measurement time interval T increases, the errors in the estimates of each filter increase, with the lower-order filters performing the worst.
A simulation with T = 0.2 and a low process noise of Q = 0.002I5×5 is also considered for each filter. From Table 25, it can be seen that the particle filter has degraded performance. This problem can be easily corrected by artificially inflating the state covariance at every time step as described in Ref. . A second simulation with a slightly larger process noise of Q = 0.005I5×5 is considered with the same measurement time interval of T = 0.2 s. As expected, in this case, the particle filter outperforms other filters, while higher-order CUT filters are able to perform better than the UKF and the CKF. This inconsistency in the performance of particle filter as a function of process noise is a well-recognized fact and common among sequential Monte Carlo class of filters, which require special remedial actions as discussed in Ref. .
Helicopter Hovering Problem: Uncertainty Quantification.
For simulation's sake, the uncertainty in the first four parameters is assumed to be uniformly distributed within following bounds plbd = [−0.0488, 0.0013, 0.126, −3.3535] and pubd = [−0.0026, 0.0247, 2.394, −0.1765].
Two methods, GLgn and CUT-U rules, are used to propagate the statistical moments up to third-order for the state variables over a period of 15 s. 100,000 MC runs are assumed to provide the ground truth and are used to calculate the relative error. The 2-norms of all the moments at regular time intervals of 0.1 s up to third-order are calculated individually. The 2-norm of the relative error over time is taken for each order of moment individually and the results are shown in Table 26. For each method, N corresponds to the number of points. It can be seen that CUT8-U and GLgn5, being of similar order, achieve the lowest error, but CUT8-U points achieve this with just one-fourth of the number of points used by GLgn5.
Control Under Parametric Uncertainty.
In this example, the rest-to-rest motion of a three mass spring damper system as shown in Fig. 12 is considered. The system is assumed to have uncertainties in both the spring constants and the damping constants, i.e., ρ = [k1, k2, c1, c2]. The uncertainty in ρ is considered to be modeled by a Gaussian PDF . The control input is applied to the first mass. Our objective is to find the control input such that the final residual energy is minimized
Due to uncertainty in the model parameters, the state vector and the cost function J are also uncertain. Hence, we seek a robust controller, which minimizes the maximum residual energy at final time due to uncertainty in system parameters. In Ref. , a polynomial chaos-based approach is described which minimizes the higher-order moments of residual energy to achieve this. This approach utilizes the conventional Gaussian quadrature methods to compute the higher-order moments of the final time residual energy and can be computationally expensive. In this work, we utilize the proposed CUT methodology to reduce the computation burden of this approach. We also compare the performance of this approach with conventional min–max controller.
The optimization problems in Eqs. (81) and (82) share the same set of constraints (83). The various parameters assumed in this simulation are , c = 1, x0 = [0, 0, 0, 0, 0, 0]T and xref(tf) = [1, 1, 1, 0, 0, 0]T. The control u is modeled as a bang–bang profile (of amplitude c = 1) with six unknown switch times. The min–max cost function in Eq. (81) is computed at a set of chosen points ρi within the parameter space with corresponding weight p(ρi). In this simulation, GH5 with 625 points and weights are used in the min–max problem. The first two (M = 2) and four (M = 4) moments of J are minimized in Eq. (82) with the relative weights being (c1 = 0.67, c2 = 0.33) and (c1 = 0.4, c2 = 0.3, c3 = 0.2, c4 = 0.1), respectively. (In particular, the weights were selected as .) The moments here are computed by only 161 CUT8-G points. The results of the optimization problems are shown in Fig. 13, where the optimal/suboptimal solutions for each method are used to evaluate the empirical distribution of J from 100,000 runs starting at random points from the parameter distribution p(ρ). From Fig. 13, it can be observed that the min–max and second-order moment cost function achieve better performance (low mean) but have wider distribution with heavy tail implying lower robustness. On the other hand, the higher-order moments such as the sixth-order moment cost can result in a PDF that is much slimmer, indicating better robustness, but achieves this at the cost of a higher mean indicating lower performance.
Discussion on Computational Complexity and Limitations of Polynomial System of Equations.
The growth of CUT is indeed exponential but with base 2. Notice that this is the minimal exponential growth with integer base that one can have. The proposed CUT points for the dimensions and orders described are much lower (by several factors) than equivalent Gaussian quadrature. Tables 18–20 show the number of points for each method. Further, all the weights are positive. The polynomial system is solved by first solving for the weights symbolically in terms of the other scaling variables ri and then substituted into the remaining moment constraint polynomial equations. This reduces the number of equations but increases the degree of the others that are a function of ri. The order can be further reduced by simple substitutions. The resultant system of equations is easier to solve by bertini  while the original system is often computationally intractable. As one goes to higher dimensions and order, the number of equations increases and the process of symbolically solving for the weights becomes intractable. Hence, the entire polynomial system of equations has to be solved by the polynomial solver, which is often numerically intractable. Further, the nature of the polynomial roots (real and positive) is generally known only after the system is solved. In the case of the uniform PDF, the problem of having all the points within the support of the uniform PDF poses an additional challenge. The axes were judiciously chosen sequentially in order to keep the number of points as low as possible and yet have a computationally feasible solution with all the desired characteristics. Nevertheless, once solved, the points can be saved and directly used. Further, the CUT points are fewer than many standard quadrature rules and offer tremendous potential in application of uncertainty quantification as made evident in numerical examples of Sec. 5. In Ref. , the problem of quantifying the uncertainty of the satellite motion was considered. The CUT8-G points proved to be very efficient in computing the higher-order moments when compared to the Gauss–Hermite and Monte Carlo methods. With just 741 points, the CUT points were able to outperform the 200,000 MC points to compute the moments up to fourth-order. In Refs.  and , the CUT methodology is successfully applied to compute a probabilistic spatial-temporal estimate of ash presence during the April 2010 eruption of the Eyjafjallajökull volcano in Iceland with the help of only 161 simulations of transport and diffusion model.
The examples related to filtering and precision control illustrated in this paper do not include the comparison to other quadratures methods as the required number of points for the same accuracy are presented in Tables 18–20. These tables reveal the relative computational burden of each quadrature method.
The curse of dimensionality, which is inherent in the conventional Gauss-quadrature methods, has motivated numerous researchers to develop nonproduct quadrature rules to embed computational efficiency in evaluating multidimensional expectation integrals. Often, these approaches result in negative weights associated with the quadrature points, which have been shown to result in inaccurate evaluation of the expectation integrals. Although negative weights associated with a particular quadrature approach may not be a problem for a specific application, one cannot be guaranteed of their accuracy a priori in evaluating multidimensional expectational integrals. The unscented transformation, a popular second-order quadrature rule for Gaussian random variables, has a linear growth in quadrature points with respect to the dimension of the uncertain space. Motivated by the unscented-transforms and labeled conjugate unscented transformation, this paper outlines a systematic approach to develop nonproduct quadrature points of desired order with positive weights to evaluate expectation integrals involving Gaussian and uniform density functions. Analytical solutions have been developed for fourth-order quadrature rules for all dimensions to evaluate expectation integrals involving Gaussian as well as uniform density functions. For sixth- and eighth-order quadrature rules, although solutions have been developed only for certain dimensions, they offer tremendous advantage in speed and accuracy for uncertainty quantification in many engineering applications. The developed CUT quadrature rules provide minimal number of quadrature points with positive weights as compared to widely used quadrature methods (e.g., Gauss quadrature and sparse Gauss quadrature methods). Numerous examples representing static mapping involving polynomial and transcendental functions are used to illustrate the efficacy of the developed methodology and provides the evidence in support of the fact that the proposed CUT methodology provides a significant reduction in function evaluations to compute multidimension expectation integrals in comparison to widely used quadrature methods. Application such as uncertainty quantification, nonlinear filtering, and robust control design are also used to illustrate the potential suite of applications that can benefit from the proposed CUT methodology.
Division of Civil, Mechanical and Manufacturing Innovation (Award Nos. CMMI-0908403 and CMMI-1054759).
The solutions of the MCE (36)–(41) and (46)–(51) of order 6 are given in Table 27. The solutions listed are accurate to the tenth decimal point. Solving the MCE accurately to include more decimal places would provide a better result in numerical evaluations. These solutions when used along with Tables 11 and 12 can generate the sigma point sets that are sixth-order accurate for dimensions 2 ≤ n ≤ 9.
The solution for the CUT8 sigma points has been computed by first deriving the MCE of order 8 for each dimension individually by using the CUT8 sigma point selection Table 16 and numerically solving them as described in Sec. 4.3. The solutions of CUT8 for each dimension that are given in Table 28 are tenth decimal accurate.