Regresion: Eigvals of bigfloat symmetric matrices does not work anymore

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

2 Likes

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
2 Likes

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