Error in matrix^power if eigenvalues are complex

[alyst]

The M^p code for the generic case (p::Number) assumes that eig() returns real eigenvalues that could be compared with 0. That gives an error if eigenvalues are complex. E.g. (julia 0.4.5):

> cplx_eigvals_mtx = [3.0 -2.0; 4.0 -1.0];

> eigvals(cplx_eigvals_mtx)
2-element Array{Complex{Float64},1}:
 1.0+2.0im
 1.0-2.0im

> cplx_eigvals_mtx^1.5
ERROR: MethodError: `isless` has no method matching isless(::Complex{Float64}, ::Int64)
Closest candidates are:
  isless(::AbstractFloat, ::Real)
  isless(::Real, ::Real)
  isless(::Char, ::Integer)
 in bitcache_lt at broadcast.jl:406
 in .< at broadcast.jl:422
 in ^ at linalg/dense.jl:180

[andreasnoack]

If we revert this, I think we'd get rid of the error.

The matrix power function could use some love. I don't think it is advisable to compute the eigen decomposition in the nonsymmetric case. We should do Algorithm 1 in http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.452.6407&rep=rep1&type=pdf which I think we are already doing for some of the other matrix function. That would also fix this issue.

[alyst]

@andreasnoack Yes, if we replace the line in question (dense.jl:180) by

eltype(v) <: Complex || (any(v.<0) && (v = complex(v)))

the error should go away.

But using "Algorithm 1" would be much better. IIUC, M^p should not contain complex entries if M.>0. It happens currently because of the numerical errors in eig().

[alyst]

See also JuliaLang/julia#5840, JuliaLang/julia#12584

[andreasnoack]

@alyst It would be great if you could prepare a PR.

[alyst]

@andreasnoack There's already JuliaLang/julia#12584, which implements state-of-the-art Schur–Padé algorithm (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.297.2021&rep=rep1&type=pdf). I hope the author is around and that PR could be merged soon.

[andreasnoack]

True. I'd forgotten that PR. Hopefully, @mfasi can be persuaded to finish it. That would be a big improvement over what we have now.

[GunnarFarneback]

I don't think it is advisable to compute the eigen decomposition in the nonsymmetric case.

I agree.

julia> A=[1 0;1 1]
2×2 Array{Int64,2}:
 1  0
 1  1

julia> A^7
2×2 Array{Int64,2}:
 1  0
 7  1

julia> A^7.00001
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> expm(7.00001 * logm(A))
2×2 Array{Float64,2}:
 1.0      0.0
 7.00001  1.0

[mfasi]

@alyst @andreasnoack I am still around and willing to work on the PR. I seem to recollect that the patch was ready to be merged, when I stopped working on it. I will rebase the branch and see what happens.

[alyst]

@mfasi Thanks!

[andreasnoack]

@mfasi That would be really great.

[stevengj]

Fixed by JuliaLang/julia#21184.