Fix for eigen not working properly for 2x2 Hermitian SMatrices

[guberger]

Addressing issue #956 on the numerical instability of the eigen function for Hermitian 2x2 SMatrix. The previous implementation also failed when A is not Hermitian but Hermitian(A) is passed as argument to eigen. For instance:

julia> A = Hermitian(@SMatrix [1.0+1.0im 1.0; 1.0 2.0-1.0im])      
2×2 Hermitian{ComplexF64, SMatrix{2, 2, ComplexF64, 4}}:
 1.0+0.0im  1.0+0.0im
 1.0-0.0im  2.0+0.0im
julia> B = Matrix(A)
2×2 Matrix{ComplexF64}:
 1.0+0.0im  1.0+0.0im
 1.0-0.0im  2.0+0.0im
julia> eigen(A)
Eigen{ComplexF64, Float64, SMatrix{2, 2, ComplexF64, 4}, SVector{2, Float64}}
values:
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 2.0
vectors:
2×2 SMatrix{2, 2, ComplexF64, 4} with indices SOneTo(2)×SOneTo(2):
 -0.57735+0.57735im       0.0+0.707107im
  0.57735-0.0im      0.707107-0.0im
julia> eigen(B)
Eigen{ComplexF64, Float64, Matrix{ComplexF64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 0.38196601125010515
 2.618033988749895
vectors:
2×2 Matrix{ComplexF64}:
 -0.850651+0.0im  0.525731+0.0im
  0.525731+0.0im  0.850651+0.0im

The new implementation seems also faster as suggested by the following testing code:

using LinearAlgebra
using StaticArrays

@inline function _eig2(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T},
        permute, scale) where {T <: Real}
    a = A.data
    TA = eltype(A)

    @inbounds if A.uplo == 'U'
        if !iszero(a[3]) # A is not diagonal
            t_half = (real(a[1]) + real(a[4])) / 2
            s_half = (real(a[1]) - real(a[4])) / 2
            tmp2 = s_half^2 + abs2(a[3]) # expansion of (tr(A)^2 - 4*det(A)) / 4
            # use abs2 to ensure that it is real and > 0 (not sure about x'x being always > 0 if a[3] ≠ 0)
            tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
            vals = SVector(t_half - tmp, t_half + tmp)

            v11 = s_half - tmp
            n1 = sqrt(abs2(v11) + abs2(a[3]))
            v11 = v11 / n1
            v12 = a[3]' / n1

            v21 = s_half + tmp
            n2 = sqrt(abs2(v21) + abs2(a[3]))
            v21 = v21 / n2
            v22 = a[3]' / n2

            vecs = @SMatrix [ v11  v21 ;
                              v12  v22 ]

            return Eigen(vals, vecs)
        end
    else # A.uplo == 'L'
        if !iszero(a[2]) # A is not diagonal
            t_half = (real(a[1]) + real(a[4])) / 2
            s_half = (real(a[1]) - real(a[4])) / 2 
            tmp2 = s_half^2 + abs2(a[2]) # expansion of (tr(A)^2 - 4*det(A)) / 4
            # use abs2 to ensure that it is real and > 0 (not sure about x'x being always > 0 if a[3] ≠ 0)
            tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
            vals = SVector(t_half - tmp, t_half + tmp)

            v11 = s_half - tmp
            n1 = sqrt(abs2(v11) + abs2(a[3]))
            v11 = v11 / n1
            v12 = a[2] / n1

            v21 = s_half + tmp
            n2 = sqrt(abs2(v21) + abs2(a[3]))
            v21 = v21 / n2
            v22 = a[2] / n2

            vecs = @SMatrix [ v11  v21 ;
                              v12  v22 ]

            return Eigen(vals, vecs)
        end
    end

    # A must be diagonal if we reached this point; treatment of uplo 'L' and 'U' is then identical
    A11 = real(a[1])
    A22 = real(a[4])
    if A11 < A22
        vals = SVector(A11, A22)
        vecs = @SMatrix [convert(TA, 1) convert(TA, 0);
                         convert(TA, 0) convert(TA, 1)]
    else # A22 <= A11
        vals = SVector(A22, A11)
        vecs = @SMatrix [convert(TA, 0) convert(TA, 1);
                         convert(TA, 1) convert(TA, 0)]
    end
    return Eigen(vals,vecs)
end

function test_eigs(nsample)
    As_list = [SMatrix{2,2}(rand(Complex{Float64}, 2, 2)) for i = 1:nsample]
    A_list = map(As -> Hermitian(As + As'), As_list)

    local eig

    for A in A_list
        eig1 = StaticArrays._eig(Size(A), A, true, true)
        eig2 = _eig2(Size(A), A, true, true)
        @assert norm(eig1.values - eig2.values) < 1e-8
        @assert norm(eig1.vectors - eig2.vectors) < 1e-8
    end

    @time for A in A_list
        eig = StaticArrays._eig(Size(A), A, true, true)
    end

    @time for A in A_list
        eig = _eig2(Size(A), A, true, true)
    end

    return eig
end

test_eigs(100_000);

which outputs:

  0.010825 seconds
  0.001010 second

Version info:

Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS =

[guberger]

I also updated the function _eigvals for 2x2 Hermitian SMatrices with the same approach. Again the performance tests suggest that the new implementation is favorable:

using LinearAlgebra
using StaticArrays

@inline function _eig2(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T},
        permute, scale) where {T <: Real}
    a = A.data
    TA = eltype(A)

    @inbounds if A.uplo == 'U'
        if !iszero(a[3]) # A is not diagonal
            t_half = (real(a[1]) + real(a[4])) / 2
            s_half = (real(a[1]) - real(a[4])) / 2
            tmp2 = s_half^2 + abs2(a[3]) # expansion of (tr(A)^2 - 4*det(A)) / 4
            # use abs2 to ensure that it is real and > 0 (not sure about x'x being always > 0 if a[3] ≠ 0)
            tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
            vals = SVector(t_half - tmp, t_half + tmp)

            v11 = s_half - tmp
            n1 = sqrt(abs2(v11) + abs2(a[3]))
            v11 = v11 / n1
            v12 = a[3]' / n1

            v21 = s_half + tmp
            n2 = sqrt(abs2(v21) + abs2(a[3]))
            v21 = v21 / n2
            v22 = a[3]' / n2

            vecs = @SMatrix [ v11  v21 ;
                              v12  v22 ]

            return Eigen(vals, vecs)
        end
    else # A.uplo == 'L'
        if !iszero(a[2]) # A is not diagonal
            t_half = (real(a[1]) + real(a[4])) / 2
            s_half = (real(a[1]) - real(a[4])) / 2 
            tmp2 = s_half^2 + abs2(a[2]) # expansion of (tr(A)^2 - 4*det(A)) / 4
            # use abs2 to ensure that it is real and > 0 (not sure about x'x being always > 0 if a[3] ≠ 0)
            tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
            vals = SVector(t_half - tmp, t_half + tmp)

            v11 = s_half - tmp
            n1 = sqrt(abs2(v11) + abs2(a[3]))
            v11 = v11 / n1
            v12 = a[2] / n1

            v21 = s_half + tmp
            n2 = sqrt(abs2(v21) + abs2(a[3]))
            v21 = v21 / n2
            v22 = a[2] / n2

            vecs = @SMatrix [ v11  v21 ;
                              v12  v22 ]

            return Eigen(vals, vecs)
        end
    end

    # A must be diagonal if we reached this point; treatment of uplo 'L' and 'U' is then identical
    A11 = real(a[1])
    A22 = real(a[4])
    if A11 < A22
        vals = SVector(A11, A22)
        vecs = @SMatrix [convert(TA, 1) convert(TA, 0);
                         convert(TA, 0) convert(TA, 1)]
    else # A22 <= A11
        vals = SVector(A22, A11)
        vecs = @SMatrix [convert(TA, 0) convert(TA, 1);
                         convert(TA, 1) convert(TA, 0)]
    end
    return Eigen(vals,vecs)
end

@inline function _eigvals2(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T},
        permute, scale) where {T <: Real}
    a = A.data
    @inbounds t_half = (real(a[1]) + real(a[4])) / 2
    @inbounds s_half = (real(a[1]) - real(a[4])) / 2
    @inbounds tmp2 = A.uplo == 'U' ? s_half^2 + abs2(a[3]) : s_half^2 + abs2(a[2]) # expansion of (tr(A)^2 - 4*det(A)) / 4
    tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
    return SVector(t_half - tmp, t_half + tmp)
end

function test_eigs(nsample)
    As_list = [SMatrix{2,2}(rand(Complex{Float64}, 2, 2)) for i = 1:nsample]
    A_list = map(As -> Hermitian(As + As'), As_list)

    for A in A_list
        eig1 = StaticArrays._eig(Size(A), A, true, true)
        eig2 = _eig2(Size(A), A, true, true)
        @assert norm(eig1.values - eig2.values) < 1e-8
        @assert norm(eig1.vectors - eig2.vectors) < 1e-8
        eigvals1 = StaticArrays._eigvals(Size(A), A, true, true)
        eigvals2 = _eigvals2(Size(A), A, true, true)
        @assert norm(eigvals1 - eigvals2) < 1e-8
    end

    local eig

    @time for A in A_list
        eig = _eig2(Size(A), A, true, true)
    end

    @time for A in A_list
        eig = StaticArrays._eig(Size(A), A, true, true)
    end

    local eigvals

    @time for A in A_list
        eigvals = _eig2(Size(A), A, true, true)
    end

    @time for A in A_list
        eigvals = StaticArrays._eig(Size(A), A, true, true)
    end

    return eig, eigvals
end

test_eigs(100_000);

outputs

  0.000779 seconds
  0.011829 seconds
  0.000925 seconds
  0.012771 seconds

[mschauer]

So there are these hard cases and we have to be careful not just moving an insufficient large table cloth on a table from one corner to the other and back. So let's add a test and we could look at #694 , Ferrite-FEM/Tensors.jl#173, #524 to check if there are other relevant test cases.