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 edge cases in SymTridiagonal when ev has an extra element (+, -, iszero, isone, etc.) #42472

Merged
merged 3 commits into from
Oct 8, 2021

Conversation

mcognetta
Copy link
Contributor

@mcognetta mcognetta commented Oct 2, 2021

In SymTridiagonal, ev can have the same number of elements as dv, but the last element in ev is typically ignored. This can, however, cause some issues for things that just check all elements of ev. For example, iszero, isone, etc can all fail due to a particular choice of the final element of ev, even though that element is not used.

This PR fixes the issue for iszero, isone, isdiag, istriu, and istril. See also #41089.

Update:

This bug is also present in some constructors (see JuliaLang/LinearAlgebra.jl#874) and in some binary operations (+/- when one operand is SymTridiagonal and the other being Tridiagonal, Bidiagonal, etc when ev has the same number of elements as ev.

This PR was updated to include fixes for these cases as well.


Head:

julia> dv = rand(3)
3-element Vector{Float64}:
 0.26763492657215415
 0.24104518570135736
 0.3485385690926436

julia> ev = zeros(3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> ev[end] = 1
1

julia> S = SymTridiagonal(dv, ev)
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 0.267635  0.0         
 0.0       0.241045  0.0
          0.0       0.348539

julia> isdiag(S)
false

julia> S = SymTridiagonal(ones(3), ev)
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 1.0  0.0    
 0.0  1.0  0.0
     0.0  1.0

julia> isone(S)
false

julia> S = SymTridiagonal(zeros(3), ev)
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 0.0  0.0    
 0.0  0.0  0.0
     0.0  0.0

julia> iszero(S)
false

This PR:

julia> dv = rand(3)
3-element Vector{Float64}:
 0.6174155214734693
 0.1508335261390339
 0.16092821432071902

julia> ev = zeros(3)
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> ev[end] = 1
1

julia> S = SymTridiagonal(dv, ev)
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 0.617416  0.0         
 0.0       0.150834  0.0
          0.0       0.160928

julia> isdiag(S)
true

julia> S = SymTridiagonal(ones(3), ev)
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 1.0  0.0    
 0.0  1.0  0.0
     0.0  1.0

julia> isone(S)
true

julia> S = SymTridiagonal(zeros(3), ev)
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 0.0  0.0    
 0.0  0.0  0.0
     0.0  0.0

julia> iszero(S)
true

@mcognetta
Copy link
Contributor Author

Seems like there are a few more that need to be added here:

julia> SymTridiagonal(rand(3), rand(2)) + SymTridiagonal(rand(3), rand(3))
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(2),), b has dims (Base.OneTo(3),), mismatch at 1")
Stacktrace:
 [1] promote_shape
   @ ./indices.jl:178 [inlined]
 [2] promote_shape
   @ ./indices.jl:169 [inlined]
 [3] +(A::Vector{Float64}, Bs::Vector{Float64})
   @ Base ./arraymath.jl:45
 [4] +(A::SymTridiagonal{Float64, Vector{Float64}}, B::SymTridiagonal{Float64, Vector{Float64}})
   @ LinearAlgebra ~/github/julia-dev/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/tridiag.jl:210
 [5] top-level scope
   @ REPL[2]:1

julia> SymTridiagonal(rand(3), rand(2)) - SymTridiagonal(rand(3), rand(3))
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(2),), b has dims (Base.OneTo(3),), mismatch at 1")
Stacktrace:
 [1] promote_shape
   @ ./indices.jl:178 [inlined]
 [2] promote_shape
   @ ./indices.jl:169 [inlined]
 [3] -(A::Vector{Float64}, B::Vector{Float64})
   @ Base ./arraymath.jl:38
 [4] -(A::SymTridiagonal{Float64, Vector{Float64}}, B::SymTridiagonal{Float64, Vector{Float64}})
   @ LinearAlgebra ~/github/julia-dev/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/tridiag.jl:211
 [5] top-level scope
   @ REPL[3]:1

@mcognetta mcognetta changed the title Fix edge cases in SymTridiagonal when ev has an extra element (iszero, isone, etc.) Fix edge cases in SymTridiagonal when ev has an extra element (+, -, iszero, isone, etc.) Oct 3, 2021
@mcognetta mcognetta force-pushed the symtridiag_ev_last_element branch from 89aff82 to 6da5108 Compare October 3, 2021 00:53
@mcognetta mcognetta force-pushed the symtridiag_ev_last_element branch from 67953dd to d199de0 Compare October 3, 2021 04:21
@kshyatt kshyatt added the linear algebra Linear algebra label Oct 4, 2021
@kshyatt kshyatt requested a review from andreasnoack October 4, 2021 15:07
Copy link
Member

@andreasnoack andreasnoack left a comment

Choose a reason for hiding this comment

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

LGTM. There is a small conflict that needs to be resolved.

@mcognetta
Copy link
Contributor Author

LGTM. There is a small conflict that needs to be resolved.

Should be good to go now. Thanks for the review.

@dkarrasch dkarrasch merged commit 146de38 into JuliaLang:master Oct 8, 2021
LilithHafner pushed a commit to LilithHafner/julia that referenced this pull request Feb 22, 2022
LilithHafner pushed a commit to LilithHafner/julia that referenced this pull request Mar 8, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants