Generic maximum/minimum function for complex numbers

Question

In julia, one can find (supposedly) efficient implementations of the min/minimum and max/maximum over collections of real numbers. As these concepts are not uniquely defined for complex numbers, I was wondering if a parametrized version of these functions was already implemented somewhere.

I am currently sorting elements of the array of interest, then taking the last value, which is as far as I know, much more costly than finding the value with the maximum absolute value (or something else).

This is mostly to reproduce the Matlab behavior of the max function over complex arrays.

Here is my current code

a = rand(ComplexF64,4)
b = sort(a,by  = (x) -> abs(x))
c = b[end]

The probable function call would look like

c = maximum/minimum(a,by=real/imag/abs/phase)

EDIT Some performance tests in Julia 1.5.3 with the provided solutions

function maxby0(f,iter)
    b = sort(iter,by  = (x) -> f(x))
    c = b[end]
end

function maxby1(f, iter)
    reduce(iter) do x, y
        f(x) > f(y) ? x : y
    end
end

function maxby2(f, iter; default = zero(eltype(iter)))
    isempty(iter) && return default
    res, rest = Iterators.peel(iter)
    fa = f(res)
    for x in rest
        fx = f(x)
        if fx > fa
            res = x
            fa = fx
        end
    end

    return res
end

compmax(CArray) = CArray[ (abs.(CArray) .== maximum(abs.(CArray))) .& (angle.(CArray) .== maximum( angle.(CArray))) ][1]


Main.isless(u1::ComplexF64, u2::ComplexF64) = abs2(u1) < abs2(u2)

function maxby5(arr)
    arr_max = arr[argmax(map(abs, arr))]
end

a = rand(ComplexF64,10)

using BenchmarkTools
@btime maxby0(abs,$a)
@btime maxby1(abs,$a)
@btime maxby2(abs,$a)
@btime compmax($a)
@btime maximum($a)
@btime maxby5($a)

Output for a vector of length 10:

>841.653 ns (1 allocation: 240 bytes)
>214.797 ns (0 allocations: 0 bytes)
>118.961 ns (0 allocations: 0 bytes)
>Execution fails
>20.340 ns (0 allocations: 0 bytes)
>144.444 ns (1 allocation: 160 bytes)

Output for a vector of length 1000:

>315.100 μs (1 allocation: 15.75 KiB)
>25.299 μs (0 allocations: 0 bytes)
>12.899 μs (0 allocations: 0 bytes)
>Execution fails
>1.520 μs (0 allocations: 0 bytes)
>14.199 μs (1 allocation: 7.94 KiB)

Output for a vector of length 1000 (with all comparisons made with abs2):

>35.399 μs (1 allocation: 15.75 KiB)
>3.075 μs (0 allocations: 0 bytes)
>1.460 μs (0 allocations: 0 bytes)
>Execution fails
>1.520 μs (0 allocations: 0 bytes)
>2.211 μs (1 allocation: 7.94 KiB)

Some remarks :

Sorting clearly (and as expected) slows the operations
Using abs2 saves a lot of performance (expected as well)

To conclude :

As a built-in function will provide this in 1.7, I will avoid using the additional Main.isless method, though it is all things considered the most performing one, to not modify the core of my julia
The maxby1 and maxby2 allocate nothing
The maxby1 feels more readable

the winner is therefore Andrej Oskin!

EDIT n°2 a new benchmark using the corrected compmax implementation

julia> @btime maxby0(abs2,$a)
  36.799 μs (1 allocation: 15.75 KiB)
julia> @btime maxby1(abs2,$a)
  3.062 μs (0 allocations: 0 bytes)
julia> @btime maxby2(abs2,$a)
  1.160 μs (0 allocations: 0 bytes)
julia> @btime compmax($a)
  26.899 μs (9 allocations: 12.77 KiB)
julia> @btime maximum($a)
  1.520 μs (0 allocations: 0 bytes)
julia> @btime maxby5(abs2,$a)
2.500 μs (4 allocations: 8.00 KiB)

Andrej Oskin Andrej Oskin · Accepted Answer · 2021-02-18T14:22:55

In Julia 1.7 you can use argmax

julia> a = rand(ComplexF64,4)
4-element Vector{ComplexF64}:
  0.3443509906876845 + 0.49984979589871426im
  0.1658370274750809 + 0.47815764287341156im
  0.4084798173736195 + 0.9688268736875587im
 0.47476987432458806 + 0.13651720575229853im

julia> argmax(abs2, a)
0.4084798173736195 + 0.9688268736875587im

Since it will take some time to get to 1.7, you can use the following analog

maxby(f, iter) = reduce(iter) do x, y
                   f(x) > f(y) ? x : y
                 end
julia> maxby(abs2, a)
0.4084798173736195 + 0.9688268736875587im

UPD: slightly more efficient way to find such maximum is to use something like

function maxby(f, iter; default = zero(eltype(iter)))
    isempty(iter) && return default
    res, rest = Iterators.peel(iter)
    fa = f(res)
    for x in rest
        fx = f(x)
        if fx > fa
            res = x
            fa = fx
        end
    end

    return res
end

Generic maximum/minimum function for complex numbers

4 Answers