Partial derivative of a diagonal matrix w.r.t a vector I am trying to find the second partial derivative of the function
$Y=diag\boldsymbol(S)\mathbb P diag \boldsymbol(\beta)diag^{-1}(\mathbb P^{T}\boldsymbol S +\mathbb P^{T}\boldsymbol E + \mathbb P^{T}\boldsymbol I +\mathbb P^{T}\boldsymbol R )\mathbb P^{T}\boldsymbol I$
where $\mathbb P=[(p_{ij})]_{n\times n},\boldsymbol S, \boldsymbol E, \boldsymbol I$ and $\boldsymbol R $ are $n\times 1$ vectors
with respect to $\boldsymbol S$. I am not sure of how to approach this. I was thinking of using the product rule whereby I will take
$$ X=\underbrace{diag\boldsymbol(S)\mathbb P diag \boldsymbol(\beta)}_{Y}\underbrace{diag^{-1}(\mathbb P^{T}\boldsymbol S +\mathbb P^{T}\boldsymbol E + \mathbb P^{T}\boldsymbol I +\mathbb P^{T}\boldsymbol R )\mathbb P^{T}\boldsymbol I}_{Z}$$ so that $$ \frac{\partial X}{\partial \boldsymbol S}=\frac{\partial Y}{\partial \boldsymbol S}Z+Y\frac{\partial Z}{\partial\boldsymbol S}$$
This is proving difficult to achieve as I am not sure if what I am thinking of doing is correct. For example, I was thinking of taking $D=diag(\mathbb P^{T}\boldsymbol S +\mathbb P^{T}\boldsymbol E + \mathbb P^{T}\boldsymbol I +\mathbb P^{T}\boldsymbol R )$ so that $$\frac{\partial Z}{\partial S}=-D\frac{\partial D}{\partial S}D^{-1}\mathbb P^{T}\boldsymbol I$$ where $$\frac{\partial D}{\partial S}=\mathbb P^{T}$$  Is this process correct? Because I am imagigning that this method (for the derivative of the inverse of a matrix) can only be applied when one is finding the derivative w.r.t a matrix and not w.r.t a vector. Could somebody help me obtain $\frac{\partial X}{\partial\boldsymbol S}$ and subsequently $\frac{\partial^{2} X}{\partial\boldsymbol S^{2}}$, $\frac{\partial^{2} X}{\partial\boldsymbol I\partial \boldsymbol E}$. I will appreciate it very much.
 A: $
\def\e{\varepsilon}
\def\l{\left}
\def\r{\right}
\def\lr#1{\l(#1\r)}
\def\d#1{\operatorname{diag}\lr{#1}\,}
\def\D#1{\operatorname{Diag}\lr{#1}\,}
\def\v#1{\operatorname{vec}\lr{#1}\,}
\def\o{{\tt1}}
\def\p{{\partial}}
\def\grad#1#2{\frac{\p #1}{\p #2}}
\def\hess#1#2#3{\frac{\p^2 #1}{\p #2\,\p #3^T}}
\def\c#1{\color{red}{#1}}
\def\E{{\cal E}}
$Notation is going to be very important in answering this question. So first, let's use a convention in which uppercase letters denote matrices and lowercase letters vectors.
Let's further stipulate that an uppercase letter denotes the diagonal matrix generated by the corresponding lowercase vector, e.g.
$$\eqalign{
S = \D{s}, \quad E = \D{e}, \quad etc \\
}$$
We'll also reserve $I$ to denote the identity matrix.
So let's rename your variables as follows
$$({\mathbb P},\beta,S,E,I,R,X) \to (P,b,s,e,a,r,x)$$
and for typing convenience, introduce a new vector
$$w = P^T(s+e+a+r)\quad\implies\quad dw = P^Tds$$
The following commutativity relationship will be key
$$\D{a}b=\D{b}a$$
Write the function, then calculate the differential and thence the gradient
(with respect to $s$).
$$\eqalign{
x &= SPBW^{-1}P^Ta \\
dx &= dS\,PBW^{-1}P^Ta + SPB\,dW^{-1}P^Ta \\
   &= dS\,PBW^{-1}P^Ta + SPB\c{\lr{-W^{-2}dW}}P^Ta \\
   &= \D{ds}PBW^{-1}P^Ta - SPBW^{-2}\D{P^Tds}P^Ta \\
   &= \D{PBW^{-1}P^Ta}ds - SPBW^{-2}\D{P^Ta}P^Tds \\
\grad{x}{s}
   &= \D{PBW^{-1}P^Ta} - SPBW^{-2}\D{P^Ta}P^T \\
}$$
This is a matrix-valued gradient, so the higher order derivatives will be tensors. I'm not sure how you want to handle that.
There are lots of options. You can use tensors, or component-wise index notation, or use Kronecker products to flatten the matrix gradient into a long vector.
$$\\$$

UPDATE
The calculation of the gradient wrt $a$ is similar, and even a bit easier
$$\eqalign{
x &= SPBW^{-1}P^Ta \\
dx &= SPBW^{-1}P^Tda + SPB\,dW^{-1}P^Ta \\
   &= SPBW^{-1}P^Tda + SPB\c{\lr{-W^{-2}dW}}P^Ta \\
   &= SPBW^{-1}P^Tda - SPBW^{-2}\D{P^Tda}P^Ta \\
   &= SPBW^{-1}P^Tda - SPBW^{-2}\D{P^Ta}P^Tda \\
\grad{x}{a}
   &= SPBW^{-1}P^T - SPBW^{-2}\D{P^Ta}P^T \\
}$$
Now calculate the differential of $\,G=\lr{\grad{x}{a}}\,$
with respect to $e$ as the first step in the cross-hessian
$$\eqalign{
dG &= SPB\,\c{ \lr{dW^{-1}}  }\,P^T - SPB\,\D{P^Ta}\c{ \lr{dW^{-2}} }P^T \\
  &= SPB\,\c{ \lr{-W^{-2}dW} }\,P^T - SPB\,\D{P^Ta}\c{ \lr{-2\,W^{-3}dW} }P^T \\
  &= 2\,SPB\,\D{P^Ta}W^{-3}\,\D{P^Tde}\,P^T - SPBW^{-2}\,\D{P^Tde}\,P^T \\
}$$
At this point we need the vectorization function
$(\v{G}\!)$
the Khatri-Rao product $(\boxtimes)$ defined in term of all-ones vectors $(\o)$ and the Kronecker $(\otimes)$ and Hadamard $(\odot)$ products, plus the obscure $\c{\rm relationship}$
$$\eqalign{
&\c{ {\rm vec}(F\,\D{g}\,H) = \left(H^T\boxtimes F\right)g } \\
&H^T\boxtimes F = (H^T\otimes{\o_m})\odot({\o_p}\otimes F)
  \;\in {\mathbb R}^{(mp)\times n} \\
&F \in {\mathbb R}^{m\times n},\quad 
g \in {\mathbb R}^{n},\quad 
H \in {\mathbb R}^{n\times p},\quad
\o_m \in {\mathbb R}^{m} \\
}$$
This yields
$$\eqalign{
dg &= \bigg(
P\boxtimes\Big( 2\,SPB\,\D{P^Ta}W^{-3} - SPBW^{-2} \Big)
\bigg)\,P^T\,de \\
\grad{g}{e} &= \bigg(
P\boxtimes\Big( 2\,SPB\,\D{P^Ta}W^{-3} - SPBW^{-2} \Big)
\bigg)\,P^T \\
}$$
So that's the vectorized version of $\lr{\hess{x}{e}{a}}$
At this point my questions are:
 $\;$ Why do you need complicated tensor quantities like these?
 $\;$ What will they be used for?
