python - numpyでのベクトル化（二乗）マハラノビス距離

Question

X (nxd)、Y (mxd)、正定値 L (dxd) があります。D_ij が (X_i - Y_i) * L * (X_i - Y_i).T である D を計算したい。n と m は約 250 です。d は約 10^4 です。

を使用できますscipy.spatial.distance.cdistが、これは非常に遅いです。

scipy.spatial.distance.cdist(X, Y, metric='mahalanobis', VI=L)

この質問に対する Dougal の回答を見て、試してみました

    diff = X[np.newaxis, :, :] - Y[:, np.newaxis, :]
    D = np.einsum('jik,kl,jil->ij', diff, L, diff)

これも非常に遅いです。

この計算をベクトル化するより効率的な方法はありますか?