add in-place kron #31069

Roger-luo · 2019-02-14T07:05:29Z

This is useful to avoid extra allocation when A, B is large in kron(A, B). Since the implementation itself is exactly the same with original kron, I guess it's better to have it directly in stdlibs

stdlib/LinearAlgebra/src/dense.jl

martinholters · 2019-02-14T12:58:37Z

stdlib/LinearAlgebra/src/diagonal.jl

 end

-function kron(A::AbstractMatrix{T}, B::Diagonal{S}) where {T<:Number, S<:Number}
+@inline function kron!(C::AbstractMatrix, A::AbstractMatrix, B::Diagonal)


This should also have a @propagate_inbounds, I suppose?

ararslan · 2019-02-14T18:29:53Z

stdlib/LinearAlgebra/src/dense.jl

+
+Like [`kron`](@ref), but stores the results in `C`.
+
+!!! tip


Is this a recognized Documenter thing? I might go with !!! info to be safe.

Yes, tip will be a green block, but info is good as well.

Oh nice, I didn't know that. Either is fine then.

stdlib/LinearAlgebra/src/dense.jl

dkarrasch · 2020-05-06T07:14:29Z

It would be nice to get this finished. @Roger-luo Could you please rebase this onto current master and resolve the merge conflicts? One thing I noticed is that there are occurrences of @inline @Base.propagate_inbounds. If I'm not mistaken, then the second macro already implies that the function is inlined, at least according to its docstring, so the first one would be redundant.

stdlib/SparseArrays/src/linalg.jl

dkarrasch

Just a few comments that I noticed. Will take a more careful look later.

dkarrasch · 2020-05-06T12:29:18Z

stdlib/SparseArrays/src/linalg.jl

-            ptr_range = (1:lB) .+ (colptrC[col]-1)
-            colptrC[col+1] = colptrC[col] + lA*lB
+            ptr_range = (1:lB) .+ (C.colptr[col]-1)
+            C.colptr[col+1] = C.colptr[col] + lA*lB


Wait, there is no C defined in this function...? Or did you intend to rename this to kron! and add a first C argument here? Independent from that, I think one is no longer supposed to refer to colptr via the dot notation, but rather via getcolptr(C), in case you'll need it. Otherwise, CI will throw and fail.

dkarrasch · 2020-05-06T12:31:23Z

stdlib/SparseArrays/src/linalg.jl

+    return z
+end
+
+kron!(C::SparseMatrixCSC, A::SparseMatrixCSC, x::SparseVector) = kron!(C, A, SparseMatrixCSC(x))


Here and below, I think all the SparseMatrixCSC can be relaxed to AbstractSparseMatrixCSC.

dkarrasch · 2020-05-06T12:32:36Z

stdlib/LinearAlgebra/src/dense.jl

+    Bounds checking can be disabled by `@inbounds`, but you need to take care of the shape
+    of `C`, `A`, `B` yourself.
+"""
+@inline Base.@propagate_inbounds function kron!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)


Suggested change

@inline Base.@propagate_inbounds function kron!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

Base.@propagate_inbounds function kron!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

should inline already, same below.

Roger-luo · 2020-05-07T11:02:06Z

Hi sorry, I found the kron seems to changed a lot since last year, thus I rewrite this PR and made force push to this branch (so I don't need to rebase I guess) Let me know if this looks good now.

I also polished the Base.@propagate_inbounds part of the kron!, but I'm actually not sure this is right way to do this, since other in place functions such as generic_matmul! etc. actually just leave the bounds checks inside the implementation. Maybe I should just remove all the Base.@propagate_inbounds in this PR?

elisno · 2020-05-07T23:00:04Z

This is impressive work!

I've been looking for ways to repeatedly compute large Kronecker sums with a fixed structure. kron! will definitely help with this task.

Below are some benchmarks showing the speedup kron! offers when calculating a single term in a Kronecker sum.

Below are some benchmarks for calculating a single term in a Kronecker sum.
kron!(C::AbstractMatrix, A::AbstractMatrix, B::Diagonal) gives a 36x speedup for dense products!

Dense

using BenchmarkTools
using LinearAlgebra
using SparseArrays

for N in 2 .^ [0,1,2,3,4,5,6,7]
    A = rand(N,N); B = Diagonal(oneunit(A)); C = kron(A,B);
    println("Out-of-place N=$N")
    @btime kron($A,$B)
    println("In-place N=$N")
    @btime kron!($C,$A,$B)
    println("")
end

Click to expand timings

Out-of-place N=1
  45.244 ns (1 allocation: 96 bytes)
In-place N=1
  15.453 ns (0 allocations: 0 bytes)

Out-of-place N=2
  60.918 ns (1 allocation: 208 bytes)
In-place N=2
  25.499 ns (0 allocations: 0 bytes)

Out-of-place N=4
  569.928 ns (1 allocation: 2.13 KiB)
In-place N=4
  87.182 ns (0 allocations: 0 bytes)

Out-of-place N=8
  4.079 μs (2 allocations: 32.08 KiB)
In-place N=8
  617.090 ns (0 allocations: 0 bytes)

Out-of-place N=16
  64.967 μs (2 allocations: 512.08 KiB)
In-place N=16
  13.821 μs (0 allocations: 0 bytes)

Out-of-place N=32
  1.149 ms (2 allocations: 8.00 MiB)
In-place N=32
  162.917 μs (0 allocations: 0 bytes)

Out-of-place N=64
  53.302 ms (2 allocations: 128.00 MiB)
In-place N=64
  2.081 ms (0 allocations: 0 bytes)

Out-of-place N=128
  754.902 ms (2 allocations: 2.00 GiB)
In-place N=128
  20.843 ms (0 allocations: 0 bytes)

Sparse

p = 1.0
for N in 2 .^ [0,1,2,3,4,5,6,7]
    A = sprand(N,N,p); B = oneunit(A); C = kron(A,B);
    println("Out-of-place N=$N")
    @btime kron($A,$B)
    println("In-place N=$N")
    @btime kron!($C,$A,$B)
    println("")
end

Click to expand timings

Out-of-place N=1
  132.470 ns (3 allocations: 288 bytes)
In-place N=1
  27.849 ns (0 allocations: 0 bytes)

Out-of-place N=2
  200.434 ns (3 allocations: 416 bytes)
In-place N=2
  83.813 ns (0 allocations: 0 bytes)

Out-of-place N=4
  659.133 ns (3 allocations: 1.44 KiB)
In-place N=4
  435.384 ns (0 allocations: 0 bytes)

Out-of-place N=8
  4.371 μs (3 allocations: 8.86 KiB)
In-place N=8
  2.997 μs (0 allocations: 0 bytes)

Out-of-place N=16
  29.972 μs (5 allocations: 66.34 KiB)
In-place N=16
  22.710 μs (0 allocations: 0 bytes)

Out-of-place N=32
  226.747 μs (5 allocations: 520.34 KiB)
In-place N=32
  182.200 μs (0 allocations: 0 bytes)

Out-of-place N=64
  1.797 ms (6 allocations: 4.03 MiB)
In-place N=64
  1.420 ms (0 allocations: 0 bytes)

Out-of-place N=128
  14.451 ms (6 allocations: 32.13 MiB)
In-place N=128
  11.239 ms (0 allocations: 0 bytes)

Roger-luo · 2020-05-08T06:11:16Z

Hi, @elisno I'm glad this helps.

Regarding large Kronecker sum, I'm not sure about your problem, but if your matrices are more special (e.g general permutation matrices) we have an even faster implementation (as a summation of quantum circuits) in Yao. You may want to check that out.

MichielStock · 2020-05-08T08:39:15Z

I would love to use this in my Kronecker package, as suggested in this issue. What would be the best approach to having access to this before it is shipped in stable Julia version?

dkarrasch · 2020-05-08T08:42:41Z

I would love to use this in my Kronecker package, as suggested in this issue. What would be the best approach to having access to this before it is shipped in stable Julia version?

First this needs to get merged, and then maybe someone could add those parts of this PR that are a new feature to Compat.jl? Or you have a if VERSION < "" ... end block in which you define what you need.

Roger-luo · 2020-05-09T12:09:23Z

bump. I think I've fixed all the tests. any chance to review this? @dkarrasch

dkarrasch

This is very nice work, indeed, and consistent with all the *-mul! business. I only have a couple of non-functional comments, so this should be quick. I'd only like to ask not to overwrite the PR by push-forced, so that I can quickly review afterwards.

dkarrasch · 2020-05-09T14:56:34Z

base/operators.jl

+function kron! end
+


What is this good for? For stack traces or something like that?

this is because kron is defined in Base (https://github.com/JuliaLang/julia/blob/master/base/operators.jl#L533), thus all other stdlibs need to import Base: kron, I think it'd be more consistent to define the generic function in Base instead of LinearAlgebra.

stdlib/LinearAlgebra/src/dense.jl

stdlib/LinearAlgebra/src/diagonal.jl

stdlib/SparseArrays/src/linalg.jl

stdlib/LinearAlgebra/src/diagonal.jl

dkarrasch · 2020-05-10T13:52:35Z

Some changes in strings/search probably unintentionally leaked into this PR, and there is a trailing whitespace issue.

base/strings/search.jl

Roger-luo · 2020-05-11T21:06:57Z

the failed test seems unrelated to this PR. I've reverted the extra commit from this PR.

dkarrasch · 2020-05-12T08:53:47Z

Let's try again.

Roger-luo · 2020-05-12T11:28:18Z

looks good now, all tests pass.

dkarrasch

This LGTM. I think we may want to announce this in NEWS.md, since this is of interest to external packages as well. I'll leave this open for a few days, hoping that another pair of eyes takes a look.

dkarrasch · 2020-05-14T09:31:50Z

@Roger-luo Could you please add an item to NEWS.md, announcing this feature? You may need to rebase the PR, but I'm not sure, to add to the Julia v1.6 version of NEWS.md.

add sparse arrays

Roger-luo · 2020-05-14T19:06:44Z

The failed test on MacOS seems unrelated.

Roger-luo · 2020-05-18T20:49:17Z

bump, anyone else could take a look at this PR maybe?

martinholters requested changes Feb 14, 2019

View reviewed changes

stdlib/LinearAlgebra/src/dense.jl Outdated Show resolved Hide resolved

Roger-luo force-pushed the in-place-kron branch from b12ea18 to 62b45de Compare February 14, 2019 10:01

martinholters approved these changes Feb 14, 2019

View reviewed changes

martinholters reviewed Feb 14, 2019

View reviewed changes

ararslan requested changes Feb 14, 2019

View reviewed changes

ararslan requested a review from andreasnoack February 14, 2019 18:35

elisno mentioned this pull request May 6, 2020

Collect GeneralizedKroneckerProducts in-place MichielStock/Kronecker.jl#61

Open

Roger-luo changed the title ~~add in-place kron~~ WIP: add in-place kron May 6, 2020

dkarrasch reviewed May 6, 2020

View reviewed changes

stdlib/SparseArrays/src/linalg.jl Outdated Show resolved Hide resolved

dkarrasch requested changes May 6, 2020

View reviewed changes

Roger-luo force-pushed the in-place-kron branch from edb7d49 to 6f049a7 Compare May 7, 2020 10:59

Roger-luo changed the title ~~WIP: add in-place kron~~ add in-place kron May 7, 2020

Roger-luo requested review from dkarrasch and ararslan May 7, 2020 13:14

dkarrasch requested changes May 9, 2020

View reviewed changes

dkarrasch reviewed May 10, 2020

View reviewed changes

stdlib/LinearAlgebra/src/diagonal.jl Outdated Show resolved Hide resolved

dkarrasch reviewed May 11, 2020

View reviewed changes

base/strings/search.jl Outdated Show resolved Hide resolved

dkarrasch closed this May 12, 2020

dkarrasch reopened this May 12, 2020

dkarrasch approved these changes May 12, 2020

View reviewed changes

Roger-luo force-pushed the in-place-kron branch 4 times, most recently from 5d08842 to 23bfba1 Compare May 14, 2020 11:10

add inplace kron

de7413d

add sparse arrays

Roger-luo force-pushed the in-place-kron branch from 23bfba1 to de7413d Compare May 14, 2020 11:12

Roger-luo closed this May 14, 2020

Roger-luo reopened this May 14, 2020

dkarrasch merged commit cebd4fa into JuliaLang:master May 19, 2020

elisno mentioned this pull request May 29, 2020

Sparse Kronecker sums QuantumBFS/LuxurySparse.jl#49

Open

simeonschaub pushed a commit to simeonschaub/julia that referenced this pull request Aug 11, 2020

add inplace kron (JuliaLang#31069)

231ce2e

vtjnash mentioned this pull request Apr 13, 2021

kron between Diagonal and non-numeric matrix gives output with undefs #40439

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add in-place kron #31069

add in-place kron #31069

Roger-luo commented Feb 14, 2019

martinholters Feb 14, 2019

Roger-luo Feb 14, 2019

ararslan Feb 14, 2019

Roger-luo Feb 14, 2019

ararslan Feb 14, 2019

dkarrasch commented May 6, 2020

dkarrasch left a comment

dkarrasch May 6, 2020

dkarrasch May 6, 2020

dkarrasch May 6, 2020

Roger-luo commented May 7, 2020

elisno commented May 7, 2020

Roger-luo commented May 8, 2020 •

edited

Loading

MichielStock commented May 8, 2020

dkarrasch commented May 8, 2020

Roger-luo commented May 9, 2020

dkarrasch left a comment

dkarrasch May 9, 2020

Roger-luo May 9, 2020

dkarrasch commented May 10, 2020

Roger-luo commented May 11, 2020

dkarrasch commented May 12, 2020

Roger-luo commented May 12, 2020

dkarrasch left a comment

dkarrasch commented May 14, 2020

Roger-luo commented May 14, 2020

Roger-luo commented May 18, 2020

	@inline Base.@propagate_inbounds function kron!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
	Base.@propagate_inbounds function kron!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

add in-place kron #31069

add in-place kron #31069

Conversation

Roger-luo commented Feb 14, 2019

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 May 6, 2020

dkarrasch 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

Roger-luo commented May 7, 2020

elisno commented May 7, 2020

Dense

Sparse

Roger-luo commented May 8, 2020 • edited Loading

MichielStock commented May 8, 2020

dkarrasch commented May 8, 2020

Roger-luo commented May 9, 2020

dkarrasch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch commented May 10, 2020

Roger-luo commented May 11, 2020

dkarrasch commented May 12, 2020

Roger-luo commented May 12, 2020

dkarrasch left a comment

Choose a reason for hiding this comment

dkarrasch commented May 14, 2020

Roger-luo commented May 14, 2020

Roger-luo commented May 18, 2020

Roger-luo commented May 8, 2020 •

edited

Loading