add Kronecker products for LinearMaps

[dkarrasch]

I realized that part of @JeffFessler's motivation for having LinearMaps <: AbstractArray is to use them in combination with Kronecker.jl. I thought I'll lower the pressure and provide a Kronecker product for LinearMaps without the need for getindex. 😛

The implementation uses the well-known vec trick, and since the Kronecker product is associative, one only needs to take care of the product of two linear maps. I couldn't find any trick for "higher-order" Kronecker products. Note that this is all lazy, and all linear map types are, as usual, inferred. This works out of the box with wrapped maps of dense/sparse/structured matrices, or function maps, etc.

[coveralls]

Coverage decreased (-0.3%) to 94.574% when pulling fd10ea2 on kronecker into 95cc89c on master.

[codecov]

Codecov Report

Merging #61 into master will decrease coverage by 0.58%.
The diff coverage is 89.28%.

@@            Coverage Diff             @@
##           master      #61      +/-   ##
==========================================
- Coverage   98.83%   98.24%   -0.59%     
==========================================
  Files           8        9       +1     
  Lines         428      456      +28     
==========================================
+ Hits          423      448      +25     
- Misses          5        8       +3

Impacted Files	Coverage Δ
src/LinearMaps.jl	100% <ø> (ø)	⬆️
src/kronecker.jl	89.28% <89.28%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a024a4d...641a770. Read the comment docs.

[codecov]

Codecov Report

Merging #61 into master will decrease coverage by 1.06%.
The diff coverage is 94.62%.

@@            Coverage Diff             @@
##           master      #61      +/-   ##
==========================================
- Coverage   97.49%   96.42%   -1.07%     
==========================================
  Files           8       10       +2     
  Lines         479      616     +137     
==========================================
+ Hits          467      594     +127     
- Misses         12       22      +10

Impacted Files	Coverage Δ
src/wrappedmap.jl	100% <100%> (ø)	⬆️
src/composition.jl	98.41% <100%> (+0.02%)	⬆️
src/uniformscalingmap.jl	96.66% <100%> (+0.3%)	⬆️
src/LinearMaps.jl	87.27% <57.14%> (-11.34%)	⬇️
src/conversion.jl	95.34% <95.34%> (ø)
src/kronecker.jl	98.05% <98.05%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 95cc89c...fd10ea2. Read the comment docs.

[JeffFessler]

Code looks fine to me.
One could easily add a method for the product of two Kronecker maps,
returning yet another Kron. map.

[dkarrasch]

One could easily add a method for the product of two Kronecker maps,
returning yet another Kron. map.

I saw that in Wikipedia, but I thought that works only when the matrix products exists, but they don't have to...?

[dkarrasch]

Hooray, I managed to add promotion of matrices (as long as there is one LinearMap among the first 8 arguments of kron. That took me a great while, to figure out the metaprogramming issues. I'm not sure I understood what turned out to be correct, but anyway.

[JeffFessler]

What I have in mind would inherit the dimension checks already there in making a CompositeMap namely
https://github.com/Jutho/LinearMaps.jl/blob/master/src/composition.jl#L74

From wikipedia, the "mixed-product property:"
(A x B) * (C x D) = (AC x BD) as long as A*C is valid and B*D is valid.

Without thinking about types (as usual for me :) ) I think it is almost a one liner (lazy rough sketch; should use if dim checks instead of try):

*(A <: KroneckerMap, B <: KroneckerMap) =
begin
 try
  kron(A.maps[1]*B.maps[1], A.maps[2]*B.maps[2])
 catch
  CompositeMap{T}(tuple(A, B))
end

Ah, but type instability... Again my use case for this is in the preamble before iterating so I don't think that instability would matter.
If you don't want to overload * then perhaps a separate kron_mul method or such.

[JeffFessler]

Impressive metaprogramming but it is hard for me to follow. A couple comments would really help!
But I really don't understand how this works with the multiply with vector.
The vector multiply as just A,B = L.maps so what happens with 3 or more maps?
I don't know how the vec trick generalizes to beyond 2, though perhaps it can with some type of recursion. Again, maybe a comment in the code to explain would help?

[dkarrasch]

In terms of numerical effort, I assume it would be beneficial to use that mixed-property? Could be that this would be a second "good" use case for the union splitting mechanism.

The promotion thing is really "simple". All it does is to take the matrix arguments and wraps them in a LinearMap. Finally, it simply calls kron again on the same list (up to to lmap-promotion). Ultimately, how is kron(A, B, C) handled? Well, by Julia Base. It iteratively calls kron(kron(A, B), C), but this time, everytime kronis called, it is called with only two arguments, each aLinearMap, so the method from two lines above gets called. The vec trick applies consecutively, first to (A x B) x C(whateverA x Bis at that point), so it'stranspose(C)X(A x B), and when you multiply A x Bwith somex`, you again apply the vec trick.

I'll add some comments to the code.

[dkarrasch]

I just did a quick check, and in fact it's beneficial to use the mixed-trick. I'll put it in the PR, then we discuss the stability issue.

[JeffFessler]

looks great!

[dkarrasch]

I thought I'll post a little comparison here for general entertainment. I'm sure I have used LazyArrays.jl incorrectly and I'm happy to fix the code below, but here is what I tested:

using BenchmarkTools, Test
using LinearAlgebra, SparseArrays
using LinearMaps, Kronecker, LazyArrays
n = 100
x = rand(n*n)
y = zeros(n*n)

@info "LinearMaps.jl"
J = LinearMaps.UniformScalingMap(true, n)
E = spdiagm(-1 => trues(n-1))
D = E + E' - 2I
A1 = kron(D, J) + kron(J, D)
@btime $A1 * $x

@info "Kronecker.jl"
J = I(n)
A2 = kronecker(D, J)
A2´ = kronecker(J, D)
@btime $A2 * $x + $A2´ * $x

@info "LazyArrays.jl"
A3 = Kron(D, J)
A3´ = Kron(J, D)
@btime $A3 * $x + $A3´ * $x

@info "Base"
A4 = kron(D, J) + kron(J, D)
@btime $A4 * $x

@test sparse(A1) == A4

and here is what I got:

[ Info: LinearMaps.jl
  956.586 μs (726 allocations: 521.30 KiB)
[ Info: Kronecker.jl
  5.115 ms (64 allocations: 1.35 MiB)
[ Info: LazyArrays.jl
  5.771 s (6 allocations: 234.61 KiB)
[ Info: Base
  55.835 μs (2 allocations: 78.20 KiB)

This is the discrete Laplace operator. Note that with Kronecker.jl and LazyArrays.jl the multiplication needs to be done in two steps: if you add the two Kronecker products, you'll get a dense matrix, and, of course, performance is horrible. Base is so fast because the Kronecker product of sparse matrices is sparse again, but Kronecker.jl has shown that for dense matrices performance is much worse compared to kronecker.

One thing I noted is that there doesn't seem to be a "lazy" version of a UniformScaling with known size. It seems like you need to construct a sparse or diagonal matrix, with obvious redundant memory consumption. @Jutho So we may wish to consider exporting our UniformScalingMap, for instance, for constructions like the above.

[JeffFessler]

Nice timing results! I didn't know about LazyArrays.jl before.
On the plus side, LazyArrays has getindex and setindex and subtypes of AbstractArray :-)
The utf8 in A2´ looks so much like A2' that it confused me at first btw!

[dkarrasch]

Yeah, I didn't know which other symbol to use so I typed some random symbol. 😛 Well, and yes, the array type behavior is what distinguishes those packages from ours here.

Some day you'll have to explain me why you can't live without getindex and the AbstractArray subtype. Until that day, I'll try to provide every other feature that you need here. I'm thinking about adding a BlockDiagonalMap type, the corresponding thing to what I saw in MIRT.jl.

[dkarrasch]

Regarding the minimal storage but known-size uniform scaling object, I just found this discussion here: JuliaLang/julia#30298

The use of I(n) in Kronecker products seems to be indeed a very classic case. But I just realized that since UniformScalings are not AbstractMatrixes, we could safely define

LinearMap(J::UniformScaling, n::Int) = UniformScalingMap(J.λ, n)

and then we don't need to export UniformScalingMaps.

[JeffFessler]

Thumbs up for Kronecker sum!

[dkarrasch]

Thanks, and mat-vec-multiplication with Kronecker products is also faster now. I reuse intermediate vectors (called temp1 and temp2) when applying the vec trick. The Kronecker sum requires a bit of work on conversion of operators to matrices. We should add a couple of specialized conversion rules, especially when operators are made up of matrices. But with that, I'm going to wait until #66 and #68 are finished and merged.

[dkarrasch]

This is still WIP, basically it needs some more tests, but I wanted to quickly share some impressive results, obtained on Julia v1.3 with the 5-arg mul!.

using BenchmarkTools, Test
using LinearAlgebra, SparseArrays
using LinearMaps, Kronecker#, LazyArrays
n = 100
x = rand(n*n)
y = zeros(n*n)

@info "LinearMaps.jl"
J = LinearMap(I, n);
E = spdiagm(-1 => trues(n-1));
D = E + E' - 2I;
L = kron(D, J) + kron(J, D);
@btime $L * $x; # 89.378 μs (27 allocations: 313.48 KiB)

@info "Kronecker.jl"
J = Kronecker.Eye(n);
A2 = kronecker(D, J);
A2´ = kronecker(J, D);
@btime ($y = $A2 * $x; $y.+= $A2´ * $x); # 89.393 μs (17 allocations: 313.20 KiB)

@info "Base"
J = I(n);
A4 = kron(D, J) + kron(J, D);
@btime $A4 * $x; # 57.114 μs (2 allocations: 78.20 KiB)

We are getting really close to the super-sparse matrix multiplication. Interestingly, the operator is nothing but a Kronecker sum.

J = LinearMap(I, n);
E = spdiagm(-1 => trues(n-1)); D = E + E' - 2I; DL = LinearMap(D);
Δ₂ = LinearMaps.kronsum(DL, DL);
Δ₃ = kroneckersum(D, D);
@btime $Δ₂ * $x; # 67.537 μs (7 allocations: 78.41 KiB)
@btime $Δ₃ * $x; # 67.327 μs (7 allocations: 78.41 KiB)

So, we are able to fully leverage all LinearAlgebra functionality in the case of matrix-based linear maps!

[dkarrasch]

This has grown big, as usual, but I think it's ready for review. For Kronecker products of matrix-based linear maps to work fast, there were quite some adjustments needed for MatrixMaps and UniformScalings. One interesting aspect is that one sometimes cannot avoid intermediate matrix results. Obviously, when these are computed from matrices, you don't want to generate these by applying them to the canonical basis vectors. As demonstrated above, there is also quite some benefit from the 5-arg mul! implementation. I'd suggest to finalize and merge #66 before this one, and afterwards make a minor release.

[dkarrasch]

Brief update: I made the Kronecker (sum) power a type parameter, so this way everything is inferable. This puts a little overhead onto the compiler, but I guess this is not going to be used with many different powers.

[Jutho]

Hi @dkarrasch , could you remind me again which open PRs are finalised and ready for review; I keep forgetting and was very busy the last weeks.

[dkarrasch]

It would be good to get this one done here.

[dkarrasch]

Sorry for messing around. When I'm bored, I'm starting to do micro-optimization. 😂 I've now optimized dimension checks, and included the option to skip them (as usual) via @inbounds. These checks always refer to whether the lengths of y and x are consistent with that of A (but not to one-based indexing). In many cases, this is guaranteed upfront, and there is no need to re-check this everytime. Accordingly, we may then wish to propagate that knowledge. I hope I have carefully distinguished when we need Base.@propagate_inbounds and when @inline.

Okay, I promise I'll leave this PR in peace for review. 🤞 But actually, the recent changes are just standard refactors of the previous code, no design changes or anything like that.

[dkarrasch]

gentle bump

[Jutho]

I tried to complete the review now but I am too tired; this time I will try not to forget.

[Jutho]

Generally looks good; I think this can be merged, maybe some small comments for consideration.