Function inner() in the
stokes package
function (M)
{
ktensor(spray(expand.grid(seq_len(nrow(M)), seq_len(ncol(M))),
c(M)))
}To cite the stokes package in publications, please use
Hankin (2022b); this function monograph
discusses inner(). Spivak (1965), in a memorable passage,
states:
The reader is already familiar with certain tensors, aside from members of \(V^*\). The first example is the inner product \(\left\langle{,}\right\rangle\in{\mathcal J}^2(\mathbb{R}^n)\). On the grounds that any good mathematical commodity is worth generalizing, we define an inner product on \(V\) to be a 2-tensor \(T\) such that \(T\) is symmetric, that is \(T(v,w)=T(w,v)\) for \(v,w\in V\) and such that \(T\) is positive-definite, that is, \(T(v,v) > 0\) if \(v\neq 0\). We distinguish \(\left\langle{,}\right\rangle\) as the usual inner product on \(\mathbb{R}^n\).
- Michael Spivak, 1969 (Calculus on Manifolds, Perseus books). Page 77
Function inner() returns the inner product corresponding
to a matrix \(M\). Spivak’s definition
requires \(M\) to be positive-definite,
but that is not necessary in the package. The inner product of two
vectors \(\mathbf{x}\) and \(\mathbf{y}\) is usually written \(\left\langle\mathbf{x},\mathbf{y}\right\rangle\)
or \(\mathbf{x}\cdot\mathbf{y}\), but
the most general form would be \(\mathbf{x}^TM\mathbf{y}\). Noting that
inner products are multilinear, that is \(\left\langle\mathbf{x},a\mathbf{y}+b\mathbf{z}\right\rangle=a\left\langle\mathbf{x},\mathbf{y}\right\rangle
+ b\left\langle\mathbf{x},\mathbf{z}\right\rangle\) and \(\left\langle
a\mathbf{x} +
b\mathbf{y},\mathbf{z}\right\rangle=a\left\langle\mathbf{x},\mathbf{z}\right\rangle
+ b\left\langle\mathbf{y},\mathbf{z}\right\rangle\) we see that
the inner product is indeed a multilinear map, that is, a tensor.
We can start with the simplest inner product, the identity matrix:
## A linear map from V^2 to R with V=R^7:
## val
## 6 6 = 1
## 7 7 = 1
## 5 5 = 1
## 3 3 = 1
## 2 2 = 1
## 4 4 = 1
## 1 1 = 1Note how the rows of the tensor appear in arbitrary order, as per
disordR discipline (Hankin 2022a). Verify:
x <- rnorm(7)
y <- rnorm(7)
V <- cbind(x,y)
LHS <- sum(x*y)
RHS <- as.function(inner(diag(7)))(V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)## LHS RHS1 RHS2 RHS3 RHS4 RHS5 RHS6 RHS7
## 5.503805 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
## diff1 diff2 diff3 diff4 diff5 diff6 diff7
## 4.503805 4.503805 4.503805 4.503805 4.503805 4.503805 4.503805Above, we see agreement between \(\mathbf{x}\cdot\mathbf{y}\) calculated
directly [LHS] and using inner()
[RHS]. A more stringent test would be to use a general
matrix:
M <- matrix(rnorm(49),7,7)
f <- as.function(inner(M))
LHS <- quad3.form(M,x,y)
RHS <- f(V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)## Warning in LHS - RHS: Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.## LHS RHS1 RHS2 RHS3 RHS4 RHS5
## -3.41065997 0.35872890 -2.22390027 -0.17262350 0.61824329 -0.05487747
## RHS6 RHS7 RHS8 RHS9 RHS10 RHS11
## -0.79533912 -0.37670272 0.26613736 -1.56378205 -0.83204330 0.50360797
## RHS12 RHS13 RHS14 RHS15 RHS16 RHS17
## -0.41151083 -0.94064916 0.37739565 -0.22732869 0.13333636 -1.16657055
## RHS18 RHS19 RHS20 RHS21 RHS22 RHS23
## -1.23753842 1.08576936 -0.42951311 1.15653700 0.43568330 -1.06559058
## RHS24 RHS25 RHS26 RHS27 RHS28 RHS29
## 0.25222345 0.25014132 -0.05710677 -0.29921512 2.44136463 -0.69095384
## RHS30 RHS31 RHS32 RHS33 RHS34 RHS35
## 0.56074609 -1.28459935 0.83204713 1.15191175 -0.54288826 -0.89192113
## RHS36 RHS37 RHS38 RHS39 RHS40 RHS41
## 0.72675075 -0.27934628 -1.26361438 -0.22426789 -0.23570656 0.80418951
## RHS42 RHS43 RHS44 RHS45 RHS46 RHS47
## 1.23830410 -0.01104548 -0.43331032 -0.64947165 0.99216037 0.04672617
## RHS48 RHS49 diff1 diff2 diff3 diff4
## -0.45278397 1.75790309 -3.76938887 -1.18675970 -3.23803647 -4.02890327
## diff5 diff6 diff7 diff8 diff9 diff10
## -3.35578250 -2.61532086 -3.03395726 -3.67679734 -1.84687792 -2.57861668
## diff11 diff12 diff13 diff14 diff15 diff16
## -3.91426795 -2.99914914 -2.47001081 -3.78805562 -3.18333128 -3.54399634
## diff17 diff18 diff19 diff20 diff21 diff22
## -2.24408943 -2.17312155 -4.49642934 -2.98114687 -4.56719697 -3.84634327
## diff23 diff24 diff25 diff26 diff27 diff28
## -2.34506939 -3.66288342 -3.66080130 -3.35355320 -3.11144486 -5.85202460
## diff29 diff30 diff31 diff32 diff33 diff34
## -2.71970613 -3.97140607 -2.12606062 -4.24270710 -4.56257173 -2.86777172
## diff35 diff36 diff37 diff38 diff39 diff40
## -2.51873885 -4.13741072 -3.13131369 -2.14704559 -3.18639209 -3.17495342
## diff41 diff42 diff43 diff44 diff45 diff46
## -4.21484948 -4.64896408 -3.39961450 -2.97734966 -2.76118833 -4.40282034
## diff47 diff48 diff49
## -3.45738615 -2.95787600 -5.16856306(function quadform::quad3.form() evaluates matrix
products efficiently; quad3.form(M,x,y) returns \(x^TMy\)). Above we see agreement, to within
numerical precision, of the dot product calculated two different ways:
LHS uses quad3.form() and RHS
uses inner(). Of course, we would expect
inner() to be a homomorphism:
M1 <- matrix(rnorm(49),7,7)
M2 <- matrix(rnorm(49),7,7)
g <- as.function(inner(M1+M2))
LHS <- quad3.form(M1+M2,x,y)
RHS <- g(V)
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)## Warning in LHS - RHS: Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.## LHS RHS1 RHS2 RHS3 RHS4 RHS5
## -5.41825307 1.66576786 0.75968355 -1.53293259 1.46734526 -0.04994462
## RHS6 RHS7 RHS8 RHS9 RHS10 RHS11
## -1.92741050 -0.39863150 1.46828030 -0.89621131 0.07688249 -0.13907752
## RHS12 RHS13 RHS14 RHS15 RHS16 RHS17
## -0.78749337 -0.65440110 0.15720910 -0.21915465 0.13854635 -1.20462248
## RHS18 RHS19 RHS20 RHS21 RHS22 RHS23
## 0.58402986 -0.01279603 0.04114776 -1.04996406 -0.75365798 -0.49024403
## RHS24 RHS25 RHS26 RHS27 RHS28 RHS29
## -1.43791924 1.23933413 -1.18786147 -0.09044245 0.69348038 1.72386316
## RHS30 RHS31 RHS32 RHS33 RHS34 RHS35
## 0.38863162 1.38006516 0.33626834 -1.39881284 -0.80081306 -0.37134753
## RHS36 RHS37 RHS38 RHS39 RHS40 RHS41
## 0.57840081 0.29912694 0.08188940 0.56008395 -2.67048744 -0.45132582
## RHS42 RHS43 RHS44 RHS45 RHS46 RHS47
## -1.50748509 -0.78737465 2.76575539 0.29644390 0.27039260 -0.64942291
## RHS48 RHS49 diff1 diff2 diff3 diff4
## -0.13198581 -0.22310848 -7.08402093 -6.17793662 -3.88532048 -6.88559833
## diff5 diff6 diff7 diff8 diff9 diff10
## -5.36830845 -3.49084257 -5.01962157 -6.88653337 -4.52204176 -5.49513556
## diff11 diff12 diff13 diff14 diff15 diff16
## -5.27917555 -4.63075970 -4.76385197 -5.57546217 -5.19909842 -5.55679942
## diff17 diff18 diff19 diff20 diff21 diff22
## -4.21363059 -6.00228293 -5.40545704 -5.45940083 -4.36828901 -4.66459509
## diff23 diff24 diff25 diff26 diff27 diff28
## -4.92800904 -3.98033383 -6.65758720 -4.23039160 -5.32781062 -6.11173345
## diff29 diff30 diff31 diff32 diff33 diff34
## -7.14211623 -5.80688469 -6.79831823 -5.75452141 -4.01944023 -4.61744001
## diff35 diff36 diff37 diff38 diff39 diff40
## -5.04690554 -5.99665388 -5.71738001 -5.50014247 -5.97833702 -2.74776563
## diff41 diff42 diff43 diff44 diff45 diff46
## -4.96692725 -3.91076798 -4.63087842 -8.18400846 -5.71469697 -5.68864567
## diff47 diff48 diff49
## -4.76883016 -5.28626726 -5.19514459Above we see numerical verification of the fact that \(I(M_1+M_2)=I(M_1)+I(M_2)\), by evaluation
at \(\mathbf{x},\mathbf{y}\), again
with LHS using direct matrix algebra and RHS
using inner(). Now, if the matrix is symmetric the
corresponding inner product should also be symmetric:
h <- as.function(inner(M1 + t(M1))) # send inner() a symmetric matrix
LHS <- h(V)
RHS <- h(V[,2:1])
c(LHS=LHS,RHS=RHS,diff=LHS-RHS)## LHS1 LHS2 LHS3 LHS4 LHS5 LHS6
## 3.13428561 1.40593168 -1.13697795 0.64200696 0.09316061 -0.89816379
## LHS7 LHS8 LHS9 LHS10 LHS11 LHS12
## 1.55400137 2.09751303 0.23943528 0.67250666 0.67250666 -0.79581232
## LHS13 LHS14 LHS15 LHS16 LHS17 LHS18
## -1.94457367 -0.79581232 0.47743224 0.51467675 0.37321138 0.47743224
## LHS19 LHS20 LHS21 LHS22 LHS23 LHS24
## 2.09751303 0.19989734 -0.89816379 -0.50732099 1.85227783 0.12450988
## LHS25 LHS26 LHS27 LHS28 LHS29 LHS30
## -1.54290567 -0.53814886 -0.23165064 -0.13010223 -1.13697795 0.79861571
## LHS31 LHS32 LHS33 LHS34 LHS35 LHS36
## 0.12450988 3.13428561 -1.84990868 0.19989734 -1.84990868 1.40593168
## LHS37 LHS38 LHS39 LHS40 LHS41 LHS42
## 1.85227783 0.79861571 0.64200696 -2.87517248 -1.56107844 -0.55155606
## LHS43 LHS44 LHS45 LHS46 LHS47 LHS48
## -1.54290567 0.37321138 1.55400137 -0.53814886 -1.56107844 -0.50732099
## LHS49 RHS1 RHS2 RHS3 RHS4 RHS5
## -0.13010223 3.13428561 1.40593168 -1.13697795 0.64200696 0.09316061
## RHS6 RHS7 RHS8 RHS9 RHS10 RHS11
## -0.89816379 1.55400137 2.09751303 0.23943528 0.67250666 0.67250666
## RHS12 RHS13 RHS14 RHS15 RHS16 RHS17
## -0.79581232 -1.94457367 -0.79581232 0.47743224 0.51467675 0.37321138
## RHS18 RHS19 RHS20 RHS21 RHS22 RHS23
## 0.47743224 2.09751303 0.19989734 -0.89816379 -0.50732099 1.85227783
## RHS24 RHS25 RHS26 RHS27 RHS28 RHS29
## 0.12450988 -1.54290567 -0.53814886 -0.23165064 -0.13010223 -1.13697795
## RHS30 RHS31 RHS32 RHS33 RHS34 RHS35
## 0.79861571 0.12450988 3.13428561 -1.84990868 0.19989734 -1.84990868
## RHS36 RHS37 RHS38 RHS39 RHS40 RHS41
## 1.40593168 1.85227783 0.79861571 0.64200696 -2.87517248 -1.56107844
## RHS42 RHS43 RHS44 RHS45 RHS46 RHS47
## -0.55155606 -1.54290567 0.37321138 1.55400137 -0.53814886 -1.56107844
## RHS48 RHS49 diff1 diff2 diff3 diff4
## -0.50732099 -0.13010223 0.00000000 0.00000000 0.00000000 0.00000000
## diff5 diff6 diff7 diff8 diff9 diff10
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff11 diff12 diff13 diff14 diff15 diff16
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff17 diff18 diff19 diff20 diff21 diff22
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff23 diff24 diff25 diff26 diff27 diff28
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff29 diff30 diff31 diff32 diff33 diff34
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff35 diff36 diff37 diff38 diff39 diff40
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff41 diff42 diff43 diff44 diff45 diff46
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## diff47 diff48 diff49
## 0.00000000 0.00000000 0.00000000Above we see that \(\mathbf{x}^TM\mathbf{y} = \mathbf{y}^TM\mathbf{x}\). Further, a positive-definite matrix should return a positive quadratic form:
M3 <- crossprod(matrix(rnorm(56),8,7)) # 7x7 pos-def matrix
as.function(inner(M3))(kronecker(rnorm(7),t(c(1,1))))>0 # should be TRUE## [1] TRUE FALSE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE
## [25] TRUE FALSE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE FALSE FALSE
## [37] FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE TRUE
## [49] TRUEAbove we see the second line evaluating \(\mathbf{x}^TM\mathbf{x}\) with \(M\) positive-definite, and correctly returning a non-negative value.
Alternating forms
The inner product on an antisymmetric matrix should be alternating:
jj <- matrix(rpois(49,lambda=3.2),7,7)
M <- jj-t(jj) # M is antisymmetric
f <- as.function(inner(M))
LHS <- f(V)
RHS <- -f(V[,2:1]) # NB negative as we are checking for an alternating form
c(LHS=LHS,RHS=RHS,diff=LHS-RHS) ## LHS1 LHS2 LHS3 LHS4 LHS5 LHS6 LHS7 LHS8 LHS9 LHS10 LHS11
## 3 3 2 2 2 2 4 -2 3 -2 1
## LHS12 LHS13 LHS14 LHS15 LHS16 LHS17 LHS18 LHS19 LHS20 LHS21 LHS22
## -1 -2 -1 1 -2 -1 -4 -5 -2 2 -2
## LHS23 LHS24 LHS25 LHS26 LHS27 LHS28 LHS29 LHS30 LHS31 LHS32 RHS1
## 4 -3 -3 1 -3 5 -3 -4 3 2 -3
## RHS2 RHS3 RHS4 RHS5 RHS6 RHS7 RHS8 RHS9 RHS10 RHS11 RHS12
## -3 -2 -2 -2 -2 -4 2 -3 2 -1 1
## RHS13 RHS14 RHS15 RHS16 RHS17 RHS18 RHS19 RHS20 RHS21 RHS22 RHS23
## 2 1 -1 2 1 4 5 2 -2 2 -4
## RHS24 RHS25 RHS26 RHS27 RHS28 RHS29 RHS30 RHS31 RHS32 diff1 diff2
## 3 3 -1 3 -5 3 4 -3 -2 6 6
## diff3 diff4 diff5 diff6 diff7 diff8 diff9 diff10 diff11 diff12 diff13
## 4 4 4 4 8 -4 6 -4 2 -2 -4
## diff14 diff15 diff16 diff17 diff18 diff19 diff20 diff21 diff22 diff23 diff24
## -2 2 -4 -2 -8 -10 -4 4 -4 8 -6
## diff25 diff26 diff27 diff28 diff29 diff30 diff31 diff32
## -6 2 -6 10 -6 -8 6 4Above we see that \(\mathbf{x}^TM\mathbf{y} = -\mathbf{y}^TM\mathbf{x}\) where \(M\) is antisymmetric.
References
disordR Package. Https://arxiv.org/abs/2210.03856; arXiv. https://doi.org/10.48550/ARXIV.2210.03856.