Opened 8 years ago

Closed 8 years ago

#13196 closed enhancement (fixed)

GL(n, GF(q)).random_element() is way too slow for what it does

Reported by: Bouillaguet Owned by: joyner
Priority: major Milestone: sage-5.2
Component: performance Keywords: matrix group, finite field, random element
Cc: malb, dimpase Merged in: sage-5.2.beta1
Authors: Charles Bouillaguet, Javier López Peña Reviewers: Dima Pasechnik, Charles Bouillaguet
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by jdemeyer)

GL(32, GF(2)).random_element() takes more than 10s, which is ridiculous.

Proposed solution : use simple rejection sampling (i.e. generate a random matrix, check if it invertible, if not try again). The result is uniformly distributed amongst invertible matrices.

Before patch:

sage: %time GL(64, GF(2)).random_element() 
5 loops, best of 3: 15.5 s per loop

After patch :

sage: %timeit GL(64, GF(2)).random_element()
625 loops, best of 3: 2ms per loop

Apply gln_gfq_random.patch and trac_13196.improvement.patch

Attachments (2)

gln_gfq_random.patch (1.6 KB) - added by Bouillaguet 8 years ago.
trac_13196.improvement.patch (863 bytes) - added by jlopez 8 years ago.
Slight improvement

Download all attachments as: .zip

Change History (26)

comment:1 follow-up: Changed 8 years ago by dimpase

OK, cool, where is the patch? :)

comment:2 in reply to: ↑ 1 Changed 8 years ago by Bouillaguet

Replying to dimpase:

OK, cool, where is the patch? :)

Sorry, it took me a while as I made a mistake (patched sage/groups/matrix_gps/matrix_group.py first instead of sage/groups/matrix_gps/general_linear.py, which is very wrong. The doctest caught it though....)

comment:3 Changed 8 years ago by Bouillaguet

  • Status changed from new to needs_review

comment:4 follow-up: Changed 8 years ago by dimpase

  • Status changed from needs_review to needs_work

Could you add in the docstring that the probaility to get an invertible element is about 1/3, and so we expect to get the answer in about 3 iterations ?

comment:5 in reply to: ↑ 4 ; follow-up: Changed 8 years ago by Bouillaguet

Replying to dimpase:

Could you add in the docstring that the probaility to get an invertible element is about 1/3, and so we expect to get the answer in about 3 iterations ?

This is not true. The expected number of iterations is asymptotically 1 + 1/q + 2/q^2 + 3/q^3 + 5/q^4 + 7/q^5 + ....

For GF(2) you have 3.5 expected iterations, but for GF(127) you have 1.008 iterations..... So that asymptotically, the number of iterations is O(1), and the complexity of the procedure is O(n^3), regardless of the size of the field.

Is that OK with you ?

comment:6 in reply to: ↑ 5 Changed 8 years ago by dimpase

Replying to Bouillaguet:

Replying to dimpase:

Could you add in the docstring that the probaility to get an invertible element is about 1/3, and so we expect to get the answer in about 3 iterations ?

This is not true. The expected number of iterations is asymptotically 1 + 1/q + 2/q^2 + 3/q^3 + 5/q^4 + 7/q^5 + ....

For GF(2) you have 3.5 expected iterations, but for GF(127) you have 1.008 iterations..... So that asymptotically, the number of iterations is O(1), and the complexity of the procedure is O(n^3), regardless of the size of the field.

Is that OK with you ?

sure. I meant to say "at most 3 iterations".

comment:7 Changed 8 years ago by Bouillaguet

I corrected the docstring :)

comment:8 Changed 8 years ago by dimpase

  • Status changed from needs_work to positive_review

Good, positive review.

comment:9 Changed 8 years ago by Bouillaguet

Modified patch to correct my own stupidity and improve performance a bit

sage: %timeit GL(64, GF(2)).random_element()
625 loops, best of 3: 752 µs per loop

@dimpase: staying positively reviewed?

comment:10 Changed 8 years ago by dimpase

  • Status changed from positive_review to needs_work

patchbot says (and I checked that manually) that the patch (the latest, or previous version, does not matter) causes doctest failures in sage/algebras/group_algebra_new.py and in sage/algebras/group_algebra.py. You can find the log by clicking on the patchbot icon (the, currently, yellow) disc near the top.

Changed 8 years ago by Bouillaguet

comment:11 follow-up: Changed 8 years ago by Bouillaguet

  • Description modified (diff)
  • Status changed from needs_work to needs_review

Patch updated. The bug was that my function returned a matrix, instead of returning an element of the group (basically, it returned g.matrix() instead of returning g). Fix: coerce the result into the group again. It takes a bit of time, and there is likely a workaround, but this is better than nothing.

comment:12 in reply to: ↑ 11 Changed 8 years ago by dimpase

  • Status changed from needs_review to positive_review

Replying to Bouillaguet:

Patch updated. The bug was that my function returned a matrix, instead of returning an element of the group (basically, it returned g.matrix() instead of returning g). Fix: coerce the result into the group again. It takes a bit of time, and there is likely a workaround, but this is better than nothing.

OK!

Changed 8 years ago by jlopez

Slight improvement

comment:13 follow-up: Changed 8 years ago by jlopez

I know this is positively reviewed already, but may I suggest a small improvement? The MatrixSpace.random_element() method creates a zero matrix and then randomizes every entry by calling the M.randomize() method. We can avoid recreating the zero matrix in every iteration by calling directly the randomize() method. For large matrices this gives me a 5-10% speed improvement over your current patch. I don't think the improvement is worth opening a whole new ticket, so can we just apply this on top of the current patch?

Last edited 8 years ago by jlopez (previous) (diff)

comment:14 in reply to: ↑ 13 Changed 8 years ago by Bouillaguet

Replying to jlopez:

+1

comment:15 Changed 8 years ago by jdemeyer

  • Status changed from positive_review to needs_work

If you add a patch, it obviously must be reviewed again.

comment:16 Changed 8 years ago by jdemeyer

  • Component changed from group theory to performance
  • Priority changed from trivial to major
  • Status changed from needs_work to needs_review

comment:17 Changed 8 years ago by jdemeyer

Also, specify which patch(es) should be applied and fill in the Reviewer.

comment:18 Changed 8 years ago by jlopez

  • Description modified (diff)

comment:19 Changed 8 years ago by jlopez

  • Authors changed from Charles Bouillaguet to Charles Bouillaguet, Javier López Peña
  • Reviewers set to Dima Pasechnik

Modified the description and author/reviewer. Since Dima already reviewed Charles' patch, I guess only the improvement patch would need to be reviewed.

patchbot: apply gln_gfq_random.patch trac_13196.improvement.patch

comment:20 Changed 8 years ago by jlopez

Charles, perhaps you can review my patch yourself?

comment:21 Changed 8 years ago by Bouillaguet

@jlopze: sure, give me a couple minutes

comment:22 Changed 8 years ago by Bouillaguet

  • Reviewers changed from Dima Pasechnik to Dima Pasechnik, Charles Bouillaguet
  • Status changed from needs_review to positive_review

Good to go !

comment:23 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:24 Changed 8 years ago by jdemeyer

  • Merged in set to sage-5.2.beta1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.