Full/square SVD and dimension mismatch

I noticed that

A = rand(2, 1)
U, S, V = svd(A; thin=false)
U*diagm(S)*V' ≈ A

produces the error (Julia 0.6.2):

ERROR: DimensionMismatch("A has dimensions (2,2) but B has dimensions (1,1)")
Stacktrace:
 [1] gemm_wrapper!(::Array{Float64,2}, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at ./linalg/matmul.jl:345
 [2] *(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at ./operators.jl:424

I expected S to be of size (2,) and V to be of size (1,2). This holds for any MxN matrix A when M /= N. The behaviour is the same with eigfact of course.

Is this fine or should I report it as a bug?

On the classic full SVD of A, m\times n, you’ll have orthogonal U and V of sizes m\times m and n \times n respectively and \Sigma m \times n. You can’t have V with more columns than that.
Note that diagm(S) won’t give you a rectangular matrix, but you can use spdiagm.

U * spdiagm(S, 0, 2, 1) * V'
1 Like

Of course, now it makes sense. Thanks for the clarification!