spdiags() function not working as expected in Python

Question

In Matlab/Octave, spdiags([-8037.500 50.000 -12.500], 0:2, 1, 51) gives following output:

(1, 1) -> -8037.5
(1, 2) ->  50
(1, 3) -> -12.500

However, when I use the following in Python, it does not produce a similar result as in Matlab/Octave:

import numpy as np
import scipy as sp
data = array([[-8037.5],
       [   50. ],
       [  -12.5]])
sp.sparse.spdiags(data, np.r_[0:2 + 1].T, 1, 51).toarray()

Python's spdiags() produce following output, which is missing the 50 and -12.5 terms at 1st and 2nd indices:

array([[-8037.5,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ]])

I took a look at this answer to a similar question, but I am not sure where I am going wrong.

Edit:

I am trying to build a matrix A that is made of A_diag1, A_diag2, and A_diag3 as shown below. I have defined A_diag1 and A_diag3 as suggested in the answer.

import numpy as np
import scipy as sp
A_diag1 = np.tile(np.array([-8037.500, 50, -12.5]), (3,1))
A_diag2 = np.reshape(np.repeat([1250, -18505, 1250], 49), (3, 49))
A_diag3 = np.tile(np.array([12.5, -50, 8037.500]), (3,1))
A = np.concatenate((sp.sparse.spdiags(A_diag1, np.r_[0:2 + 1], 1, 51).toarray(), \
              sp.sparse.spdiags(A_diag2, np.r_[0:2 + 1], 49, 51).toarray(), \
              sp.sparse.spdiags(A_diag3, np.r_[48:50 + 1], 1, 51).toarray()), axis=0)

However, five highlighted cells in last 3 rows and columns of A are showing up as zeros/singular as shown in the snapshot below. I expect those highlighted cells, currently showing as zeros, to be non-zero. [You can just copy and paste the above piece of code to reproduce A matrix from which the snapshot shown below is taken.]

EDIT2:

Following code that uses sp.sparse.diags() works as expected. Unlike sp.sparse.spdiags, input argument for the shape of the result (array dimensions) when using sp.sparse.diags() must be in a list.

import numpy as np
import scipy as sp
A_diag1 = np.array([[-8037.500], [50], [-12.5]])
A_diag2 = np.reshape(np.repeat([1250, -18505, 1250], 49), (3, 49))
A_diag3 = np.array([[12.5], [-50], [8037.500]])
A = np.concatenate((sp.sparse.diags(A_diag1, np.arange(0, 2 + 1), [1, 51]).toarray(), \
sp.sparse.diags(A_diag2, np.arange(0, 2 + 1), [49, 51]).toarray(), \
sp.sparse.diags(A_diag3, np.arange(48, 50 + 1), [1, 51]).toarray()), axis=0)

Answer 1

This makes a sparse matrix (51,1), with a value down each row:

In [5]: sparse.spdiags(data,[0,-1,-2],51,1)
Out[5]: 
<51x1 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements (3 diagonals) in DIAgonal format>
In [6]: print(_)
  (0, 0)    -8037.5
  (1, 0)    50.0
  (2, 0)    -12.5

Notice that the spdiags definition:

data : array_like matrix diagonals stored row-wise

Sparse diagonal format stores its data in a matrix, part of which can be 'off-the-screen'. So it's a little tricky to use. I usually create matrices with the coo style of input.

In [27]: M =sparse.spdiags(data,[0,-1,-2],3,3)
In [28]: M.A
Out[28]: 
array([[-8037.5,     0. ,     0. ],
       [   50. ,     0. ,     0. ],
       [  -12.5,     0. ,     0. ]])
In [29]: M.data
Out[29]: 
array([[-8037.5],
       [   50. ],
       [  -12.5]])
In [30]: M.offsets
Out[30]: array([ 0, -1, -2], dtype=int32)

What you want is its transpose (maybe)

In [32]: Mt = M.T
In [33]: Mt.A
Out[33]: 
array([[-8037.5,    50. ,   -12.5],
       [    0. ,     0. ,     0. ],
       [    0. ,     0. ,     0. ]])
In [34]: Mt.data
Out[34]: 
array([[-8037.5,     0. ,     0. ],
       [    0. ,    50. ,     0. ],
       [    0. ,     0. ,   -12.5]])
In [35]: Mt.offsets
Out[35]: array([0, 1, 2], dtype=int32)

So we can recreate Mt with:

sparse.spdiags(Mt.data, Mt.offsets, 3,3)

---

If I save the Octave matrix and load it I get:

In [40]: loadmat('diags')
Out[40]: 
{'__globals__': [],
 '__header__': b'MATLAB 5.0 MAT-file, written by Octave 4.0.0, 2017-10-19 01:24:58 UTC',
 '__version__': '1.0',
 'x': <1x51 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Compressed Sparse Column format>}
In [42]: X=_['x']
In [43]: print(X)
  (0, 0)    -8037.5
  (0, 1)    50.0
  (0, 2)    -12.5

And if I convert it to the dia format I get something like Mt:

In [48]: sparse.dia_matrix(X)
Out[48]: 
<1x51 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements (3 diagonals) in DIAgonal format>
In [49]: print(_)
  (0, 0)    -8037.5
  (0, 1)    50.0
  (0, 2)    -12.5
In [50]: _.data, _.offsets
Out[50]: 
(array([[-8037.5,     0. ,     0. ],
        [    0. ,    50. ,     0. ],
        [    0. ,     0. ,   -12.5]]), array([0, 1, 2]))

The sparse.diags function might be more intuitive:

In [92]: sparse.diags(data, [0,1,2],(1,3))
Out[92]: 
<1x3 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements (3 diagonals) in DIAgonal format>
In [93]: _.A
Out[93]: array([[-8037.5,    50. ,   -12.5]])
In [94]: print(__)
  (0, 0)    -8037.5
  (0, 1)    50.0
  (0, 2)    -12.5

---

In [56]: sp1 = sparse.spdiags(A_diag1, np.r_[0:2 + 1], 1, 51)
In [57]: sp2 = sparse.spdiags(A_diag2, np.r_[0:2 + 1], 49, 51)
In [58]: sp3 = sparse.spdiags(A_diag3, np.r_[48:50 + 1], 1, 51)

(the r_ expressions could also be np.arange(0,3) and np.arange(48,51))

These can joined with sparse.vstack (which combines the coo format attributes)

    In [69]: B = sparse.vstack((sp1,sp2,sp3))
    In [72]: B
    Out[72]: 
    <51x51 sparse matrix of type '<class 'numpy.float64'>'
        with 147 stored elements in COOrdinate format>

In [75]: B.tocsr()[45:, 46:].A
Out[75]: 
array([[  1250.,      0.,      0.,      0.,      0.],
       [-18505.,   1250.,      0.,      0.,      0.],
       [  1250., -18505.,   1250.,      0.,      0.],
       [     0.,   1250., -18505.,      0.,      0.],
       [     0.,      0.,   1250.,      0.,      0.],
       [     0.,      0.,      0.,      0.,      0.]])

matches your snapshot. (I still need to figure out what you are trying to create).

sparse.spdiags(data, diags, m, n) is just another way of calling sparse.dia_matrix((data, diags), shape=(m,n))

Going back to sparse.diags, if you want 3 diagonals, each filled with a value from data we can use:

In [111]: B = sparse.diags(data,[0,1,2],(51,51))
In [112]: B
Out[112]: 
<51x51 sparse matrix of type '<class 'numpy.float64'>'
    with 150 stored elements (3 diagonals) in DIAgonal format>

In [114]: B.tocsr()[:5,:5].A
Out[114]: 
array([[-8037.5,    50. ,   -12.5,     0. ,     0. ],
       [    0. , -8037.5,    50. ,   -12.5,     0. ],
       [    0. ,     0. , -8037.5,    50. ,   -12.5],
       [    0. ,     0. ,     0. , -8037.5,    50. ],
       [    0. ,     0. ,     0. ,     0. , -8037.5]])

In [115]: B.tocsr()[45:, 46:].A
Out[115]: 
array([[   50. ,   -12.5,     0. ,     0. ,     0. ],
       [-8037.5,    50. ,   -12.5,     0. ,     0. ],
       [    0. , -8037.5,    50. ,   -12.5,     0. ],
       [    0. ,     0. , -8037.5,    50. ,   -12.5],
       [    0. ,     0. ,     0. , -8037.5,    50. ],
       [    0. ,     0. ,     0. ,     0. , -8037.5]])

So the sp1 would have to look like

In [117]: B.tocsr()[0,:].todia().data
Out[117]: 
array([[-8037.5,     0. ,     0. ],
       [    0. ,    50. ,     0. ],
       [    0. ,     0. ,   -12.5]])

Answer 2

I have no explanation for your observation (not much of a matlab-user; but i can confirm that octave is doing it like you said), but following scipy's example-usage, you can achieve this result using:

import numpy as np
import scipy.sparse as sp

data = np.tile(np.array([-8037.5, 50., -12.5]), (3,1))
x = sp.spdiags(data, np.arange(3), 1, 51)

print(x)

Output:

(0, 0)        -8037.5
(0, 1)        50.0
(0, 2)        -12.5

The tile-step builds:

[[-8037.5    50.    -12.5]
 [-8037.5    50.    -12.5]
 [-8037.5    50.    -12.5]]

and of course everything is 0-indexing based.