Conclusion
Adapting the Invariant Coordinate Selection method to compositional data can, as it is often the case when adapting general multivariate data analysis methods, be done through applying the classical ICS method after an \(\mathop{\mathrm{ilr}}\) transformation.
Obviously, this does not provide a complete answer to the problem, since when transforming back through \(\mathop{\mathrm{ilr}}^{-1}\) we lose most of the properties of ICS, namely that the directions of the outliers are the first coordinates.
Moreover, the ICS transformation matrix \(\mathop{\mathrm{H}}\) depends on the choice of contrast matrix, even with an adjustment such as discussed at the end of The effects of changing the contrast matrix.
A first way to avoid this dependence would be to choose a contrast matrix, for instance one that simultaneously diagonalizes the scatter pair, or one that transforms the orthogonal projection on the \(\mathop{\mathrm{ilr}}\)-space into a subcomposition on the simplex, using balances (as in Section 4.8 of the book by Pawlowsky-Glahn, Egozcue, and Tolosana-Delgado (2015)).
In order to find a more intrinsic method, one could look for a link between the eigenvalues of the scatter-difference matrix \({\hat{\mathop{\mathrm{S}}}_1}^{-1} \hat{\mathop{\mathrm{S}}}_2\) and a notion of kurtosis according to the axes of the simplex (which are not orthogonal), instead of the approach that has been developed here which relies on the kurtosis according to the \(\mathop{\mathrm{ilr}}\) coordinates.
Another approach would be to create a good notion of correlation between simplex coordinates (see the article by Quinn et al. (2017) which explains why correlation is not appropriate for relative data and suggests more pertinent methods).
Naturally, results like Theorem 3 in the article by Tyler et al. (2009) (proving that under a Gaussian mixture model, the direction of the outliers can be found on the first or last Invariant Coordinates) could be adapted to the \(\mathop{\mathrm{ilr}}\) coordinates, but one could hope (if a more intrinsic method is found) to find similar results directly in simplex coordinates.