fast slicing of matrix

[davidp]

Is there a fast way to create a submatrix?

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


sage: version()
'Sage Version 4.0.alpha0, Release Date: 2009-05-15'
sage: G = graphs.GridGraph([100,100])
sage: L = G.laplacian_matrix()
sage: L
10000 x 10000 sparse matrix over Integer Ring
sage: time M = L[1:9999,1:9999]
CPU times: user 24.93 s, sys: 0.04 s, total: 24.97 s
Wall time: 25.27 s


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

I am just interested in deleting a single row and column of the matrix
(not necessarily the first).

Thanks,
Dave

[William Stein]

There is no fast way to do that right now. One could easily add code
to SAGE_ROOT/devel/sage/sage/matrix/matrix_integer_sparse.pyx that
would provide blazingly fast deletion of a row, and reasonably fast
deletion of a column. Of course it would be better to implement
arbitrary slicing in some optimized way in matrix_integer_sparse.pyx.
I hope somebody does so.

William

[Jason Grout]

Just looking at the generic code, it seems that it goes through each and
every index position in the slice, setting the new matrix entry to the
old one. This is obviously the wrong thing to do for sparse matrices,
and can likely trivially be made faster. I think all you may have to do
is override the matrix_from_rows_and_columns for sparse matrices.

Jason

[davidp]

I tried adding a delete_row method to matrix_integer_sparse.pyx but
stopped after getting the error message:

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Error converting Pyrex file to C:
------------------------------------------------------------
...
add_mpz_vector_init(&M._matrix[i], &self._matrix[i], &
(<Matrix_integer_sparse>right)._matrix[i], mul)
mpz_clear(mul)
return M

# added by David Perkinson
cpdef ModuleElement _delete_row_(self, Py_ssize_t row_number):
^
------------------------------------------------------------

/usr/local/sage-devel/devel/sage-devel/sage/matrix/
matrix_integer_sparse.pyx:216:10: C method '_delete_row_' not
previously declared in definition part of extension type
Error running command, failed with status 256.
sage: There was an error installing modified sage library code.

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

I might have to actually read the Cython documentation!

Dave

[William Stein]

On Tue, Jun 2, 2009 at 7:22 PM, davidp <dav...@reed.edu> wrote:
>
> I tried adding a delete_row method to matrix_integer_sparse.pyx but
> stopped after getting the error message:

Just def your method -- do not cpdef it. You can still use Cython
code in a cpdef's method and it will be just as fast.

Make sure that your "delete row" method returns a *new* matrix by the way...

William

--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

[Jason Grout]

William Stein wrote:
> On Tue, Jun 2, 2009 at 7:22 PM, davidp <dav...@reed.edu> wrote:
>> I tried adding a delete_row method to matrix_integer_sparse.pyx but
>> stopped after getting the error message:
>
> Just def your method -- do not cpdef it. You can still use Cython
> code in a cpdef's method and it will be just as fast.
>
> Make sure that your "delete row" method returns a *new* matrix by the way...

You might fold your functionality into matrix_from_rows_and_columns.
That way slicing will automatically use it.

Maybe you could intelligently guess from the rows and columns requested
whether it is faster to just copy the matrix, then delete some
rows/columns, or to build the new matrix up. For example, if I'm asking
for every row and column in the matrix except one, almost surely it is
faster to copy the dictionary and delete the appropriate row/column keys.

Then again, I can see the advantage of just having a delete_rows or
delete_columns method; it's definitely clear what you are doing.

Jason

[davidp]

I can't figure out how to create a new sparse matrix (see below).

Sorry,
Dave

M = Matrix_integer_sparse.__new__(Matrix_integer_sparse,
sage.matrix.matrix_space.MatrixSpace(ZZ, self._nrows-1, self._ncols,
sparse=True), None, None, None)


-------------------
...
def _delete_row_(self, Py_ssize_t row_number):
cdef Py_ssize_t i
cdef mpz_vector* self_row, *M_row
cdef Matrix_integer_sparse M

M = Matrix_integer_sparse.__new__(Matrix_integer_sparse,
sage.matrix.matrix_space.MatrixSpace(ZZ, self._nrows-1, self._ncols,
sparse=True), None, None, None)

^
------------------------------------------------------------

/usr/local/sage-devel/devel/sage-devel/sage/matrix/
matrix_integer_sparse.pyx:221:101: Compiler crash in
AnalyseExpressionsTransform

ModuleNode.body = StatListNode(matrix_integer_sparse.pyx:23:0)
StatListNode.stats[13] = StatListNode(matrix_integer_sparse.pyx:43:5)
StatListNode.stats[0] = CClassDefNode(matrix_integer_sparse.pyx:43:5,
as_name = u'Matrix_integer_sparse',
base_class_module = u'matrix_sparse',
base_class_name = u'Matrix_sparse',
class_name = u'Matrix_integer_sparse',
module_name = '',
visibility = 'private')
CClassDefNode.body = StatListNode(matrix_integer_sparse.pyx:55:4)
StatListNode.stats[9] = DefNode(matrix_integer_sparse.pyx:216:4,
name = u'_delete_row_',
num_required_args = 2,
reqd_kw_flags_cname = '0')
File 'Nodes.py', line 360, in analyse_expressions: StatListNode
(matrix_integer_sparse.pyx:217:8)
File 'Nodes.py', line 2835, in analyse_expressions:
SingleAssignmentNode(matrix_integer_sparse.pyx:221:41)
File 'Nodes.py', line 2928, in analyse_types: SingleAssignmentNode
(matrix_integer_sparse.pyx:221:41)
File 'ExprNodes.py', line 2281, in analyse_types: SimpleCallNode
(matrix_integer_sparse.pyx:221:41)
File 'ExprNodes.py', line 3105, in analyse_types: TupleNode
(matrix_integer_sparse.pyx:221:41,
is_sequence_constructor = 1)
File 'ExprNodes.py', line 2972, in analyse_types: TupleNode
(matrix_integer_sparse.pyx:221:41,
is_sequence_constructor = 1)
File 'ExprNodes.py', line 2522, in analyse_types: GeneralCallNode
(matrix_integer_sparse.pyx:221:101)

Compiler crash traceback from this point on:
File "/usr/local/sage-devel/local/lib/python2.5/site-packages/Cython/
Compiler/ExprNodes.py", line 2532, in analyse_types
if hasattr(self.function, 'entry') and not
self.function.entry.as_variable:
AttributeError: 'NoneType' object has no attribute 'as_variable'

Error running command, failed with status 256.
sage: There was an error installing modified sage library code.

[William Stein]

On Tue, Jun 2, 2009 at 8:45 PM, davidp <dav...@reed.edu> wrote:
>
> I can't figure out how to create a new sparse matrix (see below).
>
> Sorry,
> Dave

You can make a new matrix with the analogous parent but a different
shape using the new_matrix method.

sage: a = random_matrix(ZZ,3,4,sparse=True)
sage: a.new_matrix(2,4)
[0 0 0 0]
[0 0 0 0]

In your cython code you would write something like

cdef Matrix_integer_sparse M = a.new_matrix(2,4)

William

[davidp]

I might be going down the wrong path, but is seems like I would want
to add a

copy_mpz_vector_init

function to vector_integer_sparse_c.pxi that would make a copy of an
mpz_vector. Is that overkill?

My first idea was to use add_mpz_vector_init(sum, v, w, mul), which
sets sum equal to v + mul*w, but let mul = 0. However
add_mpz_vector_init returns the zero vector when mul = 0 !

Dave

[William Stein]

On Tue, Jun 2, 2009 at 10:10 PM, davidp <dav...@reed.edu> wrote:
>
> I might be going down the wrong path, but is seems like I would want
> to add a
>
>         copy_mpz_vector_init
>
> function to  vector_integer_sparse_c.pxi that would make a copy of an
> mpz_vector.  Is that overkill?

That sounds like a good idea to me.

William

[davidp]

Success!

sage: G = graphs.GridGraph([100,100])
sage: L = G.laplacian_matrix()
sage: L
10000 x 10000 sparse matrix over Integer Ring

sage: sage: time M = L[1:9999]
CPU times: user 21.98 s, sys: 0.01 s, total: 21.99 s
Wall time: 22.93 s
sage: sage: time N = L._delete_row_(0)
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01 s
sage: sage: timeit('N = L._delete_row_(0)')
25 loops, best of 3: 10.9 ms per loop

I am travelling for the next several days but will follow up on this
soon.

Thanks a lot!