Adds more Kronecker product specialisations

[zygmuntszpak]

Defines kron between two StaticVector types as well as between a StaticVector and a StaticMatrix type. The definitions also include cases for transposed StaticVector types (i.e. a StaticVector wrapped by a RowVectortype).

There are three broken tests which are due to the fact that taking the transpose of a SizedArray returns a SMatrix. I've made a comment in the relevant part of the test. Not sure of this behaviour of transpose is intentional or not?

[andyferris]

That's a good question. Our story with similar_type should probably be more nuanced... we might be able to avoid some pain by putting our "reasonable" size limits there.

[andyferris]

Dead code

[zygmuntszpak]

Ooops

[andyferris]

Probably more direct to reshape b into a matrix... but I don't think we have a method for that?

[zygmuntszpak]

Are you suggesting that we could reshape b into a 1 x 3 matrix and then dispatch on the kron that operates on two static matrices?

It seems that there is no method currently defined for that.
For example, with

s = @SVector [1, 2,  3];
reshape(s',(1,3))

1×3 Base.ReshapedArray{Int64,2,RowVector{Int64,StaticArrays.SArray{Tuple{3},Int64,1,3}},Tuple{}}:
 1  2  3

I get this Base.ReshapedArray type instead of the required StaticMatrix.

[coveralls]

Coverage increased (+0.2%) to 92.331% when pulling 12e4bef on zygmuntszpak:svec_kron into 3ee2ec2 on JuliaArrays:master.

[andyferris]

What's the purpose of the outermost array?

I'm curious about how good (performant) vcat is with four or more arguments.

[zygmuntszpak]

Is vcat generally slower than hcat? I was trying to imitate what was already done for the kron(StaticMatrix, StaticMatrix).

[zygmuntszpak]

Now that you pointed it out I realize that the outer array serves no purpose. Its an artifact of my attempt at molding the code that was already there for the Kronecker product between two matrices. I've removed it from other parts of the code as well.

[andyferris]

Regarding vcat, I mean we have a fast version of vcat(a, b) but not necessarily of vcat(a, b, c, d, e, f, g) (we should just check the generated code using @code_native looks relatively optimal over a variety of sizes).

The pattern elsewhere is mostly to construct the final array in one go (and use similar_type), rather than make temporary values and concatenate them together.

[codecov-io]

Codecov Report

Merging #407 into master will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #407   +/-   ##
=======================================
  Coverage   92.15%   92.15%           
=======================================
  Files          37       37           
  Lines        2843     2843           
=======================================
  Hits         2620     2620           
  Misses        223      223

---

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 3ee2ec2...12e4bef. Read the comment docs.

[zygmuntszpak]

I've removed the dead code and the superfluous outermost array. Would you like me to try to avoid vcat and use hcat instead?

[andyferris]

If were to do something about vcat I would consider just constructing the entire output array using similar_type and a tuple with all the computed elements, like done in many places throughout the codebase.

[zygmuntszpak]

Hi Andy,

I've been trying to figure out what you mean but to no avail. It is safe to assume I know nothing about meta programming, and am still learning aspects of Julia.

The only way that I can see how one can construct the entire output array of the Kronecker product that will avoid the use vcat is with a comprehension.
For example,

@inline kron(a::StaticMatrix, b::StaticMatrix) = _kron(_length_limit, Size(a), Size(b), a, b)
@generated function _kron(::Length{length_limit}, ::Size{SA}, ::Size{SB}, a, b) where {length_limit,SA,SB}
    outsize = SA .* SB
    if prod(outsize) > length_limit
        return :( SizedMatrix{$(outsize[1]),$(outsize[2])}( kron(drop_sdims(a), drop_sdims(b)) ) )
    end
    M = [ a[ia,ja] * b[ib,jb] for ia in 1:SA[1],  ja in 1:SA[2],  ib in 1:SB[1],  jb in 1:SB[2]]
    return quote
        @_inline_meta
        @inbounds return similar_type($(M),Size($(outsize)))($(M))  
    end
end

However, this give the error

ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type StaticArrays.SArray{Tuple{2,2},Int64,2,4}
This may have arisen from a call to the constructor StaticArrays.SArray{Tuple{2,2},Int64,2,4}(...),
since type constructors fall back to convert methods.

with the following input

A = @SMatrix [1 2; 3 4]
kron(A,A)

On the other hand, if I test the code in my own script (outside of the generated function) it seems to work:

M = [ A[ia,ja] * A[ib,jb] for ia in 1:2,  ja in 1:2,  ib in 1:2,  jb in 1:2]
N = similar_type(M,Size((4,4)))(M)
4×4 StaticArrays.SArray{Tuple{4,4},Int64,2,16}:
 1   3  2   4
 3   9  6  12
 2   6  4   8
 4  12  8  16

I also don't understand how precisely one avoids allocations. If I use a comprehension then will it not be allocating on the heap automatically? On the other hand, the existing code also uses a comprehension of rows which seems not to allocate... I guess I don't understand the design principle behind what makes this package work.

I'd appreciate if you could give me a more detailed example of how you would like me to modify the current implementation. Also, what is the reason for the error in the generated function but not in my own script?

Best wishes,

Zygmunt

[andyferris]

Apologies for brevity, but try something like @inbounds return similar_type($a, Size($(outsize)))(tuple($(M...))) (however I suspect you need to swap the comprehension order ja and ib to get the elements arranged the right way?)

[zygmuntszpak]

Yes, I need to rethink the ordering and reshaping of the comprehension but unfortunately I still get an error even after I make your suggested amendment.

ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type StaticArrays.SArray{Tuple{2,2},Int64,2,4}
This may have arisen from a call to the constructor StaticArrays.SArray{Tuple{2,2},Int64,2,4}(...),
since type constructors fall back to convert methods.
Stacktrace:
 [1] getindex(::Type{StaticArrays.SArray{Tuple{2,2},Int64,2,4}}, ::Int64, ::Int64) at ./array.jl:208
 [2] (::StaticArrays.##285#286{DataType,DataType})(::NTuple{4,Int64}) at ./<missing>:0
 [3] collect(::Base.Generator{Base.Iterators.Prod{UnitRange{Int64},Base.Iterators.Prod{UnitRange{Int64},Base.Iterators.Prod2{UnitRange{Int64},UnitRange{Int64}}}},StaticArrays.##285#286{DataType,DataType}}) at ./array.jl:475
 [4] _kron(...) at /home/zygmunt/.julia/v0.6/StaticArrays/src/linalg.jl:364
 [5] kron(::StaticArrays.SArray{Tuple{2,2},Int64,2,4}, ::StaticArrays.SArray{Tuple{2,2},Int64,2,4}) at /home/zygmunt/.julia/v0.6/StaticArrays/src/linalg.jl:355
 [6] macro expansion at /home/zygmunt/.julia/v0.6/Atom/src/repl.jl:118 [inlined]
 [7] anonymous at ./<missing>:?

That's with the following return

    return quote
        @_inline_meta
        @inbounds return similar_type($a, Size($(outsize)))(tuple($(M...))) 
    end

The stack trace points to this line (with the current incorrectly ordered comprehension)

M = [ a[ia,ja] * b[ib,jb] for ia in 1:SA[1],  ja in 1:SA[2],  ib in 1:SB[1],  jb in 1:SB[2]]

[zygmuntszpak]

The correct order of the comprehension is as follows:

M = [ A[ia,ja] * B[ib,jb] for ib in 1:SB[1],ia in 1:SA[1],jb in 1:SB[2],ja in 1:SA[2]]

if one wants to use reshape(M,outsize).

A quick test...

A = @SMatrix [1 2; 3 4]
B = @SMatrix [5 6; 7 8]
P = [1 2; 3 4]
Q = [5 6; 7 8]
M = [ A[ia,ja] * B[ib,jb] for ib in 1:SB[1],ia in 1:SA[1],jb in 1:SB[2],ja in 1:SA[2]]
N = similar_type(M,Size((4,4)))(M)
@test N == kron(P,Q)
Test Passed

A = @SMatrix [5 6; 7 8]
B = @SMatrix [1 2; 3 4]
M = [ A[ia,ja] * B[ib,jb] for ib in 1:SB[1],ia in 1:SA[1],jb in 1:SB[2],ja in 1:SA[2]]
N = similar_type(M,Size((4,4)))(M)
@test N == kron(Q,P)
Test Passed

[zygmuntszpak]

I've posted a question on the Discourse forum to see if anyone has any suggestions for resolving the error.

[zygmuntszpak]

As per your request, for each specialization the entire output array is constructed using similar_type and a tuple with all the computed elements. In so doing, I no longer return the transpose of anything and hence the previous test_broken cases are no longer an issue.

[andyferris]

Great, thanks. Love the tests, and the codegen style matches the rest of the package.

While I mentioned that in the future I'd love to consider moving some logic which chooses between SizedArray or not to similar_type, it's still important to gate this here because just like for matrix multiplication, the generated code size grows faster than linear in the input sizes. (In fact, taking that to the extreme, the output size also grows too fast too... at some point it's better just to drop the SizedArray and use Array).