I have a collection of arrays
K0,K1,K2....Kn
defined over a 1D array z
I want the following symmetric matrix in the fastest way possible without using for loop.
[np.trapz(K0*K0,z) np.trapz(K0*K1,z) np.trapz(K0*K2,z) np.trapz(K0*K3,z)...]
[ . np.trapz(K1*K1,z) np.trapz(K1*K2,z) np.trapz(K1*K3,z)...]
A = [ . . np.trapz(K2*K2,z) np.trapz(K2*K3,z)...]
[ . . . np.trapz(K3*K3,z)...]
[ . . . . ]
Below is the fastest I could manage (still not fast enough for large n... n>10000
).
I store those set of K
s in a combined array called KK
KK = []
for i in range(n):
KK.append(Ki)
KK = np.array(KK)
A = np.zeros((n,n))
for i in range(n):
A[i,i:] = A[i:,i] = np.trapz((KK[i]*KK[i:]),z)
What is a faster way to do it? I don't care how inelegant or non-pythonic the solution is. I just want to ramp up the speed.
CodePudding user response:
You are using the properties of symmetric, matrix making it very efficient. One way to speed up is to use Numba
import numpy as np
import numba as nb
@nb.njit(cache=True, nogil=True, parallel=True)
def fun(KK,z,n):
A = np.zeros((n,n))
for i in nb.prange(n):
A[i,i:] = A[i:,i] = np.trapz((KK[i]*KK[i:]),z)
return A
Old answers
np.trapz(KK.T[:,:,None]@KK.T[:,None,:],z,axis=0) # using matrix multiplication
np.trapz(np.einsum('ik,jk->ijk',KK,KK),z,axis=2) # Using einsum