16
$\begingroup$

I'm dealing with derivatives of scalar functions of matrices and wondering if Mathematica can help me here.

The standard approach of expanding it in terms of components is cumbersome. As an motivating example, I want to minimize the following function, where $X$ is a matrix

$$f(X) = \text{tr}(X'X)$$

I can use matrix differential calculus to derive that one step of gradient descent to minimize this function is:

$$X^* = X - 2 X$$

On other hand, suppose $X$ is square, and my function is $$g(X)=\text{tr}(X^2)$$

Now, a single step of gradient descent looks as follows $$X^* = X - 2 X'$$

This can get complicated to do by hand, as an example, Exercise 1 of Magnus 9.10 asks to show that gradient descent step of the following function

$$h(X) = \det A X B$$

is the following

$$X^* = X-\det(AXB)(A'(B'X'A')^{-1}B')'$$

now take $$h^*(X) = \text{tr}(AX'BXC)$$ formula to do a single step of gradient descent is $$ \text{probably something simple} $$

Is there a way to get help deriving/checking expressions above in Mathematica?

(note, I'm using "gradient descent step" instead of "derivative" because there are multiple notations for derivative which differ in shape, but reformulating it as gradient descent removes ambiguity)

$\endgroup$
  • 2
    $\begingroup$ There's a name I haven't seen in a while! Welcome back. :-) $\endgroup$ – Mr.Wizard Feb 26 '17 at 19:57
  • 3
    $\begingroup$ I missed Mathematica, glad to see it's still improving and the community is active! 5 years at Google made me too focused on practical stuff, now I'm at a research non-profit and can explore things again :) $\endgroup$ – Yaroslav Bulatov Feb 26 '17 at 20:01
  • $\begingroup$ I look forward to seeing more posts from you in that case. $\endgroup$ – Mr.Wizard Feb 26 '17 at 20:03
  • $\begingroup$ For the second one, the identity det(A) = exp(Tr(log(A))) might be useful. $\endgroup$ – Leonid Shifrin Feb 26 '17 at 20:33
  • $\begingroup$ Are you assuming X is square? $\endgroup$ – bill s Feb 26 '17 at 21:40
3
$\begingroup$

It is certainly fairly easy to check these relations for specific sizes of array:

X = Array[x, {7, 11}];
Map[D[Tr[Transpose[X].X], #] &, X, {2}] == 2 X
(* True *)

Y = Array[y, {17, 17}];
Map[D[Tr[Y.Y], #] &, Y, {2}] == 2 Transpose[Y]
(* True *)

Generalisation to the determinant expression should be relatively simple...

EDIT

Without being too rigorous about it, we can use the following relations to find some matrix derivatives (I assume in the following that all arguments to Tr are square )

With[{m = 4, n = 2, p = 3},
 A = Array[a, {m, n}];
 X = Array[x, {p, n}];
 B = Array[b, {p, p}];
 F = Array[c, {n, m}];
 Y = Array[y, {n, p}];]


Tr[B] == Tr[Transpose[B]]
(* True *)

Tr[X.Y] == Tr[Y.X]
(* True *)

This allows us to derive (for example)

Tr[A.Transpose[X].B.X.F] == Tr[B.X.F.A.Transpose[X]] == Tr[Transpose[B].X.Transpose[A].Transpose[F].Transpose[X]] // Simplify
(* True *)

and guess the form of the derivative

Map[D[Tr[A.Transpose[X].B.X.F], #] &, X, {2}] == 
  B.X.F.A + Transpose[B].X.Transpose[A].Transpose[F] // Expand
(* True *)
$\endgroup$
  • $\begingroup$ Thanks, that useful. I guess I still have to do the hand-algebra to get the actual value to check against, for instance now I'm looking at getting derivative of tr(AX'BXC) $\endgroup$ – Yaroslav Bulatov Feb 26 '17 at 22:19
11
$\begingroup$

In this post I discuss a function MatrixD which attempts to take a matrix derivative following the guidelines given in the The Matrix Cookbook.

I still want to take advantage of the normal partial derivative function D, but I need to override the default handling of matrix functions. The basic approach is the following:

  1. Add Inverse, Tr, Transpose, Det, MatrixPower, and MatrixFunction to the list of "ExcludedFunctions" in SystemOptions[]. The default handling of derivatives of these symbols needs to be avoided (D[Tr[X], X] evaluating to Tr'[X] is not useful).

  2. The matrix derivative of a matrix is not 1, but rather is given by the following: $$\frac{\partial M_{\text{ij}}}{\partial M_{\text{kl}}}=\mathbb{J}_{\text{ijkl}}$$ where $\mathbb{J}_{\text{ijkl}}\equiv \delta _{\text{ik}} \delta _{\text{jl}}$. Since D[X, X] invariably returns 1, I will instead use the following construct:

    MatrixD /: D[X, MatrixD, NonConstants->{X}] := $SingleEntryMatrix
    

    MatrixD here acts as a variable, and the NonConstants option prevents the derivative from being 0. This DownValue for D (as well as all of the others) are hidden inside the implementation of MatrixD.

  3. The result of MatrixD will need to be simplified further. For this purpose, I wrote a MatrixReduce function. Some of the simplifications performed by MatrixReduce depend on whether the matrix variable is invertible or not. This is controlled by the "Invertible" option of MatrixD and MatrixReduce.

  4. Scalar functions of a matrix (e.g., Log[X]) are not supported. Instead, one must use MatrixFunction (or MatrixPower), e.g., MatrixFunction[Log, X] instead of Log[X].

  5. MatrixFunction has a simple differentiation rule when it is inside of a Tr/Det, but much more complicated rules if it is outside of a Tr/Det. Hence, derivatives of MatrixFunction are only supported when they occur inside of a Tr/Det wrapper.

  6. For test purposes, I also created a TestMatrixD function. This function tests the MatrixD output by performing the matrix derivative using ordinary D rules as well as with the MatrixD, and comparing the output after substituting random real matrices for the variables and matrix constants.

  7. There is at least one bug in the MatrixD implementation, when the argument to one of the matrix functions Tr, Det, Inverse and Transpose is actually a scalar and not a matrix (e.g., Tr[Det[X]]). There may be more.

The code is a little lengthy, so I put it on GitHub. You can download it at:

https://github.com/carlwoll/MatrixD/releases

The source code can be viewed there as well. Download the file "MatrixD-1.0.paclet" and then install it with:

PacletInstall[file]

Load the package with:

<<MatrixD`

Now, for your first couple elementary examples:

MatrixD[Tr[X\[Transpose].X], X]//TeXForm

$2 X$

MatrixD[Tr[X.X], X]//TeXForm

$2 X^T$

Your slightly more complicated Det example:

MatrixD[Det[A.X.B], X, "Invertible"->False]//TeXForm

$\left| A.X.B\right| A^T.\left(B^T.X^T.A^T\right)^{-1}.B^T$

And the output for your last example:

MatrixD[Tr[A.X.B.X\[Transpose].C], X, "Invertible"->False]//TeXForm

$A^T.C^T.X.B^T+C.A.X.B$

I'm planning on adding support for vectors as well as matrices in the future.

$\endgroup$
  • $\begingroup$ Looks neat. You might want to add MatrixExp[] and MatrixLog[] to the list of "ExcludedFunctions" as well. $\endgroup$ – J. M. will be back soon Mar 29 '17 at 8:23
  • $\begingroup$ This is probably a long shot, but: do you think your package can support taking Gâteaux/Fréchet derivatives someday? (These are more complicated, so I'd understand if you think they're out of scope.) $\endgroup$ – J. M. will be back soon Mar 29 '17 at 13:19
  • $\begingroup$ Thanks a lot, that's gives a much easier way to verify calculations! BTW, one potential extension is to allow dealing with vectorized forms because that lets one do Hessians. IE, suppose we are trying to fit Y=W.X using least squares, what's is the Hessian of optimization problem? Since the parameter is a matrix, to see Hessian, the parameter needs to be vectorized. Following Theorem 1 of 10.6 of Magnus, the answer is XX'\otimes I . But this gets more complicated if my function is something like Y=X1X2X3 which happens for pcwise linear nn's $\endgroup$ – Yaroslav Bulatov Mar 29 '17 at 16:04
  • $\begingroup$ BTW, here's an example of computing Hessians -- math.stackexchange.com/questions/2206158/… . There's a pretty nice expression for hessian of J=Tr(e'e) where e=X1.X2....Xn and Xi's are variables, which can be used to minimize J using Newton's method (worked out by hand as hessBlock in here ) The diagonal block with respect to Xi is nice, it's (A.A')\otimes(B.B') where A is product of all matrices before Xi in e and B', all matrices after Xi in e $\endgroup$ – Yaroslav Bulatov Apr 25 '17 at 19:38
  • $\begingroup$ This is not answer but rather question to Carl Woll. Please check out this website: matrixcalculus.org/matrixCalculus Maybe you could extend your package to take care scalar and vector cases as well? What about block matrices? $\endgroup$ – Kun Deng Jun 21 '17 at 3:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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