A DFT-based approximate eigenvalue and singular value decomposition of polynomial matrices

In this article, we address the problem of singular value decomposition of polynomial matrices and eigenvalue decomposition of para-Hermitian matrices. Discrete Fourier transform enables us to propose a new algorithm based on uniform sampling of polynomial matrices in frequency domain. This formulation of polynomial matrix decomposition allows for controlling spectral properties of the decomposition. We set up a nonlinear quadratic minimization for phase alignment of decomposition at each frequency sample, which leads to a compact order approximation of decomposed matrices. Compact order approximation of decomposed matrices makes it suitable in filterbank and multiple-input multiple-output (MIMO) precoding applications or any application dealing with realization of polynomial matrices as transfer function of MIMO systems. Numerical examples demonstrate the versatility of the proposed algorithm provided by relaxation of paraunitary constraint, and its configurability to select different properties.


Introduction
Polynomial matrices have been used for a long time for modeling and realization of multiple-input multipleoutput (MIMO) systems in the context of control theory [1]. Nowadays, polynomial matrices have a wide spectrum of applications in MIMO communications [2][3][4][5][6], source separation [7], and broadband array processing [8]. They also have a dominant role in development of multirate filterbanks [9].
More recently, there have been much interest in polynomial matrix decomposition such as QR decomposition [10][11][12], eigenvalue decomposition (EVD) [13,14], and singular value decomposition (SVD) [5,11]. Lambert [15] has utilized Discrete Fourier transform (DFT) domain to change the problem of polynomial EVD to pointwise EVD. Since EVD is obtained at each frequency separately, eigenvectors are known at each frequency up to a scaling factor. Therefore, this method requires many frequency samples to avoid abrupt changes in adjacent eigenvectors.
Although, many methods of designing principle component filterbanks have been developed that are equivalent to EVD of pseudo circulant polynomial matrices [16,17], the next pioneering work on polynomial matrix EVD is presented by McWhirter et al. [13]. They use an extension of Jacobi algorithm known as SBR2 for EVD of para-Hermitian polynomial matrices which guarantees exact paraunitarity of eigenvector matrix. Since final goal of SBR2 algorithm is to have strong decorrelation, the decomposition does not necessarily satisfy spectral majorization property. SBR2 algorithm has also been modified for QR decomposition and SVD [10,11].
Jacobi-type algorithms are not the only proposed methods for polynomial matrix decomposition. Another iterative method for spectrally majorized EVD is presented in [14] which is based on the maximization of zeroth-order diagonal energy. Spectral majorization property of this algorithm is verified via simulation. Followed by the work of [6], a DFT-based approximation of polynomial SVD is also proposed in [18] which uses model order truncation by phase optimization.
In this article, we present polynomial EVD and SVD based on DFT formulation. It transforms the problem of polynomial matrix decomposition to the problem http://asp.eurasipjournals.com/content/2013/1/93 of, pointwise in frequency, constant matrix decomposition. At first it seems that applying inverse DFT on the decomposed matrices leads to polynomial EVD and SVD of the corresponding polynomial matrix. However, we will show later in this article that in order to have compact order decomposition, phase alignment of decomposed constant matrices in DFT domain results in polynomial matrices with considerably lower order. For this reason, a quadratic nonlinear minimization problem is set up to minimize the decomposition error for a given finite order constraint. Consequently, the required number of frequency samples and computational complexity of decomposition reduce dramatically. The algorithm provides compact order matrices as an approximation of polynomial matrix decomposition for an arbitrary polynomial order. This is suitable in MIMO communications and filterbank applications, where we deal with realization of MIMO linear time invariant systems. Moreover, formulation of polynomial EVD and SVD in DFT domain enables us to select the property of decomposition. We show that if eigenvalues (singular values) intersect at some frequencies in frequency domain, smooth decomposition, and spectrally majorized decomposition are distinct. The proposed algorithm is able to reach to either of these properties.
The remainder of this article is organized as follows. The relation between polynomial matrix decomposition and DFT matrix decomposition is formulated in Section 2. In Section 3, two important spectral properties of decomposition, namely spectral majorization and smooth decomposition, are provided using appropriate arrangement of singular values (eigenvalues) and corresponding singular vectors (eigenvectors). The equality of polynomial matrix and dft matrix decomposed matrices decompositions are guaranteed via the finite duration constraint, which is investigated in Section 4. The finite duration constraint imposes the phase angles of singular vector (eigenvector) to minimize a nonlinear quadratic function. A solution for this problem is proposed in Section 5. Section 6 presents the results of some computer simulations which are considered to demonstrate performance of the proposed decomposition algorithm.

Notation
Some notational conventions are as follows: constant values, vectors, and matrices are in regular character lower case, lower case over-arrow, and upper case, respectively. Coefficients of polynomial (scalar, vector, and matrix) are with indeterminate variable n in the square brackets. Any polynomial (scalar, vector, and matrix) is distinguished by bold character and indeterminate variable z in the parenthesis and its DFT by bold character and indeterminate variable k in the brackets.

Problem formulation
Denote a p × q polynomial matrix A(z) such that each element of A(z) is a polynomial. Equivalently, we can indicate this type of matrix by coefficient matrix A[ n], Define the effective degree of A(z) as N max − N min (or the length of A[ n] as N max − N min + 1). The polynomial matrix multiplication of a p × q matrix A(z) and a q × t matrix B(z) is defined as We can obtain the coefficient matrix of product by matrix convolution of A[ n] and B [ n], that is defined as where * denotes the linear convolution operator.
Denote para-conjugate of a polynomial matrix as in which, * as a subscript denotes the complex conjugate of coefficients in the polynomial matrix A(z).
A matrix is said to be para-Hermitian ifÃ(z) = A(z) or equivalently A[ n] = A H [ −n]. We call a polynomial matrix paraunitary ifŨ(z)U(z) = I, where I is a q × q identity matrix.
Thin EVD of a p × p para-Hermitian polynomial matrix A(z) is of the form and thin SVD of a p × q arbitrary polynomial matrix is of the form, where U(z) and V(z) are p × r and q × r paraunitary matrices, respectively. (z) and (z) represent r × r diagonal matrices where r is the rank of A(z). We can equivalently write EVD of a para-Hermitian matrix and SVD of a polynomial matrix in coefficient matrix form In general, EVD and SVD of a finite-order polynomial matrix are not finite order. As an example, suppose EVD of para-Hermitian polynomial matrix Eigenvalues and eigenvectors of the polynomial matrix in (6) are neither of finite order nor rational The same results can be found for polynomial QR decomposition in [12].
We mainly explain the proposed algorithm for polynomial SVD, yet wherever it seems necessary we explain the result for both decomposition.
The decomposition in (3) can also be approximated by samples of discrete-time Fourier transform, yields a decomposition off the form Such a decomposition can be obtained by taking the K- where w K = exp(−j2π/K). DFT formulation plays an important role in decomposition of polynomial matrices because it replaces the problem of polynomial SVD that involves many protracted steps with K conventional SVD that are pointwise in frequency. It also enables us to control spectral properties of the decomposition. However, it causes two inherent drawbacks: 1. Regardless of what is the trajectory of polynomial singular values in frequency domain, conventional SVD order singular values irrespectively of the ordering in neighboring frequency samples. 2. In frequency domain, samples of polynomial singular vectors are known up to a scalar complex exponential by using the SVD at each frequency sample, which yields to discontinuous variation between neighboring frequency samples.
The first issue is directly dealt with the spectral properties of the decomposition. In Section 3, we would explain why arranging singular values in decreasing order yields to approximate spectral majorization, while smooth decomposition requires rearrangement of singular values and their corresponding singular vectors.
For the second issue, suppose conventional SVD of an arbitrary constant matrix A. If the pair u and v are the left and right singular vectors corresponding to a non-zero singular value, for an arbitrary scalar phase angle θ, the pair e jθ u and e jθ v are also left and right singular vectors corresponding to the same singular value. Although this non-uniqueness is trivial in conventional SVD, it plays a crucial role in polynomial SVD. When we perform SVD at each frequency of DFT matrix as in (7), these nonuniquenesses in phase exist at each frequency regardless of other frequency samples. Denote the ith column vector of the desired matrices U(z) and V(z). Then all the vectors of the form have the chance to appear as the ith column of U [ k] and V [ k], and ith diagonal element of [ k], respectively. Moreover, in many applications, specially those which are related to MIMO precoding, we can relax constraints of the problem by letting singular values to be complex (see applications of polynomial SVD in [4,18]) Given this situation, singular values have not all their conventional meaning. For instance, the greatest singular value is conventionally 2-norm of the corresponding matrix, which is not true for complex singular values. The process of compensating singular vectors for these phases is what we call phase alignment and is developed in Section 4. Based on what was mentioned above, Algorithm 1 gives the descriptive pseudo code for DFT-based SVD. Modifications of the algorithm for EVD of para-Hermitian matrices are straightforward. If at each frequency sample all singular values are in decreasing order, REARRANGE function (which is described in Algorithm 2) is only required for smooth decomposition, otherwise for spectral majorization, no further arrangement is required. For the phase alignment, first we need to compute phase angles which is indicated in the algorithm by DOGLEG function and is described in Algorithm 3.
In many filterbank applications which are dealt with principle components filterbank, spectral majorization and strong decorrelation are both required [16]. Since smooth decomposition leads to more compact decomposition, in cases that the only objective is strong decorrelation, exploiting smooth decomposition is reasonable. The DFT-based approach of polynomial matrix decomposition is capable of decomposing a matrix with either of these properties with small modification.
Polynomial EVD of a para-Hermitian matrix is said to have spectral majorization property if [13,16] Note that, eigenvalues corresponding to para-Hermitian matrices are real in all frequencies.
If we let singular values to be complex, we can replace absolute value of singular values in the definition.
A polynomial matrix have no discontinuity in frequency domain, hence we modify definition of smooth decomposition presented in [19] to fit with our problem and avoid unnecessary discussions.
Polynomial EVD (SVD) of a matrix is said to possess smooth decomposition if eigenvectors (singular vectors) have no discontinuity in frequency domain, that is where u il is the lth element of u i . If eigenvalues (singular values) of a polynomial matrix intersect at some frequencies, the spectral majorization and smooth decomposition are not simultaneously realizable. As an example, suppose A(z) is a polynomial matrix with u 1 (z) and u 2 (z) are eigenvectors corresponding to distinct eigenvalues λ 1 (z) and λ 2 (z), respectively. Lets assume u 1 (e jω ) and u 2 (e jω ) have no discontinuity in frequency domain, and λ 1 (e jω ) and λ 2 (e jω ) intersect at some frequencies. Denote Obviously, u 1 (e jω ) and u 2 (e jω ) are eigenvectors corresponding to distinct eigenvalues λ 1 (e jω ) and λ 2 (e jω ), respectively. Note that, λ 1 (e jω ) ≥ λ 2 (e jω ) for all frequencies, which means λ 1 (e jω ) and λ 2 (e jω ) are spectrally majorized. However, u 1 (e jω ) and u 2 (e jω ) are discontinuous at intersection frequencies of λ 1 (e jω ) and λ 2 (e jω ), which implies that they are not smooth anymore. In this situation, although λ 1 (e jω ), λ 2 (e jω ), u 1 (e jω ), and u 2 (e jω ) are not even analytic, we can approximate them with finite order polynomials.
If a decomposition has spectral majorization, its eigenvalues (singular values) are of decreasing order in all frequencies. Therefore, they are in decreasing order in any arbitrary frequency sample set, including DFT frequencies. Obviously the converse is only approximately true. Hence, for polynomial EVD to possess spectral majorization approximately, it suffices to arrange sampled eigenvalues (singular values) of (7) in decreasing order. Since we only justify spectral majorization at DFT frequency samples, the resulting EVD (SVD) may possess the property only approximately. Similar results can be seen in [14,20].
To have smooth singular vectors, we propose an algorithm based on inner product of consecutive frequency samples of singular vectors. We can accumulate smoothing requirement in (11) for all r elements as Let B be the upper bound of norm of derivative and {·} be the real value of a complex value.
For an arbitrary ω we have that is, for a smooth singular vector u i (e j(ω+ ω) ) u i (e jω ) can be made to be as close to unity as desired by making ω sufficiently small. In our problem u i (e jω ) is sampled uniformly with ω = 2π K . Since EVD is performed at each frequency sample independently, u i [ k] and u i [ k + 1] are not necessarily two consecutive frequency samples of a smooth eigenvector. Therefore, we should rearrange eigenvalues and eigenvectors to yield smooth decomposition. This can be done for each sample of eigenvector u i [ k] by seeking for the eigenvector of successor sample } is not possible before phase alignment. Due to (15), for sufficiently small ω, two consecutive samples of a smooth singular vector can be as close as desired and we can approximate which allows us to use inner product of u (12) and (13)

Finite duration constraint
Phase alignment is critical to have compact order decomposition. Another aspect of this fact is revealed in the coefficient's domain perspective of (7). In this domain, the multiplication is replaced by circular convolution in which is the circular convolution operator and ((n)) K denotes n module K. http://asp.eurasipjournals.com/content/2013/1/93 Polynomial SVD corresponds to linear convolution in the coefficients domain, however the decomposition obtained from DFT corresponds to circular convolution. Recalling from discrete-time signal processing, it is well known that we can equivalently utilize circular convolution instead of linear convolution if convoluted signals are zero-padded adequately. That is, for x 1 [ n] and x 2 [ 2] are two signals with the length of N 1 and N 2 , respectively, apply zero padding such that zero padded signals have the length N 1 + N 2 − 1 [21]. Hence, if the last M − 1 coefficients of U[ n], [ n], and V [ n], are zero, the following results are hold: Therefore, the problem is to obtain the phase set { θ i [ k] } and correcting the singular vectors using (9). The phase set { θ i [ k] } should be such that the resulting coefficients satisfy (17). Without for n = M, M + 1, . . . , K − 1, in which K ≥ 2M − 1. If these conditions are satisfied, circular convolution can be used instead of linear convolution.
Since the available matrix of singular vectors at each frequency is U [ k], inserting (9) in (18) for each singular vector separately leads to for n = M, M + 1, . . . , K − 1. Without loss of generality, let θ i [ 0] = 0. In a more compact form we can express these (K − M)-folded equations in matrix form in which For polynomial EVD, Equation (20) An additional degree of freedom is obtained by letting singular values to be complex. However, an straightforward solution which yield to singular values and singular vectors of order M is complicated. Instead, we impose the finite duration constraint only two singular vectors If K ≥ 2M − 1, then the last M − 1 coefficients of resulting polynomial vectors are zero. Therefore, according to (17), U(z) and V(z) are paraunitary. On the other hand, if K ≥ 2M + N max − N min − 1, circular convolution relation of coefficient M. Another option which yields to more accurate results is by calculatingŨ(z)A(z)V(z) and replacing off-diagonal elements with zero.
Next, we provide a minimization approach to determine the unknown set {θ i [ k] }.

Gradient descent solution
In general, there may exist no phase vector θ which satisfies (20). Even when there exists a phase vector that satisfies the finite duration constraint, the solution is not straightforward. For these reasons, we can view (20) as a minimization problem [6]. We use energy of the highest order coefficients (the coefficients that we equate to zero in (18)) as the objective to the minimization problem An alternative minimization technique as a solution for this phase optimization problem is proposed in [6], which we describe it in this section.
Throughout this section, we focus on solving θ i = arg min J( θ i ) as a least square solution for a single singular vector u i [ k], so we drop the subscript "i" from the quantity θ i and use F and f , instead of F M ( u i ) and f M ( u i ) to simplify the notation. The objective J( θ) is intentionally presented as a function of θ to emphasize the fact that our problem is classified as an unconstrained optimization.
We exploit the trusted region strategy for the problem (23). By utilizing the information about the gradient vector and Hessian matrix in each step, trusted region strategy constructs a model function m k which have a similar behavior close to the current point θ (k) . The model m k is usually defined as the second-order Taylor series expansion (or its approximation) of J( θ + ϕ) around θ , that is where ∇J( θ) and ∇ 2 J( θ) are the gradient vector and the Hessian matrix corresponding to J( θ), respectively. The model m k is designed to be a good approximation of J( θ) near the current point and is not trustworthy on regions far from the current point. Consequently, the restriction in minimization of m k on a region around θ (k) is crucial, that is where R is the trusted region radius. The decision about shrinking of the trusted region is determined by comparing the actual reduction inobjective function and predicted reduction. Given a step ϕ, the ratio is used as a criterion to indicate if the trusted region is small enough. Among methods which approximate the solution of the constrained minimization (24) dogleg procedure is the only one which leads to analytical approximation. It also promises to achieve at least as much reduction in m k as is possible by Cauchy point (the minimizer of m k along the steepest descent direction −∇J( θ), subject to the trusted region) [22]. However, this procedure requires Hessian matrix (or an approximation of it) to be positive definite.

Hessian matrix modification
The gradient vector and Hessian matrix corresponding to J( θ) are as follows where X( θ) is a diagonal matrix with the kth diagonal element exp(−jθ[ k] ) and k = 0, 1, . . . , K − 1.
In general, Hessian matrix in (26) does not promise to be always positive definite. Therefore, we should modify Hessian matrix to yield a positive definite approximation.
We provide a simple modification which brings some desirable features by omitting the second term from the Hessian matrix and diagonal loading to guarantee positive definiteness The term 2 X( θ)F H FX( θ) H is positive semi-definite and in many situations, it is much more significant than the second term of Hessian matrix in (26). Hence, with diagonal loading αI (I is with conformable size and α is very small), the modified Hessian matrix guarantees (27) to be positive definite and provides the desired properties in contrast to use the exact Hessian matrix.

Dogleg method
Dogleg method starts with the unconstrained minimization of (24) When the trusted region radius is so large that φ H ≤ R, it is the exact solution of (24) and we select it as the dogleg method answer. On the other hand, for small R the solution of (24) is −R g/ g . For intermediate values of R, the optimal solution lies on a curved trajectory between these two points [22]. http://asp.eurasipjournals.com/content/2013/1/93

Algorithm 3 Trusted region dogleg algorithm
Obtain θ 0 from (33) The dogleg method approximates this trajectory by a path consisting of two line segments. The first line segment starts from the origin to the unconstrained minimized point along the steepest descent direction The second line segment starts from φ g to φ h . These two line segments form an approximate trajectory which its intersection with the sphere φ = R is the approximate solution of (24) when φ h > R.

Alternating minimization
Another solution of (23) is provided by converting the problem of multivariate minimization to a sequence of single-variate minimization problem via alternating minimization [6]. In each iteration, a series of single-variate minimization is performed, while other parameters are held unchanged. Each Iteration consists of K − 1 steps, which at each step one parameter θ[ k] is updated. Suppose we are at step k of ith iteration. At this step k − 1 first parameters were updated in the current iteration, and K − k − 2 last parameters were updated from the previous iteration. These parameters are held fixed, while θ [ k] is minimized at the current step, The cost function is guaranteed to be non-incremental at each step; however, this method is also converges to a local minima which highly depend on the initial guess of the algorithm. For solving (30) it is suffices to make the kth element of gradient vector in (26) equal to zero. Suppose the calculation are performed for phase alignment of where ∠t Fortunately, Equation (31) has a closed form solution However, only the second case of (32) has positive second partial deviation. Therefore, the global minima of (30) is

Initial guess
All algorithms of unconstrained minimization require to be supplied by a starting point, which we denoted by θ 0 .
To avoid getting stuck in local minima, we should select a good initial guess. This can be accomplished by minimizing a different but similar cost function denoted by J ( θ) Solving J ( θ) yields to a simple initial guess Based on what have been mentioned in this section, a pseudo-code description of the trusted region dogleg algorithm is given by Algorithm 3. In this algorithm, http://asp.eurasipjournals.com/content/2013/1/93 we start with the initial guess of (33) and a trusted region radius upper boundR. Then we continue the trusted region minimization procedure as described in this section.

Simulation results
In this section, we present some examples to demonstrate the performance of the proposed algorithm. For the first example, our algorithm is applied to a polynomial matrix example from [11] Frequency behavior of singular values can be seen in Figure 1. There is no intersection of singular values, so the setup of the algorithm either for spectral majorization or frequency smoothness leads to identical decomposition.
For having approximately positive singular values, we use (21). Define the average energy of highest order coefficients for the pair of polynomial singular vectors u i and v i as E u,v i = J( θ i )/(K − M) (we expect energy of highest order coefficients to be zero or at least minimized). A plot of E i versus iteration for each pair of singular vectors is depicted in Figure 2. The decomposition length is M = 9 (order is 8) and we use K = 2M + (N max − N min ) = 20 number of DFT points.
As it is seen, the use of dogleg method with approximate Hessian matrix leads to a fast convergence in contrast with using alternative minimization and Cauchy-point (which is always selected along the gradient direction). Of course we should consider that due to matrix inversion, computational complexity of Dogleg method is O(K 3 ) while computational complexity of alternative minimization and Cauchy point is O(K 2 ). The final value of average highest order coefficient for three pair of singular vectors are 5.54 × 10 −5 , 3.5 × 10 −3 , and 0.43, respectively. The first singular vector satisfies finite duration constraint almost exactly. The second singular vector fairly satisfies this constraint. However, highest order coefficients of last singular vector, possess considerable amount of energy, that seems to cause decomposition error.
Denote the relative error of the decomposition as in which · F is the extension of Frobenius norm for polynomial matrices and is defined by Since in our optimization procedure we only seek for finite duration approximation, U(z) and V(z) are only approximately paraunitary. Therefore, we also define relative error of paraunitarity as An upper bound for E U can be obtained as which means as average energy on K − M highest order goes to zero, E U diminishes. The relative error of this decomposition is E A = 1.18 × 10 −2 while the error of U(z) and V(z) are E U = 3.3 × 10 −2 and E V = 3.08 × 10 −2 , respectively. The paraunitarity error is relatively high in contrast with decomposition error. This is due to the difference between the first two singular values and the last singular value.
A plot of relative errors E A , E U , and E V for various amount of M is shown in Figure 3. The number of frequency samples is fixed at K = 2M + 2(N max − N min ).
The number of frequency samples K is an optional choice, however as discussed in Section 4, it should satisfy K ≥ 2M + N max − N min − 1. In order to demonstrate the effect of number of frequency samples on the decomposition error, a plot of relative error versus different amount of K is depicted in Figure 4. Increasing the number of frequency samples does not lead to reduction of relative error. Moreover, it increases computational burden. Therefore, a value near 2M + (N max − N min ) − 1 is a reasonable choice for the number of frequency samples. Now, lets relax the problem by allowing singular values to be complex and using (22). A plot of E u i and E v i versus iteration for each pair of singular vectors is depicted in respectively. This is while these values for right singular vectors are 1.12 × 10 −10 , 1.4 × 10 −3 , and 8.7 −4 , respectively.
Note that the average energy of highest order coefficients for the third pair of singular vectors alleviate meaningfully. Figure 1 shows that the third singular value goes to zero and then returns to positive values. If we constrain singular values to be positive, a phase jump of π radian, is imposed to one of third singular vectors near the frequency which singular vector goes to zero. However, by letting singular values to be complex, the zero crossing occur which requires no discontinuity of singular vectors.
The relative error of this decomposition is E A = 4.9 × 10 −3 while the error of U(z) and V(z) are E U = 2.5 × 10 −3 and E V = 3.5 × 10 −3 , respectively. In contrast with constraining singular values to be positive, having complex singular values decrease decomposition and paraunitarity error significantly.
Plots of relative errors E A , E U , and E V for various amount of M and K are shown in Figures 6 and 7, respectively. Letting singular values be complex causes significant reduction of all relative errors. As it was mentioned, Figure 7 shows that increasing K from 2M+N max − N min − 1 causes no improvement in relative errors while it makes additional computational burden.
McWhirter and coauthors [11] have reported the relative error of decomposition. Provided that paraunitary matrices U(z) and V(z) are of order 33, the relative error of their algorithm is 0.0469. This is while our algorithm only requires paraunitary matrices of order 3 for relative error of 0.035 with positive singular values and relative error of 2.45 × 10 −6 with complex singular values. In addition, This large difference is not caused by iteration numbers because we compare results while all algorithms relatively converges, and with continuation of iterations trivial improvement are obtained. The main reason lies on different constraints of the solution presented in [11] in contrast to our proposed method. While they impose paraunitary constraint onŨ(z)A(z)V(z) to yield a diagonalized (z), we impose the finite duration constraint and obtain approximation of U(z) and V(z) with fair fitting to the decomposed matrices at each frequency samples. Therefore, we can consider this method as a finite duration polynomial regression of matrices which is obtained by uniformly sampling U(z) and V(z) on the unit circle in z-plane.
As a second example, consider EVD of the following para-Hermitian matrix The exact smooth EVD of this matrix is of finite order Frequency behavior of eigenvalues can be seen in Figure 8. Since eigenvalues intersect at two frequencies, smooth decomposition and spectrally majorized decomposition result two distinct solutions.
To perform smooth decomposition, we need to track and rearrange eigenvectors to avoid any discontinuity using Algorithm 2. The resulting |c u ij [ k] | are shown in Figure 9 for k = 0, 1, . . . , K − 1. Using these |c u ij [ k] | the Algorithm 2 swap first and second eigenvalues and eigenvectors for k = 12 : 32 which results in continuity of eigenvalues and eigenvectors. Now that all eigenvalues and eigenvectors are rearranged in DFT domain, it's time for phase alignment of eigenvectors. A plot of E i versus iteration for M = 3 and smooth decomposition is depicted in Figure 10. It is predictable that dogleg algorithm converges rapidly while the alternative minimization and Cauchy point has a long way to converge. Since the energy of highest order coefficients of eigenvectors are trifling, using the proposed method for smooth decomposition results in very high accuracy, as seem in the figures. Relative error of smooth decomposition versus M is shown in Figure 11.
While using frequency smooth EVD of (35) leads to relative error below 10 −5 for M ≥ 3 with a few number of iterations, Spectrally majorized EVD requires a lot more polynomial order to reach a reasonable relative error.
Unlike smooth decomposition which requires rearrangement of eigenvalues and eigenvectors, spectral majorization requires only to sort eigenvalues at each frequency sample in decreasing order. Most of conventional EVD algorithm sort eigenvalues in decreasing order, which we should only align eigenvector phases using 3. A plot of E i versus iteration for M = 20 and spectrally majorized decomposition is depicted in Figure 12.
Due to an abrupt change in eigenvectors at the intersection frequency of eigenvalues, increasing the decomposition order leads to a slow decay of relative error. Figure 13 shows the relative error as a function of M.
To see the difference between smooth and spectrally majorized decomposition results, eigenvalues of spectrally majorized decomposition is shown in Figure 14, which is comparable with Figure 8 which corresponds to eigenvalues of smooth decomposition. Therefore, a low order polynomial is required using smooth decomposition and much higher polynomial order for spectrally majorized decomposition. Even with M = 20 the decomposition have relatively high error.

Conclusion
An algorithm for polynomial EVD and SVD based on DFT formulation has been presented. One of the advantages of the DFT formulation is that it enables us to control properties of decomposition. Among these properties, we introduce how to setup the decomposition to achieve spectrally majorization and frequency smoothness. We have shown, if singular values (eigenvalues) intersect at some frequency, then simultaneous achievement of spectral majorization and smooth decomposition is not possible. In this situation, setting up the decomposition to possess spectral majorization requires considerably higher order polynomial decomposition and more computational complexity. Highest order polynomial coefficients of singular vectors (eigenvectors) are utilized as square error to obtain a compact decomposition based on phase alignment of frequency samples. The algorithm has the flexibility to compute a decomposition with approximately positive singular values, and a more relaxed decomposition with complex singular values. A solution for this nonlinear quadratic problem is proposed via Newton's method. Since we apply an approximate Hessian matrix to assist the Newton optimization, a fast convergence is achieved. The algorithm capability to control the order of polynomial elements of decomposed matrices and to select properties of decomposition, make the proposed method as a good choice for filterbank and MIMO precoding applications. Finally, performance of the proposed algorithm under different conditions is demonstrated via simulations. Simulation results reveal superior decomposition accuracy in contrast with coefficient domain algorithms due to relaxation of paraunitarity.