DomainError for FactorAnalysis

[stelmo]

Hi,

I am trying to run factor analysis on a matrix of size 2457x19 (all float64's + quite sparse). When calling fit(FactorAnalysis, X) I get the error:

ERROR: DomainError:
sqrt will only return a complex result if called with a complex argument. Try sqrt(complex(x)).
Stacktrace:
[1] sqrt(::Float64) at .\math.jl:425
[2] #facm#61(::Int64, ::Float64, ::Int64, ::Float64, ::Function, ::Array{Float64,2}, ::Array{Float64,1}, ::Int64) at C:\Users\St. Elmo.julia\v0.6\MultivariateStats\src\fa.jl:118
[3] (::MultivariateStats.#kw##facm)(::Array{Any,1}, ::MultivariateStats.#facm, ::Array{Float64,2}, ::Array{Float64,1}, ::Int64) at .<missing>:0
[4] #fit#70(::Symbol, ::Int64, ::Void, ::Float64, ::Int64, ::Float64, ::Function, ::Type{MultivariateStats.FactorAnalysis}, ::Array{Float64,2}) at C:\Users\St. Elmo.julia\v0.6\MultivariateStats\src\fa.jl:172
[5] fit(::Type{MultivariateStats.FactorAnalysis}, ::Array{Float64,2}) at C:\Users\St. Elmo.julia\v0.6\MultivariateStats\src\fa.jl:164

I'm not doing anything fancy before running the method, so there might be a bug in the actual code. I could also supply my input matrix if you want... What can I do to fix this?

[wildart]

Either the covariance matrix is not positive definite or there is an accumulated floating-point error. If you could provide data sample, I'll look into it.

[stelmo]

Thanks for getting back to me! I just checked, the covariance matrix seems to be positive definite, so it might be accumulated floating point errors. I have attached an example matrix as a jld file but Github doesn't allow uploads of that file type so I added .zip to the file name i.e. don't try to extract it, just remove the .zip from the file name. Thanks!
example_matrix.jld.zip
.

[wildart]

It that a 19-dimensioanal data? If it is the case I think you incorrectly passed your data to a fit function.

All functions in this package accept column-major order data. So if you pass your matrix as it is (size 2457x19), you'll get 19 records of 2457 dims which definitely results in a non-positive definite covariance matrix.

If it is not the case just transpose your input.

[stelmo]

Hmm I'm not sure I understand. My input data has 19 samples (as columns) with 2457 features (as rows) and I checked for positive definiteness by doing eigvals(cov(X)) which returned a set of all positive eigenvalues (the smallest one was ~0.002).

[wildart]

This package uses data in a column-major order means which means that in a data matrix, each column corresponds to a sample/observation, while each row corresponds to a feature (variable or attribute or dimension).

If the matrix size is 2457x19 then you have a 19 2457-dimensional samples. By default, cov function computes covariance along the matrix dimension 1, cov(X, dims=1), which results in the covariance matrix of size 19x19, because row is first matrix dimension and we compute/fold along this dimension. This is fine if our samples are rows in the data matrix, which is not in our case.

So, you call cov(X, dims=2) for data in column-major order and get proper 2457x2457 covariance matrix, which in your case in not positive definite.

julia> count(eigvals(cov(X, dims = 2)) .< 0)
1223

[stelmo]

Oh dang, I misunderstood! Thanks for your help!!