To cite this work or the weyl
package in publications
please use Hankin (2022b). 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 (Coutinho 1997). That is, given operators \({\mathcal O}_1,{\mathcal O}_2\) and scalar
\(a\) we define the following
compounds:
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 Wolf (1975) 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 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:
library(weyl)
7*d + 4*x*d^3*x
## A member of the Weyl algebra:
## x d val
## 0 1 = 7
## 1 2 = 12
## 2 3 = 4
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:
options(polyform=TRUE)
7*d + 4*x*d^3*x
## A member of the Weyl algebra:
## +7*d +12*x*d^2 +4*x^2*d^3
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:
(d1 <- d*x + 2*d^3)
## A member of the Weyl algebra:
## 1 +x*d +2*d^3
(d2 <- 3 + 7*d -5*x^2*d^2)
## A member of the Weyl algebra:
## 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:
d1*d2
## A member of the Weyl algebra:
## 3 +7*x*d^2 +3*x*d -15*x^2*d^2 -60*x*d^4 -5*x^3*d^3 +7*d +14*d^4 -54*d^3
## -10*x^2*d^5
Standard R semantics operate, and it is possible to work with more complicated expressions:
options(polyform=TRUE)
(d1^2 + d2) * (d2 - 3*d1)
## A member of the Weyl algebra:
## -276*x*d^7 +28*d^7 +5*x^2*d^2 -732*d^6 -636*x*d^4 +28*d^4 -414*d^3
## -63*x^2*d^3 +7*d -20*x^3*d^6 -24*d^9 -70*x*d^2 -20*x^2*d^8 +77*x^3*d^3
## +20*x^4*d^4 -21*x*d +49*d^2 -198*x^2*d^5 +28*x*d^5
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]
x <- weyl(cbind(0,1))
D <- weyl(cbind(1,0))
x^2*D*x*D^2
## A member of the Weyl algebra:
## 4 +x^3*d^3 +8*x^2*d^2 +14*x*d
Above, we see agreement between weyl
and Mathematica,
although the terms are presented in a different order.
The package supports arbitrary multivariate Weyl algebras. Consider:
options(polyform=FALSE) # revert to default print method
set.seed(0)
x <- rweyl()
x
## A member of the Weyl algebra:
## x y z dx dy dz val
## 2 0 1 2 0 1 = 3
## 0 1 2 2 0 1 = 2
## 1 0 2 1 0 1 = 1
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:
x <- rweyl(n=1,d=2)
y <- rweyl(n=2,d=2)
z <- rweyl(n=3,d=2)
options(polyform=TRUE)
x*(y*z)
## A member of the Weyl algebra:
## +6*x*y^2*dx^2*dy +36*x*y^3*dx*dy^3 +2*x^2*y^4*dx*dy^4
## +4*x^2*y^2*dx*dy^5 +36*x*y^2*dx*dy^2 +x*y^4*dx*dy^3 +4*x*y^3*dx*dy^2
## +x^2*y^4*dx^2*dy^3 +12*x^2*y^2*dx*dy^2 +2*x^2*y^2*dx^2*dy^4
## +2*x*y^2*dx*dy^4 +2*x*y^2*dx*dy +2*x^2*y^2*dx^2*dy +4*x^2*y^3*dx^2*dy^2
## +3*x*y^4*dx^2*dy^3 +12*x^2*y^3*dx*dy^3 +6*x*y^4*dx*dy^4
## +12*x*y^3*dx^2*dy^2
(x*y)*z
## A member of the Weyl algebra:
## +2*x*y^2*dx*dy^4 +2*x^2*y^2*dx^2*dy^4 +3*x*y^4*dx^2*dy^3
## +x^2*y^4*dx^2*dy^3 +12*x^2*y^2*dx*dy^2 +x*y^4*dx*dy^3
## +2*x^2*y^4*dx*dy^4 +4*x^2*y^3*dx^2*dy^2 +2*x^2*y^2*dx^2*dy
## +2*x*y^2*dx*dy +4*x*y^3*dx*dy^2 +12*x^2*y^3*dx*dy^3 +6*x*y^4*dx*dy^4
## +12*x*y^3*dx^2*dy^2 +6*x*y^2*dx^2*dy +36*x*y^3*dx*dy^3
## +36*x*y^2*dx*dy^2 +4*x^2*y^2*dx*dy^5
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:
x*(y*z) - (x*y)*z
## A member of the Weyl algebra:
## the NULL multinomial of arity 4
The package can deal with arbitrarily high dimensional Weyl algebras. For example:
(x9 <- rweyl(dim=9))
## A member of the Weyl algebra:
## +3*x1*x2^2*x3^2*x4*x6*x7^2*x8^2*x9^2*d2^2*d7^2*d8*d9
## +2*x1^2*x3*x4*x5*x6*x8*x9*d3^2*d4*d6*d7^2*d8*d9^2
## +x3*x4^2*x5*x6*x7*x8^2*x9^2*d1^2*d3*d6^2*d8^2*d9
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:
options(polyform=TRUE)
x9
## A member of the Weyl algebra:
## +3*x1*x2^2*x3^2*x4*x6*x7^2*x8^2*x9^2*d2^2*d7^2*d8*d9
## +2*x1^2*x3*x4*x5*x6*x8*x9*d3^2*d4*d6*d7^2*d8*d9^2
## +x3*x4^2*x5*x6*x7*x8^2*x9^2*d1^2*d3*d6^2*d8^2*d9
options(polyform=FALSE) # revert to default
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:
f <- rweyl()
D <- as.der(f) # D(x) = xf-fx
Then
d1 <- rweyl()
d2 <- rweyl()
D(d1*d2) == d1*D(d2) + D(d1)*d2
## [1] TRUE
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()
:
weyl_prod_helper3
## function (a, b, c, d)
## {
## f <- function(r) {
## factorial(r) * choose(b, r) * choose(c, r)
## }
## ind <- numeric(0)
## val <- numeric(0)
## for (r in 0:b) {
## ind <- rbind(ind, c(a + c - r, b + d - r))
## val <- c(val, f(r))
## }
## spray(ind, val, addrepeats = TRUE)
## }
## <bytecode: 0x626da4c6e3d8>
## <environment: namespace:weyl>
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.
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:
`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)\):
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
## A member of the Weyl algebra:
## e d val
## 1 0 = 1
## 1 1 = 2
## 1 2 = 1
d^2*e
## A member of the Weyl algebra:
## e d val
## 1 0 = 1
## 1 1 = 2
## 1 2 = 1
By way of verification:
d^5*e == e*(1+d)^5
## [1] TRUE
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]
options(polyform = TRUE)
d*e*d^2*e
## A member of the Weyl algebra:
## +2*e^2 +5*e^2*d +4*e^2*d^2 +e^2*d^3
We can manipulate more complicated expressions too. Suppose we want to evaluate \((1+e^2\partial)(1-5e^3\partial^3)\):
o1 <- weyl(spray(cbind(2,1)))
o2 <- weyl(spray(cbind(3,3)))
options(polyform = FALSE)
(1+o1)*(1-5*o2)
## A member of the Weyl algebra:
## e d val
## 5 3 = -15
## 3 3 = -5
## 5 4 = -5
## 2 1 = 1
## 0 0 = 1
And of course we can display the result in symbolic form:
options(polyform = TRUE)
(1+o1)*(1-5*o2)
## A member of the Weyl algebra:
## 1 -15*e^5*d^3 -5*e^3*d^3 -5*e^5*d^4 +e^2*d
disordR
packageThe 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 (Hankin 2022a). 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 Hankin (2022a). First we
create a moderately complicated weyl
object:
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)
## A member of the Weyl algebra:
## x y dx dy val
## 2 2 4 0 = 9
## 2 2 3 0 = 18
## 0 2 2 2 = 4
## 2 2 2 0 = 9
## 1 2 3 1 = 12
## 1 2 2 0 = 6
## 1 2 3 0 = 6
## 1 2 2 1 = 6
## 0 2 2 1 = 4
The coefficients of W
may be extracted:
coeffs(W)
## A disord object with hash ef3b76da15a19ac6dd3ba83e2ec6b436a0f975f6 and elements
## [1] 9 18 4 9 12 6 6 6 4
## (in some order)
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:
try(coeffs(W)[1] , silent=FALSE)
## Error in .local(x, i, j = j, ..., drop) :
## if using a regular index to extract, must extract each element once and once only (or none of them)
However, it is perfectly well defined to ask “give all coefficients greater than 6”:
o <- coeffs(W)
o[o>6]
## A disord object with hash 908bdb13a2f83dd05aa5c66195a4fe0b570d42f1 and elements
## [1] 9 18 9 12
## (in some order)
Extraction works as expected. Using recent improvements in the
disordR
package, we take all coefficients less than 7 and
add 100 to them:
coeffs(W)[coeffs(W)<7] <- coeffs(W)[coeffs(W)<7] + 100
W
## x y dx dy val
## 0 2 2 1 = 104
## 1 2 2 1 = 106
## 1 2 3 0 = 106
## 1 2 2 0 = 106
## 1 2 3 1 = 12
## 2 2 2 0 = 9
## 0 2 2 2 = 104
## 2 2 3 0 = 18
## 2 2 4 0 = 9
disordR
Package.” https://arxiv.org/abs/2210.03856; arXiv. https://doi.org/10.48550/ARXIV.2210.03856.