What justifies this calculation of the derivative of a matrix function? In Andrew Ng's machine learning course, he uses this formula:
$\nabla_A tr(ABA^TC) = CAB + C^TAB^T$
and he does a quick proof which is shown below:
$\nabla_A tr(ABA^TC) \\
= \nabla_A tr(f(A)A^TC) \\
= \nabla_{\circ} tr(f(\circ)A^TC) + \nabla_{\circ}tr(f(A)\circ^T C)\\
=(A^TC)^Tf'(\circ) + (\nabla_{\circ^T}tr(f(A)\circ^T C)^T \\ 
= C^TAB^T + (\nabla_{\circ^T}tr(\circ^T)Cf(A))^T \\
=C^TAB^T + ((Cf(A))^T)^T \\
= C^TAB^T + CAB$
The proof seems very dense without any comments and I'm having trouble understanding it. What exactly happened from second to third equality? 
 A: There is a subtle but heavy abuse of the notation that renders many of the steps confusing.  Let's address this issue by going back to the definitions of matrix multiplication, transposition, traces, and derivatives.  For those wishing to omit the explanations, just jump to the last section "Putting It All Together" to see how short and simple a rigorous demonstration can be.

Notation and Concepts
Dimensions
For the expression $ABA^\prime C$ to make sense when $A$ is an $m\times n$ matrix, $B$ must be a (square) $n\times n$ matrix and $C$ must be an $m\times p$ matrix, whence the product is an $m\times p$ matrix.  In order to take the trace (which is the sum of diagonal elements, $\operatorname{Tr}(X)=\sum_i X_{ii}$), then $p=m$, making $C$ a square matrix.
Derivatives
The notation "$\nabla_A$" appears to refer to the derivative of an expression with respect to $A$.  Ordinarily, differentiation is an operation performed on functions $f:\mathbb{R}^N\to\mathbb{R}^M$.  The derivative at a point $x\in \mathbb{R}^N$ is a linear transformation $Df(x):\mathbb{R}^N\to\mathbb{R}^M$.  Upon choosing bases for these vector spaces, such a transformation can be represented as an $M\times N$ matrix.  That is not the case here!
Matrices as vectors
Instead, $A$ is being considered as an element of $\mathbb{R}^{mn}$: its coefficients are being unrolled (usually either row by row or column by column) into a vector of length $N=mn$.  The function $f(A)=\operatorname{Tr}(ABA^\prime C)$ has real values, whence $M=1$.  Consequently, $Df(x)$ must be a $1\times mn$ matrix: it's a row vector representing a linear form on $\mathbb{R}^{mn}$.  Howver, the calculations in the question use a different way of representing linear forms: their coefficients are rolled back up into $m\times n$ matrices.
The trace as a linear form
Let $\omega$ be a constant $m\times n$ matrix.  Then, by definition of the trace and of matrix multiplication,
$$\eqalign{
\operatorname{Tr}(A\omega^\prime) &= \sum_{i=1}^m(A\omega^\prime)_{ii} = \sum_{i=1}^m\left(\sum_{j=1}^n A_{ij}(\omega^\prime)_{ji}\right) = \sum_{i,j} \omega_{ij}A_{ij}
}$$
This expresses the most general possible linear combination of the coefficients of $A$: $\omega$ is a matrix of the same shape as $A$ and its coefficient in row $i$ and column $j$ is the coefficient of $A_{ij}$ in the linear combination.  Because $\omega_{ij}A_{ij}=A_{ij}\omega_{ij}$, the roles of $\omega$ and $A$ may switched, giving the equivalent expression
$$\sum_{i,j} \omega_{ij}A_{ij} = \operatorname{Tr}(A\omega^\prime) = \operatorname{Tr}(\omega A^\prime).\tag{1}$$
By identifying a constant matrix $\omega$ with either of the functions $A\to \operatorname{Tr}(A \omega^\prime)$ or $A\to \operatorname{Tr}(\omega A^\prime)$, we may represent linear forms on the space of $m\times n$ matrices as $m\times n$ matrices.  (Do not confuse these with derivatives of functions from $\mathbb{R}^n$ to $\mathbb{R}^m$!)

Computing a Derivative
The definition
Derivatives of many of the matrix functions encountered in statistics are most easily and reliably computed from the definition: you don't really need to resort to complicated rules of matrix differentiation.  This definition says that $f$ is differentiable at $x$ if and only if there is a linear transformation $L$ such that
$$f(x+h) - f(x) = Lh + o(|h|)$$
for arbitrarily small displacements $h\in \mathbb{R}^N$.  The little-oh notation means that the error made in approximating the difference $f(x+h)-f(x)$ by $Lh$ is arbitrarily smaller than the size of $h$ for sufficiently small $h$.  In particular, we may always ignore errors that are proportional to $|h|^2$.
The calculation
Let's apply the definition to the function in question.  Multiplying, expanding, and ignoring the term with a product of two $h$'s in it,
$$\eqalign{
f(A+h)-f(A) &= \operatorname{Tr}((A+h)B(A+h)^\prime C) - \operatorname{Tr}(ABA^\prime C) \\
&= \operatorname{Tr}(hBA^\prime C) +\operatorname{Tr}(ABh^\prime C) + o(|h|).\tag{2}
}$$
To identify the derivative $L=Df(A)$, we must get this into the form $(1)$.  The first term on the right is already in this form, with $\omega^\prime = BA^\prime C$.  The other term on the right has the form $\operatorname{Tr}(Xh^\prime C)$ for $X=AB$.  Let's write this out:
$$\operatorname{Tr}(Xh^\prime C) = \sum_{i=1}^m\sum_{j=1}^n\sum_{k=1}^m  X_{ij} h_{kj} C_{ki} = \sum_{i,j,k}h_{kj} \left(C_{ki}X_{ij}\right) =\operatorname{Tr}((CX)h^\prime).\tag{3}$$
Recalling $X=AB$, $(2)$ can be rewritten
$$f(A+h) - f(A) = \operatorname{Tr}(h\, BA^\prime C\,) + \operatorname{Tr}(CAB\, h^\prime\,)+o(|h|).$$
It is in this sense that we may consider the derivative of $f$ at $A$ to be $$Df(A) = (BA^\prime C)^\prime + CAB = C^\prime A B^\prime + CAB,$$ because these matrices play the roles of $\omega$ in the trace formulas $(1)$.

Putting It All Together
Here, then, is a complete solution.

Let $A$ be an $m\times n$ matrix, $B$ an $n\times n$ matrix, and $C$ an $m\times m$ matrix.  Let $f(A) = \operatorname{Tr}(ABA^\prime C)$.  Let $h$ be an $m\times n$ matrix with arbitrarily small coefficients. Because (by identity $(3)$) $$\eqalign{f(A+h) - f(A) &= \operatorname{Tr}(hBA^\prime C) +\operatorname{Tr}(ABh^\prime C) + o(|h|) \\
&=\operatorname{Tr}(h(C^\prime A B^\prime)^\prime + (CAB)h^\prime) + o(|h|),}$$ $f$ is differentiable and its derivative is the linear form determined by the matrix $$C^\prime A B^\prime + CAB.$$

Because this takes only about half the work and involves only the most basic manipulations of matrices and traces (multiplication and transposition), it has to be considered a simpler--and arguably more perspicuous--demonstration of the result.  If you really want to understand the individual steps in the original demonstration, you might find it fruitful to compare them to the calculations shown here.
A: There are some special matrix notations in machine learning community, which are implicit by experts but often confuse new practitioners. Some of such notations are described briefly in page 698 of "Pattern Recognition and Machine Learning" by Christopher M. Bishop. So, your life will be easier if you are familiar with that page.
If we have a scalar function (e.g., trace) of a matrix, $f(\mathbf A)$, where $\mathbf A=(A_{ij})$, you can take derivative with regard to element $A_{ij}$: $$\frac{\partial f(\mathbf A)}{\partial A_{ij}}.$$
This is well defined in any calculus text. But in machine learning community, people assembly these simple derivatives into a matrix, or "write this
result more compactly" as said by Bishop's book:
$$\frac{\partial f(\mathbf A)}{\partial\mathbf A}=\left(\frac{\partial f(\mathbf A)}{\partial A_{ij}}\right).$$
People who knows Matrix calculus would recall that this is a so-called denominator layout notation. But as that wikipage says, machine learning community does not follow a single layout for different derivative types. So, it is wise to figure out what layout the author is following before assuming something. At least Andrew Ng is following denominator layout notation for scalar-by-matrix derivatives.
Let me start by citing simple product rule of derivative $$\bigl(f(x)g(x)\bigr)'=f'(x)\color{green}{g(x)}+\color{green}{f(x)}g'(x).$$ This rule means that by taking derivative, a single product becomes a sum of two products, with one factor being further derived but another factor being considered a constant which I marked in a relieving color of green.
Now let's think of trace of the product of two matrices. Here the two matrices are general; their product does not have to be square. Yes, this is the way things are in machine learning community. Trace is the sum of leading diagonal, but unsquared matrices have no problem having a leading diagonal (see here).
Now let's think that the elements of the above two matrices are all functions of $A_{ij}$, and next take derivative of the trace of the multiplication of two matrices w.r.t. $A_{ij}$. Only one summand of the trace is involved, right? That single summand comes from one row of the left matrix and one column of the right matrix, again a sum of products, right? Each summand of such a sum has a form of $f(A_{ij})g(A_{ij})$, so if we take the derivative, it splits into a sum of two products:
$$f'(A_{ij})\color{green}{g(A_{ij})}+\color{green}{f(A_{ij})}g'(A_{ij}),$$
where the green part is considered a constant.
Now let's expand the above to matrix form because the same thing happens for all $A_{ij}$'s. We have
$$\eqalign{
\frac{\partial\operatorname{Tr}\bigl(f(\mathbf A)g(\mathbf A)\bigr)}{\partial\mathbf A}&=\operatorname{Tr}\bigl(\frac{f(\mathbf A)}{\partial\mathbf A}\color{green}{g(\mathbf A)}+\color{green}{f(\mathbf A)}\frac{\partial g(\mathbf A)}{\partial\mathbf A}\bigr)\\
&=\frac{\partial\operatorname{Tr}\bigl(f(\mathbf A)\color{green}{g(\mathbf A)}\bigr)}{\partial\mathbf A}+\frac{\partial\operatorname{Tr}\bigl(\color{green}{f(\mathbf A)}g(\mathbf A)\bigr)}{\partial\mathbf A},
}$$
where the green parts at the right hand are regarded as constant matrix. Here I used the fact $\operatorname{Tr}(\mathbf A_1+\mathbf A_2)=\operatorname{Tr}\mathbf A_1+\operatorname{Tr}\mathbf A_2$ and linearity of derivative so $\partial\mathbf A$ can penetrate $\operatorname{Tr}$. This step requires some math thinking capability. I hope you can get through it.
Plugging $f(\mathbf A)=\mathbf A\mathbf B$ and $g(\mathbf A)=\mathbf A^T\mathbf C$, we get the second equation. The circle $\circ$ represents matrix $\mathbf A$ in Andrew Ng's derivation. The following derivation is easier because of established identities in page 698 of Bishop's book.
Now we have that the original derivative is equal to $$\frac{\partial}{\partial\mathbf A}\operatorname{Tr}(\mathbf A\mathbf B\color{green}{\mathbf A^T\mathbf C})+\frac{\partial}{\partial\mathbf A}\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\mathbf C).$$ Since matrix $\mathbf B$ and $\mathbf C$ are constant too w.r.t. elements $A_{ij}$ and hence w.r.t. matrix $\mathbf A$, we color them in green as well:
$$\frac{\partial}{\partial\mathbf A}\operatorname{Tr}(\mathbf A\color{green}{\mathbf B\mathbf A^T\mathbf C})+\frac{\partial}{\partial\mathbf A}\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C}).\tag{1}$$
Now we advance to the third equality of Andrew Ng's derivation of derivative. The first term of the above result (1) turns into $(\color{green}{\mathbf B\mathbf A^T\mathbf C})^T=\mathbf C^T\mathbf A\mathbf B^T$, according to identity (C.24) in page 698 of Bishop's book:
$$\frac{\partial}{\partial\mathbf A}\operatorname{Tr}(\mathbf A\color{green}{\mathbf B})=\color{green}{\mathbf B}^T.\tag{C.24}$$
Now you should know why I tirelessly mark constant matrix as green. The first term in the third equation of Andrew Ng's derivation is confusing in that he missed the trace. I am really a bit tired at this point. Take a break of typing, but you keep on reading.
The second term contains a transpose of $\mathbf A$. We want to take the derivative w.r.t. this transpose $\mathbf A^T$ to take advantage of identity (C.24). To that end, remember we use denominator layout to arrange partial derivatives in the resulting derivative matrix. If the variable matrix is $\mathbf A^T$, the denominator layout will result in nothing but a transpose of the matrix $\frac{\partial}{\partial\mathbf A}$. So, we have
$$\frac{\partial}{\partial\mathbf A^T}\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C})=\left(\frac{\partial}{\partial\mathbf A}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C})\right)^T,$$
which leads to the second term of the third equality.
In order to derive $\frac{\partial}{\partial\mathbf A^T}\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C})$ based on identity (C.24) where the variable matrix is the left matrix of the multiplication, we apply the cyclic property of trace to deduce $\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C})=\operatorname{Tr}(\mathbf A^T\color{green}{\mathbf C\mathbf A\mathbf B}).$ In the second term of the fourth equality, Andrew Ng wrote $Cf(A)$ out of the $\operatorname{Tr}$ operator. I think that's a typo.
Now it is easy to apply Bishop's identity (C.24) to obtain
$$\eqalign{
\frac{\partial}{\partial\mathbf A}\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C})
&=\left(\frac{\partial}{\partial\mathbf A^T}\operatorname{Tr}(\color{green}{\mathbf A\mathbf B}\mathbf A^T\color{green}{\mathbf C})\right)^T\\
&=\left(\frac{\partial}{\partial\mathbf A^T}\operatorname{Tr}(\mathbf A^T\color{green}{\mathbf C\mathbf A\mathbf B})\right)^T\\
&=\bigl((\color{green}{\mathbf C\mathbf A\mathbf B})^T\bigr)^T\\
&=\mathbf C\mathbf A\mathbf B.}$$
We have proved that the first term of (1) equals $\mathbf C^T\mathbf A\mathbf B^T$. Above we get the second term $\mathbf C\mathbf A\mathbf B$. Now we can conclude the whole proof by combining them as the final result of the derivative: $\mathbf C^T\mathbf A\mathbf B^T+\mathbf C\mathbf A\mathbf B$.
