1. ホーム
  2. python

[解決済み] ガウシアンカーネル行列をnumpyで効率的に計算するには?

2022-03-04 16:06:45

質問

def GaussianMatrix(X,sigma):
    row,col=X.shape
    GassMatrix=np.zeros(shape=(row,row))
    X=np.asarray(X)
    i=0
    for v_i in X:
        j=0
        for v_j in X:
            GassMatrix[i,j]=Gaussian(v_i.T,v_j.T,sigma)
            j+=1
        i+=1
    return GassMatrix
def Gaussian(x,z,sigma):
    return np.exp((-(np.linalg.norm(x-z)**2))/(2*sigma**2))

これは私の現在の方法です。これを行列演算で行う方法はないでしょうか?Xはデータ点です。

どのように解決するのですか?

画像のスムージングなどにガウシアンカーネルを使用したいですか?もしそうなら,以下の関数があります。 gaussian_filter() をscipyで実行します。

回答を更新しました

100%正確とは言えませんが、グリッドの各セル内の確率の塊を考慮しようとするもので、これでうまくいくはずです。各セルの中点での確率密度を使用すると、特に小さなカーネルの場合、若干精度が落ちると思います。参照 https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm を例として挙げます。

import numpy as np
import scipy.stats as st

def gkern(kernlen=21, nsig=3):
    """Returns a 2D Gaussian kernel."""

    x = np.linspace(-nsig, nsig, kernlen+1)
    kern1d = np.diff(st.norm.cdf(x))
    kern2d = np.outer(kern1d, kern1d)
    return kern2d/kern2d.sum()

リンク先の図3の例でテストしてみる。

gkern(5, 2.5)*273

与える

array([[ 1.0278445 ,  4.10018648,  6.49510362,  4.10018648,  1.0278445 ],
       [ 4.10018648, 16.35610171, 25.90969361, 16.35610171,  4.10018648],
       [ 6.49510362, 25.90969361, 41.0435344 , 25.90969361,  6.49510362],
       [ 4.10018648, 16.35610171, 25.90969361, 16.35610171,  4.10018648],
       [ 1.0278445 ,  4.10018648,  6.49510362,  4.10018648,  1.0278445 ]])

以下のオリジナル(アクセプトされた)回答は誤りです。 平方根は不要であり、区間の定義も間違っている。

import numpy as np
import scipy.stats as st

def gkern(kernlen=21, nsig=3):
    """Returns a 2D Gaussian kernel array."""

    interval = (2*nsig+1.)/(kernlen)
    x = np.linspace(-nsig-interval/2., nsig+interval/2., kernlen+1)
    kern1d = np.diff(st.norm.cdf(x))
    kernel_raw = np.sqrt(np.outer(kern1d, kern1d))
    kernel = kernel_raw/kernel_raw.sum()
    return kernel