Floating-point arithmetic

Computers hold a finite amount of precision by using (binary) scientific notation.

  • Which can lead to counter-intuitive results.
# NOTE: You should really use all.equal!
(0.1 + 0.2) == 0.3
## [1] FALSE

But what if I need more precision?

Arbitrary-precision arithmetic

  • Useful when dealing with very large numbers (e.g. cryptography)
  • But inherently slower

Note: Whenever possible, use a more accurate fixed-precision algorithm

Eigenvalues of Wishart matrices

\(W \sim W_p(n, \Sigma) \Rightarrow W = X^TX,\) \(X(i,) \sim N_p(0, \Sigma), i=1,\ldots,n.\)

  • Chiani (JMA, 2014) gave an exact algorithm for computing the CDF of the largest eigenvalue of \(W\).

Chiani’s algorithm

  1. Using iterative rules, construct a skew-symmetric matrix \(A\) of dimension \(\approx \min(n,p)\).
  2. Compute \(\sqrt{\det(A)}\).
  3. Divide by normalization constant.

Problem: For \(n,p\) large, both the determinant and the normalization constant are large.


Existing R packages

Some of the available options:

  • gmp and Rmpfr: Link to C/C++ libraries
  • Ryacas: Interfaces with a Computer Algebra System

Boost via BH

“Boost provides free peer-reviewed portable C++ source libraries.”

  • BH allows R users to work with these libraries.
#include <boost/multiprecision/cpp_dec_float.hpp>

using namespace boost::multiprecision;

// Define new type with 100-digit precision
typedef number<cpp_dec_float<100>> mp_float;


“Eigen is a C++ template library for linear algebra”

  • Templates mean that you can do linear algebra with user-defined types!
#include <Eigen/Dense>

Eigen::Matrix<mp_float,100,1> vec;

// Compute dot product
mp_float dot =;


rootWishart: n=25, p=15, fixed precision

x <- seq(0, 100, length.out = 50)
y <- rootWishart::singleWishart(x, p=15, n=25,
plot(x, y, type = 'l', ylim = c(0,1))
Non-monotonic curve

Non-monotonic curve

rootWishart: n=25, p=15, arbitrary precision

x <- seq(0, 100, length.out = 50)
y <- rootWishart::singleWishart(x, p=15, n=25,
plot(x, y, type = 'l', ylim = c(0,1))
Monotonically increasing, as expected

Monotonically increasing, as expected


  • The combination of BH and RcppEigen provides a lot of flexibility.
  • The cost is slower algorithms and the need to code in C++
  • Chiani (2016) proposed a similar algorithm for double Wishart problems, also implemented in rootWishart
  • General rule: More precise algorithms > arbitrary precision