implement Tridiagonal, Bidiagonal, and Triangular eigensolvers

[stevengj]

Currently, eig and similar throw no method eigfact when applied to these types of matrices. We should be able to provide efficient LAPACK-based solvers for these types.

For Tridiagonal, we can exploit the fact that it is already in Hessenberg form to call xHSEQR directly.

In the case of Bidiagonal matrices, the eigenvalues can be read off the diagonal with no computation, and eigenvectors can be obtained in O(m^2) time (one linear-time solve per eigenvector).

Similarly, any Triangular matrix is already in Schur form (or its transpose), allowing us to read off the eigenvalues from the diagonal and to compute the eigenvectors by calling xTREVC.

• Bidiagonal
  • eig
  • eigvals
  • eigvecs
• Tridiagonal
• SymTridiagonal
  • eig
  • eigvals
  • eigvecs
• Triangular
  • eig
  • eigvals
  • eigvecs

[JeffBezanson]

+1

[ViralBShah]

@andreasnoackjensen - In case you have some spare cycles...

[jiahao]

I've implemented eigvecs for lower Triangular based on the fact that its eigenvectors are the inverse transpose of the eigenvectors of its transpose (upper Triangular) analogue. Not sure if there is a better way to do this.

[stevengj]

Can't you just compute the left eigenvectors of its transpose?

[jiahao]

Right, I forgot that xHSEQR can compute left eigenvectors. This is now implemented in 3e255d4.

For Bidiagonal eigenvectors, thanks for reminding me about the case where the off-diagonal is 0. Since this effectively breaks up the Bidiagonal into submatrices where the current method should work, I'll modify the method to calculate elements of the eigenvectors in the pertinent subblocks.

[stevengj]

See also JuliaLang/julia#5255

[jiahao]

Apparently every time I rebase a commit referencing an issue, Github creates a new cross-reference. Oopsie.

[jiahao]

bump
@andreasnoackjensen, did you say you had code for the tridiagonal problem? If not, I'll look into it.

[andreasnoack]

Sorry I haven't made a pull request yet. I have spent some time on preparing some QR code and I plan to open a pull request with later today. When that is code I'll prepare the eigensolver code.

[jakebolewski]

@andreasnoack can this issue be closed?

[andreasnoack]

Not yet. The tridiagonal case is still to be implemented.

[stevengj]

(Sorry, I confused Tridiagonal with SymTridiagonal.)

I think that eigvals may be suboptimal for SymTridiagonal in some cases where sterf is faster than stegr. I get a factor of 2 speedup from

for (sterf,elty) in ((:dsterf_,Float64), (:ssterf_,Float32))
    #   DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
    #   using the Pal-Walker-Kahan variant of the QL or QR algorithm. 
    #   Arguments:
    #      N       (input) INTEGER -- The order of the matrix.  N >= 0.
    #      D       (input/output) DOUBLE PRECISION array, dimension (N)
    #              On entry, the n diagonal elements of the tridiagonal matrix.
    #              On exit, if INFO = 0, the eigenvalues in ascending order.
    #      E       (input/output) DOUBLE PRECISION array, dimension (N-1)
    #              On entry, the (n-1) subdiagonal elements of the tridiagonal matrix.
    #              On exit, E has been destroyed.
    #      INFO    (output) INTEGER -- exit code
    @eval function sterf!(A::SymTridiagonal{$elty})
        info = Array(BlasInt, 1)
        n = size(A, 2)
        ccall(($(blasfunc(sterf)), liblapack), Void,
              (Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}), &n, A.dv, A.ev, info)
        # @lapackerror
        return A.dv
    end
end

n = 1000
A = SymTridiagonal(rand(n), rand(n-1))
@time sterf!(copy(A))
@time eigvals(A);
@time eigvals(full(A));

gives

elapsed time: 0.024711461 seconds (16296 bytes allocated)
elapsed time: 0.041806028 seconds (234488 bytes allocated)
elapsed time: 0.138864358 seconds (56377560 bytes allocated)

I'm not sure if there is some other reason to prefer stegr. The dstevx driver routine defaults to dsterf if only (all of) the eigenvalues are desired.

[stevengj]

(From an eigenvalues perspective, it would be nicer to have a HermitianTridiagonal type; a complex SymTridiagonal matrix is kind of useless. Such a matrix is similar to a real SymTridiagonal matrix via a diagonal unitary matrix, so supporting it is straightforward in principle.)

Oh, I see, SymTridiagonal is actually supposed to be Hermitian. It is just that there appears to be a printing bug: SymTridiagonal([1,2],[3im]) prints as symmetric, not Hermitian, but converting it via full produces a Hermitian matrix. But shouldn't the diagonal be constrained to be real? And the constructor checks equality of the upper and lower diagonals. I'm confused. See #157

[andreasnoack]

Yes, the sterf algorithm is the fastest when only the values are needed. I'll change that. At some point, I'd like to add an option for algorithm, but that can wait.

[andreasnoack]

@stevengj I have now changed the default solver for SymTridiagonal when all values are asked for to sterf (via stev) on 0.4. The stegr solver has the advantage that you can ask for specific parts of the spectrum so it is still used for finding the minimum and maximum eigenvalues.

[stevengj]

Thanks. (Which commit?)

[andreasnoack]

Oops. I only pushed it to anj/pwk. I'll merge when the lights are green.

[andreasnoack]

@stevengj JuliaLang/julia@06f9420

[vtjnash]

bump. what's the status on this? it's on the v0.4.x list, but looks like it has not had any work done this year. that seems to imply it is not a current high priority?

[andreasnoack]

I think we are actually good here since the original wish for solvers for Tridiagonal should have been SymTridiagonal according to #17.