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