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?

Answer 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.

Answer 2

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.

Answer 3

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

Answer 4

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 :)