Lecture 10: Recursive Least Squares Estimation

1 1. Lecture 10: Recursive Least Squares Estimation Overview Recursive Least Squares Estimation ;. The exponentially weighted Least Squares Recursive -in-time solution Initialization of the algorithm Recursion for MSE criterion Examples: Noise canceller, Channel equalization, Echo cancellation Reference : Chapter 9 from S. Haykin- Adaptive Filtering Theory - Prentice Hall, 2002. Lecture 10 2. Recursive Least Squares Estimation Problem statement Given the set of input samples {u(1), u(2), .. , u(N )} and the set of desired response {d(1), d(2), .. , d(N )}. In the family of linear filters computing their output according to M. X. y(n) = wk u(n k), n = 0, 1, 2, .. (1). k=0. Find recursively in time the parameters {w0 (n), w1 (n), .. , wM 1 (n)} such as to minimize the sum of error Squares n X n X M. X 1. E(n) = E(w0 (n), w1 (n), .. , wM 1 (n)) = (n, i)[e(i)2 ] = (n, i)[d(i) wk (n)u(i k)]2.

2 I=i1 i=i1 k=0. where the error signal is M. X 1. e(i) = d(i) y(i) = d(i) wk (n)u(i k). k=0. and the forgetting factor or weighting factor reduces the influence of old data 0 < (n, i) 1, i = 1, 2, .. , n usually taking the form (0 < < 1). (n, i) = n i , i = 1, 2, .. , n 2. Lecture 10 3. The exponentially weighted Least Squares solution Writing the criterion with an exponential forgetting factor n X n X M. X 1. E(n) = E(w0 (n), w1 (n), .. , wM 1 (n)) = n i [e(i)2 ] = n i [d(i) wk (n)u(i k)]2. i=i1 i=i1 k=0. Make the following variable changes: . u0 (i) = n i u(i); d0 (i) = n i d(i) (2). Then the criterion rewrites n X M. X 1 n X M. X 1. E(n) = n i [d(i) wk (n)u(i k)]2 = [d0 (i) wk (n)u0 (i k)]2. i=i1 k=0 i=i1 k=0. which is the standard LS criterion, in the new variables u0 (i), d0 (i). The LS solution can be obtained as n X n n n 0 0 T 1 X 0 0 X n i T 1 X.

3 W(n) = ( u (i)u (i) ) u (i)d (i) = ( u(i)u(i) ) n i u(i)d(i) = [ (n)] 1 (n) =. i=i1 i=i1 i=i1 i=i1. where we will denote (making use of Pre-windowing assumption, that data before i = 1 is zero). n X. (n) = n i u(i)u(i)T. i=1. n X. (n) = n i u(i)d(i). i=1. Lecture 10 4. Recursive in time solution We want to find a Recursive in time way to compute w(n) = [ (n)] 1 (n). using the information already available at time n 1, w(n 1) = [ (n 1)] 1 (n 1). We therefore will rewrite the variables (n) and (n) as functions of (n 1) and (n 1). n X n 1. X. (n) = n i u(i)u(i)T = n 1 i u(i)u(i)T + u(n)u(n)T = (n 1) + u(n)u(n)T. i=1 i=1. Xn n 1. X n 1 i (n) = n i u(i)d(i) = u(i)d(i) + u(n)d(n) = (n 1) + u(n)d(n). i=1 i=1. The matrix inversion formula If A and B are M M positive definite matrices, D is a N N matrix, and C is a M N matrix which are related by A = B 1 + CD 1 C T.

4 Then A 1 = B BC(D + C T BC) 1 C T B. Proof Exercise. Lecture 10 5. Derivation of the algorithm Applying the matrix inversion formula to (n) = (n 1) + u(n)u(n)T. we obtain 1 1 1 2 1 (n 1)u(n)uT (n) 1 (n 1). (n) = (n 1) . 1 + 1 uT (n) 1 (n 1)u(n). Denoting P (n) = 1 (n). and 1 P (n 1)u(n). k(n) = = Exercise = P (n)u(n). 1 + 1 uT (n)P (n 1)u(n). we obtain P (n) = 1 P (n 1) 1 k(n)uT (n)P (n 1). We are now able to derive the main time-update equation, that of w(n). w(n) = [ (n)] 1 (n) = P (n) (n) = P (n)( (n 1) + u(n)d(n)) = P (n)( (n 1)w(n 1) + u(n)d(n)). = P (n)(( (n) u(n)u(n)T )w(n 1) + u(n)d(n)) = w(n 1) P (n)u(n)u(n)T w(n 1) + P (n)u(n)d(n). = w(n 1) + P (n)u(n)(d(n) u(n)T w(n 1)) = w(n 1) + P (n)u(n) (n) = w(n 1) + k(n) (n). Lecture 10 6. where (n) = d(n) u(n)T w(n 1). is the innovation process (apriori errors). Now we can collect all necessary equations to form the RLS.

5 Algorithm: 1 P (n 1)u(n). k(n) =. 1 + 1 uT (n)P (n 1)u(n). (n) = d(n) u(n)T w(n 1). w(n) = w(n 1) + k(n) (n). P (n) = 1 P (n 1) 1 k(n)uT (n)P (n 1). Lecture 10 7. Initialization of RLS algorithm In RLS algorithm there are two variables involved in the recursions (those with time index n 1): w(n 1), Pn 1 . We must provide initial values for these variables in order to start the recursions : w(0). If we have some apriori information about the parameters w this information will be used to initialize the algorithm. Otherwise, the typical initialization is w(0) = 0. P0. 1. Recalling the significance of P (n). n X. P (n) = 1 (n) = [ n i u(i)u(i)T ] 1. i=i1. the exact initialization of the recursions uses a small initial segment of the data u(i1 ), u(i1 + 1), .. , u(0). to compute 0. X. P (0) = 1 (0) = [ i u(i)u(i)T ] 1. i=i1. However, it is not a simple matter to select the length of data required for ensuring invertibility of (0) !

6 Lecture 10 8. 2. The approximate initialization is commonly used, it doesn't require matrix inversion: P (0) = I. There is an intuitive explanation of this initialization. The significance P (n) = 1 (n) const. E(w(n) w)(w(n). T. w). can be proven. Thus, P (n) is proportional to the covariance matrix of the parameters w(n). Since our knowledge of these parameters at n = 0 is very vague, a very high covariance matrix of the parameters is to be expected, and thus we must assign a high value to . The recommended value for is > 100 u2. For large data length, the initial values assigned at n = 0 are not important, since they are forgotten due to exponential forgetting factor . Lecture 10 9. Summary of the RLS algorithm Given data u(1), u(2), u(3), .. , u(N ) and d(1), d(2), d(3), .. , d(N ). 1. Initialize w(0) = 0, P0 = I. 2. For each time instant, n = 1.

7 , N , Compute = uT (n)P (n 1) M 2 f lops = + u M f lops T. k(n) = M f lops (n) = d(n) wT (n 1)u(n) M f lops w(n) = w(n 1) + k(n) (n) M f lops 0. P = k(n) M 2 f lops P (n) = 1 (P (n 1) P 0 ) M 2 f lops We used flop as an abbreviation for one addition (subtraction) + one multiplication (floating point operations). The overall complexity of the algorithm is O(M 2 ) operations (flops) per time iteration. Lecture 10 10. Recursion for MSE criterion Suppose we have at moment n 1 the MSE criterion value n 1. X. E(w(n 1)) = n 1 i (d(i) w(n 1)T u(i))2. i=i1. and we want to know what is the MSE of the new filter w(n). n X n X. E(w(n)) = n i (d(i) w(n)T u(i))2 = (d(i))2 w(n)T (n). i=i1 i=i1. One remarkable recursion holds: E(w(n)) = E(w(n 1)) + (n)e(n). where (n) is the apriori error (n) = d(n) w(n 1)T u(n). and e(n) is the aposteriori error e(n) = d(n) w(n)T u(n).

8 Lecture 10 11. Applications of Recursive LS filtering 1. Adaptive noise canceller Single weight, dual-input adaptive noise canceller The filter order is M = 1 thus the filter output is y(n) = w(n)T u(n) = w(n)u(n). Denoting P 1 (n) = 2 (n), the Recursive Least Squares filtering algorithm can be rearranged as follows: RLS. Normalized LMS. Given data Given data u(1), u(2), u(3), .. , u(N ) and u(1), u(2), u(3), .. , u(N ) and d(1), d(2), d(3), .. , d(N ). d(1), d(2), d(3), .. , d(N ). 1. Initialize w(0) = 0, P0 = . 1. Initialize w(0) = 0. 2. For each time instant, n = 1, .. , N. 1 2. For each time instant, n = 1, .. , N. k(n) = u(n) 1. (n 1)2 + u(n)2 k(n) = u(n). (n) = d(n) w(n 1)u(n) a + u(n)2. (n) = d(n) w(n 1)u(n). w(n) = w(n 1) + (n)k(n). w(n) = w(n 1) + (n)k(n). 2 (n) = (n 1)2 + u(n)2. The normalized LMS algorithms can be obtained from RLS algorithm replacing the time varying term (n.)

9 1)2 with the constant a. Lecture 10 12. 2. Adaptive Equalization Modelling the communication channel We assume the impulse response of the channel in the form h i 1 1 + cos( 2 . h(n) = 2 W (n 2)) , n = 1, 2, 3. 0, otherwise The filter input signal will be 3. X. u(n) = (h a)(n) = h(k)a(n k) + v(n). k=1. where v2 = Selecting the filter structure The filter has M = 11 delays units (taps). y(n) = w(n)T u(n) = w0 (n)u(n) + w1 (n)u(n 1) + .. + w10 (n)u(n 10). The channel input is delayed 7 units to provide the desired response to the equalizer. d(n) = a(n 7). Lecture 10 13. Two Recursive (adaptive) filtering algorithms are compared: Recursive Least Squares (RLS) and (LMS). RLS algorithm has higher computational requirement than LMS , but behaves much better in terms of steady state MSE and transient time. For a picture of major differences between RLS and LMS, the main Recursive equation are rewritten: RLS algorithm 1.

10 Initialize w(0) = 0, P0 = I. 2. For each time instant, n = 1, .. , N. w(n) = w(n 1) + P (n)u(n)(d(n) wT (n 1)u(n)). P (n) = +u(n)T P1(n 1)u(n) (P (n 1) P (n 1)u(n)u(n)T P (n 1)). LMS algorithm 1. Initialize w(0) = 0. 2. For each time instant, n = 1, .. , N. w(n) = w(n 1) + u(n)(d(n) wT (n 1)u(n)). Lecture 10 14. Learning curve Ee2(n) for LMS and RLS algorithm 1. 10. LMS. RLS. 0. 10. 1. 10. LMS. 2. 10. RLS. 3. 10. 0 50 100 150 200 250 300 350 400 450 500. time step n The parameters: W = ( eigenvalue spread (R) = ). RLS: = 1, = 250. LMS: = Lecture 10 15. Echo Cancellation Echo in telephone networks (PSTN = public switched telephone network);. Subjectively echo is extremely annoying, the same as a too low volume or too high noise;. The echo is produced at the connection of four-wire circuits with two-wire circuits;. The customer loop (two-wire line, maximum 60 km) is connected to the central office by four-wire lines.