We are not spoilt for choices to work with polynomials of noncommutative variables. Yalmip has undocumented noncommutative functionality for Matlab. While the rest of the toolbox is regularly updated, the noncommutative parts are neglected. NCAlgebra is a package for Mathematica, but its semidefinite programming functions do not seem to be working, and those are the ones we need. SymPy, a Python library for symbolic mathematics, has basic functionality for noncommutative variables, so I took a closer look whether it is able to generate monomials of arbitrary order. Given that this is an open source solution, I would love to see it working. The discussion below uses SymPy 0.7.5.

To start with, we have to import the library:

```
from sympy import *
init_printing()
```

Defining an arbitrary number of noncommutative variables is easy:

```
a=symbols('a0:2', commutative=False)
```

Behold the result:

```
simplify(a[0]*a[1] + a[1]*a[0])
```

Trying to generate the monomials of order two, the situation is less rosy:

```
itermonomials([a[0], a[1]], 2)
```

The polynomial functions do not support noncommutative variables. An additional problem is that the noncommutative variables are not complex:

```
a[0].assumptions0
```

The library will drop an error message if you want to define symbols that are both noncommutative and complex. We are not going to get Hermitian conjugates with these variables. A project from Google Summer of Code 2010 took care of such matters:

```
from sympy.physics.quantum.operator import HermitianOperator, Operator
from sympy.physics.quantum.dagger import Dagger
A = Operator('A')
B = Operator('B')
Dagger(A)
```

```
simplify(A*B + B*A)
```

A new problem is that obtaining an arbitrary number of operators needs extra work:

```
n = 2
X=[]
for i in range(n):
X.append(HermitianOperator('X%s' % i))
```

Sympy will not generate the correct monomials of operators either:

```
itermonomials([X[0], X[1]],2)
```

A not particularly sophisticated function will do the job for me, but it is probably not sufficiently generic:

```
def get_ncmonomials(variables, degree):
if not variables:
return [S.One]
else:
_variables = variables[:]
_variables.insert(0, 1)
ncmonomials = [S.One]
for _ in range(degree):
temp = []
for var in _variables:
for new_var in ncmonomials:
temp.append(var * new_var)
if var != 1 and not var.is_hermitian:
temp.append(Dagger(var) * new_var)
ncmonomials = unique(temp[:])
return ncmonomials
def unique(seq):
seen = {}
result = []
for item in seq:
marker = item
if marker in seen:
continue
seen[marker] = 1
result.append(item)
return result
```

It apparently works:

```
get_ncmonomials(X,2)
```

What is astonishing is just how fast it is compared to Yalmip. Matlab finds it hard to handle problems as small as 512 variables with monomials of order 2. The speedup is in the range of 100–400x. The computational complexity of generating monomials does not change, but practical problems are within range.

The next question is whether SymPy is able to do a QR decomposition of a matrix of noncommutative symbolic variables. That would be an overkill.

[…] also has functions for more advanced applications, such as Hermitian Operators for quantum physics […]