EURASIP Journal on Applied Signal Processing 2002:12, 1387–1400 c ○ 2002 Hindawi Publishing Corporation Reduced-Rank Adaptive Filtering Using Krylov Subspace

A unified view of several recently introduced reduced-rank adaptive filters is presented. As all considered methods use Krylov subspace for rank reduction, the approach taken in this work is inspired from Krylov subspace methods for iterative solutions of linear systems. The alternative interpretation so obtained is used to study the properties of each considered technique and to relate one reduced-rank method to another as well as to algorithms used in computational linear algebra. Practical issues are discussed and low-complexity versions are also included in our study. It is believed that the insight developed in this paper can be further used to improve existing reduced-rank methods according to known results in the domain of Krylov subspace methods.


INTRODUCTION
Adaptive filtering is widely used in signal processing applications such as array signal processing, equalization, and multiuser detection (see [1,2,3]). Least-square adaptive filters gained considerable attention during the last three decades and numerous algorithms were proposed [4].
The frequent problem which arises when designing an adaptive filtering system is that large observation size, and therefore, large filter length, means inevitably high computational cost, slow convergence, and poor tracking performance. However, this situation corresponds to many important practical applications such as high data rate directsequence code division multiple access (DS-CDMA) systems, radar or global positioning system (GPS) array processing. Reduced-rank adaptive filters provide a way out of this dilemma [3,5]. The basic idea behind the rank reduction is to project the observation onto a lower-dimensional subspace usually defined by a set of basis vectors. The adaptation is then performed within this subspace with a low-order filter resulting in substantial computational savings, better convergence, and tracking characteristics.
The MSWF takes its origin in a decomposition of an n-dimensional minimum mean square error (MMSE) filter into a linear combination of a full rank (n-dimensional) matched filter and a reduced-rank ((n − 1)-dimensional) MMSE filter. The latter may be further expanded into a (n − 1)-rank matched filter, (n − 2)-rank MMSE filter, and so on. A k-stage MSWF, which approximates the original MMSE filter, is obtained by taking the first k matched filters of the decomposition as basis vectors set. Like any other reduced-rank method, the MSWF projects the observation onto the subspace spanned by basis vectors and filters the result with a low-rank MMSE filter. AVF schemes can be viewed as iterative procedures. The whole process is initialized by the matched filter. At each step, the current filter is linearly combined with an "auxiliary" filter. Basis vectors set of an AVF is composed of the original matched filter and of the set of all auxiliary vectors. Finally, the CGRRF builds the basis set as sequence of vectors which are orthogonal in the metric defined by the received covariance matrix.
Interestingly enough, basis vectors of the MSWF, the AVF, and the CGRRF span the same subspace [12,13]. This subspace can be shown to be Krylov subspace which is generated by taking the powers of the covariance matrix of observations on a cross-correlation (steering) vector [5]. This direct procedure of basis set generation is used in the POR receiver [9] and in the Cayley-Hamilton receiver [14] (the latter is developed in the context of centralized detection). Therefore, in the MSWF, the AVF, the CGRRF, and the POR receiver observations are projected on the same (Krylov) subspace. Our presentation, which further develops this point, is organized as follows.
In Section 2, data model is presented. Although our numerical studies use reduced-rank methods in the context of CDMA multiuser detection, the model is sufficiently general to encompass other applications, for example, radar [15] and GPS [16] space-time adaptive signal processing. Section 3 briefly introduces the reader into the problems of adaptive Wiener filtering which motivated the development of reduced-rank methods. In Section 4, the family of reducedrank methods which are based on Krylov subspace is presented and their interrelationship is discussed. Krylov subspace methods are widely used in computational linear algebra for solving large (possibly sparse) systems of linear equations [17,18,19]. For that reason, the effort is made to relate presented algorithms to known algorithms of computational linear algebra. Specifically, forward recursion of the MSWF is shown to be an interpretation of the wellknown Lanczos algorithm and AVF are derived from the steepest descent algorithm. In Section 5, the results of numerical studies are presented and Section 6 concludes our work.

DATA MODEL
Throughout the paper, the notations * , T, and H are used to denote the conjugate, transpose, and conjugate transpose operations, respectively.
Let r(k) = r 1 (k) r 2 (k) · · · r N (k) T be the N × 1 vector consisting of N data samples observed at time instant k, which is modeled as where s(k) denotes the M × 1 vector of source signals s 1 (k), s 2 (k), . . . , s M (k), H is the N × M channel matrix and n(k) stands for the N × 1 noise vector. In the sequel, s(k) and n(k) are supposed to be zero-mean and wide-sense stationary with the respective covariance matrices E[s(k)s H (k)] = diag( 1 , 2 , . . . , M ) and E[n(k)n H (k)] = R n . The model (1) can be used, for example, to represent M narrowband sources impinging on an N-element antenna array, or in the context of a symbol-synchronous DS-CDMA system. In the latter case, s(k) is the vector of signals transmitted by M system users and the ith column of channel matrix H represent the channel signature of user i, that is, ith spreading code convolved with ith channel impulse response. It can be easily seen that all methods discussed in this work apply (with minor modifications) to more complex models, such as an asynchronous CDMA system with strong intersymbol interference (ISI) or multirate systems.

MOTIVATION FOR THE REDUCED-RANK METHODS
Consider the problem of estimation of the source signal s 1 (k) given the observation (1). General linear estimator can be written asŝ where w is an N × 1 vector (filter). The well-known full-rank Wiener filter [2] is the solution of the following linear system (normal equations): where c = E[r(k)s * 1 (k)] is the desired signal data crosscorrelation vector and R = E r(k) r H (k) is the covariance matrix of r(k). The important property of the Wiener filter is that it is the only filter that minimizes the mean square error (MSE) or, in other words, average error energy. In our notations, the MSE can be written as The Wiener filter owns its popularity not only to this property but also to its relatively simple expression as a solution of a linear system (3). However, in most practical applications such as array signal processing and CDMA multiuser detection, exact values of the covariance matrix and of the cross-correlation vector are not available. For example, in a synchronous CDMA system, such characteristics as number of CDMA users, user spreading codes, user fading, and the signal-to-noise ratio are partially or completely unknown. In radar signal processing applications, little information may be available about bearings and powers of jammer signals and so forth. Moreover, noise and signal powers as well as the channel matrix H may exhibit slow variations. 1 Therefore, we have to deal with some estimates of R and c. By way of example, the estimate of R can be estimated as where 0 < γ < 1 is the forgetting factor. As soon as exact values of R and c are replaced by the time-varying estimates R(k) and c(k), system (3) has to be resolved. Each time these estimates are updated in order to take into account the most recent samples of r(k). Taking into consideration high values of N, we usually encounter in practice, may be quite a problem from the computational viewpoint when using direct inversion of R(k). Moreover, as the system to solve has the form natural questions arise such as the convergence of w(k) to the Wiener filter w N opt , convergence speed, and the behaviour of the solution w(k) in a nonstationary environment (tracking ability or "adaptivity" of the filter). These questions can only be answered while taking into account the particular method of solving (6). Unfortunately, the answers provided by conventional adaptive filtering techniques (sample matrix inversion or the recursive least square (RLS) algorithm [2]) are often unsatisfactory for applications where the amount of training data (that is, the number of observations) is limited [10], such as radar space-time adaptive processing or CDMA multiuser detection in fast fading environment.
Reduced-rank methods, as an alternative to sample matrix inversion, provide fast and efficient (approximate) solutions to (6).

REDUCED-RANK ADAPTIVE FILTERING USING KRYLOV SUBSPACES
In this section, we develop a new insight to recently introduced reduced-rank filtering techniques such as the MSWF, which is based on Krylov subspace methods of computational linear algebra. Our main motivation is to establish a link between these two branches of recent scientific research and to see how this relationship can be used in the context of adaptive filtering.

Filter rank reduction
The reduced-rank Wiener filter in subspace D is defined as The above definition includes the full rank Wiener filter as a particular case when D = Ꮿ N . Let {q j }, j = 1, . . . , D be an orthonormal basis of D . Define Q def = q 1 q 2 · · · q D . As w D opt = Qµ for some µ ∈ Ꮿ D , (7) can be rewritten as Substituting w = Qµ into (4) yields where the transformed covariance matrix R t and the transformed signal data cross-correlation vector c t are defined as It then follows that µ D opt in (8) is the solution of Therefore, the reduced-rank Wiener filter is found by solving (11) and substituting µ D opt into (8).
Contrary to (3), (11) is a system of D linear equations. Therefore, confining the filtering operation to a lowdimensional subspace D leads to substantial gains in complexity when D N. Better convergence and tracking properties can also be expected [10,20]. For example, with rank reduction, we only need 2D observations (compared to 2N for the full rank filter) to be within 3 dB from the maximum filter-output signal-to-interference-plus-noise ratio (SINR) [21].
On the other hand, confining the Wiener filter to a lowdimensional subspace implies the loss of degrees of freedom of the filter and, therefore, this operation should increase the MMSE achieved by a reduced-rank method As for the complexity, the computational overhead due to eventual estimation of Q also has to be taken into account. Therefore, "good" choice of D (and of the rank-reduction method) is always a compromise dictated by the requirements of a given application.

Reduced-rank filtering using Krylov subspace
Definition 2. Given a square matrix A and a nonzero vector v, the subspace defined by is referred to as a Dth Krylov subspace associated with the pair (A, v) and is denoted D (A, v) [17].
This work deals with a family of reduced-rank methods for which D = D (R, c). The natural question is: what kind of reasoning leads to this particular choice for D ? To answer this question, consider the gradient of the MSE (4) Now take an arbitrary i-dimensional subspace i . Let w i opt be the reduced-rank Wiener filter in i , that is, Suppose that we seek to extend the subspace i to an (i + 1)-dimensional subspace i+1 . Since J(w) decreases most rapidly in the direction of −∇J(w), a reasonable strategy is to require that It follows from (14) that, for the condition above to be satisfied, it is sufficient for i+1 to contain the pair (c, Rw i opt ). Now let { i , i = 1, 2, . . . , D} be a chain of Krylov subspaces, that is, i = i (R, c), i = 1, 2, . . . , D. It is then easy to prove by induction that in this case, condition (16) is satisfied for each i within the range 1, . . . , D. Therefore, the Krylov subspace D (R, c) results from D steps of a sequential procedure, which (i) is initialized with the matched filter ( 1 = c); (ii) at step i, solves the reduced-rank minimization problem (15) and extends the minimization subspace i with the gradient of the cost function (MSE) taken at the point w i opt .
Remark 1. In many signal processing applications, the observation noise is modeled as being white. It can be shown that in this case the full-rank Wiener filter w N opt lies in M (R, c) where M is the number of source signals. However, M can still be high (e.g., M = N/2 or even exceeds N in heavily or over loaded CDMA systems) and therefore, it defines only an upper bound for D.
Remark 2. Other approaches leading to Krylov subspaces can be found in literature. For example, consider the polynomial decomposition of R −1 : A reduced-rank filter is obtained by truncating the righthand side of (17) to D terms and by multiplying the result by c, The coefficients {α i } are chosen in order to minimize the MSE (the Cayley-Hamilton Receiver of [14]) or to maximize the signal-to-interference ratio [22]. In [6], the MSWF is developed through the decomposition of the full rank Wiener filter into a linear combination of the matched filter c N and of the reduced-rank Wiener filter v N−1 opt in the orthogonal to c N subspace The filter v N−1 opt can be further represented as a linear combination of the matched filter c N−1 (in the subspace orthogonal to c N ) and of the Wiener filter v N−2 opt of rank N −2 (in the subspace orthogonal to span{c N , c N−1 }), and so on. 2 The vectors c i so obtained again generate the Krylov subspace.

Exact algorithms
In this section, three reduced-rank methods are introduced. The common feature of all presented techniques is that they perform exact MSE minimization in the Krylov subspace D (R, c). Hence, these exact methods are mathematically equivalent and result in the same performance. Practical issues such as complexity, structural flexibility, and the robustness to rounding errors become of major importance when preferring one exact method to another (see the discussion in Section 4.4).

The powers of R (POR) receiver
The POR receiver, proposed by Honig and Xiao in [9], can be considered as the simplest Krylov subspace method of The algorithm is summarized in Algorithm 1. It is noteworthy that for the POR receiver, [i, j]th element of R t can be written as therefore, R t is a Hankel matrix.

The multistage Wiener filter (MSWF)
The MSWF [6], shown in Figure 1 for rank D = 4, consists of two distinct iterative procedures. The first one (forward recursion) builds an orthonormal basis of the Krylov subspace D (R, c) giving the projection matrix Q = q 1 q 2 . . . q D . The second procedure (backward recursion) solves system (11) giving the transformed Wiener filter µ D opt or, equivalently, the weighting of basis vectors.

Forward recursion
An orthonormal basis of D (R, c) can be constructed, for example, by applying the general Gram-Schmidt procedure to the basis (20). However, for this particular basis and for the Hermitian-symmetric R, the Gram-Schmidt orthonormalization can be simplified. The Lanczos algorithm [17] (see Algorithm 2) represents an efficient way to compute the orthonormal basis {q i } of the Krylov subspace D (R, c).
Forward recursion results from the following interpretation of the Lanczos algorithm. Define Figure 1: MSWF (rank D = 4).

Initialization:
Algorithm 2: The Lanczos algorithm. Then p i (see Algorithm 2) can be expressed as

Initialization:
Therefore, the vector p i can be viewed as the matched filter for estimating d i−1 (k) using the reduced-rank signal r i−1 (k) [3]. Assembling together the formula for q i from the Lanczos algorithm (line 1 of Algorithm 2) and (22) and (23), we obtain the forward recursion of rank D MSWF (Algorithm 3).

Backward recursion
The orthonormal basis {q i , i = 1, 2, . . . , D} 4 has a remarkable property: the transformed covariance matrix is tridiagonal [6,18]. More precisely, where Note also that the transformed signal data cross-correlation vector becomes In backward recursion of MSWF, the above-mentioned structural properties are exploited in order to solve the system Indeed, using (24) and (26), we can write system (27) in an expanded form (µ D opt = µ 1 µ 2 . . . µ D T ) The last equation can be rewritten as where Substituting (29) into where Similarly, we may show that where Equations (29), (30), (34), and (35) are sufficient to solve system (27).
Another expression for ω i can be obtained as follows. Define Then where To show (38), it is sufficient to use definition (36) and the identities (which follow from the tridiagonality of R t ): Assembling (36), (37), and (38) gives another form of backward recursion which is shown in Algorithm 4.

The conjugate-gradient reduced-rank filter (CGRRF)
The CGRRF [7] or conjugate-gradient implementation of the MSWF [8] are inspired directly from the conjugategradient algorithm (CGA) for systems of linear equations.
For that reason, we start by a brief introduction into conjugate-gradient methods. For a more detailed presenta- Algorithm 4: Backward recursion of the rank D MSWF.
tion, the reader is referred to standard textbooks on computational linear algebra [18,19]. Consider the following general iterative procedure: with the sequences of complex coefficients c i and of vectors u i chosen according to some optimization criterion.
The criterion considered here is J(w i ), so it is natural to require that J(w i ) ≤ J(w i−1 ). Note also from (42) that w i is always in ᐁ i = span{u 1 , u 2 , . . . , u i }. The question is whether it is possible to choose c i and u i to give the reduced-rank Wiener filter in ᐁ i or not. In other words, we require that The following lemma answers this question.
Lemma 1. For the requirement (43) to be satisfied, it is sufficient that (1) u i are mutually R-conjugate, that is, (2) c i is given by where Proof. See [18].
It is easy to show that the value of the coefficient c i as given by (45) minimizes the MSE in the direction of the line ᏸ = {w i−1 + cu i }. Therefore, condition (44) guarantees that the reduced-rank Wiener filter w i opt lies on ᏸ. Different versions of the CGA result from different ways to compute the sequence of R-conjugate vectors u i [18]. The version shown in Algorithm 5 requires only one matrix-byvector multiplication per iteration. After D iterations of the algorithm, the sequence {w i opt } of D reduced-rank Wiener filters in ᐁ i is generated.  The following lemma establishes the equivalence between the CGA and other exact methods (MSWF, POR).
Proof. See [18]. Therefore, the reduced-rank Wiener filter in ᐁ i generated at the ith CGA iteration is also the reduced-rank Wiener filter in the Krylov subspace i−1 (R, c).
Basically, the CGRRF of rank D performs D CGA iterations. The CGRRF has a multistage structure, as shown in Figure 2, with the stage i computing the reduced-rank Wiener filter of the rank i and filtering the received signal to give the estimateŝ i 1 (k).

Discussion
The complexity of block implementations of the each exact reduced-rank method (in multiplications per block) is given in Table 1. In this table, T denotes the block size and it is assumed that the covariance matrix R is estimated as For the reduced-rank methods, the complexity comprises the computation of basis vectors of the Krylov subspace and of the Wiener filter in the Krylov subspace. The Hankel structure of R t has been accounted for in the complexity estimate for the POR receiver. The complexity of the block sample matrix inversion (SMI) method, which solves system (3) through the Cholesky decomposition ofR, and the complexity of RLS algorithm [2] are also given for reference. It follows from Table 1 that exact methods offer significant complexity reduction over the SMI and RLS algorithms as D N. The analysis presented in [5,23] shows that in practice D = 8 is sufficient to attain the full rank (RLS) performance for wide range of system loads (M/N) and signalto-noise ratios (SNRs) ( i ).
When some sample estimate R(k) is used, exact algorithms (POR, MSWF, CGRRF) provide the exact minimum of the "sample" MSE cost function in the "sample" Krylov subspace D R(k), c = span c, R(k)c, . . . , R D−1 (k)c .
As the minimization subspace and the cost function are common to these algorithms, they are mathematically equivalent and they ideally result in the same reduced-rank solution. However, the ways of obtaining this solution are different which may have an impact when implementing exact algorithms using finite precision arithmetic. For example, POR basis vectors t i , i = 1, 2, . . . , D are not orthonormal. For small sample sizes, these basis vectors are nearly dependent [9] and the matrix R t often becomes ill-conditioned. In [9], an adaptive rank selection technique was proposed to overcome this difficulty. Another obstacle to practical use of the POR algorithm is significant disparity in norms of t i . Indeed, according to (20), t i grows exponentially with i thus, complicating the fixed point implementation of the algorithm. Clearly, dependence of basis vectors and norm disparity can be eliminated by some kind of orthonormalisation with additional computational cost. Compared to the MSWF and the POR receiver, the CGRRF has the advantage of computing Wiener filters of all ranks ranging from 1 to D, therefore, in the CGRRF filter, outputs (symbol estimates) of different ranksŝ i 1 (k), i = 1, . . . , D are simultaneously available, as shown in Figure 2.
This property simplifies real-time filter rank selection by measuring the SINR at the output of each stage and adapting the filter rank D to achieve the given target SINR. Moreover, the CGRRF of any rank i < D is always at hand while for the MSWF and POR system (11) has to be resolved for each value of i (for the MSWF, weighting coefficients ω i resulting from backward recursion have to be recomputed).

Auxiliary-vector filters (AVF)
Consider the system of normal equations (3). Note that the Wiener filter, which minimizes the MSE given by (4), also minimizes the function Starting from an arbitrary filter w 0 , the Newton-Raphson [19] iteration gives the exact solution to (3) Substituting R = ∇ 2 J 0 (w 0 ) into the above expression yields where v 1 is the solution of It is clear that to solve the original MMSE minimization problem, we have to provide the solution to (53). Suppose that the following approximation is used instead: where the unit norm auxiliary vector g 1 and constant c 1 are defined as Substituting (54) into (51) gives It can be shown that c 1 as given by (56) minimizes J 0 (w) along the line ᏸ = {w 0 + cg 1 }. The "approximate" Newton-Raphson iteration is therefore an iteration of the method of steepest descent [18,19]. Clearly, we may continue to iterate as where Initialization: Substituting ∇J 0 (w) = Rw−c in (60) and setting w 0 = c/ c leads to the algorithm given in Algorithm 6, denoted here as AVF.
It should be noted that our derivation is based on an unconstrained minimization of the MSE. In [11], the filter output energy (w i ) H Rw i is minimized under the constraint (w i ) H c = 1. The resulting algorithm, denoted here as the minimum variance AVF (MVAVF), computes auxiliary vectors as where J (w) = w H Rw. Clearly, the term (I − cc H ) in (63) guarantees the orthogonality between g i and c and leads to the constant response in the direction of c, Also, we may impose the orthogonality between g i , These constraints lead to The algorithm which employs g i , as defined above, is denoted here as the constrained AVF (CAVF), and its minimum variance counterpart, originally derived in [10], will be referred to as the CMVAVF (see Algorithm 7). It is easy to see that, for all variants of the AVF, w i ∈ span{c, g 1 , g 2 , . . . , g i } = i (R, c). Moreover, for the CMVAVF, g i = q i , i = 1, 2, . . . , D, where q i (Lanczos vectors) are generated by the algorithm of Algorithm 2.Therefore, basis vectors of the CMVAVF and of the MSWF coincide. The latter fact, however, does not imply the equivalence of an AVF with D auxiliary vectors to an exact method of rank D+1. The main difference resides in the manner the reducedrank filter solution is obtained. Indeed, exact methods perform the MSE minimization over the whole Krylov subspace span{c, g 1 , g 2 , . . . , g D } thus, giving the exact reducedrank Wiener filter, while an AVF algorithm with D auxiliary vectors performs D successive line-search optimizations in span{w i−1 , g i }. Initialization: Algorithm 7: Minimum variance AVF algorithms.
The complexity of the four auxiliary vector algorithms (in block implementation) is given in Table 2 in multiplications per block of size T. For the algorithms compared in this table, the number of auxiliary vectors n AV = D − 1.

Recursive algorithms
Various recursive implementations can be derived from the algorithms introduced in the previous section. Note that each algorithm considered previously is initialized with the matched filter (w 0 = c). After D iterations, we obtain the reduced-rank Wiener filter w D 1 in D (R, c). We may initialize the algorithm with w D 1 and run another D iterations, thus giving the reduced-rank Wiener filter w D 2 in D (R, w D 1 ). Next D iterations will minimize the MSE in D (R, w D 2 ) and so on. Another kind of recursion can be used in adaptive sampleby-sample processing when the reduced-rank method at time instant n is initialized with the filter obtained at time (n − 1). An example of this technique is the "CG1" variant of the CGA described in [24].

Approximate sample-by-sample implementations
Sample-by-sample implementation of the MSWF (adaptive residual correlation or the adaptive ResCor algorithm [25]) can be derived as follows. Recall the expressions for the matched filters p i (Section 4.3.2), We may use a stochastic approximation to compute the expectation in (67), for example, as where 0 < γ < 1 is the forgetting factor. Note that this update

Backward recursion
Initialization: Algorithm 8: Summary of the adaptive ResCor algorithm.
does not require costly matrix vector multiplications. The same trick can be applied to the expectation which appears in backward recursion of the MSWF Algorithm 9: Summary of adaptive CGRRF.  where ξ i (k) is updated as The resulting algorithm (A-ResCor) is summarized in Algorithm 8. Low-complexity version of the CGRRF is given in Algorithm 9. Similarly, approximate adaptive implementations of the AVF can be obtained [12,26]. In Table 3, the complexity of sample-by-sample algorithms is given in multiplications per sample (the complexity of the LMS algorithm is given for comparison).

COMPUTER SIMULATIONS
The considered techniques were tested in a scenario which models joint (multiuser) detection in time division CDMA (TD-CDMA), universal mobile telecommunications system (UMTS), terrestrial radio access (UTRA), time division duplex (TDD) mode [27]. The parameters of a simulated DS-CDMA system are summarized in Table 4. Each of M = 4 system users transmits an i.i.d. sequence of QPSK-modulated symbols. The vehicular environment (Type A) with a mobile speed of 120 km/h was modeled [28,29]. At this speed, the coherence time of the channel is about 500 symbols. The channel is modeled as a tapped delay line with a fixed delay. Each tap coefficient is generated using Jakes fading model [30]. The delays and relative powers of the taps are given in Table 5. The channel noise is modeled as being white.
The receiver employs two diversity antennas, that is, the signals impinging on the antennas are assumed to be uncorrelated. This is reflected in the simulation setup by modelling the path from the transmitter to each receiver antenna by an independently fading, point-to-point channel model. Let r 1 (t) and r 2 (t) denote the respective antenna element outputs. After the matched filtering and chip-rate sampling, the combined received signal vector r(k) = r T 1 (k) r T 2 (k) T is described by (1) with the observation size N = 2S f = 32. Note that due to the users' motion, the channel matrix H is time-varying. The channel signature h 1 of the desired user (first column of H in (1)) is supposed to be known at the receiver. Exact and SMI algorithms, compared here, use (5) to estimate the covariance matrix of the received signal with the forgetting factor γ = 0.995. This value of γ is also used in approximate adaptive implementations. Initially, R(0) = δI where the diagonal loading parameter δ ensures that R(k) is inversible for small k.
In Figures 3 and 4, the convergence and tracking performance of the exact reduced-rank methods (ERRM, which include the POR, MSWF, and CGRRF) and of the conventional SMI technique is explored. In these figures, the bit error rate (BER) of the first user, averaged over the independent realizations of noise, transmitted symbols, fading waveforms, and random binary spreading codes is plotted versus the number of the received observations (k). The BER attained by the RAKE receiver is given for comparison. For Figure 4, the first user has the SNR of 9 dB, while the interfering users have the SNR of 21 dB, and for Figure 3, everybody is at 16 dB.
The curves of Figure 3 illustrate two distinct zones: the zone of initial convergence (first 150 symbols periods), followed by the zone of tracking. The first zone is characterized by minor accumulated changes in the propagation channel so that the performance depends mainly on the algorithm's convergence rate. It can be seen that the reduced-rank filters converge at least as rapidly as the full rank SMI filter.
In the second zone, the channel has been significantly driven off its initial value and tracking ability becomes dominant. We see that for moderate multiple-access interference (MAI) levels (Figure 3), tracking ability is dramatically improved by the rank reduction with lower filter ranks performing the best. For strong MAI (Figure 4), accurate modelling of the interference (which is better obtained with higher ranks) prevails over the tracking performance, and the full rank solution gives the lowest BER. The latter observation also applies to the next simulation example.
In Figures 5 and 6, we compare the steady-state performance of ERRM and SMI. As a performance measure, we used the BER of the desired user over 10 6 symbol periods and averaged over random binary spreading codes. In Figure 5, this BER is plotted versus the SNR (common to all system users), while in Figure 6, only SNR of the interfering users (Interference-to-Noise Ratio, or INR) changes with the SNR of the desired user fixed at 9 dB.
It can be concluded that for low SNR (INR), ERRM of low ranks (and, in particular, the RAKE which corresponds to D = 1) perform better. In this case, white noise dominates the MAI and the gain achieved through multiuser detection is small compared to the loss resulting from the misadjustment noise. The situation changes for high SNR (INR) when more filter degrees of freedom are required in order to cope with the strong MAI. It follows from the presented figures that although reduced-rank filters are not near-far resistant, they can provide significant gains over the SMI and RAKE receivers over wide ranges of MAI levels. computed as where w(k) is the filter estimate and R I+N is the interference plus noise covariance matrix. In this experiment, the desired user is at 9 dB SNR, while the interferers are at 21 dB. The performance attained by the full rank and RAKE filters is also given for reference. Interestingly enough, for the given simulation setup algorithms which use nonorthogonal sets of basis (auxiliary) vectors (A-CGRRF, A-MVAVF) result in better performance. Figures 5 and 6 suggest that adaptive rank selection (as a function of MAI level) can be used in order to improve the near-far resistance of reduced-rank filters. This solution is demonstrated in our next experiment, where the rank of CGRRF is adapted in order to maintain constant filter output SINR.
Specifically, each 250 symbol periods the instantaneous SINR at the output of the Dth stage of CGRRF is estimated according to (71) with R I+N replaced by its estimate R I+N (k) = R(k) − 1 h 1 (k)h H 1 (k). 5 The instantaneous SINR so obtained is time-averaged using the forgetting factor of 0.95. The value of rank (D) is then either increased or decreased by 1 in order to keep time-averaged SINR within the range 3 ± 1 dB.
In Figures 8 and 9, one realization of the time-averaged SINR and of the rank D versus time is shown. The propagation channel was fixed in this experiment in order to simplify the analysis. Initially, there are 16 users in the system with the desired user having the SNR of 8 dB and interferers at 14 dB. Starting rank value is 2. Over the first 1000 samples, the rank can be seen to converge and then stabilize at D = 5. At time k = 4000, six interfering users quit the system and D decreases to 3. Finally, at k = 6000, three of interferers reenter the communication and the filter rank again grows to 4. Figure 8 shows that the reception quality is kept reasonably well within the required limits. Another advantage of this (or similar) adaptive rank selection techniques is the efficient utilization of the available processing power. Filter rank (D) Figure 9: Rank D versus time for the CGRRF with adaptive rank selection.

CONCLUSIONS
A family of recently introduced least-square adaptive filtering techniques has been studied. All of these algorithms are shown to project and filter the received observation in the same low-dimensional subspace. This (Krylov) subspace is generated by taking the powers of the received covariance matrix on the cross-correlation (steering) vector. Consequently, considered techniques are related to Krylov subspace methods for linear systems and can be divided in four distinct groups: exact, auxiliary-vector, recursive, and approximate sample-by-sample methods. Numerical studies and complexity figures compare favorably the presented methods to conventional SMI and matched filtering techniques.