# mpoly v1.1.1

0

0th

Percentile

## Symbolic Computation and More with Multivariate Polynomials

Symbolic computing with multivariate polynomials in R.

# mpoly

## Specifying polynomials

mpoly is a simple collection of tools to help deal with multivariate polynomials symbolically and functionally in R. Polynomials are defined with the mp() function:

library("mpoly")
# Registered S3 methods overwritten by 'ggplot2':
#   method         from
#   [.quosures     rlang
#   c.quosures     rlang
#   print.quosures rlang
mp("x + y")
# x  +  y

mp("(x + 4 y)^2 (x - .25)")
# x^3  -  0.25 x^2  +  8 x^2 y  -  2 x y  +  16 x y^2  -  4 y^2


Term orders are available with the reorder function:

(p <- mp("(x + y)^2 (1 + x)"))
# x^3  +  x^2  +  2 x^2 y  +  2 x y  +  x y^2  +  y^2

reorder(p, varorder = c("y","x"), order = "lex")
# y^2 x  +  y^2  +  2 y x^2  +  2 y x  +  x^3  +  x^2

reorder(p, varorder = c("x","y"), order = "glex")
# x^3  +  2 x^2 y  +  x y^2  +  x^2  +  2 x y  +  y^2


Vectors of polynomials (mpolyList’s) can be specified in the same way:

mp(c("(x+y)^2", "z"))
# x^2  +  2 x y  +  y^2
# z


## Polynomial parts

You can extract pieces of polynoimals using the standard [ operator, which works on its terms:

p
# x^3

p[1:3]
# x^3  +  x^2  +  2 x^2 y

p[-1]
# x^2  +  2 x^2 y  +  2 x y  +  x y^2  +  y^2


There are also many other functions that can be used to piece apart polynomials, for example the leading term (default lex order):

LT(p)
# x^3

LC(p)
#  1

LM(p)
# x^3


You can also extract information about exponents:

exponents(p)
# []
# x y
# 3 0
#
# []
# x y
# 2 0
#
# []
# x y
# 2 1
#
# []
# x y
# 1 1
#
# []
# x y
# 1 2
#
# []
# x y
# 0 2

multideg(p)
# x y
# 3 0

totaldeg(p)
#  3

monomials(p)
# x^3
# x^2
# 2 x^2 y
# 2 x y
# x y^2
# y^2


## Polynomial arithmetic

Arithmetic is defined for both polynomials (+, -, * and ^)…

p1 <- mp("x + y")

p2 <- mp("x - y")

p1 + p2
# 2 x

p1 - p2
# 2 y

p1 * p2
# x^2  -  y^2

p1^2
# x^2  +  2 x y  +  y^2


… and vectors of polynomials:

(ps1 <- mp(c("x", "y")))
# x
# y

(ps2 <- mp(c("2 x", "y + z")))
# 2 x
# y  +  z

ps1 + ps2
# 3 x
# 2 y  +  z

ps1 - ps2
# -1 x
# -1 z

ps1 * ps2
# 2 x^2
# y^2  +  y z


## Some calculus

You can compute derivatives easily:

p <- mp("x + x y + x y^2")

deriv(p, "y")
# x  +  2 x y

# y^2  +  y  +  1
# 2 y x  +  x


## Function coercion

You can turn polynomials and vectors of polynomials into functions you can evaluate with as.function(). Here’s a basic example using a single multivariate polynomial:

f <- as.function(mp("x + 2 y")) # makes a function with a vector argument
# f(.) with . = (x, y)

f(c(1,1))
#  3

f <- as.function(mp("x + 2 y"), vector = FALSE) # makes a function with all arguments
# f(x, y)

f(1, 1)
#  3


Here’s a basic example with a vector of multivariate polynomials:

(p <- mp(c("x", "2 y")))
# x
# 2 y

f <- as.function(p)
# f(.) with . = (x, y)

f(c(1,1))
#  1 2

f <- as.function(p, vector = FALSE)
# f(x, y)

f(1, 1)
#  1 2


Whether you’re working with a single multivariate polynomial or a vector of them (mpolyList), if it/they are actually univariate polynomial(s), the resulting function is vectorized. Here’s an example with a single univariate polynomial.

f <- as.function(mp("x^2"))
# f(.) with . = x

f(1:3)
#  1 4 9

(mat <- matrix(1:4, 2))
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

f(mat) # it's vectorized properly over arrays
#      [,1] [,2]
# [1,]    1    9
# [2,]    4   16


Here’s an example with a vector of univariate polynomials:

(p <- mp(c("t", "t^2")))
# t
# t^2

f <- as.function(p)
f(1)
#  1 1

f(1:3)
#      [,1] [,2]
# [1,]    1    1
# [2,]    2    4
# [3,]    3    9


You can use this to visualize a univariate polynomials like this:

library("tidyverse"); theme_set(theme_classic())

f <- as.function(mp("(x-2) x (x+2)"))
# f(.) with . = x
x <- seq(-2.5, 2.5, .1)

qplot(x, f(x), geom = "line") For multivariate polynomials, it’s a little more complicated:

f <- as.function(mp("x^2 - y^2"))
# f(.) with . = (x, y)
s <- seq(-2.5, 2.5, .1)
df <- expand.grid(x = s, y = s)
df$f <- apply(df, 1, f) qplot(x, y, data = df, geom = "raster", fill = f) Using tidyverse-style coding (install tidyverse packages with install.packages("tidyverse")), this looks a bit cleaner: f <- as.function(mp("x^2 - y^2"), vector = FALSE) # f(x, y) seq(-2.5, 2.5, .1) %>% list("x" = ., "y" = .) %>% cross_df() %>% mutate(f = f(x, y)) %>% ggplot(aes(x, y, fill = f)) + geom_raster() ## Algebraic geometry Grobner bases are no longer implemented in mpoly; they’re now in m2r. # polys <- mp(c("t^4 - x", "t^3 - y", "t^2 - z")) # grobner(polys)  Homogenization and dehomogenization: (p <- mp("x + 2 x y + y - z^3")) # x + 2 x y + y - z^3 (hp <- homogenize(p)) # x t^2 + 2 x y t + y t^2 - z^3 dehomogenize(hp, "t") # x + 2 x y + y - z^3 homogeneous_components(p) # x + y # 2 x y # -1 z^3  ## Special polynomials mpoly can make use of many pieces of the polynom and orthopolynom packages with as.mpoly() methods. In particular, many special polynomials are available. #### Chebyshev polynomials You can construct Chebyshev polynomials as follows: chebyshev(1) # Loading required package: polynom # # Attaching package: 'polynom' # The following object is masked from 'package:mpoly': # # LCM # x chebyshev(2) # -1 + 2 x^2 chebyshev(0:5) # 1 # x # 2 x^2 - 1 # 4 x^3 - 3 x # 8 x^4 - 8 x^2 + 1 # 16 x^5 - 20 x^3 + 5 x  And you can visualize them: s <- seq(-1, 1, length.out = 201); N <- 5 (chebPolys <- chebyshev(0:N)) # 1 # x # 2 x^2 - 1 # 4 x^3 - 3 x # 8 x^4 - 8 x^2 + 1 # 16 x^5 - 20 x^3 + 5 x df <- as.function(chebPolys)(s) %>% cbind(s, .) %>% as.data.frame() names(df) <- c("x", paste0("T_", 0:N)) mdf <- df %>% gather(degree, value, -x) qplot(x, value, data = mdf, geom = "path", color = degree) #### Jacobi polynomials s <- seq(-1, 1, length.out = 201); N <- 5 (jacPolys <- jacobi(0:N, 2, 2)) # 1 # 5 x # 17.5 x^2 - 2.5 # 52.5 x^3 - 17.5 x # 144.375 x^4 - 78.75 x^2 + 4.375 # 375.375 x^5 - 288.75 x^3 + 39.375 x df <- as.function(jacPolys)(s) %>% cbind(s, .) %>% as.data.frame names(df) <- c("x", paste0("P_", 0:N)) mdf <- df %>% gather(degree, value, -x) qplot(x, value, data = mdf, geom = "path", color = degree) + coord_cartesian(ylim = c(-25, 25)) #### Legendre polynomials s <- seq(-1, 1, length.out = 201); N <- 5 (legPolys <- legendre(0:N)) # 1 # x # 1.5 x^2 - 0.5 # 2.5 x^3 - 1.5 x # 4.375 x^4 - 3.75 x^2 + 0.375 # 7.875 x^5 - 8.75 x^3 + 1.875 x df <- as.function(legPolys)(s) %>% cbind(s, .) %>% as.data.frame names(df) <- c("x", paste0("P_", 0:N)) mdf <- df %>% gather(degree, value, -x) qplot(x, value, data = mdf, geom = "path", color = degree) #### Hermite polynomials s <- seq(-3, 3, length.out = 201); N <- 5 (hermPolys <- hermite(0:N)) # 1 # x # x^2 - 1 # x^3 - 3 x # x^4 - 6 x^2 + 3 # x^5 - 10 x^3 + 15 x df <- as.function(hermPolys)(s) %>% cbind(s, .) %>% as.data.frame names(df) <- c("x", paste0("He_", 0:N)) mdf <- df %>% gather(degree, value, -x) qplot(x, value, data = mdf, geom = "path", color = degree) #### (Generalized) Laguerre polynomials s <- seq(-5, 20, length.out = 201); N <- 5 (lagPolys <- laguerre(0:N)) # 1 # -1 x + 1 # 0.5 x^2 - 2 x + 1 # -0.1666667 x^3 + 1.5 x^2 - 3 x + 1 # 0.04166667 x^4 - 0.6666667 x^3 + 3 x^2 - 4 x + 1 # -0.008333333 x^5 + 0.2083333 x^4 - 1.666667 x^3 + 5 x^2 - 5 x + 1 df <- as.function(lagPolys)(s) %>% cbind(s, .) %>% as.data.frame names(df) <- c("x", paste0("L_", 0:N)) mdf <- df %>% gather(degree, value, -x) qplot(x, value, data = mdf, geom = "path", color = degree) + coord_cartesian(ylim = c(-25, 25)) #### Bernstein polynomials Bernstein polynomials are not in polynom or orthopolynom but are available in mpoly with bernstein(): bernstein(0:4, 4) # x^4 - 4 x^3 + 6 x^2 - 4 x + 1 # -4 x^4 + 12 x^3 - 12 x^2 + 4 x # 6 x^4 - 12 x^3 + 6 x^2 # -4 x^4 + 4 x^3 # x^4 s <- seq(0, 1, length.out = 101) N <- 5 # number of bernstein polynomials to plot (bernPolys <- bernstein(0:N, N)) # -1 x^5 + 5 x^4 - 10 x^3 + 10 x^2 - 5 x + 1 # 5 x^5 - 20 x^4 + 30 x^3 - 20 x^2 + 5 x # -10 x^5 + 30 x^4 - 30 x^3 + 10 x^2 # 10 x^5 - 20 x^4 + 10 x^3 # -5 x^5 + 5 x^4 # x^5 df <- as.function(bernPolys)(s) %>% cbind(s, .) %>% as.data.frame names(df) <- c("x", paste0("B_", 0:N)) mdf <- df %>% gather(degree, value, -x) qplot(x, value, data = mdf, geom = "path", color = degree) You can use the bernstein_approx() function to compute the Bernstein polynomial approximation to a function. Here’s an approximation to the standard normal density: p <- bernstein_approx(dnorm, 15, -1.25, 1.25) round(p, 4) # -0.1624 x^2 + 0.0262 x^4 - 0.002 x^6 + 0.0001 x^8 + 0.3796 x <- seq(-3, 3, length.out = 101) df <- data.frame( x = rep(x, 2), y = c(dnorm(x), as.function(p)(x)), which = rep(c("actual", "approx"), each = 101) ) # f(.) with . = x qplot(x, y, data = df, geom = "path", color = which) ## Bezier polynomials and curves You can construct Bezier polynomials for a given collection of points with bezier(): points <- data.frame(x = c(-1,-2,2,1), y = c(0,1,1,0)) (bezPolys <- bezier(points)) # -10 t^3 + 15 t^2 - 3 t - 1 # -3 t^2 + 3 t  And viewing them is just as easy: df <- as.function(bezPolys)(s) %>% as.data.frame ggplot(aes(x = x, y = y), data = df) + geom_point(data = points, color = "red", size = 4) + geom_path(data = points, color = "red", linetype = 2) + geom_path(size = 2) Weighting is available also: points <- data.frame(x = c(1,-2,2,-1), y = c(0,1,1,0)) (bezPolys <- bezier(points)) # -14 t^3 + 21 t^2 - 9 t + 1 # -3 t^2 + 3 t df <- as.function(bezPolys, weights = c(1,5,5,1))(s) %>% as.data.frame ggplot(aes(x = x, y = y), data = df) + geom_point(data = points, color = "red", size = 4) + geom_path(data = points, color = "red", linetype = 2) + geom_path(size = 2) To make the evaluation of the Bezier polynomials stable, as.function() has a special method for Bezier polynomials that makes use of de Casteljau’s algorithm. This allows bezier() to be used as a smoother: s <- seq(0, 1, length.out = 201) df <- as.function(bezier(cars))(s) %>% as.data.frame qplot(speed, dist, data = cars) + geom_path(data = df, color = "red") ## Other stuff I’m starting to put in methods for some other R functions: set.seed(1) n <- 101 df <- data.frame(x = seq(-5, 5, length.out = n)) df$y <- with(df, -x^2 + 2*x - 3 + rnorm(n, 0, 2))

mod <- lm(y ~ x + I(x^2), data = df)
(p <- mod %>% as.mpoly %>% round)
# 1.983 x  -  1.01 x^2  -  2.709
qplot(x, y, data = df) +
stat_function(fun = as.function(p), colour = 'red')
# f(.) with . = x s <- seq(-5, 5, length.out = n)
df <- expand.grid(x = s, y = s) %>%
mutate(z = x^2 - y^2 + 3*x*y + rnorm(n^2, 0, 3))

(mod <- lm(z ~ poly(x, y, degree = 2, raw = TRUE), data = df))
#
# Call:
# lm(formula = z ~ poly(x, y, degree = 2, raw = TRUE), data = df)
#
# Coefficients:
#                           (Intercept)
#                             -0.070512
# poly(x, y, degree = 2, raw = TRUE)1.0
#                             -0.004841
# poly(x, y, degree = 2, raw = TRUE)2.0
#                              1.005307
# poly(x, y, degree = 2, raw = TRUE)0.1
#                              0.001334
# poly(x, y, degree = 2, raw = TRUE)1.1
#                              3.003755
# poly(x, y, degree = 2, raw = TRUE)0.2
#                             -0.999536
as.mpoly(mod)
# -0.004840798 x  +  1.005307 x^2  +  0.001334122 y  +  3.003755 x y  -  0.9995356 y^2  -  0.07051218


## Installation

• From CRAN: install.packages("mpoly")

• From Github (dev version):

# install.packages("devtools")
devtools::install_github("dkahle/mpoly")


## Acknowledgements

This material is based upon work partially supported by the National Science Foundation under Grant No. 1622449.

## Functions in mpoly

 Name Description insert Insert an element into a vector. legendre Legendre polynomials lissajous Lissajous polynomials is.wholenumber Test whether an object is a whole number deriv.mpoly Compute partial derivatives of a multivariate polynomial. mpoly Multivariate polynomials in R. print.mpolyList Pretty printing of a list of multivariate polynomials. components Polynomial components as.mpoly Convert an object to an mpoly reorder.mpoly Reorder a multivariate polynomial. mpoly-equal Determine whether two multivariate polynomials are equal. hermite Hermite polynomials jacobi Jacobi polynomials terms.mpoly Extract the terms of a multivariate polynomial. swap Swap polynomial indeterminates laguerre Generalized Laguerre polynomials mpolyArithmetic Arithmetic with multivariate polynomials permutations Determine all permutations of a set. mpolyList Define a collection of multivariate polynomials. plug Switch indeterminates in a polynomial homogenize Homogenize a polynomial mpoly-defunct Defunct mpoly functions mp Define a multivariate polynomial. solve_unipoly Solve a univariate mpoly with polyroot round.mpoly Round the coefficients of a polynomial eq_mp Convert an equation to a polynomial predicates mpoly predicate functions gradient Compute gradient of a multivariate polynomial. mpolyListArithmetic Element-wise arithmetic with vectors of multivariate polynomials. partitions Enumerate the partitions of an integer tuples Determine all n-tuples using the elements of a set. vars Determine the variables in a mpoly object. print.mpoly Pretty printing of multivariate polynomials. LCM Compute the least common multiple of two numbers. burst Enumerate integer r-vectors summing to n bezier Bezier polynomials bernstein Bernstein polynomials chebyshev Chebyshev polynomials as.function.mpolyList Change a vector of multivariate polynomials into a function. as.function.mpoly Change a multivariate polynomial into a function. bernstein-approx Bernstein polynomial approximation bezier_function Bezier function No Results!