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?
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.