I am trying to transform the ODE $\partial_t \rho = \mathcal{L}(\rho)$ for a $N\times N$ matrix $\rho$ $$\rho = \frac{1}{N}I_N + \sum_{j} d_j s_j$$ to a more suitable basis $s_j$ (the generalized Pauli matrices) using the Hilbert-Schmidt product. The transformed ODE reads $\partial_t \vec d = M \vec d + \vec d_\infty$, where $\vec d$ is a vector now and the matrix $M$ has the entries $$ M_{ij} = \mathrm{Tr}\{\mathcal{L}(s_j) s_i\}.$$ Ignoring the vector $\vec d_\infty$, the solution of this ODE is $\vec d(t) = \exp (M t) \vec d(0)$. On the other hand, the solution of the original ODE would be $\rho (t) = \exp (\mathcal{L} t) \rho (0)$, which can be transformed to $$d_i(t) = \mathrm{Tr}\{ \exp (\mathcal{L} t) \rho (0)s_i\}.$$ As a crosscheck, I would like to confirm that both solutions are equivalent, i.e., $$\mathrm{Tr}\{ \exp (\mathcal{L} t) \rho (0)s_i\} = \sum_j \exp (Mt)_{ij} \mathrm{Tr}\{ \rho (0) s_j\}.$$
Any ideas how to accomplish this?
For details about the application in particular see my question in physics.SE: https://physics.stackexchange.com/questions/390458/adjoint-representation-in-liouville-von-neumann-equation