The `mvp` package: fast multivariate polynomials in R

In mvp: Fast Symbolic Multivariate Polynomials

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)

\hfill{width=10%}

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]. I present some timing results separately, in inst/timings.Rmd. 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. An mvp object is a map from terms to their coefficients; 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-discernible 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 substitute for the other variables:

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

(in this case, the default behaviour is to return the the resulting polynomial coerced to a scalar). We can suppress the coercion using the lose argument:

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

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")

Vectorised substitution

Function subvec()allows one to substitute variables for numeric values using vectorised idiom:

p <- rmvp(6,2,2,letters[1:3])
p
subvec(p,a=1,b=2,c=1:5)   # supply a named list of vectors

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)
trunc(X,3)         # truncate, retain only terms with total power <= 3
trunc1(X,x=3)    # truncate, retain only terms with  power of x <= 3
onevarpow(X,x=3) # retain only 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

Function series() will decompose an mvp object into a power series in a single variable:

p <- as.mvp("a^2 x b + x^2 a b + b c x^2 + a b c + c^6 x")
p
series(p,'x')

This works nicely with subs() if we wish to take a power series about x-v, where v is any mvp object. For example:

p %>% subs(x="xmv+a+b") %>% series("xmv") 

is a series in powers of x-a-b. We may perform a consistency check by a second substitution, returning us to the original expression:

p == p %>% subs(x="xmv+a+b") %>% subs(xmv="x-a-b")

If function series() is given a variable name ending in _m_foo, where foo is any variable name, then this is typeset as (x-foo). For example:

as.mvp('x^3 + x*a') %>% subs(x="x_m_a + a") %>% series("x_m_a")

So above we see the expansion of $x^2+ax$ in powers of $x-a$. If we want to see the expansion of a mvp in terms of a more complicated expression then it is better to use a nonce variable v:

as.mvp('x^2 + x*a+b^3') %>% subs(x="x_m_v + a^2+b") %>% series("x_m_v")

where it is understood that $v=a+b^2$. Function taylor() is a convenience wrapper that does some of the above in one step:

p <- as.mvp("1+x-x*y+a")^2
taylor(p,'x','a')

But it's not as good as I expected it to be and frankly it's overkill.

Extraction

Given a multivariate polynomial, one often needs to extract certain terms. Because the terms of an mvp object have an implementation-dependent order, this can be difficult. But we can use function onevarpow():

P <- as.mvp("1 + z + y^2 + x*z^2 +  x*y")^4
onevarpow(P,x=1,y=2)

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")

The disordR package

It is possible to examine the coefficients of an mvp object:

a <- as.mvp("5 + 8*x^2*y - 13*y*x^2 + 11*z - 3*x*yz")
a
coeffs(a)

Above, note that the result of coeffs() is a disord object, defined in the disordR package. The order of the elements unspecified as the STL map class holds the keys and values in an implementation-specific order. This device stops the user from illegal operations on the coefficients. For example, suppose we had another mvp object, b:

b <- a*2
b
coeffs(a) + coeffs(b)

above, we get an error because the coefficients of a and b are possibly stored in a different order and therefore vector addition makes no sense. However, we can operate on coefficients of a single mvp object at will:

coeffs(a) > 0
coeffs(a) + coeffs(a)^4

Extraction also works but subject to standard disordR idiom restrictions:

coeffs(a)[coeffs(a) > 0]

But "mixing" objects is forbidden:

coeffs(a)[coeffs(b) > 0]

Extraction methods work, again subject to disordR restrictions:

coeffs(a)[coeffs(a)<0] <- coeffs(a)[coeffs(a)<0] + 1000 # add 1000 to every negative coefficient
a

In cases like this where the replacement object is complicated, using magrittr would simplify the idiom and reduce the opportunity for error:

library("magrittr")
b
coeffs(b)[coeffs(b)%%3==1] %<>% `+`(100)  # add 100 to every element equal to 1 modulo 3
b

One good use for this is to "zap" small elements:

x <- as.mvp("1 - 0.11*x + 0.005*x*y")^2
x

Then we can zap as follows:

cx <- coeffs(x)
cx[abs(cx) < 0.01] <- 0
coeffs(x) <- cx
x

(I should write a method for zapsmall() that does this)

Multivariate generating functions

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)

References

mvp documentation built on March 31, 2023, 5:43 p.m.