In short, my problem is to compute $\frac{d(X^tAX)}{dX}$; where both $A$ and $X$ are matrices.
I have to maximize a negative-log likelihood function $L$
$$L = \frac{1}{2}\ln(|\Sigma|)+\frac{1}{2}\varepsilon^t\Sigma^{-1}\varepsilon;$$
where $\Sigma$ is the covariance matrix, $\varepsilon$ is a column vector of residuals (in my case) and $t$ denotes the transponse. The probelm is that $\Sigma$ is a function of other matrices
$$\Sigma = J^tCJ$$
where, both $J$ and $C$ are matrices. The matrix $J$ is again a function of a vector $\lambda$. I have to maximize the function $L$ w.r.t. vector $\lambda$. I tried chain rule to solve this problem as follows
$$\frac{dL}{d\lambda}=0\implies\frac{dJ}{d\lambda}\frac{d\Sigma}{dJ}\frac{dL}{d\Sigma}=0.$$
In the above equation, $\frac{dJ}{d\lambda}$ and $\frac{d\Sigma}{dJ}$ becomes a tensor. So, I am no longer able to write these quantities on paper. Also, taking derivative of $L$ w.r.t. an individual element of $ \lambda$ does not solve the problem. There are few online resources which suggest using $vec$ operator to deal with tensors; but they heavily use Kronecker product etc. which I have not been able to understand very well, because most of the online material is very opaque.
Can someone please point me to the solution? If someone can refer a good text dealing with a similar problem, that would be great.