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

added missing matrix exponentiation methods #29782

Merged
merged 6 commits into from
Oct 30, 2018

Conversation

ExpandingMan
Copy link
Contributor

@ExpandingMan ExpandingMan commented Oct 23, 2018

Fixes JuliaLang/LinearAlgebra.jl#579.

I added some tests for existing methods exp. exp is missing methods for Rational and Irrational; I put in commented-out tests for those. I'm not sure if inserting speculative tests like this is an acceptable practice, I can take out those comments if not.

@@ -89,6 +89,7 @@ catalan
# loop over types to prevent ambiguities for ^(::Number, x)
for T in (AbstractIrrational, Rational, Integer, Number, Complex)
Base.:^(::Irrational{:ℯ}, x::T) = exp(x)
Base.:^(::Irrational{:ℯ}, x::AbstractMatrix{T}) = exp(x)
Copy link
Member

@Keno Keno Oct 23, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The loop says it's for ambiguities. Shouldn't we instead have a general
Base.:^(::Irrational{:ℯ}, x) = exp(x) fallback?

Copy link
Contributor Author

@ExpandingMan ExpandingMan Oct 23, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know, wouldn't that still allow for the ambiguities of ^(::Number, x)? I would have thought it can create an ambiguity if we have the fallback, even if we have the specific methods.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In that case, maybe just add AbstractMatrix to the thing that T loops over?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be defined in LinearAlgebra.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't want to add AbstractMatrix as T to the loop because I thought it would be safer (in regard to method ambiguities) to ensure that only matrices with these element types have exponentials defined for them. Having said that, I just realized that this should have read AbstractMatrix{<:T} rather than what I have here.

Let me know whether you think this should indeed be fore general AbstractMatrix and I'll move it to LinearAlgebra.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since exp(::AbstractMatrix) is defined in LinearAlgebra this should too.

@ExpandingMan
Copy link
Contributor Author

Ok, wasn't quite sure where to put this. Most of the methods are in dense.jl but that didn't necessarily seem appropriate. If you want me to move it please let me know.

By the way, I don't see a exp(::AbstractMatrix) method, though there are some StridedArray methods in dense.jl.

The tests seem a bit awkward where I put them, but I didn't know where else they would go.

@fredrikekre fredrikekre added the linear algebra Linear algebra label Oct 24, 2018
@@ -340,6 +340,9 @@ rdiv!(A, B)
copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

# matrix exponentials for ℯ
Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would move this just below

exp(A::StridedMatrix{<:Union{Integer,Complex{<:Integer}}}) = exp!(float.(A))

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are you sure about that? the reason I didn't put it in there is because it is only for dense, it could work for sparse array types as well.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But it doesn't, this is the generic exp so I would put it here.

@@ -401,6 +401,10 @@ end
[4.000000000000000 -1.414213562373094 -1.414213562373095
-1.414213562373095 4.999999999999996 -0.000000000000000
0 -0.000000000000002 3.000000000000000])

# alias using ℯ, should be exact
@test ℯ^(fill(2, (4,4))) == exp(fill(2, (4,4)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would also move these test out the loop over the elty, maybe to a separate sub-testset.

@dalum
Copy link
Contributor

dalum commented Oct 24, 2018

I've had something similar on my list and was actually planning to open a PR for this just now. 😄 I think it would be really nice for consistency to also have Base.:^(x::Number, A::AbstractMatrix) = exp(log(x)*A), though. What do you think? The one specialized on \euler should still be there to avoid the overhead allocation of multiplying by 1.

@ExpandingMan
Copy link
Contributor Author

Yeah, good point, I don't know how I forgot about the general method (i.e. for Number).

Ok, the tests I wrote yesterday were awful, I've cleaned them up a bit so that they test more eltypes and they are in their own little section.

Are you really sure you want these methods in dense.jl? They are in no way specific to dense matrices, so I thought it would be confusing to put them there. If that's still what you want just let me know and I'll move it.

@ExpandingMan
Copy link
Contributor Author

Not sure what error travis had on MacOS. I don't know if it's just bad luck, but I feel like every time I make a PR I wind up with bogus failures on MacOS.

@ExpandingMan
Copy link
Contributor Author

Is any action needed here? I'm guessing this is probably not a real error on MacOS (as this has happened multiple times before and linux passes), can we re-trigger somehow?

@stevengj
Copy link
Member

Seems to be passing now?

@ExpandingMan
Copy link
Contributor Author

Ok, I don't know what the deal was with that thing. The MacOS seems to have had bogus failures every time I've ever made a PR.

@@ -489,6 +489,10 @@ julia> exp(A)
exp(A::StridedMatrix{<:BlasFloat}) = exp!(copy(A))
exp(A::StridedMatrix{<:Union{Integer,Complex{<:Integer}}}) = exp!(float.(A))

Base.:^(b::Number, A::AbstractMatrix) = exp(log(b)*A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could actually call exp!(log(b)*A) since the multiplication returns a new array.

@StefanKarpinski StefanKarpinski merged commit de278a5 into JuliaLang:master Oct 30, 2018
mortenpi added a commit to mortenpi/julia that referenced this pull request Dec 1, 2018
fredrikekre added a commit that referenced this pull request Dec 5, 2018
changes between Julia 1.0 and 1.1, including:

- Custom .css-style for compat admonitions.

- Information about compat annotations to CONTRIBUTING.md.

- NEWS.md entry for PRs #30090, #30035, #30022, #29978,
  #29969, #29858, #29845, #29754, #29638, #29636, #29615,
  #29600, #29506, #29469, #29316, #29259, #29178, #29153,
  #29033, #28902, #28761, #28745, #28708, #28696, #29997,
  #28790, #29092, #29108, #29782

- Compat annotation for PRs #30090, #30013, #29978,
  #29890, #29858, #29827, #29754, #29679, #29636, #29623,
  #29600, #29440, #29316, #29259, #29178, #29157, #29153,
  #29033, #28902, #28878, #28761, #28708, #28156, #29733,
  #29670, #29997, #28790, #29092, #29108, #29782, #25278

- Documentation for broadcasting CartesianIndices (#30230).
- Documentation for Base.julia_cmd().
- Documentation for colon constructor of CartesianIndices (#29440).
- Documentation for ^(::Matrix, ::Number) and ^(::Number, ::Matrix).

- Run NEWS-update.jl.

Co-authored-by: Morten Piibeleht <[email protected]>
Co-authored-by: Fredrik Ekre <[email protected]>
fredrikekre added a commit that referenced this pull request Dec 5, 2018
changes between Julia 1.0 and 1.1, including:

- Custom .css-style for compat admonitions.

- Information about compat annotations to CONTRIBUTING.md.

- NEWS.md entry for PRs #30090, #30035, #30022, #29978,
  #29969, #29858, #29845, #29754, #29638, #29636, #29615,
  #29600, #29506, #29469, #29316, #29259, #29178, #29153,
  #29033, #28902, #28761, #28745, #28708, #28696, #29997,
  #28790, #29092, #29108, #29782

- Compat annotation for PRs #30090, #30013, #29978,
  #29890, #29858, #29827, #29754, #29679, #29636, #29623,
  #29600, #29440, #29316, #29259, #29178, #29157, #29153,
  #29033, #28902, #28878, #28761, #28708, #28156, #29733,
  #29670, #29997, #28790, #29092, #29108, #29782, #25278

- Documentation for broadcasting CartesianIndices (#30230).
- Documentation for Base.julia_cmd().
- Documentation for colon constructor of CartesianIndices (#29440).
- Documentation for ^(::Matrix, ::Number) and ^(::Number, ::Matrix).

- Run NEWS-update.jl.

Co-authored-by: Morten Piibeleht <[email protected]>
Co-authored-by: Fredrik Ekre <[email protected]>
fredrikekre added a commit that referenced this pull request Dec 5, 2018
Addition of NEWS and compat admonitions for important changes between Julia 1.0 and 1.1, including:

- Custom .css-style for compat admonitions.

- Information about compat annotations to CONTRIBUTING.md.

- NEWS.md entry for PRs #30090, #30035, #30022, #29978,
  #29969, #29858, #29845, #29754, #29638, #29636, #29615,
  #29600, #29506, #29469, #29316, #29259, #29178, #29153,
  #29033, #28902, #28761, #28745, #28708, #28696, #29997,
  #28790, #29092, #29108, #29782

- Compat annotation for PRs #30090, #30013, #29978,
  #29890, #29858, #29827, #29754, #29679, #29636, #29623,
  #29600, #29440, #29316, #29259, #29178, #29157, #29153,
  #29033, #28902, #28878, #28761, #28708, #28156, #29733,
  #29670, #29997, #28790, #29092, #29108, #29782, #25278

- Documentation for broadcasting CartesianIndices (#30230).
- Documentation for Base.julia_cmd().
- Documentation for colon constructor of CartesianIndices (#29440).
- Documentation for ^(::Matrix, ::Number) and ^(::Number, ::Matrix).

- Run NEWS-update.jl.


Co-authored-by: Morten Piibeleht <[email protected]>
Co-authored-by: Fredrik Ekre <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

missing matrix exponentiation methods for \euler
6 participants