Efficient Matrix construction for a weighted Euclidean distance

I have M points in 2-dimensional Euclidean space, and have stored them in an array X of size M x 2.

I have constructed a cost matrix whereby element ij is the distance d(X[i, :], X[j, :]). The distance function I am using is the standard Euclidean distance weighted by an inverse of the matrix D. i.e d(x,y)= < D^{-1}(x-y) , x-y >. I would like to know if there is a more efficient way of doing this, note I have practically avoided for loops.

import numpy as np Dinv = np.linalg.inv(D) def cost(X, Dinv): Msq = len(X) ** 2 mesh = [] for i in range(2): # separate each coordinate axis xmesh = np.meshgrid(X[:, i], X[:, i]) # meshgrid each axis xmesh = xmesh[1] - xmesh[0] # create the difference matrix xmesh = xmesh.reshape(Msq) # reshape into vector mesh.append(xmesh) # save/append into list meshv = np.vstack((mesh[0], mesh[1])).T # recombined coordinate axis # apply D^{-1} Dx = np.einsum("ij,kj->ki", Dinv, meshv) return np.sum(Dx * meshv, axis=1) # dot the elements