How can I implement Newton-Raphson's method with a function of one vector and one matrix? I have a function $f(\mathbf{u}, \Sigma)$ where $\mathbf{u}$ is a $p \times 1$ vector and $\Sigma$ is a $p \times p$ real symmetric matrix (positive semi-definite).
I somehow successfully computed the partial derivatives $\frac{\partial f}{\partial \mathbf{u}}$ and $\frac{\partial f}{\partial \Sigma}$.
In this case, how do I optimize the function $f$ using Newton-Raphson's method?
=====Details=======
$y = X\mathbf{u} + \frac{1}{2}\operatorname{diag}(X\Sigma X^{T})$
$e = exp(y)$
$f = \mathbf{1}^{T} e$
where $exp$ is component-wise and $\mathbf{1}$ is a vector with ones. $X$ is not symmetric.

I assume that $X\in M_{p,p}$, $diag(V)=[v_{1,1},\cdots,v_{p,p}]^T$ and $y$ is a vector.
Consider $g:(u,\Sigma)\in \mathbb{R}^P\times \mathbb{R}^{p^2}\rightarrow (\frac{\partial f}{\partial u}(u,\Sigma),\frac{\partial f}{\partial \Sigma}(u,\Sigma))\in \mathbb{R}^P\times \mathbb{R}^{p^2}$. To do that, we stack the matrix $\Sigma$ row  by row. For example, when $p=2$, one has $(u,\Sigma)=(u_1,u_2,\Sigma_{11},\Sigma_{12},\Sigma_{21},\Sigma_{2,2})\in\mathbb{R}^6$ and $(\frac{\partial f}{\partial u},\frac{\partial f}{\partial \Sigma})=(\frac{\partial f}{\partial u_1},\frac{\partial f}{\partial u_2},\frac{\partial f}{\partial \Sigma_{1,1}},\frac{\partial f}{\partial \Sigma_{1,2}},\frac{\partial f}{\partial \Sigma_{2,1}},\frac{\partial f}{\partial \Sigma_{2,2}})\in\mathbb{R}^6$.
We want to solve the system of $p+p^2$ equations with $p+p^2$ unknowns $g=0$; clearly, we can use the Newton-Raphson's method:
$(u_{n+1},\Sigma_{n+1})=(u_n,\Sigma_n)-(Dg(u_n,\Sigma_n))^{-1}(g(u_n,\Sigma_n))$.
