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

fix equality of QRCompactWY #41363

Merged
merged 4 commits into from
Jun 28, 2021
Merged

fix equality of QRCompactWY #41363

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6679a39
fix equality of QRCompactWY
simeonschaub Jun 14, 2021
72e3b41
fix for matrix sizes >36
simeonschaub Jun 27, 2021
1c7dcbc
fix missings test
simeonschaub Jun 27, 2021
835f76b
Merge remote-tracking branch 'origin/master' into sds/fix_qr_eq
simeonschaub Jun 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions stdlib/LinearAlgebra/src/qr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,16 @@ Base.iterate(S::QRCompactWY) = (S.Q, Val(:R))
Base.iterate(S::QRCompactWY, ::Val{:R}) = (S.R, Val(:done))
Base.iterate(S::QRCompactWY, ::Val{:done}) = nothing

function Base.hash(F::QRCompactWY, h::UInt)
return hash(F.factors, hash(UpperTriangular(F.T), hash(QRCompactWY, h)))
Copy link
Member

Choose a reason for hiding this comment

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

This is not correct. Generally, the T matrix will not be square but wide where the number of rows is a block size, see

julia/stdlib/LinearAlgebra/src/qr.jl

Line 249 in 1e3ae23

qr!(A::StridedMatrix{<:BlasFloat}, ::NoPivot; blocksize=36) =
, and the number of columns is min(size(A)...). For matrices smaller than the block size T will be square but for larger matrices, it will consist of several triangular blocks side by side, so the hashing here will be slightly more complicated and depend on the block size.

Copy link
Member Author

@simeonschaub simeonschaub Jun 27, 2021 •
edited
Loading

Choose a reason for hiding this comment

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

Ah, I see! (Random appreciation for the new sparse matrix printing:

julia> using LinearAlgebra, SparseArrays

julia> a = rand(100, 100);

julia> sparse(qr(a).T .== qr(a).T)
36×100 SparseMatrixCSC{Bool, Int64} with 1767 stored entries:
⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿
⠀⠀⠀⠀⠀⠂⠛⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿
⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⢀⠐⠙⢿⣿⣿⣿⣿
⠀⠀⠐⠀⠀⠀⠀⠀⠀⢀⢙⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠁⠀⡀⠀⠙⢿⣿⣿
⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠄⠀⠙⢿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⡀⠀⠀⢀⠀⠀⠙⢿
⠀⡀⠀⠀⠀⠀⠀⠀⠂⠒⠒⠀⠀⠀⠙⢿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣿⣿⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⡀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⣈⡀⠀⠀⠀⠀⠀⠀⠙⢿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠂⠀⢀⠀

)

end
function Base.:(==)(A::QRCompactWY, B::QRCompactWY)
return A.factors == B.factors && UpperTriangular(A.T) == UpperTriangular(B.T)
end
function Base.isequal(A::QRCompactWY, B::QRCompactWY)
return isequal(A.factors, B.factors) && isequal(UpperTriangular(A.T), UpperTriangular(B.T))
end

"""
QRPivoted <: Factorization

Expand Down
45 changes: 45 additions & 0 deletions stdlib/LinearAlgebra/test/factorization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestFactorization
using Test, LinearAlgebra

@testset "equality for factorizations - $f" for f in Any[
bunchkaufman,
cholesky,
x -> cholesky(x, Val(true)),
eigen,
hessenberg,
lq,
lu,
qr,
x -> qr(x, ColumnNorm()),
svd,
schur,
]
A = randn(3, 3)
A = A * A' # ensure A is pos. def. and symmetric
F, G = f(A), f(A)

@test F == G
@test isequal(F, G)
@test hash(F) == hash(G)

f === hessenberg && continue

# change all arrays in F to have eltype Float32
F = typeof(F).name.wrapper(Base.mapany(1:nfields(F)) do i
x = getfield(F, i)
return x isa AbstractArray{Float64} ? Float32.(x) : x
end...)
# round all arrays in G to the nearest Float64 representable as Float32
G = typeof(G).name.wrapper(Base.mapany(1:nfields(G)) do i
x = getfield(G, i)
return x isa AbstractArray{Float64} ? Float64.(Float32.(x)) : x
end...)

@test F == G broken=!(f === eigen || f === qr)
@test isequal(F, G) broken=!(f === eigen || f === qr)
@test hash(F) == hash(G)
end

end
1 change: 1 addition & 0 deletions stdlib/LinearAlgebra/test/testgroups
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,3 +25,4 @@ givens
structuredbroadcast
addmul
ldlt
factorization