Skip to content

Preserve structured matrices for zero-absorbing operations#1549

Open
max-vassili3v wants to merge 8 commits intoJuliaLang:masterfrom
max-vassili3v:sparse-broadcasting
Open

Preserve structured matrices for zero-absorbing operations#1549
max-vassili3v wants to merge 8 commits intoJuliaLang:masterfrom
max-vassili3v:sparse-broadcasting

Conversation

@max-vassili3v
Copy link

@max-vassili3v max-vassili3v commented Feb 3, 2026

This PR fixes #1543 (comment).

It overloads fzero such that broadcasts involving structured matrices and zero-absorbing operations (*, dot, etc.) will return true for fzeropreserving. Logic was also added for left-absorbing operations (/, rem, etc.) where we attempt to preserve structure if the structured matrix is on the left and as with the scalar case, fall back to a dense matrix if division by zero occurs. There is also a check in the case of SymTridiagonal, where broadcasts with zero-absorbing operations can avoid materializing a dense matrix, but might need to become a Tridiagonal.

Unit tests have been added to cover the left-absorbing functionality, and some existing ones have been changed. Specifically, typeof(broadcast(*, s, fV, fA, X)) should now be whatever StructuredMatrix X is, except for Unit***erTriangular which becomes ***erTriangular and SymTridiagonal which becomes Tridiagonal for N > 1.

@codecov
Copy link

codecov bot commented Feb 3, 2026

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.07%. Comparing base (9c9d374) to head (ccc662c).

Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1549      +/-   ##
==========================================
+ Coverage   94.03%   94.07%   +0.03%     
==========================================
  Files          35       35              
  Lines       15986    16001      +15     
==========================================
+ Hits        15033    15053      +20     
+ Misses        953      948       -5     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@dlfivefifty
Copy link
Contributor

Is the untested function definition actually needed?

# broadcasts with SymTridiagonal and arrays with more than one entry will generally break symmetry.
# This is useful for knowing whether to materialize a Tridiagonal for zero-preserving functions.
function issymmetrybreaking(bc::Broadcasted{StructuredMatrixStyle{SymTridiagonal}})
any(x -> !(x isa SymTridiagonal || all(y -> length(y) == 1, axes(x))), bc.args)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there much of an advantage in preserving the structure for 1x1 matrices here, given that this comes at the price of type stability?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made it so that 1x1 SymTridiagonals always become Tridiagonals, and have changed the unit tests accordingly

@jishnub
Copy link
Member

jishnub commented Feb 7, 2026

I noticed that the following seems to have broken on this PR but is working on master:

julia> using LinearAlgebra

julia> D = Diagonal(1:4)
4×4 Diagonal{Int64, UnitRange{Int64}}:
 1      
   2    
     3  
       4

julia> (1:4) .* D .* (1:4)'
ERROR: MethodError: no method matching diagview(::UnitRange{Int64}, ::Int64)
The function `diagview` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  diagview(::Adjoint, ::Integer)
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/adjtrans.jl:579
  diagview(::Transpose, ::Integer)
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/adjtrans.jl:578
  diagview(::AbstractMatrix, ::Integer)
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/dense.jl:358
  ...

Stacktrace:
 [1] diagview(A::Adjoint{Int64, UnitRange{Int64}}, k::Int64)
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/adjtrans.jl:579
 [2] _preprocess_broadcasted(::Type{Diagonal}, d::Adjoint{Int64, UnitRange{Int64}})
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:286
 [3] preprocess_broadcasted(::Type{Diagonal}, A::Adjoint{Int64, UnitRange{Int64}})
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:276
 [4] (::LinearAlgebra.var"#preprocess_broadcasted##0#preprocess_broadcasted##1"{Diagonal})(x::Adjoint{Int64, UnitRange{Int64}})
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:278
 [5] map
   @ ./tuple.jl:358 [inlined]
 [6] preprocess_broadcasted(::Type{…}, bc::Base.Broadcast.Broadcasted{…})
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:278
 [7] copy(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{…}, Tuple{…}, typeof(*), Tuple{…}})
   @ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:291
 [8] materialize(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{…}, Nothing, typeof(*), Tuple{…}})
   @ Base.Broadcast ./broadcast.jl:893
 [9] top-level scope
   @ REPL[2]:1
Some type information was truncated. Use `show(err)` to see complete types.

This isn't the fault of this PR, so I'll try to address this separately.

@jishnub jishnub reopened this Feb 12, 2026
@jishnub
Copy link
Member

jishnub commented Feb 12, 2026

Ah sorry didn't mean to close this

jishnub added a commit to max-vassili3v/LinearAlgebra.jl that referenced this pull request Feb 21, 2026
Should fix
JuliaLang#1549 (comment).

After this, the following works:
```julia
julia> a = (1:4)'
1×4 adjoint(::UnitRange{Int64}) with eltype Int64:
 1  2  3  4

julia> diagview(a)
1-element view(::UnitRange{Int64}, 1:2:1) with eltype Int64:
 1
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Broadcasting * with Diagonal should return a Diagonal (similarly with Bidiagonal/Tridiagonal)

3 participants