Opened 10 months ago
Closed 10 months ago
#29866 closed defect (fixed)
closest_vector for IntegerLattice is broken
Reported by:  dimpase  Owned by:  

Priority:  blocker  Milestone:  sage9.2 
Component:  number theory  Keywords:  round 
Cc:  malb, slelievre, chapoton  Merged in:  
Authors:  Samuel Lelièvre, Dima Pasechnik  Reviewers:  Matthias Koeppe, Samuel Lelièvre 
Report Upstream:  N/A  Work issues:  
Branch:  99f0df7 (Commits, GitHub, GitLab)  Commit:  99f0df74fa4b4b1c836b8b850ceeb516ef012d22 
Dependencies:  Stopgaps: 
Description (last modified by )
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 sagedevel.
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 sagedevel 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.
Change History (39)
comment:1 Changed 10 months ago by
 Description modified (diff)
comment:2 Changed 10 months ago by
 Cc slelievre added
 Description modified (diff)
comment:3 Changed 10 months ago by
comment:4 Changed 10 months ago by
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.
comment:5 Changed 10 months ago by
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
comment:6 Changed 10 months ago by
 Description modified (diff)
 Priority changed from critical to blocker
 Summary changed from closest_vector for IntegerLattice is broken to closest_vector for IntegerLattice and membership test are broken
comment:7 Changed 10 months ago by
 Description modified (diff)
comment:8 Changed 10 months ago by
 Description modified (diff)
comment:9 Changed 10 months ago by
 Description modified (diff)
comment:10 Changed 10 months ago by
closest vector got broken somewhere between Sage 8.4 (works) and 8.5 (does not). Membership test already broken in 8.2.
comment:11 Changed 10 months ago by
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.
This vector is different from what Sage outputs, namely
(423434671769860, 18882583298608161305225815339)
and may be verified to be in L by calling libgap, say, or just directly:
sage: vector((423434681577544, 18882583298608161305224294058)) in L True
Note a comment on https://github.com/fplll/fpylll/issues/124 which says that the basis should be LLLreduced for this to work correctly.
comment:12 Changed 10 months ago by
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...)
comment:13 Changed 10 months ago by
 Description modified (diff)
comment:14 Changed 10 months ago by
 Description modified (diff)
 Summary changed from closest_vector for IntegerLattice and membership test are broken to closest_vector for IntegerLattice is broken
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.
comment:15 Changed 10 months ago by
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)
comment:16 Changed 10 months ago by
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)
comment:17 Changed 10 months ago by
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.
comment:18 Changed 10 months ago by
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, ab, b, a, ba)] 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=x1, xmax=x+1, ymin=y1, ymax=y+1, figsize=3)
comment:19 Changed 10 months ago by
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
comment:20 Changed 10 months ago by
I'll try bisecting Sage 8.4 and 8.5. Might take a while...
comment:21 Changed 10 months ago by
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]
comment:22 Changed 10 months ago by
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
comment:23 Changed 10 months ago by
yes, this fixes this error:

src/sage/modules/diamond_cutting.py
a b def diamond_cut(V, GM, C, verbose=False): 154 154 (A vertex at (2), A vertex at (0)) 155 155 """ 156 156 # coerce to floats 157 from sage.misc.functional import round 157 158 GM = GM.n() 158 159 C = float(C) 159 160 if verbose: … … def diamond_cut(V, GM, C, verbose=False): 223 224 cut_count += 1 224 225 if verbose: 225 226 print("\n%d) Cut using normal vector %s" % (cut_count, hv)) 226 hv = [QQ(elmt.n(digits=6)) for elmt in hv] 227 hv = [QQ(round(elmt,6)) for elmt in hv] 228 # hv = [QQ(elmt.n(digits=6)) for elmt in hv] 227 229 inequalities.append(plane_inequality(hv)) 228 230 229 231 if verbose:
comment:24 Changed 10 months ago by
 Branch set to u/dimpase/modules/corrround
 Cc chapoton added
 Commit set to d734de1617b0ab38a60d7bc2e27f8e77107bede4
 Status changed from new to needs_review
I wonder how many more of these round()
calls were butchered this way, or this was an oneoff event.
New commits:
d734de1  the correct precision of 6 digits after dot

comment:25 Changed 10 months ago by
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?
comment:26 Changed 10 months ago by
 Commit changed from d734de1617b0ab38a60d7bc2e27f8e77107bede4 to 1faa6aceeb6799ea31517c9b84b3e4164153e78d
Branch pushed to git repo; I updated commit sha1. New commits:
1faa6ac  get rid of conversion to floats and round()

comment:27 Changed 10 months ago by
this should be robust now
comment:28 followup: ↓ 29 Changed 10 months ago by
 Description modified (diff)
 Keywords round added
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?
comment:29 in reply to: ↑ 28 ; followup: ↓ 30 Changed 10 months ago by
Replying to slelievre:
Could it be that the conversion to
RDF
orRR
?
for IntegerLattice this looks irrelevant.
probably it was a dodgy old workaround, from prePPL days.
Should we add the example in the ticket description as a doctest?
OK, I can do this.
comment:30 in reply to: ↑ 29 ; followup: ↓ 32 Changed 10 months ago by
Replying to dimpase:
Replying to slelievre:
Could it be that the conversion to
RDF
orRR
?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 noninteger vector?
comment:31 Changed 10 months ago by
 Commit changed from 1faa6aceeb6799ea31517c9b84b3e4164153e78d to 685ac9bdb3b9451c97707add3118d63b1006f366
Branch pushed to git repo; I updated commit sha1. New commits:
685ac9b  adding test from trac #29866

comment:32 in reply to: ↑ 30 Changed 10 months ago by
Replying to slelievre:
Replying to dimpase:
Replying to slelievre:
Could it be that the conversion to
RDF
orRR
?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 noninteger vector?
the code we're changing only deals with the lattice, not with the vector we are finding a closest lattice point for.
comment:33 Changed 10 months ago by
Right! One last thing, QQ is now imported but not used, can you remove the import?
comment:34 Changed 10 months ago by
 Commit changed from 685ac9bdb3b9451c97707add3118d63b1006f366 to 99f0df74fa4b4b1c836b8b850ceeb516ef012d22
Branch pushed to git repo; I updated commit sha1. New commits:
99f0df7  remove unused import

comment:35 Changed 10 months ago by
done
comment:36 Changed 10 months ago by
please review
comment:37 Changed 10 months ago by
 Reviewers set to Matthias Koeppe
 Status changed from needs_review to positive_review
comment:38 Changed 10 months ago by
 Reviewers changed from Matthias Koeppe to Matthias Koeppe, Samuel Lelièvre
Thanks!
comment:39 Changed 10 months ago by
 Branch changed from u/dimpase/modules/corrround to 99f0df74fa4b4b1c836b8b850ceeb516ef012d22
 Resolution set to fixed
 Status changed from positive_review to closed
this might be just a bug in membership testing
prints
False
. (cf. continutiation of the thread on sagedevel in ticket description)