Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Avoid explicitly evaluating zero elements in certain structured matrix broadcast operations #985

Closed
jishnub opened this issue Feb 13, 2023 · 3 comments · Fixed by JuliaLang/julia#53909
Closed

Avoid explicitly evaluating zero elements in certain structured matrix broadcast operations #985

jishnub opened this issue Feb 13, 2023 · 3 comments · Fixed by JuliaLang/julia#53909

Comments

@jishnub
Copy link
Collaborator

jishnub commented Feb 13, 2023

For block-banded structured matrix types, the zero elements may not be well-defined, e.g:

julia> D = Diagonal([Float64[1 2; 3 4], Float64[5 6; 7 8]])
2×2 Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}:
 [1.0 2.0; 3.0 4.0]          ⋅         
         ⋅           [5.0 6.0; 7.0 8.0]

julia> D .+ D
ERROR: MethodError: no method matching zero(::Type{Matrix{Float64}})

Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period
   @ Dates ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/Dates/src/periods.jl:51
  zero(::AbstractIrrational)
   @ Base irrationals.jl:151
  zero(::Diagonal)
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/special.jl:324
  ...

Stacktrace:
 [1] fzero(S::Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}})
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/structuredbroadcast.jl:146
 [2] map
   @ ./tuple.jl:290 [inlined]
 [3] fzero(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, typeof(+), Tuple{Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}, Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}}})
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/structuredbroadcast.jl:149
 [4] fzeropreserving(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, typeof(+), Tuple{Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}, Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}}})
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/structuredbroadcast.jl:137
 [5] similar(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, typeof(+), Tuple{Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}, Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}}}, #unused#::Type{Matrix{Float64}})
   @ LinearAlgebra ~/packages/julias/julia-latest/share/julia/stdlib/v1.10/LinearAlgebra/src/structuredbroadcast.jl:155
 [6] copy(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, typeof(+), Tuple{Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}, Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}}})
   @ Base.Broadcast ./broadcast.jl:912
 [7] materialize(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Diagonal}, Nothing, typeof(+), Tuple{Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}, Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}}})
   @ Base.Broadcast ./broadcast.jl:887
 [8] top-level scope
   @ REPL[3]:1

However, in this case, the result may be obtained without any reference to the zeros.

julia> D + D
2×2 Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}:
 [2.0 4.0; 6.0 8.0]          ⋅         
         ⋅           [10.0 12.0; 14.0 16.0]

I wonder if it might be possible to evaluate the result using broadcasting without explicitly evaluating the zero elements?
@aravindh-krishnamoorthy
Copy link
Contributor

Hello, could you please advise on why this may be advantageous?

@jishnub
Copy link
Collaborator Author

jishnub commented May 3, 2023

This is more about getting edge cases working so that there are no surprises when writing generic code. The specific example presented above may not mean much

@jishnub
Copy link
Collaborator Author

jishnub commented Mar 31, 2024

Actually, turns out that this is needed for BlockArrays.jl

This was referenced Mar 31, 2024
Merged
Closed
jishnub referenced this issue in JuliaLang/julia Apr 7, 2024
@jishnub
Support broadcasting over structured block matrices (#53909)
243ebc3 
Fix https://github.com/JuliaLang/julia/issues/48664

After this, broadcasting over structured block matrices with
matrix-valued elements works:
```julia
julia> D = Diagonal([[1 2; 3 4], [5 6; 7 8]])
2×2 Diagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [1 2; 3 4]      ⋅     
     ⋅       [5 6; 7 8]

julia> D .+ D
2×2 Diagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [2 4; 6 8]      ⋅     
     ⋅       [10 12; 14 16]

julia> cos.(D)
2×2 Matrix{Matrix{Float64}}:
 [0.855423 -0.110876; -0.166315 0.689109]  [1.0 0.0; 0.0 1.0]
 [1.0 0.0; 0.0 1.0]                        [0.928384 -0.069963; -0.0816235 0.893403]
```
Such operations show up when using `BlockArrays`.

The implementation is a bit hacky as it uses `0I` as the zero element in
`fzero`, which isn't really the correct zero if the blocks are
rectangular. Nonetheless, this works, as `fzero` is only used to
determine if the structure is preserved.
KristofferC referenced this issue Nov 14, 2024
@jishnub
Support broadcasting over structured block matrices (#53909)
d70c149 
Fix https://github.com/JuliaLang/julia/issues/48664

After this, broadcasting over structured block matrices with
matrix-valued elements works:
```julia
julia> D = Diagonal([[1 2; 3 4], [5 6; 7 8]])
2×2 Diagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [1 2; 3 4]      ⋅     
     ⋅       [5 6; 7 8]

julia> D .+ D
2×2 Diagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [2 4; 6 8]      ⋅     
     ⋅       [10 12; 14 16]

julia> cos.(D)
2×2 Matrix{Matrix{Float64}}:
 [0.855423 -0.110876; -0.166315 0.689109]  [1.0 0.0; 0.0 1.0]
 [1.0 0.0; 0.0 1.0]                        [0.928384 -0.069963; -0.0816235 0.893403]
```
Such operations show up when using `BlockArrays`.

The implementation is a bit hacky as it uses `0I` as the zero element in
`fzero`, which isn't really the correct zero if the blocks are
rectangular. Nonetheless, this works, as `fzero` is only used to
determine if the structure is preserved.
@KristofferC KristofferC transferred this issue from JuliaLang/julia Nov 26, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Support broadcasting over structured block matrices JuliaLang/julia
2 participants
@jishnub @aravindh-krishnamoorthy