4

In my work I often need to aggregate and expand matrices of various quantities, and I am looking for the most efficient ways to do these actions. E.g. I'll have an NxN matrix that I want to aggregate from NxN into PxP where P < N. This is done using a correspondence between the larger dimensions and the smaller dimensions. Usually, P will be around 100 or so.

For example, I'll have a hypothetical 4x4 matrix like this (though in practice, my matrices will be much larger, around 1000x1000)

m=np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])

>>> m
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12],
       [13, 14, 15, 16]])

and a correspondence like this (schematically):

0 -> 0
1 -> 1
2 -> 0
3 -> 1

that I usually store in a dictionary. This means that indices 0 and 2 (for rows and columns) both get allocated to new index 0 and indices 1 and 3 (for rows and columns) both get allocated to new index 1. The matrix could be anything at all, but the correspondence is always many-to-one when I want to compress.

If the input matrix is A and the output matrix is B, then cell B[0, 0] would be the sum of A[0, 0] + A[0, 2] + A[2, 0] + A[2, 2] because new index 0 is made up of original indices 0 and 2.

The aggregation process here would lead to:

array([[ 1+3+9+11,  2+4+10+12 ],
       [ 5+7+13+15, 6+8+14+16 ]])
= array([[ 24, 28 ],
         [ 40, 44 ]])

I can do this by making an empty matrix of the right size and looping over all 4x4=16 cells of the initial matrix and accumulating in nested loops, but this seems to be inefficient and the vectorised nature of numpy is always emphasised by people. I have also done it by using np.ix_ to make sets of indices and use m[row_indices, col_indices].sum(), but I am wondering what the most efficient numpy-like way to do it is.

Conversely, what is the sensible and efficient way to expand a matrix using the correspondence the other way? For example with the same correspondence but in reverse I would go from:

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

to

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

where the values simply get replicated into the new cells.

In my attempts so far for the aggregation, I have used approaches with pandas methods with groupby on index and columns and then extracting the final matrix with, e.g. df.values. However, I don't know the equivalent way to expand a matrix, without using a lot of things like unstack and join and so on. And I see people often say that using pandas is not time-efficient.

Edit 1: I was asked in a comment about exactly how the aggregation should be done. This is how it would be done if I were using nested loops and a dictionary lookup between the original dimensions and the new dimensions:

>>> m=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])  
>>> mnew=np.zeros((2,2))  
>>> big2small={0:0, 1:1, 2:0, 3:1}  
>>> for i in range(4):
...     inew = big2small[i]
...     for j in range(4):
...         jnew = big2small[j]
...         mnew[inew, jnew] += m[i, j]
...
>>> mnew
array([[24., 28.],
       [40., 44.]])

Edit 2: Another comment asked for the aggregation example towards the start to be made more explicit, so I have done so.

4
  • Not really clear what you are looking for. Can you explain how the aggregation happens? For your second question, you can solve it with np.tile(arr, [2,2]) but I don't understand your requirements so it may not be correct. Sep 5, 2022 at 13:35
  • Hello. Because comments don't seem to allow multiline code, or I don't know how to do it, I have added an example to the end of my question to show I how would do it "the long way".
    – JohnFrum
    Sep 5, 2022 at 13:43
  • 1
    Hi. Interesting question, but hard to understand precisely. In your example, could you please add an intermediate step? That is, rewrite array([[ 24, 28 ], [ 40, 44 ]]) but as array([[ 1+3+9+11, etc ], [ etc,etc ]])
    – Stef
    Sep 13, 2022 at 13:03
  • I have done as requested.
    – JohnFrum
    Sep 13, 2022 at 14:02

1 Answer 1

3

Assuming you don't your indices don't have a regular structure I would do it try sparse matrices.

import scipy.sparse as ss
import numpy as np
# your current array of indices
g=np.array([[0,0],[1,1],[2,0],[3,1]])
# a sparse matrix of (data=ones, (row_ind=g[:,0], col_ind=g[:,1]))
# it is one for every pair (g[i,0], g[i,1]), zero elsewhere
u=ss.csr_matrix((np.ones(len(g)), (g[:,0], g[:,1])))

Aggregate

m=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
u.T @ m @ u

Expand

m2 = np.array([[1,2],[3,4]])
u @ m2 @ u.T
6
  • This is pretty great. You are right that my aggregation indices wouldn't have a nice regular structure, in general.
    – JohnFrum
    Sep 5, 2022 at 14:16
  • I should have pointed out in my question that the 4x4 is just an example and it will often be, say, 1000x1000 going down to 100x100. Would the matrix multiplications here be inefficient compared to other methods that there might be? I am assuming that my brute-force nested loops approach would be among the worst approaches.
    – JohnFrum
    Sep 5, 2022 at 14:24
  • 2
    This should be efficient, given the one to many mapping it will be bounded by O(len(m)**2), it scales as well as a copy of the matrix.
    – Bob
    Sep 5, 2022 at 14:41
  • 1
    I benchmarked dense and sparse solutions with a 1% filled matrix (1000,1000) -> (100,100) [ colab notebook] Sep 5, 2022 at 15:09
  • 1
    I've done some testing of my own with my actual code and this works brilliantly. It would never have occurred to me to use sparse matrices to do both the expansion and aggregation, and it is (as far as I can tell) perfect for my needs.
    – JohnFrum
    Sep 5, 2022 at 16:35

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

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

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