Right, thank you! If you have any favorite references for matrix calculus, I would love to know.

For your full solution, I think we can write this in matrix form where $\frac{\partial{f}}{\partial{A_{ij}}}=2\sum_{k=1}^N{x_k^Ty_k}$ and $x_k=(M^{k-1})^T(M^Nv-w)$ and $y_k=M^{N-k}v$. This makes the connection to the scalar case a little clearer perhaps. Both rank-1, and it looks like in the latter case the derivative finds interactions between different time steps of the Markov sequence, through $M^{k-1}M^{N-k}$...? So the full solution (but do confirm) is: $2\sum_{k=1}^N{(M^Nv-w)^TM^{k-1}M^{N-k}v}$.

Thanks Professor! I've numerically confirmed the online calculator and your scalar perturbation results. I discovered the online calculator is performing elementwise exponentiation, explaining the discrepancy. I'm working now to convert your elementwise derivative to a matrix expression for computation. It seems a tensor product might help. Context: I'm working on gradient descent for Markov dynamical systems. If you have any references or suggestions to help me proceed with similar calculations, I would be very grateful. I'm hoping to compute second-order (Newton's method) as well.

I see, this makes sense. Thank you. This calculator is interesting-- I'm not able to reproduce the scalar perturbation but it produces a nice-looking solution for a general derivative with input norm2((A-B*C).^N*v-w)^2. I'm curious to see if the answer there is interpretable to you. It's difficult for me to form an intuition about it unfortunately.

Thanks Professor Beenakker. I'm a little confused why the gradient here is scalar, when I expected a gradient in terms of the elements of A? That is, I expect the gradient here to tell me the change in f(A) with respect to each element of A. What am I missing?

fixed, thank you