Generate algorithmically all grid points inside a hypercube

$\def\R{\mathbb{R}}\def\Z{\mathbb{Z}}\def\n#1{\|#1\|_\infty}$The problem comes directly from computational mathematics, and can be stated as follows:

Given a regular matrix $M\in\R^{d\times d}$, find effectively all vectors $v\in\Z^d$ such that $\n{Mv}\leq1$, where $\n{Mv}$ is the maximal component of the vector in modulus.

I give an algorithm below, but it's very inefficient, since it has much more steps than the number of grid points. It is still bearable for $d=5$ in my case, but it fails completely for $d=6$, I don't have that much time; moreover, I'd like to make some work on higher dimensions as well.

My current algorithm is based on the following (disclaimer: contains more math): First, compute $\n{M^{-1}}$. Then, we observe that $\n{v}\leq\n{M^{-1}}\n{Mv}\leq\n{M^{-1}}$. This means that we can compute $L=\operatorname{floor} \n{M^{-1}}$ and then try all vectors $v\in\lbrace-L,-L+1,\dots,L-1,L\}^d$; there is exactly $(2L+1)^d$ of them. (And here lies the problem: if $d=6$ and $L=18$, I get $2.5{\rm\small E}9$ iterations, which is orders of magnitude larger than the number of points I think there are.)

• Have you looked into formulating your problem as an integer linear program and then counting solutions? I think there are methods available in IP solvers (e.g. CPLEX) to do this, but I don't know how fast they are in practice. One method for doing so is known as Barvinok's Algorithm. – Cornelius Brand Dec 3 '14 at 15:10
• @C.Brand Thanks, I'll look into it (however, it's been ages since I last seen any LP/IP). By any chance, have you got an idea whether this is present in SageMath? – yo' Dec 4 '14 at 0:12
• As to your concerns, the original formulation of the problem is already (almost) an integer program; the only thing you will have to take care of is the use of the (non-linear) $\infty$-norm. But I have no idea if an algorithm for counting IP solutions is available in SageMath. – Cornelius Brand Dec 4 '14 at 0:16

Here's another way of looking at the problem: You have a lattice generated by the columns of $M$. Use the Lenstra–Lenstra–Lovász (LLL) algorithm to obtain a reduced basis of this lattice. If you replace $M$ by a new matrix formed by the output of LLL, then the columns of $M$ will still generate the same lattice, but the basis vectors will be closer to being orthogonal to one another, and the entries of $M^{-1}$ should have smaller magnitude.

From there, it would also help to bound each component of $v$ separately: i.e., you can bound the $i$th component $|v_i|$ by $\sum_{j=1}^d |(M^{-1})_{ij}|$. (By the way, the bound $\|v\|_\infty \leq \|M^{-1}\|$ is not correct; we need to use the sum of the elements on each row, not the maximum.)

For values of $d$ up to about 30, the LLL algorithm will finish practically instantly. Asymptotically, it takes $O(d^6)$, so it will slow down for very large $d$, but at the same time the number of points we need to check grows exponentially in $d$, so the run-time of the LLL is not really the bottleneck. On the other hand, the savings in the number of points needing to be checked can be enormous. I wrote some GAP code to generate a random regular (stochastic) matrix $M$ and compare the bounds on the components of $v$ that we obtain using the original basis, compared to the LLL-reduced basis (By the way, we do not need to assume that the matrix is regular; I made this restriction only because this was the case in your application):

d:=8;
M:=IdentityMat(d);
for i in [1..d] do
for j in [1..d] do
M[i][j]:=Random([-10^8..10^8]);
od;
M[i]:=M[i]/Sum(M[i]);
od;
L:=LLLReducedBasis(M).basis;
MM:=M^-1*1.0;
LL:=L^-1*1.0;
for i in [1..d] do
for j in [1..d] do
MM[i][j] := MM[i][j]*SignFloat(MM[i][j]);
LL[i][j] := LL[i][j]*SignFloat(LL[i][j]);
od;
od;

Print("Bounds for original basis: ");
ones:=[1..d]*0+1;
v:=MM*ones;
for i in [1..d] do
v[i]:=Int(Floor(v[i]));
Print(v[i]);
Print(" ");
od;
Print("\n(");
Print(Product(v*2+1));
Print(" points to check)\n");

Print("Bounds for LLL basis: ");
v:=LL*ones;
for i in [1..d] do
v[i]:=Int(Floor(v[i]));
Print(v[i]);
Print(" ");
od;
Print("\n(");
Print(Product(v*2+1));
Print(" points to check)\n");

The following output (based on the default random seed, with $d=8$) is not atypical:

Bounds for original basis: 9 23 24 4 23 16 23 4
(258370076349 points to check)
Bounds for LLL basis: 3 3 2 2 3 4 2 3
(2701125 points to check)

Edit: This problem is a special case of the general problem of enumerating lattice points in convex polytopes, which it turns out is a well-studied problem, and there are more efficient algorithms than the one described above. See this this paper for a survey.

• I thought of getting a good base of the grid first, but I'm not sure how efficient it is in high dimensions? (Note: The problem is not currently interesting for me since we managed to get around it completely, but I'll try to keep an eye on this question still.) – yo' Feb 26 '15 at 20:09
• I added some more to my answer, with some information about the efficiency of LLL, as well as example code and output. – Brent Kerby Feb 26 '15 at 22:15