Outline
- Discussion
- Using Differentials
- Using the NonCommutativeMultiply package NCAlgebra
- Using MatrixD
- Using xAct
Discussion (code sections below)
There is an ambiguity that has not been mentioned. The shape of the matrix derivative depends on the choice of the convention as explained in this wikipedia page. The ambiguity does not matter once the convention is understood or one considers differentials (considered in the first code example below) rather than derivatives.
Moreover, the formulas can depend on whether the matrix has symmetry. Indeed, a matrix with symmetry has fewer components and thus fewer variables to include among the derivatives. For example, for a generic matrix, the partial derivative with respect to $ M_{\text{kl}}$ keeping all other matrix components fixed has the following fundamental formula:
$$\frac{\partial M_{\text{ij}}}{\partial M_{\text{kl}}} \bigg\rvert_{M_{\text{rh}},\; \text{r}\neq\text{k},\text{h}\neq\text{l}}=\mathbb{J}_{\text{ijkl}}$$
where $\mathbb{J}_{\text{ijkl}}\equiv \delta _{\text{ik}} \delta _{\text{jl}}$.
But if $M$ is symmetric, then one can not take $\frac{\partial M_{\text{ij}}}{\partial M_{\text{kl}}} \bigg\rvert_{M_{\text{rh}},\; \text{r}\neq\text{k},\text{h}\neq\text{l}}$ as in the set {$\text{r}\neq\text{k},\text{h}\neq\text{l}$} there is $M_{\text{lk}}$ which can not be held fixed while $M_{\text{kl}}$ varies as they are the same variable.
This issue with symmetric matrices caused confusion even as late as 2021 in this math stack exchange and was also the reason for this 2019 paper that was updated in 2020. My summary of the many discussions I found was that although the partial derivatives of $M_{\text{lk}}$ and $M_{\text{kl}}$ can not be taken independently, a trick that works is to pretend they can be taken independently and take an average such that $\mathbb{J}_{\text{ijkl}}\equiv \frac{1}{2}\;( \delta _{\text{ik}} \delta _{\text{jl}}+\delta _{\text{il}} \delta _{\text{jk}} )$.
That said, depending on the final result, it might not matter which $\mathbb{J}_{\text{ijkl}}$ we choose.
In the following, we will not consider matrices with symmetry.
Using differentials
The code below shows explicitly how to obtain matrix derivatives in the restricted class of scalar functions defined by a Trace as in the example of the original question. However, it is restricted to this class although including Det might be possible after a few modifications using $\text{det}(A)=\exp(\text{tr}(\log(A)))$. For more general-purpose functions see the following code sections that offer packages.
We recall the definition of the differential of a scalar function $f$ evaluated on a matrix $A$.
$$Df_A[H]=\langle H, \frac{\partial f}{\partial A} \rangle \underset{H\rightarrow 0}{=} f(A+H)-f(A) +O(H^2)$$
For a scalar function:
$$\langle H,\frac{\partial f}{\partial A} \rangle=\text{Tr}\left[H^{T}.\frac{\partial f}{\partial A}\right]$$
An aside
Notice however that one could have equally defined $\frac{\partial f}{\partial A}$ such that (alternative definition)
$$\langle H,\frac{\partial f}{\partial A} \rangle=\text{Tr}\left[ \frac{\partial f}{\partial A}.H\right]$$
In that last case $\frac{\partial f}{\partial A}$ is the transpose of the first definition. That definition has the advantage of removing the Transpose from the definition and is formally a better generalization of the case where A is a $n\times1$ column vector. Indeed, without a scalar product, the natural quantity to consider is the total derivative of the scalar function. If a scalar product is present then one can consider instead the gradient which is dual to total derivative. In Euclidean geometry, the two are related by a transpose see the same link mentioned before. However, the difference between the two is more important in non-euclidean geometry. As an example of how correctly identifying the two can even change an algorithm see this link on the natural gradient as the mathematically "right" way to use gradient descent in neural networks.
For our purpose, we will consider the first definition as the second definition leads to a transpose of $A$ in the final result.
The first step of the method consists of replacing $A$ with $A+t H$ and Taylor expanding about $t=0$ in order to obtain the differential. This is implemented in the code below:
$Assumptions = {t ∈ Reals}
DfA = Tr[Transpose[Transpose@X . A] . Transpose@X . A] /.
A -> A + t*H // TensorExpand // ExpandAll //
ReplaceAll[Tr[Plus[a__]] :> Tr /@ a] //
ReplaceAll[Tr[t^p_.*a_] :> t^p*Tr[a]] // D[#, t] & //
ReplaceAll[t -> 0]
out: (* Tr[Transpose[Transpose[X, {2, 1}] . A, {2, 1}].Transpose[X, {2, 1}] . H] + Tr[Transpose[Transpose[X, {2, 1}] . H, {2, 1}].Transpose[X, {2, 1}] . A] *)
To simplify the above expression and identify the derivative, we replace Transpose and dot with a custom transpose and dot. These custom functions have properties that will allow the expression to be simplified. Moreover, these functions will place the expression in a form that allows us to obtain the derivative. We also define an identity matrix that acts as an identity element for dot. The definitions and properties of these functions are given in the code below (notice that H is hard-coded into some of the definitions below):
totranspose = Transpose[a_, {2, 1}] :> transpose[a];
todot = Dot -> dot;
transpose[dot[a_ , b_]] := dot[transpose[b], transpose[a]]
transpose[transpose[a_]] := a
transpose[id] ^:= id
id /: dot[s___, id, g___] := dot[s, g]
dot /: Tr[dot[a___, H, b___]] := Tr[dot[transpose@H, If[{a} === {}, id, transpose@dot@a], If[{b} === {}, id, transpose@dot@b]]]
dot /: Tr[dot[a__, transpose@H, b___]] := Tr[dot[transpose[H], dot@b, If[{a} === {}, id, dot@a]]]
SetAttributes[dot, {Flat, OneIdentity}]
Then changing Dot to dot and Transpose to transpose we obtain:
DfA //. {todot, totranspose}
out: 2 Tr[dot[transpose[H], X, transpose[X], A]]
Finally, we define rules to replace transpose with Transpose and dot with Dot and to extract the derivative:
toDot = dot -> Dot;
toTranspose = transpose -> Transpose;
obtainDerivative = a_. Tr[Transpose@H . b_ ] :> a*b;
This leads to the final result:
DfA //. {todot, totranspose} //. {toDot,toTranspose} /. obtainDerivative // TeXForm
$$2 X.X^T.A $$
Using the NonCommutativeMultiply package NCAlgebra
See @dtn's answer and also this answer. Notice that @dtn's answer using the NCAlgebra package gives the transpose of the result above.
Using MatrixD
The answer on this page provides a function for matrix derivatives. After installing the package one can use the following:
MatrixD[Tr[Transpose[[email protected]][email protected]], A] // TeXForm
$$2 X.X^T.A$$
If I understood the author's answer in the link mentioned, the package uses the component-wise formula for derivatives that were mentioned in the discussion section above $\frac{\partial M_{\text{ij}}}{\partial M_{\text{kl}}} \bigg\rvert_{M_{\text{rh}},\; \text{r}\neq\text{k},\text{h}\neq\text{l}}=\mathbb{J}_{\text{ijkl}} = \delta _{\text{ik}} \delta _{\text{jl}}$
Using xAct
The author of this answer writes that xAct can only perform matrix derivatives for scalar functions (and in their opinion that is the only way that it makes sense). They explain that if the function is not scalar then one can contract with arbitrary tensors to obtain a scalar function and then identify the derivative.
One of the most popular tensor packages is xAct. In this package, calculations are done using Einstein notation which is basically sums without the sum symbol for our purpose.
The initial expression in terms of Transpose and Tr needs to be translated into a sum. I did this by hand but perhaps there is an easy way to do this automatically.
In the usual sum notation
$$ \text{tr} \left[(\vec X^T A)^T (\vec X^T A)\right]= \sum_{i,j,k}{x_i x_j A_{ik} A_{jk}}$$
In Einstein notation:
$$ \text{tr} \left[(\vec X^T A)^T (\vec X^T A)\right]=x_i x_j A_k^{i} A^{jk}$$
In xAct notation:
$ \text{tr} \left[(\vec X^T A)^T (\vec X^T A)\right]=$ x[-i] x[-j] A[i,k] A[j,-k]
Tensors and the manifold they "live" on needs to be defined in xAct. As the $A$ matrix is not square I defined the manifold for $A$ as a product manifold between $\mathbb{R}^m$ and $\mathbb{R}^n$. After installing xAct from the xact website one can write the evaluate the following definitions :
<< xAct`xTensor`
DefConstantSymbol[n]; DefConstantSymbol[m];
DefManifold[Mn, n, {an, bn, cn, en, fn}];
DefManifold[Mm, d2, {am, bm, cm, em, fm}]
DefManifold[M, {Mn, Mm}, {a, b, c, e, f}];
DefMetric[1, metric[-a, -b], cov];
DefTensor[A[a, b], M]; DefTensor[X[an], Mn];
$PrePrint = ScreenDollarIndices;
The matrix derivative can be found using VarD from the xAct package (the output in Mathematica is formatted to show up and down indices like in Einstein notation):
VarD[A[e, f]][X[-a]*X[-b]*A[a, c]*A[b, -c]] // ContractMetric
out: (* 2 A[a, -f] X[-a] X[-e] *)
