Fast way of computing covariance matrix of nonstationary kernel in Python

Question

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!

Answer 1

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)