add matrix exponentials #46
Conversation
' dexit exit eqq ! !
|
Would be good with a test as well :) |
|
It would be good with a test to see that the functionality is working properly and doesn't regress in the future. |
|
Added testing, fixed spacing and import issues! |
|
Thanks, it looks like we need to update the CI infrastructure here from Travis to Github Actions. I will do that and then CI should run on this PR. |
|
Is this PR forgotten? I think this seems ready to merge. |
src/Quaternion.jl
Outdated
| a = map(t->t.s, A) | ||
| b = map(t->t.v1, A) | ||
| c = map(t->t.v2, A) | ||
| d = map(t->t.v3, A) | ||
|
|
||
| return [a+im*b c+im*d;-c+im*d a-im*b] |
There was a problem hiding this comment.
Currently this allocates a lot of intermediate arrays. What about something like this to allocate just 1?
| a = map(t->t.s, A) | |
| b = map(t->t.v1, A) | |
| c = map(t->t.v2, A) | |
| d = map(t->t.v3, A) | |
| return [a+im*b c+im*d;-c+im*d a-im*b] | |
| n = size(A, 1) | |
| Ac = Matrix{Complex{T}}(undef, 2n, 2n) | |
| Ac[1:n, 1:n] .= complex.(getfield.(A, :s), getfield.(A, :v1)) | |
| Ac[1:n, n+1:2n] .= complex.(getfield.(A, :v2), getfield.(A, :v3)) | |
| Ac[n+1:2n, 1:n] .= .- conj.(view(Ac, 1:n, n+1:2n)) | |
| Ac[n+1:2n, n+1:2n] .= conj.(view(Ac, 1:n, 1:n)) | |
| return Ac |
| id = Matrix{Quaternion{Float64}}(I, 3, 3) | ||
| d = [Quaternion(randn(4)...) for i in 1:3] | ||
| expd = exp.(d) | ||
| D = exp(Array(Diagonal(d))) |
There was a problem hiding this comment.
Can you include a test for a dense matrix also? Perhaps you could check it using evalpoly.
Co-authored-by: Seth Axen <seth@sethaxen.com>
Co-authored-by: Seth Axen <seth@sethaxen.com>
Codecov Report
@@ Coverage Diff @@
## main #46 +/- ##
=======================================
Coverage ? 41.52%
=======================================
Files ? 3
Lines ? 354
Branches ? 0
=======================================
Hits ? 147
Misses ? 207
Partials ? 0 📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
|
While there could be reasons for a specialized implementation like this (e.g., it looks like this could take advantage of existing complex matrix exponentiation, which is reasonably performant), it appears this PR adds support exclusively for This issue could be addressed more generically by instead (or additionally) addressing JuliaLang #51008. |
Quaternion matrix exponential added using a quick detour to complex representation.