sparse vectors join the higher order function party #19690

Sacha0 · 2016-12-23T00:16:16Z

This pull request extends sparse map[!]/broadcast[!] to sparse vectors and sparse vector/matrix combinations.

Specifically, this pull request's first commit: (1) introduces a common interface to SparseVector and SparseMatrixCSC sufficient for the purposes of map[!]/broadcast[!]; (2) rewrites the existing sparse map[!]/broadcast[!] code against that interface; and (3) places the relevant code in the new SparseArrays.HigherOrderFns submodule which lives in base/sparse/higherorderfns.jl, loaded after definition of both SparseVector and SparseMatrixCSC.

This pull request's second commit: (1) condenses and systematizes existing tests for generic sparse map[!]/broadcast[!]; (2) extends those tests to sparse vectors and sparse vector/matrix combinations; and (3) collects those tests, as well older tests of sparse broadcast[!], into the new test file test/sparse/higherorderfns.jl corresponding to base/sparse/higherorderfns.jl.

Performance is more or less the same (previous generic sparse map[!]/broadcast[!] methods versus those in this pull request), faster in some cases where I fixed a type stability regression.

Best!

(Tangentially, this code (and likely an appreciable fraction of sparsevector.jl/sparsematrix.jl) could be further unified and simplified by revising SparseVector/SparseMatrixCSC with a common interface (cleaner than that in this pull request). If that idea sparks interest in this thread, I'll write a mini-julep at some point.)

stevengj · 2016-12-23T02:29:38Z

@nanosoldier runbenchmarks(ALL, vs = ":master")

nanosoldier · 2016-12-23T05:44:24Z

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

tkelman · 2016-12-23T06:33:31Z

test/sparse/higherorderfns.jl

+    @test all(A[1,:] .* 3 .== AF[1,:] .* 3)
+    #@test A[:,1] .* 3 == AF[:,1] .* 3
+    @test all(A[:,1] .* 3 .== AF[:,1] .* 3)
+    #TODO: simple comparation with == returns false because the left side is a (two-dimensional) SparseMatrixCSC


this todo is really old, predates SparseVector being in base

(Is this a request to clear this TODO in this pull request? Or a request to clear it later? Thanks!)

I was about to do this in a local branch, but it would cause a conflict here

Happy to sort that however you think best --- let me know. Thanks!

it's just uncommenting, and deleting the all and todo lines

Done. Thanks!

stevengj · 2016-12-23T12:30:18Z

Benchmarks look good; the floatexp benchmarks seem noisy as usual.

stevengj · 2016-12-23T13:07:32Z

Why do we still have broadcast(::typeof(+), ...) methods for SparseVector?

stevengj · 2016-12-23T13:08:43Z

Do we still need the methods specialized on a pair of sparse matrices?

stevengj · 2016-12-23T13:10:03Z

It still doesn't handle mixtures of sparse vectors/matrices and scalars, correct?

Sacha0 · 2016-12-23T21:51:04Z

Why do we still have broadcast(::typeof(+), ...) methods for SparseVector?

base/sparse/sparsevectors.jl contains a few not-exported map-like methods and several children [1] (e.g. the methods you mention). This pull request should enable removal of some of that code without regressions and more with potentially negligible regressions. That cleanup seems a substantial project in itself (with potentially controversial points, and largely orthogonal to this work in practice). So in the interest of keeping this pull request confined / uncontroversial / readily mergeable, I left that work separate.

[1] ... mostly defined against AbstractSparseVector, though most of those methods are SparseVector-specific?

Do we still need the methods specialized on a pair of sparse matrices?

Yes: The generic methods still appreciably lag the specialized methods in some cases. (A few performance improvements I have in mind might close the gap enough to nix some of the specialized methods. Post feature freeze work though.)

It still doesn't handle mixtures of sparse vectors/matrices and scalars, correct?

Correct: #19641 / merging #19667 blocks my solution for handling combinations including broadcast scalars. Best!

tkelman · 2016-12-24T20:49:22Z

test/sparse/higherorderfns.jl

+            @test broadcast(+, X, Y) == sparse(broadcast(+, fX, fY))
+            @test broadcast(*, X, Y) == sparse(broadcast(*, fX, fY))
+            @test broadcast(f, X, Y) == sparse(broadcast(f, fX, fY))
+            try Base.Broadcast.broadcast_indices(spzeros((shapeX .- 1)...), Y)


split this over multiple lines. if the intention is that the try should always throw, then add a @test false after so it doesn't silently pass if it happens not to throw

The intention is for the @test_throws in the catch body to only evaluate when the statement in the try body throws. Is there a better way to express this? Thanks!

(In other words, the test logic is: Check whether the arguments are shape-broadcast-incompatible. If so, make certain the method being @test_throws'd throws.)

Is there another way of determining shape compatibility than checking for something throwing? Calling some internal boolean function maybe? Otherwise seems like you might not be checking carefully enough which inputs do or do not trigger the test

Is there another way of determining shape compatibility than checking for something throwing? Calling some internal boolean function maybe?

To the best of my knowledge none such exists. (~~I should instead call check_broadcast_indices for clarity, but that similarly throws or returns nothing.~~ On second thought no, broadcast_indices is correct as check_broadcast_indices requires the broadcast shape (which you retrieve from broadcast_indices) as a argument.) Thoughts? Thanks!

(Split over multiple lines. Thanks!)

hm maybe we should make that check function return a boolean instead of nothing-or-throw?

Seems @timholy designed that mechanism --- perhaps ping him? (Orthogonal to this pull request in any case?) Thanks!

tkelman · 2016-12-24T20:49:54Z

test/sparse/sparse.jl

@@ -1674,168 +1596,9 @@ end
 # 19304
 @inferred hcat(sparse(rand(2,1)), eye(2,2))

-# Test that broadcast[!](f, [C::SparseMatrixCSC], A::SparseMatrixCSC, B::SparseMatrixCSC)
-# returns the correct (densely populated) result when f(zero(eltype(A)), zero(eltype(B))) != 0


which set of the new tests replaces this?

The broadcast[!] implementation specialized for pairs of (input) sparse vectors/matrices testset replaces those tests.

tkelman · 2016-12-24T20:54:25Z

test/sparse/sparse.jl

-    X = sparse(broadcast(f, Array(A), Array(B), Array(C)))
-    broadcast!(f, copy(X), A, B, C) # warmup for @allocated
-    @test_broken (@allocated broadcast!(f, X, A, B, C)) == 0
-    # this last test allocates 16 bytes in the entry point for broadcast!, but none of the


~~figure out what was happening here?~~ nevermind

Sacha0 · 2016-12-26T18:34:24Z

Absent objections or requests for time, I plan to merge this Wednesday morning PDT. Best!

stevengj · 2016-12-26T19:13:09Z

You can now remove the workaround in this test, yes, and just do @test exact_equal(x .* x, abs2(x)), or abs2.(x)?

Sacha0 · 2016-12-26T21:18:52Z

You can now remove the workaround in this test, yes, and just do @test exact_equal(x .* x, abs2(x)), or abs2.(x)?

Happily yes, and done. (The abs2.(x) change is in #18564.) Thanks!

(If you or @tkelman would be comfortable merging this PR prior to Wednesday, please feel welcome.)

tkelman · 2016-12-26T21:27:22Z

I have not finished reviewing this yet, but would like to do so today or tomorrow.

edit: got through it all, locally undoing some of the rearrangements as https://github.com/JuliaLang/julia/compare/d8a57182f49787055b3a864cc86401988be4bcdf...tkelman:tk/sparsevechof?w=1 just for the purposes of reviewing made it a little easier to compare

tkelman · 2016-12-28T01:03:23Z

test/sparse/higherorderfns.jl

+            @test broadcast(*, X, Y, Z) == sparse(broadcast(*, fX, fY, fZ))
+            @test broadcast(f, X, Y, Z) == sparse(broadcast(f, fX, fY, fZ))
+            try Base.Broadcast.broadcast_indices(spzeros((shapeX .- 1)...), Y, Z)
+            catch @test_throws DimensionMismatch broadcast(+, spzeros((shapeX .- 1)...), Y, Z) end


this try-catch (and the one a few lines lower) would also read more clearly if on multiple lines

though I still think it would be preferable to have this serve as an independent verification of how broadcast_indices is expected to behave for the differently-sized inputs, rather than relying on it throwing - if it were to be refactored to no longer throw, then this test would silently stop doing anything

this try-catch (and the one a few lines lower) would also read more clearly if on multiple lines

Thanks for catching these other two instances! Split into multiple lines on push.

though I still think it would be preferable to have this serve as an independent verification of how broadcast_indices is expected to behave for the differently-sized inputs, rather than relying on it throwing - if it were to be refactored to no longer throw, then this test would silently stop doing anything

Agreed. Not being certain how best to achieve that strengthening at this time, can we revisit that point later? Thanks!

leave a todo comment

tkelman · 2016-12-28T01:17:46Z

test/sparse/higherorderfns.jl

+            # entry point for broadcast! with --track-allocation=all, but that first line
+            # almost certainly should not allocate. so not certain what's going on.
+            @test broadcast!(f, Q, X, Y, Z) == sparse(broadcast!(f, fQ, fX, fY, fZ))
+            # --> test shape checks for both braodcast and broadcast! entry points


"braodcast" typo

Fixed on push. Thanks!

tkelman · 2016-12-28T01:21:05Z

test/sparse/sparse.jl

-        # separate horizontal and vertical expansion
-        @test broadcast(op, A, B, C) == sparse(broadcast(op, fA, fB, fC))
-        @test broadcast!(op, X, A, B, C) == sparse(broadcast!(op, fX, fA, fB, fC))
-        # simultaneous horizontal and vertical expansion


does this require 4 inputs or is it covered by the new loop with only 3?

Three input arguments suffice to test the relevant code paths, and indeed the testset titled "broadcast[!] implementation capable of handling >2 (input) sparse vectors/matrices" covers the relevant code. Thanks!

tkelman · 2016-12-28T10:47:50Z

base/sparse/higherorderfns.jl

+# (5) Define _map_[not]zeropres! specialized for a pair of (input) sparse vectors/matrices.
+# (6) Define general _map_[not]zeropres! capable of handling >2 (input) sparse vectors/matrices.
+# (7) Define _broadcast_[not]zeropres! specialized for a pair of (input) sparse vectors/matrices.
+# (8) Define general _broadcast_[not]zeropres! capabel of handling >2 (input) sparse vectors/matrices.


"capabel" typo

Good catch! Fixed on push. Thanks!

tkelman · 2016-12-28T11:06:22Z

base/sparse/higherorderfns.jl

+end
+# helper functions for these methods and some of those below
+@inline _densecoloffsets(A::SparseVector) = 0
+@inline _densecoloffsets(A::SparseMatrixCSC) = 0:A.m:(A.m*A.n - 1)


why isn't the end point (A.m * (A.n-1)) here?

Facepalm-worthy unintentional code obfuscation is why. Your suggestion is much better. Fixed on push. Thanks!

A.m * (0:A.n-1) might be even better, not sure

tkelman · 2016-12-28T11:27:26Z

base/sparse/sparsematrix.jl

-    @inbounds for j in 1:C.n
-        C.colptr[j] = Ck
-        ks = _startindforbccol_all(j, expandshorzs, As)
-        stopks = _stopindforbccol_all(j + 1, expandshorzs, As)


was the j+1 here a mistake or did you change the definition of this helper function?

Changed the helper function's definition IIRC, but my memory is vague on that point. Thanks!

…e vector/matrix combinations. Extend generic sparse map[!]/broadcast[!] to sparse vectors and sparse vector/matrix combinations. Do so by introducing a common interface to SparseVector and SparseMatrixCSC for the purposes of map[!]/broadcast[!], and rewriting sparse map[!]/broadcast[!] against that interface. Relocate that code to a separate file/module base/sparse/higherorderfns.jl/SparseArrays.HigherOrderFns, loaded after definition of both SparseVector and SparseMatrixCSC.

Condense and systematize existing tests for generic sparse map[!]/broadcast[!], and extend to sparse vectors and vector/matrix combinations. Relocate new test code to a separate file test/sparse/higherorderfns.jl corresponding to base/sparse/higherorderfns.jl. (Test/sparse/sparsevector.jl is hypothetically confined to SparseVectors, and test/sparse/sparse.jl mostly dedicated to SparseMatrixCSCs.) Move older tests of sparse broadcast[!] into that new file as well.

Sacha0 · 2016-12-30T20:18:04Z

Absent objections I plan to merge this once CI approves. Best!

…ned [ci skip].

Sacha0 · 2016-12-31T03:44:25Z

(AV i686 failure unrelated.)

Sacha0 mentioned this pull request Dec 23, 2016

Fix #19561 (sparse map/broadcast where the output eltype is not a concrete subtype of Number) #19589

Merged

Sacha0 added sparse Sparse arrays broadcast Applying a function over a collection labels Dec 23, 2016

Sacha0 added this to the 0.6.0 milestone Dec 23, 2016

tkelman reviewed Dec 23, 2016

View reviewed changes

Sacha0 force-pushed the sparsevechof branch from 0c539a9 to d0b450d Compare December 24, 2016 00:01

tkelman reviewed Dec 24, 2016

View reviewed changes

Sacha0 force-pushed the sparsevechof branch from d0b450d to 2f07a81 Compare December 24, 2016 21:24

This was referenced Dec 25, 2016

Deprecate vectorized two-argument complex in favor of compact broadcast syntax #19712

Merged

Deprecate vectorized min and max over pairs of sparse matrices #19713

Closed

Sacha0 force-pushed the sparsevechof branch from 2f07a81 to 15f625d Compare December 26, 2016 06:44

Sacha0 force-pushed the sparsevechof branch from 15f625d to 7bfd0b9 Compare December 26, 2016 21:16

Sacha0 mentioned this pull request Dec 27, 2016

broadcast[!] over combinations of scalars and sparse vectors/matrices #19724

Merged

Sacha0 force-pushed the sparsevechof branch from 7bfd0b9 to 89ea1a8 Compare December 27, 2016 22:25

tkelman reviewed Dec 28, 2016

View reviewed changes

Sacha0 added 2 commits December 30, 2016 12:03

Sacha0 force-pushed the sparsevechof branch from 89ea1a8 to 913f637 Compare December 30, 2016 20:17

Mark tests for sparse broadcast shape checks that should be strengthe…

2786db2

…ned [ci skip].

Sacha0 merged commit 1ed4545 into JuliaLang:master Dec 31, 2016

Sacha0 deleted the sparsevechof branch December 31, 2016 03:43

tkelman mentioned this pull request Dec 31, 2016

map lacks specializations for sparse matrices/vectors #19363

Closed

Sacha0 mentioned this pull request Dec 31, 2016

Deprecate manually vectorized methods in favor of dot syntax, monolithic edition #19791

Merged

sparse vectors join the higher order function party #19690

sparse vectors join the higher order function party #19690

Conversation

Sacha0 commented Dec 23, 2016 • edited Loading

stevengj commented Dec 23, 2016

nanosoldier commented Dec 23, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkelman Dec 23, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj commented Dec 23, 2016

stevengj commented Dec 23, 2016

stevengj commented Dec 23, 2016

stevengj commented Dec 23, 2016

Sacha0 commented Dec 23, 2016

Choose a reason for hiding this comment

Sacha0 Dec 24, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Sacha0 Dec 24, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkelman Dec 24, 2016 • edited Loading

Choose a reason for hiding this comment

Sacha0 commented Dec 26, 2016

stevengj commented Dec 26, 2016 • edited Loading

Sacha0 commented Dec 26, 2016 • edited by yuyichao Loading

tkelman commented Dec 26, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Sacha0 commented Dec 30, 2016

Sacha0 commented Dec 31, 2016

Sacha0 commented Dec 23, 2016 •

edited

Loading

tkelman Dec 23, 2016 •

edited

Loading

Sacha0 Dec 24, 2016 •

edited

Loading

Sacha0 Dec 24, 2016 •

edited

Loading

tkelman Dec 24, 2016 •

edited

Loading

stevengj commented Dec 26, 2016 •

edited

Loading

Sacha0 commented Dec 26, 2016 •

edited by yuyichao

Loading

tkelman commented Dec 26, 2016 •

edited

Loading