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.