Generate algorithmically all grid points inside a hypercube

Question

$\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.)

Answer 1

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.