In RobinHankin/Brobdingnag: Very Large Numbers in R

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 3.5, message = FALSE, warning = FALSE)
options(width = 80, tibble.width = Inf)

knitr::include_graphics(system.file("help/figures/Brobdingnag.png", package = "Brobdingnag"))

To cite the Brobdingnag package in publications, please use @hankin2007. R package Brobdingnag has basic functionality for matrices. It includes matrix multiplication and addition, but determinants and matrix inverses are not implemented. First load the package:

library("Brobdingnag")

The standard way to create a Brobdgingnagian matrix (a brobmat) is to use function brobmat() which takes arguments similar to matrix() and returns a matrix of entries created with brob():

M1 <- brobmat(-10:13,4,6)
colnames(M1) <- state.abb[1:6]
M1

Above, note that all entries of M1 are greater than zero; M1[1,1], for example, is exp(-10), or $e^{-10}\simeq 4.54\times 10^{-5}$. For negative matrix entries, function brobmat() takes a Boolean argument positive that specifies the sign:

M2 <- brobmat(
c(1,104,-66,45,1e40,-2e40,1e-200,232.2),2,4,
positive=c(T,F,T,T,T,F,T,T))
M2

Standard matrix arithmetic is implemented, thus:

rownames(M2) <- c("a","b")
colnames(M2) <- month.abb[1:4]
M2
M2[2,3] <- 0
M2
M2+1000

We can also do matrix multiplication, although it is slow:

M2 %*% M1

Numerical verification: matrix multiplication

We will verify matrix multiplication by carrying out the same operation in two different ways. First, create two largish Brobdingnagian matrices:

nrows <- 11
ncols <- 18
M3 <- brobmat(rnorm(nrows*ncols),nrows,ncols,positive=sample(c(T,F),nrows*ncols,replace=T))
M4 <- brobmat(rnorm(nrows*ncols),ncols,nrows,positive=sample(c(T,F),nrows*ncols,replace=T))
M3[1:3,1:3]

Now calculate the matrix product by coercing to numeric matrices and using base R matrix multiplication:

p1 <- as.matrix(M3) %*% as.matrix(M4)

Secondly, we use Brobdingnagian matrix multiplication, and then coercing to numeric:

p2 <- as.matrix(M3 %*% M4)

The difference:

max(abs(p1-p2))

is small. Now the other way:

q1 <- M3 %*% M4
q2 <- as.brobmat(as.matrix(M3) %*% as.matrix(M4))
max(abs(as.brob(q1-q2)))

Above we see that the difference of exp(-30.267) (about $7\times 10^{-14}$), is small.

Numerical verification: integration with the `cubature` package

The matrix functionality of the Brobdingnag package was originally written to leverage the functionality of the cubature package. Here I give some numerical verification for this.

Suppose we wish to evaluate

[ \int_{x=0}^{x=4}(x^2-4)\,dx ]

using numerical methods. See how the integrand includes positive and negative values; the theoretical value is $\frac{16}{3}=5.33\ldots$. The cubature idiom for this would be

library("cubature")

f.numeric <- function(x){x^2 - 4}

out.num <- cubature::hcubature(f = f.numeric, lowerLimit = 0, upperLimit = 4, vectorInterface = TRUE)
out.num

and the Brobdingnagian equivalent would be

f.brob <- function(x) {
    x <- as.brob(x[1, ])
    as.matrix( brobmat(x^2 - 4, ncol = length(x)))
}

out.brob <- cubature::hcubature(f = f.brob, lowerLimit = 0, upperLimit = 4, vectorInterface = TRUE)
out.brob