incorrect code for Matrix^Number

[stevengj]

There are several issues with the ^(A::AbstractMatrix, p::Number) method, illustrated by this example:

julia> B = [ 3  2
            -5 -3 ]
2×2 Array{Int64,2}:
  3   2
 -5  -3

julia> B^2.0
2×2 Array{Int64,2}:
 -1   0
  0  -1

julia> B^2.01
ERROR: MethodError: no method matching isless(::Complex{Float64}, ::Int64)
Closest candidates are:
  isless(::AbstractFloat, ::Real) at operators.jl:98
  isless(::Real, ::Real) at operators.jl:266
Stacktrace:
 [1] (::Base.LinAlg.##17#18)(::Complex{Float64}) at ./<missing>:0
 [2] broadcast_t(::Function, ::Type{Any}, ::Tuple{Base.OneTo{Int64}}, ::CartesianRange{CartesianIndex{1}}, ::Array{Complex{Float64},1}) at ./broadcast.jl:255
 [3] broadcast_c at ./broadcast.jl:298 [inlined]
 [4] broadcast(::Function, ::Array{Complex{Float64},1}) at ./broadcast.jl:415
 [5] ^(::Array{Int64,2}, ::Float64) at ./linalg/dense.jl:320

• The promotion is wrong for Matrix{Int}^Float64, since it doesn't promote the matrix if the power is integer-valued. The result is type-unstable and, more importantly, potentially wrong (since integers might overflow for large powers).

• It computes the eigenvalues v and checks whether they are negative via any(v.<0), but this check is nonsensical if the eigenvalues are complex.

• Because it relies on the eigenvalue decomposition, it is potentially very inaccurate if the matrix is close to defective (and can fail entirely if it is defective, of course). It would be better to use something like our expm algorithm here.

This is a pretty basic routine, so I was a little surprised to find it to be so buggy.

cc @jiahao and @andreasnoack (who have both touched this code).

[andreasnoack]

Interesting that the MethodError hasn't been reported earlier. The method has had that problem for a while. The real fix is probably JuliaLang/julia#12584 but I don't know if @mfasi has time to finish it.

[stevengj]

Looks like the MethodError was reported in #340.

[iamnapo]

Hello, is anyone working on this?

If not I'd want to try and fix it. I think that @stevengj is right. expm is better that svd. I'm confident that I will be able to provide a PR in the next couple of days.

[JeffBezanson]

@iamnapo Thanks, that would be great. You might want to look at JuliaLang/julia#12584 to see if anything from there should be used.

[andreasnoack]

A rebase went completely wrong in JuliaLang/julia#12584 so the branch is pretty messy. I've cherry picked the relevant commits in https://github.com/JuliaLang/julia/tree/anj/powm so it is probably a good idea to use that branch if you want to work on this issue. It looks like most of the work is done.

[iamnapo]

By the looks of it, logm produces incorrect results too.

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

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

julia> B = [3 2; -5 -3]
2×2 Array{Int64,2}:
  3   2
 -5  -3

julia> expm(logm(B))
2×2 Array{Float64,2}:
 -5.33599  -5.66504
 17.5616   18.4572

(In MATLAB expm(logm(B) is equal to B)

So, I think that if we fix logm, then A^p can simply return expm(p * logm(A)).
Do you agree? Should I turn my attention to logm ?

[stevengj]

That seems like an expensive way to compute A^p

[iamnapo]

@stevengj You are obviously correct. What I meant was that by fixing 1 function:logm, we'll get correct functionality on 2 functions immediately.

[andreasnoack]

@iamnapo If you want to work on ^ then please start out from the branch I linked to above. Fixing logm is a separate issue.

[JeffBezanson]

This is taking longer than anticipated, so moving out of 0.6.0 milestone.