Opened 9 months ago

Closed 8 months ago

#29866 closed defect (fixed)

closest_vector for IntegerLattice is broken

Reported by: dimpase Owned by:
Priority: blocker Milestone: sage-9.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:

Status badges

Description (last modified by slelievre)

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.

Change History (39)

comment:1 Changed 9 months ago by mkoeppe

  • Description modified (diff)

comment:2 Changed 9 months ago by slelievre

  • Cc slelievre added
  • Description modified (diff)

comment:3 Changed 9 months ago by 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)

Last edited 9 months ago by dimpase (previous) (diff)

comment:4 Changed 9 months ago by 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.

Last edited 9 months ago by gh-asd00012334 (previous) (diff)

comment:5 Changed 9 months ago by 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

comment:6 Changed 9 months ago by dimpase

  • 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 9 months ago by dimpase

  • Description modified (diff)

comment:8 Changed 9 months ago by dimpase

  • Description modified (diff)

comment:9 Changed 9 months ago by dimpase

  • Description modified (diff)

comment:10 Changed 9 months ago by dimpase

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

Last edited 9 months ago by dimpase (previous) (diff)

comment:11 Changed 9 months ago by 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.

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 LLL-reduced for this to work correctly.

Last edited 9 months ago by dimpase (previous) (diff)

comment:12 Changed 9 months ago by 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...)

Last edited 9 months ago by dimpase (previous) (diff)

comment:13 Changed 9 months ago by dimpase

  • Description modified (diff)

comment:14 Changed 9 months ago by dimpase

  • 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 9 months ago by 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)

comment:16 Changed 9 months ago by 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)

comment:17 Changed 9 months ago by 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.

comment:18 Changed 9 months ago by 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)
Last edited 9 months ago by slelievre (previous) (diff)

comment:19 Changed 9 months ago by 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

comment:20 Changed 9 months ago by dimpase

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

comment:21 Changed 9 months ago by 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]

comment:22 Changed 9 months ago by 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

comment:23 Changed 9 months ago by dimpase

yes, this fixes this error:

  • src/sage/modules/diamond_cutting.py

    a b def diamond_cut(V, GM, C, verbose=False): 
    154154        (A vertex at (2), A vertex at (0))
    155155    """
    156156    # coerce to floats
     157    from sage.misc.functional import round
    157158    GM = GM.n()
    158159    C = float(C)
    159160    if verbose:
    def diamond_cut(V, GM, C, verbose=False): 
    223224                cut_count += 1
    224225                if verbose:
    225226                    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]
    227229                inequalities.append(plane_inequality(hv))
    228230
    229231    if verbose:

comment:24 Changed 9 months ago by dimpase

  • Authors set to Samuel Lelièvre, Dima Pasechnik
  • 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 one-off event.


New commits:

d734de1the correct precision of 6 digits after dot

comment:25 Changed 9 months ago by 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 get it broken again?

Version 1, edited 9 months ago by dimpase (previous) (next) (diff)

comment:26 Changed 9 months ago by git

  • Commit changed from d734de1617b0ab38a60d7bc2e27f8e77107bede4 to 1faa6aceeb6799ea31517c9b84b3e4164153e78d

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

1faa6acget rid of conversion to floats and round()

comment:27 Changed 9 months ago by dimpase

this should be robust now

comment:28 follow-up: Changed 9 months ago by slelievre

  • 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 ; follow-up: Changed 9 months ago by 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.

comment:30 in reply to: ↑ 29 ; follow-up: Changed 9 months ago by 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?

comment:31 Changed 9 months ago by git

  • Commit changed from 1faa6aceeb6799ea31517c9b84b3e4164153e78d to 685ac9bdb3b9451c97707add3118d63b1006f366

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

685ac9badding test from trac #29866

comment:32 in reply to: ↑ 30 Changed 9 months ago by 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.

comment:33 Changed 9 months ago by slelievre

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

comment:34 Changed 9 months ago by git

  • Commit changed from 685ac9bdb3b9451c97707add3118d63b1006f366 to 99f0df74fa4b4b1c836b8b850ceeb516ef012d22

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

99f0df7remove unused import

comment:35 Changed 9 months ago by dimpase

done

comment:36 Changed 9 months ago by dimpase

please review

comment:37 Changed 9 months ago by mkoeppe

  • Reviewers set to Matthias Koeppe
  • Status changed from needs_review to positive_review

comment:38 Changed 9 months ago by slelievre

  • Reviewers changed from Matthias Koeppe to Matthias Koeppe, Samuel Lelièvre

Thanks!

comment:39 Changed 8 months ago by vbraun

  • Branch changed from u/dimpase/modules/corrround to 99f0df74fa4b4b1c836b8b850ceeb516ef012d22
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.