The fast recursive row-Householder subspace tracking algorithm

ARTICLE IN PRESS Signal Processing 89 (2009) 2514–2528

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

The fast recursive row-Householder subspace tracking algorithm Peter Strobach ¨hrnbach, Germany AST-Consulting Inc., Bahnsteig 6, 94133 Ro

a r t i c l e i n f o

abstract

Article history: Received 17 January 2009 Received in revised form 14 April 2009 Accepted 14 April 2009 Available online 23 April 2009

We introduce a new sequential algorithm for tracking the principal subspace and, optionally, the r dominant eigenvalues and associated eigenvectors of an exponentially updated covariance matrix of dimension N  N, where N4r. The method is based on an updated orthonormal-square (QS) decomposition using the row-Householder reduction. This new subspace tracker reaches a dominant complexity of only 3Nr multiplications per time update for tracking the principal subspace, which is the lower bound in dominant complexity for an algorithm of this kind. The new method is completely reflection based. An updating of inverse matrices is not used. & 2009 Elsevier B.V. All rights reserved.

Keywords: Subspace tracking Orthogonal iteration Eigenvalues Eigenvectors QS-factorization Householder PAST FAPI LORAF BiSVD

1. Introduction Subspace tracking is an important pre-processing step in modern adaptive systems. The goal is the recursive estimation of the principal subspace or even the r dominant or leading eigenvalues and the associated eigenvectors of a time-recursively updated covariance matrix UðtÞ of dimension N  N,

UðtÞ ¼ aUðt  1Þ þ zðtÞzT ðtÞ,

(1)

where zðtÞ is a data snapshot vector of dimension N4r and a is a positive exponential forgetting factor close to one. Subspace trackers are usually principal subspace trackers. More efforts are necessary to track the dominant eigenvalues and eigenvectors directly. Estimates of dominant eigenvectors can also be obtained by rotation of the orthonormal basis of the tracked principal subspace. E-mail address: [email protected] URL: http://www.ast-consulting.net 0165-1684/$ - see front matter & 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.04.012

A review of historical methods in subspace tracking can be found in [1]. In the early 1990s, it turned out that the principal subspace tracking problem can be solved with a dominant complexity of only OðNrÞ computations per time update. These algorithms were called fast subspace trackers. In this article, we are dealing with algorithms of this class only. One of the earlier representatives of these algorithms is the PAST (Projection Approximation Subspace Tracker), introduced by Yang in 1995 [2]. This algorithm has a characteristic structure that is similar to old recursive least squares (RLS) concepts of adaptive filtering [3]. It is based on a heavy use of the matrix inversion lemma [3] and is constructed around a time-updating of an inverse power matrix. This has the clear disadvantage that this updated inverse power matrix will ‘‘explode’’ in the case of fading signals, a case often encountered in realistic sequential online applications. The situation is literally the same as with classical RLS adaptive filters. Nevertheless, researchers in this field have silently ignored this obvious deficiency of the PAST concept and

ARTICLE IN PRESS P. Strobach / Signal Processing 89 (2009) 2514–2528

concentrated on the remedy of another problem encountered with PAST, namely the orthonormality errors that occur in typical PAST subspace estimates. A number of refined algorithms appeared, such as OPAST (Orthonormal PAST) [4], NIC (Novel Information Criterion) [5], and NP (Natural Power) [6] subspace trackers. The final state of development along this baseline of updating an inverse power matrix has recently been marked with the advent of the so-called Fast Approximated Power Iteration (FAPI) subspace tracker [7]. Our theoretical analysis reveals that this method appears optimal in terms of the fast recursive orthogonal iteration principle with linear subspace propagation model, because a ‘‘projection approximation’’ (like in PAST) is not employed in FAPI anymore. FAPI is certainly optimal in terms of the operations count, because it reaches the 3Nr lower bound of dominant complexity per time update. Nevertheless, the method still has the inherent drawback of the inverse matrix update, because it is basically an algorithm of the PASTtype, based on a heavy use of the matrix inversion lemma. From the standpoint of numerical mathematics and optimization, it is not difficult to see that a much more favorable solution to the subspace tracking problem exists in terms of appropriate orthonormal decompositions and their associated updating algorithms. An early approach of this kind is available in terms of the LORAF (Low-Rank Adaptive Filter) subspace trackers [8]. These algorithms are based on a recursively updated QR-factorization. They produce estimates of the dominant eigenvectors and corresponding eigenvalues directly, hence they could be named ‘‘eigensubspace’’ trackers as opposed to the conventional ‘‘principal subspace’’ trackers of the fast category. The LORAF3 algorithm was probably the first fast OðNrÞ eigensubspace tracker. The updating structure of LORAF3 is of the kind ‘‘QR ¼ QR þ rank one’’, a problem that has been well understood since the mid 1970s [9]. The method is merely a game of orthonormal matrices, hence the obtained eigensubspace estimates are clearly orthonormal up to machine accuracy. Moreover, no inverse matrices are updated anymore. The method has been extended and applied in several other problem areas. For instance SVD subspace tracking [10], and SVD subspace tracking using a sliding window in the serialized data case [12]. The application to adaptive ESPRIT frequency estimation has been described in [13]. Instrumental variable concepts of this kind were described in [14]. The tracking of generalized eigenvalues and eigenvectors using LORAF-type algorithms was the subject of [15]. The method generalizes to the multiple invariance case. Multiple invariance subspace trackers were described in [16], while [17–19] deal with SVD-based ESPRIT and adaptive multidimensional source localization problems. The SVD subspace trackers of [11], on the other hand, are principal subspace trackers, because they do not produce eigenvector or singular vector estimates directly. Recent investigations have shown that the LORAF concept is redundant in complexity. The simple reason is that QR is not the appropriate decomposition if only the principal subspace is sought. Too many computations are invested in maintaining the triangular structure of the R-matrix, although the particular triangular shape is

2515

nowhere needed in the algorithm as long as a direct tracking of the eigenvectors is not an issue. For tracking the principal subspace, the only requirement is that the decomposition produces an orthonormal Q -matrix. Hence we can leave the associated R-matrix completely undetermined in shape. Consequently, the complete concept can be formulated equivalently on a more general orthonormal-square (QS) decomposition, where S is simply an r  r matrix without any special shape constraints. This QS-decomposition can be updated very economically using a row-Householder reduction. This key idea presents a significant breakthrough in LORAFtype algorithms and leads directly to the new fast rowHouseholder subspace tracking algorithm, a method that is now more than four times as fast as the original LORAF algorithms. This paper is organized as follows: in Section 2, we formulate the basic theoretical background of fast recursive orthogonal iteration incorporating a linear subspace propagation model. This is the baseline of all fast subspace trackers discussed in this paper. In Section 3, we discuss the family of PAST-type subspace trackers. In Section 4, we review the concept behind the LORAF-type subspace trackers and show that this concept is independent of the underlying type of orthogonal decomposition. Then we show that the new QS-decomposition with rowHouseholder reduction gives us a striking advantage in the complexity of the resulting algorithm. Section 5 introduces the background that is needed for adaptive low-rank approximation, eigenvalue and eigenvector tracking using the row-Householder subspace tracker. Some deeper theoretical insights are gained from these considerations. In Section 6, we discuss some instructive experimental results. Section 7 concludes this paper.

2. The recursive orthogonal iteration principle All subspace trackers of this paper are based on orthogonal iteration. The orthogonal iteration principle for computing eigenvectors and eigenvalues of matrices is now well known [20,21]. One of the first obvious applications of this principle in subspace tracking traces back to the Owsley algorithm [22], introduced in 1978. This algorithm simply computes the time update of the covariance matrix according to (1), and then applies a single orthogonal iteration in each time step as follows: AðtÞ ¼ UðtÞQ ðt  1Þ, AðtÞ ¼ Q ðtÞSðtÞ : orthonormal factorization,

(2a) (2b)

where AðtÞ is an auxiliary matrix of dimension N  r. In the original paper [22], Owsley used a QR-factorization in (2b). In this case, Q ðtÞ is the matrix of estimated dominant eigenvectors, and SðtÞ is an upper-right triangular r  r-matrix with the estimated eigenvalues appearing in descending order of magnitude on the main diagonal of SðtÞ. If, however, only a basis of the principal subspace is sought, then we can leave the shape of the square r  r S-matrix completely undetermined, i.e., it can be a full square matrix. Then (2b) reduces to a

ARTICLE IN PRESS 2516

P. Strobach / Signal Processing 89 (2009) 2514–2528

QS-decomposition, which can be updated much more efficiently than the more restrictive QR-decomposition. For purposes of fast tracking, one further defines hðtÞ ¼ Q T ðt  1ÞzðtÞ.

(3)

This step is common to all fast subspace trackers. Here Q T ðt  1Þ serves as a data compressor that maps the information of interest into a subspace of usually much smaller dimension. We substitute (1) into (2a) to obtain T

AðtÞ ¼ aUðt  1ÞQ ðt  1Þ þ zðtÞh ðtÞ.

(4)

The recursion is not yet complete without a reasonable subspace propagation model. In [8], we introduced the following state space form: Q ðtÞ ¼ Q ðt  1ÞHðtÞ þ DðtÞ.

complexity to OðNrÞ complexity, but the tracked basis vectors in Q ðtÞ will no longer be the tracked eigenvectors. They become just an arbitrarily rotated variant of the estimated eigenvector set. The corresponding subspace trackers are more precisely classified as principal subspace trackers since they produce no direct estimates of the eigenvectors. This is not a disadvantage, because in many applications, the estimated eigenvectors are not needed directly. Moreover, there exists an easy to compute orthonormal r  r rotor that transforms Q ðtÞ back into the set of estimated eigenvectors, if desired. This is treated in detail in Section 5. 3. PAST-type subspace trackers

(5)

Clearly, the subspace transition matrix HðtÞ plays the role of a state transition matrix and DðtÞ is an innovations matrix that satisfies Q T ðt  1ÞDðtÞ ¼ 0.

(6)

Hence we can also write

HðtÞ ¼ Q T ðt  1ÞQ ðtÞ.

(7)

This model pays a key role in all further considerations. We substitute a one time step delayed version of (5) into (4). This yields the following recurrence for the auxiliary A-matrix:

But before we arrive at this point, it can be quite instructive to have a closer look at the concepts behind PAST-type subspace trackers. In this section, we review the classical PAST algorithm of [2] and discuss the approximations behind this algorithm. Then we show that the novel FAPI algorithm of [7] can be viewed as a PAST-algorithm without any ‘‘projection approximation’’ at all. All PAST-type subspace trackers are principal subspace trackers, i.e., they do not track the eigenvectors directly. A pre-multiplication of (10) with Q T ðt  1Þ yields the following recurrence, which can be considered as the basis of all PAST-type subspace trackers:

T

AðtÞ ¼ aAðt  1ÞHðt  1Þ þ aUðt  1ÞDðt  1Þ þ zðtÞh ðtÞ. (8) This recurrence is not yet of much practical value, because of the complexity of Uðt  1ÞDðt  1Þ. As a matter of fact, this term has commonly been neglected in all approaches of fast orthogonal iteration subspace tracking that have become known in the subspace tracking literature until now. It will be convenient, for the moment, to neglect this term as well. But we return on this subject and provide a deeper theoretical discussion in Section 5. This theoretical study will conclude with an amazing result. For the moment, we neglect the term and continue with a simplified OðNr 2 Þ recurrence for the auxiliary A-matrix, as usual:

HðtÞSðtÞ ¼ aSðt  1ÞHðt  1Þ þ hðtÞhT ðtÞ.

(11)

3.1. A review of the classical past algorithm The classical PAST algorithm has been originally derived as the solution to a curious optimization problem [2]. However, it presents no difficulty to show that PAST is nothing but a simplified form of an orthogonal iteration subspace tracker. In the classical PAST algorithm, the following two simplifications are made:

Hðt  1Þ:¼I : a posteriori projection approximation

(12)

and

T

AðtÞ ¼ aAðt  1ÞHðt  1Þ þ zðtÞh ðtÞ.

(9)

In [8], this was the basic recurrence of subspace tracker LORAF2. The problem comes with the orthonormal factorization of AðtÞ as seen in (2b). A batch orthonormal factorization cannot be carried out with less than OðNr 2 Þ operations. Hence it is clear that the goal should be a direct updating of the orthonormal factorization of the A-matrix according to T

Q ðtÞSðtÞ ¼ aQ ðt  1ÞSðt  1ÞHðt  1Þ þ zðtÞh ðtÞ.

(10)

For the special case of a QR-factorization, it requires again OðNr 2 Þ operations to update this decomposition, because the Y-matrix appears as an unpleasant perturbation that destroys the ‘‘old’’ triangular matrix. The solution to this gordian knot problem of subspace tracking lies in the application of a full square r  r S-matrix in the decomposition. This will allow us to move from OðNr 2 Þ

HðtÞ:¼I : a priori projection approximation.

(13)

This yields the basic recurrence for the S-matrix, that is underlying the PAST algorithm: T

SðtÞ ¼ aSðt  1Þ þ hðtÞh ðtÞ.

(14)

In the same fashion, simplify (10) to obtain T

Q ðtÞSðtÞ ¼ aQ ðt  1ÞSðt  1Þ þ zðtÞh ðtÞ.

(15)

Clearly, a recurrence for the orthonormal Q -matrix is now obtained by a post-multiplication of (15) with S1 ðtÞ: T

Q ðtÞ ¼ aQ ðt  1ÞSðt  1ÞS1 ðtÞ þ zðtÞh ðtÞS1 ðtÞ.

(16)

Obviously, PAST will be a game of inverse S-matrices. Define PðtÞ ¼ S1 ðtÞ,

(17)

ARTICLE IN PRESS P. Strobach / Signal Processing 89 (2009) 2514–2528

to see that T

PðtÞ ¼ ½aSðt  1Þ þ hðtÞh ðtÞ1 .

(18)

This is a case for the classical matrix inversion lemma [3]: ðD þ BCÞ1 ¼ D1  D1 BðI þ CD1 BÞ1 CD1 .

(19)

Set D ¼ aSðt  1Þ. Then D1 ¼ ð1=aÞPðt  1Þ. Moreover, set T B ¼ hðtÞ and C ¼ h ðtÞ to obtain

2517

consequences resulting from an application of the a priori projection approximation (13), because this simplification violates the underlying subspace propagation model (5) for the actual Q -estimate. The result is a significant loss of orthonormality of the subspace basis matrix Q ðtÞ. 3.2. A look at the FAPI-algorithm

We can see that S and P are both perfectly symmetric. Substitute (20) into (16):

An exact form of a PAST-type subspace tracker is meanwhile available in terms of the FAPI algorithm [7]. This algorithm can be derived along the same baseline as shown in the case of PAST, but the steps will be considerably more involved, because no Y-matrix is simplified anymore. Hence the starting point in the FAPI concept is the unmodified recurrence (11). We should have a look at this concept. It is based on two basic recurrences and an important insight. First of all, rewrite (11) in transposed form as follows:

1 Q ðtÞ ¼ aQ ðt  1ÞSðt  1Þ ½Pðt  1Þ  fðtÞgT ðtÞ

ðHðtÞSðtÞÞT ¼ aðSðt  1ÞHðt  1ÞÞT þ hðtÞh ðtÞ.

PðtÞ ¼

1

a

½Pðt  1Þ  fðtÞgT ðtÞ,

(20)

where gðtÞ ¼ Pðt  1ÞhðtÞ, 1 fðtÞ ¼ gðtÞ. T a þ h ðtÞgðtÞ

(21a) (21b)

T

a

1 T þ zðtÞh ðtÞ ½Pðt  1Þ  fðtÞgT ðtÞ.

(22)

a

Clearly, Sðt  1ÞPðt  1Þ ¼ I. Observe that as a consequence of (21b), fðtÞ and gðtÞ are linearly dependent. Therefore, we T can write: fðtÞgT ðtÞ ¼ gðtÞf ðtÞ. Moreover, (21a) can also be written as Sðt  1ÞgðtÞ ¼ hðtÞ. Finally, (21b) can be rewritT ten in the form gðtÞ ¼ ða þ h ðtÞgðtÞÞfðtÞ. Substitute all these into (22) and verify that T

T

Q ðtÞ ¼ Q ðt  1Þ  Q ðt  1ÞhðtÞf ðtÞ þ zðtÞf ðtÞ.

(23)

Introduce the complement of the orthogonal projection of zðtÞ onto the subspace spanned by Q ðt  1Þ, z? ðtÞ ¼ zðtÞ  Q ðt  1ÞhðtÞ,

(24)

to obtain the final time update of the Q -matrix according to PAST as follows: T

Q ðtÞ ¼ Q ðt  1Þ þ z? ðtÞf ðtÞ.

(25)

Equation-array (26a)–(26f) is a quasicode listing of this algorithm that is exactly the PAST-algorithm of [2]. Initialize this algorithm with Q ð0Þ:¼random orthonormal and Pð0Þ:¼s1 I, where s is a very small positive constant to account for the singularity of Sð0Þ. We see, this looks very much like classical recursive least squares: FOR t ¼ 1; 2; 3; . . . hðtÞ ¼ Q T ðt  1ÞzðtÞ

(26a)

gðtÞ ¼ Pðt  1ÞhðtÞ 1 fðtÞ ¼ gðtÞ a þ hT ðtÞgðtÞ 1 PðtÞ ¼ ½Pðt  1Þ  fðtÞgT ðtÞ

(26b)

(26d)

z? ðtÞ ¼ zðtÞ  Q ðt  1ÞhðtÞ

(26e)

a

T

Q ðtÞ ¼ Q ðt  1Þ þ z? ðtÞf ðtÞ

(26c)

(26f)

ENDFOR. Experience has shown that the a posteriori projection approximation (12) has only a minor, or even negligible effect on the practical performance of an orthogonal iteration subspace tracker. Much more severe are the

(27)

To aid an algorithm of the PAST-type, introduce a generalized inverse power matrix as follows: PðtÞ ¼ ðSðtÞHðtÞÞT .

(28)

Then T

ðHðtÞSðtÞÞT ¼ ½aP1 ðt  1Þ þ hðtÞh ðtÞ1 .

(29)

Apply the matrix inversion lemma (19) in the form T D ¼ aP1 ðt  1Þ, B ¼ hðtÞ and C ¼ h ðtÞ to find ðHðtÞSðtÞÞT ¼

1

a

T

½I  fðtÞh ðtÞ Pðt  1Þ,

(30)

where gðtÞ ¼ Pðt  1ÞhðtÞ, 1 fðtÞ ¼ gðtÞ. a þ hT ðtÞgðtÞ

(31a) (31b)

The problem is the reverse ordering of the matrix pair SðtÞHðtÞ in subsequent time steps. However, this problem can be relaxed, if the following relationship is used: PðtÞ ¼ ST ðtÞHT ðtÞ ¼ HT ðtÞðHðtÞSðtÞÞT HT ðtÞ.

(32)

Consequently, (30) can be transformed into the desired recurrence for the auxiliary inverse PðtÞ as follows: PðtÞ ¼

1

a

HT ðtÞ½I  fðtÞhT ðtÞ Pðt  1ÞHT ðtÞ.

(33)

This constitutes the first basic recurrence in the FAPI algorithm. Incorporating a posteriori (12) and a priori (13) projection approximations, one can easily see that (33) becomes identical to the PAST recurrence for the inverse power matrix (26d). Next we are searching for a recurrence that enables an exact updating of the Q -matrix in time. Recall again (10) and postmultiply by S1 ðtÞ to obtain T

Q ðtÞ ¼ aQ ðt  1ÞSðt  1ÞHðt  1ÞS1 ðtÞ þ zðtÞh ðtÞS1 ðtÞ. (34) We can verify the similarity to (16) in the case of PAST. Unlike in PAST, Hðt  1Þ complicates the expression.

ARTICLE IN PRESS 2518

P. Strobach / Signal Processing 89 (2009) 2514–2528

quadratic equation:

Rearrange (11) in the following form:

aSðt  1ÞHðt  1ÞS1 ðtÞ ¼ HðtÞ  hðtÞhT ðtÞS1 ðtÞ.

(35)

Now (34) can be rewritten as follows: T

T

Q ðtÞ ¼ Q ðt  1Þ½HðtÞ  hðtÞh ðtÞS1 ðtÞ þ zðtÞh ðtÞS1 ðtÞ T

¼ Q ðt  1ÞHðtÞ þ z? ðtÞh ðtÞS1 ðtÞ.

(36)

This is already the desired recurrence for the Q -matrix. A quick inspection reveals that this update obeys perfectly the desired exact state space structure of the subspace propagation model (5) with innovations matrix

DðtÞ ¼ z? ðtÞhT ðtÞS1 ðtÞ,

(37)

while the corresponding simplified update (26f) in the classical PAST algorithm did not satisfy the exact subspace propagation model. Moreover, one can even show that (36) can be posed in the following form: T

(38)

Q ðtÞ ¼ ½Q ðt  1Þ þ z? ðtÞf ðtÞHðtÞ.

Clearly, (36) and (38) are equivalent formulations of the same update, if T

T

f ðtÞHðtÞSðtÞ ¼ h ðtÞ.

(39)

To verify this, simply premultiply basic equation (11) by T f ðtÞ and use (31a), (31b) and (28) to see that 1

T

f ðtÞHðtÞSðtÞ ¼

T

a þ h ðtÞgðtÞ

gT ðtÞ½aSðt  1ÞHðt  1Þ

T

¼

þ hðtÞh ðtÞ 1

a þ hT ðtÞgðtÞ

T

þ gT ðtÞhðtÞh ðtÞ 1 T ½a þ gT ðtÞhðtÞh ðtÞ ¼ a þ hT ðtÞgðtÞ T

T

Q ðtÞQ ðtÞ ¼ I ¼ H ðtÞ½Q ðt  1Þ þ

T

¼ HT ðtÞ½I þ fðtÞzT? ðtÞz? ðtÞf ðtÞHðtÞ

(42)

Rearrange (41) and use the matrix inversion lemma (19) T with D ¼ I, B ¼ ZðtÞfðtÞ, and C ¼ f ðtÞ to obtain T

T

1

HðtÞH ðtÞ ¼ ½I þ ZðtÞfðtÞf ðtÞ ¼I

ZðtÞ T

1 þ ZðtÞf ðtÞfðtÞ

T

fðtÞf ðtÞ.

(43)

The key insight is that this equation holds for a HðtÞ that is of the form ‘‘identity—symmetric rank one’’:

HðtÞ ¼ I  rfðtÞf T ðtÞ.

LORAF-type subspace trackers are a straight implementation of (10). The goal is the representation of (10) in terms of an orthonormal decomposition, which is sequentially updated in time. Hereby, we are working directly on the S-matrices, not on their inverses. The updated S-matrix is a strongly diagonal dominant matrix, but is still a power matrix and not yet an estimated dominant eigenvalue matrix. Generally, we conclude that LORAFtype subspace trackers are S-domain techniques, while PAST-type subspace trackers are S1 -domain techniques. We know from adaptive filtering [23,24], that S-domain techniques are generally preferable over S1-domain techniques, if the complexity is comparable. This has motivated the development of the fast recursive rowHouseholder subspace tracker, an S-domain technique with a dominant complexity of 3Nr operations per time step. This algorithm will be presented in Section 4.3.

(46)

where (47)

Q ðtÞSðtÞ ¼ aQ ðt  1ÞSðt  1ÞHðt  1Þ (41)

where ZðtÞ ¼ zT? ðtÞz? ðtÞ.

4. LORAF-type subspace trackers

is a unit-norm variant of z? ðtÞ. Substitute (46) into (10) to obtain

½Q ðt  1Þ þ z? ðtÞf ðtÞHðtÞ T

and the FAPI algorithm is essentially complete. The special structure of HðtÞ according to (44) ensures the desired OðNrÞ complexity of the algorithm. A variant of this algorithm with minimal dominant complexity 3Nr operations per time step is summarized in Table III of [7].

z? ðtÞ ¼ Z 1=2 ðtÞz? ðtÞ

fðtÞzT? ðtÞ T

¼ HT ðtÞ½I þ ZðtÞfðtÞf ðtÞHðtÞ,

(45)

zðtÞ ¼ Z 1=2 ðtÞz? ðtÞ þ Q ðt  1ÞhðtÞ, (40)

T

¼ 0,

Begin with the orthonormal decomposition of the snapshot vector zðtÞ as follows:

The purpose of (38) is that we can now write T

T

1 þ ZðtÞf ðtÞfðtÞ

4.1. General structure of S-domain techniques

½agT ðtÞPT ðt  1Þ

¼ h ðtÞ.

ZðtÞ

T

f ðtÞfðtÞr2  2r þ

(44)

Hence it is not difficult to verify that the parameter r is given as one of the two solutions to the following

T

þ ½Z 1=2 ðtÞz? ðtÞ þ Q ðt  1ÞhðtÞh ðtÞ ¼ ½Q ðt  1Þ z? ðtÞ 2 3 aSðt  1ÞHðt  1Þ þ hðtÞhT ðtÞ 4 5.  T Z 1=2 ðtÞh ðtÞ

(48)

Determine an orthonormal rotor or reflector matrix GðtÞ of dimension ðr þ 1Þ  ðr þ 1Þ with property GT ðtÞGðtÞ ¼ I so that the old augmented and updated S-matrix at time t  1 is transformed into a new S-matrix at time t with appended zero row vector: " #   aSðt  1ÞHðt  1Þ þ hðtÞhT ðtÞ SðtÞ . (49) ¼ GðtÞ T 00 Z 1=2 ðtÞh ðtÞ Consequently, the Q -matrix is updated as follows: ½Q ðtÞ zq ðtÞ ¼ ½Q ðt  1Þ z? ðtÞGT ðtÞ,

(50)

ARTICLE IN PRESS P. Strobach / Signal Processing 89 (2009) 2514–2528

An updating structure of this kind was first seriously studied in the area of subspace tracking and low-rank adaptive filtering by Strobach in [8]. In these classical LORAF algorithms, GðtÞ was a sequence of Givens rotations [21] and SðtÞ was structurally constrained to an upperright triangular matrix. These algorithms track the eigenvectors and eigenvalues directly, just like Owsley’s classical algorithm. It is not difficult to see that in a concept of this kind, a significant amount of computations must be invested in the maintenance of the triangular shape of the S-matrix during recursion (49). This carries over to a significant complication of the orthonormal basis update (50), which is in fact an estimated eigenvector update here. However, the triangular structure of the S-matrix is nowhere needed in the algorithm, as long as the issue is principal subspace tracking only. All that is needed in the principal subspace tracking case is the zero bottom-row vector in the updated S-matrix on the left side of (49). Hence, a minimum complexity algorithm for principal subspace tracking should not expend more computations than those needed to construct this bottom row of zero elements.

4.2. The row-Householder reduction The optimal solution to meet these requirements is the row-Householder reduction. An excellent reference in the area of row-Householder transformations and their application is [25]. The goal in the row-Householder reduction is the application of a Householder matrix from the left to annihilate a row vector in a matrix of dimension ðr þ 1Þ  r. Hence the prototype problem is the design of a Householder reflector H, T

H ¼ I  2vv , where " # v , v¼

j

(51)

   Y X ¼H T , 00 x

(53)

where both X and Y are square r  r matrices. Clearly, the unique solution to this problem is a H with a bottom-row vector that lies in the null space of the appended-X matrix in (53). Hence the conditional equation that determines v is given by the bottom row of (53): " # X T ½0    0 ¼ ð½0    01  2j½v jÞ xT ¼ xT  2jvT X  2j2 xT .

(54)

Rearranging terms in (54) yields a conditional equation for a scaled variant b of v as follows: XT b ¼ x!b,

b¼

2j v. 1  2j2

(55)

(56)

A final point is the determination of j. This parameter is determined by the unit-norm constraint of v, or equivalently vT v þ j2 ¼ 1.

(57)

Note that from (56), we can write vT v ¼

ð1  2j2 Þ2 T b b. 4j2

(58)

Eliminate vT v by combination of (57) and (58): v T v þ j2 ¼

ð1  2j2 Þ2 T b b þ j2 ¼ 1, 4j2

(59)

or equivalently T

ð1  2j2 Þ2 b b þ 4j4 ¼ 4j2 . This conditional equation has the two solutions: ! 1 1 j21=2 ¼ 1  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . T 2 b bþ1

(60)

(61)

In batch applications, additional considerations are made in order to determine which one of these two solutions should be preferred from a standpoint of numerical precision. In sequential processing, there will be no time for such considerations. We always pursue the solution with the positive sign. This yields the following overall procedure for determining the parameters v and j of the desired row-Householder reflector: XT b ¼ x!b

(62a)

b ¼ 4ðbT b þ 1Þ

(62b)

1 2

with vT v ¼ 1 and HH ¼ I that annihilates the bottom-row vector xT in a problem of the kind 

where

1

j2 ¼ þ pffiffiffi (52)

2519

v¼

b

1  2j2 b 2j

(62c) (62d)

4.3. The fast recursive row-Householder subspace tracker It is now an easy task to derive the fast recursive rowHouseholder subspace tracker, a LORAF-type or S-domain technique that reaches the lower bound in dominant complexity for an algorithm of this kind. This method is based on the S and Q -updates as shown in (49) and (50), but now we are using the row-Householder reflector HðtÞ in place of the G-rotor: 2 3 " # aSðt  1ÞHðt  1Þ þ hðtÞhT ðtÞ SðtÞ 5, (63) ¼ HðtÞ4 T 00 Z 1=2 ðtÞh ðtÞ ½Q ðtÞ zq ðtÞ ¼ ½Q ðt  1Þ z? ðtÞHðtÞ.

(64)

Householder methods for updating the Q -matrix according to (64) are not new. In [26], a gradient-based algorithm was presented for adapting the H-reflector. Here, we concentrate on (63). Clearly, the H-reflector must be determined so as to annihilate the bottom row of the

ARTICLE IN PRESS 2520

P. Strobach / Signal Processing 89 (2009) 2514–2528

updated and appended S-matrix in (63). This is just a prototype problem of the kind (53). Hence, we set T

XðtÞ ¼ aSðt  1ÞHðt  1Þ þ hðtÞh ðtÞ,

(65a)

xðtÞ ¼ Z 1=2 ðtÞhðtÞ,

(65b)

and determine the reflector parameters according to (62a)–(62d). Then (63) can be rewritten as " # # " #   " vðtÞ XðtÞ XðtÞ SðtÞ T ¼ ½v  2 . (66) ðtÞ j ðtÞ jðtÞ xT ðtÞ xT ðtÞ 00

and (42), and calculate ZðtÞ ¼ zT? ðtÞz? ðtÞ T

¼ ½zT ðtÞ  h ðtÞQ T ðt  1Þ½zðtÞ  Q ðt  1ÞhðtÞ T

T

¼ z ðtÞzðtÞ  h ðtÞhðtÞ.

Finally, we should look at the Q -update (71a) and (71b) from the perspective of the subspace propagation model (5). Note that the following partitioned form of HðtÞ can be established: "

Denote " T

T

p ðtÞ ¼ ½v ðtÞ jðtÞ

XðtÞ xT ðtÞ

#

HðtÞ ¼ I  2 .

(67)

SðtÞ ¼ XðtÞ  2vðtÞpT ðtÞ, T

(68a) T

½0    0 ¼ x ðtÞ  2jðtÞp ðtÞ.

Z 1=2 ðtÞ pðtÞ ¼ hðtÞ. 2jðtÞ

½Q ðtÞ zq ðtÞ ¼ ½Q ðt  1Þ z? ðtÞ  2½Q ðt  1Þ z? ðtÞ " # vðtÞ  ½vT ðtÞ jðtÞ jðtÞ ¼ ½Q ðt  1Þ z? ðtÞ  2½Q ðt  1ÞvðtÞ þ jðtÞz? ðtÞ½vT ðtÞ jðtÞ. Extract the following relations: eðtÞ ¼ jðtÞz? ðtÞ þ Q ðt  1ÞvðtÞ,

(71a)

Q ðtÞ ¼ Q ðt  1Þ  2eðtÞvT ðtÞ.

(71b)

The appended column zq ðtÞ is uninteresting and is not computed explicitly. An important step in speeding up the algorithm is the elimination of the explicit calculation of z? ðtÞ. Recall that (46) can also be written as

Q ðt  1ÞhðtÞ þ Q ðt  1ÞvðtÞ ¼ jðtÞZ 1=2 ðtÞzðtÞ  Q ðt  1Þ (73)

This, in turn, yields

eðtÞ ¼ jðtÞZ

1=2

ðtÞzðtÞ  Q ðt  1ÞwðtÞ.

½vT ðtÞ jðtÞ

I  2vðtÞvT ðtÞ

2jðtÞvðtÞ

2jðtÞvT ðtÞ

1  2j2 ðtÞ #

dðtÞ

HðtÞ T

d ðtÞ

1  2j2 ðtÞ

#

.

(76)

Clearly, we satisfy the subspace propagation model (5) with

HðtÞ ¼ I  2vðtÞvT ðtÞ,

(77a)

dðtÞ ¼ 2jðtÞvðtÞ,

(77b)

DðtÞ ¼ z? ðtÞdT ðtÞ.

(78)

Observe the similarity of the forms obtained for the subspace state transition matrix HðtÞ in case of the FAPI algorithm (44), and in case of the fast row-Householder subspace tracker (77a). In the same fashion, compare the innovation matrices as obtained in the FAPI-case (37), and in the fast row-Householder case (78). We see that the key in these algorithms is the special structure of HðtÞ and the fact that the innovation represented by DðtÞ is only a rankone layer. As a consequence of the special structure of the Ymatrix, we can speed up the computation of the term Sðt  1ÞHðt  1Þ a little by noting that Sðt  1ÞHðt  1Þ ¼ Sðt  1Þ  2Sðt  1Þvðt  1ÞvT ðt  1Þ. (79) Hence

eðtÞ ¼ jðtÞZ 1=2 ðtÞzðtÞ  jðtÞZ 1=2 ðtÞ

wðtÞ ¼ jðtÞZ 1=2 ðtÞhðtÞ  vðtÞ,

jðtÞ

(72)

Now use (72) to eliminate z? ðtÞ in (71a). This yields

½jðtÞZ 1=2 ðtÞhðtÞ  vðtÞ.

#

and innovations matrix (70)

In the same fashion, evaluate (64) to obtain

z? ðtÞ ¼ Z 1=2 ðtÞzðtÞ  Z 1=2 ðtÞQ ðt  1ÞhðtÞ.

¼

vðtÞ

(69)

Substituting this into (68a) yields the final update of the S-matrix: Z 1=2 ðtÞ T vðtÞh ðtÞ. jðtÞ

"

(68b)

Using xðtÞ ¼ Z 1=2 ðtÞhðtÞ and pðtÞ ¼ ð1=2jðtÞÞxðtÞ, we find

SðtÞ ¼ XðtÞ 

" ¼

Then, we can write

(75)

uðt  1Þ ¼ Sðt  1Þvðt  1Þ,

(80a)

Sðt  1ÞHðt  1Þ ¼ Sðt  1Þ  2uðt  1ÞvT ðt  1Þ.

(80b)

Equation-array (81a)–(81n) is a quasicode listing of this fast recursive row-Householder subspace tracking algorithm. Initialize this algorithm with Q ð0Þ:¼ random orthonormal and Sð0Þ:¼sI, where s is a very small positive constant to account for the singularity of Sð0Þ. Additionally, set vð0Þ:¼½0    0T .

(74a) (74b)

Moreover, it remains to compute ZðtÞ without explicit quantification of z? ðtÞ. For this purpose, recall (24)

FOR t ¼ 1; 2; 3; . . . hðtÞ ¼ Q T ðt  1ÞzðtÞ

(81a)

T

(81b)

T

ZðtÞ ¼ z ðtÞzðtÞ  h ðtÞhðtÞ

ARTICLE IN PRESS P. Strobach / Signal Processing 89 (2009) 2514–2528

uðt  1Þ ¼ Sðt  1Þvðt  1Þ

(81c) T

T

XðtÞ ¼ aSðt  1Þ  2auðt  1Þv ðt  1Þ þ hðtÞh ðtÞ

(81d)

XT ðtÞbðtÞ ¼ Z 1=2 ðtÞhðtÞ!bðtÞ

(81e)

bðtÞ ¼ 4ðbT ðtÞbðtÞ þ 1Þ

(81f)

1 1 j ðtÞ ¼ þ pffiffiffiffiffiffiffiffi 2 bðtÞ

(81g)

2

gðtÞ ¼ dðtÞ ¼

1  2j2 ðtÞ 2jðtÞ jðtÞ

(81h) (81i)

Z 1=2 ðtÞ vðtÞ ¼ gðtÞbðtÞ 1 T SðtÞ ¼ XðtÞ  vðtÞh ðtÞ dðtÞ wðtÞ ¼ dðtÞhðtÞ  vðtÞ

(81j) (81k) (81l)

eðtÞ ¼ dðtÞzðtÞ  Q ðt  1ÞwðtÞ

(81m)

Q ðtÞ ¼ Q ðt  1Þ  2eðtÞvT ðtÞ

(81n)

ENDFOR. Under infinite precision conditions, the performance of this fast row-Householder subspace tracker would be identical to LORAF2 of [8]. Moreover, a fast implementation of LORAF3 can be obtained by setting vðt  1Þ ¼ ½0    0T throughout the algorithm. This corresponds to the application of the a posteriori projection approximation criterion (12). In this case, (81c) can be omitted and (81d) can be simplified accordingly. Besides the discussed projection approximations, other approximative algorithms of the recursive orthogonal iteration class appeared in the literature. For instance, the so-called ‘‘Data Projection Method’’ (DPM) and its fast implementation (FDPM) [27]. This algorithm produces orthonormal, but severely degraded subspace estimates. The approximation underlying this algorithm is a rough simplification of the S-matrix. Several experimental results shown in [27] are proven wrong. We are not recommending any further approximations beyond the discussed a posteriori projection approximation. The operations count of algorithm (81a)–(81n) is dominated by three distinct Nr computation layers:

 The input data compaction step according to (81a).  The extraction of the innovation in the actual data 

snapshot according to (81m). The updating of the Q -matrix according to (81n).

These three Nr layers constitute a set of principal computations. That is, these layers are decoupled. None of these three computation layers can be substituted by any combination of the other two computation layers in the set. Therefore, it is postulated that the operations count for a fast orthogonal iteration subspace tracker providing strictly orthonormal principal subspace basis estimates will be lower bounded by a dominant complexity of 3Nr operations per time update.

2521

4.4. The ‘‘inverse’’ row-Householder subspace tracker We should not conclude this section without a little algebraic game. The challenge is the conversion of the S-domain row-Householder subspace tracker of (81a)–(81n) into a perfectly equivalent ‘‘inverse’’ algorithm in the S1 -domain, effectively a row-Householder subspace tracker that operates on inverse power matrices, just like FAPI. To aid into the development of such an ‘‘inverse’’ rowHouseholder subspace tracker, we adopt the definition of the generalized inverse power matrix from the FAPImethod, as shown in (28). The first step in the development will be the removal of the small internal system of linear equations (81e). We recall (65a) and write T

XT ðtÞ ¼ aHðt  1ÞST ðt  1Þ þ hðtÞh ðtÞ T

¼ aP1 ðt  1Þ þ hðtÞh ðtÞ.

(82)

Clearly, now, (81e) can be solved for the desired b-vector as follows: bðtÞ ¼ Z 1=2 ðtÞXT ðtÞhðtÞ T

¼ Z 1=2 ðtÞ½aP1 ðt  1Þ þ hðtÞh ðtÞ1 hðtÞ.

(83)

This is again a case for the classical matrix inversion T lemma (19). Set D ¼ aP1 ðt  1Þ, B ¼ hðtÞ and C ¼ h ðtÞ. After some algebra, we obtain the striking result: bðtÞ ¼ Z 1=2 ðtÞfðtÞ,

(84)

where fðtÞ was introduced in (31b) in the context of the FAPI-algorithm. Interestingly, the b-vector in the fast row-Householder algorithm is perfectly aligned with the f -vector in the FAPI-algorithm. We can draw a red line of such instructive connections throughout the development of our ‘‘inverse’’ row-Householder algorithm. For instance, discover that according to (81j) and (84), we must have vðtÞ ¼ gðtÞbðtÞ ¼ gðtÞZ 1=2 ðtÞfðtÞ.

(85)

Apparently, vðtÞ and fðtÞ must be linearly dependent, because they are just scaled versions of each other in the construction of the Y-matrix according to (44) in the FAPI-case and (77a) in the row-Householder case. Introduce an auxiliary variable ZðtÞ as follows:

ZðtÞ ¼

1

gðtÞZ 1=2 ðtÞ

.

(86)

Then write fðtÞ ¼ ZðtÞvðtÞ.

(87)

A final step in the algorithm is the updating of the inverse power matrix using the Badeau–David–Richardrecursion (33), where the Y-matrix is symmetric by definition according to (77a), and the inverse is easily developed as 2 vðtÞvT ðtÞ 1  2vT ðtÞvðtÞ 2 ¼Iþ vðtÞvT ðtÞ 2j2 ðtÞ  1

H1 ðtÞ ¼ I þ

¼ I  GðtÞvðtÞvT ðtÞ,

(88)

ARTICLE IN PRESS 2522

P. Strobach / Signal Processing 89 (2009) 2514–2528

where we used vT ðtÞvðtÞ ¼ 1  2j2 ðtÞ and

GðtÞ ¼ 

2 1 . ¼ 2j2 ðtÞ  1 jðtÞgðtÞ

(89)



This yields the P-update according to (33) as follows, where we skipped the ðtÞ-index of the vectors and scalars temporarily, for convenience: PðtÞ ¼

1

a

T

½ðI  2vvT ÞðI  Zvh Þ Pðt  1ÞðI  GvvT Þ.

(90)

Note that this formula is highly redundant in complexity. An effective computation is now developed along the following line of thought: define an auxiliary vector a as follows: a ¼ 2v þ Zð2j2  1Þh.

(91)

Then, we can write PðtÞ ¼

1

a

½ðI  vaT Þ Pðt  1ÞðI  GvvT Þ.

(92)

As a minor miracle, it turns out that a ¼ 2w,

(93)

where w was given in (81l). Hence define the transformations



algorithm to obtain an efficient eigenvalue and eigenvector tracker. This will be the subject of the next section. In the inverse form (101a)–(101r) of the algorithm, exception handling becomes a bit trickier, because the updating operates on a generalized inverse. Notice that this algorithm is essentially a classical RLS adaptive filter and shares many properties with the algorithms of this class. One may argue that the Oðr 3 Þ complexity of the internal system of linear equation (81e) is a stumbling block for the S-domain variant of the algorithm. Note, however, that the recursivity that is inherent in this algorithm, is also present in the S-domain variant. It can be exploited, for instance, by QR-factorizations of the X- and S-matrices in the algorithm. It is easily seen that a time recursive updating of these factorizations amounts not more than a recurrent application of ‘‘QR ¼ QR þ rank one’’ updates. This is particularly easy to accomplish in the square r  r case, because in a problem of the kind Q ðtÞRðtÞ ¼ Q ðt  1ÞRðt  1Þ þ uðtÞvT ðtÞ, we just have to find a vector aðtÞ that satisfies

c ¼ Pðt  1Þv,

(94a)

uðtÞ ¼ Q ðt  1ÞaðtÞ.

q ¼ 2PT ðt  1Þw,

(94b)

Clearly aðtÞ ¼ Q T ðt  1ÞuðtÞ,

and introduce the filter y ¼ G½c þ ðqT vÞv,

(95)

PðtÞ ¼

1

a

½Pðt  1Þ þ vðtÞqT ðtÞ  yðtÞvT ðtÞ,

(96)

where we re-introduced the ðtÞ-index. Hereby, the ‘‘inverse’’ row-Householder subspace tracker is complete and is casted into equation-array (101a)–(101r). Like PAST (26a)–(26f), this algorithm is initialized with Q ð0Þ:¼random orthonormal and Pð0Þ:¼s1 I, where s is a very small positive constant to account for the singularity of Sð0Þ. The following remarks can be established:

 The algorithm (101a)–(101r) can be viewed as a

 

variant of an FAPI subspace tracker. While the original FAPI subspace tracker [7, Table III] operates on the f -vector, the ‘‘inverse’’ row-Householder algorithm of (101a)–(101r) operates on the v-vector of the underlying Householder reflector. The operations count of algorithm (101a)–(101r) is practically the same as for the FAPI-algorithm in Table III of [7]. A structural comparison with the S-domain method in (81a)–(81n) reveals that the S-domain method has the attraction that it produces directly the power estimates of the underlying natural modes in the signal on the main diagonal of the S-matrix. Hence, it is relatively easy to handle nasty exceptions, such as rank-drop and fading signals within this form of the algorithm. All exception processing will be concentrated on the solution of the small internal system of linear equations (81e). Moreover, it is relatively easy to extend this

(98)

(99)

and therefore Q ðtÞRðtÞ ¼ Q ðt  1Þ½Rðt  1Þ þ aðtÞvT ðtÞ

to obtain the final recursion as follows:

(97)

(100)

is a problem in which the QR-factorization can be easily restored using a pair of internal rotors or reflectors [28]. Another option exists in terms LUfactorizations and fast Bennett LU-factor updating as advocated in the quasicode of Fig. 2.1 in [29]. This way the complexity of the small internal system of linear equations can be brought down to Oðr 2 Þ without inverse matrix updating. FOR t ¼ 1; 2; 3; . . . hðtÞ ¼ Q T ðt  1ÞzðtÞ

(101a)

T

ZðtÞ ¼ z ðtÞzðtÞ  h ðtÞhðtÞ

(101b)

gðtÞ ¼ Pðt  1ÞhðtÞ 1 gðtÞ fðtÞ ¼ a þ hT ðtÞgðtÞ

(101c)

T

(101d)

bðtÞ ¼ Z 1=2 ðtÞfðtÞ

(101e)

bðtÞ ¼ 4ðbT ðtÞbðtÞ þ 1Þ

(101f)

1 1 j2 ðtÞ ¼ þ pffiffiffiffiffiffiffiffi 2 bðtÞ

(101g)

gðtÞ ¼ dðtÞ ¼

1  2j2 ðtÞ 2jðtÞ jðtÞ

Z 1=2 ðtÞ 1 GðtÞ ¼ jðtÞgðtÞ vðtÞ ¼ gðtÞbðtÞ

(101h) (101i) (101j) (101k)

ARTICLE IN PRESS P. Strobach / Signal Processing 89 (2009) 2514–2528

wðtÞ ¼ dðtÞhðtÞ  vðtÞ

(101l)

cðtÞ ¼ Pðt  1ÞvðtÞ

(101m)

2523

Compare this expression with expression (103). Clearly, we can formally introduce dominant eigenvalue and eigenvector estimates as follows:

qðtÞ ¼ 2PT ðt  1ÞwðtÞ

(101n)

yðtÞ ¼ GðtÞ½cðtÞ þ ðqT ðtÞvðtÞÞvðtÞ 1 PðtÞ ¼ ½Pðt  1Þ þ vðtÞqT ðtÞ  yðtÞvT ðtÞ

(101o)

^ r ðtÞV ^ r ðtÞ ¼ V ^ r ðtÞK ^ ðtÞ. U r

(101p)

^ r ðtÞ will be the eigenvalue matrix of ðSðtÞHðtÞÞ. Obviously, K

eðtÞ ¼ dðtÞzðtÞ  Q ðt  1ÞwðtÞ

(101q)

Q ðtÞ ¼ Q ðt  1Þ  2eðtÞvT ðtÞ

(101r)

Consequently, we construct an orthogonal iteration around this matrix product as follows:

a

ENDFOR.

(110)

WðtÞ ¼ ðSðtÞHðtÞÞUðt  1Þ,

(111a)

WðtÞ ¼ UðtÞRðtÞ : QR-decomposition.

(111b)

Since UðtÞUT ðtÞ ¼ UT ðtÞUðtÞ ¼ I and Hu ðtÞ ¼ UT ðt  1ÞUðtÞ, it is not difficult to demonstrate that

5. Additional results The main subject of this section is an inspection of the term Uðt  1ÞDðt  1Þ that has been neglected in the transition from recurrence (8) to recurrence (9). This inspection requires that we have a closer look at the Eigenvalue Decomposition (EVD) of UðtÞ. This EVD can be expressed as follows:

UðtÞ ¼ Vr ðtÞKr ðtÞVTr ðtÞ þ VNr ðtÞKNr ðtÞVTNr ðtÞ,

(102)

where

Ur ðtÞ ¼ Vr ðtÞKr ðtÞVTr ðtÞ

(103)

is the dominant part of this EVD. Clearly, Ur ðtÞ can be obtained by right-multiplication of UðtÞ with a projection matrix that projects vectors orthogonally onto the dominant (rank-r) subspace spanned by Vr ðtÞ:

Ur ðtÞ ¼ UðtÞVr ðtÞVTr ðtÞ.

(104)

On the other hand, our subspace trackers produce, in each time step, a very good approximation to this projection matrix, namely Q ðtÞQ T ðtÞ  Vr ðtÞVTr ðtÞ.

(105)

Therefore, one might be interested in computing a low^ r ðtÞ via the following projection: rank approximant U ^ r ðtÞ ¼ UðtÞQ ðtÞQ T ðtÞ. U

(106)

Unfortunately, this does not seem to be tractable directly. From (2a) to (2b), however, we can easily obtain the projection onto the one time step delayed subspace as follows:

UðtÞQ ðt  1ÞQ T ðt  1Þ ¼ Q ðtÞSðtÞQ T ðt  1Þ.

(107)

In practice, the exponential forgetting factor a will be very close to a value of 1. As a consequence, the dominant principal angle between subsequent subspaces will be very small. In such cases, the following ‘‘projection approximation’’ can be established: T

T

T

T

UðtÞQ ðtÞQ ðtÞ  UðtÞQ ðt  1ÞQ ðt  1ÞQ ðtÞQ ðtÞ.

(108)

ðSðtÞHðtÞÞ ¼ UðtÞRðtÞHu ðtÞUT ðtÞ.

(112)

Substitute this into (109) and compare with (110) to find ^ r ðtÞ ¼ Q ðtÞUðtÞðRðtÞHu ðtÞÞUT ðtÞQ T ðtÞ U T

^ r ðtÞV ^ ðtÞ. ^ r ðtÞK ¼V r

(113)

The following expressions for dominant eigenvalue and eigenvector estimates are extracted: ^ r ðtÞ ¼ diagfRðtÞHu ðtÞg, K

(114a)

^ r ðtÞ ¼ Q ðtÞUðtÞ. V

(114b)

This routine for tracking the dominant eigenvalues and eigenvectors of UðtÞ is casted into equation-array (115a)–(115e). Initialize with Uð0Þ ¼ I. The dominant ^ r ðtÞ are usually not formed explicitly. eigenvectors V WðtÞ ¼ ðSðtÞHðtÞÞUðt  1Þ

(115a)

WðtÞ ¼ UðtÞRðtÞ : QR-decomposition

(115b)

Hu ðtÞ ¼ UT ðt  1ÞUðtÞ

(115c)

^ r ðtÞ ¼ diagfRðtÞHu ðtÞg K

(115d)

^ r ðtÞ ¼ Q ðtÞUðtÞ V

(115e)

This routine can be operated inside the row-Householder subspace tracker, for instance between recursions (81k) and (81l). This method for tracking the dominant eigenvalues is not the only option. As mentioned earlier, the exponential forgetting factor a will be very close to a value of 1 in most applications. In these cases, we observe that HðtÞ  I and hence kvðtÞk2  0. This in turn will cause (recall (81k)), that SðtÞHðtÞ  SðtÞ and SðtÞ  XðtÞ  XT ðtÞ, because both SðtÞ and XðtÞ are almost perfectly symmetric matrices. Hence the eigenvalues of SðtÞHðtÞ and XT ðtÞ will be almost perfectly identical. Consequently, we may construct an orthogonal iteration around XT ðtÞ: WðtÞ ¼ XT ðtÞUðt  1Þ,

(116a)

WðtÞ ¼ UðtÞRðtÞ : QR-decomposition.

(116b)

Hence, we re-define our low-rank approximant of UðtÞ as follows:

In each time step, the following decomposition is obtained:

^ r ðtÞ ¼ UðtÞQ ðt  1ÞQ T ðt  1ÞQ ðtÞQ T ðtÞ U

XT ðtÞ ¼ UðtÞRðtÞUT ðt  1Þ.

T

¼ Q ðtÞSðtÞQ ðt  1ÞQ ðtÞQ ðtÞ ¼ Q ðtÞðSðtÞHðtÞÞQ T ðtÞ.

(117)

T

(109)

This has its attractions. Recall the key system of linear equations (81e). Obviously, we may substitute XT ðtÞ by

ARTICLE IN PRESS 2524

P. Strobach / Signal Processing 89 (2009) 2514–2528

of dðtÞ according to (77b). This yields

decomposition (117) to obtain UðtÞRðtÞUT ðt  1ÞbðtÞ ¼ Z 1=2 ðtÞhðtÞ,

(118)

T

Q T ðtÞDðtÞ ¼ Q T ðtÞz? ðtÞd ðtÞ ¼  2jðtÞZ

or equivalently

ðtÞQ T ðtÞ½zðtÞ  Q ðt  1ÞhðtÞvT ðtÞ

¼  2jðtÞZ 1=2 ðtÞ½Q T ðtÞzðtÞ

RðtÞUT ðt  1ÞbðtÞ ¼ Z 1=2 ðtÞUT ðtÞhðtÞ.

(119)

Define rðtÞ ¼ Z

1=2

 HT ðtÞhðtÞvT ðtÞ.

(125) T

1=2

T

ðtÞU ðtÞhðtÞ,

T

sðtÞ ¼ U ðt  1ÞbðtÞ,

(120a) (120b)

to see that the desired parameter vector bðtÞ can now be computed via simple back-substitution as follows: RðtÞsðtÞ ¼ rðtÞ ! back-substitution ! sðtÞ,

(121a)

bðtÞ ¼ Uðt  1ÞsðtÞ.

(121b)

Hereby, we killed two flies with one blow, because we combined orthogonal iteration and linear system solver in one operation package. The quasicode of this method is casted into equation-array (122a)–(122g): WðtÞ ¼ XT ðtÞUðt  1Þ WðtÞ ¼ UðtÞRðtÞ : QR-decomposition

The key expression here is Q ðtÞzðtÞ. Obviously, this represents a data compression step using the actual subspace estimate Q ðtÞ. This is new. The ‘‘classical’’ orthogonal iteration subspace trackers use only the ‘‘old’’ subspace estimate Q ðt  1Þ for input data compression according to (recall (3)) hðtÞ ¼ Q T ðt  1ÞzðtÞ. There exists a necessary refinement step in finally projecting zðtÞ orthogonally onto the subspace spanned by Q ðtÞ. This can be accomplished using a low-rank approximant like (109) that is based on a consecutive projection onto the subspaces spanned by Q ðt  1Þ and Q ðtÞ. The term Q T ðtÞzðtÞ can be computed efficiently. To see this, recall the subspace propagation model (5) in transposed form. We obtain

(122a)

Q T ðtÞzðtÞ ¼ HT ðtÞhðtÞ  2jðtÞvðtÞzT? ðtÞzðtÞ.

(122b)

Now it is relatively easy to verify that

rðtÞ ¼ Z 1=2 ðtÞUT ðtÞhðtÞ

(122c)

RðtÞsðtÞ ¼ rðtÞ ! back-substitution ! sðtÞ

(122d)

zT? ðtÞzðtÞ ¼ Z 1=2 ðtÞ.

bðtÞ ¼ Uðt  1ÞsðtÞ ^ r ðtÞ ¼ diagfRðtÞg K

(122e) (122f)

^ r ðtÞ ¼ Q ðtÞUðtÞ V

(122g)

This routine can be used, for instance, in place of the ordinary linear system solver (81e). Of course, the matrix ^ r ðtÞ must not be computed of tracked eigenvectors V explicitly. We are now well equipped for the desired inspection of the previously neglected term Uðt  1ÞDðt  1Þ. Discover that this term contains a computable component. It is the component that represents the signal terms in Uðt  1Þ. This component should not be ignored. It can be extracted by low-rank approximation according to (109). Replace ^ r ðt  1Þ according Uðt  1Þ by its low-rank approximant U to (109). Hereby, we should be able to extract the signalrelated components in Uðt  1ÞDðt  1Þ. This should lead us to an even more exact variant of a row-Householder subspace tracker. Begin with the following approximation:

Q ðtÞDðtÞ ¼ ?.

Hence the following compact expression is obtained: Q T ðtÞzðtÞ ¼ HT ðtÞhðtÞ  2jðtÞZ 1=2 ðtÞvðtÞ.

(128) T

Substitute (128) into (125) to see that the terms H ðtÞhðtÞ cancel out perfectly and Z 1=2 ðtÞZ 1=2 ðtÞ ¼ 1. Hence we obtain Q T ðtÞDðtÞ ¼ 4j2 ðtÞvðtÞvT ðtÞ.

(129)

Consequently, (123) attains the following form: ^ r ðt  1ÞDðt  1Þ ¼ 4j2 ðt  1ÞQ ðt  1ÞSðt  1Þ U Hðt  1Þvðt  1ÞvT ðt  1Þ.

(130)

Now recall (8)–(10), where we simply neglected the term Uðt  1ÞDðt  1Þ. A refined form of fast orthogonal iteration is now constituted on the basis of the following recurrence: ^ r ðt  1ÞDðt  1Þ þ zðtÞhT ðtÞ. AðtÞ ¼ aAðt  1ÞHðt  1Þ þ aU (131)

Q ðtÞSðtÞ ¼ aQ ðt  1ÞSðt  1ÞHðt  1Þ þ 4aj2 ðt  1ÞQ ðt  1Þ T

Sðt  1ÞHðt  1Þvðt  1ÞvT ðt  1Þ þ zðtÞh ðtÞ

(123)

This problem is solved as soon as we find a compact expression for the term Q T ðt  1ÞDðt  1Þ. A lot of index space can be saved, if we consider this problem at time instant ðtÞ. Hence the question we wish to answer is the following: T

(127)

Recall (130) and write

^ r ðt  1ÞDðt  1Þ Uðt  1ÞDðt  1Þ  U ¼ Q ðt  1ÞSðt  1ÞHðt  1Þ Q T ðt  1ÞDðt  1Þ.

(126)

(124)

Recall the definition of DðtÞ according to (78), the definition of z? ðtÞ according to (72), and the definition

¼ aQ ðt  1ÞSðt  1ÞHðt  1Þ½I þ 4j2 ðt  1Þ T

vðt  1ÞvT ðt  1Þ þ zðtÞh ðtÞ. A

final

concern

is

the

(132) evaluation

of

Hðt  1Þ½I þ 4j2 ðt  1Þvðt  1ÞvT ðt  1Þ. Use Hðt  1Þ ¼ I  2vðt  1ÞvT ðt  1Þ and vT ðt  1Þvðt  1Þ ¼ 1  j2 ðt  1Þ to demonstrate that

Hðt  1Þ½I þ 4j2 ðt  1Þvðt  1ÞvT ðt  1Þ ¼ I  2cðt  1Þvðt  1ÞvT ðt  1Þ,

(133)

ARTICLE IN PRESS P. Strobach / Signal Processing 89 (2009) 2514–2528

where

follows:

cðt  1Þ ¼ 1 þ 2j2 ðt  1Þ  4j4 ðt  1Þ.

(134)

Substitute (133) into (132). This yields Q ðtÞSðtÞ ¼ aQ ðt  1ÞSðt  1Þ½I  2cðt  1Þvðt  1ÞvT ðt  1Þ T

T

þ Z 1=2 ðtÞz? ðtÞh ðtÞ þ Q ðt  1ÞhðtÞh ðtÞ,

(135)

where we used again the orthonormal decomposition of zðtÞ according to (46). Finally, employ uðt  1Þ ¼ Sðt  1Þvðt  1Þ according to (80a) and verify that (135) can be expressed in the following standard form for rotational or reflectional updating: Q ðtÞSðtÞ ¼ ½Q ðt  1Þ z? ðtÞ 2 3 aSðt  1Þ  2acðt  1Þuðt  1ÞvT ðt  1Þ þ hðtÞhT ðtÞ 4 5.  T Z 1=2 ðtÞh ðtÞ

(136) Compare this refined form with form (48), as obtained by simple neglection of Uðt  1ÞDðt  1Þ. Clearly, (48) can be rewritten as follows (recall (77a) and (80a)): Q ðtÞSðtÞ ¼ ½Q ðt  1Þ z? ðtÞ 2 3 aSðt  1Þ  2auðt  1ÞvT ðt  1Þ þ hðtÞhT ðtÞ 5. 4 T Z 1=2 ðtÞh ðtÞ (137) This result deserves a deeper discussion:

 Obviously, the refined approach according to (136) and 





2525

the ‘‘classical’’ approach according to (137) differ only by the scalar ‘‘gain’’ factor cðt  1Þ. Under usual application conditions with a close to a value of 1, our previous considerations apply. Recall that this case is governed by kvðt  1Þk2  0, and consequently, in view of the unit-norm condition of Householder vectors (57), by jðt  1Þ  1. Hence we recall the definition of cðt  1Þ according to (134) and conclude that under these circumstances, we will have: cðt  1Þ  1. This condition is indeed almost exactly fulfilled in practice. Therefore, we can simply work with a fixed value of c ¼ 1. This is altogether a striking result, because it means that we obtain the optimal solution just by inverting the sign in the correction term associated with Hðt  1Þ in the ‘‘classical’’ approach that was regarded optimal until recently. In other words, this means that instead of working with Hðt  1Þ ¼ I  2vðt  1ÞvT ðt  1Þ, we  ðt  1Þ ¼ Iþ 2vðt  1ÞvT ðt  1Þ. should now work with Y Most ironically, this really means that a classical method like LORAF3, which assumes Hðt  1Þ ¼ I (or c ¼ 0), is closer to the true optimal solution than methods like LORAF2 or FAPI, which are based on Hðt  1Þ ¼ I  2vðt  1ÞvT ðt  1Þ (or c ¼ þ1).

In our fast row-Householder subspace tracker, only recursion (81d) will be affected. We should implement this recursion in sign-changed form as

T

XðtÞ ¼ aSðt  1Þ þ 2auðt  1ÞvT ðt  1Þ þ hðtÞh ðtÞ,

(138)

or should even use the exact c-factor as given in (134). In case of the S1 -domain or PAST-type trackers, the situation appears seriously complicated. Notice that now, the basic recurrence (11) underlying all the currently existing algorithms of the PAST-type will not hold any longer. This basic recurrence should be replaced by the modified form:  ðt  1Þ þ hðtÞhT ðtÞ, HðtÞSðtÞ ¼ aSðt  1ÞH

(139)

where  ðt  1Þ ¼ I þ 2vðt  1ÞvT ðt  1Þ. H

(140)

It is immediately clear that in this modified case, the Badeau–David–Richard-recursion (33) will not hold any longer and it is not known at this time, if there will be any tractable PAST-type solution for this modified case. We leave this problem as an interesting challenge for the PAST-community. Fortunately, the subspace estimates are relatively insensitive against modifications of the c-parameter. Such a statement does not hold for the estimated powers, and in particular, the tracked eigenvalues. In the experimental section, we demonstrate that a use of the wrong (previously ‘‘optimal’’) c ¼ þ1 causes a considerable bias in the estimated eigenvalues. This bias vanishes completely, as soon as the new optimal solution with a fixed c ¼ 1 according to (138)–(140) is employed. The subspace and eigenvector estimates show little or no dependency on the specific choice of the c-parameter. These deeper theoretical insights justify our previous recommendation of the LORAF3 concept. With a parameter choice of c ¼ 0, this method is the golden intermediate between right and wrong. It saves some calculations and can be safely used as long as subspace tracking is the only issue. In view of these aspects, it appears literally ridiculous that several authors reported a degraded performance of LORAF3. See, for instance [7,27]. These authors were trapped by a widely distributed unapproved simulation software package that contained the LORAF3 algorithm implemented with coding errors. As soon as exact eigenvalue or mode considerations become an issue, we should switch to the c ¼ 1 method. A PAST-type implementation of this method, however, does not exist at this time. 6. Computer experiments The fast recursive row-Householder subspace tracker has been implemented and tested in a standard single precision F90 simulation environment. We show experimental results for this algorithm (81a)–(81n). The eigenvalues were tracked according to (122a)–(122g). We used this set of recursions as a substitute for the linear system solver (81e) in cases where eigenvalue tracking became an issue.

ARTICLE IN PRESS 2526

P. Strobach / Signal Processing 89 (2009) 2514–2528

The following deviation measure is computed:

6.1. Data model and simulation criteria The record length for each trial run is 10 000 time steps. The problem order is N ¼ 32. The exponential forgetting factor is fixed to a value of a ¼ 0:996. In each time step, a data snapshot vector 3 2 z1 ðtÞ 7 6 6 z2 ðtÞ 7 7 6 (141) zðtÞ ¼ 6 . 7 6 .. 7 5 4 zN ðtÞ is generated, where zk ðtÞ ¼ sk ðtÞ þ nk ðtÞ;

1pkpN.

(142)

eQQ ðtÞ ¼ 10 log10 ðtrfETQQ ðtÞEQQ ðtÞgÞ ðdBÞ.

(2) The deviation from orthonormality criterion: The following basis grammian is computed: FQQ ðtÞ ¼ Q T ðtÞQ ðtÞ  I.

f QQ ðtÞ ¼ 10 log10 ðtrfFTQQ ðtÞFQQ ðtÞgÞ ðdBÞ

T

^ ðtÞVr ðtÞ. CðtÞ ¼ V r

(143)

cj ðtÞ ¼

sk ðtÞ ¼ a21 cosðot21 t þ os21 k þ o021 Þ

Ps , Pn

(145)

s2k ðtÞ,

(146a)

n2k ðtÞ.

(146b)

with N 10 000 X X k¼1

Pn ¼

t¼1

N 10 000 X X k¼1

 M   1X 180 a cosðjC j;j ðtÞjÞ ð Þ;   M i¼1 p

1pjpr.

(153)

(144)

We operated the subspace tracker with r ¼ 4, which is the necessary subspace dimension or rank to accommodate a signal of this kind in the stationary case. The signal-tonoise ratio (SNR) was set to a value of 3 dB where

Ps ¼

(152)

The main diagonal elements of CðtÞ can be regarded as the cosines of the angles of associated true and tracked eigenvectors. We show the averaged tracks of the absolute values of these angles. These average angle deviation tracks are computed as follows:

At time step t ¼ 4001, the model changes abruptly to

SNR ¼ 10 log10

(151)

is a measure of deviation from orthonormality. (3) The average eigenvector angle mismatch criterion: A grammian of true and tracked eigenvector matrices is formed:

sk ðtÞ ¼ a11 cosðot11 t þ os11 k þ o011 Þ

þ a22 cosðot22 t þ os22 k þ o022 Þ.

(150)

Then

Here sk ðtÞ denotes the signal, and nk ðtÞ denotes a stationary white Gaussian noise component. In the interval 1ptp4000, the signal component is given by the following sum of two cosine signals:

þ a12 cosðot12 t þ os12 k þ o012 Þ.

(149)

6.2. Results Fig. 1 shows the eQQ ðtÞ curves for each of the M ¼ 10 trial runs plotted in one diagram. These tracks were computed with a row-Householder subspace tracker of the new c ¼ 1 type. (We used (138) instead of (81d) in the algorithm.) We observed, however, that there occurs no visible difference in the result, if we choose c arbitrary in the interval 1pcp þ 1. From the pure subspace tracking standpoint, this result suggests that we use c ¼ 0 and omit (81c), and simplify (81d) accordingly. This is then a fast form of LORAF3, which produces literally

t¼1

The specific parameter values are given as follows:

a12 ¼ 1:6; a21 ¼ 2:0; a22 ¼ 1:0;

ot11 ¼ 2:2; ot12 ¼ 2:8; ot21 ¼ 2:7; ot22 ¼ 2:3;

os11 ¼ 0:5, os12 ¼ 0:9, os21 ¼ 1:1, os22 ¼ 0:8.

10 5

(147)

The parameters o011 , o012 , o021 , and o022 are random initial phase angles which are fixed for each trial run. Each test comprises M ¼ 10 independent trial runs. In each trial run, the following criteria were monitored: (1) The deviation of true and tracked subspaces criterion: A true EVD of the exponentially updated UðtÞ is computed as a reference in each time step. Recall (1) and (102). A difference of subspace projection matrices is formed: EQQ ðtÞ ¼ Vr ðtÞVTr ðtÞ  Q ðtÞQ T ðtÞ.

(148)

subspace error [dB]

a11 ¼ 1:4;

0 -5 -10 -15 -20 -25 -30 0

2000

4000 6000 discrete time

8000

10000

Fig. 1. Deviation of true and tracked subspaces criterion for 10 independent trial runs plotted in one diagram.

ARTICLE IN PRESS

-90

12000

-100

10000

-110

eigenvalues

orthonormality error [dB]

P. Strobach / Signal Processing 89 (2009) 2514–2528

-120 -130

2527

8000 6000 4000

-140

2000

-150

0

-160 0

2000

4000 6000 discrete time

8000

0

10000

Fig. 2. Deviation from orthonormality criterion for 10 independent trial runs plotted in one diagram.

2000

4000 6000 discrete time

8000

10000

Fig. 4. r ¼ 4 true eigenvalues and the corresponding tracked eigenvalues displayed as a function of time for the case c ¼ þ1. Average tracks of 10 independent trial runs.

12000

eigenvalues

8000 6000 4000 2000 0 0

2000

4000

6000

8000

10000

discrete time Fig. 3. r ¼ 4 true eigenvalues and the corresponding tracked eigenvalues displayed as a function of time for the case c ¼ 1. Average tracks of 10 independent trial runs.

identical results in terms of this deviation of true and tracked subspaces criterion. The corresponding deviation from orthonormality f QQ ðtÞ curves are shown in Fig. 2. The orthonormality error according to this criterion is always less than 120 dB. Expressed in degrees, this means that the deviation from orthogonality is typically only a few microdegrees and is practically ideal. Here, again, we observed little influence of the c parameter on this orthonormality performance, but a minor improvement has been observed for the new choice of c ¼ 1. Fig. 3 shows the corresponding tracks of the r ¼ 4 true, and the corresponding tracked eigenvalues for the c ¼ 1 case. These eigenvalue tracks (computed using (122f)), appear perfectly unbiased. Clearly, the deeper reason is that with the low-rank approximation of Uðt  1Þ in the formerly neglected term Uðt  1ÞDðt  1Þ, we achieved a practically perfect mapping of the complete signal energy into the reduced subspace of dimension r. In Fig. 4, we repeated the eigenvalue tracking experiment for the c ¼ þ1 case. Then the algorithm produces

eigenvector angle mismatch [degrees]

90 10000

80 70 60 50 40 30 20 10 0 0

2000

4000 6000 discrete time

8000

10000

Fig. 5. Average eigenvector angle mismatch criterion displayed for the r ¼ 4 eigenvectors for the case c ¼ 1 as a function of time.

results like LORAF2 or FAPI. We can see that the tracked eigenvalues appear considerably biased (underestimated). This is not surprising, because all the signal energy contained in the term Uðt  1ÞDðt  1Þ is ignored by this approach. The effect is more pronounced at low SNRs. Fig. 5 finally shows the average eigenvector angle mismatch criterion (153) displayed for the r ¼ 4 eigenvectors tracked according to (122g). We can see that these test data represent no simple scenario. In general, eigenvectors are difficult to estimate at low SNRs with clustered eigenvalues. They are also difficult to estimate for the smallest eigenvalues in the tracked set. Here, we observed little influence of the c-factor on the overall quality of the tracked eigenvectors.

7. Conclusions It was the purpose of this paper to shed some light on the class of fast recursive orthogonal iteration subspace tracking algorithms. These trackers can be classified in

ARTICLE IN PRESS 2528

P. Strobach / Signal Processing 89 (2009) 2514–2528

S1 -domain (or PAST-type) and S-domain (or LORAF-type) techniques. The PAST-type trackers exhibit the structure of classical RLS adaptive filters. The LORAF-type trackers are purely rotation or reflection based and avoid any inverse power matrix updating. Hence, they appear more favorable in view of many practical considerations, because with the advent of the row-Householder approach, the disadvantage of computational inefficiency has been completely relaxed. Finally, we observed that we can even incorporate a previously neglected term in the algorithm, with the effect, that the usual eigenvalue bias disappears. Until now, this new and most exact version of a fast orthogonal iteration subspace tracker is only implementable in the new fast row-Householder algorithm framework. References [1] P. Comon, G.H. Golub, Tracking a few extreme singular values and vectors in signal processing, Proceedings of the IEEE, 78(8) (August 1990). [2] B. Yang, Projection approximation subspace tracking, IEEE Transactions on Signal Processing 43 (January 1995) 95–107. [3] P. Strobach, Linear Prediction Theory: A Mathematical Basis for Adaptive Systems, Springer Series in Information Sciences, vol. 21, 1990. [4] K. Abed-Meraim, A. Chkeif, Y. Hua, Fast orthonormal PAST algorithm, IEEE Signal Processing Letters 7 (3) (March 2000) 60–62. [5] Y. Miao, Y. Hua, Fast subspace tracking and neural network learning by a novel information criterion, IEEE Transactions on Signal Processing 46 (7) (July 1998) 1967–1979. [6] Y. Hua, Y. Xiang, T. Chen, K. Abed-Meriam, Y. Miao, A new look at the power method for fast subspace tracking, Digital Signal Processing 9 (2) (October 1999) 297–314. [7] R. Badeau, B. David, G. Richard, Fast approximated power iteration subspace tracking, IEEE Transactions on Signal Processing 53 (8) (August 2005) 2931–2941. [8] P. Strobach, Low rank adaptive filters, IEEE Transactions on Signal Processing 44 (12) (December 1996) 2932–2947. [9] J.W. Daniel, W.B. Gragg, L. Kaufman, G.W. Stewart, Reorthogonalization and stable algorithms for updating the Gram–Schmidt QR factorization, Mathematics of Computation 30 (136) (October 1976) 772–795. [10] P. Strobach, Bi-iteration SVD subspace tracking algorithms, IEEE Transactions on Signal Processing 45 (5) (May 1997) 1222–1240. [11] S. Ouyang, Y. Hua, Bi-iterative least-square method for subspace tracking, IEEE Transactions on Signal Processing 53 (8) (August 2005) 2984–2996.

[12] P. Strobach, Square Hankel SVD subspace tracking algorithms, Signal Processing 57 (1) (February 1997) 1–18. [13] P. Strobach, Fast recursive subspace adaptive ESPRIT algorithms, IEEE Transactions on Signal Processing 46 (9) (September 1998) 2413–2430. [14] P. Strobach, Bi-iteration recursive instrumental variable subspace tracking and adaptive filtering, IEEE Transactions on Signal Processing 46 (10) (October 1998) 2708–2725. [15] P. Strobach, Fast orthogonal iteration adaptive algorithms for the generalized symmetric eigenproblem, IEEE Transactions on Signal Processing 46 (12) (December 1998) 3345–3359. [16] P. Strobach, Equirotational stack parametrization in subspace estimation and tracking, IEEE Transactions on Signal Processing 48 (3) (March 2000) 712–722. [17] P. Strobach, Bi-iteration multiple invariance subspace tracking and adaptive ESPRIT, IEEE Transactions on Signal Processing 48 (2) (February 2000) 442–456. [18] P. Strobach, Two-dimensional equirotational stack subspace fitting with applications to uniform rectangular arrays and ESPRIT, IEEE Transactions on Signal Processing 48 (7) (July 2000) 1902–1914. [19] P. Strobach, Total least squares phased averaging and 3-D ESPRIT for joint azimuth-elevation-carrier estimation, IEEE Transactions on Signal Processing 49 (1) (January 2001) 54–62. [20] H. Rutishauser, Simultaneous iteration method for symmetric matrices, Numerische Mathematik 13 (1969) 4–13. [21] G.H. Golub, C.F. Van Loan, Matrix Computations, second ed., John Hopkins University Press, Baltimore, MD, 1989. [22] N.L. Owsley, Adaptive data orthogonalization, in: Proceedings of the IEEE ICASSP-78, Tampa, FL, 1978, pp. 109–112. [23] S. Haykin, Adaptive Filter Theory, Prentice-Hall Information and Systems Sciences Series, fourth ed., 2001. [24] J. Cioffi, The fast Householder filters RLS adaptive filter, in: IEEE International Conference on Acoustics, Speech, and Signal Processing, ICASSP-90, April 1990, pp. 1619–1622. [25] A.W. Bojanczyk, J.G. Nagy, R.J. Plemmons, Block RLS using row Householder reflections, Linear Algebra and its Applications 188 (1993) 31–62. [26] S.C. Douglas, Numerically robust adaptive subspace tracking using Householder transformations, in: Proceedings of the IEEE Sensor Array Multichannel Signal Processing Workshop, 2000, pp. 499–503. [27] X.G. Doukopoulos, G.V. Moustakides, The fast data projection method for stable subspace tracking, in: 13th European Signal Proceedings Conference EUSIPCO-2005, Antalya, Turkey, September 2005. [28] P.E. Gill, G.H. Golub, W. Murray, M.A. Saunders, Methods for updating matrix factorizations, Mathematics of Computation 28 (1974) 505–535. [29] P. Stange, A. Griewank, M. Bollho¨fer, On the efficient update of rectangular LU-factorizations subject to low-rank modifications, Electronic Transactions on Numerical Analysis 26 (2007) 161–177.