-
Notifications
You must be signed in to change notification settings - Fork 149
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
cat SMatrices returns Arrays #856
Comments
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? |
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! |
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]
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 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 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 |
@andyferris, thank you for your kind and elaborated response! Indeed, the However, it seems that 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? |
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 |
Going back to the original issue, one small problem with these
At least it works on Julia master. |
Oh, okay. I was hoping it was because of the error message 😄 .
I am writing a macro that returns the diffusion function
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
Hi @mateuszbaran, thank you for your kind response. This means that in Julia 1.6 you are going to choose returning a 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 😆 ? |
I think it would be reasonable but there is no decision yet.
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. |
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? |
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. |
The original issue was specifically about The last two posts are about |
Hi!
I was trying to compute a static block matrix using static sub matrices and saw that the result is not a
StaticArray
.Also, trying simpler examples yield to
Array
s as well:The result is different if we just use
hcat
orvcat
.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
andvcat
then?Also, in the slack channel, @mcabbott suggested using generated functions:
Please, notice that the sub matrices might have different sizes as well.
Any ideas?
Thanks!
The text was updated successfully, but these errors were encountered: