Add solve methods for rectangular A

Author: cgarling

src/StaticArrays.jl

@@ -120,6 +120,7 @@ include("matrix_multiply.jl")
 include("lu.jl")
 include("det.jl")
 include("inv.jl")
+include("qr.jl")
 include("solve.jl")
 include("eigen.jl")
 include("expm.jl")
@@ -128,7 +129,6 @@ include("lyap.jl")
 include("triangular.jl")
 include("cholesky.jl")
 include("svd.jl")
-include("qr.jl")
 include("deque.jl")
 include("flatten.jl")
 include("io.jl")

src/solve.jl

@@ -1,4 +1,5 @@
 @inline (\)(a::StaticMatrix, b::StaticVecOrMat) = _solve(Size(a), Size(b), a, b)
+@inline (\)(Q::QR, b::StaticVecOrMat) = Q.R \ (Q.Q' * b)
 
 @inline function _solve(::Size{(1,1)}, ::Size{(1,)}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticVector{<:Any, Tb}) where {Ta, Tb}
     @inbounds return similar_type(b, typeof(a[1] \ b[1]))(a[1] \ b[1])
@@ -55,15 +56,30 @@ end
             throw(DimensionMismatch("Left and right hand side first dimensions do not match in backdivide (got sizes $Sa and $Sb)"))
         end
     end
-    if prod(Sa) ≤ 14*14 && Sa[1] == Sa[2]
+    if prod(Sa) ≤ 14*14
         # TODO: Consider triangular special cases as in Base?
-        quote
-            @_inline_meta
-            LUp = lu(a)
-            LUp.U \ (LUp.L \ $(length(Sb) > 1 ? :(b[LUp.p,:]) : :(b[LUp.p])))
+        if Sa[1] == Sa[2]
+            quote
+                @_inline_meta
+                LUp = lu(a)
+                LUp.U \ (LUp.L \ $(length(Sb) > 1 ? :(b[LUp.p,:]) : :(b[LUp.p])))
+            end
+        else
+
+            quote
+                @_inline_meta
+                q = qr(a)
+                y = q.Q' * b
+                if Sa[1] > Sa[2]
+                    R₁ = SMatrix{Sa[2], Sa[2]}(q.R[SOneTo(Sa[2]), SOneTo(Sa[2])])
+                    # return inv(R₁) * y
+                    return R₁ \ y
+                else
+                    return q.R' * ((q.R * q.R') \ y)
+                    # return pinv(q.R) * y
+                end
+            end
         end
-        # TODO: Could also use static QR here if `a` is nonsquare.
-        # Requires that we implement \(::StaticArrays.QR,::StaticVecOrMat)
     else
         # Fall back to LinearAlgebra, but carry across the statically known size.
         outsize = length(Sb) == 1 ? Size(Sa[2]) : Size(Sa[2],Sb[end])

test/solve.jl

@@ -21,11 +21,17 @@ using StaticArrays, Test, LinearAlgebra
     # So try all of these
     @testset "Mixed static/dynamic" begin
         v = @SVector([0.2,0.3])
+        # Square matrices
         for m in (@SMatrix([1.0 0; 0 1.0]), @SMatrix([1.0 0; 1.0 1.0]),
                   @SMatrix([1.0 1.0; 0 1.0]), @SMatrix([1.0 0.5; 0.25 1.0]))
-            # TODO: include @SMatrix([1.0 0.0 0.0; 1.0 2.0 0.5]), need qr methods
             @test m \ v ≈ Array(m) \ v ≈ m \ Array(v) ≈ Array(m) \ Array(v)
         end
+        # Rectangular matrices
+        for m in (@SMatrix([1.0 0.0 0.0; 1.0 2.0 0.5]), @SMatrix([1.0 2.0 0.5; 0.0 0.0 1.0]),
+                  @SMatrix([0.0 0.0 1.0; 1.0 2.0 0.5]), @SMatrix([1.0 2.0 0.5; 1.0 0.0 0.0]))
+            @test m \ v ≈ Array(m) \ v ≈ Array(m) \ Array(v)
+            @test_throws MethodError m \ Array(v) # TODO: requires adjoint(::QR) method
+        end
     end
 end