Task view: Numerical Mathematics

Last updated on 2019-06-02 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 statistics.

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 approximation.

  • 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 derivative. 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 algebra
  • Package rmumps provides a wrapper for the MUMPS library, solving large linear systems of equations applying a parallel sparse direct solver
  • 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 convergence.
  • The packages geigen and QZ compute generalized eigenvalues and -vectors for pairs of matrices, and QZ (generalized Schur) decompositions.
  • 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.
  • matrixcalc contains a collection of functions for matrix calculations, special matrices, and tests for matrix properties, e.g., (semi-)positive definiteness.
  • Package onion contains routines for manipulating quaternions and octonians (normed division algebras over the real numbers); quaternions can be useful for handling rotations in three-dimensional space.
  • 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.

Special Functions

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 contributed packages.

  • 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 many more.
  • Airy and Bessel functions, for real and complex numbers, are also computed in package Bessel, with approximations for large arguments.
  • 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.
  • 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.
  • 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 display features.
  • 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 variables.
  • 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.

Differentiation and Integration

D() and deriv() in base R compute derivatives of simple expressions symbolically. Function 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 extent.

  • 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 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. integral2() 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 function 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 quadrature formulas.
  • 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.
  • 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 interpolation, and spline() for cubic (Hermite) spline interpolation, while smooth.spline() performs cubic spline approximation. Base package splines creates periodic interpolation splines in function periodicSpline().

  • Interpolation of irregularly spaced data is possible with the akima package: aspline() for univariate data, bicubic() or 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 cubic interpolation, pchip() for piecewise cubic Hermite interpolation, and sgolay() for Savitzky-Golay smoothing.
  • Package pracma provides barycentric Lagrange interpolation (in 1 and 2 dimensions) in barylag() resp. barylag2d(), 1-dim. akima in akimaInterp(), and interpolation and approximation of data with rational functions, i.e. in the presence of singularities, in ratinterp() and rationalfit().
  • 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 akima::interp and tripack::tri.mesh functions.)
  • Package chebpol contains methods for creating multivariate Chebyshev and multilinear interpolation on regular grids, e.g. the Floater-Hormann barycenter method, or polyharmonic splines for scattered data.
  • 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 spaced data.
  • 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 two packages:

  • Package rootSolve includes function multiroot() 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 quadratic zeros).
  • For solving nonlinear systems of equations the BB package provides Barzilai-Borwein spectral methods in sane(), including a derivative-free variant in dfsane(), and 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 equations with BB or nleqslv.
  • Package FixedPoint provides algorithms for finding fixed point vectors. These algorithms include Anderson acceleration, epsilon extrapolation methods, and minimal polynomial methods.

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.
  • RcppAlgos provides flexible functions for generating combinations or permutations of a vector with or without constraints. The extension package bigIntegerAlgos 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 (OEIS).

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.
  • Package rSymPy accesses the symbolic algebra system 'SymPy' (written in Python) from R. It supports arbitrary precision computations, linear algebra and calculus, solving equations, discrete mathematics, and much more.

Python Interfaces

Python, through its modules 'NumPy', 'SciPy', 'Matplotlib', 'SymPy', and 'pandas', has elaborate and efficient numerical and graphical tools available.

  • 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.
  • R package rPython permits calls from R to Python, while RPy (with Python module 'rpy2') interfaces R from Python. SnakeCharmR is a fork of 'rPython' with several fixes and improvements.
  • 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.
  • findpython is a package designed to find an acceptable Python binary in the path, incl. minimum version or required modules.
  • '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 very structured integration of R with Python resp. Julia.
  • Note that SageMath is a free open source mathematics system based on Python, allowing to run R functions, but also providing access to Maxima, GAP, FLINT, and many more math programs. SageMath can be downloaded or used through a Web interface at CoCalc.

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.

  • The matlab emulation package contains about 30 simple functions, replicating MATLAB functions, using the respective MATLAB names and being implemented in pure R. (See also pracma for many more mathematical functions designed with MATLAB in mind.)
  • Package R.matlab provides tools to read and write MAT files, which is the MATLAB data format. It also enables a one-directional interface with a MATLAB process, sending and retrieving objects through a TCP/IP 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.

  • The Julia interface of the XRJulia package by John Chambers provides direct analogues to Julia function calls. A 'juliaExamples' package is available on Github.
  • JuliaCall provides seamless integration between R and Julia. Using the high-level interface, the user can call any Julia function just like an R function with automatic type conversion.

The commercial programs SAS and Mathematica do have facilities to call R functions. Here is another Computer Algebra System (CAS) in Pure Mathematics that can be called from R.

  • 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.


ADPF — 0.0.1

Use Least Squares Polynomial Regression and Statistical Testing to Improve Savitzky-Golay

akima — 0.6-2

Interpolation of Irregularly and Regularly Spaced Data

arrangements — 1.1.5

Fast Generators and Iterators for Permutations, Combinations and Partitions

BB — 2014.10-1

Solving and Optimizing Large-Scale Nonlinear Systems

Bessel — 0.6-0

Computations and Approximations for Bessel Functions

bigIntegerAlgos — 0.1.2

R Tool for Factoring Big Integers

Brobdingnag — 1.2-6

Very Large Numbers in R

chebpol — 2.1-1

Multivariate Interpolation

combinat — 0.0-8

combinatorics utilities

conicfit — 1.0.4

Algorithms for Fitting Circles, Ellipses and Conics Based on the Work by Prof. Nikolai Chernov

contfrac — 1.1-12

Continued Fractions

cubature — 2.0.3

Adaptive Multivariate Integration over Hypercubes

Deriv — 3.9.0

Symbolic Differentiation

eigeninv — 2011.8-1

Generates (dense) matrices that have a given set of eigenvalues

elliptic — 1.4-0

Weierstrass and Jacobi Elliptic Functions

expint — 0.1-5

Exponential Integral and Incomplete Gamma Function

expm — 0.999-4

Matrix Exponential, Log, 'etc'

fastGHQuad — 1.0

Fast 'Rcpp' Implementation of Gauss-Hermite Quadrature

feather — 0.3.5

R Bindings to the Feather 'API'

features — 2015.12-1

Feature Extraction for Discretely-Sampled Functional Data

findpython — 1.0.5

Functions to Find an Acceptable Python Binary

FixedPoint — 0.6.1

Algorithms for Finding Fixed Point Vectors of Functions

fourierin — 0.2.4

Computes Numeric Fourier Integrals

freegroup — 1.1-0

The Free Group

gaussquad — 1.0-2

Collection of functions for Gaussian quadrature

geigen — 2.3

Calculate Generalized Eigenvalues, the Generalized Schur Decomposition and the Generalized Singular Value Decomposition of a Matrix Pair with Lapack

gmp — 0.5-13.5

Multiple Precision Arithmetic

gsl — 2.1-6

Wrapper for the Gnu Scientific Library

hypergeo — 1.2-13

The Gauss Hypergeometric Function

interp — 1.0-32

Interpolation Methods

irlba — 2.3.3

Fast Truncated Singular Value Decomposition and Principal Components Analysis for Large Dense and Sparse Matrices

JuliaCall — 0.16.6

Seamless Integration Between R and 'Julia'

ktsolve — 1.2

Configurable Function for Solving Families of Nonlinear Equations

lamW — 1.3.0

Lambert-W Function

logOfGamma — 0.0.1

Natural Logarithms of the Gamma Function for Large Values

m2r — 1.0.0

Macaulay2 in R

magic — 1.5-9

Create and Investigate Magic Squares

MASS — 7.3-51.4

Support Functions and Datasets for Venables and Ripley's MASS

matlab — 1.0.2

MATLAB emulation package

Matrix — 1.2-17

Sparse and Dense Matrix Classes and Methods

matrixcalc — 1.0-3

Collection of functions for matrix calculations

MonoPoly — 0.3-10

Functions to Fit Monotone Polynomials

mpoly — 1.1.0

Symbolic Computation and More with Multivariate Polynomials

multipol — 1.0-7

Multivariate Polynomials

mvp — 1.0-8

Fast Symbolic Multivariate Polynomials

mvQuad — 1.0-6

Methods for Multivariate Quadrature

nleqslv — 3.3.2

Solve Systems of Nonlinear Equations

numbers — 0.7-1

Number-Theoretic Functions

numDeriv — 2016.8-1.1

Accurate Numerical Derivatives

onion — 1.2-7

Octonions and Quaternions

optR — 1.2.5

Optimization Toolbox for Solving Linear Systems

orthopolynom — 1.0-5

Collection of functions for orthogonal and orthonormal polynomials

Pade — 0.1-4

Padé Approximant Coefficients

partitions — 1.9-19

Additive Partitions of Integers

permutations — 1.0-4

The Symmetric Group: Permutations of a Finite Set

polyCub — 0.7.1

Cubature over Polygonal Domains

polynom — 1.4-0

A Collection of Functions to Implement a Class for Univariate Polynomial Manipulations

PolynomF — 2.0-2

Polynomials in R

pracma — 2.2.5

Practical Numerical Math Functions

PRIMME — 3.0-0

Eigenvalues and Singular Values and Vectors from Large Matrices

PythonInR — 0.1-7

Use 'Python' from Within 'R'

QZ — 0.1-7

Generalized Eigenvalues and QZ Decomposition

R.matlab — 3.6.2

Read and Write MAT Files and Call MATLAB from Within R

rARPACK — 0.11-0

Solvers for Large Scale Eigenvalue and SVD Problems

Rcpp — 1.0.2

Seamless R and C++ Integration

RcppAlgos — 2.3.4

High Performance Tools for Combinatorics and Computational Mathematics

RcppArmadillo — 0.9.700.2.0

'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library

RcppEigen —

'Rcpp' Integration for the 'Eigen' Templated Linear Algebra Library

reticulate — 1.13

Interface to 'Python'

Rlinsolve — 0.3.0

Iterative Solvers for (Sparse) Linear System of Equations

Rmpfr — 0.7-2

R MPFR - Multiple Precision Floating-Point Reliable

rmumps — 5.2.1-5

Wrapper for MUMPS Library

RootsExtremaInflections — 1.2.1

Finds Roots, Extrema and Inflection Points of a Curve

rootSolve — 1.7

Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations

rPython — 0.0-6

Package Allowing R to Call Python

Rserve — 1.7-3.1

Binary R server

RSpectra — 0.15-0

Solvers for Large-Scale Eigenvalue and SVD Problems

rSymPy — 0.2-1.2

R Interface to SymPy Computer Algebra System

Ryacas — 1.0.0

R Interface to the 'Yacas' Computer Algebra System

schumaker — 1.1

Schumaker Shape-Preserving Spline

signal — 0.7-6

Signal Processing

SimplicialCubature — 1.2

Integration of Functions Over Simplices

SnakeCharmR —

R and Python Integration

SparseGrid — 0.8.2

Sparse grid integration in R

SparseM — 1.77

Sparse Linear Algebra

SphericalCubature — 1.4

Numerical Integration over Spheres and Balls in n-Dimensions; Multivariate Polar Coordinates

statmod — 1.4.32

Statistical Modeling

stinepack — 1.4

Stineman, a Consistently Well Behaved Method of Interpolation

ssvd — 1.0

Sparse SVD

svd — 0.5

Interfaces to Various State-of-Art SVD and Eigensolvers

tripack — 1.3-8

Triangulation of Irregularly Spaced Data

VeryLargeIntegers — 0.1.6

Store and Operate with Arbitrarily Large Integers

XR — 0.7.2

A Structure for Interfaces from R

XRJulia — 0.9.0

Structured Interface to Julia

XRPython — 0.8

Structured Interface to 'Python'

Zseq — 0.2.0

Integer Sequence Generator

Task view list