  We develop a new approach for estimating the expected values of nonlinear
functions applied to multivariate random variables with arbitrary
distributions. Rather than assuming a particular distribution, we assume that
we are only given the first four moments of the distribution. The goal is to
summarize the distribution using a small number of quadrature nodes which are
called $\sigma$-points. We achieve this by choosing nodes and weights in order
to match the specified moments of the distribution. The classical scaled
unscented transform (SUT) matches the mean and covariance of a distribution. In
this paper, introduce the higher order unscented transform (HOUT) which also
matches any given skewness and kurtosis tensors. It turns out that the key to
matching the higher moments is the rank-1 tensor decomposition. While the
minimal rank-1 decomposition is NP-complete, we present a practical algorithm
for computing a non-minimal rank-1 decomposition and prove convergence in
linear time. We then show how to combine the rank-1 decompositions of the
moments in order to form the $\sigma$-points and weights of the HOUT. By
passing the $\sigma$-points through a nonlinear function and applying our
quadrature rule we can estimate the moments of the output distribution. We
prove that the HOUT is exact on arbitrary polynomials up to fourth order.
Finally, we numerically compare the HOUT to the SUT on nonlinear functions
applied to non-Gaussian random variables including an application to
forecasting and uncertainty quantification for chaotic dynamics.
