Python example for scipy dia_matrix

Question

Running this code

d0  = np.ones(N)
dp1 = np.ones(N - 1)
dm1 = np.ones(N - 1)

diag = [[d0],[dp1],[dm1]]
offsets = [0,1,-1]

A = dia_matrix( (diag,offsets), shape=(N,N), dtype=float)

I get the following error:

File "/usr/local/lib/python2.7/dist-packages/scipy/sparse/dia.py", line 109, in __init__
self.data = np.atleast_2d(np.array(arg1[0], dtype=dtype, copy=copy))

ValueError: setting an array element with a sequence.

I can't understand what I'm doing wrong! Can someone provide me a correct example to do what I'm trying to do?

Answer 1

When the first argument of dia_matrix has the form (data, offsets), data is expected to be a 2-d array, with each row of data holding a diagonal of the matrix. Since data is a rectangular matrix, some of the elements in data are ignored. Subdiagonals are "left aligned", and superdiagonals are "right aligned". (Specifically, the mapping between data and the sparse matrix A is data[i,j] == A[j - offsets[i], j].) For example, consider the following that will be used to create a 5x5 matrix:

In [28]: data
Out[28]: 
array([[ 1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10],
       [11, 12, 13, 14, 15]])

In [29]: offsets
Out[29]: [0, 1, -1]

data contains three diagonals. The since offset[0] is 0, row 0 of data contains the main diagonal. All 5 elements in this row are used in the matrix. offset[1] is 1, so the data in data[1] becomes the first superdiagonal. Only the values [7, 8, 9, 10] will be used; the first value, 6, is ignored. Similarly, the third row of data gives the first subdiagonal, and only the values [11, 12, 13, 14] are used.

In [30]: a = dia_matrix((data, offsets), shape=(5, 5))

In [31]: a.A
Out[31]: 
array([[ 1,  7,  0,  0,  0],
       [11,  2,  8,  0,  0],
       [ 0, 12,  3,  9,  0],
       [ 0,  0, 13,  4, 10],
       [ 0,  0,  0, 14,  5]])

Your example can be rewritten as follows:

In [32]: N = 5

In [33]: data = np.ones((3, 5))

In [34]: A = dia_matrix((data, offsets), shape=(N, N), dtype=float)

In [35]: A.A
Out[35]: 
array([[ 1.,  1.,  0.,  0.,  0.],
       [ 1.,  1.,  1.,  0.,  0.],
       [ 0.,  1.,  1.,  1.,  0.],
       [ 0.,  0.,  1.,  1.,  1.],
       [ 0.,  0.,  0.,  1.,  1.]])

The dia_matrix docstring has another example.

Alternatively, you can use scipy.sparse.diags to create the matrix. This is useful if you already have code that generates the "correctly" sized diagonals. With diags, you don't have to create the rectangular data matrix. For example,

In [104]: from scipy.sparse import diags

In [105]: d0 = ones(n)

In [106]: dp1 = np.ones(N - 1)

In [107]: dm1 = np.ones(N - 1)

In [108]: d = [d0, dp1, dm1]

In [109]: B = diags(d, offsets, dtype=float)

In [110]: B.A
Out[110]: 
array([[ 1.,  1.,  0.,  0.,  0.],
       [ 1.,  1.,  1.,  0.,  0.],
       [ 0.,  1.,  1.,  1.,  0.],
       [ 0.,  0.,  1.,  1.,  1.],
       [ 0.,  0.,  0.,  1.,  1.]])