Consider the function $f: \mathbb{R}^N \to \mathbb{R}$ with
$$ f(r) = ||Ap - b||_2^{2} = p^{\top}A^{\top}Ap - 2p^{\top}A^{\top}b + b^{\top}b $$
where
$A \in \mathbb{R}^{N \times N}$ and $p,b \in \mathbb{R}^N$
The elements of the matrix $A$ are linear functions of the variables $r = (r_1, \ldots, r_N)$, i.e. $a_{i,j} : \mathbb{R}^{N} \to \mathbb{R}$ and $A = (a_{i,j}(r_1, \ldots, r_N))_{i = 1, \ldots, N, j = 1, \ldots, N}$.
I know the partial derivatives of each matrix element, i.e. $ \frac{\partial a_{i,j}}{\partial r_k}$ is known.
Is there any chance for a closed form of $\nabla f$? Any hints or help is really appreciated!
PS: This question is similar to this one. But there the matrix depends only on one variable.