Is there an elegant way to take a derivative of this linear algebra problem? I'm a biologist self-learning linear algebra for some applications in my work.  Taking derivatives in non-linear-algebra contexts is quite clear to me, as I can just chain-rule my way through problems... but doing so in a linear-algebra context is a bit of a mystery to me.
I'm trying to take the derivative of the following equation $$y = \vec{d}^TP^TP\vec{\delta},$$ where $$\vec{d} = \left[ \begin{array}\\ d_1 \\ d_2 \\ d_3 \end{array} \right]$$ and $$\vec{\delta} = \left[ \begin{array}\\ o_1^2 -d_1^2 \\ o_2^2 - d_2^2 \\ o_3^2 - d_3^2 \end{array} \right].$$
$P$ is just a 3x3 matrix of constants.
I'd like to find an elegant solution to finding the derivative $y^\prime \left( \vec{d} \right)$, resulting in the Jacobian.
The challenging part is that the vector $\vec{\delta}$ is a composite vector holding $\vec{\delta} = \vec{o} - \vec{d}^{\circ2}$ (not sure on my notation here.. but $\vec{d}^{\circ2}$ is supposed to indicate all the elements of $d$ are squared).
So far, the only way I'm aware of doing this is by expanding everything out into one really large linear formula for the resulting single value, then taking the derivative of that.
My end goal here is to implement this in code... so unnecessarily large formulas are something I'm trying to avoid to keep my code readable.
Is there a more elegant solution to this?
 A: $\DeclareMathOperator{\diag}{diag}$
In what follows $d, x, y$ and $\delta$ denote column vectors.
The function $y(d) = d^T P^T P \delta$ can be rewritten as the dot product $\langle Pd, P\delta \rangle$. This is the composition of three functions:

*

*$y_1 (d) = \begin{bmatrix}d\\ \delta\end{bmatrix}$ which has derivative $Dy_1(d) = \begin{bmatrix}
I \\
-2\diag(d) \end{bmatrix}$, where $\diag(d)$ is the diagonal matrix with $(d_1, d_2, d_3)$ along the diagonal and $I$ is the $3 \times 3$ identity matrix.

*$y_2(x, y) = \begin{bmatrix}Px \\ Py \end{bmatrix}$ which has the derivative $Dy_2(x, y) = \begin{bmatrix} P & 0 \\ 0 & P \end{bmatrix}$.

*$y_3 (x, y) = \langle x, y \rangle$ which has the derivative $Dy_3(x, y) = (y^T, x^T)$.

Putting the three things together using the chain rule gives
$$ \begin{align*}
Dy(d) &= Dy_3(y_2 \circ y_1(d)) Dy_2(y_1(d)) Dy_1(d) \\
&= Dy_3(Pd, P \delta) \begin{bmatrix} P & 0 \\ 0 & P \end{bmatrix} \begin{bmatrix} I  \\ -2\diag(d) \end{bmatrix} \\
&= ((Pd)^T, (P\delta)^T) \begin{bmatrix} P & 0 \\ 0 & P \end{bmatrix} \begin{bmatrix} I  \\ -2\diag(d) \end{bmatrix} \\
&= d^TP^TP - 2 \delta^T P^T P \diag(d).
\end{align*}
$$
A: The elementwise/Hadamard product of two vectors is commutative
$$a\circ b = b\circ a$$
and is equivalent to changing one of the vectors into a diagonal matrix and then performing normal matrix mulitplication, e.g.
$$\eqalign{
 A &= {\rm Diag}(a) \quad\implies\quad &a\circ b = Ab \\
 B &= {\rm Diag}(b) \quad\implies\quad &b\circ a = Ba \\
}$$
The all-ones vector ${\tt1}$ is the identity element for Hadamard multiplication since the corresponding diagonal matrix is the identity matrix. This has an interesting implication:
$$\eqalign{
a\circ{\tt1} &= {\tt1}\circ a \\
A{\tt1} &= Ia = a \;\;(!) \\
}$$
Applying this transformation to the current problem yields
$$\eqalign{
 X &= {\rm Diag}(d),\quad W={\rm Diag}(o) \\
 y &= {\tt1}^TXP^TP(W^2-X^2){\tt1} \\
dy &= {\tt1}^T\,dX\,P^TP(W^2-X^2){\tt1} + {\tt1}^TXP^TP(-2X\,dX){\tt1} \\
 &= {\tt1}^T(W^2-X^2)P^TP\,dX{\tt1} - 2\cdot{\tt1}^TXP^TPX\,dX{\tt1} \\
 &= {\tt1}^T\Big((W^2-X^2)P^TP - 2\,XP^TPX\Big)\,dx \\
\frac{\partial y}{\partial x}
 &= {\tt1}^T\Big((W^2-X^2)P^TP - 2\,XP^TPX\Big) \\
}$$
which is a row vector.
If you prefer your gradient to be a column vector, then transpose the result
$$\eqalign{
\frac{\partial y}{\partial x}
 &= \Big(P^TP(W^2-X^2) - 2\,XP^TPX\Big){\tt1} \\
}$$
Or in terms of the original variables
$$\eqalign{
\frac{\partial y}{\partial x}
 &= P^TP\delta - 2\,{\rm Diag}(d)P^TPd \\
}$$
