Antiderivative of a matrix Assume an $M$ by $N$ matrix $Y$ whose entries are given by
$$Y_{i,j} = \mbox{tr} \left( C_{i,j}^T X \right)$$ where $\{C_{i,j}\}_{i,j=1}^{M,N}$ and $X$ are $M$ by $N$ matrices. What is the antiderivative of $Y$ with respect to $X$?

Replace the matrix $Y$ by $F$ (since $y$ is used for something different in the update section), and write it as the double-dot product between the 4th order tensor ${\cal C}$ and the matrix $X$, i.e.
$$F= {\cal C}:X \quad\implies F_{ij}= \sum_{pq} {\cal C}_{ijpq} X_{pq}$$
This can be flattened into a matrix-vector relationship 
$$x={\rm vec}(X),\quad C={\rm vecvec}({\cal C}),\quad f=C\cdot x={\rm vec}(F)$$
Your question is analogous to Hooke's Law in a linearly elastic medium.
Consider the (scalar) stress-strain energy and its differential.
$$\eqalign{
\varepsilon &= \tfrac{1}{2}x\cdot C\cdot x \\
d\varepsilon &= (C\cdot x)\cdot dx = f\cdot dx \\
}$$
Integrating yields
$$\eqalign{
\int f\cdot dx
 &= \int d\varepsilon = \varepsilon &= \tfrac{1}{2}x\cdot C\cdot x \\
 &= \int F:dX &= \tfrac{1}{2}X:{\cal C}:X \\
}$$
For this analysis to hold, the stiffness tensor must be symmetric, i.e.
$\;{\cal C}:X = X:{\cal C}$
Update
When $C$ is constant but not symmetric, you can pull it out of the integral and write 
$$\eqalign{
\int f\cdot dx = C^T:\int x\,dx^T \\
}$$
Consider the matrix integral on the RHS in $2\times 2$ case.
$$\eqalign{
M=\int x\,dx^T 
&= \int\left[\matrix{x\\y}\right] \left[\matrix{dx&dy}\right] \\
&= \left[\matrix{\int x\,dx&\int x\,dy\\\int y\,dx&\int y\,dy}\right]
=\left[\matrix{\frac{1}{2}x^2&\int x\,dy\\\int y\,dx&\frac{1}{2}y^2}\right]\\
}$$
The problem is that those off-diagonal terms are path-dependent, i.e. you cannot integrate them until you specify the path of integration $-$ the functions $y=g(x)$ and $x=h(y)$.
In physics such path-dependent line integrals indicate a non-conservative field. You can still perform the integration, but you must specify a path, and if the path is not simple then you must resort to numerical techniques.
However, if $C$ is symmetric then you have a conservative field. In this case, there exists a potential $(\varepsilon)$ and the integral is simply
$\;\int\nabla\varepsilon\cdot dx=\int d\varepsilon = \varepsilon\;$ the value of which is the same for all paths, i.e. it's path-independent.
Note that 
$\,\Big(\int x\,dy + \int y\,dx\Big) = \int d(xy) = xy\;$
 and therefore
$$\left(M+M^T\right) = \left[\matrix{x^2&xy\\xy&y^2}\right] \;=\; xx^T$$
which is another route to solving the symmetric case.
