Denote $F$ a matrix of size A x C, $G$ a matrix of size A x B, $H$ a matrix of size B x C and $F = G H$ . For a real number variable $x$ , it is known $\frac{\partial F}{\partial x} = \frac{\partial G}{\partial x} H + G \frac{\partial H}{\partial x}$. Lay out a collection of variables $x_{i,j}$ by entries of a matrix $X$ of size M x N. Then layout block-wise the $i$ in row and $j$ in column block with respect the variable $x_{i,j}$ also the definition of $\frac{\partial F}{\partial X}$ whose entry blocks are all A x C matrix:

$\frac{\partial F}{\partial X}\equiv\begin{bmatrix}\frac{\partial F}{\partial x_{1,1}}& \frac{\partial F}{\partial x_{1,2}}& \cdots & \frac{\partial F}{\partial x_{1,N}}\\ \frac{\partial F}{\partial x_{2,1}}& \frac{\partial F}{\partial x_{2,2}} & \cdots & \frac{\partial F}{\partial x_{2,N}} \\ \vdots & \vdots & \ddots & \vdots \\\frac{\partial F}{\partial x_{M,1}}& \frac{\partial F}{\partial x_{M,2}}& \cdots &\frac{\partial F}{\partial x_{M,N}}\end{bmatrix}$

Obviously, $\left(\frac{\partial F}{\partial X}\right)^T = \frac{\partial F^T}{\partial X^T}$ . Also by the mentioned identity for single variable, it leads to:

$\begin{bmatrix}\frac{\partial F}{\partial x_{1,1}} & \frac{\partial F}{\partial x_{1,2}} & \cdots & \frac{\partial F}{\partial x_{1,N}} \\ \frac{\partial F}{\partial x_{2,1}} & \frac{\partial F}{\partial x_{2,2}} & \cdots & \frac{\partial F}{\partial x_{2,N}} \\ \vdots &\vdots &\ddots& \vdots \\ \frac{\partial F}{\partial x_{M,1}} & \frac{\partial F}{\partial x_{M,2}} & \cdots & \frac{\partial F}{\partial x_{M,N}}\end{bmatrix} =\begin{bmatrix}\frac{\partial G}{\partial x_{1,1}} H + G \frac{\partial H}{\partial x_{1,1}} & \frac{\partial G}{\partial x_{1,2}} H + G \frac{\partial H}{\partial x_{1,2}} & \cdots &\frac{partial G}{\partial x_{1,N}} H + G \frac{\partial H}{\partial x_{1,N}} \\ \frac{\partial G}{\partial x_{2,1}} H + G \frac{\partial H}{\partial x_{2,1}} & \frac{\partial G}{\partial x_{2,2}} H + G \frac{\partial H}{\partial x_{2,2}} & \cdots & \frac{\partial G}{\partial x_{2,N}} H + G \frac{\partial H}{\partial x_{2,N}} \\ \vdots &\vdots &\ddots& \vdots \\\frac{\partial G}{\partial x_{M,1}} H + G \frac{\partial H}{\partial x_{M,1}} &\frac{\partial G}{\partial x_{M,2}} H + G \frac{\partial H}{\partial x_{M,2}} & \cdots & \frac{\partial G}{\partial x_{M,N}} H + G \frac{\partial H}{\partial x_{M,N}}\end{bmatrix}$

$=\begin{bmatrix}\frac{\partial G}{\partial x_{1,1}} H & \frac{\partial G}{\partial x_{1,2}} H & \cdots & \frac{\partial G}{\partial x_{1,N}} H \\ \frac{\partial G}{\partial x_{2,1}} H & \frac{\partial G}{\partial x_{2,2}} H & \cdots & \frac{\partial G}{\partial x_{2,N}} H \\\vdots &\vdots& \ddots& \vdots\\ \frac{\partial G}{\partial x_{M,1}} H & \frac{\partial G}{\partial x_{M,2}} H & \cdots & \frac{\partial G}{\partial x_{M,N}} H\end{bmatrix} + \begin{bmatrix}G \frac{\partial H}{\partial x_{1,1}} & G \frac{\partial H}{\partial x_{1,2}} & \cdots & G \frac{\partial H}{\partial x_{1,N}} \\ G \frac{\partial H}{\partial x_{2,1}} & G \frac{\partial H}{\partial x_{2,2}}& \cdots &G \frac{\partial H}{\partial x_{2,N}}\\ \vdots &\vdots &\ddots& \vdots\\ G \frac{\partial H}{\partial x_{M,1}} & G \frac{\partial H}{\partial x_{M,2}} & \cdots & G \frac{\partial H}{\partial x_{M,N}}\end{bmatrix}$

$= \begin{bmatrix}\frac{\partial G}{\partial x_{1,1}} & \frac{\partial G}{\partial x_{1,2}} &\cdots& \frac{\partial G}{\partial x_{1,N}} \\ \frac{\partial G}{\partial x_{2,1}} & \frac{\partial G}{\partial x_{2,2}}& \cdots& \frac{\partial G}{\partial x_{2,N}}\\ \vdots &\vdots& \ddots& \vdots\\ \frac{\partial G}{\partial x_{M,1}} & \frac{\partial G}{\partial x_{M,2}} &\cdots& \frac{\partial G}{\partial x_{M,N}}\end{bmatrix} \begin{bmatrix}H&&&\\&H&&\\&&\ddots&\\&&&H\end{bmatrix} + \begin{bmatrix}G&&&\\&G&&\\&&\ddots&\\&&&G\end{bmatrix} \begin{bmatrix}\frac{\partial H}{\partial x_{1,1}} & \frac{\partial H}{\partial x_{1,2}} &\cdots& \frac{\partial H}{\partial x_{1,N}} \\ \frac{\partial H}{\partial x_{2,1}} & \frac{\partial H}{\partial x_{2,2}} & \cdots & \frac{\partial H}{\partial x_{2,N}}\\ \vdots& \vdots& \ddots& \vdots\\ \frac{\partial H}{\partial x_{M,1}}& \frac{\partial H}{\partial x_{M,2}} & \cdots & \frac{\partial H}{\partial x_{M,N}}\end{bmatrix}$

Given an integer $a$ and a matrix $S$ , define $I_{a,S} \equiv \begin{bmatrix}S&&&\\&S&&\\&&\ddots&\\&&&S\end{bmatrix}$ which has $a$ number of $S$ along the diagonal blocks and all other entry blocks are zero,

Then $\frac{\partial F}{\partial X} = \frac{\partial G}{\partial X} I_{N,H} + I_{M,G} \frac{\partial H}{\partial X}$

$I_{a,S}$ is the extension of the concept of scalar product. Suppose a matrix of size M x N, right-multiply a scalar $s$ is actually right-multiply the matrix $I_{N,s}$ and left-multiply a scalar $s$ is actually left-multiply the matrix $I_{M,s}$ . Assuming matrix are all compatible size below, some facts about $I_{a,S}$ :

$\begin{aligned}I_{1,S}=S\\(I_{a,S})^T=(I_{a,S})^T\\I_{a,S_1 S_2}=I_{a,S_1}I_{a,S_2}\\I_{a,S_1 + S_2}=I_{a,S_1}+I_{a,S_2}\\I_{M,1}=I_M\\I_{M,I_N}=I_{M N}\end{aligned}$

Let $X \equiv \begin{bmatrix}x_1\\ \vdots \\x_M\end{bmatrix}$, $Y \equiv \begin{bmatrix}y_1&\cdots&y_M\end{bmatrix}$, $e_i$ be the M x 1 matrix whose $i$-th row is 1 and other entries are 0, $w_i$ be the 1 x M matrix whose $i$-th column is 1 and other entries are 0.

$\begin{aligned}\frac{\partial X^T}{\partial X} = I_M\\\frac{\partial Y^T}{\partial Y} = I_M\\\frac{\partial X}{\partial X} =\begin{bmatrix}e_1\\e_2\\\vdots\\e_M\end{bmatrix}\\\frac{\partial Y}{\partial Y} = \begin{bmatrix}w_1&\cdots&w_M\end{bmatrix}\\I_{M,X^T} \begin{bmatrix}e_1\\e_2\\\vdots\\e_M\end{bmatrix}= X\\\begin{bmatrix}w_1&\cdots&w_M\end{bmatrix} I_{M,X} =X^T\end{aligned}$

Example

$G = A X$ where $A$ is a constant 4 x 3 matrix and $X$ is 3 x 1 matrix of variables and therefore $G^T G$ is a real number function of $X$ and denoted by $L(X)$ . Then the gradient of $L(X)$ is defined as a 3 x 1 matrix $\nabla L(X) \equiv \frac{\partial L}{\partial X}$ (while it is somewhere defined as the 1 x 3 matrix $\frac{\partial L}{\partial X^T}$)

$\begin{aligned}\nabla L(X) = \frac{\partial G^T G}{\partial X} = \frac{\partial X^T A^T A X}{\partial X} = \frac{\partial X^T A^T A}{\partial X} I_{1,X} + I_{3,X^T A^T A} \frac{\partial X}{\partial X} \\= \left(\frac{\partial X^T}{\partial X} I_{1,A^T A} + 0 \right) X + I_{3,X^T A^T A} \frac{\partial X}{\partial X} = A^T A X + I_{3,X^T A^T A} \begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix}\end{aligned}$

Let $A$'s 3 column vectors be $C_1$ and $C_2$ and $C_3$ therefore $A = \begin{bmatrix}C_1&C_2&C_3\end{bmatrix}$ then $\begin{aligned}A^T A X + I_{3,X^T A^T A} \begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix} = \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X + I_{3,X^T A^T} I_{3,A} \begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix} \\= \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X + \begin{bmatrix}X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_1\\X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_2\\ X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_3\end{bmatrix}\end{aligned}$

Because   $X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_i$ is a number, it is the same as its transpose $C_i^T \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X$

therefore

$\begin{bmatrix}X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_1\\X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_2\\ X^T \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} C_3\end{bmatrix}=\begin{bmatrix}C_1^T \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X\\C_2^T \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X\\C_3^T \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X\end{bmatrix}= \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X$

Therefore the answer is $2 \begin{bmatrix}C_1^T\\C_2^T\\C_3^T\end{bmatrix} \begin{bmatrix}C_1&C_2&C_3\end{bmatrix} X = 2 A^T A X$

When the column vectors of $A$ are orthonormal, $A^T A = I_3$ ,so $A^T A X + I_{3,X^T A^T A} \begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix}$ becomes $X + I_{3,X^T} \begin{bmatrix}e_1\\e_2\\e_3\end{bmatrix} = X + X = 2X$ . Actually $L$ is $x_1^2 + x_2^2 + x_3^2$ and $\frac{\partial L}{\partial X} = \begin{bmatrix}2 x_1\\2 x_2\\2 x_3\end{bmatrix}$ followed by direct calculation of the definition.

Multiple variables integration

Let $A(Y)$ be a function of $R^M$ to $R$, $Y=F(X)$ be a $R^M$ to $R^M$ change of variables. Layout as

$\begin{bmatrix}y_1\\\vdots\\y_M\end{bmatrix} \equiv Y = F(X) \equiv \begin{bmatrix}F_1(X)\\\vdots\\F_M(X)\end{bmatrix}$ and $X \equiv \begin{bmatrix}x_1&\cdots&x_M\end{bmatrix}$

Then $\int A(y_1,\cdots,y_M) d y_1\cdots d y_M = \int A(F_1(X),\cdots,F_M(X) det(\frac{\partial F}{\partial X})d x_1 \cdots d x_M$

which is the typical change of variable of integration of one variable when $M$ is 1

Example. Calculate $\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{- \frac{y_1^2 + y_2^2}{2}} d y_1 d y_2$

Let $\begin{bmatrix}y_1\\y_2\end{bmatrix} \equiv Y = F(X) \equiv \begin{bmatrix}x_1 \cos(x_2)\\x_1 \sin(x_2)\end{bmatrix}$

and $X \equiv \begin{bmatrix}x_1&x_2\end{bmatrix}$ . Then $\frac{\partial F}{\partial X} =\begin{bmatrix}\cos(x_2)&- x_1 \sin(x_2)\\\sin(x_2)&x_1 \cos(x_2)\end{bmatrix}$ and therefore $det(\frac{\partial F}{\partial X}) = x_1$

$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{- \frac{y_1^2 + y_2^2}{2}} d y_1 d y_2 = \int_0^{2 \pi} \int_0^{\infty} e^-{\frac{x_1^2}{2}} x_1 d x_1 d x_2 = 2\pi$

As a consequence, $\int_{-\infty}^{\infty} e^{- \frac{Z^2}{2}} d Z = \sqrt{2 \pi}$ aka the density of standard normal random variable is $\frac{1}{\sqrt{2 \pi}} e^{- \frac{Z^2}{2}}$

Let $f$ be the probability density of $N$ random variables with zero mean. Layout $X = \begin{bmatrix}X_1\\\vdots\\X_N\end{bmatrix}$. Then its covariance matrix is $\int_{-\infty}^\infty X X^T f(X) d X_1 d X_2 \cdots d X_N$ aka $E(X X^T)$ aka $COV_X$ . Let $Y = R X$ where $R$ is a $N$ x $N$ matrix. $COV_Y = E(Y Y^T) = E(R X X^T R^T) = R E(X X^T) R^T = R COV_X R^T$. So if with a series of $k$ row operations $R = R_k \cdots R_1$ as well as the correspondent column operations $R^T$ on $COV_X$ leading to $COV_Y = I$ , with $A \equiv R ^{- 1} = R_1^{-1} \cdots R_k^{-1}$ , then a change of variables of $X$ defined as $X = A Y$ will have covariance matrix $COV_X = A A^T$ because $I = R COV_X R^T$ . Any symmetric matrix can be operated with this row-column operations, once a diagonal entry is not positive, this symmetric matrix fails to be a legit covariance matrix. Also $det(COV_X) = det(A) det(A^T) = det(A)^2$ so $det(A)= \sqrt{det(COV_X)}$

Demonstrate the row-column operation on $COV = \begin{bmatrix}2&4&-2\\4&10&2\\-2&2&40\end{bmatrix}$ . $COV vs A$:

$\begin{aligned}\begin{bmatrix}2&4&-2\\4&10&2\\-2&2&40\end{bmatrix} vs \begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix}\\\begin{bmatrix}2&0&0\\0&2&6\\0&6&38\end{bmatrix} vs \begin{bmatrix}1&0&0\\2&1&0\\-1&0&1\end{bmatrix}\\\begin{bmatrix}2&0&0\\0&2&0\\0&0&20\end{bmatrix} vs \begin{bmatrix}1&0&0\\2&1&0\\-1&3&1\end{bmatrix}\\\begin{bmatrix}1&0&0\\0&2&0\\0&0&20\end{bmatrix} vs \begin{bmatrix}\sqrt{2}&0&0\\2\sqrt{2}&1&0\\-\sqrt{2}&3&1\end{bmatrix}\\\begin{bmatrix}1&0&0\\0&1&0\\0&0&20\end{bmatrix} vs \begin{bmatrix}\sqrt{2}&0&0\\2 \sqrt{2}&\sqrt{2}&0\\-\sqrt{2}&3\sqrt{2}&1\end{bmatrix}\\\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix} vs \begin{bmatrix}\sqrt{2}&0&0\\2 \sqrt{2}&\sqrt{2}&0\\-\sqrt{2}&3\sqrt{2}&\sqrt{20}\end{bmatrix}\end{aligned}$

So the $A= \begin{bmatrix}\sqrt{2}&0&0\\2 \sqrt{2}&\sqrt{2}&0\\-\sqrt{2}&3\sqrt{2}&\sqrt{20}\end{bmatrix}$ and $A A^T = COV$:

$\begin{bmatrix}\sqrt{2}&0&0\\2 \sqrt{2}&\sqrt{2}&0\\-\sqrt{2}&3\sqrt{2}&\sqrt{20}\end{bmatrix} \begin{bmatrix}\sqrt{2}&2\sqrt{2}&-\sqrt{2}\\0&\sqrt{2}&3\sqrt{2}\\0&0&\sqrt{20}\end{bmatrix}= \begin{bmatrix}2&4&-2\\4&10&2\\-2&2&40\end{bmatrix}$

Example. Calculate $\int_{-\infty}^\infty e^{-\frac{1}{2}X^T (COV_X)^{-1} X} d X_1 d X_2\cdots d X_N$

Let $A$ be a matrix such that $COV_X = A A^T$ by the above procedures. Then a change of variables leads to

$\begin{aligned}\int_{-\infty}^\infty e^{-\frac{1}{2}X^T (COV_X)^{-1} X} d X_1 d X_2\cdots d X_N = \int_{-\infty}^\infty e^{-\frac{1}{2} Y^T A^T (A A^T)^{-1} A Y} det(A) d Y_1 d Y_2 \cdots d Y_N \\= \int_{-\infty}^\infty e^{-\frac{1}{2} Y^T Y} det(A) d Y_1 d Y_2 \cdots d Y_N = \sqrt{(2 \pi)^N det(COV_X)}\end{aligned}$

Meaning, let $Y_1 \cdots Y_N$ be iid standard normal distribution, then random variables $X = A Y$ will have $COV_X = A A^T$ and its density is $\frac{1}{\sqrt{(2 \pi)^N det(COV_X)}} e^{-\frac{1}{2} X^T (COV_X)^{-1} X}$ where $A$ can be found by the procedure mentioned above.