I am trying to simplify the following expression --$$(\vec{A} \times \vec{B})\cdot\underline{\underline{S}}\cdot(\vec{C} \times \vec{D}) = (A_jB_k\epsilon_{jki})S_{im}(\epsilon_{mnp}C_nD_p)$$ $S$ is a symmetric $3\times 3$ matrix, $\epsilon$ is the antisymmetric tensor object.
Can one get rid of the $\epsilon$-s somehow using the symmetry of the matrix $S$?
This doesn't use the symmetry of $S$, but you can use the following identity:
$$\epsilon_{ijk}\epsilon_{mnp}=\delta_{im}\delta_{jn}\delta_{kp}+\delta_{in}\delta_{jp}\delta_{km}+\delta_{ip}\delta_{jm}\delta_{kn}-\delta_{im}\delta_{jp}\delta_{kn}-\delta_{ip}\delta_{jn}\delta_{km}-\delta_{in}\delta_{jm}\delta_{kp}$$
to get
$$(A \times B)\cdot S\cdot(C \times D) = \text{tr}(S)\left[(A\cdot C)(B\cdot D)-(A\cdot D)(B\cdot C)\right]-(A\cdot C)(B\cdot S\cdot D)-(A\cdot S\cdot C)(B\cdot D)+(A\cdot D)(B\cdot S\cdot C)+(A\cdot S\cdot D)(B\cdot C).$$
Have fun with that!
(Note that the identity given above is just a more fundamental version of the (more commonly taught) identity
$$\epsilon_{ijk}\epsilon_{inp}=\delta_{jn}\delta_{kp}-\delta_{nk}\delta_{jp}$$ which you can derive from the above identity by contracting on $i$ and $m$)