# Ticket #2355(closed enhancement: fixed)

Opened 5 years ago

## [with patch; with positive review pending (?)] Write a clearer submatrix implementation

Reported by: Owned by: dfdeshom dfdeshom minor sage-2.11 linear algebra

### Description

The current implementation of the submatrix command could be easier to use given that matrix_from_rows_and_columns is nicely suited for this task; ie, this should just work

```A.submatrix([1..2],[0..1])
```

if A is (at least) a 3x2 matrix

## Change History

### comment:1 follow-up: ↓ 2 Changed 5 years ago by was

I vote for at least seriously considering using the same notation as numpy for submatrices, whatever that is, by extending getitem...

### comment:2 in reply to: ↑ 1 Changed 5 years ago by dfdeshom

I vote for at least seriously considering using the same notation as numpy for submatrices, whatever that is, by extending getitem...

Agreed. Looking at the relevant section of the code, it looks like it should not be too hard

### comment:3 Changed 5 years ago by dfdeshom

Here's a patch that makes the following possible:

```            sage: m=[(1, -2, -1, -1,9), (1, 8, 6, 2,2), (1, 1, -1, 1,4), (-1, 2, -2, -1,4)];M= matrix(m)
sage: M
[ 1 -2 -1 -1  9]
[ 1  8  6  2  2]
[ 1  1 -1  1  4]
[-1  2 -2 -1  4]

Get The 2 x 2 submatrix of M, starting at row index and column
index 1
sage: M[1:3,1:3]
[ 8  6]
[ 1 -1]

Get the 2 x 3 submatrix of M starting at row index and column
index 1:
sage: M[1:3,[1..3]]
[ 8  6  2]
[ 1 -1  1]

Get the second column of M:
sage: M[1:,0]
[ 1]
[ 1]
[-1]

sage: M[range(2),:]
[ 1 -2 -1 -1  9]
[ 1  8  6  2  2]
sage: M[range(2),4]
[9]
[2]
sage: M[range(3),range(5)]
[ 1 -2 -1 -1  9]
[ 1  8  6  2  2]
[ 1  1 -1  1  4]

sage: M[3,range(5)]
[-1  2 -2 -1  4]
sage: M[3,:]
[-1  2 -2 -1  4]
sage: M[3,4]
4

```

### comment:4 follow-up: ↓ 5 Changed 5 years ago by was

I have to register the obvious concern. You've replaced what I wrote:

```554	 	            if PyTuple_Size(key) != 2:
555	 	                raise IndexError, "index must be an integer or pair of integers"
556	 	            i = <object> PyTuple_GET_ITEM(key, 0)
557	 	            j = <object> PyTuple_GET_ITEM(key, 1)
```

by

``` 	587	            s1 = key[0]
588	            s2 = key[1]
589
590	            if isinstance(s1,sage.rings.integer.Integer) and \
591	                   isinstance(s2,sage.rings.integer.Integer):
592	                self.check_bounds(s1, s2)
593	                return self.get_unsafe(s1, s2)
```

If A is a matrix doing A[i,j] is a sort of "critical speed operation", so I'm concerned about speed. Thoughts? Benchmarks?

### comment:5 in reply to: ↑ 4 Changed 5 years ago by dfdeshom

The new method is slower, but not by much (surprisingly for me):

```# sage-main
sage: %timeit M[3,4]
10000 loops, best of 3: 69.5 µs per loop

#sage-getitem
sage: %timeit M[3,4]
10000 loops, best of 3: 71.2 µs per loop

#sage-main
sage: %timeit M[3:]
10000 loops, best of 3: 68.3 µs per loop

#sage-getitem
sage: %timeit M[3:]
10000 loops, best of 3: 69.5 µs per loop
```
```#sage-getitem only
sage: %timeit M[:4,range(4)]
10000 loops, best of 3: 87.8 µs per loop

sage: %timeit M[3:,:4]
10000 loops, best of 3: 78.8 µs per loop

sage: %timeit M[3:,0]
10000 loops, best of 3: 74.1 µs per loop

```

matrix_from_rows_and_columns could also be made a little faster. I plan to post another ticket/patch for it.

### comment:6 Changed 5 years ago by dfdeshom

• Status changed from new to assigned

### comment:7 Changed 5 years ago by was

• Summary changed from Write a clearer submatrix implementation to [with patch; needs review] Write a clearer submatrix implementation

### comment:8 Changed 5 years ago by dfdeshom

• Type changed from defect to enhancement

### comment:9 Changed 5 years ago by dfdeshom

• Milestone changed from sage-2.11 to sage-2.10.4

### comment:10 Changed 5 years ago by jsp

Can't give this patch a positive review. After I applied it to a brand new sage-2.10.3 a lot of code seems to be broken.

I do not quite understand what is going on.

Jaap

### comment:11 Changed 5 years ago by dfdeshom

• Summary changed from [with patch; needs review] Write a clearer submatrix implementation to [with patch; negative review] Write a clearer submatrix implementation

My changes seem to have broken matrix.__reduce___. I'll look into it.

### comment:12 Changed 5 years ago by dfdeshom

• Summary changed from [with patch; negative review] Write a clearer submatrix implementation to [with patch; needs review] Write a clearer submatrix implementation

Found the bug. All doctests for matrix0,1,2 pass now.

### comment:13 Changed 5 years ago by jsp

• Summary changed from [with patch; needs review] Write a clearer submatrix implementation to [with patch; with negative review] Write a clearer submatrix implementation

Applyin this bundle I still got an error:

```----------------------------------------------------------------------
The following tests failed:

sage -t  devel/sage-main/sage/matrix/matrix_integer_dense.pyx
Total time for all tests: 113.4 seconds

```
```sage -t  devel/sage-main/sage/matrix/matrix_integer_dense.pyx**********************************************************************
File "matrix_integer_dense.pyx", line 2961:
Exception raised:
Traceback (most recent call last):
File "/home/jaap/work/sage-2.10.3/local/lib/python2.5/doctest.py", line 1212, in __run
compileflags, 1) in test.globs
File "<doctest __main__.example_52[1]>", line 1, in <module>
File "matrix_integer_dense.pyx", line 3013, in sage.matrix.matrix_integer_dense.Matrix_integer_dense._add_row_and_maintain_echelon_form
if b % a == 0:
File "matrix0.pyx", line 1927, in sage.matrix.matrix0.Matrix.__mod__
v[i] = v[i] % p
File "integer.pyx", line 1335, in sage.rings.integer.Integer.__mod__
File "integer.pyx", line 3747, in sage.rings.integer.integer
File "integer.pyx", line 356, in sage.rings.integer.Integer.__init__
TypeError: unable to coerce element to an integer
**********************************************************************
File "matrix_integer_dense.pyx", line 2970:
Exception raised:
Traceback (most recent call last):
File "/home/jaap/work/sage-2.10.3/local/lib/python2.5/doctest.py", line 1212, in __run
compileflags, 1) in test.globs
File "<doctest __main__.example_52[4]>", line 1, in <module>
File "matrix_integer_dense.pyx", line 3013, in sage.matrix.matrix_integer_dense.Matrix_integer_dense._add_row_and_maintain_echelon_form
if b % a == 0:
File "matrix0.pyx", line 1927, in sage.matrix.matrix0.Matrix.__mod__
v[i] = v[i] % p
File "integer.pyx", line 1335, in sage.rings.integer.Integer.__mod__
File "integer.pyx", line 3747, in sage.rings.integer.integer
File "integer.pyx", line 356, in sage.rings.integer.Integer.__init__
TypeError: unable to coerce element to an integer
**********************************************************************
2 of   5 in __main__.example_52
***Test Failed*** 2 failures.
For whitespace errors, see the file .doctest_matrix_integer_dense.pyx
[3.2 s]

```

### comment:14 Changed 5 years ago by dfdeshom

• Summary changed from [with patch; with negative review] Write a clearer submatrix implementation to [with patch; needs review] Write a clearer submatrix implementation

Thanks I've updated the bundle to take care of these doctest failures. Give it one more try if you can.

### comment:15 Changed 5 years ago by jason

What should I apply to review this? The patch, the bundle, or both?

Just the bundle

### comment:17 Changed 5 years ago by jsp

• Summary changed from [with patch; needs review] Write a clearer submatrix implementation to [with bundle; negative review] Write a clearer submatrix implementation

Here we are again:

```sage -t  devel/sage-main/sage/matrix/matrix0.pyx            **********************************************************************
File "matrix0.pyx", line 553:
sage: M[1:,0]
Expected:
[ 1]
[ 1]
[-1]
Got:
1
**********************************************************************
File "matrix0.pyx", line 566:
sage: M[range(2),4]
Expected:
[9]
[2]
Got:
9
**********************************************************************
2 of  23 in __main__.example_21
***Test Failed*** 2 failures.
For whitespace errors, see the file .doctest_matrix0.pyx
[2.3 s]
sage -t  devel/sage-main/sage/matrix/matrix1.pyx
[3.9 s]
sage -t  devel/sage-main/sage/matrix/action.pyx
[1.8 s]
sage -t  devel/sage-main/sage/matrix/matrix_generic_sparse.pyx
[1.6 s]

----------------------------------------------------------------------
The following tests failed:

sage -t  devel/sage-main/sage/matrix/matrix0.pyx
Total time for all tests: 121.8 seconds

```

### comment:18 follow-ups: ↓ 19 ↓ 20 Changed 5 years ago by dfdeshom

• Summary changed from [with bundle; negative review] Write a clearer submatrix implementation to [with bundle; needs review] Write a clearer submatrix implementation

A little embarrassing, isn't it? The bundle 2355.hg has been updated and should pass all doctests in sage/matrix

### comment:19 in reply to: ↑ 18 ; follow-up: ↓ 21 Changed 5 years ago by jsp

• Summary changed from [with bundle; needs review] Write a clearer submatrix implementation to [with bundle; with negative review] Write a clearer submatrix implementation

A little embarrassing, isn't it? The bundle 2355.hg has been updated and should pass all doctests in sage/matrix

This patch once again broke a lot of code!

I tried it on a brand new sage-2.10.4.rc0

Just a snapshot:

```      File "/home/jaap/work/sage-2.10.3/local/lib/python2.5/site-packages/sage/algebras/free_algebra.py", line 229, in _repr_
self.__ngens, self.gens(), self.base_ring())
File "sage_object.pyx", line 92, in sage.structure.sage_object.SageObject.__repr__
File "/home/jaap/work/sage-2.10.3/local/lib/python2.5/site-packages/sage/algebras/free_algebra_element.py", line 74, in _repr_
x = repr_lincomb(mons, cffs).replace("*1 "," ")
File "/home/jaap/work/sage-2.10.3/local/lib/python2.5/site-packages/sage/misc/misc.py", line 508, in repr_lincomb
b = str(symbols[i])
File "sage_object.pyx", line 92, in sage.structure.sage_object.SageObject.__repr__
File "/home/jaap/work/sage-2.10.3/local/lib/python2.5/site-packages/sage/monoids/free_monoid_element.py", line 103, in _repr_
x = self.parent().variable_names()
File "parent_gens.pyx", line 375, in sage.structure.parent_gens.ParentWithGens.variable_names
RuntimeError: maximum recursion depth exceeded in cmp
**********************************************************************

```

Testing seems to hang on

```sage -t  devel/sage-main/sage/matrix/matrix_symbolic_dense.pyx
[1]+  Stopped                 ./sage -t devel/sage-main/sage/matrix/
[jaap@paix sage-2.10.4.rc0]\$

```

### comment:20 in reply to: ↑ 18 Changed 5 years ago by mabshoff

The bundle 2355.hg has been updated and should pass all doctests in sage/matrix

Didier,

please post only Mercurial patches. hg export tip > foo.patch works perfectly well. Bundles for single commits are a serious annoyance.

Cheers,

Michael

### comment:21 in reply to: ↑ 19 Changed 5 years ago by dfdeshom

• Summary changed from [with bundle; with negative review] Write a clearer submatrix implementation to [with bundle; not ready for review] Write a clearer submatrix implementation
• Milestone changed from sage-2.10.4 to sage-2.11

I tried it on a brand new sage-2.10.4.rc0

All tests pass on 2.10.3 here. 2.10.4rc0 obviously broke something so I'll wait until 2.10.4 is officially released to see what broke what.

### comment:22 Changed 5 years ago by dfdeshom

• Summary changed from [with bundle; not ready for review] Write a clearer submatrix implementation to [with patch; needs review] Write a clearer submatrix implementation
• Milestone changed from sage-3.0 to sage-2.11

Jaap, 2.4.10 with this patch passes all tests I care about (I ran sage -testall and only the matplotlib tests fail). This is on sage.math. I've reposted the patch file, it is named 2355.patch. Please try it out and if you still have failures, let me know which files/tests break what, your platform, etc.

### comment:23 Changed 5 years ago by jsp

• Summary changed from [with patch; needs review] Write a clearer submatrix implementation to [with patch; with positive review] Write a clearer submatrix implementation

The patch works for me now!

Cheers,

Jaap

### comment:24 Changed 5 years ago by mabshoff

• Summary changed from [with patch; with positive review] Write a clearer submatrix implementation to [with patch; with positive review pending performance question] Write a clearer submatrix implementation

Hi,

while the patch passes doctests for me with 2.10.4 [modulo some known problems] I am still reluctant to apply this until the performance concern that was raised has been addressed.

Cheers,

Michael

### comment:25 Changed 5 years ago by dfdeshom

Hi,

while the patch passes doctests for me with 2.10.4 [modulo some known problems] I am still reluctant to apply this until the performance concern that was raised has been addressed.

Cheers,

Michael

I've lost some speed while working on this. I've optimized the code a bit. Overall, my version is still slower because I'm doing more type-checking and handling more functionality (I think that's the reason):

```# M is the same as above

# element-by-element: slower by about 0.1 micro-seconds ~ 1.2x slower

#2355
sage: %timeit M[2,3]
1000000 loops, best of 3: 395 ns per loop

#main
sage: %timeit M[2,3]
1000000 loops, best of 3: 304 ns per loop

# single slices: slower by < 1 micro-seconds ~ 1.01x slower

#2355
sage: %timeit M[:4]
10000 loops, best of 3: 48 Âµs per loop

#main
sage: %timeit M[:4]
10000 loops, best of 3: 48.5 Âµs per loop

# Getting a whole row: faster by < 1 micro-seconds
#2355
sage: %timeit M[3]
10000 loops, best of 3: 72.2 µs per loop

#main
sage: %timeit M[3]
10000 loops, best of 3: 72.9 µs per loop
```

I'm not posting speed comparisons using the other cases (ie, M[:3,2:]) since sage doesn't handle them right now.

If people think this speed loss is unacceptable, I could try to fold the extra functionality into the submatrix method instead which would leave __getitem__ idem.

2355.patch is updated. Warning: I've added a new file: ext/python_slice.pxi to have fast access to slice objects so the patch is a bit heavy and rebuilds everything ( nearly everything depends on python.pxi). Passes all tests in matrix/.

### comment:26 Changed 5 years ago by gfurnish

• Summary changed from [with patch; with positive review pending performance question] Write a clearer submatrix implementation to [with patch; needs new review for performance] Write a clearer submatrix implementation
```# M is the same as above

# element-by-element:

#2355
sage: %timeit M[2,3]
1000000 loops, best of 3: 111 ns per loop

#main
sage: %timeit M[2,3]
1000000 loops, best of 3: 111 ns per loop

# single slices:

#2355
sage: %timeit M[:4]
10000 loops, best of 3: 28.4 Âµs per loop

#main
sage: %timeit M[:4]
10000 loops, best of 3: 29.6 Âµs per loop

# Getting a whole row:
#2355
sage: %timeit M[3]
10000 loops, best of 3: 42.3 Âµs per loop

#main
sage: %timeit M[3]
10000 loops, best of 3: 43.7 Âµs per loop
```

I note these were the best times. In my tests for (1) the new source averaged faster then main.

### Changed 5 years ago by gfurnish

Apply on top of 2355.patch

### comment:27 Changed 5 years ago by gfurnish

I ran this on sage math, and while 1 was on average faster (although best time was slower), both 2 and 3 were slower. I'm chalking this up to differences between the Core 2 and Opteron architecture. getitem is not a cdef function, so I don't see why we should be worrying over a 1Âµs slowdown on some older systems. Also, 1 is the common case, so I don't see a good reason to worry about the slower cases of 2 and 3. Note this test was done with "-O2" and the new build system.

### comment:28 Changed 5 years ago by dfdeshom

• Summary changed from [with patch; needs new review for performance] Write a clearer submatrix implementation to [with patch; with positive review pending (?)] Write a clearer submatrix implementation

### comment:29 Changed 5 years ago by mabshoff

Looks good to me. My concerns regarding performance have been addressed. It merges cleanly against my 2.11.alpha2 tree. If this doctests ok I will merge.

Cheers,

Michael

### comment:30 Changed 5 years ago by mabshoff

• Status changed from assigned to closed
• Resolution set to fixed

Merged in Sage 2.11.alpha2

Note: See TracTickets for help on using tickets.