Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Matrix Powers are type-unstable: sometime Real sometimes Complex #35667

Open
oxinabox opened this issue Apr 30, 2020 · 1 comment
Open

Matrix Powers are type-unstable: sometime Real sometimes Complex #35667

oxinabox opened this issue Apr 30, 2020 · 1 comment
Labels

Comments

@oxinabox
Copy link
Contributor

There are all kinds cases where for some matrix powers the result will be a real matrix,
and for others a complex matrix.
Such as here

function ^(A::Symmetric{<:Real}, p::Real)
isinteger(p) && return integerpow(A, p)
F = eigen(A)
if all-> λ 0, F.values)
return Symmetric((F.vectors * Diagonal((F.values).^p)) * F.vectors')
else
return Symmetric((F.vectors * Diagonal((complex(F.values)).^p)) * F.vectors')
end
end

For other functions like log it is required that you give suitbale types.
i.e. log(-1) errors, you need to do log(complex(-1))
We should be demanding the same thing from matrix powers.
Give a complex input if you want a complex output.

This was introduced back in #21184


Turnout this is really annoying for testing AD.
Making it error would not make it less annoying to test but rather remove the need to test it in the first place.

It breaks FiniteDifferences.jl (JuliaDiff/FiniteDifferences.jl#80) because FIniteDIfferences.jl really can'y handle when the number of inputs needing to be perterbed changes depending on the value of the input.
If it is a complex output there are twice as many things needed to perturbe.

It also really complicates [Zygote's tests](e.g https://github.com/FluxML/Zygote.jl/blob/ac4f1a0727d860b31197a336a02d04b33cb21219/src/lib/array.jl#L607)

@sethaxen
Copy link
Contributor

An example:

julia> using LinearAlgebra, Test

julia> x = [1.0 2.0; 3.0 4.0]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> @inferred x^1 # integer power is fine
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> @inferred x^1.0 # float power is type-unstable
ERROR: return type Array{Float64,2} does not match inferred return type Any
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] top-level scope at REPL[23]:1

julia> @inferred x^complex(1.0) # complex power is type-stable
2×2 Array{Complex{Float64},2}:
 1.0-2.35922e-16im  2.0+5.55112e-17im
 3.0-8.32667e-17im  4.0-3.33067e-16im

Another example is eigendecomposition:

julia> @inferred eigen(x)
ERROR: return type Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}} does not match inferred return type Union{Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}, Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}}
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] top-level scope at REPL[21]:1

@oxinabox oxinabox changed the title Matrix Powers are type-unstable: sometime Real somtimes Complex Matrix Powers are type-unstable: sometime Real sometimes Complex May 9, 2020
@brenhinkeller brenhinkeller added the domain:linear algebra Linear algebra label Nov 20, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants