Fix incorrect dispatch for Cholesky equality

[penelopeysm]

== on Cholesky is not correctly being dispatched when the two operands have exactly the same type:

julia> x = LowerTriangular([0.0 1.0; 0.0 0.0])
2×2 LowerTriangular{Float64, Matrix{Float64}}:
 0.0   ⋅
 0.0  0.0

julia> y = LowerTriangular([0.0 0.0; 0.0 0.0])
2×2 LowerTriangular{Float64, Matrix{Float64}}:
 0.0   ⋅
 0.0  0.0

julia> x == y
true

julia> Cholesky(x) == Cholesky(y)
false

julia> @which Cholesky(x) == Cholesky(y)
==(F::T, G::T) where T<:Factorization
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/factorization.jl:116

That's because that method

LinearAlgebra.jl/src/factorization.jl

Line 116 in 2c3fe9b

Base.:(==)( F::T, G::T) where {T<:Factorization} = all(f -> getfield(F, f) == getfield(G, f), 1:nfields(F))

is more specific than the intended method

LinearAlgebra.jl/src/cholesky.jl

Lines 683 to 686 in 2c3fe9b

	function Base.:(==)(C1::Cholesky, C2::Cholesky)
	C1.uplo == C2.uplo || return false
	C1.uplo == 'L' ? (C1.L == C2.L) : (C1.U == C2.U)
	end

This PR fixes it by adding an extra, more specific method.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.86%. Comparing base (2c3fe9b) to head (fa3cd0f).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1404      +/-   ##
==========================================
- Coverage   93.87%   93.86%   -0.02%     
==========================================
  Files          34       34              
  Lines       15830    15830              
==========================================
- Hits        14861    14859       -2     
- Misses        969      971       +2     

☔ 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.

[dkarrasch]

Thanks for identifying the issue and the cause, and for proposing a fix. I have an alternative suggestion and another thought:

There are a few more instances that probably don't work as intended (or at least worth checking):

julia> methods(==, (Factorization, Factorization))
# 5 methods for generic function "==" from Base:
 [1] ==(A::LinearAlgebra.QRCompactWY, B::LinearAlgebra.QRCompactWY)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:160
 [2] ==(C1::C, C2::D) where {C<:Cholesky, D<:Cholesky}
     @ LinearAlgebra REPL[2]:1
 [3] ==(A::Eigen, B::Eigen)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:664
 [4] ==(F::T, G::T) where T<:Factorization
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/factorization.jl:116
 [5] ==(x, y)
     @ Base.jl:207

So, Eigen and QRCompactWY are further candidates.

[penelopeysm]

So, Eigen and QRCompactWY are further candidates.

Thanks for pulling that up! I took a look. I think Eigen is fine because it doesn't have special equality behaviour (it just compares the fields, and its fields don't carry any redundant information, so it seems to me it doesn't matter which method it dispatches to).

As for QRCompactWY, I have to confess I am completely baffled, because it dispatches to the correct method:

julia> m = [1.0 2.0; 3.0 4.0]
2×2 Matrix{Float64}:
 1.0  2.0
 3.0  4.0

julia> x = y = qr(m)
LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}
Q factor: 2×2 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}
R factor:
2×2 Matrix{Float64}:
 -3.16228  -4.42719
  0.0      -0.632456

julia> @which x == y
==(A::LinearAlgebra.QRCompactWY, B::LinearAlgebra.QRCompactWY)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:160

in fact Eigen dispatches to the right method too:

julia> p = q = eigen(m)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 -0.3722813232690143
  5.372281323269014
vectors:
2×2 Matrix{Float64}:
 -0.824565  -0.415974
  0.565767  -0.909377

julia> @which p == q
==(A::Eigen, B::Eigen)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/eigen.jl:664

whereas this doesn't happen with Cholesky (pre-this PR, of course):

julia> l = LowerTriangular(zeros(2, 2))
2×2 LowerTriangular{Float64, Matrix{Float64}}:
 0.0   ⋅
 0.0  0.0

julia> a = b = Cholesky(l)
Cholesky{Float64, Matrix{Float64}}
L factor:
2×2 LowerTriangular{Float64, Matrix{Float64}}:
 0.0   ⋅
 0.0  0.0

julia> @which a == b
==(F::T, G::T) where T<:Factorization
     @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/factorization.jl:116

Do you have any idea why?

[dkarrasch]

Ah, that's because the double-Cholesky method hasn't been released yet: c68619a

It will be included in v1.12, so this PR can be closed (sorry for that).

[dkarrasch]

I found out after realizing that on v1.11 there is no such method in the list of methods.

[penelopeysm]

I was just looking into that! I did see the commit was from a year ago and I thought that it would have made it in by now, but clearly not.

Thanks for reviewing and sorry to have taken up your time :)

[penelopeysm]

(Just out of curiosity, since I've never gone as far upstream as the standard library: Do changes in LinearAlgebra.jl only get merged into new minor versions of Julia, i.e., new Julia patch releases don't get bugfixes from standard libraries?)

[dkarrasch]

That's actually a good point. That commit might be considered a bugfix, and in that case should be backported to the LTS and the current version. @jishnub What do you think?

BTW, you can see at the top if some commit is contained in a release. For the relevant commit, it just says main or master, otherwise it would mention some version (perhaps try older linalg commits on the main julia repo to see the difference).

Thanks for reviewing and sorry to have taken up your time :)

No worries at all!