De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

Sacha0 · 2016-07-04T01:09:43Z

This PR's first commit revises the implementation of vectorized unary functions over SparseMatrixCSCs, leveraging higher order functions to displace a set of macros and evals. Regarding performance, see below.

This PR's second commit then transitions those vectorized unary functions over SparseMatrixCSCs to compact broadcast syntax, accordingly revises the associated tests, expands those tests, and then adds deprecations for the vectorized calls.

Question: Both the existing code and that in this PR separate unary functions into three classes (with disjoint logic): (1) unary functions that map zeros to zeros and may map nonzeros to zeros; (2) unary functions that map zeros to zeros and nonzeros to nonzeros; and (3) unary functions that map both zeros and nonzeros to nonzeros. In the first class, when a nonzero in the original sparse matrix miraculously maps to exact zero under, e.g., sin, the implementation expunges that entry. This behavior seems like a holdover from the aggressive-stored-zero-expunging era, now inconsistent with behavior elsewhere. Perhaps classes one and two should be merged? Merging those two classes would enable some code reduction and might increase the performance of operations presently falling under class one (by eliminating a branch in an inner loop, and enabling @simd decoration of that inner loop). Thoughts?

Perf
With the exception of real and imag, for functions in the first class the existing and new implementations perform similarly or slightly favor the new implementation. The new implementation seems to resolve a type instability or so with real and imag, significantly improving performance: Where real is the existing implementation and ho_real this PR's, this bench

using BenchmarkTools

smallN = 10^1; smallA = sprand(smallN, smallN, 5/smallN); smallC = smallA + im*smallA;
largeN = 10^6; largeA = sprand(largeN, largeN, 10/largeN); largeC = largeA + im*largeA;

presmat = smallC;
benchpres = @benchmarkable real(Float64, presmat); tune!(benchpres);
benchnew = @benchmarkable ho_real(Float64, presmat); tune!(benchnew);
medianpres = median(run(benchpres));
mediannew = median(run(benchnew));
println(ratio(mediannew, medianpres))

for presmat = smallC yields

BenchmarkTools.TrialRatio:
  time:             0.07120253164556962
  gctime:           1.0
  memory:           0.5189873417721519
  allocs:           0.06557377049180328
  time tolerance:   5.00%
  memory tolerance: 1.00%

and for presmat = largeC yields

BenchmarkTools.TrialRatio:
  time:             0.16305659839888223
  gctime:           0.16709694727350574
  memory:           0.3442608362244807
  allocs:           3.4976803383995736e-7
  time tolerance:   5.00%
  memory tolerance: 1.00%

For functions in the second class, this PR's implementation exhibits significantly better performance on small matrices and somewhat better on large matrices, e.g. for the bench above with presmat = smallA, comparing expm1 and ho_expm1 yields

BenchmarkTools.TrialRatio:
  time:             0.18331303288672351
  gctime:           1.0
  memory:           0.7735849056603774
  allocs:           0.4
  time tolerance:   5.00%
  memory tolerance: 1.00%

and for presmat = largeA

BenchmarkTools.TrialRatio:
  time:             0.8084212063475197
  gctime:           0.9765840983181242
  memory:           0.9999974302672074
  allocs:           0.4375
  time tolerance:   5.00%
  memory tolerance: 1.00%

This PR's implementation also appears to resolve a type instability or so with abs and abs2, e.g. comparing abs2 and ho_abs2 for presmat = smallC the above bench yields

BenchmarkTools.TrialRatio:
  time:             0.05240027045300879
  gctime:           1.0
  memory:           0.5189873417721519
  allocs:           0.06557377049180328
  time tolerance:   5.00%
  memory tolerance: 1.00%

and for presmat = largeC

BenchmarkTools.TrialRatio:
  time:             0.11257829481710147
  gctime:           0.1572109282310488
  memory:           0.3442608362244807
  allocs:           3.4976803383995736e-7
  time tolerance:   5.00%
  memory tolerance: 1.00%

For functions in the third class, performance differences are within the noise. Best!

tkelman · 2016-07-04T01:51:24Z

base/deprecated.jl

+eval(Multimedia, quote
+export @MIME
+macro MIME(s)
+    Base.warn_once("@MIME(\"\") is deprecated, use MIME\"\" instead.")


wrong conflict resolution here

Good catch. Fixed? Thanks!

Sacha0 · 2016-07-05T18:19:37Z

Messed things up a bit while addressing comments. Should be in order again now. Thanks!

tkelman · 2016-07-18T11:07:31Z

needs rebase

Sacha0 · 2016-07-19T21:00:11Z

needs rebase

Thanks for the ping! Rebased.

tkelman · 2016-08-03T03:58:33Z

base/deprecated.jl

@@ -793,6 +793,19 @@ function transpose(x)
    return x
 end

+# Deprecate vectorized unary functions over sparse matrices in favor of compact broadcast syntax (#17265).


we won't be deprecating these in 0.5, move after the "to be deprecated in 0.6"

Sacha0 · 2016-08-05T18:35:03Z

Rebased for 0.6. Thanks!

tkelman · 2016-08-05T18:58:04Z

base/sparse/sparsematrix.jl

+# Operations that map both zeros and nonzeros to zeros, yielding a dense matrix
+"""
+Takes unary function `f` that maps both zeros and nonzeros to nonzeros, and returns a new
+`Matrix{TvB}` `B` effectively generated by applying `f` to every entry in `A`.


why "effectively" ?

Due to the B = fill(f(zero(Tv)), size(A)), this method calls f only once rather than once for each zero in A. Would make a difference if e.g. f has side effects.

worth writing that down in the docstring?

Good call. Clarified in docstring. Thanks!

tkelman · 2016-08-06T07:57:24Z

ready to merge?

ViralBShah · 2016-08-06T10:35:27Z

Looks like a good candidate for backporting as well for 0.5.x.

tkelman · 2016-08-06T11:14:32Z

No it isn't, it's a whole bunch of deprecations.

Sacha0 · 2016-08-07T18:15:38Z

ready to merge?

I assume this query was not directed at me, but in case it was: Yes, this PR should be ready to merge, though it might interact with #17302 and this PR should be easier to rebase than #17302. Thanks and best!

Sacha0 · 2016-08-10T16:43:53Z

Rebased resolving deprecation conflict. Best!

tkelman · 2016-08-25T16:45:22Z

@Sacha0 sorry we frequently do a bad job of merging your (always impeccably well-done) PR's before they hit conflicts. Rebase?

Sacha0 · 2016-08-25T21:11:34Z

No worries! Your time is much better spent elsewhere at the moment. Also thanks for the kind words! Related question: Though master is technically 0.6-dev at this point, I imagine that prior to the actual 0.5 release it's advantageous to keep master as close as possible to 0.5.x? That is, I imagine holding noncritical PRs targeting 0.6 till 0.5 is out the door makes life easier? If so, apologies for bothering you with noncritical PRs targeting 0.6 recently, and I'll hold bumps / non-bugfix work till 0.5 is out. Thanks and best!

tkelman · 2016-08-25T21:16:27Z

It's fine actually, as long as things aren't a huge risk of causing non-trivial conflicts with bugfixes that need backporting. 0.5 is most likely only one more RC away from final, and we need to get moving on 0.6 just as much as we need to get 0.5 final done.

…higher order functions and multiple dispatch to displace eval. Fixes some apparent type instabilities.

…act broadcast syntax, accordingly revise and expand the associated tests, and add deprecations for the vectorized syntax.

Sacha0 · 2016-08-27T16:24:44Z

Thanks for reviewing / merging!

stevengj · 2016-08-31T14:39:31Z

I wonder if we should have a more general function that does broadcast on sparse matrices:

function broadcast{T<:Number}(f::Function, A::AbstractSparseArray{T})
    fzero = f(zero(T))
    if fzero == zero(fzero)
       ... broadcast to sparse array, operating f only on nonzero elements ...
    else
       .... broadcast to dense array, though as an optimization we could just use fzero for the zero elements ... or throw an error?
    end
end

This assumes f is pure, though.

tkelman · 2016-08-31T14:55:19Z

https://github.com/JuliaLang/julia/issues/7010

Sacha0 · 2016-08-31T16:58:45Z

Would the sketch above be type unstable? Would a version parametric on the function type recover type stability? Best!

stevengj · 2016-08-31T17:02:22Z

@Sacha0, the problem is that it assumes f is pure; I think it should be possible to make it type stable.

I guess we could use traits for this, but they wouldn't help much for anonymous functions that are constructed during loop fusion.

tkelman reviewed Jul 4, 2016
View reviewed changes

Sacha0 force-pushed the unarybcastspmat branch from b3b8350 to fafacd9 Compare July 4, 2016 02:20

ViralBShah added the sparse Sparse arrays label Jul 4, 2016

Sacha0 force-pushed the unarybcastspmat branch from fafacd9 to 57ea2f9 Compare July 5, 2016 18:18

Sacha0 mentioned this pull request Jul 6, 2016

Deprecate functions vectorized via @vectorize_(1|2)arg in favor of compact broadcast syntax #17302

Merged

4 tasks

Sacha0 force-pushed the unarybcastspmat branch from 57ea2f9 to 08eb7b6 Compare July 19, 2016 20:59

tkelman reviewed Aug 3, 2016
View reviewed changes

Sacha0 force-pushed the unarybcastspmat branch from 08eb7b6 to e2042d7 Compare August 5, 2016 18:25

tkelman reviewed Aug 5, 2016
View reviewed changes

Sacha0 force-pushed the unarybcastspmat branch from e2042d7 to d380f1a Compare August 5, 2016 20:59

Sacha0 force-pushed the unarybcastspmat branch from d380f1a to 57b032b Compare August 10, 2016 16:43

tkelman added this to the 0.6.0 milestone Aug 25, 2016

Sacha0 added 2 commits August 25, 2016 14:25

Rewrite vectorized unary functions over SparseMatrixCSCs, leveraging …

4ea088a

…higher order functions and multiple dispatch to displace eval. Fixes some apparent type instabilities.

Transition vectorized unary functions over SparseMatrixCSCs to comp…

bd7da66

…act broadcast syntax, accordingly revise and expand the associated tests, and add deprecations for the vectorized syntax.

Sacha0 force-pushed the unarybcastspmat branch from 57b032b to bd7da66 Compare August 25, 2016 21:27

tkelman merged commit 379c248 into JuliaLang:master Aug 27, 2016

Sacha0 deleted the unarybcastspmat branch August 27, 2016 16:24

Sacha0 mentioned this pull request Aug 31, 2016

Classes of unary operations in broadcast over sparse arrays #18309

Closed

Sacha0 mentioned this pull request Sep 20, 2016

Deprecate vectorized round methods in favor of compact broadcast syntax #18590

Closed

This was referenced Dec 4, 2016

huge performance hit in 0.6: negative of a sparse matrix #19503

Closed

Fix #19503 (provide a unary minus method specialized for sparse matrices) #19530

Merged

stevengj mentioned this pull request May 26, 2017

real(::SparseMatrixCSC) should not be deprecated #22065

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

Sacha0 commented Jul 4, 2016 •

edited

Loading

tkelman Jul 4, 2016

Sacha0 Jul 4, 2016

Sacha0 commented Jul 5, 2016

tkelman commented Jul 18, 2016

Sacha0 commented Jul 19, 2016

tkelman Aug 3, 2016

Sacha0 commented Aug 5, 2016

tkelman Aug 5, 2016

Sacha0 Aug 5, 2016

tkelman Aug 5, 2016

Sacha0 Aug 5, 2016

tkelman commented Aug 6, 2016

ViralBShah commented Aug 6, 2016

tkelman commented Aug 6, 2016

Sacha0 commented Aug 7, 2016

Sacha0 commented Aug 10, 2016

tkelman commented Aug 25, 2016

Sacha0 commented Aug 25, 2016

tkelman commented Aug 25, 2016

Sacha0 commented Aug 27, 2016

stevengj commented Aug 31, 2016 •

edited

Loading

tkelman commented Aug 31, 2016

Sacha0 commented Aug 31, 2016

stevengj commented Aug 31, 2016 •

edited

Loading

De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

De-eval-ify vectorized unary functions over SparseMatrixCSCs, and then transition them to compact broadcast syntax #17265

Conversation

Sacha0 commented Jul 4, 2016 • edited Loading

tkelman Jul 4, 2016

Choose a reason for hiding this comment

Sacha0 Jul 4, 2016

Choose a reason for hiding this comment

Sacha0 commented Jul 5, 2016

tkelman commented Jul 18, 2016

Sacha0 commented Jul 19, 2016

tkelman Aug 3, 2016

Choose a reason for hiding this comment

Sacha0 commented Aug 5, 2016

tkelman Aug 5, 2016

Choose a reason for hiding this comment

Sacha0 Aug 5, 2016

Choose a reason for hiding this comment

tkelman Aug 5, 2016

Choose a reason for hiding this comment

Sacha0 Aug 5, 2016

Choose a reason for hiding this comment

tkelman commented Aug 6, 2016

ViralBShah commented Aug 6, 2016

tkelman commented Aug 6, 2016

Sacha0 commented Aug 7, 2016

Sacha0 commented Aug 10, 2016

tkelman commented Aug 25, 2016

Sacha0 commented Aug 25, 2016

tkelman commented Aug 25, 2016

Sacha0 commented Aug 27, 2016

stevengj commented Aug 31, 2016 • edited Loading

tkelman commented Aug 31, 2016

Sacha0 commented Aug 31, 2016

stevengj commented Aug 31, 2016 • edited Loading

Sacha0 commented Jul 4, 2016 •

edited

Loading

stevengj commented Aug 31, 2016 •

edited

Loading

stevengj commented Aug 31, 2016 •

edited

Loading