From your example, it seems that A \in \mathbb{R}^{n\times M} is a “wide” matrix, i.e. more columns than rows (M \ge n). In this case, maybe use the LQ factorization?
That is, lq(A) forms a factorization A=LQ, where L \mathbb{R}^{n\times n} and Q \in \mathbb{R}^{M\times n} has orthonormal rows (QQ^T = I), equivalent to the QR factorization of A^T. Given this, \Sigma = AA^T = LQQ^T L^T = LL^T, and \Sigma^{-1} = L^{-T} L^{-1}. Hence your distance is simply d(x^i) = \Vert L^{-1} (x^i - \mu) \Vert.
In Julia:
L = LowerTriangular(lq(A).L)
dmahal = norm(L \ (x - μ))
Or, using your example code, I get:
julia> norm.(eachcol(LowerTriangular(lq(AX).L) \ AX)).^2 ≈ dmahal
true
(Note that you can apply L \ to every column of a matrix at once, which can be much more efficient than solving column-by-column, especially when the matrix is large.)
In general, the standard approaches to avoiding A^T A or AA^T computations (which square condition numbers) often involve QR or LQ factorization (respectively), or the SVD.
PS. I think the Mahalanobis distance is technically just the norm(...), not the norm squared as in your code.