mpoly: Multivariate Polynomials in R

The mpoly package is a general purpose collection of tools for symbolic computing with multivariate polynomials in R. In addition to basic arithmetic, mpoly can take derivatives of polyno- mials, compute Grobner bases of collections of polynomials, and convert polynomials into a functional form to be evaluated. Among other things, it is hoped that mpoly will provide an R-based foundation for the computational needs of algebraic statisticians.


Introduction
At a basic level, R is not designed for symbolic computing.Unlike computer algebra systems founded on rational arithmetic and representations like Mathematica (Wolfram Research, 2012) and Maple (Maplesoft, 2012), R's dependence on floating-point arithmetic naturally lends itself to its purposedata analysis, numerical optimization, large(ish) data and so on.Nevertheless, the advent of algebraic statistics and its contributions to asymptotic theory in statistical models, experimental design, multiway contingency tables, and disclosure limitation has increased the need for R -as a practical tool for applied statistics -to be able to do some relatively basic operations and routines with multivariate polynomials (Diaconis and Sturmfels, 1998;Drton et al., 2009).mpoly is designed to try to fill this void entirely within R (Kahle, 2012).This article is intended to demonstrate its capabilities to a discipline/specialty-neutral audience.
A bit more specifically, the mpoly package is a general purpose collection of tools for symbolic computing with multivariate polynomials in R. While new, it is not the first package proposed to deal with multivariate polynomials.We therefore begin with a discussion of the current package intended to fulfill this need, multipol, in order to motivate the need for and capabilities of mpoly (Hankin, 2010).
Background: the multipol package R currently has three packages which deal with polynomialspolynom, PolynomF, and multipol.The first two (the second being a reboot of the first) deal exclusively with univariate polynomials and as such are not suitable for the general purpose needs cited in the introduction (Venables et al., 2009;Venables, 2010).The third, multipol, implements multivariate polynomials in a natural way suitable for some traditional applications (e.g.simple interactions in generalized linear models) but exhibits two limitations which impede its more general use and suggest a fresh re-envisioning of multivariate polynomials in R is needed (Hankin, 2008).The first limitation has to do with the impossibility of polynomial arithmetic, and the second has to do with storing sparse polynomials.
multipol is an implementation of multivariate polynomials based on the S3 class object "multipol", an unnamed multidimensional array.Consequently, multipol is fundamentally dependent on arrays as its basic data structure.This is evident in the construction of a "multipol" object: > library("multipol") Loading required package: abind In normal notation, a moment's consideration reveals that the polynomial defined above is 1 + 2x + 3y + 4xy + 5y 2 + 6xy 2 + 7z + 8xz + 9yz + 10xyz + 11y 2 z + 12xy 2 z in the ring R[x, y, z].(A similar pretty printing can be enabled by setting options(showchars = TRUE).)Thus, the elements of the array are the coefficients on the terms appearing in the polynomial.From here, the limitations are immediately accessible.
1. Notice that the definition of a multivariate polynomial in the example above makes no reference to the variable names x and y-they are determined when the print method is called.In other words, a "multipol" object has no notion of variables whatsoever, the variables are only generated when asked for via the print method for class "multipol".This has profound implications on polynomial arithmetic in multipol.In particular, arithmetic is only possible if two polynomials have the same number of dimensions (variables), and is only meaningful if those dimensions represent the same variables in the same order, the latter of which is never checked because there are no names specified which can be checked.
For example, suppose one wants to add or multiply the polynomials > options(showchars = TRUE) > ( a <-multipol( as.array( c("x^ " = , "x^1" = 1) ))) [1] 1*x^1 > ( b <-multipol( as.array( c("y^ " = , "y^1" = 1) ))) The result would seem to indicate that a + b = x + y = 2x, which is clearly in error.The location of the problem can easily be found however since we printed the polynomials at each step-regardless of our attempt to force (in a simple way) the labeling of the variables, we achieve the same incorrect result.
2. The array representation does not allow for sparse polynomials.For our purposes, a sparse (multivariate) polynomial of multidegree d is one which has "few" terms c d x d with multidegree d ≤ d (element-wise) with nonzero coefficients.As an example, consider the polynomial (3) Since multipol represents multivariate polynomials with arrays, the representation requires a 26 dimensional array (a dimension for each variable) each with three levels (e.g.a^ , a^1, a^2).This amounts to an array with 3 26 = 2, 541, 865, 828, 329 cells, all but 26 of which are nonzero-a whopping 20.33TB of space if the coefficients are stored in a double-precision floating point format.
This section has introduced two shortcomings of multipol which significantly limit its potential use for most practical applications.It should be noted, however, that multipol was not intended for use with "large" polynomials.It was built with enumerative combinatorics in mind, and its inefficiency with larger polynomials was readily acknowledged in Hankin (2008), where it is suggested that perhaps a sparse array representation (as is available in Mathematica) or outsourcing to C++ routines The R Journal Vol.5/1, June ISSN 2073-4859 could alleviate the burdens of the array representation.Its inefficiency for different purposes should therefore not be considered a flaw so much as outside the package's scope.Nevertheless, the arithmetic and storage problems caused by the array representation are fatal limitations for any problem which deals with large polynomial rings (i.e., ones with many variables).Enter mpoly.
mpoly -the basics mpoly is a complete reenvisioning of how multivariate polynomials and symbolic computing with multivariate polynomials should be implemented in R. Unlike multipol, mpoly uses as its most basic data structure the list.This fundamental change allows us to dispense of both issues with multipol discussed in the previous subsection.
In mpoly, an "mpoly" object is an S3 class object which is a list.The elements of the list are each named numeric vectors, with unique names including coef.Naturally, each element of an "mpoly" object corresponds to a term in the polynomial.The names of the elements correspond to the variables, and the values to the degrees; coef is the coefficient of the term.Constructing an "mpoly" object in this way is very straightforward using the constructor mpoly.> library("mpoly") Unlike multivariate polynomials in multipol, those in mpoly not only have variable names but also an intrinsic variable order which is taken to be the unique names of elements of the "mpoly" object minus coef.1This can be seen with the vars function.

> vars(poly) [1] "x" "y"
Viewing a multivariate polynomial as a list is a cumbersome task.To make things easier, a print method for "mpoly" objects exists and is dispatched when the object is queried by itself.

> poly 1 + 2 x^1 + 3 x^2 + 4 y^5 + 5 x y
One of the important considerations in polynomial algebra is the ordering of the terms of a multivariate polynomial.Notice the order of the terms presented in the printed version of the "mpoly" object; it is the order in which the terms are coded in the "mpoly" object itself.This can be changed in either of two ways.
> poly 1 + 2 x^1 + 3 x^2 + 4 y^5 + 5 x y > (poly2 <-reorder(poly, order = "lex")) using variable ordering -x, y Defining polynomials: mpoly vs mp The major workhorse of the package is the constructor mpoly itself.In particular, polynomial reduction (combining of like terms) and regularization (combining of coefficients and like variables within terms) are both performed when the multivariate polynomials are constructed with mpoly.
While the constructor mpoly is nice, it is inconvenient to have to keep specifying lists in order to define polynomials.The mp function was designed for this purpose and is intended to make defining multivariate polynomials quick and easy by taking them in as character strings and parsing them into "mpoly" objects.
> ( p <-mp("1 x + 2 y 3 + x^2 5 y") ) This parsing is a nontrivial process and depends heavily on the specification of the polynomial in the string.The mp function must first determine the variables that the user is specifying (which must be separated by spaces for disambiguation) and then construct the list to send to mpoly to construct the object.Because it is passed through mpoly, the "mpoly" object returned by mp is reduced and regularized.
> a^2 1 + 4 x + 6 y + 2 x y + 4 x^2 + 16 x^2 y + 9 y^2 + 24 x y^2 Here also each of the results are reduced and regularized.While the computations are done entirely in R, they are quite efficient; each of the above calculations is virtually instantaneous.
But mpoly does not stop there.The basic operations of the differential calculus, partial differentiation and gradients, are also available to the user and are efficient.A deriv method exists for "mpoly" objects which can be dispatched,3 and the gradient function is built on deriv to compute gradients.
> deriv(b, "x") 8 y z + 4 y + 6 z + 2 > gradient(b) 8 y z + 4 y + 6 z + 2 8 x z + 4 x + 7 z + 3 8 x y + 6 x + 7 y + 5 The gradient is a good example of another class object in the mpoly package, the "mpolyList"."mpolyList" objects are simply lists of "mpoly" objects and are used to hold vectors of multivariate polynomials.They can be easily specified using the mp function on a vector of character strings.
In algebraic geometry, a Gröbner basis of a polynomial ideal (or collection of polynomials generating an ideal) is a collection of polynomials which generates the same ideal and exhibits other nice properties, such as zero remainders for polynomials in the ideal and unique remainders for those outside it.They are foundational to the study of polynomial ideals and their associated varieties.In addition to differentiation, mpoly can also compute Gröbner bases of collections of multivariate polynomials ("mpolyList") by passing the proper objects in the proper syntax to the rSymPy package which has an implementation of Buchberger's algorithm (Grothendieck and Bellosta, 2010). 4The computations are performed by a Java based Python implementation and are quite fast once the Java Virtual Machine (JVM) has been initialized.The Gröbner basis is then returned as an "mpolyList" object.
> polys <-mp(c("x^2 -2 y^2", "x y -3")) > grobner(polys, varorder = c("x", "y")) 3 x -2 y^3 -9 + 2 y^4 > grobner(polys, varorder = c("x", "y"), order = "grlex") -3 x + 2 y^3 x^2 -2 y^2 -3 + x y The task of computing Gröbner bases is the only job mpoly outsources from R. This is because implementations of Gröbner bases algorithms are (1) complex and highly prone to error, (2) scarce and therefore difficult to emulate, and (3) highly iterative; thus it was thought best to leave the routine to more expert programmers and frameworks outside R. Unfortunately, there is currently no Gröbner walk algorithm available to convert a Gröbner basis in one monomial order to a Gröbner basis in another, a technique often used to quickly compute Gröbner bases in more difficult orders (e.g.lexicographic) from "easier" ones (e.g.graded reverse lexicographic), so there are still a number of improvements which can be made.

Evaluating polynomials
Apart from being interesting algebraic objects, polynomials can of course be thought of as functions.
> mpoly <-mp("x + 3 x y + z^2 x") > > f <-as.function(mpoly)f(.) with .= (x, y, z) > f(1:3) "mpolyList" objects can be converted into functions in the same way.Note that in both cases the user is messaged with the appropriate syntax for the resulting function.

Conclusion and future directions
While other functionality is also provided (such as a terms method, a function to compute least common multiples, and some combinatorial functions), the presented ones are the main contributions of the mpoly package.Put together, they provide a package which is not only user-friendly, efficient, and useful for polynomial algebra, but also a solid foundation upon which further developments can be made to make polynomials more accessible in R. In particular, it provides a nice starting point for any future package dealing with algebraic statistics.
There are, however, several improvements which can be made.First, the primary constructor function mpoly can be made significantly faster by outsourcing its job to C++, as previously suggested (Hankin, 2008).Second, improvements can be made to speed up some of the arithmetic, for example addition-chain exponentiation for polynomial exponents (as opposed to iterative multiplication).Last, improvements can be made into the flexibility of the facilitated constructor mp.In particular, parenthetical expressions are not currently handled but would be of great practical use.

Figure 2 :
Figure 2: The Rosenbrock function plotted by evaluating an "mpoly" object via as.function and apply with vector arguments.