- Email: [email protected]
Sliding window adaptive SVD using the unsymmetric householder partial compressor
ARTICLE IN PRESS Signal Processing 90 (2010) 352–362
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Sliding window adaptive SVD using the unsymmetric householder partial compressor 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 6 January 2009 Received in revised form 14 April 2009 Accepted 6 July 2009 Available online 14 July 2009
A fast algorithm for tracking the rank-r SVD-approximant Q ðtÞPðtÞUT ðtÞ of a sliding window data matrix XðtÞ of dimension L N is introduced, where PðtÞ is a square-root power matrix of dimension r r with rominfL; Ng. This algorithm is based on the unsymmetric Householder partial compressor and uses a reorthonormalizing Householder transformation for downdating. The concept is numerically self-stabilizing and requires no leakage. The dominant complexity is 4Lr þ 3Nr multiplications per time update which is the lower bound in complexity for an algorithm of this kind. Applications occur in the area of adaptive array processing and other forms of adaptive processing in finite duration subspaces. A complete algorithm summary is provided. Computer simulations illustrate the operation of the algorithm. & 2009 Elsevier B.V. All rights reserved.
Keywords: Sliding window subspace tracking Householder Reorthonormalization Adaptive SVD
1. Introduction In this article, we present an algorithm for recursive rank-r SVD approximation of N-channel data matrices X of dimension L N, which are shifted in time according to the following sliding rectangular window rule: " # " # XðtÞ xT ðtÞ ¼ . (1) T x ðt LÞ Xðt 1Þ In each time step, an old data snapshot vector xðt LÞ leaves the window, while a new snapshot vector xðtÞ is appended on the top. This scheme appears frequently in adaptive multichannel/array processing. We present a sliding window subspace tracker that computes, in each time step, a rank-r approximant XðtÞ Q ðtÞPðtÞUT ðtÞ,
(2)
where Q ðtÞ is an orthonormal L r left singular subspace basis, UðtÞ is an orthonormal N r right singular subspace basis, and PðtÞ is an r r square-root power matrix that is 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.07.002
orthonormally or (in the complex case) unitary similar to the leading singular values in the SVD of XðtÞ. Note that PðtÞ is a dense unstructured square matrix. No shape restrictions are imposed on PðtÞ. Approximants of the kind (2) feature an adaptive implementation of various pertinent processing schemes for multichannel and array data. A special algorithm of this kind is the square Hankel SVD subspace tracker of [1]. This algorithm is restricted to time series subspace tracking. Solutions to the general rectangular sliding window subspace tracking problem appeared recently in [2,3]. However, these algorithms are redundant in complexity. A heuristically accelerated variant is proposed in [3]. In [4], we adopt the basic approach suggested in [3], and present a fast and exact solution using some essential mathematical methods known from fast Householder subspace tracking [5–7]. Additionally, we introduce a reorthonormalizing Householder transformation for downdating, a fundamentally new concept for downdating in finite duration QR-type algorithms. The resulting algorithm in [4] is already maximally fast and reaches the 4Lr þ 3Nr complexity lower bound for an algorithm of this kind. However, the method is numerically unstable because the updating recursion for the
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
internal square-root power matrix forms an undamped loop. Undamped loops occur frequently in sliding-window algorithms. They run unstable in floating-point implementations due to round-off error accumulation. A simple remedy is the incorporation of a small leakage in the recursion. This effectively stabilizes the algorithm. The incorporation of leakage has the same effect as exponential forgetting. The resulting window is no longer a true rectangular window, but becomes an exponential window with truncated tail. There are many applications, where this modification can be accepted. Sometimes, it is even desirable in the sense of a ‘‘trim’’ that assigns a larger weight to the incoming data snapshots. On the other hand, there are some pertinent applications in array processing, for instance in ESPRIT spectral estimation [12], where one has to work with overlapping sub-windows. The underlying assumption in such schemes is that the basic tracker window is perfectly rectangular. In such cases, leakagestabilized algorithms cannot be applied. Therefore, we introduce a completely self-stabilizing algorithm that preserves the perfectly rectangular shape of the temporal window. The principle behind this algorithm is an unsymmetric variant of the Householder partial compressor based on the Householder partial compressor theorem introduced in [8]. This scheme extends the estimated subspace by an additional dimension, in which the updating residual is accommodated. The algorithm uses this degree of freedom to wipe out all round-off noise from the subspace into this rank-one residual subspace. This way the rectangular shape of the window is preserved and stability can be maintained without the application of leakage. This paper is organized as follows: In Section 2, we introduce the unsymmetric Householder partial compressor and its application in the updating of the QPUapproximant. Section 3 shows how downdating can be accomplished using a reorthonormalizing Householder transformation. Section 4 summarizes the complete algorithm and presents a particularly streamlined version of the algorithm that is suitable for a recursive LRfactorized representation of the square-root power matrix PðtÞ. In this form, the algorithm has a subdominant complexity of Oðr 2 Þ. Section 5 shows some simulation results in order to illustrate the operation of the algorithm. Section 6 summarizes the conclusions. 2. Updating the QPU-approximant
353
augmented form: " Xðt 1Þ ¼
#
xT ðtÞ T
Q ðt 1ÞPðt 1ÞU ðt 1Þ
.
(4)
The orthonormal decomposition of the incoming snapshot vector xðtÞ with respect to the subspace spanned by Uðt 1Þ is required. The necessary formulas are given below (note this is the same procedure as in conventional subspace tracking with exponential forgetting): hðtÞ ¼ UT ðt 1ÞxðtÞ,
(5a)
x? ðtÞ ¼ xðtÞ Uðt 1ÞhðtÞ,
(5b) T
XðtÞ ¼ xT? ðtÞx? ðtÞ ¼ xT ðtÞxðtÞ h ðtÞhðtÞ, x? ðtÞ ¼ X
1=2
(5c) (5d)
ðtÞx? ðtÞ,
xðtÞ ¼ X 1=2 ðtÞx? ðtÞ þ Uðt 1ÞhðtÞ.
(5e)
Recall that (5a) is the obligatory input data compression step that is common to all fast subspace trackers. Eq. (5b) produces the complement of the orthogonal projection of xðtÞ onto the subspace spanned by the column-vectors of Uðt 1Þ. This constitutes the necessary decomposition of the input data vector into a component that can be represented by a linear combination of the information that is already available in Uðt 1Þ, and an innovation which is represented by x? ðtÞ. The squared norm of this innovation is computed in (5c). Formula (5e) finally represents the desired orthonormal decomposition of xðtÞ with respect to the subspace spanned by Uðt 1Þ. Now it is not difficult to verify that (4) can also be expressed in a factorized form as follows: 3 32 1=2 2 T 1 00 X ðtÞ h ðtÞ # 7" 76 60 7 76 0 xT? ðtÞ 6 7 7 6 Xðt 1Þ ¼ 6 . .. T 7 76 6 .. 4 . Q ðt 1Þ 54 . Pðt 1Þ 5 U ðt 1Þ 0
0 (6)
2.1. The unsymmetric Householder partial compressor We will now demonstrate that the above augmented ‘‘old’’ QPU-approximant (6) can be updated using two Householder matrices according to the Householder compressor principle developed in [8]. For this purpose, introduce the following two Householder matrices: " # 1 2j21 2j1 vT1 , (7) H1 ¼ 2j1 v1 I 2v1 vT1 with
In this section, we develop the necessary details required in the updating of a rank-r QPU-approximant using the unsymmetric Householder partial compressor. Recall the sliding window updating scheme (1). The rank-r QPU-approximant at time step t 1 is defined as follows (recall (2)): ^ 1Þ ¼ Q ðt 1ÞPðt 1ÞUT ðt 1Þ. Xðt
(3)
We wish to update all components of this decomposition in time upon the arrival of a new data vector xðtÞ at time step t. To accomplish this, consider the following
vT1 v1 þ j21 ¼ 1,
(8a)
H1 H1 ¼ I,
(8b)
and " H2 ¼
1 2j22
2j2 vT2
2j2 v2
I 2v2 vT2
# ,
(9)
with vT2 v2 þ j22 ¼ 1,
(10a)
H2 H2 ¼ I.
(10b)
ARTICLE IN PRESS 354
P. Strobach / Signal Processing 90 (2010) 352–362
The parameters of these Householder matrices can be determined so that the augmented QPU-approximant of the previous time step is updated as follows: 3 2 3 2 T 1 00 X 1=2 ðtÞ h ðtÞ 7 6 7 6 7 6 0 7 60 7 6 7 6 7 Xðt 1Þ ¼ 6 . 7H1 H1 6 7 6 . 6 .. Q ðt 1Þ 7 . 6 . Pðt 1Þ 7 5 4 5 4 0
2
H2 H2 4
3
xT? ðtÞ UT ðt 1Þ 2 1=2
6 6 6 ¼ ½qðtÞ Q ðtÞ6 6 6 4
lmin
0
5 00
3
7" 7 uT ðtÞ # 7 7 7 UT ðtÞ ZðtÞ 7 5
0 .. . 0
T
1=2
¼ Q ðtÞZðtÞU ðtÞ þ lmin qðtÞuT ðtÞ.
(11)
Observe that the two Householder matrices in this updating scheme must be adjusted so that the energy in the augmented P-matrix of the previous time step is perfectly mapped into the r r interior matrix ZðtÞ. The inevitable residual that remains after this compression 1=2 step is the minor singular value lmin of Pðt 1Þ in the upper-left corner of the augmented matrix after compression. With the appearance of this particular updating scheme, the discussion about optimality in fast subspace tracking becomes largely obsolete, because a subspace tracker based on this scheme will obviously attain the highest possible performance, since it is based on the highest attainable compression in each time step. Verify again (11) to see that the remaining residual layer l1=2 qðtÞuT ðtÞ, which is popped off in each time step, will min 1=2 always be of minimum energy, because the factor lmin is the minimum attainable residual in a problem of this kind, and the fact that tracefqðtÞuT ðtÞg ¼ 1 is always guaranteed as a consequence of the unconditional unitnorm property of the vectors qðtÞ and uðtÞ. All that remains is to determine the rules for computing the parameters of the two Householder compressor matrices that act in this scheme. Since this is just the unsymmetric variant of the Householder compressor that has been introduced already in [8], we can easily extend these results to the unsymmetric case. For this purpose, introduce the following compact notation: 3 2 1=2 T h ðtÞ X ðtÞ 7 6 7 6 0 7. (12) Pðt 1Þ ¼ 6 . 7 6 . 4 . Pðt 1Þ 5 0 Now the core compression problem of (11) can be expressed as follows: 3 2 1=2 lmin 0 0 7 6 7 6 0 7 ¼ H1 Pðt 1ÞH2 . 6 . (13) 7 6 . 4 . ZðtÞ 5 0
Introduce the SVD of Pðt 1Þ in a partitioned form 3 2 1=2 lmin 0 0 7" vT # 6 7 s 6 0 7 T , Pðt 1Þ ¼ ½us Us 6 7 6 .. 4 . K1=2 5 Vs
(14)
s
0
l1=2 min
where is the minor (smallest) singular value of Pðt 1Þ, us and vs are the corresponding left and right minor singular vectors, Ks1=2 is the diagonal matrix of the leading (larger) singular values, and Us and Vs are the corresponding left and right singular vectors. Now introduce two arbitrary orthonormal rotors Hu and Hv of dimension r r. Apparently, we can write 3 2 1=2 lmin 00 " # 7 6 T us 7 6 0 7. 6 (15) Pðt 1Þ½v V H ¼ s s v 7 6 .. T 1=2 HTu UTs 4 . Hu Ks Hv 5 0 The important insight at this point is the following: there always exists a pair of such rotors Hu and Hv so that the partially rotated singular vector matrices on the right and left sides of Pðt 1Þ in (15) become Householder matrices. Consequently, we can write 2 3 2 3 uTs 1 2j21 2j1 vT1 4 5 ¼ H1 ¼ 4 5, (16a) 2j1 v1 I 2v1 vT1 HTu UTs 2 3 1 2j22 2j2 vT2 4 5. ½vs Vs Hv ¼ H2 ¼ (16b) 2j2 v2 I 2v2 vT2 From these equations, it becomes apparent that these Householder matrices are already uniquely determined by the minor singular vectors us and vs . The Y-matrices are never computed explicitly. The conditional equations for the desired Householder matrix parameters are directly obtained from the top row of (16a) and the leftmost column of (16b), respectively: 3 " # 2 1 2j21 u0 4 5, ¼ (17a) us ¼ uss 2j1 v1 3 " # 2 1 2j22 v0 4 5. ¼ vs ¼ (17b) 2j2 v2 vss A quick calculation yields rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 u0 1 ; v1 ¼ j1 ¼ uss , 2 2j1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 v0 1 ; v2 ¼ j2 ¼ vss , 2 2j2
(18a) (18b)
and the necessary parameters of the Householder partial compressors are uniquely determined. As a by-product, we may observe that the SVD of ZðtÞ is given as follows: ZðtÞ ¼ HTu K1=2 s Hv ,
(19)
but this result is of no interest in this context. The minor singular vectors of Pðt 1Þ are tracked using a variant of bi-inverse iteration. Inverse iteration, in general, is a very important concept in the area of
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
eigenvalue and singular value computation [16]. Biinverse iteration comprises the following iterations (where we omitted the time index for the moment): for
m ¼ 1; 2; :::; m 1
am ¼ P
T
½uðtÞ UðtÞ ¼ ½x? ðtÞ Uðt 1ÞH2 .
(24)
UðtÞ ¼ 2j2 x? ðtÞvT2 þ Uðt 1ÞðI 2v2 vT2 Þ
vm
¼ 2j2 X 1=2 ðtÞ½xðtÞ Uðt 1ÞhðtÞvT2
T
um ¼ ðbm bm Þ1=2 bm endfor
Next consider the right singular subspace update. According to (11), we can write
Now we can exploit the structure of H2 as before and additionally, we eliminate the explicit calculation of x? ðtÞ by substituting (5d) together with (5b) into (24). This yields
um1
vm ¼ ðaTm am Þ1=2 am bm ¼ P
355
(20)
þ Uðt 1ÞðI 2v2 vT2 Þ
(25)
or more compactly: The practical tracker is based on an initial LR-(lower-left upper-right)-factorization of Pðt 1Þ: Pðt 1Þ ¼ LR.
(21)
In each time step, we initialize u0 ¼ us ðt 1Þ and iterate: for
m ¼ 1; 2; . . . ; m Ld ¼ um1 ! f orward substitution ! d Ram ¼ d ! back substitution ! am
fðtÞ ¼ Uðt 1Þðv2 j2 X 1=2 ðtÞhðtÞÞ,
(26a)
UðtÞ ¼ Uðt 1Þ 2ðj2 X 1=2 ðtÞxðtÞ þ fðtÞÞvT2 ,
(26b)
where we introduced an auxiliary vector fðtÞ. Finally, we arrive at the necessary update for the Zmatrix according to (11). We take into account the special structure of the two Householder matrices used in this update and obtain
c ¼ 4j1 ðj2 X 1=2 ðtÞ þ hT ðtÞv2 Þ,
vm ¼ ðaTm am Þ1=2 am T
ZðtÞ ¼ ðI 2v1 vT1 ÞPðt 1ÞðI 2v2 vT2 Þ
T
þ v1 ðcvT2 2j1 h ðtÞÞ,
R e ¼ vm ! f orward substitution ! e
T
L bm ¼ e ! back substitution ! bm T
um ¼ ðbm bm Þ1=2 bm endfor
(22)
Experiments revealed that m ¼ 2 iterations are usually sufficient for convergence. The desired minor singular vectors are obtained as us ðtÞ ¼ um and vs ðtÞ ¼ vm . At t ¼ 0, this minor singular vector tracker is initialized with a random unit-norm vector us ð0Þ. In Section 4, we show that the LR-factorization (21) must not be computed completely anew in each time step. We formulate a completely recursive variant of the algorithm, where the L- and R-factors are updated separately; thus leading to an overall subdominant operations count of only Oðr 2 Þ. Hence the final algorithm will also be optimal in terms of its subdominant operations count since no Oðr 3 Þ computations will be required anymore. 2.2. Final updating recursions Concluding this section, we summarize the resulting updating recursions for the two subspace estimates and the square-root power matrix. Recall the basic update equation (11) and the associated Householder matrices (7) and (9). We can see that qðtÞ is an uninteresting byproduct in the updating process and can be omitted. Taking into account the structure of Householder matrix (7), the following compact expression for Q ðtÞ is obtained: " # 2j1 vT1 Q ðtÞ ¼ . (23) Q ðt 1ÞðI 2v1 vT1 Þ
(27a)
(27b)
where c is again an auxiliary quantity. Hereby, the updating step is complete. The residual layer 1=2 lmin qðtÞuT ðtÞ is popped off the updated approximant (this means the residual layer is not computed explicitly) and we finally obtain the updated rank-r approximant in the following form: ^ 1Þ ¼ Q ðtÞZðtÞUT ðtÞ. Xðt
(28)
Recall that this updated approximant is of extended dimension ðL þ 1Þ N. 3. Downdating the QPU-approximant This section is devoted to the downdating part of the algorithm. Downdating is quite tricky, because it employs some unconventional ideas which are not familiar from ordinary subspace tracking with exponential forgetting. Moreover, the underlying idea is also not familiar from related downdating problems in sliding window recursive least-squares (RLS) adaptive filtering, where a hypernormal Householder transformation was used for downdating [9,10]. We propose a different concept based on a reorthonormalizing Householder transformation here. First of all, notice that only the Q Z-part in (28) requires downdating. The U-subspace requires no downdating. In the updating step, the vertical dimension of the Q subspace has grown by one. Now the new left singular subspace basis must be constructed from the upper L rows of Q ðtÞ. For this purpose, introduce the following partition: W (29) Q ðtÞ ¼ T , b
ARTICLE IN PRESS 356
P. Strobach / Signal Processing 90 (2010) 352–362 T
where b is the bottom row vector of Q ðtÞ. Introduce a generalized form of a Householder-type matrix as follows:
C ¼ I bbbT .
(30)
Define T
B¼b b
(31)
and introduce a scaled version of b as follows: c ¼ B1=2 b.
(32)
Now C can be expressed in true Householder style with a scaled unit-norm rank-one layer:
C ¼ I bbbT ¼ I bBccT ¼ I rccT .
(33)
Clearly, this matrix turns into a conventional Householder matrix for the special case of r ¼ 2. We will next demonstrate that the free parameter r in the Householder-type matrix (33) can be tuned so that the upper submatrix W in (29) is transformed into a matrix with perfectly orthonormal columns. As soon as this is achieved, the remaining bottom row in the updated problem can be dropped, and the sliding window downdating transformation is complete. Notice that this procedure is fundamentally different from classical sliding window recursive least-squares (RLS) downdating using hypernormal Householder transformations as in [9,10]. In contrast to our G-reflector (33), a hypernormal Householder matrix (such as used in [9]) has the following form: H ¼ U 2vvT ,
(34)
and hence our G-reflector of (33) cannot be regarded as a hypernormal Householder matrix in the classical sense anymore. We use a variable scalar parameter r instead. This parameter is tuned so that the submatrix W in (29) is transformed into a matrix with perfectly orthonormal columns. A preliminary study has revealed that we could develop alternative forms of sliding window RLS algorithms on this basis as well. Now let us see how we can benefit from the G-concept in the downdating of the Q Z-part of (28). Using (29) and (30), we can write W 1 Q ðtÞZðtÞUT ðtÞ ¼ ZðtÞUT ðtÞ. (39) T CC b We shall demonstrate that there exists a special choice for the parameter r so that the application of C from the right perfectly re-orthonormalizes W. In this case, downdating is accomplished with Q ðtÞ ¼ WC ! Q T ðtÞQ ðtÞ ¼ I. On the other hand, the P-matrix is constructed via PðtÞ ¼ C1 ZðtÞ,
(41)
and the downdating step is complete. All that remains is to determine the parameter r so that the orthonormality condition (40) is fulfilled. For this purpose, recall that the updated Q -matrix of the previous time step has orthonormal columns: T
Q ðtÞQ ðtÞ ¼ I.
(42)
Moreover, a look at partitioning scheme (29) reveals that we can write T
where U is a diagonal matrix that contains þ1’s and 1’s on the main diagonal in any arbitrary order. We must have at least one 1 on the main diagonal, otherwise H reduces to an ordinary Householder matrix. In any event, it is clear that such a U satisfies:
UU ¼ I and UUU ¼ U. If, for instance, we adjust 3 2 0 6 .. 7 6 I . 7 7, U¼6 7 6 0 5 4 0 0 1
(35)
T
Q ðtÞQ ðtÞ ¼ WT W þ bb ¼ I.
(43)
Consequently, T
WT W ¼ I bb .
(44)
Now since in (40) we assumed that Q ðtÞ ¼ WC with Q T ðtÞQ ðtÞ ¼ I, we can also write
CWT WC ¼ CðI bbT ÞC ¼ I.
(36)
(40)
(45)
Substitute the C of (30) into (45) to obtain the desired conditional equation for the b-parameter that converts the G-reflector into a perfect orthonormalizer for the W-matrix: T
T
T
ðI bbb ÞðI bb ÞðI bbb Þ ¼ I.
(46)
then we could use this concept to develop an algorithm that discards data from the tail of the sliding window, if we only constrain v so that
A little algebraic exercise reveals that the b’s which satisfy (46), also satisfy the following quadratic equation:
vT Uv ¼ 1,
(37)
b2 b þ
(38)
which has the solution 1 1 b1;2 ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi , B 1B
because then: HUH ¼ U.
This is the underlying idea in classical Householder sliding window RLS algorithms as advocated in [9]. On the other hand, our concept based on the Greflector is fundamentally different from the classical hypernormal approach in that we are working with U ¼ I
2 B
1 ¼ 0, BðB 1Þ
(47)
(48)
or equivalently (using (33)) 1
r1;2 ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi . 1B
(49)
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
In most Householder problems, it makes little difference which one of the two solutions is actually pursued. But it turned out that this is clearly not the case here. Although both the ‘þ’, as well as the ‘’ sign solution allow a downdating according to (40), there will still be a significant difference in the practical result obtained for the two cases. The ‘þ’ sign solution produces a Q ðtÞ with a much higher orthonormality error than the corresponding ‘’ sign solution. Computer experiments have revealed that the orthonormality error may grow by a significant amount in some cases when the ‘þ’ sign solution for the r-parameter is used. In order to understand the reasons for this effect, we should have a look at the special case of a vanishing T bottom row vector b ¼ ½0 0 in Q ðtÞ according to (29). Clearly, since Q ðtÞ is always bounded to orthonormal columns, it follows that for a vanishing bottom row vector W will also show perfectly orthonormal columns. Hence no downdating transformation will be required anymore because the downdating solution is directly given by Q ðtÞ ¼ W. This in turn means (recall (40)) that we can assume C ¼ I or r ¼ 0 in this particular case, because no reorthonormalization will be required anymore in this T case. Now since b ¼ ½0 0 causes B ¼ 0 according to (31), we can check in (49), which one of the two choices for the sign in (49) produces the desired r ¼ 0 for a given B ¼ 0. Clearly, we can see that this is perfectly the case for the ‘’ sign. But what happens in the case of the ‘þ’ sign solution? Again, we verify from (49) that in this case, we obtain r ¼ 2 and for this choice, C according to (33) turns into a conventional Householder matrix. Obviously, this is formally correct because the application of any ordinary Householder matrix to W will not destroy the already orthonormal character of W. Nevertheless, a completely unnecessary transformation is carried out. Experimental results revealed that this has a somewhat deteriorating effect on the numerical orthonormality of the resulting Q ðtÞ, because it generates more round-off error than the ‘‘null-transform’’ solution associated with the favourable ‘’ sign. Therefore, our choice is always: 1
r ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi . 1B
(50)
With this important conclusion, the downdating operation is almost complete. All that remains is to determine the corresponding inverse transform matrix that is required in the downdating of the Z-matrix according to (41). This is again an easy exercise, if one knows how to work with matrices of the type ‘identity minus rank-one’. In this case, it is clear that we can assume
C1 ¼ I gccT .
(51)
The necessary link between g and r is then determined by evaluating
CC1 ¼ ðI rccT ÞðI gccT Þ ¼ I.
(52)
357
For safety reasons, we check the special case B ¼ 0, which gives us g ¼ 0, the desired solution for an inactive reflector. Hereby the algorithm is complete. 4. Implementation A raw pseudocode of the algorithm is listed in Table 1. These formulas are probably not in the best shape for a practical implementation. Observe that we can substitute (27b) into (51), (41) and hereby eliminate the auxiliary matrix ZðtÞ. In the same fashion, merge (23), (29) with (33), (40) to a direct time recursion for Q ðtÞ. These calculations are rather lengthy and elementary. They are therefore omitted. Table 2 shows the resulting ‘‘streamlined’’ variant of the algorithm. Also note that we computed the vector b much earlier in the algorithm. This makes it possible to bypass the explicit calculation of Q ðtÞ. Both updating followed by downdating transformations are directly applied to Q ðt 1Þ. Therefore, this matrix undergoes a transformation of the form ‘identity minus rank-two’. Next consider the resulting recursion for the square-root power matrix PðtÞ. We can see that this Table 1 Raw pseudocode of the sliding window SVD subspace tracker based on the unsymmetric Householder partial compressor.
Initialize :
8 Q ð0Þ:¼random orthonormal > > < Uð0Þ:¼random orthonormal > > : Pð0Þ:¼sI
FOR t ¼ 1; 2; 3; . . . 2
hðtÞ ¼ UT ðt 1ÞxðtÞ 6 6 XðtÞ ¼ xT ðtÞxðtÞ hT ðtÞhðtÞ 6 3 2 1=2 6 T h ðtÞ X ðtÞ 6 6 7 6 6 0 7 6 6 Pðt 1Þ ¼ 6 7 : build matrix 6 .. 7 6 6 4 . Pðt 1Þ 5 6 6 6 0 6 6 6 Pðt 1Þ ! us ðtÞ; vs ðtÞ : track minor singular vectors " # 6 6 u0 6 : partition 6 us ðtÞ ¼ uss 6 6 " # 6 6 v0 6 : partition 6 vs ðtÞ ¼ v ss 6 6 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 1 u0 1 6 6 j1 ¼ ; v1 ¼ uss 6 2 2j1 6 r ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 6 1 v0 1 6 j2 ¼ ; v2 ¼ vss 6 2 2j2 6 6 T 1=2 6 c ¼ 4j1 ðj2 X ðtÞ þ h ðtÞv2 Þ 6 6 6 ZðtÞ ¼ ðI 2v1 vT1 ÞPðt 1ÞðI 2v2 vT2 Þ þ v1 ðcvT2 2j1 hT ðtÞÞ 6 " # 6 2j1 vT1 W 6 6 Q ðtÞ ¼ : partition ¼ T T 6 Q ðt 1ÞðI 2v v Þ b 1 1 6 6 T 6 6B ¼ b b 6 6 c ¼ B1=2 b 6 6 1 6 r ¼ 1 pffiffiffiffiffiffiffiffiffiffiffi ffi 6 6 1B 6 pffiffiffiffiffiffiffiffiffiffiffiffi 6g ¼ 1 1 B 6 6 6 Q ðtÞ ¼ WðI rccT Þ 6 6 PðtÞ ¼ ðI gccT ÞZðtÞ 6 6 6 fðtÞ ¼ Uðt 1Þðv2 j2 X 1=2 ðtÞhðtÞÞ 4 UðtÞ ¼ Uðt 1Þ 2ðj2 X 1=2 xðtÞ þ fðtÞÞvT2
ð5aÞ ð5cÞ
ð12Þ
ð21Þ; ð22Þ ð17aÞ
ð17bÞ ð18aÞ ð18bÞ ð27aÞ ð27bÞ ð23Þ; ð29Þ ð31Þ ð32Þ ð50Þ ð53Þ ð33Þ; ð40Þ ð51Þ; ð41Þ ð26aÞ ð26bÞ
The solution is easily found to be pffiffiffiffiffiffiffiffiffiffiffiffi
g ¼ 1 1 B.
(53)
s is a small positive start-up constant. Equations numbered as they appear in the text.
ARTICLE IN PRESS 358
P. Strobach / Signal Processing 90 (2010) 352–362
Table 2 Streamlined pseudocode of the sliding window SVD subspace tracker based on the unsymmetric Householder partial compressor. 8 Q ð0Þ:¼random orthonormal > > < Uð0Þ:¼random orthonormal Initialize : > > : Pð0Þ:¼sI
of Pðt 1Þ as follows: 3 2 T h ðtÞ X 1=2 ðtÞ 7 6 7 6 0 7 6 7 6 Pðt 1Þ ¼ 6 7 .. 7 6 . Pðt 1Þ 5 4 0
FOR t ¼ 1; 2; 3; ::: 2
2
T
hðtÞ ¼ U ðt 1ÞxðtÞ 6 6 XðtÞ ¼ xT ðtÞxðtÞ hT ðtÞhðtÞ 6 2 3 6 T 6 X 1=2 ðtÞ h ðtÞ 6 6 7 6 6 0 7 6 7 6 Pðt 1Þ ¼ 6 6 6 Pðt 1Þ 7 4 5 6 6 0 6 6 6 Pðt 1Þ ! u ðtÞ; v ðtÞ : track minor singular vectors s 6 " #s " # 6 6 6 u ðtÞ ¼ u0 ; v ðtÞ ¼ v0 : partition s 6 s uss vss 6 6 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 6 1 u0 1 6 j1 ¼ ; v1 ¼ uss 6 2j1 2 6 r ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 6 1 v0 1 6 j2 ¼ ; v2 ¼ vss 6 2j2 2 6 6 T 6 b ¼ ½0 0 1Q ðt 1ÞðI 2v1 vT1 Þ 6 6 6 B ¼ bT b 6 6 c ¼ B1=2 b 6 6 1 6 6 r ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi 6 1B 6 p ffiffiffiffiffiffiffiffiffiffiffiffi 6 6g ¼ 1 1 B 6 6 c ¼ 4j ðj X 1=2 ðtÞ þ hT ðtÞv Þ 6 2 1 2 6 6 Z1 ¼ vT1 c 6 6 w ¼ c 2Z v 6 1 1 6 6 y ¼ PT ðt 1Þv 1 6 1 6 6 y2 ¼ Pðt 1Þv2 6 6 y ¼ PT ðt 1Þw 6 3 6 6 Z2 ¼ yT v2 1 6 6 6 Z3 ¼ yT3 v2 6 6 z1 ¼ 2y1 ð4Z þ cÞv2 þ 2j hðtÞ 2 1 6 6 6 z2 ¼ gy3 gð2Z3 Z1 cÞv2 2gZ1 j1 hðtÞ 6 6 PðtÞ ¼ Pðt 1Þ v1 zT czT 2y vT 2 2 1 2 6 " # 6 6 Q ðtÞ 2j1 ðvT1 rZ1 cT Þ 6 ¼ 6 T T v r wc Þ Q ðt 1ÞðI 2v 1 1 6 6 6 6 fðtÞ ¼ Uðt 1Þðv2 j2 X 1=2 ðtÞhðtÞÞ 4 UðtÞ ¼ Uðt 1Þ 2ðj2 X
1=2
xðtÞ þ fðtÞÞvT2
6 60 6 ¼6. 6 .. 4
ð5aÞ ð5cÞ
ð12Þ
ð21Þ; ð22Þ ð17a; bÞ ð18aÞ ð18bÞ
ð31Þ ð32Þ ð50Þ
1
0
32 1=2 X ðtÞ 76 76 0 76 76 6 .. 7 Lðt 1Þ 56 . 4 0 00
T
h ðtÞ
3
7 7 7 7. 7 Rðt 1Þ 7 5
(54)
These LR-factors of Pðt 1Þ are then used for computing the minor singular vectors according to (22). Now the overall subdominant operations count becomes quite apparent. We account approximately 6r 2 multiplications for the updating of the LR-factors of Pðt 1Þ. Assuming two iterations of (22), we account 4r 2 additional multiplications for the computation of the desired minor singular vectors. Finally 3r 2 multiplications are required for the computation of the three auxiliary y-vectors in the algorithm. This results in an overall subdominant operations count of approximately 13r2 multiplications per time update for this algorithm.
ð53Þ ð27aÞ
5. Simulation results The SVD subspace tracker of Table 2 has been tested extensively. We show results and comparisons with the sliding window Bi-SVD subspace tracker of [4]. As a test signal, we used the following vector sequence: xðtÞ ¼ sðtÞ þ nðtÞ,
(55)
where nðtÞ is an N-dimensional vector of zero-mean white Gaussian noise and ð26aÞ ð26bÞ
sT ðtÞ ¼ ½s1 ðtÞ; s2 ðtÞ; . . . ; sN ðtÞ
is a signal vector consisting of four independent cosine components as follows:
s is a small positive start-up constant. Equations numbered as they appear in the text.
matrix propagates from one time step to the next via a rank-three update. We brought it into this form because now, there exists the nice opportunity to work with an LRrepresentation of the P-matrix. LR-representations can be updated quite nicely. A fast Bennett-type algorithm for rank-one updating of an LR-decomposition is described in [11]. A Fortran-code of this algorithm is reproduced in [6]. This code comprises 12 lines of code and computes the rank-one update of an r r LR-decomposition with only 2r 2 multiplications and divisions. Hence we can update an LR-representation of the P-matrix with roughly 6r 2 multiplications. Once the P-matrix is available in LR-factorized form, we can easily establish the corresponding LR-factorization
(56)
sk ðtÞ ¼
pjm ðtÞ cos k þ om t ; 180 m¼1 4 X
k ¼ 1; 2; . . . ; N,
(57)
with o1 ¼ 0:1, o2 ¼ 0:3, o3 ¼ 0:2, and o4 ¼ 0:7. The time-varying spatial angles jm ðtÞ; m ¼ 1; 2; 3; 4 are given by the ‘jump scenario’ displayed in Fig. 1. See [4] for a detailed specification of this test signal. Results are shown for the following parameter settings: L ¼ 80, N ¼ 160, and a leakage factor of l ¼ 0:003 for the Bi-SVD tracker. All results shown here were obtained at a signal-to-noise ratio (SNR) of 0 dB. This means that we adjusted EfsT ðtÞsðtÞg ¼ EfnT ðtÞnðtÞg. Fig. 2 shows the first 10 channels of the vector signal sðtÞ, and Fig. 3 shows the corresponding 10 channels of the data sequence xðtÞ. We used an exact SVD of XðtÞ implemented in double precision as a reference. The algorithms were implemented in standard single precision real Fortran. Subspaces are compared in terms of the
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
Here Q ðtÞ is always a subspace produced by a tracker under test and VðtÞ is the corresponding reference subspace given in terms of the leading singular vectors of the reference SVD. As a second criterion, we monitored the orthonormality of the basis vector sets obtained by the trackers. The orthonormality criterion is
90 80 70 60 50 40 30
f QQ ðtÞ ¼ 10 log10 ðtrfFTQQ ðtÞFQQ ðtÞgÞ ½dB,
20
where
10
FQQ ðtÞ ¼ Q T ðtÞQ ðtÞ I.
0 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
10 signal channels
Fig. 1. Angle parameters j1 ðtÞ–j4 ðtÞ of the test signal ‘jump scenario’.
(60)
(61)
Again Q ðtÞ denotes any subspace basis, as obtained from the trackers under test. The results are shown in Figs. 4–10. Fig. 4 is the Q subspace error of the algorithm in Table 2 in terms of criterion eQQ ðtÞ. The corresponding result of the Bi-SVD tracker is the Q B -subspace error displayed in Fig. 5. Fig. 6 shows the orthonormality error of the tracker in Table 2 in terms of the criterion f QQ ðtÞ. The corresponding orthonormality error of the Q B -subspace of the Bi-SVD tracker
10
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Q-subspace error (dB)
angles (degrees)
359
5 0 -5 -10
Fig. 2. Ten channels s1 ðtÞ2s10 ðtÞ of the signal vector sðtÞ.
-15 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
10 data channels
Fig. 4. Q-subspace error according to criterion eQQ ðtÞ for the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs plotted in one diagram.
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time Fig. 3. Ten channels x1 ðtÞ2x10 ðtÞ of the data vector xðtÞ.
following criterion eQQ ðtÞ:
5 0 -5 -10 -15
eQQ ðtÞ ¼ 10 log10 ðtrfETQQ ðtÞEQQ ðtÞgÞ ½dB,
(58)
where T
QB-subspace error (dB)
10
T
EQQ ðtÞ ¼ VðtÞV ðtÞ Q ðtÞQ ðtÞ.
(59)
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 5. Q B -subspace error according to criterion eQQ ðtÞ for the Householder Bi-SVD subspace tracker of [4] for the case L ¼ 80; N ¼ 160 and a leakage factor of l ¼ 0:003. Ten independent trial runs.
ARTICLE IN PRESS
0 -20 -40 -60 -80 -100 -120 -140 -160 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 6. Q-subspace orthonormality error according to criterion f QQ ðtÞ for the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
is shown in Fig. 7. The same results are shown for the right singular subspaces. Fig. 8 is the U-subspace error of the algorithm in Table 2 and Fig. 9 is the corresponding Q A subspace error of the Bi-SVD subspace tracker. Fig. 10 shows the orthonormality error of the U-subspace or Q Asubspace estimates, respectively. These curves are practically identical for both trackers and therefore only a single curve is shown. The following observations can be made:
There are only minor differences that can be observed in the results produced by the two trackers under comparison. There exists no other tracker in the open literature that could reach the lower bound in dominant complexity while maintaining these results. This is simply because the theory how such trackers can be derived without heuristical simplifications was not known before the invention of these Householder methods. Recall that only the tracker of this paper (Table 2) based on the compressor theorem offers an exact rectangular window in temporal direction. The Bi-SVD tracker of [4], however, is essentially an exponential subspace tracker with truncated tail. There are certain applications in array processing, where an exact rectangular window is indispensable. In these cases, the Bi-SVD tracker cannot be applied. A comparison of the left singular subspace errors (Figs. 4 and 5) suggests a slightly worse result of the BiSVD tracker. The reason is that the left singular subspace of the Bi-SVD tracker is already notably curved by the exponential decay caused by the l ¼ 0:003 leakage. Also fairly identical are the corresponding orthonormality error tracks (Figs. 6 and 7) of the two algorithms under comparison. We could show (see the experimental results in [4]) that the Bi-SVD tracker runs unstable for a too small value of l. This will first result in an unbounded growth of the corresponding orthonormality error of the algorithm. Thus the leakage in a Bi-SVD subspace tracker must be chosen significantly
QB-subspace orthonormality error (dB)
P. Strobach / Signal Processing 90 (2010) 352–362
0 -20 -40 -60 -80 -100 -120 -140 -160 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 7. Q B -subspace orthonormality error according to criterion f QQ ðtÞ for the Householder Bi-SVD subspace tracker of [4] for the case L ¼ 80; N ¼ 160 and a leakage factor of l ¼ 0:003. Ten independent trial runs.
10 U-subspace error (dB)
Q-subspace orthonormality error (dB)
360
5 0 -5 -10 -15 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 8. U-subspace error according to criterion eQQ ðtÞ for the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
large. As a rule of thumb, we can say: the shorter the time window and the lower the SNR, the higher must be the value of l. Any orthonormality error below 100 dB can be regarded as ideal orthonormality and will prospectively meet all practical requirements. A double-precision implementation of the algorithms will significantly decrease the orthonormality error to values below 280 dB which is regarded totally unnecessary in practice. Other results of the algorithms, such as the subspace error, will not be affected by double precision arithmetic. This indicates that these algorithms are numerically robust, as expected from Householder-type computations on the data level, without resorting to correlations, as is the case here. The right singular subspace error produced by the BiSVD tracker (Fig. 9) looks slightly better than the comparable result of the tracker of Table 2 (Fig. 8). In this case, we should remember that the Bi-SVD can benefit from an additional exponential decay in the time window which reduces the effective window
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
QA-subspace error (dB)
10 5 0 -5 -10 -15 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
U-subspace orthonormality error (dB)
Fig. 9. Q A -subspace error according to criterion eQQ ðtÞ for the Householder Bi-SVD subspace tracker for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
0 -20 -40 -60 -80 -100 -120 -140 -160 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 10. U-subspace orthonormality error according to criterion f QQ ðtÞ of the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
length, hence resulting in a slightly lower subspace error for the right singular subspace. In all cases, we see that these trackers adapt extremely rapidly. There is no chance to achieve a comparable performance with a conventional growing window subspace tracker with exponential forgetting. In view of these results, sliding window subspace trackers will most likely become attractive in modern signal and array processing, because with the advent of the results in [4] and in this paper, the theory and the operation of such trackers is now completely understood. Fig. 10 finally shows the orthonormality error of the right singular subspace estimates. This orthonormality error is practically identical for the two trackers under comparison. But we can observe that the orthonormality error of the right singular subspace is again significantly lower than the orthonormality error or the left singular subspace. This has nothing to do with the particular dimensions L and N as chosen in this experiment. It is merely a consequence of the type of processing that is used to create these subspace
361
estimates. In the case of the right singular subspace, we have only a single Householder transformation applied to the subspace basis matrices in each time step. This is the minimum amount of computation that is necessary to track a subspace and creates the minimum amount of roundoff error. Therefore, we obtain the best possible orthonormality results here. This orthonormality is as good as the orthonormality obtained from a batch MGS QR-factorization in each time step. An increase in orthonormality error must be accepted for the left singular subspace bases because there, two Householder transformations are applied in the updating and downdating procedure. Moreover, the downdating recursion assumes that Q ðtÞ is perfectly orthonormal in each time step. But this is not the case in practice. This deviation from orthonormality of Q ðtÞ will carry over in orthonormality errors of the basis Q ðtÞ of the next time step. It is hence also the subtle nature of the downdating recursion that causes a relative increase in roundoff error in the left singular subspace basis. On the other hand, we observe that the subspace error is much smaller in the left singular subspace estimates (Figs. 4 and 5) when compared to the subspace errors of the right singular subspace (Figs. 8 and 9). The reason is again rooted in the more subtle updating/ downdating process of the left singular subspace tracking part of the algorithm. Not only the updating recursion, but also the downdating (reorthonormalization) transformation will decrease the subspace error. This is because the reorthonormalizing Householder reflector that is applied in the downdating step is mainly determined by signal components and not by noise. Therefore, this downdating reflector acts like a matched filter and further decreases the subspace error significantly. Of course, this pronounced property of the trackers will have consequences on future developments using these trackers, because algorithms will be designed so that they can benefit from the inherently higher accuracy of the left singular subspaces of such trackers. We conclude these tests with a little experiment showing the reconstruction capability of the trackers. A ‘‘toppinning’’ low-rank adaptive reconstruction filter is given by T s^ ðtÞ ¼ ½1 0 0Q ðtÞPðtÞUT ðtÞ.
(62)
Clearly, this is just an extraction of the top row of the approximant (2), which can be regarded as an estimate of the actual signal sðtÞ. Fig. 11 shows the first 10 channels of this reconstructed signal components or signal estimates and Fig. 12 shows the corresponding reconstruction error eðtÞ ¼ sðtÞ s^ ðtÞ. Fairly identical results are obtained for both the Bi-SVD tracker and the tracker of Table 2. However, it presents no difficulty to guess that we could significantly increase the reconstruction performance if multiple invariance techniques are used, because we have narrow-band signals here. It is clear that we have
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
10 reconstructed signal channels
362
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
10-channel reconstruction error
Fig. 11. Reconstructed signal s^ 1 ðtÞ2s^ 10 ðtÞ using the top pinning low-rank adaptive filter.
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 12. Corresponding reconstruction error tracks e1 ðtÞ2e10 ðtÞ for the top pinning low-rank adaptive filter.
sufficient experience in multiple invariance techniques [13,14]. The development of sliding window methods of this kind is currently in progress [15]. Note, however, that we can only use the subspace tracker of Table 2 for this purpose, because multiple invariance techniques in the temporal domain require that the window is exactly rectangular. 6. Conclusions With the recent advent of Householder methods in subspace tracking, we are now in a position that gives us a complete understanding of the solid theory required for a
competitive handling of the sliding window subspace tracking problem. Sliding window subspace trackers are of particular interest, because they allow a much faster tracking when opposed to classical exponential trackers, and also because they feature directly the adaptive implementation of multiple invariance techniques, which appear in direct connection to corresponding 2-D ESPRIT estimators, which in turn can be implemented very effectively in adaptive form now. This application area, however, is reserved for the Householder compressor theorem based tracker of this paper, because only this tracker can be operated without leakage, and therefore offers a perfectly rectangular time window. References [1] P. Strobach, Square Hankel SVD subspace tracking algorithms, Signal Processing 57 (1) (February 1997) 1–18. [2] T.M. Toolan, D.W. Tufts, Efficient and accurate rectangular window subspace tracking, in Proceedings of the IEEE Workshop on Sensor Array and Multichannel Processing, July 2006, pp. 12–14. [3] R. Badeau, G. Richard, B. David, Sliding window adaptive SVD algorithms, IEEE Trans. Signal Processing 52 (1) (January 2004) 1–10. [4] P. Strobach, The QS-Householder sliding window bIsvd subspace tracker, IEEE Transactions on Signal Processing (April 2009), accepted for publication. [5] P. Strobach, The fast recursive row-Householder subspace tracking algorithm, Signal Processing (2009), in press, doi:10.1016/j.sigpro. 2009.04.012. [6] P. Strobach, Square-root Householder subspace tracking, Numer. Math. 113 (2009) 89–121. [7] P. Strobach, The fast Householder Bi-SVD subspace tracking algorithm, Signal Processing 88 (11) (November 2008) 2651–2661. [8] P. Strobach, The Householder compressor theorem and its application in subspace tracking, Signal Processing 89 (5) (May 2009) 857–875. [9] C.M. Rader, A.O. Steinhardt, Hyperbolic Householder transformations, IEEE Trans. Acoust. Speech Signal Process. 34 (6) (December 1986) 1589–1602. [10] A.W. Bojanczyk, A.O. Steinhardt, Stabilized hyperbolic Householder transformations, in: Proceedings of the ICASSP-89, Glasgow, Scotland, 1989, pp. 1278–1281. [11] P. Stange, A. Griewank, M. Bollho¨fer, On the efficient update of rectangular LU-factorizations subject to low-rank modifications, Electron. Trans. Numer. Anal. 26 (2007) 161–177. [12] R. Roy, A. Paulraj, T. Kailath, ESPRIT—a subspace rotation approach to estimation of parameters of cisoids in noise, IEEE Trans. Acoust. Speech Signal Process. 34 (10) (October 1986) 1340–1342. [13] P. Strobach, Two-dimensional equirotational stack subspace fitting with an application to uniform rectangular arrays and ESPRIT, IEEE Trans. Signal Process. 48 (7) (July 2000) 1902–1914. [14] P. Strobach, Total least squares phased averaging and 3-D ESPRIT for joint azimuth-elevation-carrier estimation, IEEE Trans. Signal Process. 49 (1) (January 2001) 54–62. [15] P. Strobach, Frontiers in low-rank multichannel signal and array processing, AST-Consulting Inc., Internal Report (classified), December 2008. [16] G.H. Golub, C.F. Van Loan, Matrix Computations, second ed., John Hopkins University Press, Baltimore, MD, 1989.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Sliding window adaptive SVD using the unsymmetric householder partial compressor 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 6 January 2009 Received in revised form 14 April 2009 Accepted 6 July 2009 Available online 14 July 2009
A fast algorithm for tracking the rank-r SVD-approximant Q ðtÞPðtÞUT ðtÞ of a sliding window data matrix XðtÞ of dimension L N is introduced, where PðtÞ is a square-root power matrix of dimension r r with rominfL; Ng. This algorithm is based on the unsymmetric Householder partial compressor and uses a reorthonormalizing Householder transformation for downdating. The concept is numerically self-stabilizing and requires no leakage. The dominant complexity is 4Lr þ 3Nr multiplications per time update which is the lower bound in complexity for an algorithm of this kind. Applications occur in the area of adaptive array processing and other forms of adaptive processing in finite duration subspaces. A complete algorithm summary is provided. Computer simulations illustrate the operation of the algorithm. & 2009 Elsevier B.V. All rights reserved.
Keywords: Sliding window subspace tracking Householder Reorthonormalization Adaptive SVD
1. Introduction In this article, we present an algorithm for recursive rank-r SVD approximation of N-channel data matrices X of dimension L N, which are shifted in time according to the following sliding rectangular window rule: " # " # XðtÞ xT ðtÞ ¼ . (1) T x ðt LÞ Xðt 1Þ In each time step, an old data snapshot vector xðt LÞ leaves the window, while a new snapshot vector xðtÞ is appended on the top. This scheme appears frequently in adaptive multichannel/array processing. We present a sliding window subspace tracker that computes, in each time step, a rank-r approximant XðtÞ Q ðtÞPðtÞUT ðtÞ,
(2)
where Q ðtÞ is an orthonormal L r left singular subspace basis, UðtÞ is an orthonormal N r right singular subspace basis, and PðtÞ is an r r square-root power matrix that is 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.07.002
orthonormally or (in the complex case) unitary similar to the leading singular values in the SVD of XðtÞ. Note that PðtÞ is a dense unstructured square matrix. No shape restrictions are imposed on PðtÞ. Approximants of the kind (2) feature an adaptive implementation of various pertinent processing schemes for multichannel and array data. A special algorithm of this kind is the square Hankel SVD subspace tracker of [1]. This algorithm is restricted to time series subspace tracking. Solutions to the general rectangular sliding window subspace tracking problem appeared recently in [2,3]. However, these algorithms are redundant in complexity. A heuristically accelerated variant is proposed in [3]. In [4], we adopt the basic approach suggested in [3], and present a fast and exact solution using some essential mathematical methods known from fast Householder subspace tracking [5–7]. Additionally, we introduce a reorthonormalizing Householder transformation for downdating, a fundamentally new concept for downdating in finite duration QR-type algorithms. The resulting algorithm in [4] is already maximally fast and reaches the 4Lr þ 3Nr complexity lower bound for an algorithm of this kind. However, the method is numerically unstable because the updating recursion for the
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
internal square-root power matrix forms an undamped loop. Undamped loops occur frequently in sliding-window algorithms. They run unstable in floating-point implementations due to round-off error accumulation. A simple remedy is the incorporation of a small leakage in the recursion. This effectively stabilizes the algorithm. The incorporation of leakage has the same effect as exponential forgetting. The resulting window is no longer a true rectangular window, but becomes an exponential window with truncated tail. There are many applications, where this modification can be accepted. Sometimes, it is even desirable in the sense of a ‘‘trim’’ that assigns a larger weight to the incoming data snapshots. On the other hand, there are some pertinent applications in array processing, for instance in ESPRIT spectral estimation [12], where one has to work with overlapping sub-windows. The underlying assumption in such schemes is that the basic tracker window is perfectly rectangular. In such cases, leakagestabilized algorithms cannot be applied. Therefore, we introduce a completely self-stabilizing algorithm that preserves the perfectly rectangular shape of the temporal window. The principle behind this algorithm is an unsymmetric variant of the Householder partial compressor based on the Householder partial compressor theorem introduced in [8]. This scheme extends the estimated subspace by an additional dimension, in which the updating residual is accommodated. The algorithm uses this degree of freedom to wipe out all round-off noise from the subspace into this rank-one residual subspace. This way the rectangular shape of the window is preserved and stability can be maintained without the application of leakage. This paper is organized as follows: In Section 2, we introduce the unsymmetric Householder partial compressor and its application in the updating of the QPUapproximant. Section 3 shows how downdating can be accomplished using a reorthonormalizing Householder transformation. Section 4 summarizes the complete algorithm and presents a particularly streamlined version of the algorithm that is suitable for a recursive LRfactorized representation of the square-root power matrix PðtÞ. In this form, the algorithm has a subdominant complexity of Oðr 2 Þ. Section 5 shows some simulation results in order to illustrate the operation of the algorithm. Section 6 summarizes the conclusions. 2. Updating the QPU-approximant
353
augmented form: " Xðt 1Þ ¼
#
xT ðtÞ T
Q ðt 1ÞPðt 1ÞU ðt 1Þ
.
(4)
The orthonormal decomposition of the incoming snapshot vector xðtÞ with respect to the subspace spanned by Uðt 1Þ is required. The necessary formulas are given below (note this is the same procedure as in conventional subspace tracking with exponential forgetting): hðtÞ ¼ UT ðt 1ÞxðtÞ,
(5a)
x? ðtÞ ¼ xðtÞ Uðt 1ÞhðtÞ,
(5b) T
XðtÞ ¼ xT? ðtÞx? ðtÞ ¼ xT ðtÞxðtÞ h ðtÞhðtÞ, x? ðtÞ ¼ X
1=2
(5c) (5d)
ðtÞx? ðtÞ,
xðtÞ ¼ X 1=2 ðtÞx? ðtÞ þ Uðt 1ÞhðtÞ.
(5e)
Recall that (5a) is the obligatory input data compression step that is common to all fast subspace trackers. Eq. (5b) produces the complement of the orthogonal projection of xðtÞ onto the subspace spanned by the column-vectors of Uðt 1Þ. This constitutes the necessary decomposition of the input data vector into a component that can be represented by a linear combination of the information that is already available in Uðt 1Þ, and an innovation which is represented by x? ðtÞ. The squared norm of this innovation is computed in (5c). Formula (5e) finally represents the desired orthonormal decomposition of xðtÞ with respect to the subspace spanned by Uðt 1Þ. Now it is not difficult to verify that (4) can also be expressed in a factorized form as follows: 3 32 1=2 2 T 1 00 X ðtÞ h ðtÞ # 7" 76 60 7 76 0 xT? ðtÞ 6 7 7 6 Xðt 1Þ ¼ 6 . .. T 7 76 6 .. 4 . Q ðt 1Þ 54 . Pðt 1Þ 5 U ðt 1Þ 0
0 (6)
2.1. The unsymmetric Householder partial compressor We will now demonstrate that the above augmented ‘‘old’’ QPU-approximant (6) can be updated using two Householder matrices according to the Householder compressor principle developed in [8]. For this purpose, introduce the following two Householder matrices: " # 1 2j21 2j1 vT1 , (7) H1 ¼ 2j1 v1 I 2v1 vT1 with
In this section, we develop the necessary details required in the updating of a rank-r QPU-approximant using the unsymmetric Householder partial compressor. Recall the sliding window updating scheme (1). The rank-r QPU-approximant at time step t 1 is defined as follows (recall (2)): ^ 1Þ ¼ Q ðt 1ÞPðt 1ÞUT ðt 1Þ. Xðt
(3)
We wish to update all components of this decomposition in time upon the arrival of a new data vector xðtÞ at time step t. To accomplish this, consider the following
vT1 v1 þ j21 ¼ 1,
(8a)
H1 H1 ¼ I,
(8b)
and " H2 ¼
1 2j22
2j2 vT2
2j2 v2
I 2v2 vT2
# ,
(9)
with vT2 v2 þ j22 ¼ 1,
(10a)
H2 H2 ¼ I.
(10b)
ARTICLE IN PRESS 354
P. Strobach / Signal Processing 90 (2010) 352–362
The parameters of these Householder matrices can be determined so that the augmented QPU-approximant of the previous time step is updated as follows: 3 2 3 2 T 1 00 X 1=2 ðtÞ h ðtÞ 7 6 7 6 7 6 0 7 60 7 6 7 6 7 Xðt 1Þ ¼ 6 . 7H1 H1 6 7 6 . 6 .. Q ðt 1Þ 7 . 6 . Pðt 1Þ 7 5 4 5 4 0
2
H2 H2 4
3
xT? ðtÞ UT ðt 1Þ 2 1=2
6 6 6 ¼ ½qðtÞ Q ðtÞ6 6 6 4
lmin
0
5 00
3
7" 7 uT ðtÞ # 7 7 7 UT ðtÞ ZðtÞ 7 5
0 .. . 0
T
1=2
¼ Q ðtÞZðtÞU ðtÞ þ lmin qðtÞuT ðtÞ.
(11)
Observe that the two Householder matrices in this updating scheme must be adjusted so that the energy in the augmented P-matrix of the previous time step is perfectly mapped into the r r interior matrix ZðtÞ. The inevitable residual that remains after this compression 1=2 step is the minor singular value lmin of Pðt 1Þ in the upper-left corner of the augmented matrix after compression. With the appearance of this particular updating scheme, the discussion about optimality in fast subspace tracking becomes largely obsolete, because a subspace tracker based on this scheme will obviously attain the highest possible performance, since it is based on the highest attainable compression in each time step. Verify again (11) to see that the remaining residual layer l1=2 qðtÞuT ðtÞ, which is popped off in each time step, will min 1=2 always be of minimum energy, because the factor lmin is the minimum attainable residual in a problem of this kind, and the fact that tracefqðtÞuT ðtÞg ¼ 1 is always guaranteed as a consequence of the unconditional unitnorm property of the vectors qðtÞ and uðtÞ. All that remains is to determine the rules for computing the parameters of the two Householder compressor matrices that act in this scheme. Since this is just the unsymmetric variant of the Householder compressor that has been introduced already in [8], we can easily extend these results to the unsymmetric case. For this purpose, introduce the following compact notation: 3 2 1=2 T h ðtÞ X ðtÞ 7 6 7 6 0 7. (12) Pðt 1Þ ¼ 6 . 7 6 . 4 . Pðt 1Þ 5 0 Now the core compression problem of (11) can be expressed as follows: 3 2 1=2 lmin 0 0 7 6 7 6 0 7 ¼ H1 Pðt 1ÞH2 . 6 . (13) 7 6 . 4 . ZðtÞ 5 0
Introduce the SVD of Pðt 1Þ in a partitioned form 3 2 1=2 lmin 0 0 7" vT # 6 7 s 6 0 7 T , Pðt 1Þ ¼ ½us Us 6 7 6 .. 4 . K1=2 5 Vs
(14)
s
0
l1=2 min
where is the minor (smallest) singular value of Pðt 1Þ, us and vs are the corresponding left and right minor singular vectors, Ks1=2 is the diagonal matrix of the leading (larger) singular values, and Us and Vs are the corresponding left and right singular vectors. Now introduce two arbitrary orthonormal rotors Hu and Hv of dimension r r. Apparently, we can write 3 2 1=2 lmin 00 " # 7 6 T us 7 6 0 7. 6 (15) Pðt 1Þ½v V H ¼ s s v 7 6 .. T 1=2 HTu UTs 4 . Hu Ks Hv 5 0 The important insight at this point is the following: there always exists a pair of such rotors Hu and Hv so that the partially rotated singular vector matrices on the right and left sides of Pðt 1Þ in (15) become Householder matrices. Consequently, we can write 2 3 2 3 uTs 1 2j21 2j1 vT1 4 5 ¼ H1 ¼ 4 5, (16a) 2j1 v1 I 2v1 vT1 HTu UTs 2 3 1 2j22 2j2 vT2 4 5. ½vs Vs Hv ¼ H2 ¼ (16b) 2j2 v2 I 2v2 vT2 From these equations, it becomes apparent that these Householder matrices are already uniquely determined by the minor singular vectors us and vs . The Y-matrices are never computed explicitly. The conditional equations for the desired Householder matrix parameters are directly obtained from the top row of (16a) and the leftmost column of (16b), respectively: 3 " # 2 1 2j21 u0 4 5, ¼ (17a) us ¼ uss 2j1 v1 3 " # 2 1 2j22 v0 4 5. ¼ vs ¼ (17b) 2j2 v2 vss A quick calculation yields rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 u0 1 ; v1 ¼ j1 ¼ uss , 2 2j1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 v0 1 ; v2 ¼ j2 ¼ vss , 2 2j2
(18a) (18b)
and the necessary parameters of the Householder partial compressors are uniquely determined. As a by-product, we may observe that the SVD of ZðtÞ is given as follows: ZðtÞ ¼ HTu K1=2 s Hv ,
(19)
but this result is of no interest in this context. The minor singular vectors of Pðt 1Þ are tracked using a variant of bi-inverse iteration. Inverse iteration, in general, is a very important concept in the area of
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
eigenvalue and singular value computation [16]. Biinverse iteration comprises the following iterations (where we omitted the time index for the moment): for
m ¼ 1; 2; :::; m 1
am ¼ P
T
½uðtÞ UðtÞ ¼ ½x? ðtÞ Uðt 1ÞH2 .
(24)
UðtÞ ¼ 2j2 x? ðtÞvT2 þ Uðt 1ÞðI 2v2 vT2 Þ
vm
¼ 2j2 X 1=2 ðtÞ½xðtÞ Uðt 1ÞhðtÞvT2
T
um ¼ ðbm bm Þ1=2 bm endfor
Next consider the right singular subspace update. According to (11), we can write
Now we can exploit the structure of H2 as before and additionally, we eliminate the explicit calculation of x? ðtÞ by substituting (5d) together with (5b) into (24). This yields
um1
vm ¼ ðaTm am Þ1=2 am bm ¼ P
355
(20)
þ Uðt 1ÞðI 2v2 vT2 Þ
(25)
or more compactly: The practical tracker is based on an initial LR-(lower-left upper-right)-factorization of Pðt 1Þ: Pðt 1Þ ¼ LR.
(21)
In each time step, we initialize u0 ¼ us ðt 1Þ and iterate: for
m ¼ 1; 2; . . . ; m Ld ¼ um1 ! f orward substitution ! d Ram ¼ d ! back substitution ! am
fðtÞ ¼ Uðt 1Þðv2 j2 X 1=2 ðtÞhðtÞÞ,
(26a)
UðtÞ ¼ Uðt 1Þ 2ðj2 X 1=2 ðtÞxðtÞ þ fðtÞÞvT2 ,
(26b)
where we introduced an auxiliary vector fðtÞ. Finally, we arrive at the necessary update for the Zmatrix according to (11). We take into account the special structure of the two Householder matrices used in this update and obtain
c ¼ 4j1 ðj2 X 1=2 ðtÞ þ hT ðtÞv2 Þ,
vm ¼ ðaTm am Þ1=2 am T
ZðtÞ ¼ ðI 2v1 vT1 ÞPðt 1ÞðI 2v2 vT2 Þ
T
þ v1 ðcvT2 2j1 h ðtÞÞ,
R e ¼ vm ! f orward substitution ! e
T
L bm ¼ e ! back substitution ! bm T
um ¼ ðbm bm Þ1=2 bm endfor
(22)
Experiments revealed that m ¼ 2 iterations are usually sufficient for convergence. The desired minor singular vectors are obtained as us ðtÞ ¼ um and vs ðtÞ ¼ vm . At t ¼ 0, this minor singular vector tracker is initialized with a random unit-norm vector us ð0Þ. In Section 4, we show that the LR-factorization (21) must not be computed completely anew in each time step. We formulate a completely recursive variant of the algorithm, where the L- and R-factors are updated separately; thus leading to an overall subdominant operations count of only Oðr 2 Þ. Hence the final algorithm will also be optimal in terms of its subdominant operations count since no Oðr 3 Þ computations will be required anymore. 2.2. Final updating recursions Concluding this section, we summarize the resulting updating recursions for the two subspace estimates and the square-root power matrix. Recall the basic update equation (11) and the associated Householder matrices (7) and (9). We can see that qðtÞ is an uninteresting byproduct in the updating process and can be omitted. Taking into account the structure of Householder matrix (7), the following compact expression for Q ðtÞ is obtained: " # 2j1 vT1 Q ðtÞ ¼ . (23) Q ðt 1ÞðI 2v1 vT1 Þ
(27a)
(27b)
where c is again an auxiliary quantity. Hereby, the updating step is complete. The residual layer 1=2 lmin qðtÞuT ðtÞ is popped off the updated approximant (this means the residual layer is not computed explicitly) and we finally obtain the updated rank-r approximant in the following form: ^ 1Þ ¼ Q ðtÞZðtÞUT ðtÞ. Xðt
(28)
Recall that this updated approximant is of extended dimension ðL þ 1Þ N. 3. Downdating the QPU-approximant This section is devoted to the downdating part of the algorithm. Downdating is quite tricky, because it employs some unconventional ideas which are not familiar from ordinary subspace tracking with exponential forgetting. Moreover, the underlying idea is also not familiar from related downdating problems in sliding window recursive least-squares (RLS) adaptive filtering, where a hypernormal Householder transformation was used for downdating [9,10]. We propose a different concept based on a reorthonormalizing Householder transformation here. First of all, notice that only the Q Z-part in (28) requires downdating. The U-subspace requires no downdating. In the updating step, the vertical dimension of the Q subspace has grown by one. Now the new left singular subspace basis must be constructed from the upper L rows of Q ðtÞ. For this purpose, introduce the following partition: W (29) Q ðtÞ ¼ T , b
ARTICLE IN PRESS 356
P. Strobach / Signal Processing 90 (2010) 352–362 T
where b is the bottom row vector of Q ðtÞ. Introduce a generalized form of a Householder-type matrix as follows:
C ¼ I bbbT .
(30)
Define T
B¼b b
(31)
and introduce a scaled version of b as follows: c ¼ B1=2 b.
(32)
Now C can be expressed in true Householder style with a scaled unit-norm rank-one layer:
C ¼ I bbbT ¼ I bBccT ¼ I rccT .
(33)
Clearly, this matrix turns into a conventional Householder matrix for the special case of r ¼ 2. We will next demonstrate that the free parameter r in the Householder-type matrix (33) can be tuned so that the upper submatrix W in (29) is transformed into a matrix with perfectly orthonormal columns. As soon as this is achieved, the remaining bottom row in the updated problem can be dropped, and the sliding window downdating transformation is complete. Notice that this procedure is fundamentally different from classical sliding window recursive least-squares (RLS) downdating using hypernormal Householder transformations as in [9,10]. In contrast to our G-reflector (33), a hypernormal Householder matrix (such as used in [9]) has the following form: H ¼ U 2vvT ,
(34)
and hence our G-reflector of (33) cannot be regarded as a hypernormal Householder matrix in the classical sense anymore. We use a variable scalar parameter r instead. This parameter is tuned so that the submatrix W in (29) is transformed into a matrix with perfectly orthonormal columns. A preliminary study has revealed that we could develop alternative forms of sliding window RLS algorithms on this basis as well. Now let us see how we can benefit from the G-concept in the downdating of the Q Z-part of (28). Using (29) and (30), we can write W 1 Q ðtÞZðtÞUT ðtÞ ¼ ZðtÞUT ðtÞ. (39) T CC b We shall demonstrate that there exists a special choice for the parameter r so that the application of C from the right perfectly re-orthonormalizes W. In this case, downdating is accomplished with Q ðtÞ ¼ WC ! Q T ðtÞQ ðtÞ ¼ I. On the other hand, the P-matrix is constructed via PðtÞ ¼ C1 ZðtÞ,
(41)
and the downdating step is complete. All that remains is to determine the parameter r so that the orthonormality condition (40) is fulfilled. For this purpose, recall that the updated Q -matrix of the previous time step has orthonormal columns: T
Q ðtÞQ ðtÞ ¼ I.
(42)
Moreover, a look at partitioning scheme (29) reveals that we can write T
where U is a diagonal matrix that contains þ1’s and 1’s on the main diagonal in any arbitrary order. We must have at least one 1 on the main diagonal, otherwise H reduces to an ordinary Householder matrix. In any event, it is clear that such a U satisfies:
UU ¼ I and UUU ¼ U. If, for instance, we adjust 3 2 0 6 .. 7 6 I . 7 7, U¼6 7 6 0 5 4 0 0 1
(35)
T
Q ðtÞQ ðtÞ ¼ WT W þ bb ¼ I.
(43)
Consequently, T
WT W ¼ I bb .
(44)
Now since in (40) we assumed that Q ðtÞ ¼ WC with Q T ðtÞQ ðtÞ ¼ I, we can also write
CWT WC ¼ CðI bbT ÞC ¼ I.
(36)
(40)
(45)
Substitute the C of (30) into (45) to obtain the desired conditional equation for the b-parameter that converts the G-reflector into a perfect orthonormalizer for the W-matrix: T
T
T
ðI bbb ÞðI bb ÞðI bbb Þ ¼ I.
(46)
then we could use this concept to develop an algorithm that discards data from the tail of the sliding window, if we only constrain v so that
A little algebraic exercise reveals that the b’s which satisfy (46), also satisfy the following quadratic equation:
vT Uv ¼ 1,
(37)
b2 b þ
(38)
which has the solution 1 1 b1;2 ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi , B 1B
because then: HUH ¼ U.
This is the underlying idea in classical Householder sliding window RLS algorithms as advocated in [9]. On the other hand, our concept based on the Greflector is fundamentally different from the classical hypernormal approach in that we are working with U ¼ I
2 B
1 ¼ 0, BðB 1Þ
(47)
(48)
or equivalently (using (33)) 1
r1;2 ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi . 1B
(49)
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
In most Householder problems, it makes little difference which one of the two solutions is actually pursued. But it turned out that this is clearly not the case here. Although both the ‘þ’, as well as the ‘’ sign solution allow a downdating according to (40), there will still be a significant difference in the practical result obtained for the two cases. The ‘þ’ sign solution produces a Q ðtÞ with a much higher orthonormality error than the corresponding ‘’ sign solution. Computer experiments have revealed that the orthonormality error may grow by a significant amount in some cases when the ‘þ’ sign solution for the r-parameter is used. In order to understand the reasons for this effect, we should have a look at the special case of a vanishing T bottom row vector b ¼ ½0 0 in Q ðtÞ according to (29). Clearly, since Q ðtÞ is always bounded to orthonormal columns, it follows that for a vanishing bottom row vector W will also show perfectly orthonormal columns. Hence no downdating transformation will be required anymore because the downdating solution is directly given by Q ðtÞ ¼ W. This in turn means (recall (40)) that we can assume C ¼ I or r ¼ 0 in this particular case, because no reorthonormalization will be required anymore in this T case. Now since b ¼ ½0 0 causes B ¼ 0 according to (31), we can check in (49), which one of the two choices for the sign in (49) produces the desired r ¼ 0 for a given B ¼ 0. Clearly, we can see that this is perfectly the case for the ‘’ sign. But what happens in the case of the ‘þ’ sign solution? Again, we verify from (49) that in this case, we obtain r ¼ 2 and for this choice, C according to (33) turns into a conventional Householder matrix. Obviously, this is formally correct because the application of any ordinary Householder matrix to W will not destroy the already orthonormal character of W. Nevertheless, a completely unnecessary transformation is carried out. Experimental results revealed that this has a somewhat deteriorating effect on the numerical orthonormality of the resulting Q ðtÞ, because it generates more round-off error than the ‘‘null-transform’’ solution associated with the favourable ‘’ sign. Therefore, our choice is always: 1
r ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi . 1B
(50)
With this important conclusion, the downdating operation is almost complete. All that remains is to determine the corresponding inverse transform matrix that is required in the downdating of the Z-matrix according to (41). This is again an easy exercise, if one knows how to work with matrices of the type ‘identity minus rank-one’. In this case, it is clear that we can assume
C1 ¼ I gccT .
(51)
The necessary link between g and r is then determined by evaluating
CC1 ¼ ðI rccT ÞðI gccT Þ ¼ I.
(52)
357
For safety reasons, we check the special case B ¼ 0, which gives us g ¼ 0, the desired solution for an inactive reflector. Hereby the algorithm is complete. 4. Implementation A raw pseudocode of the algorithm is listed in Table 1. These formulas are probably not in the best shape for a practical implementation. Observe that we can substitute (27b) into (51), (41) and hereby eliminate the auxiliary matrix ZðtÞ. In the same fashion, merge (23), (29) with (33), (40) to a direct time recursion for Q ðtÞ. These calculations are rather lengthy and elementary. They are therefore omitted. Table 2 shows the resulting ‘‘streamlined’’ variant of the algorithm. Also note that we computed the vector b much earlier in the algorithm. This makes it possible to bypass the explicit calculation of Q ðtÞ. Both updating followed by downdating transformations are directly applied to Q ðt 1Þ. Therefore, this matrix undergoes a transformation of the form ‘identity minus rank-two’. Next consider the resulting recursion for the square-root power matrix PðtÞ. We can see that this Table 1 Raw pseudocode of the sliding window SVD subspace tracker based on the unsymmetric Householder partial compressor.
Initialize :
8 Q ð0Þ:¼random orthonormal > > < Uð0Þ:¼random orthonormal > > : Pð0Þ:¼sI
FOR t ¼ 1; 2; 3; . . . 2
hðtÞ ¼ UT ðt 1ÞxðtÞ 6 6 XðtÞ ¼ xT ðtÞxðtÞ hT ðtÞhðtÞ 6 3 2 1=2 6 T h ðtÞ X ðtÞ 6 6 7 6 6 0 7 6 6 Pðt 1Þ ¼ 6 7 : build matrix 6 .. 7 6 6 4 . Pðt 1Þ 5 6 6 6 0 6 6 6 Pðt 1Þ ! us ðtÞ; vs ðtÞ : track minor singular vectors " # 6 6 u0 6 : partition 6 us ðtÞ ¼ uss 6 6 " # 6 6 v0 6 : partition 6 vs ðtÞ ¼ v ss 6 6 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 1 u0 1 6 6 j1 ¼ ; v1 ¼ uss 6 2 2j1 6 r ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 6 1 v0 1 6 j2 ¼ ; v2 ¼ vss 6 2 2j2 6 6 T 1=2 6 c ¼ 4j1 ðj2 X ðtÞ þ h ðtÞv2 Þ 6 6 6 ZðtÞ ¼ ðI 2v1 vT1 ÞPðt 1ÞðI 2v2 vT2 Þ þ v1 ðcvT2 2j1 hT ðtÞÞ 6 " # 6 2j1 vT1 W 6 6 Q ðtÞ ¼ : partition ¼ T T 6 Q ðt 1ÞðI 2v v Þ b 1 1 6 6 T 6 6B ¼ b b 6 6 c ¼ B1=2 b 6 6 1 6 r ¼ 1 pffiffiffiffiffiffiffiffiffiffiffi ffi 6 6 1B 6 pffiffiffiffiffiffiffiffiffiffiffiffi 6g ¼ 1 1 B 6 6 6 Q ðtÞ ¼ WðI rccT Þ 6 6 PðtÞ ¼ ðI gccT ÞZðtÞ 6 6 6 fðtÞ ¼ Uðt 1Þðv2 j2 X 1=2 ðtÞhðtÞÞ 4 UðtÞ ¼ Uðt 1Þ 2ðj2 X 1=2 xðtÞ þ fðtÞÞvT2
ð5aÞ ð5cÞ
ð12Þ
ð21Þ; ð22Þ ð17aÞ
ð17bÞ ð18aÞ ð18bÞ ð27aÞ ð27bÞ ð23Þ; ð29Þ ð31Þ ð32Þ ð50Þ ð53Þ ð33Þ; ð40Þ ð51Þ; ð41Þ ð26aÞ ð26bÞ
The solution is easily found to be pffiffiffiffiffiffiffiffiffiffiffiffi
g ¼ 1 1 B.
(53)
s is a small positive start-up constant. Equations numbered as they appear in the text.
ARTICLE IN PRESS 358
P. Strobach / Signal Processing 90 (2010) 352–362
Table 2 Streamlined pseudocode of the sliding window SVD subspace tracker based on the unsymmetric Householder partial compressor. 8 Q ð0Þ:¼random orthonormal > > < Uð0Þ:¼random orthonormal Initialize : > > : Pð0Þ:¼sI
of Pðt 1Þ as follows: 3 2 T h ðtÞ X 1=2 ðtÞ 7 6 7 6 0 7 6 7 6 Pðt 1Þ ¼ 6 7 .. 7 6 . Pðt 1Þ 5 4 0
FOR t ¼ 1; 2; 3; ::: 2
2
T
hðtÞ ¼ U ðt 1ÞxðtÞ 6 6 XðtÞ ¼ xT ðtÞxðtÞ hT ðtÞhðtÞ 6 2 3 6 T 6 X 1=2 ðtÞ h ðtÞ 6 6 7 6 6 0 7 6 7 6 Pðt 1Þ ¼ 6 6 6 Pðt 1Þ 7 4 5 6 6 0 6 6 6 Pðt 1Þ ! u ðtÞ; v ðtÞ : track minor singular vectors s 6 " #s " # 6 6 6 u ðtÞ ¼ u0 ; v ðtÞ ¼ v0 : partition s 6 s uss vss 6 6 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 6 1 u0 1 6 j1 ¼ ; v1 ¼ uss 6 2j1 2 6 r ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 6 1 v0 1 6 j2 ¼ ; v2 ¼ vss 6 2j2 2 6 6 T 6 b ¼ ½0 0 1Q ðt 1ÞðI 2v1 vT1 Þ 6 6 6 B ¼ bT b 6 6 c ¼ B1=2 b 6 6 1 6 6 r ¼ 1 pffiffiffiffiffiffiffiffiffiffiffiffi 6 1B 6 p ffiffiffiffiffiffiffiffiffiffiffiffi 6 6g ¼ 1 1 B 6 6 c ¼ 4j ðj X 1=2 ðtÞ þ hT ðtÞv Þ 6 2 1 2 6 6 Z1 ¼ vT1 c 6 6 w ¼ c 2Z v 6 1 1 6 6 y ¼ PT ðt 1Þv 1 6 1 6 6 y2 ¼ Pðt 1Þv2 6 6 y ¼ PT ðt 1Þw 6 3 6 6 Z2 ¼ yT v2 1 6 6 6 Z3 ¼ yT3 v2 6 6 z1 ¼ 2y1 ð4Z þ cÞv2 þ 2j hðtÞ 2 1 6 6 6 z2 ¼ gy3 gð2Z3 Z1 cÞv2 2gZ1 j1 hðtÞ 6 6 PðtÞ ¼ Pðt 1Þ v1 zT czT 2y vT 2 2 1 2 6 " # 6 6 Q ðtÞ 2j1 ðvT1 rZ1 cT Þ 6 ¼ 6 T T v r wc Þ Q ðt 1ÞðI 2v 1 1 6 6 6 6 fðtÞ ¼ Uðt 1Þðv2 j2 X 1=2 ðtÞhðtÞÞ 4 UðtÞ ¼ Uðt 1Þ 2ðj2 X
1=2
xðtÞ þ fðtÞÞvT2
6 60 6 ¼6. 6 .. 4
ð5aÞ ð5cÞ
ð12Þ
ð21Þ; ð22Þ ð17a; bÞ ð18aÞ ð18bÞ
ð31Þ ð32Þ ð50Þ
1
0
32 1=2 X ðtÞ 76 76 0 76 76 6 .. 7 Lðt 1Þ 56 . 4 0 00
T
h ðtÞ
3
7 7 7 7. 7 Rðt 1Þ 7 5
(54)
These LR-factors of Pðt 1Þ are then used for computing the minor singular vectors according to (22). Now the overall subdominant operations count becomes quite apparent. We account approximately 6r 2 multiplications for the updating of the LR-factors of Pðt 1Þ. Assuming two iterations of (22), we account 4r 2 additional multiplications for the computation of the desired minor singular vectors. Finally 3r 2 multiplications are required for the computation of the three auxiliary y-vectors in the algorithm. This results in an overall subdominant operations count of approximately 13r2 multiplications per time update for this algorithm.
ð53Þ ð27aÞ
5. Simulation results The SVD subspace tracker of Table 2 has been tested extensively. We show results and comparisons with the sliding window Bi-SVD subspace tracker of [4]. As a test signal, we used the following vector sequence: xðtÞ ¼ sðtÞ þ nðtÞ,
(55)
where nðtÞ is an N-dimensional vector of zero-mean white Gaussian noise and ð26aÞ ð26bÞ
sT ðtÞ ¼ ½s1 ðtÞ; s2 ðtÞ; . . . ; sN ðtÞ
is a signal vector consisting of four independent cosine components as follows:
s is a small positive start-up constant. Equations numbered as they appear in the text.
matrix propagates from one time step to the next via a rank-three update. We brought it into this form because now, there exists the nice opportunity to work with an LRrepresentation of the P-matrix. LR-representations can be updated quite nicely. A fast Bennett-type algorithm for rank-one updating of an LR-decomposition is described in [11]. A Fortran-code of this algorithm is reproduced in [6]. This code comprises 12 lines of code and computes the rank-one update of an r r LR-decomposition with only 2r 2 multiplications and divisions. Hence we can update an LR-representation of the P-matrix with roughly 6r 2 multiplications. Once the P-matrix is available in LR-factorized form, we can easily establish the corresponding LR-factorization
(56)
sk ðtÞ ¼
pjm ðtÞ cos k þ om t ; 180 m¼1 4 X
k ¼ 1; 2; . . . ; N,
(57)
with o1 ¼ 0:1, o2 ¼ 0:3, o3 ¼ 0:2, and o4 ¼ 0:7. The time-varying spatial angles jm ðtÞ; m ¼ 1; 2; 3; 4 are given by the ‘jump scenario’ displayed in Fig. 1. See [4] for a detailed specification of this test signal. Results are shown for the following parameter settings: L ¼ 80, N ¼ 160, and a leakage factor of l ¼ 0:003 for the Bi-SVD tracker. All results shown here were obtained at a signal-to-noise ratio (SNR) of 0 dB. This means that we adjusted EfsT ðtÞsðtÞg ¼ EfnT ðtÞnðtÞg. Fig. 2 shows the first 10 channels of the vector signal sðtÞ, and Fig. 3 shows the corresponding 10 channels of the data sequence xðtÞ. We used an exact SVD of XðtÞ implemented in double precision as a reference. The algorithms were implemented in standard single precision real Fortran. Subspaces are compared in terms of the
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
Here Q ðtÞ is always a subspace produced by a tracker under test and VðtÞ is the corresponding reference subspace given in terms of the leading singular vectors of the reference SVD. As a second criterion, we monitored the orthonormality of the basis vector sets obtained by the trackers. The orthonormality criterion is
90 80 70 60 50 40 30
f QQ ðtÞ ¼ 10 log10 ðtrfFTQQ ðtÞFQQ ðtÞgÞ ½dB,
20
where
10
FQQ ðtÞ ¼ Q T ðtÞQ ðtÞ I.
0 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
10 signal channels
Fig. 1. Angle parameters j1 ðtÞ–j4 ðtÞ of the test signal ‘jump scenario’.
(60)
(61)
Again Q ðtÞ denotes any subspace basis, as obtained from the trackers under test. The results are shown in Figs. 4–10. Fig. 4 is the Q subspace error of the algorithm in Table 2 in terms of criterion eQQ ðtÞ. The corresponding result of the Bi-SVD tracker is the Q B -subspace error displayed in Fig. 5. Fig. 6 shows the orthonormality error of the tracker in Table 2 in terms of the criterion f QQ ðtÞ. The corresponding orthonormality error of the Q B -subspace of the Bi-SVD tracker
10
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Q-subspace error (dB)
angles (degrees)
359
5 0 -5 -10
Fig. 2. Ten channels s1 ðtÞ2s10 ðtÞ of the signal vector sðtÞ.
-15 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
10 data channels
Fig. 4. Q-subspace error according to criterion eQQ ðtÞ for the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs plotted in one diagram.
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time Fig. 3. Ten channels x1 ðtÞ2x10 ðtÞ of the data vector xðtÞ.
following criterion eQQ ðtÞ:
5 0 -5 -10 -15
eQQ ðtÞ ¼ 10 log10 ðtrfETQQ ðtÞEQQ ðtÞgÞ ½dB,
(58)
where T
QB-subspace error (dB)
10
T
EQQ ðtÞ ¼ VðtÞV ðtÞ Q ðtÞQ ðtÞ.
(59)
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 5. Q B -subspace error according to criterion eQQ ðtÞ for the Householder Bi-SVD subspace tracker of [4] for the case L ¼ 80; N ¼ 160 and a leakage factor of l ¼ 0:003. Ten independent trial runs.
ARTICLE IN PRESS
0 -20 -40 -60 -80 -100 -120 -140 -160 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 6. Q-subspace orthonormality error according to criterion f QQ ðtÞ for the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
is shown in Fig. 7. The same results are shown for the right singular subspaces. Fig. 8 is the U-subspace error of the algorithm in Table 2 and Fig. 9 is the corresponding Q A subspace error of the Bi-SVD subspace tracker. Fig. 10 shows the orthonormality error of the U-subspace or Q Asubspace estimates, respectively. These curves are practically identical for both trackers and therefore only a single curve is shown. The following observations can be made:
There are only minor differences that can be observed in the results produced by the two trackers under comparison. There exists no other tracker in the open literature that could reach the lower bound in dominant complexity while maintaining these results. This is simply because the theory how such trackers can be derived without heuristical simplifications was not known before the invention of these Householder methods. Recall that only the tracker of this paper (Table 2) based on the compressor theorem offers an exact rectangular window in temporal direction. The Bi-SVD tracker of [4], however, is essentially an exponential subspace tracker with truncated tail. There are certain applications in array processing, where an exact rectangular window is indispensable. In these cases, the Bi-SVD tracker cannot be applied. A comparison of the left singular subspace errors (Figs. 4 and 5) suggests a slightly worse result of the BiSVD tracker. The reason is that the left singular subspace of the Bi-SVD tracker is already notably curved by the exponential decay caused by the l ¼ 0:003 leakage. Also fairly identical are the corresponding orthonormality error tracks (Figs. 6 and 7) of the two algorithms under comparison. We could show (see the experimental results in [4]) that the Bi-SVD tracker runs unstable for a too small value of l. This will first result in an unbounded growth of the corresponding orthonormality error of the algorithm. Thus the leakage in a Bi-SVD subspace tracker must be chosen significantly
QB-subspace orthonormality error (dB)
P. Strobach / Signal Processing 90 (2010) 352–362
0 -20 -40 -60 -80 -100 -120 -140 -160 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 7. Q B -subspace orthonormality error according to criterion f QQ ðtÞ for the Householder Bi-SVD subspace tracker of [4] for the case L ¼ 80; N ¼ 160 and a leakage factor of l ¼ 0:003. Ten independent trial runs.
10 U-subspace error (dB)
Q-subspace orthonormality error (dB)
360
5 0 -5 -10 -15 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 8. U-subspace error according to criterion eQQ ðtÞ for the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
large. As a rule of thumb, we can say: the shorter the time window and the lower the SNR, the higher must be the value of l. Any orthonormality error below 100 dB can be regarded as ideal orthonormality and will prospectively meet all practical requirements. A double-precision implementation of the algorithms will significantly decrease the orthonormality error to values below 280 dB which is regarded totally unnecessary in practice. Other results of the algorithms, such as the subspace error, will not be affected by double precision arithmetic. This indicates that these algorithms are numerically robust, as expected from Householder-type computations on the data level, without resorting to correlations, as is the case here. The right singular subspace error produced by the BiSVD tracker (Fig. 9) looks slightly better than the comparable result of the tracker of Table 2 (Fig. 8). In this case, we should remember that the Bi-SVD can benefit from an additional exponential decay in the time window which reduces the effective window
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
QA-subspace error (dB)
10 5 0 -5 -10 -15 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
U-subspace orthonormality error (dB)
Fig. 9. Q A -subspace error according to criterion eQQ ðtÞ for the Householder Bi-SVD subspace tracker for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
0 -20 -40 -60 -80 -100 -120 -140 -160 0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 10. U-subspace orthonormality error according to criterion f QQ ðtÞ of the algorithm in Table 2 for the case L ¼ 80; N ¼ 160. Ten independent trial runs.
length, hence resulting in a slightly lower subspace error for the right singular subspace. In all cases, we see that these trackers adapt extremely rapidly. There is no chance to achieve a comparable performance with a conventional growing window subspace tracker with exponential forgetting. In view of these results, sliding window subspace trackers will most likely become attractive in modern signal and array processing, because with the advent of the results in [4] and in this paper, the theory and the operation of such trackers is now completely understood. Fig. 10 finally shows the orthonormality error of the right singular subspace estimates. This orthonormality error is practically identical for the two trackers under comparison. But we can observe that the orthonormality error of the right singular subspace is again significantly lower than the orthonormality error or the left singular subspace. This has nothing to do with the particular dimensions L and N as chosen in this experiment. It is merely a consequence of the type of processing that is used to create these subspace
361
estimates. In the case of the right singular subspace, we have only a single Householder transformation applied to the subspace basis matrices in each time step. This is the minimum amount of computation that is necessary to track a subspace and creates the minimum amount of roundoff error. Therefore, we obtain the best possible orthonormality results here. This orthonormality is as good as the orthonormality obtained from a batch MGS QR-factorization in each time step. An increase in orthonormality error must be accepted for the left singular subspace bases because there, two Householder transformations are applied in the updating and downdating procedure. Moreover, the downdating recursion assumes that Q ðtÞ is perfectly orthonormal in each time step. But this is not the case in practice. This deviation from orthonormality of Q ðtÞ will carry over in orthonormality errors of the basis Q ðtÞ of the next time step. It is hence also the subtle nature of the downdating recursion that causes a relative increase in roundoff error in the left singular subspace basis. On the other hand, we observe that the subspace error is much smaller in the left singular subspace estimates (Figs. 4 and 5) when compared to the subspace errors of the right singular subspace (Figs. 8 and 9). The reason is again rooted in the more subtle updating/ downdating process of the left singular subspace tracking part of the algorithm. Not only the updating recursion, but also the downdating (reorthonormalization) transformation will decrease the subspace error. This is because the reorthonormalizing Householder reflector that is applied in the downdating step is mainly determined by signal components and not by noise. Therefore, this downdating reflector acts like a matched filter and further decreases the subspace error significantly. Of course, this pronounced property of the trackers will have consequences on future developments using these trackers, because algorithms will be designed so that they can benefit from the inherently higher accuracy of the left singular subspaces of such trackers. We conclude these tests with a little experiment showing the reconstruction capability of the trackers. A ‘‘toppinning’’ low-rank adaptive reconstruction filter is given by T s^ ðtÞ ¼ ½1 0 0Q ðtÞPðtÞUT ðtÞ.
(62)
Clearly, this is just an extraction of the top row of the approximant (2), which can be regarded as an estimate of the actual signal sðtÞ. Fig. 11 shows the first 10 channels of this reconstructed signal components or signal estimates and Fig. 12 shows the corresponding reconstruction error eðtÞ ¼ sðtÞ s^ ðtÞ. Fairly identical results are obtained for both the Bi-SVD tracker and the tracker of Table 2. However, it presents no difficulty to guess that we could significantly increase the reconstruction performance if multiple invariance techniques are used, because we have narrow-band signals here. It is clear that we have
ARTICLE IN PRESS P. Strobach / Signal Processing 90 (2010) 352–362
10 reconstructed signal channels
362
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
10-channel reconstruction error
Fig. 11. Reconstructed signal s^ 1 ðtÞ2s^ 10 ðtÞ using the top pinning low-rank adaptive filter.
0
1000 2000 3000 4000 5000 6000 7000 8000 discrete time
Fig. 12. Corresponding reconstruction error tracks e1 ðtÞ2e10 ðtÞ for the top pinning low-rank adaptive filter.
sufficient experience in multiple invariance techniques [13,14]. The development of sliding window methods of this kind is currently in progress [15]. Note, however, that we can only use the subspace tracker of Table 2 for this purpose, because multiple invariance techniques in the temporal domain require that the window is exactly rectangular. 6. Conclusions With the recent advent of Householder methods in subspace tracking, we are now in a position that gives us a complete understanding of the solid theory required for a
competitive handling of the sliding window subspace tracking problem. Sliding window subspace trackers are of particular interest, because they allow a much faster tracking when opposed to classical exponential trackers, and also because they feature directly the adaptive implementation of multiple invariance techniques, which appear in direct connection to corresponding 2-D ESPRIT estimators, which in turn can be implemented very effectively in adaptive form now. This application area, however, is reserved for the Householder compressor theorem based tracker of this paper, because only this tracker can be operated without leakage, and therefore offers a perfectly rectangular time window. References [1] P. Strobach, Square Hankel SVD subspace tracking algorithms, Signal Processing 57 (1) (February 1997) 1–18. [2] T.M. Toolan, D.W. Tufts, Efficient and accurate rectangular window subspace tracking, in Proceedings of the IEEE Workshop on Sensor Array and Multichannel Processing, July 2006, pp. 12–14. [3] R. Badeau, G. Richard, B. David, Sliding window adaptive SVD algorithms, IEEE Trans. Signal Processing 52 (1) (January 2004) 1–10. [4] P. Strobach, The QS-Householder sliding window bIsvd subspace tracker, IEEE Transactions on Signal Processing (April 2009), accepted for publication. [5] P. Strobach, The fast recursive row-Householder subspace tracking algorithm, Signal Processing (2009), in press, doi:10.1016/j.sigpro. 2009.04.012. [6] P. Strobach, Square-root Householder subspace tracking, Numer. Math. 113 (2009) 89–121. [7] P. Strobach, The fast Householder Bi-SVD subspace tracking algorithm, Signal Processing 88 (11) (November 2008) 2651–2661. [8] P. Strobach, The Householder compressor theorem and its application in subspace tracking, Signal Processing 89 (5) (May 2009) 857–875. [9] C.M. Rader, A.O. Steinhardt, Hyperbolic Householder transformations, IEEE Trans. Acoust. Speech Signal Process. 34 (6) (December 1986) 1589–1602. [10] A.W. Bojanczyk, A.O. Steinhardt, Stabilized hyperbolic Householder transformations, in: Proceedings of the ICASSP-89, Glasgow, Scotland, 1989, pp. 1278–1281. [11] P. Stange, A. Griewank, M. Bollho¨fer, On the efficient update of rectangular LU-factorizations subject to low-rank modifications, Electron. Trans. Numer. Anal. 26 (2007) 161–177. [12] R. Roy, A. Paulraj, T. Kailath, ESPRIT—a subspace rotation approach to estimation of parameters of cisoids in noise, IEEE Trans. Acoust. Speech Signal Process. 34 (10) (October 1986) 1340–1342. [13] P. Strobach, Two-dimensional equirotational stack subspace fitting with an application to uniform rectangular arrays and ESPRIT, IEEE Trans. Signal Process. 48 (7) (July 2000) 1902–1914. [14] P. Strobach, Total least squares phased averaging and 3-D ESPRIT for joint azimuth-elevation-carrier estimation, IEEE Trans. Signal Process. 49 (1) (January 2001) 54–62. [15] P. Strobach, Frontiers in low-rank multichannel signal and array processing, AST-Consulting Inc., Internal Report (classified), December 2008. [16] G.H. Golub, C.F. Van Loan, Matrix Computations, second ed., John Hopkins University Press, Baltimore, MD, 1989.