closest_vector for IntegerLattice is broken

[dimpase]

In this example we compute v, the closest vector to the vector u that lies in the integer lattice L. Problem: v is not in L.

sage: from sage.modules.free_module_integer import IntegerLattice
sage: M = matrix(ZZ, [[20957228, -4966110], [9411844, 19625639]])
sage: L = IntegerLattice(M)
sage: u = vector([-423434678248195, -18882583298608161305227077482])
sage: v = L.closest_vector(u)
sage: v in L
False
sage: MGAP = libgap(M)
sage: vGAP = libgap(v)
sage: libgap.SolutionIntMat(MGAP, vGAP)  # independent test that v is not in L
fail

This is a simplified version of ​Taylor Huang's example posted on sage-devel.

This used to work in Sage 8.4, and was broken in 8.5.beta3 by the change from round(x, 6) to x.n(digits=6) in #26648.

The sage-devel thread also reports the membership test to be broken:

from sage.modules.free_module_integer import IntegerLattice
coef = Matrix([-44429982080874270968379672793605458,
    98931650854481334735580708522902113])
bMat = Matrix([[20957228, -4966110], [9411844, 19625639]])
L = IntegerLattice(bMat)
coef*bMat in L  # should be true, but returns False
bMatGAP = libgap(bMat)
v = libgap(coef*bMat)[0]
sol = libgap.SolutionIntMat(bMatGAP, v); sol  # prints
# [ -423434671769860, -18882583298608161305225815339 ]
vector(sol.sage())*bMat == vector(coef*bMat)  # sanity check - prints True

but this is merely a type "error", in the sense that coef is not a vector, but a 1x2 matrix (a different type, and automatic conversion does not happen here). Not sure whether the latter warrants a fix.

[dimpase]

this might be just a bug in membership testing

from sage.modules.free_module_integer import IntegerLattice
coef = Matrix([-44429982080874270968379672793605458, 98931650854481334735580708522902113])
bMat = Matrix([[20957228, -4966110],[ 9411844, 19625639]])
L = IntegerLattice(bMat)
coef*bMat in L

prints False. (cf. continutiation of the thread on sage-devel in ticket description)

[gh-asd00012334]

Hi, I'm the person who came up with these examples. I can be sure that it is not only about membership testing. I've checked my examples with Smith normal form and indeed some of the closest vectors returned is not in L. So both membership and closest_vector are broken.

[dimpase]

we confirm that closest vector is broken, nothing to do with membership test bug

sage: MGAP=libgap(M)
sage: vGAP=libgap(v)
sage: libgap.SolutionIntMat(MGAP,vGAP)
fail

[dimpase]

closest vector got broken somewhere between Sage 8.4 (works) and 8.5 (does not). Membership test already broken in 8.2.

[dimpase]

the heavy lifting is done by fplll/fpylll.

Python 3.7.7 (default, May 30 2020, 01:27:43) 
[GCC 8.3.1 20190518] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from fpylll import CVP, IntegerMatrix
>>> A = IntegerMatrix.from_matrix([[20957228, -4966110], [9411844, 19625639]])
>>> t=(-423434678248195, -18882583298608161305227077482)
>>> v0 = CVP.closest_vector(A, t)
>>> v0
(-423434681577544, -18882583298608161305224294058)

this is the same vector as computed by Sage code in the ticket description, oops. So the bug may be narrowed to fp(y)lll.

Note a comment on ​https://github.com/fplll/fpylll/issues/124 which says that the basis should be LLL-reduced for this to work correctly.

[dimpase]

Thus, in view of comment:11 (corrected) the closest_vector bug is a Sage interface to fp(y)lll bug.

rather, Sage has a different implementation --- which seems to be unchanged since 2015 (i.e. much before Sage 8.5...)

[dimpase]

The containment is "broken", as the 1x2 matrix coef*bMat fails to be coersed to a vector. Namely, with the example in the

sage: vector(coef)*bMat in L
True
sage: vector(coef*bMat) in L
True
sage: coef*bMat in L
False
sage: M=matrix([1,1])
sage: M*bMat in L
False

so this is not really a bug, it's just a conversion issue, absence of automatic conversion.

[slelievre]

This seems suspicious.

sage: from sage.modules.free_module_integer import IntegerLattice
sage: a = vector(ZZ, [20957228, -4966110])
sage: b = vector(ZZ, [9411844, 19625639])
sage: L = IntegerLattice([a, b])
sage: R = L.voronoi_relevant_vectors()
sage: S = [-b, -a, b - a, b, a, a - b]
sage: print([r - s for r, s in zip(R, S)])
[(0, 0), (0, 0), (0, -1), (0, 0), (0, 0), (0, 1)]
sage: print(R[-1]); print(a - b)
(11545384, -24591748)
(11545384, -24591749)

[dimpase]

note that on Sage 8.4 (where closest_vector() still worked) this is

...
sage: print([r - s for r, s in zip(R, S)])
[(-2133540, 44217388), (11545384, -24591749), (-9411844, -19625639), (2133540, -44217388), (0, 0), (-2133540, 44217388)]
sage: print(R[-1]); print(a - b)
(9411844, 19625639)
(11545384, -24591749)

[dimpase]

If fact, R computed in comment:15 is the same as the one in comment:16 (with Sage 8.4), but in a different order.

[slelievre]

The voronoi cell is slightly off. It should be exactly P below.

sage: V = L.voronoi_cell()
sage: ieqs = [(x^2 + y^2, 2*x, 2*y) for x, y in (-b, -a, a-b, b, a, b-a)]
sage: P = Polyhedron(ieqs=ieqs)
sage: VP = V.plot(fill='cyan', alpha=0.5) + P.plot(fill='magenta', alpha=0.5)
sage: for v in V.vertices():
....:     x, y = v.vector()
....:     VP.show(xmin=x-1, xmax=x+1, ymin=y-1, ymax=y+1, figsize=3)

[dimpase]

indeed, on Sage 9.2.beta1

sage: V.vertices()
(A vertex at (-483879947381762629327/57254902852288, 315998472693125499957/28627451426144),
 A vertex at (483879947381762629327/57254902852288, -315998472693125499957/28627451426144),
 A vertex at (5728193175888805065973/458039243775532, 1390651404779390505975/229019621887766),
 A vertex at (-708599510920502622251/229019617181844, 1552002451289776620941/114509808590922),
 A vertex at (708599510920502622251/229019617181844, -1552002451289776620941/114509808590922),
 A vertex at (-5728193175888805065973/458039243775532, -1390651404779390505975/229019621887766))

but on Sage 8.4

sage: V.vertices()
(A vertex at (5728193175888805065973/458039243775532, 1390651404779390505975/229019621887766),
 A vertex at (3871039688862599879323/458039243775532, -2527988039232444116235/229019621887766),
 A vertex at (1417199267595526864965/458039243775532, -3104005018306403526499/229019621887766),
 A vertex at (-5728193175888805065973/458039243775532, -1390651404779390505975/229019621887766),
 A vertex at (-3871039688862599879323/458039243775532, 2527988039232444116235/229019621887766),
 A vertex at (-1417199267595526864965/458039243775532, 3104005018306403526499/229019621887766))

quite a bit different

[dimpase]

I'll try bisecting Sage 8.4 and 8.5. Might take a while...

[slelievre]

Could it be #26648 with this change to the function diamond_cut in src/sage/modules/diamond_cutting.py?

-                hv = [QQ(round(elmt, 6)) for elmt in hv]
+                hv = [QQ(elmt.n(digits=6)) for elmt in hv]

[dimpase]

good catch - this change is certainly bad. cf.

sage: r=5555.55555555555555555
sage: round(r,6)
5555.555556
sage: r.n(digits=6)
5555.56

why it's always 6 that suffices is another question, though

[dimpase]

yes, this fixes this error:

src/sage/modules/diamond_cutting.py

@@ -154 +154 @@
         (A vertex at (2), A vertex at (0))
     """
     # coerce to floats
+    from sage.misc.functional import round
     GM = GM.n()
     C = float(C)
     if verbose:
@@ -223 +224 @@
                 cut_count += 1
                 if verbose:
                     print("\n%d) Cut using normal vector %s" % (cut_count, hv))
-                hv = [QQ(elmt.n(digits=6)) for elmt in hv]
+                hv = [QQ(round(elmt,6)) for elmt in hv]
+                # hv = [QQ(elmt.n(digits=6)) for elmt in hv]
                 inequalities.append(plane_inequality(hv))
 
     if verbose:

[dimpase]

I wonder how many more of these round() calls were butchered this way, or this was an one-off event.

---

New commits:

​d734de1

the correct precision of 6 digits after dot

[dimpase]

The line

hv = [QQ(round(elmt, 6)) for elmt in hv]

still bothers me a lot. Why 6 there? Is it just waiting for an example with bigger numbers to break it again?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1faa6ac

get rid of conversion to floats and round()

[dimpase]

this should be robust now

[slelievre]

Could it be that the conversion to QQ was there for inexact vectors, say vectors over RDF or RR?

Should we add the example in the ticket description as a doctest?

[dimpase]

Replying to slelievre:

Could it be that the conversion to QQ was there for inexact vectors, say vectors over RDF or RR?

for IntegerLattice this looks irrelevant.

probably it was a dodgy old workaround, from pre-PPL days.

Should we add the example in the ticket description as a doctest?

OK, I can do this.

[slelievre]

Replying to dimpase:

Replying to slelievre:

Could it be that the conversion to QQ was there for inexact vectors, say vectors over RDF or RR?

for IntegerLattice this looks irrelevant.

Even though we are working with an integer lattice, doesn't this offer to find the closest lattice point to a non-integer vector?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​685ac9b

adding test from trac #29866

[dimpase]

Replying to slelievre:

Replying to dimpase:

Replying to slelievre:

Could it be that the conversion to QQ was there for inexact vectors, say vectors over RDF or RR?

for IntegerLattice this looks irrelevant.

Even though we are working with an integer lattice, doesn't this offer to find the closest lattice point to a non-integer vector?

the code we're changing only deals with the lattice, not with the vector we are finding a closest lattice point for.

[slelievre]

Right! One last thing, QQ is now imported but not used, can you remove the import?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​99f0df7

remove unused import

[dimpase]

done

[dimpase]

please review

[slelievre]

Thanks!