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 % $$ 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}) % $$ 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
numpy.einsum
. (I do not have time for more than this hint or checking whether it actually applies here.) $\endgroup$ – Wrzlprmft Jun 2 '17 at 13:09(np.dot(a[:,:,i],b[:,:,i].T)
does not depend onv
, so could be pre-computed. And, of course (if you useC
-ordering of indices), the indexi
is the fastest running index ina[:,:,i]
orb[:,:,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 fora
andb
as ina = np.ones((g,N,v), dtype = 'complex')
(and change the code accordingly). $\endgroup$ – Walter Jun 5 '17 at 8:14