[FIXED] Fastest way to create a symmetric matrix in python, with elements as below

Issue

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 Ks 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.

Solution

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

Answered By - Murali