Matrices and Derivatives: Solving Some Problems

Written on January 29, 2019

NOTE: This post is automatically translated from Farsi using GPT 5.5

The main use of differentiation is finding extremum points; but the way derivatives are used to reach an extremum can be either through iterative optimization methods or by formally solving equations. To see an example of each of these cases, first we estimate the value of the correlation matrix from a number of samples drawn from a Gaussian distribution, and then we compute the gradient of an MLP with respect to its weights.

Estimating the correlation matrix

The formula for the multivariate Gaussian distribution is what appears below:

$$p(x|\mu ,\mathrm{\Sigma })=\text{det}(2\pi \mathrm{\Sigma }{)}^{-\frac{1}{2}}\text{exp}(-\frac{1}{2}(x-\mu {)}^T{\mathrm{\Sigma }}^{-1}(x-\mu ))$$

Well now, if we have samples $x_1$ through $x_n$ from this distribution, what values of $\mu$ and $\mathrm{\Sigma }$ maximize the probability of these samples occurring?

Assuming the samples are independent of one another, the probability of all of them occurring becomes:

$$\begin{array}{cc}& p(x_1,x_2,\dots x_n|\mu ,\mathrm{\Sigma })=\prod p(x_i|\mu ,\mathrm{\Sigma })\\ & \text{log}(p(x_1,x_2,\dots x_n|\mu ,\mathrm{\Sigma }))=\sum \text{log}(p(x_i|\mu ,\mathrm{\Sigma }))\end{array}$$

When we can do with the log-probability the same thing we can do with the probability, why not go for it? Especially for the Gaussian distribution, where the exponential function plays a fundamental role. Well, by substituting the probability distribution relation into the logarithm formulas, we have:

$$\text{log}(p(x_1,x_2,\dots x_n|\mu ,\mathrm{\Sigma }))=\sum [\text{log}(\text{det}(2\pi \mathrm{\Sigma }{)}^{-\frac{1}{2}})-\frac{1}{2}(x_i-\mu {)}^T{\mathrm{\Sigma }}^{-1}(x_i-\mu )]$$

If we simplify it even further, it becomes:

$$\begin{array}{cc}y& =\text{log}(p(x_1,x_2,\dots x_n|\mu ,\mathrm{\Sigma }))\\ & =-\frac{n}{2}\text{log}\left(2^m{\pi }^m\right)-\frac{n}{2}\text{log}(\text{det}(\mathrm{\Sigma }))-\frac{1}{2}\sum (x_i-\mu {)}^T{\mathrm{\Sigma }}^{-1}(x_i-\mu )\end{array}$$

For this simplification, I first took the determinant expression, which had nothing to do with the summation index, out of the summation, which naturally means it gets multiplied by the number of times the sum is taken (that is, $n$). Then I pulled $-\frac{1}{2}$ out of the logarithm, where it appears as a multiplier. In the next step, I simplified $\text{det}(2\pi \mathrm{\Sigma })$ using the relation $\text{det}(kA)=k^n\text{det}(A)$. Finally, I moved the coefficient $-\frac{1}{2}$ from inside the summation to the outside, which is possible because of the distributive property of multiplication (have you noticed that $n$ is the number of samples and $m$ is the number of elements in each sample?)

Now all we need to do is take two derivatives and set them equal to zero. One with respect to $\mu$ (a vector) and the other with respect to $\mathrm{\Sigma }$ (a matrix). We start with the vector first:

$$\begin{array}{cc}& \frac{\partial y}{\partial \mu }=-\frac{1}{2}({\mathrm{\Sigma }}^{-1}+{\mathrm{\Sigma }}^{-T})((\sum x_i)-n\mu )=0\\ & \Rightarrow \mu =\frac{1}{n}\sum x_i\end{array}$$

It was computed that quickly! Do you know why I paid no attention to ${\mathrm{\Sigma }}^{-1}$ and its counterpart ${\mathrm{\Sigma }}^{-T}$? Because according to the formula, these matrices are invertible, and their product with a vector becomes zero in only one case, namely when that vector itself is zero; that is, $(\sum x_i)-n\mu =0$. But the second derivative, which is also more interesting:

$$\begin{array}{cc}& \frac{\partial y}{\partial \mathrm{\Sigma }}=\frac{\partial y}{\partial {\mathrm{\Sigma }}^{-1}}\frac{\partial {\mathrm{\Sigma }}^{-1}}{\partial \mathrm{\Sigma }}\\ & \frac{\partial y}{\partial {\mathrm{\Sigma }}^{-1}}=-\frac{n}{2}{\text{vec}}^T({\mathrm{\Sigma }}^T)-\frac{1}{2}\sum ((x_i-\mu {)}^T\otimes (x_i-\mu {)}^T)\\ & \frac{\partial {\mathrm{\Sigma }}^{-1}}{\partial \mathrm{\Sigma }}=-({\mathrm{\Sigma }}^{-T}\otimes {\mathrm{\Sigma }}^{-1})\\ & \Downarrow \\ & \frac{\partial y}{\partial \mathrm{\Sigma }}=\begin{array}{cc}& \frac{n}{2}{\text{vec}}^T({\mathrm{\Sigma }}^T)({\mathrm{\Sigma }}^{-T}\otimes {\mathrm{\Sigma }}^{-1})\\ & +\frac{1}{2}\sum ((x_i-\mu {)}^T\otimes (x_i-\mu {)}^T)({\mathrm{\Sigma }}^{-T}\otimes {\mathrm{\Sigma }}^{-1})\end{array}\end{array}$$

Do not be frightened by the lengthiness of the formula above; there is nothing inside it. It is enough to use the following relation for simplification:

$$\text{vec}(AXB)=(B^T\otimes A)\text{vec}(X)$$

Or its equivalent:

$${\text{vec}}^T(AXB)={\text{vec}}^T(X)(B\otimes A^T)$$

By substitution, we arrive at the following result:

$$\frac{\partial y}{\partial \mathrm{\Sigma }}=\frac{n}{2}{\text{vec}}^T({\mathrm{\Sigma }}^{-T})+\frac{1}{2}[\sum ((x_i-\mu {)}^T\otimes (x_i-\mu {)}^T)]({\mathrm{\Sigma }}^{-T}\otimes {\mathrm{\Sigma }}^{-1})=0$$

Or, in other words:

$${\text{vec}}^T({\mathrm{\Sigma }}^{-T})({\mathrm{\Sigma }}^{-T}\otimes {\mathrm{\Sigma }}^{-1}{)}^{-1}=\frac{1}{n}\sum ((x_i-\mu {)}^T\otimes (x_i-\mu {)}^T)=0$$

Again, the same relation $\text{vec}(AXB)$. In the first post, as I said, this relation is very useful in differentiation:

$${\text{vec}}^T({\mathrm{\Sigma }}^{-T})({\mathrm{\Sigma }}^{-T}\otimes {\mathrm{\Sigma }}^{-1}{)}^{-1}={\text{vec}}^T({\mathrm{\Sigma }}^{-T})({\mathrm{\Sigma }}^T\otimes \mathrm{\Sigma })={\text{vec}}^T({\mathrm{\Sigma }}^T{\mathrm{\Sigma }}^{-T}{\mathrm{\Sigma }}^T)={\text{vec}}^T({\mathrm{\Sigma }}^T)$$

Eureka, eureka:

$${\text{vec}}^T({\mathrm{\Sigma }}^T)=\frac{1}{n}\sum ((x_i-\mu {)}^T\otimes (x_i-\mu {)}^T)$$

Again, $\text{vec}(AXB)$! If $X$ is a scalar number, $A$ is a row vector ($n\times 1$), and $B$ is a column vector, what relation is obtained?

$$\begin{array}{cc}& \text{vec}(\alpha ab)=\text{vec}(a\alpha b)=(b^T\otimes a)\text{vec}(\alpha )=\alpha (b^T\otimes a)\\ & {\text{vec}}^T(\alpha ab)={\text{vec}}^T(a\alpha b)={\text{vec}}^T(\alpha )(b\otimes a^T)=\alpha (b\otimes a^T)\end{array}$$

Isn’t it beautiful? Accordingly:

$$\begin{array}{cc}& (x_i-\mu {)}^T\otimes (x_i-\mu {)}^T={\text{vec}}^T((x_i-\mu )(x_i-\mu {)}^T)\\ & {\text{vec}}^T({\mathrm{\Sigma }}^T)=\frac{1}{n}\sum {\text{vec}}^T((x_i-\mu )(x_i-\mu {)}^T)\\ & {\text{vec}}^T({\mathrm{\Sigma }}^T)={\text{vec}}^T(\frac{1}{n}\sum ((x_i-\mu )(x_i-\mu {)}^T))\\ & \Rightarrow \mathrm{\Sigma }={\mathrm{\Sigma }}^T=\frac{1}{n}\sum ((x_i-\mu )(x_i-\mu {)}^T)\end{array}$$

Here you should ask why $\mathrm{\Sigma }={\mathrm{\Sigma }}^T$? The answer is that when ${\mathrm{\Sigma }}^T=\frac{1}{n}\sum ((x_i-\mu )(x_i-\mu {)}^T)$, its transpose is equal to itself.

Derivative with respect to the weights of an MLP

The layers of an MLP network can be represented as follows:

$$y_{n+1}=\varphi (W_n{y}_n+b_n)$$

Depending on how many layers our MLP has, when we want to train it, we compare the output of the $m$-th layer with the desired results and compute the error. Suppose we consider the least mean squares error criterion:

$$j(W_0,b_0,W_1,b_1,\dots W_{m-1},b_{m-1})=\frac{1}{2}\Vert y_m-y{\Vert }^2=\frac{1}{2}(y_m-y{)}^T(y_m-y)$$

We use the machinery we built for differentiation to compute $\frac{\partial j}{\partial W_i}$. Computing the derivative with respect to the biases is completely similar:

$$\frac{\partial j}{\partial W_i}=\frac{\partial j}{\partial y_m}\left(\prod _{m-1}^i\frac{\partial y_{j+1}}{\partial y_j}\right)\frac{\partial y_i}{\partial W_i}$$

Now it is enough to compute each component.

$$\begin{array}{cc}& \frac{\partial j}{\partial y_m}=(y_m-y{)}^T\\ & \frac{\partial y_{j+1}}{\partial y_j}=\text{diag}\left(\text{vec}({\varphi }^{\prime }(W_j{y}_j+b_j))\right)W_j\\ & \frac{\partial y_i}{\partial W_i}=(y_i^T\otimes I)\end{array}$$

I leave it to you to figure out how I have computed these using the material I presented in the previous three posts (this is a good exercise in itself.) But the discussion becomes sweeter from here on. In order to continue, I need to introduce a new relation.

$$a\odot b=\text{diag}(a)b\Rightarrow (a\odot b{)}^T=b^T\text{diag}(a)$$

The symbol $\odot$ represents Hadamard, or element-wise, multiplication. On the other hand, when $W_j{y}_j+b_j$ is a vector, if ${\varphi }^{\prime }$ is also applied to it, the output is still a vector, so the whole derivative above can be written like this:

$$\frac{\partial j}{\partial W_i}=(y_m-y{)}^T(\prod \text{diag}({\varphi }^{\prime }(W_j{y}_j+b_j))W_j)(y_i^T\otimes I)$$

It is not a bad idea to see what that $\prod$ in the middle actually does to the error value ($y_m-y$). We start with the error itself:

$$e_j^T=(y_m-y{)}^T\text{diag}({\varphi }^{\prime }(W_j{y}_j+b_j))W_j={((y_m-y)\odot {\varphi }^{\prime }(W_j{y}_j+b_j))}^T{W}_j$$

It is as if each element of the error is weighted by the value of the derivative of the function $\varphi$ at that same element, and then the resulting vector is multiplied from the left by the weight matrix of the previous layer. This is exactly backpropagation of error. Each element of the error is connected, through the corresponding weight, to one of the nodes in the lower layer. The same operation can be continued:

$$e_{j-1}^T=e_j^T\text{diag}({\varphi }^{\prime }(W_{j-1}{y}_{j-1}+b_{j-1}))W_{j-1}={(e_j\odot {\varphi }^{\prime }(W_j{y}_j+b_j))}^T{W}_{j-1}$$

Again, the error value of the upper layer is weighted by the derivative of $\varphi$ in the lower layer and propagated to the layer below. If we proceed in the same fashion, we reach the error of the desired layer, namely $e_i$:

$$\frac{\partial j}{\partial W_i}=e_i^T(y_i\otimes I)={\text{vec}}^T(e_i{y}_i^T)$$

It is a good exercise to see why the relation above holds. Examining what form $\text{vec}(e_i)$ has in this case is the starting point. Anyway.

Ultimately, what is this derivative value used for? We use it in optimization and in moving toward the trained value of the weights. In gradient-based optimization methods, the gradient of the error with respect to the parameters is always added with a negative ratio. This is the basis of all gradient-based optimization methods. But here, the gradient with respect to $W_i$ does not have dimensions similar to $W_i$; how should we add the two?

If we had gone to tensors from the beginning, this problem would not have arisen. The derivative value would also have dimensions equal to $W_i$ itself (of course because $j$ is a scalar number.) Even now, nothing special has happened. Do you remember that in order to define this kind of differentiation, which we have reviewed together over these four posts, we agreed from the beginning to represent every matrix with an ordered pair $(\text{vec}(A),(m,n))$. Now we have the derivative, and we also have the dimensions of the matrices $W_i$. It is enough to re-represent the computed derivative with the dimensions of the original matrix. The number of elements of the derivative is, by definition, the same.

The interesting point in the story is that up to now, as I write this post, without exception, rearranging the derivative for me has been equivalent to looking inside the $\text{vec}$ operator. Look at this very MLP example:

$$\frac{\partial j}{\partial W_i}=e_i^T(y_i\otimes I)={\text{vec}}^T(e_i{y}_i^T)$$

The value $e_i{y}_i^T$ has dimensions exactly equal to the original matrix!

Moral of the story

In the end, I offer these results as a summary:

1-There is no need for tensor computations to take a derivative with respect to a matrix 2-Taking derivatives is good mental exercise 3-In principle, most of the time you do not need to take derivatives with respect to matrices yourself!