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

[Bouillaguet]

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​

[dimpase]

OK, cool, where is the patch? :)

[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....)

[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 ?

[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 ?

[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".

[Bouillaguet]

I corrected the docstring :)

[dimpase]

Good, positive review.

[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?

[dimpase]

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.

[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.

[dimpase]

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!

[jlopez]

Slight improvement

[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?

[Bouillaguet]

Replying to jlopez:

+1

[jdemeyer]

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

[jdemeyer]

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

[jlopez]

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

[jlopez]

Charles, perhaps you can review my patch yourself?

[Bouillaguet]

@jlopze: sure, give me a couple minutes

[Bouillaguet]

Good to go !