^(::Matrix{Float64},::Float64) not type stable

[afniedermayer]

@code_warntype [1. 1.;1. 0.]^2. shows return type Any, same holds for Complex types, i.e. the test

@testset "type stability of ^(::Matrix{Float64},::Float64)" for elty in (Float64, Float32)
    @test Base.return_types(^, (Matrix{elty}, elty)) == [Matrix{elty}]
    @test Base.return_types(^, (Matrix{Complex{elty}}, elty)) == [Matrix{elty}]
end

fails (both on Julia 0.6 and 0.7). This issue is related to but different from JuliaLang/julia#23366 (Float64 vs Int64).

[fredrikekre]

I have thought about this before and I guess it could be reasonable to not rewrap in Hermitian here
https://github.com/JuliaLang/julia/blob/953bfc3da6d3e3bbf70ac65f045558e2ff5f3677/base/linalg/symmetric.jl#L574-L590
and same for the Symmetric case. That should fix it but not sure we wanna do that?

[afniedermayer]

After looking at this in more detail, it looks like a more general issue to me than just Hermitian and Symmetric matrices, the result of ^ might be either Real or Complex depending on the values of the parameters:

julia> [1. 1.;0. 1.]^.5
2×2 Array{Float64,2}:
 1.0  0.5
 0.0  1.0

julia> [1. 1.;0. -1.]^.5
2×2 Array{Complex{Float64},2}:
 1.0+0.0im  0.5-0.5im
 0.0+0.0im  0.0+1.0im

This is inconsistent with how scalars are handled:

julia> 1.^.5
1.0

julia> (-1.)^.5
ERROR: DomainError:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 [1] nan_dom_err at .\math.jl:300 [inlined]
 [2] ^(::Float64, ::Float64) at .\math.jl:699

i.e., you always get a Float64 (or a DomainError).

[afniedermayer]

@fredrikekre , by "rewrapping" do you mean that the input matrix A is of type Symmetric and the return type is also of type Symmetric (sometimes)? It looks like there is not only "rewrapping" but also "wrapping" going on, i.e. the input type doesn't have to be Symmetric to sometimes get a Symmetric return type:

julia> [1. 0.;0. 1.]^2.
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> [1. 0.;0. 1.]^.5
2×2 Symmetric{Float64,Array{Float64,2}}:
 1.0  0.0
 0.0  1.0

[afniedermayer]

This seems almost a dup of #21 (except for symmetry). If the design decision was made that sqrtm(A) should not be type stable, then it's consistent for A^.5 to be type unstable as well. However, it's still confusing that sqrtm and ^.5 behave differently with respect to symmetry:

julia> sqrtm([1. 0;0 1])
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

but

julia> [1. 0.;0. 1.]^.5
2×2 Symmetric{Float64,Array{Float64,2}}:
 1.0  0.0
 0.0  1.0

I guess it would make more sense that both sqrtm and ^0.5 simply return Array{Float64,2}.

[fredrikekre]

Yes, that is what I was suggesting; that we should not re-wrap into Symmetric or Hermitian. We should still use the specialized algorithms of course, but just not re-wrap just before returning.

[andreasnoack]

Part of the problem here is that, compared to the scalar case, it is much more costly to convert to complex before the computation so requiring the user to do that would be unfortunate. The alternative is to always return a complex result. I'd be in favor of that because most non-symmetric matrices will have a complex solution anyway.

[oscardssmith]

Fixed on master.

[KristofferC]

Would be good with a test.

[dkarrasch]

Wait, how is this solved? I get

julia> @code_warntype [1. 1.;1. 0.]^2.
Variables
  #self#::Core.Const(^)
  A::Matrix{Float64}
  p::Float64
  @_4::Union{Nothing, Tuple{Int64, Int64}}
  retmat::Matrix{Float64}
  TT::Type{Float64}
  n::Int64
  i::Int64

Body::Any
...

on master.