3 ms·

Can the DifferentialEquations package actually solve quaternion equations correctly? Quaternions are non-commutative and I'd hazard to guess that some solver me

by tgb 6y ago

Can the DifferentialEquations package actually solve quaternion equations correctly? Quaternions are non-commutative and I'd hazard to guess that some solver methods would assume commutativity at some point.

spacedome 6y ago
My experience with various linear solvers (which a diffeq solver typicaly relies on) is that if they assume commutivity anywhere at all, which they often do, they almost definitely will not work for quaternions. Even if a derivation of the algorithm can be done with non-commutivity in mind, the implementation typically is not. Many NLA solvers are based on orthogonality transforms, which do not translate directly to quaternions, and even solvers which use only inner products and mat-vec multipication like BiCG do not really work as-is for quaternion matrices. Luckily dense linear solves with LU are fine! One can usually use the complex matrix expansion for quaternions and solve that instead without much difficulty though, and for computing eigenvalues this is perhaps even preferable, as it gives a canonical representation. I have not tried to use DifferentialEquations.jl for quaternion problems, so maybe they have figured some of it out, but it is non-trivial.
- ChrisRackauckas 6y ago
  ODE solvers only use LU
  - spacedome 6y ago
    This is good to know, thanks Chris. I mostly solve sparse PDEs, being able to always use LU makes everything much simpler, especially with quaternions.
LolWolf 6y ago
I would say it's probably worth converting your problem to its linear-algebraic representation and then using the solvers on that. They're quite fast there, and it's probably fairly easy to make a representation map and its inverse work well. (I will say that I don't usually deal with quaternions other than in a linear-algebraic form when doing some robotics stuff, so maybe this is more of a pain than warranted!)
- spacedome 6y ago
  I found this approach to work well for linear systems, here is some rough code I used for the representation map (note the jmag/kmag functions are part of my implementation, not sure what the equivalent is for Quaternions.jl). function cmatrix(Q::AbstractVecOrMat{Quaternion{T}}) where {T} [complex.( real.(Q), imag.(Q)) complex.( jmag.(Q), kmag.(Q)); complex.(-jmag.(Q), kmag.(Q)) complex.( real.(Q), -imag.(Q))] end function qmatrix(C::AbstractMatrix{Complex{T}}) where {T} n, m = size(C) quat.(C[1:n÷2, 1:m÷2], C[1:n÷2, m÷2+1:m]) end

by tgb 6y ago

spacedome 6y ago

My experience with various linear solvers (which a diffeq solver typicaly relies on) is that if they assume commutivity anywhere at all, which they often do, they almost definitely will not work for quaternions. Even if a derivation of the algorithm can be done with non-commutivity in mind, the implementation typically is not. Many NLA solvers are based on orthogonality transforms, which do not translate directly to quaternions, and even solvers which use only inner products and mat-vec multipication like BiCG do not really work as-is for quaternion matrices. Luckily dense linear solves with LU are fine! One can usually use the complex matrix expansion for quaternions and solve that instead without much difficulty though, and for computing eigenvalues this is perhaps even preferable, as it gives a canonical representation. I have not tried to use DifferentialEquations.jl for quaternion problems, so maybe they have figured some of it out, but it is non-trivial.

ChrisRackauckas 6y ago
ODE solvers only use LU
- spacedome 6y ago
  This is good to know, thanks Chris. I mostly solve sparse PDEs, being able to always use LU makes everything much simpler, especially with quaternions.

LolWolf 6y ago

I would say it's probably worth converting your problem to its linear-algebraic representation and then using the solvers on that. They're quite fast there, and it's probably fairly easy to make a representation map and its inverse work well. (I will say that I don't usually deal with quaternions other than in a linear-algebraic form when doing some robotics stuff, so maybe this is more of a pain than warranted!)

spacedome 6y ago
I found this approach to work well for linear systems, here is some rough code I used for the representation map (note the jmag/kmag functions are part of my implementation, not sure what the equivalent is for Quaternions.jl). function cmatrix(Q::AbstractVecOrMat{Quaternion{T}}) where {T} [complex.( real.(Q), imag.(Q)) complex.( jmag.(Q), kmag.(Q)); complex.(-jmag.(Q), kmag.(Q)) complex.( real.(Q), -imag.(Q))] end function qmatrix(C::AbstractMatrix{Complex{T}}) where {T} n, m = size(C) quat.(C[1:n÷2, 1:m÷2], C[1:n÷2, m÷2+1:m]) end