See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/323995682

N-dimensional Rotation Matrix Generation Algorithm

Article · March 2018

CITATIONS READS

4 175

1 author:

Ognyan Zhelezov

3 PUBLICATIONS   4 CITATIONS   

SEE PROFILE

All content following this page was uploaded by Ognyan Zhelezov on 25 March 2018.

The user has requested enhancement of the downloaded file.

 N-dimensional Rotation Matrix Generation Algorithm
Ognyan Ivanov Zhelezov1

1 Dep. of Development, firm Data Consulting, Sofia, 1000, Bulgaria

Abstract This article presents a new algorithm for generation of N-dimensional rotation matrix M, which rotates given
N-dimensional vector X to the direction of given vector Y which has the same dimension. Algorithm, named
N-dimensional Rotation Matrix Generation Algorithm (NRMG) includes rotation of given vectors X and Y to the direction
of coordinate axis x1 using two-dimensional rotations. Matrix M is obtained as multiplication of matrix M X and inverse of
matrix MY, which rotates given vectors to the direction of axis x1. Also examined is the possibility to perform parallel cal-
culations of two-dimensional rotations.
Keywords Mathematics of computing , Mathematical analysis , Numerical analysis , Computations on matrices

M  P1.G(1,2, ).P (2)

1. Introduction
where P is matrix, which transforms initial basis B to basis
High-dimensional spaces frequently occur in mathemat- B2.
ics and the sciences, in example N-dimensional feature
space, which presents input signals of neural network or Another way to generate rotation matrix M is to use
collection of N-dimensional parameters for multidimen- Householder Reflection [4], [10], [17], [18]. If X and Y are
sional data analysis. Rotation is one of rigid transformations vectors with the same norm X  Y , there exists an or-
in geometrical space, which preserves length of vectors and thogonal symmetric matrix P such that Y  P. X where
can be presented using matrix operation like Y = M.X , P  I  W .W T and W   X  Y  / X  Y . Matrix P is ma-
where X and Y are input and output vector respectively, and trix of reflection (not a rotation) because det P = -1, (which
M is rotation matrix. This article proposes an gives the name to method) that’s why to obtain matrix of
N-dimensional rotation matrix generation algorithm for rotation M, for which det M = 1, have to be performed two
given input and output vector. subsequent reflections. Matrix of rotation can be obtained as
multiplication of two matrix of reflection P1 and P2 as M =
P1.P2 .
2. Definition of the Task
This article presents a new algorithm for generation of
Let's say that we have two N-dimensional vectors X and Y,
N-dimensional rotation matrix, which rotates given vector X
having the same dimension, X, Y  RN. We want to obtain a
to the direction of given vector Y. Algorithm, named
rotation matrix M that satisfies the equation.
N-dimensional Rotation Matrix Generation Algorithm
~ (NRMG) includes the following sequence of operations:
Y  M.X (1)
~ 1) Obtaining rotation matrix MX, which rotates given
where Y has the same norm as X and the same direction as 
~
  ~
Y, i.e. Y  X and cos Y, Y  1 .
vector X to the direction of axis x1 .
2) Obtaining rotation matrix MY, which rotates given
One of possibilities to generate rotation matrix M includes 
vector Y to the direction of axis x1 .
the following sequence of operations:
3) Obtaining rotation matrix M as multiplication of MX by
inverse matrix of MY given as M  My 1.Mx .
1) Obtain two orthogonal vectors in the plane, in which lies
the two given vectors using Gramm-Schmidt procedure
[2], [3], [16].
Below this three operations will be described subsequently
2) Enlarge this 2-dimensional basis to N- dimensional basis
B2, in which given vectors X, Y are transformed to X
and Y . 3. Rotation of Given Vector X to the
3) Create matrix of Givens rotation G(1, 2, θ) for RN, which Direction of Axis x1
makes rotation of vector X in x1 x2 -plane to the di-
rection of vector Y . Rotation of given vector X to the direction of one of co-

4) Obtaining rotation matrix M as ordinate axes (e.g. axis x1 ) can be performed by subsequent
multiplications by Givens matrices [1] as follow:
 1 is to set to zero coordinate xk 1 . It is easy to find that
X ( N 1)   Gk , k  1, k . X  M X . X  [rX ,0,0,0]T (3) equation xk 1  0 and equations (5) are satisfied simul-
k  N 1 taneously when sin( k ) and cos( k ) are calculated
Givens matrices G(k, k+1, θk), k = N-1, n-2, … 1 are de- using formulas:
fined as follows [1], [2]:
x k 1
G(k, k+1, θk) = sin(  k )   ,
x 2k  x 2k 1
1 0  0 0  0 0
xk (6)
0 1  0 0  0 0 cos( k )  if x 2k  x 2k 1  0
x 2k  x 2k 1
       
0 0  Ck  Sk  0 0 sin(  k )  0, cos( k )  1 if x 2k  x 2k 1  0

0 0  Sk Ck  0 0 Thus searched matrix MX can be calculated as multiplica-

        tion of Givens matrices
GN  1, N  2,  N 1 , GN  2, N  3,  N  2 ,, G1,2, 1 
0 0  0 0  1 0 as follows:
0 0  0 0  0 1
N 1
(4) M X   GN  k , N  k  1, N  k  (7)
k 1

where θk is the angle of rotation and coefficients

Ck=cos( θk) and Sk=sin(θk) appears at the intersections of where coefficients of rotation sin( k ) and cos( k )
k-th and k+1 -th rows and columns. are calculated using formulas (6). Calculation of angles of
Every one multiplication of vector by Givens matrix two-dimensional rotations k , k = N -1, N - 2, …, 1
G(k, k+1, θk) performs rotation of its projection in coordi- really is not needed. If x 2k  x 2k 1  0 then corresponding
nate plane (xk, xk+1), which changes values only of vector’s Givens matrix is equal to Identity matrix G(k, k+1, θk) = I –
coordinates xk , xk 1 to xk , xk 1 .This multiplication can no rotation.

be presented as multiplication of coordinates by sub matrix Schema for rotation of vector X to the direction of axis x1
A(k, k+1, θk) as follows: using two-dimensional rotations can be presented as fol-
lows:

xk xk cos k   sin  k  xk
 Ak , k  1,  k .  .
xk 1 xk 1 sin  k  cos k  xk 1
(5)

Taking in consideration that multiplication with Givens

matrix G(k, k+1, θk) changes only values of vector’s coor- 
dinates xk , xk 1 (to xk , xk 1 ), schema of multiplication Figure 2 Schema for rotation of vector X to the direction of axis x1
by Givens matrix G(k, k+1, θk) can be presented as operator
for two-dimensional rotation as follows:
Every one block for two-dimensional rotation (as it was
given on Fig. 3) presents rotation of two-dimensional vector
xk , xk 1 T to the direction of axis x k k=1,2,…N-1.
This two-dimensional rotation is the base operation of pro-
posed algorithm.

Figure 3 Block of base operation - rotation of two-dimensional vector

Figure 1 Schema of two-dimensional rotation, performed by Givens ma-
xk , xk 1 T to the direction of axis

xk
trix G(k, k+1, θk)

As it can be seen from (3) and Fig. 2, rotation of


The target of multiplication by Givens matrix G(k, k+1, θk) N-dimensional vector to the direction of axis x1 needs
 X 1  x1 , x2 , ... , xN / 2 ,0,,0T
execution of N-1 subsequent base operations. Thus if one
base operation has execution time Tb , then rotation of giv-

X 2  0,0,0, xN / 2 1 , xN / 2  2 ,, xN 
T
en N-dimensional vector X to the direction of axis x1
will take execution time (N-1).Tb . As it is described below,
this time can be reduced using parallel execution of base X  X1 X 2 (9)
operations.

Rotation of vector X to the direction of axis x1 by mul-
tiplication with matrix of rotation M will give resultant
4. Description of NRMG Algorithm vector X12(N/2) which is the sum of rotated to the direction of

NRMG algorithm for generation of N-dimensional rota- axis x1 vectors X1 and X2 as follows:
tion matrix M, which rotates given vector X to the direction
of given vector Y consist the following operations:
1) Obtaining rotation matrix MX, which rotates given X12 ( N / 2)  M.X  M.X1  X 2  X1( N / 2)  X 2( N / 2)

vector X to the direction of axis x1 N 1
2) Obtaining rotation matrix MY, which rotates given

  GN  k, N  k  1,  N  k . X1
vector Y to the direction of axis x1 . k 1
3) Obtaining rotation matrix M, which rotates given vec- N 1
tor X to the direction of given vector Y as multiplica-   GN  k, N  k  1,  N  k . X 2
k 1
tion of matrix MX and inverse matrix of MY , as fol-
lows: (10)
For zero coordinates of X1 and X2 corresponding matrices
of base operations G(k, k+1, θk) are equal to identity matrix I
M  M Y1.M X (8) (i.e. no rotation). This shows that coordinate planes, in which
are rotated projections of X1 are different from those, in
which are rotated X2. As a resultant vector we get:
It is important to note that NRMG algorithm do not need
vectors X and Y to have the same norm (as it is needed for
Householder Reflection i.e.). Multiplication of given vector N / 2 1
X by matrix of rotation M will give resultant vector Y ,
~ X12 ( N / 2)   G N / 2  k, N / 2  k  1,  N / 2  k . X1
k 1
which will have norm of vector X, but direction of vector Y. (11)
N / 2 1
Time complexity of proposed algorithm depends of algo-   G N  k, N  k  1,  N  k . X 2
rithm, used for calculation of coefficients sin(k ) and k 1
cos(k ) of two-dimensional base operations. This calcu-
lation depends of practical realization and is not a subject of where in second sum symbol for angle or rotation φ has been
this article. replaced with θ for convenience.
NRMG algorithm can be used for different reversible Vector X12(N/2) has only two non-zero coordinates - x1 and
calculations. An example that uses NRMG algorithm for xN/2+1 , so it lie in coordinate plane ( x1, x N / 2 1 ) . Vector

transformation of images appears below (6). X(N/2) to the direction of axis x1 can be obtained by only one
two-dimensional rotation of vector X12(N/2) in this coordinate
plane as follows:
5. Increasing Performance Using Paral-
lel Execution of Two-Dimensional Rota-
tions X( N / 2)  G1, N / 2  1,  N / 2 .X12( N / 2)  rX ,0,,0,,0

Parallel execution of two-dimensional rotations (base op- (12)

erations of algorithm) is one way to increase calculation
performance. This possibility is proposed in [7],[8] in rela- Schema below presents obtaining of vector X(N/2) =[rX,
tion of matrix decomposition. As this parallelism gives sig- 0, … , 0]T. using rotation of vector X12(N/2) to the direction
nificant decrease of subsequent calculations, below will be 
of axis x1 .
presented one realization of this approach for proposed al- As it can be seen from (8)÷(10) and from schema on Fig.
gorithm without pretensions of novelty. 4, after one division of given vector X the number of sub-
sequent base operations is decreased to N/2 regardless that
Let’s have N-dimensional vector X  x1, x2 , ... , xn T . the number of base operations remain the same – N-1. It is
Vector X can be presented as a sum of two vectors X1 and easy to find that separation of given vector X to more then
X2, every one of which has half coordinates of value zero two parts will give more decreasing of subsequent base op-
and half coordinates, that have the same values as the coor- erations.
dinates of given vector X, as follows:
 vector with dimension greater then 2. In example if
N=7 then difference between schema of AR and those,
given above, will be only this, that in first stage will
miss the last base operation and x7 will go straight to
base operation on second stage.
2) The number NR of stages (in which parallel execution
of base operations of AR have to be performed) is
equal to log2N whereas number of all base opera-
tions is N-1 – the same as when all rotations are sub-
sequent.
3) AR really increases calculation performance in case of
parallel execution of base operations in stages. If base
operations in stages are executed consecutively, the
only advantage of AR is decreasing of accumulation of
calculation errors (due to rounding e.g).

Figure 4 Schema for rotation of vector X, separated into two parts, to the

direction of axis x1

The limit of subsequent divisions is reached when every

one part has only two elements. In example if given vector
has dimension N=8 it can be separated to 4 parts of two
elements. Rotation of every one part of vector (as
two-dimensional vector) to the direction of one of coordi-
nates gives resultant vector, that has N/2 non-zero elements.
The same way this sub-vector of dimension N/2 can be sep-
arated to N/4 parts of two elements and rotation of every
one of this parts gives resultant vector with N/4 non-zero
elements. It’s easy to find that if N=2P, P (set of all nat-
ural numbers) then the number of subsequent rotations,
needed to obtain vector with only one non-zero element, is
equal to log2N.

Figure 5 Schema for accelerated rotation of 8-dimensional vector to the


direction of axis x1

In example if vector X has dimension N=8, the number of

needed subsequent rotations will be log28 = 3. Schema for
Figure 6 Matrices of stages for AR of 8-dimensional vector
accelerated rotation of 8-dimensional vector to the direction

of axis x1 will look as on figure above
It is important to note the following: For high dimensional vectors AR proposes significant de-
1) Accelerated rotation
 (AR) of given vector to the direc- creasing of subsequently executed base operations. In ex-
tion of axis x1 can be applied not only for vectors, ample if dimension of vectors is N=1024, the number of
that have dimension N=2P, PN, but for every one subsequently executed base operations is 1024-1, whereas
AR algorithm uses number of stages (which defines the They are matrices for parallel execution of two-dimensional
number of subsequent base operations) NStages = log21024 rotations in corresponding coordinate planes as follows:
= 10.
It is important to note, that rotation is reversible linear
 G 1  n  1.2k ,1  2.n  1.2k 1, n, k  (13)
N / 2k
operation, which mean that AR is reversible linear operation MS k 
too. Reversibility of rotation, in particularly AR, is the “key n 1
stone” of proposed NRMG algorithm.
When AR is used, rotation matrix MX and MY can be pre- In NRMG that uses AR (which is not obligate), matrices
sented as a multiplication of matrices MS1, MS2, … MSNR, of stages denote parallel execution of base operations.
which correspond to the stages of rotation, shown in Fig. 5. Schema, given below shows rotation of given
In example for vector with dimension N=8 matrices of 8-dimensional vector X to the direction of vector Y and
stages MS3.MS2.MS1 will look as in Fig. 6. where C1,1, S1,1, obtaining rotation matrix M, which performs this rotation.
C1,2, S1,2, C1,3, S1,3, C2,1, S2,1 ,C2,2, S2,2, C3,1, S1,3 are denota-
   
Every one of matrices MX and MY is calculated as multipli-
tion of cos i, j and sin i, j , index i is the number cation of log2N matrices of stages for which coefficients
of stage, and j is the number of base operation. Ci-j and Si-j are calculated using (6).
As it can be seen, matrices of stages MSk, k=1,
2, …,log2N (except the last one) are not Givens matrices.

Figure 7 Schema for rotation of 8-dimensional vector X to the direction of vector Y using NRMG algorithm and AR.

of every one pixel of which is equal to the brightness of

corresponding pixel of image I2, scaled with coeffi-
~
cient Y / X . If X and Y have the same norm, vectors Y
6. Reversible Transformation of Images and Y will be identical.
Using NRMG Algorithm Has been created Matlab program to test proposed NRMG
algorithm. Program consist function for accelerated rotation
fnAR and code, which uses this function to obtain matrix M,
Down we will give one example of using NRMG algo- which rotates given vector X to the direction of given vector
rithm for reversible image transformation. Y.
Let’s have two monochrome raster images I1 and I2
which have the same number of pixels N. As it is known
function R = fnAR(X)
[14], raster images are presented as dot matrix data structure,
N = length(X); %X have to be row vector (transposed)
so the two given images can be presented as two
R= eye(N); %Initial rotation matrix = Identity matrix
N-dimensional vectors X and Y, every element of which
step = 1; %Initial step
presents brightness of one of pixels. Using NMRG algo-
while(step<N) %Loop to create matrices of stages
rithm can be obtained rotation matrix M, which rotates
A= eye(N);
vector X to the direction of vector Y. Since rotation do not
~ n=1;
change the norm of the vector, resultant vector Y will
while(n<=N-step)
have the same norm as X. As a result, rotation will trans-
r2 = X(n)*X(n) + X(n+step)*X(n+step);
form image I1 to one presentation of image I2 , brightness
if r2 > 0
 r = sqrt(r2); are identical at least with precision of 10-6. It’s important to
pcos = X(n)/r; note, that vectors X and Y have the same norm.
psin = -X(n+step)/r;
% Base 2-dimensional rotation
A(n, n) = pcos; 7. Conclusion
A(n, n+step) = -psin;
A(n+step, n) = psin; In this article, we presents a new algorithm for generation
A(n+step, n+step) = pcos; of N-dimensional rotation matrix M, which rotates given
end N-dimensional vector X to the direction of given vector Y
n=n+2*step;
% Move to the next base operation which has the same dimension. Also examined is the possi-
end; bility to perform parallel calculation of base operations -
step = step*2; two-dimensional rotations. Algorithm can be used for dif-
X=(A*X')'; ferent reversible calculations. It is included one example of
R= A*R; % Multiply R by current matrix of stage A image transformation that uses NRMG algorithm.
end; Proposed NRMG algorithm proves that
end; 1) Every one matrix of rotation M, which rotates given
vector X to the direction of given vector Y can be pre-
sented as M  M Y1.M X where MX rotates vector X to
Example 1 Code of Matlab function for accelerated rotation of vector X to
the direction of axis x1. Function returns matrix of rotation MX
the direction of one of coordinate axes (e.g. axis x1) and
MY rotates vector Y to the direction of the same coor-
Code, that uses this function to obtain matrix M, which
dinate axis.
rotates given vector X to the direction of given vector Y is
2) Every one rotation of N-dimensional vector can be per-
given below:
formed by 2(N-1) two-dimensional rotations in N-1 co-
Mx = fnAR(X); ordinate planes.
My = fnAR(Y); As an advantage of NRMG algorithm to known solutions,
M = My'*Mx; %Obtaining rotation matrix M (shortly described in p. 2), can be pointed a possibility to
Z = M*X'; %Obtain vector Z to the direction of Y realise it as Discrete Linear System, elements in which per-
if round(Z'*100000)/100000==Y forms base operations. This possibility is a result of the lo-
disp('Z and Y are identical'); cality of data, used for base operations and independence of
end base operations inside one stage of rotation.
Example 2 Code of Matlab, which creates matrix of rotation using fnAR
function

For test data has been used the following two images:

REFERENCES
[1] H. G. Golub, J. M. Ortega, 1993, Scientific Computing and
Introduction with Parallel Computing, Academic Press, Inc.,
San Diego.
[2] G. A. Korn, T. M. Korn, 1961, Mathematical Handbook for
Scientists and Engineers (1st ed.), New York: McGraw-Hill.
pp. 55–79.
Figure 8 Test images [3] H. Friedberg, A. Insel, L. Spence, 1997, Linear Algebra (3rd
ed), Prentice Hall.
Images can be presented as 64-dimensional vectors as
follows (written using Matlab Language syntax for row [4] D. A. Harville, 1997, Matrix Algebra from Statistician’s
vectors): Perspective, Softcover.
[5] B. Buchberger, 1985, Multidimensional Systems Theory -
X =[0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 Progress Directions and Open Problems in Multidimensional
00001111000010010000100100000000 Systems, Reidel Publishing Company.
0 0]; [6] G. H. Golub, C. F. Van Loan, , 1996, Matrix Computations,
Y =[0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 4rd edition. Johns Horkins University Press, Baltimore
00001111000010000000100000000000
0 0 ]; [7] M. Cosnard, Y. Robert. Complexity of parallel QR factoriza-
tion. J. ACM 33, 4 (August 1986), 712-723.
DOI=10.1145/6490.214102, 1986.
Test program, given above, obtain vector Z to the direction
of vector Y using NRMG algorithm, and compare vectors Z [8] A. H. Sameh, D. J. Kuck. On Stable Parallel Linear System
and Y. Program displays message “Z and Y are identical”, Solvers. J. ACM 25, 1 (January 1978), 81-91. DOI=
http://dx.doi.org/10.1145/322047.322054, 1978.
which shows that two vectors (and corresponding images)
[9] N. J. Higham, 1996, Accuracy and Stability of Numerical
Algorithms, SIAM, Philadelfia.
 [10] G. W.Steward, 1976, The economical storage of plane rota-
tions, Numer. Math, 25, 2 1976, 137-139
[11] Matlock, H., and Reese, L.C., 1960, Generalized solutions for
laterally loaded piles., Journal of Soil Mechanics and Foun-
dation, 86(5), 63–91.
[12] N. K. Bose, 1985, Multidimensional Systems Theory: Pro-
gress, Directions, and Open Problems. D.Reidel Publishing
Co., Dordrecht, The Netherland.
[13] A. S. Householder, 1958, “Unitary triangularization of a
nonsimetric matrix”, J. ACM 5, 339-342, 1958.
[14] E. Anderson, 2000, “Discontinuous Plane Rotations and the
Symmetric Eigenvalue Problem” LAPACK Working Note
150, University of Tennessee, UT-CS-00-454, December 4,
2000. page 2.
[15] J. D. Foley, A. van Dam, S. K. Feiner, J. F. Hughes, 1995,
“Computer Graphics, Principles and Practice”, 2nd edition in
C, Addison-Wesley, ISBN 0-201-84840-6.
[16] G. Strang, 2006, Linear Algebra and its Applications,
Thomson Learning Ink, pages 69-135, ISBN 0-03-010567.
[17] St. Roman, Advanced Linear Algebra, second ed., 2005
Springer-Verlag, New York. pages 59-85, ISBN:
978-1-4757-2180-5
[18] J. E. Gentle, 2007, Matrix Algebra: Theory, Computations,
and Applications in Statistics, Springer, page 180, ISBN
978-0-387-70872-0
[19] I. R. Shafarevich, A. Remizov, 2013, Linear Algebra and
Geometry, Springer, pages 133-160, ISBN
978-3-642-30993-9

View publication stats