Adding three distributions: Spherical, OrthogonalMatrices and NormalSingularValues (working names)

[velochy]

Over the long weekend, I dug into BPCA topic again and implemented Nirwan et al "Rotation Invariant Householder Parameterization for Bayesian PCA", which basically involved three distributions of potentially wider interest.

1. A distribution over direction vectors, i.e. a spherically symmetric distribution.
   This basically means putting a prior on the euclidean norm of the vector. There is a longer discussion on what prior to put on it in the stan forums, but I'm thinking of implementing it generically so it can be given as a parameter.

2. A uniform distribution over all (semi-)orthogonal matrices i.e. Stiefel manifold. This is done following Mezzadri "How to generate random matrices from the classical compact groups" and basically boils down to sampling random rotation vectors and using them as bases for sequential Householder transformations.

3. A distribution of the singular values of a matrix of i.i.d. standard normal distributions. This one is adapted from James, Lee "Concise Probability Distributions of Eigenvalues of Real-Valued Wishart Matrices"

While all three work for my current use case, there is definitely a lot of cleanup work to be done to make them generic. Most notably:

• Spherical distribution needs to take in a dist fn as a parameter to make it more general
• OrthogonalMatrix distribution should be vectorized further to allow for size parameter
• OrthogonalMatrix needs to be built into a proper distribution (more on that below)
• Adding tests and documentation

Before I take that on, though, there is some discussion to be had.

a) Do all three of the distributions make sense to add?
They are somewhat esoteric, but my gut feeling tells me they do have other uses than the BPCA I needed them for. Spherical is probably the easiest to find uses for of the three. Stan discussion linked above got sparked from finding a better geometry than atan2 for Von Mises distribution .

b) What would be good names for them?
NormalSingularValues is the only one I'm happy with out of the three right now, but it could also be eigenvalues of Wishart (in which case NSV would be a simple square root of that).
Orthogonal could be called Stiefel, but that is an esoteric term not that many people are familiar with. RandomOrthogonalMatrix is more descriptive, but it also covers Semi-Orthogonal matrices, for instance.
And Spherical also probably needs some clarifying words next to it but I can't even propose good ones.

c) How should I implement the OrthogonalMatrix as a proper distribution? It basically requires K draws from Spherical as its only source of randomness, so one option would be to allow spherical to take in a list of vector lengths and then output a single vector that is the concatenation, and then do OrthogonalMatrix as just a transformation on top of that (analogous to how Chi is implemented on top of ChiSquared). But this would add some rather odd complexity into the Spherical, so I'm wondering if there are any other good options.

In short - I need someone to sanity check all of this and help me think it through on the abstract level. @zaxtax, are you up for doing a third PR with me in a row?

[zaxtax]

Sure thing! I also made sure your last commits made it into the 0.2.5 release ^_^

I also suggest posting about this on the pymc discourse and see what others have to say.

[velochy]

Sure thing! I also made sure your last commits made it into the 0.2.5 release ^_^

I noticed. Thank you for that :)

I also suggest posting about this on the pymc discourse and see what others have to say.

Good idea. Created a post here

[ricardoV94]

Pretty cool, any chance we also get a pymc-exmaples notebook/extension to the PCA model?

[velochy]

Pretty cool, any chance we also get a pymc-exmaples notebook/extension to the PCA model?

I can take that on once this is merged (as it's hard to write without having all these distributions ;) ).
I learned a lot in revisiting BPCA topic and would be happy to clean up my experiments into a notebook.

[zaxtax]

I think it would be good to have example notebooks that motivate these distributions. We have a notebooks section in the repo just for that reason!

[velochy]

Its in this same repo?
In this case yes, I can add it to the PR once we settle on the names for the variables.

[zaxtax]

Yep! I meant for this repo, at least to start.