Functions contract() and contract_elementary() in the stokes package

contract
function (K, v, lose = TRUE) 
{
    if (is.vector(v)) {
        out <- Reduce("+", Map("*", apply(index(K), 1, contract_elementary, 
            v), elements(coeffs(K))))
    }
    else {
        stopifnot(is.matrix(v))
        out <- K
        for (i in seq_len(ncol(v))) {
            out <- contract(out, v[, i, drop = TRUE], lose = FALSE)
        }
    }
    if (lose) {
        out <- lose(out)
    }
    return(disordR::drop(out))
}
contract_elementary
function (o, v) 
{
    out <- zeroform(length(o) - 1)
    for (i in seq_along(o)) {
        out <- out + (-1)^(i + 1) * v[o[i]] * as.kform(rbind(o[-i]), 
            lose = FALSE)
    }
    return(out)
}

To cite the stokes package in publications, please use Hankin (2022). Given a k-form ϕ: Vk → ℝ and a vector v ∈ V, the contraction ϕv of ϕ and v is a k − 1-form with

ϕv(v1, …, vk − 1) = ϕ(v, v1, …, vk − 1)

provided k > 1; if k = 1 we specify ϕv = ϕ(v). If Spivak (1965) does discuss this, I have forgotten it. Function contract_elementary() is a low-level helper function that translates elementary k-forms with coefficient 1 (in the form of an integer vector corresponding to one row of an index matrix) into its contraction with v; function contract() is the user-friendly front end. We will start with some simple examples. I will use phi and ϕ to represent the same object.

(phi <- as.kform(1:5))
## An alternating linear map from V^5 to R with V=R^5:
##                val
##  1 2 3 4 5  =    1

Thus k = 5 and we have ϕ = dx1 ∧ dx2 ∧ dx3 ∧ dx4 ∧ dx5. We have that ϕ is a linear alternating map with

$$\phi\left(\begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix} \right)=1$$.

The contraction of ϕ with any vector v is thus a 4-form mapping V4 to the reals with ϕv(v1, v2, v3, v4) = ϕ(v, v1, v2, v3, v4). Taking the simplest case first, if v = (1, 0, 0, 0, 0) then

v <- c(1,0,0,0,0)
contract(phi,v)
## An alternating linear map from V^4 to R with V=R^5:
##              val
##  2 3 4 5  =    1

that is, a linear alternating map from V4 to the reals with

$$\phi_\mathbf{v}\left( \begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=1$$.

(the contraction has the first argument of ϕ understood to be v = (1, 0, 0, 0, 0)). Now consider w = (0, 1, 0, 0, 0):

w <- c(0,1,0,0,0)
contract(phi,w)
## An alternating linear map from V^4 to R with V=R^5:
##              val
##  1 3 4 5  =   -1

$$\phi_\mathbf{w}\left( \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=1 \qquad\mbox{or}\qquad \phi_\mathbf{w}\left( \begin{bmatrix}1\\0\\0\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\1\\0\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\1\\0\end{bmatrix}, \begin{bmatrix}0\\0\\0\\0\\1\end{bmatrix}\right)=-1$$.

Contraction is linear, so we may use more complicated vectors:

contract(phi,c(1,3,0,0,0))
## An alternating linear map from V^4 to R with V=R^5:
##              val
##  2 3 4 5  =    1
##  1 3 4 5  =   -3
contract(phi,1:5)
## An alternating linear map from V^4 to R with V=R^5:
##              val
##  1 2 3 4  =    5
##  1 2 4 5  =    3
##  2 3 4 5  =    1
##  1 3 4 5  =   -2
##  1 2 3 5  =   -4

We can check numerically that the contraction is linear in its vector argument: ϕav + bw = aϕv + bϕw.

a <- 1.23
b <- -0.435
v <- 1:5
w <- c(-3, 2.2, 1.1, 2.1, 1.8)

contract(phi,a*v + b*w) == a*contract(phi,v) + b*contract(phi,w)
## [1] TRUE

We also have linearity in the alternating form: (aϕ + bψ)v = aϕv + bψv.

(phi <- rform(2,5))
## An alternating linear map from V^5 to R with V=R^7:
##                val
##  2 3 4 5 7  =   -2
##  1 3 4 6 7  =    1
(psi <- rform(2,5))
## An alternating linear map from V^5 to R with V=R^7:
##                val
##  1 2 3 6 7  =    2
##  2 3 5 6 7  =    1
a <- 7
b <- 13
v <- 1:7
contract(a*phi + b*psi,v) == a*contract(phi,v) + b*contract(psi,v)
## [1] TRUE

Contraction of products

Weintraub (2014) gives us the following theorem, for any k-form ϕ and l-form ψ:

(ϕ ∧ ψ)v = ϕvψ + (−1)kϕ ∧ ψv.

We can verify this numerically with k = 4, l = 5:

phi <- rform(terms=5,k=3,n=9)
psi <- rform(terms=9,k=4,n=9)
v <- sample(1:100,9)
contract(phi^psi,v) ==  contract(phi,v) ^ psi - phi ^ contract(psi,v)
## [1] TRUE

The theorem is verified. We note in passing that the object itself is quite complicated:

summary(contract(phi^psi,v))
## A kform object with 47 terms.  Summary of coefficients: 
## 
## a disord object with hash d6cdb7213e60a9d847f1752a839f18b1de98bc57 
## 
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2943.00  -516.00    48.00    44.47   768.00  2625.00 
## 
## 
## Representative selection of index and coefficients:
## 
## An alternating linear map from V^6 to R with V=R^9:
##                    val
##  1 2 4 6 8 9  =    390
##  1 2 3 5 6 7  =    420
##  1 2 3 4 6 8  =   -840
##  2 5 6 7 8 9  =   1605
##  1 2 3 6 7 9  =    355
##  1 2 3 6 8 9  =  -1200

We may also switch ϕ and ψ, remembering to change the sign:

contract(psi^phi,v) ==  contract(psi,v) ^ phi + psi ^ contract(phi,v)
## [1] TRUE

Repeated contraction

It is of course possible to contract a contraction. If ϕ is a k-form, then (ϕv)w is a k − 2 form with

(ϕu)v(w1, …, wk − 2) = ϕ(u, v, w1, …, wk − 2)

And this is straightforward to realise in the package:

(phi <- rform(2,5))
## An alternating linear map from V^5 to R with V=R^7:
##                val
##  1 4 5 6 7  =   -2
##  2 4 5 6 7  =   -1
u <- c(1,3,2,4,5,4,6)
v <- c(8,6,5,3,4,3,2)
contract(contract(phi,u),v)
## An alternating linear map from V^3 to R with V=R^7:
##             val
##  2 5 6  =    10
##  1 4 7  =     2
##  4 5 7  =    73
##  1 4 5  =    20
##  2 4 7  =     1
##  1 5 6  =    20
##  2 4 6  =   -14
##  1 6 7  =    -2
##  2 6 7  =    -1
##  2 4 5  =    10
##  5 6 7  =    73
##  4 6 7  =   -90
##  1 4 6  =   -28
##  4 5 6  =  -122

But contract() allows us to perform both contractions in one operation: if we pass a matrix M to contract() then this is interpreted as repeated contraction with the columns of M:

M <- cbind(u,v)
contract(contract(phi,u),v) == contract(phi,M)
## [1] TRUE

We can verify directly that the system works as intended. The lines below strip successively more columns from argument V and contract with them:

(o <- kform(spray(t(replicate(2, sample(9,4))), runif(2))))
## An alternating linear map from V^4 to R with V=R^9:
##                     val
##  3 7 8 9  =  -0.1482116
##  1 5 6 7  =   0.4314737
V <- matrix(rnorm(36),ncol=4)
jj <- c(
   as.function(o)(V),
   as.function(contract(o,V[,1,drop=TRUE]))(V[,-1]), # scalar
   as.function(contract(o,V[,1:2]))(V[,-(1:2),drop=FALSE]),
   as.function(contract(o,V[,1:3]))(V[,-(1:3),drop=FALSE]),
   as.function(contract(o,V[,1:4],lose=FALSE))(V[,-(1:4),drop=FALSE])
)
print(jj)
## [1] -0.4992204 -0.4992204 -0.4992204 -0.4992204 -0.4992204
max(jj) - min(jj) # zero to numerical precision
## [1] 2.775558e-16

and above we see agreement to within numerical precision. If we pass three columns to contract() the result is a 0-form:

contract(o,V)
## [1] -0.4992204

In the above, the result is coerced to a scalar which is returned in the form of a disord object; in order to work with a formal 0-form (which is represented in the package as a spray with a zero-column index matrix) we can use the lost=FALSE argument:

contract(o,V,lose=FALSE)
## An alternating linear map from V^0 to R with V=R^0:
##             val
##   =  -0.4992204

thus returning a 0-form. If we iteratively contract a k-dimensional k-form, we return the determinant, and this may be verified as follows:

o <- as.kform(1:5)
V <- matrix(rnorm(25),5,5)
LHS <- det(V)
RHS <- contract(o,V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)
##          LHS          RHS         diff 
## 6.355108e+00 6.355108e+00 1.776357e-15

Above we see agreement to within numerical error.

Contraction from first principles

Suppose we wish to contract ϕ = dxi1 ∧ ⋯ ∧ dxik with vector v = (v1e1, …, vkek). Thus we seek ϕv with ϕv(v1, …, vk − 1) = dxi1 ∧ ⋯ ∧ dxik(v, v1, …, vk − 1). Writing v = v1e1 + ⋯ + ek, we have

where we have exploited linearity. To evaluate this it is easiest and most efficient to express ϕ as $\bigwedge_{j=1}^kdx^{i_j}$ and cycle through the index j, and use various properties of wedge products:

(above, a hat indicates a term’s being omitted). From this, we see that l ∉ L → ϕ = 0 (where L = {i1, …ik} is the index set of ϕ), for all the dxp terms kill el. On the other hand, if l ∈ L we have

Worked example using contract_elementary()

Function contract_elementary() is a bare-bones low-level no-frills helper function that returns ϕv for ϕ an elementary form of the form dxi1 ∧ ⋯ ∧ dxik. Suppose we wish to contract ϕ = dx1 ∧ dx2 ∧ dx5 with vector v = (1, 2, 10, 11, 71)T.

Thus we seek ϕv with ϕv(v1, v2) = dx1 ∧ dx2 ∧ dx5(v, v1, v2). Writing v = v1e1 + ⋯ + v5e5 we have

(above, the zero terms are because the vectors e3 and e4 are killed by dx1 ∧ dx2 ∧ dx5). We can see that the way to evaluate the contraction is to go through the terms of ϕ [that is, dx1, dx2, and dx5] in turn, and sum the resulting expressions:

o <- c(1,2,5)
v <- c(1,2,10,11,71)
(
(-1)^(1+1) * as.kform(o[-1])*v[o[1]] + 
(-1)^(2+1) * as.kform(o[-2])*v[o[2]] +
(-1)^(3+1) * as.kform(o[-3])*v[o[3]]
)
## An alternating linear map from V^2 to R with V=R^5:
##          val
##  1 5  =   -2
##  2 5  =    1
##  1 2  =   71

This is performed more succinctly by contract_elementary():

contract_elementary(o,v)
## An alternating linear map from V^2 to R with V=R^5:
##          val
##  1 5  =   -2
##  2 5  =    1
##  1 2  =   71

The “meat” of contract()

Given a vector v, and kform object K, the meat of contract() would be

Reduce("+", Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))

I will show this in operation with simple but nontrivial arguments.

(K <- as.kform(spray(matrix(c(1,2,3,6,2,4,5,7),2,4,byrow=TRUE),1:2)))
## An alternating linear map from V^4 to R with V=R^7:
##              val
##  2 4 5 7  =    2
##  1 2 3 6  =    1
v <- 1:7

Then the inside bit would be

apply(index(K), 1, contract_elementary, v)
## [[1]]
## An alternating linear map from V^3 to R with V=R^7:
##            val
##  2 4 7  =    5
##  4 5 7  =    2
##  2 5 7  =   -4
##  2 4 5  =   -7
## 
## [[2]]
## An alternating linear map from V^3 to R with V=R^6:
##            val
##  1 2 6  =    3
##  2 3 6  =    1
##  1 3 6  =   -2
##  1 2 3  =   -6

Above we see a two-element list, one for each elementary term of K. We then use base R’s Map() function to multiply each one by the appropriate coefficient:

Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K)))
## [[1]]
## An alternating linear map from V^3 to R with V=R^7:
##            val
##  2 4 5  =  -14
##  2 5 7  =   -8
##  4 5 7  =    4
##  2 4 7  =   10
## 
## [[2]]
## An alternating linear map from V^3 to R with V=R^6:
##            val
##  1 2 3  =   -6
##  1 3 6  =   -2
##  2 3 6  =    1
##  1 2 6  =    3

And finally use Reduce() to sum the terms:

Reduce("+",Map("*", apply(index(K), 1, contract_elementary, v), elements(coeffs(K))))
## An alternating linear map from V^3 to R with V=R^7:
##            val
##  2 4 7  =   10
##  4 5 7  =    4
##  2 5 7  =   -8
##  1 2 3  =   -6
##  2 4 5  =  -14
##  1 3 6  =   -2
##  2 3 6  =    1
##  1 2 6  =    3

However, it might be conceptually easier to use magrittr pipes to give an equivalent definition:

K                                %>%
index                              %>%
apply(1,contract_elementary,v)       %>%
Map("*", ., K %>% coeffs %>% elements) %>%
Reduce("+",.)
## An alternating linear map from V^3 to R with V=R^7:
##            val
##  2 4 7  =   10
##  4 5 7  =    4
##  2 5 7  =   -8
##  1 2 3  =   -6
##  2 4 5  =  -14
##  1 3 6  =   -2
##  2 3 6  =    1
##  1 2 6  =    3

Well it might be clearer to Hadley but frankly YMMV.

References

Hankin, R. K. S. 2022. “Stokes’s Theorem in R.” arXiv. https://doi.org/10.48550/ARXIV.2210.17008.
Spivak, M. 1965. Calculus on Manifolds. Addison-Wesley.
Weintraub, S. T. 2014. Differential Forms: Theory and Practice. Elsevier.