17

Suppose you have either two arrays:

index = [1, 2, 3]
counts = [2, 3, 2]

or a singular array

arr = [1, 1, 2, 2, 2, 3, 3]

How can I efficiently construct the matrix

[
  [1, 1, 0, 0, 0, 0, 0],
  [1, 1, 0, 0, 0, 0, 0],
  [0, 0, 2, 2, 2, 0, 0],
  [0, 0, 2, 2, 2, 0, 0],
  [0, 0, 2, 2, 2, 0, 0],
  [0, 0, 0, 0, 0, 3, 3],
  [0, 0, 0, 0, 0, 3, 3]
]

with NumPy?

I know that

square = np.zeros((7, 7))
np.fill_diagnol(square, arr) # see arr above

produces

[
  [1, 0, 0, 0, 0, 0, 0],
  [0, 1, 0, 0, 0, 0, 0],
  [0, 0, 2, 0, 0, 0, 0],
  [0, 0, 0, 2, 0, 0, 0],
  [0, 0, 0, 0, 2, 0, 0],
  [0, 0, 0, 0, 0, 3, 0],
  [0, 0, 0, 0, 0, 0, 3]
]

How do I "expand" the diagonal by n where n is counts[index-1] for the values specified by index[I]

tmp = np.array((arr * N)).reshape((len(arr), len(arr)) 
np.floor( (tmp + tmp.T) / 2 ) # <-- this is closer


array([[1., 1., 1., 1., 1., 2., 2.],
       [1., 1., 1., 1., 1., 2., 2.],
       [1., 1., 2., 2., 2., 2., 2.],
       [1., 1., 2., 2., 2., 2., 2.],
       [1., 1., 2., 2., 2., 2., 2.],
       [2., 2., 2., 2., 2., 3., 3.],
       [2., 2., 2., 2., 2., 3., 3.]])

This gets what I want, but probably doesn't scale that well?

riffled = list(zip(index, counts))
riffled
# [(1, 2), (2, 3), (3, 2)]
a = np.zeros((len(arr), len(arr))) # 7, 7 square
last = 0 # <-- keep track of current sub square
for i, c in riffled:
    a[last:last+c, last:last+c] = np.ones((c, c)) * i 
    last += c # <-- shift square

yield

array([[1., 1., 0., 0., 0., 0., 0.],
       [1., 1., 0., 0., 0., 0., 0.],
       [0., 0., 2., 2., 2., 0., 0.],
       [0., 0., 2., 2., 2., 0., 0.],
       [0., 0., 2., 2., 2., 0., 0.],
       [0., 0., 0., 0., 0., 3., 3.],
       [0., 0., 0., 0., 0., 3., 3.]])
2
  • 4
    np.equal.outer(arr, arr) * arr Oct 26, 2021 at 17:32
  • @user3483203 also works! thank you
    – SumNeuron
    Oct 26, 2021 at 17:36

3 Answers 3

10

You can use scipy.linalg.block_diag to make that work:

import numpy as np
import scipy.linalg as linalg


a = 1*np.ones((2,2))
b = 2*np.ones((3,3))
c = 3*np.ones((2,2))

superBlock = linalg.block_diag(a,b,c)

print(superBlock)

#returns
#[[1. 1. 0. 0. 0. 0. 0.]
# [1. 1. 0. 0. 0. 0. 0.]
# [0. 0. 2. 2. 2. 0. 0.]
# [0. 0. 2. 2. 2. 0. 0.]
# [0. 0. 2. 2. 2. 0. 0.]
# [0. 0. 0. 0. 0. 3. 3.]
# [0. 0. 0. 0. 0. 3. 3.]]

if you want to get there from a list of values and a list of counts you can do this:

values = [1,2,3]
counts = [2,3,2]

mats = []
for v,c in zip(values,counts):
    thisMatrix = v*np.ones((c,c))
    mats.append( thisMatrix )


superBlock = linalg.block_diag(*mats)


print(superBlock)
4
  • 1
    So mats = [np.full((c, c), v) for v, c in zip(values, counts)]? Oct 26, 2021 at 17:42
  • This is a good approach. Oct 26, 2021 at 17:42
  • 2
    This is perfect if the user only needs to work with dense matrices after this creation step. However if all of the data can be sparsely described like the example, there could also be some benefit in using the sparse versions: scipy.sparse.block_diag docs.scipy.org/doc/scipy/reference/generated/… Oct 26, 2021 at 17:56
  • @MikeHolcomb want to propose an edit? that could be great to have a better anwer (or make your own i'll upvote it) Oct 26, 2021 at 17:57
4

Here is a generic solution.

starting from the index/count:

index = [1, 2, 1]
counts = [2, 3, 2]

arr = np.repeat(index, counts)
arr2 = np.repeat(range(len(index)), counts)
np.where(arr2 == arr2[:, None], arr, 0)

output:

array([[1, 1, 0, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 0, 0, 0, 1, 1],
       [0, 0, 0, 0, 0, 1, 1]])

starting from the array version:

arr = np.array([1, 1, 2, 2, 2, 1, 2])

arr2 = np.cumsum(np.diff(arr,prepend=np.nan) != 0)
np.where(arr2 == arr2[:, None], arr, 0)

output:

array([[1, 1, 0, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 0, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 0, 2]])
0
3

Try broadcasting:

idx = np.repeat(np.arange(len(counts)), counts)
np.where(idx==idx[:,None], arr, 0)
# or
# arr * (idx==idx[:,None])

Output;

array([[1, 1, 0, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 2, 2, 2, 0, 0],
       [0, 0, 0, 0, 0, 3, 3],
       [0, 0, 0, 0, 0, 3, 3]])
7
  • 3
    That works! Note to readers arr in question is a Python list so for this to work, convert it to numpy via np.array(arr)
    – SumNeuron
    Oct 26, 2021 at 17:33
  • 1
    @SumNeuron. You can just do np.reshape(arr, (-1, 1)) == arr to avoid an explicit conversion. It's cheaper to do an explicit conversion though. Oct 26, 2021 at 17:34
  • 1
    Note that this only works if the blocks in arr are different (e.g., arr = np.array([1, 1, 2, 2, 2, 1, 1]) wouldn't work)
    – mozway
    Oct 26, 2021 at 17:35
  • @mozway just out of curiosity do you have an adjustment for that use case?
    – SumNeuron
    Oct 26, 2021 at 17:36
  • 1
    @mozway see updated answer. We simply don't use npindex as a reference. Oct 26, 2021 at 17:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.