I have two points in an N-dimensional space like this: $$\vec r_1 = (x_1,y_1,z_1,w_1,....) \text{ and } \vec r_2 = (x_2, y_2, z_2, w_2, ...)$$
I want to calculate the gradient and Hessian of the Euclidean distance between these two points, against all of the coordinates.
So, distance (d) is $$d = \sqrt{(x_1-x_2)^2 + (y_1 -y_2)^2+ \;...}$$
Now, the gradient is quite simple to calculate: $$\vec{\nabla} d = (\frac{\partial d}{\partial x_1},\frac{\partial d}{\partial y_1},\frac{\partial d}{\partial z_1},...,\frac{\partial d}{\partial x_2},\frac{\partial d}{\partial y_2},....)$$
The first part of this would be $1/d(\vec r_1 - \vec r_2)$ and the second would be $1/d(\vec r_2 - \vec r_1)$. Then I would concatenate these two arrays. This is quite simple to do with python (numpy).
However, I am unable to find a simple matrix/vector based expression for the Hessian matrix. In the Hessian matrix, there would be cross derivative terms ($\frac{\partial^2 d}{\partial x_1 \partial x_2}$).
How do I calculate the Hessian matrix using matrix or vector operations (that would be easy to do on code)?
(I need gradient and Hessian against all the coordinates, not just coordinates of one point)