New Method of Givens Rotations for Triangularization of Square Matrices

This paper describes a new method of QR-decomposition of square nonsingular matrices (real or complex) by the Givens rotations through the unitary discrete heap transforms. This transforms can be defined by a different path, or the order of processing components of input data, which leads to different realizations of the QR-decomposition. The heap transforms are fast, because of a simple form of decomposition of their matrices. The direct calculation of the N-point discrete heap transform requires no more than 5(N − 1) multiplications, 2(N − 1) additions, plus 3(N − 1) trigonometric operations. The QR-decomposition of the square matrix N × N uses about 4/3N3


Introduction
multiplications and N(N − 1)/2 square roots.This number varies depending on the path of the heap transform, and it is shown that "the optimal path" allows for significant reduction of number of operations in QR-decomposition.The heap transform and its matrix can be described analytically, and therefore, this transform can also be applied to the QR-decomposition of complex matrices.

QR-Factorization, Givens Rotations, Householder Reflections, Heap Transform
In linear algebra, methods of QR-decomposition (or factorization) of a nonsingular matrix into a unitary matrix and a triangular matrix are well known in mathematics [1]- [3].QR-decomposition is used in many applications in computing and data analysis.This is the problem of solution of a linear system of equations written in matrix form as = Ax y .The solution 1 − = x A y can be found after the factorization of the matrix = A QR , where Q is an orthogonal matrix ( ) cients of the matrix R are positive.In this case, 1 − = x R Qy , or 1 − = x R PQy in the rank-revealing QR algo- rithm, when the diagonal elements of the matrix R are permuted P in the non-increasing order.There are several methods for computing the QR-decomposition, such as the Gramm-Schmidt process and method of Cholesky factorization.We here mention two other methods: the Householder transformations (known also as Householder reflections) and the Givens rotations.In the second method, each rotation zeros one element in the subdiagonal of the matrix.Therefore, a sequence of ( ) plane rotations is required for reduction of a square matrix ( ) N N × to triangular form.The Givens rotations require a large number of arithmetical operations, including multiplications and ( ) square roots [4].The method of Householder transforms is the most applied method for QR-decomposition, which reduces the number of square roots to at most ( ) 2 1 N − and uses about 3 4 3 N multiplications [5]- [8].In this paper, a new look on the application of Givens rotations to the QR-decomposition problem is described, which is similar to the method of Householder transformations.The concept of the discrete heap transform is applied, which has been introduced in digital signal processing to generate the signal-induced unitary transforms by Grigoryan [9]- [12].Both cases of real and complex matrices are considered and examples of performing the QR-decomposition of matrices are given.We also illustrate the importance of the path of the heap transform in such decomposition.The traditional way of consequently performing the rotations of data in natural order 1,2,3,... is not the best way, or path, in calculating the QR-decomposition.We briefly describe other more effective paths in QR-decomposition by the heap transforms and give a comparison with the known method of the Householder transformation.

Transforms with Decision Equations
Let ( ) , , f x y ϕ and ( ) , , g x y ϕ be functions of three variables; ϕ is referred to as the rotation parameter such as the angle, and x and y as the coordinates of the point ( ) , x y on the plane.These variables may have other meanings.The function ( ) , , g x y ϕ is parameterized and it is assumed that for a specified set of numbers a and for each point ( ) , x y on the plane (or its chosen subset) the equation ( ) has a unique solution with respect to ϕ .We denote the solution of this equation by ( ) , , r x y a ϕ = which is derived (or whose matrix is calculated) from the given decision equations by → is unitary.These transformations are basic stones to build unitary transforms.Figure 1 shows the graph of such a basic transformation in part a.The transformation T ϕ is defined by given x and y , namely, the angle ϕ is calculated first and then the transform ( ) , T x y ϕ is calculated.The second output is shown with the dash arrow, since the value of this output equals a which is known from given angu- lar Equation (1).The graph of this transform applied to an input ( ) Example 1: (Elementary rotation).Given a real number a , we consider the following functions: The basic transformation is defined as a rotation of the point ( ) where the rotation angle ϕ is calculated by ( ) The angular equation puts a constrain on parameter a, namely, it is required that 2 2 2 a x y ≤ + .

Discrete Heap Transforms
The composition of the N-dimensional discrete heap transform by a given vector-generator ( )  is performed by the sequential calculation of basic two-dimensional transformations.In matrix form, this composition can be written as follows: − , changes only two components of the input.We assume that the components of an input ( ) z − .This is a natural path P , and in general, such a path can be taken in many different ways.It is a very important characteristic of the heap transform and the right selection of the path leads to an effective application of the transform in calculating the QR-decomposition.We consider the case, when all basic transformations k T are parameterized by angles k ϕ .The special selection of a set of parameters is initiated by the vector-generator x through the decision equations with a given set of constants . The generator is processed first, and during this process all required angles are calculated.For a given set A, we define the following "heap" of elements: ( ) where angles k ϕ are calculated by ( ) that the first component of the vector participates in calculations of all ( ) As an example, Figure 2 shows the signal-flow graph of determination of the five-point transform by a vector x .
The N-point transformation T composing from k T in the space of N-dimensional vectors z is defined as , , , , , , , , ., , , , where ( ) The following should be mentioned about the heap transform.The transform is performed in a space of N-dimensional vectors z , but all an- gles k ϕ are found and the transformation T is composed after solving the decision equations relative to a given vector-generator x .The transform of x results in a vector with the constant components of the set A, plus ( ) , , , , .
x  The transformation T is called the N-point discrete x -signal-induced heap transformation (DsiHT), and the vector x is the generator of this transformation [6].The first component ( ) is referred to as the heap of the transform.

Transform Collecting Energy
We focus on the special case of the DsiHT, namely on transformations that collect the energy of vectors in one location.The angular equations for such transformations are defined by the set − .This condition leads to the fact that the first basis functions of such transforms are the vector- generators themselves.Other basis functions are defined based on the correlation data of components of these vectors.Matrices of these transforms are orthogonal.The transforms have simple forms of decomposition that lead to calculation of the N-point transforms with no more than ( ) operations are also required for calculating elementary trigonometric functions.In the linear space of N-dimension vectors, we construct such an orthogonal transform, T , whose matrix is defined by T E = x e where E denotes the energy of the vector and e is the unit vector ( ) The term "energy" is referred to as the norm of the vector, i.e., To define the desired transform T , we consider the following method of energy transferring.Let ( ) ′ be a 2-D vector to be rotated into a vector ( ) The matrix of the transformation, ( ) T ϕ , which is the Givens rotation, is defined by the matrix equation 0 0 The angle ϕ and component 0 y are calculated as ( ) + .When processing the N-dimensional vector x , the first pair of components, ( ) by Equation ( 2), as the point ( ) x x is rotated to the point ( ) 0 ,0 y on the horizontal line.Then, on the next kth step of calculations, when 1 k ≥ , the similar rotations are accomplished over vectors ( ) ( ) , where ( ) denotes a new value of the first component on the (k − 1)th step, and The value of the first component 0 y is renewed consequently as where the angles are defined by Here, m T are respectively the matrices of the Givens rotations ( ) − .As a result, the whole energy of the input signal is collected consequently and, then, transferred to the first component ( ) form of a vector ( ) , are calculated by the given vector ( ) and then used when transforming the vector

Analytical Expressions
It is important to note, that the basic transformations, ( ) − , that compose the N-point DsiHT, can be performed without calculation of angles and trigonometric functions.Analytical formulas can be derived for calculation of components of the heap transform ( ) , ,  , , , , .
The matrix of the transform can also be calculated analytically.To show this fact, we introduce the following notations which represent respectively the partial cross-correlation of z with the vector-generator x and energy of x : − .The components of the heap transform on the kth iteration can be expressed by correla- tion data as On the final step, the value of the first component is defined by ( ) ( ) ( ) x which is the cor- relation coefficient of the input signal z with the normalized signal-generator x .In the particular = z x case, we obtain ( ) ( ) x and the rest of coefficients ( ) ( ) of the matrix of the N-point DsiHT can be obtained from equations in (5).The mth row of the matrix of the transform is defined by applying the unit vector ( )   with 1 on the mth position, where integer − .Therefore, the coefficients can be calculated as

Transforms with Decision Equations
In this section, we give a brief analysis of the matrix decomposition by the heap transformations and the comparison with the known methods which are based on rotations in two dimensions.It should be noted first that the discrete heap transformations differ from the well-known Householder transformation [2], whose matrix ( ) ww , where the normalized vector w is calculated  is a given vector.For example, when ( )

T
These two transformations result in the same vector, ( )

Hx Tx
, where 5.4772 = x .Both transformations have the same first basis function which equals the normalized vector-generator.
The heap and Householder transformations differ much when considered with the large dimensions.For instance, when 7 N = and the vector-generator ( ) and the vector We now illustrate the difference of the DsiHT and Householder transforms applied on signals.Figure 4 shows the example with the 257-point signal ( ) in part a, along with the ( ) x t -induced heap and Householder transforms of ( ) y t respectively in b and c, and the difference of these two transforms in d.

Matrix Decomposition
We now consider the application of heap transforms for calculating the determinant of a nonsingular square matrix, as well as its decomposition.The method described in this section is based on the same idea as the method of Householder transforms, also known as the Householder reflections.Although the heap transformation and Householder transformation are different, it is evident that these two concepts should result in the same (or similar) decomposition, but by different ways.
Let X be a real square matrix (N × N) with ( ) . We denote coefficients of this matrix as Let vector a be the first column of this matrix, ( ) . We denote by a T the matrix of the heap transformation T a generated by this vector.This matrix contains an upper triangular submatrix with ( ) The product of two matrices 1 = a X T X results in a matrix ( ) . Therefore, the matrix 1 X has the following form: where we denote by 1 Y the ( ) 1,1 minor of this product.We can repeat the described above process for the submatrix 1 Y , by considering its first vector-column ( ) as a generator for the ( ) 1 N − - point heap tranformation T b , whose matrix we denote by b T .The product of two matrices ( ) On this step, the matrix X is transformed to the matrix 2 X which has the following form: where 2 Y is the (2,2) minor of this product and has the size ( ) ( )

X T X T T T X
is transformed into a matrix whose the first column is zero, except in the first row.The process of such transformations is repeated until we obtain the upper triangular matrix: since all heap transformations have determinant one.Thus, the of a square nonsingular matrix can be calculated by using the heap transforms.We denote by R the product of the following matrices ( ) ( )( )  )( ) ( )  We obtain the following decompositions of the matrix X and its inverse: 4 24 41 The decomposition of X by the heap transforms, 2 ′ = X R X , results in the unitary matrix 6

X
The coefficients of the matrix ′ R are given with precision of eight digits after the point.The decomposition of the matrix Thus, by means of the heap transforms which are composed by the Givens rotations, we achieve the same decomposition as the Householder reflections.It should be noted, that the first row-vector of the matrix R equals to the normalized first column-vector a of the original matrix X .For the example considered above, we obtain

Arithmetic Complexity of the Decomposition
The x -generated heap transform of a vector z can be calculated by using analytical Equations ( 5) and ( 6), If we use this estimation to calculate the number of multiplications for factorization of a square nonsingular matrix by heap transforms, we obtain the following number: ( )

Strong Heap Transform
The point, the heap, where the energy of the vector-generator is transferred is not required to be at the first location of the signal.The DsiHT can be defined by equation T E = x e for any unit vector ( )   .The path P along which we compose the heap is important, not a location of energy integration, when defining the discrete heap transform.As an example, we consider the basic transforms that sequentially process the last pair of components of the signal.This is the case, when the path P of the heap transformation induced by a vector-signal x is defined as ( ) − .We call this heap transformation the strong heap transformation.The energy from each component   is taken away and given to the last component

T
It is interesting to note that the energy of the signal is transferred sequentially to the first component in numbers 5, 5.3852, and 5.4772.This heap of numbers is smaller than the numbers 2.2361, 3.7417, and 5.4772, when the four-point DsiHT is composed along the ordinary path ( ) , by Equations ( 3) and ( 4).That can be explained because of small values of the first components of the signal x .
The method of triangularization of a square nonsingular matrix X by the strong heap transformations is de- scribed similarly to the method described in Section 4 for the ordinary path.The only difference is that on each kth step of calculations, where ( ) − , the generator of the strong heap transformation, which is zeroing the kth column of the renewed matrix k X , is taken from the top part of the column.To illustrate this method, we consider the matrix

X
The first column ( ) . The product of matrices a T and X equals 1 0 0.5930 1.7049 0 0.2774 0.5547 .
3.7417 5.3452 5.0780 Continuing zeroing the first two elements of the second column of this matrix, by using the strong heap transformation defined by the vector ( ) 0.5930,0.2774′ , we obtain the following decomposition of the matrix X by the product of the unitary matrix ′ R and the lower triangular matrix The determinant of X can be calculated as ( ) ( ) . Thus, we obtain the decomposition of the square matrix X by the product of the unitary matrix ′ R and lower triangular matrix 2 X .It should be noted for the comparison that when using the ordinary path, ( ) − , the following decompositions of the matrix X holds: It is evident that there is a relation between this decomposition and the decomposition of Equation ( 8).The first and third columns of the matrices ′ R have been rearranged, as well as the first and third rows of the ma- trices 2 X .In addition, the signs of the second column of ′ R and second row of 2 X have been changed.Many operations can be saved in the different steps of calculation of the QR-decomposition, if we learn how to select an optimal path in the heap transform.The optimality relates to minimization of the computation complexity of the QR-decomposition.
The calculations require six different square roots.Each component of the vector, except the first one, is processed only ones.We now consider another path 1 P in the heap transformation, which is shown in Figure 5 where the first column corresponds to the generator.
In this path, the transformation rounds three times the components of the input and intermediate data, and on each round a certain number of heaps are calculated, which then are transferred to one heap.During the first round, three different pairs of components ( ) , x x , and ( ) , x x are processed separately and three heaps of value 2 each are calculated, as shown in the second column.Component 3 x is not used in this stage since the dimension of the input is 7, odd number.In the second round, the first three obtained outputs together with 3 x are used to calculate two heaps with values of 3 and 4 .In the last round, these two heaps are rotated to collect the whole energy 7 in one heap.The calculation of the heap transform by the path 1 P is thus described by six Givens rotations as follows: ) , 2, 2,1, 0, 0, 0 step 4 : 2, 2, 2, 1,0,0,0 3, 2, 2,0,0,0,0 step 5 : 3, 2, 2,0,0,0,0 3, 4,0,0,0,0,0 step 6 : 3, 4,0,0,0,0,0 7,0,0,0,0 ,0,0 .
The underlined numbers in the input vectors show the pairs which are used in rotations in each step of calculation.The same number, six, of steps are used to accomplish the heap transform ( ) ( ) , however there are only four different square roots are calculated, 2, 3, 4 , and 7 .We can also define such a path for the general 4 N ≥ case, by composing first the transformation on the first round as and D is the empty operator if N is even, and the identity operator if N is odd.The transformations , k n T are the Givens rotations of two components of a vector with numbers k and m , when ( ) − .On the next round, the transform over the first half of outputs is composed si- milarly, Example 5: Let x be the vector ( ) This matrix is symmetric with respect to the middle column, and it has one triangle in the center, which is fully filled by zeros.We also can see two not fully filled by zeros triangles in both sides of the middle column.There are total 22 zeros in these triangles, more than in the matrix of the heap transformation defined by the path P. Indeed, the matrix of this transformation equals with 15 zero coefficients.Since the matrix 1 P H is sparse the application of the heap transformations with the path 1 P leads to the effective QR-decomposition of non-singular matrices.Indeed, in each step of the decom- position, the use of the heap transformation by the path 1 P requires fewer operations, when compared with the heap transformation by the path P , for large dimensions 3 N > .The matrix of the first transformation is sparser than the second one.For instance, in the (5 × 5) example, total 14 multiplications can be saved when using the paths 1 P in the QR-decomposition of the matrix ′ = X U T instead of P path.

Complex DsiHT
In this section, we describe briefly the case when heap transforms are generated by complex vectors and applied on complex vectors.For that we stand on Equation ( 7) which allows for defining the complex version of the heap transformation.Indeed, we can apply these equations wherein the correlation data are calculated as follows: Here, the bar denotes the operation of the complex conjugate.The transformation defined by these equations is unitary and is called the complex heap transformation.The basic two-dimensional transformation T, which is defined by a complex vector ( ) x x ′ , and then, is applied to a complex input ( ) , z z ′ is calculated in matrix form as follows: ( ) ( ) , we obtain the real transform , where 7.7460 x = .The matrix * x T is the complex conjugate to x T and its first row is the normalized generator, i.e., ( ) 1 3i, 2 i,3 4i, 4 2i 7.7460 + + + − .The complex heap transformations can be used for triangularization of the square complex matrix in a way similar to the real case described above.As an example, we consider the following matrix: 1 i 2 3i 3 4i 2 3i 3 2i 2 2i .3 i 4 3i 4 2i

X
The method of heap transforms results in the following decomposition of

Conclusion
A novel decomposition of nonsingular real and complex matrices by fast heap transforms has been described, when the paths P of the transforms are natural.The path of the heap transform is an important characteristic of the transform, and other paths P can be found, that may result in more sparse matrices than the path P does.The application of such paths will lead to effective matrix decomposition, as shown on examples with the strong heap transform.We have considered only the case of nonsingular matrices of square size.However, the calculations presented above can be applied to the case of non square matrices, as well.
system of decision equations.The value of ϕ is calculated from the second equation which is called the angular equation.Then the value of 0 y is calculated from the given input ( ) , x y and ϕ .It is also assumed that the two-dimensional transformation

Figure 1 .
Figure 1.Graphs of the basic transform T ϕ defined (a) by a point ( ) , x y and applied (b) on an input ( ) 0 1 , z z .

Figure 2 .
Figure 2. Signal-flow graph of determination of the five-point transform by a vector

Figure 3 Figure 3 .
Figure 3. Signal-flow graph of calculation of the four-point heap transform of the vector z.
The orthogonal matrix of the heap transformation generated by the vector x is not full filled and has the zero upper triangular submatrix with

Figure 4 .
Figure 4. (a) The original ( )x t and signal ( ) y t with a random noise, (b) the 257point discrete heap transform of ( ) y t , (c) the 257-point Householder transform of ( )y t , On the next step, we consider the first vector c of the submatrix 2 Y and construct the matrix c T of heap transformation generated by c .Then, the (3,3) minor 3 Y of the product of the matrices transformation is unitary, and R is unitary.The inverse matrix 1 − R can thus be represented by

X
are the upper triangular matrices, and R is unitary.Example 2: Consider the decomposition of the following matrix (3 × 3):

Figure 5 .
Figure 5. Signal-flow graph of the seven-point heap transformation .
that this transformation differs from the known definition of the complex Givens rotation[4], whose matrix is calculated as number with norm one.Example 6: Let 4 N = and let the complex vector x to be heap transformation generated by this vector has the following matrix: 0.1291 0.3873i 0.2582 0.1291i 0.3873 0.5164i 0.5164 0 where R is the unita- ry matrix 0.2000 0.2000i 0.4000 0.6000i 0.6000 0 The first seven values of the Householder and heap transforms are the following: : Hy 22.6276, −1.0608, −1.0285, −0.9915, −0.9503, −0.9052, −0.8568, : Ty 22.6276, 0.0190, 0.0366, 0.0566, 0.0788, 0.1031, 0.1291.The big difference of these transforms occurs in the first ninety components.
. The signal-flow graph of the heap transformation generated by x is given in Figure5.The orthogonal matrix of this transformation equals