One thing to point out is that because the Paulis form a basis, you can actually represent any $2^n \times 2^n$ matrix in terms of a sum of tensors of Paulis, i.e., members of the $n$-qubit Pauli group. That is, you can write any $8 \times 8$ matrix $A$ as a sum of the form
$$A =\sum_{i,j,k}h_{ijk}\ \sigma_i\otimes\sigma_j\otimes\sigma_k$$
where $h_{ijk}$ are the coefficients in the Pauli basis. The answer to this related question describes how to solve for these coefficients.
Now, once you have expressed $A$ in the Pauli basis, you can use the Q# operations you mention to implement the evolution. Some more background on implementing this in Q# is available at this link. The general idea is that once you have expressed the matrix $A$ in the Pauli basis, you can now use something like a Trotter–Suzuki expansion to approximately express the exponential $e^{iA}$ as a product of exponentials of Paulis, which can then in turn be implemented on a quantum computer (and also via built-in Q# tools such as Exp
and PauliEvolutionSet
).
Why doesn't a language like Q# include a built-in library for implementing a matrix exponential $e^{iA}$ for some general matrix $A$? Essentially because such an operation is, in general, extremely inefficient to implement on a quantum computer. To understand why, note that for a general $n$-qubit unitary, there are $4^n$ coefficients required to represent it in a basis like the Pauli basis, which means that your resulting circuit depth will be on the order of $4^n$ -- far too deep to be practical for anything beyond very small systems.
The exception is the case where the matrix $A$ has "sparsity" in some representation -- for example, if only a constant number of the $4^n$ coefficients in the Pauli basis are non-zero. In that case, the circuit resulting from a Trotter-Suzuki decomposition would have only constant depth, rather than going as $4^n$.