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

Fix multiplication of AbstractQs #46237

Merged
merged 3 commits into from
Aug 22, 2022
Merged

Fix multiplication of AbstractQs #46237

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f73de47
Fix multiplication of `AbstractQ`s
dkarrasch Aug 1, 2022
16b1976
fix typo
dkarrasch Aug 1, 2022
f1bf60f
expand tests
dkarrasch Aug 1, 2022
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
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/src/qr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -380,9 +380,9 @@ norm solution.

Multiplication with respect to either full/square or non-full/square `Q` is allowed, i.e. both `F.Q*F.R`
and `F.Q*A` are supported. A `Q` matrix can be converted into a regular matrix with
[`Matrix`](@ref). This operation returns the "thin" Q factor, i.e., if `A` is `m`×`n` with `m>=n`, then
[`Matrix`](@ref). This operation returns the "thin" Q factor, i.e., if `A` is `m`×`n` with `m>=n`, then
`Matrix(F.Q)` yields an `m`×`n` matrix with orthonormal columns. To retrieve the "full" Q factor, an
`m`×`m` orthogonal matrix, use `F.Q*Matrix(I,m,m)`. If `m<=n`, then `Matrix(F.Q)` yields an `m`×`m`
`m`×`m` orthogonal matrix, use `F.Q*I`. If `m<=n`, then `Matrix(F.Q)` yields an `m`×`m`
orthogonal matrix.

The block size for QR decomposition can be specified by keyword argument
Expand Down
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/src/special.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -347,10 +347,10 @@ end
*(A::Diagonal, Q::AbstractQ) = _qrmul(A, Q)
*(A::Diagonal, Q::Adjoint{<:Any,<:AbstractQ}) = _qrmul(A, Q)

*(Q::AbstractQ, B::AbstractQ) = _qlmul(Q, B)
*(Q::Adjoint{<:Any,<:AbstractQ}, B::AbstractQ) = _qrmul(Q, B)
*(Q::AbstractQ, B::Adjoint{<:Any,<:AbstractQ}) = _qlmul(Q, B)
*(Q::Adjoint{<:Any,<:AbstractQ}, B::Adjoint{<:Any,<:AbstractQ}) = _qrmul(Q, B)
*(Q::AbstractQ, B::AbstractQ) = Q * (B * I)
*(Q::Adjoint{<:Any,<:AbstractQ}, B::AbstractQ) = Q * (B * I)
*(Q::AbstractQ, B::Adjoint{<:Any,<:AbstractQ}) = Q * (B * I)
*(Q::Adjoint{<:Any,<:AbstractQ}, B::Adjoint{<:Any,<:AbstractQ}) = Q * (B * I)

# fill[stored]! methods
fillstored!(A::Diagonal, x) = (fill!(A.diag, x); A)
Expand Down
33 changes: 23 additions & 10 deletions stdlib/LinearAlgebra/test/special.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,16 +215,29 @@ end
atri = typ(a)
matri = Matrix(atri)
b = rand(n,n)
qrb = qr(b, ColumnNorm())
@test atri * qrb.Q ≈ matri * qrb.Q ≈ rmul!(copy(atri), qrb.Q)
@test atri * qrb.Q' ≈ matri * qrb.Q' ≈ rmul!(copy(atri), qrb.Q')
@test qrb.Q * atri ≈ qrb.Q * matri ≈ lmul!(qrb.Q, copy(atri))
@test qrb.Q' * atri ≈ qrb.Q' * matri ≈ lmul!(qrb.Q', copy(atri))
qrb = qr(b, NoPivot())
@test atri * qrb.Q ≈ matri * qrb.Q ≈ rmul!(copy(atri), qrb.Q)
@test atri * qrb.Q' ≈ matri * qrb.Q' ≈ rmul!(copy(atri), qrb.Q')
@test qrb.Q * atri ≈ qrb.Q * matri ≈ lmul!(qrb.Q, copy(atri))
@test qrb.Q' * atri ≈ qrb.Q' * matri ≈ lmul!(qrb.Q', copy(atri))
for pivot in (ColumnNorm(), NoPivot())
qrb = qr(b, pivot)
@test atri * qrb.Q ≈ matri * qrb.Q ≈ rmul!(copy(atri), qrb.Q)
@test atri * qrb.Q' ≈ matri * qrb.Q' ≈ rmul!(copy(atri), qrb.Q')
@test qrb.Q * atri ≈ qrb.Q * matri ≈ lmul!(qrb.Q, copy(atri))
@test qrb.Q' * atri ≈ qrb.Q' * matri ≈ lmul!(qrb.Q', copy(atri))
end
end
end

@testset "Multiplication of Qs" begin
for pivot in (ColumnNorm(), NoPivot()), A in (rand(5, 3), rand(5, 5), rand(3, 5))
Q = qr(A, pivot).Q
m = size(A, 1)
C = Matrix{Float64}(undef, (m, m))
@test Q*Q ≈ (Q*I) * (Q*I) ≈ mul!(C, Q, Q)
@test size(Q*Q) == (m, m)
@test Q'Q ≈ (Q'*I) * (Q*I) ≈ mul!(C, Q', Q)
@test size(Q'Q) == (m, m)
@test Q*Q' ≈ (Q*I) * (Q'*I) ≈ mul!(C, Q, Q')
@test size(Q*Q') == (m, m)
@test Q'Q' ≈ (Q'*I) * (Q'*I) ≈ mul!(C, Q', Q')
@test size(Q'Q') == (m, m)
end
end

Expand Down