Package ‘hgm’
March 27, 2014
Type Package
Depends deSolve
Title Holonomic gradient method and gradient descent
Version 1.6Date 2014-03-25
Author Nobuki Takayama, Tamio Koyama, Tomonari Sei
Maintainer Nobuki Takayama <takayama@math.kobe-u.ac.jp>
Description The holonomic gradient method (HGM, hgm) gives a way to evaluate normalization
constants of unnormalized probability distributions by utilizing holonomic
systems of differential equations. The holonomic gradient descent (HGD, hgd) gives a methodto find maximal likelihood estimates by utilizing the HGM.
License GPL-2
LazyLoad yes
URL http://www.openxm.org

R topics documented:
hgm-package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
hgm.ncorthant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
hgm.ncso3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
hgm.pwishart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
hgm.Rhgm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Index

23
4
5
7
91

2

hgm-package

hgm-package

HGM

Description
The holonomic gradient method (HGM, hgm) gives a way to evaluate normalization constants ofunnormalized probability distributions by utilizing holonomic systems of differential equations.
The holonomic gradient descent (HGD, hgd) gives a method to find maximal likelihood estimatesby utilizing the HGM.
Details
Package:Type:
License:LazyLoad:hgm
PackageGPL-2
yesThe HGM and HGD are proposed in the paper below. This method based on the fact that a broadclass of normalization constants of unnormalized probability distributions belongs to the class ofholonomic functions, which are solutions of holonomic systems of linear partial differential equations.
NoteThis package includes a small subset of the Gnu scientific library codes (http://www.gnu.org/software/gsl/). Then, it might cause a conflict with the package gsl (see gsl-package).
References• Hiromasa Nakayama, Kenta Nishiyama, Masayuki Noro, Katsuyoshi Ohara, Tomonari Sei,
Nobuki Takayama, Akimichi Takemura, Holonomic Gradient Descent and its Application toFisher-Bingham Integral, Advances in Applied Mathematics 47 (2011), 639–658, http://dx.doi.org/10.1016/j.aam.2011.03.001• http://www.openxm.org
See Also
hgm.ncorthant hgm.ncso3 hgm.pwishart hgm.Rhgm
Examples## Not run:
example(hgm.ncorthant)example(hgm.ncso3)
example(hgm.pwishart)example(hgm.Rhgm)

hgm.ncorthant3## End(Not run)hgm.ncorthant

The function hgm.ncorthant evaluates the orthant probability.

Description
The function hgm.ncorthant evaluates the orthant probability, which is the normalization constantof the multivariate normal distribution restrcted to the first orthant.
Usage
hgm.ncorthant(x,y,rk_step_size=1e-3)
Arguments
xyThis function evaluates the orthant probability for the m dimensional multivariate normal distribution whose m by m covariance matrix and the mean vector ofsize m are x and y respectively.

rk_step_size

The step size for the Runge-Kutta method to apply the HGM.

Details
The function hgm.ncorthant evaluates the orthant probability, which is the normalization constantof the m-dimensional multivariate normal distribution restrcted to the first orthant. It uses the holonomic gradient method (HGM) to evalute it. The rank of the system of differential equations for theHGM is 2^m.
Value
The output is the orthant probalibity.
Author(s)Tamio Koyama
References
Tamio Koyama, Akimichi Takemura, Calculation of orthant probabilities by the holonomic gradientmethod, http://arxiv.org/abs/1211.682.

4

hgm.ncso3

Examples## =====================================================## Example 1. Computing the orthant probability## =====================================================
x<-matrix(c(1, 0.5, 0.5, 1),nrow=2)y<-c(1.0, 0.5)
hgm.ncorthant(x,y)

The function hgm.ncso3 evaluates the normalization constant for theFisher distribution on SO(3).

hgm.ncso3

Description
The function hgm.ncso3 evaluates the normalization constant for the Fisher distribution on SO(3).
Usage
hgm.ncso3(a,b,c,t0=0.0,q=1,deg=0)
Arguments
abc

This function evaluates the normalization constant for the parameter Theta=diag(theta_ii)
of the Fisher distribution on SO(3). The variables a,b,c stand for the parameterstheta_11, theta_22, theta_33 respectively.

t0

It is the initial point to evaluate the series. If it is set to 0.0, a default value isused.

q

If it is 1, then the program works in a quiet mode.

degIt gives the approximation degree of the power series approximation of the normalization constant near the origin. If it is 0, a default value is used.

Details
The normalization constant c(Theta) of the Fisher distribution on SO(3) is defined by integral(exp(trace( transpose(Theta) X)) ) where X is the integration variable and runs over S0(3) and Theta
is a 3 x 3 matrix parameter. A general HGM algorithm to evaluate the normalization constantis given in the reference below. We use the Corollary 1 and the series expansion in 3.2 for theevaluation.
ValueThe output is c(Theta).

hgm.pwishart

5Author(s)Nobuki Takayama
References
Tomonari Sei, Hiroki Shibata, Akimichi Takemura, Katsuyoshi Ohara, Nobuki Takayama, Properties and applications of Fisher distribution on the rotation group, Journal of Multivariate Analysis,
116 (2013), 440–455, http://dx.doi.org/10.1016/j.jmva.2013.01.010
Examples## =====================================================## Example 1. Computing normalization constant of the Fisher distribution on SO(3)## =====================================================
hgm.ncso3(1,2,3)## =====================================================## Example 2. Asteroid data in the paper## =====================================================
hgm.ncso3(19.6,0.831,-0.671)

hgm.pwishart

The function hgm.pwishart evaluates the cumulative distribution function of random wishart matrix.

Description
The function hgm.pwishart evaluates the cumulative distribution function of random wishart matrixof size m times m.
Usage
hgm.pwishart(m,n,beta,q0,approxdeg,h,dp,q,mode,method,err)
Argumentsm

The dimension of the Wishart matrix.

nThe degree of freedome (a parameter of the Wishart distribution)beta

The eigenvalues of the inverse of the covariant matrix /2 (a parameter of theWishart distribution). The beta is equal to inverse(sigma)/2.

q0

The point to evaluate the matrix hypergeometric series. q0>0

approxdeg

Zonal polynomials upto the approxdeg are calculated to evaluate values near theorigin. A zonal polynomial is determined by a given partition (k1,...,km). We
call the sum k1+...+km the degree.

hA (small) step size for the Runge-Kutta method. h>0.

dp

Sampling interval of solutions by the Runge-Kutta method.

6hgm.pwishart
qThe second value y[0] of this function is the Prob(L1 < q) where L1 is the first
eigenvalue of the Wishart matrix.

mode

When mode=c(1,0,0), it returns the evaluation of the matrix hypergeometric series and its derivatives at x0. When mode=c(1,1,(m^2+1)*p), intermediate values of P(L1 < x) with respect to p-steps of x are also returned. Sampling interval
is controled by dp.

methoda-rk4 is the default value. When method="a-rk4", the adaptive Runge-Kutta
method is used. Steps are automatically adjusted by err.

err

When err=c(e1,e2), e1 is the absolute error and e2 is the relative error. As longas NaN is not returned, it is recommended to set to err=c(0.0, 1e-10), because
initial values are usually very small.

automatic

automatic=1 is the default value. If it is 1, the degree of the series approximation will be increased until |(F(i)-F(i-1))/F(i-1)| < assigned_series_error whereF(i) is the degree i approximation of the hypergeometric series with matrix argument. Step sizes for the Runge-Kutta method are also set automatically fromthe assigned_series_error if it is 1.

assigned_series_error
assigned_series_error=0.00001 is the default value.
verboseverbose=0 is the default value. If it is 1, then steps of automatic degree updates
and several parameters are output to stdout and stderr.

Details
It is evaluated by the Koev-Edelman algorithm when x is near the origin and by the HGM whenx is far from the origin. We can obtain more accurate result when the variables h is smaller, x0 isrelevant value (not very big, not very small), and the approxdeg is more larger. A heuristic method
to set parameters x0, h, approxdeg properly is to make x larger and to check if the y[0] approaches
to 1.
Value
The output is x, y[0], ..., y[2^m] in the default mode, y[0] is the value of the cumulative distributionfunction P(L1 < x) at x. y[1],...,y[2^m] are some derivatives. See the reference below.
Author(s)Nobuki Takayama
References
H.Hashiguchi, Y.Numata, N.Takayama, A.Takemura, Holonomic gradient method for the distribution function of the largest root of a Wishart matrix, Journal of Multivariate Analysis, 117, (2013)296-312, http://dx.doi.org/10.1016/j.jmva.2013.03.011,

hgm.Rhgm

7Examples## =====================================================## Example 1.
## =====================================================
hgm.pwishart(m=3,n=5,beta=c(1,2,3),q=10)## =====================================================## Example 2.
## =====================================================b<-hgm.pwishart(m=4,n=10,beta=c(1,2,3,4),q0=1,q=10,approxdeg=20,mode=c(1,1,(16+1)*100));c<-matrix(b,ncol=16+1,byrow=1);#plot(c)hgm.Rhgm

The function hgm.Rhgm performs the holonomic gradient method(HGM) for a given Pfaffian system and an initial value vector.

Description
The function hgm.Rhgm performs the holonomic gradient method (HGM) for a given Pfaffian system and an initial value vector with the deSolve package in R.
Usage
hgm.Rhgm(th0, G0, th1, dG.fun, times=NULL, fn.params=NULL)
Argumentsth0

A d-dimensional vector which is an initial point of the parameter vector th(theta).

G0

A r-dimensional vector which is the initial value of the vector G of the normalizing constant and its derivatives.

th1

A d-dimensional vector which is the target point of th.

dG.fundG.fun is the “right hand sides” of the Pfaffian system. It is a d*r-dimensional
array.

timesa vector; times in [0,1] at which explicit estimates for G are desired. If time =NULL, the set 0,1 is used, and only the final value is returned.

fn.params

fn.params: a list of parameters passed to the function dG.fun. If fn.params =NULL, no parameter is passed to dG.fun.

Details
The function hgm.Rhgm computes the value of a holonomic function at a given point, using HGM.
This is a “Step 3” function (see the reference below), which can be used for an arbitrary input, inthe HGM framework. Efficient “Step 3” functions are given for some distributions in this package.
The Pfaffian system assumed is d G_j / d th_i = (dG.fun(th, G))_i,j

8hgm.Rhgm
The inputs of hgm.Rhgm are the initial point th0, initial value G0, final point th1, and Pfaffiansystem dG.fun. The output is the final value G1.
If the argument ‘times’ is specified, the function returns a matrix, where the first column denotes
time, the following d-vector denotes th, and the remaining r-vector denotes G.

ValueThe output is the value of G at th1. The first element of G is the normalizing constant.
Author(s)Tomonari Sei
References
http://www.math.kobe-u.ac.jp/OpenXM/Math/hgm/ref-hgm.html
Examples# Example 1.
# A demo program; von Mises--Fisher on S^{3-1}G.exact = function(th){ # exact value by built-in functionc( sinh(th[1])/th[1], cosh(th[1])/th[1] - sinh(th[1])/th[1]^2 )}dG.fun = function(th, G, fn.params=NULL){ # PfaffiandG = array(0, c(1, 2))
sh = G[1] * th[1]ch = G[2] * th[1] + G[1]dG[1,1] = G[2] # Pfaffian eq’sdG[1,2] = sh/th[1] - 2*ch/th[1]^2 + 2*sh/th[1]^3dG}
th0 = 0.5
th1 = 15
G0 = G.exact(th0)G0G1 = hgm.Rhgm(th0, G0, th1, dG.fun)
G1
G1.exact = G.exact(th1)G1.exact## Example 2.
#hgm.Rhgm.demo1()# HGMIndex∗Topic Cumulative distributionfunction of random wishart
matrixhgm.pwishart, 5∗Topic Fisher distribution on SO(3)hgm.ncso3, 4∗Topic HGD
hgm-package, 2∗Topic HGM
hgm-package, 2
hgm.ncorthant, 3hgm.ncso3, 4hgm.pwishart, 5hgm.Rhgm, 7∗Topic Holonomic gradient method
hgm.ncorthant, 3
hgm.ncso3, 4hgm.pwishart, 5
hgm.Rhgm, 7∗Topic Normalization constant
hgm.ncorthant, 3
hgm.ncso3, 4
hgm.Rhgm, 7∗Topic Orthant probability
hgm.ncorthant, 3∗Topic holonomic gradient descenthgm-package, 2∗Topic holonomic gradient method
hgm-package, 2∗Topic package
hgm-package, 2
HGM (hgm-package), 2
hgm (hgm-package), 2
hgm-package, 2
hgm.ncorthant, 2, 3
hgm.ncso3, 2, 4
hgm.pwishart, 2, 5
hgm.Rhgm, 2, 7

9