Add a 2-arg flavor, argmax(f, domain)

[drewrobson]

With the rename to argmax, and dims now a kwarg, it seems very tempting to me to define a 2-arg flavor, argmax(f, itr), where itr is an iterable subset of the domain of f. This would return an x in itr that maximizes f(x), which is quite literally the mathematical definition of argmax.

(originally proposed here)

Edit: Previous wording may have been confusing, so I have attempted to improve it.

[drewrobson]

If we have #27612, then this is as simple as defining the following:

Base.argmax(f, itr) = Base.argmax(f(x) for x in itr)

(but I haven't checked what to do about dims yet)

[ararslan]

A more correct definition would be to adapt the current definitions of findmin and findmax to accept a function argument, then simply adjust the current argmax definition to

Base.argmax(f, x::AbstractArray; dims=:) = findmax(f, x, dims=dims)[2]
Base.argmax(f, itr) = findmax(f, itr)[2]

with the addition of

Base.argmax(x::AbstractArray; dims=:) = findmax(identity, x, dims=dims)[2]
Base.argmax(itr) = findmax(identity, itr)[2]

and likewise for argmin. Whether generators are supported is somewhat orthogonal here.

[drewrobson]

Thanks very much for taking a look at this. I think you're right. Would the compiler be able to inline the identity function so there is no cost to this change in the most commonly used case?

[ararslan]

Yep. This is the same pattern we use in a lot of places where we define methods that accept function arguments.

[Nosferican]

Are there any other functions that could be prime candidates to get the function argument improvement? Besides, argmin, argmax... extrema?

[drewrobson]

Wait, I think there are some subtleties in reconciling @ararslan's nice suggestion with my original intent. Let me try to explain, and please forgive me if I confuse things further (e.g. I'm not entirely sure how to talk about the indices of an AbstractArray, but I gave it a shot!)

One of the difficulties in thinking about this is that there are really two flavors of argmax, the one defined in array.jl and the one defined in reducedim.jl.

In reducedim.jl, argmax defines an operation on A::AbstractArray. Here, I think the domain is implicitly CartesianIndices(axes(A)) and the function is implicitly getindex(A, i), where i ∈ CartesianIndices(axes(A)). The way I was thinking about it, argmax(A) is a convenient shorthand for something like argmax(i -> getindex(A, i), CartesianIndices(axes(A))) and I am proposing to generalize this to allow argmax(f, itr) for any function f and collection itr where this makes sense. In the special case that itr::AbstractArray then this would be implemented in reducedim.jl and support the dims kwarg and otherwise we have a fallback function in array.jl as long as keys(itr) is defined.

I think @ararslan's suggestion is distinct from this. If I understand correctly, for A::AbstractArray the domain would still implicitly be CartesianIndices(axes(A)) and the result of the new 2-arg flavor argmax(f, A) would not be an element of A (as proposed above), but rather an index of A. But please, correct me if I have misunderstood!

[ararslan]

We may be talking past each other, but I'll try to clarify.

For any iterable x, argmax(x) always returns an index into x, not an element of x. That can be seen more clearly from its maiden name, indmax. To ensure consistency with how the function is supposed to work, argmax(f, x) should return the index i at which f(x[i]) is maximized, not the element x[i]. To get the element, just index this back into x, e.g. x[argmax(f, x)].

[jekbradbury]

I think @ararslan’s is the standard Julian extension of a 1-argument function on collections to the case with an additional function argument. But @drewrobson is proposing the traditional 2-arg mathematical argmax, which is less consistent with the existing argmax but does coexist with it in math notation-land?

[drewrobson]

Yes exactly:

I completely agree that with the former name, it was clear that indmax should apply to things that can be indexed, and it should return the index that maximizes that thing. I'm also 100 % with you that if we want to stick with that original meaning even after the rename to argmax, then your suggestion is an undoubtably useful way to add functionality that is consistent with the original meaning of indmax.

With the rename to argmax, I am raising the (perhaps heretical!) suggestion that we should change what we mean by this function. To me, argmax now applies most naturally to a function, and it should return the arg that maximizes that function. As a special case, we would preserve the existing 1-arg flavor for things that can be indexed, with the understanding that getindex is the function to be maximized.

In summary, I think this is a backwards compatible and satisfying reinterpretation of what this function should do, but I get that it feels like a radical proposal!

[ararslan]

Ohhh okay, I see now, thanks @jekbradbury and @drewrobson for clarifying.

An unfortunate consequence of the renaming of indmax to argmax is indeed that the latter conflicts a bit with the identically named function from math when it comes to applying a function to the given collection. If we're to add a function argument to argmax I think we should ensure that the behavior is consistent with that of other Julia functions that accept function arguments, in particular that g(f, x) is equivalent to g(map(f, x)). I realize that's rather unsatisfying in this case given the more general name, but I really think making the cases behave differently will cause more confusion in the long run.

[drewrobson]

Yes, avoiding confusion is important. But then, perhaps argmax should only be defined in the 2-arg form proposed above, and then indmax should be defined as indmax(A) = argmax(i -> getindex(A, i), CartesianIndices(axes(A))) (or however you write that). Those names seem closest to what each function does, but I don't know if you'd want to split it up like that.

Edit: @Nosferican's post has reminded me that in addition to the 2-arg form proposed above, we would also keep the 1-arg argmax defined in array.jl for any itr for which keys/pairs is defined. Which will include generators, if #27612 happens!

[Nosferican]

• argmax is returning the index rather than the definition due to a historical accident (renaming of indmax to argmax) which should be fixed ASAP
• argmax works only on pairs-implemented collections due to the implementation which restricts the function to only a subset of all cases included in the mathematical definition (should be fixed)
• A 2-args method for passing a non identity function should be added
  My very naive implementation complies with the above,

iter = Set(1:3)
function my_argmax(f::Function, iter)
    next = iterate(iter)
    x = next[1]
    y = f(x)
    while next !== nothing
        (i, state) = next
        yᵢ = f(i)
        if y < yᵢ
            x = i[1]
        end
        next = iterate(iter, state)
    end
    return x
end

My ignorant self doesn't know what is the appropriate abstract type for an iterate.

[drewrobson]

@Nosferican, note that with #27612 this becomes a 1-liner. In other words, first define

Base.keys(g::Base.Generator) = g.iter

and then define

argmax(f::Function, itr) = argmax(f(x) for x in itr)

[ararslan]

argmax is returning the index rather than the definition due to a historical accident (renaming of indmax to argmax) which should be fixed ASAP

Can you elaborate on what you mean by that?

[Nosferican]

An unfortunate consequence of the renaming of indmax to argmax

I may have misinterpreted the above as historical background on why argmax is returning the index rather than the value (it seems that indmax was deprecated in favor of argmax). Per the mathematical definition, the code should return the first element of findmax and not the second (value not the index). It seems to be a classic case of renaming the function to indicate a different concept, but keeping the implementation with the original concept. I would consider this a bug. Even the docstring for argmax uses the definition of indmax instead of one describing the mathematical function.

[drewrobson]

@Nosferican, I think you're talking about the mathematical definition of max, not argmax. The debate here is whether argmax should return an arg or an index, not whether it should return the max value!

[Nosferican]

Shouldn't be multi-taking. Sorry for the confusion (I suffered it the worst); indmax does return the value and position. The original point before getting lost somewhere, was that it shouldn't be the index of the solution, but the actual value in the domain.

itr = [1, 4, 3]
argmax(abs2, itr) == 4 # not 2

rather than the current behavior which is indmax: argmax(itr) == 2. That way you can support iterate which do not have pairs and it follows the actual mathematical definition: which argument in the domain gives the maximum value in the range.

[drewrobson]

There we go. Yes, that is consistent with my original proposal above.

[ararslan]

I'd be interested to hear @stevengj's thoughts on this.

[drewrobson]

And perhaps @JeffBezanson and @StefanKarpinski as well? Stefan encouraged me to look into this on Discourse and I know Jeff has seen my other issue at least. Thanks, and sorry to rock the boat! I mean well.

[Nosferican]

I think the bug was an oversight when renaming the function from indmax to argmax (two different mathematical concepts). @drewrobson, thanks for bringing that up. Wouldn't have wanted to have this wrong in a release.

[stevengj]

I don't see a problem with the single-argument argmax(a) returning the key of the maximum; it seems pretty natural to view an index container a as a map from keys k to values a[k].

[drewrobson]

I like that. That is consistent with what I was trying to say here:

"As a special case, we would preserve the existing 1-arg flavor for things that can be indexed, with the understanding that getindex is the function to be maximized."

But the remaining question is, what should be the result of the 2-arg form:

itr = [1, 4, 3]
argmax(abs2, itr)

I think @Nosferican and I are arguing that it should be 4, the arg / domain value in itr that maximizes abs2, whereas if we want argmax(abs2, itr) to be the same as argmax(map(abs2, itr)) then it should be 2, the index in itr that maximizes abs2.(itr).

Edit: If we go down the latter road then I think we'd want separate verbs indmax and argmax to avoid the confusion that @ararslan warned about. indmax(f, itr) would do what @ararslan proposed and argmax(f, itr) would do the mathematical operation proposed in this issue.

[Nosferican]

indmax should have not been deprecated for argmax as these are two different concepts. argmax as a mathematical concept is defined for sets. The actual definition for argmax doesn't deal with indices at all because it is a subset function. It subsets a set to include only elements with mappings to the maximum of the range.

obj = Set(-3:3)
argmax(abs2, obj) == Set([-3, 3])

Long-hand form,

function argmax(f::Function, domain::Any)
    output = [ (x, f(x)) for x ∈ domain ] # Generates the domain and mapping
    frontier = maximum(last.(output)) # Obtains the frontier
    output = filter(x -> isequal(last(x), frontier), output) # subsets domain
    output = Set(first.(output)) # return the elements which meet the criteria
    return output
end

Refresher: Wikipedia

[drewrobson]

@Nosferican is right to bring up the issue of ties. I have to sheepishly admit I was trying to avoid that because for efficiency reasons I wanted to do the same thing as indmax and just return an answer for one maximum. Honestly, if this constructs a Set I'd probably go back to using indmax and then indexing into my domain / collection. Hmm.

[Nosferican]

A compromise I would be willing to accept is for ordered collections to return a Vector

function argmax(f::Function, domain::AbstractVector)
    range_of_domain = f.(domain)
    frontier = maximum(range_of_domain)
    return domain[range_of_domain .== frontier]
end

For Dict, the proper way would be call argmax(f, keys(dict)) or argmax(f, values(dict)) to avoid ambiguity.

[drewrobson]

I really like the current behavior of argmax(d::Dict), which returns one key that maximizes d. I suspect in most real life use cases, the efficiency of doing a single pass and keeping the best result is what a programmer wants. Admittedly, this is a slight departure from the rigorous mathematical definition, but code has to live in the real world! And even in math, in casual usage I'm pretty sure I've seen "let x = argmax_{x} f(x)" when technically one should have written "let x ∈ argmax_{x} f(x)", because sometimes one isn't picky about which x as long as it maximizes f(x) and so the fact that there may be more than one to pick from is just a nuisance.

[Nosferican]

When in doubt, drop a keyword argument... If one only cares about a single element you could implement the easy/efficient way and the faithful mathematical function when that is what you want.
You can probably still get the short and fast behavior by calling itr[last(findmax(f, itr))] (assuming findmax is getting the two-args flavor.

[drewrobson]

I guess, but together with dims you end up with things like (N-1) dimensional Sets. We're well past the use case that I had in mind when I created this issue. I dunno, hopefully others will weigh in tomorrow!

[tkf]

Using isless simplifies a lot for supporting NaN and missing for maximum-like reductions. As we need to handle missing anyway, I don't think < simplifies minimum-like reductions.

Ultimately, I think it'd be better to create a function (say) isgreater and move the complicated logic of min to it. If we use isgreater from min and findmin, the consistency is automatic.

missing probably has to be handled explicitly

Shouldn't we aim for "implicit" handling? It's consistent with the policy to be reluctant to add missing overload for Base and stdlib functions.

Perhaps a case where people should just explicitly use skipmissing if they need to deal with missing values?

We already have minimum([-Inf, NaN, missing]) === missing. So I don't think we can change how it's handled?

[cmcaine]

@StefanKarpinski Thinking about this some more, I agree with @tkf that < doesn't give what we want on its own.

We want NaN to poison, but it just doesn't dominate all other values (NaN is not less than any other value, but other values aren't less than NaN either), so they will not be consistently chosen by any rule of the type "not less than any other value".

I also belatedly realised that floats are not partially ordered because ≤ isn't reflexive for NaN. Woops.

Here's a quick worked example of why < on its own isn't sufficient:

# In findmax somewhere:

if fm < fx
    fm, m = fx, x
end
# true if fm is less than fx
#       good
# false if fm >= fx
#       good
# false if fm is NaN
#       good
# false if fx is NaN
#       not good: we won't switch from a non NaN to a NaN

---

The existing min and max functions all define a preference order for poison values. I think it's a shame that they mostly use isless for that because it only works nicely in one direction and you don't need a total ordering to have poisoning happen, you only need to define an ordering of the poisons relative to normal values and, if you care about which poison you receive, an ordering for the poisons.

I think rather than introduce an isgreater function, it might be better to define poisoning lessthan and greaterthan operations (or a poisoning order HOF).

is_unorderable(x) = !isa(x ≤ x, Bool) || !(x ≤ x)

poison_lt(a, b) = is_unorderable(a) || is_unorderable(b) ? isless(a, b) : a < b
poison_gt(a, b) = is_unorderable(a) || is_unorderable(b) ? isless(a, b) : a > b

# or a higher order function, whatever
poisoning_order(op) = (a, b) -> is_unorderable(a) || is_unorderable(b) ? isless(a, b) : op(a, b)

That way, partial ordering is permitted and we don't introduce a new canonical total ordering into Julia. We could then use these in all of the min and max functions.

If that's too big a change or otherwise undesirable, then I think we should define isgreater. If we allow ourselves to order normally unorderable values in an arbitrary but consistent fashion (as with isless) and put missing last, then I think this definition is fine and satisfies our need for all the weird values to end up at the end.

isgreater(a, b) = is_unorderable(a) || is_unorderable(b) ? isless(a, b) : isless(b, a)

Anyway, depending on what people are willing to accept, the implementation is either this:

findmin(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmin, domain)
_rf_findmin((fm, m), (fx, x)) = !poison_gt(fm, fx) ? (fm, m) : (fx, x)

findmax(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmax, domain)
_rf_findmax((fm, m), (fx, x)) = !poison_lt(fm, fx) ? (fm, m) : (fx, x)

or this:

findmin(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmin, domain)
_rf_findmin((fm, m), (fx, x)) = isless(fm, fx) ? (fm, m) : (fx, x)

findmax(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmax, domain)
_rf_findmax((fm, m), (fx, x)) = isgreater(fm, fx) ? (fm, m) : (fx, x)

Edit: these definitions of findmax return identical results except for types like Sets, for those the first definitions find the first value that is not smaller than all other results. The definition with isgreater just gives an error.

[tkf]

OK, I agree that there is a way to use < (so the discussion on partial/total ordering matters). But why prefer partial ordering? I think important questions are:

1. Do we want maximum, argmax, and findmax to be consistent?
2. Do we want to return the first maximal element from maximum? (It does not seem to be consistent with the definition in Greatest and least elements - Wikipedia.)

I see that findmax docstring is mentioning that it returns the first maximal element. However, there is no multiple maximal element support because it uses isless in the implementation. (OK, this is technically true if you run findmax on floats and check the identity with ===. But other than this case, I don't think partial ordering is supported.)

[cmcaine]

1. Do we want maximum, argmax, and findmax to be consistent?

I vote yes. If we support partial order in one, we should support it in all of the others. So, perhaps that is beyond the scope of this issue and I should just submit a PR that adds the 2-arg versions of arg* and find* and another that changes all of the min/max functions to support types with partial orders.

2. Do we want to return the first maximal element from maximum? (It does not seem to be consistent with the definition in Greatest and least elements - Wikipedia.)

From that wikipedia article:

The greatest element of a partially ordered subset must not be confused with maximal elements of the set, which are elements that are not smaller than any other elements. A set can have several maximal elements without having a greatest element. However, if it has a greatest element, it can't have any other maximal element.

I intended to follow that definition of maximal and I think I did?

Why prefer partial ordering?

Supporting partial orders makes no difference to the existing supported uses (types with total orderings) and provides new functionality that seems reasonable for types that only have partial orderings.

I don't have a use case for this, though. Maybe it's something that nobody actually needs.

[cmcaine]

I think the partial vs total order thing comes down to this, do you think this should throw an error? If not, what do you think it should return?

min(Set(1), Set(2), Set(1:3))

If it should return something, then so too should all the other functions we're talking about.

[tkf]

I intended to follow that definition of maximal and I think I did?

Yes, but I think that's the problem; i.e., maximum returns a maximal.

In other words, a maximal is not (necessary) the maximum in the definition mentioned in Wikipedia (which I think is pretty standard). So, if you take the definition literally, returning a maximal from a function named maximum is confusing.

[cmcaine]

Sorry for the misunderstanding! I wondered if we should define a maximal function, but I thought it might become confusing.

In any case, if there is a maximum or greatest element, then it is also the only maximal. If there are multiple maximals, then we can (and already have, for findmax) just say that we'll give the first one.

I think a line like "If there are multiple maximal elements, then the first one is returned" is pretty clear and better than these functions just erroring on partially ordered types, but I can see your perspective too: we are no longer guaranteeing that the element is the maximum element (though it sort of wasn't before because of NaN and missing anyway).

My intention is to submit a PR that adds the 2-arg variants of find* and arg* and ask for triage on whether these functions should error or return the first maximal for types that can only be partially ordered. Does that sound good to you?

[tkf]

I think a line like "If there are multiple maximal elements, then the first one is returned" is pretty clear and better than these functions just erroring on partially ordered types

Yes, I can see that a function like this may be useful sometimes. At the same time, I think having a value guaranteed to be the maximum is also useful, too (e.g., using it as a sentinel value).

it sort of wasn't before because of NaN and missing anyway

Current behavior of findmax for NaN is the correct maximum w.r.t isless. This is why it is so easy to implement findmax when using isless. This is not the case for findmin, though (which is why I mentioned isgreater).

By the way, another problem for <-based approach is:

julia> maximum([-0.0, 0.0])
0.0

julia> foldl([-0.0, 0.0]) do a, b
           a < b ? b : a
       end
-0.0

[cmcaine]

Current behavior of findmax for NaN is the correct maximum w.r.t isless. This is why it is so easy to implement findmax when using isless. This is not the case for findmin, though (which is why I mentioned isgreater).

isless doesn't guarantee in it's documentation that NaN will be the maximum, or that unorderable elements are bigger than orderable elements, those are extra guarantees that are specified for some of the min/max functions. We could just rewrite the docs to guarantee that, though.

Maybe the original sin here is that I believed the docstring for <. That says < implements a partial order, but that's just not true (1. It's not reflexive even for integers; 2. <= is a partial order for lots of stuff, but not for Floats because of NaN).

I think that that probably merits a documentation change, and, if we wanted to take partial orders seriously then there should probably be a new function for the default partial order relation on a type (and we'd need to decide what it should do for unorderables and missings: error or put them at the end).

[float zeroes are a problem]

Thanks for pointing that out.-0.0 == 0.0, so that's fine in a partial order, but it's plausible that someone is relying on the sign of the result and that would make this a breaking change, contrary to my statement above :( It would be fairly easy to fix this case, but I'm not sure that a partial-ish order is a very good idea.

I basically just need to finalise the docstrings for this and check that isgreater isn't hilariously wrong and I'll submit the PR for review, probably tomorrow. I'll be editing the docstrings of isless, isgreater, findmax/min, argmax/min, <, and maybe, but probably not minimum, min, maximum, max.

[StefanKarpinski]

We want NaN to poison, but it just doesn't dominate all other values (NaN is not less than any other value, but other values aren't less than NaN either), so they will not be consistently chosen by any rule of the type "not less than any other value".

I haven't read all the rest here, but that's exactly why the condition to update the max is that the max is less than the new value: because NaN < 0.0 and NaN > 0.0 are both false. I demonstrated the correct NaN poisoning behavior already with that logic. I also demonstrated the correct behavior with partially ordered types like sets.

[cmcaine]

The issue is that if the current max is say 3.0, then we won't switch to NaN because 3.0 < NaN is false.

[cmcaine]

If floats were partially ordered it would be fine, but they're not, NaNs are unordered and will only be picked as the max if they're the first item in the iterable.

[tkf]

isless doesn't guarantee in it's documentation that NaN will be the maximum, or that unorderable elements are bigger than orderable elements, those are extra guarantees that are specified for some of the min/max functions.

That's a fair point. isless on concrete floats are defined in Base so I'd say it's OK to exploit the implementation details inside Base. But users can define their own AbstractFloat subtypes. So I agree it's better to document ordering of NaN for floats w.r.t isless. It's documented for missing so I think it's reasonable.

[float zeroes are a problem]

Thanks for pointing that out.-0.0 == 0.0, so that's fine in a partial order, but it's plausible that someone is relying on the sign of the result and that would make this a breaking change

Yeah, it's true that this behavior is not in the docstring but it doesn't mean changing it is OK. (In a way "true" SemVer is a fantasy...)

It is already documented that isequal(-0.0, 0.0) is false and isequal and isless are compatible. So, documenting isless(-0.0, 0.0) === !isless(0.0, -0.0) === true seems to not hurt anything. Once it is documented that maximum-like reductions are using isless as-is, I suppose there is no ambiguity in the specification any more for Union{T,Missing} where T is one of concrete AbstractFloat subtypes defined in Base?

[cmcaine]

I think it's defined somewhere or other that negative zero isless positive zero. If I can't find where then I'll add it to the docstring in implementer notes or something.

[cmcaine]

I finally submitted that PR. It's here: #35316. 🎉