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 #48664

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

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

jishnub opened this issue Feb 13, 2023 · 3 comments · Fixed by #53909
Labels
domain:linear algebra Linear algebra

Comments

@jishnub
Copy link
Contributor

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?
@jishnub jishnub added the domain:linear algebra Linear algebra label Feb 13, 2023
@aravindh-krishnamoorthy
Copy link
Contributor

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

@jishnub
Copy link
Contributor 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
Contributor 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 added a commit that referenced this issue Apr 7, 2024
@jishnub
Support broadcasting over structured block matrices (#53909)
243ebc3 
Fix #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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
domain:linear algebra Linear algebra
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