15
$\begingroup$

Let's say I want to calculate the following scalar-by-matrix derivative

$$\frac{\partial}{\partial A} \text{tr} \left[(\vec X^T A)^T (\vec X^T A)\right],$$

with $\vec X$ and $A$ being a $n \times 1$ and a $n \times m$ matrix, respectively. Is there a way in Mathematica to get the result

$$2 \vec X (\vec X^T A)$$

without explicitly defining (for instance)

n=3
m=2
A=Array[a,{n,m}]
X=Array[x,{n,1}]

and calculating

D[Tr[Transpose[Transpose[X].A].Transpose[X].A],{A}]

? The problem with this approach is that the Mathematica result cannot be easily cast back into an human-readable form like

$$2 \vec X (\vec X^T A).$$

$\endgroup$

3 Answers 3

7
$\begingroup$

Your expression simplifies to this $$\vec X \vec X^T A + A^T \vec X \vec X^T$$ using just these rules

Unprotect[D, Transpose, Dot];
(*Derivative rules*)
D[Tr[A_], X_] := Tr[D[A, X]]
D[Transpose[A_], X_] := D[A, X]\[Transpose]
D[A_ .B_, X_] := D[A, X].B + A.D[B, X]

(*Tranpose rules*)
0\[Transpose] := 0
1\[Transpose] := 1
(A_\[Transpose])\[Transpose] := A
(A_ .B_)\[Transpose] := B\[Transpose].A\[Transpose]

(*Dot rules*)
Dot[_, 0] := 0
Dot[0, _] := 0
Dot[1, A_] := A
Dot[A_, 1] := A
Protect[D, Transpose, Dot];

If it's correct then I'm quite sure these two parts are different (transpose of one another), so you can't add them to get coefficient 2.

$\endgroup$
7
  • 2
    $\begingroup$ Thank you very much for your answer. Are you sure that the first rule is correct? I think that the matrix derivative of the trace (i.e. a scalar) should give a matrix, whereas the trace of the tensor D[A,X] will give a scalar. You would have different objects on both sides of D[Tr[A_], X_] := Tr[D[A, X]]. $\endgroup$
    – cryo111
    Apr 11, 2014 at 12:55
  • 1
    $\begingroup$ @cryo111 Wouldn't a matrix derivative of a matrix be a fourth rank tensor? Therefore its trace should be a matrix, not a scalar. $\endgroup$
    – swish
    Apr 11, 2014 at 18:06
  • 1
    $\begingroup$ The Mathematica help says that "For a higher-rank tensor, Tr gives the sum of elements with equal indices" (i.e., a scalar). $\endgroup$
    – cryo111
    Apr 11, 2014 at 19:25
  • 2
    $\begingroup$ But yes, a matrix-by-matrix derivative should give a fourth-rank tensor. $\endgroup$
    – cryo111
    Apr 11, 2014 at 19:31
  • 2
    $\begingroup$ When I use your approach, I get $\text{tr}(\vec X \vec X^T A + A^T \vec X \vec X^T)$ as a result. Apart from the additional trace (that should not be there), $\vec X \vec X^T A$ is a $n \times m$ matrix, whereas $A^T \vec X \vec X^T$ is a $m \times n$ matrix. Therefore, these cannot be added. The result $2 \vec X \vec X^T A$ is correct and cross-checked with the explicit example given above (n=3;m=2). $\endgroup$
    – cryo111
    Apr 14, 2014 at 14:54
5
$\begingroup$

I tried the NCAlgebra package:

<< NC`;
<< NCAlgebra`;
SNC[X, A]
NCGrad[tr[tp[tp[X] ** A] ** (tp[X] ** A)], A]

And got this result:

2 A^T ** X ** X^T

This is a result similar to that which can be obtained by hand.

$\endgroup$
3
$\begingroup$

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] *)

$\endgroup$
1
  • $\begingroup$ Note that in the first section using differentials, H is hard coded in the TagSetDelayed lines for dot. Thus, the variable H used for obtaining the differential plays a special role in the code. Therefore, the code would not work if H is replaced by another variable. This is not important if no matrices are called H in the notebook. $\endgroup$ Oct 6, 2022 at 4:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.