RFC: Methods for eigen on complex hermitian tridiagonal matrices

[MasonProtter]

Background:

We have a type SymTridiagonal which is an unnatural abomination type representing symmetric tridiagonal matrices, but we have no such type for hermitian tridiagonal matrices which are much more natural objects in linear algebra.

This leads to an unfortunate situation where trying to get the eigenvalues of a real SymTridiagonal matrix is very fast, but if you have complex elements and a hermitian structure, then you end up hitting the generic fallback and performance will greatly suffer.

The upside is that for any hermitian tridiagonal matrix A, there exists a diagonal unitary matrix S such that S' * A * S is a real symmetric tridiagonal matrix. B = S' * A * S will have the same main diagonal as A, and the off diagonals will simply be the absolute values of the off-diagonals of A.

This is something known by experts, and sounds plausible to people with some experience, but I think isn't particularly common knowledge amongst wider users of linear algebra, so if there's an easy way we can do this for them automatically, that'd be good.

This PR

I wasn't really sure what would be the best approach would be to put this transformation into LinearAlgebra, so I figured I'd just do the most obvious thing, put it out there and then request advice for how to integrate it better into the system.

Ideally, this would work with regular Tridiagonal{<:Complex} matrices by just checking if they're hermitian, but I was a bit wary of the resulting dispatch gymnastics, so I opted to only operate on Hermitian{Complex{T}, <:Tridiagonal} for now.

The approach is that for eigen, I generate the similarity transform S, and the then calculate the eigenvalues and eigenvectors of

B = SymTridiagonal(real.(A.d), abs.(A.dl))

and the plain eigenvalues, as well as the eigenvectors which are transformed by S.

For eigvals, that transformation is not necessary, so I simply do eigvals(SymTridiagonal(real.(d), abs.(dl)); kwargs...) which saves us a matmul relative to eigen

Any advice or comments on the implementation / approach before I start writing up some tests and such?

Performance comparisons

Before:

#+begin_src julia
let N = 500
    D = rand(N) |> complex
    E = rand(ComplexF64, N-1)

    A = Tridiagonal(E, D, adjoint.(E))
    @btime eigvals(Hermitian($A))
    @btime eigen(Hermitian($A))
end;
#+end_src

#+RESULTS:
:   16.222 ms (25 allocations: 4.21 MiB)
:   37.915 ms (29 allocations: 11.84 MiB)

and after

#+begin_src julia
let N = 500
    D = rand(N) |> complex
    E = rand(ComplexF64, N-1)

    A = Tridiagonal(E, D, adjoint.(E))
    @btime eigvals(Hermitian($A))
    @btime eigen(Hermitian($A))
end;
#+end_src

#+RESULTS:
:   2.758 ms (15 allocations: 23.89 KiB)
:   17.495 ms (39 allocations: 5.87 MiB)

---

Edit: I forgot to mention in the original post, this algorithm was sketched out in https://web.archive.org/web/20201028202318/http:https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=3592# , and it was brought to my attention by @stevengj with this Discourse post: https://discourse.julialang.org/t/what-is-the-point-of-the-symmetric-type/40773/5

[antoine-levitt]

Can't you dispatch to the lapack hermitian routine directly?

[MasonProtter]

Yeah, that might be preferable. I wonder if it's any faster than what I did. I preferred this approach because then I didn't have to worry about any of the different options which might be supplied to eigen and just let the real symmetric solver take care of them.

[dkarrasch]

Can't you dispatch to the lapack hermitian routine directly?

There doesn't seem to be such a direct hermitian routine. LAPACK suggests that users reduce a matrix of interest (probably dense hermitian?) to tridiagonal form first and then use one of the standard routines. For the time being, I think merging this PR as is is a good idea. If there is sufficient interest and room for further improvement, then we could do this later?

[antoine-levitt]

There doesn't seem to be such a direct hermitian routine

Oh wow that's crazy but yeah that appears to be correct. So OK let's just merge this. Someone really needs to rewrite lapack at some point...

[MasonProtter]

@dkarrasch @antoine-levitt any thoughts on if I should try and add an intermediate method for regular eigen(::Tridiagonal) which checks if it's Hermitian and if so, dispatches to this method? Or is it not worth it?

My concern is that trying to intercept that case, and then if it isn't hermitian, carry on as before is kinda hard. I guess I'd need to do an invoke, but that could cause problems.

[brenhinkeller]

It seems like this is already a significant performance improvement!

[danielwe]

Would this warrant a new HermTridiagonal type? It would provide a nice symmetry between Symmetric/Hermitian and SymTridiagonal/HermTridiagonal, which might be helpful for discoverability and useability of specialized methods. (Never mind it might be the first LinearAlgebra structured matrix type that does not represent a LAPACK-supported structure.)

Related reading: #43983.

There doesn't seem to be such a direct hermitian routine

As far as I understand, LAPACK's support for real symmetric tridiagonal matrices exists not so much for the interest in such matrices in their own right, but because they are an intermediate representation in the eigendecomposition of hermitian matrices, much like Hessenberg matrices for QR decomposition. Thus SymTridiagonal exists and is Sym rather than Herm because LinearAlgebra is in large part designed as a LAPACK wrapper. Somewhat tangential, but @andreasnoack's comment here gives an idea about why these types exist and work the way they do: #30137 (comment)

The way I think about Bidiagonal, SymTridiagonal, and Tridiagional is that they are the building blocks for SVD, the symmetric eigenproblem, and a special case of LU respectively and that it shouldn't really be expected that all linear algebra methods work for these.

[KlausC]

Hi @MasonProtter , I agree with this PR; was missing the HermTridiagonal as well and can live with your proposed API.

I will suggest some improvements in implementation, which are attached to the source code.

[stevengj]

LGTM.

Needs tests, though.

[MasonProtter]

Is there anything in particular you'd like to see tested? This path is already in the test suite fwiw, and compared to the default dense result, both for eigenvectors and eigenvalues

[danielwe]

To elaborate on the argument for defining HermTridiagonal: it seems inconsistent that optimized methods are defined for Hermitian{<:Complex,<:Tridiagonal} and SymTridiagonal{<:Real}, but not Symmetric{<:Real,<:Tridiagonal}.

[stevengj]

Where is this tested? I don't see the Hermitian constructor used in https://github.com/JuliaLang/julia/blob/c1e1d5c451f1f565f1c1fe529f333931633d1a71/stdlib/LinearAlgebra/test/tridiag.jl

[MasonProtter]

Ah yeah, I see how that's confusing. For some reason, this test lives in test/eigen.jl instead of test/tridag.jl or test/symmetriceigen.jl:

julia/stdlib/LinearAlgebra/test/eigen.jl

Line 48 in c1e1d5c

for T in (Tridiagonal(a), Hermitian(Tridiagonal(a)))

The only reason I knew it was being tested was that I was getting test failures before because I wasn't checking if the Hermitian was U or L

[MasonProtter]

To elaborate on the argument for defining HermTridiagonal: it seems inconsistent that optimized methods are defined for Hermitian{<:Complex,<:Tridiagonal} and SymTridiagonal{<:Real}, but not Symmetric{<:Real,<:Tridiagonal}.

While it'd be nice, adding a whole new HermTridiagonal type seems like a much bigger project that I'm not sure I am currently willing to take on. I was hoping to just make this quite important, but not widespread optimized routine available to people who might need it.

I don't think it's really that unusual though to want to specialized on Hermitian{<:Complex,<:Tridiagonal}.

[danielwe]

I don't think it's really that unusual though to want to specialized on Hermitian{<:Complex,<:Tridiagonal}.

It's certainly not unusual to want to diagonalize a hermitian tridiagonal matrix, I've often wanted that myself. It's also reasonable that Hermitian{<:Complex,<:Tridiagonal} hits the optimized method. I'm just thinking about consistency in the zoo of structured matrix types and their behavior, the lack of which has always been a bit of a stumbling block with LinearAlgebra.jl.

Perhaps this PR should also add methods for Hermitian{<:Real,<:Tridiagonal} and Symmetric{<:Real,<:Tridiagonal}? Those would be simple forwarding methods like what you're doing for eigvals, except you don't even have to use abs. I think this would help a lot with consistency and reliability in dispatching to the best method available.

[danielwe]

I've often wanted that myself

Btw., emphasizing this and that I'm really excited about having this transformation in LinearAlgebra. In the past, I've always resorted to making a dense copy. Thank you for making this happen!

[dkarrasch]

This appears to be NEWS-worthy, doesn't it? And then we can merge, I think.

As for the real cases, we could have another PR, and then discuss whether we should have

function Symmetric(A::Tridiagonal, uplo::Char)
    checkuplo(uplo)
    uplo == 'U' ? SymTridiagonal(A.d, A.dv) : SymTridiagonal(A.d, A.ev)
end

or whether we want specialized eig* methods for Hermitian{<:Real,<:Tridiagonal} and Symmetric{<:Real,<:Tridiagonal}. To me, at least, it's not immediately clear what to do.

[MasonProtter]

Okay, I added the NEWS.md entry, should be good to go now!