高次元データでカーネル密度推定を行った際に発生するエラーの原因を知りたい

Question

前提・実現したいこと

いつもお世話になっております。
Pythonで密度推定を行っているのですが、
データの次元数を増やすとエラーが発生します。
このエラーなぜ発生するのか原因が理解できず、ご質問させていただきました。
（発生原因をご教示いただけますと幸いです）

発生している問題・エラーメッセージ

次元数が4次元以上になると以下のようなエラーが発生することがあります。
発生割合は5割程度です。（次元数が増えるほどエラー発生確率も増加）

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-11-ce4c335d8dd1> in <module>
      2 xy = np.vstack([a1,a2,a3,a4,a5])
      3 kernel = stats.gaussian_kde(xy)
----> 4 z_est = kernel.evaluate(xy)

~\program\anaconda\lib\site-packages\scipy\stats\kde.py in evaluate(self, points)
    244         result = zeros((m,), dtype=float)
    245 
--> 246         whitening = linalg.cholesky(self.inv_cov)
    247         scaled_dataset = dot(whitening, self.dataset)
    248         scaled_points = dot(whitening, points)

~\program\anaconda\lib\site-packages\scipy\linalg\decomp_cholesky.py in cholesky(a, lower, overwrite_a, check_finite)
     89     """
     90     c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True,
---> 91                          check_finite=check_finite)
     92     return c
     93 

~\program\anaconda\lib\site-packages\scipy\linalg\decomp_cholesky.py in _cholesky(a, lower, overwrite_a, clean, check_finite)
     38     if info > 0:
     39         raise LinAlgError("%d-th leading minor of the array is not positive "
---> 40                           "definite" % info)
     41     if info < 0:
     42         raise ValueError('LAPACK reported an illegal value in {}-th argument'

LinAlgError: 1-th leading minor of the array is not positive definite

該当のソースコード

# import
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
%matplotlib inline

# データ作成（テキトーに生成）
N1 = np.random.normal(size=400)
N2 = np.random.normal(scale=0.5, size=400)
N3 = np.random.normal(scale=0.5, size=400)
N4 = np.random.normal(scale=0.5, size=400)
N5 = np.random.normal(scale=0.5, size=400)
N6 = np.random.normal(scale=0.5, size=400)

a1 = N1+1*N2
a2 = N1-1*N2
a3 = N1+1*N3
a4 = N1-1*N3
a5 = N1+1*N4
a6 = N1-1*N4
a7 = N1+1*N5
a8 = N1-1*N5
a9 = N1+1*N6
a0 = N1-1*N6

# 密度推定（a4以上にするとエラーが発生）
xy = np.vstack([a1,a2,a3])
kernel = stats.gaussian_kde(xy)
z_est = kernel.evaluate(xy)

# 可視化
x = a1
y = a2
plt.scatter(x, y, c=z_est)

試したこと

3次元まではエラーは発生しませんでした。（30回程度試した）
4次元以上では5割ほどの確率でエラーが発生しました。
（次元数が増えるほどエラー発生確率も増加）

補足情報（FW/ツールのバージョンなど）

python3.7
anaconda corei5

以上、ご助力いただけますと幸いです。

Answer 1

自己回答になりますが、
１：この問題を解決することは困難（evaluate時のコレスキー分解による問題）
２：sklearn以外のkernel densityを計算できるライブラリーでは、色々問題があって高次元では算出できない（webではできると書いている人が多いが、3次元以上をまともに調査した結果は無い）

ということで、結局sklaernのkerneldensityをCVさせてbandwidthを決定して使っています。
ただ、sklearnはscipyとかと計算法が違うので若干違いは出ます
（違いが出ても、どちらもどうせ真の値からはずれるので問題ないです）

たぶんこれ以外に現状、あっさり解決する方法はなさそうです。
誰かの参考になれば幸いです。