I have tried to compute the 1st order approximation using Taylor's expansion of the Mahalanobis distance:
$f(\mathbf{X})=\mathbf{a^TXa}$, where $\mathbf{a}\in \mathbb{R}^N$.
The function maps $\mathbf{X}\in \mathbb{R}^{N\times N}$ to $\mathbb{R}$. As the function $f$ is first order differentiable,
$\partial f/\partial\mathbf{X}=\mathbf{aa^T}$,
so the 1st order approximation at $\mathbf{X=X}_0$ would be like:
$f(\mathbf{X})\approx f(\mathbf{X}_0)+ \nabla f(\mathbf{X}_0)(\mathbf{X-X}_0) = \mathbf{a^TX_0a}+\mathbf{aa^T}(\mathbf{X-X}_0)$
However, the first part of the approximation $\mathbf{a^TX_0a}$ is a scalar, but the second part of the approximation $\mathbf{aa^T}(\mathbf{X-X}_0)$ is a matrix. This is definitely incorrect.
Anyone knows how to compute this in the right way? Thanks.
The reason for using the Frobenius product, won't quite fit into a comment.
Let's say you have a 3rd order tensor $F_{ijk}$ which is a function of a 2nd order tensor $X_{pq}$
Then the derivative of $F$ must be a 5th order tensor $\frac {\partial F_{ijk}} {\partial X_{pq}} = Q_{ijkpq}$
If you wish to use this quantity in differentials or Taylor expansions, then it needs to have this nice property $$\eqalign{ dF_{ijk} &= \sum_{p} \sum_{q} Q_{ijkpq}\,\,dX_{pq} \cr &= Q_{ijkpq}\,\,dX_{pq} \cr } $$ The last line employs the widely used "index notation" from physics and engineering. In index notation, it is clear that you are contracting over 2 indices. That is one of its many strengths.
Writing this same expression in matrix notation requires some way to indicate the double contraction, e.g. $\,\,dF = Q:dX$
In general, you need to contract over all of the indices in the independent variable. For example, if $X$ were a 3rd order tensor, then you would use a triple dot product to calculate $dF$.