$ \def\o{{\tt1}}\def\p{\partial} \def\L{\left}\def\R{\right} \def\LR#1{\L(#1\R)} \def\BR#1{\Big(#1\Big)} \def\bR#1{\big(#1\big)} \def\vecc#1{\operatorname{vec}\LR{#1}} \def\diag#1{\operatorname{diag}\LR{#1}} \def\Sym#1{\operatorname{Sym}\!\BR{#1}} \def\sym#1{\operatorname{Sym}\LR{#1}} \def\Diag#1{\operatorname{Diag}\LR{#1}} \def\trace#1{\operatorname{Tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\c#1{\color{red}{#1}} $ En primer lugar, defina algunas variables nuevas $$\eqalign{ B &= \tfrac 12D \\ y_k &= Mx_k \\ X &= \m{x_1&x_2&\ldots&x_n}\;\in{\mathbb R}^{m\times n} \\ Y &= \m{y_1&y_2&\,\ldots&\,y_n}\;= MX \qiq \c{dY = dM\,X} \\ J &=X\oslash X \qiq J_{ij}=\o\quad\big({\rm all\!-\!ones\;matrix}\big) \\ X &= J\odot X \\ {\cal D}_Y &= {\rm Diag}\BR{\!\vecc{Y}}\,\in{\mathbb R}^{mn\times mn} \\ }$$ donde $\{\odot,\oslash\}$ denotan multiplicación y división elemental/Hadamard.
Entonces reescribe el problema en notación matricial pura y calcular su diferencial. $$\eqalign{ 2B &= (X\odot Y)^TJ + J^T(X\odot Y) - X^TY - Y^TX \\ &= 2\,\Sym{J^T(X\odot{Y}) - X^T{Y}} \\ dB &= \Sym{J^T(X\odot{\c{dY}}) - X^T{\c{dY}}} \\ &= \Sym{J^T(X\odot\LR{\c{dM\,X}}) - X^T{\c{dM\,X}}} \\ }$$ donde $\;\sym{A} \doteq \tfrac 12\LR{A+A^T}$
En este punto, se puede calcular un gradiente por componentes $$\eqalign{ \grad{B}{M_{ij}} &= \Sym{J^T(X\odot\LR{E_{ij}\,X}) - X^T{E_{ij}\,X}} \\ }$$ donde $E_{ij}$ es una matriz de una sola entrada cuya $(i,j)$ elemento es $\o$ y todos los demás son cero.
En términos de la variable original $$\eqalign{ \grad{B}{M_{ij}} &= \frac 12 \LR{\grad{D}{M_{ij}}} \\ }$$ El problema con el gradiente solicitado $\LR{\grad DM}$ es que es un cuarto orden tensor que no puede escribirse en notación matricial estándar.
Actualización
Has añadido una restricción a tu pregunta original, es decir $$M=A^TA \qiq \c{dM=dA^TA+A^TdA}$$ El procedimiento para hallar el nuevo gradiente es sencillo. Se escribe la diferencial, se cambia la variable independiente de $P\to A$ y, a continuación, aislar el gradiente.
He aquí el cálculo del gradiente por componentes $$\eqalign{ dB &= \Sym{J^T(X\odot\LR{\c{dM\,X}}) - X^T\LR{\c{dM\,X}}} \\ &= \Sym{J^T(X\odot\LR{\c{dA^TAX+A^TdA\,X}} - X^T\LR{\c{dA^TAX+A^TdA\,X}}} \\ &= 2 \Sym{J^T\bR{\!\LR{AX}\odot\LR{dA\,X}\!} - X^TA^TdA\,X} \\ \grad{B}{A_{ij}} &= 2 \Sym{J^T\bR{\!\LR{AX}\odot\LR{E_{ij}X}\!}-X^TA^TE_{ij}X} \\ }$$ El documento enlazado es bastante bueno si lo comparamos con la mayoría Aprendizaje automático pero las matemáticas aún están lejos demasiado "a mano". Está dirigido a un público que solo quiere codificar la expresión del gradiente directamente en sus programas Python sin molestarse en comprobar las matemáticas.
Derivar el resultado de su gradiente utilizando la notación matricial sería un esfuerzo hercúleo. Incluso convertir su resultado en notación matricial sería muy difícil, porque la suma anidada viola la convención de suma habitual, ya que contiene una expresión con cuatro $i$ -índices y cuatro $j$ -índices. La notación convencional de índices sólo permite dos de cualquier índice.
Si quieres que alguien te descifre ese embrollo, entonces, como sugirió Steph, deberías publicar una nueva pregunta $-$ Dejaré que responda a esa pregunta.
