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)

hpaulj · Accepted Answer · 2017-10-19T02:07:49.787

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]])

Thanks. Can you please take a look at why matrix `A`, as described under the edit, is not giving correct output as expected? — user11, Oct 19 '17 at 01:26
`sparse.diags` is probably the function that you want to use, not `spdiags`. — hpaulj, Oct 19 '17 at 02:08

sascha · Answer 2 · 2017-10-18T21:51:02.603

1

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.

edited Oct 18 '17 at 21:51

answered Oct 18 '17 at 21:42

sascha

32,238
6
68
110

Can you suggest how to make it work when added at the end like following because it doesn't work when concatenated at the end? `A = np.concatenate(your_solved_matrix, \ sp.sparse.spdiags(diagA_innerPts, np.r_[0:2 + 1], 49, 51).toarray(), \ your_solved_matrix ), axis=0)` – user11 Oct 18 '17 at 21:56
Can't follow you exactly, but you are calling numpy's concatenate (for dense numpy-arrays), while you probably want scipy's ```sparse.hstack``` or ```sparse.vstack``` Be careful when mixing those two very different types. – sascha Oct 18 '17 at 22:02
Please let me know if you still need some more clarification. Thanks. – user11 Oct 19 '17 at 01:11
@user11 To be honest... the general design of your matrix is unclear to me. Try to be more formal within your edit. I'm not sure if that's something for this question or a new one... I won't look at it today, maybe tomorrow. But especially the cropped example-output is not enough. – sascha Oct 19 '17 at 01:13
Please see the edit now. I have added more clarification and you can just reproduce the matrix `A` by copying and running those few lines of code under the edit. – user11 Oct 19 '17 at 01:25

spdiags() function not working as expected in Python

2 Answers2