"Matrix" multiplication of Dirac delta distributions? Forgive me if I am not very rigorous, part of the following question is to hopefully clear up any inconsistencies coming from my lack of rigor. In any case I will try to make things clear enough that you at least understand what I am trying to ask.
I would like to consider a type of matrix multiplication between distributions $M(x,y)$ and $N(x,y)$, with $x,y\in S$, defined as
$$\big(M*N\big)(x,y)=\int_Sdz\,M(x,z)N(z,y)$$
More specifically, I'd like to consider the matrix multiplication between $f(x)\frac{\partial}{\partial x}\delta (x-y)$ and $g(x)\frac{\partial}{\partial x}\delta (x-y)$, with $S=\mathbb{R}$. Assume that $f$ and $g$ are differentiable on $\mathbb{R}$. I seem to be getting different answers depending on how I evaluate the integral:
$$\int_{-\infty}^\infty dz\,\Big(f(x)\frac{\partial}{\partial x}\delta (x-z)\Big)g(z)\frac{\partial}{\partial z}\delta (z-y)=f(x)\frac{\partial}{\partial x}\Big(\int_{-\infty}^{\infty}dz\,\delta (x-z)g(z)\frac{\partial}{\partial z}\delta (z-y)\Big)=f(x)\frac{\partial}{\partial x}\Big(g(x)\frac{\partial}{\partial x}\delta(x-y)\Big). \tag{1}$$
The right hand side is obtained by using the definition of the first dirac delta distribution in the integrand. If I instead apply the definition to the second dirac delta distribution, I find
$$\int_{-\infty}^\infty \,dz\Big(f(x)\frac{\partial}{\partial x}\delta (x-z)\Big)g(z)\frac{\partial}{\partial z}\delta (z-y)=-f(x)\int_{-\infty}^{\infty}dz\,\frac{\partial}{\partial z}\Big(g(z)\frac{\partial}{\partial x}\delta (x-z)\Big)\delta (z-y)=-f(x)\frac{\partial}{\partial y}\Big(g(y)\frac{\partial}{\partial x}\delta (x-y)\Big). \tag{2}$$
These two expression don't appear to be equal. Are they equal? And if not, what went wrong?
 A: Here's my 2 a.m. take on the integral:
$$
\int_{-\infty}^{\infty} dz \, \left( f(x)\frac{\partial}{\partial x}\delta (x-z) \right) \left( g(z)\frac{\partial}{\partial z}\delta(z-y) \right)
\\
= \int_{-\infty}^{\infty} dz \, f(x) \, \delta'(x-z) \, g(z) \, \delta'(z-y)
\\
= f(x) \, \int_{-\infty}^{\infty} dz \, \delta'(x-z) \, g(z) \, \delta'(z-y)
\\
= -f(x) \, \int_{-\infty}^{\infty} dz \, \frac{\partial}{\partial z}\left( \delta'(x-z) \, g(z) \right) \, \delta(z-y)
\\
= -f(x) \, \int_{-\infty}^{\infty} dz \, \left( -\delta''(x-z) \, g(z) + \delta'(x-z)\,g'(z) \right) \, \delta(z-y)
\\
= -f(x) \left( -\delta''(x-y) \, g(y) + \delta'(x-y) \, g'(y) \right)
\\
= f(x) \, \delta''(x-y) \, g(y) - f(x) \, \delta'(x-y) \, g'(y)
$$
A: There is a lot that can be said about "integral operators" of the form $Tf(x)=\int_{\mathbb R^n} k(x-y)\,f(y)\;dy$, where $k$ is meets a suitable smoothness condition away from $0$, and its Fourier transform (e.g.) has suitable decay. These are "singular integral operators of convolution type", for a search-able key-phrase.
With suitable tweaks of the conditions on the kernel, these give bounded operators $L^p\to L^p$, and has other good/interesting features.
The prototype is the Hilbert transform, with $k(x)=1/x$. This already illustrates the point that the "integral" defining the transform cannot quite be a literal integral, but must be a principal value version of the integral. The simplest $n$-dimensional analogue is the Riesz transforms, with $k(x)=x_i/|x|^{n+1}$. On $\mathbb R^2$, for $m=1,2,3,\ldots$, $k((x,y))=(x\pm iy)^m/(x^2+y^2)^{{m\over 2}+1}$ gives another family of this sort.
For kernels $K(x,y)$ not of the form $k(x-y)$, the key-phrase "Calderon-Zygmund kernel/operator" is a good search-phrase.
A: Hint: One may show that the RHS of OP's eqs. (1) & (2) are equal by using the following 2 identities
$$ \left(\frac{\partial}{\partial x}+\frac{\partial}{\partial y}\right)\delta(x-y)~=~0, \tag{A}$$
and
$$ \left(g(x)-g(y)\right)\delta(x-y)~=~0, \tag{B}$$
where $g$ is a test function.
