Use echelonize instead of echelon_form in a few places

[tscrim]

We do not need to create another copy of the augmented matrix to keep memory usage down (including putting it in a cache) and because it is faster.

We do this change in __invert__, right_kernel_matrix, and _solve_right_nonsingular_square.

We also can take advantage of the sparse structure and our knowledge of the augmented matrix [I|A^-1] when reconstructing the matrix.

Component: performance

Author: Travis Scrimshaw

Branch/Commit: cfbb776

Reviewer: Vincent Delecroix

Issue created by migration from https://trac.sagemath.org/ticket/24122

[tscrim]

Branch: public/linear_algebra/echelonize_in_invert-24122

[tscrim]

Author: Travis Scrimshaw

[tscrim]

comment:1

Bad branch not based off develop. Fixed.

[tscrim]

Commit: d37d983

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

3990da4

Do not create another matrix in __invert__.

[sagetrac-git]

Changed commit from d37d983 to 3990da4

[sagetrac-git]

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

294150f

Do different things for sparse and dense matrices from the augmented matrix.

[sagetrac-git]

Changed commit from 3990da4 to 294150f

[tscrim]

comment:4

This last change nets me ~8% speedup:

sage: %timeit ~(matrix({(200-i,i): 1/i for i in range(1,200,1)}, sparse=True) + 1)
10 loops, best of 3: 159 ms per loop
sage: %timeit ~(matrix({(100-i,i): 1/i for i in range(1,100,1)}, sparse=True) + 1)
10 loops, best of 3: 28.6 ms per loop
sage: %timeit ~(matrix({(20-i,i): 1/i for i in range(1,20,1)}, sparse=True) + 1)
The slowest run took 24.10 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 1.27 ms per loop

versus

sage: %timeit ~(matrix({(200-i,i): 1/i for i in range(1,200,1)}, sparse=True) + 1)
10 loops, best of 3: 167 ms per loop
sage: %timeit ~(matrix({(100-i,i): 1/i for i in range(1,100,1)}, sparse=True) + 1)
10 loops, best of 3: 30.7 ms per loop
sage: %timeit ~(matrix({(20-i,i): 1/i for i in range(1,20,1)}, sparse=True) + 1)
1000 loops, best of 3: 1.36 ms per loop

[tscrim]

comment:5

Also versus develop, which this test is not so good for showing how the first change would improve things as the matrices are relatively tiny:

sage: %timeit ~(matrix({(200-i,i): 1/i for i in range(1,200,1)}, sparse=True) + 1)
10 loops, best of 3: 171 ms per loop
sage: %timeit ~(matrix({(100-i,i): 1/i for i in range(1,100,1)}, sparse=True) + 1)
10 loops, best of 3: 30.2 ms per loop
sage: %timeit ~(matrix({(20-i,i): 1/i for i in range(1,20,1)}, sparse=True) + 1)
The slowest run took 11.65 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 1.3 ms per loop

[sagetrac-git]

Changed commit from 294150f to d6e5a09

[sagetrac-git]

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

d6e5a09

Use echelonize in two other places.

[sagetrac-git]

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

cfbb776

Merge branch 'develop' into public/linear_algebra/echelonize_in_invert-24122

[sagetrac-git]

Changed commit from d6e5a09 to cfbb776

[videlec]

comment:9

This is a bit weird

+        for i in range(nrows):
+            del data[i,i]
+        data = {(r,c-nrows): data[r,c] for (r,c) in data}

You are modifying data and then override it. If you can not do the modification inplace, the following looks simpler

data = {(r,c-nrows): data[r,c] for (r,c) in data if c >= nrows}

[tscrim]

comment:10

Yes, it looks more simple, but it actually takes longer because it has to do the if check on every nonzero object. Whereas the version I have just removes objects from the dict. I suspect this is more of a micro improvement, but that could be a lot of extra checks for a really big sparse matrix (at least, I have a 62001 x 62001 sparse matrix with 1756951 non-zero entries).

Ideally I would like to modify the data in place and just update the indices, but hash tables are not good for doing that. One of these days, someone (me) should implement another data structure for sparse matrices...

[videlec]

Reviewer: Vincent Delecroix

[videlec]

comment:11

Oh, I see.

Note that for matrices over ZZ or QQ, there is a custom C datatype which is an array of sparse vectors. This is not ideal but already faster than dictionaries in most contexts.

[tscrim]

comment:13

Thanks for the review.

[vbraun]

Changed branch from public/linear_algebra/echelonize_in_invert-24122 to cfbb776