Compute efficiently a 1D function relying on a 2D convolution

Question

Let $X = [0,1]$, $h$ the Gaussian function (i.e. $\forall x \in X, h(x) = e^{-\frac{x}{2}}$) and $p \in L^2(X^2)$

I would like to compute numerically the following function :

$$ \forall x \in X, \quad f(x) = \int_{X^2} p(s_1,s_2) h(s_1 - x) h(s_2 - x) \,\mathrm{d}(s_1,s_2) $$

meaning that i have, a discretized grid for $X$, a matrix of all values of $p$ evaluated on (the grid of) $X^2$ and I would like to compute $f$ at every point of the grid. I already did:

• a numerical discretization of the integral
• a computation of the 2D convolution between the matrix $p$ and the outer product of $h$ (evaluated on the grid). Indeed in some sense $f(x) = (p \ast h \cdot h )(x,x)$ (it's just formal computation, excuse me for this kind of horror!) so taking the diagonal of the computed 2D convolution gives me the $f$ evaluated at the grid

My question is the following: would you have an idea of the operation to be implemented in order not to calculate the whole convolution matrix but just its diagonal? Would you have other ideas to calculate this function?

I specify that I code on Python and that I would like to be able to use the libraries already implemented, so if you have ideas that use SciPy functions it would be great!

Answer 1

I didn't manage to find a suitable formulation for the diagonal. But doing the convolution of the column then the convolution of the rows was enough computation time gain for me

h_vec = gaussienne(X)
convol_row = scipy.signal.convolve(p, h_vec, 'same').T
adj = np.diag(scipy.signal.convolve(convol_row, h_vec, 'same'))