Opened 12 years ago

Last modified 8 years ago

#10733 needs_work enhancement

Faster echelon form code for matrix_modn_sparse

Reported by: Gonzalo Tornaría Owned by: William Stein
Priority: major Milestone: sage-pending
Component: linear algebra Keywords:
Cc: Merged in:
Authors: Gonzalo Tornaría, William Stein Reviewers: William Stein, Gonzalo Tornaría
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #10734 Stopgaps:

Status badges

Description (last modified by Jeroen Demeyer)

In the method _echelon_in_place_classical of matrix_modn_sparse there's this comment:

TODO: Implement switching to a dense method when the matrix gets
      dense.

That's what we will do in this ticket, by switching to linbox echelon form code when rows start to get non-dense.

Attachments (6)

trac_10733.01-echelon_better_verbose (2.5 KB) - added by Gonzalo Tornaría 11 years ago.
Improve verbose logging of echelon_in_place_classical
trac_10733.02-echelon_make_faster (4.2 KB) - added by Gonzalo Tornaría 11 years ago.
Make matrix modn sparse echelon form way faster by switching to dense linbox at a good time when the matrix gets too dense
trac_10733.03-documentation (2.4 KB) - added by Gonzalo Tornaría 11 years ago.
Add a test for matrix_modn_sparse._echelon_in_place_classical
trac_10733.04-referee_refactor_to_support_p2.2.patch (1.3 KB) - added by Gonzalo Tornaría 11 years ago.
refactor code to support p=2 case (part 2)
trac_10733.05-doctest_mod2.patch (1.8 KB) - added by Gonzalo Tornaría 11 years ago.
Add testing for echelonize of mod2 sparse matrices
trac_10733.06-fastpath_for_strictly_dense_method.patch (1.7 KB) - added by Gonzalo Tornaría 11 years ago.
Implement a fast path for strictly using a dense method (could be useful for mod2)

Download all attachments as: .zip

Change History (32)

comment:1 Changed 12 years ago by Gonzalo Tornaría

NOTE: The patches in this ticket require #10734 to be applied first.

comment:2 Changed 12 years ago by Gonzalo Tornaría

Status: newneeds_review

comment:3 Changed 12 years ago by William Stein

Status: needs_reviewneeds_work

There are no doctests to illustrate the new functionality added by this patch.

comment:4 Changed 12 years ago by Gonzalo Tornaría

Status: needs_workneeds_review

The patch doesn't add any new functionality.

The function didn't have a doctest before.

I've attached a patch (to be applied in top of the other two) which adds a test for the sparse echelon form. The test takes about 20 seconds, so it's labeled as long.

What the test does is to construct a 1000x1000 sparse matrix with density=1%, and compute the echelon form using different values for "switch_density", and also using the dense matrix echelon form method. Of course everything should give the same echelon form. I'm not sure if it makes sense to test with smaller matrices.

In any case, feel free to add more tests or modify the current one in any way you see fits.

comment:5 Changed 12 years ago by Gonzalo Tornaría

I replaced the patch by a new version. In the new version, I copied the test twice, so it's always run for 100x100 matrices (5% density) and optionally (long) run for 1000x1000 matrices (1% density).

The test for 100x100 matrices takes about 0.1 s, so it's ok. I think it's still a useful test since it runs the algorithm switching to dense at several different columns.

The test for 1000x1000 takes about 25 s (I changed the modulus from 97 to 997) and it's still optional (long). Is that too long for a "long" test?

comment:6 Changed 12 years ago by Gonzalo Tornaría

PS: 500x500 with 1% density takes about 5s, so we may want to replace the 1000x1000 test by a 500x500 test (and still make it long time optional).

comment:7 Changed 12 years ago by William Stein

Description: modified (diff)
Milestone: sage-4.7

comment:8 Changed 12 years ago by William Stein

REVIEW:

Annoyingly, there are issues when the modulus is 2, e.g.:

sage -t  sage/coding/linear_code.py
**********************************************************************
File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/devel/sage-main/sage/coding/linear_code.py", line 1523:
    sage: C.is_permutation_automorphism(g)
Exception raised:
    Traceback (most recent call last):
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/bin/ncadoctest.py", line 1231, in run_one_test
        self.run_one_example(test, example, filename, compileflags)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/bin/sagedoctest.py", line 38, in run_one_example
        OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/bin/ncadoctest.py", line 1172, in run_one_example
        compileflags, 1) in test.globs
      File "<doctest __main__.example_31[10]>", line 1, in <module>
        C.is_permutation_automorphism(g)###line 1523:
    sage: C.is_permutation_automorphism(g)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/lib/python/site-packages/sage/coding/linear_code.py", line 1532, in is_permutation_automorphism
        HGm = H*g.matrix()
      File "element.pyx", line 2282, in sage.structure.element.Matrix.__mul__ (sage/structure/element.c:15874)
      File "coerce.pyx", line 709, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:6368)
      File "action.pyx", line 144, in sage.matrix.action.MatrixMatrixAction._call_ (sage/matrix/action.c:2818)
      File "matrix_modn_sparse.pyx", line 683, in sage.matrix.matrix_modn_sparse.Matrix_modn_sparse.dense_matrix (sage/matrix/matrix_modn_sparse.c:11767)
    TypeError: Cannot convert sage.matrix.matrix_mod2_dense.Matrix_mod2_dense to sage.matrix.matrix_modn_dense.Matrix_modn_dense
**********************************************************************

and

sage -t  sage/homology/cell_complex.py
**********************************************************************
File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/devel/sage-main/sage/homology/cell_complex.py", line 456:
    sage: P.homology(base_ring=GF(2))
Exception raised:
    Traceback (most recent call last):
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/bin/ncadoctest.py", line 1231, in run_one_test
        self.run_one_example(test, example, filename, compileflags)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/bin/sagedoctest.py", line 38, in run_one_example
        OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/bin/ncadoctest.py", line 1172, in run_one_example
        compileflags, 1) in test.globs
      File "<doctest __main__.example_15[4]>", line 1, in <module>
        P.homology(base_ring=GF(Integer(2)))###line 456:
    sage: P.homology(base_ring=GF(2))
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/lib/python/site-packages/sage/homology/cell_complex.py", line 552, in homology
        dimensions=dims, **kwds)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/lib/python/site-packages/sage/homology/delta_complex.py", line 662, in chain_complex
        return ChainComplex(data=differentials, degree=-1, **kwds)
      File "/tmp/wstein/sage-4.7.alpha3-sage.math.washington.edu-x86_64-Linux/local/lib/python/site-packages/sage/homology/chain_complex.py", line 462, in __init__
        raise TypeError, "The differentials d_{%s} and d_{%s} are not compatible: their product is not defined." % (n, n+degree)
    TypeError: The differentials d_{1} and d_{0} are not compatible: their product is not defined.

The result of running the whole test suite:


Here is another style comment. In Cython now doing

for c from 0 <= c < self._ncols:

is *identical* to doing:

for c in range(self._ncols):

Otherwise, everything looks great, and this is going to be *extremeley* useful for my research!!

comment:9 Changed 12 years ago by William Stein

Status: needs_reviewneeds_work

comment:10 Changed 12 years ago by William Stein

Oh, here are the files with failures:

 
----------------------------------------------------------------------

The following tests failed:

	sage -t  devel/sage-main/sage/modular/modsym/ambient.py # 4 doctests failed
	sage -t  devel/sage-main/sage/matrix/matrix_space.py # 1 doctests failed
	sage -t  devel/sage-main/sage/coding/linear_code.py # 2 doctests failed
	sage -t  devel/sage-main/sage/rings/polynomial/pbori.pyx # 1 doctests failed
	sage -t  devel/sage-main/sage/homology/cell_complex.py # 2 doctests failed
	sage -t  devel/sage-main/sage/homology/examples.py # 2 doctests failed
	sage -t  devel/sage-main/sage/homology/delta_complex.py # 3 doctests failed
----------------------------------------------------------------------
Timings have been updated.
Total time for all tests: 328.6 seconds

It's all p=2 weirdness, I suspect.

comment:11 Changed 12 years ago by William Stein

Comment: I haven't run the full test suite yet -- will report back later...

comment:12 Changed 12 years ago by William Stein

Status: needs_workneeds_review

Test suite run and it fully passes for me. Gonzalo (say) needs to review patch 4 and then this can get a positive review.

comment:13 Changed 12 years ago by Gonzalo Tornaría

Status: needs_reviewpositive_review

I reviewed the new patch and ran the test suite, which passes for me as well.

The problematic method was Matrix_modn_sparse.set_block_unsafe(..., Matrix_modn_dense m) which is meant to set a block of a sparse matrix given by a dense matrix. But for mod 2, the dense matrices are implemented by a different class (Matrix_mod2_dense) so calls to set_block_unsafe will fail. What the reviewer patch does is to factor out a common superclass for both Matrix_mod_dense which can be used in the implementation of set_block_unsafe.

I noticed that computing echelon form for a sparse mod2 matrix is much slower than converting the matrix to dense and using echelon form (which uses M4RI and is quite fast). Thus, it could make sense to reimplement sparse mod2 echelon form via conversion to dense. However, note that conversion from mod2 dense to mod2 sparse is quite slow currently.

In spite of the above comment, I'm giving positive review to patch 4 based on:

  • the code works properly for mod2 sparse matrices, so no functionality regression.
  • it runs definitely faster than what we had before for sparse matrices, so no performance regression.

Since this ticket gives significant improvements for larger modulus, it's important that we merge it even if it's not the best we could do for the modulus 2.

comment:14 Changed 12 years ago by Gonzalo Tornaría

For the record, this is what I think could be used for mod2 sparse via dense:

   dense = self.dense_matrix()
   dense.echelonize()
   self.set_block_unsafe(0, 0, dense)

The last step does the conversion from dense to sparse, and it's way faster than calling dense.matrix_sparse(), but it's still slow.

comment:15 Changed 12 years ago by Gonzalo Tornaría

I added two more patches to the sequence:

  • trac_10733.05: adds a doctest for the case of mod2, because it is using different code it's worth testing separately
  • trac_10733.06: adds a special case to the method for the case when "switch_density" is 0.

The latter basically implements what I suggested in my last comment as part of the _echelon_in_place_classical() method. As is now, the default for "switch_density" is always 0.1 so this code would not be normally used. But if one explicitly uses "switch_density=0", then the method will run faster after this patch. And for the case of mod 2, it seems that it could indeed be faster to use switch_density=0 instead of the default 0.1.

I guess these two patches need to go back to review, although trac doesn't give me the option of changing the status to needs_review...

comment:16 Changed 12 years ago by Gonzalo Tornaría

Note that the "switch_density=0" case is already tested by the doctests so there's no need to add another doctest after patch 06.

comment:17 Changed 12 years ago by William Stein

Positive review (on the two extra little patches added, hence everything)!

comment:18 Changed 12 years ago by Jeroen Demeyer

Milestone: sage-4.7sage-4.7.1

comment:19 Changed 12 years ago by Jeroen Demeyer

Milestone: sage-4.7.1sage-feature

Moving the milestone since this depends on a non-positively-reviewed ticket.

comment:20 Changed 12 years ago by Jeroen Demeyer

Dependencies: #10734
Description: modified (diff)

comment:21 Changed 11 years ago by Jeroen Demeyer

Milestone: sage-featuresage-wait

Changed 11 years ago by Gonzalo Tornaría

Improve verbose logging of echelon_in_place_classical

Changed 11 years ago by Gonzalo Tornaría

Make matrix modn sparse echelon form way faster by switching to dense linbox at a good time when the matrix gets too dense

Changed 11 years ago by Gonzalo Tornaría

Attachment: trac_10733.03-documentation added

Add a test for matrix_modn_sparse._echelon_in_place_classical

Changed 11 years ago by Gonzalo Tornaría

refactor code to support p=2 case (part 2)

Changed 11 years ago by Gonzalo Tornaría

Add testing for echelonize of mod2 sparse matrices

Changed 11 years ago by Gonzalo Tornaría

Implement a fast path for strictly using a dense method (could be useful for mod2)

comment:22 Changed 11 years ago by Gonzalo Tornaría

I've reuploaded all the patches with the right format (see #10734).

Also, I moved part of the referee patch to #10734, otherwise some doctests would fail if that ticket is merged without the current one.

There's no new code at all, so I assume the positive review is still valid.


OTOH, I screwed up filenames because (a) I made a mistake uploading one patch (b) I don't have permissions to delete incorrectly uploaded patches (c) I don't have permissions to replace the referee patch.

Whoever has the rights, what needs to be done is:

  1. delete attachment trac_10733.2.03-documentation (uploaded by mistake)
  2. delete the original referee patch trac_10733.04-referee_refactor_to_support_p2.patch, and replace it by the new trac_10733.04-referee_refactor_to_support_p2.2.patch

After those two changes, the six patches are applied one after other in numerical order.

comment:23 Changed 11 years ago by Gonzalo Tornaría

apply trac_10733.01-echelon_better_verbose, trac_10733.02-echelon_make_faster, trac_10733.03-documentation, trac_10733.04-referee_refactor_to_support_p2.2.patch, trac_10733.05-doctest_mod2.patch, trac_10733.06-fastpath_for_strictly_dense_method.patch

comment:24 Changed 10 years ago by Jeroen Demeyer

Please fill in your real name as Reviewer.

comment:25 Changed 9 years ago by Volker Braun

Authors: Gonzalo TornariaGonzalo Tornaría, William Stein
Reviewers: William Stein, Gonzalo Tornaría

comment:26 Changed 8 years ago by Volker Braun

Status: positive_reviewneeds_work

I'm setting this to needs works since its not a git branch.

Note: See TracTickets for help on using tickets.