[FIXED] Numpy scalable diagonal matrices

Issue

Assuming I have the variables:

A = 3
B = 2
C = 1

How can i transform them into diagonal matrices in the following form:

np.diag([1, 1, 1, 0, 0, 0])
Out[0]: 
array([[1, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0]])

np.diag([0,0,0,1,1,0])
Out[1]: 
array([[0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0],
       [0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 0]])

np.diag([0,0,0,0,0,1])
Out[2]: 
array([[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, 1]])

I would like this to be scalable, so for instance with 4 variables a = 500, b = 20, c = 300, d = 200 the size of the matrix will be 500 + 20 + 300 + 200 = 1020. What is the easiest way to do this?

Solution

You can achieve even better performance by just allocating the array once, then setting the values all at once by specifying the indices. The indices are fortunately easy to obtain.

import numpy as np

a = [3, 2, 1] # Put your values in a list
s = np.sum(a)
m = np.zeros((len(a), s, s), dtype=int) # Initialize array once
indices = (np.repeat(range(len(a)), a), *np.diag_indices(s, 2)) # Get indices
m[indices] = 1 # Set the diagonals at once
return m

Output:

[[[1 0 0 0 0 0]
  [0 1 0 0 0 0]
  [0 0 1 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 1 0 0]
  [0 0 0 0 1 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 1]]]

Comparing to @Ben Grossmann's answer, with A=3000, B=2000, C=1000 and 100 repeats:

def A():
    '''My solution'''
    a = [3000, 2000, 1000] # Put your values in a list
    s = np.sum(a)
    m = np.zeros((len(a), s, s), dtype=int) # Initialize array once
    indices = (np.repeat(range(len(a)), a), *np.diag_indices(s, 2)) # Get indices
    m[indices] = 1 # Set the diagonals at once
    return m

def B():
    '''Bens solution'''
    A = 3000
    B = 2000
    C = 1000

    n_list = [A,B,C]
    ab_list = np.cumsum([0] + n_list)
    ran = np.arange(ab_list[-1])
    return [np.diag(((a <= ran) & (ran < b)).astype('int')) for a,b in zip(ab_list[:-1], ab_list[1:])]

print(f'Timings:')
timeA = timeit.timeit(A, number=100)
timeB = timeit.timeit(B, number=100)
ratio = timeA / timeB
print(f'This solution: {timeA} seconds')
print(f'Current accepted answer: {timeB} seconds')
if ratio < 1:
    print(f'This solution is {1 / ratio} times faster than Bens solution')
else:
    print(f'Bens solution is {ratio} times faster than this solution')

Output:

Timings:
This solution: 1.6834218999993027 seconds
Current accepted answer: 5.096610300000066 seconds
This solution is 3.027529997086397 times faster than Bens solution

EDIT: Changed the "indices" algorithm to use np.repeat instead of np.concatenate.

Answered By - Naphat Amundsen