# python – How to calculate a Gaussian kernel matrix efficiently in numpy?

## python – How to calculate a Gaussian kernel matrix efficiently in numpy?

I myself used the accepted answer for my image processing, but I find it (and the other answers) too dependent on other modules. Therefore, 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)
gauss = np.exp(-0.5 * np.square(ax) / np.square(sig))
kernel = np.outer(gauss, gauss)
return kernel / np.sum(kernel)
```

Edit: Changed arange to linspace to handle even side lengths

Edit: Use separability for faster computation, thank you Yves Daoust.

Do you want to use the Gaussian kernel for e.g. image smoothing? If so, theres a function `gaussian_filter()`

in scipy:

**Updated answer**

This should work – while its still not 100% accurate, it attempts to account for the probability mass within each cell of the grid. I think that using the probability density at the midpoint of each cell is slightly less accurate, especially for small kernels. See https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm for an example.

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

Testing it on the example in Figure 3 from 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 accepted is wrong**

The square root is unnecessary, and the definition of the interval 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
```

#### python – How to calculate a Gaussian kernel matrix efficiently in numpy?

Im trying to improve on FuzzyDucks answer here. I think this approach is shorter and easier to understand. Here Im using `signal.scipy.gaussian`

to get the 2D 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
```

Plotting it using `matplotlib.pyplot`

:

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