[SciPy-User] preserving sparse matrix type under operations

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[SciPy-User] preserving sparse matrix type under operations

Gideon Simpson-3
I was wondering if there’s a way to (cheaply) retain sparsity type of matrices which are of the same sparsity pattern?  For instance, if I have two DIA formatted matrices with the same banding, and I want to compute

A + s * B, where s is a scalar, it seems to convert everything over to CSR.  Is there a way to avoid this?

-gideon


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: preserving sparse matrix type under operations

Juan Nunez-Iglesias
Hey Gideon,

I believe this behaviour exists to prevent having to implement different adding operations for all possible combinations of sparse matrices. You can see in the source code that it’s doing exactly what you suspect:

In [24]: type(A)
Out[24]: scipy.sparse.dia.dia_matrix

In [25]: %psource A.__add__
    def __add__(self, other):   # self + other
        return self.tocsr().__add__(other)

If the banding is indeed the same, you can just add the `.data` attributes together:

In [16]: A = sparse.spdiags([[0, 1, 2], [0, 1, 2.]], [0, 1], 3, 3)
In [17]: B = sparse.spdiags([[9, 2, 1], [5, 5, 5.]], [0, 1], 3, 3)
In [18]: s = 0.5

In [19]: C = A.copy()
In [20]: C.data += s * B.data

In [21]: C.A
Out[21]:
array([[ 4.5,  3.5,  0. ],
       [ 0. ,  2. ,  4.5],
       [ 0. ,  0. ,  2.5]])

If you are going to be doing a lot of this and want to preserve nice linear algebra syntax, you can monkey-patch as follows

In [28]: def add_diag(self, other):
    ...:     if (isinstance(other, sparse.dia.dia_matrix) and
    ...:             np.all(self.offsets == other.offsets)):
    ...:         result = self.copy()
    ...:         result.data += other.data
    ...:         return result
    ...:     else:
    ...:         return self.tocsr().__add__(other)
    ...:

In [29]: sparse.dia.dia_matrix.__add__ = add_diag

In [30]: type(A + s * B)
Out[30]: scipy.sparse.dia.dia_matrix

Hope that helps!

Juan.

On 3 Apr 2017, 1:05 PM -0400, Gideon Simpson <[hidden email]>, wrote:
I was wondering if there’s a way to (cheaply) retain sparsity type of matrices which are of the same sparsity pattern?  For instance, if I have two DIA formatted matrices with the same banding, and I want to compute

A + s * B, where s is a scalar, it seems to convert everything over to CSR.  Is there a way to avoid this?

-gideon

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Loading...