Implemented homogeneous option to Quaternion.to_rotation_matrix

[evbernardes]

Brief description of what is fixed or changed

Added a bool parameter called homogeneous defaulting to False.

When homogeneous == True, a slightly different but equivalent definition of the rotation matrix is given.

Other comments

This is another one in a series of small quality-of-life contributions I'm making to Quaternion.

The "homogeneous" version of the rotation matrix can be slightly better for symbolic manipulation (its formula can be seen here).

To demonstrate this, the test function I implemented shows that the homogeneous version is equivalent to the original non-homogeneous version after being passed to simplify.

As an interesting observation, the homogeneous version also has the property of being equivalent to a rotation + a scaling when the quaternion norm is ignored. If the norm is ignored in the non-homogeneous version, the rotation does not correctly represent anything.

Release Notes

• algebras
  • Implemented homogeneous option to Quaternion.to_rotation_matrix

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [41d90958]       [d505ea97]
     <sympy-1.11.1^0>                 
-         984±3μs          630±1μs     0.64  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-        2.79±0ms      1.18±0.01ms     0.42  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.58±0.01ms         1.72±0ms     0.31  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[sympy-bot]

✅

Hi, I am the SymPy bot (v169). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• algebras
  • Implemented homogeneous option to Quaternion.to_rotation_matrix (#24382 by @evbernardes)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.12.

Click here to see the pull request description that was parsed.

#### Brief description of what is fixed or changed
Added a bool parameter called `homogeneous` defaulting to `False`.

When `homogeneous == True`, a slightly different but equivalent definition of the rotation matrix is given.

#### Other comments
This is another one in a series of small quality-of-life contributions I'm making to `Quaternion`.

The "homogeneous" version of the rotation matrix can be slightly better for symbolic manipulation (its formula can be seen [here](https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation)). 

To demonstrate this, the test function I implemented shows that the homogeneous version is equivalent to the original non-homogeneous version after being passed to `simplify`.

As an interesting observation, the homogeneous version also has the property of being equivalent to a rotation + a scaling when the quaternion norm is ignored. If the norm is ignored in the non-homogeneous version, the rotation does not correctly represent anything.

#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* algebras
  * Implemented homogeneous option to Quaternion.to_rotation_matrix
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[smichr]

Either the parameter description or the docstring proper should probably have a description of what this option does.

[evbernardes]

Either the parameter description or the docstring proper should probably have a description of what this option does.

You're right, I completely forgot about is somehow. I added a line to explain this, but maybe it needs some more work!

[smichr]

(its formula can be seen here)

Can you help me find it on the page? I don't find the word "homogeneous" there. Are you using something that is common nomenclature when working with quaternions?

Also, there is a table of conventions near the discussion of JPL convention. Which method does SymPy use?

[evbernardes]

Can you help me find it on the page?

The formula I implemented can be found in Conversion to and from the matrix representation -> From a quaternion to an orthogonal matrix

I don't find the word "homogeneous" there. Are you using something that is common nomenclature when working with quaternions?

I have seen the word homogeneous before but I can't find it again, maybe we could change the name. I might have read it in some journal article.

Also, there is a table of conventions near the discussion of JPL convention. Which method does SymPy use?

In this website (Derivation of Equations section) we can see how, during the derivation, first we get a matrix in which every element is the product between two quaternion elements, and then we get a matrix similar to the one defined originally in SymPy when supposing that a**2 + b**2 + c**2 + d**2 = 1. The problem with this is that we get that 1 term in the diagonal.

This formula (the original form SymPy) makes sense for numerical purposes because it has less operations, but I don't think it makes much sense to use it for symbolic calculations!

Maybe instead of the a parameter called homogeneous that is False by default, maybe we could instead use a different parameter, True by default, that uses the original formula when True. Don't know how to call it though, I'd have to think a bit about a good parameter name.

[evbernardes]

@smichr As an extra reasoning for using the formula I presented here, declaring the following:

from sympy.abc import a, b, c, d, x, y, z
q = Quaternion(a, b, c, d)
v = Quaternion(0, x, y, z)

And then using the product to find the rotation (supposing |q| = 1):

v_rotated = Matrix((q * v * q.conjugate()).args[1:])

We can symbolically differentiate v_rotated w.r.t. to each variable to find the rotation matrix:

R = Matrix([[
    v_rotated.diff(x),
    v_rotated.diff(y),
    v_rotated.diff(z)]])

Once again, this gives the same formula as the one I'm proposing (minus the normalization):

$$\left[ \begin{matrix} % or pmatrix or bmatrix or Bmatrix or ... a^{2} + b^{2} - c^{2} - d^{2} & - 2 a d + 2 b c & 2 a c + 2 b d \\\ 2 a d + 2 b c & a^{2} - b^{2} + c^{2} - d^{2} & - 2 a b + 2 c d \\ - 2 a c + 2 b d & 2 a b + 2 c d & a^{2} - b^{2} - c^{2} + d^{2} \end{matrix} \right]$$

So it would make the symbolic formula more consistent with the calculation of a rotation by the product with q and q.conjugate(). Maybe we could change the parameter homogeneous to simplified, to imply what I said before, that it's the same as passing the result to the simplify function.

EDIT.: This can also be used to create a code generator for the rotation matrix, that can be copy-pasted to the rotation matrix function implementation:

from sympy.abc import a, b, c, d, x, y, z
q = Quaternion(a, b, c, d)
v = Quaternion(0, x, y, z)

v_rotated = Matrix((q * v * q.conjugate()).args[1:])

R = Matrix([[
    v_rotated.diff(x),
    v_rotated.diff(y),
    v_rotated.diff(z)]])

def replace(line):
    line_new = line.replace("a","q.a")
    line_new = line_new.replace("b","q.b")
    line_new = line_new.replace("c","q.c")
    line_new = line_new.replace("d","q.d")
    return line_new

for i, j in product(range(3), repeat=2):
    if i==j:
        line = f'm{i}{j} = s*({R[i,j]})'
    else:
        line = f'm{i}{j} = 2*s*({R[i,j]/2})'
    print(replace(line))

Generates:

m00 = s*(q.a**2 + q.b**2 - q.c**2 - q.d**2)
m01 = 2*s*(-q.a*q.d + q.b*q.c)
m02 = 2*s*(q.a*q.c + q.b*q.d)
m10 = 2*s*(q.a*q.d + q.b*q.c)
m11 = s*(q.a**2 - q.b**2 + q.c**2 - q.d**2)
m12 = 2*s*(-q.a*q.b + q.c*q.d)
m20 = 2*s*(-q.a*q.c + q.b*q.d)
m21 = 2*s*(q.a*q.b + q.c*q.d)
m22 = s*(q.a**2 - q.b**2 - q.c**2 + q.d**2)

[smichr]

So this option is just giving an equivalent form of the expression: the rhs below.

That would be what simplify(lhs) would give but since it may not actually be the same in general (e.g. if b = 1/(x + 1), for example) I am hesitant to make it simplify=True.

On the page you reference they say that the non-homogeneous form is "more efficient calculation in which the quaternion does not need to be unit normalized". Perhaps the flag could be called normal=False since the elements of the matrix will be expressed over a denominator (and that is what the normal() method would do). The docstring could say that normal=True gives an expression that may be more efficient for symbolic calculations but less so for direct evaluation.

[evbernardes]

That would be what simplify(lhs) would give but since it may not actually be the same in general (e.g. if b = 1/(x + 1), for example) I am hesitant to make it simplify=True.

This is actually the same in general (except, of course, if the norm is zero)!

By simply multiplying and dividing that lone 1 by $a^2 + b^2 + c^2 + d^2$ and assembling the elements, you arrive at the form on the right. It wouldn't be the same if one of the forms was not being normalized at the end, but they are.

I didn't really understand the b = 1/(1+x) example though, what do you mean by that?

[evbernardes]

Just changed the parameter from homogeneous to normal.

As a last step, I'd personally prefer that normal should be True by default! But that's a different discussion :)

[smichr]

what do you mean by that?

[21] simplifies to [22] and this is not the same as [23]

(But for that matter, the normal of [21] is not [23], either.)

But if there is no denominator on any of the arguments then the normal form is consistent with what normal does.

I am not opposed to changing the default value since there is no guarantee of the format of the result, only the equivalency of it. And since it is the same for numbers and better for symbolic expressions, I don't think there would be any complaint. We can see if anyone else has an opinion.

[evbernardes]

Ah alright, you meant in the context of sympy simplifications! I see what you mean!

Alright! In the meantime, I think this is finished! So feel free to merge it or tell me to change the default value whenever you feel like it :)

[evbernardes]

@smichr If you're curious, I found where I had read the word homogeneous, it was in a different page.

[smichr]

Thanks for the ref to the use of "homogeneous". Given some precedence of the term use, I think it makes sense to change normal back to homogeneous. And let's go ahead and try making the default True.

[evbernardes]

Given some precedence of the term use, I think it makes sense to change normal back to homogeneous. And let's go ahead and try making the default True.

I'm not sure about the name anymore tbh, because it can be misinterpreted as a homogeneous 4x4 matrix...

I agree with making it by default though!

[smichr]

ok, if you want to change the function parameter and docstring we can see what fails (if anything) with the change.

[evbernardes]

Switched it back to homogeneous and is now True by default.