# triple cross prouct of tensor

Im trying to compute a triple cross product of vectors a,b, and c in real space and integrate over the entire space. The result is a term in the hamiltonian for an electronic system so there are indices for the band number too (you can think of a and b as gradient of wavefunctions in real space and of c as a vector field): $$h_{nm} = \int \boldsymbol{a}_n(x)\times\boldsymbol{c}(x)\times\boldsymbol{b}_m(x)\,d x %![equation](https://i.stack.imgur.com/VeGe3.gif)$$ where $n$ and $m$ are indices from i.e. 0 to 50. I can write the triple product as $$\boldsymbol{A}\times\boldsymbol{B}\times\boldsymbol{C} = \boldsymbol{B}(\boldsymbol{A}\cdot\boldsymbol{C})-\boldsymbol{C}(\boldsymbol{A}\cdot\boldsymbol{B}) %![eq](https://i.stack.imgur.com/aiEle.gif)$$ which should be faster.

So for the first term in python I'm doing

import numpy as np

N = 50   #  band index
v = 3    # direction x,y,z
g = 30000 # real space grid points

a = np.ones((N,v,g), dtype = 'complex')
b = np.ones((N,v,g), dtype = 'complex')
c = np.ones((v,g), dtype = 'complex')
H_vnn = np.zeros((N,v,N), dtype = 'complex')

for i in range(g):
H_vnn += np.array([c[v,i]*(np.dot(a[:,:,i],b[:,:,i].T)) for v in range(3)])


which, when doesnt run into memory error, takes forever. Im pretty sure this is not the most efficient way of doing it, but the only alternative I can think of is looping over the g index which seems to be even more expansive. Thank to everybody and particular to those who can help me

• where does the integration come in? It looks you're just summing over $n$ and $m$. So should your first equation not just be a double sum? – Walter Jun 2 '17 at 12:39
• Ops you are right; i have written the wrong code. Instead of the double for loop, is one single loop: for i in range(g): – al_j Jun 2 '17 at 12:42
• It may be worth taking a look at numpy.einsum. (I do not have time for more than this hint or checking whether it actually applies here.) – Wrzlprmft Jun 2 '17 at 13:09
• In your code, the sub-expression (np.dot(a[:,:,i],b[:,:,i].T) does not depend on v, so could be pre-computed. And, of course (if you use C-ordering of indices), the index i is the fastest running index in a[:,:,i] or b[:,:,i], so your dot product is over the slower indices, implying extremely inefficient memory access (cache misses all the time). Thus suggests to re-order the indices for a and b as in a = np.ones((g,N,v), dtype = 'complex') (and change the code accordingly). – Walter Jun 5 '17 at 8:14

Try this

import numpy as np

N = 50   #  band index
v = 3    # direction x,y,z
g = 30000 # real space grid points

a = np.ones((g,N,v), dtype = 'complex')
b = np.ones((g,N,v), dtype = 'complex')
c = np.ones((g,v), dtype = 'complex')
H_vnn = np.zeros((N,v,N), dtype = 'complex')

for i in range(g):

• How much then? (btw, I appreciate an up-vote or acceptance tick) Also, are you sure the .T does actually what you want? – Walter Jun 6 '17 at 17:34