Skip to content

Commit

Permalink
Integer matrix exponentiation in schurpow (#51992)
Browse files Browse the repository at this point in the history
This improves type-inference by avoiding recursion, as the `A^p` method
calls `schurpow` if `p` is a float. After this, the result of `schurpow`
is inferred as a small `Union`:

```julia
julia> @inferred Union{Matrix{ComplexF64}, Matrix{Float64}} LinearAlgebra.schurpow([1.0 2; 3 4], 2.0)
2×2 Matrix{Float64}:
  7.0  10.0
 15.0  22.0
```

One concern here might be that for large `p`, the `A^Int(p)` computation
might be expensive by repeated multiplication, as opposed to
diagonalization. However, this may only be the case for really large
`p`, which may not be commonly encountered.

I've added a test, but I'm unsure if `schurpow` is deemed to be an
internal function, and this test is unwise. Unfortunately, the return
type of `A^p` isn't concretely inferred yet as there are too many
possible types that are returned, so I couldn't test for that.
  • Loading branch information
jishnub authored Dec 4, 2023
1 parent fb5f7ad commit b69398a
Show file tree
Hide file tree
Showing 2 changed files with 7 additions and 4 deletions.
4 changes: 2 additions & 2 deletions stdlib/LinearAlgebra/src/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -490,7 +490,7 @@ end
function schurpow(A::AbstractMatrix, p)
if istriu(A)
# Integer part
retmat = A ^ floor(p)
retmat = A ^ floor(Integer, p)
# Real part
if p - floor(p) == 0.5
# special case: A^0.5 === sqrt(A)
Expand All @@ -501,7 +501,7 @@ function schurpow(A::AbstractMatrix, p)
else
S,Q,d = Schur{Complex}(schur(A))
# Integer part
R = S ^ floor(p)
R = S ^ floor(Integer, p)
# Real part
if p - floor(p) == 0.5
# special case: A^0.5 === sqrt(A)
Expand Down
7 changes: 5 additions & 2 deletions stdlib/LinearAlgebra/test/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1028,8 +1028,8 @@ end
@test lyap(1.0+2.0im, 3.0+4.0im) == -1.5 - 2.0im
end

@testset "Matrix to real power" for elty in (Float64, ComplexF64)
# Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014
@testset "$elty Matrix to real power" for elty in (Float64, ComplexF64)
# Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014
#Aa : only positive real eigenvalues
Aa = convert(Matrix{elty}, [5 4 2 1; 0 1 -1 -1; -1 -1 3 0; 1 1 -1 2])

Expand Down Expand Up @@ -1065,6 +1065,9 @@ end
@test (A^(2/3))*(A^(1/3)) A
@test (A^im)^(-im) A
end

Tschurpow = Union{Matrix{real(elty)}, Matrix{complex(elty)}}
@test (@inferred Tschurpow LinearAlgebra.schurpow(Aa, 2.0)) Aa^2
end

@testset "diagonal integer matrix to real power" begin
Expand Down

0 comments on commit b69398a

Please sign in to comment.