In RobinHankin/freealg: The Free Algebra

knitr::opts_chunk$set(echo = TRUE)
set.seed(0)

Numerical verification of the freealg package

![](`r system.file("help/figures/freealg.png", package = "freealg")`){width=10%}

Here we verify that the freealg package is consistent.

First, we observe that $n\times n$ matrices are an algebra over the reals, and therefore should obey all the identities of a free algebra. We will evaluate an expression in two different ways (one using the package and the other using base R matrix algebra) and verify that the results agree. We will evaluate an expression on three indeterminates, A,B,C, here random $5\times 5$ matrices:

library("freealg")
 A <- matrix(rnorm(25),5,5)
 B <- matrix(rnorm(25),5,5)
 C <- matrix(rnorm(25),5,5)

Now we will take a random-ish element of the free algebra and call it x:

x <- structure(list(indices = list(1L, 1:2, c(1L, 2L, 1L, 1L), 2L, 
    3L), coeffs = c(9, 6, 1, 9, 3)), class = "freealg")
x

We will evaluate x^2 in two different ways. Firstly, direct matrix multiplication:

xm <- 9*A + 6*A%*%B + A%*%B%*%A%*%A + 9*B + 3*C  # object "x" adapted for matrix multiplication
(xcubed_matrix <- xm %*% xm %*% xm)

And secondly using freealg:

x^2

We can turn the above into an expression (rather clumsily) using sed. Cut-and-paste the above output into a file called e.txt [if using the text above, strip out the # symbols and the first line saying "free algebra..."] and then the following command

cat e.txt|tr 'a-z' 'A-Z' | sed -e ':loop' -e 's/\([A-Z]\)\([A-Z]\)/\1%*%\2/g' -e 't loop'

gives this:

xcubed_freealg <- (
+ 729*A%*%A%*%A + 486*A%*%A%*%A%*%B + 81*A%*%A%*%A%*%B%*%A%*%A + 729*A%*%A%*%B + 486*A%*%A%*%B%*%A + 81*A%*%A%*%B%*%A%*%A%*%A + 54*A%*%A%*%B%*%A%*%A%*%A%*%B
+ 9*A%*%A%*%B%*%A%*%A%*%A%*%B%*%A%*%A + 81*A%*%A%*%B%*%A%*%A%*%B + 27*A%*%A%*%B%*%A%*%A%*%C + 324*A%*%A%*%B%*%A%*%B + 54*A%*%A%*%B%*%A%*%B%*%A%*%A + 486*A%*%A%*%B%*%B +
162*A%*%A%*%B%*%C + 243*A%*%A%*%C + 729*A%*%B%*%A + 486*A%*%B%*%A%*%A + 81*A%*%B%*%A%*%A%*%A%*%A + 54*A%*%B%*%A%*%A%*%A%*%A%*%B + 9*A%*%B%*%A%*%A%*%A%*%A%*%B%*%A%*%A
+ 81*A%*%B%*%A%*%A%*%A%*%B + 54*A%*%B%*%A%*%A%*%A%*%B%*%A + 9*A%*%B%*%A%*%A%*%A%*%B%*%A%*%A%*%A + 6*A%*%B%*%A%*%A%*%A%*%B%*%A%*%A%*%A%*%B + 1*A%*%B%*%A%*%A%*%A%*%B%*%A%*%A%*%A%*%B%*%A%*%A +
9*A%*%B%*%A%*%A%*%A%*%B%*%A%*%A%*%B + 3*A%*%B%*%A%*%A%*%A%*%B%*%A%*%A%*%C + 36*A%*%B%*%A%*%A%*%A%*%B%*%A%*%B + 6*A%*%B%*%A%*%A%*%A%*%B%*%A%*%B%*%A%*%A + 54*A%*%B%*%A%*%A%*%A%*%B%*%B +
18*A%*%B%*%A%*%A%*%A%*%B%*%C + 27*A%*%B%*%A%*%A%*%A%*%C + 324*A%*%B%*%A%*%A%*%B + 81*A%*%B%*%A%*%A%*%B%*%A + 54*A%*%B%*%A%*%A%*%B%*%A%*%A + 54*A%*%B%*%A%*%A%*%B%*%A%*%B +
9*A%*%B%*%A%*%A%*%B%*%A%*%B%*%A%*%A + 81*A%*%B%*%A%*%A%*%B%*%B + 27*A%*%B%*%A%*%A%*%B%*%C + 27*A%*%B%*%A%*%A%*%C%*%A + 18*A%*%B%*%A%*%A%*%C%*%A%*%B + 3*A%*%B%*%A%*%A%*%C%*%A%*%B%*%A%*%A +
27*A%*%B%*%A%*%A%*%C%*%B + 9*A%*%B%*%A%*%A%*%C%*%C + 972*A%*%B%*%A%*%B + 324*A%*%B%*%A%*%B%*%A + 81*A%*%B%*%A%*%B%*%A%*%A + 54*A%*%B%*%A%*%B%*%A%*%A%*%A +
36*A%*%B%*%A%*%B%*%A%*%A%*%A%*%B + 6*A%*%B%*%A%*%B%*%A%*%A%*%A%*%B%*%A%*%A + 54*A%*%B%*%A%*%B%*%A%*%A%*%B + 18*A%*%B%*%A%*%B%*%A%*%A%*%C + 216*A%*%B%*%A%*%B%*%A%*%B + 36*A%*%B%*%A%*%B%*%A%*%B%*%A%*%A
+ 324*A%*%B%*%A%*%B%*%B + 108*A%*%B%*%A%*%B%*%C + 162*A%*%B%*%A%*%C + 729*A%*%B%*%B + 486*A%*%B%*%B%*%A + 324*A%*%B%*%B%*%A%*%B +
54*A%*%B%*%B%*%A%*%B%*%A%*%A + 486*A%*%B%*%B%*%B + 162*A%*%B%*%B%*%C + 243*A%*%B%*%C + 162*A%*%B%*%C%*%A + 108*A%*%B%*%C%*%A%*%B + 18*A%*%B%*%C%*%A%*%B%*%A%*%A
+ 162*A%*%B%*%C%*%B + 54*A%*%B%*%C%*%C + 243*A%*%C%*%A + 162*A%*%C%*%A%*%B + 27*A%*%C%*%A%*%B%*%A%*%A + 243*A%*%C%*%B + 81*A%*%C%*%C +
729*B%*%A%*%A + 486*B%*%A%*%A%*%B + 81*B%*%A%*%A%*%B%*%A%*%A + 729*B%*%A%*%B + 486*B%*%A%*%B%*%A + 81*B%*%A%*%B%*%A%*%A%*%A + 54*B%*%A%*%B%*%A%*%A%*%A%*%B +
9*B%*%A%*%B%*%A%*%A%*%A%*%B%*%A%*%A + 81*B%*%A%*%B%*%A%*%A%*%B + 27*B%*%A%*%B%*%A%*%A%*%C + 324*B%*%A%*%B%*%A%*%B + 54*B%*%A%*%B%*%A%*%B%*%A%*%A + 486*B%*%A%*%B%*%B +
162*B%*%A%*%B%*%C + 243*B%*%A%*%C + 729*B%*%B%*%A + 486*B%*%B%*%A%*%B + 81*B%*%B%*%A%*%B%*%A%*%A + 729*B%*%B%*%B + 243*B%*%B%*%C +
243*B%*%C%*%A + 162*B%*%C%*%A%*%B + 27*B%*%C%*%A%*%B%*%A%*%A + 243*B%*%C%*%B + 81*B%*%C%*%C + 243*C%*%A%*%A + 162*C%*%A%*%A%*%B +
27*C%*%A%*%A%*%B%*%A%*%A + 243*C%*%A%*%B + 162*C%*%A%*%B%*%A + 27*C%*%A%*%B%*%A%*%A%*%A + 18*C%*%A%*%B%*%A%*%A%*%A%*%B + 3*C%*%A%*%B%*%A%*%A%*%A%*%B%*%A%*%A +
27*C%*%A%*%B%*%A%*%A%*%B + 9*C%*%A%*%B%*%A%*%A%*%C + 108*C%*%A%*%B%*%A%*%B + 18*C%*%A%*%B%*%A%*%B%*%A%*%A + 162*C%*%A%*%B%*%B + 54*C%*%A%*%B%*%C + 81*C%*%A%*%C +
243*C%*%B%*%A + 162*C%*%B%*%A%*%B + 27*C%*%B%*%A%*%B%*%A%*%A + 243*C%*%B%*%B + 81*C%*%B%*%C + 81*C%*%C%*%A + 54*C%*%C%*%A%*%B + 9*C%*%C%*%A%*%B%*%A%*%A
+ 81*C%*%C%*%B + 27*C%*%C%*%C
)

and then

xcubed_freealg - xcubed_matrix

is zero, to within numerical tolerances.

Calculus

We will verify that the deriv() function behaves correctly.

o <- as.freealg("4*aaabacaa + 7*ccab -9*babaa")
do <- deriv(o,1)  # differentiate WRT a
o
do

Now evaluate o in matrix-speak:

f <- function(A,B,C){4*A%*%A%*%A%*%B%*%A%*%C%*%A%*%A + 7*C%*%C%*%A%*%B -9*B%*%A%*%B%*%A%*%A}
small <- 1e-5
dA <- small*matrix(rnorm(25),5,5)
fdash <- function(A,B,C,dA){(     
   + 4*A%*%A%*%A%*%B%*%A%*%C%*%A%*%(dA) + 4*A%*%A%*%A%*%B%*%A%*%C%*%(dA)%*%A
   + 4*A%*%A%*%A%*%B%*%(dA)%*%C%*%A%*%A + 4*A%*%A%*%(dA)%*%B%*%A%*%C%*%A%*%A
   + 4*A%*%(dA)%*%A%*%B%*%A%*%C%*%A%*%A
   - 9*B%*%A%*%B%*%A%*%(dA) - 9*B%*%A%*%B%*%(dA)%*%A - 9*B%*%(dA)%*%B%*%A%*%A
   + 7*C%*%C%*%(dA)%*%B + 4*(dA)%*%A%*%A%*%B%*%A%*%C%*%A%*%A
)}
approximate <- fdash(A,B,C,dA)
exact <- f(A+dA,B,C) - f(A,B,C)
approximate
exact
approximate-exact

Above, we see the difference between the approximate and the exact is small, as might be expected. We can improve the accuracy by evaluating the approximation at an intermediate point:

approximate <- fdash(A+dA/2,B,C,dA)
exact <- f(A+dA,B,C) - f(A,B,C)
approximate-exact

and we see an improved accuracy.

RobinHankin/freealg documentation built on Dec. 24, 2024, 3:16 a.m.