8
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Cornelius Brand Dec 3 '14 at 15:10
  • $\begingroup$ @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? $\endgroup$ – yo' Dec 4 '14 at 0:12
  • $\begingroup$ 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. $\endgroup$ – Cornelius Brand Dec 4 '14 at 0:16
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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.) $\endgroup$ – yo' Feb 26 '15 at 20:09
  • $\begingroup$ I added some more to my answer, with some information about the efficiency of LLL, as well as example code and output. $\endgroup$ – Brent Kerby Feb 26 '15 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.