From e271af6e730b9154da36dde6fccab6e865cfc097 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 15:42:46 +0200 Subject: [PATCH 01/32] Add orthomax factor rotations --- src/MultivariateStats.jl | 164 ++++++++++++++++++----------------- src/facrot.jl | 179 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 265 insertions(+), 78 deletions(-) create mode 100644 src/facrot.jl diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index e3ec71f..02df790 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -2,7 +2,7 @@ module MultivariateStats using LinearAlgebra using StatsBase: SimpleCovariance, CovarianceEstimator import Statistics: mean, var, cov, covm - import Base: length, size, show, dump + import Base: length, size, show, dump, eltype import StatsBase: fit, predict, ConvergenceException import SparseArrays import LinearAlgebra: eigvals @@ -10,109 +10,116 @@ module MultivariateStats export ## common - evaluate, # evaluate discriminant function values (imported from Base) - predict, # use a model to predict responses (imported from StatsBase) - fit, # fit a model to data (imported from StatsBase) - centralize, # subtract a mean vector from each column - decentralize, # add a mean vector to each column - indim, # the input dimension of a model - outdim, # the output dimension of a model - projection, # the projection matrix - reconstruct, # reconstruct the input (approximately) given the output - transform, # apply a model to transform a vector or a matrix - - # lreg - llsq, # Linear Least Square regression - ridge, # Ridge regression - - # whiten - Whitening, # Type: Whitening transformation - - invsqrtm, # Compute inverse of matrix square root, i.e. inv(sqrtm(A)) - invsqrtm!, # Compute inverse of matrix square root inplace - cov_whitening, # Compute a whitening transform based on covariance - cov_whitening!, # Compute a whitening transform based on covariance (input will be overwritten) - invsqrtm, # Compute C^{-1/2}, i.e. inv(sqrtm(C)) + evaluate, # evaluate discriminant function values (imported from Base) + predict, # use a model to predict responses (imported from StatsBase) + fit, # fit a model to data (imported from StatsBase) + centralize, # subtract a mean vector from each column + decentralize, # add a mean vector to each column + indim, # the input dimension of a model + outdim, # the output dimension of a model + projection, # the projection matrix + reconstruct, # reconstruct the input (approximately) given the output + transform, # apply a model to transform a vector or a matrix + + ## lreg + llsq, # Linear Least Square regression + ridge, # Ridge regression + + ## whiten + Whitening, # Type: Whitening transformation + + invsqrtm, # Compute inverse of matrix square root, i.e. inv(sqrtm(A)) + invsqrtm!, # Compute inverse of matrix square root inplace + cov_whitening, # Compute a whitening transform based on covariance + cov_whitening!, # Compute a whitening transform based on covariance (input will be overwritten) + invsqrtm, # Compute C^{-1/2}, i.e. inv(sqrtm(C)) ## pca - PCA, # Type: Principal Component Analysis model + PCA, # Type: Principal Component Analysis model - pcacov, # PCA based on covariance - pcasvd, # PCA based on singular value decomposition of input data - principalratio, # the ratio of variances preserved in the principal subspace - principalvar, # the variance along a specific principal direction - principalvars, # the variances along all principal directions + pcacov, # PCA based on covariance + pcasvd, # PCA based on singular value decomposition of input data + principalratio, # the ratio of variances preserved in the principal subspace + principalvar, # the variance along a specific principal direction + principalvars, # the variances along all principal directions - tprincipalvar, # total principal variance, i.e. sum(principalvars(M)) - tresidualvar, # total residual variance - tvar, # total variance + tprincipalvar, # total principal variance, i.e. sum(principalvars(M)) + tresidualvar, # total residual variance + tvar, # total variance ## ppca - PPCA, # Type: the Probabilistic PCA model + PPCA, # Type: the Probabilistic PCA model - ppcaml, # Maximum likelihood probabilistic PCA - ppcaem, # EM algorithm for probabilistic PCA - bayespca, # Bayesian PCA - loadings, # factor loadings matrix + ppcaml, # Maximum likelihood probabilistic PCA + ppcaem, # EM algorithm for probabilistic PCA + bayespca, # Bayesian PCA + loadings, # factor loadings matrix ## kpca - KernelPCA, # Type: the Kernel PCA model + KernelPCA, # Type: the Kernel PCA model ## cca - CCA, # Type: Correlation Component Analysis model + CCA, # Type: Correlation Component Analysis model - ccacov, # CCA based on covariances - ccasvd, # CCA based on singular value decomposition of input data + ccacov, # CCA based on covariances + ccasvd, # CCA based on singular value decomposition of input data - xindim, # input dimension of X - yindim, # input dimension of Y - xmean, # sample mean of X - ymean, # sample mean of Y - xprojection, # projection matrix for X - yprojection, # projection matrix for Y - xtransform, # transform for X - ytransform, # transform for Y - correlations, # correlations of all projected directions + xindim, # input dimension of X + yindim, # input dimension of Y + xmean, # sample mean of X + ymean, # sample mean of Y + xprojection, # projection matrix for X + yprojection, # projection matrix for Y + xtransform, # transform for X + ytransform, # transform for Y + correlations, # correlations of all projected directions ## cmds MDS, - classical_mds, # perform classical MDS over a given distance matrix - eigvals, # eignenvalues of the transformation - stress, # stress evaluation + classical_mds, # perform classical MDS over a given distance matrix + eigvals, # eignenvalues of the transformation + stress, # stress evaluation - gram2dmat, gram2dmat!, # Gram matrix => Distance matrix - dmat2gram, dmat2gram!, # Distance matrix => Gram matrix + gram2dmat, gram2dmat!, # Gram matrix => Distance matrix + dmat2gram, dmat2gram!, # Distance matrix => Gram matrix ## lda - Discriminant, # Abstract Type: for all discriminant functionals - LinearDiscriminant, # Type: Linear Discriminant functional - MulticlassLDAStats, # Type: Statistics required for training multi-class LDA - MulticlassLDA, # Type: Multi-class LDA model - SubspaceLDA, # Type: LDA model for high-dimensional spaces - - ldacov, # Linear discriminant analysis based on covariances - - classweights, # class-specific weights - classmeans, # class-specific means - withclass_scatter, # with-class scatter matrix - betweenclass_scatter, # between-class scatter matrix - multiclass_lda_stats, # compute statistics for multiclass LDA training - multiclass_lda, # train multi-class LDA based on statistics - mclda_solve, # solve multi-class LDA projection given scatter matrices - mclda_solve!, # solve multi-class LDA projection (inputs are overriden) + Discriminant, # Abstract Type: for all discriminant functionals + LinearDiscriminant, # Type: Linear Discriminant functional + MulticlassLDAStats, # Type: Statistics required for training multi-class LDA + MulticlassLDA, # Type: Multi-class LDA model + SubspaceLDA, # Type: LDA model for high-dimensional spaces + + ldacov, # Linear discriminant analysis based on covariances + + classweights, # class-specific weights + classmeans, # class-specific means + withclass_scatter, # with-class scatter matrix + betweenclass_scatter, # between-class scatter matrix + multiclass_lda_stats, # compute statistics for multiclass LDA training + multiclass_lda, # train multi-class LDA based on statistics + mclda_solve, # solve multi-class LDA projection given scatter matrices + mclda_solve!, # solve multi-class LDA projection (inputs are overriden) ## ica - ICA, # Type: the Fast ICA model + ICA, # Type: the Fast ICA model - icagfun, # a function to get a ICA approx neg-entropy functor - fastica!, # core algorithm function for the Fast ICA + icagfun, # a function to get a ICA approx neg-entropy functor + fastica!, # core algorithm function for the Fast ICA ## fa - FactorAnalysis, # Type: the Factor Analysis model + FactorAnalysis, # Type: the Factor Analysis model - faem, # Maximum likelihood probabilistic PCA - facm # EM algorithm for probabilistic PCA + faem, # Maximum likelihood probabilistic PCA + facm, # EM algorithm for probabilistic PCA + ## facrot + FactorRotationAlgorithm, # Type: Factor rotation algorithm + Orthomax, # Type: Orthomax factor rotation algorithm + + FactorRotation, # Type: Return type for factor rotations + + rotatefactors # Alternative interface to factor rotations ## source files include("common.jl") @@ -126,5 +133,6 @@ module MultivariateStats include("lda.jl") include("ica.jl") include("fa.jl") + include("facrot.jl") end # module diff --git a/src/facrot.jl b/src/facrot.jl new file mode 100644 index 0000000..8c07fc2 --- /dev/null +++ b/src/facrot.jl @@ -0,0 +1,179 @@ +""" + FactorRotationAlgorithm + +An abstract type for factor rotation algorithms. +""" +abstract type FactorRotationAlgorithm end + +""" + Orthomax <: FactorRotationAlgorithm + +A type representing the orthomax factor rotation algorithm. + +The positive parameter `γ::Real` determines which type of rotation is performed. +For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` +is quartimax rotation, `γ = p / 2` is equamax rotation, and +`γ = n (p - 1) / (n + p - 2)` is parsimax rotation. + +The parameter `maxiter::Integer` controls the maximum number of iterations to +perform (default `1000`) and `ϵ::Real` is a small positive constant determining +convergence (default `1e-12`). +""" +struct Orthomax <: FactorRotationAlgorithm + γ::Real + miniter::Integer + maxiter::Integer + ϵ::Real + + Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin + maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) + miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) + γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) + ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) + + new(γ, miniter, maxiter, ϵ) + end +end + +function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) + summary(io, alg); println(io) + println(io, "γ = $(alg.γ)") + println(io, "miniter = $(alg.miniter)") + println(io, "maxiter = $(alg.maxiter)") + println(io, "ϵ = $(alg.ϵ)") +end + +""" + FactorRotation{T <: Real} + +The return type for factor rotations. + +`F` contains the rotated factors and `R` is the rotation matrix that +was applied to the original factors. +""" +struct FactorRotation{T <: Real, Ta <: FactorRotationAlgorithm} + F::Matrix{T} + R::Matrix{T} + + alg::Ta + + FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) +end + +eltype(::Type{FactorRotation{T, Ta}}) where {T,Ta} = T + +FactorRotation(F::Matrix{T}, R::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation{T, Ta}(F, R, alg) +FactorRotation(F::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation(F, Matrix{eltype(F)}(I, size(F, 2), size(F, 2)), alg) + +function FactorRotation(F::Matrix, R::Matrix, alg::Ta) where {Ta <: FactorRotationAlgorithm} + return FactorRotation(promote(F, R)..., alg) +end + +function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any,<:Any}) + summary(io, FR); println(io) +end + + +## Comparison to other implementations +# The implementation of varimax in R row-normlises the input matrix before +# application of the algorithm and rescales the rows afterwards. +## Reference +# 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 +function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) + n, p = size(F) + if n < 2 + return (F, Matrix{eltype(F)}(I, p, p)) + end + + # Warm up step + # Do one step. If that first step did not lead away from the identity + # matrix enough use a random orthogonal matrix as a starting point. + M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) + R = M.U * M.Vt + if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ + R, _ = qr(randn(p, p)).Q + end + + # Main iterations + d = 0 + converged = false + for i in 1:maxiter + dold = d + B = F * R + M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) + R = M.U * M.Vt + d = sum(M.S) + if abs(d - dold) / d < ϵ && i >= miniter + converged = true + break + end + end + + if !converged + @warn "orthomax did not converge in $(maxiter) iterations" + end + + (F * R, R) +end + +""" + fit(::Type{FactorRotation}, F::AbstractMatrix; ...) -> FactorRotation + +Fit a factor rotation to the matrix `F` and apply it. The algorithm used to +perform the factor rotation is by default [`Orthomax`](@ref) and can be changed +with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. +""" +function fit(::Type{FactorRotation}, F::AbstractMatrix; + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + if isa(alg, Orthomax) + F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + return FactorRotation(F, R, alg) + end +end + +""" + fit(::Type{FactorRotation}, F::FactorAnalysis; ...) -> FactorRotation + +Fit a factor rotation to the loading matrix of the [`FactorAnalysis`](@ref) +object and apply it. + +See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. +""" +function fit(::Type{FactorRotation}, F::FactorAnalysis; + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + return fit(FactorRotation, F.W; alg) +end + +## Alternative interface + +""" + rotatefactors(F::AbstractMatrix, [alg::FactorRotationAlgorithm]) -> FactorRotation + +Rotate the factors in the matrix `F`. The algorithm to be used can be passed in +via the second argument `alg`. By default [`Orthomax`](@ref) is used. +""" +function rotatefactors(F::AbstractMatrix, alg::FactorRotationAlgorithm = Orthomax()) + F, R = _rotatefactors(F, alg) + return FactorRotation(F, R, alg) +end + +# Use multiple dispatch to decide on the algorithm implementation + +function _rotatefactors(F::AbstractMatrix, alg::Orthomax) + return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) +end + +""" + rotatefactors(F::FactorAnalysis, [alg::FactorRotationAlgorithm]) -> FactorAnalysis + +Rotate the factors in the loading matrix of `F` which is of type [`FactorAnalysis`](@ref). +The algorithm to be used can be passed in via the second argument `alg`. +By default [`Orthomax`](@ref) is used. +""" +function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) + FR = rotatefactors(F.W, alg) + return FactorAnalysis(F.mean, FR.F, F.Ψ) +end From 40067c6c445c598de864ff68178f671b138ad261 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 16:09:48 +0200 Subject: [PATCH 02/32] Make backwards-compatible with Julia 1.0 --- src/facrot.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/facrot.jl b/src/facrot.jl index 8c07fc2..4d2e824 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -144,7 +144,7 @@ See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword argumen """ function fit(::Type{FactorRotation}, F::FactorAnalysis; alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - return fit(FactorRotation, F.W; alg) + return fit(FactorRotation, F.W; alg = alg) end ## Alternative interface From b7a067ea77d27c60c1a8d3ecee4adb62ed94441d Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 16:19:16 +0200 Subject: [PATCH 03/32] Add documentation for `miniter` --- src/facrot.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index 4d2e824..2e77281 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -16,8 +16,9 @@ is quartimax rotation, `γ = p / 2` is equamax rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. The parameter `maxiter::Integer` controls the maximum number of iterations to -perform (default `1000`) and `ϵ::Real` is a small positive constant determining -convergence (default `1e-12`). +perform (default `1000`), `miniter::Integer` controls the minimum number of +iterations taken before convergence is checked, and `ϵ::Real` is a small positive +constant determining convergence (default `1e-12`). """ struct Orthomax <: FactorRotationAlgorithm γ::Real From caead0178295c256f6867650e54b0fa2833a787f Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 16:48:26 +0200 Subject: [PATCH 04/32] Fixed small mistake --- src/facrot.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/facrot.jl b/src/facrot.jl index 2e77281..cb4d392 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -95,7 +95,7 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) R = M.U * M.Vt if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ - R, _ = qr(randn(p, p)).Q + R = qr(randn(p, p)).Q end # Main iterations From afbb5a47d1b08cd045271fca1faf8ac451e006d3 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 21:48:03 +0200 Subject: [PATCH 05/32] Adjusted indent --- src/facrot.jl | 132 +++++++++++++++++++++++++------------------------- 1 file changed, 66 insertions(+), 66 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index cb4d392..1b7ee8a 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -21,27 +21,27 @@ iterations taken before convergence is checked, and `ϵ::Real` is a small positi constant determining convergence (default `1e-12`). """ struct Orthomax <: FactorRotationAlgorithm - γ::Real - miniter::Integer - maxiter::Integer - ϵ::Real - - Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin - maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) - miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) - γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) - ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) - - new(γ, miniter, maxiter, ϵ) - end + γ::Real + miniter::Integer + maxiter::Integer + ϵ::Real + + Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin + maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) + miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) + γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) + ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) + + new(γ, miniter, maxiter, ϵ) + end end function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) - summary(io, alg); println(io) - println(io, "γ = $(alg.γ)") - println(io, "miniter = $(alg.miniter)") - println(io, "maxiter = $(alg.maxiter)") - println(io, "ϵ = $(alg.ϵ)") + summary(io, alg); println(io) + println(io, "γ = $(alg.γ)") + println(io, "miniter = $(alg.miniter)") + println(io, "maxiter = $(alg.maxiter)") + println(io, "ϵ = $(alg.ϵ)") end """ @@ -53,12 +53,12 @@ The return type for factor rotations. was applied to the original factors. """ struct FactorRotation{T <: Real, Ta <: FactorRotationAlgorithm} - F::Matrix{T} - R::Matrix{T} + F::Matrix{T} + R::Matrix{T} - alg::Ta + alg::Ta - FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) + FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) end eltype(::Type{FactorRotation{T, Ta}}) where {T,Ta} = T @@ -67,11 +67,11 @@ FactorRotation(F::Matrix{T}, R::Matrix{T}, alg::Ta) where {T <: Real, Ta <: Fact FactorRotation(F::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation(F, Matrix{eltype(F)}(I, size(F, 2), size(F, 2)), alg) function FactorRotation(F::Matrix, R::Matrix, alg::Ta) where {Ta <: FactorRotationAlgorithm} - return FactorRotation(promote(F, R)..., alg) + return FactorRotation(promote(F, R)..., alg) end function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any,<:Any}) - summary(io, FR); println(io) + summary(io, FR); println(io) end @@ -84,40 +84,40 @@ end # Proc. SPIE 6144, Medical Imaging 2006: Image Processing, 61441G # (10 March 2006); doi: 10.1117/12.651293 function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) - n, p = size(F) - if n < 2 - return (F, Matrix{eltype(F)}(I, p, p)) - end - - # Warm up step - # Do one step. If that first step did not lead away from the identity - # matrix enough use a random orthogonal matrix as a starting point. - M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) - R = M.U * M.Vt - if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ - R = qr(randn(p, p)).Q - end - - # Main iterations - d = 0 - converged = false - for i in 1:maxiter - dold = d - B = F * R - M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) + n, p = size(F) + if n < 2 + return (F, Matrix{eltype(F)}(I, p, p)) + end + + # Warm up step + # Do one step. If that first step did not lead away from the identity + # matrix enough use a random orthogonal matrix as a starting point. + M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) R = M.U * M.Vt - d = sum(M.S) - if abs(d - dold) / d < ϵ && i >= miniter - converged = true - break + if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ + R = qr(randn(p, p)).Q end - end - if !converged - @warn "orthomax did not converge in $(maxiter) iterations" - end + # Main iterations + d = 0 + converged = false + for i in 1:maxiter + dold = d + B = F * R + M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) + R = M.U * M.Vt + d = sum(M.S) + if abs(d - dold) / d < ϵ && i >= miniter + converged = true + break + end + end - (F * R, R) + if !converged + @warn "orthomax did not converge in $(maxiter) iterations" + end + + (F * R, R) end """ @@ -128,11 +128,11 @@ perform the factor rotation is by default [`Orthomax`](@ref) and can be changed with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. """ function fit(::Type{FactorRotation}, F::AbstractMatrix; - alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - if isa(alg, Orthomax) - F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) - return FactorRotation(F, R, alg) - end + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + if isa(alg, Orthomax) + F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + return FactorRotation(F, R, alg) + end end """ @@ -144,8 +144,8 @@ object and apply it. See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. """ function fit(::Type{FactorRotation}, F::FactorAnalysis; - alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - return fit(FactorRotation, F.W; alg = alg) + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + return fit(FactorRotation, F.W; alg = alg) end ## Alternative interface @@ -157,14 +157,14 @@ Rotate the factors in the matrix `F`. The algorithm to be used can be passed in via the second argument `alg`. By default [`Orthomax`](@ref) is used. """ function rotatefactors(F::AbstractMatrix, alg::FactorRotationAlgorithm = Orthomax()) - F, R = _rotatefactors(F, alg) - return FactorRotation(F, R, alg) + F, R = _rotatefactors(F, alg) + return FactorRotation(F, R, alg) end # Use multiple dispatch to decide on the algorithm implementation function _rotatefactors(F::AbstractMatrix, alg::Orthomax) - return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) end """ @@ -175,6 +175,6 @@ The algorithm to be used can be passed in via the second argument `alg`. By default [`Orthomax`](@ref) is used. """ function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) - FR = rotatefactors(F.W, alg) - return FactorAnalysis(F.mean, FR.F, F.Ψ) + FR = rotatefactors(F.W, alg) + return FactorAnalysis(F.mean, FR.F, F.Ψ) end From 80e7316de92ae0bfdedbd480f3d5bc996d3267ce Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 00:22:15 +0200 Subject: [PATCH 06/32] Better exception use/getters for loadings/rotation --- src/MultivariateStats.jl | 4 +++- src/facrot.jl | 45 ++++++++++++++++++++++++++-------------- 2 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 02df790..8a1be64 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -118,7 +118,9 @@ module MultivariateStats Orthomax, # Type: Orthomax factor rotation algorithm FactorRotation, # Type: Return type for factor rotations - + loadings, # rotated loading matrix + rotation, # rotation matrix applied to loadings + rotatefactors # Alternative interface to factor rotations ## source files diff --git a/src/facrot.jl b/src/facrot.jl index 1b7ee8a..7117c1c 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -11,8 +11,8 @@ abstract type FactorRotationAlgorithm end A type representing the orthomax factor rotation algorithm. The positive parameter `γ::Real` determines which type of rotation is performed. -For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` -is quartimax rotation, `γ = p / 2` is equamax rotation, and +For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` +is quartimax rotation, `γ = p / 2` is equamax rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. The parameter `maxiter::Integer` controls the maximum number of iterations to @@ -27,10 +27,10 @@ struct Orthomax <: FactorRotationAlgorithm ϵ::Real Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin - maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) - miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) - γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) - ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) + γ ≥ zero(eltype(γ)) || throw(DomainError("Orthomax: γ needs to be non-negative")) + miniter > zero(eltype(miniter)) || throw(DomainError("Orthomax: miniter needs to be positive")) + maxiter > zero(eltype(maxiter)) || throw(DomainError("Orthomax: maxiter needs to be positive")) + ϵ > zero(eltype(ϵ)) || throw(DomainError("Orthomax: ϵ needs to be positive")) new(γ, miniter, maxiter, ϵ) end @@ -74,13 +74,26 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any summary(io, FR); println(io) end +""" + loadings(FR::FactorRotation) + +Returns the loading matrix. +""" +loadings(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.F + +""" + rotation(FR::FactorRotation) + +Returns the rotation matrix. +""" +rotation(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.R ## Comparison to other implementations # The implementation of varimax in R row-normlises the input matrix before # application of the algorithm and rescales the rows afterwards. -## Reference -# Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen, "Sparse modeling of -# landmark and texture variability using the orthomax criterion," +## Reference +# 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 function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) @@ -100,6 +113,7 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) # Main iterations d = 0 + lastchange = NaN converged = false for i in 1:maxiter dold = d @@ -107,15 +121,14 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) R = M.U * M.Vt d = sum(M.S) - if abs(d - dold) / d < ϵ && i >= miniter + lastchange = abs(d - dold) / d + if lastchange < ϵ && i >= miniter converged = true break end end - if !converged - @warn "orthomax did not converge in $(maxiter) iterations" - end + converged || throw(ConvergenceException(maxiter, lastchange, ϵ)) (F * R, R) end @@ -127,7 +140,7 @@ Fit a factor rotation to the matrix `F` and apply it. The algorithm used to perform the factor rotation is by default [`Orthomax`](@ref) and can be changed with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. """ -function fit(::Type{FactorRotation}, F::AbstractMatrix; +function fit(::Type{FactorRotation}, F::AbstractMatrix; alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} if isa(alg, Orthomax) F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) @@ -143,7 +156,7 @@ object and apply it. See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. """ -function fit(::Type{FactorRotation}, F::FactorAnalysis; +function fit(::Type{FactorRotation}, F::FactorAnalysis; alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} return fit(FactorRotation, F.W; alg = alg) end @@ -171,7 +184,7 @@ end rotatefactors(F::FactorAnalysis, [alg::FactorRotationAlgorithm]) -> FactorAnalysis Rotate the factors in the loading matrix of `F` which is of type [`FactorAnalysis`](@ref). -The algorithm to be used can be passed in via the second argument `alg`. +The algorithm to be used can be passed in via the second argument `alg`. By default [`Orthomax`](@ref) is used. """ function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) From c0e2389cebb6ca86bded6f6d1307849fec52c0d6 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 00:22:36 +0200 Subject: [PATCH 07/32] Created some tests --- Project.toml | 9 +++--- test/facrot.jl | 74 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 +- 3 files changed, 81 insertions(+), 5 deletions(-) create mode 100644 test/facrot.jl diff --git a/Project.toml b/Project.toml index dc0e941..774fe7f 100644 --- a/Project.toml +++ b/Project.toml @@ -13,11 +13,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +[compat] +StatsBase = ">=0.29" + [extras] +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] - -[compat] -StatsBase = ">=0.29" +test = ["Test", "Random"] diff --git a/test/facrot.jl b/test/facrot.jl new file mode 100644 index 0000000..501b3dd --- /dev/null +++ b/test/facrot.jl @@ -0,0 +1,74 @@ +using MultivariateStats +using LinearAlgebra +using Test +import Random + +@testset "Factor Analysis" begin + + Random.seed!(37273) + + X = randn(5, 2) + + ## Varimax rotation + # Comparison with R's `varimax` function called as (after `using RCall`) + # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" + # R"F <- qm$loadings" + # R"R <- qm$rotmat" + # @rget F + # @rget R + + F = [-2.56254778 -0.05712524; + -0.37387139 0.53683094; + -1.76868624 0.67513474; + 0.08057810 0.34826689; + -0.06196338 0.20433233] + R = [ 0.96256370 0.27105556; + -0.27105556 0.96256370] + + FR = fit(FactorRotation, X) + @test isapprox(FR.F, F, rtol = √eps()) + @test isapprox(FR.R, R, rtol = √eps()) + + ## Quartimax rotation + # Comparison with R's `GPArotation::quartimax` function called as + # R"qm <- GPArotation::quartimax($X, eps = 1e-6, maxit = 10000)" + # R"F <- qm$loadings" + # R"R <- qm$Th" + # @rget F + # @rget R + + F = [-2.55665706 -0.18280898; + -0.39976499 0.51783707; + -1.79968636 0.58752613; + 0.06339043 0.35180152; + -0.07191598 0.20104540] + R = [ 0.94810241 0.31796513; + -0.31796513 0.94810241] + + FR = fit(FactorRotation, X, alg = Orthomax(γ = 0.0)) + @test isapprox(FR.F, F, rtol = 1e-6) + @test isapprox(FR.R, R, rtol = 1e-6) + + ## Equamax + X = randn(10, 5) + + # Comparsion with Matlab's `rotatefactors` + # using MATLAB + # mat"rotatefactors($X, 'Method', 'equamax')" + + F = [-0.01195028 0.05771308 -0.35271260 2.83724330 -0.68373015; + 0.07085630 2.43064454 -0.56171788 0.00250704 0.05441219; + -0.12946985 0.78708715 -0.83039068 -2.01918172 0.65545648; + 1.95142102 -1.08251779 -0.49721317 1.13050103 -1.23799305; + -0.55556677 -0.16288095 1.56941619 -0.36283530 2.46150446; + 0.21075154 2.16186107 1.16767716 -1.03674456 0.71948514; + 0.09579139 0.15765568 1.21512619 0.03140918 0.17671690; + 0.00855404 0.96679333 -1.69727390 -0.20320960 -0.75931306; + 0.01368675 0.34330631 0.00514439 -0.80566209 1.21579113; + -2.13711263 -0.36494414 -0.25451404 0.25001421 -0.07785436] + + # Equamax for n x p matrix means γ = p / 2 + FR = fit(FactorRotation, X, alg = Orthomax(γ = 5 / 2, maxiter = 20000)) + @test isapprox(FR.F, F, rtol = √eps()) + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 7350589..fc4a0fb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,8 @@ tests = ["lreg", "ica", "ppca", "kpca", - "fa"] + "fa", + "facrot"] for test in tests include(test*".jl") From 2e0bfc3eebea0c2d2fed16704e3633bb96b7d32d Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 10:24:30 +0200 Subject: [PATCH 08/32] Convenience functions for common rotations --- src/facrot.jl | 87 ++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 72 insertions(+), 15 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index 7117c1c..f327b7a 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -10,15 +10,15 @@ abstract type FactorRotationAlgorithm end A type representing the orthomax factor rotation algorithm. -The positive parameter `γ::Real` determines which type of rotation is performed. -For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` -is quartimax rotation, `γ = p / 2` is equamax rotation, and -`γ = n (p - 1) / (n + p - 2)` is parsimax rotation. +The positive parameter `γ::Real` determines which type of rotation is +performed. For a `n x p` matrix, the default `γ = 1.0` leads to varimax +rotation, `γ = 0.0` is quartimax rotation, `γ = p / 2` is equamax +rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. -The parameter `maxiter::Integer` controls the maximum number of iterations to -perform (default `1000`), `miniter::Integer` controls the minimum number of -iterations taken before convergence is checked, and `ϵ::Real` is a small positive -constant determining convergence (default `1e-12`). +The parameter `maxiter::Integer` controls the maximum number of iterations +to perform (default `1000`), `miniter::Integer` controls the minimum number +of iterations taken before convergence is checked, and `ϵ::Real` is a small +positive constant determining convergence (default `1e-12`). """ struct Orthomax <: FactorRotationAlgorithm γ::Real @@ -26,16 +26,73 @@ struct Orthomax <: FactorRotationAlgorithm maxiter::Integer ϵ::Real - Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin - γ ≥ zero(eltype(γ)) || throw(DomainError("Orthomax: γ needs to be non-negative")) - miniter > zero(eltype(miniter)) || throw(DomainError("Orthomax: miniter needs to be positive")) - maxiter > zero(eltype(maxiter)) || throw(DomainError("Orthomax: maxiter needs to be positive")) - ϵ > zero(eltype(ϵ)) || throw(DomainError("Orthomax: ϵ needs to be positive")) - - new(γ, miniter, maxiter, ϵ) + Orthomax(;γ::Union{Real, Integer} = 1.0, + miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) = begin + γ ≥ zero(eltype(γ)) || + throw(DomainError("Orthomax: γ needs to be non-negative")) + miniter > zero(eltype(miniter)) || + throw(DomainError("Orthomax: miniter needs to be positive")) + maxiter > zero(eltype(maxiter)) || + throw(DomainError("Orthomax: maxiter needs to be positive")) + ϵ > zero(eltype(ϵ)) || + throw(DomainError("Orthomax: ϵ needs to be positive")) + + new(float(γ), miniter, maxiter, ϵ) end end +""" + Varimax() -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = 1.0`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Varimax(;miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end +""" + Quartimax() -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = 0.0`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Quartimax(;miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = 0.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end +""" + Equamax(p) -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = p / 2`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Equamax(p::Integer; miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = float(p) / 2.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end +""" + Parsimax(n, p) -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = n (p - 1) / (n + p - 2)`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Parsimax(n::Integer, p::Integer; miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = float(n) * (float(p) - 1) / (float(n) + float(p) - 2), + miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end + function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) summary(io, alg); println(io) println(io, "γ = $(alg.γ)") From 855ee0e0c8e46d40a969e5626308c4f8c435dae1 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 10:24:48 +0200 Subject: [PATCH 09/32] Corrected name of test set for factor rotations --- test/facrot.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/facrot.jl b/test/facrot.jl index 501b3dd..56226f4 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -3,14 +3,15 @@ using LinearAlgebra using Test import Random -@testset "Factor Analysis" begin +@testset "Factor Rotations" begin Random.seed!(37273) X = randn(5, 2) ## Varimax rotation - # Comparison with R's `varimax` function called as (after `using RCall`) + # Comparison with R's `varimax` function called as + # using RCall # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" # R"F <- qm$loadings" # R"R <- qm$rotmat" @@ -52,7 +53,7 @@ import Random ## Equamax X = randn(10, 5) - # Comparsion with Matlab's `rotatefactors` + # Comparison with Matlab's `rotatefactors` called as # using MATLAB # mat"rotatefactors($X, 'Method', 'equamax')" From 275d34a3e2a0336295d6f4b1421e884d2b48146c Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 16:08:23 +0200 Subject: [PATCH 10/32] Convenience wrappers for common Orthomax rotations --- src/MultivariateStats.jl | 4 +++ src/facrot.jl | 3 ++ test/facrot.jl | 70 ++++++++++++++++++++++++++-------------- 3 files changed, 53 insertions(+), 24 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 8a1be64..c7bce11 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -116,6 +116,10 @@ module MultivariateStats ## facrot FactorRotationAlgorithm, # Type: Factor rotation algorithm Orthomax, # Type: Orthomax factor rotation algorithm + Varimax, # Convenience interfaces to Orthomax factor rotation algorithms + Quartimax, + Equamax, + Parsimax, FactorRotation, # Type: Return type for factor rotations loadings, # rotated loading matrix diff --git a/src/facrot.jl b/src/facrot.jl index f327b7a..235de40 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -55,6 +55,7 @@ function Varimax(;miniter::Integer = 20, ϵ::Real = 1e-12) Orthomax(γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end + """ Quartimax() -> Orthomax @@ -67,6 +68,7 @@ function Quartimax(;miniter::Integer = 20, ϵ::Real = 1e-12) Orthomax(γ = 0.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end + """ Equamax(p) -> Orthomax @@ -79,6 +81,7 @@ function Equamax(p::Integer; miniter::Integer = 20, ϵ::Real = 1e-12) Orthomax(γ = float(p) / 2.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end + """ Parsimax(n, p) -> Orthomax diff --git a/test/facrot.jl b/test/facrot.jl index 56226f4..33dcae9 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -1,5 +1,4 @@ using MultivariateStats -using LinearAlgebra using Test import Random @@ -11,6 +10,8 @@ import Random ## Varimax rotation # Comparison with R's `varimax` function called as + # using Printf + # Base.show(io::IO, f::Float64) = @printf(io, "%1.8f", f) # using RCall # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" # R"F <- qm$loadings" @@ -26,7 +27,7 @@ import Random R = [ 0.96256370 0.27105556; -0.27105556 0.96256370] - FR = fit(FactorRotation, X) + FR = fit(FactorRotation, X, alg = Varimax()) @test isapprox(FR.F, F, rtol = √eps()) @test isapprox(FR.R, R, rtol = √eps()) @@ -46,30 +47,51 @@ import Random R = [ 0.94810241 0.31796513; -0.31796513 0.94810241] - FR = fit(FactorRotation, X, alg = Orthomax(γ = 0.0)) + FR = fit(FactorRotation, X, alg = Quartimax()) @test isapprox(FR.F, F, rtol = 1e-6) @test isapprox(FR.R, R, rtol = 1e-6) - ## Equamax - X = randn(10, 5) - - # Comparison with Matlab's `rotatefactors` called as - # using MATLAB - # mat"rotatefactors($X, 'Method', 'equamax')" - - F = [-0.01195028 0.05771308 -0.35271260 2.83724330 -0.68373015; - 0.07085630 2.43064454 -0.56171788 0.00250704 0.05441219; - -0.12946985 0.78708715 -0.83039068 -2.01918172 0.65545648; - 1.95142102 -1.08251779 -0.49721317 1.13050103 -1.23799305; - -0.55556677 -0.16288095 1.56941619 -0.36283530 2.46150446; - 0.21075154 2.16186107 1.16767716 -1.03674456 0.71948514; - 0.09579139 0.15765568 1.21512619 0.03140918 0.17671690; - 0.00855404 0.96679333 -1.69727390 -0.20320960 -0.75931306; - 0.01368675 0.34330631 0.00514439 -0.80566209 1.21579113; - -2.13711263 -0.36494414 -0.25451404 0.25001421 -0.07785436] - - # Equamax for n x p matrix means γ = p / 2 - FR = fit(FactorRotation, X, alg = Orthomax(γ = 5 / 2, maxiter = 20000)) - @test isapprox(FR.F, F, rtol = √eps()) + # ## Equamax + # X = randn(10, 5) + + # # Comparison with Matlab's `rotatefactors` called as + # # using MATLAB + # # mat"rotatefactors($X, 'Method', 'equamax', 'Normalize', 'off')" + + # F = [ 0.20295895 2.07838653 -0.54122597 -0.12487534 0.98358184; + # -1.03815211 0.17239352 -1.27261240 -0.41485812 1.44036802; + # -1.32942778 0.75883602 0.53246754 -0.39153390 2.81316627; + # 0.17943030 0.10200669 -0.02036786 -2.05727725 0.36179405; + # -0.57150075 0.83976222 0.13897091 -0.58537812 1.78231347; + # -0.24047953 -0.14609212 -0.10808767 -0.31099008 0.53277677; + # 1.63020394 -0.61318538 -0.77879612 -0.44644772 -0.30331055; + # 2.92004862 0.09873379 -0.22817883 -0.00487858 -1.05015088; + # 0.36654224 0.21341754 -2.75662510 -0.02498979 -0.47774001; + # -0.28002665 -0.92746802 -0.32811579 -0.52183972 0.15056833] + + # # Equamax for n x p matrix means γ = p / 2 + # FR = fit(FactorRotation, X, alg = Equamax(size(X, 2))) + # @test isapprox(FR.F, F, rtol = √eps()) + + # ## Parsimax + + # # Comparison with Matlab's `rotatefactors` called as + # # using MATLAB + # # mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')" + + # F = [ 0.20088755 2.13833957 -0.51952937 -0.15258917 0.85486221; + # -1.04430793 0.26969187 -1.25045141 -0.47113722 1.42298177; + # -1.32771913 0.92728472 0.57264722 -0.48321345 2.74040511; + # 0.17856831 0.12930037 -0.00621592 -2.06837501 0.28192375; + # -0.57130614 0.94788128 0.16710829 -0.64273775 1.70426783; + # -0.24126103 -0.11168423 -0.10096713 -0.33074538 0.53024057; + # 1.62645119 -0.62482282 -0.79024454 -0.44334505 -0.27314675; + # 2.91932555 0.03504886 -0.25094949 0.03018747 -1.05060055; + # 0.35443644 0.20160255 -2.76135234 -0.02259016 -0.46464651; + # -0.28200304 -0.91296776 -0.32805047 -0.53293768 0.19126898] + + # # Equamax for n x p matrix means γ = p / 2 + # FR = fit(FactorRotation, X, alg = Parsimax(size(X)..., maxiter = 2000, ϵ = 1e-5)) + # @test isapprox(FR.F, F, rtol = √eps()) end \ No newline at end of file From e51ede172a26c537530fe879071dacaa92bff047 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 14 Oct 2021 10:17:32 +0200 Subject: [PATCH 11/32] added vscode --- .vscode/settings.json | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..e1e4146 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "julia.environmentPath": "/Users/hefelix/.julia/environments/dev" +} \ No newline at end of file From 9327d9baab37e26636f8c3265c35190f66587d98 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Tue, 21 Dec 2021 15:20:30 +0100 Subject: [PATCH 12/32] Remove vscode --- .vscode/settings.json | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index e1e4146..0000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,3 +0,0 @@ -{ - "julia.environmentPath": "/Users/hefelix/.julia/environments/dev" -} \ No newline at end of file From cbb0687bf0fae33521534726e7e2c9c872bfe1b6 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 2 Mar 2022 17:10:56 +0100 Subject: [PATCH 13/32] Add orthomax factor rotations --- src/MultivariateStats.jl | 39 +++++---- src/facrot.jl | 179 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 203 insertions(+), 15 deletions(-) create mode 100644 src/facrot.jl diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 58943f8..a37450a 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -40,29 +40,29 @@ module MultivariateStats invsqrtm, # Compute C^{-1/2}, i.e. inv(sqrtm(C)) ## pca - PCA, # Type: Principal Component Analysis model + PCA, # Type: Principal Component Analysis model - pcacov, # PCA based on covariance - pcasvd, # PCA based on singular value decomposition of input data - principalratio, # the ratio of variances preserved in the principal subspace - principalvar, # the variance along a specific principal direction - principalvars, # the variances along all principal directions + pcacov, # PCA based on covariance + pcasvd, # PCA based on singular value decomposition of input data + principalratio, # the ratio of variances preserved in the principal subspace + principalvar, # the variance along a specific principal direction + principalvars, # the variances along all principal directions tprincipalvar, # total principal variance, i.e. sum(principalvars(M)) tresidualvar, # total residual variance ## ppca - PPCA, # Type: the Probabilistic PCA model + PPCA, # Type: the Probabilistic PCA model ppcaml, # Maximum likelihood probabilistic PCA ppcaem, # EM algorithm for probabilistic PCA bayespca, # Bayesian PCA ## kpca - KernelPCA, # Type: the Kernel PCA model + KernelPCA, # Type: the Kernel PCA model ## cca - CCA, # Type: Correlation Component Analysis model + CCA, # Type: Correlation Component Analysis model ccacov, # CCA based on covariances ccasvd, # CCA based on singular value decomposition of input data @@ -73,8 +73,8 @@ module MultivariateStats classical_mds, # perform classical MDS over a given distance matrix stress, # stress evaluation - gram2dmat, gram2dmat!, # Gram matrix => Distance matrix - dmat2gram, dmat2gram!, # Distance matrix => Gram matrix + gram2dmat, gram2dmat!, # Gram matrix => Distance matrix + dmat2gram, dmat2gram!, # Distance matrix => Gram matrix ## lda LinearDiscriminant, # Type: Linear Discriminant functional @@ -93,15 +93,23 @@ module MultivariateStats mclda_solve, # solve multi-class LDA projection given sStatisticalModel ## ica - ICA, # Type: the Fast ICA model + ICA, # Type: the Fast ICA model fastica!, # core algorithm function for the Fast ICA ## fa - FactorAnalysis, # Type: the Factor Analysis model + FactorAnalysis, # Type: the Factor Analysis model - faem, # EM algorithm for factor analysis - facm # CM algorithm for factor analysis + faem, # Maximum likelihood probabilistic PCA + facm, # EM algorithm for probabilistic PCA + + ## facrot + FactorRotationAlgorithm, # Type: Factor rotation algorithm + Orthomax, # Type: Orthomax factor rotation algorithm + + FactorRotation, # Type: Return type for factor rotations + + rotatefactors # Alternative interface to factor rotations ## source files include("types.jl") @@ -116,6 +124,7 @@ module MultivariateStats include("lda.jl") include("ica.jl") include("fa.jl") + include("facrot.jl") ## deprecations @deprecate indim(f) size(f)[1] diff --git a/src/facrot.jl b/src/facrot.jl new file mode 100644 index 0000000..8c07fc2 --- /dev/null +++ b/src/facrot.jl @@ -0,0 +1,179 @@ +""" + FactorRotationAlgorithm + +An abstract type for factor rotation algorithms. +""" +abstract type FactorRotationAlgorithm end + +""" + Orthomax <: FactorRotationAlgorithm + +A type representing the orthomax factor rotation algorithm. + +The positive parameter `γ::Real` determines which type of rotation is performed. +For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` +is quartimax rotation, `γ = p / 2` is equamax rotation, and +`γ = n (p - 1) / (n + p - 2)` is parsimax rotation. + +The parameter `maxiter::Integer` controls the maximum number of iterations to +perform (default `1000`) and `ϵ::Real` is a small positive constant determining +convergence (default `1e-12`). +""" +struct Orthomax <: FactorRotationAlgorithm + γ::Real + miniter::Integer + maxiter::Integer + ϵ::Real + + Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin + maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) + miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) + γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) + ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) + + new(γ, miniter, maxiter, ϵ) + end +end + +function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) + summary(io, alg); println(io) + println(io, "γ = $(alg.γ)") + println(io, "miniter = $(alg.miniter)") + println(io, "maxiter = $(alg.maxiter)") + println(io, "ϵ = $(alg.ϵ)") +end + +""" + FactorRotation{T <: Real} + +The return type for factor rotations. + +`F` contains the rotated factors and `R` is the rotation matrix that +was applied to the original factors. +""" +struct FactorRotation{T <: Real, Ta <: FactorRotationAlgorithm} + F::Matrix{T} + R::Matrix{T} + + alg::Ta + + FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) +end + +eltype(::Type{FactorRotation{T, Ta}}) where {T,Ta} = T + +FactorRotation(F::Matrix{T}, R::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation{T, Ta}(F, R, alg) +FactorRotation(F::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation(F, Matrix{eltype(F)}(I, size(F, 2), size(F, 2)), alg) + +function FactorRotation(F::Matrix, R::Matrix, alg::Ta) where {Ta <: FactorRotationAlgorithm} + return FactorRotation(promote(F, R)..., alg) +end + +function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any,<:Any}) + summary(io, FR); println(io) +end + + +## Comparison to other implementations +# The implementation of varimax in R row-normlises the input matrix before +# application of the algorithm and rescales the rows afterwards. +## Reference +# 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 +function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) + n, p = size(F) + if n < 2 + return (F, Matrix{eltype(F)}(I, p, p)) + end + + # Warm up step + # Do one step. If that first step did not lead away from the identity + # matrix enough use a random orthogonal matrix as a starting point. + M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) + R = M.U * M.Vt + if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ + R, _ = qr(randn(p, p)).Q + end + + # Main iterations + d = 0 + converged = false + for i in 1:maxiter + dold = d + B = F * R + M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) + R = M.U * M.Vt + d = sum(M.S) + if abs(d - dold) / d < ϵ && i >= miniter + converged = true + break + end + end + + if !converged + @warn "orthomax did not converge in $(maxiter) iterations" + end + + (F * R, R) +end + +""" + fit(::Type{FactorRotation}, F::AbstractMatrix; ...) -> FactorRotation + +Fit a factor rotation to the matrix `F` and apply it. The algorithm used to +perform the factor rotation is by default [`Orthomax`](@ref) and can be changed +with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. +""" +function fit(::Type{FactorRotation}, F::AbstractMatrix; + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + if isa(alg, Orthomax) + F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + return FactorRotation(F, R, alg) + end +end + +""" + fit(::Type{FactorRotation}, F::FactorAnalysis; ...) -> FactorRotation + +Fit a factor rotation to the loading matrix of the [`FactorAnalysis`](@ref) +object and apply it. + +See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. +""" +function fit(::Type{FactorRotation}, F::FactorAnalysis; + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + return fit(FactorRotation, F.W; alg) +end + +## Alternative interface + +""" + rotatefactors(F::AbstractMatrix, [alg::FactorRotationAlgorithm]) -> FactorRotation + +Rotate the factors in the matrix `F`. The algorithm to be used can be passed in +via the second argument `alg`. By default [`Orthomax`](@ref) is used. +""" +function rotatefactors(F::AbstractMatrix, alg::FactorRotationAlgorithm = Orthomax()) + F, R = _rotatefactors(F, alg) + return FactorRotation(F, R, alg) +end + +# Use multiple dispatch to decide on the algorithm implementation + +function _rotatefactors(F::AbstractMatrix, alg::Orthomax) + return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) +end + +""" + rotatefactors(F::FactorAnalysis, [alg::FactorRotationAlgorithm]) -> FactorAnalysis + +Rotate the factors in the loading matrix of `F` which is of type [`FactorAnalysis`](@ref). +The algorithm to be used can be passed in via the second argument `alg`. +By default [`Orthomax`](@ref) is used. +""" +function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) + FR = rotatefactors(F.W, alg) + return FactorAnalysis(F.mean, FR.F, F.Ψ) +end From c266980f2f0d0443f69187ced48a24ddfdec1c63 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 16:09:48 +0200 Subject: [PATCH 14/32] Make backwards-compatible with Julia 1.0 --- src/facrot.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/facrot.jl b/src/facrot.jl index 8c07fc2..4d2e824 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -144,7 +144,7 @@ See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword argumen """ function fit(::Type{FactorRotation}, F::FactorAnalysis; alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - return fit(FactorRotation, F.W; alg) + return fit(FactorRotation, F.W; alg = alg) end ## Alternative interface From 764512ea1b1492e87515fcaa5cb13107702ced85 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 16:19:16 +0200 Subject: [PATCH 15/32] Add documentation for `miniter` --- src/facrot.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index 4d2e824..2e77281 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -16,8 +16,9 @@ is quartimax rotation, `γ = p / 2` is equamax rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. The parameter `maxiter::Integer` controls the maximum number of iterations to -perform (default `1000`) and `ϵ::Real` is a small positive constant determining -convergence (default `1e-12`). +perform (default `1000`), `miniter::Integer` controls the minimum number of +iterations taken before convergence is checked, and `ϵ::Real` is a small positive +constant determining convergence (default `1e-12`). """ struct Orthomax <: FactorRotationAlgorithm γ::Real From c3c8bd9cc6bd874a8555f7da9ddecbcd46c8a830 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 16:48:26 +0200 Subject: [PATCH 16/32] Fixed small mistake --- src/facrot.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/facrot.jl b/src/facrot.jl index 2e77281..cb4d392 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -95,7 +95,7 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) R = M.U * M.Vt if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ - R, _ = qr(randn(p, p)).Q + R = qr(randn(p, p)).Q end # Main iterations From 3dcfb226a1e6bce8d0cb675c35ec1338e42a15ae Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 21 Oct 2020 21:48:03 +0200 Subject: [PATCH 17/32] Adjusted indent --- src/facrot.jl | 132 +++++++++++++++++++++++++------------------------- 1 file changed, 66 insertions(+), 66 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index cb4d392..1b7ee8a 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -21,27 +21,27 @@ iterations taken before convergence is checked, and `ϵ::Real` is a small positi constant determining convergence (default `1e-12`). """ struct Orthomax <: FactorRotationAlgorithm - γ::Real - miniter::Integer - maxiter::Integer - ϵ::Real - - Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin - maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) - miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) - γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) - ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) - - new(γ, miniter, maxiter, ϵ) - end + γ::Real + miniter::Integer + maxiter::Integer + ϵ::Real + + Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin + maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) + miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) + γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) + ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) + + new(γ, miniter, maxiter, ϵ) + end end function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) - summary(io, alg); println(io) - println(io, "γ = $(alg.γ)") - println(io, "miniter = $(alg.miniter)") - println(io, "maxiter = $(alg.maxiter)") - println(io, "ϵ = $(alg.ϵ)") + summary(io, alg); println(io) + println(io, "γ = $(alg.γ)") + println(io, "miniter = $(alg.miniter)") + println(io, "maxiter = $(alg.maxiter)") + println(io, "ϵ = $(alg.ϵ)") end """ @@ -53,12 +53,12 @@ The return type for factor rotations. was applied to the original factors. """ struct FactorRotation{T <: Real, Ta <: FactorRotationAlgorithm} - F::Matrix{T} - R::Matrix{T} + F::Matrix{T} + R::Matrix{T} - alg::Ta + alg::Ta - FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) + FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) end eltype(::Type{FactorRotation{T, Ta}}) where {T,Ta} = T @@ -67,11 +67,11 @@ FactorRotation(F::Matrix{T}, R::Matrix{T}, alg::Ta) where {T <: Real, Ta <: Fact FactorRotation(F::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation(F, Matrix{eltype(F)}(I, size(F, 2), size(F, 2)), alg) function FactorRotation(F::Matrix, R::Matrix, alg::Ta) where {Ta <: FactorRotationAlgorithm} - return FactorRotation(promote(F, R)..., alg) + return FactorRotation(promote(F, R)..., alg) end function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any,<:Any}) - summary(io, FR); println(io) + summary(io, FR); println(io) end @@ -84,40 +84,40 @@ end # Proc. SPIE 6144, Medical Imaging 2006: Image Processing, 61441G # (10 March 2006); doi: 10.1117/12.651293 function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) - n, p = size(F) - if n < 2 - return (F, Matrix{eltype(F)}(I, p, p)) - end - - # Warm up step - # Do one step. If that first step did not lead away from the identity - # matrix enough use a random orthogonal matrix as a starting point. - M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) - R = M.U * M.Vt - if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ - R = qr(randn(p, p)).Q - end - - # Main iterations - d = 0 - converged = false - for i in 1:maxiter - dold = d - B = F * R - M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) + n, p = size(F) + if n < 2 + return (F, Matrix{eltype(F)}(I, p, p)) + end + + # Warm up step + # Do one step. If that first step did not lead away from the identity + # matrix enough use a random orthogonal matrix as a starting point. + M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) R = M.U * M.Vt - d = sum(M.S) - if abs(d - dold) / d < ϵ && i >= miniter - converged = true - break + if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ + R = qr(randn(p, p)).Q end - end - if !converged - @warn "orthomax did not converge in $(maxiter) iterations" - end + # Main iterations + d = 0 + converged = false + for i in 1:maxiter + dold = d + B = F * R + M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) + R = M.U * M.Vt + d = sum(M.S) + if abs(d - dold) / d < ϵ && i >= miniter + converged = true + break + end + end - (F * R, R) + if !converged + @warn "orthomax did not converge in $(maxiter) iterations" + end + + (F * R, R) end """ @@ -128,11 +128,11 @@ perform the factor rotation is by default [`Orthomax`](@ref) and can be changed with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. """ function fit(::Type{FactorRotation}, F::AbstractMatrix; - alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - if isa(alg, Orthomax) - F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) - return FactorRotation(F, R, alg) - end + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + if isa(alg, Orthomax) + F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + return FactorRotation(F, R, alg) + end end """ @@ -144,8 +144,8 @@ object and apply it. See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. """ function fit(::Type{FactorRotation}, F::FactorAnalysis; - alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - return fit(FactorRotation, F.W; alg = alg) + alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} + return fit(FactorRotation, F.W; alg = alg) end ## Alternative interface @@ -157,14 +157,14 @@ Rotate the factors in the matrix `F`. The algorithm to be used can be passed in via the second argument `alg`. By default [`Orthomax`](@ref) is used. """ function rotatefactors(F::AbstractMatrix, alg::FactorRotationAlgorithm = Orthomax()) - F, R = _rotatefactors(F, alg) - return FactorRotation(F, R, alg) + F, R = _rotatefactors(F, alg) + return FactorRotation(F, R, alg) end # Use multiple dispatch to decide on the algorithm implementation function _rotatefactors(F::AbstractMatrix, alg::Orthomax) - return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) end """ @@ -175,6 +175,6 @@ The algorithm to be used can be passed in via the second argument `alg`. By default [`Orthomax`](@ref) is used. """ function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) - FR = rotatefactors(F.W, alg) - return FactorAnalysis(F.mean, FR.F, F.Ψ) + FR = rotatefactors(F.W, alg) + return FactorAnalysis(F.mean, FR.F, F.Ψ) end From d00ceb5c95b9193e2e2f798f73719b842b6fda73 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 00:22:15 +0200 Subject: [PATCH 18/32] Better exception use/getters for loadings/rotation --- src/MultivariateStats.jl | 4 +++- src/facrot.jl | 45 ++++++++++++++++++++++++++-------------- 2 files changed, 32 insertions(+), 17 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index a37450a..2d72969 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -108,7 +108,9 @@ module MultivariateStats Orthomax, # Type: Orthomax factor rotation algorithm FactorRotation, # Type: Return type for factor rotations - + loadings, # rotated loading matrix + rotation, # rotation matrix applied to loadings + rotatefactors # Alternative interface to factor rotations ## source files diff --git a/src/facrot.jl b/src/facrot.jl index 1b7ee8a..7117c1c 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -11,8 +11,8 @@ abstract type FactorRotationAlgorithm end A type representing the orthomax factor rotation algorithm. The positive parameter `γ::Real` determines which type of rotation is performed. -For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` -is quartimax rotation, `γ = p / 2` is equamax rotation, and +For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` +is quartimax rotation, `γ = p / 2` is equamax rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. The parameter `maxiter::Integer` controls the maximum number of iterations to @@ -27,10 +27,10 @@ struct Orthomax <: FactorRotationAlgorithm ϵ::Real Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin - maxiter > zero(eltype(maxiter)) || throw(ArgumentError("Orthomax: maxiter needs to be positive")) - miniter > zero(eltype(miniter)) || throw(ArgumentError("Orthomax: miniter needs to be positive")) - γ > zero(eltype(γ)) || throw(ArgumentError("Orthomax: γ needs to be positive")) - ϵ > zero(eltype(ϵ)) || throw(ArgumentError("Orthomax: ϵ needs to be positive")) + γ ≥ zero(eltype(γ)) || throw(DomainError("Orthomax: γ needs to be non-negative")) + miniter > zero(eltype(miniter)) || throw(DomainError("Orthomax: miniter needs to be positive")) + maxiter > zero(eltype(maxiter)) || throw(DomainError("Orthomax: maxiter needs to be positive")) + ϵ > zero(eltype(ϵ)) || throw(DomainError("Orthomax: ϵ needs to be positive")) new(γ, miniter, maxiter, ϵ) end @@ -74,13 +74,26 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any summary(io, FR); println(io) end +""" + loadings(FR::FactorRotation) + +Returns the loading matrix. +""" +loadings(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.F + +""" + rotation(FR::FactorRotation) + +Returns the rotation matrix. +""" +rotation(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.R ## Comparison to other implementations # The implementation of varimax in R row-normlises the input matrix before # application of the algorithm and rescales the rows afterwards. -## Reference -# Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen, "Sparse modeling of -# landmark and texture variability using the orthomax criterion," +## Reference +# 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 function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) @@ -100,6 +113,7 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) # Main iterations d = 0 + lastchange = NaN converged = false for i in 1:maxiter dold = d @@ -107,15 +121,14 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) R = M.U * M.Vt d = sum(M.S) - if abs(d - dold) / d < ϵ && i >= miniter + lastchange = abs(d - dold) / d + if lastchange < ϵ && i >= miniter converged = true break end end - if !converged - @warn "orthomax did not converge in $(maxiter) iterations" - end + converged || throw(ConvergenceException(maxiter, lastchange, ϵ)) (F * R, R) end @@ -127,7 +140,7 @@ Fit a factor rotation to the matrix `F` and apply it. The algorithm used to perform the factor rotation is by default [`Orthomax`](@ref) and can be changed with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. """ -function fit(::Type{FactorRotation}, F::AbstractMatrix; +function fit(::Type{FactorRotation}, F::AbstractMatrix; alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} if isa(alg, Orthomax) F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) @@ -143,7 +156,7 @@ object and apply it. See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. """ -function fit(::Type{FactorRotation}, F::FactorAnalysis; +function fit(::Type{FactorRotation}, F::FactorAnalysis; alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} return fit(FactorRotation, F.W; alg = alg) end @@ -171,7 +184,7 @@ end rotatefactors(F::FactorAnalysis, [alg::FactorRotationAlgorithm]) -> FactorAnalysis Rotate the factors in the loading matrix of `F` which is of type [`FactorAnalysis`](@ref). -The algorithm to be used can be passed in via the second argument `alg`. +The algorithm to be used can be passed in via the second argument `alg`. By default [`Orthomax`](@ref) is used. """ function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) From 5f203921b7131b8282440ce05d768ff415e7b3c9 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Wed, 2 Mar 2022 17:12:16 +0100 Subject: [PATCH 19/32] Created some tests --- test/facrot.jl | 74 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 4 +-- 2 files changed, 76 insertions(+), 2 deletions(-) create mode 100644 test/facrot.jl diff --git a/test/facrot.jl b/test/facrot.jl new file mode 100644 index 0000000..501b3dd --- /dev/null +++ b/test/facrot.jl @@ -0,0 +1,74 @@ +using MultivariateStats +using LinearAlgebra +using Test +import Random + +@testset "Factor Analysis" begin + + Random.seed!(37273) + + X = randn(5, 2) + + ## Varimax rotation + # Comparison with R's `varimax` function called as (after `using RCall`) + # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" + # R"F <- qm$loadings" + # R"R <- qm$rotmat" + # @rget F + # @rget R + + F = [-2.56254778 -0.05712524; + -0.37387139 0.53683094; + -1.76868624 0.67513474; + 0.08057810 0.34826689; + -0.06196338 0.20433233] + R = [ 0.96256370 0.27105556; + -0.27105556 0.96256370] + + FR = fit(FactorRotation, X) + @test isapprox(FR.F, F, rtol = √eps()) + @test isapprox(FR.R, R, rtol = √eps()) + + ## Quartimax rotation + # Comparison with R's `GPArotation::quartimax` function called as + # R"qm <- GPArotation::quartimax($X, eps = 1e-6, maxit = 10000)" + # R"F <- qm$loadings" + # R"R <- qm$Th" + # @rget F + # @rget R + + F = [-2.55665706 -0.18280898; + -0.39976499 0.51783707; + -1.79968636 0.58752613; + 0.06339043 0.35180152; + -0.07191598 0.20104540] + R = [ 0.94810241 0.31796513; + -0.31796513 0.94810241] + + FR = fit(FactorRotation, X, alg = Orthomax(γ = 0.0)) + @test isapprox(FR.F, F, rtol = 1e-6) + @test isapprox(FR.R, R, rtol = 1e-6) + + ## Equamax + X = randn(10, 5) + + # Comparsion with Matlab's `rotatefactors` + # using MATLAB + # mat"rotatefactors($X, 'Method', 'equamax')" + + F = [-0.01195028 0.05771308 -0.35271260 2.83724330 -0.68373015; + 0.07085630 2.43064454 -0.56171788 0.00250704 0.05441219; + -0.12946985 0.78708715 -0.83039068 -2.01918172 0.65545648; + 1.95142102 -1.08251779 -0.49721317 1.13050103 -1.23799305; + -0.55556677 -0.16288095 1.56941619 -0.36283530 2.46150446; + 0.21075154 2.16186107 1.16767716 -1.03674456 0.71948514; + 0.09579139 0.15765568 1.21512619 0.03140918 0.17671690; + 0.00855404 0.96679333 -1.69727390 -0.20320960 -0.75931306; + 0.01368675 0.34330631 0.00514439 -0.80566209 1.21579113; + -2.13711263 -0.36494414 -0.25451404 0.25001421 -0.07785436] + + # Equamax for n x p matrix means γ = p / 2 + FR = fit(FactorRotation, X, alg = Orthomax(γ = 5 / 2, maxiter = 20000)) + @test isapprox(FR.F, F, rtol = √eps()) + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index a551076..fc4a0fb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,8 +8,8 @@ tests = ["lreg", "ica", "ppca", "kpca", - "fa" - ] + "fa", + "facrot"] for test in tests include(test*".jl") From 975782f157811399fa9856212c739f680efa37f2 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 10:24:30 +0200 Subject: [PATCH 20/32] Convenience functions for common rotations --- src/facrot.jl | 87 ++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 72 insertions(+), 15 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index 7117c1c..f327b7a 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -10,15 +10,15 @@ abstract type FactorRotationAlgorithm end A type representing the orthomax factor rotation algorithm. -The positive parameter `γ::Real` determines which type of rotation is performed. -For a `n x p` matrix, the default `γ = 1.0` leads to varimax rotation, `γ = 0.0` -is quartimax rotation, `γ = p / 2` is equamax rotation, and -`γ = n (p - 1) / (n + p - 2)` is parsimax rotation. +The positive parameter `γ::Real` determines which type of rotation is +performed. For a `n x p` matrix, the default `γ = 1.0` leads to varimax +rotation, `γ = 0.0` is quartimax rotation, `γ = p / 2` is equamax +rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. -The parameter `maxiter::Integer` controls the maximum number of iterations to -perform (default `1000`), `miniter::Integer` controls the minimum number of -iterations taken before convergence is checked, and `ϵ::Real` is a small positive -constant determining convergence (default `1e-12`). +The parameter `maxiter::Integer` controls the maximum number of iterations +to perform (default `1000`), `miniter::Integer` controls the minimum number +of iterations taken before convergence is checked, and `ϵ::Real` is a small +positive constant determining convergence (default `1e-12`). """ struct Orthomax <: FactorRotationAlgorithm γ::Real @@ -26,16 +26,73 @@ struct Orthomax <: FactorRotationAlgorithm maxiter::Integer ϵ::Real - Orthomax(;γ = 1.0, miniter = 20, maxiter = 1000, ϵ = 1e-12) = begin - γ ≥ zero(eltype(γ)) || throw(DomainError("Orthomax: γ needs to be non-negative")) - miniter > zero(eltype(miniter)) || throw(DomainError("Orthomax: miniter needs to be positive")) - maxiter > zero(eltype(maxiter)) || throw(DomainError("Orthomax: maxiter needs to be positive")) - ϵ > zero(eltype(ϵ)) || throw(DomainError("Orthomax: ϵ needs to be positive")) - - new(γ, miniter, maxiter, ϵ) + Orthomax(;γ::Union{Real, Integer} = 1.0, + miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) = begin + γ ≥ zero(eltype(γ)) || + throw(DomainError("Orthomax: γ needs to be non-negative")) + miniter > zero(eltype(miniter)) || + throw(DomainError("Orthomax: miniter needs to be positive")) + maxiter > zero(eltype(maxiter)) || + throw(DomainError("Orthomax: maxiter needs to be positive")) + ϵ > zero(eltype(ϵ)) || + throw(DomainError("Orthomax: ϵ needs to be positive")) + + new(float(γ), miniter, maxiter, ϵ) end end +""" + Varimax() -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = 1.0`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Varimax(;miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end +""" + Quartimax() -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = 0.0`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Quartimax(;miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = 0.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end +""" + Equamax(p) -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = p / 2`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Equamax(p::Integer; miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = float(p) / 2.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end +""" + Parsimax(n, p) -> Orthomax + +Creates an orthomax factor rotation algorithm object for a matrix of +size `n x p` with `γ = n (p - 1) / (n + p - 2)`. +Remaining keyword parameters as for [`Orthomax`](@ref). +""" +function Parsimax(n::Integer, p::Integer; miniter::Integer = 20, + maxiter::Integer = 1000, + ϵ::Real = 1e-12) + Orthomax(γ = float(n) * (float(p) - 1) / (float(n) + float(p) - 2), + miniter = miniter, maxiter = maxiter, ϵ = ϵ) +end + function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) summary(io, alg); println(io) println(io, "γ = $(alg.γ)") From 3c5cd0c3e23d3437242cdcbfbf7f4fdfe0ed527c Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 10:24:48 +0200 Subject: [PATCH 21/32] Corrected name of test set for factor rotations --- test/facrot.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/facrot.jl b/test/facrot.jl index 501b3dd..56226f4 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -3,14 +3,15 @@ using LinearAlgebra using Test import Random -@testset "Factor Analysis" begin +@testset "Factor Rotations" begin Random.seed!(37273) X = randn(5, 2) ## Varimax rotation - # Comparison with R's `varimax` function called as (after `using RCall`) + # Comparison with R's `varimax` function called as + # using RCall # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" # R"F <- qm$loadings" # R"R <- qm$rotmat" @@ -52,7 +53,7 @@ import Random ## Equamax X = randn(10, 5) - # Comparsion with Matlab's `rotatefactors` + # Comparison with Matlab's `rotatefactors` called as # using MATLAB # mat"rotatefactors($X, 'Method', 'equamax')" From 0d848410ea2ce3761894ce70ec90d2a7b4aeccc0 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Thu, 22 Oct 2020 16:08:23 +0200 Subject: [PATCH 22/32] Convenience wrappers for common Orthomax rotations --- src/MultivariateStats.jl | 4 +++ src/facrot.jl | 3 ++ test/facrot.jl | 70 ++++++++++++++++++++++++++-------------- 3 files changed, 53 insertions(+), 24 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 2d72969..7aec743 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -106,6 +106,10 @@ module MultivariateStats ## facrot FactorRotationAlgorithm, # Type: Factor rotation algorithm Orthomax, # Type: Orthomax factor rotation algorithm + Varimax, # Convenience interfaces to Orthomax factor rotation algorithms + Quartimax, + Equamax, + Parsimax, FactorRotation, # Type: Return type for factor rotations loadings, # rotated loading matrix diff --git a/src/facrot.jl b/src/facrot.jl index f327b7a..235de40 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -55,6 +55,7 @@ function Varimax(;miniter::Integer = 20, ϵ::Real = 1e-12) Orthomax(γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end + """ Quartimax() -> Orthomax @@ -67,6 +68,7 @@ function Quartimax(;miniter::Integer = 20, ϵ::Real = 1e-12) Orthomax(γ = 0.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end + """ Equamax(p) -> Orthomax @@ -79,6 +81,7 @@ function Equamax(p::Integer; miniter::Integer = 20, ϵ::Real = 1e-12) Orthomax(γ = float(p) / 2.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end + """ Parsimax(n, p) -> Orthomax diff --git a/test/facrot.jl b/test/facrot.jl index 56226f4..33dcae9 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -1,5 +1,4 @@ using MultivariateStats -using LinearAlgebra using Test import Random @@ -11,6 +10,8 @@ import Random ## Varimax rotation # Comparison with R's `varimax` function called as + # using Printf + # Base.show(io::IO, f::Float64) = @printf(io, "%1.8f", f) # using RCall # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" # R"F <- qm$loadings" @@ -26,7 +27,7 @@ import Random R = [ 0.96256370 0.27105556; -0.27105556 0.96256370] - FR = fit(FactorRotation, X) + FR = fit(FactorRotation, X, alg = Varimax()) @test isapprox(FR.F, F, rtol = √eps()) @test isapprox(FR.R, R, rtol = √eps()) @@ -46,30 +47,51 @@ import Random R = [ 0.94810241 0.31796513; -0.31796513 0.94810241] - FR = fit(FactorRotation, X, alg = Orthomax(γ = 0.0)) + FR = fit(FactorRotation, X, alg = Quartimax()) @test isapprox(FR.F, F, rtol = 1e-6) @test isapprox(FR.R, R, rtol = 1e-6) - ## Equamax - X = randn(10, 5) - - # Comparison with Matlab's `rotatefactors` called as - # using MATLAB - # mat"rotatefactors($X, 'Method', 'equamax')" - - F = [-0.01195028 0.05771308 -0.35271260 2.83724330 -0.68373015; - 0.07085630 2.43064454 -0.56171788 0.00250704 0.05441219; - -0.12946985 0.78708715 -0.83039068 -2.01918172 0.65545648; - 1.95142102 -1.08251779 -0.49721317 1.13050103 -1.23799305; - -0.55556677 -0.16288095 1.56941619 -0.36283530 2.46150446; - 0.21075154 2.16186107 1.16767716 -1.03674456 0.71948514; - 0.09579139 0.15765568 1.21512619 0.03140918 0.17671690; - 0.00855404 0.96679333 -1.69727390 -0.20320960 -0.75931306; - 0.01368675 0.34330631 0.00514439 -0.80566209 1.21579113; - -2.13711263 -0.36494414 -0.25451404 0.25001421 -0.07785436] - - # Equamax for n x p matrix means γ = p / 2 - FR = fit(FactorRotation, X, alg = Orthomax(γ = 5 / 2, maxiter = 20000)) - @test isapprox(FR.F, F, rtol = √eps()) + # ## Equamax + # X = randn(10, 5) + + # # Comparison with Matlab's `rotatefactors` called as + # # using MATLAB + # # mat"rotatefactors($X, 'Method', 'equamax', 'Normalize', 'off')" + + # F = [ 0.20295895 2.07838653 -0.54122597 -0.12487534 0.98358184; + # -1.03815211 0.17239352 -1.27261240 -0.41485812 1.44036802; + # -1.32942778 0.75883602 0.53246754 -0.39153390 2.81316627; + # 0.17943030 0.10200669 -0.02036786 -2.05727725 0.36179405; + # -0.57150075 0.83976222 0.13897091 -0.58537812 1.78231347; + # -0.24047953 -0.14609212 -0.10808767 -0.31099008 0.53277677; + # 1.63020394 -0.61318538 -0.77879612 -0.44644772 -0.30331055; + # 2.92004862 0.09873379 -0.22817883 -0.00487858 -1.05015088; + # 0.36654224 0.21341754 -2.75662510 -0.02498979 -0.47774001; + # -0.28002665 -0.92746802 -0.32811579 -0.52183972 0.15056833] + + # # Equamax for n x p matrix means γ = p / 2 + # FR = fit(FactorRotation, X, alg = Equamax(size(X, 2))) + # @test isapprox(FR.F, F, rtol = √eps()) + + # ## Parsimax + + # # Comparison with Matlab's `rotatefactors` called as + # # using MATLAB + # # mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')" + + # F = [ 0.20088755 2.13833957 -0.51952937 -0.15258917 0.85486221; + # -1.04430793 0.26969187 -1.25045141 -0.47113722 1.42298177; + # -1.32771913 0.92728472 0.57264722 -0.48321345 2.74040511; + # 0.17856831 0.12930037 -0.00621592 -2.06837501 0.28192375; + # -0.57130614 0.94788128 0.16710829 -0.64273775 1.70426783; + # -0.24126103 -0.11168423 -0.10096713 -0.33074538 0.53024057; + # 1.62645119 -0.62482282 -0.79024454 -0.44334505 -0.27314675; + # 2.91932555 0.03504886 -0.25094949 0.03018747 -1.05060055; + # 0.35443644 0.20160255 -2.76135234 -0.02259016 -0.46464651; + # -0.28200304 -0.91296776 -0.32805047 -0.53293768 0.19126898] + + # # Equamax for n x p matrix means γ = p / 2 + # FR = fit(FactorRotation, X, alg = Parsimax(size(X)..., maxiter = 2000, ϵ = 1e-5)) + # @test isapprox(FR.F, F, rtol = √eps()) end \ No newline at end of file From f468367148defa07674624f87dc1d6c4de82804d Mon Sep 17 00:00:00 2001 From: Felix Held Date: Mon, 4 Jul 2022 15:33:00 +0200 Subject: [PATCH 23/32] Commit to merge upstream changes --- src/MultivariateStats.jl | 9 ++- src/facrot.jl | 2 +- test/facrot.jl | 132 +++++++++++++++++++-------------------- 3 files changed, 73 insertions(+), 70 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 2cae1d2..c3c148b 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -3,7 +3,7 @@ module MultivariateStats using StatsBase: SimpleCovariance, CovarianceEstimator, RegressionModel, AbstractDataTransform, pairwise! import Statistics: mean, var, cov, covm, cor - import Base: length, size, show, dump + import Base: length, size, show, dump, eltype import StatsBase: fit, predict, predict!, ConvergenceException, coef, weights, dof, pairwise, r2 import SparseArrays @@ -101,11 +101,16 @@ module MultivariateStats FactorAnalysis, # Type: the Factor Analysis model faem, # EM algorithm for factor analysis - facm # CM algorithm for factor analysis + facm, # CM algorithm for factor analysis ## facrot FactorRotationAlgorithm, # Type: Factor rotation algorithm Orthomax, # Type: Orthomax factor rotation algorithm + # Convenience types for factor rotation algorithms + Varimax, # Type + Quartimax, # Type + Equamax, # Type + Parsimax, # Type FactorRotation, # Type: Return type for factor rotations diff --git a/src/facrot.jl b/src/facrot.jl index 235de40..4400c10 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -53,7 +53,7 @@ Remaining keyword parameters as for [`Orthomax`](@ref). function Varimax(;miniter::Integer = 20, maxiter::Integer = 1000, ϵ::Real = 1e-12) - Orthomax(γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) + Orthomax(; γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) end """ diff --git a/test/facrot.jl b/test/facrot.jl index 33dcae9..f8d0c7c 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -1,97 +1,95 @@ using MultivariateStats using Test -import Random +using StableRNGs @testset "Factor Rotations" begin - Random.seed!(37273) + rng = StableRNG(37273) - X = randn(5, 2) + X = randn(rng, 5, 2) ## Varimax rotation - # Comparison with R's `varimax` function called as + # Ground-truth was extracted by specifying the following + # print command + # ```` # using Printf # Base.show(io::IO, f::Float64) = @printf(io, "%1.8f", f) + # ```` + # Comparison with R's `varimax` function called as + # ```` # using RCall # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)" - # R"F <- qm$loadings" - # R"R <- qm$rotmat" + # R"F <- vm$loadings" + # R"R <- vm$rotmat" # @rget F # @rget R + # ```` - F = [-2.56254778 -0.05712524; - -0.37387139 0.53683094; - -1.76868624 0.67513474; - 0.08057810 0.34826689; - -0.06196338 0.20433233] - R = [ 0.96256370 0.27105556; - -0.27105556 0.96256370] + F = [ 0.90503610 1.24332783; + -2.32641022 0.36067316; + 0.97744296 1.36280123; + 0.85650467 0.23818243; + 0.45865489 -2.40335769] + R = [ 0.89940646 0.43711327; + -0.43711327 0.89940646] FR = fit(FactorRotation, X, alg = Varimax()) - @test isapprox(FR.F, F, rtol = √eps()) - @test isapprox(FR.R, R, rtol = √eps()) + @test isapprox(FR.F, F, rtol = √eps(Float64)) + @test isapprox(FR.R, R, rtol = √eps(Float64)) ## Quartimax rotation # Comparison with R's `GPArotation::quartimax` function called as + # ``` + # using RCall # R"qm <- GPArotation::quartimax($X, eps = 1e-6, maxit = 10000)" # R"F <- qm$loadings" # R"R <- qm$Th" # @rget F # @rget R + # ``` - F = [-2.55665706 -0.18280898; - -0.39976499 0.51783707; - -1.79968636 0.58752613; - 0.06339043 0.35180152; - -0.07191598 0.20104540] - R = [ 0.94810241 0.31796513; - -0.31796513 0.94810241] + F = [ 0.89959033 1.24727370; + -2.32796522 0.35049625; + 0.97147404 1.36706259; + 0.85545490 0.24192567; + 0.46916046 -2.40132899] + R = [ 0.89748635 0.44104222; + -0.44104222 0.89748635] FR = fit(FactorRotation, X, alg = Quartimax()) - @test isapprox(FR.F, F, rtol = 1e-6) - @test isapprox(FR.R, R, rtol = 1e-6) - - # ## Equamax - # X = randn(10, 5) - - # # Comparison with Matlab's `rotatefactors` called as - # # using MATLAB - # # mat"rotatefactors($X, 'Method', 'equamax', 'Normalize', 'off')" - - # F = [ 0.20295895 2.07838653 -0.54122597 -0.12487534 0.98358184; - # -1.03815211 0.17239352 -1.27261240 -0.41485812 1.44036802; - # -1.32942778 0.75883602 0.53246754 -0.39153390 2.81316627; - # 0.17943030 0.10200669 -0.02036786 -2.05727725 0.36179405; - # -0.57150075 0.83976222 0.13897091 -0.58537812 1.78231347; - # -0.24047953 -0.14609212 -0.10808767 -0.31099008 0.53277677; - # 1.63020394 -0.61318538 -0.77879612 -0.44644772 -0.30331055; - # 2.92004862 0.09873379 -0.22817883 -0.00487858 -1.05015088; - # 0.36654224 0.21341754 -2.75662510 -0.02498979 -0.47774001; - # -0.28002665 -0.92746802 -0.32811579 -0.52183972 0.15056833] - - # # Equamax for n x p matrix means γ = p / 2 - # FR = fit(FactorRotation, X, alg = Equamax(size(X, 2))) - # @test isapprox(FR.F, F, rtol = √eps()) - - # ## Parsimax - - # # Comparison with Matlab's `rotatefactors` called as - # # using MATLAB - # # mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')" - - # F = [ 0.20088755 2.13833957 -0.51952937 -0.15258917 0.85486221; - # -1.04430793 0.26969187 -1.25045141 -0.47113722 1.42298177; - # -1.32771913 0.92728472 0.57264722 -0.48321345 2.74040511; - # 0.17856831 0.12930037 -0.00621592 -2.06837501 0.28192375; - # -0.57130614 0.94788128 0.16710829 -0.64273775 1.70426783; - # -0.24126103 -0.11168423 -0.10096713 -0.33074538 0.53024057; - # 1.62645119 -0.62482282 -0.79024454 -0.44334505 -0.27314675; - # 2.91932555 0.03504886 -0.25094949 0.03018747 -1.05060055; - # 0.35443644 0.20160255 -2.76135234 -0.02259016 -0.46464651; - # -0.28200304 -0.91296776 -0.32805047 -0.53293768 0.19126898] - - # # Equamax for n x p matrix means γ = p / 2 - # FR = fit(FactorRotation, X, alg = Parsimax(size(X)..., maxiter = 2000, ϵ = 1e-5)) - # @test isapprox(FR.F, F, rtol = √eps()) + @test isapprox(FR.F, F, rtol = √eps(Float64)) + @test isapprox(FR.R, R, rtol = √eps(Float64)) + + ## Equamax + # Comparison with Matlab's `rotatefactors` called as + # using MATLAB + # F = mat"rotatefactors($X, 'Method', 'equamax', 'Normalize', 'off')" + + F = [ 0.90503611 1.24332782; + -2.32641022 0.36067319; + 0.97744298 1.36280122; + 0.85650467 0.23818242; + 0.45865486 -2.40335769] + + # Equamax for n x p matrix means γ = p / 2 + FR = fit(FactorRotation, X, alg = Equamax(size(X, 2))) + @test isapprox(FR.F, F, rtol = √eps(Float64)) + + + ## Parsimax + + # Comparison with Matlab's `rotatefactors` called as + # using MATLAB + # F = mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')" + + F = [ 0.90503611 1.24332782; + -2.32641022 0.36067319; + 0.97744298 1.36280122; + 0.85650467 0.23818242; + 0.45865486 -2.40335769] + + # Equamax for n x p matrix means γ = p / 2 + FR = fit(FactorRotation, X, alg = Parsimax(size(X)...)) + @test isapprox(FR.F, F, rtol = √eps(Float64)) end \ No newline at end of file From 5f9524652c23867ac5afa35d46afd0d2c28bda94 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Mon, 4 Jul 2022 17:10:37 +0200 Subject: [PATCH 24/32] Translation of GPA ortho rotations from R code --- src/MultivariateStats.jl | 2 +- src/facrot.jl | 67 ++++++++++++++++++++++++++++++++++++++++ test/facrot.jl | 6 ++++ 3 files changed, 74 insertions(+), 1 deletion(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index e35055b..6f278df 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -8,7 +8,7 @@ module MultivariateStats ConvergenceException, pairwise, pairwise!, CoefTable import Statistics: mean, var, cov, covm, cor - import Base: length, size, show + import Base: length, size, show, eltype import StatsAPI: fit, predict, coef, weights, dof, r2 import LinearAlgebra: eigvals, eigvecs diff --git a/src/facrot.jl b/src/facrot.jl index 4400c10..e8dbd04 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -193,6 +193,73 @@ function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) (F * R, R) end +# Compute value of the optimization criterion and gradient of the +# underlying quality measure +function varimax(L::AbstractMatrix) + Q = L.^2 .- mean.(eachcol(L.^2))' + (-L .* Q, -abs(sum(diag(Q' * Q))) / 4) +end + +function quartimax(L::AbstractMatrix) + Q = L.^2 + (-L .* Q, -sum(diag(Q' * Q)) / 4) +end + +## Orthogonal rotations computed by gradient projection algorithms +# This algorithm computes orthogonal rotations using a gradient +# project algorithm. +## Reference +# Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms +# and Software for Arbitrary Rotation Criteria in Factor Analysis. +# Educational and Psychological Measurement, 65, 676–696 +# doi 10.1177/0013164404272507 +function gpaortho(F::AbstractMatrix, gradval, maxiter, lsiter, ϵ) + n, p = size(F) + if n < 2 + return (F, Matrix{eltype(F)}(I, p, p)) + end + + # Setup + R = Matrix{eltype(F)}(I, p, p) + α = 1.0 + L = F * R + + Gq, f = gradval(L) + G = F' * Gq + + Gqt, ft = gradval(L) + s = 0 + for _ in 1:maxiter + M = R' * G + S = (M + M') / 2 + Gp = G - R * S + s = sqrt(sum(diag(Gp' * Gp))) + if s < ϵ + break + end + α *= 2.0 + Rt = Matrix{eltype(F)}(I, p, p) + for _ in 1:lsiter + X = R - α * Gp + UDV = svd(X) + Rt = UDV.U * UDV.Vt + L = F * Rt + Gqt, ft = gradval(L) + if (ft < (f - 0.5 * s^2 * α)) + break + end + α /= 2.0 + end + R = Rt + f = ft + G = F' * Gqt + end + + (s < ϵ) || throw(ConvergenceException(maxiter, s, ϵ)) + + (L, R) +end + """ fit(::Type{FactorRotation}, F::AbstractMatrix; ...) -> FactorRotation diff --git a/test/facrot.jl b/test/facrot.jl index f8d0c7c..b619fee 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -37,6 +37,9 @@ using StableRNGs @test isapprox(FR.F, F, rtol = √eps(Float64)) @test isapprox(FR.R, R, rtol = √eps(Float64)) + isapprox(gpaortho(X, varimax, 1000, 10, 1e-6)[1], F, rtol = 1e-6) + isapprox(gpaortho(X, varimax, 1000, 10, 1e-6)[2], R, rtol = 1e-6) + ## Quartimax rotation # Comparison with R's `GPArotation::quartimax` function called as # ``` @@ -60,6 +63,9 @@ using StableRNGs @test isapprox(FR.F, F, rtol = √eps(Float64)) @test isapprox(FR.R, R, rtol = √eps(Float64)) + isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[1], F, rtol = 1e-6) + isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[2], R, rtol = 1e-6) + ## Equamax # Comparison with Matlab's `rotatefactors` called as # using MATLAB From 588be725fc51193603418001b3846963372bbba8 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Mon, 4 Jul 2022 18:26:23 +0200 Subject: [PATCH 25/32] First translation of GPA oblique rotations from R code --- src/facrot.jl | 65 ++++++++++++++++++++++++++++++++++++++++++++++++++ test/facrot.jl | 26 ++++++++++++++++++-- 2 files changed, 89 insertions(+), 2 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index e8dbd04..c5a0994 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -260,6 +260,71 @@ function gpaortho(F::AbstractMatrix, gradval, maxiter, lsiter, ϵ) (L, R) end +function quartimin(L::AbstractMatrix) + _, p = size(L) + Q = L.^2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) + (L .* Q, sum(L.^2 .* Q) / 4) +end + +## Oblique rotations computed by gradient projection algorithms +# This algorithm computes oblique rotations using a gradient +# project algorithm. +## Reference +# Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms +# and Software for Arbitrary Rotation Criteria in Factor Analysis. +# Educational and Psychological Measurement, 65, 676–696 +# doi 10.1177/0013164404272507 +function gpaoblique(F::AbstractMatrix, gradval, maxiter, lsiter, ϵ) + n, p = size(F) + if n < 2 + return (F, Matrix{eltype(F)}(I, p, p)) + end + + # Setup + R = Matrix{eltype(F)}(I, p, p) + α = 1.0 + # L = F * inv(R)' + L = (R \ F')' + # isapprox((R \ F')', F * inv(R)') + + Gq, f = gradval(L) + # G = -(L' * Gq * inv(R))' + G = R' \ (-Gq' * L) + # isapprox(R' \ (-Gq' * L), -(L' * Gq * inv(R))') + + Gqt, ft = gradval(L) + s = 0 + for _ in 1:maxiter + Gp = G - R .* sum.(eachcol(R .* G))' + s = sqrt(sum(diag(Gp' * Gp))) + # if s < ϵ + # break + # end + α *= 2.0 + Rt = Matrix{eltype(F)}(I, p, p) + for _ in 1:lsiter + X = R - α * Gp + v = 1.0 ./ sqrt.(sum.(eachcol(X.^2))) + Rt = X .* v' + # L = F * inv(Rt)' + L = (Rt \ F')' + Gqt, ft = gradval(L) + if (ft < (f - 0.5 * s^2 * α)) + break + end + α /= 2.0 + end + R = Rt + f = ft + # G = -(L' * Gqt * inv(R))' + G = R' \ (-Gqt' * L) + end + + (s < ϵ) || throw(ConvergenceException(maxiter, s, ϵ)) + + (L, R) +end + """ fit(::Type{FactorRotation}, F::AbstractMatrix; ...) -> FactorRotation diff --git a/test/facrot.jl b/test/facrot.jl index b619fee..b22d650 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -63,8 +63,8 @@ using StableRNGs @test isapprox(FR.F, F, rtol = √eps(Float64)) @test isapprox(FR.R, R, rtol = √eps(Float64)) - isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[1], F, rtol = 1e-6) - isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[2], R, rtol = 1e-6) + isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[1], F, rtol = √eps(Float64)) + isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[2], R, rtol = √eps(Float64)) ## Equamax # Comparison with Matlab's `rotatefactors` called as @@ -98,4 +98,26 @@ using StableRNGs FR = fit(FactorRotation, X, alg = Parsimax(size(X)...)) @test isapprox(FR.F, F, rtol = √eps(Float64)) + ## Quartimin rotation + # Comparison with R's `GPArotation::quartimin` function called as + # ``` + # using RCall + # R"qm <- GPArotation::quartimin($X, eps = 1e-6, maxit = 10000)" + # R"F <- qm$loadings" + # R"R <- qm$Th" + # @rget F + # @rget R + # ``` + + F = [ 0.94295134 1.29282920 + -2.32274063 0.24011390 + 1.01896033 1.41630052 + 0.86570192 0.28326547 + 0.39161552 -2.38397278] + R = [ 0.87548617 0.41144611 + -0.48324317 0.91143409] + + isapprox(gpaoblique(X, quartimin, 1000, 20, 1e-6)[1], F, rtol = √eps(Float64)) + isapprox(gpaoblique(X, quartimin, 1000, 10, 1e-6)[2], R, rtol = √eps(Float64)) + end \ No newline at end of file From 434d484435f2327045f2a67981908c0fe0453b93 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Tue, 5 Jul 2022 14:59:24 +0200 Subject: [PATCH 26/32] Improved interface using dispatch on types --- src/MultivariateStats.jl | 24 +- src/facrot.jl | 623 ++++++++++++++++++++++----------------- test/facrot.jl | 10 +- 3 files changed, 374 insertions(+), 283 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 6f278df..96fe700 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -8,7 +8,7 @@ module MultivariateStats ConvergenceException, pairwise, pairwise!, CoefTable import Statistics: mean, var, cov, covm, cor - import Base: length, size, show, eltype + import Base: length, size, show import StatsAPI: fit, predict, coef, weights, dof, r2 import LinearAlgebra: eigvals, eigvecs @@ -108,17 +108,17 @@ module MultivariateStats facm, # CM algorithm for factor analysis ## facrot - FactorRotationAlgorithm, # Type: Factor rotation algorithm - Orthomax, # Type: Orthomax factor rotation algorithm - # Convenience types for factor rotation algorithms - Varimax, # Type - Quartimax, # Type - Equamax, # Type - Parsimax, # Type - - FactorRotation, # Type: Return type for factor rotations - - rotatefactors # Alternative interface to factor rotations + FactorRotationMethod, # Type: a factor rotation method + Orthogonal, # Type: Orthogonal factor rotation method + Oblique, # Type: Oblique factor rotation method + + FactorRotationCriterion, # Type: a factor rotation criterion + CrawfordFerguson, # Type: Crawford-Fergusion rotation criterion + Varimax, # Type: Varimax rotation criterion + Quartimax, # Type: Quartimax rotation criterion + MinimumEntropy, # Type: MinimumEntropy rotation criterion + + rotate # Rotate factors ## source files include("types.jl") diff --git a/src/facrot.jl b/src/facrot.jl index c5a0994..0684dca 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -1,314 +1,425 @@ """ - FactorRotationAlgorithm + FactorRotationMethod -An abstract type for factor rotation algorithms. +An abstract type for factor rotation methods. """ -abstract type FactorRotationAlgorithm end +abstract type FactorRotationMethod end """ - Orthomax <: FactorRotationAlgorithm + Orthogonal -A type representing the orthomax factor rotation algorithm. - -The positive parameter `γ::Real` determines which type of rotation is -performed. For a `n x p` matrix, the default `γ = 1.0` leads to varimax -rotation, `γ = 0.0` is quartimax rotation, `γ = p / 2` is equamax -rotation, and `γ = n (p - 1) / (n + p - 2)` is parsimax rotation. +A type representing orthogonal factor rotations. +""" +struct Orthogonal <: FactorRotationMethod end -The parameter `maxiter::Integer` controls the maximum number of iterations -to perform (default `1000`), `miniter::Integer` controls the minimum number -of iterations taken before convergence is checked, and `ϵ::Real` is a small -positive constant determining convergence (default `1e-12`). """ -struct Orthomax <: FactorRotationAlgorithm - γ::Real - miniter::Integer - maxiter::Integer - ϵ::Real - - Orthomax(;γ::Union{Real, Integer} = 1.0, - miniter::Integer = 20, - maxiter::Integer = 1000, - ϵ::Real = 1e-12) = begin - γ ≥ zero(eltype(γ)) || - throw(DomainError("Orthomax: γ needs to be non-negative")) - miniter > zero(eltype(miniter)) || - throw(DomainError("Orthomax: miniter needs to be positive")) - maxiter > zero(eltype(maxiter)) || - throw(DomainError("Orthomax: maxiter needs to be positive")) - ϵ > zero(eltype(ϵ)) || - throw(DomainError("Orthomax: ϵ needs to be positive")) - - new(float(γ), miniter, maxiter, ϵ) - end -end + Oblique +A type representing oblique factor rotations. """ - Varimax() -> Orthomax +struct Oblique <: FactorRotationMethod end -Creates an orthomax factor rotation algorithm object for a matrix of -size `n x p` with `γ = 1.0`. -Remaining keyword parameters as for [`Orthomax`](@ref). """ -function Varimax(;miniter::Integer = 20, - maxiter::Integer = 1000, - ϵ::Real = 1e-12) - Orthomax(; γ = 1.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) -end + FactorRotationCriterion{T <: FactorRotationMethod} +An abstract type for factor rotation criteria for a specific factor +rotation type. """ - Quartimax() -> Orthomax +abstract type FactorRotationCriterion{T <: FactorRotationMethod} end -Creates an orthomax factor rotation algorithm object for a matrix of -size `n x p` with `γ = 0.0`. -Remaining keyword parameters as for [`Orthomax`](@ref). """ -function Quartimax(;miniter::Integer = 20, - maxiter::Integer = 1000, - ϵ::Real = 1e-12) - Orthomax(γ = 0.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) + CrawfordFerguson{T} <: FactorRotationCriterion{T} + +A type representing the Crawford-Ferguson family of factor rotation criteria. +These criteria are valid for both orthogonal and oblique rotation. The method +of rotation can be set through the type parameter `T`. + +This criterion minimizes + +`Q(Λ) = (1 - κ) * tr((Λ.^2)' * Λ.^2 * N) / 4 + κ * tr((Λ.^2)' * M * Λ.^2) / 4` + +where `Λ` is a `d x p` rotated loading matrix, `N` is a `p x p` +matrix with zeros on the diagonal and ones everywhere else, `M` is an +analogous `d x d` matrix, and `κ` is a non-negative shape parameter. + +**Parameters** +- `κ` is a non-negative shape parameter. In the orthogonal setting, + some classical special cases are + - `κ = 0` is quartimax rotation + - `κ = 1 / d` is varimax rotation + - `κ = p / (2 * d)` is equimax rotation + - `κ = (p - 1) / (d + p - 2)` is parsimax rotation + - `κ = 1` is factor parsimony rotation + +**References** +- Crawford, C.B. and Ferguson, G.A. (1970). A general rotation criterion and + its use in orthogonal rotation. Psychometrika, 35, 321-332. + doi 10.1007/BF02310792 +- Browne, M.W. (2001). An overview of analytic rotation in exploratory factor + analysis. Multivariate Behavioral Research, 36, 111-150. + doi 10.1207/S15327906MBR3601_05 +""" +struct CrawfordFerguson{T} <: FactorRotationCriterion{T} + κ::Real + + CrawfordFerguson{T}(; κ::Union{Real, Integer} = 0.0) where {T <: FactorRotationMethod} = begin + κ ≥ zero(eltype(κ)) || + throw(DomainError("CrawfordFerguson: κ needs to be non-negative")) + + new(float(κ)) + end end -""" - Equamax(p) -> Orthomax +function ∇Qf(L::AbstractMatrix, C::CrawfordFerguson) + d, p = size(L) + N = ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p) + M = ones(eltype(L), d, d) - Matrix{eltype(L)}(I, d, d) + L2 = L.^2 -Creates an orthomax factor rotation algorithm object for a matrix of -size `n x p` with `γ = p / 2`. -Remaining keyword parameters as for [`Orthomax`](@ref). -""" -function Equamax(p::Integer; miniter::Integer = 20, - maxiter::Integer = 1000, - ϵ::Real = 1e-12) - Orthomax(γ = float(p) / 2.0, miniter = miniter, maxiter = maxiter, ϵ = ϵ) + return ( + (1 - C.κ) * L .* (L2 * N) + C.κ * L .* (M * L2), + (1 - C.κ) * sum(L2 .* (L2 * N)) / 4.0 + C.κ * sum(L2 .* (M * L2)) / 4.0 + ) end """ - Parsimax(n, p) -> Orthomax + Varimax <: FactorRotationCriterion{Orthogonal} -Creates an orthomax factor rotation algorithm object for a matrix of -size `n x p` with `γ = n (p - 1) / (n + p - 2)`. -Remaining keyword parameters as for [`Orthomax`](@ref). -""" -function Parsimax(n::Integer, p::Integer; miniter::Integer = 20, - maxiter::Integer = 1000, - ϵ::Real = 1e-12) - Orthomax(γ = float(n) * (float(p) - 1) / (float(n) + float(p) - 2), - miniter = miniter, maxiter = maxiter, ϵ = ϵ) -end +Convenience type for the varimax factor rotation criterion. -function show(io::IO, mime::MIME{Symbol("text/plain")}, alg::Orthomax) - summary(io, alg); println(io) - println(io, "γ = $(alg.γ)") - println(io, "miniter = $(alg.miniter)") - println(io, "maxiter = $(alg.maxiter)") - println(io, "ϵ = $(alg.ϵ)") -end +This criterion minimizes + +`Q(Λ) = -norm(Λ.^2 .- mean.(eachcol(Λ.^2))')^2 / 4.0` +where `Λ` is a `d x p` matrix of rotated loadings. + +**References** +- Kaiser, H.F. (1958). The varimax criterion for analytic rotation in factor + analysis. Psychometrika, 23, 187-200. +- Harman, H.H. (1976). Modern factor analysis (3rd. ed.). + Chicago: The University of Chicago Press. Page 290. """ - FactorRotation{T <: Real} +struct Varimax <: FactorRotationCriterion{Orthogonal} end -The return type for factor rotations. +function ∇Qf(L::AbstractMatrix, ::Varimax) + Q = L.^2 .- mean.(eachcol(L.^2))' + return (-L .* Q, -norm(Q)^2 / 4.0) +end -`F` contains the rotated factors and `R` is the rotation matrix that -was applied to the original factors. """ -struct FactorRotation{T <: Real, Ta <: FactorRotationAlgorithm} - F::Matrix{T} - R::Matrix{T} + Quartimax <: FactorRotationCriterion{Orthogonal} - alg::Ta +Convenience type for the quartimax factor rotation criterion. - FactorRotation{T, Ta}(F, R, alg) where {T <: Real, Ta <: FactorRotationAlgorithm} = new{T, Ta}(F, R, alg) -end +This criterion minimizes -eltype(::Type{FactorRotation{T, Ta}}) where {T,Ta} = T +`Q(Λ) = -norm(Λ.^2)^2 / 4.0` -FactorRotation(F::Matrix{T}, R::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation{T, Ta}(F, R, alg) -FactorRotation(F::Matrix{T}, alg::Ta) where {T <: Real, Ta <: FactorRotationAlgorithm} = FactorRotation(F, Matrix{eltype(F)}(I, size(F, 2), size(F, 2)), alg) +where `Λ` is a `d x p` matrix of rotated loadings. -function FactorRotation(F::Matrix, R::Matrix, alg::Ta) where {Ta <: FactorRotationAlgorithm} - return FactorRotation(promote(F, R)..., alg) -end +**References** +- Carroll, J.B. (1953). An analytic solution for approximating simple + structure in factor analysis. Psychometrika, 18, 23-38. +- Ferguson, G.A. (1954). The concept of parsimony in factor analysis. + Psychometrika, 19, 281-290. +- Neuhaus, J.O. & Wrigley, C. (1954). The quartimax method: An analytical + approach to orthogonal simple structure. + British Journal of Statistical Psychology, 7, 81-91. +""" +struct Quartimax <: FactorRotationCriterion{Orthogonal} end -function show(io::IO, mime::MIME{Symbol("text/plain")}, FR::FactorRotation{<:Any,<:Any}) - summary(io, FR); println(io) +function ∇Qf(L::AbstractMatrix, ::Quartimax) + Q = L.^2 + return (-L .* Q, -norm(Q)^2 / 4.0) end """ - loadings(FR::FactorRotation) + MinimumEntropy <: FactorRotationCriterion{Orthogonal} -Returns the loading matrix. -""" -loadings(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.F +A type representing a simple entropy criterion for factor rotation. + +This criterion minimizes +`Q(Λ) = -tr((Λ.^2)' * log.(Λ.^2)) / 2.0` + +where `Λ` is a `d x p` matrix of rotated loadings. + +**Note** +No oblique version of this criterion exits. + +**References** +- Jennrich, R. I. (2004). Rotation to simple loadings using component + loss functions: The orthogonal case. Psychometrika, 69, 257-273. + doi 10.1007/BF02295943 """ - rotation(FR::FactorRotation) +struct MinimumEntropy <: FactorRotationCriterion{Orthogonal} end + +function ∇Qf(L::AbstractMatrix, ::MinimumEntropy) + L2 = L.^2 + + return ( + -L .* log.(L2) - L, + -sum(L2 .* log.(L2)) / 2.0 + ) +end -Returns the rotation matrix. """ -rotation(FR::FactorRotation{T, Ta}) where {T, Ta} = FR.R - -## Comparison to other implementations -# The implementation of varimax in R row-normlises the input matrix before -# application of the algorithm and rescales the rows afterwards. -## Reference -# 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 -function orthomax(F::AbstractMatrix, γ, miniter, maxiter, ϵ) - n, p = size(F) - if n < 2 - return (F, Matrix{eltype(F)}(I, p, p)) - end + Oblimin{T} <: FactorRotationCriterion{T} + +Type for the oblimin family of factor rotation criteria. +These criteria are valid for both orthogonal and oblique rotation. The method +of rotation can be set through the type parameter `T`. + +This criterion minimizes + +`Q(Λ) = -tr((Λ.^2)' * (I - γ / d * C) * Λ.^2 * N) / 4.0` + +where `Λ` is a `d x p` matrix of rotated loadings, `I` is the `d`-dimensional +identity matrix, `C` is a `d x d` matrix with only ones, and `N` is a `p x p` +matrix with zeros on the diagonal and ones everywhere else. + +**Parameters** +- `γ` is a shape parameter. Negative values are allowed and might be useful + for oblique rotations. + + In the setting of oblique factor rotation, some special cases are + - `γ = 0` is the quartimin criterion + - `γ = 1/2` is the biquartimin criterion + - `γ = 1` is the covarimin criterion + + In the setting of orthogonal factor rotation, the oblimin family of factor + rotations is equivalent to the orthomax family of factor rotation criteria + and some special cases are + - `γ = 0` is the quartimax criterion + - `γ = 0.5` is the biquartimax criterion + - `γ = 1` is the varimax criterion + - `γ = d / 2.0` is the equamax criterion + +**References** +- Carroll, J.B. (1960). IBM 704 program for generalized analytic rotation + solution in factor analysis. Harvard University, unpublished. +- Harman, H.H. (1976). Modern factor analysis (3rd. ed.). + Chicago: The University of Chicago Press. Page 322. +- Jennrich, R.I. (1979). Admissible values of κ in direct oblimin rotation. + Psychometrika, 44, 173-177. +""" +struct Oblimin{T} <: FactorRotationCriterion{T} + γ::Real - # Warm up step - # Do one step. If that first step did not lead away from the identity - # matrix enough use a random orthogonal matrix as a starting point. - M = svd(F' * (F .^ 3 - γ / n * F * Diagonal(vec(sum(F .^ 2, dims=1))))) - R = M.U * M.Vt - if norm(R - Matrix{eltype(R)}(I, p, p)) < ϵ - R = qr(randn(p, p)).Q + Oblimin{T}(; γ::Union{Real, Integer} = 0.0) where {T <: FactorRotationMethod} = begin + new(float(γ)) end +end - # Main iterations - d = 0 - lastchange = NaN - converged = false - for i in 1:maxiter - dold = d - B = F * R - M = svd(F' * (B .^ 3 - γ / n * B * Diagonal(vec(sum(B .^ 2, dims=1))))) - R = M.U * M.Vt - d = sum(M.S) - lastchange = abs(d - dold) / d - if lastchange < ϵ && i >= miniter - converged = true - break - end +function ∇Qf(L::AbstractMatrix, C::Oblimin) + d, p = size(L) + Q = L.^2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) + if C.γ != zero(eltype(C.γ)) + Q = (Matrix{eltype(L)}(I, d, d) - C.γ / d * ones(eltype(L), d, d)) * Q end + return (L * Q, sum(L.^2 .* Q) / 4.0) +end + +""" + Quartimin <: FactorRotationCriterion{Oblique} - converged || throw(ConvergenceException(maxiter, lastchange, ϵ)) +Convenience type for the quartimin factor rotation criterion. - (F * R, R) -end +This criterion minimizes -# Compute value of the optimization criterion and gradient of the -# underlying quality measure -function varimax(L::AbstractMatrix) - Q = L.^2 .- mean.(eachcol(L.^2))' - (-L .* Q, -abs(sum(diag(Q' * Q))) / 4) -end +`Q(Λ) = -tr((Λ.^2)' * Λ.^2 * N) / 4.0` -function quartimax(L::AbstractMatrix) - Q = L.^2 - (-L .* Q, -sum(diag(Q' * Q)) / 4) +where `Λ` is a `d x p` matrix of rotated loadings, and `N` is a `p x p` +matrix with zeros on the diagonal and ones everywhere else. + +**References** +- Carroll, J.B. (1960). IBM 704 program for generalized analytic rotation + solution in factor analysis. Harvard University, unpublished. +""" +struct Quartimin <: FactorRotationCriterion{Oblique} end + +function ∇Qf(L::AbstractMatrix, ::Quartimin) + _, p = size(L) + Q = L.^2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) + return (L .* Q, sum(L.^2 .* Q) / 4.0) end -## Orthogonal rotations computed by gradient projection algorithms -# This algorithm computes orthogonal rotations using a gradient -# project algorithm. -## Reference -# Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms -# and Software for Arbitrary Rotation Criteria in Factor Analysis. -# Educational and Psychological Measurement, 65, 676–696 -# doi 10.1177/0013164404272507 -function gpaortho(F::AbstractMatrix, gradval, maxiter, lsiter, ϵ) - n, p = size(F) - if n < 2 +""" + gparotate(F, C::FactorRotationCriterion{Orthogonal}; ...) + +Compute an orthogonal rotation of a loading matrix according to the +provided rotation criterion. A gradient projection algorithm is used +to perform the computation + +**Parameters** +- `F` is the matrix of loadings to be rotated +- `C` is a factor rotation criterion +- If `normalizerows` is true, then the rows of `F` are normalized to + length 1 before rotation and the rows of the rotated loadings are + scaled back. +- If `randominit` is true, then random initial values are used. + Otherwise, the algorithm is initialized with the identity matrix. +- `maxiter` determines the maximum number of iterations +- `lsiter` determines the maximum number of iterations spent on + line search +- `ϵ` is the convergence tolerance + +**References** +- Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms + and Software for Arbitrary Rotation Criteria in Factor Analysis. + Educational and Psychological Measurement, 65, 676-696. + doi 10.1177/0013164404272507 +- Jennrich, R. I. (2001). A simple general procedure for orthogonal rotation. + Psychometrika, 66, 289-306. doi 10.1007/BF02294840 +""" +function gparotate(F::AbstractMatrix, + C::FactorRotationCriterion{Orthogonal}; + normalizerows = false, + randominit = false, + maxiter::Integer = 1000, + lsiter::Integer = 10, + ϵ::Float64 = 1.0e-6) + d, p = size(F) + if d < 2 return (F, Matrix{eltype(F)}(I, p, p)) end + w = zeros(eltype(F), d) + if normalizerows + w = norm.(eachrow(F)) + F ./= w + end + # Setup - R = Matrix{eltype(F)}(I, p, p) + if randominit + T = qr(randn(eltype{F}, p, p)).Q + else + T = Matrix{eltype(F)}(I, p, p) + end α = 1.0 - L = F * R + L = F * T - Gq, f = gradval(L) - G = F' * Gq + ∇Q, f = ∇Qf(L, C) + ∇f = F' * ∇Q - Gqt, ft = gradval(L) s = 0 for _ in 1:maxiter - M = R' * G + M = T' * ∇f S = (M + M') / 2 - Gp = G - R * S - s = sqrt(sum(diag(Gp' * Gp))) + ∇fp = ∇f - T * S + + # Check for convergence + s = norm(∇fp) if s < ϵ break end + α *= 2.0 - Rt = Matrix{eltype(F)}(I, p, p) + # Create temporaries here so they are not local to the loop + Tt = zeros(eltype(T), p, p) + ∇Qt, ft = ∇Q, f for _ in 1:lsiter - X = R - α * Gp + # Line search to project the gradient step back onto + # the Stiefel manifold + X = T - α * ∇fp UDV = svd(X) - Rt = UDV.U * UDV.Vt - L = F * Rt - Gqt, ft = gradval(L) + Tt = UDV.U * UDV.Vt + L = F * Tt + ∇Qt, ft = ∇Qf(L, C) if (ft < (f - 0.5 * s^2 * α)) break end α /= 2.0 end - R = Rt + T = Tt f = ft - G = F' * Gqt + ∇f = F' * ∇Qt end (s < ϵ) || throw(ConvergenceException(maxiter, s, ϵ)) - (L, R) -end + if normalizerows + L .*= w + end -function quartimin(L::AbstractMatrix) - _, p = size(L) - Q = L.^2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) - (L .* Q, sum(L.^2 .* Q) / 4) + (L, T) end -## Oblique rotations computed by gradient projection algorithms -# This algorithm computes oblique rotations using a gradient -# project algorithm. -## Reference -# Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms -# and Software for Arbitrary Rotation Criteria in Factor Analysis. -# Educational and Psychological Measurement, 65, 676–696 -# doi 10.1177/0013164404272507 -function gpaoblique(F::AbstractMatrix, gradval, maxiter, lsiter, ϵ) +""" + gparotate(F, C::FactorRotationCriterion{Oblique}; ...) + +Compute an oblique rotation of a loading matrix according to the +provided rotation criterion. A gradient projection algorithm is used +to perform the computation + +**Parameters** +- `F` is the matrix of loadings to be rotated +- `C` is a factor rotation criterion +- If `normalizerows` is true, then the rows of `F` are normalized to + length 1 before rotation and the rows of the rotated loadings are + scaled back. +- If `randominit` is true, then random initial values are used. + Otherwise, the algorithm is initialized with the identity matrix. +- `maxiter` determines the maximum number of iterations +- `lsiter` determines the maximum number of iterations spent on + line search +- `ϵ` is the convergence tolerance + +**References** +- Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms + and Software for Arbitrary Rotation Criteria in Factor Analysis. + Educational and Psychological Measurement, 65, 676-696. + doi 10.1177/0013164404272507 +- Jennrich, R. I. (2002). A simple general method for oblique rotation. + Psychometrika, 67, 7-19. doi 10.1007/BF02294706 +""" +function gparotate(F::AbstractMatrix, + C::FactorRotationCriterion{Oblique}; + normalizerows = false, + randominit = false, + maxiter::Integer = 1000, + lsiter::Integer = 10, + ϵ::Float64 = 1.0e-6) n, p = size(F) if n < 2 return (F, Matrix{eltype(F)}(I, p, p)) end + w = zeros(eltype(F), d) + if normalizerows + w = norm.(eachrow(F)) + F ./= w + end + # Setup - R = Matrix{eltype(F)}(I, p, p) + if randominit + T = randn(eltype{F}, p, p) + T ./= norm.(eachcol(T))' + else + T = Matrix{eltype(F)}(I, p, p) + end α = 1.0 - # L = F * inv(R)' - L = (R \ F')' - # isapprox((R \ F')', F * inv(R)') + L = (T \ F')' - Gq, f = gradval(L) - # G = -(L' * Gq * inv(R))' - G = R' \ (-Gq' * L) - # isapprox(R' \ (-Gq' * L), -(L' * Gq * inv(R))') + ∇Q, f = ∇Qf(L, C) + ∇f = T' \ (-∇Q' * L) - Gqt, ft = gradval(L) + ∇Qt, ft = ∇Qf(L, C) s = 0 for _ in 1:maxiter - Gp = G - R .* sum.(eachcol(R .* G))' - s = sqrt(sum(diag(Gp' * Gp))) - # if s < ϵ - # break - # end + ∇fp = ∇f - T .* sum.(eachcol(T .* ∇f))' + s = norm(∇fp) + if s < ϵ + break + end α *= 2.0 Rt = Matrix{eltype(F)}(I, p, p) for _ in 1:lsiter - X = R - α * Gp - v = 1.0 ./ sqrt.(sum.(eachcol(X.^2))) + X = R - α * ∇fp + v = 1.0 ./ norm.(eachcol(X)) Rt = X .* v' - # L = F * inv(Rt)' L = (Rt \ F')' - Gqt, ft = gradval(L) + ∇Qt, ft = ∇Qf(L, C) if (ft < (f - 0.5 * s^2 * α)) break end @@ -316,70 +427,44 @@ function gpaoblique(F::AbstractMatrix, gradval, maxiter, lsiter, ϵ) end R = Rt f = ft - # G = -(L' * Gqt * inv(R))' - G = R' \ (-Gqt' * L) + ∇f = R' \ (-∇Qt' * L) end (s < ϵ) || throw(ConvergenceException(maxiter, s, ϵ)) - (L, R) -end - -""" - fit(::Type{FactorRotation}, F::AbstractMatrix; ...) -> FactorRotation - -Fit a factor rotation to the matrix `F` and apply it. The algorithm used to -perform the factor rotation is by default [`Orthomax`](@ref) and can be changed -with the keyword argument `alg` which is of type `<:FactorRotationAlgorithm`. -""" -function fit(::Type{FactorRotation}, F::AbstractMatrix; - alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - if isa(alg, Orthomax) - F, R = orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) - return FactorRotation(F, R, alg) + if normalizerows + L .*= w end -end - -""" - fit(::Type{FactorRotation}, F::FactorAnalysis; ...) -> FactorRotation - -Fit a factor rotation to the loading matrix of the [`FactorAnalysis`](@ref) -object and apply it. -See [`fit(::Type{FactorRotation}, F::AbstractMatrix)`](@ref) for keyword arguments. -""" -function fit(::Type{FactorRotation}, F::FactorAnalysis; - alg::T = Orthomax()) where {T <: FactorRotationAlgorithm} - return fit(FactorRotation, F.W; alg = alg) -end - -## Alternative interface - -""" - rotatefactors(F::AbstractMatrix, [alg::FactorRotationAlgorithm]) -> FactorRotation - -Rotate the factors in the matrix `F`. The algorithm to be used can be passed in -via the second argument `alg`. By default [`Orthomax`](@ref) is used. -""" -function rotatefactors(F::AbstractMatrix, alg::FactorRotationAlgorithm = Orthomax()) - F, R = _rotatefactors(F, alg) - return FactorRotation(F, R, alg) -end - -# Use multiple dispatch to decide on the algorithm implementation - -function _rotatefactors(F::AbstractMatrix, alg::Orthomax) - return orthomax(F, alg.γ, alg.miniter, alg.maxiter, alg.ϵ) + (L, R) end """ - rotatefactors(F::FactorAnalysis, [alg::FactorRotationAlgorithm]) -> FactorAnalysis - -Rotate the factors in the loading matrix of `F` which is of type [`FactorAnalysis`](@ref). -The algorithm to be used can be passed in via the second argument `alg`. -By default [`Orthomax`](@ref) is used. + rotate(F, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) + +Rotate the factors in matrix `F` using the criterion `C`. + +**Parameters** +- `F` is the matrix of loadings to be rotated +- `C` is a factor rotation criterion +- If `normalizerows` is true, then the rows of `F` are normalized to + length 1 before rotation and the rows of the rotated loadings are + scaled back. +- If `randominit` is true, then random initial values are used. + Otherwise, the algorithm is initialized with the identity matrix. +- `maxiter` determines the maximum number of iterations +- `lsiter` determines the maximum number of iterations spent on + line search +- `ϵ` is the convergence tolerance """ -function rotatefactors(F::FactorAnalysis, alg::FactorRotationAlgorithm = Orthomax()) - FR = rotatefactors(F.W, alg) - return FactorAnalysis(F.mean, FR.F, F.Ψ) +function rotate(F::AbstractMatrix, + C::FactorRotationCriterion{T}; + normalizerows::Bool = false, + randominit::Bool = false, + maxiter::Integer = 1000, + lsiter::Integer = 10, + ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} + L, R = gparotate(F, C; normalizerows, randominit, maxiter, lsiter, ϵ) + + return (L, R) end diff --git a/test/facrot.jl b/test/facrot.jl index b22d650..05b9ffc 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -40,6 +40,10 @@ using StableRNGs isapprox(gpaortho(X, varimax, 1000, 10, 1e-6)[1], F, rtol = 1e-6) isapprox(gpaortho(X, varimax, 1000, 10, 1e-6)[2], R, rtol = 1e-6) + FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = 1.0 / 5.0)) + isapprox(FR[1], F, rtol = 1e-6) + isapprox(FR[2], F, rtol = 1e-6) + ## Quartimax rotation # Comparison with R's `GPArotation::quartimax` function called as # ``` @@ -85,9 +89,11 @@ using StableRNGs ## Parsimax # Comparison with Matlab's `rotatefactors` called as + # ``` # using MATLAB # F = mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')" - + # ``` + F = [ 0.90503611 1.24332782; -2.32641022 0.36067319; 0.97744298 1.36280122; @@ -117,7 +123,7 @@ using StableRNGs R = [ 0.87548617 0.41144611 -0.48324317 0.91143409] - isapprox(gpaoblique(X, quartimin, 1000, 20, 1e-6)[1], F, rtol = √eps(Float64)) + isapprox(gpaoblique(X, quartimin, 1000, 50, 1e-6)[1], F, rtol = √eps(Float64)) isapprox(gpaoblique(X, quartimin, 1000, 10, 1e-6)[2], R, rtol = √eps(Float64)) end \ No newline at end of file From 0211a3a3fba655f7d45f70496176ad0d2be8fafe Mon Sep 17 00:00:00 2001 From: Felix Held Date: Tue, 5 Jul 2022 19:52:36 +0200 Subject: [PATCH 27/32] Add interface for statistical models --- src/MultivariateStats.jl | 10 ++++- src/facrot.jl | 95 ++++++++++++++++++++++++++++++++-------- test/facrot.jl | 61 ++++++++++++++++---------- 3 files changed, 122 insertions(+), 44 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 96fe700..73230ed 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -113,12 +113,18 @@ module MultivariateStats Oblique, # Type: Oblique factor rotation method FactorRotationCriterion, # Type: a factor rotation criterion - CrawfordFerguson, # Type: Crawford-Fergusion rotation criterion + CrawfordFerguson, # Type: Crawford-Fergusion rotation criteria Varimax, # Type: Varimax rotation criterion Quartimax, # Type: Quartimax rotation criterion MinimumEntropy, # Type: MinimumEntropy rotation criterion + Oblimin, # Type: Oblimin rotation criteria + Quartimin, # Type: Quartimin rotation criterion - rotate # Rotate factors + FactorRotation, # Type: Result of a general factor rotation + rotation, # Extract rotation matrix + + rotate, # Rotate factors + rotate! # Rotate in-place ## source files include("types.jl") diff --git a/src/facrot.jl b/src/facrot.jl index 0684dca..bdc0835 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -47,7 +47,7 @@ analogous `d x d` matrix, and `κ` is a non-negative shape parameter. some classical special cases are - `κ = 0` is quartimax rotation - `κ = 1 / d` is varimax rotation - - `κ = p / (2 * d)` is equimax rotation + - `κ = p / (2 * d)` is equamax rotation - `κ = (p - 1) / (d + p - 2)` is parsimax rotation - `κ = 1` is factor parsimony rotation @@ -217,7 +217,7 @@ function ∇Qf(L::AbstractMatrix, C::Oblimin) if C.γ != zero(eltype(C.γ)) Q = (Matrix{eltype(L)}(I, d, d) - C.γ / d * ones(eltype(L), d, d)) * Q end - return (L * Q, sum(L.^2 .* Q) / 4.0) + return (L .* Q, sum(L.^2 .* Q) / 4.0) end """ @@ -342,7 +342,7 @@ function gparotate(F::AbstractMatrix, L .*= w end - (L, T) + (Matrix(L), Matrix(T)) end """ @@ -380,8 +380,8 @@ function gparotate(F::AbstractMatrix, maxiter::Integer = 1000, lsiter::Integer = 10, ϵ::Float64 = 1.0e-6) - n, p = size(F) - if n < 2 + d, p = size(F) + if d < 2 return (F, Matrix{eltype(F)}(I, p, p)) end @@ -404,30 +404,36 @@ function gparotate(F::AbstractMatrix, ∇Q, f = ∇Qf(L, C) ∇f = T' \ (-∇Q' * L) - ∇Qt, ft = ∇Qf(L, C) s = 0 for _ in 1:maxiter ∇fp = ∇f - T .* sum.(eachcol(T .* ∇f))' s = norm(∇fp) + + # Check for convergence if s < ϵ break end + α *= 2.0 - Rt = Matrix{eltype(F)}(I, p, p) + # Create temporaries here so they are not local to the loop + Tt = zeros(eltype(F), p, p) + ∇Qt, ft = ∇Q, f for _ in 1:lsiter - X = R - α * ∇fp + # Line search to project the gradient step back onto + # the oblique manifold + X = T - α * ∇fp v = 1.0 ./ norm.(eachcol(X)) - Rt = X .* v' - L = (Rt \ F')' + Tt = X .* v' + L = (Tt \ F')' ∇Qt, ft = ∇Qf(L, C) if (ft < (f - 0.5 * s^2 * α)) break end α /= 2.0 end - R = Rt + T = Tt f = ft - ∇f = R' \ (-∇Qt' * L) + ∇f = T' \ (-∇Qt' * L) end (s < ϵ) || throw(ConvergenceException(maxiter, s, ϵ)) @@ -436,16 +442,29 @@ function gparotate(F::AbstractMatrix, L .*= w end - (L, R) + (Matrix(L), Matrix(T)) end """ - rotate(F, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) + FactorRotation{T <: Real} -Rotate the factors in matrix `F` using the criterion `C`. +This type +""" +struct FactorRotation{T <: Real} + L::Matrix{T} + R::Matrix{T} +end + +loadings(M::FactorRotation) = M.L +rotation(M::FactorRotation) = M.R + +""" + rotate(L, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) + +Rotate the loadings in matrix `L` using the criterion `C`. **Parameters** -- `F` is the matrix of loadings to be rotated +- `L` is the matrix of loadings to be rotated - `C` is a factor rotation criterion - If `normalizerows` is true, then the rows of `F` are normalized to length 1 before rotation and the rows of the rotated loadings are @@ -457,14 +476,52 @@ Rotate the factors in matrix `F` using the criterion `C`. line search - `ϵ` is the convergence tolerance """ -function rotate(F::AbstractMatrix, +function rotate(L::AbstractMatrix, C::FactorRotationCriterion{T}; normalizerows::Bool = false, randominit::Bool = false, maxiter::Integer = 1000, lsiter::Integer = 10, ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} - L, R = gparotate(F, C; normalizerows, randominit, maxiter, lsiter, ϵ) + return FactorRotation( + gparotate(L, C; normalizerows, randominit, maxiter, lsiter, ϵ)... + ) +end + +""" + rotate(M::FactorAnalysis, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) - return (L, R) +Rotate the loadings of the Factor Analysis model `M` using the criterion `C`. +Modifies the loading matrix in `M`. Parameters as in [`rotate`](@ref). +""" +function rotate!(M::FactorAnalysis, + C::FactorRotationCriterion{T}; + normalizerows::Bool = false, + randominit::Bool = false, + maxiter::Integer = 1000, + lsiter::Integer = 10, + ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} + FR = rotate(loadings(M), C; normalizerows, randominit, maxiter, lsiter, ϵ) + M.W .= loadings(FR) + + return M end + +""" + rotate(M::PCA, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) + +Rotate the components of the PCA model `M` using the criterion `C`. +Modifies the projection matrix in `M`. Parameters as in [`rotate`](@ref). +""" +function rotate!(M::PCA, + C::FactorRotationCriterion{T}; + normalizerows::Bool = false, + randominit::Bool = false, + maxiter::Integer = 1000, + lsiter::Integer = 10, + ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} + FR = rotate(projection(M), C; normalizerows, randominit, maxiter, lsiter, ϵ) + M.proj .= loadings(FR) + + return M +end \ No newline at end of file diff --git a/test/facrot.jl b/test/facrot.jl index 05b9ffc..d362341 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -6,7 +6,9 @@ using StableRNGs rng = StableRNG(37273) - X = randn(rng, 5, 2) + d = 5 + p = 2 + X = randn(rng, d, p) ## Varimax rotation # Ground-truth was extracted by specifying the following @@ -33,16 +35,18 @@ using StableRNGs R = [ 0.89940646 0.43711327; -0.43711327 0.89940646] - FR = fit(FactorRotation, X, alg = Varimax()) - @test isapprox(FR.F, F, rtol = √eps(Float64)) - @test isapprox(FR.R, R, rtol = √eps(Float64)) - - isapprox(gpaortho(X, varimax, 1000, 10, 1e-6)[1], F, rtol = 1e-6) - isapprox(gpaortho(X, varimax, 1000, 10, 1e-6)[2], R, rtol = 1e-6) + FR = rotate(X, Varimax()) + @test isapprox(FR[1], F, rtol = 1e-6) + @test isapprox(FR[2], R, rtol = 1e-6) + # Different equivalent ways of computing varimax FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = 1.0 / 5.0)) - isapprox(FR[1], F, rtol = 1e-6) - isapprox(FR[2], F, rtol = 1e-6) + @test isapprox(FR[1], F, rtol = 1e-6) + @test isapprox(FR[2], F, rtol = 1e-6) + + FR = rotate(X, Oblimin{Orthogonal}(γ = 1.0)) + @test isapprox(FR[1], F, rtol = 1e-6) + @test isapprox(FR[2], F, rtol = 1e-6) ## Quartimax rotation # Comparison with R's `GPArotation::quartimax` function called as @@ -63,12 +67,9 @@ using StableRNGs R = [ 0.89748635 0.44104222; -0.44104222 0.89748635] - FR = fit(FactorRotation, X, alg = Quartimax()) - @test isapprox(FR.F, F, rtol = √eps(Float64)) - @test isapprox(FR.R, R, rtol = √eps(Float64)) - - isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[1], F, rtol = √eps(Float64)) - isapprox(gpaortho(X, quartimax, 1000, 10, 1e-6)[2], R, rtol = √eps(Float64)) + FR = rotate(X, Quartimax()) + @test isapprox(loadings(FR), F, rtol = √eps(Float64)) + @test isapprox(rotation(FR), R, rtol = √eps(Float64)) ## Equamax # Comparison with Matlab's `rotatefactors` called as @@ -81,9 +82,9 @@ using StableRNGs 0.85650467 0.23818242; 0.45865486 -2.40335769] - # Equamax for n x p matrix means γ = p / 2 - FR = fit(FactorRotation, X, alg = Equamax(size(X, 2))) - @test isapprox(FR.F, F, rtol = √eps(Float64)) + # Equamax for d x p matrix means orthogonal Crawford-Ferguson with κ = p / (2 * d) + FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = p / (2 * d))) + @test isapprox(loadings(FR), F, rtol = √eps(Float64)) ## Parsimax @@ -100,9 +101,9 @@ using StableRNGs 0.85650467 0.23818242; 0.45865486 -2.40335769] - # Equamax for n x p matrix means γ = p / 2 - FR = fit(FactorRotation, X, alg = Parsimax(size(X)...)) - @test isapprox(FR.F, F, rtol = √eps(Float64)) + # Parsimax for d x p matrix means orthogonal Crawford-Ferguson with κ = (p - 1) / (d + p - 2) + FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = (p - 1) / (d + p - 2))) + @test isapprox(loadings(FR), F, rtol = √eps(Float64)) ## Quartimin rotation # Comparison with R's `GPArotation::quartimin` function called as @@ -123,7 +124,21 @@ using StableRNGs R = [ 0.87548617 0.41144611 -0.48324317 0.91143409] - isapprox(gpaoblique(X, quartimin, 1000, 50, 1e-6)[1], F, rtol = √eps(Float64)) - isapprox(gpaoblique(X, quartimin, 1000, 10, 1e-6)[2], R, rtol = √eps(Float64)) + FR = rotate(X, Quartimin()) + @test isapprox(loadings(FR), F, rtol = √eps(Float64)) + @test isapprox(rotation(FR), R, rtol = √eps(Float64)) + + # Test application to factor analysis and PCA models + + X = randn(10, 5) + + M = fit(FactorAnalysis, X) + loadings(M) + rotate!(M, Varimax()) + loadings(M) + M = fit(PCA, X) + projection(M) + rotate!(M, Varimax()) + projection(M) end \ No newline at end of file From 8c206d92036cb08a0dc5ae056c0b28ed4ee8d2cb Mon Sep 17 00:00:00 2001 From: Felix Held Date: Mon, 25 Jul 2022 15:25:56 +0200 Subject: [PATCH 28/32] Adress code review --- src/MultivariateStats.jl | 2 +- src/facrot.jl | 188 +++++++++++++++++++++------------------ test/facrot.jl | 12 +-- 3 files changed, 109 insertions(+), 93 deletions(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 73230ed..41191dd 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -10,7 +10,7 @@ module MultivariateStats import Statistics: mean, var, cov, covm, cor import Base: length, size, show import StatsAPI: fit, predict, coef, weights, dof, r2 - import LinearAlgebra: eigvals, eigvecs + import LinearAlgebra: eigvals, eigvecs, rotate! export diff --git a/src/facrot.jl b/src/facrot.jl index bdc0835..f7fcec8 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -25,7 +25,7 @@ struct Oblique <: FactorRotationMethod end An abstract type for factor rotation criteria for a specific factor rotation type. """ -abstract type FactorRotationCriterion{T <: FactorRotationMethod} end +abstract type FactorRotationCriterion{T<:FactorRotationMethod} end """ CrawfordFerguson{T} <: FactorRotationCriterion{T} @@ -36,20 +36,20 @@ of rotation can be set through the type parameter `T`. This criterion minimizes -`Q(Λ) = (1 - κ) * tr((Λ.^2)' * Λ.^2 * N) / 4 + κ * tr((Λ.^2)' * M * Λ.^2) / 4` +`Q(Λ) = (1 - k) * tr((Λ.^2)' * Λ.^2 * N) / 4 + k * tr((Λ.^2)' * M * Λ.^2) / 4` where `Λ` is a `d x p` rotated loading matrix, `N` is a `p x p` matrix with zeros on the diagonal and ones everywhere else, `M` is an -analogous `d x d` matrix, and `κ` is a non-negative shape parameter. +analogous `d x d` matrix, and `k` is a non-negative shape parameter. **Parameters** -- `κ` is a non-negative shape parameter. In the orthogonal setting, +- `kappa` is a non-negative shape parameter. In the orthogonal setting, some classical special cases are - - `κ = 0` is quartimax rotation - - `κ = 1 / d` is varimax rotation - - `κ = p / (2 * d)` is equamax rotation - - `κ = (p - 1) / (d + p - 2)` is parsimax rotation - - `κ = 1` is factor parsimony rotation + - `kappa = 0` is quartimax rotation + - `kappa = 1 / d` is varimax rotation + - `kappa = p / (2 * d)` is equamax rotation + - `kappa = (p - 1) / (d + p - 2)` is parsimax rotation + - `kappa = 1` is factor parsimony rotation **References** - Crawford, C.B. and Ferguson, G.A. (1970). A general rotation criterion and @@ -60,13 +60,13 @@ analogous `d x d` matrix, and `κ` is a non-negative shape parameter. doi 10.1207/S15327906MBR3601_05 """ struct CrawfordFerguson{T} <: FactorRotationCriterion{T} - κ::Real + kappa::Real - CrawfordFerguson{T}(; κ::Union{Real, Integer} = 0.0) where {T <: FactorRotationMethod} = begin - κ ≥ zero(eltype(κ)) || - throw(DomainError("CrawfordFerguson: κ needs to be non-negative")) + CrawfordFerguson{T}(; kappa::Union{Real,Integer}=0.0) where {T<:FactorRotationMethod} = begin + kappa ≥ zero(eltype(kappa)) || + throw(DomainError("CrawfordFerguson: kappa needs to be non-negative")) - new(float(κ)) + new(float(kappa)) end end @@ -74,11 +74,11 @@ function ∇Qf(L::AbstractMatrix, C::CrawfordFerguson) d, p = size(L) N = ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p) M = ones(eltype(L), d, d) - Matrix{eltype(L)}(I, d, d) - L2 = L.^2 + L2 = L .^ 2 return ( - (1 - C.κ) * L .* (L2 * N) + C.κ * L .* (M * L2), - (1 - C.κ) * sum(L2 .* (L2 * N)) / 4.0 + C.κ * sum(L2 .* (M * L2)) / 4.0 + (1 - C.kappa) * L .* (L2 * N) + C.kappa * L .* (M * L2), + (1 - C.kappa) * sum(L2 .* (L2 * N)) / 4.0 + C.kappa * sum(L2 .* (M * L2)) / 4.0 ) end @@ -102,7 +102,7 @@ where `Λ` is a `d x p` matrix of rotated loadings. struct Varimax <: FactorRotationCriterion{Orthogonal} end function ∇Qf(L::AbstractMatrix, ::Varimax) - Q = L.^2 .- mean.(eachcol(L.^2))' + Q = L .^ 2 .- mean.(eachcol(L .^ 2))' return (-L .* Q, -norm(Q)^2 / 4.0) end @@ -129,7 +129,7 @@ where `Λ` is a `d x p` matrix of rotated loadings. struct Quartimax <: FactorRotationCriterion{Orthogonal} end function ∇Qf(L::AbstractMatrix, ::Quartimax) - Q = L.^2 + Q = L .^ 2 return (-L .* Q, -norm(Q)^2 / 4.0) end @@ -155,7 +155,7 @@ No oblique version of this criterion exits. struct MinimumEntropy <: FactorRotationCriterion{Orthogonal} end function ∇Qf(L::AbstractMatrix, ::MinimumEntropy) - L2 = L.^2 + L2 = L .^ 2 return ( -L .* log.(L2) - L, @@ -172,28 +172,29 @@ of rotation can be set through the type parameter `T`. This criterion minimizes -`Q(Λ) = -tr((Λ.^2)' * (I - γ / d * C) * Λ.^2 * N) / 4.0` +`Q(Λ) = -tr((Λ.^2)' * (I - gamma / d * C) * Λ.^2 * N) / 4.0` where `Λ` is a `d x p` matrix of rotated loadings, `I` is the `d`-dimensional -identity matrix, `C` is a `d x d` matrix with only ones, and `N` is a `p x p` -matrix with zeros on the diagonal and ones everywhere else. +identity matrix, `C` is a `d x d` matrix with only ones, `N` is a `p x p` +matrix with zeros on the diagonal and ones everywhere else, and `gamma` is a +shape parameter. **Parameters** -- `γ` is a shape parameter. Negative values are allowed and might be useful +- `gamma` is a shape parameter. Negative values are allowed and might be useful for oblique rotations. In the setting of oblique factor rotation, some special cases are - - `γ = 0` is the quartimin criterion - - `γ = 1/2` is the biquartimin criterion - - `γ = 1` is the covarimin criterion + - `gamma = 0` is the quartimin criterion + - `gamma = 1/2` is the biquartimin criterion + - `gamma = 1` is the covarimin criterion In the setting of orthogonal factor rotation, the oblimin family of factor rotations is equivalent to the orthomax family of factor rotation criteria and some special cases are - - `γ = 0` is the quartimax criterion - - `γ = 0.5` is the biquartimax criterion - - `γ = 1` is the varimax criterion - - `γ = d / 2.0` is the equamax criterion + - `gamma = 0` is the quartimax criterion + - `gamma = 0.5` is the biquartimax criterion + - `gamma = 1` is the varimax criterion + - `gamma = d / 2.0` is the equamax criterion **References** - Carroll, J.B. (1960). IBM 704 program for generalized analytic rotation @@ -206,18 +207,18 @@ matrix with zeros on the diagonal and ones everywhere else. struct Oblimin{T} <: FactorRotationCriterion{T} γ::Real - Oblimin{T}(; γ::Union{Real, Integer} = 0.0) where {T <: FactorRotationMethod} = begin + Oblimin{T}(; γ::Union{Real,Integer}=0.0) where {T<:FactorRotationMethod} = begin new(float(γ)) end end function ∇Qf(L::AbstractMatrix, C::Oblimin) d, p = size(L) - Q = L.^2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) - if C.γ != zero(eltype(C.γ)) - Q = (Matrix{eltype(L)}(I, d, d) - C.γ / d * ones(eltype(L), d, d)) * Q + Q = L .^ 2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) + if C.gamma != zero(eltype(C.gamma)) + Q = (Matrix{eltype(L)}(I, d, d) - C.gamma / d * ones(eltype(L), d, d)) * Q end - return (L .* Q, sum(L.^2 .* Q) / 4.0) + return (L .* Q, sum(L .^ 2 .* Q) / 4.0) end """ @@ -240,8 +241,8 @@ struct Quartimin <: FactorRotationCriterion{Oblique} end function ∇Qf(L::AbstractMatrix, ::Quartimin) _, p = size(L) - Q = L.^2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) - return (L .* Q, sum(L.^2 .* Q) / 4.0) + Q = L .^ 2 * (ones(eltype(L), p, p) - Matrix{eltype(L)}(I, p, p)) + return (L .* Q, sum(L .^ 2 .* Q) / 4.0) end """ @@ -273,12 +274,12 @@ to perform the computation Psychometrika, 66, 289-306. doi 10.1007/BF02294840 """ function gparotate(F::AbstractMatrix, - C::FactorRotationCriterion{Orthogonal}; - normalizerows = false, - randominit = false, - maxiter::Integer = 1000, - lsiter::Integer = 10, - ϵ::Float64 = 1.0e-6) + C::FactorRotationCriterion{Orthogonal}; + normalizerows=false, + randominit=false, + maxiter::Integer=1000, + lsiter::Integer=10, + ϵ::Real=1.0e-6) d, p = size(F) if d < 2 return (F, Matrix{eltype(F)}(I, p, p)) @@ -296,7 +297,7 @@ function gparotate(F::AbstractMatrix, else T = Matrix{eltype(F)}(I, p, p) end - α = 1.0 + α = 1 L = F * T ∇Q, f = ∇Qf(L, C) @@ -307,14 +308,14 @@ function gparotate(F::AbstractMatrix, M = T' * ∇f S = (M + M') / 2 ∇fp = ∇f - T * S - + # Check for convergence s = norm(∇fp) if s < ϵ break end - α *= 2.0 + α *= 2 # Create temporaries here so they are not local to the loop Tt = zeros(eltype(T), p, p) ∇Qt, ft = ∇Q, f @@ -326,10 +327,10 @@ function gparotate(F::AbstractMatrix, Tt = UDV.U * UDV.Vt L = F * Tt ∇Qt, ft = ∇Qf(L, C) - if (ft < (f - 0.5 * s^2 * α)) + if (ft < (f - s^2 * α / 2)) break end - α /= 2.0 + α /= 2 end T = Tt f = ft @@ -374,12 +375,12 @@ to perform the computation Psychometrika, 67, 7-19. doi 10.1007/BF02294706 """ function gparotate(F::AbstractMatrix, - C::FactorRotationCriterion{Oblique}; - normalizerows = false, - randominit = false, - maxiter::Integer = 1000, - lsiter::Integer = 10, - ϵ::Float64 = 1.0e-6) + C::FactorRotationCriterion{Oblique}; + normalizerows=false, + randominit=false, + maxiter::Integer=1000, + lsiter::Integer=10, + ϵ::Real=1.0e-6) d, p = size(F) if d < 2 return (F, Matrix{eltype(F)}(I, p, p)) @@ -398,7 +399,7 @@ function gparotate(F::AbstractMatrix, else T = Matrix{eltype(F)}(I, p, p) end - α = 1.0 + α = 1 L = (T \ F')' ∇Q, f = ∇Qf(L, C) @@ -414,7 +415,7 @@ function gparotate(F::AbstractMatrix, break end - α *= 2.0 + α *= 2 # Create temporaries here so they are not local to the loop Tt = zeros(eltype(F), p, p) ∇Qt, ft = ∇Q, f @@ -422,14 +423,14 @@ function gparotate(F::AbstractMatrix, # Line search to project the gradient step back onto # the oblique manifold X = T - α * ∇fp - v = 1.0 ./ norm.(eachcol(X)) + v = 1 ./ norm.(eachcol(X)) Tt = X .* v' L = (Tt \ F')' ∇Qt, ft = ∇Qf(L, C) - if (ft < (f - 0.5 * s^2 * α)) + if (ft < (f - s^2 * α / 2)) break end - α /= 2.0 + α /= 2 end T = Tt f = ft @@ -448,14 +449,27 @@ end """ FactorRotation{T <: Real} -This type +This type contains the result of a factor rotation. The matrix `L` contains +the rotated loadings and the matrix `R` contains the rotation matrix which +was applied to the original loadings. """ -struct FactorRotation{T <: Real} +struct FactorRotation{T<:Real} L::Matrix{T} R::Matrix{T} end +""" + loadings(M::FactorRotation) + +Return rotated loading matrix. +""" loadings(M::FactorRotation) = M.L + +""" + rotation(M::FactorRotation) + +Return rotation matrix used to obtain rotated loadings. +""" rotation(M::FactorRotation) = M.R """ @@ -466,6 +480,8 @@ Rotate the loadings in matrix `L` using the criterion `C`. **Parameters** - `L` is the matrix of loadings to be rotated - `C` is a factor rotation criterion + +**Keyword parameters** - If `normalizerows` is true, then the rows of `F` are normalized to length 1 before rotation and the rows of the rotated loadings are scaled back. @@ -474,17 +490,17 @@ Rotate the loadings in matrix `L` using the criterion `C`. - `maxiter` determines the maximum number of iterations - `lsiter` determines the maximum number of iterations spent on line search -- `ϵ` is the convergence tolerance +- `tol` is the convergence tolerance """ function rotate(L::AbstractMatrix, - C::FactorRotationCriterion{T}; - normalizerows::Bool = false, - randominit::Bool = false, - maxiter::Integer = 1000, - lsiter::Integer = 10, - ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} + C::FactorRotationCriterion{T}; + normalizerows::Bool=false, + randominit::Bool=false, + maxiter::Integer=1000, + lsiter::Integer=10, + tol::Real=1.0e-6) where {T<:FactorRotationMethod} return FactorRotation( - gparotate(L, C; normalizerows, randominit, maxiter, lsiter, ϵ)... + gparotate(L, C; normalizerows, randominit, maxiter, lsiter, tol)... ) end @@ -495,15 +511,15 @@ Rotate the loadings of the Factor Analysis model `M` using the criterion `C`. Modifies the loading matrix in `M`. Parameters as in [`rotate`](@ref). """ function rotate!(M::FactorAnalysis, - C::FactorRotationCriterion{T}; - normalizerows::Bool = false, - randominit::Bool = false, - maxiter::Integer = 1000, - lsiter::Integer = 10, - ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} - FR = rotate(loadings(M), C; normalizerows, randominit, maxiter, lsiter, ϵ) + C::FactorRotationCriterion{T}; + normalizerows::Bool=false, + randominit::Bool=false, + maxiter::Integer=1000, + lsiter::Integer=10, + tol::Real=1.0e-6) where {T<:FactorRotationMethod} + FR = rotate(loadings(M), C; normalizerows, randominit, maxiter, lsiter, tol) M.W .= loadings(FR) - + return M end @@ -514,14 +530,14 @@ Rotate the components of the PCA model `M` using the criterion `C`. Modifies the projection matrix in `M`. Parameters as in [`rotate`](@ref). """ function rotate!(M::PCA, - C::FactorRotationCriterion{T}; - normalizerows::Bool = false, - randominit::Bool = false, - maxiter::Integer = 1000, - lsiter::Integer = 10, - ϵ::Float64 = 1.0e-6) where {T <: FactorRotationMethod} - FR = rotate(projection(M), C; normalizerows, randominit, maxiter, lsiter, ϵ) + C::FactorRotationCriterion{T}; + normalizerows::Bool=false, + randominit::Bool=false, + maxiter::Integer=1000, + lsiter::Integer=10, + tol::Real=1.0e-6) where {T<:FactorRotationMethod} + FR = rotate(projection(M), C; normalizerows, randominit, maxiter, lsiter, tol) M.proj .= loadings(FR) - + return M end \ No newline at end of file diff --git a/test/facrot.jl b/test/facrot.jl index d362341..7f79d3d 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -36,17 +36,17 @@ using StableRNGs -0.43711327 0.89940646] FR = rotate(X, Varimax()) - @test isapprox(FR[1], F, rtol = 1e-6) - @test isapprox(FR[2], R, rtol = 1e-6) + @test isapprox(loadings(FR), F, rtol = 1e-6) + @test isapprox(rotation(FR), R, rtol = 1e-6) # Different equivalent ways of computing varimax FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = 1.0 / 5.0)) - @test isapprox(FR[1], F, rtol = 1e-6) - @test isapprox(FR[2], F, rtol = 1e-6) + @test isapprox(loadings(FR), F, rtol = 1e-6) + @test isapprox(rotation(FR), R, rtol = 1e-6) FR = rotate(X, Oblimin{Orthogonal}(γ = 1.0)) - @test isapprox(FR[1], F, rtol = 1e-6) - @test isapprox(FR[2], F, rtol = 1e-6) + @test isapprox(loadings(FR), F, rtol = 1e-6) + @test isapprox(rotation(FR), R, rtol = 1e-6) ## Quartimax rotation # Comparison with R's `GPArotation::quartimax` function called as From 10f6949e68802980b36a2bbc691f2189768ea83a Mon Sep 17 00:00:00 2001 From: Felix Held Date: Fri, 19 Jan 2024 11:29:41 +0100 Subject: [PATCH 29/32] latest changes --- docs/make.jl | 27 +++++++------ docs/src/facrot.md | 4 ++ src/facrot.jl | 12 +++--- test/facrot.jl | 98 +++++++++++++++++++++++----------------------- test/runtests.jl | 26 ++++++------ 5 files changed, 86 insertions(+), 81 deletions(-) create mode 100644 docs/src/facrot.md diff --git a/docs/make.jl b/docs/make.jl index ddc4db9..a95c1f6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,20 +5,21 @@ if Base.HOME_PROJECT[] !== nothing end makedocs( - sitename = "MultivariateStats.jl", - modules = [MultivariateStats], - pages = ["Home"=>"index.md", - "whiten.md", - "lreg.md", - "lda.md", - "pca.md", - "ica.md", - "cca.md", - "fa.md", - "mds.md", - "Development"=>"api.md"] + sitename="MultivariateStats.jl", + modules=[MultivariateStats], + pages=["Home" => "index.md", + "whiten.md", + "lreg.md", + "lda.md", + "pca.md", + "ica.md", + "cca.md", + "fa.md", + "mds.md", + "facrot.md", + "Development" => "api.md"] ) deploydocs( - repo = "github.com/JuliaStats/MultivariateStats.jl.git" + repo="github.com/JuliaStats/MultivariateStats.jl.git" ) diff --git a/docs/src/facrot.md b/docs/src/facrot.md new file mode 100644 index 0000000..c3771b4 --- /dev/null +++ b/docs/src/facrot.md @@ -0,0 +1,4 @@ +# Factor Rotation + +Factor rotations are used for post-processing of a factor matrix determined +from a [Factor Analysis](../src/fa.md) or [Principal Component Analysis](../src/pca.md). \ No newline at end of file diff --git a/src/facrot.jl b/src/facrot.jl index f7fcec8..7166db9 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -36,11 +36,11 @@ of rotation can be set through the type parameter `T`. This criterion minimizes -`Q(Λ) = (1 - k) * tr((Λ.^2)' * Λ.^2 * N) / 4 + k * tr((Λ.^2)' * M * Λ.^2) / 4` +`Q(Λ) = (1 - kappa) * tr((Λ.^2)' * Λ.^2 * N) / 4 + kappa * tr((Λ.^2)' * M * Λ.^2) / 4` where `Λ` is a `d x p` rotated loading matrix, `N` is a `p x p` matrix with zeros on the diagonal and ones everywhere else, `M` is an -analogous `d x d` matrix, and `k` is a non-negative shape parameter. +analogous `d x d` matrix, and `kappa` is a non-negative shape parameter. **Parameters** - `kappa` is a non-negative shape parameter. In the orthogonal setting, @@ -205,10 +205,10 @@ shape parameter. Psychometrika, 44, 173-177. """ struct Oblimin{T} <: FactorRotationCriterion{T} - γ::Real + gamma::Real - Oblimin{T}(; γ::Union{Real,Integer}=0.0) where {T<:FactorRotationMethod} = begin - new(float(γ)) + Oblimin{T}(; gamma::Union{Real,Integer}=0.0) where {T<:FactorRotationMethod} = begin + new(float(gamma)) end end @@ -500,7 +500,7 @@ function rotate(L::AbstractMatrix, lsiter::Integer=10, tol::Real=1.0e-6) where {T<:FactorRotationMethod} return FactorRotation( - gparotate(L, C; normalizerows, randominit, maxiter, lsiter, tol)... + gparotate(L, C; normalizerows, randominit, maxiter, lsiter, ϵ=tol)... ) end diff --git a/test/facrot.jl b/test/facrot.jl index 7f79d3d..dc2366b 100644 --- a/test/facrot.jl +++ b/test/facrot.jl @@ -27,26 +27,26 @@ using StableRNGs # @rget R # ```` - F = [ 0.90503610 1.24332783; - -2.32641022 0.36067316; - 0.97744296 1.36280123; - 0.85650467 0.23818243; - 0.45865489 -2.40335769] - R = [ 0.89940646 0.43711327; - -0.43711327 0.89940646] + F = [0.90503610 1.24332783 + -2.32641022 0.36067316 + 0.97744296 1.36280123 + 0.85650467 0.23818243 + 0.45865489 -2.40335769] + R = [0.89940646 0.43711327 + -0.43711327 0.89940646] FR = rotate(X, Varimax()) - @test isapprox(loadings(FR), F, rtol = 1e-6) - @test isapprox(rotation(FR), R, rtol = 1e-6) + @test isapprox(loadings(FR), F, rtol=1e-6) + @test isapprox(rotation(FR), R, rtol=1e-6) # Different equivalent ways of computing varimax - FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = 1.0 / 5.0)) - @test isapprox(loadings(FR), F, rtol = 1e-6) - @test isapprox(rotation(FR), R, rtol = 1e-6) + FR = rotate(X, CrawfordFerguson{Orthogonal}(kappa=1.0 / 5.0)) + @test isapprox(loadings(FR), F, rtol=1e-6) + @test isapprox(rotation(FR), R, rtol=1e-6) - FR = rotate(X, Oblimin{Orthogonal}(γ = 1.0)) - @test isapprox(loadings(FR), F, rtol = 1e-6) - @test isapprox(rotation(FR), R, rtol = 1e-6) + FR = rotate(X, Oblimin{Orthogonal}(gamma=1.0)) + @test isapprox(loadings(FR), F, rtol=1e-6) + @test isapprox(rotation(FR), R, rtol=1e-6) ## Quartimax rotation # Comparison with R's `GPArotation::quartimax` function called as @@ -59,32 +59,32 @@ using StableRNGs # @rget R # ``` - F = [ 0.89959033 1.24727370; - -2.32796522 0.35049625; - 0.97147404 1.36706259; - 0.85545490 0.24192567; - 0.46916046 -2.40132899] - R = [ 0.89748635 0.44104222; - -0.44104222 0.89748635] + F = [0.89959033 1.24727370 + -2.32796522 0.35049625 + 0.97147404 1.36706259 + 0.85545490 0.24192567 + 0.46916046 -2.40132899] + R = [0.89748635 0.44104222 + -0.44104222 0.89748635] FR = rotate(X, Quartimax()) - @test isapprox(loadings(FR), F, rtol = √eps(Float64)) - @test isapprox(rotation(FR), R, rtol = √eps(Float64)) + @test isapprox(loadings(FR), F, rtol=√eps(Float64)) + @test isapprox(rotation(FR), R, rtol=√eps(Float64)) ## Equamax # Comparison with Matlab's `rotatefactors` called as # using MATLAB # F = mat"rotatefactors($X, 'Method', 'equamax', 'Normalize', 'off')" - F = [ 0.90503611 1.24332782; - -2.32641022 0.36067319; - 0.97744298 1.36280122; - 0.85650467 0.23818242; - 0.45865486 -2.40335769] + F = [0.90503611 1.24332782 + -2.32641022 0.36067319 + 0.97744298 1.36280122 + 0.85650467 0.23818242 + 0.45865486 -2.40335769] # Equamax for d x p matrix means orthogonal Crawford-Ferguson with κ = p / (2 * d) - FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = p / (2 * d))) - @test isapprox(loadings(FR), F, rtol = √eps(Float64)) + FR = rotate(X, CrawfordFerguson{Orthogonal}(kappa=p / (2 * d))) + @test isapprox(loadings(FR), F, rtol=√eps(Float64)) ## Parsimax @@ -94,16 +94,16 @@ using StableRNGs # using MATLAB # F = mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')" # ``` - - F = [ 0.90503611 1.24332782; - -2.32641022 0.36067319; - 0.97744298 1.36280122; - 0.85650467 0.23818242; - 0.45865486 -2.40335769] + + F = [0.90503611 1.24332782 + -2.32641022 0.36067319 + 0.97744298 1.36280122 + 0.85650467 0.23818242 + 0.45865486 -2.40335769] # Parsimax for d x p matrix means orthogonal Crawford-Ferguson with κ = (p - 1) / (d + p - 2) - FR = rotate(X, CrawfordFerguson{Orthogonal}(κ = (p - 1) / (d + p - 2))) - @test isapprox(loadings(FR), F, rtol = √eps(Float64)) + FR = rotate(X, CrawfordFerguson{Orthogonal}(kappa=(p - 1) / (d + p - 2))) + @test isapprox(loadings(FR), F, rtol=√eps(Float64)) ## Quartimin rotation # Comparison with R's `GPArotation::quartimin` function called as @@ -116,21 +116,21 @@ using StableRNGs # @rget R # ``` - F = [ 0.94295134 1.29282920 - -2.32274063 0.24011390 - 1.01896033 1.41630052 - 0.86570192 0.28326547 - 0.39161552 -2.38397278] - R = [ 0.87548617 0.41144611 - -0.48324317 0.91143409] + F = [0.94295134 1.29282920 + -2.32274063 0.24011390 + 1.01896033 1.41630052 + 0.86570192 0.28326547 + 0.39161552 -2.38397278] + R = [0.87548617 0.41144611 + -0.48324317 0.91143409] FR = rotate(X, Quartimin()) - @test isapprox(loadings(FR), F, rtol = √eps(Float64)) - @test isapprox(rotation(FR), R, rtol = √eps(Float64)) + @test isapprox(loadings(FR), F, rtol=√eps(Float64)) + @test isapprox(rotation(FR), R, rtol=√eps(Float64)) # Test application to factor analysis and PCA models - X = randn(10, 5) + X = randn(rng, 10, 5) M = fit(FactorAnalysis, X) loadings(M) diff --git a/test/runtests.jl b/test/runtests.jl index fc4a0fb..8345cb2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,16 +1,16 @@ -tests = ["lreg", - "whiten", - "pca", - "cca", - "cmds", - "lda", - "mclda", - "ica", - "ppca", - "kpca", - "fa", - "facrot"] +tests = [ #"lreg", + # "whiten", + # "pca", + # "cca", + # "cmds", + # "lda", + # "mclda", + # "ica", + # "ppca", + # "kpca", + # "fa", + "facrot"] for test in tests - include(test*".jl") + include(test * ".jl") end From fc49ff78cc070c51329bf158d34e65af9c97b672 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Fri, 19 Jan 2024 11:42:18 +0100 Subject: [PATCH 30/32] version dependent import --- src/MultivariateStats.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/MultivariateStats.jl b/src/MultivariateStats.jl index 41191dd..ecef985 100644 --- a/src/MultivariateStats.jl +++ b/src/MultivariateStats.jl @@ -10,7 +10,11 @@ module MultivariateStats import Statistics: mean, var, cov, covm, cor import Base: length, size, show import StatsAPI: fit, predict, coef, weights, dof, r2 - import LinearAlgebra: eigvals, eigvecs, rotate! + import LinearAlgebra: eigvals, eigvecs + + @static if VERSION ≥ v"1.5" + import LinearAlgebra: rotate! + end export From 9eca2dca3021eb42b23c5459d47822aa38b9ad80 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Fri, 19 Jan 2024 11:58:37 +0100 Subject: [PATCH 31/32] convert to old style keyword arguments --- src/facrot.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/facrot.jl b/src/facrot.jl index 7166db9..764f5bb 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -500,7 +500,7 @@ function rotate(L::AbstractMatrix, lsiter::Integer=10, tol::Real=1.0e-6) where {T<:FactorRotationMethod} return FactorRotation( - gparotate(L, C; normalizerows, randominit, maxiter, lsiter, ϵ=tol)... + gparotate(L, C; normalizerows=normalizerows, randominit=randominit, maxiter=maxiter, lsiter=lsiter, ϵ=tol)... ) end @@ -517,7 +517,7 @@ function rotate!(M::FactorAnalysis, maxiter::Integer=1000, lsiter::Integer=10, tol::Real=1.0e-6) where {T<:FactorRotationMethod} - FR = rotate(loadings(M), C; normalizerows, randominit, maxiter, lsiter, tol) + FR = rotate(loadings(M), C; normalizerows=normalizerows, randominit=randominit, maxiter=maxiter, lsiter=lsiter, tol=tol) M.W .= loadings(FR) return M @@ -536,7 +536,7 @@ function rotate!(M::PCA, maxiter::Integer=1000, lsiter::Integer=10, tol::Real=1.0e-6) where {T<:FactorRotationMethod} - FR = rotate(projection(M), C; normalizerows, randominit, maxiter, lsiter, tol) + FR = rotate(projection(M), C; normalizerows=normalizerows, randominit=randominit, maxiter=maxiter, lsiter=lsiter, tol=tol) M.proj .= loadings(FR) return M From 21c6a3a70f151b6027ed6d928b0613293646b857 Mon Sep 17 00:00:00 2001 From: Felix Held Date: Fri, 19 Jan 2024 13:09:16 +0100 Subject: [PATCH 32/32] improved documentation --- docs/make.jl | 2 +- docs/src/facrot.md | 85 ++++++++++++++++++++++++++++++++++++++++++++-- docs/src/index.md | 2 +- src/facrot.jl | 10 ++++-- 4 files changed, 93 insertions(+), 6 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index a95c1f6..fdc42ec 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,8 +15,8 @@ makedocs( "ica.md", "cca.md", "fa.md", - "mds.md", "facrot.md", + "mds.md", "Development" => "api.md"] ) diff --git a/docs/src/facrot.md b/docs/src/facrot.md index c3771b4..a02b69a 100644 --- a/docs/src/facrot.md +++ b/docs/src/facrot.md @@ -1,4 +1,85 @@ # Factor Rotation -Factor rotations are used for post-processing of a factor matrix determined -from a [Factor Analysis](../src/fa.md) or [Principal Component Analysis](../src/pca.md). \ No newline at end of file +[Factor rotations](https://en.wikipedia.org/wiki/Factor_analysis#Rotation_methods) +are used for post-processing of a factor matrix obtained by +[`FactorAnalysis`](@ref) or [`PCA`](@ref). + +Let ``L`` be a matrix of factor loadings. Two common methods of factor rotation[^1] +are typically considered: +- A matrix ``\Lambda`` is an *orthogonal rotation*[^2] of ``L`` if + ``\Lambda = L R`` for some orthogonal matrix ``R``, i.e. ``R^{-1} = R^\top``. +- A matrix ``\Lambda`` is an *oblique rotation*[^3] of ``L`` if + ``\Lambda = L (R')^{-1}`` for some matrix ``R`` that has columns of length 1. + +The matrix ``R`` is chosen by minimizing a factor rotation criterion +``Q(\Lambda)``. This package uses a gradient projection algorithm[^1] +to minimize the criterion ``Q``. + +## Applying factor rotations + +This package defines the function [`rotate`](@ref) which can be used to rotate +a matrix of loadings. + +```@docs +rotate(::AbstractMatrix, ::FactorRotationCriterion{T}) where {T<:FactorRotationMethod} +FactorRotation +loadings +rotation +``` + +[`rotate!`](@ref) modifies the supplied model in place. + +```@docs +rotate!(::FactorAnalysis, ::FactorRotationCriterion{T}) where {T<:FactorRotationMethod} +rotate!(::PCA, ::FactorRotationCriterion{T}) where {T<:FactorRotationMethod} +``` + +## Factor rotation methods + +There are two types of factor rotations: *orthogonal* and *oblique*. +Which type of rotation method should be used can be chosen by setting the +appropriate subtype of [`FactorRotationMethod`](@ref). + +```@docs +FactorRotationMethod +Orthogonal +Oblique +``` + +## Factor rotation criteria + +The actual factor rotation criteria are described by subtypes of +[`FactorRotationCriterion`](@ref) and many parametrized forms of standard +rotation methods are provided. + +```@docs +FactorRotationCriterion{T} where {T<:FactorRotationMethod} +``` + +Sometimes, factor rotation criteria can be used as orthogonal and oblique methods. + +```@docs +CrawfordFerguson{T} where {T<:FactorRotationMethod} +Oblimin{T} where {T<:FactorRotationMethod} +``` + +### Orthogonal factor rotation criteria + +```@docs +Varimax +Quartimax +MinimumEntropy +``` + +### Oblique factor rotation criteria + +```@docs +Quartimin +``` + +## References + +[^1] Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms and Software for Arbitrary Rotation Criteria in Factor Analysis. Educational and Psychological Measurement, 65, 676-696. doi [10.1177/0013164404272507](https://doi.org/10.1177/0013164404272507) +[^2] Jennrich, R. I. (2001). A simple general procedure for orthogonal rotation. Psychometrika, 66, 289-306. doi [10.1007/BF02294840](https://doi.org/10.1007/BF02294840) +[^3] Jennrich, R. I. (2002). A simple general method for oblique rotation. Psychometrika, 67, 7-19. doi [10.1007/BF02294706](https://doi.org/10.1007/BF02294706) + diff --git a/docs/src/index.md b/docs/src/index.md index a45f2e8..0003acc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -11,7 +11,7 @@ end [MultivariateStats.jl](https://github.com/JuliaStats/MultivariateStats.jl) is a Julia package for multivariate statistical analysis. It provides a rich set of useful analysis techniques, such as PCA, CCA, LDA, ICA, etc. ```@contents -Pages = ["whiten.md", "lreg.md", "lda.md", "pca.md", "ica.md", "cca.md", "fa.md", "mds.md", "api.md"] +Pages = ["whiten.md", "lreg.md", "lda.md", "pca.md", "ica.md", "cca.md", "fa.md", "facrot.md", "mds.md", "api.md"] Depth = 2 ``` diff --git a/src/facrot.jl b/src/facrot.jl index 764f5bb..4f9d0ef 100644 --- a/src/facrot.jl +++ b/src/facrot.jl @@ -265,6 +265,8 @@ to perform the computation line search - `ϵ` is the convergence tolerance +Returns the rotated loading matrix and the rotation matrix as a tuple. + **References** - Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms and Software for Arbitrary Rotation Criteria in Factor Analysis. @@ -366,6 +368,8 @@ to perform the computation line search - `ϵ` is the convergence tolerance +Returns the rotated loading matrix and the rotation matrix as a tuple. + **References** - Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms and Software for Arbitrary Rotation Criteria in Factor Analysis. @@ -481,6 +485,8 @@ Rotate the loadings in matrix `L` using the criterion `C`. - `L` is the matrix of loadings to be rotated - `C` is a factor rotation criterion +Returns a [`FactorRotation`](@ref) object. + **Keyword parameters** - If `normalizerows` is true, then the rows of `F` are normalized to length 1 before rotation and the rows of the rotated loadings are @@ -507,7 +513,7 @@ end """ rotate(M::FactorAnalysis, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) -Rotate the loadings of the Factor Analysis model `M` using the criterion `C`. +Rotate the loadings of the [`FactorAnalysis`](@ref) model `M` using the criterion `C`. Modifies the loading matrix in `M`. Parameters as in [`rotate`](@ref). """ function rotate!(M::FactorAnalysis, @@ -526,7 +532,7 @@ end """ rotate(M::PCA, C::FactorRotationCriterion{T <: FactorRotationMethod}; ...) -Rotate the components of the PCA model `M` using the criterion `C`. +Rotate the components of the [`PCA`](@ref) model `M` using the criterion `C`. Modifies the projection matrix in `M`. Parameters as in [`rotate`](@ref). """ function rotate!(M::PCA,