# Function vector_cross_product() in the stokes package vector_cross_product
## function (M)
## {
##     n <- nrow(M)
##     stopifnot(n == ncol(M) + 1)
##     (-1)^n * sapply(seq_len(n), function(i) {
##         (-1)^i * det(M[-i, ])
##     })
## }
## <bytecode: 0x560d83916930>
## <environment: namespace:stokes>

# The vector cross product

Spivak (1965), p83, considers the standard vector cross product $$\mathbf{u}\times\mathbf{v}=\det\begin{pmatrix} i & j & k \\ u_1&u_2&u_3\\ v_1&v_2&v_3 \end{pmatrix}$$ and places it in a more general and rigorous context. In a memorable passage, he states:

If $$v_1,\ldots,v_{n-1}\in\mathbb{R}^n$$ and $$\phi$$ is defined by

$\phi(w)=\det\left(\begin{array}{c}v_1\\ \vdots\\ v_{n-1}\\w\end{array}\right)$

then $$\phi\in\Lambda^1\left(\mathbb{R}^n\right)$$; therefore there is a unique $$z\in\mathbb{R}^n$$ such that

$\left\langle w,z\right\rangle=\phi(w)= \det\left(\begin{array}{c}v_1\\ \vdots\\ v_{n-1}\\w\end{array}\right).$

This $$z$$ is denoted $$v_1\times\ldots\times v_{n-1}$$ and is called the cross product of $$v_1,\ldots,v_{n-1}$$.

- Michael Spivak, 1969 (Calculus on Manifolds, Perseus books). Pages 83-84

The reason that $$\mathbf{w}$$ is at the bottom rather than the top is that it ensures that the the $$n$$-tuple $$(\mathbf{v}_1,\ldots,\mathbf{v}_{n-1},\mathbf{w})$$ has positive orientation with respect to the standard basis vectors of $$\mathbb{R}^n$$. In $$\mathbb{R}^3$$ we get the standard elementary mnemonic for $$\mathbf{u}=(u_1,u_2,u_3)$$, $$\mathbf{v}=(v_1,v_2,v_3)$$:

$\mathbf{u}\times\mathbf{v}= \mathrm{det} \begin{pmatrix} i&j&k\\ u_1&u_2&u_3\\ v_1&v_2&v_3 \end{pmatrix}.$

This is (universal) shorthand for the formal definition of the cross product, although sometimes it is better to return to Spivak’s formulation and, writing $$\mathbf{w}=(w_1,w_2,w_3)$$, use the definition directly obtaining

$(\mathbf{u}\times\mathbf{v})\cdot\mathbf{w}= \mathrm{det} \begin{pmatrix} w_1&w_2&w_3\\ u_1&u_2&u_3\\ v_1&v_2&v_3 \end{pmatrix}.$

## R implementation

In the package (Hankin 2022b), R function vector_cross_product() takes a matrix with $$n$$ rows and $$n-1$$ columns: the transpose of the work above. This is because stokes (and R) convention is to interpret columns of a matrix as vectors. If we wanted to take the cross product of $$\mathbf{u}=(5,-2,1)$$ with $$\mathbf{v}=(1,2,0)$$:

(M <- cbind(c(5,-2,1),c(1,2,0)))
##      [,1] [,2]
## [1,]    5    1
## [2,]   -2    2
## [3,]    1    0
vector_cross_product(M)
##  -2  1 12

But of course we can work with higher dimensional spaces:

vector_cross_product(matrix(rnorm(30),6,5))
##  -7.892174  0.927769  1.070595 -5.334296  1.771056 -4.517277

# Verification

## Orientation

We can demonstrate that the function has the correct orientation. We need to ensure that the vectors $$\mathbf{v}_1,\ldots,\mathbf{v}_n,\mathbf{v}_1\times\cdots\times\mathbf{v}_n$$ constitute a right-handed basis:

det(cbind(M,vector_cross_product(M)))>0
##  TRUE

So it is right-handed in this case. Here is a more severe test of the right-handedness::

f <- function(n){
M <- matrix(rnorm(n^2+n),n+1,n)
det(cbind(M,vector_cross_product(M)))>0
}

all(sapply(sample(3:10,100,replace=TRUE),f))
##  TRUE

Above, we see that in each case the vectors are right-handed. We may further verify that the rules for determinants are being obeyed by taking a dot product as follows:

M <- matrix(rnorm(42),7,6)
crossprod(M,vector_cross_product(M))
##               [,1]
## [1,]  4.440892e-16
## [2,] -1.665335e-16
## [3,]  2.664535e-15
## [4,]  8.881784e-16
## [5,] -8.881784e-16
## [6,]  4.440892e-16

Writing $$M=[v_1,\ldots,v_6]$$, $$v_i\in\mathbb{R}^7$$, we see that the dot product $$v_i\cdot v_1\times\cdots\times v_6$$ as implemented by crossprod() vanishes (up to numerical precision), as the determinants in question have two identical columns.

## Immediate properties

Spivak gives the following properties:

$\mathbf{v}_{\sigma(1)}\times\cdots\times\mathbf{v}_{\sigma(n)} = \operatorname{sgn}\sigma\cdot \mathbf{v}_{1}\times\cdots\times\mathbf{v}_{n}$

$\mathbf{v}_{1} \times\cdots\times a\mathbf{v}_i \times\cdots\times \mathbf{v}_{n} = a\cdot \mathbf{v}_{1} \times\cdots\times \mathbf{v}_i \times\cdots\times \mathbf{v}_{n}$

$\mathbf{v}_{1} \times\cdots\times \left(\mathbf{v}_i+{\mathbf{v}'}_i\right) \times\cdots\times \mathbf{v}_{n} = \mathbf{v}_{1} \times\cdots\times \mathbf{v}_i \times\cdots\times \mathbf{v}_{n} + \mathbf{v}_{1} \times\cdots\times {\mathbf{v}'}_i \times\cdots\times \mathbf{v}_{n}$

For the first we use a permutation sigma from the permutations package (Hankin 2020) with a sign of $$-1$$:

M <- matrix(rnorm(30),6,5)
sigma <- as.cycle("(12)(345)")
sgn(sigma)
##  -1
Mdash <- M[,as.function(sigma)(seq_len(5))]
vector_cross_product(M) + vector_cross_product(Mdash)
##   1.665335e-16  2.664535e-15 -1.554312e-15 -6.661338e-16 -4.440892e-16
##   0.000000e+00

Above we see that the the two vector cross products add to zero (up to numerical precision), as they should because sigma is an odd permutation. For the second:

Mdash <- M
Mdash[,3] <- pi*Mdash[,3]
vector_cross_product(Mdash) - vector_cross_product(M) * pi
##  -2.220446e-16  0.000000e+00 -8.881784e-16  1.776357e-15 -2.664535e-15
##   0.000000e+00

Above we see that the second product is $$\pi$$ times the first (to numerical precision), by linearity of the vector cross product. For the third:

M1 <- M
M2 <- M
Msum <- M
v1 <- runif(6)
v2 <- runif(6)
M1[,3] <- v1
M2[,3] <- v2
Msum[,3] <- v1+v2
vector_cross_product(M1) + vector_cross_product(M2) - vector_cross_product(Msum)
##  2.220446e-16 0.000000e+00 0.000000e+00 8.881784e-16 0.000000e+00
##  1.387779e-17

Above we see that the sum of the first two products is equal to that of the third (up to numerical precision), again by linearity of the vector cross product.

## Vector products and Hodge

The cross product has a coordinate-free definition as the Hodge conjugate of the wedge product of its arguments. This is not used in function vector_cross_product() because it is computationally inefficient and (I think) prone to numerical roundoff errors. We may verify that the definitions agree, using a six-dimensional test case:

set.seed(2)
M <- matrix(rnorm(30),6,5)
(ans1 <- vector_cross_product(M))
##    4.431826  -1.966102  -3.344998  -6.853352 -11.879641   7.170485

We can see that vector_cross_product() returns an R vector. To verify that this is correct, we compare the output with the value calculated directly with the wedge product:

(jj <- as.1form(M[,1]) ^ as.1form(M[,2]) ^ as.1form(M[,3]) ^ as.1form(M[,4]) ^ as.1form(M[,5]))
## An alternating linear map from V^5 to R with V=R^6:
##                      val
##  1 2 3 4 5  =   7.170485
##  1 2 4 5 6  =   3.344998
##  1 2 3 5 6  =  -6.853352
##  1 3 4 5 6  =  -1.966103
##  2 3 4 5 6  =  -4.431826
##  1 2 3 4 6  =  11.879641
(ans2 <- hodge(jj))
## An alternating linear map from V^1 to R with V=R^6:
##               val
##  5  =  -11.879641
##  1  =    4.431826
##  2  =   -1.966103
##  4  =   -6.853352
##  3  =   -3.344998
##  6  =    7.170485

Above we see agreement between ans1 and ans2 although the elements might appear in a different order (as per disordR discipline). Actually it is possible to produce the same answer using slightly slicker idiom:

(ans3 <- hodge(Reduce(^,lapply(1:5,function(i){as.1form(M[,i])}))))
## An alternating linear map from V^1 to R with V=R^6:
##               val
##  4  =   -6.853352
##  3  =   -3.344998
##  1  =    4.431826
##  5  =  -11.879641
##  2  =   -1.966103
##  6  =    7.170485

[again note the different order in the output]. Above, we see that the output of vector_cross_product() [ans1] is an ordinary R vector, but the direct result [ans2] is a 1-form. In order to compare these, we first need to coerce ans2 to a 1-form and then subtract:

(diff <- as.1form(ans1) - ans3)
## An alternating linear map from V^1 to R with V=R^6:
##        val
##  1  =    0
##  2  =    0
##  3  =    0
##  5  =    0
##  6  =    0
coeffs(diff)
## A disord object with hash dd72472bf71aceb8700b0990953a3c9cf7326de3 and elements
##   3.552714e-15  1.332268e-15 -4.440892e-16  7.105427e-15 -4.440892e-15
## (in some order)

Above we see that ans1 and ans3 match to within numerical precision.

## Vector cross products in 3 dimensions

Taking Spivak’s definition at face value, we could define the vector cross product $$\mathbf{u}\times\mathbf{v}$$ of three-vectors $$\mathbf{u}$$ and $$\mathbf{v}$$ as a map from the tangent space to the reals, with $$\left(\mathbf{u}\times\mathbf{v}\right)(\mathbf{w})= \left(\mathbf{u}\times\mathbf{v}\right)\cdot\mathbf{w} =\left(I_\mathbf{u}\right)_\mathbf{v}(\mathbf{w})$$, where $$I$$ is the 3-volume element and subscripts refer to contraction. Package idiom for this would be:

function(u,v){contract(volume(3),cbind(u,v))}

However, note that 3D vector cross products are implemented in the package as function vcp3(), which uses different idiom:

body(vcp3)
## {
##     hodge(as.1form(u)^as.1form(v))
## }

This is preferable on the grounds that coercion to a 1-form is explicit. Suppose we wish to take the vector cross product of $$\mathbf{u}=\left(1,4,2\right)^T$$ and $$\mathbf{v}=\left(2,1,5\right)^T$$:

u <- c(1,4,2)
v <- c(2,1,5)
(p <- vcp3(u,v))  # 'p' for (cross) product
## An alternating linear map from V^1 to R with V=R^3:
##        val
##  1  =   18
##  2  =   -1
##  3  =   -7

Above, note the order of the lines is implementation-specific as per disordR discipline (Hankin 2022a), but the form itself should agree with basis vector evaluation given below. Object p is the vector cross product of $$\mathbf{u}$$ and $$\mathbf{v}$$, but is given as a one-form. We can see the mnemonic in operation by coercing p to a function and then evaluating this on the three basis vectors of $$\mathbb{R}^3$$:

ucv  <- as.function(p)
c(i=ucv(ex), j=ucv(ey), k=ucv(ez))
##  i  j  k
## 18 -1 -7

and we see agreement with the mnemonic $$\det\begin{pmatrix}i&j&k\\1&4&2\\2&1&5\end{pmatrix}$$. Further, we may evaluate the triple cross product $$(\mathbf{u}\times\mathbf{v})\cdot\mathbf{w}$$ by evaluating ucv() at $$\mathbf{w}$$.

w <- c(1,-3,2)
ucv(w)
##  7

This shows agreement with the elementary mnemonic $$\det\begin{pmatrix}1&-3&2\\1&4&2\\2&1&5\end{pmatrix}=7$$, as expected from linearity.

## Vector cross product identities

The following identities are standard results:

\begin{aligned} \mathbf{u}\times(\mathbf{v}\times\mathbf{w}) &= \mathbf{v}(\mathbf{w}\cdot\mathbf{u})-\mathbf{w}(\mathbf{u}\cdot\mathbf{v})\\ (\mathbf{u}\times\mathbf{v})\times\mathbf{w} &= \mathbf{v}(\mathbf{w}\cdot\mathbf{u})-\mathbf{u}(\mathbf{v}\cdot\mathbf{w})\\ (\mathbf{u}\times\mathbf{v})\times(\mathbf{u}\times\mathbf{w}) &= (\mathbf{u}\cdot(\mathbf{v}\times\mathbf{w}))\mathbf{u} \\ (\mathbf{u}\times\mathbf{v})\cdot(\mathbf{w}\times\mathbf{x}) &= (\mathbf{u}\cdot\mathbf{w})(\mathbf{v}\cdot\mathbf{x}) - (\mathbf{u}\cdot\mathbf{x})(\mathbf{v}\cdot\mathbf{w}) \end{aligned}

We may verify all four together:

x <- c(-6,5,7)  # u,v,w as before
c(
hodge(as.1form(u) ^ vcp3(v,w))        == as.1form(v*sum(w*u) - w*sum(u*v)),
hodge(vcp3(u,v) ^ as.1form(w))        == as.1form(v*sum(w*u) - u*sum(v*w)),
as.1form(as.function(vcp3(v,w))(u)*u) == hodge(vcp3(u,v) ^ vcp3(u,w))     ,
hodge(hodge(vcp3(u,v)) ^ vcp3(w,x))   == sum(u*w)*sum(v*x) - sum(u*x)*sum(v*w)
)         
##  TRUE TRUE TRUE TRUE