Voronoi cells

Recently, I have implemented the diamond-cutting algorithm to calculate the Voronoi cell of a lattice as described in the paper Computing the Voronoi Cell of a Lattice: The Diamond-Cutting Algorithm by Emanuele Viterbo and Ezio Biglieri.

Basically, it starts with a parallelotope defined by hyperplanes through (and normal to) the halves of the basis vectors (and their negated counterparts). Then it gradually combines basis vectors and cuts off  further half spaces, much like a diamond is cut.

The algorithm considers vectors within a certain distance from the origin. This radius can be specified explicitly to speed things up significantly; if not, a "conservative" safe choice is made.

In the current implementation, I am using Sage's existing Polyhedron class to perform the cutting. It seems reasonably fast.

Emanuele Viterbo was so friendly as to provide his implementation of the algorithm from 1995. Comparing running times, it does not seem like it is worth switching from the Polyhedron class to a custom-tailored version of cutting (i.e. intersecting polyhedrons with half spaces). For example, calculating the Voronoi cell of the following lattice:

2.45 0.00 0.00 0.00 0.00 0.00 
-0.41 2.42 0.00 0.00 0.00 0.00 
-0.41 -0.48 2.37 0.00 0.00 0.00 
-0.41 -0.48 -0.59 2.29 0.00 0.00 
-0.41 -0.48 -0.59 -0.76 2.16 0.00 
-0.41 -0.48 -0.59 -0.76 -1.08 1.87

with ball radius 3 takes 51.4 sec using Emanuele's program but only 13.5 sec with my implementation using Polyhedron on my 2.4 GHz Intel Core 2 Duo Mac. (104 cutting operations are performed in both cases.)

A nice consequence of using the Polyhedron class is that it already provides methods to deal with the resulting Voronoi cell. For example, a representation either by (in)equalities or by vertices/edges/faces can be computed easily.

I had to take special care to deal with non-quadratic basis matrices. To deal with "over-specified" lattices (e.g. redundant basis vectors as in (1,0), (2,0), (0,1)), an LLL-reduction is performed prior to diamond-cutting. (This helps to speed up the algorithm anyway.)

"Non-full" lattices (less basis vectors than dimension of embedded space) are dealt with in a special way, too:

1. To get a quadratic basis matrix, "artificial" basis vectors that represent infinity are added first (they are simply chosen to be longer than any other basis vector);
2. LLL-reduction is performed to get a quadratic basis matrix.
3. After diamond-cutting, all inequalities induced by artificial basis vectors are removed again.

This allows the following:

sage: L = Lattice([[2, 0, 0], [0, 2, 0]])
sage: V = L.voronoi_cell()
sage: V.Hrepresentation()
(An inequality (-1, 0, 0) x + 1 >= 0, An inequality (0, -1, 0) x + 1 >= 0, An inequality (1, 0, 0) x + 1 >= 0, An inequality (0, 1, 0) x + 1 >= 0)

Of course, this is still sort of a workaround—maybe there is a better solution to the problem of non-quadratic basis matrices.