Inverse of a SVD factorization of a complex matrix broken in Julia 1.3

[carstenbauer]

On Julia 1.1:

julia> using LinearAlgebra

julia> inv(svd(rand(ComplexF64,2,2)))
2×2 Array{Complex{Float64},2}:
 -1.03305+2.82321im  -1.07298-3.37716im
 0.623141-1.89883im  0.946827+1.23427im

On Julia 1.3:

julia> using LinearAlgebra

julia> inv(svd(rand(ComplexF64,2,2)))
ERROR: MethodError: no method matching eps(::Type{Complex{Float64}})
Closest candidates are:
  eps(::Dates.Time) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\Dates\src\types.jl:387
  eps(::Dates.Date) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\Dates\src\types.jl:386
  eps(::Dates.DateTime) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\Dates\src\types.jl:385
  ...
Stacktrace:
 [1] inv(::SVD{Complex{Float64},Float64,Array{Complex{Float64},2}}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\svd.jl:284
 [2] top-level scope at REPL[2]:1

The relevant eps(T) call is here.

I guess I'm to blame for this, because it seems this was introduced in JuliaLang/julia#32126. Sorry 🤦‍♂

[carstenbauer]

Goes without saying that I'll prepare a PR to fix this :)

[dkarrasch]

Can't you quickly fix it by replacing eps(T)*F.S[1] by eps(abs(F.S[1]))*one(T)?

[carstenbauer]

I was thinking eps(real(T))*F.S[1].

[StefanKarpinski]

Isn't it better to take the eps of F.S[1] instead of multiplying?

[carstenbauer]

Well, strictly speaking, this would be a slightly different truncation. I chose eps(real(T))*F.S[1] to match what we do in ldiv!.

Personally, I don't mind if we'd drop the truncation altogether. But in any case, I think truncation changes should be discussed in #652.

[carstenbauer]

On an unrelated note, shouldn't this issue get a bug label (or similar)?

[dkarrasch]

But isn't eps(real(T))*F.S[1] == eps(F.S[1]) just without multiplication? In my suggestion, I forgot that singular values are real, so away goes the abs and the *one(T).

[carstenbauer]

Maybe I'm misunderstanding - which might well be the case when talking about floating point numbers - but why would the two be equivalent? Wouldn't it depend on the value of F.S[1]?

julia> s = rand()
0.9013815392525062

julia> eps(real(typeof(s)))*s == eps(s)
false

julia> eps(real(typeof(s)))*s - eps(s)
8.912460530752368e-17

This value is smaller than eps(s) but in a comparison the difference would still matter, wouldn't it?

[dkarrasch]

But if you divide the difference by s, you get 9.887556093220253e-17, so the relative error of the two epss is below eps(Float64).

[carstenbauer]

This I get. But isn't it the absolute value that matters here? Eventually, in searchsortedlast, we are comparing the elements of F.S to our eps in absolute terms, aren't we?

Sorry for perhaps asking the obvious :)