Add orthomax factor rotations

[cyianor]

This is a start to add factor rotations to MultivariateStats.jl. The implemented algorithm is currently what is described in the review paper

Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen, 
"Sparse modeling of landmark and texture variability using the orthomax criterion," 
Proc. SPIE 6144, Medical Imaging 2006: Image Processing, 61441G (10 March 2006); 
doi: 10.1117/12.651293

with the addition of choosing a random rotation matrix if the identity matrix is not a good starting point (this is what Matlab does in rotatefactors and was suggested by @haotian127).

Currently there are two suggestions for interfaces: Either the standard fit interface or rotatefactors which is more similar to Matlabs interface. I am not a big fan of the fit interface since, contrary to the scikit-learn original where it presumably came from, is not really chainable since keyword parameters differ from one fit function to the other. However, if this is the preferred way of adding things to MultivariateStats.jl then I will not stand in the way.

I like the idea of having a different type for each rotation algorithm. This is what JuliaLang started to do with e.g. the Algorithm for SVD in the LinearAlgebra package as well. This way multiple dispatch can be used (as I currently have it in the rotatefactors interface) instead of a chain of isa() checks or checking a symbol switch such as :orthomax.

Many things are still missing

• Decide on proper interface: fit vs rotatefactors
• Helper functions? e.g. to access the rotation matrix or loadings in the returned FactorRotation object
• Is the special functionality for e.g. the FactorAnalysis object desirable?
• Add tests
• Add documentation to the Sphinx source
• Other rotations (oblique, based on GPA rotations, ...)

Addresses #67

[cyianor]

Addresses issue #67

[wildart]

4 space indentation, please.

[cyianor]

It seems the method I found in the article by Stegman et al. (2006) only works for small gamma, i.e. Varimax (gamma = 1) and Quartimax (gamma = 0) rotations. Tests comparing to R's varimax and R's GPArotation::quartimax give the correct results.

For larger gamma I cannot reproduce the results with established software. Looking at Matlab they also use a different algorithm for gamma > 1. I guess it will be most worthwhile to implement the methods from GPArotation since these seem most flexible for extension to other factor rotations than the Orthomax family. I'll continue with that in a little while when I have time.

[wildart]

Currently there are two suggestions for interfaces: Either the standard fit interface or rotatefactors which is more similar to Matlabs interface. I am not a big fan of the fit interface since, contrary to the scikit-learn original where it presumably came from, is not really chainable since keyword parameters differ from one fit function to the other. However, if this is the preferred way of adding things to MultivariateStats.jl then I will not stand in the way.

I think fit usage is not good choice for this functionality. The factor (projection) matrix is a part of the statistical model object, and hardly can be used separate. So, following an established convention of this package, there should be two layers of the exposed functionality: low-level functions that are applied directly to the factor matrix, and high-level interface that is applied to statistical method objects, e.g. PCA or FactorAnalysis.

The low-level functions should be named after the particular rotation method, and accept the factor matrix as a parameter with any additional parameters, default or kw arguments, e.g.

varimax(F, γ = 1.0, maxiter = 10, tol = 1e-6)

The high-level interface has to operate on stat. model objects with additional parameters for rotation method specification. Rotation algorithm can be encoded as a type, so it would be easy to specialize a particular implementation.

M = fit(PCA, X)
rotate!(M, Orthmax)
rotate!(M, Varimax; γ = 1.0, maxiter = 10, tol = 1e-6)

The high-level rotation function retrieves model's factor matrix, projection(M), calls low-level function on it, and performs in-place update of the model's factor matrix. Because, the stat. model object is usually immutable, for the models that allow rotations an additional method for updating factor matrix has to be added to avoid direct access to the model properties.

[cyianor]

Thanks for the feedback! Will continue to work on this when time allows.

[cyianor]

Here is finally a new version. I replaced the old algorithm with the method that is used in the GPArotation package for R, since it is easy this way to implement all kinds of factor rotations. Especially, it is possible to provide oblique as well as orthogonal rotations.

I still need to write more tests and documentation, but I wanted to double check on the interface first.

The numerical crunch is done by the two versions of the gparotate function, one for orthogonal and one for oblique rotations. Using types for the rotation criteria turned out very beneficial to avoid symbol switch branches. Each rotation criterion requires a function that computes the criterion value and a gradient. By using types, Julia automatically picks the right one and criteria that require parameters have them readily accessible in the struct. By making criteria parametric types and introducing types for Orthogonal and Oblique rotation methods it is possible to let Julia pick the correct version of gparotate, especially for cases where criteria can be used in the orthogonal as well as in the oblique setting.

One concern could be that introducing more criteria will lead to more types. Not sure if that is a concern though.

In terms of user interface, I am offering one general rotate method which can be applied to any arbitrary matrix. In addition to the rotated loadings it also returns the matrix used to perform the rotation. Then there is the rotate! method which can be applied to statistical models (currently PCA and FactorAnalysis) and which modifies the projection/loading matrix in place. By using the .= operator I didn't have to add additional functions to the models to modify these matrices.

Is this interface in line with the general package philosophy? I've seen that symbols are used in other places to select methods, which could of course be a possibility here too, but it would lead to long switch statements that need to be modified each time a new rotation criterion gets added.

[cyianor]

Just to clarify, the interface would be

M = fit(FactorAnalysis, X)
rotate!(M, CrawfordFerguson{Orthogonal}(κ = 0.0))

or

rotate!(M, Oblimin{Oblique}(ɣ = 0.0))

when applied to existing models (currently PCA and FA). If one simply has a matrix of loadings L one can also call

Lrot = rotate(L, Oblimin{Oblique}(ɣ = 0.0))

There could be convenience types like Varimax() which is a very standard rotation and is equivalent to Oblimin{Orthogonal}(ɣ = 1.0).

[wildart]

The rotate interface usage, you presented is good. Just one thing I'd like you to change. Please, provide latinised parameters for rotation constructors instead of Unicode symbols. Using Unicode for variables internally is ok, but Unicode parameters sometimes could be tricky for library user to handle, especially new ones.

[cyianor]

Thank you for the feedback! I addressed what you pointed out and made the appropriate changes. I will add the missing documentation, some more factor rotation types, and appropriate testing code soon.

[p-gw]

What is the status on getting this feature released? I am very interested and would be happy to help if there is anything missing still.

[cyianor]

I haven't managed to work on this in a long time since my projects moved in another direction.

The code to perform factor rotations is basically finished and ready to use. It should be reasonably well documented in facrot.jl.

The major road block to getting the PR out is testing and writing more documentation. Also it's been 1.5 years since I last touched it, so I guess it will have to be synchronized with an updated version of the package.

Compared to something like the R package GPArotation it also supports quite few rotations currently. But I think that is more of a "would be nice to have" thing than a "must have" thing.

[cyianor]

Looking at the code now, there might be enough tests to get this off the ground. What do you think, @wildart?

I think the one thing missing is to update the Sphinx documentation.

[EDIT] Oh I see there is a Documenter.jl documentation now also. Is it enough to update that one or is the Sphinx documentation still in use as well?

[p-gw]

Thanks for the heads-up. I ended up writing my own implementation in a separate package. I think it has somewhat more functionality than proposed here, such as 'gradient-free' methods via automatic differentiation.

If you are interested you can take a look: https://github.com/p-gw/FactorRotations.jl