I recently discovered a unit test file in the SymPy repository that demonstrates how to use SymPy for matrix calculus, specifically for finding derivatives of symbolic matrix expressions. This is of course very useful when working with optimization problems in e.g. machine-learning. The point is that SymPy can do this directly (in matrix form), yet this is not really obvious from the available documentation / content on forums.