mvp package: fastmultivariate polynomials in RTo cite themvp package in publications please useHankin (2022b). Multivariate polynomialsare interesting and useful objects, with a wide range of applications.Themvp package provides some functionality for fastmanipulation of multivariate polynomials, using the Standard Templatelibrary ofC++(Stroustrup 1997), commonly knownas theSTL. It is comparable in speed to thespray package for sparse arrays, while retaining thesymbolic capabilities of thempoly package(Kahle 2013). Ipresent some timing results separately, ininst/timings.Rmd. Themvp package uses theexcellent print and coercion methods ofmpoly. Themvp package provides improved speed overmpoly, the ability to handle negative powers, and a moresophisticated substitution mechanism.
STL map classAmap is a sorted associative container that containskey-value pairs with unique keys. It is interesting here because searchand insertion operations have logarithmic complexity. Multivariatepolynomials are considered to be the sum of a finite number ofterms, each multiplied by a coefficient. Aterm issomething like\(x^2y^3z\). We mayconsider this term to be the map
{"x" -> 2, "y" -> 3, "z" -> 1}where the map takes symbols to their (integer) power; it isunderstood that powers are nonzero. Anmvp object is a mapfrom 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. InC++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 aterm maps a string to a (signed) integer, and amvp maps terms to doubles. One reason why themap class is fast is that the order in which the keys arestored is undefined: the compiler may store them in the order which itregards as most propitious. This is not an issue for the maps consideredhere 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 thezero polynomial (which is simply an empty map).
Consider a simple multivariate polynomial\(3xy+z^3+xy^6z\) and its representation inthe following R session:
library("mvp",quietly=TRUE)library("magrittr")(p<-as.mvp("3 x y + z^3 + x y^6 z"))#> mvp object algebraically equal to#> 3 x y + x y^6 z + z^3Coercion and printing are accomplished by thempolypackage (there is no way I could improve upon Kahle’s work). Notecarefully that the printed representation of the mvp object is createdby thempoly package and the print method can rearrangeboth the terms of the polynomial (\(3xy+z^3+xy^6z = z^3+3xy+xy^6z\), forexample) and the symbols within a term (\(3xy=3yx\), for example) to display thepolynomial in a human-friendly form.
However, note carefully that such rearranging does not affect themathematical properties of the polynomial itself. In themvp package, the order of the terms is not preserved (oreven defined) in the internal representation of the object; and neitheris the order of the symbols within a single term. Although this mightsound odd, if we consider a marginally more involved situation, suchas
(M<-as.mvp("3 stoat goat^6 -4 + 7 stoatboat^3 bloat -9 float boat goat gloat^6"))#> mvp object algebraically equal to#> -4 + 7 bloat stoatboat^3 - 9 boat float gloat^6 goat + 3 goat^6 stoatdput(M)#> structure(list(names = list(character(0), c("bloat", "stoatboat"#> ), c("boat", "float", "gloat", "goat"), c("goat", "stoat")),#> power = list(integer(0), c(1L, 3L), c(1L, 1L, 6L, 1L), c(6L,#> 1L)), coeffs = c(-4, 7, -9, 3), "a42fbf7de80773710e0bf126ddb55cd6ece4a8dc"), class = "mvp")it is not clear that any human-discernible ordering is preferable toany other, and we would be better off letting the compiler decide apropitious ordering. In any event, thempoly package canspecify a print order:
print(M,order="lex",varorder=c("stoat","goat","boat","bloat","gloat","float","stoatboat"))#> mvp object algebraically equal to#> 3 stoat goat^6 - 9 goat boat gloat^6 float + 7 bloat stoatboat^3 - 4Note in passing that variable names may be any character string, notjust single letters.
The arithmetic operations*,+,- and^ work as expected:
(S1<-rmvp(5,2,2,4))#> mvp object algebraically equal to#> 3 + 5 a + 2 a b + 4 c d + c^2(S2<-rmvp(5,2,2,4))#> mvp object algebraically equal to#> 2 + 5 c + 9 dS1+ S2#> mvp object algebraically equal to#> 5 + 5 a + 2 a b + 5 c + 4 c d + c^2 + 9 dS1* S2#> mvp object algebraically equal to#> 6 + 10 a + 4 a b + 10 a b c + 18 a b d + 25 a c + 45 a d + 15 c + 8 c d + 36 c#> d^2 + 2 c^2 + 29 c^2 d + 5 c^3 + 27 dS1^2#> mvp object algebraically equal to#> 9 + 30 a + 12 a b + 16 a b c d + 4 a b c^2 + 40 a c d + 10 a c^2 + 25 a^2 + 20#> a^2 b + 4 a^2 b^2 + 24 c d + 6 c^2 + 16 c^2 d^2 + 8 c^3 d + c^4The package has two substitution functionalities. Firstly, we cansubstitute one or more variables for a numeric value. Define a mvpobject:
(S3<-as.mvp("x + 5 x^4 y + 8 y^2 x z^3"))#> mvp object algebraically equal to#> x + 8 x y^2 z^3 + 5 x^4 yAnd then we may substitute\(x=1\):
Note the natural R idiom, and that the return value is another mvpobject. We may substitute for the other variables:
(in this case, the default behaviour is to return the resultingpolynomial coerced to a scalar). We can suppress the coercion using thedrop (formerlylose) argument:
The idiom also allows one to substitute a variable for anmvp object:
Note carefully thatsubs() depends on the order ofsubstitution:
Functionsubvec()allows one to substitute variables fornumeric values using vectorised idiom:
Differentiation is implemented. First we have thederiv() method:
(S<-as.mvp("a + 5 a^5*b^2*c^8 -3 x^2 a^3 b c^3"))#> mvp object algebraically equal to#> a - 3 a^3 b c^3 x^2 + 5 a^5 b^2 c^8deriv(S, letters[1:3])#> mvp object algebraically equal to#> -27 a^2 c^2 x^2 + 400 a^4 b c^7deriv(S,rev(letters[1:3]))# should be the same.#> mvp object algebraically equal to#> -27 a^2 c^2 x^2 + 400 a^4 b c^7Also a slightly different form:aderiv(), here used toevaluate\(\frac{\partial^6S}{\partiala^3\partial b\partial c^2}\):
Again, pipes work quite nicely:
The package includes functionality to deal with Taylor and Laurentseries:
(X<-as.mvp("1+x+x^2 y")^3)#> mvp object algebraically equal to#> 1 + 3 x + 3 x^2 + 3 x^2 y + x^3 + 6 x^3 y + 3 x^4 y + 3 x^4 y^2 + 3 x^5 y^2 +#> x^6 y^3trunc(X,3)# truncate, retain only terms with total power <= 3#> mvp object algebraically equal to#> 1 + 3 x + 3 x^2 + 3 x^2 y + x^3trunc1(X,x=3)# truncate, retain only terms with power of x <= 3#> mvp object algebraically equal to#> 1 + 3 x + 3 x^2 + 3 x^2 y + x^3 + 6 x^3 yonevarpow(X,x=3)# retain only terms with power of x == 3#> mvp object algebraically equal to#> 1 + 6 y## 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)#> mvp object algebraically equal to#> 0.8912877 + 0.4534028 y - 0.4458333 y^2 - 0.07597222 y^3 + 0.03659722 y^4 +#> 0.003291667 y^5 - 0.001527778 y^6 - 0.0001984127 y^7(t2%>%subs(y=0.3))-sin(1.4)# numeric; should be small#> [1] -1.416914e-05Functionseries() will decompose anmvpobject 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#> mvp object algebraically equal to#> a b c + a b x^2 + a^2 b x + b c x^2 + c^6 xseries(p,'x')#> x^0(a b c) + x^1(a^2 b + c^6) + x^2(a b + b c)This works nicely withsubs() if we wish to take a powerseries aboutx-v, wherev is anymvp object. For example:
p%>%subs(x="xmv+a+b")%>%series("xmv")#> xmv^0(a b c + 2 a b^2 c + a b^3 + a c^6 + a^2 b c + 3 a^2 b^2 + 2#> a^3 b + b c^6 + b^3 c) + xmv^1(2 a b c + 2 a b^2 + 3 a^2 b + 2 b^2 c#> + c^6) + xmv^2(a b + b c)is a series in powers ofx-a-b. We may perform aconsistency check by a second substitution, returning us to the originalexpression:
If functionseries() is given a variable name ending in_m_foo, wherefoo is any variable name, thenthis is typeset as(x-foo). For example:
as.mvp('x^3 + x*a')%>%subs(x="x_m_a + a")%>%series("x_m_a")#> (x-a)^0(a^2 + a^3) + (x-a)^1(a + 3 a^2) + (x-a)^2(3 a) + (x-a)^3(1)So above we see the expansion of\(x^2+ax\) in powers of\(x-a\). If we want to see the expansion of amvp in terms of a more complicated expression then it is better to use anonce variablev:
as.mvp('x^2 + x*a+b^3')%>%subs(x="x_m_v + a^2+b")%>%series("x_m_v")#> (x-v)^0(a b + 2 a^2 b + a^3 + a^4 + b^2 + b^3) + (x-v)^1(a + 2 a^2#> + 2 b) + (x-v)^2(1)where it is understood that\(v=a+b^2\). Functiontaylor()is a convenience wrapper that does some of the above in one step:
p<-as.mvp("1+x-x*y+a")^2taylor(p,'x','a')#> (x-a)^0(1 + 4 a - 2 a y + 4 a^2 - 4 a^2 y + a^2 y^2) + (x-a)^1(2 +#> 4 a - 6 a y + 2 a y^2 - 2 y) + (x-a)^2(1 - 2 y + y^2)But it’s not as good as I expected it to be and frankly it’soverkill.
Given a multivariate polynomial, one often needs to extract certainterms. Because the terms of anmvp object have animplementation-dependent order, this can be difficult. But we can usefunctiononevarpow():
Themvp package handles negative powers, although theidiom is not perfect and I’m still working on it. There is theinvert() function:
(p<-as.mvp("1+x+x^2 y"))#> mvp object algebraically equal to#> 1 + x + x^2 yinvert(p)#> mvp object algebraically equal to#> 1 + x^-2 y^-1 + x^-1In the above,p is a regular multivariate polynomialwhich includes negative powers. It obeys the same arithmetic rules asother mvp objects:
disordR packageIt is possible to examine the coefficients of anmvpobject:
a<-as.mvp("5 + 8*x^2*y - 13*y*x^2 + 11*z - 3*x*yz")a#> mvp object algebraically equal to#> 5 - 3 x yz - 5 x^2 y + 11 zcoeffs(a)#> A disord object with hash b0695ec1b46bb6da92f016f17e49971e1036d786 and elements#> [1] 5 -3 -5 11#> (in some order)Above, note that the result ofcoeffs() is adisord object, defined in thedisordR package(Hankin2022a). The order of the elements unspecified as theSTL map class holds the keys and values in animplementation-specific order. This device stops the user from illegaloperations on the coefficients. For example, suppose we had anothermvp object,b:
b<- a*2b#> mvp object algebraically equal to#> 10 - 6 x yz - 10 x^2 y + 22 zcoeffs(a)+coeffs(b)#> coeffs(a) + coeffs(b)#> Error in check_matching_hash(e1, e2, match.call()):#> hash codes b0695ec1b46bb6da92f016f17e49971e1036d786 and 286de701fbfc3887e4489c2b07127f0811071f43 do not matchabove, we get an error because the coefficients ofa andb are possibly stored in a different order and thereforevector addition makes no sense. However, we can operate on coefficientsof a singlemvp object at will:
coeffs(a)>0#> A disord object with hash b0695ec1b46bb6da92f016f17e49971e1036d786 and elements#> [1] TRUE FALSE FALSE TRUE#> (in some order)coeffs(a)+coeffs(a)^4#> A disord object with hash b0695ec1b46bb6da92f016f17e49971e1036d786 and elements#> [1] 630 78 620 14652#> (in some order)Extraction also works but subject to standarddisordRidiom restrictions:
coeffs(a)[coeffs(a)>0]#> A disord object with hash 54f03094a678b0e4328942e555e5a308c29f40fe and elements#> [1] 5 11#> (in some order)But “mixing” objects is forbidden:
coeffs(a)[coeffs(b)>0]#> .local(x = x, i = i, j = j, drop = drop)#> Error in check_matching_hash(x, i, match.call()):#> hash codes b0695ec1b46bb6da92f016f17e49971e1036d786 and 286de701fbfc3887e4489c2b07127f0811071f43 do not matchExtraction methods work, again subject todisordRrestrictions:
coeffs(a)[coeffs(a)<0]<-coeffs(a)[coeffs(a)<0]+1000# add 1000 to every negative coefficienta#> mvp object algebraically equal to#> 5 + 997 x yz + 995 x^2 y + 11 zIn cases like this where the replacement object is complicated, usingmagrittr would simplify the idiom and reduce theopportunity for error:
coeffs(b)[coeffs(b)%%3==1]%<>%`+`(100)# add 100 to every element equal to 1 modulo 3b#> mvp object algebraically equal to#> 110 - 6 x yz - 10 x^2 y + 122 zOne good use for this is to “zap” small elements:
x<-as.mvp("1 - 0.11*x + 0.005*x*y")^2x#> mvp object algebraically equal to#> 1 - 0.22 x + 0.01 x y + 0.0121 x^2 - 0.0011 x^2 y + 0.000025 x^2 y^2Then we can zap as follows:
cx<-coeffs(x)cx[abs(cx)<0.01]<-0coeffs(x)<- cxx#> mvp object algebraically equal to#> 1 - 0.22 x + 0.01 x y + 0.0121 x^2(I should write a method forzapsmall() that doesthis)
We can see the generating function for a chess knight:
knight(2)#> mvp object algebraically equal to#> a^-2 b^-1 + a^-2 b + a^-1 b^-2 + a^-1 b^2 + a b^-2 + a b^2 + a^2 b^-1 + a^2 bHow many ways are there for a 4D knight to return to its startingsquare after four moves? Answer:
mvp Package.”https://arxiv.org/abs/2210.15991; arXiv.https://doi.org/10.48550/ARXIV.2210.15991.