calculus: High Dimensional Numerical and Symbolic Calculus in R
ccalculus: High Dimensional Numerical and SymbolicCalculus in R Emanuele Guidotti
University of Neuchˆatel
Abstract
The R package calculus implements C++ optimized functions for numerical and sym-bolic calculus, such as the Einstein summing convention, fast computation of the Levi-Civita symbol and generalized Kronecker delta, Taylor series expansion, multivariate Her-mite polynomials, high-order derivatives, ordinary differential equations, differential oper-ators and numerical integration in arbitrary orthogonal coordinate systems. The libraryapplies numerical methods when working with functions or symbolic programming whenworking with characters or expressions . The package handles multivariate numericalcalculus in arbitrary dimensions and coordinates and implements the symbolic counterpartof the numerical methods whenever possible, without depending on external computer al-gebra systems. Except for Rcpp , the package has no strict dependencies in order toprovide a stable self-contained toolbox that invites re-use.
Keywords : symbolic programming, finite difference, differential operators, numerical integra-tion, coordinate systems, Einstein summation, Taylor series, Hermite polynomials, R .
1. Introduction
Multivariate calculus underlies a wide range of applications in the natural and social sci-ences. In statistics, asymptotic expansion formulas for stochastic processes (Yoshida 1992)can be obtained by solving high dimensional systems of ordinary differential equations. Thetransition density of multivariate diffusions can be approximated using Hermite polynomials(A¨ıt-Sahalia 2002) or Taylor-like expansions (Li et al. R (R Core Team 2020) has shown to be a viable computing environment for implementingand applying numerical methods as a practical tool for applied statistics. However, suchmethods are seldom flexible enough to handle multivariate calculus in arbitrary dimensionsand coordinates. The package numDeriv (Gilbert and Varadhan 2019) sets the standard for https://cloud.r-project.org/web/views/NumericalMathematics.html a r X i v : . [ c s . M S ] D ec calculus : High Dimensional Numerical and Symbolic Calculus in R numerical differentiation in R , providing numerical gradients, Jacobians, and Hessians, butdoes not support higher order derivatives or differentiation of tensor-valued functions. tensorA (van den Boogaart 2020) implements the Einstein summing convention but does not supportarbitrary expressions involving more than two tensors or tensors with repeated indices. mpoly (Kahle 2013) implements univariate but not multivariate Hermite polynomials. In a similarway, pracma (Borchers 2019) supports the computation of Taylor series for univariate butnot multivariate functions. cubature (Narasimhan, Johnson, Hahn, Bouvier, and Kiˆeu 2019)provides an efficient interface for multivariate integration but limited to Cartesian coordinates.On the other hand, R is not designed for symbolic computing. Nevertheless, the advent ofalgebraic statistics and its contributions to asymptotic theory in statistical models, experi-mental design, multiway contingency tables, and disclosure limitation has increased the needfor R to be able to do some relatively basic operations and routines with multivariate symboliccalculus (Kahle 2013). Although there exist packages to interface external computer algebrasystems, R still lacks a native support that invites re-use. The package Ryacas (Andersenand Hojsgaard 2019) interfaces the computer algebra system
Yacas , while caracas (Ander-sen and Hojsgaard 2020) - based on reticulate (Ushey, Allaire, and Tang 2020) - and rSymPy (Grothendieck and Bellosta 2019) - based on rJython (Grothendieck and Bellosta 2012) - bothaccess the symbolic algebra system SymPy .This work presents the R package calculus for high dimensional numerical and symbolic cal-culus in R . The contribution is twofold. First, the package handles multivariate numericalcalculus in arbitrary dimensions and coordinates via C++ optimized functions, improving thestate-of-the-art both in terms of flexibility and efficiency. It achieves approximately the sameaccuracy for numerical differentiation as the numDeriv (Gilbert and Varadhan 2019) packagebut significantly reduces the computational time. It supports higher order derivatives and thedifferentiation of possibly tensor-valued functions. Differential operators such as the gradient,divergence, curl, and Laplacian are made available in arbitrary orthogonal coordinate systems.The Einstein summing convention supports expressions involving more than two tensors andtensors with repeated indices. Besides being more flexible, the summation proves to be fasterthan the alternative implementation found in the tensorA (van den Boogaart 2020) packagefor advanced tensor arithmetic with named indices. Unlike mpoly (Kahle 2013) and pracma (Borchers 2019), the package supports multidimensional Hermite polynomials and Taylor se-ries of multivariate functions. The package integrates seamlessly with cubature (Narasimhan et al. C and extends the numerical integration toarbitrary orthogonal coordinate systems. Second, the symbolic counterpart of the numericalmethods are implemented whenever possible to meet the growing needs for R to handle basicsymbolic operations. The package provides, among others, symbolic high order derivativesof possibly tensor-valued functions, symbolic differential operators in arbitrary orthogonalcoordinate systems, symbolic Einstein summing convention and Taylor series expansion ofmultivariate functions. This is done entirely in R , without depending on external computeralgebra systems in order to provide a self-contained toolbox that invites re-use.The remainder of the paper is organized as follows: Section 2 introduces the package and theunderlying philosophy, Section 3 and 4 provide basic utilities for vector and matrix algebra,Section 5 presents tensor algebra with particular focus on the Einstein summation, Section6 provides fast and accurate derivatives, Section 7 presents the Taylor series of possibly manuele Guidotti
2. The R package calculus
The R package calculus implements C++ optimized functions for numerical and symboliccalculus, such as the Einstein summing convention, fast computation of the Levi-Civita symboland generalized Kronecker delta, Taylor series expansion, multivariate Hermite polynomials,high-order derivatives, ordinary differential equations, differential operators and numericalintegration in arbitrary orthogonal coordinate systems.
Several unit tests are implemented via the standard framework offered by testthat (Wickham2011) and run via continuous integration on Travis CI . The package integrates seamlessly with cubature (Narasimhan et al. C . However, except for Rcpp (Eddelbuettel and Fran¸cois 2011), thepackage has no strict dependencies in order to provide a stable self-contained toolbox thatinvites re-use.
The stable release version of calculus is hosted on the Comprehensive R Archive Network(CRAN) at https://CRAN.R-project.org/package=calculus and it can be installed using:
R> install.packages("calculus")
The package provides a unified interface to work with mathematical objects in R . The libraryapplies numerical methods when working with functions or symbolic programming whenworking with characters or expressions . To describe multidimensional objects such asvectors, matrices, and tensors, the package uses the class array regardless of the dimension.This is done to prevent unwanted results due to operations among different classes such as vector for unidimensional objects or matrix for bidimensional objects. Basic elementwise operations are supported for arrays of the same dimensions with the follow-ing data types: numeric , complex , character , expression . Automatic type conversion issupported and string manipulation is performed in C++ to improve performance. A minimal https://travis-ci.com calculus : High Dimensional Numerical and Symbolic Calculus in R simplification algorithm is included to simplify operations involving zeros. Below a unidi-mensional example on the sum, difference, product, and division among the different datatypes. R> ("a + b" %prod% 1i) %sum% (0 %prod% "c") %diff% (expression(d+e) %div% 3) [1] "((a + b) * (0+1i)) - ((d + e) / 3)"
The characters are automatically wrapped in parentheses when performing basic symbolicoperations to prevent unwanted results, e.g.: a + b · c + d instead of ( a + b ) · ( c + d )To disable this behaviour the user can set options(calculus.auto.wrap = FALSE) .
3. Vector algebra
A vector can be regarded as a 1-dimensional tensor. In R , it can be regarded as a 1-dimensional array so that the methods presented for multidimensional tensors in Section 5 are availablefor vectors. The package also implements a few vector-specific utilities, such as the crossproduct. The cross product or vector product is an operation on n − n -dimensional space.The results is a n -dimensional vector that is perpendicular to the n − R : R> cross(c(1,0,0), c(0,1,0)) [1] 0 0 1
And in R : R> cross(c(1,0,0,0), c(0,1,0,0), c(0,0,0,1)) [1] 0 0 1 0
Consistently with the philosophy of the package, the same interface is provided for character vectors.
R> cross(c("a","b","c"), c("d","e","f")) manuele Guidotti [1] "((b)*((f)) + -(e)*((c))) * 1" "((a)*((f)) + -(d)*((c))) * -1"[3] "((a)*((e)) + -(d)*((b))) * 1"
4. Matrix algebra
A matrix can be regarded as a 2-dimensional tensor. In R , it can be regarded as a 2-dimensional array so that the methods presented for multidimensional tensors in Section5 are available for matrices. The package also implements a few matrix-specific utilities, suchas the symbolic determinant, inverse, and matrix product. The function mxdet computes the numerical or symbolic determinant of matrices dependingon the data type. If the elements of the matrix are of type numeric , then the determinant iscomputed via the function det available in base R . R> mxdet(matrix(1:4, nrow = 2)) [1] -2
If the elements are of type character , then the symbolic determinant is computed recursivelyin
C++ . R> mxdet(matrix(letters[1:4], nrow = 2)) [1] "a*(d) + -b*(c)"
The symbolic determinant offers a significant gain in performance when computing deter-minants for a large number of matrices. The following test compares the performance oftwo different approaches to compute the determinant of 2 numeric :compute the numeric determinant for each matrix. Method symbolic : compute the symbolicdeterminant of a 4x4-matrix and evaluate it for each matrix. R> n <- 4R> e <- letters[1:n^2]R> grid <- expand.grid(lapply(1:n^2, function(e) runif(2)))R> colnames(grid) <- eR> microbenchmark(R+ "numeric" = {R+ x <- apply(grid, 1, function(e) det(matrix(e, nrow = n)))R+ },R+ "symbolic" = {R+ x <- evaluate(c2e(mxdet(matrix(e, nrow = n))), grid)R+ }R+ ) calculus : High Dimensional Numerical and Symbolic Calculus in R Unit: millisecondsexpr min lq mean median uq max neval cldnumeric 911.6423 935.7297 967.09 964.9941 993.513 1086.197 100 bsymbolic 7.4456 8.1825 10.44 9.2427 13.423 16.848 100 a
The function mxinv computes the numerical or symbolic inverse of matrices depending on thedata type. If the elements of the matrix are of type numeric , then the inverse is computedvia the function solve available in base R . R> mxinv(matrix(1:4, byrow = TRUE, nrow = 2)) [,1] [,2][1,] -2.0 1.0[2,] 1.5 -0.5
If the elements are of type character , then the symbolic inverse is computed based on thedeterminants in Section 4.1.
R> mxinv(matrix(letters[1:4], byrow = TRUE, nrow = 2)) [,1] [,2][1,] "(d) / (a*(d) + -c*(b))" "-(b) / (a*(d) + -c*(b))"[2,] "-(c) / (a*(d) + -c*(b))" "(a) / (a*(d) + -c*(b))"
The symbolic inverse offers a gain in performance when inverting a large number of matrices,as shown by replicating the test in section 4.1 and replacing the determinant with the inverse.
R> n <- 4R> e <- letters[1:n^2]R> grid <- expand.grid(lapply(1:n^2, function(e) runif(2)))R> colnames(grid) <- eR> microbenchmark(R+ "numeric" = {R+ x <- apply(grid, 1, function(e) solve(matrix(e, nrow = n)))R+ },R+ "symbolic" = {R+ x <- evaluate(c2e(mxinv(matrix(e, nrow = n))), grid)R+ }R+ )
Unit: millisecondsexpr min lq mean median uq max neval cldnumeric 1507.0 1625.25 1693.77 1702.14 1732.92 2066.72 100 bsymbolic 192.5 224.77 279.43 249.18 342.02 469.23 100 a manuele Guidotti The matrix product can be expressed in Einstein notation as shown in Section 5, thus inher-iting the support for symbolic calculations.
R> a <- matrix(1:4, nrow = 2, byrow = TRUE)R> b <- matrix(letters[1:4], nrow = 2, byrow = TRUE)R> a %mx% b [,1] [,2][1,] "1 * (a) + 2 * (c)" "1 * (b) + 2 * (d)"[2,] "3 * (a) + 4 * (c)" "3 * (b) + 4 * (d)"
5. Tensor algebra
A tensor may be represented as a multidimensional array. Just as a vector in an n -dimensionalspace is represented by a 1-dimensional array with n components, a matrix is represented bya 2-dimensional array with n x n components, and any tensor can be represented by a d -dimensional array with n x . . . x n d components. This makes the class array available inbase R an ideal candidate to represent mathematical tensors. In particular, the class stores itsdimensions in the attribute dim that contains a vector giving the length for each dimension. R> A <- array(1:24, dim = c(2,3,4))R> attributes(A) $dim[1] 2 3 4
The package calculus reads this attribute to represent tensors in index notation, such as A ijk .In particular, the function index is used to assign indices to the dimensions of the tensor bysetting names to the attribute dim . R> index(A) <- c("i","j","k")R> attributes(A) $dimi j k2 3 4
In this way, a tensor with named indices is represented by an array with named dim . Thepackage calculus builds a set of tools to work with tensors on top of this class and providesthe implementation of the Levi-Civita symbol and Generalized Kronecker delta that oftenappears in tensor algebra. At the time of writing, the package makes no distinction betweenupper and lower indices, i.e. vectors and covectors . Levi-Civita symbol https://en.wikipedia.org/wiki/Einstein_notation calculus : High Dimensional Numerical and Symbolic Calculus in R In mathematics, particularly in linear algebra, tensor analysis, and differential geometry, theLevi-Civita symbol represents a collection of numbers; defined from the sign of a permutationof the natural numbers 1 , , . . . , n , for some positive integer n . It is named after the Italianmathematician and physicist Tullio Levi-Civita. Other names include the permutation sym-bol, antisymmetric symbol, or alternating symbol, which refer to its antisymmetric propertyand definition in terms of permutations. In the general n -dimensional case, the Levi-Civitasymbol is defined by: ε i i ...i n = +1 if ( i , i , . . . , i n ) is an even permutation of (1 , , . . . , n ) − i , i , . . . , i n ) is an odd permutation of (1 , , . . . , n )0 otherwiseThe function epsilon determines the parity of the permutation in C++ via efficient cycledecomposition and constructs the Levi-Civita symbol in arbitrary dimension. For examplethe 2-dimensional Levi-Civita symbol is given by: ε ij = +1 if ( i, j ) = (1 , − i, j ) = (2 , i = j R> epsilon(2) [,1] [,2][1,] 0 1[2,] -1 0
And in 3 dimensions: ε ijk = +1 if ( i, j, k ) is (1 , , , (2 , , , or (3 , , , − i, j, k ) is (3 , , , (1 , , , or (2 , , , i = j, or j = k, or k = i R> epsilon(3) , , 1[,1] [,2] [,3][1,] 0 0 0[2,] 0 0 1[3,] 0 -1 0, , 2 https://en.wikipedia.org/wiki/Levi-Civita_symbol manuele Guidotti [,1] [,2] [,3][1,] 0 0 -1[2,] 0 0 0[3,] 1 0 0, , 3[,1] [,2] [,3][1,] 0 1 0[2,] -1 0 0[3,] 0 0 0 Generalized Kronecker delta
The generalized Kronecker delta or multi-index Kronecker delta of order 2 p is a type ( p, p )tensor that is a completely antisymmetric in its p upper indices, and also in its p lower indices. In terms of the indices, the generalized Kronecker delta is defined as (Frankel 2011): δ µ ...µ p ν ...ν p = +1 if ( ν . . . ν p ) is an even permutation of ( µ . . . µ p ) − ν . . . ν p ) is an odd permutation of ( µ . . . µ p )0 otherwiseWhen p = 1, the definition reduces to the standard Kronecker delta that corresponds to the n x n identity matrix I ij = δ ij where i and j take the values 1 , , . . . , n . R> delta(n = 3, p = 1) [,1] [,2] [,3][1,] 1 0 0[2,] 0 1 0[3,] 0 0 1
Tensor contraction can be seen as a generalization of the trace for a square matrix. In thegeneral case, a tensor can be contracted by summing over pairs of repeated indices that sharethe same dimension. This is achieved via the
C++ optimized function contraction . Considerthe following 2x2x2 tensor:
R> x <- array(1:8, dim = c(2,2,2))R> print(x) , , 1 https://en.wikipedia.org/wiki/Kronecker_delta calculus : High Dimensional Numerical and Symbolic Calculus in R [,1] [,2][1,] 1 3[2,] 2 4, , 2[,1] [,2][1,] 5 7[2,] 6 8 The trace of the tensor T = (cid:80) i T iii is obtained with: R> contraction(x) [1] 9
The contraction on the first and third dimension T j = (cid:80) i T iji can be computed with: R> index(x) <- c("i","j","i")R> contraction(x) [1] 7 11
Finally, it is possible to preserve the dummy dimensions T ij = T iji by setting the argument drop = FALSE : R> index(x) <- c("i","j","i")R> contraction(x, drop = FALSE) [,1] [,2][1,] 1 6[2,] 3 8
In this way, it is possible to compute arbitrary contraction of tensors such as T klm = (cid:80) ij T ikiiljjm or T ijklm = T ikiiljjm to preserve the dummy dimensions. In mathematics, the Einstein notation or Einstein summation convention is a notationalconvention that implies summation over a set of repeated indices. When an index variableappears twice, it implies summation over all the values of the index . For instance the matrixproduct can be written in terms of Einstein notation as: C ij = A ik B kj ≡ (cid:88) k A ik B kj https://en.wikipedia.org/wiki/Einstein_notation manuele Guidotti D k = A ijj B iijk C j ≡ (cid:88) ij A ijj B iijk C j = (cid:88) j (cid:16)(cid:88) i A ijj B iijk (cid:17) C j is implemented as follows:1) Contract the first tensor and preserve the dummy dimensions: A ijj → A ij .2) Contract the second tensor and preserve the dummy dimensions: B iijk → B ijk .3) Permute and move the summation indices to the end: A ij → A ij , B ijk → B kij .4) Compute the elementwise product on the repeated indices: ( AB ) kij = A ij B kij .5) Sum over the summation indices that do now appear in the other tensors:( AB ) kj = (cid:88) i ( AB ) kij
6) Contract the third tensor and preserve the dummy dimensions: C j → C j .7) Permute and move the summation indices to the end: ( AB ) kj → ( AB ) kj , C j → C j .8) Compute the elementwise product on the repeated indices: ( ABC ) kj = ( AB ) kj C j .9) Sum over the summation indices that do now appear in the other tensors: D k = ( ABC ) k = (cid:88) j ( ABC ) kj
10) Iterate until all the tensors in the summation are considered.The function einstein provides a convenient way to compute general Einstein summationsamong two or more tensors, with or without repeated indices appearing in the same tensor.The function supports both numerical and symbolical calculations implemented via the usageof
C++ templates that operate with generic types and allow the function to work on thedifferent data types without being rewritten for each one. The following example illustratesa sample Einstein summation with mixed data types: D jk = A ij B ki C ii R> A <- array(1:6, dim = c(i = 2, j = 3))R> print(A) [,1] [,2] [,3][1,] 1 3 5[2,] 2 4 6
R> B <- array(1:4, dim = c(k = 2, i = 2))R> print(B) [,1] [,2][1,] 1 3[2,] 2 4 calculus : High Dimensional Numerical and Symbolic Calculus in R R> C <- array(letters[1:4], dim = c(i = 2, i = 2))R> print(C) [,1] [,2][1,] "a" "c"[2,] "b" "d"
R> einstein(A, B, C) [,1] [,2][1,] "1 * (a) + 6 * (d)" "2 * (a) + 8 * (d)"[2,] "3 * (a) + 12 * (d)" "6 * (a) + 16 * (d)"[3,] "5 * (a) + 18 * (d)" "10 * (a) + 24 * (d)"
In the particular case of Einstein summations between two numeric tensors that, after propercontraction and permutation, can be rewritten as C i ...i a ,j ...j b = A i ...i a ,k ...k n B k ...k n ,j ...j b the function implements the following scheme:1) Reshape the tensor A i ...i a ,k ...k n in the matrix A I,K where the dimension of I is theproduct of the dimensions of i . . . i a and the dimensions of K is the product of thedimensions of k . . . k n .2) Reshape the tensor B k ...k n ,j ...j b in the matrix B K,J where the dimension of K is theproduct of the dimensions of k . . . k n and the dimensions of J is the product of thedimensions of j . . . j b .3) Compute the matrix product C IJ = A IK B KJ .4) Reshape the matrix C IJ in the tensor C i ...i a ,j ...j b .In this way, it is sufficient to change the attribute dim of the arrays and the Einstein sum-mation is written in terms of a matrix product that can be computed efficiently in base R .This approach is almost twice as fast as the alternative implementation for the Einstein sum-mation in the R package tensorA for advanced tensor arithmetic with named indices (van denBoogaart 2020). R> a <- array(1:1000000, dim = c(a=2, i=5, j=100, k=50, d=20))R> b <- array(1:100000, dim = c(a=2, j=100, i=5, l=100))R>R> Ta <- tensorA::to.tensor(a)R> Tb <- tensorA::to.tensor(b)R>R> microbenchmark(R+ "calculus" = calculus::einstein(a, b),R+ "tensorA" = tensorA::einstein.tensor(Ta, Tb)R+ ) manuele Guidotti Unit: millisecondsexpr min lq mean median uq max neval cldcalculus 87.165 89.85 95.212 91.716 96.103 207.46 100 atensorA 146.204 152.58 170.172 157.642 169.796 274.36 100 b
The inner product is computed in base R for numeric arrays or via Einstein summation for character arrays: A i ...i n B i ...i n R> 1:3 %inner% letters[1:3] [1] "1 * (a) + 2 * (b) + 3 * (c)"
Dot Product
The dot product between arrays with different dimensions is computed by taking the innerproduct on the last dimensions of the two arrays. It is written in Einstein notation as: A i ...i a j ...j n B j ...j n R> matrix(1:6, byrow = TRUE, nrow = 2, ncol = 3) %dot% letters[1:3] [1] "1 * (a) + 2 * (b) + 3 * (c)" "4 * (a) + 5 * (b) + 6 * (c)"
The outer product is computed in base R for numeric arrays or via Einstein summation for character arrays: A i ...i a B j ...j b R> 1:3 %outer% letters[1:3] [,1] [,2] [,3][1,] "1 * (a)" "1 * (b)" "1 * (c)"[2,] "2 * (a)" "2 * (b)" "2 * (c)"[3,] "3 * (a)" "3 * (b)" "3 * (c)"
The package extends the generalized kronecker product available in base R with support forarrays of type character .4 calculus : High Dimensional Numerical and Symbolic Calculus in R R> 1:3 %kronecker% letters[1:3] [1] "1 * (a)" "1 * (b)" "1 * (c)" "2 * (a)" "2 * (b)" "2 * (c)" "3 * (a)"[8] "3 * (b)" "3 * (c)"
6. Derivatives
The function derivative performs high-order symbolic and numerical differentiation forgeneric tensors with respect to an arbitrary number of variables. The function behaves differ-ently depending on the arguments order , the order of differentiation, and var , the variablenames with respect to which the derivatives are computed.When multiple variables are provided and order is a single integer n , then the n -th orderderivative is computed for each element of the tensor with respect to each variable: D = ∂ ( n ) ⊗ F that is: D i,...,j,k = ∂ ( n ) k F i,...,j where F is the tensor of functions and ∂ ( n ) k denotes the n -th order partial derivative withrespect to the k -th variable.When order matches the length of var , it is assumed that the differentiation order is providedfor each variable. In this case, each element is derived n k times with respect to the k -thvariable, for each of the m variables. D i,...,j = ∂ ( n )1 · · · ∂ ( n m ) m F i,...,j The same applies when order is a named vector giving the differentiation order for eachvariable. For example, order = c(x=1, y=2) differentiates once with respect to x and twicewith respect to y . A call with order = c(x=1, y=0) is equivalent to order = c(x=1) .To compute numerical derivatives or to evaluate symbolic derivatives at a point, the functionaccepts a named vector for the argument var ; e.g. var = c(x=1, y=2) evaluates the deriva-tives in x = 1 and y = 2. For functions where the first argument is used as a parametervector, var should be a numeric vector indicating the point at which the derivatives are tobe calculated. Symbolic derivatives are computed via the D function available in base R . The function isiterated multiple times for second and higher order derivatives. manuele Guidotti f with respect to one or morevariables is approximated up to the degree O ( h p . . . h pm ) by: ∂ n ,...,n m f = ∂ ( n ) x . . . ∂ ( n m ) x m f ( x , . . . , x m ) == n ! . . . n m ! h n . . . h n m m i ( n (cid:88) j = − i ( n · · · i ( nm ) (cid:88) j m = − i ( nm ) C ( n ) j · · · C ( n m ) j m f ( x + j h , . . . , x m + j m h m )where n k is the order of differentiation with respect to the k -th variable, h are the step sizes, i are equal to i ( n ) = (cid:98) ( n + p − / (cid:99) , and the coefficients C ( n ) j are computed by solving thefollowing linear system for each n : C − i C − i +1 C − i +2 ... C − i + n +1 ... C i = ( − i ) . . . ( − . . . i ( − i ) . . . ( − . . . i ( − i ) . . . ( − . . . i ... ... ... ... ...( − i ) n +1 . . . ( − n +1 n +1 . . . i n +1 ... ... ... ... ...( − i ) i . . . ( − i i . . . i i − The summation is computed via Einstein notation by setting: C ( n ) j · · · C ( n m ) j m F j ,...,j m ≡ i ( n (cid:88) j = − i ( n · · · i ( nm ) (cid:88) j m = − i ( nm ) C ( n ) j · · · C ( n m ) j m f ( x + j h , . . . , x m + j m h m ) Symbolic derivatives of univariate functions: ∂ x sin ( x ). R> derivative(f = "sin(x)", var = "x") [1] "cos(x)"
Evaluation of symbolic and numerical derivatives: ∂ x sin ( x ) | x =0 . R> sym <- derivative(f = "sin(x)", var = c(x = 0))R> num <- derivative(f = function(x) sin(x), var = c(x = 0))
Symbolic Numeric1 1
High order symbolic and numerical derivatives: ∂ (4) x sin ( x ) | x =0 .6 calculus : High Dimensional Numerical and Symbolic Calculus in R R> sym <- derivative(f = "sin(x)", var = c(x = 0), order = 4)R> num <- derivative(f = function(x) sin(x), var = c(x = 0), order = 4)
Symbolic Numeric0.000e+00 -2.922e-11
Symbolic derivatives of multivariate functions: ∂ (1) x ∂ (2) y y sin ( x ). R> derivative(f = "y^2*sin(x)", var = c("x", "y"), order = c(1, 2)) [1] "2 * cos(x)"
Numerical derivatives of multivariate functions: ∂ (1) x ∂ (2) y y sin ( x ) | x =0 ,y =0 with degree of accu-racy O ( h ). R> f <- function(x, y) y^2*sin(x)R> derivative(f, var = c(x=0, y=0), order = c(1, 2), accuracy = 6) [1] 2
Symbolic gradient of multivariate functions: ∂ x,y x y . R> derivative("x^2*y^2", var = c("x", "y")) [,1] [,2][1,] "2 * x * y^2" "x^2 * (2 * y)"
High order derivatives of multivariate functions: ∂ (6) x,y x y . R> derivative("x^6*y^6", var = c("x", "y"), order = 6) [,1] [,2][1,] "6 * (5 * (4 * (3 * 2))) * y^6" "x^6 * (6 * (5 * (4 * (3 * 2))))"
Numerical gradient of multivariate functions: ∂ x,y x y | x =1 ,y =2 . R> f <- function(x, y) x^2*y^2R> derivative(f, var = c(x=1, y=2)) [,1] [,2][1,] 8 4
Numerical Jacobian of vector valued functions: ∂ x,y [ xy, x y ] | x =1 ,y =2 . R> f <- function(x, y) c(x*y, x^2*y^2)R> derivative(f, var = c(x=1, y=2)) manuele Guidotti [,1] [,2][1,] 2 1[2,] 8 4 Numerical Jacobian of vector valued functions where the first argument is used as a param-eter vector: ∂ X [ (cid:80) i x i , (cid:81) i x i ] | X =0 . R> f <- function(x) c(sum(x), prod(x))R> derivative(f, var = c(0, 0, 0)) [,1] [,2] [,3][1,] 1 1 1[2,] 0 0 0
The following table compares the accuracy of Richardson extrapolation implemented in thepackage numDeriv (Gilbert and Varadhan 2019) with central finite differences implementedin calculus using accuracy = 4 by default. 10 derivatives have been computed for the fourfunctions: x e x , xsin ( x ), xlog ( x ), e sin ( x ) . The table shows the mean relative error and thecorresponding standard deviation.pkg N Mean SDx exp(x) calculus 1 . × . × − . × − numDeriv 1 . × − . × − . × − x sin(x ) calculus 1 . × . × − . × − numDeriv 1 . × − . × − . × − x log(x ) calculus 1 . × . × − . × − numDeriv 1 . × − . × − . × − exp(sin(x)) calculus 1 . × − . × − . × − numDeriv 1 . × − . × − . × − Although it is known that Richardson extrapolation is usually more accurate than finitedifferences, the results of the two packages are very similar. Both packages produce accuratederivatives with relative errors close to the precision of the machine (10 − ). On the otherhand, calculus proves to be significantly faster than numDeriv for multivariate functions asshown in the following benchmarking. R> x <- rep(0, 1000)R> f <- function(x) sum(x)R> microbenchmark(R+ "calculus O(2)" = calculus::derivative(f, x, accuracy = 2), https://en.wikipedia.org/wiki/Richardson_extrapolation calculus : High Dimensional Numerical and Symbolic Calculus in R R+ "calculus O(4)" = calculus::derivative(f, x, accuracy = 4),R+ "calculus O(6)" = calculus::derivative(f, x, accuracy = 6),R+ "calculus O(8)" = calculus::derivative(f, x, accuracy = 8),R+ "numDeriv" = numDeriv::grad(f, x)R+ )
Unit: millisecondsexpr min lq mean median uq max neval cldcalculus O(2) 17.653 18.546 20.531 20.604 22.260 25.735 100 acalculus O(4) 28.807 30.599 33.137 33.416 34.406 50.392 100 bcalculus O(6) 39.758 43.676 49.144 45.171 50.534 156.399 100 ccalculus O(8) 51.556 55.662 59.491 56.580 59.303 173.515 100 dnumDeriv 102.813 105.561 115.707 108.687 110.493 227.502 100 e
7. Taylor Series
Based on the derivatives in the previous section, the function taylor provides a convenientway to compute the Taylor series of arbitrary unidimensional or multidimensional functions.The mathematical function can be specified both as a character string or as a function .Symbolic or numerical methods are applied accordingly. For univariate functions, the n -thorder Taylor approximation centered in x is given by: f ( x ) (cid:39) n (cid:88) k =0 f ( k ) ( x ) k ! ( x − x ) k where f ( k ) ( x ) denotes the k -th order derivative evaluated in x . By using multi-index nota-tion, the Taylor series is generalized to multidimensional functions with an arbitrary numberof variables: f ( x ) (cid:39) n (cid:88) | k | =0 f ( k ) ( x ) k ! ( x − x ) k where now x = ( x , . . . , x d ) is the vector of variables, k = ( k , . . . , k d ) gives the order ofdifferentiation with respect to each variable f ( k ) = ∂ ( | k | ) f∂ ( k x ··· ∂ ( kd ) xd , and: | k | = k + · · · + k d k ! = k ! · · · k d ! x k = x k · · · x k d d The summation runs for 0 ≤ | k | ≤ n and identifies the set { ( k , · · · , k d ) : k + · · · k d ≤ n } that corresponds to the partitions of the integer n . These partitions can be computed withthe function partitions that is included in the package and optimized in C++ for speed andflexibility. For example, the following call generates the partitions needed for the 2-nd orderTaylor expansion for a function of 3 variables: manuele Guidotti R> partitions(n = 2, length = 3, fill = TRUE, perm = TRUE, equal = FALSE) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10][1,] 0 0 0 1 0 0 2 0 1 1[2,] 0 0 1 0 0 2 0 1 0 1[3,] 0 1 0 0 2 0 0 1 1 0
Based on these partitions, the function taylor computes the corresponding derivatives andbuilds the Taylor series. The output is a list containing the Taylor series, the order of theexpansion, and a data.frame containing the variables, coefficients and degrees of each termin the Taylor series.
R> taylor("exp(x)", var = "x", order = 2) $f[1] "(1) * 1 + (1) * x^1 + (0.5) * x^2"$order[1] 2$termsvar coef degree0 1 1.0 01 x^1 1.0 12 x^2 0.5 2
By default, the series is centered in x = 0 but the function also supports x (cid:54) = 0, themultivariable case, and the approximation of user defined R functions . R> f <- function(x, y) log(y)*sin(x)R> taylor(f, var = c(x = 0, y = 1), order = 2) $f[1] "(0.999999999969436) * x^1*(y-1)^1"$order[1] 2$terms var coef degree0,0 1 0 00,1 (y-1)^1 0 11,0 x^1 0 10,2 (y-1)^2 0 22,0 x^2 0 21,1 x^1*(y-1)^1 1 2 calculus : High Dimensional Numerical and Symbolic Calculus in R
8. Hermite Polynomials
Hermite polynomials are obtained by differentiation of the Gaussian kernel: H ν ( x, Σ) = exp (cid:16) x i Σ ij x j (cid:17) ( − ∂ x ) ν exp (cid:16) − x i Σ ij x j (cid:17) where Σ is a d -dimensional square matrix and ν = ( ν . . . ν d ) is the vector representing theorder of differentiation for each variable x = ( x . . . x d ). In the case where Σ = 1 and x = x the formula reduces to the standard univariate Hermite polynomials: H ν ( x ) = e x ( − ν d ν dx ν e − x High order derivatives of the kernel e − x cannot performed efficiently in base R . The followingexample shows the naive calculation of d dx e − x via the function D : R> D(D(expression(exp(-x^2/2)), "x"), "x") -(exp(-x^2/2) * (2/2) - exp(-x^2/2) * (2 * x/2) * (2 * x/2))
The resulting expression is not simplified and this leads to more and more iterations of thechain rule to compute higher order derivatives. The expression grows fast and soon requireslong computational times and gigabytes of storage.
R> f <- expression(exp(-x^2/2))R> for(i in 1:14) f <- D(f, "x")R> object.size(f)
To overcome this difficulty, the function hermite implements the following scheme. First, itdifferentiates the gaussian kernel. Then, the kernel is dropped from the resulting expression.In this way, the expression becomes a polynomial of degree 1. The taylor series of order1 is computed in order to extract the coefficients of the polynomial and rewrite it compactform. The polynomial is now multiplied by the gaussian kernel and differentiated again. Thekernel is dropped so that the expression becomes a polynomial of degree 2. The taylor series of order 2 is computed and the scheme is iterated until reaching the desired degree ν . The same applies when ν = ( ν . . . ν d ) represents the multi index of multivariate Hermitepolynomials. The scheme allows to reduce the computational time and storage, return a wellformatted output, and generate recursively all the Hermite polynomials of degree ν (cid:48) where | ν (cid:48) | ≤ | ν | . The output is a list of Hermite polynomials of degree ν (cid:48) , where each polynomial isrepresented by the corresponding taylor series. In the univariate case, for ν = 2 the Hermitepolynomials H , H , and H are generated: R> hermite(order = 2) manuele Guidotti $ ` ` $ ` ` $f[1] "(1) * 1"$ ` ` $order[1] 0$ ` ` $termsvar coef degree0 1 1 0$ ` ` $ ` ` $f[1] "(1) * x^1"$ ` ` $order[1] 1$ ` ` $termsvar coef degree0 1 0 01 x^1 1 1$ ` ` $ ` ` $f[1] "(-1) * 1 + (1) * x^2"$ ` ` $order[1] 2$ ` ` $termsvar coef degree0 1 -1 01 x^1 0 12 x^2 1 2 In the multivariate case, where for simplicity Σ ij = δ ij , x = ( x , x ), and | ν | = 2: R> hermite(order = 2, sigma = diag(2), var = c("x1", "x2")) $ ` ` $ ` ` $f[1] "(1) * 1"$ ` ` $order calculus : High Dimensional Numerical and Symbolic Calculus in R [1] 0$ ` ` $termsvar coef degree0,0 1 1 0$ ` ` $ ` ` $f[1] "(1) * x2^1"$ ` ` $order[1] 1$ ` ` $termsvar coef degree0,0 1 0 00,1 x2^1 1 11,0 x1^1 0 1$ ` ` $ ` ` $f[1] "(1) * x1^1"$ ` ` $order[1] 1$ ` ` $termsvar coef degree0,0 1 0 00,1 x2^1 0 11,0 x1^1 1 1$ ` ` $ ` ` $f[1] "(-1) * 1 + (1) * x2^2"$ ` ` $order[1] 2$ ` ` $termsvar coef degree0,0 1 -1 00,1 x2^1 0 11,0 x1^1 0 1 manuele Guidotti ` ` $ ` ` $f[1] "(-1) * 1 + (1) * x1^2"$ ` ` $order[1] 2$ ` ` $termsvar coef degree0,0 1 -1 00,1 x2^1 0 11,0 x1^1 0 10,2 x2^2 0 22,0 x1^2 1 21,1 x1^1*x2^1 0 2$ ` ` $ ` ` $f[1] "(1) * x1^1*x2^1"$ ` ` $order[1] 2$ ` ` $termsvar coef degree0,0 1 0 00,1 x2^1 0 11,0 x1^1 0 10,2 x2^2 0 22,0 x1^2 0 21,1 x1^1*x2^1 1 2
9. Ordinary differential equations
The function ode provides solvers for systems of ordinary differential equations of the type: dydt = f ( t, y ) , y ( t ) = y where y is the vector of state variables. Two solvers are available: the simpler and faster Euler4 calculus : High Dimensional Numerical and Symbolic Calculus in R scheme or the more accurate 4-th order Runge-Kutta method . Although many packagesalready exist to solve ordinary differential equations in R , they usually represent the function f either with an R function - see e.g. deSolve (Soetaert, Petzoldt, and Setzer 2010), odeintr (Keitt 2017), and pracma (Borchers 2019) - or with characters - see e.g. yuima (Brouste,Fukasawa, Hino, Iacus, Kamatani, Koike, Masuda, Nomura, Ogihara, Shimuzu, Uchida, andYoshida 2014). While the representation via R functions is usually more efficient, the sym-bolic representation is easier to adopt for beginners and more flexible for advanced users tohandle systems that might have been generated via symbolic programming. The package cal-culus supports both the representations and uses hashed environments to improve symbolicevaluations. Consider the following system: ddt (cid:20) xy (cid:21) = (cid:20) xx (1 + cos(10 t )) (cid:21) , (cid:20) x y (cid:21) = (cid:20) (cid:21) The vector-valued function f representing the system can be specified as a vector of characters ,or a function returning a numeric vector, giving the values of the derivatives at time t . Theinitial conditions are set with the argument var and the time variable can be specified with timevar . R> sim <- ode(f = c("x", "x*(1+cos(10*t))"),R+ var = c(x = 1, y = 1),R+ times = seq(0, 2*pi, by = 0.001),R+ timevar = "t")
Time V a l ue xy Figure 1: Solution to the system of ordinary differential equations in Section 9 obtained withthe function ode using the Runge-Kutta method. https://en.wikipedia.org/wiki/Euler_method https://en.wikipedia.org/wiki/Runge-Kutta_methods https://cran.r-project.org/web/views/DifferentialEquations.html manuele Guidotti
10. Differential operators
Orthogonal coordinates are a special but extremely common case of curvilinear coordinateswhere the coordinate surfaces all meet at right angles. The chief advantage of non-Cartesiancoordinates is that they can be chosen to match the symmetry of the problem. For example,spherical coordinates are the most common curvilinear coordinate systems and are used inEarth sciences, cartography, quantum mechanics, relativity, and engineering. These coor-dinates may be derived from a set of Cartesian coordinates by using a transformation thatis locally invertible (a one-to-one map) at each point. This means that one can convert apoint given in a Cartesian coordinate system to its curvilinear coordinates and back. Differ-ential operators such as the gradient, divergence, curl, and Laplacian can be transformed fromone coordinate system to another via the usage of scale factors. The package implementsthese operators in Cartesian, polar, spherical, cylindrical, parabolic coordinates, and supportsarbitrary orthogonal coordinates systems defined by custom scale factors.Curvilinear coordinates( q , q , q ) Transformation fromcartesian ( x, y, z ) Scale factorsSpherical polar coordinates( r, θ, φ ) x = r sin θ cos φy = r sin θ sin φz = r cos θ h = 1 h = rh = r sin θ ————— ————— —————Cylindrical polarcoordinates ( r, φ, z ) x = r cos φy = r sin φz = z h = h = 1 h = r ————— ————— —————Parabolic coordinates( u, v, φ ) x = uv cos φy = uv sin φz = 12 ( u − v ) h = h = (cid:112) u + v h = uv ————— ————— —————Parabolic cylindricalcoordinates ( u, v, z ) x = 12 ( u − v ) y = uvz = z h = h = (cid:112) u + v h = 1 The gradient of a scalar-valued function F is the vector ( ∇ F ) i whose components are thepartial derivatives of F with respect to each variable i . In arbitrary orthogonal coordinatesystems, the gradient is expressed in terms of the scale factors h i as follows: https://en.wikipedia.org/wiki/Curvilinear_coordinates https://en.wikipedia.org/wiki/Orthogonal_coordinates calculus : High Dimensional Numerical and Symbolic Calculus in R ( ∇ F ) i = 1 h i ∂ i F The function gradient implements the symbolic and numeric gradient of functions , expressions and characters . In Cartesian coordinates: R> gradient("x*y*z", var = c("x", "y", "z")) [1] "y * z" "x * z" "x * y" and in spherical coordinates:
R> gradient("x*y*z", var = c("x","y","z"), coordinates = "spherical") [1] "1/1 * (y * z)" "1/x * (x * z)" "1/(x*sin(y)) * (x * y)"
To support arbitrary orthogonal coordinate systems, it is possible to pass custom scale factorsto the argument coordinates . For instance, the following call is equivalent to the previousexample in spherical coordinates where the scale factors are now explicitly specified:
R> gradient("x*y*z", var = c("x","y","z"), coordinates = c(1,"x","x*sin(y)")) [1] "1/(1) * (y * z)" "1/(x) * (x * z)" "1/(x*sin(y)) * (x * y)"
Numerical methods are applied when working with functions with the same sintax intro-duced for derivatives in section 6:
R> f <- function(x, y, z) x*y*zR> gradient(f, var = c(x = 1, y = pi/2, z = 0), coordinates = "spherical") [1] 0.0000 0.0000 1.5708 or in vectorized form:
R> f <- function(x) x[1]*x[2]*x[3]R> gradient(f, var = c(1, pi/2, 0), coordinates = "spherical") [1] 0.0000 0.0000 1.5708
When the function F is a tensor-valued function F d ,...,d n , the gradient is computed for eachscalar component. ( ∇ F d ,...,d n ) i = 1 h i ∂ i F d ,...,d n In particular, this reduces to the Jacobian matrix for vector-valued functions F d : manuele Guidotti R> f <- function(x) c(prod(x), sum(x))R> gradient(f, var = c(3, 2, 1)) [,1] [,2] [,3][1,] 2 3 6[2,] 1 1 1 that may be expressed in arbitrary orthogonal coordinate systems.
R> f <- function(x) c(prod(x), sum(x))R> gradient(f, var = c(3, 2, 1), coordinates = "cylindrical") [,1] [,2] [,3][1,] 2 1.00000 6[2,] 1 0.33333 1
Jacobian
The function jacobian is a wrapper for gradient that always returns the Jacobian as a matrix , even in the case of unidimensional scalar-valued functions.
R> f <- function(x) x^2R> jacobian(f, var = c(1)) [,1][1,] 2
Hessian
In Cartesian coordinates, the Hessian of a scalar-valued function F is the square matrix ofsecond-order partial derivatives: ( H ( F )) ij = ∂ ij F It might be tempting to apply the definition of the Hessian as the Jacobian of the gradientto write it in terms of the scale factors. However, this results in a Hessian matrix that is notsymmetric and ignores the distinction between vector and covectors in tensor analysis (seee.g. Masi (2007)). The generalization to arbitrary coordinate system is out of the scope ofthis paper and only Cartesian coordinates are supported:
R> f <- function(x, y, z) x*y*zR> hessian(f, var = c(x = 3, y = 2, z = 1)) [,1] [,2] [,3][1,] 1.2223e-11 1.0000e+00 2.0000e+00[2,] 1.0000e+00 2.7501e-11 3.0000e+00[3,] 2.0000e+00 3.0000e+00 1.1001e-10 calculus : High Dimensional Numerical and Symbolic Calculus in R When the function F is a tensor-valued function F d ,...,d n , the hessian is computed for eachscalar component. ( H ( F d ,...,d n )) ij = ∂ ij F d ,...,d n In this case, the function returns an array of Hessian matrices:
R> f <- function(x, y, z) c(x*y*z, x+y+z)R> h <- hessian(f, var = c(x = 3, y = 2, z = 1)) that can be extracted with the corresponding indices.
R> h[1,,] [,1] [,2] [,3][1,] 1.2223e-11 1.0000e+00 2.0000e+00[2,] 1.0000e+00 2.7501e-11 3.0000e+00[3,] 2.0000e+00 3.0000e+00 1.1001e-10
R> h[2,,] [,1] [,2] [,3][1,] -1.8334e-11 7.8835e-12 -3.6273e-12[2,] 7.8835e-12 -6.4170e-11 1.0907e-11[3,] -3.6273e-12 1.0907e-11 9.1671e-11
The divergence of a vector-valued function F i produces a scalar value ∇ · F representing thevolume density of the outward flux of the vector field from an infinitesimal volume around agiven point. In terms of scale factors, it is expressed as follows: ∇ · F = 1 J (cid:88) i ∂ i (cid:32) Jh i F i (cid:33) where J = (cid:81) i h i . When F is an array of vector-valued functions F d ,...,d n ,i , the divergence is computed for each vector:( ∇ · F ) d ,...,d n = 1 J (cid:88) i ∂ i (cid:32) Jh i F d ,...,d n ,i (cid:33) = 1 J (cid:88) i ∂ i ( J h − i ) F d ,...,d n ,i + J h − i ∂ i ( F d ,...,d n ,i )where the last equality is preferable in practice as the derivatives of the scale factor can becomputed symbolically and the computation of the derivatives of F is more efficient than thedirect computation of ∂ i (cid:0) Jh i F d ,...,d n ,i (cid:1) via finite differences. In Cartesian coordinates: https://en.wikipedia.org/wiki/Divergence manuele Guidotti R> f <- c("x^2", "y^2", "z^2")R> divergence(f, var = c("x","y","z")) [1] "2 * x + 2 * y + 2 * z"
In polar coordinates:
R> f <- c("sqrt(r)/10", "sqrt(r)")R> divergence(f, var = c("r","phi"), coordinates = "polar") [1] "(0.5 * r^-0.5/10 * r + (sqrt(r)/10)) / (1*r)"
And for tensors of vector-valued functions:
R> f <- matrix(c("x^2","y^2","z^2","x","y","z"), nrow = 2, byrow = TRUE)R> divergence(f, var = c("x","y","z")) [1] "2 * x + 2 * y + 2 * z" "1 + 1 + 1"
The same syntax holds for functions where numerical methods are automatically applied:
R> f <- function(x,y,z) matrix(c(x^2,y^2,z^2,x,y,z), nrow = 2, byrow = TRUE)R> divergence(f, var = c(x = 0, y = 0, z = 0)) [1] 0 3
The curl of a vector-valued function F i at a point is represented by a vector whose length anddirection denote the magnitude and axis of the maximum circulation. In 2 dimensions, thecurl is written in terms of the scale factors h and the Levi-Civita symbol (cid:15) as follows: ∇ × F = 1 h h (cid:88) ij (cid:15) ij ∂ i (cid:16) h j F j (cid:17) = 1 h h (cid:32) ∂ (cid:16) h F (cid:17) − ∂ (cid:16) h F (cid:17)(cid:33) In 3 dimensions: ( ∇ × F ) k = h k J (cid:88) ij (cid:15) ijk ∂ i (cid:16) h j F j (cid:17) where J = (cid:81) i h i . This suggests to implement the curl in m + 2 dimensions in such a waythat the formula reduces correctly to the previous cases: https://en.wikipedia.org/wiki/Curl_(mathematics) calculus : High Dimensional Numerical and Symbolic Calculus in R ( ∇ × F ) k ...k m = h k · · · h k m J (cid:88) ij (cid:15) ijk ...k m ∂ i (cid:16) h j F j (cid:17) And in particular, when F is an array of vector-valued functions F d ,...,d n ,i the curl is com-puted for each vector:( ∇ × F ) d ...d n ,k ...k m == h k · · · h k m J (cid:88) ij (cid:15) ijk ...k m ∂ i (cid:16) h j F d ...d n ,j (cid:17) = (cid:88) ij h i h j (cid:15) ijk ...k m ∂ i (cid:16) h j F d ...d n ,j (cid:17) = (cid:88) ij h i h j (cid:15) ijk ...k m (cid:16) ∂ i ( h j ) F d ...d n ,j + h j ∂ i ( F d ...d n ,j ) (cid:17) where the last equality is preferable in practice as the derivatives of the scale factor can becomputed symbolically and the computation of the derivatives of F is more efficient thanthe direct computation of ∂ i (cid:0) h j F d ...d n ,j (cid:1) via finite differences. In 2-dimensional Cartesiancoordinates: R> f <- c("x^3*y^2","x")R> curl(f, var = c("x","y")) [1] "(1) * 1 + (x^3 * (2 * y)) * -1"
In 3 dimensions, for an irrotational vector field:
R> f <- c("x","-y","z")R> curl(f, var = c("x","y","z")) [1] "0" "0" "0"
And for tensors of vector-valued functions:
R> f <- matrix(c("x","-y","z","x^3*y^2","x","0"), nrow = 2, byrow = TRUE)R> curl(f, var = c("x","y","z")) [,1] [,2] [,3][1,] "0" "0" "0"[2,] "0" "0" "(1) * 1 + (x^3 * (2 * y)) * -1"
The same syntax holds for functions where numerical methods are automatically appliedand for arbitrary orthogonal coordinate systems as shown in the previous sections. manuele Guidotti F , resulting in a scalar value giving the flux density of the gradient flowof a function. The Laplacian occurs in differential equations that describe many physicalphenomena, such as electric and gravitational potentials, the diffusion equation for heat andfluid flow, wave propagation, and quantum mechanics. In terms of the scale factor, theoperator is written as: ∇ F = 1 J (cid:88) i ∂ i (cid:32) Jh i ∂ i F (cid:33) where J = (cid:81) i h i . When the function F is a tensor-valued function F d ,...,d n , the laplacian is computed for each scalar component:( ∇ F ) d ...d n = 1 J (cid:88) i ∂ i (cid:32) Jh i ∂ i F d ...d n (cid:33) = 1 J (cid:88) i ∂ i (cid:16) J h − i (cid:17) ∂ i F d ...d n + J h − i ∂ i F d ...d n where the last equality is preferable in practice as the derivatives of the scale factor can becomputed symbolically and the computation of the derivatives of F is more efficient than thedirect computation of ∂ i (cid:0) Jh i ∂ i F (cid:1) via finite differences. In Cartesian coordinates: R> f <- "x^3+y^3+z^3"R> laplacian(f, var = c("x","y","z")) [1] "3 * (2 * x) + 3 * (2 * y) + 3 * (2 * z)"
And for tensors of scalar-valued functions:
R> f <- array(c("x^3+y^3+z^3", "x^2+y^2+z^2", "y^2", "z*x^2"), dim = c(2,2))R> laplacian(f, var = c("x","y","z")) [,1] [,2][1,] "3 * (2 * x) + 3 * (2 * y) + 3 * (2 * z)" "2"[2,] "2 + 2 + 2" "z * 2"
The same syntax holds for functions where numerical methods are automatically appliedand for arbitrary orthogonal coordinate systems as shown in the previous sections.
11. Integrals
The package integrates seamlessly with cubature (Narasimhan et al. C . The function integral provides the interface for multidimensionalintegrals of functions , expressions , and characters in arbitrary orthogonal coordinatesystems. If the package cubature is not installed, the package implements a naive Monte https://en.wikipedia.org/wiki/Laplace_operator calculus : High Dimensional Numerical and Symbolic Calculus in R Carlo integration by default. The function returns a list containing the value of the in-tegral as well as other information on the estimation uncertainty. The integration boundsare specified via the argument bounds : a list containing the lower and upper bound for eachvariable. If the two bounds coincide, or if a single number is specified, the correspondingvariable is not integrated and its value is fixed. For arbitrary orthogonal coordinates q . . . q n the integral is computed as: (cid:90) J · f ( q . . . q n ) dq . . . dq n where J = (cid:81) i h i is the Jacobian determinant of the transformation and is equal to the productof the scale factors h . . . h n . Univariate integral (cid:82) xdx : R> i <- integral(f = "x", bounds = list(x = c(0,1)))R> i$value [1] 0.5 that is equivalent to:
R> i <- integral(f = function(x) x, bounds = list(x = c(0,1)))R> i$value [1] 0.5
Univariate integral (cid:82) yxdx | y =2 : R> i <- integral(f = "y*x", bounds = list(x = c(0,1), y = 2))R> i$value [1] 1
Multivariate integral (cid:82) (cid:82) o yxdxdy : R> i <- integral(f = "y*x", bounds = list(x = c(0,1), y = c(0,1)))R> i$value [1] 0.25
Area of a circle (cid:82) π (cid:82) dA ( r, θ ) R> i <- integral(f = 1,R+ bounds = list(r = c(0,1), theta = c(0,2*pi)),R+ coordinates = "polar")R> i$value manuele Guidotti [1] 3.1416 Volume of a sphere (cid:82) π (cid:82) π (cid:82) dV ( r, θ, φ ) R> i <- integral(f = 1,R+ bounds = list(r = c(0,1), theta = c(0,pi), phi = c(0,2*pi)),R+ coordinates = "spherical")R> i$value [1] 4.1888
As a final example consider the electric potential in spherical coordinates V = πr arisingfrom a unitary point charge: R> V <- "1/(4*pi*r)"
The electric field is determined by the gradient of the potential E = −∇ V : R> E <- -1 %prod% gradient(V, c("r","theta","phi"), coordinates = "spherical")
Then, by Gauss’s law , the total charge enclosed within a given volume is equal to the surfaceintegral of the electric field q = (cid:82) E · dA where · denotes the scalar product between the twovectors. In spherical coordinates, this reduces to the surface integral of the radial componentof the electric field (cid:82) E r dA . The following code computes this surface integral on a spherewith fixed radius r = 1: R> i <- integral(E[1],R+ bounds = list(r = 1, theta = c(0,pi), phi = c(0,2*pi)),R+ coordinates = "spherical")R> i$value [1] 1
As expected q = (cid:82) E · dA = (cid:82) E r dA = 1, the unitary charge generating the electric potential.
12. Summary
This work has presented the calculus package for high dimensional numerical and symbolic cal-culus in R . The library applies numerical methods when working with functions or symbolicprogramming when working with characters or expressions . To describe multidimensionalobjects such as vectors, matrices, and tensors, the package uses the class array regardless ofthe dimension. This is done to prevent unwanted results due to operations among differentclasses such as vector for unidimensional objects or matrix for bidimensional objects. https://en.wikipedia.org/wiki/Electric_potential https://en.wikipedia.org/wiki/Gauss%27s_law calculus : High Dimensional Numerical and Symbolic Calculus in R The package handles multivariate numerical calculus in arbitrary dimensions and coordinatesvia
C++ optimized functions. It achieves approximately the same accuracy for numericaldifferentiation as the numDeriv (Gilbert and Varadhan 2019) package but significantly re-duces the computational time. It supports higher order derivatives and the differentiationof possibly tensor-valued functions. Differential operators such as the gradient, divergence,curl, and Laplacian are made available in arbitrary orthogonal coordinate systems. The Ein-stein summing convention supports expressions involving more than two tensors and tensorswith repeated indices. Besides being more flexible, the summation proves to be faster thanthe alternative implementation found in the tensorA (van den Boogaart 2020) package foradvanced tensor arithmetic with named indices. Unlike mpoly (Kahle 2013) and pracma (Borchers 2019), the package supports multidimensional Hermite polynomials and Taylor se-ries of multivariate functions. The package integrates seamlessly with cubature (Narasimhan et al. C and extends the numerical integration toarbitrary orthogonal coordinate systems.The symbolic counterpart of the numerical methods are implemented whenever possible tomeet the growing needs for R to handle basic symbolic operations. The package provides,among others, symbolic high order derivatives of possibly tensor-valued functions, symbolicdifferential operators such as the gradient, divergence, curl, and Laplacian in arbitrary orthog-onal coordinate systems, symbolic Einstein summing convention and Taylor series expansionof multivariate functions. This is done entirely in R , without depending on external computeralgebra systems.Except for Rcpp (Eddelbuettel and Fran¸cois 2011), the calculus package has no strict depen-dencies in order to provide a stable self-contained toolbox that invites re-use.
13. Computational details
The results in this paper were obtained using R 4.0.0 (R Core Team 2020) with the packages numDeriv tensorA cubature et al. microbenchmark calculus R itself and all packages used are available from the Comprehensive R Archive Network(CRAN) at https://CRAN.R-project.org/ . References
A¨ıt-Sahalia Y (2002). “Maximum likelihood estimation of discretely sampled diffusions: aclosed-form approximation approach.”
Econometrica , (1), 223–262.Andersen MM, Hojsgaard S (2019). “Ryacas: A computer algebra system in R.” Journal ofOpen Source Software , (42). URL https://doi.org/10.21105/joss.01763 .Andersen MM, Hojsgaard S (2020). caracas: Computer Algebra . R package version 1.0.1,URL https://CRAN.R-project.org/package=caracas .Borchers HW (2019). pracma: Practical Numerical Math Functions . R package version 2.2.9,URL https://CRAN.R-project.org/package=pracma . manuele Guidotti Journal of Statistical Software , (4), 1–51. doi:10.18637/jss.v057.i04 . URL .Eberly D (2008). “Derivative approximation by finite differences.” Magic Software, Inc .Eddelbuettel D, Fran¸cois R (2011). “Rcpp: Seamless R and C++ integration.”
Journal ofStatistical Software , (8), 1–18. URL .Frankel T (2011). The geometry of physics: an introduction . Cambridge university press.Gilbert P, Varadhan R (2019). numDeriv: Accurate Numerical Derivatives . R package version2016.8-1.1, URL https://CRAN.R-project.org/package=numDeriv .Grothendieck G, Bellosta CCJG (2019). rSymPy: R Interface to SymPy Computer AlgebraSystem . R package version 0.2-1.2, URL https://CRAN.R-project.org/package=rSymPy .Grothendieck G, Bellosta CJG (2012). rJython: R interface to Python via Jython . R packageversion 0.0-4, URL https://CRAN.R-project.org/package=rJython .Kahle D (2013). “mpoly: Multivariate Polynomials in R.”
R Journal , (1).Keitt TH (2017). odeintr: C++ ODE Solvers Compiled on-Demand . R package version 1.7.1,URL https://CRAN.R-project.org/package=odeintr .Li C, et al. (2013). “Maximum-likelihood estimation for diffusion processes via closed-formdensity expansions.” The Annals of Statistics , (3), 1350–1380.Li J, Bien J, Wells MT (2018). “rTensor: An R Package for Multidimensional Array (Tensor)Unfolding, Multiplication, and Decomposition.” Journal of Statistical Software , (1), 1–31.Masi M (2007). “On compressive radial tidal forces.” American Journal of Physics , (2),116–124.Mersmann O (2019). microbenchmark: Accurate Timing Functions . R package version 1.4-7,URL https://CRAN.R-project.org/package=microbenchmark .Narasimhan B, Johnson SG, Hahn T, Bouvier A, Kiˆeu K (2019). cubature: Adaptive Multivari-ate Integration over Hypercubes . R package version 2.0.4, URL https://CRAN.R-project.org/package=cubature .R Core Team (2020). R: A Language and Environment for Statistical Computing . R Foun-dation for Statistical Computing, Vienna, Austria. URL .Sidiropoulos ND, De Lathauwer L, Fu X, Huang K, Papalexakis EE, Faloutsos C (2017).“Tensor decomposition for signal processing and machine learning.”
IEEE Transactions onSignal Processing , (13), 3551–3582.Soetaert K, Petzoldt T, Setzer RW (2010). “Solving Differential Equations in R: PackagedeSolve.” Journal of Statistical Software , (9), 1–25. ISSN 1548-7660. doi:10.18637/jss.v033.i09 . URL .6 calculus : High Dimensional Numerical and Symbolic Calculus in R Ushey K, Allaire J, Tang Y (2020). reticulate: Interface to ’Python’ . R package version 1.18,URL https://CRAN.R-project.org/package=reticulate .van den Boogaart KG (2020). tensorA: Advanced Tensor Arithmetic with Named Indices . Rpackage version 0.36.2, URL https://CRAN.R-project.org/package=tensorA .Wickham H (2011). “testthat: Get Started with Testing.”
The R Journal , , 5–10. URL https://journal.r-project.org/archive/2011-1/RJournal_2011-1_Wickham.pdf .Yoshida N (1992). “Asymptotic expansion for statistics related to small diffusions.” Journalof the Japan Statistical Society, Japanese Issue , (2), 139–159. Affiliation: