LinearAlgebra.eigen fails on Transposed / Adjoint matrices (Julia 1.0.0)

[kragol]

The eigen / eigen! functions of the LinearAlgebra package fail with the Adjoint and Transposed types.

Given a square matrix A (e.g. a 10×10 Array{Float64,2}), both the following statements fail:
E = eigen(A')
E = eigen(transpose(A))
with errors
ERROR: MethodError: no method matching eigen(::Adjoint{Float64,Array{Float64,2}})
and
ERROR: MethodError: no method matching eigen(::Transpose{Float64,Array{Float64,2}})
respectively.

A simple fix could be to define the method
LinearAlgebra.eigen(x::Union{Transpose{T, A}, Adjoint{T, A}}) where {T, A<:AbstractMatrix} = eigen(copy(x)) (credit: Twan Koolen on discourse).

However, it could be more efficient to use LAPACK functions directly computing left eigenvectors (i.e. eigenvectors of the transposed) without making another copy of the matrix. It looks like the wrapper LinearAlgebra.LAPACK.geevx! has such a feature. In eigen.jl, the lower level LAPACK call (line 40 as of tag v1.0.0) is
A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
According to the documentation,
A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'V', 'N', 'N', A)
should compute left eigenvectors instead.

I am very new to Julia (I have been mostly a Matlab programmer so far), so I am not quite sure what is a clean way to insert such a fix into eigen.jl. At least, duplicating the eigen! code in eigen.jl does not seem to be a very good idea. Also, adding this may not cover all eigen methods for special matrices. (The eigen! method in eigen.jl applies to StridedMatrix{T}, which I assume excludes at least sparse matrices.)

Note: this issue is similar to Issue #523 regarding svdfact (now svd) which was fixed in PR JuliaLang/julia#27916 . The fix for svd will not work here though.

[fredrikekre]

It seems like you have already looked into the source code for this. Your suggestions make sense to me, feel free to submit a pull request!

[haampie]

Looks like a great idea! Note that there is still an edge case: the transpose (not adjoint) of a complex-valued matrix: transpose(rand(ComplexF64, 10, 10)). I don't think LAPACK supports it, but you could still compute the left eigenvectors and eigenvalues, and then afterwards replace each value with its complex conjugate.

[kragol]

I have been experimenting with some solutions. It's easy enough to modify the code for eigen!(::StridedMatrix{T};...) to extract left eigenvectors. I wouldn't say my solution is very elegant as I end up with 3 copies of pretty much the same code for StridedMatrix{T}, Adjoint{T,StridedMatrix{T}} and Transpose{T,StridedMatrix{T}}, with only minor variations like taking a conjugate here or there. There must be a way to define these functions without writing the same code 3 times (using macros or some auxiliary functions?) but I am not sure what is the best approach in Julia so that it does not impact performance and keeps the code readable. My understanding is that macros should have zero performance impact but I still need to get used to using them.

I still need to go through the dozen specialized eigen! methods to make sure they all work with Adjoint/Transpose types. Now the big question is: is it worth it?. On the one hand, duplicating all eigen! methods seems like it will increase the maintenance burden on the code with only a limited performance benefit. On the other hand, writing a generic eigen! method for Adjoint/Transpose types that makes a copy of the matrix pretty much ruins the purpose of the eigen! methods.

Besides, it seems that there will be pretty much zero benefit for the more common eigen methods since they make a copy of the matrix in the first place. For these, it seems that, instead of calling specialized eigen! methods for Adjoint/Transpose types, it would be much cleaner to write a single method like

LinearAlgebra.eigen(x::Union{Transpose{T, A}, Adjoint{T, A}}) where {T, A<:AbstractMatrix} = eigen(copy(x))

(ignoring optional keyword arguments here)

I will keep experimenting and do some benchmarks to assess the possible performance benefits of avoiding copying the matrix when possible. We will see from there if it is worth going for a complex solution involving many additional methods.

Finally, I have only been looking at eigen(A...) methods so far. The same should be done for eigen(A, B...) methods for the generalized spectral decomposition problem. This will take a while...

Would it make sense to do a first PR fixing the eigen/ eigen! methods using a simple copy strategy (as above)? This issue has not been tagged as a bug, but it certainly looks like one from a user perspective. Adding specialized eigen! methods for better performance could be done later once the bug is fixed.

[jleugeri]

Any news on this?
IMO something as simple as eigen(A') should not throw a nondescript error message - not for a language that wants to displace Matlab.

eigen(copy(A')) seems to be an acceptable workaround, but it looks stupid in user code and unless you know that

1. A' = Adjoint(A) is a lazy wrapper around A (for which no method is implemented)
2. copy doesn't just copy, but also turns the Adjoint into an Array

people won't come up with that workaround themselves when they see that error message or understand it, when they see it in code.
Since people shouldn't have to dig through github issues to find a fix for a "backend"-problem as that, could we just go with the solution suggested by @kragol for now, with one minor change of using eigen! instead of eigen, since we're making a copy anyways? IE:

LinearAlgebra.eigen(x::Union{Transpose{T, A}, Adjoint{T, A}}) where {T, A<:AbstractMatrix} = eigen!(copy(x))

Performance issues could then be dealt with later via updates to LinearAlgebra.
Once everyone hard-codes the workaround eigen!(copy(A')) into their code (as I'll have to do for now), that'll be impossible!

If @kragol doesn't call dibs, I could create a PR for that one-liner.

[kragol]

@jleugeri Feel free to create such a PR. I did not do it because I thought that was a really poor solution and I was looking for something a bit better (see my previous post). Maybe I was wrong though, and adding a poor man's fix would be the right thing to do now.

While digging through the code, I realised, that fixing that eigen bit properly would require fixing quite a few functions in the eigen suite, which was a bit daunting for a newb like me. At that point I looked at other functions in the LinearAlgebra package and noticed that many actually did not work with Transpose and / or Adjoint types. My conclusion was that the Adjoint / Transpose types were basically half-baked (they were like 3 months old at that point) and that a larger effort was necessary to fix the whole package. So, I gave up on julia at the time, thinking that despite the 1.0 mark, it was not ready for my daily work and I should better wait for the language and its ecosystem to mature.

As you can see, there has not been much progress regarding these issues in the last 6 months. As of now, it looks like several LinearAlgebra functions are still broken with Adjoint/Transpose types (e.g. cholesky, lq, schur), not to mention issues with other matrix types like Bidiagonal, Tridiagonal, SymTridiagonal, etc (see open github issues) . While there might have been some improvements since August, I still don't feel like it is ready to displace matlab. (although to be really fair I believe that the main reason why that won't happen right now is more like the whole matlab user experience is still much more polished and accessible, especially when it comes to plotting stuff)

Edit: I spoke too fast. Thanks @andreasnoack for addressing the issue (I would not have been able to devise such a fix). While you are looking at it, you may want to also fix the eigen(A,B) methods for generalized eigenvalue decomposition or even some other functions like cholesky, lq or schur.

[jleugeri]

Splendid! Thanks @kragol for looking at this and thanks @andreasnoack for the swift response. I didn't even have time yet to look at this and you already have a better PR than what I'd have ham-fisted in there. ( Side-remark regarding the polished user experience of Matlab: yes, until people start using strings & structured arrays or MEX files and all hell breaks loose 😖 )