A convenient way is to note that each entry could also be written as :

with above notation the computation could be much easier and:

def G_einsum(M, gamma, R, Z, nu=1.5, length_scale=1.0):

        # --- 1) Determine dimensions ---
        n, d = R.shape  # R, Z each are (n, d)
    
        # --- 2) Build the big array X of shifted points [γ * Z[i] + R[ell]] ---
        X = np.empty((n * n, d))
        # Use tqdm to show progress while constructing X
        for i in tqdm(range(n), desc="Building X array"):
            for ell in range(n):
                X[i * n + ell, :] = gamma * Z[i, :] + R[ell, :]
    
        # --- 3) Use scikit-learn's Matern kernel to compute the (n*n) x (n*n) kernel matrix ---
        kernel = Matern(length_scale=length_scale, nu=nu)
        K_big = kernel(X, X)  # shape: (n*n, n*n)
    
        # --- 4) Reshape K_big into a 4D tensor (n, n, n, n) ---
        K_4D = K_big.reshape(n, n, n, n)
    
        # --- 5) Contract with M on both sides using Einstein summation ---
        M_flat = M.ravel()  # Ensure shape is (n,)
        G = np.einsum('iljm,l,m->ij', K_4D, M_flat, M_flat)
    
        return G