rL RL RL rT RR - IEEE Xplore

0 downloads 0 Views 135KB Size Report
r(0) rT. L. rL. RL. = RL rb;L. rT b;L r(0). (2). Manuscript received March 9, 2005; revised August 10, 2005. The associate editor coordinating the review of this ...
2362

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 6, JUNE 2006

Computation of the Condition Number of a Nonsingular Symmetric Toeplitz Matrix With the Levinson–Durbin Algorithm Jacob Benesty, Senior Member, IEEE, and Tomas Gänsler, Member, IEEE

Index Terms—Condition number, Levinson–Durbin, linear prediction, Toeplitz, Wiener–Hopf.

I. INTRODUCTION Consider a zero-mean stationary real signal x(n), where n = 0; 1; . . ., is the time index. The corresponding L 2 L covariance matrix is L

=

( ) ( ) r(0) r(1) r(1) r(0)

xL n x L n T

.. .

... ...

.. .

..

r(L 0 1) r(L 0 2)

.

...

r(L 0 1) r(L 0 2) .. . r(0)

where E f1g is the mathematical expectation, superscript transpose of a vector or a matrix, and

( ) = [ x (n )

xL n

x(n 0 1)

(1)

T

( ) = [ x (n)

x(n 0 1)

...

=E = =

( )

r(0) rL

RL rb;L

r(0)

+ 01 b

L 0 01 bTL

L

0 01 b

bT L

L01 L

L

:

L

(3)

where bL

= R01 r = [b 1 L

b;L

...

bL;2

L;

bL;L ]T

(4)

is the backward predictor of length L

= r(0) 0 r b = r(0) 0 r a is the prediction error energy, and a = J b

L

T b;L T L

L

(5)

L

L L L is the forward predictor with JL being the co-identity matrix. Equation (3) is important and will be used later for a fast computation of the condition number. The inverse of the covariance matrix plays a fundamental role in many important signal processing problems: Wiener–Hopf equations for identification and prediction, Yule–Walker equations in autoregressive (AR) modeling, etc, [2]–[4]. Therefore, it is essential to determine how this matrix is conditioned, especially when the data is perturbed (e.g., presence of noise). It is well known that the solution of a linear system, where RL+1 has to be inverted, is very sensitive to noise when the condition number is high [5]. This correspondence derives the condition numbers based on the Frobenius norm and gives an efficient way to compute them for a symmetric Toeplitz matrix. Some of the ideas used here are borrowed from [6] and [7]. This work is also an extension of the pioneering work by Cybenko [8].

Usually, the condition number is computed by using the 2-norm matrix. However, in the context of Toeplitz matrices, it is more convenient to use the Frobenius norm, as explained below. The covariance matrix RL+1 = R is symmetric, positive, and assumed to be nonsingular. It can be diagonalized as follows: T

Q RQ

=3

(6)

where T

( )

= QQ = I 3 = diag [1 ; 2 ; . . . ;  +1 ] T

Q Q

L

rT L

RL rT b;L

01 RL

01 =

RL+1

xL+1 n x L+1 n T

T

x(n 0 L) ]T ;

the covariance matrix of size (L + 1) 2 (L + 1) is RL+1

T

II. CONDITION NUMBER OF THE COVARIANCE MATRIX T

is a vector containing the L most recent samples of the signal x(n). Matrix RL is symmetric and Toeplitz [1]; we further assume that it is nonsingular. For a vector of length L + 1 xL+1 n

rb;L

= [ r(1) r(2) . . . r(L) ] ; = [ r(L) r(L 0 1) . . . r(1) ]

denotes

x(n 0 L + 1) ]

...

rL

By using the Schur complements, it is easy to invert RL+1 , as follows:

Abstract—One well-known and widely used concept in signal processing is the optimal Wiener filtering, where a linear system (Wiener–Hopf equations) has to be solved. The symmetric Toeplitz matrix that naturally appears in this system is the covariance matrix. If this matrix is ill-conditioned and the data is perturbed, the accuracy of the solution will suffer a lot if the linear system is solved directly. One way to improve the accuracy is to regularize the covariance matrix. However, this regularization depends on the condition number: the higher the condition number, the larger the regularization. Therefore, it is important to be able to estimate this condition number in an efficient way, in order to use this information for improving the quality of the solution. Many other problems require the knowledge of this condition number for different reasons. Therefore, it is of great interest to find a practical algorithm to determine this condition number, which is the focus of this correspondence.

R =E

where

(2)

and 0 < 1 is

 2  1 1 1   R

Manuscript received March 9, 2005; revised August 10, 2005. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Ta-Hsin Li. J. Benesty is with the Université du Québec, INRS-EMT, Montréal, QC H5A 1K6, Canada (e-mail: [email protected]). T. Gänsler is with MH Acoustics, Summit, NJ 07901 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TSP.2006.873494

(7) (8)

+1 . By definition, the square root of R

L

1=2 = Q31=2 QT :

(9)

The condition number of a matrix R is [5]

[R] = kRkkR01 k

(10)

where k1k can be any matrix norm. Note that  [R] depends on the underlying norm and subscripts will be used to distinguish the different condition numbers.

1053-587X/$20.00 © 2006 IEEE

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 6, JUNE 2006

2363

Consider the Frobenius norm

RF=

k

Moreover, by using the two inequalities

RT R]

ftr[

k

=

1 2

g

:

L

(11)

l=0 L

We have

R=

1 2

R]

= ftr[

F

L+1 =

l=1

g

=

1 2

=

l=0

1 2

l

(12)

1 2

R0 ]

L+1 =

Hence, the condition number of

R=

F

1 2

R=

R=

g

=

1 2

=

F

(13)

R



L + 1:

(14)

R

R

R

R

1 2

1

(L + 1)2



R

R=



1 2

2 [R] 

c2 



1 2

 1 [R ]   2 [R ] :

1

L+1

L+1 : 1 0 1 R  1=1 and tr[R]  L+1 , we have 2

Since tr

R=

1 2

R]tr[R0 ] 1

tr[

=



R]

tr[

1

22

obtain

R] tr

tr [

R01



(L + 1)

R]

tr[

1



R=

1 2



(L + 1)2



(L + 1)=1 ,

(L + 1)

1 2

L+1 1

(20)

:

R=

1 2



F

R=

1 2



(L + 1)2

R=

1 2

:

(22)

R] : (25)

2

1 2

1 2

In this section, we need to efficiently compute the two norms

RL=

1 2 +1

2

and forward. Indeed

F

RL0 =

2

1 2 +1

F

. The calculation of the first one is straight-

RL=

1 2 +1

2

F

RL

= tr [

+1 ]

= (L + 1)r (0):

(26)

Expression (26) requires one multiplication only. Consider the ma01 where its diagonal elements are gL+1;ii , i = trix L+1 = L +1 1; 2; . . . ; L + 1. It is clear from (3) that the last diagonal component 01 . The Lth diagonal element of of L+1 is gL+1;(L+1)(L+1) = L 0 1 0 1 2 is g =

+

b L+1 L+1;LL L01 L L;L . Continuing the same process, we easily find

R

L

gL+1;ii = i0011 +

l=i

2

l01 bl;i

(27)

with 0 = r(0). Therefore, from (27), we deduce that

RL0 =

1 2 +1

2

F

aa

bb

GL

= tr [

L

l=0 L l=0

+1 ]

l01

1+

bTl bl

l01

1+

alT al

(28)

with 0T 0 = T0 0 = 0. Finally, the condition number is

F2

RL=

1 2 +1

= (L + 1)r (0)

L l=0

l01

1+

alT al

:

(29)

By using the Toeplitz structure, the Levinson–Durbin algorithm 01 L , in O(L2 ) opersolves the linear prediction equation L = L ations, instead of O(L3 ). This algorithm computes all predictors l , l = 1; 2; . . . ; L, and this is exactly what we need to compute (29). Expression (29) shows also a very nice link between the condition number and the predictors of all orders. The proposed algorithm, which has roughly the same complexity as the Levinson–Durbin algorithm, is summarized in Table I. Note that a very efficient algorithm

a

(21)

(L +1)F [



III. FAST COMPUTATION OF THE CONDITION NUMBER

we

Therefore, we deduce that

2

1 2

R

=

2

R=

R=

2 but F [ ] 6= F

1 2

(18)

thus

F

R=

=

(19)

(24)

1 2

hence

F [R1=2 ]  2 [R1=2 ]: Also, since tr[R]  (L + 1)L+1 and tr R01

l=0

1 2

2

(17)

L+1 1



2l

1 2

(16)

We now show the same principle for the F- and 2-norm matrices and for 1=2 . We recall that

R

L

R = . According to expressions (22) and (25), F R = and F R = are a good measure of the condition number of matrices R = and R, respectively. Basically, there is no difference in the trend of the condition numbers of R and R = . In other words, if R = is ill-conditioned (respectively, well-conditioned), so is R. In Section III, we will show how to compute F R = by using the Levinson–Durbin algorithm.

G R= : (15) G For example, for the 1- and 2-norm matrices and for R, we can show G [5] R=

(23)

(L + 1)

F [R]  F2



R

c1 



1 2

(The inequality in the previous expression is easy to show by using the Cauchy–Schwartz inequality.) In this correspondence, we choose to 1=2 (rather than F [ ]), because as it will be shown in work on F Section III it is easy to derive efficient algorithms to estimate its value. As far as we know, it does not seem obvious to estimate efficiently F [ ] . 1=2 is large, then 1=2 is said to be an ill-conditioned maIf  trix. Note that this is a norm-dependent property. However, according to 1=2 1=2 [5], any two condition numbers  and  are equivalent in that constants c1 and c2 can be found for which

R

l=0

2

l

2l

1 2

:

associated with k1kF is 1 2

1 2

Note that 2 [ ] =

1 2

1

R0 =

F

R=

2 L+1 F

l=1 l

1 2

1 2

=

1

= ftr[

F



> 0, 8 l, it is easy to show that

1

and

R0 =

where l

L

2

l

R r

a

2364

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 6, JUNE 2006

R

R

(L + 1)F [ L+1 ], for different values of L + 1. F [ L+1 ] can be easily calculated from MATLAB with the function “cond(R,’fro’).” As 2 1=2 and F [ L+1 ] we can see, the two condition numbers F L+1 have exactly the same trend, which confirms that the first one can be used as a good measure of the condition of L+1 . Moreover, Fig. 1 1=2 shows that expression (25) is always verified. Also, F L+1 can be computed with MATLAB, using the function “cond(R0.5, ’fro’),” which we did. The two different ways to compute it give virtually the same results. The only difference is that our approach is much faster.

TABLE I COMPUTATION OF THE CONDITION NUMBER WITH THE LEVINSON–DURBIN ALGORITHM

R R

R

R

V. CONCLUSION Being able to estimate the condition number of the covariance matrix of a random signal x(n), with a reasonable amount of operations, can be very useful in many applications. In this correspondence, we have shown on how to compute, in an efficient way, the condition number of the square root of a symmetric Toeplitz matrix with the Frobenius norm. We have shown that the square of this condition number can be used as a good measure of the condition of L+1 . We have also presented a nice link between the condition number and the predictors of all orders. Simulations confirm that the proposed algorithm works well.

R

ACKNOWLEDGMENT The authors would like to express their sincere thanks to the four anonymous reviewers whose positive and constructive comments helped to enhance the quality of this correspondence. REFERENCES

Fig. 1. Condition number in decibels as a function of L + 1 (from 20 to 1000). (a)  [ ], (b)  , and (c) (L + 1) [ ].

R

R

R

TR

01 g was recently proposed by Dias and Leitão [9] to compute trf (where is a Toeplitz matrix, this form is a much general form than the one used in this section) with the Trench algorithm. Using these techniques here, we can reduce even more the complexity [to O(L ln L)] for the estimation of the overall algorithm.

T

IV. SIMULATIONS In this section, we show on one example the behavior of the proposed condition number. The signal x(n) is a 5-s (zero-mean) noisy speech signal sampled at 8 kHz. The correlation coefficients r(l) of the symmetric Toeplitz matrix L+1 were estimated, using the whole sequence of the speech, with the MATLAB function “xcorr.” The condition numbers were calculated for different sizes of the covariance matrix (from 20 2 20 to 1000 2 1000, every 20l). We implemented 1=2 the condition number 2F L+1 with the proposed method (see Table I). Fig. 1 compares this condition number with F [ L+1 ] and

R

R

R

[1] R. M. Gray, “Toeplitz and Circulant matrices: A review,” Stanford Univ., Stanford, CA, Internal Rep., 2002. [2] M. G. Bellanger, Adaptive Digital Filters and Signal Analysis. New York: Marcel Dekker, 1987. [3] B. Widrow and S. D. Stearns, Adaptive Signal Processing. Englewood Cliffs, N.J: Prentice-Hall, 1985. [4] S. Haykin, Adaptive Filter Theory, 4th ed. Upper Saddle River, N.J: Prentice-Hall, 2002. [5] G. H. Golub and C. F. Van Loan, Matrix Computations. Baltimore, MD: The Johns Hopkins Univ. Press, 1996. [6] J. Benesty and T. Gänsler, “New insights into the RLS algorithm,” EURASIP J. Appl. Signal Process., vol. 2004, pp. 331–339, Mar. 2004. , “A recursive estimation of the condition number in the RLS al[7] gorithm,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing (ICASSP), vol. 4, Mar. 2005, pp. iv/25–iv/28. [8] G. Cybenko, “The numerical stability of the Levinson–Durbin algorithm for Toeplitz systems and equations,” SIAM J. Sci. Statistical Comp., vol. 1, pp. 303–310, 1980. g [9] J. M. B. Dias and J. M. N. Leitão, “Efficient computation of trf for Toeplitz matrices,” IEEE Signal Process. Lett., vol. 9, no. 2, pp. 54–56, Feb. 2002.

TR