Polynomials & Integration

Introduction

Aww hell, might as well finish this "Intro Series" off properly and briefly explore a few other areas within numerical computing with applications to applied statistical analysis. Notably, this series wouldn't be complete without at least discussing fitting polynomials and numerical integration. However, I find the code for these subjects a bit tedious so instead of a step-by-step walkthrough, I want to use this as an opportunity to play around with Rcpp and RcppArmadillo.

This piece is structured as follows,

• Polynomial Basics
• Rcpp: LaGrange Polynomials
• Rcpp: Trapezoidal Rule
• RcppArmadillo: Gauss Legendre Quadrature
• RcppArmadillo: Cubic Spline
• Links

Although some variant of the package has existed for longer, Rcpp has now existed, and continued to grow, in its current form since 2008 Seamless R and C++ Integration with Rcpp, Eddelbuettel (2013) with the aim of making R-C++ linkage as seamless as possible. As compared to the R C API, see here, Rcpp is relatively easy and straightforward to use. Moreover, there are now a large number of tutorials posted across the web, numerous SO questions and answers, and even a Rcpp-devel mailing list. Given the volume of introductory material available, I plan to omit basics and instead give a few simple examples of translating R code to Rcpp. Finally, instead of more, uh, "naive" implementations, I'm going to stick to translating functions from the pracma package as these algorithms are often direct translations of Matlab code from numerical computing text-books Github link....

Rcpp

Polynomial Basics

To begin, let's review polynomial interpolation. As you are well aware, the basic idea is to fit a line through the \(x,y\) plane such that \(p_n(x_i) = y_i \), with coefficients \(p_n(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_n x^n\). For each additional coefficient up to \(a_n\), we add an additional polynomial term, \(a_0 \ , a_1 \ , \dots , a_n \) such that the terms are equal to \(n\). Then, each equation, for each value of \(x\), is set equal to its corresponding \(y\) value, represented by a Vandermonde matrix See for a brief overview..., \[ V = \left( \begin{array}{cccccc} 1 & x_0 & x_0^2 & . & . & x_0^n \\ 1 & x_1 & x_1^2 & . & . & x_1^n \\ 1 & x_2 & x_2^2 & . & . & x_2^n \\ . & . & . & . & . & . \\ . & . & . & . & . & . \\ 1 & x_n & x_n^2 & . & . & x_n^n \\ \end{array} \right)\ \]

library(pracma)
library(polynom)

# We'll use these simply x,y vectors throughout..

y <- cos(1:10)
x <- 1:10

(X <- pracma::vander(x))

##             [,1]      [,2]     [,3]    [,4]   [,5]  [,6] [,7] [,8] [,9]
##  [1,]          1         1        1       1      1     1    1    1    1
##  [2,]        512       256      128      64     32    16    8    4    2
##  [3,]      19683      6561     2187     729    243    81   27    9    3
##  [4,]     262144     65536    16384    4096   1024   256   64   16    4
##  [5,]    1953125    390625    78125   15625   3125   625  125   25    5
##  [6,]   10077696   1679616   279936   46656   7776  1296  216   36    6
##  [7,]   40353607   5764801   823543  117649  16807  2401  343   49    7
##  [8,]  134217728  16777216  2097152  262144  32768  4096  512   64    8
##  [9,]  387420489  43046721  4782969  531441  59049  6561  729   81    9
## [10,] 1000000000 100000000 10000000 1000000 100000 10000 1000  100   10
##       [,10]
##  [1,]     1
##  [2,]     1
##  [3,]     1
##  [4,]     1
##  [5,]     1
##  [6,]     1
##  [7,]     1
##  [8,]     1
##  [9,]     1
## [10,]     1

Then, pick your poison to solve the system of equations and you have your interpolating polynomial.

solve(X,y)

##  [1]  1.332046e-06 -5.491529e-05  8.333860e-04 -5.090844e-03  2.495563e-03
##  [6]  9.043252e-02 -2.331367e-01 -1.430786e-02 -4.872132e-01  1.186343e+00

# Or, alternatively...

polynom::poly.calc(x,y)

## 1.186343 - 0.4872132*x - 0.01430786*x^2 - 0.2331367*x^3 + 0.09043252*x^4 +
## 0.002495563*x^5 - 0.005090844*x^6 + 0.000833386*x^7 - 5.491529e-05*x^8 +
## 1.332046e-06*x^9

LaGrange Polynomials

However, this is quite inefficient, so instead of translating this to Rcpp, we'll instead look at LaGrange polynomials, a more clever approximation approach. Here, the method, is represented by the following formula where \[f(x_k) = P(x_k),\ \sum_{k=0}^n\] and, \[L_{n,k}(x) = \frac{(x-x_0)(x-x_1)\dots(x-x_{k-1})(x-x_{k-1})\dots(x-x_n)}{(x_k - x_0)(x_k - x_1)\dots(x_k - x_{x-1})(x_k - x_{x+1})\dots(x_k - x_n)} = \prod_{i\neq k,\;i = 0}^n \frac{(x-x_i)}{(x_k - x_i)}\] Although this is formula appears complicated, note that for all points where \(i \ne k \), the function takes a value of zero and terms drop out These polynomials are also orthogonal which as we will see later, is a an important property . This type of indicator function is referred to as Kronecker's Delta, \(\delta_{ik} = \begin{cases} 1, & \text{if } i=k,\\ 0, & \text{if } i\neq k. \end{cases}\). Now, let's consider the following implementation via pracma::lagrangeInterp,

pracma::lagrangeInterp

## function (x, y, xs)
## {
##     stopifnot(is.numeric(x), is.numeric(y))
##     if (!is.numeric(xs))
##         stop("Argument 'xs' must be empty or a numeric vector.")
##     x <- c(x)
##     y <- c(y)
##     n <- length(x)
##     if (length(y) != n)
##         stop("Vectors 'x' and 'y' must be of the same length.")
##     A <- matrix(0, n, n)
##     A[, 1] <- y
##     for (i in 2:n) {
##         A[i:n, i] <- (A[i:n, i - 1] - A[i - 1, i - 1])/(x[i:n] -
##             x[i - 1])
##     }
##     ys <- A[n, n] * (xs - x[n - 1]) + A[n - 1, n - 1]
##     for (i in 2:(n - 1)) {
##         ys <- ys * (xs - x[n - i]) + A[n - i, n - i]
##     }
##     return(ys)
## }
## <bytecode: 0x000000000fabb438>
## <environment: namespace:pracma>

lagrangeInterp(x,y,5.5)

## [1] 0.7085453

As you can see, the code pretty closely follows the formula and utilizes a double-loop to iterate over all \(i,k\). Using arbitrary values, you can see the function mostly does a good job. Note, however, that this technique is for 'interpolation', not 'extrapolation' as it does a poor job outside the range of the data.

cos(1:15)[15]

## [1] -0.7596879

lagrangeInterp(x,y,15)

## [1] 547.3912

In terms of my translation to Rcpp, the only major change I made was to eliminate the matrix subsetting and instead use vectors ⊕ I also remove error checking... to simplify the code.

library(Rcpp)

Rcpp::cppFunction('double LaGrangeRcpp(NumericVector x, NumericVector y, double z) {

                  int n = x.size();
                  double P = 0;
                  NumericVector L(n,1.0);


                  for(int i = 0; i < n; i++) {
                   for(int j = 0; j < n; j++) {
                    if(i != j) {
                     L[i] = L[i] * ((z-x[j])/(x[i]-x[j]));

                    }
                   }

                   P =  (L[i] * y[i]) + P;
                  }

                  return P;

                  }')


LaGrangeRcpp(x,y,5.5)

## [1] 0.7085453

Clearly, this is considerably more simple than using the C API with the NumericVector declaration abstracting away having to deal with SEXPs and what not. Moreover, it is considerably faster than the R version, most likely due to the loops Note that a) these benchmarks are ad-hoc at best and b) they are something of a straw-man comparison as again, pracma is not designed for speed, but rather more educational/accuracy... .

library(microbenchmark)

microbenchmark(lagrangeInterp(x,y,5.5))

## Unit: microseconds
##                       expr    min     lq     mean median     uq    max
##  lagrangeInterp(x, y, 5.5) 21.151 21.517 22.77079 21.881 22.246 58.348
##  neval
##    100

microbenchmark(LaGrangeRcpp(x,y,5.5))

## Unit: nanoseconds
##                     expr min  lq     mean median   uq    max neval
##  LaGrangeRcpp(x, y, 5.5) 730 730 11086.72   1095 1095 994828   100

Trapezoid Rule

Now, what if we want to find the interval between two points for the above data pretending that those data points are the only thing we know. One basic, but effective method is the trapezoid rule, which subdivides the desired interval into \(n\) smaller intervals, approximate polynomials for said subintervals, and then integrate and sum the areas. The area of a trapezoid is, \(A = \frac{a+b}{2}h\), while the formula for the trapezoid rule is, \[\int_{a}^b f(x) = h\left [\frac{1}{2}f(x_0) + \sum_{i=1}^{n-1}f(x_i) + \frac{1}{2}f(x_n)\right] \] Below is a obligatory pictoral representation as well as code for pracma::trapz,

pracma::trapz

## function (x, y)
## {
##     if (missing(y)) {
##         if (length(x) == 0)
##             return(0)
##         y <- x
##         x <- seq(along = x)
##     }
##     if (length(x) == 0 && length(y) == 0)
##         return(0)
##     if (!(is.numeric(x) || is.complex(x)) || !(is.numeric(y) ||
##         is.complex(y)))
##         stop("Arguments 'x' and 'y' must be real or complex vectors.")
##     m <- length(x)
##     if (length(y) != m)
##         stop("Arguments 'x', 'y' must be vectors of the same length.")
##     if (m <= 1)
##         return(0)
##     xp <- c(x, x[m:1])
##     yp <- c(numeric(m), y[m:1])
##     n <- 2 * m
##     p1 <- sum(xp[1:(n - 1)] * yp[2:n]) + xp[n] * yp[1]
##     p2 <- sum(xp[2:n] * yp[1:(n - 1)]) + xp[1] * yp[n]
##     return(0.5 * (p1 - p2))
## }
## <bytecode: 0x000000000f2b22e8>
## <environment: namespace:pracma>

Again, for translation purposes, I omit type/length checks and also preallocate space as opposed to growing the vectors \(xp\) and \(yp\). Pay special attention to the numerous Rcpp index/access functions (e.g., x.begin()) which do differ substantially from R.

Rcpp::cppFunction('double TrapRcpp(NumericVector x, NumericVector y) {

  int m = x.length(); // C++ OOP uses frequent method chaining

  NumericVector xp(m * 2);
  NumericVector yp(m * 2);

  // Using basic standard template library functions here...
  std::copy(x.begin(),x.end(),xp.begin());
  std::reverse_copy(x.begin(),x.end(),xp.begin()+m);
  std::reverse_copy(y.begin(),y.end(),yp.begin()+m);

  int n = m * 2;

  // Creating index vectors..
  IntegerVector index1 = Range(0,n-2);
  IntegerVector index2 = Range(1,n-1);

  // Notice the vectorized sum function...
  double p1 = sum(xp[index1] * yp[index2]) + xp[n-1] * yp[0];
  double p2 = sum(xp[index2] * yp[index1]) + xp[0] * yp[n-1];

  double out = 0.5 * (p1 - p2);

  return(out);}')



TrapRcpp(x,y)

## [1] -1.268063

pracma::trapz(x,y)

## [1] -1.268063

microbenchmark::microbenchmark(pracma::trapz(x,y))

## Unit: microseconds
##                 expr   min     lq     mean median      uq    max neval
##  pracma::trapz(x, y) 9.846 10.211 10.87884 10.576 10.9405 43.397   100

microbenchmark::microbenchmark(TrapRcpp(x,y))

## Unit: microseconds
##            expr   min    lq    mean median    uq     max neval
##  TrapRcpp(x, y) 1.459 1.824 9.96352  1.824 1.825 806.292   100

How accurate is our numerical approximation? To find out, we'll calculate area again this time using the indefinite integral. Since the function in question, \(cos\), is trivial to integrate or reverse differentiate, so to is the process to find the exact area...

# Exact
sin(y[10]) - sin(y[1])

## [1] -1.258418

# Error
TrapRcpp(x,y) - (sin(y[10]) - sin(y[1]))

## [1] -0.009644797

What if we want a more accurate approximation? One way is to take more points, another is to use a slightly different method like Simpson's Rule Adaptive Simpsons' Method is often the workhorses of numerical computing packages... , but instead, I want to cover Gaussian Quadrature using Rcpp Armadillo.

RcppArmadillo

Gaussian Quadrature

Whereas other methods, like the trapezoidal rule, use evenly spaced intervals, Gaussian Quadrature instead allows both the intervals, or nodes, and the weights to vary. The idea, then, is to find both the weights and nodes which maximize the degree of exactness of the quadrature rule The term quadrature means finding the area of a geometric shape... , \(\int_{a}^bf(x)w(x)dx \approx \sum_{i=0}^n w_if(x_i)\) where \(w_i = \int_{a}^b\ell_i(x)w(x)\), which will be exact for all polynomials degree \(2n-1\) or less. Note that the significance of \(2n-1\) is that there are \(n\) nodes and \(n\) weights.

How is this exact? Full exposition is beyond the scope (and my ability) here, but the answer relates to a series of recurring orthogonal polynomials "The three-term recurrence relation satisfied by orthogonal polynomials is arguably the single most important piece of information for the constructive and computational use of orthogonal polynomials. Apart from its obvious use in generating values of orthogonal polynomials and their derivatives, both within and without the spectrum of the measure, knowledge of the recursion coefficients (i) allows the zeros of orthogonal polynomials to be readily computed as eigenvalues of a symmetric tridiagonal matrix, and with them the all-important Gaussian quadrature rule, (ii) yields immediately the normalization coefficients needed to pass from monic orthogonal to orthonormal polynomials, (iii) opens access to polynomials of the second kind and related continued fractions, and (iv) allows an efficient evaluation of expansions in orthogonal polynomials."
Orthogonal Polynomials - Gautschi (2004) , \(\int_{a}^b w(x) p_n (x) p_k(x) \propto \delta_{nk}\), that can be used to find to the roots of the polynomial, \(P_N (x_i) = 0\), which happen to be the optimal nodes for Gaussian Quadrature. Accordingly, the quadrature sum, \(\sum_{i=1}^N q(x_i)p_N(x_i) = 0\), gives the exact value for the integral, \(\int_{a}^b w(x)q(x)p_N(x)dx=0\), where, \( q(x)\lt N \).

There are a number of different polynomials used in Gaussian Quadrature (e.g., Hermite, LaGuerre, etc.), but the one we are interested in is the Legendre polynomials characterized by a) the first term is a monic polynomial (e.g, \(x^3 = 1 \cdot x^3\)) and b) for those polynomials with degree less than \(n\), the following holds, \(\int_{-1}^{1}P_{n}(x) \cdot P(x)dx = 0 \).

Notably, the roots of these polynomials are distinct, are all within the interval \([-1,1]\), and as noted previously, are ultimately used to represent the nodes for Gauss Legendre Quadrature See this handout.... Although the first several Legendre polynomials and roots are generally given in texts introducing the subject, \(P_{0}(x) = 1, P_{1}(x) = x, P_{2}(x) = X^2 - \frac{1}{3}, \dots\), modern algorithms don't rely on pre-computed values, but instead solve a tri-diagonal eigenvalue problem for the Jacobi matrix, \[ \textbf{J} = \left[ \begin{array}{ccccc} \alpha_0 & \beta_1 & & & \\ \beta_1 & \alpha_1 & \beta_2 & & \\ & \ddots & \ddots & \ddots & \\ & & \beta_{n-2} & \alpha_{n-2} & \beta_{n-1} \\ & & & \beta_{n-1} & \alpha_{n-1} \end{array} \right] \] , where the resulting eigenvalues are the nodes for evaluation and the first components of the eigenvectors which are used to create the weights after squaring and normalizing. For the special case of the Gauss Legendre Golub-Welsch Algorithm Golub-Welsch Algorithm, 1967 , \(\alpha = 0\) and \(\beta_k = \frac{k}{\sqrt(2k-1)(2k-1)}\). Also because limiting ourselves to the interval \([-1,1]\), isn't that practical, it is necessary to scale the interval as follows, \[ \int_{a}^{b} f(x)dx = \frac{b-a}{2} \int_{-1}^{1}f (\frac{b-a}{2}x + \frac{a+b}{2}x)dx \] and then, \[ \int_{a}^{b} f(x)dx = \frac{b-a}{2} \sum_{i=1}^n w_i f(\frac{b-a}{2}x + \frac{a+b}{2}x)dx \] The material above covered a lot of ground quickly, but code from, pracma::gaussLegendre() should help clarify, or, at least procedurally...

pracma::gaussLegendre

## function (n, a, b)
## {
##     stopifnot(is.numeric(a), length(a) == 1, is.numeric(b), length(b) ==
##         1, is.numeric(n), length(n) == 1, n >= 2)
##     i <- seq(1, n - 1, by = 1)
##     d <- i/sqrt(4 * i^2 - 1)
##     E <- eigen(Diag(d, 1) + Diag(d, -1), symmetric = TRUE)
##     L <- E$values
##     V <- E$vectors
##     inds <- order(L)
##     x <- L[inds]
##     V <- t(V[, inds])
##     w <- 2 * V[, 1]^2
##     x <- 0.5 * ((b - a) * x + a + b)
##     w <- -0.5 * (a - b) * w
##     return(list(x = x, w = w))
## }
## <bytecode: 0x0000000014725538>
## <environment: namespace:pracma>

Now, moving this to Rcpp Armadillo is fairly straightforward. Armadillo purposefully uses Matlab-like syntax Armadillo Syntax Table to make translation simple. Given that pracma also heavily uses Matlab function names, our job is even easier...

GLsrc <- 'List GLRcpp(int n, double a, double b){

// You might see const references in production code often...

const arma::vec& i = arma::linspace(1, n-1, n-1);

arma::mat J(n,n,arma::fill::zeros);

// Beta diagonal from Golub-Welsch

const arma::vec& d = i / (arma::sqrt(4 * arma::square(i) - 1));

// Setting off-diagonal elements

J.diag(1) = d;
J.diag(-1) = d;

// Initialize matrix and vector for eigenvalues/vectors

arma::vec L;
arma::mat V;

arma::eig_sym(L, V, J);

// Only need first row...

arma::vec w = 2 * arma::vectorise(arma::square(V.row(0)));

arma::vec x = .5 * ((b - a) * L + a + b);
w = -.5 * (a - b) * w;

return List::create(
Named("x") = x,
Named("w") = w);

}

'

Rcpp::cppFunction(GLsrc, depends="RcppArmadillo")

head(pracma::gaussLegendre(100,y[1],y[10])$x)

## [1] 0.5401049 0.5392622 0.5377471 0.5355607 0.5327050 0.5291829

head(pracma::gaussLegendre(100,y[1],y[10])$w)

## [1] -0.0005066678 -0.0011789457 -0.0018510682 -0.0025214686 -0.0031894227
## [6] -0.0038542657

head(GLRcpp(100,y[1],y[10])$x)

##           [,1]
## [1,] 0.5401049
## [2,] 0.5392622
## [3,] 0.5377471
## [4,] 0.5355607
## [5,] 0.5327050
## [6,] 0.5291829

head(GLRcpp(100,y[1],y[10])$w)

##               [,1]
## [1,] -0.0005066678
## [2,] -0.0011789457
## [3,] -0.0018510682
## [4,] -0.0025214686
## [5,] -0.0031894227
## [6,] -0.0038542657

# Computing area
gl <- GLRcpp(100,y[1],y[10])
sum(gl$w * cos(gl$x))

## [1] -1.258418

# Speed

microbenchmark::microbenchmark(pracma::gaussLegendre(100,y[1],y[10]))

## Unit: milliseconds
##                                     expr      min       lq     mean
##  pracma::gaussLegendre(100, y[1], y[10]) 1.198679 1.297506 1.844624
##    median       uq      max neval
##  1.395603 1.903227 26.24548   100

microbenchmark::microbenchmark(GLRcpp(100,y[1],y[10]))

## Unit: milliseconds
##                      expr      min       lq     mean   median       uq
##  GLRcpp(100, y[1], y[10]) 1.624981 1.695911 1.829103 1.762099 1.861107
##      max neval
##  2.56766   100

Somewhat to my surprise, there isn't much improvement in terms of speed I profiled a number of different parts of the algorithm, but couldn't pinpoint much, but undoubtedly there is a better way to do it... . However, considering the function mostly relies on R primatives, doesn't utilize a for-loop, and uses the highly optimized LAPACK eigenvalue decomposition, perhaps this isn't that surprising...

Cubic Splines

To finish up, let's briefly look at a translation of a simple natural cubic spline algorithm to Rcpp Armadillo. As the name implies, a (interpolating) cubic spline fits a cubic polynomial between each interval.

Ultimately, we want to end up with the following... \[S(x) = \begin{cases} S_{0}\left( x\right) =a_{0}+b_{0}\left( x-x_{i}\right) +c_{0}\left( x-x_{0}\right) ^{2}+d_{0}\left( x-x_{0}\right) ^{3} & \text{if }x_{0}\leq x\leq x_{1} \\ S_{1}\left( x\right) =a_{1}+b_{1}\left( x-x_{1}\right) +c_{1}\left( x-x_{1}\right) ^{2}+d_{1}\left( x-x_{1}\right) ^{3} & \text{if }x_{1}\leq x\leq x_{2} \\ : & : \\ S_{n-1}\left( x\right) =a_{n-1}+b_{n-1}\left( x-x_{n-1}\right) +c_{n-1}\left( x-x_{n-1}\right) ^{2}+d_{n-1}\left( x-x_{n-1}\right) ^{3} & \text{if }x_{n-1}\leq x\leq x_{n}\\ \end{cases} \]

To achieve this, we impose the following conditions,

• \(S_{i}(x_{i}) = f(x_{i})\) and \(S_{i+1}(x_{i+1}) = f(x_{i+1})\) for each \(i=0,1,\dots,n-1\). Or, in other words, the spline goes through each observation and is continuous for all interior points.
• Moreover, the first derivatives of all interior points are equal, (e.g., rate of change is the same), \(S'_{i+1}(x_{i+1}) = S'_{i}(x_{i+1})\) for each \(i=0,1,\dots,n-2\)
• Likewise, so are the second derivatives at all interior points, \(S''_{i+1}(x_{i+1}) = S''_{i}(x_{i+1})\) for each \(i=0,1,\dots,n-2\)

And finally, that either the second derivatives at end points are set to zero (natural), \(S''(x_{0}) = S''(x_{n}) = 0\), or that the first derivatives at the end points are fixed to a particular value (clamped), \(S'(x_{0}) = S'(x_{n})\) and \(S'(x_{n}) = f'(x_{n})\).

Here, we'll go with common requirement that the second derivatives at the endpoints are set to zero, which is a natural cubic spline. Everything else is just constructing the system of equations. In particular, analytical integration and alegraic manipulation A short, but excellent exposition gives us the following, \[h_i = x_{i+1} - x_i, i = 0,1,\dots,n-1\] \[ S= \left[ \begin{array}{ccccccccc} 1 & 0 & 0 & 0 & 0 & ... & 0 & 0 & 0 \\ h_{0} & 2( h_{0}+h_{1}) & h_{1} & 0 & 0 & ... & 0 & 0 & 0 \\ 0 & h_{1} & 2( h_{1}+h_{2}) & h_{2} & 0 & ... & 0 & 0 & 0 \\ : & : & : & : & : & ... & 0 & : & : \\ 0 & 0 & 0 & 0 & 0 & ... & h_{n-2} & 2( h_{n-2}+h_{n-1}) & h_{n-1} \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & 0 & 1 \end{array}\right] \] \[ \vec{v}=\left[ \begin{array}{l} 0 \\ 3\left[ \dfrac{1}{h_{1}}\left( f(x_{2}) - f(x_{1})\right) -\dfrac{1}{h_{0}}\left( f(x_{1}) - f(x_{0})\right) \right] \\ 3\left[ \dfrac{1}{h_{2}}\left( f(x_{3}) - f(x_{2})\right) -\dfrac{1}{h_{1}}\left( f(x_{2}) - f(x_{1})\right) \right] \\ : \\ 3\left[ \dfrac{1}{h_{n-1}}\left( f(x_{n})-f(x_{n-1})\right) -\dfrac{1}{h_{n-2}}% \left( f(x_{n-1}) - f(x_{n-2})\right) \right] \\ 0% \end{array}% \right] \] Accordingly, one then has a system of equations where one solves \(A\vec{c}=\vec{v}\) and then after obtaining a solution for \(c\) Note that this form is slightly different than a lot of other expositions where one solves for the second derivatives as opposed to the first... , finds both \(b\) and \(d\), \[\begin{aligned} b_{i} &=( a_{i+1}-a_{i}) - ( 2c_{i}+c_{i+1}),\ i=0,1,\dots,n-1 \\ d_{i} &=( c_{i+1}-c_{i}),\ i=0,1,\dots,n-1 \end{aligned}\] The pracma version gives the option to specify the derivatives at the endpoints and/or interpolate a value. Also, near the end, there class construction code consistent with Matlab.

pracma::cubicspline

## function (x, y, xi = NULL, endp2nd = FALSE, der = c(0, 0))
## {
##     n <- length(x)
##     h <- x[2:n] - x[1:(n - 1)]
##     e <- 2 * c(h[1], h[1:(n - 2)] + h[2:(n - 1)], h[n - 1])
##     A <- Diag(e) + Diag(h, -1) + Diag(h, 1)
##     d <- (y[2:n] - y[1:(n - 1)])/h
##     rhs <- 3 * (d[2:(n - 1)] - d[1:(n - 2)])
##     der0 <- der[1]
##     dern <- der[2]
##     if (endp2nd) {
##         A[1, 1] <- 2 * h[1]
##         A[1, 2] <- h[1]
##         A[n, n] <- 2 * h[n - 1]
##         A[n - 1, n - 2] <- h[n - 1]
##         rhs <- c(3 * (d[1] - der0), rhs, 3 * (dern - d[n - 1]))
##     }
##     else {
##         A[1, ] <- 0
##         A[1, 1] <- 1
##         A[n, ] <- 0
##         A[n, n] <- 1
##         rhs <- c(der0, rhs, dern)
##     }
##     S <- zeros(n, 4)
##     S[, 3] <- solve(A, rhs)
##     for (m in 1:(n - 1)) {
##         S[m, 4] = (S[m + 1, 3] - S[m, 3])/3/h[m]
##         S[m, 2] = d[m] - h[m]/3 * (S[m + 1, 3] + 2 * S[m, 3])
##         S[m, 1] = y[m]
##     }
##     S <- S[1:(n - 1), 4:1]
##     pp <- mkpp(x, S)
##     if (is.null(xi)) {
##         return(pp)
##     }
##     else {
##         yi <- ppval(pp, xi)
##         return(yi)
##     }
## }
## <bytecode: 0x0000000014940b10>
## <environment: namespace:pracma>

I'm going to simplify things by getting rid of everything else except the code which gives us the spline coefficients.

cubicsrc <- ' List CubicSplineRcpp(arma::vec x,arma::vec y) {

int n = x.n_elem;
arma::vec h = x.subvec(1, n-1) - x.subvec(0, n-2);
arma::vec e(n);

e.head(1) = h.head(1); e.tail(1) = h.tail(1);
e.subvec(1,n-2) = h.subvec(0,n-3) + h.subvec(1,n-2);
e = 2 * e;

arma::mat A = arma::diagmat(e) + arma::diagmat(h,-1) + arma::diagmat(h,1);
arma::vec d = (y.subvec(1,n-1) - y.subvec(0,n-2)) / h ;

arma::vec rhs(n);
rhs.subvec(1,n-2) = 3 * (d.subvec(1,n-2) - d.subvec(0,n-3));
A.row(0).zeros(); A.row(n-1).zeros();
A.at(0,0) = 1; A.at(n-1,n-1) = 1; // Setting up matrix...

arma::mat  S(n,4);
S.col(2) = solve(A, rhs); // solve for derivatives

// Now, go back and find other terms

S.col(3).rows(0,n-2) = (S.col(2).rows(1,n-1) - S.col(2).rows(0,n-2)) / 3;
arma::vec tmp1 = h / 3; // Removed loop and Vectorized
arma::vec tmp2 = S.col(2).rows(1,n-1) + 2 * S.col(2).rows(0,n-2);
arma::vec tmp3 = tmp1 % tmp2;

S.col(1).rows(0,n-2) = d - tmp3;
S.col(0).rows(0,n-2) = y.subvec(0,n-2);

S.shed_row(n-1); S = fliplr(S);

return List::create(
Named("S") = S);}'

Rcpp::cppFunction(cubicsrc,depends="RcppArmadillo")

cubicspline(x,y)$coefs

##              [,1]        [,2]       [,3]       [,4]
##  [1,]  0.04910445  0.00000000 -1.0055536  0.5403023
##  [2,]  0.13708122  0.14731336 -0.8582402 -0.4161468
##  [3,] -0.06983828  0.55855702 -0.1523699 -0.9899925
##  [4,] -0.16696569  0.34904217  0.7552293 -0.6536436
##  [5,] -0.12405359 -0.15185490  0.9524166  0.2836622
##  [6,]  0.04120163 -0.52401568  0.2765460  0.9601703
##  [7,]  0.14888895 -0.40041079 -0.6478805  0.7539023
##  [8,]  0.19014888  0.04625607 -1.0020352 -0.1455000
##  [9,] -0.20556757  0.61670270 -0.3390764 -0.9111303

CubicSplineRcpp(x,y)

## $S
##              [,1]        [,2]       [,3]       [,4]
##  [1,]  0.04910445  0.00000000 -1.0055536  0.5403023
##  [2,]  0.13708122  0.14731336 -0.8582402 -0.4161468
##  [3,] -0.06983828  0.55855702 -0.1523699 -0.9899925
##  [4,] -0.16696569  0.34904217  0.7552293 -0.6536436
##  [5,] -0.12405359 -0.15185490  0.9524166  0.2836622
##  [6,]  0.04120163 -0.52401568  0.2765460  0.9601703
##  [7,]  0.14888895 -0.40041079 -0.6478805  0.7539023
##  [8,]  0.19014888  0.04625607 -1.0020352 -0.1455000
##  [9,] -0.20556757  0.61670270 -0.3390764 -0.9111303

For this function, there is definitely a speed-up with RcppArmadillo.

microbenchmark::microbenchmark(cubicspline(x,y)$coefs)

## Unit: microseconds
##                     expr    min     lq     mean median     uq     max
##  cubicspline(x, y)$coefs 68.924 70.383 76.85909 71.841 75.123 175.773
##  neval
##    100

microbenchmark::microbenchmark(CubicSplineRcpp(x,y))

## Unit: microseconds
##                   expr   min    lq     mean median     uq     max neval
##  CubicSplineRcpp(x, y) 9.117 9.482 17.84412  9.482 9.6645 831.454   100

Now, go forth and animate or something Cubic splines are frequently used in computer animation/graphics... !

Links

To close out, I'll provide links to a number of resources related to what I covered,

Rcpp Links

• Rcpp Official Documentation
• Rcpp Book (unaffiliated)
• Rcpp Gallery
• General Rcpp Overview
• Eddelbuettel Rcpp Slides
• Advanced R: Rcpp Chapter
• Rcpp FAQ
• Quick Reference Guide

RcppArmadillo

• RcppArmadillo Page
• RcppArmadillo
• Armadillo Cpp Homepage
• Armadillo 2010 Paper

Polynomial Interpolation

• Basic Polynomial Handout
• LaGrange Polynomial Handout 1
• LaGrange Polynomial Handout 2
• LaGrange Polynomial Handout 3

Natural Cubic Splines

• Cubic Spline Handout 1 (helpful!)
• Cubic Spline Handout 2
• Cubic Spline Handout 3
• Cubic Spline Handout 4

Gaussian Quadrature

• Golub & Welsch, 1967
• Gauss. Quad. Handout 1
• Gauss. Quad. Handout 2
• Gauss. Quad. Handout 3
• Gauss. Quad. Handout 4
• Gauss. Quad. Handout 5
• Gauss. Quad. Handout 6

Appendix

Cubic Spline Construction

The necessary cubic spline equations I rely on Cubic Spline Handout 1 (above) and Numerical Analysis (Burden & Faires, 2011) here... can be derived directly from the required conditions.

First, we know that the first polynomial term, \(a_i\), must equal \(f(x_i)\), due to the fact \(S(x_i) = f(x_i)\). Thus, \[f(x_i) = a_i = S_i(x_i)\] Moreover, because the spline should interpolate the next consecutive point, \(S_{n-1}(x_n) = y_n\), it leads to, \[f(x_n) = f(x_{n-1}) + b_{n-1}(x_n - x_{n-1}) + c_{n-1}(x_n - x_{n-1})^2 + d_{n-1}(x_n - x_{n-1})^3 \] Based on the requirement that \(S_i(x_{i+1}) = S_{i+1}(x_{i+1}) \), or continuity for all interior points, there are \(n-1\) equations of the following form, \[f(x_{i+1}) = f(x_i) + b_{i}(x_{i+1} - x_{i}) + c_{i}(x_{i+1} - x_{i})^2 + d_{i}(x_{i+1} - x_{i})^3 \] Differentiating the spline equation once results in, \(S_i'(x) = b_i + 2c_i(x-x_i) + 3d_i(x-x_i)^2\), which along with the requirement, \(S_i'(x_{i+1}) = S_{i+1}'(x_{i+1})\), gives another \(n-1\) of, \[b_{i+1} = b_i + 2c_i(x-x_i) + 3d_i(x-x_i)^2\] The requirement that, \(S_i''(x_{i+1}) = S_{i+1}''(x_{i+1})\), and differentiating again then gives \(n-1\) equations of, \[2c_{i+1} = 2c_i + 6d_i(x_{i+1}-x_i)\] Finally, the natural spline requires \(S_0''(x_0) = 0, \ S_{n-1}''(x_n) = 0\), so it follows that, \[ \begin{aligned} &2c_0 = 0 \\ &6d_{n-1}(x_n-x_{n-1}) + 2c_{n-1} = 0 \end{aligned} \]