Add matrix exponential for 1x1 and 2x2, and fallback to Base (fixes #159) #160

Godisemo · 2017-05-10T17:18:41Z

I have only tried this in julia v0.5, which is why the PR is to the julia-0.5 branch. You could probably cherry pick it to master.

Note: On my computer this gives a speedup of over 300x

a = randn(SMatrix{2,2})
b = Array(a)

@benchmark expm(a)
# BenchmarkTools.Trial:
#  memory estimate:  48 bytes
#  allocs estimate:  1
#  --------------
#  minimum time:     57.861 ns (0.00% GC)
#  median time:      58.580 ns (0.00% GC)
#  mean time:        65.692 ns (8.97% GC)
#  maximum time:     3.145 μs (96.47% GC)
#  --------------
#  samples:          10000
#  evals/sample:     983
@benchmark expm(b)
# BenchmarkTools.Trial:
#  memory estimate:  5.67 KiB
#  allocs estimate:  76
#  --------------
#  minimum time:     15.294 μs (0.00% GC)
#  median time:      18.114 μs (0.00% GC)
#  mean time:        19.030 μs (4.38% GC)
#  maximum time:     4.279 ms (98.87% GC)
#  --------------
#  samples:          10000
#  evals/sample:     1

andyferris

Thanks @Godisemo! This looks nice, and is appreciated.

I think that the functionality here is OK to merge (there are some necessary fixes in the code comments).

In the future it would be nice to leverage eig or other decompositions to do this... e.g. since we can do eig for symmetric/hermitian 3x3 matrices then we extend expm for symmetric/hermitian 3x3 matrices. (Arguably we could define all the functions which work on in independent eigenvalues like expm, logm, etc in terms of eig, and add specializations to that over time). To be clear, I don't think this is necessary to be part of this PR, though.

andyferris · 2017-05-10T23:28:40Z

src/expm.jl

+@inline expm(A::StaticMatrix) = _expm(Size(A), A)
+
+@generated function _expm(::Size{(1,1)}, A::StaticMatrix)
+    @assert size(A) == (1,1)


Probably overkill, given the above.

you mean the assertion?

andyferris · 2017-05-10T23:45:42Z

src/expm.jl

+
+@generated function _expm(::Size{(1,1)}, A::StaticMatrix)
+    @assert size(A) == (1,1)
+    T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)


In Julia v0.6, you can only perform "hyper"-pure operations inside the function generator. Since the eltype might not be defined before static arrays is loaded, this just isn't "pure enough" for the assumptions made by the Julia 0.6 compiler.

However, it's safe to move this into the body of the function (inside the quote) - where possible inference will perform type calculations at compile time.

Also - why sqrt(one) not exp(zero)? (zero is slightly more flexible from a group vs. algebra perspective, i.e. in theory eltype(A) might be a generator that supports addition and zero only.)

Moving

T = promote_type(typeof(exp(zero(eltype(A)))), Float32) newtype = similar_type(A,T)

into the body removes the need for using a generated functions in the first place, right?

actually, yes! I think the eig and inv code works like that.

andyferris · 2017-05-10T23:46:01Z

src/expm.jl

+@generated function _expm(::Size{(1,1)}, A::StaticMatrix)
+    @assert size(A) == (1,1)
+    T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)
+    newtype = similar_type(A,T)


Also move this to body - it's "free"

andyferris · 2017-05-10T23:46:13Z

src/expm.jl

+@generated function _expm{S<:Real}(::Size{(2,2)}, A::StaticMatrix{S})
+    @assert size(A) == (2,2)
+    T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)
+    newtype = similar_type(A,T)


Move these two lines to the body.

andyferris · 2017-05-10T23:47:43Z

src/expm.jl

+        v = (a-d)^2 + 4*b*c
+
+        # if v == 0
+        z1::$T = 1.0


These are unnecessary since if statements don't introduce a new scoping block.

The reason for this was to ensure type stability. Moved it to the else-block instead.

andyferris · 2017-05-10T23:48:15Z

src/expm.jl

+
+@generated function _expm(::Size, A::StaticArray)
+  T = promote_type(typeof(sqrt(one(eltype(A)))), Float32)
+  newtype = similar_type(A,T)


Move these two lines to the generated body

andyferris · 2017-05-10T23:53:15Z

src/expm.jl

+
+  quote
+      $(Expr(:meta, :inline))
+      ($newtype)(Base.expm(Array(A)))


Generally, I've been using SizedArray for this purpose with eig, inv, etc. The static size of the array is not forgotten, the data doesn't need to be copied to the stack unnecessarily, and the user gets feedback that a fully "native" static arrays method was not used.

However, it does violate similar_type, which doesn't seem that desirable...

andyferris · 2017-05-11T00:43:11Z

src/expm.jl

+@inline expm(A::StaticMatrix) = _expm(Size(A), A)
+
+@inline function _expm(::Size{(1,1)}, A::StaticMatrix)
+    T = promote_type(typeof(exp(zero(eltype(A)))), Float32)


Why are we promoting with Float32 here?

I probably copied from the cholesky source

andyferris · 2017-05-11T00:43:41Z

src/expm.jl

+
+@inline function _expm(s::Size, A::StaticArray)
+    T = promote_type(typeof(exp(zero(eltype(A)))), Float32)
+    newtype = similar_type(A,T)


newtype unused

andyferris · 2017-05-11T00:44:38Z

src/expm.jl

+        z1 = cos(z / 2)
+        z2 = sin(z / 2) / z
+    else # if v == 0
+        z1 = T(1.0)


~~zero(T)~~ one(T)

Or maybe not. I dunno.

since there is no half(T), I would prefer not to use one(T)

andyferris · 2017-05-11T00:47:58Z

src/expm.jl

+    (newtype)((exp(A[1]), ))
+end
+
+@inline function _expm{S<:Real}(::Size{(2,2)}, A::StaticMatrix{S})


StaticMatrix has new type parameters, StaticMatrix{M,N,T}.

Try _expm(::Size{(2,2)}, A::StaticMatrix{<:Any, <:Any, <:Real})

Then it won't work on 0.5 though. How to handle backwards compatibility?

How to handle backwards compatibility?

We don't :) They have diverged.

It's not an ideal situation, but that's where we are.

I think there is a branch somewhere where we can merge backports. It's just more work.

Ok, up until the last commit, it worked on julia 0.5. I think I might have messed up when changing target branch for this PR to master.

I see... there seems to be some conflicts.

One way of doing this one PR to julia-0.5 and one PR to master. The difference will be trivial.

yeah, I changed back to the julia-0.5 branch. Will copy the files over to master and create a new PR. You will have to revert the last commit though and pull 1f2ff2b

This reverts commit f8e1de2.

Godisemo · 2017-05-22T11:45:40Z

Is this PR and #161 not ready for merging?

andyferris · 2017-05-23T00:08:10Z

Travis is showing errors...

Godisemo · 2017-05-23T00:12:55Z

You can clearly see that the errors are caused by previous commits... The tests that fail are related to Binary IO and write. As for #161 I don't see how that error is related to these changes either.

andyferris · 2017-05-23T01:59:06Z

I'm totally confused.

I re-ran travis... what exactly are these errors about?

c42f · 2017-05-23T03:03:38Z

Yes, these errors are really weird. Looks like I forgot a Compat dependency for take! when doing the IO tests... but as to why the tests pass locally I'm really confused.

coveralls · 2017-05-23T05:24:44Z

Changes Unknown when pulling 72fff23 on Godisemo:matrix-exponential into ** on JuliaArrays:julia-0.5**.

c42f · 2017-05-23T05:26:46Z

After merging #174 the test failures here are now only spurious failures on nightly - tests for julia-0.5 pass so I'm going to merge this.

Add matrix exponential for 1x1 and 2x2, and fallback to Base

963b530

Godisemo mentioned this pull request May 10, 2017

Matrix exponential #159

Closed

Godisemo changed the title ~~Add matrix exponential for 1x1 and 2x2, and fallback to Base (issue #159)~~ Add matrix exponential for 1x1 and 2x2, and fallback to Base (fixes #159) May 10, 2017

Godisemo added 2 commits May 10, 2017 19:49

Remove use of complex numbers in expm

86a565f

Rename type parameter

6820f2d

andyferris requested changes May 10, 2017

View reviewed changes

Fix code review comments

2071dff

andyferris reviewed May 11, 2017

View reviewed changes

Godisemo added 3 commits May 11, 2017 02:48

Remove unused code

27a718f

Don't promote with float32

1f2ff2b

Update type parameters of staticmatrix

f8e1de2

Godisemo changed the base branch from julia-0.5 to master May 11, 2017 01:08

Godisemo changed the base branch from master to julia-0.5 May 11, 2017 01:15

Godisemo mentioned this pull request May 11, 2017

Add matrix exponential for 1x1 and 2x2, and fallback to Base #161

Merged

Revert "Update type parameters of staticmatrix"

72fff23

This reverts commit f8e1de2.

c42f mentioned this pull request May 23, 2017

Use Compat for take!() #174

Merged

c42f merged commit c353641 into JuliaArrays:julia-0.5 May 23, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add matrix exponential for 1x1 and 2x2, and fallback to Base (fixes #159) #160

Add matrix exponential for 1x1 and 2x2, and fallback to Base (fixes #159) #160

Godisemo commented May 10, 2017 •

edited

Loading

andyferris left a comment

andyferris May 10, 2017

Godisemo May 11, 2017

andyferris May 11, 2017

andyferris May 10, 2017

Godisemo May 11, 2017

andyferris May 11, 2017

andyferris May 10, 2017

andyferris May 10, 2017

andyferris May 10, 2017

Godisemo May 11, 2017

andyferris May 10, 2017

andyferris May 10, 2017

andyferris May 11, 2017

Godisemo May 11, 2017

andyferris May 11, 2017

andyferris May 11, 2017 •

edited

Loading

Godisemo May 11, 2017

andyferris May 11, 2017

Godisemo May 11, 2017

andyferris May 11, 2017

Godisemo May 11, 2017

andyferris May 11, 2017

andyferris May 11, 2017

Godisemo May 11, 2017

andyferris May 11, 2017

Godisemo May 11, 2017

Godisemo commented May 22, 2017

andyferris commented May 23, 2017

Godisemo commented May 23, 2017 •

edited

Loading

andyferris commented May 23, 2017

c42f commented May 23, 2017

coveralls commented May 23, 2017

c42f commented May 23, 2017

Add matrix exponential for 1x1 and 2x2, and fallback to Base (fixes #159) #160

Add matrix exponential for 1x1 and 2x2, and fallback to Base (fixes #159) #160

Conversation

Godisemo commented May 10, 2017 • edited Loading

andyferris left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andyferris May 11, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Godisemo commented May 22, 2017

andyferris commented May 23, 2017

Godisemo commented May 23, 2017 • edited Loading

andyferris commented May 23, 2017

c42f commented May 23, 2017

coveralls commented May 23, 2017

c42f commented May 23, 2017

Godisemo commented May 10, 2017 •

edited

Loading

andyferris May 11, 2017 •

edited

Loading

Godisemo commented May 23, 2017 •

edited

Loading