In RobinHankin/weyl: The Weyl Algebra

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:

library(weyl)

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:

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:

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

d1*d2

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

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]

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:

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:

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:

x*(y*z) - (x*y)*z

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

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:

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:

f <- rweyl()
D <- as.der(f)  # D(x) = xf-fx

Then

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

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:

`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
d^2*e

By way of verification:

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]

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)$:

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:

options(polyform = TRUE)
(1+o1)*(1-5*o2)

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:

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:

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:

try(coeffs(W)[1] , silent=FALSE)

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

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:

coeffs(W)[coeffs(W)<7] <- coeffs(W)[coeffs(W)<7] + 100
W

References

RobinHankin/weyl documentation built on April 14, 2025, 11:49 a.m.