How to efficiently calculate the gaussian kernel matrix in numpy?

Question

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(xz)**2))/(2*sigma**2)) 

This is my current path. Is there a way to use matrix operation for this? X are data points.

Answer 1

Do you want to use the Gaussian kernel, for example, to smooth images? If so, there is a gaussian_filter() function in scipy:

Updated Answer

This should work - although it is still not 100% accurate, it is trying to account for the mass of probability in each grid cell. I think that using the probability density at the midpoint of each cell is slightly less accurate, especially for small cores. See Example https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm .

 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() 

We test it using the example in Figure 3 at the link:

 gkern(5, 2.5)*273 

gives

 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 ]]) 

The original (accepted) answer below is incorrect. The square root is not needed, and the interval definition is incorrect.

 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 

Answer 2

You can simply gaussian filter a simple 2D dirac function , the result is a filter function that is used:

 import numpy as np import scipy.ndimage.filters as fi def gkern2(kernlen=21, nsig=3): """Returns a 2D Gaussian kernel array.""" # create nxn zeros inp = np.zeros((kernlen, kernlen)) # set element at the middle to one, a dirac delta inp[kernlen//2, kernlen//2] = 1 # gaussian-smooth the dirac, resulting in a gaussian filter mask return fi.gaussian_filter(inp, nsig) 

Answer 3

I myself used the accepted answer for image processing, but I find it (and other answers) too dependent on other modules. So here is my compact solution:

 import numpy as np def gkern(l=5, sig=1.): """\ creates gaussian kernel with side length l and a sigma of sig """ ax = np.linspace(-(l - 1) / 2., (l - 1) / 2., l) xx, yy = np.meshgrid(ax, ax) kernel = np.exp(-0.5 * (np.square(xx) + np.square(yy)) / np.square(sig)) return kernel / np.sum(kernel) 

Edit: changed arange value to linspace to handle equal sides

Answer 4

I am trying to improve FuzzyDuck answer here. I think this approach is shorter and easier to understand. Here I use signal.scipy.gaussian to get a two-dimensional Gaussian kernel.

 import numpy as np from scipy import signal def gkern(kernlen=21, std=3): """Returns a 2D Gaussian kernel array.""" gkern1d = signal.gaussian(kernlen, std=std).reshape(kernlen, 1) gkern2d = np.outer(gkern1d, gkern1d) return gkern2d 

matplotlib.pyplot graphics using matplotlib.pyplot :

 import matplotlib.pyplot as plt plt.imshow(gkern(21), interpolation='none') 

Answer 5

The two-dimensional Gaussian kernel matrix can be computed using numpy broadcast,

 def gaussian_kernel(size=21, sigma=3): """Returns a 2D Gaussian kernel. Parameters ---------- size : float, the kernel size (will be square) sigma : float, the sigma Gaussian parameter Returns ------- out : array, shape = (size, size) an array with the centered gaussian kernel """ x = np.linspace(- (size // 2), size // 2) x /= np.sqrt(2)*sigma x2 = x**2 kernel = np.exp(- x2[:, None] - x2[None, :]) return kernel / kernel.sum() 

For small kernel sizes, this should be fast enough.

Note: this makes it easier to change the sigma parameter with respect to the received response.

Answer 6

linalg.norm takes an axis parameter. With a little experiment, I found that I can calculate the rate for all combinations of rows with

 np.linalg.norm(x[None,:,:]-x[:,None,:],axis=2) 

It extends x into a 3d array of all the differences and accepts the norm in the last dimension.

Therefore, I can apply this to your code by adding the axis parameter to your Gaussian :

 def Gaussian(x,z,sigma,axis=None): return np.exp((-(np.linalg.norm(xz, axis=axis)**2))/(2*sigma**2)) x=np.arange(12).reshape(3,4) GaussianMatrix(x,1) 

produces

 array([[ 1.00000000e+00, 1.26641655e-14, 2.57220937e-56], [ 1.26641655e-14, 1.00000000e+00, 1.26641655e-14], [ 2.57220937e-56, 1.26641655e-14, 1.00000000e+00]]) 

Matching:

 Gaussian(x[None,:,:],x[:,None,:],1,axis=2) array([[ 1.00000000e+00, 1.26641655e-14, 2.57220937e-56], [ 1.26641655e-14, 1.00000000e+00, 1.26641655e-14], [ 2.57220937e-56, 1.26641655e-14, 1.00000000e+00]]) 

Answer 7

Based on the answer of Teddy Hartanto. You can simply compute your own one-dimensional Gaussian functions, and then use np.outer to compute the two-dimensional. Very fast and efficient way.

Using the code below you can also use different Sigmas for each dimension

 import numpy as np def generate_gaussian_mask(shape, sigma, sigma_y=None): if sigma_y==None: sigma_y=sigma rows, cols = shape def get_gaussian_fct(size, sigma): fct_gaus_x = np.linspace(0,size,size) fct_gaus_x = fct_gaus_x-size/2 fct_gaus_x = fct_gaus_x**2 fct_gaus_x = fct_gaus_x/(2*sigma**2) fct_gaus_x = np.exp(-fct_gaus_x) return fct_gaus_x mask = np.outer(get_gaussian_fct(rows,sigma), get_gaussian_fct(cols,sigma_y)) return mask 

Answer 8

I tried to use only NumPy. Here is the code

 def get_gauss_kernel(size=3,sigma=1): center=(int)(size/2) kernel=np.zeros((size,size)) for i in range(size): for j in range(size): diff=np.sqrt((i-center)**2+(j-center)**2) kernel[i,j]=np.exp(-(diff**2)/(2*sigma**2)) return kernel/np.sum(kernel) 

You can visualize the result using:

 plt.imshow(get_gauss_kernel(5,1))