The ```mvp``` package: fast multivariate polynomials R

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 3.5, message = FALSE, warning = FALSE)
options(width = 80, tibble.width = Inf)

Introduction

The mvp package provides some functionality for fast manipulation of multivariate polynomials, using the Standard Template library of C++, commonly known as the STL. It is comparable in speed to the spray package for sparse arrays, while retaining the symbolic capabilities of the mpoly package [@kahle2013]. The mvp package uses the excellent print and coercion methods of mpoly. The mvp package provides improved speed over mpoly, the ability to handle negative powers, and a more sophisticated substitution mechanism.

The STL map class

A map is a sorted associative container that contains key-value pairs with unique keys. It is interesting here because search and insertion operations have logarithmic complexity. Multivariate polynomials are considered to be the sum of a finite number of terms, each multiplied by a coefficient. A term is something like $x^2y^3z$. We may consider this term to be the map

{"x" -> 2, "y" -> 3, "z" -> 1}

where the map takes symbols to their (integer) power; it is understood that powers are nonzero. A mvp object is a map from terms to their coeffients; thus $7xy^2 -3x^2yz^5$ would be

{{"x" -> 1, "y" -> 2} -> 7, {"x" -> 2, 'y" -> 1, "z" ->5} -> -3}

and we understand that coefficients are nonzero. In C++ the declarations would be

typedef vector <signed int> mypowers;  
typedef vector <string> mynames;  

typedef map <string, signed int> term; 
typedef map <term, double> mvp; 

Thus a term maps a string to a (signed) integer, and a mvp maps terms to doubles. One reason why the map class is fast is that the order in which the keys are stored is undefined: the compiler may store them in the order which it regards as most propitious. This is not an issue for the maps considered here as addition and multiplication are commutative and associative.

Note also that constant terms are handled with no difficulty (constants are simply maps from the empty map to its value), as is the zero polynomial (which is simply an empty map).

The package in use

Consider a simple multivariate polynomial $3xy+z^3+xy^6z$ and its representation in the following R session:

library("mvp",quietly=TRUE)
(p <- as.mvp("3 x y + z^3 + x y^6 z"))

Coercion and printing are accomplished by the mpoly package (there is no way I could improve upon Kahle's work). Note carefully that the printed representation of the mvp object is created by the mpoly package and the print method can rearrange both the terms of the polynomial ($3xy+z^3+xy^6z = z^3+3xy+xy^6z$, for example) and the symbols within a term ($3xy=3yx$, for example) to display the polynomial in a human-friendly form.

However, note carefully that such rearranging does not affect the mathematical properties of the polynomial itself. In the mvp package, the order of the terms is not preserved (or even defined) in the internal representation of the object; and neither is the order of the symbols within a single term. Although this might sound odd, if we consider a marginally more involved situation, such as

(M <- as.mvp("3 stoat goat^6 -4 + 7 stoatboat^3 bloat -9 float boat goat gloat^6"))
dput(M)

it is not clear that any human-discernable ordering is preferable to any other, and we would be better off letting the compiler decide a propitious ordering. In any event, the mpoly package can specify a print order:

print(M,order="lex", varorder=c("stoat","goat","boat","bloat","gloat","float","stoatboat"))

Arithmetic operations

The arithmetic operations *, +, - and ^ work as expected:

(S1 <- rmvp(5,2,2,4))
(S2 <- rmvp(5,2,2,4))
S1 + S2
S1 * S2
S1^2

Substitution

The package has two substitution functionalities. Firstly, we can substitute one or more variables for a numeric value. Define a mvp object:

(S3 <- as.mvp("x + 5 x^4 y + 8 y^2 x z^3"))

And then we may substitute $x=1$:

subs(S3, x = 1)

Note the natural R idiom, and that the return value is another mvp object. We may subsitute for the other variables:

subs(S3, x = 1, y = 2, z = 3)

(in this case, the default behaviour is to return a "dropped" version of the resulting polynomial, that is, coerced to a scalar). The idiom also allows one to substitute a variable for an mvp object:

subs(as.mvp("a+b+c"), a="x^6")

Note carefully that subs() depends on the order of substitution:

subs(as.mvp("a+b+c"), a="x^6",x="1+a")
subs(as.mvp("a+b+c"), x="1+a",a="x^6")

Pipes

Substitution works well with pipes:

as.mvp("a+b") %>% subs(a="a^2+b^2") %>% subs(b="x^6")

Differentiation

Differentiation is implemented. First we have the deriv() method:

(S <- as.mvp("a + 5 a^5*b^2*c^8 -3 x^2 a^3 b c^3"))
deriv(S, letters[1:3])
deriv(S, rev(letters[1:3]))  # should be the same.

Also a slightly different form: aderiv(), here used to evaluate $\frac{\partial^6S}{\partial a^3\partial b\partial c^2}$:

aderiv(S, a = 3, b = 1, c = 2)

Again, pipes work quite nicely:

S %<>% aderiv(a=1,b=2) %>% subs(c="x^4") %>% `+`(as.mvp("o^99"))
S

Taylor series

The package includes functionality to deal with Taylor and Laurent series:

(X <- as.mvp("1+x+x^2 y")^3)
taylor(X,3)        # terms with total power <= 3
taylor(X,3,"x")    # terms with power of x <= 3
onevarpow(X,3,"x") # terms with power of x == 3
## second order taylor expansion of f(x)=sin(x+y) for x=1.1, about x=1:
sinxpy <- horner("x+y",c(0,1,0,-1/6,0,+1/120,0,-1/5040))  # sin(x+y)
dx <- as.mvp("dx")
t2 <- sinxpy  + aderiv(sinxpy,x=1)*dx + aderiv(sinxpy,x=2)*dx^2/2
(t2 %<>% subs(x=1,dx=0.1))  # (Taylor expansion of sin(y+1.1), left in symbolic form)
(t2 %>% subs(y=0.3))  - sin(1.4)  # numeric; should be small

Negative powers

The mvp package handles negative powers, although the idiom is not perfect and I'm still working on it. There is the invert() function:

(p <- as.mvp("1+x+x^2 y"))
invert(p)

In the above, p is a regular multivariate polynomial which includes negative powers. It obeys the same arithmetic rules as other mvp objects:

p + as.mvp("z^6")

We can see the generating function for a chess knight:

knight(2)

How many ways are there for a 4D knight to return to its starting square after four moves? Answer:

constant(knight(4)^4)

Some timings

I will show some timings using a particularly favourable example that exploits the symbolic nature of the mvp package.

library("spray")
library("mpoly")
n <- 100
k <- kahle(n, r = 3, p = 1:3, symbols = paste0("x", sprintf("%03d", 1:n)))

In the above, polynomial k has 500 terms of the form $xy^2z^3$. Coercing k to spray form would need a $500\times 500$ matrix for the indices, almost every element of which would be zero. This makes the spray package slower:

library("microbenchmark")

spray_k <- mvp_to_spray(k)
microbenchmark( k^2, spray_k^2 )

In the above, the first line uses mvp functionality, and the second line uses spray functionality. The speedup increases for larger polynomials.

References



Try the mvp package in your browser

Any scripts or data that you put into this service are public.

mvp documentation built on Feb. 4, 2019, 5:05 p.m.