Re: [math-fun] Robert Smith's vector problem...
I didn't understand your "#" operation. Could you please try again? At 10:46 AM 6/26/2014, Robert Smith wrote:
Thanks for the info.
I agree that a quaternion multiply is not sufficient, unless you define quaternion-vector "multiplication" by the operator # above.
Robert
On Thu, Jun 26, 2014 at 6:07 AM, Henry Baker <hbaker1@pipeline.com> wrote:
I did provide you with a single quaternion s.t.
Q A (1/Q) = A' Q B (1/Q) = B'
However, one normally (pun intended!) does the computation like this:
(Q A conjugate(Q))/N(Q) = A' (Q B conjugate(Q))/N(Q) = B'
where N(Q) = Q conjugate(Q)
because it is cheaper/easier/more precise to normalize at the end.
If you compute the _unit_ quaternion Q/sqrt(N(Q)), then in this case,
(1/Q) = conjugate(Q)
There is no way to do the rotation with only one quaternion multiply, because each quaternion multiply only computes a _reflection_, and you need 2 reflections to make a rotation.
See Coxeter's wonderful paper about quaternion rotations & reflections here:
www.math.utah.edu/~ptrapa/math-library/coxeter/Coxeter-Quaternions-and-reflections-AMM-1946.pdf
Coxeter shows that quaternions can even do all the rotations in 4D, but there you need a different quaternion for the pre- and post- multiplication in order to get enough degrees of freedom.
At 12:13 AM 6/26/2014, Robert Smith wrote:
Last one for the night. (It's midnight here.)
I got the above math to work numerically as well, but now I want to see if it's possible to boil that down into a single quaternion Q, as opposed to a function, so that for
Q # v = Q * v * conj(Q),
we have the simultaneously satisfied
Q # A = A' Q # B = B'.
I am introducing all of these symbols (#, ., *, implicit) just to be absolutely clear. I guess it's the programmer in me.
Do you think there's a nice expression for this? I'm betting on "no", but you've surprised me above.
Robert
On Wed, Jun 25, 2014 at 11:20 PM, Robert Smith <quad@symbo1ics.com> wrote:
Sorry for the noise. I also meant
Let R = (A' - A) x (B' - B) Let U = R x A Let U' = R x A'
Robert
On Wed, Jun 25, 2014 at 10:22 PM, Robert Smith <quad@symbo1ics.com> wrote:
Sorry, I guess what I said doesn't quite make sense in the last bit. When I said Q * A, I meant more rotate(Q, A) = A', where
rotate(Q, A) = Q * A * conj(Q).
Robert
On Wed, Jun 25, 2014 at 6:06 PM, Robert Smith <quad@symbo1ics.com> wrote:
Henry:
You've provided a function which maps a quaternion
v = r + xi + yj + zk
to a new quaternion
v' = r' + x'i + y'j + z'k.
If I followed along correctly, the expression can be rewritten as follows:
A -> A' B -> B'
Let S = (A' - A) x (B' - B) Let U = R x A Let U = R x A'
Then the mapping sought is
v -> (|U|^2 - U' * U) * v * (|U|^2 - U * U') --------------------------------------- 2(U . U)(U . (U + U'))
which is equal to
v -> (|U|^2 - U' * U) * v * (|U|^2 - U * U') --------------------------------------- 2 |U|^2 (|U|^2 + U . U')
where implicit multiplication is scalar-by-scalar, the dot is vector dot product (which in all these cases, if U is interpreted as a quaternion, then it has a realpart of 0), and the asterisk * is quaternion multiply associating left to right.
If the above interpretation is right, how can we recover a quaternion Q such that
Q * A = A' Q * B = B'
where a vector [x y z] is the quaternion xi + yj + ck? Just use v = 1?
I haven't verified this correct yet because I'm trying to ensure all the types check out. :)
Cheers,
Robert
On Wed, Jun 25, 2014 at 3:50 PM, Robert Smith <quad@symbo1ics.com> wrote: > Fantastic. I didn't expect it to be as nice as that, for certain! > > Robert > > On Fri, Jun 13, 2014 at 9:54 AM, Henry Baker <hbaker1@pipeline.com> wrote: >> I have confirmed via a number of numerical experiments that >> >> v->(|SxA|^2-(SxA')(SxA))v(|SxA|^2-(SxA)(SxA'))/Nq >> >> where S=(A'-A)x(B'-B) >> >> does indeed transform A->A' and B->B'. >> >> Here's a slightly simpler form for >> >> Nq = 2*|SxA|^2*(|SxA|^2+(SxA).(SxA')) >> = 2*[(SxA).(SxA)]*[(SxA).(Sx(A+A'))] >> >> At 09:24 PM 6/10/2014, Henry Baker wrote: >>>Another solution based on vectors & quaternions w/o matrices. >>> >>>If A,B are unit vectors (|A|=|B|=1), then rotating them >>>will preserve length (|A'|=|B'|=1). >>> >>>We need to warm up a bit before solving the full problem; >>>we first consider mapping A->A' without worrying about B. >>> >>>The most straightforward rotation that maps A->A' is >>>the quaternion q=sqrt((-A')A). Here (-A')A is >>>the quaternion product of the quaternion inverse of A' >>>and A. (For unit vector quaternions like A', (A')^-1=-A'.) >>>We use the standard mapping of vectors into >>>quaternions: A -> 0+A; i.e., the scalar part of (0+A) >>>is 0, and the vector part of (0+A) is A itself. >>> >>>The mapping >>> >>>v -> q v (q)^-1 = (q v q*)/(qq*) >>> >>>will transform A->A'. q* is conjugate(q). >>> >>>(The second form is preferred, since we don't have to >>>normalize q, and we don't have to compute absolute values.) >>> >>>The "sqrt()" above is the _quaternion_ square root, which >>>can be computed: >>> >>>sqrt(q) = (|q|+q)/sqrt(2*(|q|+S(q))) >>> >>>(Note that the denominator is a real scalar quantity.) >>> >>>Here S(q) is the "scalar part" of q, or >>> >>>S(q) = (q+q*)/2 >>> >>>But since we don't need a _normalized_ q, we can use a >>>cheap sqrt: >>> >>>elcheaposqrt(q) = |q|+q (!) >>> >>>Proof: >>> >>>(|q|+q)^2 = |q|^2+2|q|q+q^2 >>> = |q|^2+2|q|q+(Sq+Vq)^2 >>> = |q|^2+2|q|q+Sq^2+2SqVq+Vq^2 >>> = |q|^2+2|q|q+Sq^2+2SqVq-|Vq|^2 >>> = |q|^2+2|q|q+Sq^2+2SqVq-|Vq|^2-Sq^2+Sq^2 >>> = |q|^2+2|q|q+2Sq^2+2SqVq-|q|^2 >>> = 2|q|q+2Sq^2+2SqVq >>> = 2|q|q+2Sq(Sq+Vq) >>> = 2|q|q+2Sq q >>> = 2(|q|+Sq) q >>> >>>QED >>> >>>In particular, the function that takes A->A' is >>> >>>v -> (1-A'A)v(1-AA')/2/(1+A.A') >>> >>>(Here A.A' is the dot product of A & A'.) >>> >>>We now attack your actual question: you want a >>>rotation that _simultaneously_ transports >>>A->A' and B->B'. >>> >>>As Dan pointed out, the axis of rotation that we >>>want is (A'-A)x(B'-B), but we can't use (A'-A).(B'-B) >>>as the scalar part of a quaternion, because that >>>would rotate the projection of A into the projection >>>of B (or equivalently, the projection of A' into the >>>projection of B') -- not what we want. >>> >>>Let S=(A'-A)x(B'-B) be our axiS. If we are careful, >>>we won't have to normalize S into a unit vector. We >>>also no longer need B,B' since we only needed B,B' >>>to establish the axis S. >>> >>>If we compute SxA and SxA', we note that they are in >>>a plane orthogonal to S, and during the rotation that >>>maps A->A' and B->B', we will also map SxA into SxA'. >>>So the axis of this rotation is S itself. Also, >>>obviously, |SxA|=|SxA'|. >>> >>>So, we can now utilize the method of our warmup >>>exercise to map SxA into SxA' via the mapping: >>> >>>v->(|SxA|^2-(SxA')(SxA))v(|SxA|^2-(SxA)(SxA'))/Nq >>> >>>and Nq is the "absolute value squared" of >>>(|SxA|^2-(SxA')(SxA)) or of (|SxA|^2-(SxA)(SxA')). >>> >>>This is the mapping that we seek that will >>>simultaneously map A->A' and B->B'. >>> >>>Thus, although we talked about square roots and >>>absolute values in our development, we find that >>>the final result requires no square roots or >>>absolute values. >>> >>>I'd be surprised if there is a simpler expression >>>for computing this quaternion. >>> >>>At 04:03 PM 6/9/2014, Henry Baker wrote: >>>>We can take our matrix product below, which is >>>>an orthogonal matrix. Let's call it M. >>>> >>>>Then we can easily recover the quaternion from >>>>the orthogonal matrix M using the inverse _Cayley >>>>Transform_: >>>> >>>>S = (I-M)(I+M)^-1 >>>> >>>>The matrix I+M is invertible, because we have _rigid_ >>>>rotations (i.e., no reflections, hence no eigenvalues >>>>of -1). >>>> >>>>S is a _skew-symmetric_ matrix: >>>> >>>>[ 0 z -y] >>>>[-z 0 x] >>>>[ y -x 0] = S >>>> >>>>Then our (non-normalized) quaternion is [1,x,y,z]. >>>> >>>>http://en.wikipedia.org/wiki/Cayley_transform >>>> >>>>At 02:06 PM 6/7/2014, Henry Baker wrote: >>>>>We can easily do it with vectors & matrices: >>>>> >>>>>We are given A,B,A',B', where |A|=|B|=|A'|=|B'|. >>>>>R(A)=A', R(B)=B'. >>>>> >>>>>Because the rotation R is rigid, we know that >>>>> >>>>>R(AxB) = A'xB' = R(A)xR(B) >>>>> >>>>>Assuming that B is not a multiple of A, we can >>>>>form the 3x3 matrices R=matrix(A,B,AxB) and >>>>>R'=matrix(A',B',A'xB'). >>>>> >>>>>Our answer is then the matrix product >>>>> >>>>>R^-1 R' = >>>>> >>>>>matrix(A,B,AxB)^-1 matrix(A',B',A'xB') >>>>> >>>>>You are now free to convert this _orthogonal_ >>>>>matrix into a quaternion at your leisure. >>>>> >>>>>(Although AxB, A'xB' are not unit vectors, >>>>>hence R,R' are not orthogonal matrices, >>>>>nevertheless, R^-1 R' _is_ orthogonal.) >>>>> >>>>>(Contrary to Knuth, I did actually try this >>>>>on some numbers & it worked!) >>>>> >>>>>At 04:25 AM 6/7/2014, Dan Asimov wrote: >>>>>>In almost all cases, V := (A-A')x(B-B') will be a vector in the direction of the axis of rotation. Knowing V makes it easy to project say A and A' onto the perpendicular plane to C to determine the angle T. >>>>>> >>>>>>The remaining cases should be easy to exclude or deal with. >>>>>> >>>>>>--Dan >>>>>> >>>>>>On Jun 7, 2014, at 1:32 AM, Robert Smith <quad@symbo1ics.com> wrote: >>>>>> >>>>>>> Let A and B be unit vectors in R^3. Suppose they are rotated about >>>>>>> some vector V by an angle T, resulting in A' and B' respectively. What >>>>>>> are V and T? >>>>>>> >>>>>>> I set up a quadratic system using quaternions and got a result that >>>>>>> was 3 million terms large. Am I missing something?
Most formally speaking, for a quaternion Q and a vector v in R^3, and for * meaning quaternion multiplication, Q # v = Q * (v . <i, j, k>) * Q^{-1} where i, j, and k are the pure unit quaternions. Originally I said conj(Q) instead of Q^{-1} which is only correct some of the time. Robert On Thu, Jun 26, 2014 at 11:23 AM, Henry Baker <hbaker1@pipeline.com> wrote:
I didn't understand your "#" operation. Could you please try again?
At 10:46 AM 6/26/2014, Robert Smith wrote:
Thanks for the info.
I agree that a quaternion multiply is not sufficient, unless you define quaternion-vector "multiplication" by the operator # above.
Robert
On Thu, Jun 26, 2014 at 6:07 AM, Henry Baker <hbaker1@pipeline.com> wrote:
I did provide you with a single quaternion s.t.
Q A (1/Q) = A' Q B (1/Q) = B'
However, one normally (pun intended!) does the computation like this:
(Q A conjugate(Q))/N(Q) = A' (Q B conjugate(Q))/N(Q) = B'
where N(Q) = Q conjugate(Q)
because it is cheaper/easier/more precise to normalize at the end.
If you compute the _unit_ quaternion Q/sqrt(N(Q)), then in this case,
(1/Q) = conjugate(Q)
There is no way to do the rotation with only one quaternion multiply, because each quaternion multiply only computes a _reflection_, and you need 2 reflections to make a rotation.
See Coxeter's wonderful paper about quaternion rotations & reflections here:
www.math.utah.edu/~ptrapa/math-library/coxeter/Coxeter-Quaternions-and-reflections-AMM-1946.pdf
Coxeter shows that quaternions can even do all the rotations in 4D, but there you need a different quaternion for the pre- and post- multiplication in order to get enough degrees of freedom.
At 12:13 AM 6/26/2014, Robert Smith wrote:
Last one for the night. (It's midnight here.)
I got the above math to work numerically as well, but now I want to see if it's possible to boil that down into a single quaternion Q, as opposed to a function, so that for
Q # v = Q * v * conj(Q),
we have the simultaneously satisfied
Q # A = A' Q # B = B'.
I am introducing all of these symbols (#, ., *, implicit) just to be absolutely clear. I guess it's the programmer in me.
Do you think there's a nice expression for this? I'm betting on "no", but you've surprised me above.
Robert
On Wed, Jun 25, 2014 at 11:20 PM, Robert Smith <quad@symbo1ics.com> wrote:
Sorry for the noise. I also meant
Let R = (A' - A) x (B' - B) Let U = R x A Let U' = R x A'
Robert
On Wed, Jun 25, 2014 at 10:22 PM, Robert Smith <quad@symbo1ics.com> wrote:
Sorry, I guess what I said doesn't quite make sense in the last bit. When I said Q * A, I meant more rotate(Q, A) = A', where
rotate(Q, A) = Q * A * conj(Q).
Robert
On Wed, Jun 25, 2014 at 6:06 PM, Robert Smith <quad@symbo1ics.com> wrote: > Henry: > > You've provided a function which maps a quaternion > > v = r + xi + yj + zk > > to a new quaternion > > v' = r' + x'i + y'j + z'k. > > If I followed along correctly, the expression can be rewritten as follows: > > A -> A' > B -> B' > > Let S = (A' - A) x (B' - B) > Let U = R x A > Let U = R x A' > > Then the mapping sought is > > v -> > (|U|^2 - U' * U) * v * (|U|^2 - U * U') > --------------------------------------- > 2(U . U)(U . (U + U')) > > which is equal to > > v -> > (|U|^2 - U' * U) * v * (|U|^2 - U * U') > --------------------------------------- > 2 |U|^2 (|U|^2 + U . U') > > where implicit multiplication is scalar-by-scalar, the dot is vector > dot product (which in all these cases, if U is interpreted as a > quaternion, then it has a realpart of 0), and the asterisk * is > quaternion multiply associating left to right. > > If the above interpretation is right, how can we recover a quaternion > Q such that > > Q * A = A' > Q * B = B' > > where a vector [x y z] is the quaternion xi + yj + ck? Just use v = 1? > > I haven't verified this correct yet because I'm trying to ensure all > the types check out. :) > > Cheers, > > Robert > > On Wed, Jun 25, 2014 at 3:50 PM, Robert Smith <quad@symbo1ics.com> wrote: >> Fantastic. I didn't expect it to be as nice as that, for certain! >> >> Robert >> >> On Fri, Jun 13, 2014 at 9:54 AM, Henry Baker <hbaker1@pipeline.com> wrote: >>> I have confirmed via a number of numerical experiments that >>> >>> v->(|SxA|^2-(SxA')(SxA))v(|SxA|^2-(SxA)(SxA'))/Nq >>> >>> where S=(A'-A)x(B'-B) >>> >>> does indeed transform A->A' and B->B'. >>> >>> Here's a slightly simpler form for >>> >>> Nq = 2*|SxA|^2*(|SxA|^2+(SxA).(SxA')) >>> = 2*[(SxA).(SxA)]*[(SxA).(Sx(A+A'))] >>> >>> At 09:24 PM 6/10/2014, Henry Baker wrote: >>>>Another solution based on vectors & quaternions w/o matrices. >>>> >>>>If A,B are unit vectors (|A|=|B|=1), then rotating them >>>>will preserve length (|A'|=|B'|=1). >>>> >>>>We need to warm up a bit before solving the full problem; >>>>we first consider mapping A->A' without worrying about B. >>>> >>>>The most straightforward rotation that maps A->A' is >>>>the quaternion q=sqrt((-A')A). Here (-A')A is >>>>the quaternion product of the quaternion inverse of A' >>>>and A. (For unit vector quaternions like A', (A')^-1=-A'.) >>>>We use the standard mapping of vectors into >>>>quaternions: A -> 0+A; i.e., the scalar part of (0+A) >>>>is 0, and the vector part of (0+A) is A itself. >>>> >>>>The mapping >>>> >>>>v -> q v (q)^-1 = (q v q*)/(qq*) >>>> >>>>will transform A->A'. q* is conjugate(q). >>>> >>>>(The second form is preferred, since we don't have to >>>>normalize q, and we don't have to compute absolute values.) >>>> >>>>The "sqrt()" above is the _quaternion_ square root, which >>>>can be computed: >>>> >>>>sqrt(q) = (|q|+q)/sqrt(2*(|q|+S(q))) >>>> >>>>(Note that the denominator is a real scalar quantity.) >>>> >>>>Here S(q) is the "scalar part" of q, or >>>> >>>>S(q) = (q+q*)/2 >>>> >>>>But since we don't need a _normalized_ q, we can use a >>>>cheap sqrt: >>>> >>>>elcheaposqrt(q) = |q|+q (!) >>>> >>>>Proof: >>>> >>>>(|q|+q)^2 = |q|^2+2|q|q+q^2 >>>> = |q|^2+2|q|q+(Sq+Vq)^2 >>>> = |q|^2+2|q|q+Sq^2+2SqVq+Vq^2 >>>> = |q|^2+2|q|q+Sq^2+2SqVq-|Vq|^2 >>>> = |q|^2+2|q|q+Sq^2+2SqVq-|Vq|^2-Sq^2+Sq^2 >>>> = |q|^2+2|q|q+2Sq^2+2SqVq-|q|^2 >>>> = 2|q|q+2Sq^2+2SqVq >>>> = 2|q|q+2Sq(Sq+Vq) >>>> = 2|q|q+2Sq q >>>> = 2(|q|+Sq) q >>>> >>>>QED >>>> >>>>In particular, the function that takes A->A' is >>>> >>>>v -> (1-A'A)v(1-AA')/2/(1+A.A') >>>> >>>>(Here A.A' is the dot product of A & A'.) >>>> >>>>We now attack your actual question: you want a >>>>rotation that _simultaneously_ transports >>>>A->A' and B->B'. >>>> >>>>As Dan pointed out, the axis of rotation that we >>>>want is (A'-A)x(B'-B), but we can't use (A'-A).(B'-B) >>>>as the scalar part of a quaternion, because that >>>>would rotate the projection of A into the projection >>>>of B (or equivalently, the projection of A' into the >>>>projection of B') -- not what we want. >>>> >>>>Let S=(A'-A)x(B'-B) be our axiS. If we are careful, >>>>we won't have to normalize S into a unit vector. We >>>>also no longer need B,B' since we only needed B,B' >>>>to establish the axis S. >>>> >>>>If we compute SxA and SxA', we note that they are in >>>>a plane orthogonal to S, and during the rotation that >>>>maps A->A' and B->B', we will also map SxA into SxA'. >>>>So the axis of this rotation is S itself. Also, >>>>obviously, |SxA|=|SxA'|. >>>> >>>>So, we can now utilize the method of our warmup >>>>exercise to map SxA into SxA' via the mapping: >>>> >>>>v->(|SxA|^2-(SxA')(SxA))v(|SxA|^2-(SxA)(SxA'))/Nq >>>> >>>>and Nq is the "absolute value squared" of >>>>(|SxA|^2-(SxA')(SxA)) or of (|SxA|^2-(SxA)(SxA')). >>>> >>>>This is the mapping that we seek that will >>>>simultaneously map A->A' and B->B'. >>>> >>>>Thus, although we talked about square roots and >>>>absolute values in our development, we find that >>>>the final result requires no square roots or >>>>absolute values. >>>> >>>>I'd be surprised if there is a simpler expression >>>>for computing this quaternion. >>>> >>>>At 04:03 PM 6/9/2014, Henry Baker wrote: >>>>>We can take our matrix product below, which is >>>>>an orthogonal matrix. Let's call it M. >>>>> >>>>>Then we can easily recover the quaternion from >>>>>the orthogonal matrix M using the inverse _Cayley >>>>>Transform_: >>>>> >>>>>S = (I-M)(I+M)^-1 >>>>> >>>>>The matrix I+M is invertible, because we have _rigid_ >>>>>rotations (i.e., no reflections, hence no eigenvalues >>>>>of -1). >>>>> >>>>>S is a _skew-symmetric_ matrix: >>>>> >>>>>[ 0 z -y] >>>>>[-z 0 x] >>>>>[ y -x 0] = S >>>>> >>>>>Then our (non-normalized) quaternion is [1,x,y,z]. >>>>> >>>>>http://en.wikipedia.org/wiki/Cayley_transform >>>>> >>>>>At 02:06 PM 6/7/2014, Henry Baker wrote: >>>>>>We can easily do it with vectors & matrices: >>>>>> >>>>>>We are given A,B,A',B', where |A|=|B|=|A'|=|B'|. >>>>>>R(A)=A', R(B)=B'. >>>>>> >>>>>>Because the rotation R is rigid, we know that >>>>>> >>>>>>R(AxB) = A'xB' = R(A)xR(B) >>>>>> >>>>>>Assuming that B is not a multiple of A, we can >>>>>>form the 3x3 matrices R=matrix(A,B,AxB) and >>>>>>R'=matrix(A',B',A'xB'). >>>>>> >>>>>>Our answer is then the matrix product >>>>>> >>>>>>R^-1 R' = >>>>>> >>>>>>matrix(A,B,AxB)^-1 matrix(A',B',A'xB') >>>>>> >>>>>>You are now free to convert this _orthogonal_ >>>>>>matrix into a quaternion at your leisure. >>>>>> >>>>>>(Although AxB, A'xB' are not unit vectors, >>>>>>hence R,R' are not orthogonal matrices, >>>>>>nevertheless, R^-1 R' _is_ orthogonal.) >>>>>> >>>>>>(Contrary to Knuth, I did actually try this >>>>>>on some numbers & it worked!) >>>>>> >>>>>>At 04:25 AM 6/7/2014, Dan Asimov wrote: >>>>>>>In almost all cases, V := (A-A')x(B-B') will be a vector in the direction of the axis of rotation. Knowing V makes it easy to project say A and A' onto the perpendicular plane to C to determine the angle T. >>>>>>> >>>>>>>The remaining cases should be easy to exclude or deal with. >>>>>>> >>>>>>>--Dan >>>>>>> >>>>>>>On Jun 7, 2014, at 1:32 AM, Robert Smith <quad@symbo1ics.com> wrote: >>>>>>> >>>>>>>> Let A and B be unit vectors in R^3. Suppose they are rotated about >>>>>>>> some vector V by an angle T, resulting in A' and B' respectively. What >>>>>>>> are V and T? >>>>>>>> >>>>>>>> I set up a quadratic system using quaternions and got a result that >>>>>>>> was 3 million terms large. Am I missing something?
participants (2)
-
Henry Baker -
Robert Smith