This is another test of Joy Christian’s one-page paper: http://arxiv.org/pdf/1103.1879v1.pdf

This time I implement Albert Jan Wonnink’s “fix” to the central derivation in the analysis, namely by switching the “handedness” convention for cross product depending on the sign of the hidden variable.

I use the R package “onion” to implement the quaternions, through which I am able to define the various bivector algebras used by Joy Christian in his paper.

First of all we load the library “onion” and we define a function for sampling points on S2, uniformly at random. We are going to test Christian’s model by picking two measurement directions a, b, completely at random

library(onion)


S2 <- function(N) {    
    z <- runif(N, -1, 1)
    t <- runif(N, 0, 2*pi)
    r <- sqrt(1 - z^2)
    x <- r * cos(t)
    y <- r * sin(t)
    rbind(x, y, z)
}

set.seed(1234)
a <- S2(1)
b <- S2(1)

a
##         [,1]
## x -0.4564959
## y -0.4412610
## z -0.7725932
b
##         [,1]
## x -0.6970031
## y -0.6829515
## z  0.2185495

Next we will create the bivector basis elements and test that they have the desired properties; Christian’s formula (3).

Note that in “onion”, H1, Hi, Hj, Hk stand for the quaternion (real) unit and the three quaternionic roots of minus one. To get the bivector algebra, we simply have to multiply Hi, Hk and Hk each by minus 1.

## test 

beta1 <- -Hi
beta2 <- -Hj
beta3 <- -Hk

(beta1 * beta2 == -beta3)
## [1] TRUE
(beta2 * beta3 == -beta1)
## [1] TRUE
(beta3 * beta1 == -beta2)
## [1] TRUE
(beta2 * beta1 == beta3)
## [1] TRUE
(beta3 * beta2 == beta1)
## [1] TRUE
(beta1 * beta3 == beta2)
## [1] TRUE

Now we will also test the two derived algebras generated through choice of lambda = +1 or -1, in particular, testing Christian’s formulas (1) and (2).

My “abeta” is Christian’s \(\sum a_j\beta_j\) (note: he uses Einstein summation convention in his paper; here I write out the summation).

My “abetaL” is Christian’s \(\sum a_k\beta_k(\lambda)\).

My “bbeta” is Christian’s \(\sum b_k\beta_k\).

My “bbetaL” is Christian’s \(\sum b_j\beta_j(\lambda)\).

I do not explicitly test (4) but the reader is welcome to do that also. Christian’s new bivector basis elements \(\beta_j(\lambda)\) are equal to my “betaL1”, “betaL2”, “betaL3.”

lambda <- 1
betaL1 <- beta1 * lambda
betaL2 <- beta2 * lambda
betaL3 <- beta3 * lambda
abeta <- a[1]*beta1 +a[2]*beta2 +a[3]*beta3
abetaL <- a[1]*betaL1 +a[2]*betaL2 +a[3]*betaL3
bbeta <- b[1]*beta1 +b[2]*beta2 +b[3]*beta3
bbetaL <- b[1]*betaL1 +b[2]*betaL2 +b[3]*betaL3

all.equal((-abeta) * abetaL, lambda*H1)
## [1] TRUE
all.equal((bbetaL) * bbeta, -lambda*H1)
## [1] TRUE
lambda <- -1
betaL1 <- beta1 * lambda
betaL2 <- beta2 * lambda
betaL3 <- beta3 * lambda
abeta <- a[1]*beta1 +a[2]*beta2 +a[3]*beta3
abetaL <- a[1]*betaL1 +a[2]*betaL2 +a[3]*betaL3
bbeta <- b[1]*beta1 +b[2]*beta2 +b[3]*beta3
bbetaL <- b[1]*betaL1 +b[2]*betaL2 +b[3]*betaL3

all.equal((-abeta) * abetaL, lambda*H1)
## [1] TRUE
all.equal((bbetaL) * bbeta, -lambda*H1)
## [1] TRUE

Now we will test formulas (5) through (7). I only do 1000 samples, but the reader is welcome to test with a much larger number. For that purpose it would be useful to “vectorise” the code; i.e., do not do an explict “for” loop, but instead create a long vector of values of lambda and use R’s built-in vector and matrix algebra.

Here, I want to the code to be as close as possible to the formulas in the paper.

Nrep <- 1000
S <- 0 * H1
for (i in 1:1000){
  lambda <- sample(c(-1, 1), 1)
  betaL1 <- beta1 * lambda
  betaL2 <- beta2 * lambda
  betaL3 <- beta3 * lambda
  abeta <- a[1]*beta1 +a[2]*beta2 +a[3]*beta3
  abetaL <- a[1]*betaL1 +a[2]*betaL2 +a[3]*betaL3
  bbeta <- b[1]*beta1 +b[2]*beta2 +b[3]*beta3
  bbetaL <- b[1]*betaL1 +b[2]*betaL2 +b[3]*betaL3
  
  A <- (-abeta) * abetaL
  B <- bbetaL * bbeta
  if (lambda == 1) S <- S + A * B else S <- S + B * A
}

(abeta)^{-1} * (S/Nrep) * bbeta   ## E(a, b)
##           [1]
## Re -0.4506891
## i  -0.6240811
## j   0.6382668
## k   0.0042043

We still did not get a real number, and certainly not the real number a.b

What did we get? I calculate the bivectorial product a b and the real number a.b, so we can see what has happened.

abeta * bbeta
##           [1]
## Re -0.4506891
## i  -0.6240811
## j   0.6382668
## k   0.0042043
sum(a * b)
## [1] 0.4506891

This didn’t work. The reason it didn’t work is because of the initial definitions of A and B, which force A and B to equal +/-1 and be opposite in sign to one another, whatever convention we use to define bivector product!