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!