`vector_cross_product()`

and `vcp3()`

in the
`stokes`

package```
function (M)
{
stopifnot(is.matrix(M))
n <- nrow(M)
stopifnot(n == ncol(M) + 1)
(-1)^n * sapply(seq_len(n), function(i) {
(-1)^i * det(M[-i, , drop = FALSE])
})
}
```

`function(u,v){hodge(as.1form(u)^as.1form(v))}`

To cite the `stokes`

package in publications, please use
Hankin (2022b). The *vector cross
product* of two vectors \(\mathbf{u},\mathbf{v}\in\mathbb{R}^3\),
denoted \(\mathbf{u}\times\mathbf{v}\),
is defined in elementary mechanics as \(|\mathbf{u}||\mathbf{v}|\sin(\theta)\,\mathbf{n}\),
where \(\theta\) is the angle between
\(\mathbf{u}\) and \(\mathbf{v}\), and \(\mathbf{n}\) is the unit vector
perpendicular to \(\mathbf{u}\) and
\(\mathbf{v}\) such that \((\mathbf{u},\mathbf{v},\mathbf{n})\) is
positively oriented. Vector cross products find wide applications in
physics, engineering, and computer science. Spivak (1965) 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\cdots\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 for \(\mathbf{w}\) being at the bottom rather than the top is to ensure that 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}. \]

In the `stokes`

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)\):

```
## [,1] [,2]
## [1,] 5 1
## [2,] -2 2
## [3,] 1 0
```

`## [1] -2 1 12`

But of course we can work with higher dimensional spaces:

`## [1] 4.715354 -1.003152 -12.051733 3.459023 -12.902338 -6.943296`

In the case \(n=2\) the vector cross
product becomes a unary operator of a *single* vector \(\left(u,v\right)\in\mathbb{R}^2\),
returning its argument rotated counterclockwise by \(\pi/2\); this case is discussed at the end,
along with \(n=1\).

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:

`## [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:

```
## [,1]
## [1,] 8.881784e-16
## [2,] -1.346145e-15
## [3,] 2.664535e-15
## [4,] 3.552714e-15
## [5,] 8.881784e-16
## [6,] 4.787837e-15
```

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.

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\):

`## [1] -1`

```
## [1] -1.776357e-15 -2.220446e-16 0.000000e+00 -1.776357e-15 1.776357e-15
## [6] -3.330669e-16
```

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:

```
## [1] 3.552714e-15 -8.326673e-16 -1.776357e-15 3.552714e-15 -3.552714e-15
## [6] 2.220446e-16
```

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] -1.776357e-15 4.440892e-16 0.000000e+00 -3.552714e-15 1.776357e-15
## [6] 2.775558e-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.

The cross product has a coordinate-free definition as the Hodge conjugate of the wedge product of its arguments. In \(d\) dimensions:

\[ \mathbf{v}_1\times\cdots\times\mathbf{v}_{d-1}={\star}\left( \mathbf{v}_1\wedge\cdots\wedge\mathbf{v}_{d-1}\right). \]

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:

`## [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:

```
## 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
```

```
## 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:

```
## 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:

```
## An alternating linear map from V^1 to R with V=R^6:
## val
## 1 = 0
## 2 = 0
## 3 = 0
## 5 = 0
## 6 = 0
```

```
## 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.

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:

However, note that 3D vector cross products are implemented in the
package as function `vcp3()`

, which uses different idiom:

`## function(u,v){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\):

```
## 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\):

```
## 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}\).

`## [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.

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`

Function `vector_cross_product()`

takes a matrix with
\(n\) rows and \(n-1\) columns. Here I consider the cases
\(n=2\) and \(n=1\). Firstly, \(n=2\). Going back to Spivak’s definition,
we see that the cross product is a unary operation which takes a single
vector \(\mathbf{v}\in\mathbb{R}^2\);
we might write \({\times}(\mathbf{v})={\times}(v_1,v_2)\).
Formally we define

\[ \phi(w)= \mathrm{det} \begin{pmatrix} v_1&v_2\\ w_1&w_2 \end{pmatrix} \]

and seek a vector \(\mathbf{z}=(z_1,z_2)\in\mathbb{R}^2\) such that \(\left\langle\mathbf{w},\mathbf{z}\right\rangle=\phi(\mathbf{w})\). Thus \(\phi(\mathbf{w})=v_1w_2-v_2w_1\) and we see \(\mathbf{z}=(-v_2,v_1)\). Numerically:

`## [1] -7 4`

We see that the vector cross product of a single vector \(\mathbf{v}\in\mathbb{R}^2\) is vector \(\mathbf{v}\) rotated by \(\pi/2\) counterclockwise; the dot product of \(\mathbf{v}\) with \({\times}{\left(\mathbf{v}\right)}\) is zero.

Now we try the even more peculiar case \(n=1\), corresponding to a matrix with one
row and *zero* columns. Formally, the cross product is a
*nullary* operation which takes zero vectors \(\mathbf{v}\in\mathbb{R}^1\) and returns a
“vector” \(\mathbf{z}\in\mathbb{R}^1\).
The vector cross product in the case \(n=1\) is thus a scalar. Again following
Spivak we see that \(\phi\) is a map
from \(\mathbb{R}^1\) to the reals,
with \(\phi(w_1)=\det(w_1)=w_1\); we
then seek \(z_1\in\mathbb{R}\) such
that \(\phi(w_1)=\left\langle
w_1,z_1\right\rangle\); so \(w_1z_1=w_1\) and then \(z_1=1\). Matrices with zero columns and one
row are easily created in R:

```
##
## [1,]
```

`## structure(logical(0), dim = 1:0)`

Function `vector_cross_product()`

takes such an
argument:

`## [1] 1`

thus returning scalar 1 as intended. Examining the body of
`vector_cross_product`

at the head of the document we see
that the function boils down to returning the determinant of

`## <0 x 0 matrix>`

The determinant of a zero-by-zero matrix is equal to 1 [any zero by zero matrix maps \(\left\lbrace 0\right\rbrace\) to itself and is thus the identity map, which has by definition a determinant of 1]. Numerically:

`## [1] 1`

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.