refactor local model code to rely on symmetry explicitly, using a decomposition

[tpapp]

Currently we use B=J'J which is positive semidefinite by construction, but that information is implicitly hardcoded into a lot of functions. The relevant code should be revamped:

1. using J=QR, store R for B=R'R. (double check algebra), in a wrapper type, eg called SymmetricProduct,
2. remove _factorize, just implement \ using the decomposition above and benchmark,
3. reimplement ellipsoidal_norm using this type.
4. expose checking if B is positive definite using the decomposition.
5. add AbstractMatrix(...) for this type (for use in hvcat/GES)

Then, later on, if a non-symmetric local model is considered, either defaults will be used or we can code that nicely, without making implicit assumptions about B like we currently do now.