inconsistent Hermitian-ness of SymTridiagonal{Complex}

[stevengj]

Julia is inconsistent about whether SymTridiagonal is Hermitian or merely symmetric for complex types.

The comment at the top of tridiag.jl says that it represents "Hermitian tridiagonal matrices". The manual says that it is a "real symmetric tridiagonal matrix", which is clearly wrong since complex elements are allowed.

If I define A = SymTridiagonal([1,2],[3im]), it prints as:

2x2 Base.LinAlg.SymTridiagonal{Complex{Int64}}:
 1+0im  0+3im
 0+3im  2+0im

but full(A) gives

2x2 Array{Complex{Int64},2}:
 1+0im  0-3im
 0+3im  2+0im

and worse, A * x and full(A) * x disagree.

The SymTridiagonal(A::AbstractMatrix) constructor requires that A be symmetric, not Hermitian, and transpose is defined to be the identity while ctranspose conjugates it (i.e. it treats the matrix as symmetric, not Hermitian).

My feeling is that making it Hermitian would be much more useful. Factorization right now is of the L*D*L.' form, but of course we could use L*D*L' for Hermitian matrices. Cholesky could be supported. And efficient eigensolvers exist because complex-Hermitian matrices are trivially (via a diagonal unitary scaling) similar to real-symmetric matrices.

I guess there is a place for complex-symmetric tridiagonal matrices too, but I don't see that we gain much in that case vs. just using Tridiagonal.

cc: @andreasnoack, @jiahao

[jiahao]

This problem might be a good motivation to get started on JuliaLang/julia#8240. We could:

• typealias Tridiagonal{T} Banded{T,1,1, :col}
• deprecate SymTridiagonal for Hermitian{Tridiagonal} and Symmetric{Tridiagonal}.

[andreasnoack]

@stevengj I wonder if it is only a bug in full and that the rest of the code is handling symmetric matrices correctly (and a wrong comment). I agree that Hermitian feels more natural than complex symmetric, but what about restricting the type to real numbers? I guess that this type is mainly for Hermitian eigenvalue problems and LAPACK reduces general Hermitian problems directly to real symmetric tridiagonal problems. Are you aware of cases where it is more natural to formulate a Hermitian tridiagonal eigenproblem with complex off diagonals?

@jiahao There are some considerations about extra memory allocation here.

1. Our present Tridiagonal allocates an extra super diagonal which is used when factorizing the matrix. This is necessary because of pivoting.
2. Changing SymTridiagonal to Hermitian{Tridiagonal} introduces an uncessary super diagonal. This is kind of similar to the unused triangle of Hermitian{Matrix}.

[stevengj]

@andreasnoack, a simple example is for a 1d 2nd-order PDE with Bloch-periodic boundary conditions (where one side = other side × complex phase) whoops, no, that is not tridiagonal because the boundary conditions introduce two (complex) entries in the upper-right and lower-left corners. This can be transformed back to real-symmetric with a scaling, but the complex formulation is much more convenient.

I agree that the short-term fix is just to change full and the comment, and possibly to restrict it to real types for now.

[stevengj]

I'm not super-concerned about the memory required by the extra diagonal, since it is only a constant factor....it's hard for me to imagine a case in which a tridiagonal solver is memory-bound.

[andreasnoack]

I have now fixed the convert method called by full and changed the comment on both 0.3 and 0.4. The eigvals! and eigfact! methods have also been restricted to real element types for now. I believe these changed should fix the bug part of this issue and I have therefore removed that label.

We should start experimenting with some of the ideas proposed by @jiahao to find the right way of expressing a Hermitian tridiagonal matrix for 0.4.

I'm not super-concerned about the memory required by the extra diagonal, since it is only a constant factor....it's hard for me to imagine a case in which a tridiagonal solver is memory-bound.

No, I think you are right. It will probably not matter for perfomance.

[stevengj]

Okay closing this in favor of JuliaLang/julia#8240 now that the behavior is consistent.