serac
serac copied to clipboard
Implement useful tensor operations
There are some core tensor operations that are needed, as well as some that would be nice to have.
In my opinion, these are the needs. They're all for 2D and 3D 2nd order tensors:
- Spectral decomposition for symmetric matrices
- Matrix exponential for symmetric matrices and matrix log function for SPD matrices
- Matrix square root for symmetric SPD matrices
- Polar decomposition for real matrices with positive determinant
Once the spectral decomposition is in place, we could use it to write convenience functions for the others.
Nice to have:
- Extend the linear solve function to work on matrices of right hand sides (ie
A * x = B, whereBis a matrix of RHS vectors). This would make it possible to avoid most calls to the matrix inverse. - A representation of symmetric tensors, at least for 2nd order tensors. This would automatically guard some algorithms that only make sense on symmetric matrices and would save a little storage.
- Some means of invoking linear algebra libraries on 4th order tensors, especially if they have at least major symmetry. Computing eigenvalues of elastic tensors is the main use case I'm thinking about. Being able to cast the 4-th order tensors as 9x9 matrices (or 6x6 and symmetric 6x6 when there are major/minor symmetries) is handy.
- More flexible tensor contractions. It's hard to give all of the possible contractions meaningful names, so it would be nice to be able to express things in index notation. Numpy's
einsum(...)function is one possible model (ie,C_ijml, F_km -> B_ijkl) - Convenience function for ranges, like Python's
a = range(5). Sam posted some ideas here