Hacker News Reader

Comment by tgb | Hacker News Reader

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.

View on HN

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.

ODE solvers only use LU

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.

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!)

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

Reply on news.ycombinator.com