Unexpected error from sparse.spdiags()
Question: In Python 3 I am trying to run the following line of code to get a particular
sparse matrix.

`sparse.spdiags(np.concatenate((-np.ones((9,1)), np.ones((9,1))), axis=1), [0,
1], 9, 10)`

This gives the following error message:

    
    
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "/usr/lib/python3/dist-packages/scipy/sparse/construct.py", line 61, in spdiags
        return dia_matrix((data, diags), shape=(m,n)).asformat(format)
      File "/usr/lib/python3/dist-packages/scipy/sparse/dia.py", line 138, in __init__
        % (self.data.shape[0], len(self.offsets)))
    ValueError: number of diagonals (9) does not match the number of offsets (2)
    

Running what I understand to be the equivalent code in Octave seems to get me
a sparse matrix.

    
    
    spdiags([-ones(9,1) ones(9,1)],[0 1],9,10)
    
    Compressed Column Sparse (rows = 9, cols = 10, nnz = 18 [20%])
    
    (1, 1) -> -1
    (1, 2) ->  1
    (2, 2) -> -1
    (2, 3) ->  1
    (3, 3) -> -1
    (3, 4) ->  1
    (4, 4) -> -1
    (4, 5) ->  1
    (5, 5) -> -1
    (5, 6) ->  1
    (6, 6) -> -1
    (6, 7) ->  1
    (7, 7) -> -1
    (7, 8) ->  1
    (8, 8) -> -1
    (8, 9) ->  1
    (9, 9) -> -1
    (9, 10) ->  1
    

Any ideas on why they are behaving differently, and how to fix it?

**ADDITION**

I'm having an additional problem with Scipy.sparse's output vs Octave's.

PYTHON

    
    
    >>> sparse.spdiags(np.concatenate((-np.ones((9,1)),np.ones((9,1))), axis=1).T, [0,1],9,10).A
    array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])
    
    >>> sparse.spdiags(np.concatenate((-np.ones((9,1)),np.ones((9,1))), axis=1).T, [0,1],9,10).A.shape                                                      
    (9, 10)
    
    >>> sparse.spdiags(np.vstack([-np.ones(9),np.ones(9)]), [0,1],9,10).A                         
    array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])
    
    >>> sparse.spdiags(np.vstack([-np.ones(9),np.ones(9)]), [0,1],9,10).A.shape
    (9, 10) 
    
    >>> sparse.spdiags(np.ones(9)*[[-1],[1]], [0,1],9,10).A
    array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])
    
    >>> sparse.spdiags(np.ones(9)*[[-1],[1]], [0,1],9,10).A.shape
    (9, 10) 
    

OCTAVE

    
    
    >full(spdiags([-ones(9,1) ones(9,1)],[0 1],9,10))
    ans =
    
      -1   1   0   0   0   0   0   0   0   0
       0  -1   1   0   0   0   0   0   0   0
       0   0  -1   1   0   0   0   0   0   0
       0   0   0  -1   1   0   0   0   0   0
       0   0   0   0  -1   1   0   0   0   0
       0   0   0   0   0  -1   1   0   0   0
       0   0   0   0   0   0  -1   1   0   0
       0   0   0   0   0   0   0  -1   1   0
       0   0   0   0   0   0   0   0  -1   1
    
    >size(full(spdiags([-ones(9,1) ones(9,1)],[0 1],9,10)))
    ans =
    
    9   10
    

Why does scipy and Octave not give the same value at the last column of the
last row?


Answer: Your `concatenate` produces a (9,2) matrix:

    
    
    In [310]: np.concatenate((-np.ones((9,1)), np.ones((9,1))), axis=1)
    Out[310]: 
    array([[-1.,  1.],
           [-1.,  1.],
           [-1.,  1.],
           [-1.,  1.],
           [-1.,  1.],
           [-1.,  1.],
           [-1.,  1.],
           [-1.,  1.],
           [-1.,  1.]])
    
    In [311]: _.shape
    Out[311]: (9, 2)
    

`spdiags` doc describes this `data` parameter as `matrix diagonals stored row-
wise`. That is, each row of the matrix corresponds to a diagonal. 9 rows, but
only 2 values in `[0,1]`.

This is an important difference, that I alluded to in my previous answer,
though maybe I didn't stress it enough.

If you want 2 diagonals, you need give it a `(2,9)` array, such as the
transpose of this matrix:

    
    
    In [317]: sparse.spdiags(np.concatenate((-np.ones((9,1)),
        np.ones((9,1))), axis=1).T, [0,1],9,10)
    Out[317]: 
    <9x10 sparse matrix of type '<class 'numpy.float64'>'
        with 18 stored elements (2 diagonals) in DIAgonal format>
    

You could also construct the diagonals with:

    
    
    In [321]: np.concatenate([-np.ones((1,9)), np.ones((1,9))],axis=0)
    Out[321]: 
    array([[-1., -1., -1., -1., -1., -1., -1., -1., -1.],
           [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]])
    

Or `np.vstack([-np.ones(9),np.ones(9)])` or `np.ones(9)*[[-1],[1]]`.

* * *

Look at my previous answer, but change the final shape, to make more columns
than rows):

    
    
    octave:17> reshape (1:12, 4, 3)
    ans =
        1    5    9
        2    6   10
        3    7   11
        4    8   12
    
    octave:18> full(spdiags (reshape (1:12, 4, 3), [-1 0 1], 4,5))
    ans =
        5    9    0    0    0
        1    6   10    0    0
        0    2    7   11    0
        0    0    3    8   12
    
    In [327]: np.arange(1,13).reshape(3,4)
    Out[327]: 
    array([[ 1,  2,  3,  4],
           [ 5,  6,  7,  8],
           [ 9, 10, 11, 12]])
    
    In [328]: sparse.spdiags(np.arange(1,13).reshape(3,4), [-1, 0, 1], 4,5).A
    Out[328]: 
    array([[ 5, 10,  0,  0,  0],
           [ 1,  6, 11,  0,  0],
           [ 0,  2,  7, 12,  0],
           [ 0,  0,  3,  8,  0]])
    

In Octave (and presumably MATLAB) the +1 diagonal starts with 9, ends with 12,
the full last column from the input matrix. Look at [2,6,10] - in a right
angle arrangment.

The scipy diagonal starts with 10, ends with an added 0. The `9` is invisible
in the nonexistent row above. Look at [2,6,10] - in one column.

They are both consistent - in their own way. So at least when there are more
columns than rows, you'll need to take the difference into account when
creating the input matrix.

Another `scipy` function takes the ambiguity out, by expecting the correct
number of elements for each diagonal (as a list of lists):

    
    
    In [337]: sparse.diags([[1,2,3],[5,6,7,8],[9,10,11,12]],[-1,0,1],(4,5),dtype=int).A
    Out[337]: 
    array([[ 5,  9,  0,  0,  0],
           [ 1,  6, 10,  0,  0],
           [ 0,  2,  7, 11,  0],
           [ 0,  0,  3,  8, 12]])
    

Note that I had to omit the `4`.

Look at the `tocoo` method in `scipy/sparse/dia.py` to see more about how a
`dia_matrix` maps the diagonals data on to the sparse coordinates (`coo`
format).

