---
title: "The `inner()` function in the `stokes` package"
author: "Robin K. S. Hankin"
output: html_vignette
bibliography: stokes.bib
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{inner}
  %\usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
set.seed(0)
library("stokes")
library("quadform")    # needed for quad3.form()
knitr::opts_chunk$set(echo = TRUE)
options(rmarkdown.html_vignette.check_title = FALSE)
knit_print.function <- function(x, ...){dput(x)}
registerS3method(
  "knit_print", "function", knit_print.function,
  envir = asNamespace("knitr")
)
```

```{r out.width='20%', out.extra='style="float:right; padding:10px"',echo=FALSE}
knitr::include_graphics(system.file("help/figures/stokes.png", package = "stokes"))
```

```{r, label=showAlt,comment=""}
inner
```


To cite the `stokes` package in publications, please use
@hankin2022_stokes.  @spivak1965, in a memorable passage, states:

<div class="warning" style='padding:0.1em; background-color:#E9D8FD; color:#69337A'>
<span>
<p>

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 <strong>inner product</strong> on $V$ to be a 2-tensor $T$
such that $T$ is <strong>symmetric</strong>, that is $T(v,w)=T(w,v)$
for $v,w\in V$ and such that $T$ is
<strong>positive-definite</strong>, that is, $T(v,v) > 0$ if $v\neq
0$.  We distinguish $\left\langle{,}\right\rangle$ as the
<strong>usual inner product</strong> on $\mathbb{R}^n$.

</p>
<p style='margin-bottom:1em; margin-right:1em; text-align:right; font-family:Georgia'> <b>- Michael Spivak, 1969</b> <i>(Calculus on Manifolds, Perseus books).  Page 77</i>
</p></span>
</div>

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:

```{r}
inner(diag(7))
```

Note how the rows of the tensor appear in arbitrary order, as per
`disordR` discipline [@hankin2022_disordR].  Verify:

```{r}
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)
```

Above, 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:

```{r}
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)
```

(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:

```{r}
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)
```

Above 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:

```{r}
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)
```

Above we see that $\mathbf{x}^TM\mathbf{y} = \mathbf{y}^TM\mathbf{x}$.
Further, a positive-definite matrix should return a positive quadratic
form:

```{r}
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
```

Above 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:

```{r}
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) 
```

Above we see that $\mathbf{x}^TM\mathbf{y} = -\mathbf{y}^TM\mathbf{x}$
where $M$ is antisymmetric.

# References