The closed-loop dynamics of a linear optimal controller are simple but have interesting properties. From a starting state $\mathbf{v}(0)$ the dynamic can be iterated to reach a final state $\mathbf{v}(N)$ as in


$$
\mathbf{v}(N) = (A - BC)^N\mathbf{v}(0).
$$


I would like to optimize the transition parameter matrix $A$ through gradient descent with rate $\eta$:


$$
A' = A - \eta\frac{df}{dA}
$$


with loss function $f$ is the error between a generic resultant state $(A - BC)^N\mathbf{v}$ and a target state $\mathbf{w}$. $\newcommand\norm[1]{\lVert#1\rVert}$I'm struggling to find the derivative $\frac{\partial{f(A)}}{\partial{A}}$ of this scalar loss function (L2 norm)
$$
f(A) = \norm{(A - BC)^N\mathbf{v} - \mathbf{w}}\_2^2
$$
where $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{m\times n}$, and $\mathbf{v},\mathbf{w}\in\mathbb{R}^{n}$.


Using chain rule, I can write this as
$$
\frac{\partial{f(X)}}{\partial{X}} = 
\frac{\partial{f(X)}}{\partial{g(X)}}
\frac{\partial{g(X)}}{\partial{h(X)}}
\frac{\partial{h(X)}}{\partial{X}}
$$
where
\begin{align\*}
f(X) & {}= g(h(X)) \\
g(X) & {}= \norm{X\mathbf{v} - \mathbf{w}}\_2^2 \\
h(X) & {}= (X-BC)^N.
\end{align\*}


The problem is, I'm not clear on how to take these these intermediate derivatives.



Gradient Descent for Markov Dynamics
Note that $A$ only appears in the combination $M=A-BC$, so the derivative with respect to $A$ equals the derivative with respect to $M$; The function $f(A)$ is given by
$$f(A)=||(A - BC)^Nv - w||\_2^2=(M^Nv)^T M^Nv+w^Tw-2w^TM^Nv.$$
For a simple case, let me first consider a scalar perturbation, $f(A+\epsilon I)=f(A)+\epsilon df/dA$, with derivative
$$\frac{df}{dA} =2N(M^Nv-w)^TM^{N-1}v.$$


Next consider the derivative with respect to a given matrix element $A\_{ij}$ of $A$. The expressions are more lengthy, basically each matrix $M$ gives a separate term so we have a sum $\sum\_{k=1}^N$ instead of the factor $N$:
$$\frac{\partial f}{\partial A\_{ij}}= 2\sum\_{k=1}^N\sum\_{p,q=1}^n (M^Nv-w)\_q (M^{k-1})\_{qi} (M^{N-k})\_{jp} v\_p.$$
Note that the scalar perturbation is the trace of the matrix of elementwise perturbations, $df/dA=\sum\_{i=1}^n \partial f/\partial A\_{ii}$.


