2.6 Subspace Tracking Algorithms

It is seen from Section 2.5 that the linear multiuser detectors are obtained once the signal subspace components are identified. The classic approach to subspace estimation is through batch eigenvalue decomposition (ED) of the sample autocorrelation matrix or batch singular value decomposition (SVD) of the data matrix, both of which are computationally too expensive for adaptive applications. Modern subspace tracking algorithms are recursive in nature and update the subspace in a sample-by-sample fashion. An adaptive blind multiuser detector can be based on subspace tracking by sequentially estimating the signal subspace components and forming the closed-form detector based on these estimates. Specifically, suppose that at time ( i - 1), the estimated signal subspace rank is K [ i - 1] and the components are U _s [ i - 1], L _s [ i - 1], and s ² [ i - 1]. Then at time i , the adaptive detector performs the following steps to update the detector and to detect the data.

Algorithm 2.6: [Blind adaptive linear MMSE detector based on subspace tracking ”synchronous CDMA]

• Update the signal subspace: Using a particular signal subspace tracking algorithm, update the signal subspace rank K [ i ] and the subspace components U _s [ i ], L _s [ i ], and s ² [ i ].

• Form the detector and perform detection:

  

  
  

Various subspace tracking algorithms are described in the literature (e.g., [43, 87, 97, 406, 412, 461, 493, 586]). Here we present two low-complexity subspace tracking algorithms: the PASTd algorithm [586] and the more recently developed NAHJ algorithm [412].

2.6.1 PASTd Algorithm

Let be a random vector with autocorrelation matrix C _r = E { r [ i ] r [ i ] ^H }. Consider the scalar function

Equation 2.150

with a matrix argument ( r < N ). It can be shown that [586]:

• W is a stationary point of if and only if W = U _r Q , where contains any r distinct eigenvectors of C _r and is any unitary matrix.

• All stationary points of are saddle points except when U _r contains the r dominant eigenvectors of C _r . In that case, attains the global minimum.

Therefore, for r = 1, the solution of minimizing is given by the most dominant eigenvector of C _r . In real applications, only sample vectors r [ i ] are available. Replacing (2.150) with the exponentially weighted sums yields

Equation 2.151

where 0 < b < 1 is a forgetting factor. The key issue of the PASTd (projection approximation subspace tracking with deflation) approach is to approximate W [ i ] ^H r [ n ] in (2.151), the unknown projection of r [ n ] onto the columns of W [ i ], by y [ n ] = W [ n “ 1] ^H r [ n ], which can be calculated for 1 n i at time i . This results in a modified cost function

Equation 2.152

The recursive least-squares (RLS) algorithm can then be used to solve for the W [ i ] that minimizes the exponentially weighted least-squares criterion (2.152).

The PASTd algorithm for tracking the eigenvalues and eigenvectors of the signal subspace is based on the deflation technique and can be described as follows . For r = 1, the most dominant eigenvector is updated by minimizing in (2.152). Then the projection of the current data vector r [ i ] onto this eigenvector is removed from r [ i ] itself. Now the second most dominant eigenvector becomes the most dominant one in the updated data vector and it can be extracted similarly. This procedure is applied repeatedly until all the K eigenvectors are estimated sequentially.

Based on the estimated eigenvalues, using information theoretic criteria such as the Akaike information criterion (AIC) or the minimum description length (MDL) criterion [557], the rank of the signal subspace, or equivalently, the number of active users in the channel, can be estimated adaptively as well [585]. The quantities AIC and MDL are defined as follows:

Equation 2.153

Equation 2.154

where M is the number of data samples used in the estimation. When an exponentially weighted window with forgetting factor b is applied to the data, the equivalent number of data samples is M = 1/(1 - b ). a ( k ) in the definitions above is defined as

Equation 2.155

The AIC (respectively, MDL) estimate of subspace rank is given by the value of k that minimizes the quantity (2.153) [respectively, (2.154)]. Finally, the PASTd algorithm for both rank and signal subspace tracking is summarized in Table 2.1. The computational complexity of this algorithm is operations per update. The convergence dynamics of the PASTd algorithm are studied in [587]. It is shown there that with a forgetting factor b = 1, under mild conditions, this algorithm converges globally and almost surely to the signal eigenvectors and eigenvalues.

Simulation Examples

In what follows we provide two simulation examples to illustrate the performance of the subspace blind adaptive detector employing the PASTd algorithm.

Example 1: Performance Comparison Between Subspace and MOE Blind Detectors This example compares the performance of the subspace-based blind MMSE detector with the performance of the MOE blind adaptive detector. It assumes a real-valued synchronous CDMA system with a processing gain N = 31 and six users ( K = 6). The desired user is user 1. There are four 10-dB multiple-access interferers (MAIs) and one 20-dB MAI (i.e., = 10 for k = 2, ..., 5 and = 100 for k = 6). The performance measure is the output SINR. For the PASTd subspace tracking algorithm, we found that with a random initialization, the convergence is fairly slow. Therefore, in the simulations, the initial estimates of the eigencomponents of the signal subspace are obtained by applying an SVD to the first 50 data vectors. The PASTd algorithm is then employed thereafter for tracking the signal subspace. The time-averaged output SINR versus number of iterations is plotted in Fig. 2.10.

Figure 2.10. Performance comparison between a subspace-based blind linear MMSE multiuser detector and an RLS MOE blind adaptive detector. The processing gain N = 31. There are four 10-dB MAIs and one 20-dB MAI in the channel, all relative to the desired user's signal. The signature sequence of the desired user is a m -sequence, whereas the signature sequences of the MAIs are randomly generated. The signal-to-ambient noise ratio after despreading is 20 dB. The forgetting factor used in both algorithms is 0.995. The data plotted are averages over 100 simulations.

Table 2.1. PASTd Algorithm [585, 586] for Tracking the Rank and Signal Subspace Components of Received Signal r [ i ].

Updating the eigenvalues and eigenvectors of signal subspace
	x 1 [ i ]	=	r [ i ]
FOR	k = 1: K i -1		DO
	y k [ i ]	=	u k [ i -1] H x k [ i ]
	l k [ i ]	=	b l k [ i -1] + y k [ i ] 2
	u k [ i ]	=	u k [ i -1] + ( x k [ i ] - u k [ i -1] y k [ i ] ) y k [ i ] * / l k [ i ]
	x k +1 [ i ]	=	x k [ i ] - u k [ i ] y k [ i ]
END
	s 2 [ i ]	=	b s 2 [ i -1] + x K i-1 + 1 [ i ] 2 /( N-K i -1 )
Updating the rank of signal subspace K i [a]
FOR	k = 1: K i -1		DO
	a ( k )	=
	AIC( k )	=	( N - k ) ln [ a ( k )] /(1 - b )+ k (2 N - k )
END
	K i	=	arg min k N -1 AIC( k ) + 1
IF	K i < K i -1		THEN
ELSE IF	K i > K i -1		THEN
	u K i [ i ]	=	x K i -1 +1 ( i )/ x K i -1 +1 [ i ]
	l K i [ i ]	=	s 2 [ i ]
END

^[a] The rank estimation is based on the Akaike information criterion.

As a comparison, the simulated performance of the recursive least-squares (RLS) version of the MOE blind adaptive detector is also shown in Fig. 2.10. It has been shown in [389] that the steady-state SINR of this algorithm is given by

Equation 2.156

where SINR ^* is the output SINR value of the exact linear MMSE detector, and is the forgetting factor). Hence the steady-state SINR performance of this algorithm is upper bounded by 1/ d . This behavior can be observed in Fig. 2.10. Although an analytical expression for the steady-state SINR of the subspace-based blind adaptive detector is very difficult to obtain, as the dynamics of the subspace tracking algorithms are fairly complicated, it is seen from Fig. 2.10 that with the same forgetting factor, the subspace blind adaptive detector well outperforms the RLS MOE detector. Moreover, the RLS MOE detector has a computational complexity of per update, whereas the complexity per update of the subspace detector is

Example 2: Tracking Performance in a Dynamic Environment This example illustrates the performance of the subspace blind adaptive detector in a dynamic multiple-access channel, where interferers may enter or exit the channel. The simulation starts with six 10-dB MAIs in the channel; at t = 2000, a 20-dB MAI enters the channel; at t = 4000, the 20-dB MAI and three of the 10-dB MAIs exit the channel. The performance of the proposed detector is plotted in Fig. 2.11. It is seen that this subspace-based blind adaptive multiuser detector can adapt fairly rapidly to the dynamic channel traffic.

Figure 2.11. Performance of a subspace-based blind linear MMSE multiuser detector in a dynamic multiple-access channel where interferers may enter or exit the channel. At t = 0, there are six 10-dB MAIs in the channel; at t = 2000, a 20-dB MAI enters the channel; at t = 4000, the 20-dB MAI and three of the 10-dB MAIs exit the channel. The processing gain N = 31. The signal-to-noise ratio after despreading is 20 dB. The forgetting factor is 0.995. The data plotted are averages over 100 simulations.

2.6.2 QR-Jacobi Methods

QR-Jacobi methods constitute a family of SVD-based subspace tracking algorithms that rely extensively on Givens rotations during the updating process. This reduces complexity and has the advantage of maintaining the orthogonality of matrices. Members of this family include the algorithms presented in [340, 398, 461].

Let

Equation 2.157

denote an N x ( i + 1) matrix whose columns contain the exponentially windowed first i snapshots of the received signal. The sample autocorrelation matrix of Y [ i ] and its eigendecomposition are given by

Equation 2.158

Alternatively, the SVD of the data matrix Y [ i ] is given by

Equation 2.159

where in both (2.158) and (2.159) the columns of U _s [ i ] contain the eigenvectors that span the signal subspace and contains the square roots of the corresponding eigenvalues. Generally speaking, SVD-based subspace tracking algorithms attempt to track the SVD of a data matrix of growing dimension, defined recursively as

The matrix V [ i ] need not be tracked. Furthermore, since the noise subspace does not need to be calculated for the blind multiuser detection algorithm, we do not need to track U _n [ i ]. This allows us to reduce complexity using noise averaging [224]. Since calculating the SVD from scratch at each iteration is time consuming and expensive, the issue then is how to best use the new measurement vector, r [ i + 1], to update the decomposition in (2.159).

Noise-averaged QR-Jacobi algorithms begin with a Householder transformation that rotates the noise eigenvectors such that the projection of the new measurement vector r [ i + 1] onto the noise subspace is parallel to the first noise vector, which we denote by u _n . Specifically, let

Equation 2.160

Equation 2.161

where g = r [ i + 1] - U _s [ i ] r _s . Then we may write the modified factorization

Equation 2.162

where represents the subspace of U _n [ i ] that is orthogonal to u _n . The second step in QR-Jacobi methods, sometimes called the QR step , involves the use of Givens rotations to zero each entry of the measurement vector's projection onto the signal subspace. We refer the reader to [176] for details concerning the use of Givens matrices for this purpose. The QR step replaces the last row in the middle matrix in the decomposition in (2.162) with zeros. These are row-type transformations involving premultiplication of the middle matrix with a sequence of orthogonal matrices. We do not need to accumulate these transformations in V [ i ] since it does not need to be tracked.

The next step, the diagonalization step, involves at least one set each of column and row-type rotations to further concentrate the energy in the middle matrix along its diagonal. Sometimes called the refinement step , this is where many of the existing algorithms begin to diverge. The algorithm in [398], for example, performs two fixed sets of rotations in the diagonalization step but leaves the middle matrix in upper triangular form and does not attempt a true diagonalization. This is particularly efficient for applications that do not require a full set of eigenvalues but is not useful here. The algorithm in [398], on the other hand, attempts to optimize the choice of rotations to achieve the best diagonalization possible.

2.6.3 NAHJ Subspace Tracking

The algorithm we present here was developed in [411, 412]. It is a member of the QR-Jacobi family in the sense that it uses Givens rotations during the updating process. However, this algorithm avoids the QR step entirely. Instead of working with the SVD-type decomposition in (2.159), we work with the eigendecomposition of the form

Equation 2.163

where S ² [ i ] is Hermitian and almost diagonal. This is simply the eigendecomposition (2.158) except that we have relaxed the assumption that L [ i ] is perfectly diagonal. At each iteration we use a Householder transformation and a vector outer product to update S ² [ i ] directly. We then use a single set of two-sided Givens rotations to partially diagonalize the resulting Hermitian matrix. There is no need for a separate QR step. Essentially, the diagonalization process used in this algorithm is a partial implementation of the well-known symmetric Jacobi SVD algorithm [176] (not to be confused with the family of QR-Jacobi update algorithms). This algorithm is used to find the eigenstructure of a general fixed symmetric matrix and is known to generate more accurate eigenvalues and eigenvectors than the symmetric QR SVD algorithm but with a higher computational complexity [313]. However, we do not perform the full sweep of K ( K - 1)/2 rotations required for the symmetric Jacobi algorithm, only a carefully selected set of about K rotations. This is sufficient because the matrix that we wish to diagonalize already has much of its energy concentrated along the diagonal. This is a situation that the Jacobi algorithm can take advantage of but which the QR algorithm cannot. The Jacobi algorithm also has an inherent parallelism, which the QR algorithm does not. Table 2.2 contains a summary of this algorithm, which we term NAHJ (noise-averaged Hermitian “Jacobi) subspace tracking.

Table 2.2. NAHJ Subspace Tracking Algorithm

Given: , and U _s [ i - 1]

1. Calculate r _s , u _n , and g according to (2.160) and (2.161).

2. Dropping the indices, generate the modified factorization:

   

   
   

3. Let Y _s be the K + 2 principal submatrix of the matrix sum in step 2. Apply a sequence of r + 1 Givens rotations to Y _s to produce

4. Set equal to the K + 1 principal submatrix of Y _a .

5. Let U _s [ i ] be composed of the first K+1 columns of

6. Reaverage the noise power:

   where

7. Let L _s [ i ] be the diagonal matrix whose diagonal is equal to the first K elements of the diagonal of Y _a .

Algorithm

The first step in NAHJ subspace tracking is the Householder transformation mentioned previously. The second step involves generating a modified factorization that maintains the equality

Equation 2.164

Step 3 requires that we apply K + 1 Givens rotations to partially diagonalize Y _s . Ideally, we would apply these rotations to those off-diagonal elements having the largest magnitudes . However, since the off-diagonal maxima can be located anywhere in Y _s , finding the optimal set of rotations requires an search for each rotation. This leads to an complexity algorithm. To maintain low complexity we have implemented a suboptimal alternative that is simple, yet effective. Let be the vector whose outer product is used in the modified factorization of step 2. Suppose that i , 1 i K + 2, is the index of the element in z that has the largest magnitude. The set of elements we choose to annihilate with the Givens rotations is given by . Of course, if ( Y _s ) _{i , j} is annihilated, so is ( Y _s ) _j,i . This choice of rotations is not optimal; in fact, since we retain the off-diagonal information from the previous iteration we cannot even be sure that we annihilate the off-diagonal element in Y _s with the largest magnitude. Nevertheless, we see that the technique is very simple and is somewhat heuristically pleasing. Ultimately, performance is the measure of merit, and simulations show that it performs very well. The total computational complexity of the NAHJ subspace tracking algorithm is per update.

To adapt to changes in the size of the signal subspace (number of users), the tracking algorithm must be rank-adaptive. As before, both the AIC and the MDL criteria can be used for this purpose. To use this algorithm, we must track at least one extra eigenvalue “eigenvector pair ”hence the appearance of K + 1 in Table 2.2.

Simulation Example

This example compares the performance of the subspace blind adaptive multiuser detector using the NAHJ subspace tracking algorithm with that of the LMS MOE blind adaptive multiuser detector. It assumes a synchronous CDMA system with seven users ( K = 7), each employing a gold sequence of length 15 ( N = 15). The desired user is user 1. There are two 0-dB and four 10-dB interferers. The performance measure is the output SINR. The performance is shown in Fig. 2.12. It is seen that the subspace blind detector significantly outperforms the LMS MOE blind detector, in terms of both convergence rate and steady-state SINR. Further applications of the NAHJ subspace tracking algorithm are found in later chapters (cf. Sections 2.7.4, 3.5.2, 5.5.4, and 5.6.3). Simulations show that the NAHJ subspace tracking algorithm substantially outperforms the PASTd algorithm, especially in multipath environments. Note that both algorithms have a computational complexity of .

Figure 2.12. Performance comparison between a subspace blind adaptive multiuser detector using the NAHJ subspace tracking algorithm and an LMS MOE blind adaptive multiuser detector.