Polynomial Library in BsplineQuantReg

Introduction

The BsplineQuantReg package provides a comprehensive set of tools for polynomial manipulation, implemented in pure R with a consistent interface. All polynomial functions follow the decreasing power order convention, where a polynomial \(p(x) = a_0 + a_1 x + a_2 x^2\) is represented as c(a_2, a_1, a_0).

This vignette covers: - Basic polynomial operations (addition, multiplication, evaluation, differentiation) - Piecewise polynomial (PP) form - Callable and non-callable PP objects - Extracting parameters from PP objects

Basic Polynomial Operations

Polynomial Multiplication

# (1 + x) * (1 + x) = 1 + 2x + x^2
p1 <- c(1, 1)  # 1 + x
p2 <- c(1, 1)  # 1 + x
product <- polymul(p1, p2)
print(product)  # c(1, 2, 1) → 1 + 2x + x^2
## [1] 1 2 1
# (x^2 + 2x + 1) * (x - 1) = x^3 + x^2 - x - 1
p3 <- c(1, 2, 1)  # x^2 + 2x + 1
p4 <- c(1, -1)    # x - 1
polymul(p3, p4)
## [1]  1  1 -1 -1

Polynomial Addition

# (1 + x) + (1 - x) = 2
polyadd(c(1, 1), c(1, -1))
## [1] 2 0
# (x^2 + 1) + (x + 1) = x^2 + x + 2
polyadd(c(1, 0, 1), c(1, 1))
## [1] 1 1 2

Polynomial Evaluation

# P(x) = 1 + x + x^2
p <- c(1, 1, 1)  # x^2 + x + 1
poly_eval(p, c(0, 1, 2))  # returns c(1, 3, 7)
## [1] 1 3 7
# Evaluate at many points
x <- seq(-2, 2, length.out = 10)
y <- poly_eval(c(1, 0, -1), x)  # 1 - x^2
plot(x, y, type = "l", main = "1 - x^2")

Polynomial Differentiation

# P(x) = x^2 → P'(x) = 2x
p <- c(1, 0, 0)  # x^2
polyderiv(p, der = 1)  # c(2, 0) → 2x
## [1] 2 0
# Second derivative: P''(x) = 2
polyderiv(p, der = 2)  # c(2)
## [1] 2
# Higher order derivatives
p <- c(1, 2, 3, 4)  # 4x^3 + 3x^2 + 2x + 1
polyderiv(p, der = 2)  # 24x + 6
## [1] 6 4

Polynomial Reduction

# Remove leading zeros
reduce_pol(c(0, 0, 1, 2, 1))  # c(1, 2, 1) → x^2 + 2x + 1
## [1] 1 2 1
reduce_pol(c(0, 0, 0, 5))     # c(5) → constant polynomial
## [1] 5

Taylor Basis Change

# Convert (x-1)^2 to expansion around 0
# (x-1)^2 = x^2 - 2x + 1
coeffs_a <- c(1, 0, 0)  # (x-1)^2 in basis centered at a=1
change_polynomial_base_taylor(coeffs_a, a = 1, b = 0)
## [1]  1 -2  1
# Returns c(1, -2, 1) → x^2 - 2x + 1

Piecewise Polynomial (PP) Form

The PP form represents a function defined piecewise by polynomials on different intervals.

Creating a PP Object

# Create a piecewise polynomial with two intervals
# Interval 1: x^2 on [0, 1]
# Interval 2: 2x - 1 on [1, 2]
coeff <- matrix(c(
  1, 0, 0,   # x^2
  0, 2, -1   # 2x - 1
), nrow = 2, byrow = TRUE)

knots <- c(0, 1, 2)

pp <- makpp(coeff, knots)
print(pp)
## Piecewise Polynomial (PP) (non-callable)
## ================================
##   $degree: 2 
##   $knot: 0 1 2 
##  coefficients dimension: 2 x 3 
## 
##   $coeff:
##     Intervals 1 : 1, 0, 0 
##     Intervals 2 : 0, 2, -1 
##  Usage: pp_eval(pp, x_values)
# Evaluate the PP form
x <- seq(0, 2, length.out = 100)
y <- evalpp(pp, x)

plot(x, y, type = "l", main = "Piecewise Polynomial")
abline(v = knots, col = "red", lty = 2)

Non-Callable PP Objects

# A non-callable PP object is a list with components:
# - coeff: matrix of polynomial coefficients
# - knot: knot positions
# - degree: polynomial degree

pp <- makpp(coeff, knots, callable = FALSE)
class(pp)  # "non_callable_pp"
## [1] "non_callable_pp"
# Access components
pp$coeff
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2   -1
pp$knot
## [1] 0 1 2
pp$degree
## [1] 2

Callable PP Objects

A callable PP object can be evaluated directly like a function.

# Create a callable PP
pp_call <- makpp(coeff, knots, callable = TRUE)
class(pp_call)  # "callable_pp" "function"
## [1] "callable_pp" "function"
# Evaluate directly
x <- seq(0, 2, length.out = 10)
y <- pp_call(x)

# Print shows information
print(pp_call)
## Callable Piecewise Polynomial (PP) Object
## ==========================================
##   Degree: 2 
##   Intervals: 2 
##   Knots ( 3 ) : 0, 1, 2 
##  coefficients dimension: 2 x 3 
## 
##   $coeff:
##     Intervals 1 : 1, 0, 0 
##     Intervals 2 : 0, 2, -1 
##  Usage: pp(x_values) or evalpp(pp, x_values)
# Can be used in plots
plot(pp_call, xlim = c(0, 2))

Create a callable PP

pp_call <- makpp(coeff, knots, callable = TRUE) class(pp_call) # “callable_pp” “function”

Evaluate directly

x <- seq(0, 2, length.out = 10) y <- pp_call(x)

Can be used in plots

plot(pp_call, xlim = c(0, 2))

Extrapolating values

When evalpp() is called with x values outside the knot range (x < knots[1] or x > knots[length(knots)]), it automatically performs extrapolation by extension using the border polynomial piece (the first or last polynomial in the PP structure). For example, if the knot vector is c(0, 1, 2), points with x < 0 will be evaluated using the polynomial from the first interval [0, 1], while points with x > 2 will use the polynomial from the last interval [1, 2]. This is equivalent to “extending” the first and last piece of the spline. A warning message is printed to alert the user when extrapolation occurs.

PP with two intervals: x^2 on [0,1], 2x-1 on [1,2]

coeff <- matrix(c(1, 0, 0, 0, 2, -1), nrow = 2, byrow = TRUE) knots <- c(0, 1, 2) pp <- makpp(coeff, knots, callable = TRUE)

Evaluate inside and outside the knot range

x_test <- c(-0.5, 0.5, 1.5, 2.5) y_test <- pp(x_test) data.frame(x = x_test, y = y_test) ### The values at -0.5 and 2.5 are extrapolated


### Extracting Parameters


``` r
### For callable PP objects, use get_parameters()
params <- get_parameters(pp_call)
print(params$degree)
## [1] 2
print(params$knot)
## [1] 0 1 2
print(params$coeff)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2   -1
### For non-callable PP, access directly
pp=makpp(pp,callable=FALSE,verbose=TRUE)
## Already a 'non_callable_pp'
## 
pp$degree
## [1] 2
pp$knot
## [1] 0 1 2
pp$coeff
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2   -1

Advanced: PP Form Manipulation

Extracting a Single Interval Polynomial

### Get polynomial for interval i
pp=makpp(pp,callable=FALSE)
interval_poly <- pp$coeff[1, ]  ### First interval
print(interval_poly)
## [1] 1 0 0
### Evaluate on a subinterval
x_sub <- seq(0, 1, length.out = 20)
y_sub <- poly_eval(interval_poly, x_sub)

Combining PP Objects

# Create two PP objects and combine them
# This is useful for constructing complex piecewise functions
coeff1 <- matrix(c(1, 0), nrow = 1)
coeff2 <- matrix(c(0, 1), nrow = 1)

pp1 <- makpp(coeff1, c(0, 1))
pp2 <- makpp(coeff2, c(1, 2))

### Combine coefficients and knots
combined_coeff <- rbind(coeff1, coeff2)
combined_knots <- c(0, 1, 2)
pp_combined <- makpp(combined_coeff, combined_knots)

### Evaluate combined function
x <- seq(0, 2, length.out = 100)
y_combined <- evalpp(pp_combined, x)
plot(x, y_combined, type = "l")
abline(v = c(0, 1, 2), col = "red", lty = 2)

Visualizing PP with Basis Functions

### Create a PP with multiple pieces
knots <- seq(0, 1, length.out = 5)
degree <- 2
n_pieces <- length(knots) - 1

### Generate random coefficients
coeff <- matrix(rnorm(n_pieces * (degree + 1)), nrow = n_pieces)

pp <- makpp(coeff, knots, callable = TRUE)

### Evaluate and plot
x <- seq(0, 1, length.out = 200)
y <- pp(x)

plot(x, y, type = "l", lwd = 2,
     main = "Random Piecewise Quadratic")
abline(v = knots, col = "red", lty = 2)
grid()

Comparison: PP vs B-spline

### Create a B-spline and convert to PP form
sn <- c(0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- runif(basis$n_splines)

### Convert to PP
pp_from_bspline <- Bsplinetopp(basis)

### Both represent the same function
x <- seq(0, 1, length.out = 100)
y_bspline <- spline_eval(basis, x)
y_pp <- evalpp(pp_from_bspline, x)

### They should match
max(abs(y_bspline - y_pp))
## [1] 43.43534

Summary

Feature Function Description
Polynomial multiplication polymul() Multiply two polynomials
Polynomial addition polyadd() Add two polynomials
Polynomial evaluation poly_eval() Evaluate at points
Polynomial derivative polyderiv() Differentiate polynomial
Taylor shift change_polynomial_base_taylor() Change expansion point
PP creation makpp() Create PP form
PP evaluation evalpp() Evaluate PP at points
Callable PP makpp(callable = TRUE) PP as function
Extract params get_parameters() Extract from callable

All functions are implemented in pure R and follow matlab consistent conventions, making them robust and easy to use for B-spline manipulation and quantile regression.

Displaying polynomial in a human readable form

Introduction

The show_poly() function displays a polynomial as a human-readable mathematical equation. It supports both canonical and local bases, and can convert between different expansion points using Taylor’s formula.

A polynomial \(P(x) = \sum_{k=0}^d c_k (x-a)^k\) represented by coefficients c(c_d, c_{d-1}, ..., c_0) (decreasing power order) can be displayed in any basis \((x-b)^k\).

Basic Usage

Canonical Basis (Default)

By default, show_poly() displays the polynomial in the canonical basis \(1, x, x^2, ...\):

# P(x) = 3x^3 - 2x^2 + x - 5
p <- c(3, -2, 1, -5)
show_poly(p)
## [1] "3*x^3-2*x^2+1*x^1-5"
# "3*x^3 - 2*x^2 + 1*x - 5"

Working with Local Bases

Display in Local Basis

To display a polynomial in the local basis \((x-a)^k\), use the b parameter:

# P(x) = 3(x-2)^3 - 2(x-2)^2 + (x-2) - 5
p <- c(3, -2, 1, -5)
show_poly(p, b = 2)
## [1] "3*(x-2)^3+16*(x-2)^2+29*(x-2)^1+13"
# "3*(x-2)^3 - 2*(x-2)^2 + 1*(x-2) - 5"

Converting Between Bases

The a parameter specifies the basis in which the coefficients are given. The function automatically converts using change_polynomial_base_taylor():

# Coefficients are in basis (x-2)^k, display in canonical basis
p <- c(3, -2, 1, -5)  # In basis (x-2)^k
show_poly(p, a = 2, b = 0)
## [1] "3*x^3-20*x^2+45*x^1-39"
# This expands to: 3x^3 - 18x^2 + 37x - 27

Examples with Different Bases

Positive and Negative Centers

# Center at positive value
p <- c(1, 0, 2)  # (x-3)^2 + 2
show_poly(p, b = 3)
## [1] "1*(x-3)^2+6*(x-3)^1+11"
# "(x-3)^2 + 2"

# Center at negative value
p <- c(1, 0, 2)  # (x+2)^2 + 2
show_poly(p, b = -2)
## [1] "1*(x+2)^2-4*(x+2)^1+6"
# "(x+2)^2 + 2"

Controlling Output Format

Number of Digits

p <- c(1/3, -2/7, 1/5)
show_poly(p, digits = 2)
## [1] "0.33*x^2-0.29*x^1+0.2"
# "0.33*x^2 - 0.29*x + 0.2"

show_poly(p, digits = 4)
## [1] "0.3333*x^2-0.2857*x^1+0.2"
# "0.3333*x^2 - 0.2857*x + 0.2"

Working with B-splines

When working with B-splines, show_poly() helps visualize the polynomial pieces:

# Create a B-spline and display its polynomial form
sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75)

# Convert to PP
pp <- Bsplinetopp(basis, callable = FALSE)

# Display each piece
cat("Piece 1: ", show_poly(pp$coeff[1, ], b = 0), "\n")
## Piece 1:  -133.3333*x^3+91.6667*x^2-15*x^1+1
cat("Piece 2: ", show_poly(pp$coeff[2, ], b = 0), "\n")
## Piece 2:  28.4203*x^3-28.3333*x^2+4*x^1+1.15
cat("Piece 3: ", show_poly(pp$coeff[3, ], b = 0), "\n")
## Piece 3:  43.0325*x^3-2.7551*x^2-5.3265*x^1+0.5673

Integration with show_pp()

The show_poly() function is used internally by show_pp() to display piecewise polynomials:

# Create a PP with two intervals
knot <- c(0, 0.5, 1)
coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE)
pp <- makpp(coeff, knot)

# Display in local basis
show_pp(pp, local = TRUE, digits = 3)
##      [,1]                  [,2]                       
## [1,] "  [0.0000, 0.5000] " "1*x^2+2*x^1+0.5"          
## [2,] "  [0.5000, 1.0000] " "0*(x-0.5)^2+1*(x-0.5)^1-1"
# Create a PP with two intervals
knot <- c(0, 0.5, 1)
coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE)
pp <- makpp(coeff, knot)

# Display in canonical basis (1, x, x², ...)
show_pp(pp, local = FALSE)
##      [,1]                  [,2]             
## [1,] "  [0.0000, 0.5000] " "1*x^2+2*x^1+0.5"
## [2,] "  [0.5000, 1.0000] " "0*x^2+1*x^1-1.5"
# [0.000, 0.500] 0.5x^2 + 2x + 1
# [0.500, 1.000] -1x^2 + 1x + 0

# Display in local basis ((x-a)^i)
show_pp(pp, local = TRUE)
##      [,1]                  [,2]                       
## [1,] "  [0.0000, 0.5000] " "1*x^2+2*x^1+0.5"          
## [2,] "  [0.5000, 1.0000] " "0*(x-0.5)^2+1*(x-0.5)^1-1"
# [0.000, 0.500] 0.5(x-0)^2 + 2(x-0) + 1
# [0.500, 1.000] -1(x-0.5)^2 + 1(x-0.5) + 0

For spline objects, the conversion to PP is automatic:

# Create a B-spline and display its polynomial form
sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75)
basis<-make_spline(basis)

show_pp(basis, local = TRUE, verbose = TRUE)
## Piecewise Polynomial (PP) (non-callable)
## ================================
##   $degree: 3 
##   $knot: 0 0.3 0.6 1 
##  coefficients dimension: 3 x 4 
## 
##   $coeff:
##     Intervals 1 : -133.3333, 91.6667, -15, 1 
##     Intervals 2 : 28.4203, -28.3333, 4, 1.15 
##     Intervals 3 : 43.0325, -2.7551, -5.3265, 0.5673 
##  Usage: pp_eval(pp, x_values)
##      [,1]                 
## [1,] "  [0.0000, 0.3000] "
## [2,] "  [0.3000, 0.6000] "
## [3,] "  [0.6000, 1.0000] "
##      [,2]                                                        
## [1,] "-133.3333*x^3+91.6667*x^2-15*x^1+1"                        
## [2,] "28.4203*(x-0.3)^3-28.3333*(x-0.3)^2+4*(x-0.3)^1+1.15"      
## [3,] "43.0325*(x-0.6)^3-2.7551*(x-0.6)^2-5.3265*(x-0.6)^1+0.5673"
# PP Information:
# Degree: 3
# Knots: 0, 0.3, 0.6, 1
# Number of intervals: 3
# 
# [0.000, 0.300] 1(x-0)^3 + 2.5(x-0)^2 + 1.5(x-0) + 1
# [0.300, 0.600] -4.63(x-0.3)^3 + ... 
# [0.600, 1.000] ...

Customization Options The display functions offer several customization options:

local: Control the basis (TRUE for local, FALSE for canonical)

digits: Number of significant digits (default: 4)

verbose: Display additional PP information

Summary show_pp()

Parameter Description Default
obj Polynomial coefficients (decreasing power order) Required
a Base of input coefficients \((x-a)^k\) 0
b Base for output display \((x-b)^k\) 0
digits Number of significant digits 4
verbose Print additional information FALSE