Define det(Q) for {QR,LQ}PackedQ and QRCompactWYQ (#32887)

* Define det(Q) = ±1 for {QR,LQ}PackedQ

* Handle general {QR,LQ}PackedQ in det(Q)

* Inline x'y to remove O(n) allocations in det(Q)

* Revert "Inline x'y to remove O(n) allocations in det(Q)"

This reverts commit 940d16c.

* Revert "Handle general {QR,LQ}PackedQ in det(Q)"

This reverts commit f152120.

* Implement det(Q) for non-real case

* Implement det(Q::QRCompactWYQ)

* Fix non-real Q; handle iszero(τ) explicitly

* Fix det(Q::LQPackedQ); reflectors are implicitly adjoint

* Return ±one(eltype(Q)) instead of ±1::Int

* Support blocked QRCompactWYQ

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.τ))

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)

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

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