# Finding the rotation matrix in n-dimensions

Suppose that we know two real vectors with n components, which are linked by some arbitrary transformation/scaling/rotation/shearing...

Now, I think that it is possible to know which is the scaling matrix and the rotation matrix. For example, the scaling matrix would be a diagonal matrix with n entries representing the n scaling factors.

On the other side, I can normalize the two vectors and then compute the rotation matrix between the two, isn't it?

1) How can I retrieve the rotation matrix? (if it is possible, obviously)

2) I need such matrix in order to use it in a computational mechanics context. Then, how would you find it in an operative way?

I have been reading somewhere that in multiple dimensions we can apply rotations that are combination of rotations in n-2 hyperplanes. Is this a convenient way to compute such matrix?

Thank you very much!!

EDIT: Let us put this way, in order to understand better what I should do...

Vector $x= [ 2, 4, 5, 3, 6 ]^T$ and vector $y= [ 6, 2, 0, 1, 7 ]^T$. I would like to find the rotation matrix that aligns vector x to vector y. First of all, I understood that one need to find the base of the orthogonal complement to the two vectors, i.e. find the hyperplane containing such two vectors. Consequently, one first compute the null space:

• Reducing rows: $$A= \left[ \begin{array}{ccccc} 1 & 0 & -1/2 & -1/10 & 4/5 \\ 0 & 1 & 3/2 & 4/5 & 11/10 \\ \end{array} \right]$$

• Vectors of orthogonal space (linked to $x_3, x_4, x_5$) $$v_1^{\perp}= \left[ \begin{array}{c} 1/2 \\ -3/2 \\ 1 \\ 0 \\ 0 \\ \end{array} \right], v_2^{\perp}= \left[ \begin{array}{c} 1/10 \\ -4/5 \\ 0 \\ 1 \\ 0 \\ \end{array} \right], v_3^{\perp}= \left[ \begin{array}{c} -4/5 \\ -11/10 \\ 0 \\ 0 \\ 1 \\ \end{array} \right]$$

• Then one can do Gram-Schmidt over the two groups of vectors, the one representing the plane containing the two initial vectors and the one representing the orthogonal subspace. The resulting matrix presents in the column vectors the basis of the space:

$$E= \left[ \begin{array}{ccccc} 0.7255 & -0.0117 & 0.2673 & -0.0716 & -0.6301 \\ 0 & 0.4429 & -0.8018 & -0.2409 & -0.3209 \\ -0.3627 & 0.6701 & 0.5345 & -0.3255 & -0.1663 \\ -0.0725 & 0.3555 & 0 & 0.9115 & -0.1937 \\ 0.5804 & 0.4778 & 0 & 0 & 0.6594 \\ \end{array} \right]$$

The last matrix is correct and represent a base because obviously the scalar product between two columns i,j is equal to the Kronecker delta $\delta_{i,j}$.

Now, I have information about almost everything, I can compute the angle between the two vectors saying that

$$\cos(\theta)=\frac{x \cdot y}{\left|| x \right|| \: \left|| y \right||}$$

but now how do I construct the rotation matrix?

One way to do this is to find two orthonormal vectors in the plane generated by your two vectors, and then extend it to an orthonormal basis of $\mathbb R^n$. Then with respect to this basis, consider the rotation by angle $\theta$ in the plan generated by the first two vectors, and the identity on the space generated by the rest of the orthonormal basis. Use Gram-Schmidt to find the orthonormal basis.

As you said in a previous comment, you cannot rotate around an axis except in 3D. Rather you need to rotate about an $n-2$-dimensional subspace.

So suppose you want to rotate $x$ to $y$, and you happen to know they are the same norm. Let $u = x/|x|$, and $v = (y-(u.y)u)/|y-(u.y)u|$. Then $P=uu^T + vv^T$ is a projection onto the space generated by $x$ and $y$, and $Q=I-uu^T-vv^T$ is the projection onto the $n-2$-dimensional complemented subspace. So the "rotation" part just has to take place on the range of $P$. That is, $z \mapsto (z.u,z.v)$ is a isomorphic isometry of the range of $P$ to $\mathbb R^2$. Do the rotation on $\mathbb R^2$. Then map this back to $\mathbb R^n$ by $[a,b]^T \mapsto au + bv$.

So the whole rotation is $$I - uu^T -vv^T + [u\ v]R_\theta[u\ v]^T$$ where $$R_\theta = \begin{bmatrix}\cos(\theta) & -\sin(\theta) \\ \sin(\theta & \cos(\theta) \end{bmatrix}$$ where $\cos(\theta) = x\cdot y/(|x||y|)$.

(By the way, if a reflection matrix is OK, then consider Householder transformations instead.)

Here is some MATLAB code (but tested using octave so maybe it has syntax errors in matlab).

x=[2,4,5,3,6]';

y=[6,2,0,1,7]';

u=x/norm(x);

v=y-u'*y*u;

v=v/norm(v);

cost=x'*y/norm(x)/norm(y);

sint=sqrt(1-cost^2);

R = eye(length(x))-u*u'-v*v' + [u v]* [cost -sint;sint cost] *[u v]';

fprintf('testing\n');

should_be_identity = R*R'

should_be_y = R*x

• I lost your reasoning when you introduce the isomorphic isometry of the range of P to $R^2$. I understand that any vector multiplied by P is projected to the plane given by the two vectors, so since the rotation will happen in this plane, the only components that matter lays in such plane, the others (for which Q is the "corresponding" part) will stay the same. But then, what do you mean by "Do the rotation on $R^2$"? Are you saying that I could use a classical "sine/cosine" matrix to change components on the plane? then, how would you map back in the n-dimensional space? – L. B. Dec 9 '13 at 14:18
• Yes, do the rotation on $\mathbb R^$, just as you say. And the map back to $\mathbb R^n$ is given by the last line of the second to last paragraph. – Stephen Montgomery-Smith Dec 9 '13 at 14:36
• I made the rotation more explicit. – Stephen Montgomery-Smith Dec 9 '13 at 14:43
• And I gave you MATLAB code. – Stephen Montgomery-Smith Dec 9 '13 at 15:02
• Ok, now it is way clearer, I really appreciate it. Some doubts: if the whole rotation is happening in $\mathbb{R}^2$ why don't we use only the $[u v]R_{\theta}[u v]^T$ part? As you said the rest of the values on the other axes will be unchanged, isn't it? You suggested also Householder transformations: those are reflections by an hyperplane.. how would you use them to do a rotation? the only thing that comes to my mind now would be finding the "bisectrix" hyperplane and to reflect with that... – L. B. Dec 9 '13 at 17:50

Are you familiar with quaternions? Clifford algebra offers a good $N$-dimensional generalization to quaternions.

Let me take your example. You have two vectors $x$ and $y$:

\begin{align*} x &= 2e_1 + 4e_2 + 5e_3 + 3e_4 + 6 e_5 \\ y &= 6 e_1 + 2 e_2 + 0 e_3 + 1e_4 + 7 e_5\end{align*}

With quaternions, one would compute a cross product to find the rotation axis. Using clifford algebra, one directly computes the plane of rotation using the wedge product of vectors, as one would do in exterior algebra. (A general term for such geometric primitives is blade. Vectors are just 1-dimensional blades; here, we will compute a 2-dimensional blade.)

The wedge product is

\begin{align*}x \wedge y = (4 - 24) e_{12} &+ (0 -30) e_{13} + (2 - 18) e_{14} + (14 - 36) e_{15} \\ +(0 - 10) e_{23} &+ (4 - 6) e_{24} + (28 - 12) e_{25} \\ + (5 - 0) e_{34} &+ (35 - 0) e_{35} \\ + (21 - 6) e_{45} & \end{align*}

You should appreciate how direct this is compared to finding the null space of a matrix; when you're finding the null space, you're finding the orthogonal complement to this blade; the wedge product allows us to proceed directly instead.

Anyway, to proceed as one would with quaternions, we must made the blade unit. This is similar to the vector case--simply sum the squares of the coefficients and divide by the square root. I will not carry out this arithmetic here.

Once one has the normalized blade $\hat B$, one can construct the versor $R = \exp(-\hat B \theta/2)$. The minus sign is conventional in clifford algebra, to correspond to counterclockwise rotation. This is different from how quaternions usually do.

The versor $R$ is then used to construct a rotation map $\underline R$ that takes in any vector $a$ and produces its rotated counterpart:

$$\underline R(a) = Ra R^{\dagger} = e^{-\hat B \theta/2} a e^{\hat B \theta/2}$$

The quantity $\hat B$ obeys $\hat B^2 = -1$, and so this generates the trigonometric functions through power series, just as in quaternions.

The products here are the geometric product of clifford algebra. This can be computed in terms of basis vectors like so: for any $i, j = 1, 2, \ldots, N$, set $e_i e_i = +1$, and when $i \neq j$, $e_i e_j = - e_i e_j$.

For instance, say that $R$ takes the form $\alpha - \beta e_{12}$ and $a = a^1 e_1 + a^2 e_2$. The computation would then look like

\begin{align*}\underline R(a) &= (\alpha - \beta e_1 e_2)(a^1 e_1 + a^2 e_2)(\alpha + \beta e_1 e_2) \\ &= (\alpha a - \beta a^1 e_1 e_2 e_1 - \beta a^2 e_1 e_2 e_2)(\alpha + \beta e_1 e_2) \\ &= (\alpha a + \beta a^1 e_2 - \beta a^2 e_1) (\alpha + \beta e_1 e_2) \\ &= \ldots \end{align*}

(I won't finish this calculation; I just hoped to illustrate how geometric products are computed.)

This is the most direct rotation that connects the two directions. There could be other, subsequent rotations that rotate other vectors orthogonal to these while leaving these two intact, so there are several rotation maps that map these vectors as you specified.

Once you have the rotation map $\underline R$, you can extract the matrix components by plugging in basis vectors and taking dot products. $\underline R_{12} = e_1 \cdot \underline R(e_2)$, for instance. But if you actually go to all the trouble to write something that can compute this rotation in terms of clifford algebra, it makes one wonder why you would even bother with a matrix. There is simply no need.

A special use of this clifford algebra formulation is that, using the geometric product, it can be used to rotate higher-dimensional blade, not just vectors. Indeed, for a general blade $K$, the formula is just

$$\underline R(K) = R K R^\dagger$$

This is a special property of rotations. If you have a direct expression for the blade $K$, you need not bother decomposing it into vectors or rotate those vectors individually to recover the image of $K$. This is handy.

• I really appreciate your point of view of the problem from an general algebraic approach as the Clifford Algebra.. Intuitively, I understand the steps, however I never studied those kind of algebras, but it seems very interesting for its neat and powerful approach... could you give me a possible reference, an introductory book or similar for such topic? – L. B. Dec 9 '13 at 23:48
• I knew about using quaternions to represent 3D rotations. But I had never come across using Clifford algebras to represent $N$-D rotations. The nice thing is that the dimension of the "rank 2" clifford subalgebra is $N(N-1)/2$, which I know to be correct. It is rather like (perhaps identical) to alternating forms on $\mathbb R^N$. I think this is very clever. – Stephen Montgomery-Smith Dec 9 '13 at 23:57
• So the unit quaternions provide a double cover of $SO(3)$. What happens in higher dimensions? Is it still a double cover? Oh, but maybe there is no group structure in the higher dimensional Clifford algebras. – Stephen Montgomery-Smith Dec 9 '13 at 23:58
• @L.B. I've mostly read books that lean toward the physics side of things as far as clifford (or geometric, as it's sometimes called) algebras go. Alan MacDonald has some very introductory books on the subject. Something more advanced would be Hestenes and Sobczyk, or perhaps Gull, Doran, and Lasenby. – Muphrid Dec 10 '13 at 0:15
• @StephenMontgomery-Smith Yes! The 2-blades (or bivectors) are at the heart of alternating forms. As you might've noticed, any rotation can be completely described by its versor (or "rotor") $R$, as $R^\dagger$ contains no other information. The product of two rotors is another rotor, and $-R$ represents the same rotation as $R$, so there is a double covering in all dimensions. – Muphrid Dec 10 '13 at 0:17

To find the rotation matrix:

1. Find the plane $P$ containing the two vectors and the origin (always possible).

2. Compute the angle between the two vectors. Call this angle $\alpha.$

3. Rotate by $\alpha$ around the axis which is the normal vector to $P.$

EDIT For more detail see: http://mathforum.org/library/drmath/view/74532.html

• Such procedure is very abstract. However, if you think it a little better, it is neither so trivial nor operative. Find the plane P means "find the span of the two vectors". Compute the angle means "find the cosine with scalar product". Finally rotate by alpha around the axis does not help since the formulas that I found till now gives such rotation matrix for 2D and 3D, while here I would like a solution for a n-dimensional problem... – L. B. Dec 8 '13 at 22:00
• It is quite operative. However, if you want more help, look at the edit. – Igor Rivin Dec 8 '13 at 22:17
• Oh, and thanks for the downvote! – Igor Rivin Dec 9 '13 at 0:43
• Sorry for the downvote, but it wasn't me. I saw your edit, but still I do not have a clear idea of what procedure would you use to recover the rotation matrix. I added an example in the first post in order to make things more understandable. – L. B. Dec 9 '13 at 14:14
• If in dimensions higher than 3 you cannot rotate around an axis, how could you rotate around the normal vector of the plane? What you describe is a simple axis-angle to rotation matrix conversion, which only works in 3D as I know. – plasmacel Aug 16 '16 at 16:13

This is not technically an answer to your question, but it may be helpful to some and is a bit longer than a comment. If you just want a way to rotate vectors, computing the rotation matrix is not efficient when the dimension of the vector space is large. The following is much faster; all operations are linear complexity in the dimension of your vector space.

Suppose v and w are your (unit) vectors. To compute the rotation that takes v to w (note: (,) denotes inner product):

1) If v and w are not orthogonal, let w' = w - (w,v) v / |w - (w,v) v|. Then w' is a unit vector in the v,w plane perpendicular to w. Otherwise, let w'=w.

2) Let x be the vector to be rotated. Let x_v = (x,v) and let x_w' = (x,w'). Then x projects orthogonally onto x_v v + x_w' w' in the vw plane. Apply the desired rotation in this two element basis to obtain a vector a v + b w'.

3) The solution is then the "deprojection" x - (x_v v + x_w' w') + a v + b w'.

Here's a super easy way to find a rotation matrix $R$ that rotates $u$ to align with $v$ in $N$ dimensions. One feature of this approach is that $R$ does not rotate vectors outside the subspace of $u$ and $v$.

Another way to generate rotation matrix M, which rotates given n-dimensional vector X to the direction of given n-dimensional vectot Y is to use NRMG algorithm, that consist the following operations:

1) Obtaining rotation matrix Mx, which rotates given vector X to the direction of axis $x_1$

2) Obtaining rotation matrix My, which rotates given vector Y to the direction of axis $x_1$

3) Obtaining rotation matrix M, which rotates given vector X to the direction of given vector Y as multiplication of matrix Mx and inverse matrix of My $$M=M_y^{-1}M_x=M_y^TM_x$$

Rotation of n-dimensional vector to the direction of one of coordinate axes (i.e. to axis $x_1$) can be performed by subsequent multiplications by $n-1$ Givens matrises $G(n-k,n-k+1,\theta_{n-k}), k=1,...,n-1$ , therefore corresponding rotation matrix can be obtained as follows $$M_x =\prod_{k=1}^{n-1} G(n-k,n-k+1,\theta_{n-k})$$ Coefficients $sin(\theta_{n-k})$ and $cos(\theta_{n-k})$ of Givens matrises have to be calculated using formulas: $$\begin{cases} \sin(\theta_{n-k})=-\frac{x_{n-k+1}}{\sqrt{x_{n-k}^2+x_{n-k+1}^2}}, \cos(\theta_{n-k}) = \frac{x_{n-k}}{\sqrt{x_{n-k}^2+x_{n-k+1}^2}}, \quad if\quad x_{n-k}^2+x_{n-k+1}^2>0 \\ \sin(\theta_{n-k})=0, \cos(\theta_{n-k})=1, \quad if\quad x_{n-k}^2+x_{n-k+1}^2=0 \end{cases}$$ where $x_{n-k}$ and $x_{n-k+1}$ are components of given vector. More detailed explanation of algorithm and it's Matlab implementation you can find on N-dimensional Rotation Matrix Generation Algorithm