## 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;

## RcppEigen

“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);

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

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))

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

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))

## Conclusion

• 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