title: How to vectorize two nested for loops using numpy?
body: <p>I would like to vectorize part of my code below. Specifically the two for loops that iterate from a <code>range(0,n)</code>. If anything I could just vectorize the <code>j</code> loop because when I include it in my operations the run time exceeds 10 minutes.</p>
  <p>I was thinking that for the exponential part that is called wavea and waveb, to do something with <code>np.reshape</code> and then subtract these so I can just input them on the exponential function but that will give me a matrix and I am not sure how to treat it.</p>
  <pre class="lang-py prettyprint-override"><code>
  epsilon = 0.15
  def delta(x):
      
      # Definining the delta dirac function as a Lorentzian distribution 
      # with a width epsilon
      
      return (1/np.pi)*(epsilon)/(epsilon**2+x**2)
  
  </code></pre>
  <pre class="lang-py prettyprint-override"><code># Doing it inefficiently at first
  
  im = complex(0,1)
  n = 244      
  kdots = 1575      
  
  
  def quantity_abs_squared(Vecs, Evalue, Eigen, rx,ry):
      holder = 0
      
      for k in range(0,kdots):
          
          for i in range(0,n):
  
                  ind_G2a = G2*hexind[i,0] #G1, G2 are fixed (1,2) vectors
                  ind_G1a = G1*hexind[i,1] #hexind is a list of (1,2) tuples of length n
  
                  vectora = Vecs[k][:,i]    #Vecs has shape (1575, 244, 244) 
  
                  wavea = np.exp(im*(ind_G1a[0] + ind_G2a[0])*rx + im*(ind_G1a[1] + ind_G2a[1])*ry)
  
  
              for j in range(0,n):
  
                  ind_G2b = G2*hexind[j,0]
                  ind_G1b = G1*hexind[j,1]
  
                  # Extracting the eigenvector for the eigenvalues
  
                  vectorb = Vecs[k][:,j]
  
                  waveb = np.exp(im*(ind_G1b[0] + ind_G2b[0])*rx + im*(ind_G1b[1] + ind_G2b[1])*ry)
  
                  coeffs = np.dot(vectora, np.conjugate(vectorb))
  
                  exps = np.dot(wavea, np.conjugate(waveb))
  
                  # For every kdots, there are n number of Eigenvalues
                  # For every kdots, there are n number of Eigenvectors of length n
  
                  holder = holder + coeffs*exps*delta((Egrid - Eigen[k][i]))
  
      quant = holder
      
      return (2*2)*quant/kdots
  </code></pre>
  <p>EDIT:</p>
  <p>I have an idea but I am not sure how to bring it to fruition. The most important thing to get rid of is the combinations of the array <code>hexind</code> which is (n,2) in shape and holds the numbers that I need to multiply the <code>G1</code> and <code>G2</code> vectors with. Well what if I do something like this:</p>
  <pre class="lang-py prettyprint-override"><code>
  differences = hexind.reshape([n, 1, 2])-hexind.reshape([1,n,2])
  
  </code></pre>
  <p><code>differences</code> is now a (n,n,2) matrix that holds the index needed for the <code>G1</code> and <code>G2</code> part BUT also gives me the difference that I can input inside of the exponential. Again I am not sure how to loop through this matrix without for loops.</p>
  <p>EDIT #2:</p>
  <p>Typing out a minimal reproducible example. Let us say index is the following array:</p>
  <pre class="lang-py prettyprint-override"><code>
  G1 = np.random.rand((1,2))
  G2 = np.random.rand((1,2))
  
  index = np. array([[ 0,  0],
         [ 1,  1],
         [ 1,  0],
         [ 0, -1],
         [-1, -1],
         [-1,  0],
         [ 0,  1]])
  </code></pre>
  <p>I then change this index in the following way:</p>
  <pre class="lang-py prettyprint-override"><code>
  hexind2 = np.repeat(index,2, axis=0)   # repeat two times for every row
  hexind = np.tile(hexind2,(2,1))
  
  </code></pre>
  <p>Then we can just create some random numpy arrays for the eigenvalues and eigenvectors:</p>
  <pre class="lang-py prettyprint-override"><code>kdots = 100
  n = 28
  
  E = np.random.rand(100, n)
  vec = np.random.rand(100,n,n) + 1j*np.random.rand(100,n,n)
  
  </code></pre>
  <p>Finally I call the function <code>quantity_abs_squared</code> as such:</p>
  <pre class="lang-py prettyprint-override"><code>Eval = 300
  
  rx = np.linspace(-20, 20)
  ry = np.linspace(-20, 20)
  
  Rx,Ry = np.meshgrid(rx,ry)
  
  result = np.zeros((len(Rx), len(Ry)))
  
  
  result = quantity_abs_squared(vec, Eval, E, Rx,Ry) 
  
  </code></pre>