nonsquare (underdetermined) LQ solves (#34350)

* nonsquare (underdetermined) LQ solves

* fix NEWS item

* fix NEWS item

* narrow method signature to eliminate ambiguity

* fix docs -- under/overdetermined terms were swapped

* Update stdlib/LinearAlgebra/test/lq.jl

Co-Authored-By: Stefan Karpinski <stefan@karpinski.org>

Co-authored-by: Stefan Karpinski <stefan@karpinski.org>

Author: 2 people

NEWS.md

@@ -36,6 +36,7 @@ Standard library changes
 #### LinearAlgebra
 * The BLAS submodule now supports the level-2 BLAS subroutine `hpmv!` ([#34211]).
 * `normalize` now supports multidimensional arrays ([#34239])
+* `lq` factorizations can now be used to compute the minimum-norm solution to under-determined systems ([#34350]).
 
 #### Markdown
 

stdlib/LinearAlgebra/src/lq.jl

@@ -80,7 +80,9 @@ orthogonal/unitary component via `S.Q`, such that `A ≈ S.L*S.Q`.
 
 Iterating the decomposition produces the components `S.L` and `S.Q`.
 
-The LQ decomposition is the QR decomposition of `transpose(A)`.
+The LQ decomposition is the QR decomposition of `transpose(A)`, and it is useful
+in order to compute the minimum-norm solution `lq(A) \\ b` to an underdetermined
+system of equations (`A` has more columns than rows, but has full row rank).
 
 # Examples
 ```jldoctest
@@ -174,11 +176,13 @@ end
 
 
 ## Multiplication by LQ
-lmul!(A::LQ, B::StridedVecOrMat) =
-    lmul!(LowerTriangular(A.L), lmul!(A.Q, B))
+function lmul!(A::LQ, B::StridedVecOrMat)
+    lmul!(LowerTriangular(A.L), view(lmul!(A.Q, B), 1:size(A,1), axes(B,2)))
+    return B
+end
 function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
     TAB = promote_type(TA, TB)
-    lmul!(Factorization{TAB}(A), copy_oftype(B, TAB))
+    _cut_B(lmul!(Factorization{TAB}(A), copy_oftype(B, TAB)), 1:size(A,1))
 end
 
 ## Multiplication by Q
@@ -303,36 +307,33 @@ _rightappdimmismatch(rowsorcols) =
         "(the factorization's originating matrix's number of rows)")))
 
 
-function (\)(A::LQ{TA}, b::StridedVector{Tb}) where {TA,Tb}
-    S = promote_type(TA,Tb)
-    m = checksquare(A)
-    m == length(b) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has length $(length(b))"))
-    AA = Factorization{S}(A)
-    x = ldiv!(AA, copy_oftype(b, S))
-    return x
-end
-function (\)(A::LQ{TA},B::StridedMatrix{TB}) where {TA,TB}
+function (\)(A::LQ{TA},B::StridedVecOrMat{TB}) where {TA,TB}
     S = promote_type(TA,TB)
-    m = checksquare(A)
-    m == size(B,1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(B,1)) rows"))
+    m, n = size(A)
+    m ≤ n || throw(DimensionMismatch("LQ solver does not support overdetermined systems (more rows than columns)"))
+    m == size(B,1) || throw(DimensionMismatch("Both inputs should have the same number of rows"))
     AA = Factorization{S}(A)
-    X = ldiv!(AA, copy_oftype(B, S))
-    return X
+    X = _zeros(S, B, n)
+    X[1:size(B, 1), :] = B
+    return ldiv!(AA, X)
 end
 # With a real lhs and complex rhs with the same precision, we can reinterpret
 # the complex rhs as a real rhs with twice the number of columns
 function (\)(F::LQ{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal
     require_one_based_indexing(B)
-    c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2))
-    x = ldiv!(F, c2r)
-    return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(x, div(length(x), 2), 2))))),
+    X = zeros(T, size(F,2), 2*size(B,2))
+    X[1:size(B,1), 1:size(B,2)] .= real.(B)
+    X[1:size(B,1), size(B,2)+1:size(X,2)] .= imag.(B)
+    ldiv!(F, X)
+    return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(X, div(length(X), 2), 2))))),
                            isa(B, AbstractVector) ? (size(F,2),) : (size(F,2), size(B,2)))
 end
 
 
 function ldiv!(A::LQ{T}, B::StridedVecOrMat{T}) where T
-    lmul!(adjoint(A.Q), ldiv!(LowerTriangular(A.L),B))
-    return B
+    require_one_based_indexing(B)
+    ldiv!(LowerTriangular(A.L), view(B, 1:size(A,1), axes(B,2)))
+    return lmul!(adjoint(A.Q), B)
 end
 
 

stdlib/LinearAlgebra/test/lq.jl

@@ -5,48 +5,42 @@ module TestLQ
 using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!
 
-n = 10
-
-# Split n into 2 parts for tests needing two matrices
-n1 = div(n, 2)
-n2 = 2*n1
+m = 10
 
 Random.seed!(1234321)
 
-areal = randn(n,n)/2
-aimg  = randn(n,n)/2
-a2real = randn(n,n)/2
-a2img  = randn(n,n)/2
-breal = randn(n,2)/2
-bimg  = randn(n,2)/2
+asquare = randn(ComplexF64, m, m) / 2
+awide = randn(ComplexF64, m, m+3) / 2
+bcomplex = randn(ComplexF64, m, 2) / 2
 
 # helper functions to unambiguously recover explicit forms of an LQPackedQ
 squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); lmul!(Q, Matrix{eltype(Q)}(I, n, n)))
 rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
 
-@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
-    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
-    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
-    asym = a' + a                 # symmetric indefinite
-    apd  = a' * a                 # symmetric positive-definite
+@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64), n in (m, size(awide, 2))
+    adata = m == n ? asquare : awide
+    a = convert(Matrix{eltya}, eltya <: Complex ? adata : real(adata))
     ε = εa = eps(abs(float(one(eltya))))
+    n1 = n ÷ 2
+
+    α = rand(eltya)
+    aα = fill(α,1,1)
+    @test lq(α).L*lq(α).Q ≈ lq(aα).L*lq(aα).Q
+    @test abs(lq(α).Q[1,1]) ≈ one(eltya)
 
     @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
-        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
+        b = eltyb == Int ? rand(1:5, m, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? bcomplex : real(bcomplex))
         εb = eps(abs(float(one(eltyb))))
         ε = max(εa,εb)
 
-        α = rand(eltya)
-        aα = fill(α,1,1)
-        @test lq(α).L*lq(α).Q ≈ lq(aα).L*lq(aα).Q
-        @test abs(lq(α).Q[1,1]) ≈ one(eltya)
         tab = promote_type(eltya,eltyb)
 
-        for i = 1:2
-            let a = i == 1 ? a : view(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : view(b, 1:n - 1), n = i == 1 ? n : n - 1
+        @testset for isview in (false,true)
+            let a = isview ? view(a, 1:m - 1, 1:n - 1) : a, b = isview ? view(b, 1:m - 1) : b, m = m - isview, n = n - isview
                 lqa   = lq(a)
+                x = lqa\b
                 l,q   = lqa.L, lqa.Q
-                qra   = qr(a)
+                qra   = qr(a, Val(true))
                 @testset "Basic ops" begin
                     @test size(lqa,1) == size(a,1)
                     @test size(lqa,3) == 1
@@ -69,30 +63,31 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
                     @test Matrix{eltya}(q) isa Matrix{eltya}
                 end
                 @testset "Binary ops" begin
-                    @test a*(lqa\b) ≈ b atol=3000ε
-                    @test lqa*b ≈ qra.Q*qra.R*b atol=3000ε
-                    @test (sq = size(q.factors, 2); *(Matrix{eltyb}(I, sq, sq), adjoint(q))*squareQ(q)) ≈ Matrix(I, n, n) atol=5000ε
+                    @test a*x ≈ b rtol=3000ε
+                    @test x ≈ qra \ b rtol=3000ε
+                    @test lqa*x ≈ a*x rtol=3000ε
+                    @test (sq = size(q.factors, 2); *(Matrix{eltyb}(I, sq, sq), adjoint(q))*squareQ(q)) ≈ Matrix(I, n, n) rtol=5000ε
                     if eltya != Int
                         @test Matrix{eltyb}(I, n, n)*q ≈ convert(AbstractMatrix{tab},q)
                     end
-                    @test q*b ≈ squareQ(q)*b atol=100ε
-                    @test transpose(q)*b ≈ transpose(squareQ(q))*b atol=100ε
-                    @test q'*b ≈ squareQ(q)'*b atol=100ε
-                    @test a*q ≈ a*squareQ(q) atol=100ε
-                    @test a*transpose(q) ≈ a*transpose(squareQ(q)) atol=100ε
-                    @test a*q' ≈ a*squareQ(q)' atol=100ε
-                    @test a'*q ≈ a'*squareQ(q) atol=100ε
-                    @test a'*q' ≈ a'*squareQ(q)' atol=100ε
-                    @test_throws DimensionMismatch q*b[1:n1 + 1]
-                    @test_throws DimensionMismatch adjoint(q) * Matrix{eltya}(undef,n+2,n+2)
-                    @test_throws DimensionMismatch Matrix{eltyb}(undef,n+2,n+2)*q
+                    @test q*x ≈ squareQ(q)*x rtol=100ε
+                    @test transpose(q)*x ≈ transpose(squareQ(q))*x rtol=100ε
+                    @test q'*x ≈ squareQ(q)'*x rtol=100ε
+                    @test a*q ≈ a*squareQ(q) rtol=100ε
+                    @test a*transpose(q) ≈ a*transpose(squareQ(q)) rtol=100ε
+                    @test a*q' ≈ a*squareQ(q)' rtol=100ε
+                    @test q*a'≈ squareQ(q)*a' rtol=100ε
+                    @test q'*a' ≈ squareQ(q)'*a' rtol=100ε
+                    @test_throws DimensionMismatch q*x[1:n1 + 1]
+                    @test_throws DimensionMismatch adjoint(q) * Matrix{eltya}(undef,m+2,m+2)
+                    @test_throws DimensionMismatch Matrix{eltyb}(undef,m+2,m+2)*q
                     if isa(a, DenseArray) && isa(b, DenseArray)
                         # use this to test 2nd branch in mult code
                         pad_a = vcat(I, a)
-                        pad_b = hcat(I, b)
-                        @test pad_a*q ≈ pad_a*squareQ(q) atol=100ε
-                        @test transpose(q)*pad_b ≈ transpose(squareQ(q))*pad_b atol=100ε
-                        @test q'*pad_b ≈ squareQ(q)'*pad_b atol=100ε
+                        pad_x = hcat(I, x)
+                        @test pad_a*q ≈ pad_a*squareQ(q) rtol=100ε
+                        @test transpose(q)*pad_x ≈ transpose(squareQ(q))*pad_x rtol=100ε
+                        @test q'*pad_x ≈ squareQ(q)'*pad_x rtol=100ε
                     end
                 end
             end