Stream: helpdesk (published)

Topic: ✔ Hermitian tridiagonal arrays

Mason Protter (Apr 26 2023 at 20:04):

Do we have any eigensolvers laying around that are specialized for hermitian tridiagonal arrays?

Mason Protter (Apr 26 2023 at 20:05):

We have eigen(::SymTridiagonal) but SymTridiagonal only seems to support being symmetric but not Hermitian which is a really weird choice because symmetric, non-hermitian matrices are quite rare in linear algebra

Mason Protter (Apr 26 2023 at 22:37):

Okay, so using https://web.archive.org/web/20201028202318/http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=3592#p9189 it seems the solution is to use a similarity transform.

real_form(A::Hermitian{T, <:Tridiagonal}) where {T <: Real} = (;T = A,  S = Diagonal(ones(T, size(A, 1))), )
function real_form(A::Hermitian{T, <:Tridiagonal}) where {T}
    (; dl, d) = parent(A)
    N = length(d)
    if N == 1
        return (; T = SymTridiagonal(real.(d), real.(dl)), S = Diagonal([one(T)]))
    end
    S = Vector{ComplexF64}(undef, N)
    E = dl
    Er = abs.(E)
    S[1] = 1
    S[2] = Er[1] == 0 ? one(T) : S[1] * E[1]/Er[1]
    for i ∈ 2:N-1
        S[i+1] = Er[i] == 0 ? one(T) : S[i] * E[i]/Er[i]
    end
    (; T = SymTridiagonal(real.(d), Er), S = Diagonal(S))
end

This will produce a real SymTridiagonal matrix T which has the same eigenvalues as the hermitian matrix A, as well as a diagonal matrix S such that S * T * S' ≈ A, so that if U are the eigenvectors of T, then

S * U * Diagonal(eigvals(T)) * U' * S' ≈ A

Notification Bot (Apr 26 2023 at 22:38):

Mason Protter has marked this topic as resolved.

Last updated: Nov 06 2024 at 04:40 UTC