Optimizing assignment into an array from various arrays - NumPy

Question

I have four square matrices with dimension 3Nx3N, called A, B, C and D.

I want to combine them in a single matrix. The code with for loops is the following:

import numpy
N = 3
A = numpy.random.random((3*N, 3*N))
B = numpy.random.random((3*N, 3*N))
C = numpy.random.random((3*N, 3*N))
D = numpy.random.random((3*N, 3*N))

final = numpy.zeros((6*N, 6*N))

for i in range(N):
    for j in range(N):
        for k in range(3):
            for l in range(3):
                final[6*i + k][6*j + l] = A[3*i+k][3*j+l]
                final[6*i + k + 3][6*j + l + 3] = B[3*i+k][3*j+l]
                final[6*i + k + 3][6*j + l] = C[3*i+k][3*j+l]
                final[6*i + k][6*j + l + 3] = D[3*i+k][3*j+l]

Is it possible to write the previous for loops in a numpythonic way?

(off-topic) | I was thinking to myself, the only time I've heard someone say "numpythonic" was some other person on SO. But then I checked and that was you too :) — miradulo
– miradulo, Commented Jan 24, 2017 at 16:28
Yes, I'm fighting with different ways of make for loops in a Numpythonic way. The eye is a mistake, It should be zero. — nunodsousa
– nunodsousa, Commented Jan 24, 2017 at 16:30
@Guillaume Jacquenot N is not 3, is any number! But now I understood, is just to be an executable code! Cool — nunodsousa
– nunodsousa, Commented Jan 24, 2017 at 16:35

Divakar · Accepted Answer · 2017-01-24 17:12:42Z

Great problem for practicing array-slicing into multi-dimensional tensors/arrays!

We will initialize the output array as a multi-dimensional 6D array and simply slice it and assign the four arrays being reshaped as 4D arrays. The intention is avoid any stacking/concatenating as those would be expensive specially when working with large arrays by instead working with reshaping of input arrays, which would be merely views.

Here's the implementation -

out = np.zeros((N,2,3,N,2,3),dtype=A.dtype)
out[:,0,:,:,0,:] = A.reshape(N,3,N,3)
out[:,0,:,:,1,:] = D.reshape(N,3,N,3)
out[:,1,:,:,0,:] = C.reshape(N,3,N,3)
out[:,1,:,:,1,:] = B.reshape(N,3,N,3)
out.shape = (6*N,6*N)

Just to explain a bit more, we had :

            |------------------------ Axes for selecting A, B, C, D
np.zeros((N,2,3,N,2,3),dtype=A.dtype)
                  |------------------------- Axes for selecting A, B, C, D

Thus, those two axes (second and fifth) of lengths (2x2) = 4 were used to select between the four inputs.

Runtime test

Approaches -

def original_app(A, B, C, D):
    final = np.zeros((6*N,6*N),dtype=A.dtype)
    for i in range(N):
        for j in range(N):
            for k in range(3):
                for l in range(3):
                    final[6*i + k][6*j + l] = A[3*i+k][3*j+l]
                    final[6*i + k + 3][6*j + l + 3] = B[3*i+k][3*j+l]
                    final[6*i + k + 3][6*j + l] = C[3*i+k][3*j+l]
                    final[6*i + k][6*j + l + 3] = D[3*i+k][3*j+l]
    return final

def slicing_app(A, B, C, D):
    out = np.zeros((N,2,3,N,2,3),dtype=A.dtype)
    out[:,0,:,:,0,:] = A.reshape(N,3,N,3)
    out[:,0,:,:,1,:] = D.reshape(N,3,N,3)
    out[:,1,:,:,0,:] = C.reshape(N,3,N,3)
    out[:,1,:,:,1,:] = B.reshape(N,3,N,3)
    return out.reshape(6*N,6*N)

Timings and verification -

In [147]: # Setup input arrays
     ...: N = 200
     ...: A = np.random.randint(11,99,(3*N,3*N))
     ...: B = np.random.randint(11,99,(3*N,3*N))
     ...: C = np.random.randint(11,99,(3*N,3*N))
     ...: D = np.random.randint(11,99,(3*N,3*N))
     ...: 

In [148]: np.allclose(slicing_app(A, B, C, D), original_app(A, B, C, D))
Out[148]: True

In [149]: %timeit original_app(A, B, C, D)
1 loops, best of 3: 1.63 s per loop

In [150]: %timeit slicing_app(A, B, C, D)
100 loops, best of 3: 9.26 ms per loop

The command numpy.allclose(out, final) returns True: this ensures that it is what we want to perform
@GuillaumeJacquenot Thanks for the confirmation! Added runtime test alongwith the verification part too.

hpaulj · Accepted Answer · 2017-01-24 19:32:01Z

I'll start with a couple of generic observations

For numpy arrays we normally use the [ , ] syntax rather than [][] final[6*i + k][6*j + l] final[6*i + k, 6*j + l]

For new arrays built from others we often use things like reshape and slicing so that we can then add them together as blocks rather than with iterative loops

For an simple example, to take successive differences:

 y = x[1:] - x[:-1]

Regarding the title, 'matrix creation' is clearer. 'load' has more of the sense of reading data from a file, as in np.loadtxt.

=================

So with N=1,

In [171]: A=np.arange(0,9).reshape(3,3)
In [172]: B=np.arange(10,19).reshape(3,3)
In [173]: C=np.arange(20,29).reshape(3,3)
In [174]: D=np.arange(30,39).reshape(3,3)

In [178]: final
Out[178]: 
array([[ 0,  1,  2, 30, 31, 32],
       [ 3,  4,  5, 33, 34, 35],
       [ 6,  7,  8, 36, 37, 38],
       [20, 21, 22, 10, 11, 12],
       [23, 24, 25, 13, 14, 15],
       [26, 27, 28, 16, 17, 18]])

Which can be created with one call to bmat:

In [183]: np.bmat([[A,D],[C,B]]).A
Out[183]: 
array([[ 0,  1,  2, 30, 31, 32],
       [ 3,  4,  5, 33, 34, 35],
       [ 6,  7,  8, 36, 37, 38],
       [20, 21, 22, 10, 11, 12],
       [23, 24, 25, 13, 14, 15],
       [26, 27, 28, 16, 17, 18]])

bmat uses a mix of hstack and vstack. It also produces a np.matrix, hence the need for .A. @Divakar's solution is bound to be faster.

This does not match with N=3. 3x3 blocks are out of order. But expanding the array to 6d (as Divakar does), and swapping some axes, puts the sub blocks into the the right order.

For N=3:

In [57]: block=np.bmat([[A,D],[C,B]])
In [58]: b1=block.A.reshape(2,3,3,2,3,3)
In [59]: b2=b1.transpose(1,0,2,4,3,5)
In [60]: b3=b2.reshape(18,18)
In [61]: np.allclose(b3,final)
Out[61]: True

In quick time tests (N=3), my approach is about half the speed of slicing_app.

As a matter of curiosity, bmat works with a string input: np.bmat('A,D;C,B'). That's because np.matrix was trying, years ago, to give a MATLAB feel.

epattaro · Accepted Answer · 2017-01-24 16:53:00Z

1

you can just concat em

concat A and B horizontally concat C and D horizontally

concat the conjunction of AB with the conjucntion of CD vertically

example:

AB = numpy.concatenate([A,B],1)
CD = numpy.concatenate([C,D],1)
ABCD = numpy.concatenate([AB,CD],0)

i hope that helps :)

answered Jan 24, 2017 at 16:53

epattaro

2,4981 gold badge21 silver badges32 bronze badges

3 Comments

Guillaume Jacquenot Over a year ago

The result is not the one expected: numpy.allclose(ABCD, final) returns false

epattaro Over a year ago

did you try just to print ABCD? it shld be the desired matrix. you can (should) check the values manually printing the matrix to see if it corresponds to what is expected. try with minor 2x2 matrixes as an example/test

hpaulj Over a year ago

You are on the right track. My answer with bmat does this sort of concatenate. But the OP's layout is more complicated, requiring some reshaping and axis swapping.

B. M. · Accepted Answer · 2017-01-25 10:16:05Z

1

Yet another way to do that, using view_as_blocks :

from skimage.util import view_as_blocks

def by_blocks():
    final = numpy.empty((6*N,6*N))
    a,b,c,d,f= [view_as_blocks(X,(3,3)) for X in [A,B,C,D,final]]
    f[0::2,0::2]=a
    f[1::2,1::2]=b
    f[1::2,0::2]=c
    f[0::2,1::2]=d
    return final

You just have to think block by block, letting view_as_blocks manage strides and shapes for you. It is as fast as other numpy solutions.

edited Jan 25, 2017 at 10:16

answered Jan 24, 2017 at 21:19

B. M.

18.7k2 gold badges40 silver badges56 bronze badges

Collectives™ on Stack Overflow

Optimizing assignment into an array from various arrays - NumPy

4 Answers 4

2 Comments

Comments

3 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

4 Answers 4

2 Comments

Comments

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related