R-only QR and LQ decompositions

.gitignore

@@ -4,3 +4,4 @@
 .*.swp
 Manifest.toml
 docs/build/
+.vscode

ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl

@@ -66,19 +66,16 @@ function qr_householder!(
 
     m, n = size(A)
     minmn = min(m, n)
-    computeR = length(R) > 0
 
     # compute QR
     Q̃, R̃ = qr!(A)
-    lmul!(Q̃, one!(Q))
+    isempty(Q) || lmul!(Q̃, one!(Q))
 
-    if computeR
+    if !isempty(R)
         copyto!(view(R, 1:minmn, :), R̃)
         zero!(view(R, (minmn + 1):size(R, 1), :))
-        positive && gaugefix!(qr_householder!, Q, R, diagview(R̃))
-    elseif positive
-        gaugefix!(qr_householder!, Q, nothing, diagview(R̃))
     end
+    positive && gaugefix!(qr_householder!, Q, R, diagview(R̃))
 
     return Q, R
 end

src/common/gauge.jl

@@ -15,9 +15,11 @@ _largest(x::Complex, y::Complex) = abs2(x) < abs2(y) ? y : x
 
 function gaugefix!(::typeof(qr_householder!), Q, R, Rd)
     ax = Base.OneTo(length(Rd))
-    Qf = view(Q, axes(Q, 1), ax)
-    Qf .*= sign_safe.(transpose(Rd))
-    if !isnothing(R)
+    if !isempty(Q)
+        Qf = view(Q, axes(Q, 1), ax)
+        Qf .*= sign_safe.(transpose(Rd))
+    end
+    if !isempty(R)
         Rf = view(R, ax, axes(R, 2))
         Rf .*= conj.(sign_safe.(Rd))
     end
@@ -26,9 +28,11 @@ end
 
 function gaugefix!(::typeof(lq_householder!), L, Q, Ld)
     ax = Base.OneTo(length(Ld))
-    Qf = view(Q, ax, axes(Q, 2))
-    Qf .*= sign_safe.(Ld)
-    if !isnothing(L)
+    if !isempty(Q)
+        Qf = view(Q, ax, axes(Q, 2))
+        Qf .*= sign_safe.(Ld)
+    end
+    if !isempty(L)
         Lf = view(L, axes(L, 1), ax)
         Lf .*= conj.(sign_safe.(transpose(Ld)))
     end

src/implementations/lq.jl

@@ -12,7 +12,7 @@ function check_input(::typeof(lq_full!), A::AbstractMatrix, LQ, ::AbstractAlgori
     @assert L isa AbstractMatrix && Q isa AbstractMatrix
     isempty(L) || @check_size(L, (m, n))
     @check_scalar(L, A)
-    @check_size(Q, (n, n))
+    isempty(Q) || @check_size(Q, (n, n))
     @check_scalar(Q, A)
     return nothing
 end
@@ -23,7 +23,7 @@ function check_input(::typeof(lq_compact!), A::AbstractMatrix, LQ, ::AbstractAlg
     @assert L isa AbstractMatrix && Q isa AbstractMatrix
     isempty(L) || @check_size(L, (m, minmn))
     @check_scalar(L, A)
-    @check_size(Q, (minmn, n))
+    isempty(Q) || @check_size(Q, (minmn, n))
     @check_scalar(Q, A)
     return nothing
 end
@@ -42,7 +42,7 @@ function check_input(::typeof(lq_full!), A::AbstractMatrix, (L, Q), ::DiagonalAl
     @assert Q isa Diagonal && L isa Diagonal
     isempty(L) || @check_size(L, (m, n))
     @check_scalar(L, A)
-    @check_size(Q, (n, n))
+    isempty(Q) || @check_size(Q, (n, n))
     @check_scalar(Q, A)
     return nothing
 end
@@ -148,7 +148,7 @@ function lq_householder!(
 
     m, n = size(A)
     minmn = min(m, n)
-    computeL = length(L) > 0
+    computeL = !isempty(L)
     inplaceQ = Q === A
 
     if inplaceQ
@@ -162,25 +162,25 @@ function lq_householder!(
     if blocksize > 1
         # L doubles as workspace for T, so Q is constructed before L is extracted
         mb = min(minmn, blocksize)
-        if computeL # first use L as space for T
-            A, T = gelqt!(driver, A, view(L, 1:mb, 1:minmn))
-        else
-            A, T = gelqt!(driver, A, similar(A, mb, minmn))
-        end
-        Q = gemlqt!(driver, 'R', 'N', A, T, one!(Q))
+        T = computeL ? view(L, 1:mb, 1:minmn) : similar(A, mb, minmn)
+        A, T = gelqt!(driver, A, T)
+        isempty(Q) || (Q = gemlqt!(driver, 'R', 'N', A, T, one!(Q)))
         computeL && _lq_copyL!(L, A)
-        positive && gaugefix!(lq_householder!, computeL ? L : nothing, Q, diagview(A))
+        positive && gaugefix!(lq_householder!, L, Q, diagview(A))
     else
         A, τ = gelqf!(driver, A)
         computeL && _lq_copyL!(L, A)
-        Lf = computeL ? L : nothing
         if inplaceQ
-            Ld = positive ? copy(diagview(A)) : nothing # unglq! destroys the diagonal of A
-            unglq!(driver, A, τ) # Q === A, so no need to rebind Q
-            positive && gaugefix!(lq_householder!, Lf, Q, Ld)
+            if positive
+                Ld = copy(diagview(A)) # unglq! destroys the diagonal of A
+                unglq!(driver, A, τ) # Q === A, so no need to rebind Q
+                gaugefix!(lq_householder!, L, Q, Ld)
+            else
+                unglq!(driver, A, τ)
+            end
         else
-            Q = unmlq!(driver, 'R', 'N', A, τ, one!(Q))
-            positive && gaugefix!(lq_householder!, Lf, Q, diagview(A))
+            isempty(Q) || (Q = unmlq!(driver, 'R', 'N', A, τ, one!(Q)))
+            positive && gaugefix!(lq_householder!, L, Q, diagview(A))
         end
     end
 
@@ -220,6 +220,7 @@ function lq_householder!(
             L[i, j] = A[i, j]
         end
     end
+    isempty(Q) && return L, Q
     # build Q
     one!(Q)
     @inbounds for i in minmn:-1:1
@@ -324,13 +325,13 @@ function lq_via_qr!(
     Qt = (A === Q) ? At : similar(Q')
     Lt = similar(L')
     n = size(A, 2)
-    if size(Q) == (n, n)
+    if size(Q) == (n, n) || size(L, 2) == n
         Qt, Lt = qr_full!(At, (Qt, Lt), qr_alg)
     else
         Qt, Lt = qr_compact!(At, (Qt, Lt), qr_alg)
     end
-    adjoint!(Q, Qt)
-    !isempty(L) && adjoint!(L, Lt)
+    isempty(Q) || adjoint!(Q, Qt)
+    isempty(L) || adjoint!(L, Lt)
     return L, Q
 end
 
@@ -348,38 +349,28 @@ end
 function lq_full!(A::AbstractMatrix, LQ, alg::DiagonalAlgorithm)
     check_input(lq_full!, A, LQ, alg)
     L, Q = LQ
-    lq_diagonal!(A, L, Q; alg.kwargs...)
+    _lq_diagonal!(A, L, Q; alg.kwargs...)
     return L, Q
 end
 function lq_compact!(A::AbstractMatrix, LQ, alg::DiagonalAlgorithm)
     check_input(lq_compact!, A, LQ, alg)
     L, Q = LQ
-    lq_diagonal!(A, L, Q; alg.kwargs...)
+    _lq_diagonal!(A, L, Q; alg.kwargs...)
     return L, Q
 end
 function lq_null!(A::AbstractMatrix, N, alg::DiagonalAlgorithm)
     check_input(lq_null!, A, N, alg)
-    return lq_null_diagonal!(A, N; alg.kwargs...)
+    return _lq_null_diagonal!(A, N; alg.kwargs...)
 end
 
-function lq_diagonal!(
+function _lq_diagonal!(
         A::AbstractMatrix, L::AbstractMatrix, Q::AbstractMatrix; positive::Bool = true
     )
-    # note: Ad and Qd might share memory here so order of operations is important
-    Ad = diagview(A)
-    Ld = diagview(L)
-    Qd = diagview(Q)
-    if positive
-        @. Ld = abs(Ad)
-        @. Qd = sign_safe(Ad)
-    else
-        Ld .= Ad
-        one!(Q)
-    end
+    Q, L = _qr_diagonal!(A, Q, L; positive)
     return L, Q
 end
 
-lq_null_diagonal!(A::AbstractMatrix, N; positive::Bool = true) = N
+_lq_null_diagonal!(A::AbstractMatrix, N; positive::Bool = true) = N
 
 # Deprecations
 # ------------

src/implementations/qr.jl

@@ -10,7 +10,7 @@ function check_input(::typeof(qr_full!), A::AbstractMatrix, QR, ::AbstractAlgori
     m, n = size(A)
     Q, R = QR
     @assert Q isa AbstractMatrix && R isa AbstractMatrix
-    @check_size(Q, (m, m))
+    isempty(Q) || @check_size(Q, (m, m))
     @check_scalar(Q, A)
     isempty(R) || @check_size(R, (m, n))
     @check_scalar(R, A)
@@ -21,7 +21,7 @@ function check_input(::typeof(qr_compact!), A::AbstractMatrix, QR, ::AbstractAlg
     minmn = min(m, n)
     Q, R = QR
     @assert Q isa AbstractMatrix && R isa AbstractMatrix
-    @check_size(Q, (m, minmn))
+    isempty(Q) || @check_size(Q, (m, minmn))
     @check_scalar(Q, A)
     isempty(R) || @check_size(R, (minmn, n))
     @check_scalar(R, A)
@@ -40,7 +40,7 @@ function check_input(::typeof(qr_full!), A::AbstractMatrix, (Q, R), alg::Diagona
     m, n = size(A)
     @assert m == n && isdiag(A)
     @assert Q isa Diagonal && R isa Diagonal
-    @check_size(Q, (m, n))
+    isempty(Q) || @check_size(Q, (m, n))
     @check_scalar(Q, A)
     isempty(R) || @check_size(R, (m, n))
     @check_scalar(R, A)
@@ -170,7 +170,7 @@ function qr_householder!(
 
     m, n = size(A)
     minmn = min(m, n)
-    computeR = length(R) > 0
+    computeR = !isempty(R)
     inplaceQ = Q === A
 
     if inplaceQ
@@ -185,14 +185,11 @@ function qr_householder!(
     if blocksize > 1
         # R doubles as workspace for T, so Q is constructed before R is extracted
         nb = min(minmn, blocksize)
-        if computeR # first use R as space for T
-            A, T = geqrt!(driver, A, view(R, 1:nb, 1:minmn))
-        else
-            A, T = geqrt!(driver, A, similar(A, nb, minmn))
-        end
-        Q = gemqrt!(driver, 'L', 'N', A, T, one!(Q))
+        T = computeR ? view(R, 1:nb, 1:minmn) : similar(A, nb, minmn)
+        A, T = geqrt!(driver, A, T)
+        isempty(Q) || (Q = gemqrt!(driver, 'L', 'N', A, T, one!(Q)))
         computeR && _qr_copyR!(R, A)
-        positive && gaugefix!(qr_householder!, Q, computeR ? R : nothing, diagview(A))
+        positive && gaugefix!(qr_householder!, Q, R, diagview(A))
     else
         if pivoted
             A, τ, jpvt = geqp3!(driver, A)
@@ -201,14 +198,18 @@ function qr_householder!(
             A, τ = geqrf!(driver, A)
             computeR && _qr_copyR!(R, A)
         end
-        Rf = computeR ? R : nothing # gaugefix! rescales rows, which commutes with the pivoting
+        # Gauge fixing rescales rows, which commutes with pivoting.
         if inplaceQ
-            Rd = positive ? copy(diagview(A)) : nothing # ungqr! destroys the diagonal of A
-            ungqr!(driver, A, τ) # Q === A, so no need to rebind Q
-            positive && gaugefix!(qr_householder!, Q, Rf, Rd)
+            if positive
+                Rd = copy(diagview(A)) # ungqr! destroys the diagonal of A
+                ungqr!(driver, A, τ) # Q === A, so no need to rebind Q
+                gaugefix!(qr_householder!, Q, R, Rd)
+            else
+                ungqr!(driver, A, τ)
+            end
         else
-            Q = _qr_buildQ!(driver, Q, A, τ, minmn)
-            positive && gaugefix!(qr_householder!, Q, Rf, diagview(A))
+            isempty(Q) || (Q = _qr_buildQ!(driver, Q, A, τ, minmn))
+            positive && gaugefix!(qr_householder!, Q, R, diagview(A))
         end
     end
 
@@ -248,6 +249,7 @@ function qr_householder!(
             R[i, j] = A[i, j]
         end
     end
+    isempty(Q) && return Q, R
     # build Q
     one!(Q)
     @inbounds for j in minmn:-1:1
@@ -330,22 +332,22 @@ end
 function qr_full!(A::AbstractMatrix, QR, alg::DiagonalAlgorithm)
     check_input(qr_full!, A, QR, alg)
     Q, R = QR
-    _diagonal_qr!(A, Q, R; alg.kwargs...)
+    _qr_diagonal!(A, Q, R; alg.kwargs...)
     return Q, R
 end
 function qr_compact!(A::AbstractMatrix, QR, alg::DiagonalAlgorithm)
     check_input(qr_compact!, A, QR, alg)
     Q, R = QR
-    _diagonal_qr!(A, Q, R; alg.kwargs...)
+    _qr_diagonal!(A, Q, R; alg.kwargs...)
     return Q, R
 end
 function qr_null!(A::AbstractMatrix, N, alg::DiagonalAlgorithm)
     check_input(qr_null!, A, N, alg)
-    _diagonal_qr_null!(A, N; alg.kwargs...)
+    _qr_null_diagonal!(A, N; alg.kwargs...)
     return N
 end
 
-function _diagonal_qr!(
+function _qr_diagonal!(
         A::AbstractMatrix, Q::AbstractMatrix, R::AbstractMatrix; positive::Bool = true
     )
     # note: Ad and Qd might share memory here so order of operations is important
@@ -354,15 +356,15 @@ function _diagonal_qr!(
     Rd = diagview(R)
     if positive
         @. Rd = abs(Ad)
-        @. Qd = sign_safe(Ad)
+        isempty(Q) || (@. Qd = sign_safe(Ad))
     else
         Rd .= Ad
         one!(Q)
     end
     return Q, R
 end
 
-_diagonal_qr_null!(A::AbstractMatrix, N; positive::Bool = true) = N
+_qr_null_diagonal!(A::AbstractMatrix, N; positive::Bool = true) = N
 
 # Deprecations
 # ------------

src/interface/lq.jl

@@ -10,6 +10,9 @@ Compute the full LQ decomposition of the rectangular matrix `A`, such that `A =
 where `Q` is a square unitary matrix with the same number of rows as `A` and `L` is a
 lower triangular matrix with the same size as `A`.
 
+Pass an empty `Q` output to `lq_full!` to compute only `L` without constructing `Q`.
+Conversely, an empty `L` output requests only `Q`.
+
 !!! note
     The bang method `lq_full!` optionally accepts the output structure and
     possibly destroys the input matrix `A`. Always use the return value of the function
@@ -29,6 +32,9 @@ Compute the compact LQ decomposition of the rectangular matrix `A` of size `(m,n
 such that `A = L * Q` where the matrix `Q` of size `(min(m,n), n)` has orthogonal rows
 spanning the image of `A'`, and the matrix `L` of size `(m, min(m,n))` is lower triangular.
 
+Pass an empty `Q` output to `lq_compact!` to compute only `L` without constructing `Q`.
+Conversely, an empty `L` output requests only `Q`.
+
 !!! note
     The bang method `lq_compact!` optionally accepts the output structure and
     possibly destroys the input matrix `A`. Always use the return value of the function

src/interface/orthnull.jl

@@ -76,6 +76,10 @@ deduce the decomposition kind with [`left_orth_alg(alg)`](@ref left_orth_alg).
     destroys the input matrix `A`. Always use the return value of the function as it may
     not always be possible to use the provided `CV` as output.
 
+For QR-based decompositions, pass an empty `V` output to `left_orth!` to compute only `C` without constructing `V`.
+Conversely, an empty `C` output requests only `V`.
+The returned basis matrix remains empty, and its construction is skipped.
+
 See also [`right_orth(!)`](@ref right_orth), [`left_null(!)`](@ref left_null) and
 [`right_null(!)`](@ref right_null).
 """
@@ -157,6 +161,9 @@ deduce the decomposition kind with [`right_orth_alg`](@ref).
     destroys the input matrix `A`. Always use the return value of the function as it may not
     always be possible to use the provided `CVᴴ` as output.
 
+For LQ-based decompositions, pass an empty `Vᴴ` output to `right_orth!` to compute only `C` without constructing `Vᴴ`.
+Conversely, an empty `C` output requests only `Vᴴ`.
+
 See also [`left_orth(!)`](@ref left_orth), [`left_null(!)`](@ref left_null) and
 [`right_null(!)`](@ref right_null).
 """

src/interface/qr.jl

@@ -10,6 +10,9 @@ Compute the full QR decomposition of the rectangular matrix `A`, such that `A =
 where `Q` is a square unitary matrix with the same number of rows as `A` and `R` is an
 upper triangular matrix with the same size as `A`.
 
+Pass an empty `Q` output to `qr_full!` to compute only `R` without constructing `Q`.
+Conversely, an empty `R` output requests only `Q`.
+
 !!! note
     The bang method `qr_full!` optionally accepts the output structure and
     possibly destroys the input matrix `A`. Always use the return value of the function
@@ -30,6 +33,9 @@ such that `A = Q * R` where the isometric matrix `Q` of size `(m, min(m,n))` has
 orthogonal columns spanning the image of `A`, and the matrix `R` of size `(min(m,n), n)`
 is upper triangular.
 
+Pass an empty `Q` output to `qr_compact!` to compute only `R` without constructing `Q`.
+Conversely, an empty `R` output requests only `Q`.
+
 !!! note
     The bang method `qr_compact!` optionally accepts the output structure and
     possibly destroys the input matrix `A`. Always use the return value of the function