Merge pull request #231 from JuliaArrays/matrix-factorizations

Matrix factorizations

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,
+             lu, 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
@@ -93,6 +94,8 @@ include("eigen.jl")
 include("expm.jl")
 include("sqrtm.jl")
 include("cholesky.jl")
+include("svd.jl")
+include("lu.jl")
 include("deque.jl")
 include("io.jl")
 

src/linalg.jl

@@ -323,3 +323,4 @@ 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))
+

src/lu.jl

@@ -0,0 +1,30 @@
+# LU decomposition
+function lu(A::StaticMatrix, pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{true})
+    L,U,p = _lu(Size(A), A, pivot)
+    (L,U,p)
+end
+
+# For the square version, return explicit lower and upper triangular matrices.
+# We would do this for the rectangular case too, but Base doesn't support that.
+function lu(A::StaticMatrix{N,N}, pivot::Union{Type{Val{false}},Type{Val{true}}}=Val{true}) where {N}
+    L,U,p = _lu(Size(A), A, pivot)
+    (LowerTriangular(L), UpperTriangular(U), p)
+end
+
+
+@inline function _lu(::Size{S}, A::StaticMatrix, pivot) where {S}
+    # For now, just call through to Base.
+    # TODO: statically sized LU without allocations!
+    f = lufact(Matrix(A), pivot)
+    T = eltype(A)
+    # Trick to get the output eltype - can't rely on the result of f[:L] as
+    # it's not type inferrable.
+    T2 = typeof((one(T)*zero(T) + zero(T))/one(T))
+    L = similar_type(A, T2, Size(Size(A)[1], diagsize(A)))(f[:L])
+    U = similar_type(A, T2, Size(diagsize(A), Size(A)[2]))(f[:U])
+    p = similar_type(A, Int, Size(Size(A)[1]))(f[:p])
+    (L,U,p)
+end
+
+# Base.lufact() interface is fairly inherently type unstable.  Punt on
+# implementing that, for now...

src/svd.jl

@@ -0,0 +1,37 @@
+# Singular Value Decomposition
+
+# 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")
+
+function svdvals(A::StaticMatrix)
+    sv = svdvals(Matrix(A))
+    # We should be using `T2=eltype(sv)`, but it's not inferrable for complex
+    # eltypes.  See https://github.com/JuliaLang/julia/pull/22443
+    T = eltype(A)
+    T2 = promote_type(Float32, real(typeof(one(T)/norm(one(T)))))
+    similar_type(A, T2, 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, eltype(f.S), 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

src/traits.jl

@@ -123,3 +123,7 @@ end
 @noinline function _throw_size_mismatch(as...)
     throw(DimensionMismatch("Sizes $(map(_size, as)) of input arrays do not match"))
 end
+
+# Return the "diagonal size" of a matrix - the minimum of the two dimensions
+diagsize(A::StaticMatrix) = diagsize(Size(A))
+@pure diagsize(::Size{S}) where {S} = min(S...)

test/linalg.jl

@@ -1,3 +1,5 @@
+using StaticArrays, Base.Test
+
 @testset "Linear algebra" begin
 
     @testset "SVector as a (mathematical) vector space" begin
@@ -138,7 +140,7 @@
         @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/lu.jl

@@ -0,0 +1,21 @@
+using StaticArrays, Base.Test
+
+@testset "LU decomposition" begin
+    # Square case
+    m22 = @SMatrix [1 2; 3 4]
+    @test @inferred(lu(m22)) isa Tuple{LowerTriangular{Float64,SMatrix{2,2,Float64,4}}, UpperTriangular{Float64,SMatrix{2,2,Float64,4}}, SVector{2,Int}}
+    @test lu(m22)[1]::LowerTriangular{<:Any,<:StaticMatrix} ≊ lu(Matrix(m22))[1]
+    @test lu(m22)[2]::UpperTriangular{<:Any,<:StaticMatrix} ≊ lu(Matrix(m22))[2]
+    @test lu(m22)[3]::StaticVector ≊ lu(Matrix(m22))[3]
+
+    # Rectangular case
+    m23 = @SMatrix Float64[3 9 4; 6 6 2]
+    @test @inferred(lu(m23)) isa Tuple{SMatrix{2,2,Float64,4}, SMatrix{2,3,Float64,6}, SVector{2,Int}}
+    @test lu(m23)[1] ≊ lu(Matrix(m23))[1]
+    @test lu(m23)[2] ≊ lu(Matrix(m23))[2]
+    @test lu(m23)[3] ≊ lu(Matrix(m23))[3]
+
+    @test lu(m23')[1] ≊ lu(Matrix(m23'))[1]
+    @test lu(m23')[2] ≊ lu(Matrix(m23'))[2]
+    @test lu(m23')[3] ≊ lu(Matrix(m23'))[3]
+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")
@@ -31,6 +33,7 @@ using Base.Test
     include("chol.jl")
     include("deque.jl")
     include("io.jl")
+    include("svd.jl")
 
     include("fixed_size_arrays.jl")
 end

test/svd.jl

@@ -0,0 +1,47 @@
+using StaticArrays, Base.Test
+
+@testset "SVD factorization" begin
+    m3 = @SMatrix Float64[3 9 4; 6 6 2; 3 7 9]
+    m3c = Complex128.(m3)
+    m23 = @SMatrix Float64[3 9 4; 6 6 2]
+
+    @testset "svd" begin
+        @testinf svdvals(@SMatrix [2 0; 0 0])::StaticVector ≊ [2, 0]
+        @testinf svdvals((@SMatrix [2 -2; 1 1]) / sqrt(2)) ≊ [2, 1]
+
+        @testinf svdvals(m3) ≊ svdvals(Matrix(m3))
+        @testinf svdvals(m3c) isa SVector{3,Float64}
+
+        @testinf svd(m3)[1]::StaticMatrix ≊ svd(Matrix(m3))[1]
+        @testinf svd(m3)[2]::StaticVector ≊ svd(Matrix(m3))[2]
+        @testinf svd(m3)[3]::StaticMatrix ≊ 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, 0.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)
+
+        @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]
+
+        @testinf svdfact(m3c).U  ≊ svdfact(Matrix(m3c))[:U]
+        @testinf svdfact(m3c).S  ≊ svdfact(Matrix(m3c))[:S]
+        @testinf svdfact(m3c).Vt ≊ svdfact(Matrix(m3c))[:Vt]
+
+        @testinf svdfact(m3c).U  isa SMatrix{3,3,Complex128}
+        @testinf svdfact(m3c).S  isa SVector{3,Float64}
+        @testinf svdfact(m3c).Vt isa SMatrix{3,3,Complex128}
+    end
+end

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