Altfunction (S, give_kform = FALSE)
{
if (is.kform(S)) {
return(S)
}
if (give_kform) {
return(kform(S)/factorial(arity(S)))
}
out <- kill_trivial_rows(S)
if (nrow(index(out)) == 0) {
return(S * 0)
}
ktensor(include_perms(consolidate(out))/factorial(ncol(index(out))))
}To cite the stokes package in publications, please use
Hankin (2022). Spivak (1965), in a memorable passage,
states:
A k-tensor \omega\in{\mathcal J}(V) is called alternating if
\omega(v_1,\ldots,v_i,\ldots,v_j,\ldots,v_k)= -\omega(v_1,\ldots,v_j,\ldots,v_i,\ldots,v_k)\qquad\mbox{for all $v_1,\ldots,v_k\in V$.}
\ldots The set of all alternating k-tensors is clearly a subspace \Lambda^k(V) of {\mathcal J}^k(V). Since it requires considerable work to produce the determinant, it is not surprising that alternating k-tensors are difficult to write down. There is, however, a uniform way of expressing all of them. Recall that the sign of a permutation \sigma, denoted \operatorname{sgn}\sigma, is +1 if \sigma is even and -1 if \sigma is odd. If T\in{\mathcal J}^k(V), we define \operatorname{Alt}(T) by
\operatorname{Alt}(T){\left(v_1,\ldots,v_k\right)}= \frac{1}{k!}\sum_{\sigma\in S_k}\mathrm{sgn}(\sigma)\cdot T{\left(v_{\sigma(1)},\ldots,v_{\sigma(k)}\right)}
where S_k is the set of all permutations of numbers 1 to k.
- Michael Spivak, 1969 (Calculus on Manifolds, Perseus books). Page 78
For example, if k=3 we have
\operatorname{Alt}(T)\left(v_1,v_2,v_3\right)= \frac{1}{6}\left(\begin{array}{cc} +T{\left(v_1,v_2,v_3\right)}& -T{\left(v_1,v_3,v_2\right)}\cr -T{\left(v_2,v_1,v_3\right)}& +T{\left(v_2,v_3,v_1\right)}\cr +T{\left(v_3,v_1,v_2\right)}& -T{\left(v_3,v_2,v_1\right)} \end{array} \right)
In the stokes package, function Alt()
performs this operation on a given k-tensor T.
Idiom is straightforward:
S <- as.ktensor(rbind(c(1,7,8)))*6 # the "6" is to stop rounding error
S## A linear map from V^3 to R with V=R^8:
## val
## 1 7 8 = 6Above, object S=6\phi_1\otimes\phi_7\otimes\phi_8, where
\phi_i\colon\mathbb{R}^8\longrightarrow\mathbb{R},
with \phi_i(e_j)=\delta_{ij} [e_j are basis vectors in \mathbb{R}^8]. We may calculate \operatorname{Alt}(S) in the package using
Alt():
Alt(S)## A linear map from V^3 to R with V=R^8:
## val
## 8 7 1 = -1
## 8 1 7 = 1
## 7 8 1 = 1
## 7 1 8 = -1
## 1 8 7 = -1
## 1 7 8 = 1Above we see that \operatorname{Alt}(S)=\phi_1\otimes\phi_7\otimes\phi_8-\phi_1\otimes\phi_8\otimes\phi_7\pm\cdots
-\phi_8\otimes\phi_7\otimes\phi_1, and we can see the pattern
clearly, observing in passing that the order of the rows in the R object
is arbitrary (as per disordR discipline). Now,
Alt(S) is an alternating form, which is easy to verify, by
making it act on an odd permutation of a set of vectors:
V <- matrix(rnorm(24),ncol=3)
c(as.function(Alt(S))(V),as.function(Alt(S))(V[,c(2,1,3)]))## [1] -0.1645455 0.1645455Observing that \operatorname{Alt} is linear, we might apply it to a more complicated tensor, here with two terms:
(S <- as.ktensor(rbind(c(1,2,4),c(2,2,3)),c(12,1000)))## A linear map from V^3 to R with V=R^4:
## val
## 2 2 3 = 1000
## 1 2 4 = 12Above we have S=12\phi_1\otimes\phi_2\otimes\phi_4+1000\phi_2\otimes\phi_2\otimes\phi_3. Calculating \operatorname{Alt}(S):
S## A linear map from V^3 to R with V=R^4:
## val
## 2 2 3 = 1000
## 1 2 4 = 12
Alt(S)## A linear map from V^3 to R with V=R^4:
## val
## 4 2 1 = -2
## 4 1 2 = 2
## 2 4 1 = 2
## 2 1 4 = -2
## 1 4 2 = -2
## 1 2 4 = 2Note that the 2 2 3 term (with coefficient 1000)
disappears in Alt(S): the even permutations (being
positive) cancel out the odd permutations (being negative) in pairs.
This must be the case for an alternating map by definition. We can see
why this cancellation occurs by considering T=\phi_2\otimes\phi_2\otimes\phi_3, a linear
map corresponding to the [2 2 3] row of S:
T(v_1,v_2,v_3) = (v_1\cdot e_2)(v_2\cdot e_2)(v_1\cdot e_3)x
(where x\in\mathbb{R} is any real number and “\cdot” denotes the dot product), and so
T(v_2,v_1,v_3) = (v_2\cdot e_2)(v_1\cdot e_2)(v_3\cdot e_3)x = T(v_1,v_2,v_3).
Thus if we require T to be
alternating [that is, T(v_2,v_1,v_3) =
-T(v_1,v_2,v_3)], we must have x=0, which is why the [2 2 3]
term with coefficient 1000 vanishes (such terms are killed in the
function by applying kill_trivial_rows()). We can check
that terms with repeated entries are correctly discarded by taking the
\operatorname{Alt} of a tensor all of
whose entries include at least one repeat:
S <- as.ktensor(matrix(c(
3,2,1,1,
1,4,1,4,
1,1,2,3,
7,7,4,7,
1,2,3,3
),ncol=4,byrow=TRUE),1:5)
S## A linear map from V^4 to R with V=R^7:
## val
## 1 2 3 3 = 5
## 7 7 4 7 = 4
## 1 1 2 3 = 3
## 1 4 1 4 = 2
## 3 2 1 1 = 1
Alt(S)## The zero linear map from V^4 to R with V=R^n:
## empty sparse array with 4 columnsWe should verify that Alt() does indeed return an
alternating tensor for complicated tensors (this is trivial
algebraically because \operatorname{Alt}(\cdot) is linear, but it
is always good to check):
S <- rtensor(k=5,n=9)
S## A linear map from V^5 to R with V=R^9:
## val
## 6 1 1 1 2 = 9
## 8 6 7 6 6 = 8
## 7 6 8 7 3 = 7
## 8 2 7 7 7 = 5
## 9 8 6 9 9 = 4
## 6 6 8 9 1 = 3
## 9 2 6 4 7 = 6
## 7 3 3 8 6 = 2
## 9 7 3 4 5 = 1Then we swap columns of V, using both even and odd permutations, to verify that \operatorname{Alt}(S) is in fact alternating:
V_even <- V[,c(1,2,5,3,4)] # an even permutation
V_odd <- V[,c(2,1,5,3,4)] # an odd permutation
V_rep <- V[,c(2,1,5,2,4)] # not a permutation
c(as.function(AS)(V),as.function(AS)(V_even)) # should be identical (even permutation)## [1] 0.226959 0.226959
c(as.function(AS)(V),as.function(AS)(V_odd)) # should differ in sign only (odd permutation)## [1] 0.226959 -0.226959
as.function(AS)(V_rep) # should be zero## [1] -1.01644e-20Further properties of Alt()
In his theorem 4.3, Spivak proves the following statements:
- If T is a tensor, then \operatorname{Alt}(T) is alternating
- If \omega is alternating [he writes \omega\in\Lambda^k(V)] then \operatorname{Alt}(\omega)=\omega
- If T is a tensor, then \operatorname{Alt}(\operatorname{Alt}(T))=\operatorname{Alt}(T).
We have demonstrated the first point above. For the second, we need
to construct a tensor that is alternating, and then show that
Alt() does not change it:
P <- as.ktensor(1+diag(2),c(-7,7))
P## A linear map from V^2 to R with V=R^2:
## val
## 1 2 = 7
## 2 1 = -7
P == Alt(P)## [1] TRUEThe third point, idempotence is also easy:
## [1] TRUEWedge product
(Main article: wedge product). Spivak defines the wedge product as follows. Given alternating forms \omega\in\Lambda^k(V) and \eta\in\Lambda^l(V) we have
\omega\wedge\eta=\frac{(k+l)!}{k!l!}\operatorname{Alt}(\omega\otimes\eta)
So for example:
omega <- as.ktensor(2+diag(2),c(-7,7))
eta <- Alt(ktensor(6*spray(matrix(c(1,2,3,1,4,7,4,5,6),3,3,byrow=TRUE),1:3)))
omega## A linear map from V^2 to R with V=R^3:
## val
## 2 3 = 7
## 3 2 = -7
eta## A linear map from V^3 to R with V=R^7:
## val
## 1 4 7 = 2
## 3 2 1 = -1
## 4 6 5 = -3
## 1 7 4 = -2
## 4 5 6 = 3
## 7 1 4 = 2
## 4 1 7 = -2
## 4 7 1 = 2
## 6 5 4 = -3
## 6 4 5 = 3
## 5 4 6 = -3
## 1 3 2 = -1
## 2 1 3 = -1
## 1 2 3 = 1
## 2 3 1 = 1
## 7 4 1 = -2
## 3 1 2 = 1
## 5 6 4 = 3Above we see that omega is alternating by construction,
and eta is alternating by virtue of Alt().
Thus tensor \omega\wedge\eta is defined
as per Spivak’s definition, and we may calculate it directly:
## An alternating linear map from V^5 to R with V=R^7:
## val
## 2 3 4 5 6 = 2.1
## 1 2 3 4 7 = 1.4
f <- as.function(Alt(omega %X% eta))Observe that the tensor \omega\otimes\eta (which would be returned if
the default argument give_kform=FALSE was sent to
Alt()) is quite long, having 2\cdot 5!=240 nonzero components, which is
why it is not printed in full. We may verify that f() is
alternating by evaluating it on a randomly chosen point in \left(\mathbb{R}^7\right)^5:
## [1] -0.852457 0.852457Above we see a verification that f() is fact
alternating: writing the matrix V in terms of its five
vectors as V=\left[v_1,v_2,v_3,v_4,v_5\right], where
v_i\in\mathbb{R}^7, we see that f{\left(v_1,v_2,v_3,v_4,v_5\right)}=-f{\left(v_2,v_1,v_3,v_4,v_5\right)}.
Further further properties of Alt()
Spivak goes on to prove the following three statements. If S,T are tensors and \omega,\eta,\theta are alternating tensors of arity k,l,m respectively, then
- If \operatorname{Alt}(S)=0, then \operatorname{Alt}(S\otimes T)=\operatorname{Alt}(T\otimes S)=0.
- \operatorname{Alt}(\operatorname{Alt}(\omega\otimes\eta)\otimes\theta) =\operatorname{Alt}(\omega\otimes\eta\otimes\theta)= \operatorname{Alt}(\omega\otimes\operatorname{Alt}(\eta\otimes\theta))
- (\omega\wedge\eta)\wedge\theta = \omega\wedge(\eta\wedge\theta) = \frac{(k+l+m)!}{k!l!m!}\operatorname{Alt}(\omega\otimes\eta\otimes\theta).
Taking the points in turn. Firstly \operatorname{Alt}(S\otimes T)=\operatorname{Alt}(T\otimes S)=0:
(S <- as.ktensor(rbind(c(1,2,3,3),c(1,1,2,3)),1000:1001)) ## A linear map from V^4 to R with V=R^3:
## val
## 1 1 2 3 = 1001
## 1 2 3 3 = 1000
Alt(S) # each row of S includes repeats## The zero linear map from V^4 to R with V=R^n:
## empty sparse array with 4 columns## [1] TRUE TRUEsecondly, \operatorname{Alt}(\operatorname{Alt}(\omega\otimes\eta)\otimes\theta) =\operatorname{Alt}(\omega\otimes\eta\otimes\theta)= \operatorname{Alt}(\omega\otimes\operatorname{Alt}(\eta\otimes\theta)):
omega <- Alt(as.ktensor(rbind(1:3),6))
eta <- Alt(as.ktensor(rbind(4:5),60))
theta <- Alt(as.ktensor(rbind(6:7),14))
omega## A linear map from V^3 to R with V=R^3:
## val
## 3 2 1 = -1
## 3 1 2 = 1
## 2 3 1 = 1
## 2 1 3 = -1
## 1 3 2 = -1
## 1 2 3 = 1
eta## A linear map from V^2 to R with V=R^5:
## val
## 5 4 = -30
## 4 5 = 30
theta## A linear map from V^2 to R with V=R^7:
## val
## 7 6 = -7
## 6 7 = 7
f1 <- as.function(Alt(Alt(omega %X% eta) %X% theta))
f2 <- as.function(Alt(omega %X% eta %X% theta))
f3 <- as.function(Alt(omega %X% Alt(eta %X% theta)))
V <- matrix(rnorm(9*14),ncol=9)
c(f1(V),f2(V),f3(V))## [1] -3.154794 -3.154794 -3.154794Verifying the third identity (\omega\wedge\eta)\wedge\theta = \omega\wedge(\eta\wedge\theta) = \frac{(k+l+m)!}{k!l!m!}\operatorname{Alt}(\omega\otimes\eta\otimes\theta) needs us to coerce from a k-form to a k-tensor:
omega <- rform(2,2,19)
eta <- rform(3,2,19)
theta <- rform(2,2,19)
a1 <- as.ktensor(omega ^ (eta ^ theta))
a2 <- as.ktensor((omega ^ eta) ^ theta)
a3 <- Alt(as.ktensor(omega) %X% as.ktensor(eta) %X% as.ktensor(theta))*90
c(is.zero(a1-a2),is.zero(a1-a3),is.zero(a2-a3))## [1] TRUE TRUE TRUEArgument give_kform
Function Alt() takes a Boolean argument
give_kform. We have been using Alt() with
give_kform taking its default value of FALSE,
which means that it returns an object of class ktensor.
However, an alternating form can be much more efficiently represented as
an object of class kform, and this is returned if
give_kform is TRUE. Here I verify that the two
options return identical objects:
(rand_tensor <- rtensor(k=5,n=9)*120)## A linear map from V^5 to R with V=R^9:
## val
## 7 3 6 4 7 = 120
## 8 1 3 4 4 = 240
## 2 6 4 5 7 = 360
## 8 4 2 8 8 = 840
## 7 1 3 5 1 = 480
## 1 6 4 7 6 = 600
## 6 7 3 3 4 = 1080
## 9 2 3 5 3 = 720
## 6 5 6 5 4 = 960## A ktensor object with 120 terms. Summary of coefficients:
##
## a disord object with hash 866c13390582e2f7e562be3e046abda44be30079
##
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3 -3 0 0 3 3
##
##
## Representative selection of index and coefficients:
##
## A linear map from V^5 to R with V=R^7:
## val
## 7 6 5 4 2 = 3
## 2 6 4 7 5 = -3
## 7 4 2 6 5 = 3
## 2 7 6 4 5 = -3
## 4 6 2 5 7 = -3
## 6 5 2 4 7 = -3Above, we see that S1 is a rather extensive object, with
120 terms. However, if argument give_kform = TRUE is passed
to Alt() we get a kform object which is much
more succinct:
(SA1 <- Alt(rand_tensor,give_kform=TRUE))## An alternating linear map from V^5 to R with V=R^7:
## val
## 2 4 5 6 7 = 3Verification that objects S1 and SA1 are
the same object:
V <- matrix(rnorm(45),ncol=5)
LHS <- as.function(S1)(V)
RHS <- as.function(SA1)(V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)## LHS RHS diff
## 2.386928e+01 2.386928e+01 -1.065814e-14Above we see agreement to within numerical error.