search

JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.47k stars 5.46k forks source link

Eigen fails for real-symmetric matrix #40260

Closed timholy closed 1 month ago

timholy commented 3 years ago 
using LinearAlgebra

Q11, Q22, Q33 = -0.0003713258198415712, -3.259560133917638e-5, -0.00031149545839469857
Q1e, Q2e, Q3e, Qee = 0.00011266199942972238, 0.00022565487110446598, 0.0001394626512764405, 0.4901233123942332
Q = [Q11 zeros(1, 8) Q1e; 0 Q22 zeros(1, 7) Q2e; 0 0 Q33 zeros(1, 6) Q3e; 0 0 0 Q33 zeros(1, 5) Q3e; 0 0 0 0 Q33 zeros(1, 4) Q3e; 0 0 0 0 0 Q33 0 0 0 Q3e; 0 0 0 0 0 0 Q33 0 0 Q3e; 0 0 0 0 0 0 0 Q33 0 Q3e; zeros(1, 8) Q33 Q3e; Q1e Q2e fill(Q3e, 1, 7) Qee]
@show issymmetric(Q)
eigen(Q)
julia> include("eigenbug.jl")
issymmetric(Q) = true
ERROR: LoadError: LAPACKException(1)
Stacktrace:
 [1] chklapackerror(ret::Int64)
   @ LinearAlgebra.LAPACK ~/src/julia-1/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:38
 [2] syevr!(jobz::Char, range::Char, uplo::Char, A::Matrix{Float64}, vl::Float64, vu::Float64, il::Int64, iu::Int64, abstol::Float64)
   @ LinearAlgebra.LAPACK ~/src/julia-1/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/lapack.jl:5109
 [3] #eigen!#95
   @ ~/src/julia-1/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:675 [inlined]
 [4] eigen!(A::Matrix{Float64}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby))
   @ LinearAlgebra ~/src/julia-1/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:151
 [5] eigen(A::Matrix{Float64}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby))
   @ LinearAlgebra ~/src/julia-1/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:237
 [6] eigen(A::Matrix{Float64})
   @ LinearAlgebra ~/src/julia-1/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/eigen.jl:235
 [7] top-level scope
   @ /tmp/eigenbug.jl:7
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:444
 [9] top-level scope
   @ REPL[1]:1
in expression starting at /tmp/eigenbug.jl:7

It fails during the second iteration, but I've not looked more deeply than that.

jaakkor2 commented 3 years ago

Issue on Linux, seems to work on Mac and Windows

julia> include("eigenbug.jl")
issymmetric(Q) = true
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
10-element Vector{Float64}:
 -0.0003713518258088031
 -0.0003117730492286244
 -0.00031149545839469857
 -0.00031149545839469857
 -0.00031149545839467965
 -0.00031149545839467965
 -0.00031149545839467965
 -0.00031149545839467965
 -3.269937555021306e-5
  0.49012371976524544
vectors:
10×10 Matrix{Float64}:
  0.999999      0.00142322    3.81639e-17   6.17562e-16  …  -6.93889e-17  -2.30371e-15  -0.000153003  0.00022969
  0.000153763  -0.000608082   1.17961e-16   8.32667e-17     -2.77556e-17   2.77556e-16   1.0          0.000460373
  0.000537827  -0.377964      0.883255     -0.0957135       -0.128067     -0.119113     -0.000230047  0.000284365
  0.000537827  -0.377964     -0.2046        0.578961        -0.19561       0.0835241    -0.000230047  0.000284365
  0.000537827  -0.377964     -0.231285      0.11919          0.0802217    -0.0313724    -0.000230047  0.000284365
  0.000537827  -0.377964      0.020105      0.180654     …   0.813838     -0.200094     -0.000230047  0.000284365
  0.000537827  -0.377964     -0.00647631   -0.164934        -0.0345753     0.86541      -0.000230047  0.000284365
  0.000537827  -0.377964     -0.32484      -0.750971        -0.0111665    -0.228941     -0.000230047  0.000284365
  0.000537827  -0.377964     -0.136159      0.132814        -0.524641     -0.369414     -0.000230047  0.000284365
 -0.000230832   0.000752311   0.0           0.0              3.36131e-17  -1.24272e-16  -0.00045988   1.0
dkarrasch commented 3 years ago

Fails for me on Mac with the same error as in the OP on Julia v1.0 and v1.6.

jaakkor2 commented 3 years ago

My old Mac is running 10.11.6 OS X El Capitan, and eigen worked.

Julia Version 1.6.0
Commit f9720dc2eb (2021-03-24 12:55 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin19.6.0)
  CPU: Intel(R) Core(TM)2 Duo CPU     P7350  @ 2.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, penryn)
timholy commented 3 years ago

It may more likely be CPU-dependent. For me,

julia> versioninfo()
Julia Version 1.6.1-pre.0
Commit f25d9885f6* (2021-03-25 03:13 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_CPU_THREADS = 4
pablosanjose commented 3 years ago

This reminds me of https://github.com/Reference-LAPACK/lapack/issues/475. That was also CPU dependent, and rather unpredictable.

The issue there was for a badly broken iteration procedure in the generalized eigen(A, B) methods (or schur(A, B), ultimately). Iit got fixed upstream, but only for the generalized case. Maybe the standard case also suffers from something similar?

pablosanjose commented 3 years ago

By the way, much like in Reference-LAPACK/lapack#475, this fails with OpenBLAS for me, but does not fail when using Julia 1.5.3 + MKL.jl

johnnychen94 commented 3 years ago

It might be related to the high decimal accuracy and float point errors it introduces, a workaround is to use fewer decimal points... 

julia> @show issymmetric(Float32.(Q))
issymmetric(Float32.(Q)) = true
true

julia> eigen(Float32.(Q))
Eigen{Float32, Float32, Matrix{Float32}, Vector{Float32}}
values:
julia> eigen(round.(Q; digits=64))
ERROR: LAPACKException(1)
Stacktrace:
 [1] chklapackerror(ret::Int64)
   @ LinearAlgebra.LAPACK /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/lapack.jl:43
 [2] syevr!(jobz::Char, range::Char, uplo::Char, A::Matrix{Float64}, vl::Float64, vu::Float64, il::Int64, iu::Int64, abstol::Float64)
   @ LinearAlgebra.LAPACK /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/lapack.jl:5120
 [3] #eigen!#169
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetriceigen.jl:4 [inlined]
 [4] eigen!(A::Matrix{Float64}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby))
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:151
 [5] eigen(A::Matrix{Float64}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby))
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:237
 [6] eigen(A::Matrix{Float64})
   @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/eigen.jl:235
 [7] top-level scope
   @ REPL[16]:1

julia> eigen(round.(Q; digits=63)) # works :)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
andreasnoack commented 3 years ago

I can reproduce this. The symmetric tridiagonal matrix looks fine

julia> diag(LAPACK.hetrd!('L', copy(Q))[1], 0)
10-element Vector{Float64}:
 -0.0003713258198415712
  0.4901233123942331
 -0.00023557898323602305
 -0.00010851207649784496
 -0.00031149545839469846
 -0.00031149545839469857
 -0.00031149545839469857
 -0.00031149545839469857
 -0.00031149545839469857
 -0.0003114954583946985

julia> diag(LAPACK.hetrd!('L', copy(Q))[1], -1)
9-element Vector{Float64}:
 -0.00011266199942972238
  0.00043251466860733215
 -0.00012413614650615042
 -2.065162908528577e-19
  1.489831575676201e-19
 -1.470481574344822e-20
 -2.4468077391424153e-20
  4.793334368155311e-20
  4.8633506824568657e-20

but maybe it's problem for the dqds algorithm that many of the eigenvalues seem to already have converged. I suspect it's either a bug or a limitation in the algorithm. It might be good to expose the QR algorithm as well and maybe default to using it since it should generally be more robust and indeed it works fine here

julia> eigvals(Q)
10-element Vector{Float64}:
 -0.0003713518258088077
 -0.0003117730492286513
 -0.00031149545839469873
 -0.0003114954583946987
 -0.0003114954583946986
 -0.0003114954583946985
 -0.0003114954583946985
 -0.0003114954583946984
 -3.26993755502334e-5
  0.49012371976524544

(eigvals is based on the square root free QR)

dkarrasch commented 8 months ago

This still fails on my Mac machine on v1.6, but works fine on v1.9 and v1.10. Does the error persist on other machines or can this be closed?

giordano commented 1 month ago

I could reproduce the issue on AMD Ryzen 9 3950X 16-Core Processor with Julia v1.6, v1.8, v1.9, but not v1.7, v1.10, v1.11 and nightly (5230d27de95), so I'm assuming this was fixed on recent versions of OpenBLAS and I'm going to close this ticket.