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$?
 A: 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.
A: Rephrasing, we have the following problem.

Given matrix $\mathrm A \in \mathbb R^{n \times n}$, find the scalar field $f : \mathbb R^n \to \mathbb R$ such that $\nabla_{\rm{x}} f ( \rm{x} ) = \rm{A} \rm{x}$.

If matrix $\rm{A}$ is not symmetric, I do not believe function $f$ exists.

$2 \times 2$ case
$$\begin{aligned} \partial_1 f (x_1, x_2) &= a_{11} x_1 + a_{12} x_2\\ \partial_2 f (x_1, x_2) &= a_{21} x_1 + a_{22} x_2\end{aligned}$$
Integrating the first equation,
$$f (x_1, x_2) = \frac12 a_{11} x_1^2 + a_{12} x_1 x_2 + g (x_2)$$
Taking the partial derivative with respect to $x_2$ and using the 2nd equation,
$$g' (x_2) = a_{22} x_2 + (\underbrace{a_{21} - a_{12}}_{=0}) x_1 = a_{22} x_2$$
Note that if $a_{21} \neq a_{12}$ then $g$ is also a function of $x_1$, which is nonsense.
