Polynomials & Integration
A (slightly less) Naive Introduction to Algorithms: Part III
BB
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
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} \]
Comments
2019-07-11 23:42:00
http://mewkid.net/buy-amoxicillin/ - Amoxicillin Without Prescription <a href="http://mewkid.net/buy-amoxicillin/">Buy Amoxicillin</a> dmz.mixm.raw-r.org.gjw.ji http://mewkid.net/buy-amoxicillin/
2019-07-12 00:04:00
http://mewkid.net/buy-amoxicillin/ - Amoxicillin - Prix.achetercommander.fr <a href="http://mewkid.net/buy-amoxicillin/">Buy Amoxicillin Online</a> zww.hdxv.raw-r.org.fdg.zr http://mewkid.net/buy-amoxicillin/
2020-06-16 11:02:00
Online Pharmacy No Prescription Canada https://cheapcialisll.com/ - where to buy cialis cialis prostate surgery <a href=https://cheapcialisll.com/#>buy cialis online india</a> Achat Viagra 50mg Telephone
2020-07-06 14:20:00
Amoxil Acide Clavulanique SeneFewWeesk https://ascialis.com/# - Cialis Autown Acheter Xenical Sans Ordonnance sahshibers <a href=https://ascialis.com/#>Cialis</a> Ralgaith Ciprofloxacin Alternative
2020-08-25 01:59:00
Healthyman.Co SeneFewWeesk https://asocialiser.com/ - tadalafil cialis from india Autown Viagra Generique Prix sahshibers <a href=https://asocialiser.com/#>cialis 20mg</a> Ralgaith Propecia En Hombres
2020-09-19 14:01:00
cialis 10 mg 4 comprimidos precio SeneFewWeesk https://artsocialist.com/ - Cialis Autown Buy Generic Viagra With Pay Pal sahshibers <a href=https://artsocialist.com/#>Cialis</a> Ralgaith Levitra Generika Sicher Kaufen
2020-12-08 07:34:00
Propecia Comparativa SeneFewWeesk <a href=https://xbuycheapcialiss.com/>buying cialis generic</a> Autown What Do You Use Amoxicillin For
2021-01-10 00:38:00
Pills information leaflet. Generic Name. <a href="https://viagra4u.top">can i purchase generic viagra for sale</a> in Canada. Some information about medicine. Get now.
<a href=https://amp.en.vaskar.co.in/translate/1?to=ru&from=en&source=Medicament%20information%20leaflet.%20Brand%20names.%20%3Ca%20href%3D%22https%3A%2F%2Fprednisone4u.top%22%3Ecan%20you%20get%20prednisone%20tablets%3C%2Fa%3E%20in%20Canada.%20Some%20about%20meds.%20Read%20here.%20%0D%0A%3Ca%20href%3Dhttp%3A%2F%2Feshop.goodmat.cz%2Findex.php%3Faction%3Ddetail%26id%3DID1482%3EAll%20trends%20of%20meds.%3C%2Fa%3E%20%3Ca%20href%3Dhttp%3A%2F%2Fkireevsk-crb-zdrav.ru%2Fforums%2Ftopic%2Feverything-about-medicine%2F%3EEverything%20about%20medicine.%3C%2Fa%3E%20%3Ca%20href%3Dhttp%3A%2F%2Fjapaninfo.jp%2Frecuruitforum%2F%3Fmod%3Ddocument%26uid%3D136%23kboard-comments-136%3EEverything%20trends%20of%20medicament.%3C%2Fa%3E%20%204f703ca%20&result=%D0%98%D0%BD%D1%84%D0%BE%D1%80%D0%BC%D0%B0%D1%86%D0%B8%D0%BE%D0%BD%D0%BD%D1%8B%D0%B9%20%D0%BB%D0%B8%D1%81%D1%82%D0%BE%D0%BA%20%D0%BE%20%D0%BB%D0%B5%D0%BA%D0%B0%D1%80%D1%81%D1%82%D0%B2%D0%B0%D1%85.%20%D0%A4%D0%B8%D1%80%D0%BC%D0%B5%D0%BD%D0%BD%D1%8B%D0%B5%20%D0%BD%D0%B0%D0%B8%D0%BC%D0%B5%D0%BD%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D1%8F.%20%3Ca%20href%3D%22https%3A%2F%2Fprednisone4u.top%22%20%3E%20%D0%B2%D1%8B%20%D0%BC%D0%BE%D0%B6%D0%B5%D1%82%D0%B5%20%D0%BF%D0%BE%D0%BB%D1%83%D1%87%D0%B8%D1%82%D1%8C%20%D1%82%D0%B0%D0%B1%D0%BB%D0%B5%D1%82%D0%BA%D0%B8%20%D0%BF%D1%80%D0%B5%D0%B4%D0%BD%D0%B8%D0%B7%D0%BE%D0%BD%D0%B0%3C%2Fa%3E%20%D0%B2%20%D0%9A%D0%B0%D0%BD%D0%B0%D0%B4%D0%B5.%20%D0%9A%D0%BE%D0%B5-%D1%87%D1%82%D0%BE%20%D0%BE%20%D0%BB%D0%B5%D0%BA%D0%B0%D1%80%D1%81%D1%82%D0%B2%D0%B0%D1%85.%20%D0%A7%D0%B8%D1%82%D0%B0%D0%B9%D1%82%D0%B5%20%D0%B7%D0%B4%D0%B5%D1%81%D1%8C.%20%3Ca%20href%3Dhttp%3A%2F%20%2F%20eshop.goodmat.cz%20%2F%20index.php%3Faction%3Ddetail%26id%3DID1482%3E%D0%B2%D1%81%D0%B5%20%D1%82%D0%B5%D0%BD%D0%B4%D0%B5%D0%BD%D1%86%D0%B8%D0%B8%20%D0%BB%D0%B5%D0%BA%D0%B0%D1%80%D1%81%D1%82%D0%B2.%3C%20%2F%20a%3E%20%3Ca%20href%3Dhttp%3A%2F%2Fkireevsk-crb-zdrav.ru%2F%D1%84%D0%BE%D1%80%D1%83%D0%BC%D1%8B%2F%D1%82%D0%B5%D0%BC%D0%B0%20%2F%20%D0%B2%D1%81%D0%B5-%D0%BE-%D0%BC%D0%B5%D0%B4%D0%B8%D1%86%D0%B8%D0%BD%D0%B5%20%2F%20%3E%D0%B2%D1%81%D0%B5%20%D0%BE%20%D0%BC%D0%B5%D0%B4%D0%B8%D1%86%D0%B8%D0%BD%D0%B5.%3C%20%2F%20a%3E%20%3Ca%20href%3Dhttp%3A%2F%2Fjapaninfo.jp%20%2F%20recuruitforum%2F%3F%20mod%3Ddocument%26uid%3D136%23kboard-comments-136%3E%D0%B2%D1%81%D0%B5%20%D1%82%D0%B5%D0%BD%D0%B4%D0%B5%D0%BD%D1%86%D0%B8%D0%B8%20%D0%BC%D0%B5%D0%B4%D0%B8%D1%86%D0%B8%D0%BD%D1%8B.%3C%20%2F%20a%3E%204f703ca>Best trends of medication.</a> <a href=https://dinhduongnet.com/san-pham/sua-anka-well-iq1-400gram/#comment-67466>Best information about drugs.</a> <a href=http://bpo.gov.mn/content/185>Everything news about drug.</a> 8_791c8
2021-01-10 16:59:00
Medicines information for patients. What side effects? <a href="https://prednisone4u.top">can you buy cheap prednisone for sale</a> in USA. Best trends of meds. Get now.
<a href=http://shop.khunjib.com/product/1540/%E0%B9%80%E0%B8%8B%E0%B8%95%E0%B8%AA%E0%B8%9B%E0%B8%AD%E0%B8%A3%E0%B9%8C%E0%B8%95%E0%B9%80%E0%B8%81%E0%B8%B4%E0%B8%A5%E0%B9%8C%E0%B8%A5%E0%B9%81%E0%B8%95%E0%B9%88%E0%B8%87%E0%B9%81%E0%B8%96%E0%B8%9A%E0%B8%82%E0%B9%89%E0%B8%B2%E0%B8%87%20%E0%B8%AA%E0%B8%B5%E0%B8%94%E0%B8%B3.html>Everything news about meds.</a> <a href=https://actcycle.jp/spot/r459.html#comment-60029>Everything news about meds.</a> <a href=https://libertyarmenia.am/best-religious-sights-in-armenia/#comment-32>All what you want to know about medicines.</a> 9d6c91a