Static wrappers for svdvals(), svd() and svdfact()

A bunch of wrappers to ensure that the shape of static matrices is
preserved when using various SVD factorization related functions.

Author: c42f

src/StaticArrays.jl

@@ -9,6 +9,7 @@ import Base: getindex, setindex!, size, similar, vec, show,
              mapreduce, broadcast, broadcast!, conj, transpose, ctranspose,
              hcat, vcat, ones, zeros, eye, one, cross, vecdot, reshape, fill,
              fill!, det, inv, eig, eigvals, expm, sqrtm, trace, vecnorm, norm, dot, diagm, diag,
+             svd, svdvals!, svdfact,
              sum, diff, prod, count, any, all, minimum,
              maximum, extrema, mean, copy, rand, randn, randexp, rand!, randn!,
              randexp!, normalize, normalize!, read, read!, write

src/linalg.jl

@@ -323,3 +323,57 @@ end
 
 @inline Base.LinAlg.Symmetric(A::StaticMatrix, uplo::Char='U') = (Base.LinAlg.checksquare(A);Symmetric{eltype(A),typeof(A)}(A, uplo))
 @inline Base.LinAlg.Hermitian(A::StaticMatrix, uplo::Char='U') = (Base.LinAlg.checksquare(A);Hermitian{eltype(A),typeof(A)}(A, uplo))
+
+
+#-------------------------------------------------------------------------------
+# Wrappers for various matrix decompositions which correctly propagate the
+# StaticMatrix type in the returned factorization.
+
+
+#--------------------------------------------------
+# lu
+# TODO
+
+
+#--------------------------------------------------
+# qr
+# TODO
+
+#--------------------------------------------------
+# SVD
+# We need our own SVD factorization struct, as Base.LinAlg.SVD assumes
+# Base.Vector for `S`, and that the `U` and `Vt` have the same
+struct SVD{T,TU,TS,TVt} <: Factorization{T}
+    U::TU
+    S::TS
+    Vt::TVt
+end
+SVD(U::AbstractArray{T}, S::AbstractVector, Vt::AbstractArray{T}) where {T} = SVD{T,typeof(U),typeof(S),typeof(Vt)}(U, S, Vt)
+
+getindex(::SVD, ::Symbol) = error("In order to avoid type instability, StaticArrays.SVD doesn't support indexing the output of svdfact with a symbol.  Instead, you can access the fields of the factorization directly as f.U, f.S, and f.Vt")
+
+# Return the "diagonal size" of a matrix - the minimum of the two dimensions
+@generated function diagsize(A::StaticMatrix{N,M}) where {N,M}
+    :($(min(N,M)))
+end
+
+function svdvals!(A::StaticMatrix)
+    sv = svdvals!(Matrix(A))
+    similar_type(A, eltype(sv), Size(diagsize(A)))(sv)
+end
+
+function svdfact(A::StaticMatrix)
+    # "Thin" SVD only for now.
+    f = svdfact(Matrix(A))
+    U = similar_type(A, eltype(f.U), Size(Size(A)[1], diagsize(A)))(f.U)
+    S = similar_type(A, Size(diagsize(A)))(f.S)
+    Vt = similar_type(A, eltype(f.Vt), Size(diagsize(A), Size(A)[2]))(f.Vt)
+    SVD(U,S,Vt)
+end
+
+function svd(A::StaticMatrix)
+    # Need our own version of `svd()`, as `Base` passes the `thin` argument
+    # which makes the resulting dimensions uninferrable.
+    f = svdfact(A)
+    (f.U, f.S, f.Vt')
+end

test/linalg.jl

@@ -1,3 +1,5 @@
+using StaticArrays, Base.Test
+
 @testset "Linear algebra" begin
 
     @testset "SVector as a (mathematical) vector space" begin
@@ -138,11 +140,44 @@
         @test trace(@SMatrix [1.0 2.0; 3.0 4.0]) === 5.0
         @test_throws DimensionMismatch trace(@SMatrix rand(5,4))
     end
-    
+
     @testset "size zero" begin
         @test vecdot(SVector{0, Float64}(()), SVector{0, Float64}(())) === 0.
         @test vecnorm(SVector{0, Float64}(())) === 0.
         @test vecnorm(SVector{0, Float64}(()), 1) === 0.
         @test trace(SMatrix{0,0,Float64}(())) === 0.
     end
+
+    @testset "svd" begin
+        @testinf svdvals(@SMatrix [2 0; 0 0]) ≊ [2, 0]
+        @testinf svdvals((@SMatrix [2 -2; 1 1]) / sqrt(2)) ≊ [2, 1]
+
+        m3 = @SMatrix Float64[3 9 4; 6 6 2; 3 7 9]
+        @testinf svdvals(m3) ≈ svdvals(Matrix(m3))
+
+        @testinf svd(m3)[1] ≈ svd(Matrix(m3))[1]
+        @testinf svd(m3)[2] ≈ svd(Matrix(m3))[2]
+        @testinf svd(m3)[3] ≈ svd(Matrix(m3))[3]
+    end
+
+    @testset "svdfact" begin
+        @test_throws ErrorException svdfact(@SMatrix [1 0; 0 1])[:U]
+
+        @testinf svdfact(@SMatrix [2 0; 0 0]).U === eye(SMatrix{2,2})
+        @testinf svdfact(@SMatrix [2 0; 0 0]).S === SVector(2, 0)
+        @testinf svdfact(@SMatrix [2 0; 0 0]).Vt === eye(SMatrix{2,2})
+
+        @testinf svdfact((@SMatrix [2 -2; 1 1]) / sqrt(2)).U ≊ [-1 0; 0 1]
+        @testinf svdfact((@SMatrix [2 -2; 1 1]) / sqrt(2)).S ≊ [2, 1]
+        @testinf svdfact((@SMatrix [2 -2; 1 1]) / sqrt(2)).Vt ≊ [-1 1; 1 1]/sqrt(2)
+
+        m23 = @SMatrix Float64[3 9 4; 6 6 2]
+        @testinf svdfact(m23).U  ≊ svdfact(Matrix(m23))[:U]
+        @testinf svdfact(m23).S  ≊ svdfact(Matrix(m23))[:S]
+        @testinf svdfact(m23).Vt ≊ svdfact(Matrix(m23))[:Vt]
+
+        @testinf svdfact(m23').U  ≊ svdfact(Matrix(m23'))[:U]
+        @testinf svdfact(m23').S  ≊ svdfact(Matrix(m23'))[:S]
+        @testinf svdfact(m23').Vt ≊ svdfact(Matrix(m23'))[:Vt]
+    end
 end

test/runtests.jl

@@ -1,6 +1,8 @@
 using StaticArrays
 using Base.Test
 
+include("testutil.jl")
+
 @testset "StaticArrays" begin
     include("SVector.jl")
     include("MVector.jl")

test/testutil.jl

@@ -0,0 +1,36 @@
+using Base.Test, StaticArrays
+
+"""
+    x ≊ y
+
+Inexact equality comparison. Like `≈` this calls `isapprox`, but with a
+tighter tolerance of `rtol=10*eps()`.  Input with "\\approxeq".
+"""
+≊(x,y) = isapprox(x, y, rtol=10*eps())
+
+"""
+    @testinf a op b
+
+Test that the type of the first argument `a` is inferred, and that `a op b` is
+true.  For example, the following are equivalent:
+
+    @testinf SVector(1,2) + SVector(1,2) == SVector(2,4)
+    @test @inferred(SVector(1,2) + SVector(1,2)) == SVector(2,4)
+"""
+macro testinf(ex)
+    @assert ex.head == :call
+    infarg = ex.args[2]
+    if !(infarg isa Expr) || infarg.head != :call
+        # Workaround for an oddity in @inferred
+        infarg = :(identity($infarg))
+    end
+    ex.args[2] = :(@inferred($infarg))
+    esc(:(@test $ex))
+end
+
+@testset "@testinf" begin
+    @testinf [1,2] == [1,2]
+    x = [1,2]
+    @testinf x == [1,2]
+    @testinf (@SVector [1,2]) == (@SVector [1,2])
+end