cov(A, mean=0) gives InexactError

Question

I am getting InexactError() when calculating the covariance of rows in a matrix when setting mean=0:

julia> A = [1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
5x3 Array{Int64,2}:
  1  -1  -1
 -1   1   1
 -1   1  -1
  1  -1  -1
  1  -1   1

julia> cov(A) # Works
3x3 Array{Float64,2}:
  1.2  -1.2  -0.2
 -1.2   1.2   0.2
 -0.2   0.2   1.2

julia> cov(A, mean=[0. 0. 0.]) # Works
3x3 Array{Float64,2}:
  1.25  -1.25  -0.25
 -1.25   1.25   0.25
 -0.25   0.25   1.25

julia> cov(A, mean=0) # Expecting same output as above
ERROR: InexactError()
 in generic_scale! at linalg/generic.jl:18
 in covzm at statistics.jl:253
 in cov at statistics.jl:286

Docs for Base.cov say:

mean: allow users to supply mean values that are known. By default, it is set to nothing, which indicates that the mean(s) are unknown, and the function will compute the mean. Users can use mean=0 to indicate that the input data are centered, and hence there’s no need to subtract the mean.

The Julia source (statistics.jl, line 286) clearly tests if mean=0 has been passed:

function cov(x::AbstractMatrix; vardim::Int=1, corrected::Bool=true, mean=nothing)
    mean == 0 ? covzm(x; vardim=vardim, corrected=corrected) :
    mean == nothing ? covm(x, _vmean(x, vardim); vardim=vardim, corrected=corrected) :
    isa(mean, AbstractArray) ? covm(x, mean; vardim=vardim, corrected=corrected) :
    error("Invalid value of mean.")
end

Am I doing something wrong or is this a bug?

Answer 1

I think this should be considered a bug. After having stepped through the code, the error in question occurs on line 18 of generic.jl, in particular the following function.

function generic_scale!(C::AbstractArray, X::AbstractArray, s::Number)
    length(C) == length(X) || error("C must be the same length as X")
    for i = 1:length(X)
        @inbounds C[i] = X[i]*s
    end
    C
end

where the 3 variable form is invoked by the 2 variable form like this:

generic_scale!(X::AbstractArray, s::Number) = generic_scale!(X, X, s)

At this point, from looking at the source, s will have value 0.25 and X will have the result of the following

x = unscaled_covzm(A, 1)
3x3 Array{Int64,2}:
  5  -5  -1
 -5   5   1
 -1   1   5

Cutting and pasting the above definition into the REPL allow one to test generic_scale directly

julia> generic_scale!(x,x,0.25)
ERROR: InexactError()
 in generic_scale! at none:4

And the reason is X is a Matrix of Int64 which causes convert to throw an inexact error because 1.25 can't be represented by an int. If instead X were Float64 then there is no error:

julia> x = unscaled_covzm(A, 1)*1.0
3x3 Array{Float64,2}:
  5.0  -5.0  -1.0
 -5.0   5.0   1.0
 -1.0   1.0   5.0

julia> generic_scale!(x,x,0.25)
3x3 Array{Float64,2}:
  1.25  -1.25  -0.25
 -1.25   1.25   0.25
 -0.25   0.25   1.25

And so if you define a Float64 version of A from your example

julia> B=A*1.0
5x3 Array{Float64,2}:
  1.0  -1.0  -1.0
 -1.0   1.0   1.0
 -1.0   1.0  -1.0
  1.0  -1.0  -1.0
  1.0  -1.0   1.0

julia> cov(B,mean=0)
3x3 Array{Float64,2}:
  1.25  -1.25  -0.25
 -1.25   1.25   0.25
 -0.25   0.25   1.25

Answer 2

Maybe it's Not the right answer. I just want to say, if you replace the 'mean=0' with 'mean=[0.]' in you lines, Julia get the following Answer. :)

julia> cov(A,mean=[0.])
3x3 Array{Float64,2}:
  1.25  -1.25  -0.25
 -1.25   1.25   0.25
 -0.25   0.25   1.25