rank(Any[1 2; 3 4]) fails as one(Any) is not defined

[tpapp]

julia> using LinearAlgebra

julia> rank(Any[1 2; 3 4])
ERROR: MethodError: no method matching one(::Type{Any})
Closest candidates are:
  one(::Type{Union{Missing, T}}) where T at missing.jl:105
  one(::Type{Missing}) at missing.jl:103
  one(::BitArray{2}) at bitarray.jl:400
  ...
Stacktrace:
 [1] one(::Type{Any}) at ./missing.jl:105
 [2] rank(::Array{Any,2}) at /home/tamas/src/julia-git/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/generic.jl:918
 [3] top-level scope at REPL[43]:1

julia> versioninfo()
Julia Version 1.3.0-rc4.1
Commit 8c4656b97a (2019-10-15 14:08 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_CMDSTAN_HOME = /home/tamas/src/cmdstan/

[KlausC]

Why would you expect linear algebra operations to work on matrices with non-numeric element types?
You can always do

julia> rank(Matrix{Number}(Any[1 2; 3 4]))
2

[tpapp]

I would expect them to work on all supertypes, as long as the actual elements are numbers. It actually works fine for Real, Number, or even Union{Missing,Float64} and similar.

Of course I would not expect the actual method to be fast or optimal, but having a fallback is reasonable.

I am just not sure what the best way to fix this is, which is why I opened a discussion. I am happy to make a PR then.

[StefanKarpinski]

We do generally try to make things work if the elements are of usable types, it's tough though.

[tpapp]

Is there a common idiom in Base for trying to "shrink" to a useful type and just giving up if can't? Eg

_shrink_type(iter) = mapreduce(typeof, promote_type, iter)

"Try to convert argument to a concrete element type. Second value is `false` iff failed."
function shrink_type(A::Array{T}) where T
    if isconcretetype(T)
        A, true
    else
        S = _shrink_type(A)
        if isconcretetype(S)
            S.(A), true
        else
            A, false
        end
    end
end

This should be a no-op for most cases with nice concrete types, and would still do something meaningful in other cases.

[andreasnoack]

I don't think we can make this work in general for the linear algebra routines. There are too many places where we need to know something about the element type. E.g. often tolerances are computes as function of the element type.

Is there a common idiom in Base for trying to "shrink" to a useful type and just giving up if can't?

Often people broadcast with the identity

julia> identity.(Any[1 2; 3 4])
2×2 Array{Int64,2}:
 1  2
 3  4

[tpapp]

I have reviewed the code in LinearAlgebra and I believe that a lot of use cases could be fixed by something like

_safe_one(::Type{T}) where {T <: Number} = one(T)
_safe_one(_) = 1 # fallback

and using it instead of one. This way,

1. non-numeric container element types will work (or at least, this will not be the issue preventing it),
2. uses of eps(float(one(...))) etc would default to eps(Float64), which is a sensible choice.

I am happy to do a PR, even an experimental one, but thought I would ask first if there is an a priori reason for rejecting this approach.

[andreasnoack]

Then we'd have to fundamentally change how many of the routines work since they rely on the element type of the container. E.g., the original question was related to rank and there we compute a cutoff for the singular values based on the element type. As you mention, the _safe_one fallback definition then the eps becomes eps(Float64) for an Any array insteady of failing. That's just wrong if the array contains BigFloats or Float32s. What should be the default value for the cutoff for an array that contains a Float32, BigFloat, and an Int?

Things are already a bit broken for abstractly typed matrices when they don't fail, e.g.

julia> svdvals(Matrix{Number}([1 big(2); 3 4]))
2-element Array{Float64,1}:
 5.464985704219043
 0.36596619062625746

[tpapp]

First, it isn't "wrong" in the sense that default cutoffs based on eps are in most cases already arbitrary. Eg there is no theoretical reason behind rtol = min(size(A)...)*eps(real(float(one(eltype(A)))) for rank, it just happens to be convenient in practice.

An array with element types BigFloat or Float32 would be fine.

What would break is that for some functions f, f(M) and f(Matrix{T}(M)) may not give the same result (numerically) if T is not a concrete type. But AFAIK this is not a guarantee that the API offers even for concrete T.

I would also argue that svdvals above is not broken. It is not a reasonable expectation that LinearAlgebra methods should make and effort to guess a container type when the user fails to provide one. They should just provide a reasonable fallback; if the user wants a particular intermediate type for caculation or a specific tolerance, they should just specify that explicitly.