Fast way of computing covariance matrix of nonstationary kernel in Python Suppose I have symmetric positive definite covariance function $k:\mathbb{R}\times\mathbb{R}\rightarrow \mathbb{R}$ that is non-stationary (i.e. $k(x,y) \neq g(|x-y|)$ for any function $g$). Is there a fast way in Python given design points $(x_1,\ldots,x_n$) to calculate its covariance matrix $(k(x_i,x_j))_{i,j}$? 

If the covariance function is stationary then we can compute the whole matrix at once using numpy's matrix operations and avoid slow Python loops - e.g. in this. 

Currently my implementation is:

dim = len(X)
kern_mat = np.zeros((dim,dim))
   for i in range(dim):
      for j in range(i+1):
         kern_mat[i,j] = kern(X[i],X[j])
         kern_mat[j,i] = kern_mat[i,j]

Any help with speedups or otherwise is appreciated!
I still would apply numpy's covaranice function using numpy.apply_along_axis 

import numpy as np

x = np.array([[0, 2], [1, 1], [2, 0]]).T
np.apply_along_axis(func1d=np.cov, arr=x, axis=0)