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.

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

Fix size-1 StructuredMatrix's broadcast. #54190

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

N5N3
Copy link
Member

@N5N3 N5N3 commented Apr 22, 2024

  1. size-1 StructuredMatrix should behave like scalar during broadcast. Thus their fzero should return the only element.

  2. But for simple broadcast with only one StructuredMatrix, we can keep stability as the structure is "preserved" even for size-1 case. Thus count_structedmatrix is added to check that.

  3. count_structedmatrix is fused to keep Bidiagonal stability.

fix JuliaLang/LinearAlgebra.jl#1063
replace #54067

1. size-1 StructuredMatrix should behave like scalar during broadcast. Thus their `fzero` should return the only element.
(fix #54087)

2. But for simple broadcast with only one StructuredMatrix, we can keep stability as the structure is "preserved" even for size-1 case. Thus `count_structedmatrix` is added to check that.

3. `count_structedmatrix` is fused to keep `Bidiagonal` stability.
(replace JuliaLang#54067)
@N5N3 N5N3 added linear algebra Linear algebra broadcast Applying a function over a collection labels Apr 22, 2024
@N5N3 N5N3 requested a review from jishnub April 22, 2024 15:24
@jishnub
Copy link
Contributor

jishnub commented Apr 23, 2024

The impact of the type-instability is small on performance

julia> D = Diagonal(ones(10000));

julia> @btime $D .+ $D;
  7.359 μs (5 allocations: 78.22 KiB)

julia> @btime $D + $D;
  7.330 μs (3 allocations: 78.19 KiB)

but I wonder if this might lead to cascading type-instabilities downstream if one tries to store the result in a struct/array. The single-element structured matrix is a rather niche use case as such. I wonder if this should just be deprecated, with a suggestion to use a Matrix instead?

@N5N3
Copy link
Member Author

N5N3 commented Apr 23, 2024

The single-element structured matrix is a rather niche use case as such.

I agree the size-1 case has no realistic usage, but throwing for inputs which are valid for general fallback looks too broken.

but I wonder if this might lead to cascading type-instabilities downstream

I guess there's no need to worry too much for this.
The stability of StructuredMatrix's broadcast is also fragile on master, e.g.

julia> bad(T) = T .* 2
bad (generic function with 1 method)

julia> good(T) = T .* 2.
good (generic function with 1 method)

julia> T = Tridiagonal(1:1,1:2,1:1)
2×2 Tridiagonal{Int64, UnitRange{Int64}}:
 1  1
 1  2

julia> @code_warntype good(T)
MethodInstance for good(::Tridiagonal{Int64, UnitRange{Int64}})
  from good(T) @ Main REPL[2]:1
Arguments
  #self#::Core.Const(good)
  T::Tridiagonal{Int64, UnitRange{Int64}}
Body::Tridiagonal{Float64, Vector{Float64}}
1%1 = Main.:*::Core.Const(*)
│   %2 = Base.broadcasted(%1, T, 2.0)::Core.PartialStruct(Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Tridiagonal}, Nothing, typeof(*), Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Float64}}, Any[Core.Const(LinearAlgebra.StructuredMatrixStyle{Tridiagonal}()), Core.Const(*), Core.PartialStruct(Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Float64}, Any[Tridiagonal{Int64, UnitRange{Int64}}, Core.Const(2.0)]), Nothing])
│   %3 = Base.materialize(%2)::Tridiagonal{Float64, Vector{Float64}}
└──      return %3


julia> @code_warntype bad(T)
MethodInstance for bad(::Tridiagonal{Int64, UnitRange{Int64}})
  from bad(T) @ Main REPL[1]:1
Arguments
  #self#::Core.Const(bad)
  T::Tridiagonal{Int64, UnitRange{Int64}}
Body::Union{Tridiagonal{Int64, Vector{Int64}}, Matrix{Int64}}
1%1 = Main.:*::Core.Const(*)
│   %2 = Base.broadcasted(%1, T, 2)::Core.PartialStruct(Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Tridiagonal}, Nothing, typeof(*), Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Int64}}, Any[Core.Const(LinearAlgebra.StructuredMatrixStyle{Tridiagonal}()), Core.Const(*), Core.PartialStruct(Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Int64}, Any[Tridiagonal{Int64, UnitRange{Int64}}, Core.Const(2)]), Nothing])
│   %3 = Base.materialize(%2)::Union{Tridiagonal{Int64, Vector{Int64}}, Matrix{Int64}}
└──      return %3

Users can always use typeassert to restore stability where possible.

@jishnub
Copy link
Contributor

jishnub commented Apr 23, 2024

I'm not really in favor of making this type-unstable to address a rarely encountered case, but I wonder what @mbauman might have to say on this

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
broadcast Applying a function over a collection linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Structured matrices don't extend missing dimensions correctly in broadcasting
2 participants