---
title: "Quantum algebra with R: the weyl package"
author: "Robin K. S. Hankin"
bibliography: weyl.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{The weyl package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



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

To cite this work or the `weyl` package in publications please use
@hankin2022_weyl_arxiv.  In this short document I show how Weyl
algebras are implemented in the `weyl` package.  Consider the vector
space ${\mathcal A}$ of linear operators on univariate functions;
${\mathcal A}$ can be made into an algebra where multiplication
(denoted by juxtaposition) is defined as operator composition
[@coutinho1997].  That is, given operators ${\mathcal O}_1,{\mathcal
O}_2$ and scalar $a$ we define the following compounds:

* $({\mathcal O}_1{\mathcal O}_2)f={\mathcal O}_1({\mathcal O}_2f)$
* $({\mathcal O}_1+{\mathcal O}_2)f={\mathcal O}_1(f)+{\mathcal O}_2(f)$
* $(a{\mathcal O}_1)f=a{\mathcal O}_1(f)$

where $f$ is any univariate function.  Here we consider the algebra
generated by the set $\left\lbrace\partial,x\right\rbrace$ where
$\partial\colon f\longrightarrow f'$ [that is, $(\partial f)(x) =
f'(x)$] and $x\colon f\longrightarrow xf$ [that is, $(xf)(x) =
xf(x)$].  This is known as the (first) Weyl algebra.  We observe that
the Weyl algebra is not commutative: $\partial xf=(xf)'=f+xf'$ but
$x\partial f=xf'$, so $\partial x=x\partial+1$.  The algebra generated
by $\left\lbrace x,\partial\right\rbrace$ will include elements such
as $7\partial + 4\partial x\partial^3 x$, which would map $f$ to $7f'
+ 4\left(x\left(xf\right)'''\right)'$.  It can be shown that any
element of the Weyl algebra can be expressed in the standard form

\[
\sum_i a_i \partial^{p_i}x^{q_i}\]

for real $a_i$ and nonnegative integers $p_i,q_i$.  Converting a
general word to standard form is not straightforward but we have

\[
\partial x^n = x^n\partial +  nx^{n-1}\]

and

\[
\partial^n x = x\partial^n +  n\partial^{n-1}.\]

We can apply these rules recursively to find standard form for
products $(\partial^i x^j)(\partial^l x^m)$.  Alternatively we may
follow @wolf1975 and use the fact that

\[
(\partial^i x^j)(\partial^lx^m)=
\sum_{r=0}^j{j\choose r}{l\choose r}\partial^{i+l-r}x^{j+m-r}.\]

These rules can be used to show that $7\partial + 4\partial
x\partial^3 x$ can be expressed as $7\partial + 12x\partial^2 +
4x^2\partial d^3$, which is in standard form.

## The package in use

The package includes functionality to automate the above calculations.
In particular, package idiom represents the generating elements
$\partial$ and $x$ of the first Weyl algebra as R objects `d` and `x`
respectively.  These may be manipulated with standard arithmetic
operations:

```{r loadlib,echo=TRUE,print=FALSE,results="hide",message=FALSE}
library(weyl)
```		

```{r firstexample}
7*d + 4*x*d^3*x
```

Above, the result of the input is given in standard form.  We see the
terms, one per row, with coefficients in the rightmost column (viz
$7,12,4$).  Thus the first row is $7\partial$, the second is
$12x\partial^2$, and the third is $4x^2\partial^3$.  We may choose to
display the result in symbolic form rather than matrix form:

```{r firstexamplepolyform}
options(polyform=TRUE)
7*d + 4*x*d^3*x
```

which is arguably a more natural representation.  The package allows
one to use R semantics.  For example, consider $d_1=\partial x +
2\partial^3$ and $d_2=3+7\partial -5x^2\partial^2$.  Observing that
$d_1$ and $d_2$ are in standard form, package idiom to create these
operators would be:


```{r defineandprintd1d2}
(d1 <- d*x + 2*d^3)
(d2 <- 3 + 7*d  -5*x^2*d^2)
```

(object `d1` is converted to standard form automatically).  Observe
that, like the `spray` package, the order of the terms is not defined.
We may apply the usual rules of arithmetic to these objects:

```{r showmultiplication}
d1*d2
```

Standard R semantics operate, and it is possible to work with more
complicated expressions:

```{r showinpolyform}
options(polyform=TRUE)
(d1^2 + d2) * (d2 - 3*d1)
```

### Comparison with mathematica

Mathematica can deal with operators and we may compare the two
systems' results for $\partial^2x\partial x^2$:

`In[1] := D[D[x*D[x^2*f[x],x],x],x] // Expand`

`Out[1] := 4 f[x] + 14 x f'[x] + 8 x^2 f''[x] + x^3f'''[x]`

```{r}
x <- weyl(cbind(0,1))
D <- weyl(cbind(1,0))
x^2*D*x*D^2
```

Above, we see agreement between `weyl` and Mathematica, although the
terms are presented in a different order.

# Further Weyl algebras

The package supports arbitrary multivariate Weyl algebras.  Consider:

```{r showmultivariateweyl}
options(polyform=FALSE)  # revert to default print method
set.seed(0)
x <- rweyl()
x
```

Above, object `x` is a member of the operator algebra generated by
$\left\lbrace\partial_x,\partial_y,\partial_z,x,y,z\right\rbrace$.
Object `x` might be expressed as $xz^2\partial_x\partial_z +
3x^2z\partial_x^2\partial_z + 2yz^2\partial_x^2\partial_z^2$ although
as ever the rows are presented in an implementation-dependent order.
We may verify associativity of multiplication:

```{r associativityisnontrivial,cache=TRUE}
x <- rweyl(n=1,d=2)
y <- rweyl(n=2,d=2)
z <- rweyl(n=3,d=2)
options(polyform=TRUE)
x*(y*z)
(x*y)*z
```

Comparing the two results above, we see that they apparently differ.
But the apparent difference is due to the fact that the terms appear
in a different order, a feature that is not algebraically meaningful.
We may verify that the expressions are indeed algebraically identical:

```{r verifyassociativity}
x*(y*z) - (x*y)*z
```

The package can deal with arbitrarily high dimensional Weyl algebras.
For example:

```{r}
(x9 <- rweyl(dim=9))
```

Above we see a member of the ninth Weyl algebra; see how the column
headings no longer use the `x y z` notation and revert to numeric
labels.  Symbolic notation is available but can be difficult to read:

```{r}
options(polyform=TRUE)
x9
options(polyform=FALSE) # revert to default
```

# Derivations

A _derivation_ $D$ of an algebra ${\mathcal A}$ is a linear operator
that satisfies $D(d_1d_2)=d_1D(d_2) + D(d_1)d_2$, for every
$d_1,d_2\in{\mathcal A}$.  If a derivation is of the form $D(d) =
[d,f] = df-fd$ for some fixed $f\in{\mathcal A}$, we say that $D$ is
an _inner_ derivation:

\[
D(d_1d_2) =
d_1d_2f-fd_1d_2 =
d_1d_2f-d_1fd_2 + d_1fd_2-fd_1d_2 = 
d_1(d_2f-fd_2) + (d_1f-fd_1)d_2 =
d_1D(d_2) + D(d_1)d_2
\]

Dirac showed that all derivations are inner derivations for some
$f\in{\mathcal A}$.  The package supports derivations:

```{r define_derivation_D}
f <- rweyl()
D <- as.der(f)  # D(x) = xf-fx
```

Then

```{r show_D_is_a_derivation,cache=TRUE}
d1 <- rweyl()
d2 <- rweyl()
D(d1*d2) == d1*D(d2) + D(d1)*d2
```

# Low-level considerations and generalizations

In the package, the product is customisable.  In general, product
`a*b` [where `a` and `b` are `weyl` objects] is dispatched to the
following sequence of functions:

* `weyl_prod_multivariate_nrow_allcolumns()`
* `weyl_prod_multivariate_onerow_allcolumns()`
* `weyl_prod_multivariate_onerow_singlecolumn()`
* `weyl_prod_univariate_onerow()`
* `weyl_prod_helper3()` [default]

In the above, "univariate" means "generated by $\left\lbrace
x,\partial_x\right\rbrace$" [so the corresponding `spray` object has
_two_ columns]; and "multivariate" means that the algebra is generated
by more than one variable, typically something like $\left\lbrace
x,y,z,\partial_x,\partial_y,\partial_z\right\rbrace$.

The penultimate function `weyl_prod_univariate_onerow()` is sensitive
to option `prodfunc` which specifies the recurrence relation used.
This defaults to `weyl_prod_helper3()`:

```{r}
weyl_prod_helper3
```

Function `weyl_prod_helper3()` follows Wolf.  This gives the
univariate concatenation product $(\partial^a x^b)(\partial^c x^d)$ in
terms of standard generators:

\[
\partial^a x^b \partial^c x^d=\sum_{r=0}^b
r!{b\choose r}{c\choose r}
\partial^{a+c-r}x^{b+d-r}
\]

The package also includes lower-level function `weyl_prod_helper1()`
implementing $\partial^a x^b \partial^c
x^d=\partial^ax^{b-1}\partial^cx^{d+1} +
c\partial^ax^{b-1}\partial^{c-1}x^d$ (together with suitable
bottoming-out).  I expected function `weyl_prod_helper3()` to be much
faster than `weyl_prod_helper1()` but there doesn't seem to be much
difference between the two.

# Generalized commutator relations

We can exploit this package customisability by considering, instead of
$\left\lbrace x,\partial\right\rbrace$, the algebra generated by
$\left\lbrace e,\partial\right\rbrace$, where $e$ maps $f$ to $e^xf$:
if $f$ maps $x$ to $f(x)$, then $ef$ maps $x$ to $e^xf(x)$.  We see
that $\partial e-e\partial=e$.  With this, we can prove that
$\partial^ne=e(1+\partial)^n$ and $e^n\partial=e^n\partial+ne^n$ and,
thus

\[
(e^a\partial^b)(e^c\partial^d)
=e^{a+1}(1+\partial)^be^{c-1}\partial^d
=e^{a}\partial^{b-1}e^{c}\partial^{d+1}+ce^{a}\partial^{b-1}e^{c}\partial^d\]

We may implement this set in package idiom as follows:

```{r}
`weyl_e_prod` <- function(a,b,c,d){
    if(c==0){return(spray(cbind(a,b+d)))}
    if(b==0){return(spray(cbind(a+c,d)))}
    return(
    Recall(a,b-1,c,d+1) +
    c*Recall(a,b-1,c,d)  # cf: c*Recall(a,b-1,c-1,d)) for regular Weyl algebra
    )  
}
```

Then, for example, to calculate $\partial^2e=e(1+2\partial+\partial^2)$:

```{r}
options(prodfunc = weyl_e_prod) 
options(weylvars = "e")  # changes print method
d <- weyl(spray(cbind(0,1)))
e <- weyl(spray(cbind(1,0)))
d*d*e
d^2*e
```

By way of verification:

```{r}
d^5*e == e*(1+d)^5
```

which verifies that indeed $\partial^5e=e(1+\partial)^5$.  Another
verification would be to cross-check with Mathematica, here working with
$\partial e\partial^2e$:

`In[1] := D[Exp[x]*D[D[Exp[x]*f[x],x],x],x]`

`Out[1] := 2E^2x f[x] + 5E^2x f'[x] + 4E^2xf''[x] + E^2x f'''[x]`

```{r}
options(polyform = TRUE)
d*e*d^2*e
```

We can manipulate more complicated expressions too.  Suppose we want to
evaluate $(1+e^2\partial)(1-5e^3\partial^3)$:

```{r}
o1 <- weyl(spray(cbind(2,1)))
o2 <- weyl(spray(cbind(3,3)))
options(polyform = FALSE)
(1+o1)*(1-5*o2)
```

And of course we can display the result in symbolic form:

```{r}
options(polyform = TRUE)
(1+o1)*(1-5*o2)
```


```{r echo=FALSE}
options(polyform = NULL) # restore default print method
```

# Note on use of `disordR` package

The coefficients of a `weyl` object, and the rows of its `spray`
matrix, are stored in an implementation-specific order.  Thus,
extraction and replacement use the `disordR` package
[@hankin2022_disordR_arxiv].  A short example follows in the context
of the `weyl` package; a much more extensive and detailed discussion
is given in the `disordR` vignette and in @hankin2022_disordR_arxiv.
First we create a moderately complicated `weyl` object:

```{r,makearandomweyl}
options(weylvars = NULL)  # revert to default names
(W <- weyl(spray(matrix(c(0,1,1,1,1,2,1,0),2,4),2:3))^2)
```

The coefficients of `W` may be extracted:

```{r,showcoeffsinuse}
coeffs(W)
```

The object returned is a `disord` object.  There is no way to extract
(e.g.) the first coefficient, for the order of the matrix rows is not
defined.  If we try we will get an error:


```{r,showdisorderror,error=TRUE}
try(coeffs(W)[1] , silent=FALSE)
```

However, it is perfectly well defined to ask "give all coefficients
greater than 6":

```{r,ogreaterthantwo}
o <- coeffs(W)
o[o>6]
```

Extraction works as expected.  Using recent improvements in the
`disordR` package, we take all coefficients less than 7 and add 100 to
them:

```{r,showextractionworking}
coeffs(W)[coeffs(W)<7] <- coeffs(W)[coeffs(W)<7] + 100
W
```

## References