The ‘jack’ package: Jack polynomials

R-CMD-check

library(jack)
library(microbenchmark)

Schur polynomials have applications in combinatorics and zonal polynomials have applications in multivariate statistics. They are particular cases of Jack polynomials. This package allows to evaluate these polynomials. It can also compute their symbolic form.

The functions JackPol, ZonalPol, ZonalQPol and SchurPol respectively return the Jack polynomial, the zonal polynomial, the quaternionic zonal polynomial, and the Schur polynomial.

Each of these polynomials corresponds is given by a positive integer, the number of variables, and an integer partition, the lambda argument; the Jack polynomial has one more parameter, the alpha argument, a positive number.

To get an exact symbolic polynomial with JackPol, you have to supply a bigq rational number for the parameter alpha:

jpol <- JackPol(2, lambda = c(3, 1), alpha = gmp::as.bigq("2/5"))
jpol
## 98/25*x^(3, 1) + 98/25*x^(1, 3) + 28/5*x^(2, 2)

This is a qspray object, from the qspray package. Here is how you can evaluate this polynomial:

qspray::evalQspray(jpol, c("2", "3/2"))
## Big Rational ('bigq') :
## [1] 1239/10

By default, ZonalPol, ZonalQPol and SchurPol return exact symbolic polynomials.

zpol <- ZonalPol(2, lambda = c(3, 1))
zpol
## 24/7*x^(3, 1) + 24/7*x^(1, 3) + 16/7*x^(2, 2)

It is also possible to convert a qspray polynomial to a function whose evaluation is performed by the Ryacas package:

zyacas <- as.function(zpol)

You can provide the values of the variables of this function as numbers or character strings:

zyacas(2, "3/2")
## [1] "594/7"

You can even pass a variable name to this function:

zyacas("x", "x")
## [1] "(64*x^4)/7"

If you want to substitute a variable with a complex number, use a character string which represents this number, with I denoting the imaginary unit:

zyacas("2 + 2*I", "2/3")
## [1] "Complex((-2176)/63,2944/63)"

Jack polynomials with Julia

As of version 2.0.0, the Jack polynomials can be calculated with Julia. The speed is amazing:

julia <- Jack_julia()
## Starting Julia ...
x <- c(1/2, 2/3, 1, 2/3, -1, -2, 1)
lambda <- c(5, 3, 2, 2, 1)
alpha <- 3
print(
  microbenchmark(
        R = Jack(x, lambda, alpha),
    Julia = julia$Jack(x, lambda, alpha),
    times = 6L,
    unit  = "seconds"
  ),
  signif = 6L
)
## Unit: seconds
##   expr        min        lq     mean    median        uq     max neval cld
##      R 7.73008000 7.8626500 7.903070 7.9140600 7.9534400 8.04412     6   b
##  Julia 0.00338613 0.0459389 0.216757 0.0697467 0.0948473 1.01688     6  a

Jack_julia() returns a list of functions. ZonalPol, ZonalQPol and SchurPol always return an exact expression of the polynomial, i.e. with rational coefficients (integers for SchurPol). If you want an exact expression with JackPol, you have to give a rational number for the argument alpha, as a character string:

JP <- julia$JackPol(m = 2, lambda = c(3, 1), alpha = "2/5")
JP
## 98/25*x^(1, 3) + 28/5*x^(2, 2) + 98/25*x^(3, 1)

Again, Julia is faster:

n <- 5
lambda <- c(4, 3, 3)
alpha <- "2/3"
alphaq <- gmp::as.bigq(alpha)
microbenchmark(
      R = JackPol(n, lambda, alphaq),
  Julia = julia$JackPol(n, lambda, alpha),
  times = 6L
)
## Unit: milliseconds
##   expr       min        lq      mean    median        uq       max neval cld
##      R 1170.8920 1188.7068 1225.4087 1211.0391 1236.0073 1334.7680     6   b
##  Julia  464.3045  491.7447  549.3006  496.9847  546.5387  799.2461     6  a
JuliaConnectoR::stopJulia()