MIT-licensed SparseMatrixCSC fkeep! and children #14702

Sacha0 · 2016-01-17T01:03:11Z

Followup to #13001 and #14631. This pull request replaces the LGPL-licensed SparseMatrixCSC fkeep! method and children tril!, triu!, droptol!, and dropzeros[!] (noticed signature during cleanup) in base/sparse/csparse.jl with MIT-licensed versions.

These methods run in O(A.n, nnz(A)) time and require no space beyond that passed in. The new methods perform better than the old methods, particularly for large matrices; are the new implementations missing anything that would require additional time?

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.

Benchmarks:

using Benchmarks: @benchmark, SummaryStatistics

function prettytimes(bench)
    stats = SummaryStatistics(bench)
    timecenter = stats.elapsed_time_center
    timelower = get(stats.elapsed_time_lower)
    timeupper = get(stats.elapsed_time_upper)
    # based on Benchmarks.pretty_time_string
    timecenter < 1_000.0 ? (scalefactor = 1.0; units = "ns") :
        timecenter < 1_000_000.0 ? (scalefactor = 1_000.0; units = "μs") :
            timecenter < 1_000_000_000.0 ? (scalefactor = 1_000_000.0; units = "ms") :
                (scalefactor = 1_000_000_000.0; units = " s")
    @sprintf("%6.2f %s [%6.2f,%6.2f]", timecenter/scalefactor, units, timelower/scalefactor, timeupper/scalefactor)
end

smallN, smallerN = 600, 400;
smallsqrA = sprand(smallN, smallN, 0.01);
smallrectA = sprand(smallN, smallerN, 0.01);
largeN, largerN = 100000, 200000;
largesqrA = sprand(largeN, largeN, 0.001);
largerectA = sprand(largeN, largerN, 0.001);

wm, wt = 12, 25;
println(" $(lpad("method ", wm)) $(lpad("small square A", wt)), $(lpad("small rect A", wt)) | $(lpad("large square A", wt)), $(lpad("large rect A", wt)) ")
@printf("%s: %s , %s | %s, %s\n", lpad("tril!", wm),
    prettytimes(@benchmark triu!(copy(smallsqrA))),
    prettytimes(@benchmark triu!(copy(smallrectA))),
    prettytimes(@benchmark triu!(copy(largesqrA))),
    prettytimes(@benchmark triu!(copy(largerectA))) )
@printf("%s: %s , %s | %s, %s\n", lpad("triu!", wm),
    prettytimes(@benchmark triu!(copy(smallsqrA))),
    prettytimes(@benchmark triu!(copy(smallrectA))),
    prettytimes(@benchmark triu!(copy(largesqrA))),
    prettytimes(@benchmark triu!(copy(largerectA))) )
@printf("%s: %s , %s | %s, %s\n", lpad("droptol!", wm),
    prettytimes(@benchmark Base.SparseArrays.droptol!(copy(smallsqrA), 0.5)),
    prettytimes(@benchmark Base.SparseArrays.droptol!(copy(smallrectA), 0.5)),
    prettytimes(@benchmark Base.SparseArrays.droptol!(copy(largesqrA), 0.5)),
    prettytimes(@benchmark Base.SparseArrays.droptol!(copy(largerectA), 0.5)) )
@printf("%s: %s , %s | %s, %s\n", lpad("dropzeros!", wm),
    prettytimes(@benchmark Base.SparseArrays.dropzeros!(copy(smallsqrA))),
    prettytimes(@benchmark Base.SparseArrays.dropzeros!(copy(smallrectA))),
    prettytimes(@benchmark Base.SparseArrays.dropzeros!(copy(largesqrA))),
    prettytimes(@benchmark Base.SparseArrays.dropzeros!(copy(largerectA))) )

On master:

    method             small square A,              small rect A |            large square A,              large rect A
     tril!:   5.69 μs [  5.62,  5.75] ,   2.50 μs [  2.48,  2.52] |  16.45 ms [ 15.39, 17.51],  47.47 ms [ 45.94, 49.01]
     triu!:   5.66 μs [  5.61,  5.71] ,   2.46 μs [  2.45,  2.47] |  16.32 ms [ 15.33, 17.31],  47.88 ms [ 45.88, 49.89]
  droptol!:   4.64 μs [  4.61,  4.67] ,   3.09 μs [  3.05,  3.13] |  14.32 ms [  9.82, 18.83],  29.77 ms [ 16.43, 43.10]
dropzeros!:  13.10 μs [ 13.07, 13.14] ,   8.57 μs [  8.53,  8.62] |  25.05 ms [ 24.10, 25.99],  50.34 ms [ 49.09, 51.59]

On this PR's branch:

    method             small square A,              small rect A |            large square A,              large rect A
     tril!:   4.31 μs [  4.30,  4.32] ,   1.60 μs [  1.60,  1.61] |   6.51 ms [  5.26,  7.76],  18.19 ms [ 15.06, 21.32]
     triu!:   4.31 μs [  4.29,  4.33] ,   1.60 μs [  1.59,  1.60] |   6.51 ms [  5.27,  7.74],  17.90 ms [ 15.17, 20.64]
  droptol!:   4.47 μs [  4.45,  4.50] ,   2.17 μs [  2.16,  2.19] |   8.59 ms [  3.83, 13.34],  18.41 ms [  4.40, 32.42]
dropzeros!:  10.50 μs [ 10.42, 10.59] ,   6.48 μs [  6.46,  6.51] |  11.45 ms [ 10.88, 12.01],  22.98 ms [ 22.15, 23.81]

The dropzeros! benchmark isn't particularly meaningful given the absence of stored zeros in the test matrices. Thanks, and best!

Sacha0 · 2016-01-17T02:36:00Z

Travis discontent unrelated?

tkelman · 2016-01-17T03:37:09Z

yeah

fatal error during bootstrap:
(parse-error "read: unexpected ']'")
ERROR: fatal error loading system image

there must be a race condition in a makefile or the parser or something...

ViralBShah · 2016-01-17T08:07:58Z

Nice performance gain. Would you be able to add some tests as well? Doesn't seem like there are any at the moment.

Sacha0 · 2016-01-17T19:29:41Z

Would you be able to add some tests as well? Doesn't seem like there are any at the moment.

tril! and triu! tests exist in 876240c/test/sparsedir/sparse.jl#L974-L975. A droptol! test exists in 876240c/test/sparsedir/sparse.jl#L956-L957. The preceding tests implicitly test fkeep!. dropzeros! appears to lack tests.

I added a test for dropzeros!. Are the existing tests for tril!, triu!, droptol!, and thereby fkeep! adequate, or would you like additional tests? If the latter, what would you like to see? Thanks for the review!

hayd · 2016-01-24T20:08:21Z

lgtm

ViralBShah · 2016-01-24T20:43:28Z

@tkelman Merge this?

tkelman · 2016-01-25T00:41:43Z

base/sparse/sparsematrix.jl

+triu!(A::SparseMatrixCSC, k::Integer = 0) = fkeep!(A, TriuFunc(), k)
+
+immutable DroptolFunc <: Base.Func{4} end
+call{Tv,Ti}(::DroptolFunc, i::Ti, j::Ti, x::Tv, tol::Real) = abs(x) >= tol


previous implementation of droptol used > here

Fixed, thanks!

hayd · 2016-01-27T04:06:44Z

@tkelman Are more tests needed here?

tkelman · 2016-01-27T04:16:39Z

Yes, aside from the clarification about "attach" each of the differences I highlighted means a missing corner case that should be tested against (and would fail until fixed).

tkelman · 2016-01-27T22:42:34Z

base/sparse/sparsematrix.jl

+call{Tv,Ti}(::TrilFunc, i::Ti, j::Ti, x::Tv, k::Integer) = i + k >= j
+call{Tv,Ti}(::TriuFunc, i::Ti, j::Ti, x::Tv, k::Integer) = j >= i + k
+function tril!(A::SparseMatrixCSC, k::Integer = 0)
+    if k >= A.n - 1


the old boundserror conditions were, effectively, k > n || -k > m

Do you prefer the old checks? I would be happy with either. The new checks are what I came up with in vacuum on the basis of logical consistency. I imagine you saw my comment above?

Added bounds checking on k. When adding tests for said bounds checks, I noticed (by way of existing tests) that the old checks may not be as tight as the new checks and throw BoundsErrors rather than ArgumentErrors. Thoughts on the new checks? Thanks!

edit: Sorry. Ignore this. Posted prior to seeing your response below. Thanks!

Sacha0 · 2016-01-27T22:42:44Z

each of the differences I highlighted means a missing corner case that should be tested against

Added tests for these edge cases. Thanks!

tkelman · 2016-01-27T22:49:35Z

ref

julia/base/linalg/dense.jl

Line 91 in 91dcdaa

if (k > 0 && k > n) || (k < 0 && -k > m)

for the current bounds check conditions on dense tril! - if we're going to offset this check by one (or two?), we should be consistent between dense and sparse.

edit: may as well make the error message consistent as well

Sacha0 · 2016-01-27T22:51:15Z

we should be consistent between dense and sparse.

Nice catch! Shall fix.

Sacha0 · 2016-01-27T23:07:34Z

Actually, the k > n || -k > m test seems inconsistent with the associated error message? An m-by-n matrix has only n-1 superdiagonals and m-1 subdiagonals, so k == n and k == -m are already out of bounds. Hence the checks I wrote a priori. Am I missing something? If not, perhaps it makes more sense to revise the existing k > n || -k > m tests? Thanks!

tkelman · 2016-01-27T23:39:20Z

Our last few comments have been crossing paths. It does probably make sense to offset the current dense checks by one (though would need to check all the other tril[!] and triu[!] methods too), but at the very least in a separate commit. And the asymmetry of having one check inclusive and the other exclusive is maybe questionable.

Sacha0 · 2016-01-28T02:52:52Z

And the asymmetry of having one check inclusive and the other exclusive is maybe questionable.

Agreed. I would be happy to symmetrize the checks. My reasoning for the asymmetric checks was as follows: tril!(A, A.n - 1) and triu!(A, 1 - A.m) are expensive identities and hence most likely logical errors, whereas tril!(A, 1 - A.m) and triu!(A, A.n - 1) are legitimate (assuming I did not mix that up).

Our last few comments have been crossing paths. It does probably make sense to offset the current dense checks by one (though would need to check all the other tril[!] and triu[!] methods too), but at the very least in a separate commit.

Cheers, proposal for action on my part: Keep the offset checks (symmetric, asymmetric, whichever you prefer) in this pull request. Open an issue regarding consistency of tri(u|l)[!] k-checks. On consensus regarding desired behavior, create a new pull request revising either the methods from this pull request or the other, presently-existing methods in accord with consensus. Thoughts? Thanks again!

Sacha0 · 2016-01-28T02:54:49Z

Travis x86_64 unhappiness seemingly unrelated?

tkelman · 2016-01-28T03:03:19Z

are expensive identities and hence most likely logical errors

I don't think it's our job to guess what the user's intent was. It may be a silly thing to ask for, but it's valid and within the bounds of the array so I don't think an exception is called for if you ask for the entire array.

Sacha0 · 2016-01-28T05:12:29Z

Checks symmetrized. Thanks! Edit: And fixed tests.

tkelman · 2016-01-28T06:13:19Z

32 bit travis failure is a new one, looks unrelated though:

julia: /home/travis/build/JuliaLang/julia/src/gc.c:607: find_region: Assertion `maybe && "find_region failed"' failed.

    From worker 3:       * docs                 Worker 3 terminated.

ERROR (unhandled task failure): EOFError: read end of file

Sacha0 · 2016-01-28T16:48:31Z

AppVeyor i686 failures also seem unrelated?

Sacha0 · 2016-02-24T19:01:48Z

Rebased. Is this in shape to merge? Did I miss anything above? Thanks! (Prompted by #14798 (comment).)

tkelman · 2016-02-24T19:09:51Z

This probably doesn't need to use functors any more?

Sacha0 · 2016-02-24T19:20:00Z

This probably doesn't need to use functors any more?

True, and likewise with sparse! and children and qftranspose! and children :). If acceptable, I would prefer to merge this pull request and #14798 as they are (with established good performance), and later carefully clean up all those methods in one shot. Does that sound good? Thanks!

tkelman · 2016-02-24T22:53:53Z

Yeah that makes sense.

tkelman · 2016-02-24T22:54:59Z

base/sparse/sparsematrix.jl

+
+immutable TrilFunc <: Base.Func{4} end
+immutable TriuFunc <: Base.Func{4} end
+call{Tv,Ti}(::TrilFunc, i::Ti, j::Ti, x::Tv, k::Integer) = i + k >= j


oh, but call is now deprecated

I can put the functors-to-functions PR together in the not-too-distant future, certainly within the deprecation cycle. Would that be alright? Thanks!

This can still use functors, but it should do so using the non-deprecated syntax to avoid introducing warnings to the tests. The tests should currently be running with --depwarn=error on CI so I'm surprised this didn't fail.

Thanks for the pointer; I was unaware of the new syntax. Assuming I found the right issues, this should be fixed now (?). Thanks!

…tril!, triu!, droptol!, and dropzeros[!] with MIT-licensed versions. See JuliaLang#13001 and JuliaLang#14631. Also add a test for dropzeros!.

Sacha0 · 2016-02-25T18:26:17Z

The Appveyor i686 discontent (failure on 7f845ff/test/arrayops.jl#L1060) seems unrelated?

MIT-licensed SparseMatrixCSC fkeep! and children

KristofferC · 2016-02-26T06:39:41Z

Nice!

Sacha0 · 2016-02-26T16:33:47Z

@tkelman Much thanks for the thorough review here as well!

tkelman added the sparse Sparse arrays label Jan 17, 2016

Sacha0 force-pushed the fkeepvrienden branch from 0d970aa to 786f3f7 Compare January 17, 2016 19:29

Sacha0 mentioned this pull request Jan 17, 2016

Intermittent Travis build failure? #14705

Closed

jrevels added the potential benchmark Could make a good benchmark in BaseBenchmarks label Jan 24, 2016

tkelman reviewed Jan 25, 2016
View reviewed changes

tkelman added the needs tests Unit tests are required for this change label Jan 25, 2016

Sacha0 mentioned this pull request Jan 26, 2016

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

Closed

Sacha0 force-pushed the fkeepvrienden branch from 786f3f7 to c56da1d Compare January 27, 2016 22:35

tkelman reviewed Jan 27, 2016
View reviewed changes

tkelman removed the needs tests Unit tests are required for this change label Jan 27, 2016

Sacha0 force-pushed the fkeepvrienden branch from c56da1d to a783405 Compare January 27, 2016 22:45

Sacha0 mentioned this pull request Jan 27, 2016

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

Closed

Sacha0 force-pushed the fkeepvrienden branch from a783405 to fe77fc2 Compare January 28, 2016 05:11

Sacha0 force-pushed the fkeepvrienden branch 2 times, most recently from 7101000 to 58a6799 Compare January 28, 2016 05:17

Sacha0 force-pushed the fkeepvrienden branch from 58a6799 to 83352ee Compare February 24, 2016 19:00

tkelman reviewed Feb 24, 2016
View reviewed changes

Replace the LGPL-licensed SparseMatrixCSC fkeep! method and children …

f893c06

…tril!, triu!, droptol!, and dropzeros[!] with MIT-licensed versions. See JuliaLang#13001 and JuliaLang#14631. Also add a test for dropzeros!.

Sacha0 force-pushed the fkeepvrienden branch from 83352ee to f893c06 Compare February 24, 2016 23:49

tkelman added a commit that referenced this pull request Feb 26, 2016

Merge pull request #14702 from Sacha0/fkeepvrienden

31e041e

MIT-licensed SparseMatrixCSC fkeep! and children

tkelman merged commit 31e041e into JuliaLang:master Feb 26, 2016

Sacha0 deleted the fkeepvrienden branch February 26, 2016 16:33

This was referenced Jun 14, 2016

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

Merged

taking sub/superdiagonal argument checks too seriously #16959

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

MIT-licensed SparseMatrixCSC fkeep! and children #14702

MIT-licensed SparseMatrixCSC fkeep! and children #14702

Sacha0 commented Jan 17, 2016

Sacha0 commented Jan 17, 2016

tkelman commented Jan 17, 2016

ViralBShah commented Jan 17, 2016

Sacha0 commented Jan 17, 2016

hayd commented Jan 24, 2016

ViralBShah commented Jan 24, 2016

tkelman Jan 25, 2016

Sacha0 Jan 27, 2016

hayd commented Jan 27, 2016

tkelman commented Jan 27, 2016

tkelman Jan 27, 2016

Sacha0 Jan 27, 2016

Sacha0 commented Jan 27, 2016

tkelman commented Jan 27, 2016

Sacha0 commented Jan 27, 2016

Sacha0 commented Jan 27, 2016

tkelman commented Jan 27, 2016

Sacha0 commented Jan 28, 2016

Sacha0 commented Jan 28, 2016

tkelman commented Jan 28, 2016

Sacha0 commented Jan 28, 2016

tkelman commented Jan 28, 2016

Sacha0 commented Jan 28, 2016

Sacha0 commented Feb 24, 2016

tkelman commented Feb 24, 2016

Sacha0 commented Feb 24, 2016

tkelman commented Feb 24, 2016

tkelman Feb 24, 2016

Sacha0 Feb 24, 2016

tkelman Feb 24, 2016

Sacha0 Feb 24, 2016

Sacha0 commented Feb 25, 2016

KristofferC commented Feb 26, 2016

Sacha0 commented Feb 26, 2016

MIT-licensed SparseMatrixCSC fkeep! and children #14702

MIT-licensed SparseMatrixCSC fkeep! and children #14702

Conversation

Sacha0 commented Jan 17, 2016

Sacha0 commented Jan 17, 2016

tkelman commented Jan 17, 2016

ViralBShah commented Jan 17, 2016

Sacha0 commented Jan 17, 2016

hayd commented Jan 24, 2016

ViralBShah commented Jan 24, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hayd commented Jan 27, 2016

tkelman commented Jan 27, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Sacha0 commented Jan 27, 2016

tkelman commented Jan 27, 2016

Sacha0 commented Jan 27, 2016

Sacha0 commented Jan 27, 2016

tkelman commented Jan 27, 2016

Sacha0 commented Jan 28, 2016

Sacha0 commented Jan 28, 2016

tkelman commented Jan 28, 2016

Sacha0 commented Jan 28, 2016

tkelman commented Jan 28, 2016

Sacha0 commented Jan 28, 2016

Sacha0 commented Feb 24, 2016

tkelman commented Feb 24, 2016

Sacha0 commented Feb 24, 2016

tkelman 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

Sacha0 commented Feb 25, 2016

KristofferC commented Feb 26, 2016

Sacha0 commented Feb 26, 2016