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

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

\(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\).

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.

Some of the available options:

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

“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 = vec.dot(vec);

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

Non-monotonic curve

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

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