Last updated on 2020-08-22
by Hans W. Borchers
This task view on numerical mathematics lists R packages and functions
that are useful for solving numerical problems in linear algebra and
analysis. It shows that R is a viable computing environment for
implementing and applying numerical methods, also outside the realm of
The task view will not cover differential equations,
optimization problems and solvers, or packages and functions operating
on times series, because all these topics are treated extensively in
the corresponding task views DifferentialEquations,
Optimization, and TimeSeries.
All these task views together will provide a good selection of what is
available in R for the area of numerical mathematics.
The HighPerformanceComputing task view with its many
links for parallel computing may also be of interest.
The task view has been created to provide an overview of the topic.
If some packages are missing or certain topics in numerical math
should be treated in more detail, please let the maintainer know.
Numerical Linear Algebra
As statistics is based to a large extent on linear algebra, many
numerical linear algebra routines are present in R, and some only
implicitly. Examples of explicitly available functions are vector and
matrix operations, matrix (QR) decompositions, solving linear equations,
eigenvalues/-vectors, singular value decomposition, or least-squares
- The recommended package Matrix provides classes and methods
for dense and sparse matrices and operations on them, for example
Cholesky and Schur decomposition, matrix exponential, or norms and
conditional numbers for sparse matrices.
- Recommended package MASS adds generalized (Penrose)
inverses and null spaces of matrices.
- expm computes the exponential, logarithm, and square root
of square matrices, but also powers of matrices or the Frechet
expm() is to be preferred to the function
with the same name in Matrix.
- SparseM provides classes and methods for sparse matrices
and for solving linear and least-squares problems in sparse linear
- Package rmumps provides a wrapper for the MUMPS library,
solving large linear systems of equations applying a parallel sparse
- Rlinsolve is a collection of iterative solvers for sparse
linear system of equations. Stationary iterative solvers such as
Jacobi or Gauss-Seidel, as well as nonstationary (Krylov subspace)
methods are provided.
- svd provides R bindings to state-of-the-art implementations
of singular value decomposition (SVD) and eigenvalue/eigenvector
computations. Package ssvd will obtain sparse SVDs using an
iterative thresholding method, while irlba will compute
approximate singular values/vectors of large matrices.
- Package PRIMME interfaces PRIMME, a C library for computing
eigenvalues and corresponding eigenvectors of real symmetric or
complex Hermitian matrices. It can find largest, smallest, or interior
eigen-/singular values and will apply preconditioning to accelerate
- The packages geigen and QZ compute generalized
eigenvalues and -vectors for pairs of matrices, and QZ (generalized
- eigeninv generates matrices with a given set of
eigenvalues ('inverse eigenvalue problem').
- Package rARPACK, a wrapper for the ARPACK library, is
typically used to compute only a few eigenvalues/vectors, e.g., a
small number of largest eigenvalues.
- Package RSpectra interfaces the 'Spectra' library for
large-scale eigenvalue decomposition and SVD problems.
- optR uses elementary methods of linear algebra (Gauss, LU,
CGM, Cholesky) to solve linear systems.
- Package mbend for bending non-positive-definite (symmetric)
matrices to positive-definiteness, using weighted and unweighted
- matrixcalc contains a collection of functions for matrix
calculations, special matrices, and tests for matrix properties,
e.g., (semi-)positive definiteness; mainly used for teaching and research purposes
- matlib contains a collection of matrix functions for
teaching and learning matrix linear algebra as used in multivariate
statistical methods; mainly for tutorial purposes in learning matrix
algebra ideas using R.
- Package onion contains routines for manipulating
quaternions and octonians (normed division algebras over the real
numbers); quaternions can be useful for handling rotations in
- clifford provides a suite of routines for arbitrary
dimensional Clifford algebras and discusses special cases such as
Lorentz transforms or quaternion multiplication.
- Packages RcppArmadillo and RcppEigen enable the
integration of the C++ template libraries 'Armadillo' resp. 'Eigen'
for linear algebra applications written in C++ and integrated in R
using Rcpp for performance and ease of use.
Many special mathematical functions are present in R, especially
logarithms and exponentials, trigonometric and hyperbolic functions, or
Bessel and Gamma functions. Many more special functions are available in
- Package gsl provides an interface to the 'GNU Scientific
Library' that contains implementations of many special functions,
for example the Airy and Bessel functions, elliptic and exponential
integrals, the hypergeometric function, Lambert's W function, and
- Airy and Bessel functions, for real and complex numbers, are also
computed in package Bessel, with approximations for large
- Package pracma includes special functions, such as
error functions and inverses, incomplete and complex gamma function,
exponential and logarithmic integrals, Fresnel integrals, the
polygamma and the Dirichlet and Riemann zeta functions.
- The hypergeometric (and generalized hypergeometric) function, is
computed in hypergeo, including transformation formulas
and special values of the parameters.
- HypergeoMat evaluates the hypergeometric functions of a
matrix argument through a C++ implementation of Koev and Edelman's algorithm.
- Elliptic and modular functions are available in package
elliptic, including the Weierstrass P function and
Jacobi's theta functions.
There are tools for visualizing complex functions.
- Carlson evaluates Carlson elliptic and incomplete
elliptic integrals (with compex arguments).
- Package expint wraps C-functions from the GNU Scientific
Library to calculate exponential integrals and the incomplete Gamma
function, including negative values for its first argument.
- fourierin computes Fourier integrals of functions of one
and two variables using the Fast Fourier Transform.
- logOfGamma uses approximations to compute the natural
logarithms of the Gamma function for large values.
- Package lamW implements both real-valued branches of the
Lambert W function (using Rcpp).
Function polyroot() in base R determines all zeros of a polynomial,
based on the Jenkins-Traub algorithm. Linear regression function lm()
can perform polynomial fitting when using
poly() in the model
formula (with option
raw = TRUE).
- Packages PolynomF (recommended) and polynom
provide similar functionality for manipulating univariate polynomials,
like evaluating polynomials (Horner scheme), or finding their roots.
'PolynomF' generates orthogonal polynomials and provides graphical
- polyMatrix (based on 'polynom') implements basic matrix
operations and provides thus an infrastructure for the manipulation
of polynomial matrices.
- Package MonoPoly fits univariate polynomials to given data,
applying different algorithms.
- For multivariate polynomials, package multipol provides
various tools to manipulate and combine these polynomials of several
- Package mpoly facilitates symbolic manipulations on
multivariate polynomials, including basic differential calculus
operations on polynomials, plus some Groebner basis calculations.
- mvp enables fast manipulation of symbolic multivariate
polynomials, using print and coercion methods from the 'mpoly'
package, but offers speed improvements.
- Package orthopolynom consists of a collection of functions
to construct orthogonal polynomials and their recurrence relations,
among them Chebyshev, Hermite, and Legendre polynomials, as well as
spherical and ultraspherical polynomials. There are functions to
operate on these polynomials.
- Symbolic calculation and evaluation of the Jack polynomials, zonal
polynomials (appear in random matrix theory), and Schur polynomials
(appear in combinatorics) is available in package jack.
- The Free Algebra in R package freealg handles multivariate
polynomials with non-commuting indeterminates.
Differentiation and Integration
deriv() in base R compute
derivatives of simple expressions symbolically.
integrate() implements an approach for numerically
integrating univariate functions in R. It applies adaptive Gauss-Kronrod
quadrature and can handle singularities and unbounded domains to a certain
- Package Deriv provides an extended solution for symbolic
differentiation in R; the user can add custom derivative rules, and
the output for a function will be an executable function again.
- numDeriv sets the standard for numerical differentiation
in R, providing numerical gradients, Jacobians, and Hessians, computed
by simple finite differences, Richardson extrapolation, or the highly
accurate complex step approach.
- Package dual achieves automatic differentiation (for
univariate functions) by employing dual numbers; for a mathematical
function its value and its exact first derivative are returned.
- Package Non-Contradiction/autodiffr (on Github)
provides an R wrapper for the Julia packages ForwardDiff.jl and
ReverseDiff.jl to do automatic differentiation for native R functions.
- pracma contains functions for computing numerical
derivatives, including Richardson extrapolation or complex step.
fderiv() computes numerical derivatives of higher orders.
pracma has several routines for numerical integration:
adaptive Lobatto quadrature, Romberg integration, Newton-Cotes
formulas, Clenshaw-Curtis quadrature rules.
integrates functions in two dimensions, also for domains characterized
by polar coordinates or with variable interval limits.
- Package gaussquad contains a collection of functions to
perform Gaussian quadrature, among them Chebyshev, Hermite, Laguerre,
and Legendre quadrature rules, explicitly returning nodes and weights
in each case. Function
gaussquad() in package
statmod does a similar job.
- Package fastGHQuad provides a fast Rcpp-based
implementation of (adaptive) Gauss-Hermite quadrature.
- Adaptive multivariate integration over hyper-rectangles in
n-dimensional space is available in package cubature as
adaptIntegrate(), based on a C library of the
same name. The integrand functions can even be multi-valued.
vegas() includes an approach to Monte Carlo
integration based on importance sampling.
- mvQuad provides methods for generating multivariate grids
that can be used for multivariate integration. These grids will be
based on different quadrature rules such as Newton-Cotes or Gauss
- Package SparseGrid provides another approach to
multivariate integration in high-dimensional spaces. It creates sparse
n-dimensional grids that can be used as with quadrature rules.
- Package SphericalCubature employs cubature to
integrate functions over unit spheres and balls in n-dimensional
space; SimplicialCubature provides methods to integrate
functions over m-dimensional simplices in n-dimensional space.
Both packages comprise exact methods for polynomials.
- Package polyCub holds some routines for numerical
integration over polygonal domains in two dimensions.
- Package Pade calculates the numerator and denominator
coefficients of the Pade approximation, given the Taylor series
coefficients of sufficient length.
- calculus provides efficient functions for high-dimensional
numerical and symbolic calculus, including accurate higher-order
derivatives, Taylor series expansion, differential operators, and Monte-Carlo integration in orthogonal coordinate systems.
- features extracts features from functional data, such as
first and second derivatives, or curvature at critical points, while
RootsExtremaInflections finds roots, extrema and
inflection points of curves defined by discrete points.
Interpolation and Approximation
Base R provides functions
approx() for constant and linear
spline() for cubic (Hermite) spline
smooth.spline() performs cubic spline
approximation. Base package splines creates periodic interpolation
splines in function
- Interpolation of irregularly spaced data is possible with the
aspline() for univariate data,
interp() for data on a 2D
rectangular domain. (This package is distributed under ACM license and
not available for commercial use.)
- Package signal contains several filters to smooth
discrete data, notably
interp1() for linear, spline, and
pchip() for piecewise cubic Hermite
sgolay() for Savitzky-Golay
- Package pracma provides barycentric Lagrange interpolation
(in 1 and 2 dimensions) in
barylag2d(), 1-dim. akima in
and interpolation and approximation of data with rational functions,
i.e. in the presence of singularities, in
- The interp package provides bivariate data interpolation
on regular and irregular grids, either linear or using splines.
Currently the piecewise linear interpolation part is implemented.
(It is intended to provide a free replacement for the ACM licensed
- Package chebpol contains methods for creating multivariate
Chebyshev and other multilinear interpolations on regular grids, e.g.
the Floater-Hormann barycenter method, or polyharmonic splines for
- Package kubik provides (constructs, plots and evaluates)
constrained cubic Hermite splines and computes their derivatives,
indefinite integrals, as well as roots, optima, or inflection points.
- tripack for triangulation of irregularly spaced data is a
constrained two-dimensional Delaunay triangulation package providing
both triangulation and generation of Voronoi mosaics of irregular
sinterp() in package stinepack realizes
interpolation based on piecewise rational functions by applying
Stineman's algorithm. The interpolating function will be monotone in
regions where the specified points change monotonically.
Schumaker() in package schumaker implements
shape-preserving splines, guaranteed to be monotonic resp. concave
or convex if the data is monotonic, concave, or convex.
- ADPF uses least-squares polynomial regression and
statistical testing to improve Savitzky-Golay smoothing.
- Package conicfit provides several (geometric and algebraic)
algorithms for fitting circles, ellipses, and conics in general.
Root Finding and Fixed Points
uniroot(), implementing the Brent-Decker algorithm, is the
basic routine in R to find roots of univariate functions. There are
implementations of the bisection algorithm in several contributed
packages. For root finding with higher precision there is function
unirootR() in the multi-precision package Rmpfr.
And for finding roots of multivariate functions see the following
- Package rootSolve includes function
for finding roots of systems of nonlinear (and linear) equations; it
also contains an extension
uniroot.all() that attempts to
find all zeros of a univariate function in an intervall (excepting
- For solving nonlinear systems of equations the BB package
provides Barzilai-Borwein spectral methods in
including a derivative-free variant in
multi-start features with sensitivity analysis.
- Package nleqslv solves nonlinear systems of equations
using alternatively the Broyden or Newton method, supported by
strategies such as line searches or trust regions.
- ktsolve defines a common interface for solving a set of
- Package FixedPoint provides algorithms for finding fixed
point vectors. These algorithms include Anderson acceleration,
epsilon extrapolation methods, and minimal polynomial methods.
- Package daarem implements the DAAREM method for accelerating
the convergence of any smooth, monotone, slow fixed point iteration.
- Algorithms for accelerating the convergence of slow, monotone
sequences from smooth contraction mappings such as the
expectation-maximization (EM) algorithm are provided in packages
SQUAREM resp. turboEM.
Discrete Mathematics and Number Theory
Not so many functions are available for computational number theory.
Note that integers in double precision can be represented exactly up to
2^53 - 1, above that limit a multi-precision package such as
gmp is needed, see below.
- Package numbers provides functions for factorization, prime
numbers, twin primes, primitive roots, modular inverses, extended GCD,
etc. Included are some number-theoretic functions like divisor
functions or Euler's Phi function.
- contfrac contains various utilities for evaluating
continued fractions and partial convergents.
- magic creates and investigates magical squares and
hypercubes, including functions for the manipulation and analysis
of arbitrarily dimensioned arrays.
- Package freegroup provides functionality for manipulating
elements of a free group including juxtaposition, inversion,
multiplication by a scalar, power operations, and Tietze forms.
- The partitions package enumerates additive partitions
of integers, including restricted and unequal partitions.
- permutations treats permutations as invertible functions of
finite sets and includes several mathematical operations on them.
- Package combinat generates all permutations or all
combinations of a certain length of a set of elements (i.e. a vector);
it also computes binomial coefficients.
- Package arrangements provides generators and iterators for
permutations, combinations and partitions. The iterators allow users
to generate arrangements in a fast and memory efficient manner.
Permutations and combinations can be drawn with/without replacement
and support multisets.
- Package set6 implements (as R6 classes) many forms of
mathematical sets (sets, tuples, intervals) and allows for standard
operations on them (unions, products, differences).
- RcppAlgos provides flexible functions for generating
combinations or permutations of a vector with or without constraints;
the extension package RcppBigIntAlgos features a
quadratic sieve algorithm for completely factoring large integers.
- Package Zseq generates well-known integer sequences; the
'gmp' package is adopted for computing with arbitrarily large numbers.
Every function has on its help page a hyperlink to the corresponding
entry in the On-Line Encyclopedia of Integer Sequences
- Package primes provides quite fast (Rcpp) functions for
identifying and generating prime numbers. And primefactr
uses prime factorization for computations such as reducing ratios
of large factorials.
Multi-Precision Arithmetic and Symbolic Mathematics
- Multiple precision arithmetic is available in R through package
gmp that interfaces to the GMP C library. Examples are
factorization of integers, a probabilistic prime number test, or
operations on big rationals -- for which linear systems of equations
can be solved.
- Multiple precision floating point operations and functions are
provided through package Rmpfr using the MPFR and GMP
libraries. Special numbers and some special functions are included,
as well as routines for root finding, integration, and optimization
in arbitrary precision.
- Brobdingnag handles very large numbers by holding their
logarithm plus a flag indicating their sign. (An excellent vignette
explains how this is done using S4 methods.)
- VeryLargeIntegers implements a multi-precision library that
allows to store and manage arbitrarily big integers; it includes
probabilistic primality tests and factorization algorithms.
- Package Ryacas interfaces the computer algebra system
'Yacas'; it supports symbolic and arbitrary precision computations
in calculus and linear algebra.
- Packages caracas (based on 'reticulate') and
rSymPy (based on 'rjava') both access the symbolic algebra
system 'SymPy'; supported are symbolic operations in linear algebra and
calculus, such as eigenvalues, derivatives, integrals, limits, etc.,
computing special functions, or solving systems of equations.
- Package symengine provides an interface to 'SymEngine',
a C++ library for fast symbolic calculations, such as manipulating
mathematical expressions, finding exact derivatives, performing
symbolic matrix computations, or solving ordinary differential equations (numerically).
Python, through its modules 'NumPy', 'SciPy', 'Matplotlib', 'SymPy',
and 'pandas', has elaborate and efficient numerical and graphical tools
- reticulate is an R interface to Python modules, classes,
and functions. When calling Python in R data types are automatically
converted to their equivalent Python types; when values are returned
from Python to R they are converted back to R types. This package from
the RStudio team is a kind of standard for calling Python from R.
- PythonInR is another package to interact with Python from
within R. It provides Python classes for vectors, matrices and
data.frames which allow an easy conversion from R to Python and back.
- feather provides bindings to read and write feather files,
a lightweight binary data store designed for maximum speed.
This storage format can also be accessed in Python, Julia, or Scala.
- 'pyRserve' is a Python module for connecting Python to an R process
running Rserve as an RPC gateway. This R process can run on
a remote machine, variable access and function calls will be delegated
through the network.
- XRPython (and 'XRJulia') are based on John Chambers'
XR package and his "Extending R" book and allow for a
structured integration of R with Python resp. Julia.
MATLAB, Octave, Julia, and other Interfaces
Interfaces to numerical computation software such as MATLAB
(commercial) or Octave (free) will be important when solving difficult
numerical problems. Unfortunately, at the moment there is no package
allowing to call Octave functions from within R.
- The matlab emulation package contains about 30 simple
functions, replicating MATLAB functions, using the respective MATLAB
names and being implemented in pure R.
- Packages rmatio and R.matlab provides tools to
read and write MAT files (the MATLAB data format) for versions 4 and 5.
'R.matlab' also enables a one-directional interface with a MATLAB v6
process, sending and retrieving objects through a TCP connection.
Julia is "a high-level, high-performance dynamic programming language
for numerical computing", which makes it interesting for optimization
problems and other demanding scientific computations in R.
- JuliaCall provides seamless integration between R and Julia;
the user can call Julia functions just like any R function, and R
functions can be called in the Julia environment, both with reasonable
automatic type conversion.
- JuliaConnectoR provides a functionally oriented interface
for integrating Julia with R; imported Julia functions can be called
as R functions; data structures are converted automatically.
- Package XRJulia provides an interface from R to computations
in the Julia language, based on the interface structure described in
the book "Extending R" by John M. Chambers.
Java Math functions can be employed through the 'rjava' or 'rscala'
interfaces. Then package commonsMath allows calling Java JAR
files of the Apache Commons Mathematics Library, a specialized library for
all aspects of numerics, optimization, and differential equations.
SageMath is an open source
mathematics system based on Python, allowing to run R functions, but
also providing access to systems like Maxima, GAP, FLINT, and many more
math programs. SageMath can be freely used through a Web interface at
Package m2r provides a persistent interface to Macauley2,
an extended software program supporting research in algebraic geometry
and commutative algebra. Macauley2 has to be installed independently,
otherwise a Macauley2 process in the cloud will be instantiated.
Please note that commercial programs such as MATLAB, Maple, or Mathematica have facilities to call R functions.