No svd method for Hermitian matrix

[lstagner]

There is currently no method for svd with a hermitian matrix.

julia> using LinearAlgebra

julia> x = Hermitian(rand(5,5))
5×5 Hermitian{Float64,Array{Float64,2}}:
 0.0584939  0.45097   0.73698    0.945566  0.332148
 0.45097    0.91507   0.272532   0.514991  0.946693
 0.73698    0.272532  0.0604401  0.635022  0.906136
 0.945566   0.514991  0.635022   0.352346  0.5379
 0.332148   0.946693  0.906136   0.5379    0.484081

julia> svd(x)
ERROR: MethodError: no method matching svd(::Hermitian{Float64,Array{Float64,2}})
Closest candidates are:
  svd(::BitArray{2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/bitarray.jl:90
  svd(::LinearAlgebra.AbstractTriangular) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:2471
  svd(::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, ::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/svd.jl:281
  ...
Stacktrace:
 [1] top-level scope at none:0

julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 12
  JULIA_EDITOR = nvim

[dkarrasch]

For Hermitian and Symmetric matrices, the SVD is just the eigendecomposition, so probably one should write a specialized function that calls eigen, and return an SVD object composed of the Eigen result. But I'm not sure if there are other, more performant ways to do it.

[dkarrasch]

I think I have working code for that. Will make a PR for discussion/optimization later.

[dkarrasch]

Well, before making a PR, let me drop a few lines and ask for comments. I have an svd1 that leverages eigen and avoids materializing the Hermitian matrix, and an svd2 that does simply that.

using LinearAlgebra

function svd1(A::Union{Hermitian,Symmetric{<:Real}}, full::Bool=false)
    s = eigen(A)
    vals = s.values
    σ = map(x -> x >= 0 ? 1 : -1, vals)
    vals .*= σ
    I = sortperm(vals; rev=true)
    vecs = s.vectors[:,I]
    return SVD(vecs, vals[I], (vecs*Diagonal(σ[I]))')
end

function svd2(At::Union{Hermitian,Symmetric{<:Real}}, full::Bool=false)
    A = Matrix(At)
    return svd!(A)
end

Especially for large matrices, the first one seems advantageous, but even for smaller sizes svd1 seems never worse than svd2.

using Test, BenchmarkTools
for A in (Hermitian(rand(ComplexF64,5,5)), Hermitian(zeros(5,5)), Hermitian(ones(5,5)))
    MA = Matrix(A)
    F = svd1(A)
    FM = svd(MA)
    @test A ≈ F.U * Diagonal(F.S) * F.Vt
    @test all(x -> x >= 0, F.S)
    @test F.U'F.U ≈ Matrix(I, 5, 5)
    @test F.Vt*F.V ≈ Matrix(I, 5, 5)
    @test F.S ≈ FM.S
end

At = Hermitian(rand(ComplexF64,1000,1000));
@btime svd1($At) # 377.069 ms (37 allocations: 92.38 MiB)
@btime svd2($At) #  690.571 ms (15 allocations: 85.05 MiB)

I couldn't improve svdvals, which does not throw a MethodError anyways.

[StefanKarpinski]

Any thoughts on this, @kshyatt, @andreasnoack, @stevengj?

[stevengj]

Seems fine to me.

[oxinabox]

svd1 looks like it can be optimised to cut out some of those allocations.
e.g. using abs instead of product with the sign?

[dkarrasch]

Replacing vals .*= σ by vals .= abs.(vals) doesn't change the number of allocations reported by @btime. I assume that this line doesn't allocate either way. I have another version, which uses lots of inplace permutation.

function svd(A::RealHermSymComplexHerm, full::Bool=false)
    vals, vecs = eigen(A)
    σ = map(x -> x >= 0 ? 1 : -1, vals)
    vals .= abs.(vals)
    I = sortperm(vals; rev=true)
    permute!(vals, I)
    Σ = Diagonal(σ[I])
    Base.permutecols!!(vecs, I)
    return SVD(vecs, vals, (vecs*Σ)')
end

I'm not sure it's reproducibly faster, but it definitely allocates less, basically only σ, I, and σ[I], besides what needs to be allocated for the output.

[oxinabox]

I thought it should be possible to not have to allocate σ at all, by also using abs in place of (vecs*Diagonal[σ])' ?

[stevengj]

@oxinabox, you need the signs to scale the eigenvectors. You could do something like:

I = sortperm(vals; by=abs, rev=true)
permute!(vals, I)
Base.permutecols!!(vecs, I)
V = copy(vecs)
@inbounds for i = 1:length(vals)
    if vals[i] < 0
        vals[i] = -vals[i]
        for j = 1:size(vecs,1); V[j,i] = -V[j,i]; end
    end
end
return SVD(vecs, vals, V')

Even this much trouble may not be worth it, since it only reduces the allocations by a small factor — by far the biggest allocations are for constructing V anyway. There are some additional allocations that are done here for I and inside permute!, but it is probably not worth the trouble to eliminate them too.

[dkarrasch]

Thanks @stevengj and @oxinabox for your comments.

I knew one could do I = sortperm(vals; by=abs, rev=true), but together with vals[i] = -vals[i] that felt like "taking abs of vals twice", which I wanted to avoid. I could replace my vecs*Σ with

V = copy(vecs)
@inbounds for i = 1:length(σ)
    if σ[i] < 0
        for j = 1:size(vecs,1); V[j,i] = -V[j,i]; end
    end
end

to avoid "multiplication with 1s".

[oxinabox]

@stevengj thanks, I think I see now.
Yes, probably not worth it since introducing the branchs (with the if) is probably worse for performance than the extra allocations anyway.