Build diagonal matrix without using for loop

Question

I am trying to build the following matrix in Python without using a for loop:

A 
 [[ 0.1  0.2  0.   0.   0. ]
 [ 1.   2.   3.   0.   0. ]
 [ 0.   1.   2.   3.   0. ]
 [ 0.   0.   1.   2.   3. ]
 [ 0.   0.   0.   4.   5. ]]

I tried the fill_diagonal method in NumPy (see matrix B below) but it does not give me the same matrix as shown in matrix A:

B 
 [[ 1.   0.2  0.   0.   0. ]
 [ 0.   2.   0.   0.   0. ]
 [ 0.   0.   3.   0.   0. ]
 [ 0.   0.   0.   1.   0. ]
 [ 0.   0.   0.   4.   5. ]]

Here is the Python code that I used to construct the matrices:

import numpy as np
import scipy.linalg as sp   # maybe use scipy to build diagonal matrix?

#---- build diagonal square array using "for" loop
m = 5

A = np.zeros((m, m))

A[0, 0] = 0.1
A[0, 1] = 0.2

for i in range(1, m-1):
    A[i, i-1] = 1   # m-1
    A[i, i] = 2     # m
    A[i, i+1] = 3   # m+1

A[m-1, m-2] = 4
A[m-1, m-1] = 5

print('A \n', A)

#---- build diagonal square array without loop
B = np.zeros((m, m))

B[0, 0] = 0.1
B[0, 1] = 0.2

np.fill_diagonal(B, [1, 2, 3])

B[m-1, m-2] = 4
B[m-1, m-1] = 5

print('B \n', B)

So is there a way to construct a diagonal matrix like the one shown by matrix A without using a for loop?

Answer 1

There are functions for this in scipy.sparse, e.g.:

from scipy.sparse import diags

C = diags([1,2,3], [-1,0,1], shape=(5,5), dtype=float)

C = C.toarray()

C[0, 0] = 0.1
C[0, 1] = 0.2
C[-1, -2] = 4
C[-1, -1] = 5

Diagonal matrices are generally very sparse, so you could also keep it as a sparse matrix. This could even have large efficiency benefits, depending on the application.

---

The efficiency gains sparse matrices could give you depend very much on matrix size. For a 5x5 array you can't really be bothered I guess. But for larger matrices creating the array could be a lot faster with sparse matrices, illustrated by the following example with an identity matrix:

%timeit np.eye(3000)
# 100 loops, best of 3: 3.12 ms per loop

%timeit sparse.eye(3000)
# 10000 loops, best of 3: 79.5 µs per loop

But the real strength of the sparse matrix data type is shown when you need to do mathematical operations on arrays that are sparse:

%timeit np.eye(3000).dot(np.eye(3000))
# 1 loops, best of 3: 2.8 s per loop

%timeit sparse.eye(3000).dot(sparse.eye(3000))
# 1000 loops, best of 3: 1.11 ms per loop

Or when you need to work with some very large but sparse array:

np.eye(1E6)
# ValueError: array is too big.

sparse.eye(1E6)
# <1000000x1000000 sparse matrix of type '<type 'numpy.float64'>'
# with 1000000 stored elements (1 diagonals) in DIAgonal format>

Answer 2

Notice that the number of 0 is always 3 (or a constant whenever you want to have a diagonal matrix like this):

In [10]:    
import numpy as np
A1=[0.1, 0.2]
A2=[1,2,3]
A3=[4,5]
SPC=[0,0,0] #=or use np.zeros #spacing zeros
np.hstack((A1,SPC,A2,SPC,A2,SPC,A2,SPC,A3)).reshape(5,5)
Out[10]:
array([[ 0.1,  0.2,  0. ,  0. ,  0. ],
       [ 1. ,  2. ,  3. ,  0. ,  0. ],
       [ 0. ,  1. ,  2. ,  3. ,  0. ],
       [ 0. ,  0. ,  1. ,  2. ,  3. ],
       [ 0. ,  0. ,  0. ,  4. ,  5. ]])

---

In [11]:    
import itertools #A more general way of doing it
np.hstack(list(itertools.chain(*[(item, SPC) for item in [A1, A2, A2, A2, A3]]))[:-1]).reshape(5,5)
Out[11]:
array([[ 0.1,  0.2,  0. ,  0. ,  0. ],
       [ 1. ,  2. ,  3. ,  0. ,  0. ],
       [ 0. ,  1. ,  2. ,  3. ,  0. ],
       [ 0. ,  0. ,  1. ,  2. ,  3. ],
       [ 0. ,  0. ,  0. ,  4. ,  5. ]])