cat SMatrices returns Arrays

[rvignolo]

Hi!

I was trying to compute a static block matrix using static sub matrices and saw that the result is not a StaticArray.

julia> m1 = @SMatrix rand(2,3);

julia> m2 = @SMatrix rand(2,3);

julia> cat(m1, m2, dims = (1, 2))
4×6 Array{Float64,2}:
 0.611586  0.396562  0.30208    0.0       0.0       0.0
 0.574913  0.909981  0.0731041  0.0       0.0       0.0
 0.0       0.0       0.0        0.495477  0.512808  0.232148
 0.0       0.0       0.0        0.880756  0.449617  0.305318

Also, trying simpler examples yield to Arrays as well:

julia> cat(m1, m2, dims = 1)
4×3 Array{Float64,2}:
 0.611586  0.396562  0.30208
 0.574913  0.909981  0.0731041
 0.495477  0.512808  0.232148
 0.880756  0.449617  0.305318

julia> cat(m1, m2, dims = 2)
2×6 Array{Float64,2}:
 0.611586  0.396562  0.30208    0.495477  0.512808  0.232148
 0.574913  0.909981  0.0731041  0.880756  0.449617  0.305318

The result is different if we just use hcat or vcat.

julia> hcat(m1, m2)
2×6 SArray{Tuple{2,6},Float64,2,12} with indices SOneTo(2)×SOneTo(6):
 0.611586  0.396562  0.30208    0.495477  0.512808  0.232148
 0.574913  0.909981  0.0731041  0.880756  0.449617  0.305318

julia> vcat(m1, m2)
4×3 SArray{Tuple{4,3},Float64,2,12} with indices SOneTo(4)×SOneTo(3):
 0.611586  0.396562  0.30208
 0.574913  0.909981  0.0731041
 0.495477  0.512808  0.232148
 0.880756  0.449617  0.305318

However, I still need to find a method that computes a block matrix using submatrices.

I have a function that computes several static matrices in the body that need to be concatenated as a block static matrix (which is then returned). Should I create static sub-matrices filled with zeros and use hcat and vcat then?

Also, in the slack channel, @mcabbott suggested using generated functions:

I guess this might end up more efficient. I think it ought to be possible to re-use cat by making arrays of expressions like :(args[n][i,j]) and cat-ing those… The one thing I contributed to StaticArrays is like that: https://github.com/JuliaArrays/StaticArrays.jl/blob/0a7f3f7eae1fab11a50699fe33332d5099534be9/src/abstractarray.jl#L218

Please, notice that the sub matrices might have different sizes as well.

Any ideas?

Thanks!

[andyferris]

The theory is that Julia linear algebra works also with nested arrays, so you can have a static array of static arrays of numbers. It’s similar to how a mathematician might represent a block matrix on paper/blackboard/whiteboard. In the case that your matrix has a block diagonal structure, the outer array can be diagonal so that we don’t need to store all those zeros and making computation faster.

(Theoretically the structure should also make algorithms for solving linear algebra problems, eigenvalues, etc faster but in practice I’m not sure how much of that currently works).

Does that make sense?

[rvignolo]

The theory is that Julia linear algebra works also with nested arrays, so you can have a static array of static arrays of numbers. It’s similar to how a mathematician might represent a block matrix on paper/blackboard/whiteboard.

Do you mean something like:

julia> m1 = @SMatrix rand(2,3);

julia> m2 = @SMatrix rand(3,3);

julia> v = SVector(m1, m2)
2-element SArray{Tuple{2},SArray{S,Float64,2,L} where L where S<:Tuple,1,2} with indices SOneTo(2):
 [0.8836511927390833 0.7779904745859783 0.12003972967976773; 0.9281821132469268 0.10058434721711107 0.2922066338953775]
 [0.620506339300614 0.1835853176439728 0.12548038761092162; 0.07619298087618231 0.6922824755394659 0.12951291742159943; 0.03283194019384905 0.33477716333092733 0.19351822469878832]

Working with this would be great but the output of this function, i.e. the block matrix, is not used in my solvers, it is used in StochasticDiffEq.jl solvers. I would have to discuss if this is possible.

Please, let me know if I have understand your suggestion correctly. Thank you!

[andyferris]

OK, what I had in mind was this:

julia> using StaticArrays

julia> m1 = rand(SMatrix{2,2})
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
 0.403896  0.410033
 0.262563  0.15044

julia> m2 = rand(SMatrix{2,2})
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
 0.476859  0.567338
 0.723458  0.634039

julia> m = SDiagonal(m1, m2)
2×2 LinearAlgebra.Diagonal{SArray{Tuple{2,2},Float64,2,4},SArray{Tuple{2},SArray{Tuple{2,2},Float64,2,4},1,2}}:
 [0.403896 0.410033; 0.262563 0.15044]          ⋅         
         ⋅                              [0.476859 0.567338; 0.723458 0.634039]

m is an a linear operator of a vector space which happens to be block diagonal. Here's a proof - we can create a corresponding block vector v and do all the usual linear algebra operations on it:

julia> v = SVector(rand(SVector{2}), rand(SVector{2}))
2-element SArray{Tuple{2},SArray{Tuple{2},Float64,1,2},1,2} with indices SOneTo(2):
 [0.5771816379090431, 0.5291493155287397]
 [0.057315507065596405, 0.6548649874117727]

julia> v + v
2-element SArray{Tuple{2},SArray{Tuple{2},Float64,1,2},1,2} with indices SOneTo(2):
 [1.1543632758180862, 1.0582986310574793]
 [0.11463101413119281, 1.3097299748235454]

julia> 2.0 * v
2-element SArray{Tuple{2},SArray{Tuple{2},Float64,1,2},1,2} with indices SOneTo(2):
 [1.1543632758180862, 1.0582986310574793]
 [0.11463101413119281, 1.3097299748235454]

julia> m * v
2-element SArray{Tuple{2},SArray{Tuple{2},Float64,1,2},1,2} with indices SOneTo(2):
 [0.3801852387122504, 0.345769483620375]
 [0.14995451067455998, 0.4059756370380156]

So m is a linear operator. Like all linear operators, it lives in its own vector space:

julia> m + m
2×2 LinearAlgebra.Diagonal{SArray{Tuple{2,2},Float64,2,4},SArray{Tuple{2},SArray{Tuple{2,2},Float64,2,4},1,2}}:
 [0.807792 0.820065; 0.525125 0.30088]          ⋅         
         ⋅                              [0.953718 1.13468; 1.44692 1.26808]

julia> 2.0 * m
2×2 LinearAlgebra.Diagonal{SArray{Tuple{2,2},Float64,2,4},SArray{Tuple{2},SArray{Tuple{2,2},Float64,2,4},1,2}}:
 [0.522284 0.867275; 0.264996 1.01784]          ⋅         
         ⋅                              [0.202188 0.440275; 1.20041 1.13481]

Theoretically, we can do all the usual things with linear operator m like solving m \ v, squaring m, inverting m, finding its eigenvalues, etc. I think this is super cool:

julia> m \ v
2-element SArray{Tuple{2},SArray{Tuple{2},Float64,1,2},1,2} with indices SOneTo(2):
 [0.8520191806581197, 0.8179269896030508]
 [1.493174936346978, -0.42535125674651797]

It's even the right answer up to floating-point rounding error:

julia> m * (m \ v) - v
2-element SArray{Tuple{2},SArray{Tuple{2},Float64,1,2},1,2} with indices SOneTo(2):
 [0.0, 0.0]
 [2.7755575615628914e-17, 0.0]

However - the idea here isn't necessarily fully baked (yet). Take this example:

julia> m * m
Unreachable reached at 0x7f3c8cba0775

... and it crashes. Lol - @c42f, maybe we should fix that one?

@rvignolo I have no idea whether StochasticDiffEq supports nested arrays, but if you are just doing * and + and it uses randn internally, it could potentially be fine. You could always try reach out on discourse, slack or zulip to those guys to see what they think. Also - using static arrays and block diagonal structures are optimization techniques - it's always worth checking the performance of the simpler flat Array problem first and optimizing from there.

[rvignolo]

@andyferris, thank you for your kind and elaborated response! Indeed, the SDiagonal hint is really cool!

However, it seems that SDiagonal works only if the all the sub matrices have the same size:

julia> m1 = @SMatrix rand(2,3);

julia> m2 = @SMatrix rand(4,2);

julia> SDiagonal(m1, m2)
2×2 LinearAlgebra.Diagonal{SArray{S,Float64,2,L} where L where S<:Tuple,SArray{Tuple{2},SArray{S,Float64,2,L} where L where S<:Tuple,1,2}}:
Error showing value of type LinearAlgebra.Diagonal{SArray{S,Float64,2,L} where L where S<:Tuple,SArray{Tuple{2},SArray{S,Float64,2,L} where L where S<:Tuple,1,2}}:
ERROR: The size of type `SArray{S,Float64,2,L} where L where S<:Tuple` is not known.

Is this a bug?

In my case, the block matrix is formed by sub matrices that do not necessarily have equal sizes, so I would need this feature to work with your suggestions.

Any ideas?

[andyferris]

No it’s not a bug but it is a known limitation (the elements are assumed to have the same type). Have you checked out the performance using Julia’s built in Array and Diagonal? It might be plenty fast for your problem.

[mateuszbaran]

Going back to the original issue, one small problem with these cat methods is that constant propagation of keyword arguments is needed to make this type stable and it's a Julia 1.6 thing. So on Julia <1.6 we have to choose between returning Array and being type unstable.

julia> m * m
Unreachable reached at 0x7f3c8cba0775

... and it crashes. Lol - @c42f, maybe we should fix that one?

At least it works on Julia master.

[rvignolo]

No it’s not a bug but it is a known limitation (the elements are assumed to have the same type).

Oh, okay. I was hoping it was because of the error message 😄 .

Have you checked out the performance using Julia’s built in Array and Diagonal? It might be plenty fast for your problem.

I am writing a macro that returns the diffusion function g of a Stochastic Differential Equation in both versions:

1. in place (IIP) + Arrays: the function g(du, u, p, t) receives a Matrix du that is modified in place using the Array u, parameters p and current time t.
2. out of pace (OOP) + SArrays: the function g(u, p, t) returns a SMatrix that is computed using the SArray u, parameters p and current time t.

The IIP version might be plenty fast but I have seen that each subproblem (i.e. each problem defined by one block matrix) is WAY faster when I use SArrays. I was hopping that I could define the OOP function version with SArrays as well, which means that I have to return a block matrix.

Going back to the original issue, one small problem with these cat methods is that constant propagation of keyword arguments is needed to make this type stable and it's a Julia 1.6 thing.

Hi @mateuszbaran, thank you for your kind response. This means that in Julia 1.6 you are going to choose returning a SArray (and being type stable)?

At this point I would like to say that I would really like to be able to use Static Arrays for my problems. What approach should I take 😆 ?

[mateuszbaran]

Hi @mateuszbaran, thank you for your kind response. This means that in Julia 1.6 you are going to choose returning a SArray (and being type stable)?

I think it would be reasonable but there is no decision yet.

At this point I would like to say that I would really like to be able to use Static Arrays for my problems. What approach should I take ?

Depending on the size of your problem, making your own heterogeneous diagonal matrix type that stores diagonal elements in a tuple may be a reasonable solution.

[manuelbb-upb]

I was caught a bit of guard today when I tried concatenating matrices vertically and horizontally at the same time, e.g.,

a = @SMatrix(rand(2,2))
[ a ; a ]   # isa SMatrix 
[ a a ]     # isa SMatrix 

# however this is not an SMatrix, but a Matrix{Float64}
[ a a; 
  a a]

I guess its the same issue? Using [email protected] and [email protected].

[Wu-Chenyang]

I encountered the same issue when I am trying to concatenate SMatrix with SVector.

a = SVector{3}(1,2,3)
M = [a a a] # isa SMatrix

julia> [M a; a' 0]
4×4 Matrix{Int64}:
 1  1  1  1
 2  2  2  2
 3  3  3  3
 1  2  3  0

It is not hard to fix though. Just concatenate rows and columns separately.

julia> [[M a]; [a' 0]]
4×4 SMatrix{4, 4, Int64, 16} with indices SOneTo(4)×SOneTo(4):
 1  1  1  1
 2  2  2  2
 3  3  3  3
 1  2  3  0

Only that you need to be highly conscious of this issue.

[mcabbott]

The original issue was specifically about cat & block-diagonal matrices (and keyword arguments, and blockwise linear algebra).

The last two posts are about hvcat which is what's used for block-matrix syntax. It takes no keywords, but does take an integer describing the shape. Probably #1262 is the best issue.