Function vector_cross_product() in the stokes package

Robin K. S. Hankin

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)
## [1] -2  1 12

But of course we can work with higher dimensional spaces:

vector_cross_product(matrix(rnorm(30),6,5))
## [1] -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
## [1] 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))
## [1] 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] -1
Mdash <- M[,as.function(sigma)(seq_len(5))]
vector_cross_product(M) + vector_cross_product(Mdash)
## [1]  1.665335e-16  2.664535e-15 -1.554312e-15 -6.661338e-16 -4.440892e-16
## [6]  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
## [1] -2.220446e-16  0.000000e+00 -8.881784e-16  1.776357e-15 -2.664535e-15
## [6]  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)
## [1] 2.220446e-16 0.000000e+00 0.000000e+00 8.881784e-16 0.000000e+00
## [6] 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))
## [1]   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
## [1]  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)
## [1] 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)
)         
## [1] TRUE TRUE TRUE TRUE

References

Hankin, R. K. S. 2020. “Introducing the Permutations R Package.” SoftwareX 11. https://doi.org/10.1016/j.softx.2020.100453.

———. 2022a. “Disordered Vectors in R: Introducing the disordR Package.” https://arxiv.org/abs/2210.03856; arXiv. https://doi.org/10.48550/ARXIV.2210.03856.

———. 2022b. “Stokes’s Theorem in R.” arXiv. https://doi.org/10.48550/ARXIV.2210.17008.

Spivak, M. 1965. Calculus on Manifolds. Addison-Wesley.