---
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 x^{q_i}\partial^{p_i}\]
for real $a_i$ and nonnegative integers $p_i,q_i$. Standard form is
desirable operationally: it is a lot easier to evaluate
$x^i\partial^j$ [that is, differentiate $j$ times then multiply by
$x^i$] than to evaluate $\partial^kx^l$ [that is, multiply by $x^l$
then differentiate $k$ times] Converting a general word to standard
form is not straightforward but we have
\[
\partial x^n = x^n\partial + nx^{n-1}.\]
(or, using functional notation,
$\left(x^nf(x)\right)'=x^nf'(x)+nx^{n-1}f(x)$). Noting that operator
composition is associative, we may apply this rule recursively to find
standard form for products $(x^i\partial^j)(x^l\partial^m)$.
Alternatively we may follow @wolf1975 and use the fact that
\[
(x^i\partial^j)(x^l\partial^m)=
\sum_{r=0}^j{j\choose r}{l\choose r}x^{i+l-r}\partial^{j+m-r}.\]
These rules can be used to show, for example, that $7\partial +
4x\partial^3 x$ can be expressed as $7\partial + 12x\partial^2 +
4x^2\partial^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}
options(polyform=FALSE) # default print method
(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