Clearly document problems with mixing sympy / numpy

[greschd]

Matrices / arrays from sympy / numpy don't mix well when doing operations -- the resulting shapes are inconsistent and wrong.

TODO:

• Find places to enforce / assert having the same matrices, or correct matrix shapes
• Clearly document the use of the numeric parameter

[greschd]

List of issues with sympy matrices (incomplete):

• HDF5 serialization
• np.allclose in is_identity