diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index c2366a5bf432d..7b3829023b4ca 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -334,3 +334,8 @@ function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T lmul!(adjoint(A.Q), ldiv!(LowerTriangular(A.L),B)) return B end + + +# In LQ factorization, `Q` is expressed as the product of the adjoint of the +# reflectors. Thus, `det` has to be conjugated. +det(Q::LQPackedQ) = conj(_det_tau(Q.τ)) diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index a14cb658bdc1c..48949aa393284 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -917,3 +917,25 @@ end ## Lower priority: Add LQ, QL and RQ factorizations # FIXME! Should add balancing option through xgebal + + +det(Q::QRPackedQ) = _det_tau(Q.τ) + +det(Q::QRCompactWYQ) = + prod(i -> _det_tau(_diagview(Q.T[:, i:min(i + size(Q.T, 1), size(Q.T, 2))])), + 1:size(Q.T, 1):size(Q.T, 2)) + +_diagview(A) = @view A[diagind(A)] + +# Compute `det` from the number of Householder reflections. Handle +# the case `Q.τ` contains zeros. +_det_tau(τs::AbstractVector{<:Real}) = + isodd(count(!iszero, τs)) ? -one(eltype(τs)) : one(eltype(τs)) + +# In complex case, we need to compute the non-unit eigenvalue `λ = 1 - c*τ` +# (where `c = v'v`) of each Householder reflector. As we know that the +# reflector must have the determinant of 1, it must satisfy `abs2(λ) == 1`. +# Combining this with the constraint `c > 0`, it turns out that the eigenvalue +# (hence the determinant) can be computed as `λ = -sign(τ)^2`. +# See: https://github.com/JuliaLang/julia/pull/32887#issuecomment-521935716 +_det_tau(τs) = prod(τ -> iszero(τ) ? one(τ) : -sign(τ)^2, τs) diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl index be03dc0125b5a..1318124e883e1 100644 --- a/stdlib/LinearAlgebra/test/lq.jl +++ b/stdlib/LinearAlgebra/test/lq.jl @@ -191,4 +191,19 @@ end @test_throws DimensionMismatch adjoint(C) * adjoint(implicitQ) end +@testset "det(Q::LQPackedQ)" begin + @testset for n in 1:3, m in 1:3 + @testset "real" begin + _, Q = lq(randn(n, m)) + @test det(Q) ≈ det(collect(Q)) + @test abs(det(Q)) ≈ 1 + end + @testset "complex" begin + _, Q = lq(randn(ComplexF64, n, m)) + @test det(Q) ≈ det(collect(Q)) + @test abs(det(Q)) ≈ 1 + end + end +end + end # module TestLQ diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index b3337d48a166a..9fdde9df0f7e5 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -246,4 +246,26 @@ end @test c0 == c end +@testset "det(Q::Union{QRCompactWYQ, QRPackedQ})" begin + # 40 is the number larger than the default block size 36 of QRCompactWY + @testset for n in [1:3; 40], m in [1:3; 40], pivot in [false, true] + @testset "real" begin + @testset for k in 0:min(n, m, 5) + A = cat(Array(I(k)), randn(n - k, m - k); dims=(1, 2)) + Q, = qr(A, Val(pivot)) + @test det(Q) ≈ det(collect(Q)) + @test abs(det(Q)) ≈ 1 + end + end + @testset "complex" begin + @testset for k in 0:min(n, m, 5) + A = cat(Array(I(k)), randn(ComplexF64, n - k, m - k); dims=(1, 2)) + Q, = qr(A, Val(pivot)) + @test det(Q) ≈ det(collect(Q)) + @test abs(det(Q)) ≈ 1 + end + end + end +end + end # module TestQR