contract() and related functions in the stokes package

Robin K. S. Hankin

contract
## function (K, v, lose = TRUE) 
## {
##     if (is.vector(v)) {
##         out <- Reduce("+", Map("*", apply(index(K), 1, contract_elementary, 
##             v), 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))
## }
## <bytecode: 0x55abc654b598>
## <environment: namespace:stokes>
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)
## }
## <bytecode: 0x55abc65a9ce0>
## <environment: namespace:stokes>

Contractions

Given a \(k\)-form \(\phi\colon V^k\longrightarrow\mathbb{R}\) and a vector \(\mathbf{v}\in V\), the contraction \(\phi_\mathbf{v}\) of \(\phi\) and \(\mathbf{v}\) is a \(k-1\)-form with

\[ \phi_\mathbf{v}\left(\mathbf{v}^1,\ldots,\mathbf{v}^{k-1}\right) = \phi\left(\mathbf{v},\mathbf{v}^1,\ldots,\mathbf{v}^{k-1}\right) \]

provided \(k>1\); if \(k=1\) we specify \(\phi_\mathbf{v}=\phi(\mathbf{v})\). 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 \(\mathbf{v}\); function contract() is the user-friendly front end.

We will start with some simple examples. I will use phi and \(\phi\) 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 \(\phi=dx^1\wedge dx^2\wedge dx^3\wedge dx^4\wedge dx^5\). We have that \(\phi\) 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 \(\phi\) with any vector \(\mathbf{v}\) is thus a \(4\)-form mapping \(V^4\) to the reals with \(\phi_\mathbf{v}\left(\mathbf{v}^1,\mathbf{v}^2,\mathbf{v}^3,\mathbf{v}^4\right)=\phi\left(\mathbf{v},\mathbf{v}^1,\mathbf{v}^2,\mathbf{v}^3,\mathbf{v}^4\right)\). Taking the simplest case first, if \(\mathbf{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 \(V^4\) 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 \(\phi\) understood to be \(\mathbf{v}=(1,0,0,0,0)\)). Now consider \(\mathbf{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: \(\phi_{a\mathbf{v}+b\mathbf{w}}= a\phi_\mathbf{v}+b\phi_\mathbf{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\phi+b\psi)_\mathbf{v}=a\phi_\mathbf{v} + b\psi_\mathbf{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

It is of course possible to contract a contraction. If \(\phi\) is a \(k\)-form, then \(\left(\phi_\mathbf{v}\right)_\mathbf{w}\) is a \(k-2\) form with

\[ \left(\phi_\mathbf{u}\right)_\mathbf{v}\left(\mathbf{w}^1,\ldots,\mathbf{w}^{k-2}\right)=\phi\left(\mathbf{u},\mathbf{v},\mathbf{w}^1,\ldots,\mathbf{w}^{k-2}\right) \]

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
##  2 4 5 6 7  =    2
##  1 2 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
##  1 5 7  =   15
##  4 5 6  =   60
##  1 2 6  =   14
##  1 2 5  =  -10
##  1 5 6  =  -30
##  2 4 7  =   -2
##  2 6 7  =   38
##  1 2 7  =   -1
##  2 5 7  =  -29
##  1 6 7  =  -18
##  2 4 6  =   28
##  4 5 7  =  -30
##  2 4 5  =  -20
##  5 6 7  =  -48
##  2 5 6  =   26
##  4 6 7  =   36

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^8:
##                     val
##  4 6 7 8  =   0.0233312
##  1 3 4 6  =  -0.7893562
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.3797459 -0.3797459 -0.3797459 -0.3797459 -0.3797459
max(jj) - min(jj) # zero to numerical precision
## [1] 2.220446e-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.3797459

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.3797459

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,RHS,LHS-RHS)
## [1]  6.831085e-01  6.831085e-01 -2.220446e-16