Scipy sparse diags matrix design error

When using scipy.sparse.spdiags or scipy.sparse.diags, I noticed that I want me to read the error in routines, for example

scipy.sparse.spdiags([1.1,1.2,1.3],1,4,4).toarray() 

returns

 array([[ 0. , 1.2, 0. , 0. ], [ 0. , 0. , 1.3, 0. ], [ 0. , 0. , 0. , 0. ], [ 0. , 0. , 0. , 0. ]]) 

That is, for positive diagonals, it discards the first k data. One could argue that for this there is some kind of grandiose reason for programming, and I just need to fill with zeros. OK annoying, as it may be, you can use scipy.sparse.diags, which gives the correct result. However, this procedure has an error that cannot be processed.

 scipy.sparse.diags([1.1,1.2],0,(4,2)).toarray() 

gives

 array([[ 1.1, 0. ], [ 0. , 1.2], [ 0. , 0. ], [ 0. , 0. ]]) 

nice, and

 scipy.sparse.diags([1.1,1.2],-2,(4,2)).toarray() 

gives

 array([[ 0. , 0. ], [ 0. , 0. ], [ 1.1, 0. ], [ 0. , 1.2]]) 

but

 scipy.sparse.diags([1.1,1.2],-1,(4,2)).toarray() 

gives an error: ValueError: The length of the diagonal (index 0: 2 at offset -1) does not agree with the size of the matrix (4, 2). Obvious answer

 array([[ 0. , 0. ], [ 1.1, 0. ], [ 0. , 1.2], [ 0. , 0. ]]) 

and for additional random behavior we have

 scipy.sparse.diags([1.1],-1,(4,2)).toarray() 

gives

 array([[ 0. , 0. ], [ 1.1, 0. ], [ 0. , 1.1], [ 0. , 0. ]]) 

Does anyone know if there is a function to build diagonal sparse matrices that really work?

+6
source share
1 answer

Summary: spdiags works correctly even if the matrix input is not the most intuitive. diags has an error that affects some offsets in rectangular matrices. Fixed bug in scipy github.


Example for spdiags :

 >>> data = array([[1,2,3,4],[1,2,3,4],[1,2,3,4]]) >>> diags = array([0,-1,2]) >>> spdiags(data, diags, 4, 4).todense() matrix([[1, 0, 3, 0], [1, 2, 0, 4], [0, 2, 3, 0], [0, 0, 3, 4]]) 

Note that the third column of data always appears in the third column of sparse. The remaining columns also line up. But they descend where they "fall from the edge."

The input of this function is a matrix, and the input to diags is a ragged list. The diagonals of a sparse matrix have a different number of values. Therefore, the specification must conform to this in one form or another. spdiags does this by ignoring certain values, diags , accepting list input.

The error sparse.diags([1.1,1.2],-1,(4,2)) is bewildering.

spdiags equivalent works:

 In [421]: sparse.spdiags([[1.1,1.2]],-1,4,2).A Out[421]: array([[ 0. , 0. ], [ 1.1, 0. ], [ 0. , 1.2], [ 0. , 0. ]]) 

The error in this block of code is:

 for j, diagonal in enumerate(diagonals): offset = offsets[j] k = max(0, offset) length = min(m + offset, n - offset) if length <= 0: raise ValueError("Offset %d (index %d) out of bounds" % (offset, j)) try: data_arr[j, k:k+length] = diagonal except ValueError: if len(diagonal) != length and len(diagonal) != 1: raise ValueError( "Diagonal length (index %d: %d at offset %d) does not " "agree with matrix size (%d, %d)." % ( j, len(diagonal), offset, m, n)) raise 

Actual matrix constructor in diags :

 dia_matrix((data_arr, offsets), shape=(m, n)) 

This is the same constructor that uses spdiags , but without any manipulation.

 In [434]: sparse.dia_matrix(([[1.1,1.2]],-1),shape=(4,2)).A Out[434]: array([[ 0. , 0. ], [ 1.1, 0. ], [ 0. , 1.2], [ 0. , 0. ]]) 

In the dia format, the inputs are stored exactly as spdiags indicates (complete with this matrix with additional values):

 In [436]: M.data Out[436]: array([[ 1.1, 1.2]]) In [437]: M.offsets Out[437]: array([-1], dtype=int32) 

As @user2357112 points out, length = min(m + offset, n - offset is incorrect, producing 3 in the test example. Changing it to length = min(m + k, n - k) does all the cases for this (4,2) matrix But it fails with transposition: diags([1.1,1.2], 1, (2, 4))

Correction as of October 5th for this problem:

https://github.com/pv/scipy-work/commit/529cbde47121c8ed87f74fa6445c05d71353eb6c

 length = min(m + offset, n - offset, min(m,n)) 

diags([1.1,1.2], 1, (2, 4)) works with this fix.

+2
source

All Articles