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%}
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.
STL map
classA 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).
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"))
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
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")
Substitution works well with pipes:
as.mvp("a+b") %>% subs(a="a^2+b^2") %>% subs(b="x^6")
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 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
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.
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)
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")
disordR
packageIt 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)
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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.