I think I have a big data (N = 1e6 and dimension = 3 ) scenario. I require to do some matrix manipulation such as einsum, matrix inversion etc several times in my code. To give an idea I want to do something like below.

```
import numpy.random as rd

ndata, kdata = 1e6, 1e5

x = rd.normal(0,1,(ndata, kdata,3,3))

y = rd.normal(0,1,(ndata, kdata,3,3))

```

For small ndata, kdata following would be efficient and convenient approach,

```
xy =  einsum('pqrs, pqsu -> pqru', x, y )

```

Since I have large ndata and kdata above approach becomes a memory bound problem so next bet would be the dot product with nested for loop over ndata and kdata as follows:

```
xyloop1 = np.empty((ndata, kdata, 3, 3))

for j in xrange(ndata):

    for k in xrange(kdata):

        xyloop1[j,k] =  np.dot(x[j,k], y[j,k] )

```

Given what I am taught for loops are nasty in python. Also, I want to use benefits of numpy so thought block matrix approach would be preferable something like following:

```
nstep = 200
ndiv  = ndata/nstep   

kstep = 200
kdiv  = kdata/kstep   

xyloop2 = np.empty((ndata, kdata, 3, 3))

for j in xrange(ndiv):

    ji, jf = j*nstep, (j+1)*nstep     

    for k in xrange(kdiv):

        ki, kf = k*kstep, (k+1)*kstep     

        xyloop2[ji:jf,ki:kf] =  einsum('pqrs, pqsu -> pqru', x[ji:jf,ki:kf], y[ji:jf,ki:kf] )

```

Also, I need these xy or xyloop1 or xyloop2 for my further calculation. So I have to write and read it after every computation. Given the bandwidth of system I/O do you reckon best approach would be approach 3 as it means less I/O and also small number of iterations in compare to approach 2? If you have any other idea or need more info please let me know. 

I am new to the stack so please be gentle with me :). Any help will be highly appreciated. BTW I am trying to solve a mixture modelling problem for a big data. Thanks!