You can do deterministic differentiation of matrix exponentation with a variety of algorithms, e.g. ChainRules.jl implements the algorithm in Al Mohy et al. (2009). (And, in particular, if you eventually are computing a scalar-valued function, e.g. an ML “loss” function, with the matrix exponential as an intermediate step, then the chain rule can be carried efficiently through the matrix exponential.)

Part of the confusion in the math.stackexchange question you linked is how to think about what the derivative *is*, and in particular the proposed answer \frac{\partial \exp(DA)}{\partial (DA)} = DA \exp(DA) is a category error: it’s not even the right *kind* of object for a matrix derivative. What you need is a Fréchet derivative, a linear operator on matrices, which is *not* simply a matrix of the same size. See e.g. our matrix calculus course, or this short talk on the nature of the problem.