Holonomic Gradient Method Research Articles

Nonlinear Bayesian filtering via holonomic gradient method with quasi moment generating function

SummaryA symbolic‐numeric method is proposed for addressing the Bayesian filtering problems of a class of discrete‐time nonlinear stochastic systems. We first approximate the posterior probability density function to be Gaussian. The update law of the mean and variance is formulated as the evaluation of several integrals depending on certain parameters. Unlike existing methods, such as the extended Kalman filter (EKF), unscented Kalman filter (UKF), and particle filter (PF), this formulation considers the nonlinearity of system dynamics exactly. To evaluate the integrals efficiently, we introduce an integral transform motivated by the moment generating function (MGF), which we call a quasi MGF. Furthermore, the quasi MGF is compatible with the Fourier transform of differential operators. We utilize this compatibility to decrease the number of computations of Gröbner bases in the noncommutative rings of differential operators, which reduces the offline computational time. A numerical example is presented to show the efficiency of the proposed method compared to that of other existing methods such as the EKF, UKF, and PF.

Holonomic gradient method for two-way contingency tables

The holonomic gradient method gives an algorithm to efficiently and accurately evaluate normalizing constants and their derivatives. We apply the holonomic gradient method in the case of the conditional Poisson or multinomial distribution on two way contingency tables. We utilize the modular method in computer algebra for an efficient and exact evaluation, and we discuss on complexities of these algorithms and their implementation. We also discuss on a theoretical aspect of the distribution from the viewpoint of the conditional maximum likelihood estimation.

Algebraic analysis of rotation data

We develop algebraic tools for statistical inference from samples of rotation matrices. This rests on the theory of D-modules in algebraic analysis. Noncommutative Gr\"obner bases are used to design numerical algorithms for maximum likelihood estimation, building on the holonomic gradient method of Sei, Shibata, Takemura, Ohara, and Takayama. We study the Fisher model for sampling from rotation matrices, and we apply our algorithms for data from the applied sciences. On the theoretical side, we generalize the underlying equivariant D-modules from SO(3) to arbitrary Lie groups. For compact groups, our D-ideals encode the normalizing constant of the Fisher model.

Holonomic extended least angle regression

One of the main problems studied in statistics is the fitting of models. Ideally, we would like to explain a large dataset with as few parameters as possible. There have been numerous attempts at automatizing this process. Most notably, the Least Angle Regression algorithm, or LARS, is a computationally efficient algorithm that ranks the covariates of a linear model. The algorithm is further extended to a class of distributions in the generalized linear model by using properties of the manifold of exponential families as dually flat manifolds. However this extension assumes that the normalizing constant of the joint distribution of observations is easy to compute. This is often not the case, for example the normalizing constant may contain a complicated integral. We circumvent this issue if the normalizing constant satisfies a holonomic system, a system of linear partial differential equations with a finite-dimensional space of solutions. In this paper we present a modification of the holonomic gradient method and add it to the extended LARS algorithm. We call this the holonomic extended least angle regression algorithm, or HELARS. The algorithm was implemented using the statistical software R, and was tested with real and simulated datasets.

The 8-parameter Fisher–Bingham distribution on the sphere

The Fisher–Bingham distribution ( $$\mathrm {FB}_8$$ ) is an eight-parameter family of probability density functions (PDF) on the unit sphere that, under certain conditions, reduce to spherical analogues of bivariate normal PDFs. Due to difficulties in its interpretation and estimation, applications have been mainly restricted to subclasses of $$\mathrm {FB}_8$$ , such as the Kent ( $$\mathrm {FB}_5$$ ) or von Mises–Fisher (vMF) distributions. However, these subclasses often do not adequately describe directional data that are not symmetric along great circles. The normalizing constant of $$\mathrm {FB}_8$$ can be numerically integrated, and recently Kume and Sei showed that it can be computed using an adjusted holonomic gradient method. Both approaches, however, can be computationally expensive. In this paper, I show that the normalization of $$\mathrm {FB}_8$$ can be expressed as an infinite sum consisting of hypergeometric functions, similar to that of the $$\mathrm {FB}_5$$ . This allows the normalization to be computed under summation with adequate stopping conditions. I then fit the $$\mathrm {FB}_8$$ to two datasets using a maximum-likelihood approach and show its improvements over a fit with the more restrictive $$\mathrm {FB}_5$$ distribution.

On the fractional moments of a truncated centered multivariate normal distribution

In this paper, we study the fractional moments of a truncated centered multivariate normal distribution, with a focus on their computation. We develop computational methods, including ones based on the holonomic gradient method, the second-order Laplace approximation, and the Monte Carlo method. These methods enable us to compute higher order fractional moments without evaluating multiple integrals. Via numerical experiments, we investigate their performances. Some applications, including robust graphical modeling based on the alternative multivariate t-distribution, are also presented.

Holonomic properties and recurrence formula for the distribution of sample correlation coefficient

This paper presents the holonomic properties and recurrence formula for the distribution of the sample correlation coefficient. The probability density function (pdf) is holonomic. Therefore, it is computed exactly based on the holonomic gradient method (HGM). The initial values for computation are expressed in terms of Gaussian hypergeometric functions with specific parameters that can be transformed to a rational equation of gamma functions. Using the integral algorithm in the D-module theory, the cumulative distribution function (cdf) is also holonomic. It can be computed using HGM. Next, we derive the recurrence formula for the Gaussian hypergeometric function related to the degrees of freedom and apply it to exact computation of the pdf under a fixed population correlation coefficient and increasing degrees of freedom. We conclude with discussion of the quantile function of the sample correlation coefficient which satisfies a nonlinear differential equation of second order.

A-hypergeometric distributions and Newton polytopes

We give a bijection between a quotient space of the parameters and the space of moments for any A-hypergeometric distribution. An algorithmic method to compute the inverse image of the map is proposed utilizing the holonomic gradient method and an asymptotic equivalence of the map and the iterative proportional scaling. The algorithm gives a method for solving a conditional maximum likelihood estimation problem in statistics. The interplay between the theory of hypergeometric functions and statistics allows us to give some new formulas for A-hypergeometric polynomials.

Distribution of the ratio of two Wishart matrices and cumulative probability evaluation by the holonomic gradient method

We study the distribution of the ratio of two central Wishart matrices with different covariance matrices. We first derive the density function of a particular matrix form of the ratio and show that its cumulative distribution function can be expressed in terms of the hypergeometric function 2F1 of a matrix argument. Then we apply the holonomic gradient method for numerical evaluation of the hypergeometric function. This approach enables us to compute the power function of Roy’s maximum root test for testing the equality of two covariance matrices.

Approximate maximum likelihood estimation of the Bingham distribution

Maximum likelihood estimation of the Bingham distribution is difficult because the density function contains a normalization constant that cannot be computed in closed form. Given the availability of sufficient statistics, Approximate Maximum Likelihood Estimation (AMLE) is an appealing method that allows one to bypass the evaluation of the likelihood function. The impact of the input parameters of the AMLE algorithm is investigated and some methods for choosing their numerical values are suggested. Moreover, AMLE is compared to the standard approach which numerically maximizes the (approximate) likelihood obtained with the normalization constant estimated via the Holonomic Gradient Method (HGM). For the Bingham distribution on the sphere, simulation experiments and real-data applications produce similar outcomes for both methods. On the other hand, AMLE outperforms HGM when the dimension increases.

Exact ZF Analysis and Computer-Algebra-Aided Evaluation in Rank-1 LoS Rician Fading

We study zero-forcing detection (ZF) for multiple-input/multiple-output (MIMO) spatial multiplexing under transmit-correlated Rician fading for an N_R X N_T channel matrix with rank-1 line-of-sight (LoS) component. By using matrix transformations and multivariate statistics, our exact analysis yields the signal-to-noise ratio moment generating function (m.g.f.) as an infinite series of gamma distribution m.g.f.'s and analogous series for ZF performance measures, e.g., outage probability and ergodic capacity. However, their numerical convergence is inherently problematic with increasing Rician K-factor, N_R , and N_T. We circumvent this limitation as follows. First, we derive differential equations satisfied by the performance measures with a novel automated approach employing a computer-algebra tool which implements Groebner basis computation and creative telescoping. These differential equations are then solved with the holonomic gradient method (HGM) from initial conditions computed with the infinite series. We demonstrate that HGM yields more reliable performance evaluation than by infinite series alone and more expeditious than by simulation, for realistic values of K , and even for N_R and N_T relevant to large MIMO systems. We envision extending the proposed approaches for exact analysis and reliable evaluation to more general Rician fading and other transceiver methods.

Holonomic gradient method for distribution function of a weighted sum of noncentral chi-square random variables

We apply the holonomic gradient method to compute the distribution function of a weighted sum of independent noncentral chi-square random variables. It is the distribution function of the squared length of a multivariate normal random vector. We treat this distribution as an integral of the normalizing constant of the Fisher–Bingham distribution on the unit sphere and make use of the partial differential equations for the Fisher–Bingham distribution.

MIMO Zero-Forcing Performance Evaluation Using the Holonomic Gradient Method

For multiple-input multiple-output (MIMO) spatial-multiplexing transmission, zero-forcing detection (ZF) is appealing because of its low complexity. Our recent MIMO ZF performance analysis for Rician--Rayleigh fading, which is relevant in heterogeneous networks, has yielded for the ZF outage probability and ergodic capacity infinite-series expressions. Because they arose from expanding the confluent hypergeometric function $ {_1\! F_1} (\cdot, \cdot, \sigma) $ around 0, they do not converge numerically at realistically-high Rician $ K $-factor values. Therefore, herein, we seek to take advantage of the fact that $ {_1\! F_1} (\cdot, \cdot, \sigma) $ satisfies a differential equation, i.e., it is a \textit{holonomic} function. Holonomic functions can be computed by the \textit{holonomic gradient method} (HGM), i.e., by numerically solving the satisfied differential equation. Thus, we first reveal that the moment generating function (m.g.f.) and probability density function (p.d.f.) of the ZF signal-to-noise ratio (SNR) are holonomic. Then, from the differential equation for $ {_1\! F_1} (\cdot, \cdot, \sigma) $, we deduce those satisfied by the SNR m.g.f. and p.d.f., and demonstrate that the HGM helps compute the p.d.f. accurately at practically-relevant values of $ K $. Finally, numerical integration of the SNR p.d.f. produced by HGM yields accurate ZF outage probability and ergodic capacity results.

Calculation of orthant probabilities by the holonomic gradient method

We apply the holonomic gradient method (Nakayama et al. Adv Appl Math 47:639–658, 2011) to the calculation of the probabilities of the multivariate normal distribution. The holonomic gradient method applied to finding the orthant probabilities is found to be a variant of Plackett’s recurrence relation. However, an implementation of the method yields recurrence relations that are more suitable for numerical computation than is Plackett’s recurrence relation. We derive some theoretical results on the holonomic system for the orthant probabilities. These results show that multivariate normal orthant probabilities possess some remarkable properties from the viewpoint of holonomic systems. Finally, we show that the numerical performance of our method is comparable or superior to that of existing methods.

Calculating the normalising constant of the Bingham distribution on the sphere using the holonomic gradient method

In this paper we implement the holonomic gradient method to exactly compute the normalising constant of Bingham distributions. This idea is originally applied for general Fisher–Bingham distributions in Nakayama et al. (Adv. Appl. Math. 47:639–658, 2011). In this paper we explicitly apply this algorithm to show the exact calculation of the normalising constant; derive explicitly the Pfaffian system for this parametric case; implement the general approach for the maximum likelihood solution search and finally adjust the method for degenerate cases, namely when the parameter values have multiplicities.

The holonomic gradient method for the distribution function of the largest root of a Wishart matrix

We apply the holonomic gradient method introduced by Nakayama et al. (2011) [23] to the evaluation of the exact distribution function of the largest root of a Wishart matrix, which involves a hypergeometric function 1F1 of a matrix argument. Numerical evaluation of the hypergeometric function has been one of the longstanding problems in multivariate distribution theory. The holonomic gradient method offers a totally new approach, which is complementary to the infinite series expansion around the origin in terms of zonal polynomials. It allows us to move away from the origin by the use of partial differential equations satisfied by the hypergeometric function. From the numerical viewpoint we show that the method works well up to dimension 10. From the theoretical viewpoint the method offers many challenging problems both to statistics and D-module theory.