Optimized unique and union! methods

[gajomi]

This PR implements performance improvements for unique and union! and adds an additional method signature for unique. The basic outline of the different performance issues can be found in #14649. There are four separates commits in the pull request which implement

1. optimized methods for associative collections
2. optimized methods for collections with "small" eltypes, as well as an additional method unique(c,k), which returns at most k unique elements from a collection c. This is used in the implementation of the aforementioned optimizations, but I think is a useful method in its own right and worth exposing to the user.
3. optimized unique method for range types
4. optimized unique and union methods for sparse-like arrays such as Diagonal,Bidiagonal,SparseVector etc.

Each of these optimizations have different regimes in which performance improvements become apparent. A benchmark that illustrates performance improvements for unique in action can be seen by looking at the diff in this gist: https://gist.github.com/gajomi/0c184160042ddfc329ab/revisions (results similar for union!). Generally speaking, one sees improvements in walltime performance, although there are some minor performance regressions in memory allocation for some methods.

[gajomi]

NOTE: The build generates warnings about ambiguous methods along the lines of

WARNING: New definition 
    union!(Base.Set, Base.LinAlg.Tridiagonal{#T<:Any}) at linalg/tridiag.jl:495
is ambiguous with: 
    union!(Base.Set{Int8}, AbstractArray{Int8, N<:Any}) at set.jl:64.
To fix, define 
    union!(Base.Set{Int8}, Base.LinAlg.Tridiagonal{Int8})
before the new definition.

On my system the right method seems to be called anyway

julia> A = Diagonal(Int8[1,2,3,4])
4x4 Diagonal{Int8}:
 1  ⋅  ⋅  ⋅
 ⋅  2  ⋅  ⋅
 ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  4

julia> @which union!(Set(),A)
union!{T}(s::Set{T<:Any}, D::Diagonal{T}) at linalg/diagonal.jl:240

but it would be preferable to get rid of these warnings. As far as I can tell the recommended approach would require duplicating the functionality for collections with small eltypes to all the structured array types. I would like to avoid this and was hoping to find suggestions to work around it.

[timholy]

From a brief skim, this looks amazingly thorough. I'm not sure I understand the reason for all the constructs like [v[1]; zero(T); v[2:end]], can you explain?

As far as the ambiguity warning goes: dispatching on both the element type and the container type is often a formula for trouble. I'm not sure which should take precedence, but let's suppose it's the container type. Then rather than "sharing an implicit common node" (which is what the ambiguity warning is suggesting you define), I think a better approach is to cascade to differently-named methods:

## Split on container type
function union!(s::Set, x::Diagonal)
    # code to handle Diagonal
end
function union!(s::Set, Tridiagonal)
    # code to handle Tridiagonal
end
union!(s::Set, x::AbstractArray) = _union!(s, x)  # fallback method

## Split on element type
function _union!{T<:SpecialElementTypes}(s::Set{T}, x::AbstactArray{T})
    ...
end
function _union!(s::Set, x::AbstractArray)
    # the generic method
    ...
end

(Stated otherwise, from a given node in the method "tree," split on either container type OR element type, not both simultaneously.)

[gajomi]

EDIT: I gave a very confusing answer and so thought it best to write a new one

The most straightforward way to unique methods for structured matrices like Diagonal, Bidiagonal, etc, would be to loop over over the elements of the matrix accumulating unique entries in column major order until one reached an entry that belonged to the structure zero indexes. At this point one would push the zero entry into the unique if it had not yet been seen. After this point one would continue iterating through the values of the matrix, but skipping over all indices corresponding to structured zeros, since this element is already known to be in the unique result.

The main problem with the above "straightforward" method is that is essentially requires reimplementing all the logic of the generic unique method across a bunch of specialized ones. Furthermore, if one wants to benefit from optimizations for small eltypes as described above this has to be reimplemented all over again.

So the idea behind the implementation in the pull request is to form a new matrix that only contains the data of the nonzero regions in the correct (column major) order plus one or more zero elements in appropriate places so that the result of calling generic unique on this new matrix is the same as calling generic unique on the full matrix. Gettting this matrix requires the funny looking constructs mentioned above.

This perhaps best explained with an example:

julia> A = SymTridiagonal([1,2,3,4,5],[6,7,8,1])
5x5 SymTridiagonal{Int64}:
 1  6  ⋅  ⋅  ⋅
 6  2  7  ⋅  ⋅
 ⋅  7  3  8  ⋅
 ⋅  ⋅  8  4  1
 ⋅  ⋅  ⋅  1  5

julia> B = hcat([A.dv[1],A.ev[1]],
                [zero(Int64),A.dv[2]],
                [A.ev[2:end] A.dv[3:end]]')
2x5 Array{Int64,2}:
 1  0  7  8  1
 6  2  3  4  5

julia> unique(B)'
1x9 Array{Int64,2}:
 1  6  0  2  7  3  8  4  5

This approach allows for zeros to be found before the structured zero portion of the matrix, but without the need to branch upon their detection. Furthermore, in the case that the matrix has small eltype, the intermediate matrix B will have this eltype so that the final call uses the unique method for small eltypes. In this way, one may possibly use both optimizations from the same top-level call.

[gajomi]

@timholy not sure I fully grok the cascading method approach, but I'll give it a try.

[tkelman]

generally these should always be sorted as part of the data structure convention

[gajomi]

Got it, this fact was documented for SparseMatrixCSC but I didn't see it for SparseVectors. I may elect to take the liberty of adding this to comments.

[tkelman]

yes, please do

[tkelman]

could use some inline comments explaining the approach and constructions used. do these work for the case of 1x1 matrices of these types?

[nalimilan]

Maybe this should return an Array instead (i.e. call collect)? The documentation for unique isn't very clear about what kind of "array" it returns. Not sure.

[hayd]

This was the main idea of #14649.

IMO unique needn't always be an array, if the datastructure is already unique e.g. Set/Range/etc. then a no-op should be preferred.

[nalimilan]

I'm fine with that. But note that sets are applied collect, as they are not AbstractArrays.

[hayd]

I don't think unique on Sets need collect, and return an AbstractArray, that was also in 14649. Or maybe I'm misunderstanding you.

[nalimilan]

All I'm saying is that this is what happens with this PR, and that this sounds justified, as Sets are not AbstractArrays and e.g. cannot be indexed.

[DNF2]

Isn't this a concern?

_unique(x::Range) = x
length(unique(linspace(0,0,10)))
1
length(_unique(linspace(0,0,10)))
10

I was under the impression that Ranges are not in general unique.

[gajomi]

This was brought up in #14649. I had thought that the above behavior behavior was supposed to be a bug but reading back on the comments it looks like its still ambiguous. This PR has lingered on a bit long, so I modified the method to accommodate either case. WRT to the larger issue of confusion about range with zero step sizes see #15218.

[timholy]

@gajomi, I edited my post above to hopefully make the intention clearer. Here's a little diagram that might help:

"Product-specialization" is the solution that's recommended in the ambiguity warning, but over time I've found that it's rarely a good idea to follow that advice. The "cascade" model is the one I usually find more productive.

[gajomi]

Thanks everyone for the extensive comments. I have mostly fixed things up and as a natural result of the process of refactoring have slightly expanded the situations in which the optimization apply, so I'll have to think a bit about upgrading benchmarks.

[gajomi]

OK, this took quite a bit longer than I expected (there were many forks in the road I'll to revisit later), but the latest commits should address all the issues raised above.

New performance benchmarks can be found here: https://gist.github.com/gajomi/617466a1502dbf920252/revisions?diff=split

Differences from first commits:

• The various union! specializations to container types versus eltypes follows the strategy suggested by @timholy.
• The new union! implementation distinguishes between diagonal dispatch on specialized eltypes in both set/array arguments and single dispatch on the eltype of the array, which extends the optimization to new cases.
• Where appropriate Set annotation are generalized to AbstractSet.
• The test coverage was expanded to cover union! methods.

[gajomi]

It seems that a test associated with IntSet failed on most architectures:

@test union!(IntSet(UInt8[]),[collect(0x00:0xff); 0x01]) == IntSet(collect(0x00:0xff))

This gave a deprecation warning about pushing zeros into IntSets on my system, but otherwise passed. I removed the test and the associated method in latest commits.

[tkelman]

can you add a few more comments on this?

[gajomi]

Yeah I added some comments (and ended up refactoring lower branch)/

[gajomi]

This seems to be good to go (modulo some potentially needed fixups stemming from nearly orthogonal issue discussed at JuliaLang/LinearAlgebra.jl#310).

[tkelman]

what is the purpose of having the extra step here?

[gajomi]

union! for Set types needs a size hint before calling rawunion!, whereas in IntSet one gets a size hint for every push! (as the underlying container size depends on largest value).

[Keno]

I'm gonna call this one stale. As a general comment, it seemed like this PR combined quite a number of different changes. Generally it's easier to get things in if they are done if very small pieces. New API surface in particular usually is quite slow because it's hard to change later. That said, if anybody ever wants to pick this up again, you're welcome to of course.