# 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!