Allow SVD of vectors

[dkarrasch]

Closes #39047.

[andreasvarga]

Occasionally a scaling of F.R[1] and F.Vt with -1 is also necessary to ensure nonnegative singular value. Daniel Karrasch <[email protected]> schrieb am Di., 5. Jan. 2021, 21:16:

…

***@***.**** commented on this pull request. ------------------------------ In stdlib/LinearAlgebra/src/svd.jl <#39087 (comment)>: > @@ -98,7 +98,9 @@ function svd!(A::StridedMatrix{T}; full::Bool = false, alg::Algorithm = default_ end SVD(u,s,vt) end - +function svd!(A::StridedVector{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where T<:BlasFloat + svd!(reshape(A, (length(A), 1)), full = full, alg = alg) I don't think the reshape was @stevengj <https://github.com/stevengj>'s issue, rather the fact that I call a function instead of writing out the result manually. Writing out the result is easy for the default full = false, but for full = true we need the left square unitary matrix. This is what I could come up with so far (which is super-ugly!): function mysvd!(A::StridedVector{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where T if !full normA = norm(A) normalize!(A) return SVD(A, [normA], [one(T)]) else F = qr!(reshape(A, (length(A), 1))) # need to reshape, no qr!(::Vector) return SVD(F.Q*Matrix(I, length(A), length(A)), [first(F.R)], fill!(Matrix{T}(undef, 1, 1), 1)) endend I'm not sure this is worth the hassle, but maybe we find better ways to code it. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#39087 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ALJDHEG6PV26CPHP5424VZ3SYNXTJANCNFSM4VS7HVBA> .

[dkarrasch]

We need to reshape here to make the function type stable: in the other branch, the left factor is a matrix, so it must be here as well. Next, since the left factor must be of matrix type, so must be the right factor. Still, I observed a speed-up of almost factor 10 compared to calling the line below.

[andreasvarga]

I played around a little bit with using a simple Householder reflector to do the whole job for a vector.
As basis I adapted the already existing code in LinearAlgebra.reflector!. Here is my code (still needs probably some polishing):

function svd1!(A::StridedVector;full = false)
    # based on LinearAlgebra.reflector! 
    n = length(A)
    @inbounds begin
        ξ1 = A[1]
        normA = norm(A)
        if iszero(normA)
            full ? (return SVD(Matrix{eltype(A)}(I,n,n), [normA], ones(eltype(A),1,1) ) )  : 
                    (return SVD(Matrix{eltype(A)}(I,n,1), [normA], ones(eltype(A),1,1)) )
        end
        ν = copysign(normA, real(ξ1))
        ξ1 += ν
        A[1] = -ν
        for i = 2:n
            A[i] /= ξ1
        end
    end
    τ = ξ1/ν 
    A[1] = one(eltype(A))
    full ? (return SVD(I - (conj(τ)*A)*A', [normA], -ones(eltype(A),1,1)) )  : 
             (return SVD(Matrix{eltype(A)}(I,n,1)-conj(τ)*A, [normA], -ones(eltype(A),1,1)) )
end

n = 4
A = rand(n)
Asave = copy(A)
# full case
A = copy(Asave);
F = svd1!(A,full = true); F.U
F.U*F.S[1,1]*Matrix{eltype(A)}(I,4,1)*F.Vt - Asave
# thin case
A = copy(Asave);
F = svd1!(A,full = false); U1
F.U*F.S*F.Vt - Asave

With this occasion I observed that svdvals does not work on a vector:

svdvals([1])
ERROR: MethodError: no method matching svdvals(::Array{Int64,1})
Closest candidates are:
  svdvals(::LinearAlgebra.AbstractTriangular; kwargs...) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\triangular.jl:2672
  svdvals(::StridedArray{T, 2}, ::StridedArray{T, 2}) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\svd.jl:510
  svdvals(::StridedArray{TA, 2}, ::StridedArray{TB, 2}) where {TA, TB} at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\svd.jl:536
  ...
Stacktrace:
 [1] top-level scope at REPL[193]:1

[dkarrasch]

With this occasion I observed that svdvals does not work on a vector:

That's right, currently it doesn't, but it's included in this PR already.

[dkarrasch]

I have tested your reflector-based code, and it takes much much longer than the reshape-based version (as in e8285b4). I have also included your zero-norm case handling, thanks for the reminder!

[dkarrasch]

Are we okay with this, or do we have ideas how to construct the full SVD more performant than what I currently have?