Builtin argmin is slower than manual

I don’t understand how this is possible. I tried it many times, targmin always wins, and by quite a significant margin (1.5.2)

using BenchmarkTools

function targmin(a)
    m = Inf
    ind = 0
    for i in 1:length(a)
        if a[i] < m
            m = a[i]
            ind = i
        end
    end
    ind
end

a = rand(10 ^ 8);

println(argmin(a) == targmin(a))

@btime argmin(a)
@btime targmin(a)

true
382.492 ms (1 allocation: 16 bytes)
182.026 ms (1 allocation: 16 bytes)

I can replicate this on master. It would be great to get to the bottom of this — ideally with a fix, but at least opening an issue.

Note that your code assumes that eltype(a) is comparable to Float64, since m is Inf in the first iteration. That might not necessarily be the case. You’re also assuming 1-based indexing, eachindex would be better.

That said, I think the crux of the matter is that argmin does findmin, which itself calls findminmax! (seemingly very complex function)

Additionally, argmin allows for reduction across arbitrary dimensions, your targmin always reduces over the whole array.

1 Like

The core of argmin is:

It does some extra stuff to handle NaN etc and it is generic on the type of the array. But maybe it could be optimized as well (the code in question is quite old), for example, is both the m != m and ai != ai needed in every loop iteration?

Oh wow, that code is ages old!

Those checks are needed to ensure the same behaviour as min and max for NaN, I think.

Yeah, I commented on the NaN handling but it does for exampe m != m over and over in the loop, even if m hasn’t changed.

I wonder why the “change in iteration protocol” commit didn’t change to the new for loop as well and instead did it manually?

Probably because you want to peel off the first value to use as the initial minimum.

How about:

findmax(a) = _findmax(a, :)

function _findmax(a, ::Colon)
    p = pairs(a)
    y = @inbounds iterate(p)
    if y === nothing
        throw(ArgumentError("collection must be non-empty"))
    end
    (mi, m), s = y
    isnan(m) && return m, mi
    i = mi
    while true
        y = @inbounds iterate(p, s)
        y === nothing && break
        (i, ai), s = y
        isnan(ai) && return ai, i
        # Neither `m` nor `ai` can be NaN here so can use `<` instead of `isless`.
        # Edit: not true, consider 0.0 < -0-0 == false
        if m < ai
            m, mi = ai, i
        end
    end
    return m, mi
end

does anyone see any problems with that?

Edit, the isnan should be changed back to ai != ai to handle missing…
Edit Edit: actually, findmax is already broken with respect to missing

julia> findmax([missing, 1])
ERROR: TypeError: non-boolean (Missing) used in boolean context

Edit Edit Edit: My code fails to handle signed zero I think.

1 Like

You would not need to bail out early at all, isless already treats NaN the way you want

          while true
               y = iterate(p, s)
               y === nothing && break
               (i, ai), s = y
               if isless(m, ai)
                   m = ai
                   mi = i
               end
           end
julia> findmax([1, NaN,  2])
(NaN, 2)

Not sure if you should have a branch just for finding NaNs faster at all.

But early exit is nice in case you have a NaN at the beginning, or?

There are many early exits you could imagine, e.g. typemax(T) where T <: Union{Int32, Int64, Int128} etc., the NaN early exit seems to be of less use. And here luckily

julia> isless(NaN, NaN)
false

so you actually find the first NaN as promised.

2 Likes

Yeah, that’s true.

PS: This behaviour is also not very convincing to me:

julia> struct LargerThanNaN
       end

julia> Base.isless(_, ::LargerThanNaN) = true

julia> findmax([1, NaN, LargerThanNaN()])
(NaN, 2)

3 Likes

IEEE specifies that NaN is supposed to propagate to indicate an error in an earlier computation. We’re not returning NaN because it’s the maximum (or minimum), but because it flags an illegal operation in a preceding computation.

Personally, since m < NaN is false, I would never want NaN as the result of maximum (or minimum for that matter), but that’s so 2008 of me. As of 2019, we should really be specifying whether we want maximum or maximumNumber to clarify our intent as to whether we want errors to propagate or be discarded.

How fun it is that Julia predates the latest floating point standard.