[Paper Review] Holonomic gradient method for the distribution function of the largest root of a Wishart matrix

Hiroki Hashiguchi, Yasuhide Numata+2 more|arXiv (Cornell University)|Jan 2, 2012

Random Matrices and Applications22 references4 citations

TL;DR

This paper introduces the holonomic gradient method (HGM) to compute the cumulative distribution function of the largest eigenvalue of a Wishart matrix, leveraging partial differential equations (PDEs) satisfied by the confluent hypergeometric function of matrix argument ${}_1F_1$. The method enables accurate numerical evaluation beyond the convergence limitations of series expansions at the origin, achieving high accuracy for dimensions up to 10 with controlled discretization errors.

ABSTRACT

We apply the holonomic gradient method introduced by Nakayama et al.(2011) 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 numerical viewpoint we show that the method works well up to dimension 10. From theoretical viewpoint the method offers many challenging problems both to statistics and D-module theory.

Motivation & Objective

• To address the longstanding computational challenge in multivariate statistics of evaluating the cumulative distribution function of the largest root of a Wishart matrix.
• To overcome the slow convergence of infinite series expansions in zonal polynomials, especially for large argument values.
• To develop a numerically robust method that moves away from the origin by solving PDEs satisfied by the hypergeometric function ${}_1F_1$.
• To demonstrate the feasibility and accuracy of the holonomic gradient method for dimensions up to 10.
• To establish a bridge between statistical multivariate analysis and $D$-module theory through a novel computational approach.

Proposed method

• The holonomic gradient method is applied to the hypergeometric function ${}_1F_1$ of a matrix argument, which arises in the distribution of the largest eigenvalue of a Wishart matrix.
• The method uses a Pfaffian system—a system of first-order linear PDEs—derived from the holonomic system satisfied by ${}_1F_1$, enabling numerical integration from initial values.
• Initial values are computed via series expansion in zonal polynomials around the origin, providing the starting point for numerical integration.
• The PDE system is solved numerically using an ODE integrator (e.g., from the deSolve package in R), advancing the solution along a path from the origin to the desired argument value.
• The method is implemented in R for dimension $m=2$, with symbolic computation used to derive the PDE coefficients and verify holonomicity via Gröbner basis techniques.
• The solution is normalized using gamma functions and exponential terms to yield the final cumulative distribution function.

Experimental results

Research questions

• RQ1 Can the holonomic gradient method provide a numerically stable and accurate alternative to series expansions for computing the distribution of the largest eigenvalue of a Wishart matrix?
• RQ2 How well does the HGM perform in terms of accuracy and convergence speed compared to infinite series expansions, especially for large argument values?
• RQ3 What is the maximum dimension for which the HGM remains computationally feasible and accurate?
• RQ4 Can the holonomic system approach be systematically extended to higher-dimensional Wishart distributions using $D$-module theory?
• RQ5 What are the theoretical and computational challenges in constructing and solving the Pfaffian system for general $m$?

Key findings

• The holonomic gradient method achieves high accuracy in computing the cumulative distribution function of the largest eigenvalue of a Wishart matrix for dimensions up to 10.
• The method provides exact results with errors only from numerical discretization and initial value computation, unlike series expansions that suffer from slow convergence at large arguments.
• For dimension $m=2$, the method is implemented in R using the deSolve package, and numerical results confirm consistency with known analytical forms.
• The underlying system of PDEs is shown to be holonomic via Gröbner basis computation, confirming the existence of a finite-dimensional solution space.
• The method is complementary to existing approaches: it excels for small degrees of freedom where Laplace approximations fail, while series expansions are effective near the origin.
• Theoretical analysis confirms that the left ideal generated by the differential operators is holonomic, ensuring the existence of a well-defined, finite-dimensional solution space.