Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
carstenbauer opened this issue Feb 25, 2020 · 10 comments · Fixed by #34872
Closed

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

carstenbauer opened this issue Feb 25, 2020 · 10 comments · Fixed by #34872
Labels
domain:linear algebra Linear algebra kind:bug Indicates an unexpected problem or unintended behavior

Comments

@carstenbauer
Copy link
Member

carstenbauer commented Feb 25, 2020

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 #32126. Sorry 🤦‍♂

@carstenbauer
Copy link
Member Author

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

@dkarrasch
Copy link
Member

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

@carstenbauer
Copy link
Member Author

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

carstenbauer added a commit to carstenbauer/julia that referenced this issue Feb 25, 2020
@StefanKarpinski
Copy link
Sponsor Member

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

@carstenbauer
Copy link
Member Author

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 #32880.

@carstenbauer
Copy link
Member Author

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

@rfourquet rfourquet added kind:bug Indicates an unexpected problem or unintended behavior domain:linear algebra Linear algebra labels Feb 25, 2020
@dkarrasch
Copy link
Member

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
Copy link
Member Author

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
Copy link
Member

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
Copy link
Member Author

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 :)

carstenbauer added a commit to carstenbauer/julia that referenced this issue Feb 29, 2020
dkarrasch pushed a commit that referenced this issue Mar 2, 2020
KristofferC pushed a commit that referenced this issue Mar 23, 2020
ravibitsgoa pushed a commit to ravibitsgoa/julia that referenced this issue Apr 9, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
domain:linear algebra Linear algebra kind:bug Indicates an unexpected problem or unintended behavior
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants