MIT-licensed sparse() parent method and expert driver #14798

Sacha0 · 2016-01-26T03:16:52Z

Followup to #13001, #14631, and #14702. This pull request replaces the LGPL-licensed sparse parent method with an MIT-licensed version and introduces an expert driver sparse! underlying that method (see documentation for details on the expert driver).

In accord with discussion in #12605, #9928, #9906, and #6769, this PR's sparse method does not automatically purge numerical zeros. Hence this pull request comments two tests of numerical-zero purging. Rather than simply being commented, perhaps those tests should remain with dropzeros!(sparse(...)) in place of sparse(...) and a reference to this PR?

Hats off to those who contributed to the existing implementation! Matching its performance in all tested cases was challenging. Benchmarks results, benchmark code; the results should be comprehensible without reading the code. tl;dr: In these benchmarks, this PR's implementation ranges from matching the existing implementation's performance to reducing runtime by ~30%, and seems to allocate approximately Vector{Ti}(m) + Vector{Ti}(n) less storage. Notably the sparser the matrix, the better this PR's implementation's relative performance. Though solving #13400 requires a fundamentally different algorithm, #13400's example serves as an extreme illustration of the preceding point:

using Benchmarks
@benchmark sparse([1;100000000], [1;1], [1;1], 100000000, 1, Base.AddFun())
@benchmark mitsparse([1;100000000], [1;1], [1;1], 100000000, 1, Base.AddFun())

(with mitsparse this PR's implementation) yields

================ Benchmark Results ========================
     Time per evaluation: 1.94 s [1.90 s, 1.98 s]
Proportion of time in GC: 9.71% [9.01%, 10.40%]
        Memory allocated: 2.24 gb
   Number of allocations: 13 allocations
       Number of samples: 4
   Number of evaluations: 4
 Time spent benchmarking: 9.67 s

================ Benchmark Results ========================
     Time per evaluation: 735.76 ms [723.62 ms, 747.89 ms]
Proportion of time in GC: 6.50% [5.98%, 7.03%]
        Memory allocated: 762.94 mb
   Number of allocations: 11 allocations
       Number of samples: 12
   Number of evaluations: 12
 Time spent benchmarking: 9.73 s

When editing base/sparse/csparse.jl to remove the existing methods, I folded all blocks to avoid peaking at the LGPL code. So someone should check that the removal was graceful given it was blind but for method signatures.

Apologies for responding slowly elsewhere while working on this!

hayd · 2016-01-26T04:31:46Z

This is great!

Travis shows a whitespace failure.

git rebase --whitespace=fix HEAD~1

Sacha0 · 2016-01-26T04:36:08Z

Travis shows a whitespace failure.

Derp! Fixed, thanks @hayd!

KristofferC · 2016-01-26T05:23:42Z

Very impressive work!

nalimilan · 2016-01-26T09:57:29Z

base/sparse/sparsematrix.jl

+ sparse(I, J, V, [m, n, combine])
+
+Create a sparse matrix `S` of dimensions `m` x `n` such that `S[I[k], J[k]] = V[k]`. The
+ `combine` function is used to combine duplicates. If `m` and `n` are not specified, they


I think the common style does not add indentation here.

Fixed, thanks!

StefanKarpinski · 2016-01-26T15:10:08Z

This is excellent work. Let's review the technical aspects of this and not focus too much on cosmetic issues, which can be fixed after merging this.

nalimilan · 2016-01-26T15:33:30Z

This is excellent work. Let's review the technical aspects of this and not focus too much on cosmetic issues, which can be fixed after merging this.

Well, these are clearly minor, but they are probably easier to correct while fixing more significant aspects than after merging (when nobody will care enough to make another PR).

Great work, BTW, though I'm not in position of evaluating it!

StefanKarpinski · 2016-01-26T15:37:46Z

I just don't want to nitpick a great PR to death with cosmetic stuff, but yes, it can be easier to fix before merging. @Sacha0, if it's easy for you to fix the cosmetic stuff, great. If not, then we can just merge and fix after.

ViralBShah · 2016-01-26T17:51:48Z

Nicely done!

tkelman · 2016-01-26T19:51:01Z

There are test failures that obviously need fixing, so "just merge and fix cosmetic stuff later" is a bit premature.

StefanKarpinski · 2016-01-26T20:10:06Z

Yes, obviously the test failures need to get fixed before merging, but let's focus on those rather than a slew of indentation and formatting issues.

Sacha0 · 2016-01-26T21:55:49Z

I greatly appreciate the detailed review / feedback, cosmetic included :). Constructive criticism is a gift; the more and sharper, the more I learn, so much thanks @nalimilan! I will do my best to address the comments.

CI is brilliant. I completely missed the arnoldi failure when testing locally! Thanks!

Sacha0 · 2016-01-26T22:29:57Z

I can reproduce the failure locally, but curiously Base.runtests("linalg/arnoldi") succeeds where Base.runtest("linalg") fails. Somehow the former call does not force bounds checking whereas the latter does?

Sacha0 · 2016-01-26T23:16:59Z

The linalg/arnoldi test failure should be fixed. I still wonder about the difference between Base.runtests("linalg/arnoldi") and Base.runtests("linalg") mentioned above. Thanks!

tkelman · 2016-01-27T01:29:53Z

If you're using Base.runtests, then I don't think anything unusual happens with bounds checking. More likely, running different sets of tests in the same process results in different random numbers.

Sacha0 · 2016-01-27T01:55:14Z

If you're using Base.runtests, then I don't think anything unusual happens with bounds checking. More likely, running different sets of tests in the same process results in different random numbers.

The issue was an out-of-bounds access. Running with --check-bounds=yes enabled reproduction under all conditions (at the REPL, via Base.runtests("linalg/arnoldi"), and via Base.runtests("linalg")). Without --check-bounds=yes, the issue appeared only during Base.runtests("linalg"). So it seems bounds checking state must have been involved somehow? Moreover, there are no random numbers involved in the offending line? Thanks!

Edit: For easy reference, the offending line:

A = sparse([1, 1, 2, 3, 4], [2, 1, 1, 3, 1], [2.0, -1.0, 6.1, 7.0, 1.5])

tkelman · 2016-01-27T02:08:32Z

Base.runtests("linalg") runs multiple processes, so it may be that the check-bounds flag only gets propagated to child processes when running multiple tests.

tkelman · 2016-01-27T02:12:05Z

base/sparse/sparsematrix.jl

+ "$(length(V)))")))
+ end
+
+ if m == 0 || n == 0 || coolen == 0


this should throw a ~~BoundsError~~ ArgumentError if elements of I or J are outside of m-by-n

edit: my bad, sorry

Fixed in lines 421 and 440 (ArgumentError -> BoundsError), thanks!

No, if m == 0 || n == 0 this should not return successfully if the input indices are out of bounds (aka if there are any of them)

Ah, I see, thanks! Fixing now.

Fixed (I think), thanks!

I suppose the first check (!isempty(I)) suffices given the earlier check length(I) == length(J) == length(V). But perhaps this is more clear. Thoughts?

Sure, that works. I think the error message can be the same as the non-empty case though, ~~row values I[k] must satisfy 1 <= I[k] <= m~~ row indices I[k] must satisfy 1 <= I[k] <= m etc. Also always good to add more tests for this kind of corner case.

After a little thought I might advocate for the distinct, explicit error message in now: The I[k]-specific error message may be confusing where, for example, m == 2, n == 0, and (I,J,V) = (2, 1, 1.0).

The current version that starts "where ..., any entry is necessarily out-of-bounds," reads to me as too meandering for an error message. I'd do something like this

if coolen != 0 if n == 0 throw(BoundsError("column indices J[k] must satisfy 1 <= J[k] <= n")) elseif m == 0 throw(BoundsError("row indices I[k] must satisfy 1 <= I[k] <= m")) end end

edit: except it should be ArgumentError, whoops

Beautiful. Copied verbatim. Thanks!

mbauman · 2016-01-27T02:38:50Z

base/sparse/sparsematrix.jl

+ @inbounds for k in 1:coolen
+ Ik, Jk = I[k], J[k]
+ if 1 > Jk || n < Jk
+ throw(BoundsError(string("row values I[k] must satisfy 1 <= I[k] <= m")))


column indices J[k]…, no? This also covers the TODO comment below, doesn't it?

Fixed, thanks!

Seems not? Both error messages still say "row values" where I agree with @mbauman that one should say "row indices" and the other should say "column indices"

Good catch --- only the I -> J part of the comment registered. Thanks! Fixing now.

Solid now? Thanks!

Sacha0 · 2016-01-27T18:27:38Z

My hunch is that checking the lengths of the input vectors wouldn't be all that expensive compared to the actual work that the function does. But if you've tried that and it's too costly, then it may be a necessary tradeoff.

@mbauman I suspect you are correct, at least for most use cases. I would be happy going either way. My philosophy with the expert driver was to get out of the user's way as much as possible, checks included in case the user was doing something unanticipated, e.g. manipulating tiny sparse matrices in a tight loop. I've added a line to the documentation suggesting testing with --check-bounds=yes to partially address this concern. Thoughts? Thanks!

mbauman · 2016-01-27T23:01:40Z

Given that it's not exported, I'm not too worried here. I think it definitely needs those checks if it ever gets exported, but as it stands it's probably fine.

mbauman · 2016-01-27T23:02:26Z

base/sparse/sparsematrix.jl

+`length(cscrowval) >= nnz(S)` and `length(cscnzval) >= nnz(S)`; hence, if `nnz(S)` is
+unknown at the outset, passing in empty vectors of the appropriate type (`Vector{Ti}()`
+and `Vector{Tv}()` respectively) suffices, or calling the `sparse!` method
+neglecting `cscrowval` and `cscnzval`.


These arrays are currently not optional.

Fixed, thanks!

mbauman · 2016-01-27T23:18:10Z

base/sparse/sparsematrix.jl

+intermediate CSR forms and require `length(csrrowptr) >= m + 1`,
+`length(csrcolval) >= length(I)`, and `length(csrnzval >= length(I))`. Input
+array `klasttouch`, workspace for the second stage, requires `length(klasttouch) >= n`.
+Optional input arrays `csccolptr`, `cscrowval`, and `cscnzval` constitute storage for the


There's another optional here.

The CSC arrays are currently optional; see the method definitions immediately following the main sparse! definition. Thanks!

But there is a method argument style issue with the immediately following definitions. Fixing. Thanks!

tkelman · 2016-02-24T04:29:20Z

test/sparsedir/sparse.jl

@@ -303,7 +310,8 @@ mfe22 = eye(Float64, 2)
 @test_throws ArgumentError sparsevec([3,5,7],[0.1,0.0,3.2],4)

 # issue #5169
-@test nnz(sparse([1,1],[1,2],[0.0,-0.0])) == 0
+# @test nnz(sparse([1,1],[1,2],[0.0,-0.0])) == 0
+# commented following change to sparse() in #14798, also see #12605, #9928, #9906, #6769


This PR seems to have fallen off the radar, I think the only thing we were missing was running some package tests and comparing performance on real use cases.

Also for these tests where you've commented them out, I think it would be preferable to change the tested-for value, maybe leaving a comment that points to these issues as being the reason for the value to be what it is after this change.

This PR seems to have fallen off the radar, I think the only thing we were missing was running some package tests and comparing performance on real use cases.

Thanks for the bump! To clarify, was the request for package testing and benchmarking with e.g. Convex.jl directed at me? Apologies if so, I misinterpreted the request. Also if so, how should I go about package testing? Similarly, what benchmarking with Convex.jl would you like to see? Not being a Convex.jl user, I do not know what would constitute a meaningful benchmark.

Also for these tests where you've commented them out, I think it would be preferable to change the tested-for value, maybe leaving a comment that points to these issues as being the reason for the value to be what it is after this change.

Good call --- I will rebase, touch those tests up, and update the PR. Thanks again!

I can try some benchmarks with my Finite Element code.

I can try some benchmarks with my Finite Element code.

I would be delighted to see those, particularly if your code can benefit from using sparse! in place of sparse. I should have a new version up shortly --- in the process of building and testing. Thanks and best!

No it shouldn't be required for you to run such benchmarks, and convex.jl in particular may have issues on nightly. More directed at other reviewers and people who want to see this merged, we can merge this without doing much benchmarking ahead of time but having some more testing of it in advance would be good to be better informed in case there are any situations where this might end up leading to regressions.

In my code, master and this branch performed almost identically. I don't have time to test the sparse! driver now but it is great it is available!

In my code, master and this branch performed almost identically.

Cheers, thanks!

Sacha0 · 2016-02-24T18:30:38Z

Rebased, touched up tests, and added a note to sparse!'s documentation making it clear that output arrays csccolptr, cscrowval, and cscnzval can alias input arrays I, J, and V in space-constrained environments.

tkelman · 2016-02-24T19:14:24Z

base/sparse/sparsematrix.jl

+`combine` function is used to combine duplicates. If `m` and `n` are not specified, they
+are set to `maximum(I)` and `maximum(J)` respectively. If the `combine` function is not
+supplied, duplicates are added by default. All elements of `I` must satisfy
+`1 <= I[k] <= m`, and all elements of `J` must satisfy `1 <= J[k] <= n`.


there have been updates to this docstring that you need to incorporate.

make sure the signature is consistent between here and the rst docs and run make docs to make sure you aren't undoing recent changes

there have been updates to this docstring that you need to incorporate.

make sure the signature is consistent between here and the rst docs and run make docs to make sure you aren't undoing recent changes

First time playing in the doc sandbox, so apologies in advance for inevitable mistakes. I updated the signature for sparse in doc/stdlib/arrays.rst and issued make julia-genstdlib. No complaints about sparse, though complaints about other things:

UndefVarError(:build_sysimg) WARNING: Mod Base build_sysimg INFO: devdocs/sysimg.rst: no docs for build_sysimg(sysimg_path=default_sysimg_path, cpu_target="native", userimg_path=nothing; force=false) in Base WARNING: Exported method missing doc for Base.StackTraces.StackTrace WARNING: Exported method missing doc for StackFrame WARNING: Exported method missing doc for Base.LibGit2.with WARNING: Exported method missing doc for Base.Docs.doc WARNING: Exported method missing doc for @threadcall WARNING: Missing 5 exported doc strings INFO: Missing 121 unexported doc strings

Guessing that constitutes success, there being no complaints about sparse? Subsequently issued make docs, and then make -C doc html and make -C doc latex for good measure; all seemed to complete fine. Should that do the trick? Will push after clarifying the sparse! docs. Thanks again!

check the diff on the rst and make sure you're incorporating rather than overwriting the recent updates to this docstring

Thanks for the pointer! Diff seems solid?

diff --git a/doc/stdlib/arrays.rst b/doc/stdlib/arrays.rst index 5847201..ee0c291 100644 --- a/doc/stdlib/arrays.rst +++ b/doc/stdlib/arrays.rst @@ -855,11 +855,13 @@ Sparse Vectors and Matrices Sparse vectors and matrices largely support the same set of operations as their dense counterparts. The following functions are specific to sparse arrays. -.. function:: sparse(I,J,V,[m,n,combine]) +.. function:: sparse(I, J, V,[ m, n, combine]) .. Docstring generated from Julia source - Create a sparse matrix ``S`` of dimensions ``m x n`` such that ``S[I[k], J[k]] = V[k]``\ . The ``combine`` function is used to combine duplicates. If ``m`` and ``n`` are not specified, they are set to ``maximum(I)`` and ``maximum(J)`` respectively. If the ``combine`` function is not supplied, ``combine`` defaults to ``+`` unless the elements of ``V`` are Booleans in which case ``combine`` defaults to ``|``\ . All elements of ``I`` must satisfy ``1 <= I[k] <= m``\ , and all elements of ``J`` must satisfy ``1 <= J[k] <= n``\ . + Create a sparse matrix ``S`` of dimensions ``m x n`` such that ``S[I[k], J[k]] = V[k]``\ . The ``combine`` function is used to combine duplicates. If ``m`` and ``n`` are not specified, they are set to ``maximum(I)`` and ``maximum(J)`` respectively. If the ``combine`` function is not supplied, ``combine`` defaults to ``+`` unless the elements of ``V`` are Booleans in which case ``combine`` defaults to ``|``\ . All elements of ``I`` must satisfy ``1 <= I[k] <= m``\ , and all elements of ``J`` must satisfy ``1 <= J[k] <= n``\ . + + For additional documentation and an expert driver, see ``Base.SparseArrays.sparse!``\ .

yep. looks like that isn't pushed yet?

Added the doc modification below and pushed. Thanks!

tkelman · 2016-02-25T03:13:30Z

Ah, we should really document the change in behavior regarding not dropping zero values any more. And maybe it would be a good idea to add new copies of those changed tests but exercising dropzeros!.

Sacha0 · 2016-02-25T04:07:04Z

Ah, we should really document the change in behavior regarding not dropping zero values any more.

Sure, a line at the end of sparse's documentation noting it retains numerical zeros as structural nonzeros? Or something else?

And maybe it would be a good idea to add new copies of those changed tests but exercising dropzeros!

Shall add these tomorrow. Best!

tkelman · 2016-02-25T04:11:52Z

Yeah, in sparse's documentation, and I think this is also worthy of mentioning in NEWS.md.

… version. See JuliaLang#13001 and JuliaLang#14631.

Sacha0 · 2016-02-25T19:00:37Z

Sorry, I accidentally closed this and pushed a revision prior to noticing. GH will not allow me to reopen this pull request post-push, and a brief search seems to indicate there is no recourse. How should I proceed?

Concerning the push, I added a line to sparse's documentation and a mention in NEWS.md regarding the numerical zero retention change. I also added the dropzeros!-exercising test versions mentioned above. Thanks, and best!

StefanKarpinski · 2016-02-25T19:07:20Z

Yeah, looks like it can't be reopened by an admin either so you'll have to make a new PR :-\

Sacha0 · 2016-02-25T19:13:40Z

Opened anew (#15242). Apologies for the glitch!

Sacha0 force-pushed the mitsparse branch from 8e7ceb9 to 8058db3 Compare January 26, 2016 03:18

tkelman added the domain:arrays:sparse Sparse arrays label Jan 26, 2016

Sacha0 force-pushed the mitsparse branch from 8058db3 to 1c03e86 Compare January 26, 2016 04:35

nalimilan reviewed Jan 26, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch from 1c03e86 to 5b14cfb Compare January 26, 2016 23:06

tkelman reviewed Jan 27, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch from 5b14cfb to ad5dae2 Compare January 27, 2016 02:23

mbauman reviewed Jan 27, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch from ad5dae2 to fad055b Compare January 27, 2016 02:39

Sacha0 force-pushed the mitsparse branch from e232c46 to bec6318 Compare January 27, 2016 18:26

mbauman reviewed Jan 27, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch from bec6318 to 5ddfd81 Compare January 27, 2016 23:14

mbauman reviewed Jan 27, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch from 5ddfd81 to 6fa8436 Compare January 27, 2016 23:24

Sacha0 mentioned this pull request Jan 27, 2016

Remove csc_permute and ereach from base/sparse/csparse.jl #12231

Closed

tkelman reviewed Feb 24, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch 2 times, most recently from 5a8b334 to c2fb43e Compare February 24, 2016 18:27

Sacha0 mentioned this pull request Feb 24, 2016

MIT-licensed SparseMatrixCSC fkeep! and children #14702

Merged

tkelman reviewed Feb 24, 2016
View reviewed changes

Sacha0 force-pushed the mitsparse branch from c2fb43e to 3d80e68 Compare February 24, 2016 23:02

tkelman mentioned this pull request Feb 25, 2016

add keepzeros option in sparse #12608

Closed

Replace the LGPL-licensed sparse() parent method with an MIT-licensed…

4496d42

… version. See JuliaLang#13001 and JuliaLang#14631.

Sacha0 force-pushed the mitsparse branch from 3d80e68 to 4496d42 Compare February 25, 2016 18:41

Sacha0 closed this Feb 25, 2016

Sacha0 mentioned this pull request Feb 25, 2016

MIT-licensed sparse() parent method and expert driver, take two #15242

Merged

matthieugomez mentioned this pull request Mar 20, 2016

A keepzeros option for sparse(I, J, V) #12605

Closed

This was referenced Jun 14, 2016

MIT-licensed SparseMatrixCSC permute[!] and refactored [c]transpose[!] #16931

Merged

Test, document, and export dropzeros[!] #16947

Merged

MIT-licensed sparse() parent method and expert driver #14798

MIT-licensed sparse() parent method and expert driver #14798

Conversation

Sacha0 commented Jan 26, 2016 • edited Loading

hayd commented Jan 26, 2016

Sacha0 commented Jan 26, 2016

KristofferC commented Jan 26, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

StefanKarpinski commented Jan 26, 2016

nalimilan commented Jan 26, 2016

StefanKarpinski commented Jan 26, 2016

ViralBShah commented Jan 26, 2016

tkelman commented Jan 26, 2016

StefanKarpinski commented Jan 26, 2016

Sacha0 commented Jan 26, 2016

Sacha0 commented Jan 26, 2016

Sacha0 commented Jan 26, 2016

tkelman commented Jan 27, 2016

Sacha0 commented Jan 27, 2016

tkelman commented Jan 27, 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

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 Jan 27, 2016

mbauman commented Jan 27, 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

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 Feb 24, 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkelman commented Feb 25, 2016

Sacha0 commented Feb 25, 2016

tkelman commented Feb 25, 2016

Sacha0 commented Feb 25, 2016

StefanKarpinski commented Feb 25, 2016

Sacha0 commented Feb 25, 2016

Sacha0 commented Jan 26, 2016 •

edited

Loading