Regresion: Eigvals of bigfloat symmetric matrices does not work anymore

lrnv · July 2, 2021, 4:05pm

Hi,

I have troubles to get eigenvalues of Symmetric matrices. Here is an example :

julia> X = big.(rand(10,10));

julia> X = X .+ X';

julia> eigvals(X)
ERROR: UndefVarError: eigvals not defined
Stacktrace:
 [1] top-level scope
   @ REPL[7]:1

julia> using LinearAlgebra

julia> eigvals(X)
ERROR: MethodError: no method matching eigvals!(::Matrix{BigFloat})
Closest candidates are:
  eigvals!(::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\tridiag.jl:293
  eigvals!(::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}, ::UnitRange) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\tridiag.jl:296
  eigvals!(::SymTridiagonal{var"#s832", V} where {var"#s832"<:Union{Float32, Float64}, V<:AbstractVector{var"#s832"}}, ::Real, ::Real) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\tridiag.jl:301
  ...
Stacktrace:
 [1] eigvals(A::Matrix{BigFloat}; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:326
 [2] eigvals(A::Matrix{BigFloat})
   @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:326
 [3] top-level scope
   @ REPL[9]:1

julia> using GenericLinearAlgebra

julia> eigvals(X)
10-element Vector{BigFloat}:
 -1.40896494093833723000919442666358694782060277534420585719212498189431446573352
 -0.9628570629553533920427015265454329908211049396906281502199649555971476817879393
 -0.765572319270339668165400589029211855446345366148112121821247891083801835348516
 -0.4988457885526648034975679430523053256471389136676962921706146821259835194905027
 -0.04967253288626804519859817605198605470194529038496302903243869042530808557833691
  0.7331486204381298592384410093739200409781239753824975991067791372244422807346604
  0.8346942215521769622694146059645529605227057642425001766162913257146661200708872
  1.089539314412129190793027457231843700294486594865326765338587929055980425160901
  2.201640880358843793326393001708644130370018951703733240381220100406889408981475
  9.844172310877760484923259844070267149553296139666547668993512708724577352991

julia> Y = Symmetric(X)
10×10 Symmetric{BigFloat, Matrix{BigFloat}}:
 1.81304   0.688118  1.22362   1.34483   1.22881   1.32082   1.00047   0.88352    1.64804   0.619048
 0.688118  0.243196  0.816204  0.573237  0.963247  0.381922  0.862619  0.123559   0.502359  0.979202
 1.22362   0.816204  0.639077  0.421363  1.01484   0.63954   0.83318   0.446309   1.56783   1.09914
 1.34483   0.573237  0.421363  1.59545   1.44708   0.526422  0.828333  0.49427    1.2904    0.719624
 1.22881   0.963247  1.01484   1.44708   1.96836   0.9886    0.345132  0.791248   1.5455    0.984035
 1.32082   0.381922  0.63954   0.526422  0.9886    1.34552   0.953025  1.38539    1.04277   1.62095
 1.00047   0.862619  0.83318   0.828333  0.345132  0.953025  0.858865  0.851573   0.425886  1.67392
 0.88352   0.123559  0.446309  0.49427   0.791248  1.38539   0.851573  0.0927274  1.02034   0.827739
 1.64804   0.502359  1.56783   1.2904    1.5455    1.04277   0.425886  1.02034    0.734844  1.08006
 0.619048  0.979202  1.09914   0.719624  0.984035  1.62095   1.67392   0.827739   1.08006   1.72621

julia> eigvals(Y)
ERROR: MethodError: no method matching eigvals!(::Symmetric{BigFloat, Matrix{BigFloat}})
Closest candidates are:
  eigvals!(::StridedMatrix{var"#s829"} where var"#s829"<:Union{Float32, Float64}; permute, scale, sortby) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:291
  eigvals!(::StridedMatrix{var"#s830"} where var"#s830"<:Union{ComplexF32, ComplexF64}; permute, scale, sortby) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:296
  eigvals!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\eigen.jl:544
  ...
Stacktrace:
 [1] eigvals(A::Symmetric{BigFloat, Matrix{BigFloat}})
   @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\symmetric.jl:740
 [2] top-level scope
   @ REPL[13]:1

julia> Z = Hermitian(X)
10×10 Hermitian{BigFloat, Matrix{BigFloat}}:
 1.81304   0.688118  1.22362   1.34483   1.22881   1.32082   1.00047   0.88352    1.64804   0.619048
 0.688118  0.243196  0.816204  0.573237  0.963247  0.381922  0.862619  0.123559   0.502359  0.979202
 1.22362   0.816204  0.639077  0.421363  1.01484   0.63954   0.83318   0.446309   1.56783   1.09914
 1.34483   0.573237  0.421363  1.59545   1.44708   0.526422  0.828333  0.49427    1.2904    0.719624
 1.22881   0.963247  1.01484   1.44708   1.96836   0.9886    0.345132  0.791248   1.5455    0.984035
 1.32082   0.381922  0.63954   0.526422  0.9886    1.34552   0.953025  1.38539    1.04277   1.62095
 1.00047   0.862619  0.83318   0.828333  0.345132  0.953025  0.858865  0.851573   0.425886  1.67392
 0.88352   0.123559  0.446309  0.49427   0.791248  1.38539   0.851573  0.0927274  1.02034   0.827739
 1.64804   0.502359  1.56783   1.2904    1.5455    1.04277   0.425886  1.02034    0.734844  1.08006
 0.619048  0.979202  1.09914   0.719624  0.984035  1.62095   1.67392   0.827739   1.08006   1.72621

julia> eigvals(Z)
10-element Vector{BigFloat}:
 -1.40896494093833723000919442666358694782060277534420585719212498189431446573352
 -0.9628570629553533920427015265454329908211049396906281502199649555971476817879393
 -0.765572319270339668165400589029211855446345366148112121821247891083801835348516
 -0.4988457885526648034975679430523053256471389136676962921706146821259835194905027
 -0.04967253288626804519859817605198605470194529038496302903243869042530808557833691
  0.7331486204381298592384410093739200409781239753824975991067791372244422807346604
  0.8346942215521769622694146059645529605227057642425001766162913257146661200708872
  1.089539314412129190793027457231843700294486594865326765338587929055980425160901
  2.201640880358843793326393001708644130370018951703733240381220100406889408981475
  9.844172310877760484923259844070267149553296139666547668993512708724577352991

julia>

The problem is that this call is far away in my call stack through several packages, and I cannot force the packages to upgrade their matrices to Hermitian instead of Symmetric.

Which square in the disptch field do you propose I patch on my local environnement to make this work ?

I tried to add :

GenericLinearAlgebra.eigvals!(X::Symmetric{BigFloat, Matrix{BigFloat}}) = GenericLinearAlgebra.eigvals!(GenericLinearAlgebra.Hermitian(X))

But it was’nt enough. The call stack I have on my real code is

ERROR: MethodError: no method matching geigvals!(::LinearAlgebra.Symmetric{BigFloat, Matrix{BigFloat}}, ::LinearAlgebra.QRIteration)
Closest candidates are:
  geigvals!(::LinearAlgebra.SymTridiagonal{T, V} where V<:AbstractVector{T}, ::Any) where T<:AbstractFloat at C:\Users\u009192\.julia\packages\GenericSchur\i4qpR\src\symtridiag.jl:24
Stacktrace:
 [1] eigvals!(A::LinearAlgebra.Symmetric{BigFloat, Matrix{BigFloat}}, alg::LinearAlgebra.QRIteration; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ GenericSchur ~\.julia\packages\GenericSchur\i4qpR\src\GenericSchur.jl:98
 [2] eigvals!(A::LinearAlgebra.Symmetric{BigFloat, Matrix{BigFloat}}, alg::LinearAlgebra.QRIteration) (repeats 2 times)
   @ GenericSchur ~\.julia\packages\GenericSchur\i4qpR\src\GenericSchur.jl:98
 [3] eigvals(A::LinearAlgebra.Symmetric{BigFloat, Matrix{BigFloat}})
   @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\symmetric.jl:740
 [4] quadform(x::Convex.AdditionAtom, A::Matrix{BigFloat})
   @ Convex ~\.julia\packages\Convex\ubaUN\src\atoms\second_order_cone\quadform.jl:12
 [5] top-level scope
   @ REPL[57]:1

I also tried adding :

LinearAlgebra.eigvals(A::LinearAlgebra.Symmetric{BigFloat, Matrix{BigFloat}}) = GenericLinearAlgebra.eigvals!(GenericLinearAlgebra.Hermitian(A))

And still same problem, although different call stack:

ERROR: MethodError: no method matching geigvals!(::LinearAlgebra.Hermitian{BigFloat, Matrix{BigFloat}}, 
::LinearAlgebra.QRIteration)
Closest candidates are:
  geigvals!(::LinearAlgebra.SymTridiagonal{T, V} where V<:AbstractVector{T}, ::Any) where T<:AbstractFloat at C:\Users\u009192\.julia\packages\GenericSchur\i4qpR\src\symtridiag.jl:24
Stacktrace:
 [1] eigvals!(A::LinearAlgebra.Hermitian{BigFloat, Matrix{BigFloat}}, alg::LinearAlgebra.QRIteration; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ GenericSchur ~\.julia\packages\GenericSchur\i4qpR\src\GenericSchur.jl:98
 [2] eigvals!(A::LinearAlgebra.Hermitian{BigFloat, Matrix{BigFloat}}, alg::LinearAlgebra.QRIteration) (repeats 2 times)
   @ GenericSchur ~\.julia\packages\GenericSchur\i4qpR\src\GenericSchur.jl:98
 [3] eigvals(A::LinearAlgebra.Symmetric{BigFloat, Matrix{BigFloat}})
   @ Main .\REPL[58]:1
 [4] quadform(x::Convex.AdditionAtom, A::Matrix{BigFloat})
   @ Convex ~\.julia\packages\Convex\ubaUN\src\atoms\second_order_cone\quadform.jl:12
 [5] top-level scope
   @ REPL[59]:1

Here is my MWE :

using Convex
x = Variable(3)
y = rand(3)
M = big.(rand(3,3))
M = M+M'

Prob = Convex.minimize(Convex.quadform(x - y, M), numeric_type = BigFloat) # error

using LinearAlgebra
M = Symmetric(M)

Prob = Convex.minimize(Convex.quadform(x - y, M), numeric_type = BigFloat) # error

using GenericLinearAlgebra

Prob = Convex.minimize(Convex.quadform(x - y, M), numeric_type = BigFloat) # error

M = Hermitian(M)

Prob = Convex.minimize(Convex.quadform(x - y, M), numeric_type = BigFloat) # error

M = GenericLinearAlgebra.Hermitian(M)
Prob = Convex.minimize(Convex.quadform(x - y, M), numeric_type = BigFloat) # error


# What should I do ?

EDIT: By taking a look, I found that :

julia> Symmetric{BigFloat, Matrix{BigFloat}} <: Hermitian{BigFloat, Matrix{BigFloat}}
false

Why is that ?

In this file : julia/symmetric.jl at master · JuliaLang/julia · GitHub, it says :

ishermitian(A::Hermitian) = true
ishermitian(A::Symmetric{<:Real}) = true
ishermitian(A::Symmetric{<:Complex}) = isreal(A)
issymmetric(A::Hermitian{<:Real}) = true
issymmetric(A::Hermitian{<:Complex}) = isreal(A)
issymmetric(A::Symmetric) = true

So clearly, for reals, Symmetric matrices could be treated by the same code as Hermitian matrices, which is probably why there is no code specific for them in GenericLinearAlgebra.jl.

Even with all this digging around, I still cannot find the right spot in the method field to patch

Ralph_Smith · July 3, 2021, 4:04am

Your “real code” and the examples are using different packages. The “real code” encounters a twisty gap in the method tree of my GenericSchur package, which I will address in the next release. (Thanks for your accidental contribution, and sorry for the inconvenience.)

A possible workaround for now is

import LinearAlgebra.eigvals
eigvals(A::Symmetric{BigFloat}) = eigen(Hermitian(A)).values

lrnv · July 5, 2021, 8:20am

Thanks a lot for the workaround. I’ll check your changes once they are released for sure.

Topic		Replies	Views
Methods to diagonalize a Matrix{BigFloat}? Numerics	22	6590	January 26, 2021
ERROR: UndefVarError: eigvecs not defined New to Julia	3	4742	September 14, 2018
Eigen not working for arbitrary precision BigFloat Numerics eigenvalues	16	2318	September 7, 2021
Eigenvalues are approximately correct, not exactly New to Julia linearalgebra , eigenvalues	6	887	April 4, 2022
Eigen solution of Matrix{Rational{BigInt}} General Usage	4	529	January 14, 2019

Regresion: Eigvals of bigfloat symmetric matrices does not work anymore

Related topics