stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

dkarrasch · 2018-09-21T07:27:35Z

Fixes #29111. Previously, lmul! and rmul! had no special method for multiplying a sparse matrix with a diagonal matrix in-place, but fell back to the generic LinearAlgebra/Diagonal-method. Tests had been included before.

Add rmul! and lmul! of sparse matrix with Diagonal

StefanKarpinski

Thank you for the nice contribution! I'm not the best person to review this but it looks good to me.

StefanKarpinski · 2018-09-21T13:57:22Z

stdlib/SparseArrays/src/linalg.jl

+    (n == size(D, 1)) || throw(DimensionMismatch())
+    Anzval = A.nzval
+    for col = 1:n, p = A.colptr[col]:(A.colptr[col+1]-1)
+        @inbounds Anzval[p] *= D.diag[col]


It might be possible to gain a bit of speed here by manually hoisting the D.diag[col] lookup outside of the inner loop. Of course, the compiler may well be able to do that automatically. Depends on how much it knows about whether D.diag can change or move.

From my experience with the SparseArray code, this is only needed when the object one getfield into is mutable. There is a lot of hoisting out fields for SparseMatrixCSC that was put there before #16371 (see discussion in #15668) but it shouldn't be needed now. Always worth benchmarking though...

I'm sorry I can't follow. I did the change @StefanKarpinski suggested. Was that just (potentially) unnecessary or even "wrong"?

I benchmarked the previous ("my") against the latest version ("Stefan's") on these:

P = sprand(12000, 12000, 0.01) q = Diagonal(dropdims(sum(P, dims=2), dims=2))

and Stefan's won with 607.996 μs vs. 631.496 μs. So I guess it was worth it.

Hm, might be neglible. For

P = sprand(120000, 120000, 0.01)

(an order of magnitude larger) and q as before I get 104.832 ms vs. 103.934 ms.

The question is if the compiler tries to reload D.diag` every time in the loop or if it realizes that this will never change and hoist the load outside the loop.

Seems like it doesn't reload, at least it doesn't show up in performance loss. Should I, for better readability, then merge the two loops again as before, but leave the @inbounds upfront like in the last commit?

Doesn't matter that much but perhaps it makes sense to at least have the two methods in this PR be consistent with eachother.

dkarrasch · 2018-09-28T13:10:32Z

Is there anything that I can do more to get this merged? This is my first PR to Julia, so that would make me really proud. Maybe @andreasnoack may want to have a look at it? I can't seem to make much sense of the CI failures...

mcognetta · 2018-09-28T13:55:49Z

stdlib/SparseArrays/src/linalg.jl

 function lmul!(b::Number, A::SparseMatrixCSC)
    lmul!(b, A.nzval)
    return A
 end

+function rmul!(A::SparseMatrixCSC, D::Diagonal{T}) where T


T is never used in the function body so I think you can make this signature

function rmul!(A::SparseMatrixCSC, D::Diagonal). (same for lmul)

mcognetta · 2018-09-28T13:57:45Z

Do you mind giving a small case where rmul and lmul would be called? Just doing regular multiplication doesn't hit those methods.

Also the CI failures are likely not because of your codet and just need to be restarted.

dkarrasch · 2018-09-28T14:45:54Z

One use case is shown in #29111 : given a sparse "weight matrix", you may wish to turn it into a row-stochastic matrix (graph Laplacian) by row-normalizing by the sum of each row:

P = sprand(12000, 12000, 0.01)
q = Diagonal(dropdims(sum(P, dims=2), dims=2))
rmul!(P, q)

Very often, the weight matrix is huge and therefore sparse, and row-normalization corresponds to the rmul operation.

Sacha0

These changes look great! Thanks and welcome @dkarrasch! :)

Perhaps add tests exercising these new methods?

dkarrasch · 2018-09-28T19:05:16Z

Thank you! I thought tests already existed at

julia/stdlib/SparseArrays/test/sparse.jl

Line 375 in e1f42b1

@test dA * Diagonal(b) == rmul!(copy(sA), Diagonal(b))

but after a more thorough look I see that the result of rmul!´ and lmul` (a few lines below) were tested against a dense array. I'll change the above line then to

@test sA * Diagonal(b) == rmul!(copy(sA), Diagonal(b))

and similarly for lmul!, as I don't see any reason why one would want to have the result be dense in any case.

Sacha0 · 2018-09-28T19:14:23Z

Thanks for pointing the existing tests out! The existing tests look good to me. I imagine the rationale for testing against the result of an equivalent dense operation is that the relevant code paths are more likely to be disjoint, whereas e.g. sA * Diagonal(b) could quite well be implemented in terms of the mutating calls in question. Best!

dkarrasch · 2018-09-28T19:14:49Z

I was not aware of these details (or better, never thought about it), but the density/sparsity doesn't matter in the test. So even though dA * Diagonal(b) returns a dense array and rmul!(copy(sA), Diagonal(b)) a sparse one (in this PR), this is all fine as long as they give the same matrix eventually. So I think the existing tests cover both rmul! and lmul! already.

Sacha0

Thanks @dkarrasch! :)

(Absent objections, requests for time, or someone beating me to it, I plan to merge these changes tomorrow morning PT or later.)

Sacha0 · 2018-09-29T18:39:19Z

Once again thanks and welcome @dkarrasch! :)

chriscoey · 2018-09-29T19:14:43Z

Thank you thank you! This improves the speed of my code

dkarrasch · 2018-09-29T20:30:17Z

I know, mine, too! :-) But all I did was tracing back what scale! used to do in Julia v0.6, when it had a special method for sparse and diagonal. The change was really minor, after all. So, when and how is this "delivered" to the people officially? With Julia v1.0.2, since it's in a stdlib?

Sacha0 · 2018-09-30T17:42:19Z

So, when and how is this "delivered" to the people officially? With Julia v1.0.2, since it's in a stdlib?

Good question! I imagine these changes are a backport candidate. @ararslan?

dkarrasch · 2018-11-19T17:44:05Z

@ararslan, is there any chance for a backport?

ararslan · 2018-11-19T18:54:32Z

Yep should be possible to put into 1.0.3. cc @KristofferC

…al (#29296) * Add rmul! and lmul! of sparse matrix with Diagonal (cherry picked from commit 9454c8d)

dkarrasch added 2 commits September 21, 2018 09:14

Add rmul! and lmul! of sparse matrix with Diagonal

0537f6b

Merge pull request #1 from dkarrasch/dk_rmul_lmul_sp_diag

b93376c

Add rmul! and lmul! of sparse matrix with Diagonal

StefanKarpinski reviewed Sep 21, 2018

View reviewed changes

dkarrasch added 2 commits September 21, 2018 16:40

Some optimization

43af178

merge loops again for consistency

2efc646

dkarrasch changed the title ~~Add rmul! and lmul! of sparse matrix with Diagonal~~ stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal Sep 23, 2018

chriscoey mentioned this pull request Sep 25, 2018

Alfonso matlab code comparison jump-dev/Hypatia.jl#23

Closed

mcognetta reviewed Sep 28, 2018

View reviewed changes

remove unused parametrized type

2edf05f

fredrikekre added linear algebra Linear algebra sparse Sparse arrays labels Sep 28, 2018

fredrikekre requested review from andreasnoack and Sacha0 September 28, 2018 15:38

Sacha0 reviewed Sep 28, 2018

View reviewed changes

Sacha0 approved these changes Sep 28, 2018

View reviewed changes

Sacha0 merged commit 9454c8d into JuliaLang:master Sep 29, 2018

ararslan added the backport pending 1.0 label Nov 19, 2018

dkarrasch mentioned this pull request Nov 28, 2018

[SparseArrays] Respect order of mul in (l)mul!(::Diagonal,::Sparse) #30163

Merged

KristofferC pushed a commit that referenced this pull request Nov 28, 2018

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagon…

ef19ed2

…al (#29296) * Add rmul! and lmul! of sparse matrix with Diagonal (cherry picked from commit 9454c8d)

KristofferC mentioned this pull request Nov 28, 2018

Backports for 1.0.3 #30010

Merged

61 tasks

KristofferC pushed a commit that referenced this pull request Dec 12, 2018

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagon…

e77ebd2

…al (#29296) * Add rmul! and lmul! of sparse matrix with Diagonal (cherry picked from commit 9454c8d)

KristofferC added performance Must go faster and removed backport pending 1.0 labels Dec 12, 2018

KristofferC pushed a commit that referenced this pull request Feb 11, 2019

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagon…

ca50865

…al (#29296) * Add rmul! and lmul! of sparse matrix with Diagonal (cherry picked from commit 9454c8d)

KristofferC pushed a commit that referenced this pull request Feb 20, 2020

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagon…

5ad620b

…al (#29296) * Add rmul! and lmul! of sparse matrix with Diagonal (cherry picked from commit 9454c8d)

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

dkarrasch commented Sep 21, 2018

StefanKarpinski left a comment

StefanKarpinski Sep 21, 2018

KristofferC Sep 21, 2018

dkarrasch Sep 21, 2018

dkarrasch Sep 21, 2018

dkarrasch Sep 21, 2018

KristofferC Sep 21, 2018

dkarrasch Sep 21, 2018

KristofferC Sep 21, 2018

dkarrasch commented Sep 28, 2018

mcognetta Sep 28, 2018

mcognetta commented Sep 28, 2018

dkarrasch commented Sep 28, 2018

Sacha0 left a comment

dkarrasch commented Sep 28, 2018

Sacha0 commented Sep 28, 2018

dkarrasch commented Sep 28, 2018

Sacha0 left a comment

Sacha0 commented Sep 29, 2018

chriscoey commented Sep 29, 2018

dkarrasch commented Sep 29, 2018

Sacha0 commented Sep 30, 2018

dkarrasch commented Nov 19, 2018

ararslan commented Nov 19, 2018

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

stdlib/SparseArrays: add rmul! and lmul! of sparse matrix with Diagonal #29296

Conversation

dkarrasch commented Sep 21, 2018

StefanKarpinski left a 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

dkarrasch commented Sep 28, 2018

Choose a reason for hiding this comment

mcognetta commented Sep 28, 2018

dkarrasch commented Sep 28, 2018

Sacha0 left a comment

Choose a reason for hiding this comment

dkarrasch commented Sep 28, 2018

Sacha0 commented Sep 28, 2018

dkarrasch commented Sep 28, 2018

Sacha0 left a comment

Choose a reason for hiding this comment

Sacha0 commented Sep 29, 2018

chriscoey commented Sep 29, 2018

dkarrasch commented Sep 29, 2018

Sacha0 commented Sep 30, 2018

dkarrasch commented Nov 19, 2018

ararslan commented Nov 19, 2018