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!

  • $\begingroup$ Have you tried scipy.signal.convolve? $\endgroup$
    – nicoguaro
    Feb 3 at 19:47
  • $\begingroup$ yes, i was more thinking about how to formulate the problem tbf. But I think I got my answer $\endgroup$
    – Bast
    Feb 3 at 21:42
  • 2
    $\begingroup$ If you have your answer you can add it and even accept it. It might be useful to someone else in the future. $\endgroup$
    – nicoguaro
    Feb 3 at 22:36

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.